1
$\begingroup$

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: Simulation exponential Vasicek/Ornstein-Uhlenbeck

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
$\endgroup$
9
  • $\begingroup$ Why do they make OU exponential instead of using traditional OU? $\endgroup$ – develarist Jan 13 at 11:27
  • $\begingroup$ @develarist Because that is what Schwartz 1997 uses and is used a lot in other literature as well. $\endgroup$ – Tharmis Jan 13 at 11:36
  • $\begingroup$ my question is why do they use it $\endgroup$ – develarist Jan 13 at 11:37
  • 2
    $\begingroup$ @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$. $\endgroup$ – Tharmis Jan 13 at 12:02
  • 1
    $\begingroup$ Welcome Quant SE, @Tharmis! $\endgroup$ – Jan Stuller Jan 14 at 8:43
1
$\begingroup$

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

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$ – Tharmis Jan 13 at 11:39
  • $\begingroup$ 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 $\endgroup$ – Kermittfrog Jan 13 at 13:08
  • $\begingroup$ 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? $\endgroup$ – Tharmis Jan 13 at 15:03
  • $\begingroup$ I think that you can literally just copy my code and calculate exp of the resulting vector afterwards. $\endgroup$ – Kermittfrog Jan 14 at 10:57
  • $\begingroup$ 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. $\endgroup$ – Tharmis Jan 14 at 11:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.