Dynamical systems forced by shot noise as a new paradigm in the interest rate modeling

SFB 649 Discussion Paper 2010-037

Dynamical systems

forced by shot noise as a new paradigm in the interest rate modeling

Alexander L. Baranovski*

* WestLB AG

This research was supported by the Deutsche

Forschungsgemeinschaft through the SFB 649 "Economic Risk".

http://sfb649.wiwi.hu-berlin.de ISSN 1860-5664

S FB

6 4 9

E C O N O M I C

R I S K

B E R L I N

Dynamical systems forced by shot noise as a new paradigm in the interest rate modeling

Alexander L. Baranovski (WestLB AG)

Abstract. In this paper we give a generalized model of the interest rates term structure including Nelson-Siegel and Svensson structure. For that we introduce a continuous m- factor exponential-polynomial form of forward interest rates and demonstrate its considerably better performance in a fitting of the zero-coupon curves in comparison with the well known Nelson-Siegel and Svensson ones. In the sequel we transform the model into a dynamic model for interest rates by designing a switching dynamical system of the considerably reduced dimension n < m generating the forward rate curves in form a càdlàg function. A system is described by n-th order linear differential equation driven by a stochastic or chaotic shot noise. From fitted forward rates we specify the parameters of the switching system and discuss perspectives of our models to produce term-structure forecasts at both short and long horizons.

Keywords: forward interest rates, shot noise processes, switching dynamical systems, chaotic Brownian subordination, chaotic maps

JEL classification: C13, C20 and C22

Disclaimer: The ideas presented below reflect the personal view of the author and are not necessarily identical to the official methodology used at WestLB AG

Acknowledgements: The partial financial support from the Deutsche Forschungsgemeinschaft via SFB 649 “Ökonomisches Risiko” is gratefully acknowledged

1. Exponential-polynomial models of interest rates

We introduce an exponential-polynomial term structure model of interest rates by

∑

⁻

=

−

⋅

= +

=

1

1 0

