Geometric Brownian Motion

A few simple facts to get up to speed:

• GBM assumes that percentage returns — not price changes — are random, normally distributed, and independent over time.
• Prices therefore stay positive and their terminal distribution is log-normal.
• 2 parameters fully describe the model — both estimable directly from historical data:
  – the drift $\mathbf{\mu }$ (expected annual log-growth plus half-variance)
  – the volatility $\mathbf{\sigma }$ (annualised standard deviation of log returns)

The logic before diving into equations:

• Returns compound, prices don't add
  – A $10 move means something very different at SPY = 100 than at SPY = 600.
  – What is comparable across price levels is the relative move.
  – GBM makes randomness proportional to the current price: $\mathbf{dS}=\mathbf{\mu S}\mathbf{dt}+\mathbf{\sigma S}\mathbf{dW}$.
• Why log-normal
  – If each day multiplies the price by a small random gross return, the log-price is a sum of many small independent shocks.
  – Sums of independent shocks are (approximately) normal.
  – Normal log-price $\Rightarrow$ log-normal price: skewed right, floored at zero.
• Drift vs. noise
  – Over short horizons the noise term dominates: $\mathbf{\sigma }\sqrt{\mathbf{t}}$ shrinks slower than $\mathbf{\mu t}$ as $\mathbf{t}\to 0$.
  – Over long horizons drift wins.
  – This is the same $\sqrt{\mathbf{t}}$ rule that governs plain Brownian motion — GBM simply wraps it in an exponential.

Below the theoretical foundation and respective equations:

• The Stochastic Differential Equation (SDE) and its solution
  – A differential equation with a random term — so it yields a distribution of paths, not only one.
  – $\mathbf{\mu }{\mathbf{S}}_{\mathbf{t}}\mathbf{dt}$ is the drift (deterministic trend); $\mathbf{\sigma }{\mathbf{S}}_{\mathbf{t}}\mathbf{d}{\mathbf{W}}_{\mathbf{t}}$ is the diffusion (the noise).

  $$\mathbf{d}{\mathbf{S}}_{\mathbf{t}}=\mathbf{\mu }{\mathbf{S}}_{\mathbf{t}}\mathbf{dt}+\mathbf{\sigma }{\mathbf{S}}_{\mathbf{t}}\mathbf{d}{\mathbf{W}}_{\mathbf{t}}$$
• Applying Itô's lemma to $\ln {\mathbf{S}}_{\mathbf{t}}$
  – The $-\frac{1}{2}{\mathbf{\sigma }}^2$ correction comes from $(\mathbf{dW}{)}^2=\mathbf{dt}$.

  $${\mathbf{S}}_{\mathbf{t}}={\mathbf{S}}_0\exp \left[\right(\mathbf{\mu }-\frac{1}{2}{\mathbf{\sigma }}^2)\mathbf{t}+\mathbf{\sigma }{\mathbf{W}}_{\mathbf{t}}]$$
• Exact discretisation
  – Because the solution is closed-form, we can simulate without discretisation error at any step size $\mathbf{\Delta }\mathbf{t}$.

  $${\mathbf{S}}_{\mathbf{t}+\mathbf{\Delta }\mathbf{t}}={\mathbf{S}}_{\mathbf{t}}\cdot \exp \left[\right(\mathbf{\mu }-\frac{1}{2}{\mathbf{\sigma }}^2)\mathbf{\Delta }\mathbf{t}+\mathbf{\sigma }\sqrt{\mathbf{\Delta }\mathbf{t}}\mathbf{Z}],\mathbf{Z}\sim \mathcal{N}(0,1)$$

The key implementations of those small functions: one draws per-step gross returns, the other compounds them into price paths.

# Implementation

