http://www.iecn.u-nancy.fr/~etore/ ProjetOMEGA,IECN,BP239,F-54506Vand÷uvre-lès-NancyCEDEXPierre.Etore@iecn.u-nancy.fr PierreÉtoré Onrandomwalksimulationofone-dimensionaldiusionprocesseswithdiscontinuouscoecients

E l e c t r on ic

Jour n a l o f

Pr o

b a b i l i t y Vol. 11 (2006), Paper no. 9, pages 249275.

Journal URL

http://www.math.washington.edu/∼ejpecp/

On random walk simulation of one-dimensional diusion processes with discontinuous coecients

Pierre Étoré

Projet OMEGA, IECN, BP 239, F-54506 Vand÷uvre-lès-Nancy CEDEX Pierre.Etore@iecn.u-nancy.fr

http://www.iecn.u-nancy.fr/~etore/

Abstract

In this paper, we provide a scheme for simulating one-dimensional processes generated by divergence or non-divergence form operators with discontinuous coecients. We use a space bijection to transform such a process in another one that behaves locally like a Skew Brownian motion. Indeed the behavior of the Skew Brownian motion can easily be approached by an asymmetric random walk.

Key words: Monte Carlo methods, random walk, Skew Brownian motion, one-dimensional process, divergence form operator, local time.

AMS 2000 Subject Classication: Primary 660J60, 65C.

Submitted to EJP on May 11, 2005. Final version accepted on February 21, 2006.

∗Support acknowledgement: This work has been supported by the GdR MOMAS.

1 Introduction

In this paper we provide a random walk based scheme for simulating one-dimensional processes generated by operators of type

L=ρ 2∇

a∇

+b∇. (1.1)

These operators appear in the modelisation of a wide variety of diusion phenomena, for instance in uid mechanics, in ecology, in nance (see [DDG05]). Ifa,ρ, andbare measurable and bounded and ifaand ρare uniformly elliptic it can be shown thatLis the innitesimal generator of a Markov processX. Note that these operators contain the case of operators of typeL= ^a₂∆ +b∇whose interpretation in terms of Stochastic Dierential Equation is well known. In the general case there is still such an interpretation.

Indeed if we assume for simplicity thatb= 0it can be shown (see Section 4) thatX solves

X_t=X₀+ Z t

0

pa(X_s)ρ(X_s)dW_s+ Z t

0

ρ(Xs)a⁰(Xs)

2 ds+X

x∈I

a(x+)−a(x−)

a(x+) +a(x−)L^x_t(X), (1.2) whereI is the set of the points of discontinuity ofaandL^x_t(X)is the symmetric local time ofX in x. For the coecients a, ρ, and b may be discontinuous, providing a scheme to simulate trajectories of X is challenging: we cannot use the panel of Monte-Carlo methods available for smooth coecients (see [KP92]). However, some authors have recently provided schemes to simulateX in the case of coecients having some discontinuities.

In [Mar04] (see also [MT06]) M. Martinez treated the case of a coecient a having one point of discontinuity. He applied an Euler scheme after a space transformation that allows to get rid of the local time in (1.2). To estimate the speed of convergence of his method he needs a to be C⁶ outside the point of discontinuity. The initial condition has to be C⁴ almost everywhere and to satisfy other restrictive conditions.

In [LM06] (see also [Mar04]) A. Lejay and M. Martinez proposed a dierent scheme. After a piecewise constant approximation of the coecients and a dierent space transformation, they propose to use an exact simulation method of the Skew Brownian Motion (SBM). This one is based on the simulation of the exit times of a Brownian motion. In general the whole algorithm is slow and costly but allows to treat the case of coecientsa,ρandb being right continuous with left limits, and of classC¹except on countable set of points, without cluster point. Besides the initial condition can be taken inH¹, and the algorithm is well adapted to the case of coecients being at on large intervals outside their points of discontinuity.

Here, under the same hypotheses on a, ρ and b, but with the initial condition in W^1,∞∩H¹∩ C0 we propose a new scheme based mostly on random walks.

Roughly the idea is the following: assume for simplicity thatb = 0. First, like in [LM06], we replacea andρ by piecewise constantaⁿ and ρⁿ in order to obtain Xⁿ that is a good weak approximation of X and that solves

X_tⁿ=X₀ⁿ+ Z t

0

paⁿ(X_sⁿ)ρⁿ(X_sⁿ)dWs+ X

xⁿ_i∈Iⁿ

a(xⁿ_i+)−a(xⁿ_i−)

a(xⁿ_i+) +a(xⁿ_i−)L^x_tⁿⁱ(Xⁿ),

whereIⁿ is the set of the points of discontinuity ofaⁿ. Second by a proper change of scale we transform Xⁿ intoYⁿ that solves

Y_tⁿ=Y₀ⁿ+Wt+ X

xⁿ_k∈Iⁿ

β_kⁿL^k/n_t (Yⁿ),

where theβⁿ_k's are explicitely known. ThusYⁿ behaves around each k/nlike a SBM of parameter β_kⁿ (see Subsection5.2for a brief presentation of the SBM). That is, heuristically,Yⁿwhen ink/nmoves up with probability(β_kⁿ+ 1)/2and down with probability(β_kⁿ−1)/2, and behaves like a standard Brownian motion elsewhere. Thus a random walk on the grid{k/n:k∈Z}can reect the behaviour ofYⁿas was shown in [Leg85]. We use this to nally construct an approximationYbⁿ ofYⁿ.

We obtain a very easy to implement algorithm that only requires simulations of Bernoulli random variables. We estimate the speed of weak convergence of our algorithm by mixing an estimate of a weak error and an estimate of a strong error. Indeed computing the strong error of the algorithm presents diculties we were not able to overcome. On the other hand computing directly the weak error without using a strong error estimate would lead to more complicated computations without any improvement of the speed of convergence: basically our approach relies on the Donsker theorem and we cannot get better than an error inO(n^−1/2). Moreover to make such computations we should require additionnal smoothness on the data (see [Mar04]).

We nally make numerical experiments: the proposed scheme appears to be satisfying compared to the ones proposed in [Mar04] or [LM06].

Hypothesis. We make some assumptions from now till the end of the paper, for the sake of simplicity but without loss of generality.

