1 Stochastic Taylor Expansion

In this lecture, we discuss the stochastic version of the Taylor expansion to understand how
stochastic integration methods are designed. In addition, we illustrate why the Euler method
is strongly convergent with order 1/2 and is weakly convergent with order 1.

1.1 Ito-Taylor Expansion

First, let us recall how we can obtain a Taylor expansion from an integral representation for
the deterministic case. Consider the autonomous ODE,
d
X (t) = a [X (t)] .
dt
Let f be a function of X (t) , then the evolution of the function f is governed by

d ∂
f [X (t)] = a [X (t)] f [X (t)] (1)
dt ∂X
via the chain rule.

Integral Equation: By defining a linear operator,

∂
L ≡ a (X) ,
∂X
we can rewrite Eq. (1) in terms of the integral equation:
Z t
f [X (t)] = f [X (t0 )] + Lf [X (τ 1 )] dτ 1 (2)
t0

where t0 is the initial time.

Iterations: Iterating Eq. (2) leads to

Z t · Z τ1 ¸
f [X (t)] = f [X (t0 )] + L f [X (t0 )] + Lf [X (τ 2 )] dτ 2 dτ 1
t0 t0
Z t Z t Z τ1
= f [X (t0 )] + Lf [X (t0 )] dτ 1 + L2 f [X (τ 2 )] dτ 2 dτ 1
t0 t t
Z0t Z0τ1
= f [X (t0 )] + Lf [X (t0 )] (t − t0 ) + L2 f [X (τ 2 )] dτ 2 dτ 1
t0 t0

1
 If we iterate once more using Eq. (2) for f (X (τ 2 )) in the integrand of the double integral,
· Z τ2 ¸
2 2
L f [X (τ 2 )] = L f [X (t0 )] + Lf [X (τ 3 )] dτ 3
t0
Z τ2
2 2
= L f [X (t0 )] + L Lf [X (τ 3 )] dτ 3
t0
| {z }
O(L3 )

then the contribution of the first term above, which is a constant, to the double integral is
Z t Z τ1 Z t Z τ1
2 2
L f [X (t0 )] dτ 2 dτ 1 = L f [X (t0 )] dτ 2 dτ 1
t0 t0 t0 t0
1 2
= L f [X (t0 )] (t − t0 )2
2
Putting all these together, we have
1 ¡ ¢
f [X (t)] = f [X (t0 )] + Lf [X (t0 )] (t − t0 ) + L2 f [X (t0 )] (t − t0 )2 + O L3
2
which is precisely the Taylor expansion we are familiar with.

1.1.1 Ito-Taylor Expansion

Now we discuss a similar expansion for the stochastic differential equation:

dX (t) = a [X (t)] dt + b [X (t)] dW (t) . (3)
Again, for simplicity, we consider the autonomous case, i.e., a = a [X (t)] , b = b [X (t)] and
they do not depend on time explicitly.
The Ito lemma leads to
½ ¾
∂ 1 2 ∂2 ∂
df [X (t)] = a f [X (t)] + b [X (t)] f [X (t)] dt + b [X (t)] f [X (t)] dW (t)
∂X 2 ∂X 2 ∂X
(4)
Defining
∂ 1 ∂2
L0 ≡ a + b2 [X]
∂X 2 ∂X 2
∂
L1 ≡ b [X]
∂X
then Eq. (4) becomes
df [X (t)] = L0 f [X (t)] dt + L1 f [X (t)] dW (t)
i.e., Z Z
t t
0
f [X (t)] = f [X (t0 )] + L f [X (s)] ds + L1 f [X (s)] dW (s) (5)
t0 t0
In equation (5), if we choose the following f (x):

2
 1. Choose f (x) = x, then Eq. (5) becomes
Z t Z t
X (t) = X (t0 ) + a [X (s)] ds + b [X (s)] dW (s) (6)
t0 t0

2. Choose f (x) = a (x) , then Eq. (5) becomes

Z t Z t
0
a [X (t)] = a [X (t0 )] + L a [X (s)] ds + L1 a [X (s)] dW (s) (7)
t0 t0

3. Similarly, choose f (x) = b (x) , then Eq. (5) becomes

