Numerical error for SDE: Asymptotic expansion and hyperdistributions

C. R. Acad. Sci. Paris, Ser. I 336 (2003) 851–856

Probability Theory

Numerical error for SDE: Asymptotic expansion and hyperdistributions Paul Malliavin a , Anton Thalmaier b a 10, rue Saint Louis en l’Isle, 75004 Paris, France b Université d’Évry, laboratoire d’analyse et probabilité, bd François Mitterrand, 91025 Evry cedex, France

Received 13 March 2003; accepted 1 April 2003 Presented by Paul Malliavin

Abstract The principal part of the error in the Euler scheme for an SDE with smooth coefficients can be expressed as a generalized Watanabe distribution on Wiener space. To cite this article: P. Malliavin, A. Thalmaier, C. R. Acad. Sci. Paris, Ser. I 336 (2003).  2003 Académie des sciences. Published by Éditions scientifiques et médicales Elsevier SAS. All rights reserved. Résumé Erreur du schéma d’Euler pour une EDS et hyperdistribution. La partie principale de l’erreur dans l’intégration par le schéma d’Euler d’une EDS avec des coefficients réguliers est une distribution de Watanabe généralisée. Pour citer cet article : P. Malliavin, A. Thalmaier, C. R. Acad. Sci. Paris, Ser. I 336 (2003).  2003 Académie des sciences. Published by Éditions scientifiques et médicales Elsevier SAS. All rights reserved.

1. Algebraic computation of the error Asymptotic expansion depends upon the choice of a topology for measuring the size of the remainder term; either a norm topology, or on the opposite side, some weak topology may be chosen. The “weakest” among weak topologies are topologies on spaces of distributions: this will be the point of view of this Note. Our main concern is to establish asymptotic expansions “upstairs”, that is on the probability space itself. We shall deal with the Rd -valued SDE dξW (t) =

n


    Ak ξW (t) dW k + A0 ξW (t) dt;

  ξW t 0 = ξ0 ,

(1)

k=1

where W denotes n-dimensional Brownian motion on Wiener space W n and A0 , A1 , . . . , An are smooth bounded vector fields on Rd with bounded derivatives of all orders. E-mail addresses: [email protected] (P. Malliavin), [email protected] (A. Thalmaier). 1631-073X/03/$ – see front matter  2003 Académie des sciences. Published by Éditions scientifiques et médicales Elsevier SAS. All rights reserved. doi:10.1016/S1631-073X(03)00189-4

852

P. Malliavin, A. Thalmaier / C. R. Acad. Sci. Paris, Ser. I 336 (2003) 851–856

Given ε > 0, the Euler scheme of mesh ε is defined by the following recursion formula: n
       ξWε (qε) − ξWε (q − 1)ε = Ak ξWε (q − 1)ε W k (qε) − W k (q − 1)ε k=1

   + A0 ξWε (q − 1)ε ε,

  for t ∈ tε0 , t 0 ,

ξWε (t) = ξ0

where for t
0 we use the notation tε = [t/ε]ε with [a] the largest integer  a. The Euler scheme for all times is then defined as solution of the delayed SDE dξWε (t) =

n


    Ak ξWε (tε ) dW k (t) + A0 ξWε (tε ) dt

for t
t 0 ,

  ξWε | tε0 , t 0 = ξ0 .

(2)

k=1

The remainder term θε (t) := ξWε (t) − ξW (t) satisfies the SDE dθε =

n
          Ak ξW (t) + θε (t) − Ak ξW (t) dW k (t) + A0 ξW (t) + θε (t) − A0 ξW (t) dt + dχ(t),

(3)

k=1

where θε (t 0 ) = 0 and dχ(t) :=

n


ε

Rk (t) dW k (t) + ε R0 (t) dt,

ε

    Rk (t) := Ak ξWε (tε ) − Ak ξWε (t) .

k=1