(A1)b= 0.

Indeed, as explained in [LM06] Section 2, if we can treat the case L= ρ

2∇ a∇

, (1.3)

we can treat the case (1.1) for any measurable bounded bby dening the coecientsρand ain (1.3) in the following manner:

If ab :=aexp^Ψ andρb :=ρexp^−Ψ,

with Ψ(x) = Z x

0

h(y)dy andh(x) = 2 b(x) ρ(x)a(x), then ρb

2∇ ab∇

= ρ 2∇

a∇

+b∇.

(A2) Let beG= (l, r)an open bounded interval ofR. We will assume the processX starts fromx∈G and is killed when reaching{l, r}. >From a PDEs point of view this means the parabolic PDEs involving L we will study are submitted to uniform Dirichlet boundary conditions. We could treat Neumann boundary conditions (thanks to the results of [BC05] for instance) and, by localization arguments, the case of an unbounded domainG(see [LM06]). But this assumption will make the material of the paper simpler and clearer.

Outline of the paper. In Section 2 we dene precisely Divergence Form Operators (DFO) and recall some of their properties. In Section 3 we speak of Stochastic Dierential Equations involving Local Time (SDELT): we state a general change of scale formula and recall some convergence results established by J.F. Le Gall in [Leg85]. In Section 4 we link DFO and SDELT: a process generated by a DFO of coecientsaandρhaving countable discontinuities without cluster points is solution of a SDELT with coecients determined by a andρ. In Section 5 we present our scheme. In Section 6 we estimate the speed of convergence of this scheme. Section 7 is devoted to numerical experiments.

Some notations. For1≤p <∞we denote byL^p(G)the set of measurable functionsf onGsuch that kfk_p:=

Z

G

|f(x)|^pdx1/p

<∞.

Let be0< T <∞xed. For1≤p, q <∞we denote byL^q(0, T; L^p(G))the set of measurable functions f on(0, T)×Gsuch that

|||f|||_p,q :=

Z T 0

kfk^q_pdt^1/q

<∞.

For u∈ L^p(G) we denote by ^du_dx the rst derivative of u in the distributional sense. It is standard to denote byW^1,p(G)the space of functionsu∈L^p(G)such that ^du_dx ∈L^p(G), and byW₀^1,p(G)the closure of C_c^∞(G)inW^1,p(G)equipped with the norm(kuk^p_p+

^du_dx

p

p)^1/p. We denote byH¹(G)the spaceW^1,2(G), and byH¹₀(G)the spaceW^1,2₀ (G).

For u ∈ L²(0, T; L²(G)) we denote by ∂tu the distribution such that for all ϕ ∈ C_c^∞((0, T)×G), we have,h∂tu, ϕi =−RT

0

R

Gu∂tϕ. We still denote by ^du_dx the rst derivative of uwith respect to xin the distributional sense.

We will classicaly denote byk.k_∞ and|||.|||_∞,∞ the supremum norms.

For the use of probability we will denote by C0(G)the set of continuous bounded functions on G. The symbol∼=will denote equality in law.

2 On divergence form operators

For0< λ <Λ<∞let us denote byEll(λ,Λ)the set of functionsf onGthat are measurable and such that

∀x∈G, λ≤f(x)≤Λ.

Forρ∈Ell(λ,Λ)let us dene the measuremρ(dx) :=ρ⁻¹(x)dx.

For any measure m with a bounded density with respect to the Lebesgue measure we then denote by L²(G, m)the Hilbert space of functions inL²(G)equipped with the scalar product

(f, g)7−→

Z

G

f(x)g(x)m(dx).

This is done in order that the operator we dene below is symmetric onL²(G, mρ).

Denition 2.1 Letaandρbe inEll(λ,Λ)for some0< λ <Λ<∞. We call Divergence form operator of coecients aandρ, and we noteL(a, ρ), the operator(L, D(L))on L²(G, mρ)dened by

L=ρ 2

d dx

a d

dx

, and

D(L) ={u∈H¹₀(G), Lu∈L²(G)}.

Actually ifa, ρ ∈Ell(λ,Λ)the operatorL(a, ρ)has sucient properties to generate a continuous Markov process. We sum up these properties in the next theorem.

Theorem 2.1 Letaandρbe in Ell(λ,Λ)for some 0< λ <Λ<∞. Then we have:

i) The operatorL(a, ρ)on L²(G, mρ)is closed and self-adjoint, with dense domain.

ii) This operator is the innitesimal generator of a strongly continuous semigroup of contraction (St)_t≥0 onL²(G, mρ).

iii) Moreover (St)_t≥0 is a Feller semigroup. Thus L(a, ρ) is the innitesimal generator of a Markov process(Xt, t≥0).

iv) The process(X_t, t≥0) has continuous trajectories.

Proof. We give the great lines of the proof and refer the reader to [Lej00] and [Str82] for details.

We set(L, D(L)) =L(a, ρ). First it is possible to build a symmetric bilinear formE onL²(G, m_ρ)dened by

E(u, v) = Z

G

a 2

du dx

dv

dxdx, ∀(u, v)∈ D(E)×D(E), and D(E) = H¹₀(G), that veries,

E(u, v) =− hLu, vi_L2(G,m_ρ), ∀(u, v)∈D(L)×D(E).

Thus the resolvent of(L, D(L))can be built and we get i). An application of the Hille-Yosida theorem then leads to ii).

Further it is a classical result of PDEs that the semigroup(Pt)_t≥0 has a densityp(t, x, y)with respect to the measuremρ such that

u(t, x) = Z

G

p(t, x, y)f(y)ρ⁻¹(y)dy (2.1)

is a continuous version ofP_tf(x). Then forf inC₀(G),P_tf belongs toC₀(G). By the use of the maximum principle it can be shown that(P_t)_t≥0 is semi-markovian and we get iii).

Finally Aronson estimates on the densityp(t, x, y)can be used to show for example that

∀ε >0, ∀x∈G, lim

t↓0

1 t

Z

|y−x|>ε

p(t, x, y)ρ⁻¹(y)dy= 0, and thus we get iv) (see Proposition 2.9 in chapter 4 of [EK86]).

We have a consistency theorem.

