Approximations of solutions to SDEs

Further, due to the properties of the Wiener process we can simulate its values at the
selected points by

Wtn+1 = Wtn + ∆Wn , Wt0 = W0 = 0,

where ∆Wn ∼ N (0, δt) are indep.

Consequently, the stochastic integral
Z T
σs dWs
0

is the limit of the ‘partial sums’

N
X
σtn (Wtn+1 − Wtn ),
n=0

as the mesh of the partition P = (t0 , t1 , . . . , tN ) of [0, T ] approaches zero, i.e., δt → 0.

Approximations of solutions to SDEs
Further, due to the properties of the Wiener process we can simulate its values at the
selected points by

Wtn+1 = Wtn + ∆Wn , Wt0 = W0 = 0,

where ∆Wn ∼ N (0, δt) are indep.

Consequently, the stochastic integral
Z T
σs dWs
0

is the limit of the ‘partial sums’

N
X
σtn (Wtn+1 − Wtn ),
n=0

as the mesh of the partition P = (t0 , t1 , . . . , tN ) of [0, T ] approaches zero, i.e., δt → 0.

Two main questions on approximation

Our rule of thumb

Therefore, if δt is small enough, the above formula is a good approximation for the
stochastic integral.

What is the algorithm to simulate the stochastic integral?

Two main questions on approximation

Our rule of thumb

Therefore, if δt is small enough, the above formula is a good approximation for the
stochastic integral.

What is the algorithm to simulate the stochastic integral?

How to measure the error of the approximation?
Two main questions on approximation

Our rule of thumb

Therefore, if δt is small enough, the above formula is a good approximation for the
stochastic integral.

What is the algorithm to simulate the stochastic integral?

How to measure the error of the approximation?

One simulation ∼ one path ∼ one ‘ω’.

A recap from ODEs

Consider the Cauchy problem

y ′ (t) = f (t, y (t)), y (t0 ) = y0 .

Choose h for the size of every step and set tn = t0 + nh. Set tn+1 = tn + h and consider

yn+1 = yn + hf (tn , yn ).

The value of yn is an approximation of the solution to the ODE at time tn , that is,
yn ≈ y (tn ).
A recap from ODEs

Consider the Cauchy problem

y ′ (t) = f (t, y (t)), y (t0 ) = y0 .

Choose h for the size of every step and set tn = t0 + nh. Set tn+1 = tn + h and consider

yn+1 = yn + hf (tn , yn ).

The value of yn is an approximation of the solution to the ODE at time tn , that is,
yn ≈ y (tn ).

The Euler method is explicit.

Euler’s method ctd

Note that any ODE of order N can be represented as a system of first-order ODEs.
Indeed, given
y (N) (t) = f (t, y (t), y ′ (t), . . . , y (N−1) (t))

set z1 (t) = y (t), z2 (t) = y ′ (t), . . . , zN (t) = y (N−1) (t) to get

     
z1′ (t) y ′ (t) z2 (t)
 ..   ..   .. 
′

z (t) = 
 .  
 =  .  
 =  . 

′ (t) y (N−1) (t) 
zN−1 zN (t)
    

′
zN (t) (N)
y (t) f (t, z1 (t), . . . , zN (t))
The Euler–Maruyama scheme

The framework
We consider the SDE
dXt = a(t, Xt )dt + b(t, Xt )dWt

with the initial condition X0 = x.

The Euler–Maruyama scheme

For the above stochastic differential equation the scheme has the form

Xn+1 = Xn + a(tn , Xn )δt + b(tn , Xn )∆Wn ,

where
∆Wn = Wtn+1 − Wtn .
The Euler–Maruyama scheme
In order to obtain this scheme the following approximation is used
Z tn+1
b(s, Xs ) dWs ≈ b(tn , Xn )∆Wn ,
tn
Z tn+1
a(s, Xs ) ds ≈ a(tn , Xn )δt.
tn
The Euler–Maruyama scheme
In order to obtain this scheme the following approximation is used
Z tn+1
b(s, Xs ) dWs ≈ b(tn , Xn )∆Wn ,
tn
Z tn+1
a(s, Xs ) ds ≈ a(tn , Xn )δt.
tn

