Geometric Brownian motion

Prices change in percent

Arithmetic Brownian motion adds the same dollar noise whether a stock costs 5 or 500, and it can go below zero. A stock behaves differently: a 1% move is a 1% move at any price. Geometric Brownian motion puts the randomness in the percentage change: dS = mu * S * dt + sigma * S * dW. It is the model of Black–Scholes. mu is the expected return per year and sigma the volatility, the quantity quoted in the market.

import math

def gbm_step(s, mu, sigma, dt, z):
    return s * math.exp((mu - 0.5 * sigma ** 2) * dt + sigma * math.sqrt(dt) * z)

print(round(gbm_step(100.0, 0.08, 0.2, 1 / 252, 0.0), 4))     # a calm day: only the (mu - sigma^2/2) drift acts
print(round(gbm_step(100.0, 0.08, 0.2, 1 / 252, 1.0), 4))     # a one-sigma up day
print(round(gbm_step(100.0, 0.08, 0.2, 1 / 252, -50.0), 4))   # an absurd fall: still positive

The exact solution, and the minus one half sigma squared

Solving the equation gives S(t) = S(0) * exp((mu - sigma**2 / 2) * t + sigma * W(t)). Why the strange - sigma**2 / 2? Because the exponential is convex: by Jensen's inequality, noise raises the average of exp(noise), and subtracting sigma**2 / 2 exactly cancels it, so the *expected* price grows at rate mu. Typical *paths* grow more slowly, at mu - sigma**2 / 2, which is why the median stock price is below the mean. The next topic, Itô's lemma, derives this term.

import numpy as np

def gbm_path(s0, mu, sigma, dt, z):
    log_moves = (mu - 0.5 * sigma ** 2) * dt + sigma * np.sqrt(dt) * z
    return s0 * np.exp(np.concatenate([[0.0], np.cumsum(log_moves)]))

rng = np.random.default_rng(1)
z = rng.standard_normal(252)
path = gbm_path(100.0, 0.08, 0.2, 1 / 252, z)
print(path[:4].round(3), round(float(path[-1]), 3))

Many paths: the mean and the median

Simulate thousands of one-year paths and look at where they end. The *average* final price is close to S0 * exp(mu * T), the formula for the expected price. The *median* is lower, S0 * exp((mu - sigma**2/2) * T), and with a high volatility the gap is large: most paths finish below the average, pulled up by a few spectacular ones. This right-skew, a lognormal distribution, is the shape of stock prices over long horizons.

import numpy as np

rng = np.random.default_rng(5)
s0, mu, sigma, T, steps, paths = 100.0, 0.08, 0.4, 1.0, 252, 20000
dt = T / steps
z = rng.standard_normal((paths, steps))
final = s0 * np.exp(np.sum((mu - 0.5 * sigma ** 2) * dt + sigma * np.sqrt(dt) * z, axis=1))

print(round(float(final.mean()), 1), round(s0 * np.exp(mu * T), 1))                      # close to each other: the expected price
print(round(float(np.median(final)), 1), round(s0 * np.exp((mu - 0.5 * sigma ** 2) * T), 1))   # the median is lower
print(round(float(np.mean(final < s0)), 3))                                              # the share of paths that lose money

Estimating mu and sigma

From prices sampled every dt years, the log returns r = diff(log(prices)) are independent normals with mean (mu - sigma**2/2) * dt and standard deviation sigma * sqrt(dt). So sigma = std(r) / sqrt(dt) and mu = mean(r) / dt + sigma**2 / 2. The volatility estimate is good with a few hundred points; the drift is notoriously noisy, because its standard error shrinks only with the total *time*, not with the number of observations.

import numpy as np

def estimate_gbm(prices, dt):
    r = np.diff(np.log(prices))
    sigma = r.std(ddof=1) / np.sqrt(dt)
    return r.mean() / dt + 0.5 * sigma ** 2, sigma

rng = np.random.default_rng(8)
dt = 1 / 252
for years in (1, 10, 100):
    n = int(years * 252)
    log_moves = (0.10 - 0.5 * 0.25 ** 2) * dt + 0.25 * np.sqrt(dt) * rng.standard_normal(n)
    prices = 100 * np.exp(np.concatenate([[0.0], np.cumsum(log_moves)]))
    mu, sigma = estimate_gbm(prices, dt)
    print(years, "years:", round(float(mu), 3), round(float(sigma), 3))           # true mu 0.10 and sigma 0.25

• GBM: S(t) = S(0) * exp((mu - sigma**2 / 2) * t + sigma * W(t)); percentage changes are random and prices stay positive.
• The - sigma**2 / 2 cancels the convexity of exp: the expected price is S0 * exp(mu * t), but the median path grows at mu - sigma**2 / 2.
• Estimate from log returns: sigma = std(r) / sqrt(dt), mu = mean(r) / dt + sigma**2 / 2; volatility is estimated well, drift poorly.

Next up: Itô's lemma, the chain rule that explains where the - sigma**2 / 2 comes from.

Lab · what's graded

• Take one exact step with gbm_step25 pts
• Simulate a price path with gbm_path25 pts
• Average over all futures with gbm_expected25 pts
• Read the parameters back with estimate_gbm25 pts