Stochastic Differential Equations

5.1Euler–Maruyama

Most SDEs can't be solved with a formula, so computers solve them step by step: from where you are now, add the predictable part of the change, add a random kick, repeat. This recipe is the Euler–Maruyama method. To check it, we use an SDE that can be solved exactly (a model of a stock price) and compare the two, driven by the same Brownian path.

Drag the number of steps. Left: the exact solution (thick grey) and the step-by-step one (orange). Right: the average error over 500 paths for every step count.

$X_k$

The computed value after $k$ steps (the orange dots). It starts at $X_0 = 1$.

$\Delta t$

The size of one time step: $1/N$ for $N$ steps.

$a(x)$

The drift: the predictable part of the change per unit of time. Here $a(x) = \mu x$ with $\mu = 0.2$: growth of 20% per unit of time, like interest.

$b(x)$

The volatility: how strongly the random kicks act. Here $b(x) = \sigma x$ with $\sigma = 0.8$: kicks proportional to the current value, so the price moves in percentage terms.

$\Delta W_k$

The Brownian path's rise during step $k$: a bell-curve random number with variance $\Delta t$ (section 1.2).

exact solution

For this model Itô's lemma gives a formula: $X_t = e^{(\mu - \sigma^2/2)t + \sigma W_t}$. The $-\sigma^2/2$ is the Itô correction from Chapter 4.

In plain words: each step freezes the drift and the volatility at their values at the start of the step. That is Itô's left-end rule again (section 3.1), which is why this method matches Itô equations. Freezing them makes a small error per step. More, smaller steps shrink the total error, but slowly: the average error falls like $\sqrt{\Delta t}$, so you need 4× the steps to halve it. The dashed line on the right shows that rate.

Right panel: horizontal axis is the number of steps, doubling each tick. The orange dot is your current choice.

Euler–Maruyama converges to the true solution as the steps shrink, with error proportional to $\sqrt{\Delta t}$ (called "strong order ½"). It is the workhorse behind almost every simulation of an SDE, including every other chart on this page.

import numpy as np
rng = np.random.default_rng(13)
mu, sigma, paths, fine = 0.2, 0.8, 2_000, 1024
dW = rng.normal(0, np.sqrt(1 / fine), (paths, fine))
exact = np.exp((mu - sigma ** 2 / 2) + sigma * dW.sum(axis=1))   # exact GBM value at t = 1
for N in (8, 32, 128, 512):
    X = np.ones(paths)
    for incr in dW.reshape(paths, N, -1).sum(axis=2).T:   # N coarse steps, same noise
        X += mu * X / N + sigma * X * incr
    print(f"N = {N:>3}: average error {np.mean(np.abs(X - exact)):.4f}")   # about halves per 4x N

N =   8: average error 0.1595
N =  32: average error 0.0848
N = 128: average error 0.0412
N = 512: average error 0.0203

5.2A Zoo of Diffusions

Every SDE is a choice of drift $a$ and volatility $b$. Three choices appear everywhere. Pick a model and drag its two settings. All three are driven by the same fixed random kicks, so what you see changing is the model, not the luck.

$dX_t$

The tiny change in $X$ over a tiny time step $dt$ (section 4.2).

$a(X_t)\,dt$

Drift: the predictable push, which depends on where $X$ is now.

$b(X_t)\,dW_t$

Noise: a random kick whose size depends on where $X$ is now.

GBM

Geometric Brownian motion, $a = \mu x$, $b = \sigma x$. Grows and wobbles in percentage terms, never below 0. The standard model of a stock price.

OU

Ornstein–Uhlenbeck, $a = \theta(0 - x)$, $b = \sigma$. A pull back toward 0 that gets stronger the further away $X$ is, like a spring. Models velocities, interest rates and anything "mean-reverting".

CIR

Cox–Ingersoll–Ross, $a = \kappa(1 - x)$, $b = \sigma\sqrt{x}$. Pulled toward 1, and the noise fades near 0. Models interest rates and volatility, which should not go negative.

In plain words: the drift decides where the paths want to go, and the volatility decides how wildly they get there. In GBM both scale with the value itself, so paths fan out upward. In OU the spring balances the noise, and the spread stops growing at a fixed width, $\sigma/\sqrt{2\theta}$. In CIR the noise shrinks near 0, so if the pull is strong enough compared with the noise ($2\kappa \ge \sigma^2$, the Feller condition), paths never reach 0 at all.

Grey: 60 of 200 paths (one in orange). Blue bars: where all 200 end at t = 2. Black curve: the exact end-value distribution, shown for GBM and OU where it has a simple formula.

