Numerical Methods for SDEs driven

by Fractional Brownian Motion:

Exact Convergence Rates

A. Neuenkirch
Goethe-Universität Frankfurt am Main

Joint work with

I. Nourdin (Paris VI), S. Tindel (Nancy I) and J. Unterberger (Nancy I)

1
 SDEs driven by FBM

m
X (j)
(SDE) dXt = a(Xt )dt + σ (j) (Xt )dBt , t≥0
j=1

X0 = x0 ∈ Rd

where

• a : Rd → Rd drift vector

• σ = (σ (1) , . . . , σ (m) ) with σ (j) : Rd → Rd diffusion matrix

• B = (B (1) , . . . , B (m) )′ m-dimensional fractional Brownian motion with Hurst

parameter H ∈ (0, 1), i.e.
B (1) , . . . , B (m) independent centered Gaussian processes with continuous sample
paths and
(j)
E |Bt − Bs(j) |2 = |t − s|2H , s, t ≥ 0

2
 SDEs driven by FBM

m
X (j)
(SDE) dXt = a(Xt )dt + σ (j) (Xt )dBt , t≥0
j=1

X0 = x0 ∈ Rd

Here: Pathwise approach (Coutin, Qian (2002))

Existence and uniqueness under standard assumptions on a and σ for H > 1/4 using
rough paths theory (Lyons 1994, 1998)

Lin (1995); Klingenhöfer, Zähle (1999); Mikosch, Norvaiša (2000); Nualart, Răşcanu (2002); Errami,
Russo (2003); Gubinelli (2004); Zähle (2005); Hu, Nualart (2006); Nourdin, Simon (2007); ...

H = 1/2: classical Stratonovich SDE

Other approach: Wick-Itô resp. Skorohod integral equation

Biagini et al. (2002); Mishura (2004); Nourdin, Tudor (2006);...

3
 dXt = −0.5Xt dt + 2Xt dBt , X0 = 1

H = 0.5 H = 0.7

Blue: (Bt )t∈[0,1] Red: (Xt )t∈[0,1]

4
 Outline of the Talk

Problem: Approximation of X1 by

φn (B1/n , B2/n , . . . , B1 )

with φn: Rm,n → Rd , measurable

` ´
2 1/2
Mean square error: e(φn ) = E kX1 − φn (B1/n , B2/n , . . . , B1 )k

In the following:

• simple φn

• optimal φn

for

• one-dimensional case (m = d = 1) for H > 1/2

• fractional Lévy area (m = 2, d = 3, non-commutative noise) for H > 1/4

5
 Outline of the Talk

Problem: Approximation of X1 by

φn (B1/n , B2/n , . . . , B1 )

with φn: Rm,n → Rd , measurable

` ´
2 1/2
Mean square error: e(φn ) = E kX1 − φn (B1/n , B2/n , . . . , B1 )k

Remark
Exact simulation of fractional noise (Coeurjolly (2000)):

• non-equidistant: Choleski method, cost O(n2 )

• equidistant: Methods of Wood-Chan and Davis-Harte, cost O(n log(n))

(Input: n iid N (0, 1) random numbers)

6
 One-dimensional Case

(SDE) dXt = a(Xt )dt + σ(Xt )dBt , t ∈ [0, 1]

X0 = x0 ∈ R

Here: Doss-Sussmann representation

Xt = ϕ(At , Bt ), t ∈ [0, 1]

where
∂ϕ
• ∂y
(x, y)
= σ(ϕ(x, y)), ϕ(x, 0) = x
RB
• A′t = exp(− 0 t σ ′ (ϕ(At , s)) ds) a(ϕ(At , Bt )), A0 = x0

Assumptions:

(A1) H > 1/2

(A2) a ∈ C 2 , σ ∈ C 3 with bounded derivatives
(A3) a, σ bounded, σ elliptic, i.e. inf x∈R |σ(x)| > 0

7
 Euler Scheme

1
Notation: ∆= n
, tk = k · ∆, ∆k B = Btk+1 − Btk

(n)
X0 = x0 ,
(n) (n) (n) (n)
X tk+1 = X tk + a(X tk )∆ + σ(X tk )∆k B, k = 0, . . . , n − 1

Theorem N, Nourdin (2007)

˛Z ˛2 !1/2
(n) 1 ˛ 1 ˛
e(X 1 ) = E ˛˛ σ (Xs )Ds X1 ds˛˛
′
· n−2H+1 + o(n−2H+1 )
2 0