Z t Z t
0
b [X (t)] = b [X (t0 )] + L b [X (s)] ds + L1 b [X (s)] dW (s) (8)
t0 t0

Substituting Eqs. (7) and (8) into Eq. (6) leads to

Z t½ Z s1 Z s1 ¾
0 1
X (t) = X (t0 ) + a [X (t0 )] + L a [X (s2 )] ds2 + L a [X (s2 )] dW (s2 ) ds1
t0 t0 t0
(9)
Z t½ Z s1 Z s1 ¾
+ b [X (t0 )] + L0 b [X (s2 )] ds2 + 1
L b [X (s2 )] dW (s2 ) dW (s1 )
t0 t0 t0

Note that
∂ 1 ∂2 1
L0 a = a a + b2 [X] 2
a ≡ aa0 + b2 a00
∂X 2 ∂X 2
1
L0 b = ab0 + b2 b00
2
∂
L1 a = b a = ba0
∂X
L1 b = bb0

Separating the constant terms out in the integrand in Eq. (9) from the remaining terms,
which are the double integral terms:
Z t Z s1 Z t Z s1
0
R ≡ L a [X (s2 )] ds2 ds1 + L1 a [X (s2 )] dW (s2 ) ds1
t0 t0 t0 t0
Z t Z s1 Z t Z s1
0
+ L b [X (s2 )] ds2 dW (s1 ) + L1 b [X (s2 )] dW (s2 ) dW (s1 ) ,
t0 t0 t0 t0

leads to Z t Z t
X (t) = X (t0 ) + a [X (t0 )] ds1 + b [X (t0 )] dW (s1 ) + R (10)
t0 t0

3
 Note that the essence of the method is to use the substitution repeatedly to obtain
constant integrands in higher and higher order terms. For example, the last term in the
remainder, R, (which is of the lowest order in ∆t in R if ∆t ≡ t − t0 is small) is
Z t Z s1
L1 b [X (s2 )] dW (s2 ) dW (s1 )
t t
Z 0t Z 0s1 ½ Z s2 Z s2 ¾
1 0 1 1 1
= L b [X (t0 )] + L L b [X (s3 )] ds3 + L L b [X (s3 )] dW (s3 ) dW (s2 ) dW (s1 )
t0 t0 t0 t0

by selecting f = L1 b in Eq. (5). The first term in the last line of the above equation is
Z t Z s1
0
b [X (t0 )] b [X (t0 )] dW (s2 ) dW (s1 )
t0 t0

by noting L1 b = bb0 . Therefore, Eq. (10) becomes

Z t Z t
X (t) = X (t0 ) + a [X (t0 )] ds1 + b [X (t0 )] dW (s1 )
t0 t0
Z t Z s1
0
+b [X (t0 )] b [X (t0 )] dW (s2 ) dW (s1 ) + R̃ (11)
t0 t0

where R̃ is a new remainder. Eq. (11) is the Ito-Taylor expansion for the process (3).

Note that the double integral in Eq. (11) is evaluated to be

Z t Z s1
1 1
dW (s2 ) dW (s1 ) = [W (t) − W (t0 )]2 − (t − t0 ) ,
t0 t0 2 2

which can be shown as follows

Z tZ s1 Z t
dW (s2 ) dW (s1 ) = [W (s1 ) − W (t0 )] dW (s1 )
t0 t0 t
Z 0t Z t
= W (s1 ) dW (s1 ) − W (t0 ) dW (s1 )
t0 t0
Z Z t
1 t£ 2
¤
= dW (s1 ) − dt − W (t0 ) dW (s1 )
2 t0 t0
¡ ¢
∵ dW 2 = 2W dW + dt (Ito Lemma)
1£ 2 ¤
= W (t) − W 2 (t0 ) − (t − t0 ) − W (t0 ) [W (t) − W (t0 )]
2
1 1
= [W (t) − W (t0 )]2 − (t − t0 ) (12)
2 2

4
1.1.2 Numerical Integration Schemes

Once we have the Ito-Taylor expansion, we can construct numerical integration schemes for
the process (3). For the interval [ti , ti+1 ] , by choosing