Theorem 2.2 Let be 0 < λ < Λ < ∞. Let a and ρ be in Ell(λ,Λ) and (aⁿ, ρⁿ) be a sequence of Ell(λ,Λ)×Ell(λ,Λ).

Let us denote by S and X respectively the semigroup and the process generated by L(a, ρ) and by(Sⁿ) and(Xⁿ) the sequences of semigroups and processes generated by the sequence of operatorsL(aⁿ, ρⁿ). Assume that

1 aⁿ

L²(G)

−−−−*

n→∞

1

a, and 1 ρⁿ

L²(G)

−−−−*

n→∞

1 ρ. Then for anyT >0 and any f ∈L²(G)we have :

i) The functionS_tⁿf(x)converges weakly inL²(0, T; H¹₀(G))toStf(x).

ii) The continuous version of S_tⁿf(x) given by (2.1) with preplaced by pⁿ converges uniformly on each compact of(0, T)×Gto the continuous version of S_tf(x) given by (2.1).

iii)

(X_tⁿ, t≥0)−−−−→^L

n→∞ (Xt, t≥0).

Proof. See in [LM06] the proofs of Propositions 3 and 4.

3 On SDE involving local time

First we introduce a new class of coecients. For0 < λ <Λ <∞ we denote byCoeff(λ,Λ)the set of the elementsf ofEll(λ,Λ)that verify:

i) f is right continuous with left limits (r.c.l.l.).

ii) f belongs toC¹(G\ I), whereI is a countable set without cluster point.

Let us also denote byMthe space of all bounded measuresν onGsuch that|ν({x})|<1for allxinG. Denition 3.1 Let σ be inCoeff(λ,Λ) for some 0< λ <Λ<∞, and ν be in M. We call Stochastic Dierential Equation with Local Time of coecientsσandν, and we noteSde(σ, ν), the following SDE

X_t=X₀+ Z t

0

σ(X_s)dW_s+ Z

R

ν(dx)L^x_t(X), whereL^x_t(X)is the symmetric local time of the unknown processX.

In [Leg85] J.F. Le Gall studied some properties of SDEs of the typeSde(σ, ν). We will recall here some results of this work we will use in the sequel.

We will see below that σ ∈ Coeff(λ,Λ) and ν ∈ Mis a sucient condition to have a unique strong solution toSde(σ, ν). We rst x some additionnal notations.

Forf inCoeff(λ,Λ)we denote byf⁰(dx)the bounded measure corresponding to the rst derivative off in the generalized sense. We denote byf(x+) andf(x−)respectively the right and left limits off inx. We will also denote byf⁰(x)the r.c.l.l. density of the absolutely continuous part off⁰(dx)(that it is to say we take forf⁰(x)the right derivative of f in x).

Forν inMwe denote byν^c the absolutely continuous part ofν.

3.1 A change of scale formula

Let us dene the class of bijections we will use in our change of scale.

For0< λ <Λ<∞we denote byT(λ,Λ)the set of all functionsΦonGthat have a rst derivativeΦ⁰ that belongs toCoeff(λ,Λ). The assumption made on the bijection is minimal and we can then state a very general change of scale formula.

Proposition 3.1 Letσbe in Coeff(λ,Λ)for some0< λ <Λ<∞. Let ν(dx) =b(x)dx+ X

xi∈I

cx_iδx_i(dx), be in M(i.e., b is measurable and bounded, and each|cx_i|<1).

Let Φbe inT(λ⁰,Λ⁰)for some0< λ⁰<Λ⁰<∞and letJ be the set of the points of discontinuity of Φ⁰. Then the next statements are equivalent:

i) The processX solves Sde(σ, ν).

ii) The process Y := Φ(X)solves Sde(γ, µ)with

γ(y) = (σΦ⁰)◦Φ⁻¹(y), and

µ(dy) = Φ⁰b+¹₂(Φ⁰)⁰

(Φ⁰)² ◦Φ⁻¹(y)dy+ X

xi∈I∪J

βx_iδ_Φ(xi)(dy), where,

β_x= Φ⁰(x+)(1 +c_x)−Φ⁰(x−)(1−c_x) Φ⁰(x+)(1 +cx) + Φ⁰(x−)(1−cx), withc_x= 0 if x∈ J \ I.

Remark 3.1 Note thatγobviously belongs toCoeff(λ⁰⁰,Λ⁰⁰)for some0< λ⁰⁰<Λ⁰⁰<∞, and thatµis inM, so it makes sense to speak ofSde(γ, µ).

We can say that the class of SDE of typeSde(σ, ν)is stable by transformation by a bijection belonging toT(λ,Λ)for some 0< λ <Λ<∞.

Proof of Proposition 3.1. We provei)⇒ii). The converse can be proven in the same manner quite being technically more cumbersome.

By the symmetric Itô-Tanaka formula we rst get:

Y_t= Φ(X_t) = Φ(X₀) +Rt

0(σΦ⁰)(X_s)dW_s+Rt

0(σ²bΦ⁰)(X_s)ds +P

x_i∈I

Φ⁰(x_i+)+Φ⁰(x_i−)

2 cx_iL^x_tⁱ(X) +¹₂Rt

0[σ²(Φ⁰)⁰](X_s)ds+P

x_i∈J

Φ⁰(x_i+)−Φ⁰(x_i−) 2 L^x_tⁱ(X)

= Φ(X0) +Rt

0(σΦ⁰)◦Φ⁻¹(Ys)dWs

+Rt

0[σ²(Φ⁰b+¹₂(Φ⁰)⁰)]◦Φ⁻¹(Y_s)ds+P

xi∈I∪JK_xiL^x_tⁱ(X), withK_x=c_x(Φ⁰(x+) + Φ⁰(x−))/2 + (Φ⁰(x+)−Φ⁰(x−))/2.

We have then to expressL^x_t(X)in function ofL^Φ(x)_t (Y)forx∈ I ∪ J.

Using Corollary VI.1.9 of [RY91] it can be shown that L^Φ(x)

±

t (Y) = Φ⁰(x±)L^x±_t (X). (3.1)

Besides theorem VI.1.7 of [RY91] leads to (L^x+_t (X)−L^x−_t (X))/2 = cxL^x_t(X) and combining with (L^x+_t (X) +L^x−_t (X))/2 =L^x_t(X)we get

