Skip to content
All library documents

Simulating an Exponential Ornstein–Uhlenbeck Commodity Price Model

Article Quant Q&A · Author: Tharmis
Read the original

Summary

The document investigates a simulation of commodity prices using an exponential Vasicek or Ornstein–Uhlenbeck model. It presents a closed-form time-step expression for the logarithm of the price, along with a Python implementation and an unexpected upward path pattern. The response points to a mismatch in how the recurrence is applied: the stated expression evolves log price, so the state variable and the final transformation back to price must be handled consistently. The example highlights that applying logarithms and exponentials at the wrong stage changes the simulated process.

The discussion connects mean reversion in the modeled state to the behavior of simulated prices, and includes an updated attempt that evolves a log-state before exponentiating the paths. However, the answer's replacement recurrence appears to mix the log-price equation with a level-price update, so readers should verify the model definition and parameterization before relying on that code. The post provides no statistical validation or comparison against theoretical moments. Its practical lesson is about translating a stochastic differential equation into a discrete simulation without confusing the log state with the price level.

Key ideas

  • The displayed transition equation specifies the evolution of the logarithm of the price.
  • A simulation must keep its state variable consistent with the model equation and exponentiate only when converting log prices to prices.
  • Mixing a log-price transition with a level-price recurrence can produce incorrect paths.
  • The example motivates checking the model specification and theoretical moments before interpreting simulated behavior.

Tags

CommoditiesStatisticsMean reversionBacktesting
Full text
# Simulating exponential Vasicek/Ornstein-Uhlenbeck


# 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[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 by Kermittfrog (score 1)

https://quant.stackexchange.com/a/60483

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

Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)

This summary was written by Stratmill's research agent from the original; it is not a copy of the source.