t0 = ti , t = ti+1 ,
∆t = ti+1 − ti
∆Wi = W (ti+1 ) − W (ti ) ,

combining Eq. (11) and Eq. (12) yields

1 £ ¤
X (ti+1 ) = X (ti ) + a [X (ti )] ∆t + b [X (ti )] ∆Wi + b [X (ti )] b0 [X (ti )] (∆Wi )2 − ∆t + R̃.
2
(13)

1. Keeping the first three terms in Eq. (13) gives us the explicit Euler method:
h i h i
X̂i+1 = X̂i + a X̂i ∆t + b X̂i ∆Wi

2. Keeping all terms of O (∆t) gives us the Milstein method:

h i h i 1 h i h i£ ¤
X̂i+1 = X̂i + a X̂i ∆t + b X̂i ∆Wi + b X̂i b0 X̂i (∆Wi )2 − ∆t (14)
2

Note that the Milstein needs to compute the derivative b0 for the last term.

Runge-Kutta Methods Note that the Milstein method requires to compute the deriva-
tive of b. Some times it is computationally costly to compute derivatives. We can construct
Runge-Kutta schemes to avoid this, as for the deterministic case.

For ∆X = a∆t + b∆W, we have

¡ ¢
b (X + ∆X) − b (X) = b0 (X) ∆X + O (∆X)2 ,

therefore,
¡ ¢
b (X + ∆X) − b (X) = b0 (X) [a (X) ∆t + b (X) ∆W ] + O (∆t) ∵ (∆X)2 ∼ ∆t
= b0 (X) b (X) ∆W + O (∆t)

∴
1
b0 (Xi ) b (Xi ) ∼ √ [b (Xi + a (Xi ) ∆t + b (Xi ) ∆Wi ) − b (Xi )] + O (∆t)1/2
∆t
1 h ³ √ ´ i
∼ √ b Xi + a (Xi ) ∆t + b (Xi ) ∆t − b (Xi ) + O (∆t)1/2
∆t

5
Thus, we have the following RK method
 √

 X̃i = Xi + a (Xi ) ∆t + b (Xi ) ∆t


 1 h ³ ´ i£ ¤
 Xi+1 = Xi + a (Xi ) ∆t + b (Xi ) ∆Wi + √ b X̃i − b (Xi ) (∆Wi )2 − ∆t
2 ∆t
which has order-one strong convergence. (Question: how would you demonstrate this con-
vergence of order one numerically?)

Higher order RK schemes can be constructed similarly.

Clearly, if b =constant, the RK method above reduces to the Euler scheme.

1.2 Weak and Strong Convergence for Euler Scheme

Now we turn to getting some intuitive feeling why the Euler scheme has strong order 1/2
and weak order 1. For simplicity, we analyze the convergence for the geometric Brownian
motion
dS (t) = µS (t) dt + σS (t) dW (t)
where µ and σ are constant. Recall that the exact solution for this process is

S (t) = S (0) e(µ− 2 σ )t+σW (t) .

1 2
(15)
The Euler scheme for this process is
Ŝ (ti+1 ) = S (ti ) + µŜ (ti ) ∆t + σ Ŝ (ti ) ∆W (ti )

i.e., Ŝ (ti+1 ) = Ŝ (ti ) [1 + µ∆t + σ∆W (ti )]

√
where ∆W (ti ) = ∆tZi , Zi ∼ N (0, 1) . Then, obviously, the numerical solution at t = tk is
Y
k−1
Ŝ (tk ) = S (0) [1 + µ∆t + σ∆W (ti )]
i=0

1.2.1 The Error for Strong Convergence

For t0 = 0, t1 = Tn , t2 = 2T
n
,··· , tn = T, by definition, the error in strong convergence is

¯ ¯
¯ ¯
E ¯Ŝ (T ) − S (T )¯
¯n−1 ¯
¯Y ¯
¯ (µ− 12 σ 2 )T +σW (T ) ¯
= S (0) E ¯ [1 + µ∆t + σ∆W (ti )] − e ¯ (16)
¯ ¯
i=0