L^x+_t (X) = (1 +cx)L^x_t(X). (3.2)

In a similar manner we have(L^Φ(x)+_t (Y)−L^Φ(x)−_t (Y))/2 =KxL^x_t(X)and we can get KxLt(X) +L^Φ(x)_t (Y) =L^Φ(x)+_t (Y).

Then using (3.1) and (3.2) we get

L^Φ(x)_t (Y) = (Φ⁰(x+)(1 +cx)−Kx)L^x_t(X), and the formula is proved.

To prove the proposition below, Le Gall used in [Leg85] a space bijection that enters in the general setting of Proposition3.1.

Proposition 3.2 (Le Gall 1985) Let σ be in Coeff(λ,Λ) for some 0 < λ <Λ<∞ and ν be in M. There is a unique strong solution toSde(σ, ν).

We need two lemmas.

Lemma 3.1 (Le Gall 1985) Let σ be in Coeff(λ,Λ) for some 0 < λ < Λ < ∞. There is a unique strong solution toSde(σ,0).

Proof. See [Leg85].

The next lemma will play a great role for calculations in the sequel.

Lemma 3.2 Let ν be in M. There exists a function fν in Coeff(λ,Λ) (for some 0 < λ < Λ < ∞), unique up to a multiplicative constant, such that:

f_ν⁰(dx) + (fν(x+) +fν(x−))ν(dx) = 0. (3.3) Proof of Proposition 3.2. It suces to set

Φν(x) = Z x

0

fν(y)dy.

By Lemma3.2Φ_νobviously belongs toT(λ,Λ)for some0< λ <Λ<∞. By Proposition3.1and Lemma 3.2 we get thatX solvesSde(σ, ν)if and only if Y := Φ_ν(X) solvesSde (σf_ν)◦Φ⁻¹_ν ,0. By Lemma 3.1the proof is completed.

3.2 Convergence results

Le Gall proved the next consistency result for equations of the typeSde(σ, ν).

Theorem 3.1 (Le Gall 1985) Let be two sequences(σⁿ)and(νⁿ)for which there exist0< λ <Λ<∞, 0< M <∞ andδ >0 such that

(H1) σⁿ ∈Coeff(λ,Λ), ∀n∈N,

(H2) |νⁿ({x})| ≤1−δ, ∀n∈N,∀x∈G.

(H3) |νⁿ|(G)≤M, ∀n∈N,

so that each νⁿ is in M. Assume that there exist two functions σ and f in Coeff(λ⁰,Λ⁰) (for some 0< λ⁰<Λ⁰ <∞) such that

σn L¹_loc(R)

−−−−−→

n→∞ σ and fνⁿ

L¹_loc(R)

−−−−−→

n→∞ f, and set:

ν(dx) =− f⁰(dx)

f(x+) +f(x−). (3.4)

Let (Ω,F,(Ft)_t≥0,P^x) be a ltered probability space carrying an adapted Brownian motion W. On this space, for each n∈N let be Xⁿ the strong solution of Sde(σⁿ, νⁿ), and let be X the strong solution of Sde(σ, ν). Then:

E[ sup

0≤s≤t

|X_sⁿ−Xs|]−−−−→

n→∞ 0 and (X_tⁿ, t≥0)−−−−→^L

n→∞ (Xt, t≥0).

Remark 3.2 In this theoremνⁿapproachesν in the sense thatf_νntends tof_ν for theL¹_locconvergence.

Note we can haveνⁿ * ν₁, butf_νn →f_ν2 for theL¹_loc convergence, with ν₁6=ν₂ (See [Leg85] p 65 for an example). The theorem asserts thatXⁿ tends toX that solvesSde(σ, ν₂)and notSde(σ, ν₁)! Le Gall also proved a Donsker theorem for solution to SDEs of the typeSde(σ, ν)forσ≡1. Let beµ inMandY be the solution toSde(1, µ).

We dene some coecientsβ_kⁿ for allk∈Z, and alln∈N^∗, by:

1−β_kⁿ

1 +β_kⁿ = exp −2µ^c(]k n,k+ 1

n ]) Y

k n<y≤^k+1_n

1−µ({y}) 1 +µ({y})

= fµ(^k+1_n )

fµ(^k_n) . (3.5) We dene a sequence(µⁿ)of measures inMby

µⁿ=X βⁿ_kδk

n, (3.6)

and a sequence(Yⁿ)of processes such that eachYⁿ solvesSde(1, µⁿ). Finally we dene for alln∈N^∗a sequence (τ_pⁿ)_p∈Nof stopping times by,

τ₀ⁿ= 0

τ_p+1ⁿ = inf{t > τ_pⁿ :|Y_tⁿ−Y_τⁿn

p|= ¹_n}. (3.7)

We have the next theorem.

Theorem 3.2 (Le Gall 1985) In the previous contextS_pⁿ:=nY_τⁿn

p denes a sequence of random walks on the integers such that:

i)

Sⁿ₀ = 0, ∀n∈N^∗,

P[S_p+1ⁿ =k+ 1|S_pⁿ=k] = ¹₂(1 +β_kⁿ), ∀n, p∈N^∗,∀k∈Z, P[S_p+1ⁿ =k−1|S_pⁿ=k] = ¹₂(1−β_kⁿ), ∀n, p∈N^∗,∀k∈Z. ii)

The sequence of processes dened by Ye_tⁿ := (1/n)S_[nⁿ2t], where b.c stands for the integer part of a non negative real number, veries for all0< T <∞:

E|Ye_tⁿ−Yt| −−−−→

n→∞ 0, ∀t∈[0, T] and (Ye_tⁿ, t≥0)−−−−→^L

n→∞ (Yt, t≥0).

4 Link between DFO and SDELT

This link is stated by the following proposition.

Proposition 4.1 Letaandρbe inCoeff(λ,Λ)for some0< λ <Λ<∞. Let us denote byI the set of the points of discontinuity ofa. ThenL(a, ρ)is the innitesimal generator of the unique strong solution ofSde(√

aρ, ν)with,

ν(dx) =a⁰(x) 2a(x)

dx+ X

xi∈I

a(xi+)−a(xi−)