with (Ds X1 )s∈[0,1] Malliavin derivative of X1

Remark Under (A1), (A2):

Z 1
(n) a.s. 1
n2H−1 · (X1 − X 1 ) −→ σ ′ (Xs )Ds X1 ds
2 0
Asymptotic error distribution (N, Nourdin (2007))

8
 Euler Scheme

How to analyze the error of the Euler scheme?

Direct method not possible, since isometry for integrals with respect to fBm involves
Malliavin derivative of the integrand: equations / inequalities not closable!

Use Doss-Sussmann representation:

(1) transform (SDE) to (RODE) At = ϕ(Xt , −Bt )

(n) (n)
(2) transform X to corresponding approximation A , i.e.
(n) (n)
Atk = ϕ(X tk , −Btk ), k = 0, 1, . . . , n
(n)
(3) analyze A1 − A1

(4) transform back

9
 Optimal Approximation

Recall:

• Approximation of X1 by φn (B1/n , B2/n , . . . , B1 )

` ´1/2
• Mean square error: e(φn ) = E kX1 − φn (B1/n , B2/n , . . . , B1 )k2

Optimal φn ? Clearly:

φ∗n (x1 , x2 , . . . , xn ) = E(X1 | B1/n = x1 , B2/n = x2 , . . . , B1 = xn )

∗
Rate of convergence of φn ?

Remark
Known for H = 1/2:
„Z 1 «1/2
1
e(φ∗n ) = √ E |Yt |2 dt · n−1 + o(n−1 )
12 0

where (Yt )t∈[0,1] random weight function

Cameron, Clark (1980); Newton (1986,...); Cambanis, Hu (1996); Müller-Gronbach (2004)

10
 Optimal Approximation

Theorem N (2008)
′
(1) If a σ − aσ ′ = 0, then
e(φ∗n ) = 0
(2) Let
′ ′
R1 ′
R1
Yt = (a σ − aσ )(Xt ) exp( t
a (Xτ ) dτ + t
σ ′ (Xτ ) dBτ ), t ∈ [0, 1]
If Z 1
(ND) | E Yt | dt > 0,
0
then
e(φ∗n ) ≥ C(a, σ, x0 , H) · n−H−1/2

Remarks

• Y : time derivative (in mean square sense) of Malliavin derivative (Dt X1 )t∈[0,1]
• Minimal error for non-equidistant discretizations: In general also
inf{e( E( X1 | Bt1 , . . . , Btn ) ) : t1 , . . . , tn ∈ [0, 1]} ≥ C · n−H−1/2

11
 Optimal Approximation

Proof

(1) Doss-Sussmann representation

a(x0 )
a′ σ − aσ ′ = 0 =⇒ a/σ = const. =⇒ A′t = σ(At ) σ(x0)
, A0 = x0

Thus X1 = ϕ(A1 , B1 ) = ϕ̃(B1 )

(2) Linearization via fractional Wiener chaos decomposition

R1
X1 = E X1 + 0
E Dt X1 dBt + . . .
P
φ∗n (B1/n , . . . , B1 ) = E φ∗n (B1/n , . . . , B1 ) + n
i=1 ai Bi/n + . . .

Thus (by orthogonality of the decomposition)

R1 Pn 2 1/2
e(φ∗n ) ≥ (E | 0
E Dt X1 dBt − i=1 ai B i/n | )

Linearized problem:
Embeddings of reproducing kernel Hilbert spaces; results of Stein (1995), Ritter (2000)

12
 An Optimal Implementable Method

Notation: ∆ = 1/n, tk = k · ∆, ∆k B = Btk+1 − Btk

b (n) = x0
X 0

X bt(n) )∆ + σ(X
bt(n) + a(X
bt(n) = X bt(n) )∆k B + 1 σσ ′ (X
bt(n) )(∆k B)2
k+1 k k k
2 k

1 bt(n) )∆k B∆ + 1 aa′ (X bt(n) )∆2

+ (a′ σ + aσ ′ )(X
2 k
2 k

1 bt(n) )(∆k B)3 ,

+ (σ 2 σ ′′ + σ(σ ′ )2 )(X k = 0, . . . , n − 1
6 k

Remark

For H = 1/2:
(i) McShane’s method (1974) for Stratonovich SDEs

(ii) Müller-Gronbach (2004): truncated Wagner-Platen method for Itô SDEs

13
 An Optimal Implementable Method

Theorem N (2008)
„Z «1/2
(n)
p 1
b ) =
e(X |ζ(−2H)| E |Yt |2 dt · n−H−1/2 + o(n−H−1/2 )
1
0

where
′ ′
R1 ′
R1
Yt = (a σ − aσ )(Xt ) exp( t
a (Xτ ) dτ + t
σ ′ (Xτ ) dBτ ), t ∈ [0, 1]

Corollary If
Z 1
(ND) | E Yt | dt > 0,
0
b
then X
(n)
optimal, i.e. same convergence rate as conditional expectation!

14
 Multi-dimensional Case

m
X (j)
(SDE) dXt = a(Xt )dt + σ (j) (Xt )dBt , t≥0
j=1

X0 = x0 ∈ Rd

Commutative Noise:

D(i) σ (j) = D(j) σ (i) , i, j = 1, . . . , m

where
d
X (i) ∂
D(i) = σk
k=1
∂xk

In this case: Doss-Sussmann representation Xt = ϕ(At , Bt )

SDE with commutative noise “ ⇐⇒ “ one-dimensional SDE

15
 Fractional Lévy Area

0 1 0 1
1 0
B C B C
(fLA) dXt = B 0 C dBt(1) + B 1 C dBt(2) , t ∈ [0, 1]
@ A @ A
(2) (1)
− 21 Xt 1
2
X t

X0 = 0

Prototype for SDEs with non-commutative noise

Z t Z t
(1) (1) (2) (2) (3) 1 1
Xt = Bt , Xt = Bt , Xt = Bs(1) dBs(2) − Bs(2) dBs(1)
2 0 2 0

Approximation of (fLA): quadrature problem

Assumption: H > 1/4

16
 Euler Scheme

n−1 n−1
(3)(n) 1 X (1) (2) (2) 1 X (2) (1) (1)
X1 = Bi/n (B(i+1)/n − Bi/n ) − Bi/n (B(i+1)/n − Bi/n )
2 i=0 2 i=0

Theorem N, Tindel, Unterberger (2008)

8
>
> α1 (H) · n
−2H+1/2
+ o(n−2H+1/2 ) for H ∈ (1/4, 3/4)
< p p
(n) 9
e(X 1 ) = · log(n)n −1
+ o( log(n)n−1 ) for H = 3/4
>
>
128
:
α2 (H) · n−1 + o(n−1 ) for H ∈ (3/4, 1)

Remarks

• Euler scheme for general multi-dimensional SDE (Mishura, Shevchenko (2008)):

−2H+1+ε
pathwise convergence rate n for H > 1/2

• No ”diagonal” noise in (fLA)

• Critical point H = 3/4: compare with behavior of quadratic variation of fBm
Guyon, Leon (1989); Istas, Lang (1997); Corcuera et al. (2006);...

17
 Midpoint Scheme

n−1
b (3)(2n) 1 X (1) (2) (2)
X 1 = B(i+1/2)/n (B(i+1)/n − Bi/n )
2 i=0
n−1
1 X (2) (1) (1)
− B(i+1/2)/n (B(i+1)/n − Bi/n )
2 i=0

Theorem N, Tindel, Unterberger (2008)

b (2n)
e(X 1 ) = α3 (H) · n−2H+1/2 + o(n−2H+1/2 )

Questions

• Asymptotic error distributions for Euler- and Midpoint scheme?

• Midpoint scheme optimal, i.e. e(E(X1 | B1/n , . . . , B1 )) ≥ C(H) · n−2H+1/2
for (fLA)?
1
Known for H = 1/2: e(E(X1 | B1/n , . . . , B1 )) = √
8
· n−1/2
Cameron, Clark (1980); ...

18
 Summary

Exact convergence rates in the one-dimensional case for H > 1/2

• Euler scheme: n−2H+1
• McShane’s method: n−H−1/2
• Conditional expectation:
n−∞ if drift and diffusion commute
n−H−1/2 under a non-degeneracy condition on the Malliavin derivative

Exact convergence rates for the fractional Lévy area

8
>
> n −2H+1/2
if H ∈ (1/4, 3/4)
< p
• Euler scheme: log(n)n−1 if H = 3/4
>
>
:
n−1 if H ∈ (3/4, 1)

• Midpoint scheme: n−2H+1/2

• Conjecture for conditional expectation: n−2H+1/2

19