Black–Scholes
We know that the solution to

dSt = St µ dt + St σ dWt .

is
St = S0 exp(σWt + (µ − 12 σ 2 )t).
The Euler–Maruyama scheme
In order to obtain this scheme the following approximation is used
Z tn+1
b(s, Xs ) dWs ≈ b(tn , Xn )∆Wn ,
tn
Z tn+1
a(s, Xs ) ds ≈ a(tn , Xn )δt.
tn

Black–Scholes
We know that the solution to

dSt = St µ dt + St σ dWt .

is
St = S0 exp(σWt + (µ − 12 σ 2 )t).

Let us simulate it!

Pseudocode

delta_t = T/N

Let N[0..N-1] be an array of N independent draws from N (0,1)

W[0] = 0
For n = 0, 1, ..., N-1
W[n+1] = W[n] + sqrt (delta_t) * N[n]
S[0] = S0
For n = 0, 1, ..., N-1
S[n+1] = S[n] + S[n] * mu * delta_t + S[n] * sigma * (W[n+1] - W[n])

return S[N]
Approximation errors

We may vary the step size in order to get some idea about the approximation error. We
will measure the performance in two ways. We make M independent runs of the
algorithm and compute

M
1 X σ2
|S0 exp(σW m [N] + (µ − )T − S m [N]| ≈ E ST − S[N]
M 2
m=1

M
σ2 1 X m
|E(S0 exp(σWT + (µ − )T ) − S [N]| ≈ EST − ES[N]
2 M
m=1
Piecewise approximations

When writing algorithms we use the following notation:

tn = n · δt,
Xn is the computed value at tn .

In order to study convergence of numerical schema we introduceXtδt a process that

results from connecting of the points obtained by a numerical scheme with line
segments,
t − tn
Xtδt = Xn +


(Xn+1 − Xn ) t ∈ [tn , tn+1 )
tn+1 − tn
i.e.

δt is the step size used in the algorithm.

Xt is the exact solution.
Numerical stability

A numerical scheme is

strongly convergent, if

lim E(|XT − XTδt |) = 0

δt→0

weakly convergent, if
lim |Eg (XT ) − Eg (XTδt )| = 0
δt→0

for every polynomial g .

Stability

Stability generally means that the numerical solution remains bounded and behaves
similarly to the exact solution as the time step size decreases. Stability can be
categorised into different types, including:

Mean-Square Stability: The numerical solution’s mean-square (second moment)

remains bounded over time.
Asymptotic Stability: The numerical solution converges to a steady state or a
stable equilibrium as time progresses.
Almost Sure Stability: The numerical solution remains bounded with probability
one.
Stability
Stability often depends on the choice of the time step size ∆t. If a scheme is stable for
a particular range of ∆t, it is important to ensure that ∆t is chosen appropriately in
practical computations.

Mean-Square Stability: EM is not always mean-square stable, especially for stiff SDEs.
For certain SDEs, especially those with large coefficients, EM can exhibit numerical
instability unless the time step size ∆t is chosen to be very small.
Almost Sure Stability: EM can also exhibit instability in an almost sure sense,
particularly for stiff SDEs or those with certain types of coefficients.

The Milstein scheme is a more accurate approximation than EM, incorporating an additional
term to account for the derivative of the diffusion coefficient.

Mean-Square Stability: The Milstein scheme is generally more stable than EM. It tends
to exhibit better mean-square stability properties, particularly for SDEs where the
diffusion coefficient g (Xt ) varies significantly.
Almost Sure Stability: The Milstein scheme also tends to be more stable in the almost
Weak/strong convergence of EM

The Euler–Maruyama scheme is both strongly and weakly convergent.

Unfortunately, it is not unconditionally stable (unconditional stability means that
no matter how large ∆t is, the numerical method will not produce solutions that
diverge or exhibit unstable behaviour).

For example the following conditions are sufficient for stability:

Weak/strong convergence of EM

The Euler–Maruyama scheme is both strongly and weakly convergent.

Unfortunately, it is not unconditionally stable (unconditional stability means that
no matter how large ∆t is, the numerical method will not produce solutions that
diverge or exhibit unstable behaviour).

For example the following conditions are sufficient for stability:

functions a and b are C 4 with bounded first derivatives, and

Weak/strong convergence of EM

The Euler–Maruyama scheme is both strongly and weakly convergent.

Unfortunately, it is not unconditionally stable (unconditional stability means that
no matter how large ∆t is, the numerical method will not produce solutions that
diverge or exhibit unstable behaviour).

For example the following conditions are sufficient for stability:

functions a and b are C 4 with bounded first derivatives, and

they do not grow too fast with parameters.
Weak/strong convergence of EM

The Euler–Maruyama scheme is both strongly and weakly convergent.

Unfortunately, it is not unconditionally stable (unconditional stability means that
no matter how large ∆t is, the numerical method will not produce solutions that
diverge or exhibit unstable behaviour).

For example the following conditions are sufficient for stability:

functions a and b are C 4 with bounded first derivatives, and

they do not grow too fast with parameters.

Since the Euler–Maruyama scheme proved to be an approximation of the solution to the

SDE, we may ask how good this approximation is. This gives rise to the notion of
convergence order.
Convergence orders

A numerical scheme is strongly convergent with order γ, whenever

E XT − XTδt ) ⩽ KT · (δt)γ .
Convergence orders

A numerical scheme is strongly convergent with order γ, whenever

E XT − XTδt ) ⩽ KT · (δt)γ .