a(xi+) +a(xi−)δx_i(dx). (4.1) In [LM06] the authors proved the proposition above by the use of Dirichlet forms and Revuz measures.

We give here a more simple proof, based on smoothing the coecients and using the consistency theorems of the two preceeding sections.

Proof of Proposition 4.1. Asaandρare inCoeff(λ,Λ)the function√

ρais inCoeff(λ,Λ). Besides, as|a−b|/|a+b|<1for anya,binR^∗+, the measureνdened by (4.1) is inM. The existence of a unique strong solutionX toSde(√

aρ, ν)follows from Proposition3.2.

We then identify the innitesimal generator ofX. We can build two sequences(aⁿ)and(ρⁿ)of functions inCoeff(λ,Λ)∩ C^∞(G), such that

aⁿ −−−−→

n→∞ a a.e. and ρⁿ−−−−→

n→∞ ρ a.e.

For anyninNwe denote byXⁿ the process generated byL(aⁿ, ρⁿ).

On one hand, by dominated convergence the hypotheses of Theorem2.2are fullled, and we have,

(X_tⁿ, t≥0)−−−−→^L

n→∞ (Xet, t≥0), (4.2)

where the processXe is generated byL(a, ρ).

On the other hand we will show by Theorem3.1that (X_tⁿ, t≥0)−−−−→^L

n→∞ (Xt, t≥0). (4.3)

Thus taking in care (4.2) and (4.3) we will conclude that the innitesimal generator ofX isL(a, ρ). Asaⁿ andρⁿ areC^∞,(Lⁿ, D(Lⁿ)) =L(aⁿ, ρⁿ)can be written,

Lⁿ= ρⁿ 2

h aⁿ⁰ d

dx +aⁿ d² dx²

i , so it is standard to say thatXⁿ solves

X_tⁿ=x+ Z t

0

pρⁿ(X_sⁿ)aⁿ(X_sⁿ)dWs+ Z t

0

ρⁿ(X_sⁿ)aⁿ⁰(X_sⁿ)

2 ds. (4.4)

AsdhXⁿi_s=ρⁿ(X_sⁿ)aⁿ(X_sⁿ)ds, by the occupation time density formula we can rewrite (4.4) and assert thatXⁿ solves,

X_tⁿ=x+ Z t

0

pρⁿ(X_sⁿ)aⁿ(X_sⁿ)dWs+ Z

R

νn(dx)L^x_t(Xⁿ), whereνn(dx) = (aⁿ⁰(x)/2aⁿ(x))λ(dx).

Then elementary calculations show that the functionfν_n associated toνn by Lemma3.2 is of the form fν_n(x) = K/aⁿ(x) with K a real number. This obviously tends to K/a(x) =: f(x) for the L¹_loc(R)- convergence. We then determine the measure ν associated to f by (3.4). First we check that ν^c(dx) = (a⁰(x)/2(a(x))λ(dx). Second, the set{x∈G: ν({x})6= 0}is equal toI, and we have for allx∈ I,

ν({x}) =−f(x+)−f(x−)

f(x+) +f(x−) =a(x+)−a(x−) a(x+) +a(x−). So the measureν is equal to the one dened by (4.1). As it is obvious that

√ρⁿa^{n L}

1 loc(R)

−−−−−→

n→∞

√ρa,

and that the hypotheses(H1)−(H3)of Theorem3.1are fullled, we can say that (4.3) holds. The proof is completed.

5 Random walk approximation

5.1 Monte Carlo Approximation

From now the horizon 0 < T <∞ is xed. For any a, ρ ∈Coeff(λ,Λ) and any initial conditionf we denote by(P)(a, ρ, f)the parabolic PDE

(P)(a, ρ, f)















∂u(t,x)

∂t =Lu(t, x), for(t, x)∈[0, T]×G, u(t, l) =u(t, r) = 0 fort∈[0, T],

u(0, x) =f(x) forx∈G, with(L, D(L)) =L(a, ρ).

Let be0< λ <Λ<∞. From now till the end of this paper we assume that aandρare inCoeff(λ,Λ). We denote byI={xi}i∈I the set of the points of discontinuity ofa(I={0≤i≤k₁} ⊂Zis nite). We setX to be the process generated byL(a, ρ).

We seek for a probabilistic numerical method to approximate the solution of (P)(a, ρ, f). By Theorem 2.1and some standard PDEs renements we know that for allf ∈L²(G),(P)(a, ρ, f)has a unique weak solutionu(t, x)inC([0, T],L²(G, m_ρ))∩L²(0, T; H¹₀(G)). We know thatE^x[f(X_t)]is a continuous version ofu(t, x).

Our goal is to build a processXbⁿ easy to simulate and such that (Xb_tⁿ, t≥0)−−−−→^L

n→∞ (X_t, t≥0). (5.1)

ThusE^x[f(Xb_tⁿ)]→E^x[f(Xt)]for anyt∈[0, T], and the strong law of large numbers asserts that, 1

N

N

X

i=1

f(Xb_t^n,(i))−−−−→^n→∞

N→∞ u(t, x), (5.2)

where for eachi,Xb^n,(i) is a realisation of the random variableXbⁿ.

5.2 Skew Brownian Motion

The Skew Brownian Motion (SBM) of parameterβ ∈(−1,1)starting from y, which we denote byY^β,y is known to solve:

Y_t^β,y=y+Wt+βL^y_t(Y^β,y), (5.3) i.e. Y^β,y solvesSde(1, βδ0)(see [HS81]).

It was rst constructed by Itô and McKean in [IM74] (Problem 1 p115) by ipping the excursions of a reected Brownian motion with probabilityα= (β+ 1)/2. On SBM see also [Wal78] . It behaves like a Brownian motion except iny where its behaviour is pertubated, so that

P(Y_t^β,y > y) =α, ∀t >0. (5.4)

We denote by T(∆) the law of the stopping time τ = inf{t ≥ 0, |Wt| = ∆} where W is a standard Brownian motion starting at zero. For the SBM Y^β,0 and ∆ in R^∗+ we dene the stopping time τ∆ = inf{t≥0, Y_t^β,0∈ {∆,−∆}}. Our approach relies on the following lemma.

Lemma 5.1 Lety andxbe in R,∆ in R^∗+ andβ in (−1,1). Setα= (β+ 1)/2. Then i)Y^β,y+x∼=Y^β,y+x, and ii)(τ∆, Y_τ^β,0_∆ )∼=T(∆)⊗Ber(α).

