Simulación de Vasicek/Ornstein-Uhlenbeck exponencial

Question

Estoy intentando simular los precios de las materias primas utilizando el modelo exponencial Vasicek/Ornstein-Uhlenbeck de Schwartz 1997 p. 926 Ecuación (1). Estoy utilizando la solución de forma cerrada de Vega 2018 p. 5 Ecuación (9) que es:

$$\ln(X_{t})=\ln(X_{t-1})e^{-\theta \Delta t}+\left(\mu-\frac{\sigma^2}{2\theta}\right)(1-e^{-\theta \Delta t})+\sigma\sqrt{\frac{1}{2\theta}(1-e^{-2\theta \Delta t})}\epsilon_i$$

Aquí está mi código en Python:

import numpy as np
import matplotlib.pyplot as plt
np.random.seed(123)

def gen_paths(X0, theta, mu, sigma, T, num_steps, num_sims):
    dt = float(T) / num_steps
    paths = np.zeros((num_steps + 1, num_sims), np.float64)
    paths[0] = X0
    for t in range(1, num_steps + 1):
        rand = np.random.standard_normal(num_sims)
        paths[t] = np.exp(np.log(paths[t-1]) * np.exp(-theta * dt) 
           + (mu - (sigma ** 2) / (2 * theta)) * (1 - np.exp(-theta * dt)) 
           + np.sqrt((1 - np.exp(-2 * theta * dt)) * (sigma ** 2) / (2 * theta)) * rand)
    return paths

X0 = 5
theta = 0.4
mu = 5
sigma = 0.15
T = 1
num_steps = 365
num_sims = 5

paths = gen_paths(X0, theta, mu, sigma, T, num_steps, num_sims)
plt.plot(paths[:, :10])
plt.grid(True)
plt.xlabel('time steps')
plt.ylabel('index level')
plt.show()

Y este es el resultado que estoy obteniendo:

Este no es ciertamente el resultado que esperaba, esperaba que la trayectoria fluctuara alrededor de la media a largo plazo $\mu$ y no un aumento exponencial.

Pregunta: ¿He entendido mal la exponencial Vasicek/Ornstein-Uhlenbeck que significa que este resultado es correcto y esperado o hay algo mal en mi simulación?

Actualización:

Aquí está mi nueva función después de la sugerencia en la respuesta de abajo:

def gen_paths(X0, theta, mu, sigma, T, num_steps, num_sims):
    dt = float(T) / num_steps
    paths = np.zeros((num_steps + 1, num_sims), np.float64)
    paths[0] = X0
    for t in range(1, num_steps + 1):
        rand = np.random.standard_normal(num_sims)
        z = np.log(paths[t-1])
        paths[t] = z * np.exp(-theta * dt) 
           + (mu - (sigma ** 2) / (2 * theta)) * (1 - np.exp(-theta * dt)) 
           + np.sqrt((1 - np.exp(-2 * theta * dt)) * (sigma ** 2) / (2 * theta)) * rand

    paths_new = np.exp(paths)
    return paths_new

Answer 1

Creo que tienes que ajustar tu Python un poco para:

...

def gen_paths(X0, theta, mu, sigma, T, num_steps, num_sims):
    ...
    for t in range(1, num_steps + 1):
        rand = np.random.standard_normal(num_sims)
        paths[t] = paths[t-1] * np.exp(-theta * dt) 
           + (mu - (sigma ** 2) / (2 * theta)) * (1 - np.exp(-theta * dt)) 
           + np.sqrt((1 - np.exp(-2 * theta * dt)) * (sigma ** 2) / (2 * theta)) * rand
    return paths

En concreto, hay que deshacerse del exp / log cosas. Se puede tomar la exponencial de las trayectorias después de la simulación - o se incorpora la exp / log correctamente (no se muestra en mi código anterior).