Cómo determinar el orden de convergencia de Euler-Maruyama método?

Question

Para hacer esta sencilla consideremos el Geométrico Browniano Movimientos.

Mis preguntas:

• 1. ¿Cómo puedo mostrar que el de Euler-Maruyama Método es convergente el uso de GBM?
• 2. ¿Cómo puedo determinar el orden de convergencia

De acuerdo con el Análisis Numérico de la teoría de Un SDE solver tiene orden p si el valor esperado del error es de pth pedido en el tiempo el tamaño de paso de

Ahora dejando $S_t$ ser una GBM con $S_0=1,$ $\mu=0.1$ y $\sigma=0.15$ entonces $E[S_{10}] =e$

Cuando yo ejecute el solucionador de $10,000$ veces para diferentes tamaños de $dt$ me sería de esperar que la diferencia entre la media muestral y la media real, $e$ , disminuirá a medida que $dt$ se hace más pequeño. Cómo siempre, cuando ejecuto estas simulaciones esto es lo que obtengo: Esto no indica que el error converge a cero como $dt$ va a cero? Por eso es que .....

Si alguien está interesado en el código

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 = ":")

El código es en su mayoría sólo copia-pega de StackOverflow

Answer 1

Para GBM se puede mostrar por la teoría (ver Kloeden) o mostrar empíricamente como sigue. Se puede observar que la EM tiene una fuerte orden de convergencia igual a 0.5, mientras que Milstein es 1.0. No siempre el último es superior. E. g. simple EM beats simple Mistein de Heston. En este caso es buena para un uso adaptado esquemas (Giles), por ejemplo, especialmente relacionados con los saltos.

A continuación tienes que utilizar bucle interno deltaW para cada P (Higham). Este bucle es necesario para ver el efecto del cambio de la hora de paso para la convergencia en sentido fuerte. No es necesario aplicar para los débiles de error de tan sólo tenemos que la media de la solución, sin embargo, en caso fuerte por debajo de este bucle es muy deseable. Este código calcula el fuerte de los errores de los puntos finales en la muestra mth ruta de acceso para la pth del tamaño de paso. Es útil, por ejemplo, para algunos la UE opción pricings, sin embargo, para el estricto ruta que depende de las opciones que usted tiene para almacenar matrices con toda rutas, es decir, para cualquier t.
Para ver el efecto sólo tienes que copiar y pegar :)

#######+++++++++++++++++++++++++++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('