6
where Ŝ (T ) is numerically evaluated using the Euler method and S (T ) is the exact solution
(Eq.(15)).

To O (∆t)2 , Taylor expansion leads to

·µ ¶ ¸
(µ− 12 σ2 )∆t+σW (ti ) 1 2
e = 1 + µ − σ ∆t + σW (ti )
2
·µ ¶ ¸2
1 1 2
+ µ − σ ∆t + σW (ti )
2 2
·µ ¶ ¸3
1 1 2
+ µ − σ ∆t + σW (ti )
6 2
+ · · ·µ
· ¶ ¸
1 2
= 1 + µ − σ ∆t + σW (ti )
2
µ ¶
1 2 1 ¡ ¢
+ µ − σ σ∆t∆W (ti ) + σ 2 [∆W (ti )]2 N.B. (∆W )2 = ∆t
2 2
1 3
+ σ [∆W (ti )]3 + O (∆t)2
6
1
= 1 + µ∆t + σW (ti ) + µ0 σ∆t∆W (ti ) + σ 3 [∆W (ti )]3 + O (∆t)2
6
¡ ¢
where µ0 ≡ µ − 12 σ 2 . Therefore,

1 + µ∆t + σW (ti )
1
= e(µ− 2 σ )∆t+σW (ti ) − µ0 σ∆t∆W (ti ) − σ 3 [∆W (ti )]3 − O (∆t)2
1 2

6
∴
Y
n−1 Y·
n−1
1
¸
e(µ− 2 σ )∆t+σW (ti ) − µ0 σ∆t∆W (ti ) − σ 3 [∆W (ti )]3
1 2
[1 + µ∆t + σW (ti )] =
i=0 i=0
6

Y
n−1
[1 + µ∆t + σW (ti )] = e(µ− 2 σ )T +σW (T ) + nO (∆t∆W ) + nO (∆W )3 + nO (∆t)2 (17)
1 2
∴
i=0

7
Hence, the error (16) is
¯ ¯
¯ ¯
E ¯Ŝ (T ) − S (T )¯
¯ ¯
= E ¯nO (∆t∆W ) + nO (∆W )3 + nO (∆t)2 ¯
¯ ¯ µ ¶
¯T T T ¯ T
= E ¯¯ O (∆t∆W ) + O (∆W ) + 3
O (∆t) ¯¯
2
∵ n=
∆t ∆t ∆t ∆t
¯ ¯
¯ ¯
¯ ¯
¯ 1 1 1 ¯
¯ ¯
= T E ¯ O (∆t∆W ) + O (∆W )3 + O (∆t)2 ¯
¯|∆t {z } |∆t {z } |∆t {z }¯¯
¯
¯ O(∆t)1/2 O(∆t)1/2 O(∆t) ¯
= O (∆t)1/2

This demonstrates that the Euler scheme has strong convergence of order 1.

1.2.2 The Error for Weak Convergence

For weak convergence, we need to specify the function f in the error estimate:
¯ h ³ ´i ¯
¯ ¯
¯E f Ŝ (T ) − E [f (S (T ))]¯ ≤ β∆tq

For the convergence of expectation, i.e, the first moment, we choose

f (x) = x

For the geometric Brownian motion, we thus have the “weak” error
¯ h i ¯
¯ ¯
¯E Ŝ (T ) − E [S (T )]¯

¯ "n−1 # ¯
¯ Y h i¯
¯ ¯
[1 + µ∆t + σW (ti )] − E e(µ− 2 σ )T +σW (T ) ¯
1 2
= S (0) ¯E
¯ ¯
i=0

¯ h i h i¯
¯ (µ− 12 σ2 )T +σW (T ) 3 2 (µ− 12 σ2 )T +σW (T ) ¯
= S (0) ¯E e + nO (∆t∆W ) + nO (∆W ) + nO (∆t) − E e ¯

Y
n−1
(∵ [1 + µ∆t + σW (ti )]
i=0

= e(µ− 2 σ )T +σW (T ) + nO (∆t∆W ) + nO (∆W )3 + nO (∆t)2 )

1 2

¯ £ ¤¯
= S (0) ¯E nO (∆t∆W ) + nO (∆W )3 + nO (∆t)2 ¯

