Espero que todos estén teniendo un día bendecido,
Estoy trabajando en calibrar el modelo de Heston a datos observados, y uno de los pasos (como propuesto por Mikhalov, 2003) es comparar el desempeño de mi fórmula de precios de Heston con los resultados de una simulación de Monte Carlo.
Decidí usar la siguiente discretización para hacerlo:
$ S_{t+dt} = S_t * exp((r-\frac{1}{2} * v_t)*dt + \sqrt{v_t * dt}*Z_s) $ (Milstein)
$ v_{t+dt} = (\theta + (v_t - \theta)*e^{-\kappa * dt}) * exp(-\frac{1}{2}*\Gamma_t^{2} + \Gamma_t*Z_v)$
$ \Gamma_t = ln(1 + \frac{\sigma^2v_t(1-e^{-2\kappa dt})}{2\kappa (\theta + (v_t - \theta)*e^{-\kappa * dt})^2})$
Con las ecuaciones para $ v_{t+dt}$ y $\Gamma_t$ provenientes del esquema de coincidencia de momentos de Andersen (Proviene de F. Rouah - El modelo de Heston y sus extensiones en Matlab y C#).
Intenté implementar esto en python con el siguiente código:
# Parámetros utilizados solo para la implementación, generaré nuevos parámetros al verificar la precisión de mi precio
S0 =100
K = 100
v0 = 0.1
r = 0.02
kappa = 1.5768
theta = 0.0398
sigma = 0.3
rho = -0.5711
tau = 1
N = 252
M = 1000
def heston_model_sim(S0, v0, rho, kappa, theta, sigma, t, N, M):
# Inicializando otros parámetros
dt = t/N
mu = np.array([0,0])
cov = np.array([[1,rho],[rho,1]])
g0 = np.log(1 + (sigma ** 2 * v0 * (1 - np.exp(-2 * kappa * dt)))/(2 * kappa * (theta + (v0 - theta) * np.exp(- kappa * dt)) ** 2))
# Arrays para almacenar precios y varianzas
S = np.full(shape = (N+1, M), fill_value = S0)
v = np.full(shape = (N+1, M), fill_value = v0)
gamma = np.full(shape = (N+1, M), fill_value = g0)
# Muestreo de movimientos brownianos correlacionados bajo medida de riesgo neutral
Z = np.random.multivariate_normal(mu, cov, (N, M))
for i in range(1, N+1):
S[i] = S[i-1] * np.exp((r - 0.5 * v[i-1]) * dt + np.sqrt(v[i-1] * dt) * Z[i-1,:,1])
gamma[i-1] = np.log(1 + (sigma ** 2 * v[i-1] * (1 - np.exp(-2 * kappa * dt)))/(2 * kappa * (theta + (v[i-1] - theta) * np.exp(- kappa * dt)) ** 2))
v[i] = (theta + (v[i-1] - theta) * np.exp(- kappa * dt)) * np.exp(-0.5 * gamma[i-1]**2 + gamma[i-1] * Z[i-1,:,0])
Return S,v(Observación - No soy el mejor para pegar código en stackexchange pero la indexación de la declaración de retorno es correcta en mi código.)
Y este es el resultado que obtengo:
Parece que hay un (gran) problema con mi proceso de volatilidad (quizás también con el proceso de acción pero es menos llamativo) pero no puedo encontrar de dónde proviene y cómo puedo solucionarlo.
¿Ves algún error en mi código o en la discretización que estoy utilizando?
Muchas gracias de antemano
Actualización:
esto es lo que sucede cuando bajo la volatilidad inicial:
Pareciera ser más correcto, ¿pero por qué hay un problema al comenzar con una volatilidad más alta?
