# Simulating exponential Vasicek/Ornstein-Uhlenbeck

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 = 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 = 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

• Why do they make OU exponential instead of using traditional OU? – develarist Jan 13 at 11:27
• @develarist Because that is what Schwartz 1997 uses and is used a lot in other literature as well. – Tharmis Jan 13 at 11:36
• my question is why do they use it – develarist Jan 13 at 11:37
• @develarist I get that, but it does seem like the exponential does also fluctuate roughly around $\mu$, not the same as the normal OU though. Maybe I formulated it poorly, I never expected the same results as the normal O/U, I was just trying to say that based on my research the result I am getting can't be correct and my expectation was for the results to still fluctuate somewhere around $\mu$. – Tharmis Jan 13 at 12:02
• Welcome Quant SE, @Tharmis! – Jan Stuller Jan 14 at 8:43

## 1 Answer

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).

• The results from your code look a lot more likely, but am I still simulating the same stochastic process? Care to explain what exactly is wrong with the exp / log stuff or how I can fix that? I used exactly the formula from the mentioned paper, so I am not sure why in Python it should be different. – Tharmis Jan 13 at 11:39
• I think you had one closing bracket at the wrong place. If you look carefully, you can replace $ln(X_t)$ with some variable $z_t$ everywhere in your code and simply do the $exp$ing afterwards. HTH – Kermittfrog Jan 13 at 13:08
• I did exactly that (see the Update in my question) however that also doesn't work. If I don't immediatly do the exping the values in the log get negative which yields to an error. Did I misunderstand you somehow? – Tharmis Jan 13 at 15:03
• I think that you can literally just copy my code and calculate exp of the resulting vector afterwards. – Kermittfrog Jan 14 at 10:57
• That gives me values around 140, but since my $X_0$ and $\mu$ is 5 that is not the correct result. Either I am completely missing something or my results are correct and I just misunderstood the exponential O/U completely. – Tharmis Jan 14 at 11:36