Simulating exponential Vasicek/Ornstein-Uhlenbeck

Question

I am trying to simulate commodity prices using the exponential Vasicek/Ornstein-Uhlenbeck model from Schwartz 1997 p. 926 Equation (1). I am using the closed form solution from Vega 2018 p. 5 Equation (9) which is:

$$\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$$

Here is my code in 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()

And here is the result I am getting:

This is certainly not the result I expected, I expected the path to fluctuate around the long term mean $\mu$ and not a exponential rise.

Question: Did I misunderstood the exponential Vasicek/Ornstein-Uhlenbeck meaning that this result is correct and expected or is there something wrong in my simulation?

Update:

Here is my new function after the suggestion in the answer below:

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

I think you have to adjust your Python a bit to:

...

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

Specifically, you need to get rid of the exp / log stuff. You can take the exponential of the paths after simulation - or you incorporate the exp / log properly (not shown in my code above).