The constant KT depends on T and the considered SDE. The numerical scheme is
weakly convergent with order γ, whenever

Eg (XT ) − Eg (XTδt ) ⩽ KTg · (δt)γ .

Convergence orders

A numerical scheme is strongly convergent with order γ, whenever

E XT − XTδt ) ⩽ KT · (δt)γ .

The constant KT depends on T and the considered SDE. The numerical scheme is
weakly convergent with order γ, whenever

Eg (XT ) − Eg (XTδt ) ⩽ KTg · (δt)γ .

for every polynomial g . The constant KTg depends on T and g .

Convergence order explained

If a numerical scheme is convergent with order γ and we make the step k times smaller,
then the approximation error will decrease by a factor of k γ .
Thus, the order equal 1 means that if we want to decrease error 100 times, we have to
make the step 100 times smaller. The order equal to 1/2 means that if we want to
decrease error 100 times, we have to make the step 1002 = 10000 times smaller. And
the computation time grows by the same factor.
Euler-Maruyama scheme is

weakly convergent with order 1,

strongly convergent with order 0.5.
Euler-Maruyama scheme is

weakly convergent with order 1,

strongly convergent with order 0.5.

Milstein scheme (to be introduced), is

weakly convergent with order 1,

strongly convergent with order 1.
Euler-Maruyama scheme is

weakly convergent with order 1,

strongly convergent with order 0.5.

Milstein scheme (to be introduced), is

weakly convergent with order 1,

strongly convergent with order 1.

Weak convergence:

says something about the distribution at T only,

useful for pricing of European path-independent options.
Euler-Maruyama scheme is

weakly convergent with order 1,

strongly convergent with order 0.5.

Milstein scheme (to be introduced), is

weakly convergent with order 1,

strongly convergent with order 1.

Weak convergence:

says something about the distribution at T only,

useful for pricing of European path-independent options.

Strong convergence:

pathwise property, should be used if whole paths play a rôle!

good for path-dependent options.
E–M revisited

The Euler–Maruyama scheme may be successfully applied to pricing of

path-independent options (options with payoffs depending only the stock price at the
moment of exercise) i.e. of the form h(ST ).
E–M revisited

The Euler–Maruyama scheme may be successfully applied to pricing of

path-independent options (options with payoffs depending only the stock price at the
moment of exercise) i.e. of the form h(ST ).
This is caused by its poor strong convergence order. If we make the step 100 times
smaller, the approximation improves only by a factor of 10. The step cannot be too
small because of:

(computational errors) errors caused by computations in a computer with the use

of finite precision numbers (with finite number of digits),
computation time.
The Milstein scheme

Still considering
dXt = a(Xt ) dt + b(Xt ) dWt X0 = x

