Measuring Euler–Maruyama Convergence for Geometric Brownian Motion

# How to determine the order of convergence of the Euler-Maruyama method?


# How to determine the order of convergence of the Euler-Maruyama method?












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: 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)
    print(mn)
    err[j]      = abs( mn - solution)

plt.clf();
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

## Answer by Lucy Nowacki (score 2, accepted)

https://quant.stackexchange.com/a/45959

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)))
    rhs=np.log10(errors)
    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)
#Exact
ones = np.ones(M)
Xexact = S0*np.exp((r-0.5*sigma**2)*ones+sigma*W[:, -1])

#np.random.seed(1234)
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.legend(loc='best')
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)
```