def gbm_returns(mu, sigma, dt, n_steps, n_paths):
    """One-step gross returns under GBM (shape: n_steps x n_paths)."""
    z = np.random.normal(size=(n_steps, n_paths))
    return np.exp((mu - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * z)

def gbm_paths(s0, mu, sigma, dt, n_steps, n_paths):
    """Price paths of shape (n_steps + 1, n_paths), starting at s0."""
    rets = gbm_returns(mu, sigma, dt, n_steps, n_paths)
    return s0 * np.vstack([np.ones(rets.shape[1]), rets]).cumprod(axis=0)

• Calibrate $\mathbf{\mu }$ and $\mathbf{\sigma }$ from real data
  – We pull daily SPY prices, estimate annualised drift and volatility from log returns, and let the data drive the simulation.

TICKER = "SPY"
START, END = "2018-01-01", "2024-12-31"

import yfinance as yf

px = (yf.download(TICKER, start=START, end=END, auto_adjust=True, progress=False)["Close"]
        .squeeze().rename(TICKER).dropna())

log_ret = np.log(px / px.shift(1)).dropna()

s0        = float(px.iloc[-1])
sigma_hat = float(log_ret.std() * np.sqrt(252))
mu_hat    = float(log_ret.mean() * 252 + 0.5 * sigma_hat**2)  # drift of the SDE, not of log returns

print(f"{TICKER}: {len(px)} daily observations, last close = {s0:,.2f}")
print(f"mu_hat    = {mu_hat:.2%}  (annualised drift)")
print(f"sigma_hat = {sigma_hat:.2%}  (annualised volatility)")

SPY: 1760 daily observations, last close = 578.32
mu_hat    = 14.75%  (annualised drift)
sigma_hat = 19.54%  (annualised volatility)

• Simulate 1,000 5-year paths
  – Grey lines are individual paths, the dashed line is the mean of the simulated distribution, and the band spans the 5th–95th percentile.

DT      = 1 / 252
HORIZON = 252 * 5          # 5 years
N_PATHS = 1_000

paths = gbm_paths(s0, mu_hat, sigma_hat, DT, HORIZON, N_PATHS)
t_ax  = np.arange(paths.shape[0]) / 252

p5, p95 = np.percentile(paths, [5, 95], axis=1)

fig, ax = plt.subplots(figsize=(10, 5))
ax.fill_between(t_ax, p5, p95, color=AQUA, alpha=0.35, lw=0, zorder=1,
                label="5\u201395% band")
ax.plot(t_ax, paths[:, :200], color="#BBBBBB", lw=0.35, alpha=0.45, zorder=2)
ax.plot(t_ax, paths.mean(axis=1), color=INK, lw=2.2, ls="--", zorder=3,
        label="Mean path")
ax.set_title(f"{TICKER}: {N_PATHS:,} GBM paths, 5 years  "
             f"(mu={mu_hat:.1%}, sigma={sigma_hat:.1%})", **HEAD)
ax.set_xlabel("Years"); ax.set_ylabel("Simulated price"); ax.legend()
plt.tight_layout(); plt.show()

• Terminal distribution — is it log-normal?
  – Verifies the code: the simulated histogram should match the closed-form density $\ln ({\mathbf{S}}_{\mathbf{T}}/{\mathbf{S}}_0)\sim \mathcal{N}((\mathbf{\mu }-\frac{1}{2}{\mathbf{\sigma }}^2)\mathbf{T},{\mathbf{\sigma }}^2\mathbf{T})$ — if it doesn't, the ½σ² term is usually the culprit.
  – Proves the claim: section 2 asserts prices end up log-normal; this is where you see the right skew rather than take it on faith.
  – The practical lesson: the mean sits above the median, so the “average” projected outcome is not the typical one — most paths land below it.

from scipy import stats

T   = HORIZON / 252
s_T = paths[-1]

x   = np.linspace(s_T.min(), s_T.max(), 400)
pdf = stats.lognorm.pdf(x, s=sigma_hat*np.sqrt(T),
                        scale=s0*np.exp((mu_hat - 0.5*sigma_hat**2)*T))

fig, ax = plt.subplots(figsize=(10, 4))
ax.hist(s_T, bins=60, density=True, color=AQUA, alpha=0.55,
        edgecolor=INK, linewidth=0.25, label="Simulated $S_T$")
ax.plot(x, pdf, color=INK, lw=2.2, label="Theoretical log-normal")
ax.axvline(s_T.mean(),     color=INK, lw=1.8, ls="--", label=f"Mean   {s_T.mean():,.0f}")
ax.axvline(np.median(s_T), color=INK, lw=1.8,          label=f"Median {np.median(s_T):,.0f}")
ax.set_title(f"{TICKER}: terminal price distribution after {T:.0f} years", **HEAD)
ax.set_xlabel("Price"); ax.legend()
plt.tight_layout(); plt.show()

• Sanity check: switch the drift off
  – With $\mathbf{\mu }=0$ the process becomes a martingale — no expected gain — so the mean terminal price must land back at ${\mathbf{S}}_0$.
  – Zeroing one parameter isolates the rest: this catches a missing $-\frac{1}{2}{\mathbf{\sigma }}^2$ correction that a full simulation might hide.
  – How to read the result: the sample mean wobbles around ${\mathbf{S}}_0$ with standard error $\approx {\mathbf{S}}_0\sqrt{{\mathbf{e}}^{{\mathbf{\sigma }}^2\mathbf{T}}-1}/\sqrt{\mathbf{N}}$. Within ~2 SE is noise; a persistent gap in one direction that doesn't shrink as $\mathbf{N}$ grows is a bug.

paths_nd = gbm_paths(s0, 0.0, sigma_hat, DT, HORIZON, N_PATHS)
print(f"Mean terminal price (no drift): {paths_nd[-1].mean():,.2f}   vs   S0 = {s0:,.2f}")

Mean terminal price (no drift): 588.06   vs   S0 = 578.32

Strengths

• Analytically tractable — closed-form solutions (Black–Scholes) make it the natural baseline
• Positive prices, log-normal terminal distribution — matches the basic stylised fact that prices can't go negative
• Only two parameters — both estimable directly from historical log returns
• Cheap to simulate — vectorised NumPy generates millions of paths in seconds

Weaknesses & Limitations

• Constant volatility — realilty has vol clusters and spikes
• No jumps — crashes like March 2020 are far outside its reach
• Normal log-returns — real returns have fat tails
• Independent increments — momentum and mean-reversion exist

Applications in Practice

• Baseline for option pricing and Monte-Carlo risk engines
• Scenario cones for wealth projections
• Base model against which fancier models must justify their complexity

Alternatives & Extensions

• Stochastic volatility — Heston
• Jump-diffusion — Merton
• GARCH-family models — for clustered volatility