Simulating Geometric Brownian Motion (GBM)

We first assume the asset price follows a Stochastic Differential Equation (SDE) and it's generalised as GBM.

$$dS_t=\mu S_t \ dt +\sigma S_t \ dW_t$$

$$ ln(S_t)\sim N\Bigg( ln(S_0) + \big( \mu-\frac{1}{2}\sigma^2 \big) \Bigg)t,\sigma^2 t $$

and

$$ S_t=S_0\ e^{(\mu-\frac{1}{2}\sigma^2)t+\sigma W_t} $$

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Define the Parameters here
# number of simulation
M = 1000
# number of steps
n = 100
# time in year
T=1

# drift
mu = 0.1
# volatility
sigma = 0.2
# initial price
S0=100

For the Stochastic Term (the 2nd term) of our GBM, $$W_t \sim N(0, \sqrt{t})$$ Because we assume the stochastic term follows a Brownian Motion. Also, we showed that the Brownian Motion converges to a Normal Distribution.

The Variance of that N( , ) is the qudratic variation of our Brownian Process. $Var(W)=t$

$$ S_t=S_0\ e^{(\mu-\frac{1}{2}\sigma^2)t+\sigma W_t} $$

We let, $S_t=S_0 \cdot af$, and $af=e^{(\mu-\frac{1}{2}\sigma^2)t+\sigma W_t}$

dt = T/n

# simulate af, af is the random rate each period
af = np.exp(
    (mu - 1/2*np.square(sigma))*dt
    + sigma * np.random.normal(0, np.sqrt(dt), size = (n,M))
) # size = (n,M) means n row, and M collumn. Each Vector represents a simulation.

#We need the Asset Price begins from S0, so stack S_t with a row of 100.
af= np.vstack([np.ones(M), af])

# cumproduct af to get a "compounding" return, and multiply it by the Base S0
St=S0*af.cumprod(axis=0)

St.shape

(101, 1000)

Set the interval for x-axis

# Define a time interval, there are n+1 periods, coz a initial S0 were added
time = np.linspace(0,T,n+1) 

# Create an np.array in the same shape as St
# fill in value of variable "time" in each row, 
# so the shape must initally be create in (M row, n+1 collumn), and tranpose after filling values
tt = np.full(shape=(M,n+1),fill_value=time).T 

plt.plot(tt, St)
plt.xlabel('Year $(t)$')
plt.ylabel('Stock Price $(S_t)$')
plt.title(
    "Realizations of Geometric Brownian Motion\n $dS_t = \mu S_t dt + \sigma S_t dW_t$\n $S_0 = {0}, \mu = {1}, \sigma = {2}$".format(S0, mu, sigma)
)
plt.show()

P.S. We could see the distribution of ending value of S, St.

It's not normally distributed, and a right skewness is here.

Why is that?

Because it is surely not a normal distribution. $ln(S_t)$ follows a normal, so $S_t$ follows a log-normal distribution!

plt.hist(St[-1])
plt.title('$(S_t \sim log Normal)$')
plt.show()

plt.hist(np.log(St[-1]))
plt.title('$(ln(S_t) \sim log Normal)$')
plt.show()

by Fanyu Zhao