0( , ), ( , ) ( ) exp( )

) , (

l

i

i i

l f t Z f t Z t t

c Z t

f ϕ γ (1)

where z c c_p^l _l Z R^m

l ∈ ⊂

={ ,..., ⁻₋ , ,..., ₋} : ₁^{(1) ( 1)} _{1 1}

1 γ γ , t > 0 denotes time to maturity T,

∑

=

= − pi

k k i k

i t c t

1 1 )

) (

ϕ ( are polynomials of degree p_i −1 with coefficients c_k⁽ⁱ⁾ ; k =1,...,p_i (p_i ∈

{

1,2,3,...

}

) and γ_i are positive real numbers. The total number of parameters in (1) is m p l

l

i i +

=

∑

⁻

= 1

1

.

We note that the widely used Nelson-Siegel (N-S) [Nelson 1987] and Svensson (SV) [Svensson 1994] families of f(t) can be easily derived from (1). Assuming l =2 and

1 =2

p , p₂ =1, i.e. dimensionality m=4, (1) leads to N-S forward rate curve

f_NS(t)=(c₁⁽¹⁾+c₂⁽¹⁾t)⋅exp(−γ₁t)+c₂ (2) as well as the SV curve is given by

3 2 )

2 ( 1 1 )

1 ( 2 ) 1 (

1 ) exp( ) exp( )

( )

(t c c t t c t c

f_SV = + ⋅ −γ + ⋅ −γ + (3)

with l=3, p₁ =2 and p₂ =1, p₃ =1 i.e. m = 6.

Another special case of (1) is a curve of the exponentials mixture underp_i =1,∀i, i.e.

( ) ∑

⁻

( )

=

+

−

⋅

=

1

1

exp

l

i

l i i

EXP t c t c

f γ (4)

We show that a performance of a new term structure (1) in a fitting of the yields is considerably higher of the well known Nelson-Siegel and Svensson ones. By other words a today’s choice of the parameters in (1) is to be not limited by the state space

R4

Z ⊂ as for the N-S curve or Z ⊂ R⁶ in the Svensson model, i.e. Z ⊂Rⁿ,n≥6. The corresponding term structure of the bond prices will be then given by

( )

_^









−

=

∫

^{t f x Z dx}

Z t B

0

, exp

: ) ,

( at Z ⊂Rⁿ

2. ODE for the interest rate models

We establish that dynamics of the interest rates f (t) in a model (1) follows a n-th order ODE

where 1

1

1

+

=

∑

⁻

= l

i

pi

n , p_i is the multiplicity of the root γ_i, i =1,...,l of the corresponding characteristic polynomial

∑

⁻

=

−

+ −

=

1

1

) (

n

i

i n i n

D γ γn β γ (6)

The coefficients of the ODE are given by the Vieta formula

∑

<

<

<

− = − ⋅ ⋅

n

i i i

i i i j j

n

j

... j

~

~

~

2 1

2

1 ...

) 1

( γ γ γ

β (7)

with γ₁^~ =γ₂^~ =...=γ^~_pi =γ₁, ^~ ₁ ^~ ₂ ^~ ₂

2 1 1

1 γ ... γ γ

γ_p+ = _p+ = = _p+p = and so on.

Example (Nelson-Siegel).

Recall that the Nelson-Siegel model corresponds to the case l =2 and p₁ =2. Then 3

1 +1=

= p

n , γ₁^~ =γ₂^~ =γ₁ and γ₃^~ =γ₂ =0. It follows that

( )

( )

0 )

1 ( : ) 3 (

) 1 ( : ) 2 (

2 )

1 ( : ) 1 (

~ 3

~ 2

~ 1 3 0

2 1

~ 3

~ 2

~ 3

~ 1

~ 2

~ 1 2 1

1

~ 3

~ 2

~ 1 1 2

=

−

=

=

= +

+

−

=

=

−

= + +

−

=

=

γ γ γ β

γ γ γ γ γ γ γ β

γ γ

γ γ β

j j j

Thus, the corresponding differential equation is 0 2 ₁ ′′+ ₁² ′=

′′−

′ f f

f γ γ .

It is easy to check that the Nelson-Siegel curve f(t)=c₁+(c₂ +c₃t)⋅exp(−γ₁t)

is a general solution of the ODE as a combination of the two obvious particular solutions of r′=0 and r′′−2γ₁r′+γ₁²r =0, such that

2 1 1

1

) 0 ( ) 0 ( ) 2 0

( γ γ

f f f

c ′′

′ + +

= , ₂

1 1

2

) 0 ( ) 0 2 (

γ γ

f

c f ′′

′ −

−

= ,

1 3

) 0 ) (

0

( γ

f f

c ′′

′ −

−

= . (8)

Example (Svensson).

By analogy to the previous example the SV curve follows a 4^th order ODE in the form

( )^IV

(

2 1 2

) (

1² 2 1 2

)

1² 2 0 f − γ +γ f′′′+ γ + γ γ f′′−γ γ f′=

Generally ODE (5) can be presented in a matrix form

( )

^{, 0}

d t T

d tf=F γ ⋅f ≤ ≤ (9)

where f =





















−1 1 0

fn

f f

M ^,





















= γl

γ γ M

2 1

γ , F =

( ) ( ) ( )

0 1 1

0 1 0 0

0 0 1 0 0

0 0 1 0

0 0 1

β β βn₋

 

 

 

 

 

 

 

 

− − − 

 γ γ γ 

L L

M O O O M

L L L L L

L L L

and

( )

t t f d f d _i

i

i ≡ .

Eq. (9) has a vector solution ^f

( )

^{t =eF( )γ⋅t⋅f}

( )

^{0 .}

Define by

( )

tτ f

(

t Z_τ

)

r , := , (10)

the instantaneous forward rate at time t for date τ and by

(

τ

)

r

( )

t τ f

(

Z_τ

)

r 0, limt , 0,

0 =

= → (11)

the short rates, where stochastic process Z_τ with values in Z ⊂R^m contains two groups of processes ^{c1⁽¹⁾

( )

τ ^,...,c⁽_p^l_l⁻₋¹₁⁾

( )

τ ^} and {γ₁

( )

τ ,...,γ_l−1

( )

τ } such that the processes

( )

τ

j

ci depend on the stochastic processes {γ₁

( )

τ ,...,γ_l−1

( )

τ } and the short rates (initial conditions) as shown by (8) for a Nelson-Siegel term structure.

Taking into account (9), (10) the evolution of the forward rates on the date τ is given by

d

( )

d tr= F γ^τ ⋅r (12)

Assuming Z_τis given τ – evolution of r (t, τ) can be described by two different ways.

First one presents τ – evolution of r (t, τ) as a train of curves (12) for dates τ = 1, 2,…

in a form (see Fig. 1)

( )

i

( )

^{, 0}

i

d t i T t

d t^r=^{F γτ} ⋅ +^r

∑

^ξδ − ⋅ ≤ ≤ ∞ or equivalently

( ) ( )

d t T

d tr=F γ^τ ⋅ +r ξ^τδ − ⋅τ , (13) where +1



= T

τ t is a counting process ( t denotes the floor function ( largest integer smaller or equal to t ) and

(

0, 1

) (

T,

)

τ = τ+ − τ

ξ f Z f Z . (14)

We call the vector ξ_i as the stochastic amplitude of the impulse perturbation, which acts on system (13) at times t = τT, τ = 1,2,3, …. such that

(

τ ⋅T +0,τ

)

=r

(

τ⋅T −0,τ

)

+ξτ

r .

Thus r is a cadlag function, i.e. a right continuous function ^r

(

^{τ⋅T +0,τ}

)

^≡r

(

^τ⋅T,τ

)

^,

defined on Rⁿ and has a left limit.

Between kicks iT and (i +1)T a state vector is governed by the homogeneous system of linear differential equations (12) at τ = i+1.

Fig. 1 illustrates the system (13) with a jump ξ1 ≡ f

(

0,Z2

)

− f

(

T,Z1

)

at t = T

0 5 10 15 20 t

0.01 0.02 0.03 0.04 0.05 0.06 rHtL

Fig.1 Train of the two first curves with T = 10 years.

The solution to (13) is explicitly given by

( )( ( ) )

( )

_



 



 ×



 



 ⋅∑

⋅

= ₋

= =

−

−

⋅

∑

1

1

1 exp _i

i j i j

T t

F T F

e

r γ ξ

τ τ

τ γ_τ

where we denote ξ₀ ≡r

( )

0 .

The second approach is based on a simple idea to express a random variable f

(

^t,Zτ

)

by f

(

t^,Z_τfixed

)

+ζ

^{( )}

t in the interval 0≤t≤T . By other words we assume that

dynamics of the interest rates r (t, τ) can be modelled by

, 0 ,

dr= F r⋅ dt+dζ ≤ ≤t T ∀τ (15) where F is a n x n matrix with the constant coefficients ^{βi ≡βi}

(

^γfixed

)

for any τ . The model (15) generates a predicted term structure, whose exponential-polynomial shape depends on the model parameters and the initial short rate. One can show that (15) is more general and includes a class of equilibrium models such as Vasicek, CIR, lognormal models.

3. Approach I – demonstrating example: from estimating the yield curve to its dynamic modelling

First we empirically estimate processZ_τ. For that we are going to fit the default-free yield spreads, downloaded from the Reuters database. The observable time period is 23.02. 2006 – 14.01. 2008, i.e. contains q = 478 dates. In framework of the above exponential-polynomial approach we introduce a term structure model of yields curves by

( )

∫

=

t

dx Z x t f Z t Y

0

1 , : ) ,

( (16)

(16) can be easily done in a closed form. For this we need the relation [Prudnikov 1981]

( )

^!₁^{, 0,1,2,...,}

0

= +

−

= ₊

∫

^{t xke}−^{γxdx Mk t eγt} _γ^{kk k where}

^{( )}

0

₍

^!

₎

_! 1^, <⁰

= − ₊

−

∑

= ^γ

γⁱ

i k k

k i

t i k t k

M .

Substituting the model (1) into (16) leads to

( )

( )

∑

( )

∑

₌ −

−

= 



 



 −

+

− +

= _l ^l_{i k}^pi _kⁱ _k k _k t M t c

c Z t

Y 1 ¹_{1 1 1} 1!

) ,

( γ ⁽¹⁷⁾

We introduce a minimization criterion as:

( ) ( ( ) ( ) )

i Z

i

i Y t Z

N1 s , 2 min

→

−

=

∑

τ τ

τ

ρ (18)

where

{

s_i

( )

τ , i=1,2,...,N_τ

}

are quotes of the yields on the base dateτ . The cost function ρ

( )

τ is to be minimized by the appropriate choice of the m p l

l

i i +

=

∑

⁻

= 1

1

parameters of the state space Z. It is clear that in the fitting problem the following restriction

( )

Z p l

( {

N q

} )

m

l

i

i min , 1,..., dim

1

1

=

≤ +

=

=

∑

⁻

=

τ τ (19)

is to be provided.

Given the number l of the parameters γ_l, i=1,...,l and distribution of their multiplicitiesp_l,i=1,...,l such that the condition (19) holds the nonlinear regression technique for the least squares criterion (18) leads to the minimum of the cost function with the optimal parameters c c_p^l _l Z R^m

l⁻₋ , ,..., ₋}∈ ⊂ ,...,

{ ₁^{(1) ( 1)} _{1 1}

1 γ γ .

We extend the criterion (7) by

( )

pl

p l q

q ₁ ^,min^1,...,

1 →

=

∑

τ=

τ ρ

ρ (20)

with the obvious restrictions for the parameters

{

1,2,3,...

}

, ,

, p₁ p ∈Z⁺ =

l K _l . (21)

We note that (20) according to a law of large numbers/ Birkhoff ergodic theorem approaches the mean of the stochastic/ chaotic cost function withq→∞.

The problem (18)-(21) for euro swap rates has the following twofold solution

1) l =5, p_l =1, p_i =3,i=1,...,l−1 and { ₁⁽¹⁾,..., ^{( 1)}, ₁,..., ₁} ¹⁷

1 Z R

c

c _p^l _l

l⁻₋ γ γ ₋ ∈ ⊂ (m =17,

τ =60,∀τ

N ) for a “laminar” period of the observed yields quotes on the bond market:

(

355

)

07 . 07 . 17 ) 1 ( 06 . 02 .

23 τ = − τ =

2) l =5, p_l =1, p_l−1=1,p_i =3,i=1,...,l−2 and { ₁⁽¹⁾,..., ^{( 1)}, ₁,..., ₁} ¹⁵

1 Z R

c

c _p^l _l

l⁻₋ γ γ ₋ ∈ ⊂

(m=15) for a “turbulence” period ^18.07.07

(

τ =³⁵⁶

)

−^14.01.08(τ =⁴⁷⁸⁾ .

The above calculated m+l+1 parameters specify a general term structure model of interest rates by the exponential-quadratic curves (1) ( p_l =3 ) as well as a general term structure model of yields by the exponential-cubic curves (17).

Fig. 2 demonstrates a performance of the general model for the Euro swap rates by comparison of its cost function (18) with the cost functions of the conventional Nelson-Siegel and Svensson models and exponential model (4) with l = 6.

0 100 200 300 400 tdays

0.005 0.010 0.015 0.020 0.025 cost functionrHtL

Fig. 2 Comparison of the cost functions for the observable time period: 23.02. 2006 – 14.01. 2008. N-S is a green curve, SV is a black, the exponential is a blue, and the general model is a red curve.

Moreover the ratios 00074 9 . 0

00653394 .

0 ≈

=

GEN NS

ρ

ρ , 6

00074 . 0

00419495 .

0 ≈

=

GEN SV

ρ

ρ and 5

00074 . 0

0039557 .

0 ≈

=

GEN EXP

ρ

ρ (22)

quantify the performance of the general model with the derived optimal exponential-

with l = 6, respectively.

The mean of the cost function (20) for the exponential model (4) has a local minimum at l = 6 in value ρ_EXP =0.0039 and at l =5 in value ρ_GEN =0.00074 for a general model as shown in Fig. 3.

0 1 2 3 4 5 6 l-1

0.001 0.002 0.003 0.004 0.005 0.006 0.007

rHlL

Fig. 3. The means of the cost functions: blue curve – exponential and red one – general model

Fig. 4 collects all base curves fitted to the available data on the date 25.03.08 .

0%

1%

2%

3%

4%

5%

6%

7%

8%

9%

1 34 67 100 133 166 199 232 265 298 331 364 397 430 463 496 529 562 595

Maturity in months

%

GOVCHF GOVEUR GOVGBP GOVJPY GOVUSD GOVRUB SWAPCHF SWAPEUR SWAPGBP SWAPJPY SWAPUSD

Fig. 4 Fitted zero curves on the spot date 25.03.08

Repeating a fitting procedure to the another date, let’s say τ +1 = 26.03.08, we get a similar set of the base curves (16) which can be described by the system (13), where the impulse perturbation ξ_τ (14) is to be predetermined.

As an example we design a dynamical system for N-S instantaneous forward rate curve

( ) ( ( ) ( ) )

exp( ) )

,

(t c₁ c₂ c₃ t t

r τ = τ + τ + τ ⋅ −γ_τ described by the following low-dimensional ODE

(

t T

)

r r

r′′′−2γ_τ ′′+γ_τ² ′=ξ_τδ −τ⋅ (23) where

(

1

)

₂

(

1

)

₁

( )

( ₂

( )

₃

( )

) exp( )

1 c c c c T T

c _τ

τ τ τ τ τ τ γ

ξ = + + + − − + ⋅ − . (24)

according to (10) and (14).

From the output of the above fitting procedure we retrieve the time series

( ) ( ) ( )

{

c1τ ^,c2 τ ^,c3 τ ^,γ_ττ =^1,2,K^,q=⁴⁷⁸

}

for the observable time period.

Applying (24) we immediately get time series of stochastic perturbation.

Let us introduce the k-th order increments for both processes ξ_τ ,γ_τ by

τ τ

τ ξ ξ

ξ ^{k k}

k =∆ −∆

∆⁺¹ ₊₁ (25)

τ τ

τ γ γ

γ ^{k k}

k =∆ −∆

∆⁺¹ ₊₁ , (26)

where ∆⁰ ≡1, k = 0, 1, 2, ….

We are now able to do an elementary statistical analysis of both processes ξ_τ ,γ_τ.

Table 1 contains the histograms of the processes and Table 2 collects the histograms of the increments.

Histogram of γ_τ Histogram of ξ_τ

2 4 6 8 10

0 20 40 60 80 100 120 140

-0.015 -0.01 -0.005 0

0 10 20 30 40 50 60

Table 1.

Histograms of the increments ∆^kγ_τ, k =1,10 Histograms of the increments ∆^kξ_τ,k =1,10

-4 -2 0 2 4

0 50 100 150

-0.004 -0.002 0 0.002 0.004

0 50 100 150

-400 -200 0 200 400

0 20 40 60 80 100 120

-1 0 1 2

0 50 100 150

Table 2.

We note that the both processes are diffusion processes characterized by symmetrical bell-like but not-Gaussian distributions of their increments. A χ² Pearson’s test with the confidence level 0.95 rejects a hypothesis of the independence of the increments for both forced signal

ξτ and γ_τ.

The means of ∆^kγ_τand ∆^kξ_τ,∀k ≥1 are zeros. The sample variances σ∆²k_γ,σ∆²k_ξ grow exponentially with the order k as shown in Fig. in a log scale

10 20 30 40 k

-10 0 20 30 40 50 ln@s2_{_Dg}&xHkLD

Fig.5. Red line islnσ∆²k_ξ, black line is lnσ∆²k_γ

A distance between the above two lines does not remain constant with k but grows slowly with rate ln ₂ 0.00594878

2

∞

→

∆

∆ →

k k

k

dk d

ξ γ

σ

σ . It means that a variance of the ∆^kγ_τ grows a bit quicker

than a variance of ∆^kξ_τ,∀k≥1. The sample autocorrelation functions

( ) ( )

,...

2 , 1 , 0 ,

1 2 1

=

∆

− ∆

= −

∑

^{− −}

=

+

−

∆

∆ j

j k j q

C

j k q

j k

k k

k

τ

τ τ γ

γ σ γ γ

(27)

( ) ( )

,..

2 , 1 , 0 ,

1 2 1

=

∆

− ∆

= −

∑

^{− −}

=

+

−

∆

∆ j

j k j q

C

j k q

j k

k k

k

τ

τ τ ξ

ξ σ ξ ξ

(28)

are presented in Table 3.

Autocorrelations of ∆^kγ_τ, k =1,10, 20 Autocorrelations of ∆^kξ_τ, k =1,10, 20

10 20 30 40 t

-0.5 0.5 1.0 acf_xHtL

10 20 30 40 t

-1.0 -0.5 0.5 1.0 acf_xHtL

Table 3.

To estimate a mutual correlation of ξ_τ and γ_τ and their increments we introduce the Pearson product-moment correlation coefficient [Norman L. Johnson 1995]

(

^,

)

^{; 0,1,...; 0,1,...}

1 1 1

=

=

∆

⋅

− ∆

= −

∑

^{− −}

=

+

−

∆

−

∆ j k

j k k q

j R

j k q

j k

k k

k

τ

τ τ ξ

γσ ξ γ

σ (29)

depicted in Figures 6 – 8.

Fig. 6 . Geometric interpretation of the matrix R of dimension 40x40

-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8

1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 51 53 55 57 59 61

R(j,0) R(j,1) R(j,10) R(j,20)

Fig. 7 The Pearson’s correlation coefficient at j = 0,1,…,60; k = 0, 1, 10, 20

-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8

1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 51 53 55 57 59 61

R(0,k) R(1,k) R(10,k) R(20,k)

Fig. 8 The Pearson’s correlation coefficient at j = 0,1,10,20; k = 0, 1, …, 60

4. Approach II – forcing signal as a shot noise

We note that applying Vieta formula (7) at j = n the coefficient β₀ is equal to zero since 0

l =

γ . It implies that a general solution of the ODE (5) is a combination of the obvious solution of f′=0 and the particular solution η(t) of

0 ... ₁

) 3 ( 2 ) 2 ( 1 ) 1

( ⁻ + ₋ f ⁻ + ₋ f ⁻ + + f =

f ⁿ β_n ⁿ β_n ⁿ β (30)

i.e. f = ≡ +r η c.

We are specifically interested in the behaviour of the system (30) with a shot noise as a chaotic/stochastic perturbation, i.e.

( )

( )

^{, 0}

N t

i i

i

d t t t T

d t^η=^F⋅ +^η

∑

^Aδ − ≤ ≤ ⁽³¹⁾

0 1

2 n

η η η ₋

 

 

 

= 

 

 

η M ^,

1 2 1

0 1 0 0

0 0 0

0 0 1 0

0 0 1

n n

β β ₋ β ₋

 

 

 

 

= 

 

 

 

− − − 

 

F

L

M O O O L L

L L L L L

.

where tk are successive occurrence or arrival times of δ-impulses,

( )

0 0 1 ... k ...

t = <t < <t < <T

{

^{k t t}

}

t

N( )=max : _k ≤ is a counting process.

The impulse perturbation acts on system (31) at times t = tk, k = 1,2,3, …. such that

(

tk +⁰

)

=

(

tk −⁰

)

+ k

η η A _{. (32)}

In sequel we assume that η is a càdlàg function.

We introduce the positive inter-arrival times Tk such that

∑

=

− + =

=

k

i i k

k

k t T T

t

1

1

Between kicks a state vector is governed by the following homogeneous system of linear

differential equations

d

dtη= ⋅F η

and the initial condition of the system η0 ≡η

(

t0+0

)

defines an evolution of a state vector (Cauchy theorem).

4.1. One-dimensional case

We consider a special case of (31) in a form of one-dimensional ODE:

( )

( )

^{, 0}

N t

i i

i

d A t t t T

d tη=

∑

δ − ≤ ≤ ⁽³³⁾

4.1.1. Response to a Shot Noise. The generalized Wiener process

Integrating Eq. (33) we immediately get a solution

(34)

A plot of the process η is depicted in Fig.

ηk

ηk+1

η

t_k t_k+1 t T_k

Fig. 9. A solution of Eq. (1)

Thus the process η(t) (Fig. 9) is a rectangular signal with step heights ηk satisfying the relation:

1

k k Ak

η

=

η

₋ + ₍₃₅₎

Probability density function

( )

^{( )} _{( )}

( )

1

(0) 0

N t

k N t

k

t A

η η η

=

=

∑

≡ =

We first establish that the distribution function of η(t) is

( ( ) ) ( ) ( )

1

k k

k

P η t z P t P η z

∞

=

≤ =

∑

≤ ,

where P tk

( )

=P t

(

k₋1 ≤ <t tk

)

=P N t

( ( )

=k

)

.

We assume zero mean for the magnitudes Ak . It immediately implies that η(t) is a martingale with a zero mean

( ) ( ( ) ) ( )

⁰

(

⁰

) (

^|

^{( )} ) ( ) ( )

E η =E N t E A = A= as E η N t =N t E A (36)

To calculate a variance of η (t) we use a law of total variance

( ) (

^|

( ) ) (

^|

( ) )

D η =E D^ η N t ^+D E^ η N t ^ One can show that the conditional variance is

( ( ) ) ^{( ) ( )}

^{( ) 1}

^{( )}

1

| 2 ( ) ( )

N t

A n

D η N t N t D A N t n c n

−

=

= +

∑

− (37)

where c kA^{( )}=E A A

(

i i k₊

)

is the autocorrelation function (acf) and σA² ≡^D

( )

^A is the variance of the Ak.

Hence, for i.i.d. random or uncorrelated chaotic magnitudes Ak (cA(k)=0) we have

( )

_{( ) 0}

( ) ( ( ) ) ( ( ) )

²

( ) ( )

A t

c n

D η D A E N t D N t E A D A N

= = + = (38)

where

( ) { ( ) } { }

( )

t k

t t

k k k

N H t kP N t k P t t t

E t

>>

≡ =

∑

= =

∑

< is the intensity function.

For correlated random/chaotic A_k we assume that the first moment of the autocorrelation function cA(k) is finite

( )

1 A k

k c k

∞

=

< ∞

∑

(39)

and then the variance is given by

( )

2 k 2 o 1

σ = σ + . (40)

Let us introduce a new variable

k

k k

k

η

η η

ε σ

= − with E

( )

εk =⁰ and E

( )

εk² =¹ . It can be shown that ε_k converges in distribution to the standard normal law, i.e. the central limit theorem holds both with i.i.d.random [Feller] and chaotic magnitudes Ak [Chernov 1995]. In [Baranovski 2003], authors have presented the analytical expressions for the characteristic functions of the chaotic partial sums η_kof the magnitudes A_k generated by PWL onto maps and shown their fast convergence to the limit ^exp

(

^{−ω2/ 2}

)

^.

We consider a piecewise constant function Wk

( )

t on ^t∈

[ ]

^0,1 such that

(41)

where   x is the floor function (it gives the greatest integer less than or equal to x).

Then for any k

{ }

Wk induces a measure on the space of continuous functions on [0,1] . According to the invariance principle this measure converges weakly, as k→ ∞, to the Wiener process W [Chernov 1995 ] Fig. 10 depicts examples of functions

{ }

Wk for different k when the magnitudes Ak are chaotic variables generated by a tent map on [-1,1]:

1 1 2 , 1, 2,...

n n

A₊ = − A n= (42)

Fig. 10 Three realizations of the process W for k =100,300 and 10000 (red, green and blue

0 0.2 0.4 0.6 0.8 1

t -2

-1.5 -1 -0.5 0 0.5 1

WkHtL

( ) ( ) ( ) [ ]

1

1 , 0,1 , 0,1,...,

kt k t

L i

i

W t A t k L

D A L D A L

η_{ } ^{  }

=

= =

∑

∈ =

The weak invariance principle known also as the functional central limit theorem provides an approximation deterministic dynamical systems by a Brownian motion on large space and time scales.

Thus the distribution of η(t) tends to the Gaussian law with the mean (36) and variance (38).

This confirms the diffusion character of η(t). It follows that the Eq. (33) can be used for stochastic and chaotic modeling of the Wiener process.

Example 1. Valuation of the European call option.

The underlying asset of the European option is assumed to grow at the constant risk-free rate r perturbed by a stochastic/chaotic marked point process η(t). Thus an asset price is modeled as (43) Properties:

1) Markov property: the next asset price (S+dS) depends solely on today’s price 2) The next value for S is higher than the old by an amount

3) Variance of dS is

