Simulación de Movimiento Browniano Geométrico en Python: resultados extraños

Question

Estoy tratando de simular el Movimiento Browniano Geométrico en Python, sin embargo los resultados que obtengo parecen muy extraños y en mi opinión no pueden ser correctos. Mi objetivo es simular cada día de 1 año. Básicamente, utilicé dos enfoques ligeramente diferentes. Basado en mi investigación, debería ser posible usar deriva y volatilidad diaria con dt = 1 o anualizar la deriva y volatilidad y usar dt = 1/365. Calculé la deriva y volatilidad diarias basadas en datos históricos diarios para el primer enfoque y las anualicé para el segundo. En consecuencia, también se cambia T. A continuación se muestra el código en Python que obtuve parcialmente de esta respuesta (Geometric Brownian Motion simulation in Python):

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(1000)

quandl.ApiConfig.api_key = "XXXXXX"

def gen_paths(S0, r, sigma, T, M, I):
    dt = float(T) / M
    paths = np.zeros((M + 1, I), np.float64)
    paths[0] = S0
    for t in range(1, M + 1):
        rand = np.random.standard_normal(I)
        paths[t] = paths[t - 1] * np.exp((r - 0.5 * sigma ** 2) * dt + sigma * np.sqrt(dt) * rand)
    return paths

# Obtener datos
df = quandl.get(...)
value = df['Value']

# Calcular rendimientos
returns = value.pct_change(1)

# Deriva y volatilidad diarias
S0 = value[-1]              # precio inicial de la acción
r = np.mean(returns)        # deriva
sigma = np.std(returns)     # volatilidad
T = 365                     # tiempo
M = 365                     # número de pasos
I = 100000                  # número de simulaciones

paths_1 = gen_paths(S0, r, sigma, T, M, I)
print("Promedio: " + str(np.average(paths_1[-1])))

# Deriva y volatilidad anuales
r = np.mean(returns) * 365
sigma = np.std(returns) * np.sqrt(365)
T = 1.0

paths_2 = gen_paths(S0, r, sigma, T, M, I)
print("Promedio: " + str(np.average(paths_2[-1])))

# Visualizar 5 caminos simulados para cada uno
fig, axs = plt.subplots(2)
axs[0].plot(paths_1[:, :5])
axs[1].plot(paths_2[:, :5])
plt.show() 

Este es el resultado que obtengo:

Para los promedios obtengo 117183 para el enfoque con deriva y volatilidad diarias y 76145 para el enfoque con deriva y volatilidad anuales. Aquí están mis preguntas:

1. ¿Estos resultados pueden ser correctos o qué hay de malo en mi código?
2. ¿No deberían ser muy similares los promedios para ambos enfoques?
3. ¿Qué enfoque es más favorable usar?

Answer 1

¡Gran pregunta!

Resumen: Tu código y matemáticas son correctos, pero utilizas un vol y deriva demasiado altos para ser realistas. Tus simulaciones se desvanecen a cero debido al alto vol y LogNormalidad.

---

Básicamente, utilicé dos enfoques ligeramente diferentes. Según mi investigación, debería ser posible utilizar deriva y volatilidad diaria con dt = 1 o anualizar deriva y volatilidad y utilizar dt = 1/365.

La solución utilizando cantidades anuales $\{t, \sigma, r\}$ es \begin{align} S(t+\Delta t) = S(t) \exp \left((r-\tfrac{1}{2}\sigma^2)\Delta t + \sigma\sqrt{\Delta t}\cdot Z \right). \end{align}

Cree nuevas cantidades diarias $\{d, s, \nu\}$ sustituyendo con $\Delta t := d/365$, $\; \sigma := s \sqrt{365}$, $\; r:= \nu \cdot 365$, obtenemos \begin{align} S(d+\Delta d) = S(d) \exp \left((\nu-\tfrac{1}{2}s^2)\Delta d + s\sqrt{\Delta d}\cdot Z \right). \end{align}

Por lo tanto, la "fórmula" se ve igual, ya sea que uses parámetros anuales o diarios.

---

1. ¿Estos resultados pueden ser correctos o qué está mal con mi código?

Tu código y resultados son correctos. El problema es que tus parámetros son enormes. Lo que asumes que son valores diarios deberían ser anuales. Lo que asumes que son anuales son por lo tanto gigantes. La tasa diaria debería ser $\nu=\tfrac{0.04}{365}$; la volatilidad debería ser $s=\tfrac{0.4}{\sqrt{365}}.$ Las versiones anuales son entonces $\sigma=0.4$ y $r=0.04$.

Lo que es confuso, aunque esperado, es que todas tus simulaciones decaen hacia cero. Esto se debe a que tu GBM es una distribución lognormal, y cuando tu varianza se acerca a infinito (tu vol era gigante), casi toda la densidad terminal se encuentra en cero. Existe una pequeña probabilidad de que el valor terminal sea muy grande, pero esto tiene una probabilidad muy pequeña, por lo que no lo vemos. Puedes ver esto trazando la PDF y CDF con una gran varianza.

---

2. ¿No deberían ser muy similares los promedios para ambos enfoques?

Estás utilizando números aleatorios diferentes para tus dos simulaciones. Haz que los números aleatorios sean iguales eliminando tu np.random.seed(1000) en la parte superior y añadiéndolo una vez antes de cada llamada a tus funciones, como así

...
np.random.seed(1000)
paths_1 = gen_paths(S0, r, sigma, T, M, I)
...
np.random.seed(1000)
paths_2 = gen_paths(S0, r, sigma, T, M, I)

Esto asegurará que utilicen la misma aleatoriedad, y verás que tus dos gráficos son idénticos.

---

¿Qué enfoque es favorable de usar?

Anual (casi) siempre, a menos que trabajes en modelado intradiario y con una resolución muy alta y para evitar problemas numéricos con underflow. Nunca he visto realmente a nadie utilizar otra cosa que no sea la anual para el modelado y la fijación de precios del GBM estándar.

---