Proof. This follows from the construction of the SBM by Itô and McKean in [IM74].

5.3 Some possible approaches

Recently two methods have been proposed to buildXbⁿ satisfying (5.1).

In [Mar04] M. Martinez proposed to use an Euler scheme. We know by Proposition4.1 that X solves Sde(√

aρ, ν) with ν dened by (4.1). We have seen in the proof of Proposition 3.2 that if we dene Φ(x) =Rx

0 fν(y)then Y = Φ(X)solvesSde(γ,0) withγ in some Coeff(m, M). Thus an Euler scheme approximation Ybⁿ of Y can be built and by setting Xbⁿ = Φ⁻¹(Ybⁿ) we get an approximation of X. Because the coecientγis not Lipschitz ifaandρare not, evaluating the speed of convergence of such a scheme is not easy.

In [LM06], A. Lejay and M. Martinez proposed to use the SBM. They rst build a piecewise constant approximation (aⁿ, ρⁿ) of (a, ρ) in order that the process Xⁿ generated by L(aⁿ, ρⁿ) solves Sde(√

aⁿρⁿ, νⁿ)with νⁿ satisfying(νⁿ)^c = 0. Second by a proper bijection Φⁿ ∈ T(m, M), they get that Yⁿ = Φⁿ(Xⁿ) solves Sde(1, µⁿ) with µⁿ = P

βkδy_k, i.e. Yⁿ behaves locally like a SBM. Third they proposed a schemeYbⁿ forYⁿ based on Lemma5.1and simulations of exit times of the SBM and they nally setXbⁿ= (Φⁿ)⁻¹(Ybⁿ).

Our method can be seen as a variation of this last approach because it also deeply relies on getting such aYⁿ and using Lemma 5.1. But we then use random walks instead of the scheme proposed in [LM06].

5.4 The basic idea of our approach

We focus on weak convergence and propose a three-step approximation scheme diering slightly from the one proposed by Theorem3.2.

We xn∈N^∗, and1/n will be the spatial discretization step size.

STEP 1. We build(aⁿ, ρⁿ)in Coeff(λ,Λ)×Coeff(λ,Λ)such that:

i) The functions aⁿ and ρⁿ are piecewise constant. The points of discontinuity of either aⁿ and ρⁿ are included in some set Iⁿ. We assume Iⁿ = {xⁿ_k}_k∈Iⁿ for Iⁿ = {0 ≤ k ≤ k₁ⁿ} ⊂ Z nite and xⁿ_k < xⁿ_k+1, ∀k∈Iⁿ.

ii) For eachxⁿ_k ∈ Iⁿ we haveaⁿ(xⁿ_k) =a(xⁿ_k)andρⁿ(xⁿ_k) =ρ(xⁿ_k). iii) Consider the function

Φⁿ(x) =

k_n,x−1

X

k=0

xⁿ_k+1−xⁿ_k

pa(xⁿ_k)ρ(xⁿ_k)+ x−xⁿ_k

n,x

qa(xⁿ_k

n,x)ρ(xⁿ_k

n,x)

, (5.5)

where the integerk_n,xveriesxⁿ_k

n,x ≤x≤xⁿ_k

n,x+1.

The setIⁿ satisesΦⁿ(Iⁿ) = {k/n, k∈Z} ∩Φⁿ(G). From now we assume xⁿ_k is the point of Iⁿ such thatΦⁿ(xⁿ_k) =k/n.

Remark 5.1 In fact the rst thing to do is to construct the grid Iⁿ satisfying iii). It is very easy and only requires to know the coecientsaandρ(see point 1 of the algorithm in Subsection 5.5). Thenaⁿ andρⁿ can be constructed.

Remark 5.2 The setsI andIⁿ may have no common points.

We takeXⁿ to be the process generated byL(aⁿ, ρⁿ).

Remark 5.3 It can be shown by Theorem2.2thatXⁿ converges in law toX. STEP 2. By Proposition4.1the processXⁿ solvesSde(√

aⁿρⁿ, νⁿ)with νⁿ= X

xⁿ_k∈Iⁿ

aⁿ(xⁿ_k+)−aⁿ(xⁿ_k−) aⁿ(xⁿ_k+) +aⁿ(xⁿ_k−)δ_xⁿ

k.

The functionΦⁿ dened by (5.5) belongs toT(1/Λ,1/λ). The points of discontinuity ofΦⁿ⁰ are those in Iⁿ, and(Φⁿ⁰)⁰= 0, so by Proposition3.1the processYⁿ= Φⁿ(Xⁿ)solves

Y_tⁿ=Y₀ⁿ+Wt+ X

xⁿ_k∈Iⁿ

β_kⁿL^k/n_t (Yⁿ), (5.6)

where

βⁿ_k =

pa(xⁿ_k)/ρ(xⁿ_k)−q

a(xⁿ_k−1)/ρ(xⁿ_k−1) pa(xⁿ_k)/ρ(xⁿ_k) +q

a(xⁿ_k−1)/ρ(xⁿ_k−1)

. (5.7)

To write these coecients we have used the fact thataⁿandρⁿare r.c.l.l. and that for instanceaⁿ(xⁿ_k+) = a(xⁿ_k)andaⁿ(xⁿ_k−) =a(xⁿ_k−1).

Remark 5.4 We have gotYⁿ that solvesSde(1,Pβⁿ_kδ_k/n)in a dierent way than the one used by Le Gall in Theorem3.2. We now use his method to getYeⁿ that veries

E|Ye_tⁿ−Y_tⁿ| −−−−→

n→∞ 0, ∀t∈[0, T]. (5.8)

STEP 3. Like in (3.7) we dene a sequence(τ_pⁿ)_p∈Nof stopping times by, τ₀ⁿ= 0, and τ_p+1ⁿ = inf{t > τ_pⁿ:|Y_tⁿ−Y_τⁿn

p|= 1/n}.

Thanks to the uniformity of the grid{k/n, k∈Iⁿ}we have the following lemma.