We want to price a call option, i.e.

(

^,

)

^{r t}

( )

^{r t}

( ( ) ) ( )

^|

( )

k

C t K =e^{− ⋅} E S−K ⁺=e^{− ⋅}

∑

P N t =k ⋅E S−K ⁺ N t =k, (44) where K is a strike price. We calculate a conditional expectation

( ) ( ) ( ) ( ) ( ( ) )

( ) (

^{( )}

( ) ) ^{( )}

( )

0

0

0

0

, |

ln

k ln

k

r t t

x x k

K K

K

C t K E S K N t k x K P S x N t k

x K P S e x N t k x K P x r t

S

x dx

x K p r t

S x

η

η

η

∞

+ +

∞ ∞

⋅ +

∞

 

=  − = = − ∂ ≤ = =

   

− ⋅ ∂ ≤ = = − ⋅ ∂  ≤  − ⋅ =

   

− ⋅   − ⋅ 

∫

∫ ∫

∫

% %

%

(45)

where ² ₂

( )

2 , r r D A

T σ σ

= − =

% _.

( ) (

²

)

²

^{( )} (

²

^{( )}

²

)

²

^{( )}

²

^{( ) ( )}

D dS =E dS −E dS =E S dη =S D dη =S D A dH t dS rdt d

S = + η