Let Ak be the d × d matrix defined by differentiating the components of the vector field Ak with respect to the coordinate vector fields; then almost surely the derivative of the solution ξWε (t) of (2) with respect to the initial data ξ0 defines a random flow of diffeomorphisms; its Jacobian is given by the matrix-valued delayed SDE   n
    Wε Wε k Ak ξW (tε ) dW + A0 ξW (tε ) dt , t
t 0 , dJt ←t 0 = Jt ←t 0 k=1 ε JtW ←t 0

= id for t ∈ [tε0 , t 0 ]. The derivatives of ε Rk (t) may be computed in terms of the Jacobian matrix:  ε  Wε     Dτ, ε Rk (t) = 1{τ
t } Ak ξWε (tε ) JtW (4) ←τ (A ) − Ak ξWε (t) Jt ←τ (A ) ;

where

the derivatives Dτ, θε (t) =: u(t) are computed by differentiating (3). We get du −

n


Ak (ξW + θε )u dW k − A0 (ξW + θε )u dt =: dΓ ≡

k=1

n


Γk dW k + Γ0 dt,

(5)

k=1

where Γ0 , Γ1 , . . . , Γn may be computed using (4) and standard computations of derivatives along the stochastic flow to SDE (1). Then, by Itô’s formula, the following version  of the Lagrange formula (variation of constants) is established: in terms of the compensation vector field Z := nk=1 Ak Γk , we have t    (6) Jt 0 ←τ dΓ (τ ) − Z(τ ) dτ , u(t) = Jt ←t 0 t0

where the Itô stochastic integral inside the brackets has to be computed first. We introduce a parameter λ ∈ [0, 1] and define λ θ as the solution of the SDE dλ θ (t) =

n


    Ak ξW (t) + λ θ (t) − Ak ξW (t) dW k (t)

k=1

     + A0 ξW (t) + λ θ (t) − A0 ξW (t) dt + λ dχ(t),

λ

  θ t 0 = 0;

(7)

P. Malliavin, A. Thalmaier / C. R. Acad. Sci. Paris, Ser. I 336 (2003) 851–856 d λ as 0 θ (t) = 0 for all t, denoting dλ θ = λ u, we have θε = λ the following linear SDE for u:

1λ 0

853

u dλ. By differentiating (7) with respect to λ, we get

d λu = d εQ · λ u + dχ,  where d εQ = nk=1 Ak (ξW + λ θ ) dW k + A0 (ξW + λ θ ) dt.

(8)

2. Diagonal-regular functions and diagonal-distributions on Wiener space We consider the diagonal of Rr which is the 1-dimensional subspace γ = {τ1 = · · · = τr }; let Φ r be the quotient r map Rr → Rr /γ . A function U ∈ L2 (Rr ) is said to be diagonal-continuous if EΦ [U 2 ] is a continuous function, r γ p p where EΦ denotes the conditional expectation under Φ r ; the space Dr is the subspace of Dr such that derivatives up to the order r are diagonal-continuous. We define a diagonal Sobolev norm  r


j p/2 p

E Φ |Dτ ,...,τ f |2 ∞ j f γ p = E |f |p + . 1

Dr

j =1

j

L (R /γ )

We remark that solutions of SDE with infinitely bounded differentiable coefficients give rise to functionals γ p γ p belonging to Dr . A diagonal-distribution (see [8]) is a linear form on Dr such that |f, S|  c f γ Drp .

3. Estimation of the error in Sobolev norm of positive order p