8
Note that £ ¤
E [O (∆t∆W )] = 0, E O (∆W )3 = 0
this leads to the order of weak convergence for Euler scheme:
¯ h i ¯ ¯ £
¯ ¯ 2 ¤¯
¯E Ŝ (T ) − E [S (T )]¯ = S (0) ¯E nO (∆t) ¯
= S (0) nO (∆t)2
T
= O (∆t)2
∆t
= O (∆t)
which has order one.

The main difference between the strong and weak convergence lies in whether the odd
powered terms, such as ∆W , ∆W 3 , can be eliminated:
¯ ¯
¯ ¯
For strong convergence: E [|O (∆W ) + h.o.t.|] = E ¯O (∆t)1/2 + h.o.t.¯
For weak convergence : |E [O (∆W ) + h.o.t.]| = |0 + h.o.t.|
where h.o.t. is “higher order terms”.

Finally, we try to get a little more intuition of why the Euler scheme has order one of
strong convergence for the process
dS = µdt + σdW
when σ is deterministic. In this case, the Euler scheme is
Ŝ (ti+1 ) = Ŝ (ti ) + µ (ti ) ∆t + σ (ti ) ∆W (ti )
therefore at tn = T,
X
n−1 X
n−1
Ŝ (T ) = S (0) + µ (ti ) ∆t + σ (ti ) ∆W (ti ) .
i=0 i=0

Since the exact solution is

Z T Z T
S (T ) = S (0) + µ (t) dt + σ (t) dW (t) ,
0 0

we have
X
n−1 Z T
Ŝ (T ) − S (T ) = µ (ti ) ∆t − µ (t) dt
i=0 0

X
n−1 Z T
+ σ (ti ) ∆W (ti ) − σ (t) dW (t)
|i=0
{z } |0 {z }
≡B
≡A

9
the first line of which clearly has the O (∆t) error. Now we analyze the error for the second
line in the above equation. Define

X
n−1
A ≡ σ (ti ) ∆W (ti ) ,
i=0
Z T
B ≡ σ (t) dW (t)
0

Note that

1. A is a Gaussian (since a linear superposition of Gaussian random variables is a Gaussian

random variable) with the expectation
à n−1 !
X
EA = E σ (ti ) ∆W (ti ) = 0
i=0

and variance
X
n−1
VarA = σ 2 (ti ) ∆t. (18)
i=0

Therefore, the Gaussian random variable A can be expressed as

v
u n−1
W (T ) u X
A= √ t σ 2 (ti ) ∆t
T i=0

since the right-hand side of the equation above is a Gaussian random variable with
mean zero and variance the same as Eq. (18).

2. B is also a Gaussian random variable with mean and variance:

µZ T ¶
EB = E σ (t) dW (t) = 0
0
Z T
VarB = σ 2 (t) dt
0

therefore, B can be expressed as

s
Z T
W (T )
B= √ σ 2 (t) dt.
T 0

Hence, v 
u n−1 s
uX Z T
W (T ) t
A−B = √ σ 2 (ti ) ∆t − σ 2 (t) dt .
T i=0 0

10
Since Z
X
n−1 T
2
σ (ti ) ∆t − σ 2 (t) dt = O (∆t) ,
i=0 0

X
n−1 Z T
2
σ (ti ) ∆t = σ 2 (t) dt + O (∆t) .
i=0 0
RT
∴ By defining Σ ≡ 0
σ 2 (t) dt, we have

W (T ) hp √ i
A−B = √ Σ + O (∆t) − Σ
T " #
r
W (T ) √ O (∆t) √
= √ Σ 1+ − Σ
T Σ
· µ ¶ ¸
W (T ) √ 1 O (∆t) 2
√
= √ Σ 1+ + O (∆t) − Σ
T 2 Σ
W (T ) 1 O (∆t)
= √ √
T 2 Σ
Therefore,

1 W (T ) O (∆t)
Ŝ (T ) − S (T ) = O (∆t) + √ √
2 T Σ
= O (∆t)

i.e., the Euler scheme has strong convergence of order one for deterministic σ (t) .

11