( )

⁽

( )

⁰⁾

E dS =rSdt as E dη =

A pdf of η_k can be found via its characteristic function. We note that

( ) ( )

1

k , ,

k

p k

p

E i A

ψη ω ω ω ω ω

=

 

=  = Θ



∑

 ^{K (46)}

where

( ) ( )

1 2 1 1

1 1

X

( , ,..., ) exp ^k exp ^k ( ,..., ) ...

k k i i i i A k k

i i

E i A i x p x x dx dx

= =

 

Θ ω ω ω =  ⋅∑ω = ⋅∑ω ⋅ ⋅

 

∫ ∫

^{L (47)}

is a k-dimensional characteristic function of a sequence

{

A₁,K,A_k

}

having a joint pdf

(

_k

)

A x x x

p ₁, ₂,K .

For a case of i.i.d. random values Ak (46) simplifies to

( )

ω

( )

ω

ψ_η ^k

k =Θ1 (48)

where ₁( )ω e^i⋅ωxp_A( )x dx

Χ

Θ =

∫

is the characteristic function of the distribution of Ak .

Here we focus on a special case of (46) when the magnitudes Ak are generated by a chaotic mapping

(

−1

)

= _k

k A

A ϕ (49)

in an interval X.

A joint pdf does not factorize in this case and calculates as