We get the following estimates of the remainder term in the Lp -norm, resp. γ Dr -Sobolev norm; see [2,6,7] for related results. Proposition 3.1. For any p ∈ ]1, ∞[ and any integer r > 0, we have   p 1/p



√  √  E sup θε (t)Rd =O ε ; sup θε (t) γ D p = O ε . t ∈[t 0 ,T ]

t ∈[t 0 ,T ]

(9)

r

ε W Proof. Denote by εJ W t ←τ the solution to the linear equation (8) with initial condition J τ ←τ = id; as its coefficients Ak are bounded, we have uniformly with respect to ε:   2p    cp < ∞; E sup εJ W t ←τ t,τ ∈[t 0 ,T ]

2p

by the same reason the compensation vector field Z is bounded in Lp ; we get E[|ε Rk |Rd ] = O(εp ). Using (6), we have t n 
W    ε λ ε W ε k ε W ε u(t) = J t ←t 0 J t 0 ←τ Rk (τ ) dW (τ ) + J t 0 ←τ R0 (τ ) − Z(τ ) dτ . t0

k=1

Consequently Dτ,k θ satisfies an SDE of the same nature as (3); the second part of (9) is verified along the same lines. ✷ Corollary 3.1 (cf. [4]). The Euler scheme converges quasi-surely.

854

P. Malliavin, A. Thalmaier / C. R. Acad. Sci. Paris, Ser. I 336 (2003) 851–856

4. Asymptotic expansion in terms of a diagonal-distribution

λv

We localize the method of Section 3 by using a Taylor expansion at λ = 0; denoting the second derivatives by := (d/dλ) λ u we have by means of (9) θ = ◦u +

1◦ v + o(ε). 2

The question of asymptotic expansion of θ as o(ε) is therefore reduced to the asymptotic expansion of ◦ u, ◦ v which we abbreviate as u, v. We compute u from (8) for λ = 0 which realizes a computation along the path of the original diffusion. In the same way differentiating equation (8) with respect to λ and letting λ = 0, the terms d


R1k (t) :=

,p=1

 ∂ 2 Ak  ξW (t) u (t)up (t),  p ∂ξ ∂ξ

n


  Ak · R1k ξW (t)

(10)

   1 1 R1k (τ ) dW k (τ ) + JtW . 0 ←τ R0 (τ ) − Z (τ )

(11)

Z 1 (t) :=

k=1

appear, and we get t v(t) = JtW ←t 0

dτ

n


JtW 0 ←τ





k=1

t0

Theorem 4.1 (Rate of weak convergence). There exist Rd -valued functions ak , bk , c on Rd , computable in terms γ of the coefficients of (1) and its first fourth derivatives, such that for any f ∈ D3∞−0 ,  1  lim E θε (t) f = ε→0 ε



t
n n
   2    E ak ξW (τ ) Dτ,k;τ, f + bk ξW (τ ) Dτ,k f + c ξW (τ ) f dτ. t0

k,

(12)

k=1

Remark 1. This result should be compared with [5,1,3]. Proof.Assuming that diagonal distributions  S1, . . . ,  S d of third order are given, we consider the formal expression q k γ ∞−0 W q S = k (Jt ←t 0 )k S . As the coefficients of the matrix JtW and as the space of third order 0 ←t belong to D3

diagonal distributions is a module over the algebra D3∞−0 , we deduce that S is a diagonal distribution of third  ˜ := JtW order. Defining u(t) ˜ := JtW 0 ←t u(t) and v(t) 0 ←t v(t), we are reduced to find a diagonal distribution S such that γ

   1  E u(t)f ˜ = f,  S, ε→0 ε

 1  E v(t)f ˜ = 0. ε→0 ε lim

lim

To short-hand the notation, we write W (·) := JtW ˜ = u˜ 1 (t) + u˜ 2 (t) where 0 ←τ (·); then u(t) u˜ 1 (t) :=

t
n  W ε t0



t

Rk (τ ) dW (τ ), k

u˜ 2 (t) :=

k=1

t0

Using the fact that an Itô integral is a divergence we get   E f u˜ 1 (t) =



t
n   (Dτ,k f ) W ε Rk (τ ) dτ ; t0

k=1



W ε

 R0 (τ ) − Z(τ ) dτ.

P. Malliavin, A. Thalmaier / C. R. Acad. Sci. Paris, Ser. I 336 (2003) 851–856

by Itô formula applied to the delayed SDE (2) we have

τ


τ   ε     ε s Rk (τ ) = − Ak ξWε (λ) · As ξWε (τε ) dW (λ) − LAk ξWε (λ) dλ, s

τε

855

(13)

τε

where  ε   p  ∂ 2 Ak  
  ∂Ak   1
 ξWε (λ) + ξWε (λ) . LAk ξWε (λ) := As ξWε (τε ) As ξWε (τε ) A0 ξWε (τε ) 2 ∂ξ ∂ξp ∂ξ s,,p



We want to eliminate the first term in (13). To this end we remark       E (Dτ,r f ) W ε Rr (τ ) = E ENτ (Dτ,r f ) W ε Rr (τ ) ; 

τ
    Nτ ε s E E (Dτ,r f ) Ak ξWε (λ) · As ξWε (τε ) dW (λ) = 0. s

τε

By the Ocone–Karatzas formula, we get n




τ

Nτ

E

Nτ ε

[Dτ,r f ] − E

[Dτ,r f ] =

  2 ENλ Dτ,r;λ,k f dW k (λ)

k=1 τε W (·)

and with analogously as above,   E f u˜ 1 (t)  t

τ 
n n  2    W   
 W ε Dτ,r;λ,k f Ar ξWε (λ) · Ak ξWε (τε ) + dτ dλ (Dτ,k f ) LAk ξWε (λ) =E t0

r,k=1

τε

k=1

which as ε → 0 gives rise to the equivalence  n  t n

     ε     2 E f u˜ 1 (t)  E Dτ,r;τ,k dτ f W (Ar · Ak ) ξW (τ ) + (Dτ,k f ) W 0LAk ξW (τ ) . 2 t0

r,k=1

k=1

Finally t    W  W E f u˜ 2 (t) = E f · R0 − Z dτ t0

ˆ cˆ ˆ b, may be computed as before. There appear integrals along the paths of Dτ,k f , f , and we get existence of a, such that 

t
n n
    2    1  ˜ = E bˆk ξW (τ ) Dτ,k f + cˆ ξW (τ ) f dτ. lim E u(t)f aˆ k ξW (τ ) Dτ,k;τ, f + (14) ε→0 ε t0

k,

k=1

It remains to take care of v: ˜ combining (9) and (11) gives immediately v˜ = O(ε). To get the sharper estimate o(ε) we use that (10) gives R1k as a bilinear expression in u. A typical term is  t t W ∂ 2A W ∂ 2A   k  p k  k p = u u dW E B(τ ) u (τ ) dτ, B := (D h) u; E h τ,k ∂ξ  ∂ξ p ∂ξ  ∂ξ p t0

t0

856

P. Malliavin, A. Thalmaier / C. R. Acad. Sci. Paris, Ser. I 336 (2003) 851–856

applying formula (14) and taking f = B, we get  t    W ∂ 2A k  p  k  u u dW   cεBγD s  cεhγ D 2s uγ D 2s , E h 2 3 2   ∂ξ  ∂ξ p t0

and using (9), with θε (·) replaced by u, we obtain the wanted order O(ε).

✷

References [1] V. Bally, D. Talay, The law of the Euler scheme for stochastic differential equations. I. Convergence rate of the distribution function, Probab. Theory Related Fields 104 (1996) 43–60. [2] E. Gobet, Weak approximation of killed diffusion using Euler schemes, Stochastic Process. Appl. 87 (2000) 167–197. [3] A. Kohatsu-Higa, S. Ogawa, Weak rate of convergence for an Euler scheme of nonlinear SDE’s, Monte Carlo Methods Appl. 3 (1997) 327–345. [4] D. Laurent, Analyse quasi-sure du schéma d’Euler, C. R. Acad. Sci. Paris 315 (1992) 599–602. [5] D. Talay, L. Tubaro, Expansion of the global error for numerical schemes solving stochastic differential equations, Stochastic Anal. Appl. 8 (1990) 483–509. [6] D. Talay, Z. Zheng, Quantiles of the Euler scheme for diffusion processes and financial applications, Math. Finance 13 (2003) 187–199. [7] E. Teman, Analysis of error with Malliavin calculus: application to hedging, Math. Finance 13 (2003) 201–214. [8] S. Watanabe, Analysis of Wiener functionals (Malliavin calculus) and its applications to heat kernels, Ann. Probab. 15 (1987) 1–39.