Lemma 5.2 i) For allk∈Zand allp∈N,(Y_τⁿn

p+u−Y_τⁿn

p,0≤u≤τ_p+1ⁿ −τ_pⁿ)knowing that{Y_τⁿn

p =k/n}

has the same law as(Y_t^βkn,0,0≤t≤τ_1/n). ii)

∀p∈N, σⁿ_p :=n² τ_pⁿ−τ_p−1ⁿ ∼=T(1), and the σ_pⁿ's are independent.

Proof. The statement i) follows simply from point i) of Lemma5.1and the comparison between (5.3) and (5.6). From i), the strong Markov property, and point ii) of Lemma5.1, we get that(τ_pⁿ−τ_p−1ⁿ )∼=T(1/n²), and the statement ii) follows by scaling. Using again the Markov property we get the independence of theσ_pⁿ's.

Le Gall used this lemma in his proof of Theorem3.2. Indeed it is obvious by the i) of Lemma 5.2 and the ii) of Lemma5.1thatS_pⁿ:=nY_τⁿn

p satises

S₀ⁿ= 0,

P[S_p+1ⁿ =k+ 1|S_pⁿ=k] = ¹₂(1 +β_kⁿ) =:αⁿ_k, ∀p∈N^∗,∀k∈Iⁿ, P[S_p+1ⁿ =k−1|S_pⁿ=k] = ¹₂(1−β_kⁿ) = 1−αⁿ_k, ∀p∈N^∗,∀k∈Iⁿ.

(5.9)

Moreover the ii) of Lemma5.2allows to show thatYe_tⁿ:= (1/n)S_[nⁿ2t]=Y_τⁿn

[n2t] satises (5.8).

Thus the idea is to take

Xb_tⁿ:= (Φⁿ)⁻¹ 1 nSb_[nⁿ2t]

, (5.10)

where Sbⁿ is a random walk on the integers dened by (5.9). The processXbⁿ is a random walk on the grid Iⁿ. In fact this grid is made in order that Xbⁿ spends the same average time in each of its cells.

Combining remark5.3and theorem3.2we should have (5.1). To sum up this section we write our scheme in the algorithm form. In the next section we will estimate the approximation error of our scheme.

5.5 The algorithm

Note that by construction(Φⁿ)⁻¹ k/n

=xⁿ_k for allk∈Iⁿ. We dene a functionALGOin the next manner:

INPUT DATA: the coecientsaandρ, the starting pointx, the precision ordernand the nal timet. OUTPUT DATA: an approximation in lawXbⁿ ofX at time t.

1. Setxⁿ₀ ←l. whilexⁿ_k ≤r

set xⁿ_k ←p

a(xⁿ_k)ρ(xⁿ_k)(1/n) +xⁿ_k andk←k+ 1. endwhile

2. Compute theαⁿ_k = (1 +β_kⁿ)/2 withβ_kⁿ dened by (5.7).

3. Sety←Φⁿ(x). if(n y− bn yc)<0.5

set s0← bn yc. else

set s0← bn yc+ 1. endif

4. fori= 0 toi=bn²tc=:N ifxⁿ_si∈R\(l, r)

Returnxⁿ_s

i. endif

We have si=k for somek∈I^k. Simulate a realizationB ofBer(αⁿ_k). Then set si←si+B.

endfor 5. Returnxⁿ_sN.

6 Speed of convergence

In this section we will prove the following theorem.

Theorem 6.1 Assume thata, ρ∈Coeff(λ,Λ)for some 0< λ≤Λ<∞. Let be0< T <∞ andX the process generated byL(a, ρ). Forn∈N consider the processXbⁿ starting fromxdened by,

∀t∈[0, T], Xb_tⁿ=ALGO(a, ρ, x, n, t).

For allf ∈W^1,∞₀ (G)∩ C₀(G), allε >0, and allγ∈(0,1/2)there exists a constantC depending onε,γ, T,λ,Λ,G,ka⁰k_∞,kρ⁰k_∞ kfk_∞,kdf /dxk₂,kdf /dxk_∞,sup_i∈I1/(x_i+1−x_i), and the two rst moments ofT(1) such that, fornlarge enough,

sup

(t,x)∈[ε,T]×G¯

E^xf(X_t)−E^xf(Xb_tⁿ)

≤Cn^−γ. We have,

|E^xf(Xt)−E^xf(Xb_tⁿ)| ≤ |E^xf(Xt)−E^xf(X_tⁿ)|+|E^xf(X_tⁿ)−E^xf(Xb_tⁿ)|

=: e1(t, x, n) +e2(t, x, n).

(6.1)

We will estimatee1(t, x, n)by PDEs techniques ande2(t, x, n)by very simple probabilistic techniques.

6.1 Estimate of a weak error

In this subsection we prove the following proposition.

Proposition 6.1 Assumefbelongs toH¹₀(G)∩C0(G). Let beu(t, x)anduⁿ(t, x)respectively the solutions of(P)(a, ρ, f)and(P)(aⁿ, ρⁿ, f), withaⁿ andρⁿ like in Subsection 5.4, Step 1. Then for all ε >0there is a constantC₁ depending onε, T, λ, Λ,G,kfk_∞,kdf /dxk₂,ka⁰k_∞,kρ⁰k_∞, and sup_i∈I1/(x_i+1−x_i) such that forn large enough,

sup

(t,x)∈[ε,T]×G¯

|u(t, x)−uⁿ(t, x)| ≤ C1

√1 n.

As we will see in Proposition6.2, if we hadI ⊂ Iⁿ we could obtain an upper bound for|||u−uⁿ|||_∞,∞

of the formK(ka−aⁿk²_∞+kρ−ρⁿk_∞). But this is not necessary the case (see Remark5.2). However it is possible to modifyaand ρin order to rend us in a situation close to this one, and we will do that to prove Proposition6.1.

Proposition 6.2 Let bef ∈H¹₀(G)∩C0(G). Let bea₁, ρ₁, a₂, ρ₂∈Coeff(λ,Λ), andI1andI2respectively the set of points of discontinuity of a₁ andρ₁ anda₂ andρ₂. Assume I₁⊂ I₂. Let be u₁(t, x) the weak solution of (P)(a1, ρ1, f) and u2(t, x) the weak solution of (P)(a2, ρ2, f). There exists a constant Ce1

