A random walk on moving spheres approach for the simulation of Bessel hitting times S. Herrmann University of Burgundy, Dijon, France joint work with M. Deaconu (INRIA Nancy, France) February th֒ S.Herrmann (University of Burgundy) Dijon, France 1 / 16 February th֒
Bessel hitting time Outline 1 Introduction to the Bessel first hitting time 2 Method of images and curved boundaries: the particular Bessel case 3 Walk on Moving Spheres Algorithm (WoMS) S.Herrmann (University of Burgundy) Dijon, France 2 / 16 February th֒
Bessel hitting time Introduction to the first Bessel hitting time 1) A financial model linked to the Bessel process The Cox-Ingersoll-Ross (1985) model permits to describe the short interest rates in finance. � dX δ t = ( a + bX δ X δ | X δ t ) dt + c t | dB t , 0 = x 0 with δ = 4 a / c 2 , x 0 ≥ 0, a ≥ 0, b ∈ R , c > 0. Aim: to simulate easily the first hitting time of the level L by the short rate. T L = inf { t ≥ 0 : X δ t = L } . Proposition (Îto+time change) The CIR model has the same distribution as X defined by: X t = e bt Y δ � c 2 � 4 b ( 1 − e − bt ) X 0 = Y δ ( 0 ) . , where Y δ is the squared Bessel process of dimension δ = 4 a c 2 . S.Herrmann (University of Burgundy) Dijon, France 3 / 16 February th֒
Bessel hitting time The squared Bessel process of dimension δ (index ν = δ/ 2 − 1): � t � Y δ t = y 0 + δ t + 2 | Y δ s | dB s . 0 The Bessel process Z δ t satisfies: � t t = x + δ − 1 s ) − 1 ds + B t , Z δ ( Z δ 2 0 Proposition � � T L has the same distribution as − 1 b log 1 − 4 b where c 2 τ ψ � � � τ ψ = inf { t ≥ 0 : Z δ t ≥ ψ ( t ) } with ψ ( t ) = 1 − 4 b L c 2 t δ = 4 a / c 2 et Z δ 0 . 0 = X δ It suffices to simulate the first time the Bessel process hits a curved boudary ψ . S.Herrmann (University of Burgundy) Dijon, France 4 / 16 February th֒
Bessel hitting time 2. Spiking neuron models: Feller’s model Inter-neuronal communication: the information is transmitted by the fluctuations of the electric potential membrane (spike trains). Classical model: integrate and fire (LIF model) based of the hitting times for an Ornstein-Uhlenbeck process. Feller’s model (more realistic): − X t � � � dX t = τ + µ dt + σ X t − V I dB t τ : caracteristics of the neuron (time constante of the membrane), µ linked to the input signal. Finally we also need to simulate the first time a Bessel process hits a given level L . S.Herrmann (University of Burgundy) Dijon, France 5 / 16 February th֒
Bessel hitting time The Bessel process of dimension δ > − 1 is solution of the SDE � t t = x + δ − 1 s ) − 1 ds + B t , Z δ ( Z δ 2 0 where ( B t , t ≥ 0 ) is a one-dimensional Brownian motion, x ≥ 0. Study of the first hitting time: τ L = inf { t ≥ 0 : Z δ t ≥ L } , L > x . Laplace transform We denote by I ν the modified Bessel function of order ν = δ 2 − 1. Then √ √ = x − ν 2 λ ) 2 λ ) ν 1 I ν ( x ( L � � � � e − λτ L e − λτ L , x > 0 , √ = √ E x E 0 2 ν Γ( ν + 1 ) L − ν 2 λ ) 2 λ ) I ν ( L I ν ( L S.Herrmann (University of Burgundy) Dijon, France 6 / 16 February th֒
Bessel hitting time The explicit expression of the Laplace transform: an eigenvalue problem associated with the infinitesimal generator 1 d x 2 ν + 1 d G ( ν ) = � � . 2 x 2 ν + 1 dx dx Then u ( x ) = E x [ e − λτ L ] satisfies G ( ν ) u = λ u , u ( L ) = 1 . Inversion of the Laplace transform: j ν − 1 ∞ j 2 1 ν, k � J ν + 1 ( j ν, k ) e − ν, k 2 L 2 t . P ( τ L > t ) = 2 ν − 1 Γ( ν + 1 ) k = 1 where J · is the Bessel function of the first kind and j · , k are the associated sequence of positive zeros. ➽ difficult to use for the simulation of τ L . Aim of the study Construction of a new algorithm in order to simulate the hitting time without using these Bessel functions and the Euler scheme (unboundedness of the stopping time). S.Herrmann (University of Burgundy) Dijon, France 7 / 16 February th֒
Bessel case Method of images Method of images and curved boundaries: the particular Bessel case It is difficult to handle with the hitting time of the level L . Is it easier to deal with a time-dependent boundary ? Let F ( dy ) be a positive σ -finite measure and a > 0 a constant. We define u ( t , x ) = p 0 ( t , x ) − 1 � p y ( t , x ) F ( dy ) , a R + a ψ ( t ) = inf { x ∈ R : u ( t , x ) < 0 } c τ = inf { t ≥ 0 : Z δ t ≥ ψ ( t ) } b Then � ψ ( t ) P 0 ( τ > t ) = u ( t , x ) dx . 0 S.Herrmann (University of Burgundy) Dijon, France 8 / 16 February th֒
Method of images In particular... � For F ( dy ) = y 2 ν + 1 1 { y > 0 } dy , we get ψ ( t ) = 2 t log Γ( ν + 1 ) t ν + 1 2 ν and a � 1 a � ν + 1 � � 2 t log P 0 ( τ ∈ dt ) = dt . 2 at Γ( ν + 1 ) t ν + 1 2 ν The case δ = 2 that is ν = 0 leads to simple expressions: � P 0 ( τ ∈ dt ) = 1 2 t log a 2 a log a and ψ a ( t ) = t dt . t ➽ τ can be easily simulated for this particular curved boundary. ➽ the method of images is useful for simulating the initial hitting time τ L . From now on, we fix ν = 0 ( δ = 2) ➽ This situation can be generalized to δ ∈ N ∗ . S.Herrmann (University of Burgundy) Dijon, France 9 / 16 February th֒
WoMS Algorithm Walk on Moving Spheres Algorithm (WoMS) �� � ( Z ( 2 ) ( B ( 1 ) ) 2 + ( B ( 2 ) , t ≥ 0 ) ∼ ) 2 , t ≥ 0 . t t t Simulating τ L is equivalent to simulate the first Brownian exit time from the sphere of radius L in R 2 . Idea for writing the algorithm: initialization: we start at the origin X ( 0 ) at time Θ 0 = 0 at each step n : we choose the constant a such that ψ a ( t ) is smaller than the distance between X ( n − 1 ) and ∂ D L for all t < a . We simulate the Brownian exit location of the sphere of radius ψ a : X ( n ) (uniform) and the exit time (method of image) θ n . We set Θ n = Θ n − 1 + θ n We stop when the norm of X ( n ) is ǫ -close to L . Number of steps denoted by S ǫ and the approximated hitting time is given by Θ S ǫ S.Herrmann (University of Burgundy) Dijon, France 10 / 16 February th֒
WoMS Algorithm Theorem 1 (Deaconu-H.) The number of steps S ǫ is almost surely finite. There exists a constant C > 0 and ǫ 0 > 0 such that E [ S ǫ ] ≤ C | log ǫ | , ∀ ǫ ≤ ǫ 0 . Theorem 2 (Deaconu-H.) As ǫ → 0, Θ S ǫ converges in distribution towards τ L . Moreover, for any small α > 0, ǫ � � 1 − F ǫ ( t − α ) ≤ F ( t ) ≤ F ǫ ( t ) , ∀ t > 0 . √ 2 απ S.Herrmann (University of Burgundy) Dijon, France 11 / 16 February th֒
WoMS Algorithm Algorithm (A1). Let us fix a parameter 0 < γ < 1 . Initialization: set X ( 0 ) = 0 , θ 0 = 0 , Θ 0 = 0 , A 0 = γ 2 L 2 e / 2 . First step : Let ( U 1 , V 1 , W 1 ) uniformly distributed on [ 0 , 1 ] 3 . θ 1 = A 0 U 1 V 1 , Θ 1 = Θ 0 + θ 1 , � cos ( 2 π W 1 ) � X ( 1 ) ⊺ = ( X 1 ( 1 ) , X 2 ( 1 )) ⊺ = X ( 0 ) ⊺ + ψ A 0 ( θ 1 ) . sin ( 2 π W 1 ) At the end of this step we set A 1 = γ 2 ρ ( X ( 1 )) 2 e / 2 . The n -th step : While X ( n − 1 ) ∈ D ε , simulate ( U n , V n , W n ) uniformly distributed on [ 0 , 1 ] 3 and define θ n = A n − 1 U n V n , Θ n = Θ n − 1 + θ n , � cos ( 2 π W n ) � X ( n ) ⊺ = ( X 1 ( n ) , X 2 ( n )) ⊺ = X ( n − 1 ) ⊺ + ψ A n − 1 ( θ n ) . sin ( 2 π W n ) At the end of this step we set A n = γ 2 ρ ( X ( n )) 2 e / 2 . ∈ D ε the algorithm is stopped. When X ( n − 1 ) / S.Herrmann (University of Burgundy) Dijon, France 12 / 16 February th֒
WoMS Algorithm Examples For L = 2, ν = 2, γ = 0 . 9 and 100000 simulations in order to evaluate the average number of steps. 10 − 1 10 − 2 10 − 3 10 − 4 ε E [ N ε ] 4.0807 7.53902 9.50845 10.83133 10 − 5 10 − 6 10 − 7 ε E [ N ε ] 10.94468 11.30869 11.62303 Figure: Average number of steps versus ε S.Herrmann (University of Burgundy) Dijon, France 13 / 16 February th֒
WoMS Algorithm For L = 2 with ε = 10 − 3 , γ = 0 . 9, 50000 simulations for each estimation of the average number of steps { 2 , 3 , . . . , 18 } 0 0.5 1 1.5 2 ν E [ N ε ] 6.819 7.405 8.270 8.887 9.594 2.5 3 3.5 4 ν E [ N ε ] 10.256 10.542 10.995 11.096 Figure: Average number of steps versus δ = 2 ν + 2 S.Herrmann (University of Burgundy) Dijon, France 14 / 16 February th֒
WoMS Algorithm Remarks... the algorithm can be generalized to any integer dimension δ we can reach other curved boundaries ψ : decreasing functions with bounded derivatives and � α > 0 , β > 0 . ψ ( t ) = α − β t , We therefore obtain a method to simulate the CIR hitting times for integers δ . this method can also be used to simulate the Brownian exit time of various bounded domains. Open question: what can we do for δ / ∈ N ? S.Herrmann (University of Burgundy) Dijon, France 15 / 16 February th֒
Recommend
More recommend