Three one-line equations give three very different kinds of randomness: exploding (GBM), settling (OU) and positive-with-a-floor (CIR). Choosing $a$ and $b$ to match what you are modelling is most of the art of applied stochastic calculus.

import numpy as np
rng = np.random.default_rng(14)
paths, steps, T = 5_000, 400, 2.0
dt = T / steps
models = {   # name: (drift a(x), volatility b(x), start, lowest allowed value)
    "GBM": (lambda x: 0.3 * x, lambda x: 0.4 * x, 1.0, -np.inf),
    "OU": (lambda x: 1.0 * (0 - x), lambda x: 0.6, 2.0, -np.inf),
    "CIR": (lambda x: 1.0 * (1 - x), lambda x: 0.5 * np.sqrt(x), 1.0, 0.0),
}
for name, (a, b, x0, floor) in models.items():
    X = np.full(paths, x0)
    for _ in range(steps):
        X = np.maximum(X + a(X) * dt + b(X) * rng.normal(0, np.sqrt(dt), paths), floor)
    print(f"{name}: mean at t = 2 is {X.mean():.3f}, spread {X.std():.3f}")
# theory: GBM mean e^0.6 = 1.822; OU mean 2e^-2 = 0.271, spread 0.420; CIR mean 1.000

GBM: mean at t = 2 is 1.820, spread 1.147
OU: mean at t = 2 is 0.268, spread 0.415
CIR: mean at t = 2 is 1.000, spread 0.350

5.3Paths ↔ Densities (Fokker–Planck)

There are two ways to describe an SDE. You can follow individual random paths, as above. Or you can describe the whole crowd at once with a density: a curve showing how likely each position is at each time. The density obeys an ordinary (non-random) equation, the Fokker–Planck equation. Here 3,000 particles follow the OU equation from near $x = 2$, while the computer separately solves the Fokker–Planck equation. Scrub the time or press Play.

$p(t, x)$

The density: how crowded position $x$ is at time $t$. The area under the curve between two positions is the fraction of particles between them.

$\frac{\partial p}{\partial t}$

How fast the density at one fixed position is changing over time. The curly ∂ ("partial") means only one variable changes while the others are held fixed.

$\frac{\partial}{\partial x}$, $\frac{\partial^2}{\partial x^2}$

The slope, and the curvature, of a curve drawn across positions at a fixed time.

$a(x)$, $b(x)$

The drift and volatility from section 5.2. Here (OU): $a(x) = -x$, $b(x) = 1$.

In plain words: the first term moves the crowd: drift carries density along like a current, here back toward 0. The second term spreads it: noise smooths out peaks the way heat spreads through a metal bar. Density flows out of crowded, peaked regions (negative curvature) into sparse ones. Together they give one equation that predicts the shape of a crowd of 3,000 random paths without simulating any of them. The black curve (the equation) and the blue bars (the particles) should track each other at every moment.

Left: 40 of the 3,000 particles up to the current time (orange line). Right: histogram of all 3,000 particles now (blue) and the Fokker–Planck solution (black).

Random paths and a smooth equation tell the same story from two sides. This bridge between probability and differential equations runs through all of applied mathematics. The last section of this site, Feynman–Kac, crosses it in the opposite direction.

import numpy as np
rng = np.random.default_rng(15)
theta, sigma, x0, T = 1.0, 1.0, 2.0, 1.0
X = x0 + 0.15 * rng.normal(size=20_000)              # particles
for _ in range(200):                                 # 200 Euler-Maruyama steps of 0.005
    X += -theta * X * 0.005 + sigma * np.sqrt(0.005) * rng.normal(size=X.size)
x = np.linspace(-4, 5, 181)                           # Fokker-Planck on a grid
dx, dt = x[1] - x[0], 0.001
p = np.exp(-0.5 * ((x - x0) / 0.15) ** 2) / (0.15 * np.sqrt(2 * np.pi))
for _ in range(1000):                                # dp/dt = d/dx(theta x p) + 1/2 sigma^2 d2p/dx2
    flux = theta * x * p
    p[1:-1] += dt * ((flux[2:] - flux[:-2]) / (2 * dx) + 0.5 * sigma ** 2 * (p[2:] - 2 * p[1:-1] + p[:-2]) / dx ** 2)
print(f"t = 1: particle mean {X.mean():.3f}, PDE mean {np.sum(x * p) * dx:.3f}, exact {x0 * np.exp(-theta * T):.3f}")

t = 1: particle mean 0.738, PDE mean 0.735, exact 0.736