The scheme:

X0 = x
1
Xn+1 = Xn a(Xn )δt + b(Wn )∆Wn + b ′ (Xn )b(Xn )((∆Wn )2 − δt)
2
Derivation

Still, the equation

dXt = a(Xt ) dt + b(Xt ) dWt X0 = x

but this time with some smoothness assumptions on a and b.

Derivation

Still, the equation

dXt = a(Xt ) dt + b(Xt ) dWt X0 = x

but this time with some smoothness assumptions on a and b.

Let ti , ti+1 be two consecutive points in our time discretisation. By Ito’s formula
Z s Z s
′ 1 ′′
f (Xs ) = f (Xti ) + 2
((f (Xu )a(Xu ) + f (Xu )b(Xu ) ) + f ′ (Xu )b(Xu )dW
ti 2 ti

We can apply the Itô formula on the expressions a(Xs ) and b(Xs ), which are the
coefficients in our SDE.
Derivation

Rt Rs
Xti+1 = Xti + ti i+1 (a(Xti + ti (a′ (Xu )a(Xu ) + 12 a′′ (Xu )b 2 (Xu ))du
Rs
+ ti a′ (Xu )b(Xu )dWu )ds
Rt Rs
+ ti i+1 (b(Xti + ti (b ′ (Xu )a(Xu )
+ 12 b ′′ (Xu )b 2 (Xu ))du
Rs
+ ti b ′ (Xu )b(Xu )dWu )dWs .
Derivation

We want to achieve a method which converges strongly of order 1. By using a time

discretisation, the differentials dW and dt are replaced by the corresponding discrete
versions ∆W and δt.
Derivation

We want to achieve a method which converges strongly of order 1. By using a time

discretisation, the differentials dW and dt are replaced by the corresponding discrete
versions ∆W and δt.
We have

δt · δt = O((δt)2 )
∆W · δt = O((δt)3/2 )
∆W · ∆W = O(δt)
Derivation

We want to achieve a method which converges strongly of order 1. By using a time

discretisation, the differentials dW and dt are replaced by the corresponding discrete
versions ∆W and δt.
We have

δt · δt = O((δt)2 )
∆W · δt = O((δt)3/2 )
∆W · ∆W = O(δt)

Therefore: If we are up for a method which converges strongly of order 1, we can

neglect the double integrals from the previous slide, which are of type dWs · ds and
ds · ds.
Derivation

We then obtain

Rt
Xti+1 ≈ Xti + ti i+1 a(Xti )ds
Rt Rs
+ ti i+1 (b(Xti + ti b ′ (Xu )b(Xu )dWu )dWs
Rt Rs
≈ Xti + a(Xti )δt + b(Xti )∆Wi + ti i+1 ti b ′ (Xu )b(Xu )dWu dWs

The first two summands in the equation above are well known from the
Euler–Maruyama scheme. The third one is new. We approximate the third term above
by

Z ti+1 Z s Z ti+1 Z s
b ′ (Xu )b(Xu )dWu dWs ≈ b ′ (Xti )b(Xti ) dWu dWs
ti ti ti ti
Derivation

The integral on the right hand side of the last equality is well known in Continuous
Time Finance. It is
Z ti+1 Z s
′ 1
b (Xti )b(Xti ) dWu dWs = b ′ (Xti )b(Xti )((∆Wi )2 − δt)
ti ti 2
Derivation

The integral on the right hand side of the last equality is well known in Continuous
Time Finance. It is
Z ti+1 Z s
′ 1
b (Xti )b(Xti ) dWu dWs = b ′ (Xti )b(Xti )((∆Wi )2 − δt)
ti ti 2

Substituting this in our previous approximation we finally obtain the Milstein scheme.
Derivation

The integral on the right hand side of the last equality is well known in Continuous
Time Finance. It is
Z ti+1 Z s
′ 1
b (Xti )b(Xti ) dWu dWs = b ′ (Xti )b(Xti )((∆Wi )2 − δt)
ti ti 2

Substituting this in our previous approximation we finally obtain the Milstein scheme.
One can prove that the Milstein scheme converges strongly with order 1 to the solution
of the SDE.