# Geometric Brownian Motion simulation in Python: strange results

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?
• Hi, usually, volatility is defined as an annualized term. So, in your example: 40% is the annualized volatility. So if you have M=365 steps, then you have dt=1/365 and you use 0.4 as vol. If you want to simulate on a daily scale, you set dt=1 and you then have to use a volatility of sqr(365)*40%. – Kermittfrog Dec 16 '20 at 11:27

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.

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

• This makes sense, however I have one issue. I computed the daily drift and volatility based on daily historical data (I edited the code in my question to show the computations). This is the time series I am working with Time Series. Does your answer in that case imply that GBM isn't suited for my specific time series or do I need another method to compute drift and volatility? I got a suggestion that my time series should be simulated using a GBM. – Willart Dec 17 '20 at 8:41
• @Willart: well in this case, then your vol and drifts might be okey. The behaviour of the source data and your simulated data look very alike. Note that the Fee is basically dependent on the volume. – Pontus Hultkrantz Dec 17 '20 at 9:53
• Well, when I create 100000 simulated paths for 1 year using drift and volatility from the last 9 years of the historical time series I get an average value which is basically 0, probably because of your explanation above that all simulations decay to zero. I don't think I'm doing something wrong, which makes me believe that my approach of simulating this time series with a GBM is wrong, probably since it doesn't fit this time series at all. – Willart Dec 17 '20 at 11:12
• @Willart: I suppose you calculate mean and vol using returns right? Yeah, well, I think the GBM isn't too bad, but check what people use for modelling trading volume. – Pontus Hultkrantz Dec 17 '20 at 18:38
• Yes exactly, daily_return = (today - yesterday)/yesterday is what I use, which should be correct I guess? I was thinking if my method to calculate drift and volatility based on returns is correct for my specific time series or do I need to use some other method to calculate them? In theory, GBM should be a good fit since the time series should also be influenced heavily by the market price which should behave very similar as a stock. Also, I got the recommendation to simulate this time series with a GBM, but with this behaviour I don't know what to make out of it. – Willart Dec 18 '20 at 10:35