Geometric Brownian Motion simulation, terminal distributions, and Black-Scholes comparison
I post the models as I build them, the changes I make to them, and the results that did not work. Shorter than these write ups, and there are a lot more of them.
Follow me on X
Want the numbers rather than the code? The Monte Carlo calculator simulates the paths in the browser and compares the result against the closed form.
Under the risk-neutral measure, the stock price follows a Geometric Brownian Motion (GBM):
$$dS = (r - q)\, S \, dt + \sigma \, S \, dW_t$$
Applying Itô's Lemma to $\ln S$ and integrating over $[0, T]$ gives the exact solution:
$$S_T = S_0 \; \exp\!\left[\left(r - q - \tfrac{\sigma^2}{2}\right)T + \sigma\sqrt{T}\; Z\right]$$
where $Z \sim \mathcal{N}(0, 1)$.
The continuously compounded return is:
$$\ln\!\left(\frac{S_T}{S_0}\right) = \underbrace{\left(r - q - \tfrac{\sigma^2}{2}\right)T}_{\text{drift}} + \underbrace{\sigma \sqrt{T}\; Z}_{\text{diffusion}}$$
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt# ── Parameters ──
S0 = 50 # Initial stock price
K = 50 # Strike price
T = 9 / 12 # Time to maturity (years)
r = 0.04 # Risk-free interest rate
q = 0.00 # Continuous dividend yield
sigma = 0.30 # Volatility (base)
N = 10_000_000 # Number of simulationsnp.random.seed(42)
z = np.random.normal(0, 1, N)
st = S0 * np.exp((r - q - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * z)# Terminal prices at two different implied volatilities
sigma_low = 0.20
sigma_high = 0.30
S_T_low = S0 * np.exp((r - q - 0.5 * sigma_low**2) * T + sigma_low * np.sqrt(T) * z)
S_T_high = S0 * np.exp((r - q - 0.5 * sigma_high**2) * T + sigma_high * np.sqrt(T) * z)
# Both share the same risk-neutral mean
mean_both = np.mean(S_T_low)
median_low = np.median(S_T_low)
median_high = np.median(S_T_high)
fig, ax = plt.subplots(figsize=(10, 5))
ax.hist(S_T_low, bins=100, density=True, alpha=0.55, color="steelblue", label=f"sigma = {sigma_low:.0%}")
ax.hist(S_T_high, bins=100, density=True, alpha=0.45, color="indianred", label=f"sigma = {sigma_high:.0%}")
ax.axvline(mean_both, color="black", linestyle="-", linewidth=1.3, label=f"Mean = {mean_both:.2f}")
ax.axvline(median_low, color="steelblue", linestyle="--", linewidth=1.2, label=f"Median (20%) = {median_low:.2f}")
ax.axvline(median_high, color="indianred", linestyle="--", linewidth=1.2, label=f"Median (30%) = {median_high:.2f}")
ax.set_title("Terminal Price Distribution -- sigma = 20% vs 30%", fontweight="bold")
ax.set_xlabel("$S_T$")
ax.set_ylabel("Density")
ax.legend()
plt.tight_layout()
plt.show()For the closed-form alternative we benchmark against here, see the Black-Scholes notebook directly.
# Monte Carlo simulation (single-step terminal value)
# Call
payoff_call = np.maximum(st - K, 0)
premium_call = np.mean(payoff_call) * np.exp(-r * T)
se_call = np.std(payoff_call, ddof=1) / np.sqrt(N)
# Put
payoff_put = np.maximum(K - st, 0)
premium_put = np.mean(payoff_put) * np.exp(-r * T)
se_put = np.std(payoff_put, ddof=1) / np.sqrt(N)# Black-Scholes analytical prices (for comparison)
d1 = (np.log(S0 / K) + (r - q + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
bs_call = S0 * np.exp(-q * T) * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
bs_put = K * np.exp(-r * T) * norm.cdf(-d2) - S0 * np.exp(-q * T) * norm.cdf(-d1)
print(premium_call, bs_call)
print(premium_call - bs_call)## Paths matrix, now we go with two dimensions -> (n_paths, n_steps)
n_paths = 175
n_steps = int(T * 252)
dt = T / n_steps
np.random.seed(7)
Z = np.random.normal(0, 1, (n_paths, n_steps))
# Build price matrix: each row is one path
paths = np.zeros((n_paths, n_steps + 1))
paths[:, 0] = S0# Paths composition and terminal value
for t in range(n_steps):
paths[:, t + 1] = paths[:, t] * np.exp(
(r - q - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * Z[:, t]
)
S_terminal = paths[:, -1]# Classification of terminal values, call option
atm_band = 0.5
is_itm = S_terminal > (K + atm_band)
is_otm = S_terminal < (K - atm_band)
is_atm = ~is_itm & ~is_otmdays = np.arange(n_steps + 1)
plt.figure(figsize=(10, 5))
for i in np.where(is_otm)[0]:
plt.plot(days, paths[i], color="red", alpha=0.15, linewidth=0.5)
for i in np.where(is_atm)[0]:
plt.plot(days, paths[i], color="gold", alpha=0.6, linewidth=0.8)
for i in np.where(is_itm)[0]:
plt.plot(days, paths[i], color="green", alpha=0.15, linewidth=0.5)
plt.axhline(K, color="black", linestyle="--", linewidth=1.2, label=f"Strike K = {K}")
plt.ylabel("Stock Price ($)")
plt.xlabel("Days")
plt.title(f"MC simulation -- {n_paths} paths (sigma = {sigma:.0%})")
plt.legend([
plt.Line2D([], [], color="green", linewidth=2),
plt.Line2D([], [], color="gold", linewidth=2),
plt.Line2D([], [], color="red", linewidth=2),
plt.Line2D([], [], color="black", linewidth=1.2, linestyle="--"),
], [
f"ITM ({np.sum(is_itm)})",
f"ATM ({np.sum(is_atm)})",
f"OTM ({np.sum(is_otm)})",
f"Strike = {K}",
], loc="upper left")
plt.tight_layout()
plt.show()