Applications

6.1Black–Scholes by Hedging

A call option is a contract that pays $\max(S_T - K,\,0)$ in one year: the amount by which the stock price $S$ ends above an agreed price $K$, or nothing. What is a fair price for it today? Black, Scholes and Merton's answer: sell the option, then trade the stock so that your holdings always move exactly like the option. If that works, the risk cancels and the fair price is whatever it costs to run this hedge.

Here you sell the option for the Black–Scholes price and hedge it along 1,000 simulated stock paths. Drag how often you rebalance.

$S_t$

The stock price, modelled as geometric Brownian motion (section 5.2), starting at 100 with 20% volatility.

$K$

The strike: the agreed price, here 100.

$C(t, S)$

The Black–Scholes value of the option at time $t$ when the stock is at $S$ (thick grey line).

$V_t$

Your hedge portfolio: the money received for the option, plus gains from trading the stock (orange line).

$\Delta_t$

"Delta": how many shares to hold, equal to how much the option's value changes when the stock moves by 1. $N(d_1)$ is its formula, a bell-curve area like $P(W_t \ge a)$ in section 2.3.

$dV_t = \Delta_t\,dS_t$

Your portfolio changes by (shares held) × (change in stock price): an Itô integral (section 3.2), with the stock as the "price".

In plain words: by Itô's lemma (Chapter 4), the option's value changes by $\Delta \times$ (the stock's change), plus a predictable, non-random term in $dt$. Holding $\Delta$ shares cancels the random part exactly. Black and Scholes chose the price so that the predictable part cancels too. Hedge continuously and the portfolio ends up exactly equal to the option's payoff on every path. In practice you rebalance only now and then, so a small error is left, and it shrinks as you trade more often. Notice that the stock's expected growth (here 10% a year) never enters the price: the hedge removes it.

Left: one stock path. Vertical lines mark rebalancing dates (hidden when frequent). Right: final profit/loss of "sold option + hedge" across all 1,000 paths. Interest rate set to 0 for simplicity.

The more often you rebalance, the tighter the profit/loss clusters around 0, and the orange portfolio tracks the grey option value. In the limit, the option's risk is replicated away completely, so its price must be the cost of the hedge: the Black–Scholes formula. This argument earned Scholes and Merton the 1997 Nobel prize in economics (Black had died in 1995).

import numpy as np
def Phi(z):                                          # bell-curve area N(z), Abramowitz & Stegun 7.1.26
    t = 1 / (1 + 0.3275911 * np.abs(z) / np.sqrt(2))
    poly = t * (0.254829592 + t * (-0.284496736 + t * (1.421413741 + t * (-1.453152027 + t * 1.061405429))))
    erf = 1 - poly * np.exp(-z * z / 2)
    return 0.5 * (1 + np.sign(z) * erf)

S0 = K = 100.0
sigma, mu, days, paths = 0.2, 0.1, 252, 2_000

def call(S, tau):                                    # Black-Scholes value and delta, rate 0
    d1 = (np.log(S / K) + 0.5 * sigma ** 2 * tau) / (sigma * np.sqrt(tau))
    return S * Phi(d1) - K * Phi(d1 - sigma * np.sqrt(tau)), Phi(d1)

rng = np.random.default_rng(16)
steps = (mu - sigma ** 2 / 2) / days + sigma * np.sqrt(1 / days) * rng.normal(size=(paths, days))
S = S0 * np.exp(np.concatenate([np.zeros((paths, 1)), np.cumsum(steps, axis=1)], axis=1))
price = float(call(S0, 1.0)[0])
for rebalances in (4, 21, 252):
    every, V = days // rebalances, np.full(paths, price)
    for i in range(days):
        if i % every == 0:
            delta = call(S[:, i], 1 - i / days)[1]
        V += delta * (S[:, i + 1] - S[:, i])        # gains from holding delta shares
    pnl = V - np.maximum(S[:, -1] - K, 0)
    print(f"price {price:.3f}, rebalance {rebalances:>3}x a year: P&L spread {pnl.std():.3f}")

price 7.966, rebalance   4x a year: P&L spread 3.315
price 7.966, rebalance  21x a year: P&L spread 1.448
price 7.966, rebalance 252x a year: P&L spread 0.445

6.2Change of Measure (Girsanov)

These paths drift upward at rate $\lambda$ ("lambda"): $X_t = W_t + \lambda t$. There are two ways to remove a drift. You can move the paths, or you can keep the paths and change how likely each one is. Tick "reweight": each path gets a weight, shown by how strongly it is drawn. Paths that drifted up a lot become rare, and paths that stayed low become common.

$P$

The original probabilities: every simulated path equally likely.

$Q$

New probabilities for the same set of paths: a re-weighting.

$\frac{dQ}{dP}$

The weight each path gets: how much more (or less) likely it is under $Q$ than under $P$. It is called the Radon–Nikodym derivative. Weights average to 1, so they only shift likelihood between paths.

$W_1$

Where the underlying Brownian path ends at $t = 1$. The weight depends only on that.

$\lambda$

The drift rate, set by the slider.

In plain words: a path that ends high is typical when there's an upward drift, and rare when there isn't. Girsanov's theorem gives the exact weight that makes the drifting paths look statistically identical to drift-free Brownian motion: the average path goes flat, and the end values follow $N(0,1)$. Nothing about any individual path changes. Only the odds change.

Left: 80 of 2,000 paths. With reweighting on, more strongly drawn means more weight. Orange: the (weighted) average path. Right: (weighted) histogram of end values with the matching bell curve. For large drifts a handful of paths carry most of the weight, so the weighted numbers get noisier.

Changing probabilities can remove a drift exactly. In finance this is risk-neutral pricing. Real stocks drift upward, but reweighting the paths so the stock grows only at the risk-free interest rate lets you price an option as a plain weighted average of its payoff. It gives the same answer as the hedging argument of section 6.1, without building a hedge.

import numpy as np
rng = np.random.default_rng(17)
lam, paths = 1.0, 200_000
W1 = rng.normal(size=paths)                          # where each Brownian path ends at t = 1
X1 = W1 + lam                                        # the drifting paths end here
Z = np.exp(-lam * W1 - lam ** 2 / 2)                 # Girsanov weight dQ/dP
mean_q = np.mean(Z * X1)
print(f"plain mean {X1.mean():.3f} (lambda)   weighted mean {mean_q:.3f} (0)   mean weight {Z.mean():.3f} (1)")
print(f"weighted variance {np.mean(Z * X1 ** 2) - mean_q ** 2:.3f} (1)")

plain mean 1.002 (lambda)   weighted mean 0.007 (0)   mean weight 0.997 (1)
weighted variance 0.983 (1)

6.3Feynman–Kac

The colour map shows the solution $u(t, x)$ of the heat equation, the equation for how heat spreads along a bar: stronger blue means larger. The equation is fixed by its value at the final time, $u(1, x) = g(x)$: a triangular bump centred on $x = 0$. Click anywhere on the map. Random Brownian paths start from that point, and the average of $g$ at their end points is plotted on the right.

$u(t, x)$

The unknown function, with one value for each time $t$ (horizontal) and position $x$ (vertical): the colour map.

$g(x)$

The known final values at $t = 1$: $g(x) = \max(0,\,1 - \lvert x\rvert)$, a triangle of height 1 over $-1 \le x \le 1$.

$x + W_{1-t}$

Where a Brownian path started at position $x$ at time $t$ ends up at time 1, after running for the remaining time $1 - t$.

$\mathbb{E}[\;\cdot\;]$

The average over all such random paths (section 2.2).

$\frac{\partial u}{\partial t} + \frac12 \frac{\partial^2 u}{\partial x^2} = 0$

A partial differential equation (PDE): the rate of change of $u$ in time is tied to its curvature in space (section 5.3's notation). This one is the heat equation run backward from $t = 1$.

In plain words: to find the value of the PDE's solution at any point, start random walkers there, let them wander until the final time, look up $g$ where they land, and average. No grid and no solving are needed: randomness does the calculus. It also works the other way round: the average of a random quantity can be computed by solving a PDE. With a GBM in place of $W$, a call option's payoff as $g$, and discounting at the interest rate, this is the Black–Scholes formula.

Left: the solution (stronger blue = larger, from 0 to 1), with 25 of the random paths from your point in orange. Right: running average over 3,000 paths (log scale), against the PDE's value (dashed). Keyboard: Tab to the map and use the arrow keys to move the point.

Wherever you click, the average over random paths converges to the value the PDE gives there. This is the Feynman–Kac formula, the reverse of section 5.3. Fokker–Planck turns paths into a density equation. Feynman–Kac turns a PDE into an average over paths. It is how banks price complex derivatives by Monte Carlo simulation, and it ties together everything on this site: Brownian motion (Ch 1), its rough-but-regular behaviour (Ch 2), the Itô integral and lemma (Ch 3–4) and SDEs (Ch 5).

import numpy as np
from math import sqrt
rng = np.random.default_rng(18)
g = lambda x: np.maximum(0, 1 - np.abs(x))           # final condition u(1, x) = g(x)
t, x = 0.3, 0.8
monte_carlo = g(x + sqrt(1 - t) * rng.normal(size=200_000)).mean()
xs = np.linspace(-6, 6, 601)                         # heat equation, solved backward from t = 1
dx = xs[1] - xs[0]
dt = 0.4 * dx ** 2
u = g(xs)
for _ in range(round((1 - t) / dt)):
    u[1:-1] += dt * 0.5 * (u[2:] - 2 * u[1:-1] + u[:-2]) / dx ** 2
print(f"Monte Carlo {monte_carlo:.4f}   PDE {np.interp(x, xs, u):.4f}")

Monte Carlo 0.2963   PDE 0.2965