depending onT,λ,Λ,G,kfk_∞, andkdf /dxk₂, such that,

|||u₁−u₂|||_∞,∞≤ Ce₁

ka₁−a₂k²_∞+kρ₁−ρ₂k_∞ .

We need a lemma asserting some standard estimates.

Lemma 6.1 i) Let be f ∈ H¹₀(G). Let be u(t, x) the weak solution of (P)(a, ρ, f). Then ∂tu is in L²(0, T; L²(G))and more precisely,

|||∂tu|||_2,2≤ Λ 2

df dx ₂

. (6.2)

ii) Let bef ∈L²(G). Let beu(t, x) the weak solution of(P)(a, ρ, f). Then ^du_dx is inL²(0, T; L²(G)) and more precisely,

du dx _2,2

≤ 1

λkfk₂. (6.3)

Proof. i)Step 1. Assume rst that a and ρ are C^∞(G) and that f is C_c^∞(G) so that u(t, x) is itself C^∞((0, T)×G). Asu(t, x)is a weak solution of (P)(a, ρ, f), andu, Luand∂_tuareC^∞, using∂_tuas a test function and integrating by parts with respect tox, we get

Z T 0

Z

G

|∂tu|²mρ(dx)dt= Z T

0

Z

G

∂tu Lu mρ(dx)dt=−1 2

Z T 0

Z

G

adu dx

d(∂tu) dx dx dt.

Then using Fubini's theorem, interverting the partial derivatives, integrating by parts with respect tot and then again with respect tox, we get

2 Z T

0

Z

G

|∂tu|²mρ(dx)dt= 1 2

Z

G

a

du(0, x) dx

2

dx−1 2 Z

G

a

du(T, x) dx

2

dx, which leads to (6.2).

Step 2. In the general case, with a and ρ in Coeff(λ,Λ), and f in H¹₀(G), we use a regularization argument, Theorem2.2, Step 1, a compactness argument and an integration by parts with respect tot, to exhibit a functionw∈L²(0, T; L²(G))satifying:

∀ϕ∈ C_c^∞((0, T)×G), Z T

0

Z

G

wϕ=− Z T

0

Z

G

u∂tϕ, and |||w|||_2,2≤ Λ 2

df dx ₂

. That is∂tuis inL²(0, T; L²(G))and veries (6.2).

ii) Thanks to point i) and becauseLu(t, .)∈L²(G)we can write Z T

δ

h∂tu, ui_L2(G,m_ρ)dt= Z T

δ

hLu, ui_L2(G,m_ρ)dt, (6.4) for allδ >0. Asuis in C¹([δ, T],L²(G, m_ρ))dtand we have (see [Bre83]),

2h∂tu, ui_L2(G,m_ρ)= d

dtkuk²_L2(G,m_ρ), ∀t∈[δ, T], (6.5) using an integration by part with respect toxin the right hand side of (6.4), and makingδtend to0we get,

λ 2

du dx

2

2,2

≤1 2

ku(0, .)k²_L2(G,m_ρ)− ku(T, .)k²_L2(G,m_ρ)

, which leads to (6.3).

Proof of Proposition 6.2. Step 1. We introduce the following norm on C(0, T; L²(G))∩L²(0, T; H¹₀(G)):

|v|_G,T :=

sup

t∈[0,T]

kv(t, .)k²₂+

dv dx

2

2,2

1/2

. We have the following estimate:

|||v|||_∞,∞≤ K|v|_G,T, (6.6)

where the constantKdepends only onT andG(see [LSU68], II.Ÿ.3 inequality (3) p 74 withr=∞and q=∞).

We setv:=u2−u1. Our goal is now to estimate|v|G,T.

Step 2. Set(L2, D(L2)) =L(a2, ρ2). Elementary computations show that,v(t, x)is a weak solution to

∂_tv(t, x) =L₂v(t, x) +ρ₂(x) 2

d

dx (a₂(x)−a₁(x))du₁(t, x) dx

+ (ρ₂(x)

ρ1(x)−1)∂_tu₁(t, x), that is to say for allϕ∈ C_c^∞((0, T)×G)we have,

−RT 0

R

Gv∂_tϕm_ρ2(dx)dt=−RT 0

R

Ga_2dx^dvdϕ_dxdx dt +RT

0

R

G

ρ₂ 2

d

dx((a2−a1)^du_dx¹) + (^ρ_ρ²

1 −1)∂tu1

ϕmρ₂(dx)dt.

(6.7)

We takev as a test function in (6.7). Using (6.5) with v, integration by parts, ab≤(λ/2)a²+ (2/λ)b² anda1, ρ2∈Coeff(λ,Λ)we nally get,

|v|G,T ≤ κ⁰ κ

Z T 0

Z

G

(a1−a2)²du1

dx 2

dxdt+ Z T

0

Z

G

1 ρ₁ − 1

ρ₂

∂tu1vdxdt

, (6.8)

whereκ= min(1/Λ, λ/4) andκ⁰ = max(1/λ,1).

Step 3. As f ∈ C0(G)and the semigroup (S_t¹)t≥0 and (S_t²)t≥0 generated respectively by (L1, D(L1)) = L(a1, ρ1)and(L2, D(L2))are Feller, we can consider that

∀t∈[0, T], kv(t, .)k_∞≤ ku2(t, .)k_∞+ku1(t, .)k_∞≤2kfk_∞.

Thus|||v|||_∞,∞≤2kfk_∞; besidesf ∈H¹₀(G)thus using point i) of Lemma6.1and Hölder inequality we get,

RT 0

R

G(_ρ¹

1−_ρ¹

2)∂tu1vdxdt ≤ 2

1 ρ₁ −_ρ¹

2

_∞kfk_∞|||∂tu1|||_1,1

≤ _λ²2

pT|G| kρ1−ρ₂k_∞kfk_∞|||∂tu₁|||_2,2.

≤ _λ^Λ2

pT|G| kfk_∞

df dx

₂kρ1−ρ2k_∞.

(6.9)

To nish we have,

Z TZ

(a1−a2)²(du₁

dx)²dxdt≤ ka1−a2k²_∞

du₁ dx

2

,