A random walk on moving spheres approach for the simulation of - - PowerPoint PPT Presentation

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 February th֒  1 / 16

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 February th֒  2 / 16

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 δ t ) dt + c

• |X δ

t | dBt,

X δ

0 = x0

with δ = 4a/c2, x0 ≥ 0, a ≥ 0, b ∈ R, c > 0. Aim: to simulate easily the first hitting time of the level L by the short

• rate. TL = inf{t ≥ 0 : X δ

t = L}.

Proposition (Îto+time change) The CIR model has the same distribution as X defined by: X t = ebtY δc2 4b (1 − e−bt)

• ,

X 0 = Y δ(0). where Y δ is the squared Bessel process of dimension δ = 4a

c2 .

S.Herrmann (University of Burgundy) Dijon, France February th֒  3 / 16

Bessel hitting time

The squared Bessel process of dimension δ (index ν = δ/2 − 1): Y δ

t = y0 + δt + 2

t

• |Y δ

s | dBs.

The Bessel process Z δ

t satisfies:

Z δ

t = x + δ − 1

2 t (Z δ

s )−1 ds + Bt,

Proposition TL has the same distribution as − 1

b log

• 1 − 4b

c2 τψ

• where

τψ = inf{t ≥ 0 : Z δ

t ≥ ψ(t)} with ψ(t) =

• L
• 1 − 4b

c2 t

• δ = 4a/c2 et Z δ

0 = X δ 0 .

It suffices to simulate the first time the Bessel process hits a curved boudary ψ.

S.Herrmann (University of Burgundy) Dijon, France February th֒  4 / 16

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): dXt =

• − Xt

τ + µ

• dt + σ
• Xt − VI dBt

τ: 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 February th֒  5 / 16

Bessel hitting time

The Bessel process of dimension δ > −1 is solution of the SDE Z δ

t = x + δ − 1

2 t (Z δ

s )−1 ds + Bt,

where (Bt, 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

Ex

• e−λτL
• = x−ν

L−ν Iν(x √ 2λ) Iν(L √ 2λ) , x > 0, E0

• e−λτL
• =

(L √ 2λ)ν 2νΓ(ν + 1) 1 Iν(L √ 2λ)

S.Herrmann (University of Burgundy) Dijon, France February th֒  6 / 16

Bessel hitting time

The explicit expression of the Laplace transform: an eigenvalue problem associated with the infinitesimal generator G(ν) = 1 2x2ν+1 d dx

• x2ν+1 d

dx

• .

Then u(x) = Ex[e−λτL] satisfies G(ν)u = λu, u(L) = 1. Inversion of the Laplace transform: P(τL > t) = 1 2ν−1Γ(ν + 1)

∞

• k=1

jν−1

ν,k

Jν+1(jν,k) e−

j2 ν,k 2L2 t.

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 February th֒  7 / 16

Method of images

Bessel case

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) = p0(t, x) − 1 a

• R+

py(t, x)F(dy), ψ(t) = inf{x ∈ R : u(t, x) < 0} τ = inf{t ≥ 0 : Z δ

t ≥ ψ(t)}

c a b

Then P0(τ > t) = ψ(t) u(t, x) dx.

S.Herrmann (University of Burgundy) Dijon, France February th֒  8 / 16

Method of images

In particular... For F(dy) = y2ν+11{y>0}dy, we get ψ(t) =

• 2t log

a Γ(ν+1)tν+12ν and

P0(τ ∈ dt) = 1 2at

• 2t log

a Γ(ν + 1)tν+12ν ν+1 dt. The case δ = 2 that is ν = 0 leads to simple expressions: ψa(t) =

• 2t log a

t and P0(τ ∈ dt) = 1 2a log a t dt. ➽ τ 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 February th֒  9 / 16

WoMS Algorithm

Walk on Moving Spheres Algorithm (WoMS) (Z (2)

t

, t ≥ 0) ∼

• (B(1)

t

)2 + (B(2)

t

)2, t ≥ 0

• .

Simulating τL is equivalent to simulate the first Brownian exit time from the sphere of radius L in R2. 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 ∂DL 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 February th֒  10 / 16

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 −

ǫ √ 2απ

• F ǫ(t − α) ≤ F(t) ≤ F ǫ(t),

∀t > 0.

S.Herrmann (University of Burgundy) Dijon, France February th֒  11 / 16

WoMS Algorithm

Algorithm (A1). Let us fix a parameter 0 < γ < 1. Initialization: set X(0) = 0, θ0 = 0, Θ0 = 0, A0 = γ2L2e/2. First step : Let (U1, V1, W1) uniformly distributed on [0, 1]3.      θ1 = A0U1V1, Θ1 = Θ0 + θ1, X(1)⊺ = (X1(1), X2(1))⊺ = X(0)⊺ + ψA0(θ1) cos(2πW1) sin(2πW1)

• .

At the end of this step we set A1 = γ2ρ(X(1))2e/2. The n-th step : While X(n − 1) ∈ Dε, simulate (Un, Vn, Wn) uniformly distributed on [0, 1]3 and define      θn = An−1UnVn, Θn = Θn−1 + θn, X(n)⊺ = (X1(n), X2(n))⊺ = X(n − 1)⊺ + ψAn−1(θn) cos(2πWn) sin(2πWn)

• .

At the end of this step we set An = γ2ρ(X(n))2e/2. When X(n − 1) / ∈ Dε the algorithm is stopped.

S.Herrmann (University of Burgundy) Dijon, France February th֒  12 / 16

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 February th֒  13 / 16

WoMS Algorithm

For L = 2 with ε = 10−3, γ = 0.9, 50000 simulations for each estimation

• f the average number of steps {2, 3, . . . , 18}

ν 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 February th֒  14 / 16

WoMS Algorithm

Remarks... the algorithm can be generalized to any integer dimension δ we can reach other curved boundaries ψ: decreasing functions with bounded derivatives and ψ(t) =

• α − βt,

α > 0, β > 0. 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 February th֒  15 / 16

WoMS Algorithm

Algorithme (A2). Let 0 < γ < 1. Initialization: χ(0) = 0, ξ0 = 0, Ξ0 = 0, A0 = (γ2L2e/(ν + 1))ν+1 Γ(ν+1)

2

. Step n: while

• χ(n − 1) < L − ε, we choose Un uniformly distributed on

[0, 1]⌊ν⌋+2, Gn gaussian and Vn uniformly distr. on Sδ. We choose Un, Gn and Vn independent. Then        ξn =

• An−1

Γ(ν+1)2ν Un(1) . . . Un(⌊ν⌋ + 2)

• 1

ν+1 exp

• − ν−⌊ν⌋

ν+1 G 2 n

• ,

Ξn = Ξn−1 + ξn, χ(n) = χ(n − 1) + 2π1(Vn)

• χ(n − 1)ψAn−1(ξn) + ψ2

An−1(ξn).

End of the step: An =

• γ2(L −
• χ(n))2e/(ν + 1)

ν+1 Γ(ν + 1) 2 . If

• χ(n) ≥ l − ε then stop.

Output: Hitting time Ξn and the value of the Markov chain χ(n).

S.Herrmann (University of Burgundy) Dijon, France February th֒  16 / 16