1

( )

1 2 1 1 1

1

( , ,..., ) ( ) ( ( ))

k

i

A k A i

i

p x x x p x δ x ϕ x

− +

=

= ⋅

∏

− ^{, (50)}

where pA (x) is the invariant density of the map ϕ .

The goal equation (46) simplifies for piece-wise linear onto maps

( )

x

{

i

( )

x a x b xi i^, J ii^{, 1, 2,...,}m

ϕ = ϕ = + ∈ = (51)

such that ∀i^{: :}ϕ Ji →X =

(

^0,1

)

.

We collect their main probabilistic properties:

• The invariant density is uniform with 1

A=2 and variance ² 1

A 12 σ =

• The autocorrelation function is

A property (39)-(40) can be easily illustrated with the exponentially decaying acf.

We next substitute the acf (52) into (37) and get

( )

( )

2 2 2

2

1 2 1

1 1

k

k

A A

r r

k r

r r

ση = σ ⁻ −σ ⁻

+ −

This confirms (39) at large k as r^k→0.

The characteristic function can be also calculated analytically. Substituting (50) into (47) for the inner integral we have

( )

1 1 1 1

1

1

1 2

( 1) 1 ( 1)

2 1 1 1 1 1 1

1 1 ( )

( 2) 2

1

(

3 2

1

( ( )) ... ( ( )) ( ( ( ))

1 ( ) ... ( )

1 ( ( )) ...

l l

l l

l l

k m

i v x k i v x i

k i l

l J i x z

m i v z b

a k

k

l l

x b

m i v

a k

k

l l

e x x x x dx e x x dx

e x z x z dz

a

e x x x

a

ϕ

δ ϕ δ ϕ δ ϕ ϕ

δ δ ϕ

δ ϕ δ ϕ

−

⋅ − ⋅ −

+

= =

Χ =

⋅ ⋅ −

−

= Χ

⋅ ⋅ −

=

⋅ − ⋅ ⋅ − = ⋅ −

= ⋅ ⋅ − ⋅ ⋅ −

= ⋅ ⋅ − ⋅ ⋅ −

∑ ∏

∫ ∫

∑ ∫

∑ (

^{−2)( ) .x2}

)

Hence the following recurrence equation can be obtained

1 1

1 1 2 3

1

( ,..., ) 1 ( , ,..., )

l l

m i b

a

k k k k

l l l

a e a

ω ω

ω ω ^{− ⋅ ⋅} ₋ ω ω ω

=

Θ =

∑

⋅ ⋅Θ +

the solution of which is

1 1

1

1 2 1

1 1 1 1

1 1

, ,..., 1 1 1

1 1

( ,..., ) ,

p

k k

n ip

n p n l n il

k n p

k i b k

m k

a

k k n k

i i i n i n p n i

e a

a

ω

ω ω ω ω

− −

= = =

−

− − ⋅∑ ⋅∑ ⋅∏ − −

= = = =

 

 

Θ = ⋅ Θ  ⋅ + 

 

∑ ∏ ∑ ∏

⁽⁵³⁾

where

1 1

0

( ) ( ) 1

i

i x i x

A

e p x dx e dx e i

ω

ω ω

ω _Χ ω

Θ =

∫

=

∫

= − is the characteristic function of the uniform distribution. Setting ω₁ =ω₂ =K=ω_k =ω in (53) and substituting the result into (46) we first get a characteristic function and then a required pdf of η_k by use an inverse Fourier transform. In [Baranovski 2003] authors have shown a fast convergence of the characteristic function Θ_k( ,..., )ω ω of the cumulative sum

1 k

k p

p

η A

=

=

∑

^to

( )

2 2

cos exp

2 ω ση

η ω ^− ^

 , which is the characteristic function of a normal distribution with the mean η and variance σ_η².

For example, a tent map on the unit interval has the following characteristic function [Baranovski 2003]