Simulación de Movimiento Browniano Geométrico en R

Question

Usando R, me gustaría simular una trayectoria de una ecuación de movimiento browniano geométrico usando

\begin{equation*} S(t) = S(0) \exp\left(\left(\mu - \frac{\sigma^{2}}{2}\right)t + \sigma B_{t}\right), \end{equation*}

donde $(B_t)$ es el proceso de Wiener, es decir $B_t\sim N(0,t)$ para todo $t$.

Me gustaría comparar esta trayectoria con la que obtengo utilizando el esquema de Euler-Maruyama:

\begin{equation*} S(i+1) = S(i) + mu*S(i)*delta_t + sigma*S(i)*B_{t} \end{equation*}

Me gustaría reproducir el gráfico en la página 534 del paper Higham (2001) "An algorithmic introduction to numerical simulation of SDE":

Obtuve un resultado incorrecto usando el código:

rm(list=ls())
#Simulating  Geometric Brownian motion (GMB)
tau <- 1 #tiempo de vencimiento
N <- 1000 #número de subintervalos
dt <- tau/N #longitud de cada subintervalo de tiempo
time <- seq(from=0, to=tau, by=dt) #momentos de tiempo en los que simulamos el proceso
length(time) #debería ser N+1

mu <- 0.05 #parámetro 1 de GBM
sigma <- 0.9 #parámetro 2 de GBM
X0 <- 10 #condición inicial

#simular 1 trayectoria de un movimiento browniano geométrico
Z <- rnorm(N, mean = 0, sd = 1) #muestra normal estándar de N elementos
dW <- Z*sqrt(dt) #incrementos de movimiento browniano
W <- c(0, cumsum(dW)) #movimiento browniano en cada instante de tiempo con N+1 elementos

#Solución analítica
X_analytic <- numeric(N+1) #vector de ceros, N+1 elementos
X_analytic[1] <- X0 #primer elemento de X_analytic es X0. con el bucle for encontramos los otros N elementos

for(i in 2:length(X_analytic)){
  X_analytic[i] <- X_analytic[1]*exp(mu - 0.5*sigma^2*i*dt + sigma*W[i-1])
}

#graficar X contra el tiempo
plot(time, X_analytic, type = "l", main = "Trayectoria GBM con solución analítica", 
     xlab = expression("t"[i]), ylab = expression("W"[t[i]]))

#Esquema de Euler-Maruyama
X_EM <- numeric(N+1) #vector de ceros, N+1 elementos
X_EM[1] <- X0 #primer elemento de X_EM es X0. con el bucle for encontramos los otros N elementos

for(i in 2:length(X_EM)){
  X_EM[i] <- X_EM[i-1] + mu*X_EM[i-1]*dt + sigma*dW[i-1]
}

#graficar X contra el tiempo
plot(time, X_EM, type = "l", main = "Trayectoria GBM con esquema de Euler-Maruyama", 
     xlab = expression("t"[i]), ylab = expression("W"[t[i]]))

#graficar W contra el tiempo
matplot(time, cbind(X_analytic, X_EM), type = "l", main = "GBM", 
        xlab = expression("t"[i]), ylab = expression("X"[t[i]]))

Este es el resultado:

No sé cuál de los dos es incorrecto y por qué

Answer 1

El problema es que no trazas una ruta de muestra, sino que para cada punto en el tiempo $t$, simplemente trazas una posible realización de la variable aleatoria $S_t(\omega)$. Por lo tanto, no obtienes una ruta conectada.

(Solo como un detalle, necesitarías corchetes en el exponencial en tu bucle for, es decir,

X_analytic[i] <- X_analytic[1]*exp((mu - 0.5*sigma^2)*time[i] + sigma*Z[i-1]*sqrt(time[i]))

Sin embargo, para simular una ruta de muestra de un movimiento Browniano geométrico, nota que \begin{align*} S_{t_i}=S_{t_{i-1}}\cdot\exp\left(\left(\mu-\frac{1}{2}\sigma^2\right)(t_i-t_{i-1})+\sigma B_{t_i-t_{i-1}}\right) \end{align*} En tu caso, elegiste un tamaño de paso fijo $\Delta t=t_i-t_{i-1}$ para todos los $i$ de manera que \begin{align*} S_{t_i}=S_{t_{i-1}}\cdot\exp\left(\left(\mu-\frac{1}{2}\sigma^2\right)\Delta t+\sigma \sqrt{\Delta t}Z\right), \end{align*} donde $Z\sim N(0,1)$.

Por lo tanto, cambia la línea de tu bucle for a

X_analytic[i] <- X_analytic[i-1]*exp((mu - 0.5*sigma^2)*dt + sigma*Z[i-1]*sqrt(dt))

Además, es posible que desees cambiar la línea time <- seq(from=0, to=1, by=dt) a time <- seq(from=0, to=tau, by=dt) de manera que puedas realmente hacer uso de tau. Finalmente, la línea sigma <- 0.9 es un poco ambiciosa, una volatilidad del 90% es bastante alta. Si estás modelando precios de acciones, un valor de 0.1 a 0.4 es más apropiado.

---

Edición en respuesta a tu código actualizado

No necesitas el primer bucle for para calcular la solución "analítica". Simplemente utiliza

X_analytic = X0 * exp((mu-0.5*sigma^2)*time+sigma*W)

En segundo lugar, en tu aproximación de Euler, te faltó $S$ en el último término, por lo que el paso del bucle for debería poren verdad ser

X_EM[i] <- X_EM[i-1] + mu*X_EM[i-1]*dt + sigma*X_EM[i-1]*dW[i-1]