Simulating Geometric Brownian Motion (GBM) in Python
Stochastic processes are essential tools in quantitative finance for modeling the random evolution of market variables like stock prices. Geometric Brownian Motion (GBM) is arguably the most fundamental and widely used model, forming the basis of the Black-Scholes-Merton option pricing theory. This article demonstrates how to simulate GBM using Python, covering both terminal values and full price paths, and discusses its applications in finance and trading. We’ll use Monte Carlo simulation, a technique relying on repeated random sampling to obtain numerical results.
Setup and Preliminaries
Before simulating, we need to import necessary libraries and set up our environment.
Python
# -*- coding: utf-8 -*-
import math
import numpy as np
import numpy.random as npr
from pylab import plt, mpl
import scipy.stats as scs # Used for statistics function below
# Setup plotting style and font
# plt.style.use('seaborn-v0_8') # Use a specific available style if needed
mpl.rcParams['font.family'] = 'serif'
# %config InlineBackend.figure_format = 'svg' # Jupyter/Colab specific line
# Set seed for reproducibility and numpy print precision
npr.seed(100)
np.set_printoptions(precision=4)
# Helper function for printing statistics (useful for comparing methods)
def print_statistics(a1, a2):
''' Prints selected statistics for two datasets. '''
sta1 = scs.describe(a1)
sta2 = scs.describe(a2)
print('%14s %14s %14s' % ('statistic', 'data set 1', 'data set 2'))
print(45 * "-")
print('%14s %14.3f %14.3f' % ('size', sta1[0], sta2[0]))
print('%14s %14.3f %14.3f' % ('min', sta1[1][0], sta2[1][0]))
print('%14s %14.3f %14.3f' % ('max', sta1[1][1], sta2[1][1]))
print('%14s %14.3f %14.3f' % ('mean', sta1[2], sta2[2]))
print('%14s %14.3f %14.3f' % ('std', np.sqrt(sta1[3]), np.sqrt(sta2[3])))
print('%14s %14.3f %14.3f' % ('skew', sta1[4], sta2[4]))
print('%14s %14.3f %14.3f' % ('kurtosis', sta1[5], sta2[5]))
print("=== GBM Simulation Section ===")- We import
numpyfor numerical arrays,numpy.randomfor random numbers,matplotlib(pylab) for plotting, andscipy.statsfor descriptive statistics. - Setting a
seedensures our random simulations are repeatable.
The GBM Model
GBM models the price St of an asset over time. It assumes the percentage change in price follows a random walk with drift. Its Stochastic Differential Equation (SDE) is:
\[d S_t=r S_t d t+\sigma S_t d W_t\]
where \(r\) is the constant risk-free interest rate (drift), \(σ\) is the constant volatility, and \(dW_t\) represents the randomness (a Wiener process).
1. Simulating Terminal Asset Value
For some applications, like pricing simple European options, we only need the distribution of the asset price at a specific future time \(T\). The analytical solution for \(S_T\) is:
\[S_T=S_0 \exp \left(\left(r-\frac{1}{2} \sigma^2\right) T+\sigma \sqrt{T} Z\right)\]
where \(S_0\) is the initial price and \(Z∼N(0,1)\) is a standard normal random variable.
Python
# === Simulation: Random Variables (Terminal Value) ===
print("\n--- Simulating Terminal Asset Value (GBM) ---")
S0 = 100 # Initial asset price
r = 0.05 # Risk-free rate
sigma = 0.25 # Volatility
T = 2.0 # Time horizon (years)
I = 10000 # Number of simulations
# Method 1: Analytical formula for terminal value
ST1 = S0 * np.exp((r - 0.5 * sigma ** 2) * T +
sigma * math.sqrt(T) * npr.standard_normal(I))
# Plot histogram for Method 1
plt.figure(figsize=(10, 6))
plt.hist(ST1, bins=50, label='Terminal Value (Analytical Formula)')
plt.xlabel('index level')
plt.ylabel('frequency')
plt.title('Simulated Terminal Asset Value (GBM) - Method 1')
plt.legend()
plt.show()
# Method 2: Using lognormal distribution directly
ST2 = S0 * npr.lognormal((r - 0.5 * sigma ** 2) * T,
sigma * math.sqrt(T), size=I)
# Plot histogram for Method 2
plt.figure(figsize=(10, 6))
plt.hist(ST2, bins=50, label='Terminal Value (Lognormal Function)')
plt.xlabel('index level')
plt.ylabel('frequency')
plt.title('Simulated Terminal Asset Value (GBM) - Method 2')
plt.legend()
plt.show()
# Compare statistics
print("\nComparing statistics of ST1 and ST2:")
print_statistics(ST1, ST2)— Simulating Terminal Asset Value (GBM) —
Comparing statistics of ST1 and ST2: statistic data set 1 data set 2 ——————————————— size 10000.000 10000.000 min 28.230 27.105 max 414.825 409.110 mean 110.468 110.537 std 40.271 40.317 skew 1.167 1.188 kurtosis 2.533 2.697
- Explanation: We simulate \(I=10000\) possible terminal prices \(S_T\). Method 1 uses the formula directly
with standard normal random numbers (
npr.standard_normal). Method 2 leverages the fact that \(S_T\) is log-normally distributed and usesnpr.lognormal, providing the correct mean and standard deviation for the underlying normal distribution (\(\ln(ST/S0)\)). - The statistics confirm both methods yield virtually identical results, validating our understanding.
- Use Cases: Pricing European options (averaging payoff \(\max(S_T−K,0)\) or \(\max(K−S_T,0)\) over simulations), estimating terminal Value at Risk (VaR).
2. Simulating GBM Paths
For path-dependent options or dynamic analysis, we need the entire price trajectory. We discretize time \(T\) into \(M\) steps of size \(Δt=T/M\) and simulate iteratively:
\[S_{t+\Delta t}=S_t \exp \left(\left(r-\frac{1}{2} \sigma^2\right) \Delta t+\sigma \sqrt{\Delta t} Z_t\right)\]
Python
# --- Geometric Brownian Motion (GBM) Path ---
print("\n--- Geometric Brownian Motion (GBM) Path Simulation ---")
M = 50 # Number of time steps
dt = T / M # Time interval
S = np.zeros((M + 1, I)) # Array to hold paths S[time_step, path_index]
S[0] = S0 # Set initial value for all paths
for t in range(1, M + 1):
S[t] = S[t - 1] * np.exp((r - 0.5 * sigma ** 2) * dt +
sigma * math.sqrt(dt) * npr.standard_normal(I))
# Plot histogram of final values from path simulation
plt.figure(figsize=(10, 6))
plt.hist(S[-1], bins=50, label='Final Value from Path Sim')
plt.xlabel('index level')
plt.ylabel('frequency')
plt.title('Final Asset Values from GBM Path Simulation')
plt.legend()
plt.show()
# Plot first 10 simulated paths
plt.figure(figsize=(10, 6))
plt.plot(np.linspace(0, T, M+1), S[:, :10], lw=1.5) # Add time axis
plt.xlabel('time (years)')
plt.ylabel('index level')
plt.title('Sample GBM Paths')
plt.show()— Geometric Brownian Motion (GBM) Path Simulation —
- Explanation: A 2D NumPy array S stores
the paths. The loop progresses through time, updating all I
paths simultaneously at each step t using the value from
the previous step t-1 and a new set of random numbers
z. Vectorization makes this efficient. - Comparing the
final simulated values S[-1] with the terminal values
ST2 confirms the path simulation is consistent with the
direct method. - Use Cases: Pricing path-dependent
options (barriers, lookbacks, Asians), simulating dynamic hedging
strategies, multi-period VaR calculations, backtesting trading
strategies based on price evolution.
Conclusion
Geometric Brownian Motion is a foundational model in finance. Simulating GBM using Python and NumPy allows us to price derivatives, assess risk, and analyze strategies. While simple, it provides the basis for understanding more complex models. This article showed how to simulate both terminal values and full price paths efficiently using Monte Carlo methods.