Geometric Brownian Motion simulation in Python: strange results

Question

I am trying to simulate Geometric Brownian Motion in Python, however the results that I get seem very strange and in my opinion they can't be correct. My goal is to simulate each day of 1 year. Basically, I used two slightly different approaches. Based on my research, it should be either possible to use daily drift and volatility with dt = 1 or annualize drift and volatility and use dt = 1/365. I computed the daily drift and volatility based on daily historical data for the first approach and annualized them for the second. Accordingly also T is changed. Below is the code in Python which I partially got from this answer (Geometric Brownian Motion simulation in Python):

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(1000)

quandl.ApiConfig.api_key = "XXXXXX"


def gen_paths(S0, r, sigma, T, M, I):
    dt = float(T) / M
    paths = np.zeros((M + 1, I), np.float64)
    paths[0] = S0
    for t in range(1, M + 1):
        rand = np.random.standard_normal(I)
        paths[t] = paths[t - 1] * np.exp((r - 0.5 * sigma ** 2) * dt + sigma * np.sqrt(dt) * rand)
    return paths


# Get data
df = quandl.get(...)
value = df['Value']

# Compute returns
returns = value.pct_change(1)


# Daily drift and volatility
S0 = value[-1]              # initial stock price
r = np.mean(returns)        # drift
sigma = np.std(returns)     # volatility
T = 365                     # time
M = 365                     # number of steps
I = 100000                  # number of simulations

paths_1 = gen_paths(S0, r, sigma, T, M, I)
print("Average: " + str(np.average(paths_1[-1])))


# Yearly drift and volatility
r = np.mean(returns) * 365
sigma = np.std(returns) * np.sqrt(365)
T = 1.0

paths_2 = gen_paths(S0, r, sigma, T, M, I)
print("Average: " + str(np.average(paths_2[-1])))

# Visualize 5 simulated paths for each
fig, axs = plt.subplots(2)
axs[0].plot(paths_1[:, :5])
axs[1].plot(paths_2[:, :5])
plt.show() 

Here is the result that I get:

For the averages I get 117183 for the approach with daily drift and volatility and 76145 for the approach with yearly drift and volatility. Here are my questions:

1. Can these results be correct or what is wrong with my code?
2. Shouldnt the averages be very similar for both approaches?
3. Which approach is favorable to use?

Answer 1

Great question!

Abstract: Your code and math are correct, but you use too high vol and drift to be real world realistic. Your simulations decay to zero due to high vol and LogNormality.

---

Basically, I used two slightly different approaches. Based on my research, it should be either possible to use daily drift and volatility with dt = 1 or annualize drift and volatility and use dt = 1/365

The solution using yearly quantities $\{t, \sigma, r\}$ is \begin{align} S(t+\Delta t) = S(t) \exp \left((r-\tfrac{1}{2}\sigma^2)\Delta t + \sigma\sqrt{\Delta t}\cdot Z \right). \end{align}

Create new daily quantities $\{d, s, \nu\}$ by subsituting with $\Delta t := d/365$, $\; \sigma := s \sqrt{365}$, $\; r:= \nu \cdot 365$, we get \begin{align} S(d+\Delta d) = S(d) \exp \left((\nu-\tfrac{1}{2}s^2)\Delta d + s\sqrt{\Delta d}\cdot Z \right). \end{align}

Hence, the "formula" looks the same, whether you use annual or daily parameters.

---

1. Can these results be correct or what is wrong with my code?

Your code and results are correct. The thing is that your parameters are enormous. What you assume to be daily values should be annual. What you assume to be annual are therefore gigantic. Daily rate should be $\nu=\tfrac{0.04}{365}$; vol should be $s=\tfrac{0.4}{\sqrt{365}}.$ Annual versions are then $\sigma=0.4$ and $r=0.04$.

What is confusing, although expected, is that all your simulations decay towards zero. This is because your GBM is a lognormal distribution, and when your variance approaches infinity (your vol was gigantic), almost all of the terminal density lies at zero. There is a small probability that the terminal value is very large, but this has very small probabiltiy, so we don't see this. You can see this by plotting the PDF and CDF with a large variance.

---

2. Shouldnt the averages be very similar for both approaches?

You use different random numbers for your two simulations. Enforce the random numbers to be the same by remove your np.random.seed(1000) at the top and add it once before each call to your functions, like this

...
np.random.seed(1000)
paths_1 = gen_paths(S0, r, sigma, T, M, I)
...
np.random.seed(1000)
paths_2 = gen_paths(S0, r, sigma, T, M, I)

This will ensure that they use the same randomness, and you will see that your two plots are identical.

---

Which approach is favorable to use?

Annual (almost) always, unless you work on intraday modelling and really high resolution and to avoid numerical issues with underflow. I have never really seen anyone use anything other than annual for standard GBM modelling and pricing.

---