To make this simple let us consider the Geometric Brownian Motions.

My questions:

  • 1. How can I show that the Euler-Maruyama Method is convergent using GBM?
  • 2. How can I determine the order of convergence

According to Numerical Analysis theory An SDE solver has order p if the expected value of the error is of pth order in the time step size

Now letting $S_t$ be a GBM with $S_0=1,$ $\mu=0.1$ and $\sigma=0.15$ then $E[S_{10}] =e$

When I run the solver $10,000$ times for different sizes $dt$ I would expect that the difference between my sample mean and true mean, $e$ , will decrease as $dt$ gets smaller. How ever when I run these simulations this is what I get: enter image description here This does not indicate that the error converges to zero as $dt$ goes to zero? Why is that .....

If anyone is interested in the code

import matplotlib.pyplot as plt
import numpy as np

T = 10
mu = 0.1
sigma = 0.15
S0 = 1
ns = 10000

solution = S0*np.exp( mu*T ) 

dt_ = np.array([0.1,0.05,0.01,0.005])
err = np.zeros( len(dt_ )); 
for j in range ( len(dt_ )):
    dt = dt_[j]
    Sn = np.zeros( (ns) ) 
    for i in range(ns):
        N = int(round(T/dt))
        t = np.linspace(0, T, N)
        ex= np.linspace(0, T, N)
        W = np.random.standard_normal(size = N) 
        W = np.cumsum(W)*np.sqrt(dt) ### standard brownian motion ###
        X = (mu-0.5*sigma**2)*t + sigma*W 
        S = S0*np.exp(X) ### geometric brownian motion ###
        Sn[i]= S[-1]

    mn          = np.mean(Sn)
    err[j]      = abs( mn - solution)

plt.loglog(dt_,err, color ="black", label = "Error (abs)");plt.xlabel("dt",fontsize = 20); plt.ylabel("Error (abs)",fontsize = 20)
plt.loglog(dt_,err, 'o', color ="black" , label = "x")
plt.title("loglog-plot",fontsize = 30);
plt.loglog(dt_,0.05*dt_**0.5, linestyle = ":")

The code is mostly just copy-pasted from StackOverflow


1 Answer 1


For GBM you can show it by theory (see Kloeden) or show it empirically as follows. You can observe that EM has strong order of convergence equal to 0.5, whereas Milstein's is 1.0. Not always the latter is superior. E.g. simple EM beats simple Mistein for Heston. In this case is good to use adapted schemes (Giles) for instance, especially dealing with jumps.

In below you have to use inner loop deltaW for each P (Higham). This loop is necessary to see the effect of the change of the time-step for convergence in strong sense. You don't need to apply it for weak error as we need only the mean of the solution in here, however in strong case below this loop is very desirable. This code computes strong errors of the end-points in the sample mth path for the pth step-size. It is useful, e.g. for some EU option pricings, however for strict path-dependant options you have to store arrays with entire paths, i.e. for any t.
To see effect just copy and paste :)

#######+++++++++++++++++++++++++++Lucy Juicy   +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++##############
import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
import numpy.random as npr

def slope(steps, errors):
    steps = DDt.reshape(-1,1)
    A = np.hstack((np.ones((P,1)), np.log10(steps)))
    sol = np.linalg.lstsq(A,rhs, rcond=None)[0]
    q =sol[1]
    resid=np.linalg.norm(np.dot(A,sol) - rhs)
    print('The rate of convergence is: %1.6f and residual = %1.6f' % (q, resid))

def a(X):
    return r*X

def b(X):
    return sigma*X

def bd(X):
    return sigma*np.ones_like(X)

# Problem definition
M =10000 ; T=1; N = 2**9
P = 9             # Number of discretizations
dt = 1.0*T/N
r = 0.07; sigma = 0.1; S0 = 80.0; K=100

dW = np.random.normal(0.0, np.sqrt(dt), (M,N+1))
dW[ : , 0] = 0.00
W = np.cumsum(dW, axis=1)
ones = np.ones(M)
Xexact = S0*np.exp((r-0.5*sigma**2)*ones+sigma*W[:, -1])

Xemerr = np.empty((M,P))
Xmilerr = np.empty((M,P))
DDt =  np.ones(P)
for p in range(P):
    R = 2**p; L = N/R; Dt = R*dt    
    DDt[p] = Dt
    Xem = S0*ones
    Xmil = S0*ones
    Wc = W[:, ::R]
    for j in range(int(L)):
        deltaW = Wc[:, j+1]-Wc[: ,j]
        deltaW2 = Wc[:, j+1]-Wc[: ,j]
        Xem += Dt*a(Xem) + deltaW*b(Xem)
        Xmil += Dt*a(Xmil) + deltaW*b(Xmil) + 0.5*b(Xmil)*bd(Xmil)*(deltaW**2-Dt)
    Xemerr[:,p] = np.abs(Xem - Xexact)
    Xmilerr[:,p] = np.abs(Xmil - Xexact)
Xemerr_tot_strong = np.mean(Xemerr, axis=0)
Xmilerr_tot_strong = np.mean(Xmilerr, axis=0)

slope(DDt, Xemerr_tot_strong)
slope(DDt, Xmilerr_tot_strong)

fig = plt.figure(figsize=(12,6))
ax = fig.add_subplot(111)
ax.loglog(DDt,Xemerr_tot_strong, 'm-', label='EM')
ax.loglog(DDt,Xemerr_tot_strong,linestyle='None', marker='o',markersize='6', markerfacecolor='c',markeredgewidth='1.7', markeredgecolor='m')
ax.loglog(DDt,Xmilerr_tot_strong, 'r-', label='Milstain')
ax.loglog(DDt,Xmilerr_tot_strong,linestyle='None', marker='o',markersize='6', markerfacecolor='g',markeredgewidth='1.7', markeredgecolor='r')
ax.loglog(DDt,DDt, 'c--')
ax.grid(True, which='both', ls='--')
ax.set_xlabel('$\Delta t$'); ax.set_ylabel('$E| $(Price(T)) - Exact Solution(T) |$')
ax.set_title('Strong Error for Geometric Brownian Motion', fontsize=16)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.