Stochastic differential equations, Itô calculus and the Fokker–Planck equation

Before this lesson: order of convergence and time stepping (lesson 02); Taylor expansion to second order (calculus 03); the Gaussian distribution (probability 05); Monte Carlo estimates and their error bars (probability 14). Langevin dynamics and Markov chains (probability 15) help but are not required. Time: about 110 minutes, plus the notebook.

---

Why this matters #

• Diffusion models are SDEs. The forward process of a diffusion model adds noise along an OU-type SDE. The sampler integrates a reverse-time SDE or a deterministic ODE with the same marginals (lesson 13). Both come from the Fokker–Planck equation of this lesson.
• One-shot training needs the exact solution. Training a diffusion model draws a noisy sample at any time $t$ in one step, with no simulation. That is possible because the OU process has a Gaussian solution, which we derive with Itô's lemma.
• Sampler cost is a question of weak order. A generative sampler must get the distribution right, not one particular noise path. The weak order of the solver tells you how fast the error falls with the number of steps.
• Climate models use stochastic terms. Hasselmann's stochastic climate model is an OU process, and ensemble prediction systems perturb their physics with random patterns that evolve as discretised OU processes.
• Langevin samplers draw from a density $\pi$ by simulating an SDE whose stationary density is $\pi$. The zero-flux condition of the Fokker–Planck equation proves that this works.

1. Predict first #

2. The picture #

Brownian motion: kicks of size $\sqrt h$ #

Brownian motion (also called the Wiener process) $W_t$ is the simplest random path. It starts at $W_0 = 0$, it is continuous, and its increments are independent Gaussians:

So a step of length $h$ moves the path by a random amount with standard deviation $\sqrt h$. For small $h$, $\sqrt h$ is much larger than $h$: for $h = 0.01$, $\sqrt h = 0.1$. This makes the path very rough. It is continuous, but it has no slope anywhere.

Cut $[0, T]$ into $N$ steps of size $h = T/N$ and add up the increments in two ways:

• the total variation $\sum_k |\Delta W_k|$ (the length of the path, seen at resolution $h$). Each term has mean $\sqrt{2h/\pi}$, so the sum has mean $\sqrt{2NT/\pi}$. It grows without limit as $N \to \infty$;
• the quadratic variation $\sum_k (\Delta W_k)^2$. Each term has mean $h$, so the sum has mean $T$. Its variance is $N \cdot 2h^2 = 2T^2/N$, which goes to 0.

So the sum of squared increments does not vanish: it converges to exactly $T$. In a smooth function the squares of small changes are $O(h^2)$ and disappear. For Brownian motion, $(\Delta W)^2$ is of order $h$, the same order as $dt$. In the limit we write

This one rule produces the extra term of Itô's lemma below.

Interactive figure: loading…

An SDE is a limit of a simple recipe #

An SDE in one dimension is written

The function $a$ is the drift: it moves the mean. The function $b$ is the diffusion coefficient: it sets the size of the random kicks. The SDE is defined as the limit of this update, the Euler–Maruyama (EM) method:

Each step adds the mean change $a h$ and the variance $b^2 h$. The coefficients are evaluated at the left point $(X_n, t_n)$, before the noise of the step is known. This choice defines the Itô interpretation of the SDE. With $b = 0$ it is forward Euler for the ODE $x' = a(x, t)$ (lesson 02).

The chain rule gets an extra term #

Take geometric Brownian motion (GBM), $dX = \mu X\,dt + \sigma X\,dW$: the noise is proportional to $X$. If the ordinary chain rule held, $\log X$ would grow at rate $\mu$. It does not. Section 3 shows that $\log X$ grows at rate $\mu - \sigma^2/2$. The mean $E X_t = x_0 e^{\mu t}$ still grows at rate $\mu$, but the median (the typical path) grows only at rate $\mu - \sigma^2/2$. A few large paths pull the mean up.

Interactive figure: loading…

From paths to a density #

You can follow single paths, or you can follow the density $p(x, t)$ of a very large cloud of paths. The density is not random. It obeys a partial differential equation, the Fokker–Planck equation: the drift transports probability (like wind carries smoke), and the noise spreads it (like heat diffusion). An ensemble forecast with many members is a Monte Carlo approximation of this equation.

The simplest example with both effects is the Ornstein–Uhlenbeck process $dX = -\theta X\,dt + \sigma\,dW$. The drift $-\theta X$ pulls $X$ back toward 0 at rate $\theta$; the noise pushes it away. The two balance in a stationary Gaussian with variance $\sigma^2/(2\theta)$.

Interactive figure: loading…

3. Definitions and main results #

Brownian motion and the Itô integral #

Brownian motion $W_t$: $W_0 = 0$; increments over disjoint intervals are independent; $W_t - W_s \sim \mathcal N(0, t - s)$ for $s < t$; paths are continuous. It follows that $\operatorname{Cov}(W_s, W_t) = \min(s, t)$.

A process $g_t$ is non-anticipating if $g_t$ depends only on the path up to time $t$, not on the future. For such $g$, the Itô integral is the limit of left-point sums:

Because $g_{t_k}$ is independent of the next increment, two properties follow:

• mean zero: $E\int_0^T g\,dW = 0$;
• Itô isometry: $E\Big[\big(\int_0^T g\,dW\big)^2\Big] = \int_0^T E[g_t^2]\,dt$.

An Itô SDE $dX = a(X, t)\,dt + b(X, t)\,dW$ is shorthand for the integral equation $X_t = X_0 + \int_0^t a\,ds + \int_0^t b\,dW$.

Theorem 1 (Itô's lemma) #

Let $X$ solve $dX = a\,dt + b\,dW$, and let $f(x, t)$ be once differentiable in $t$ and twice in $x$. Then $Y_t = f(X_t, t)$ solves

Why this is true (intuition). Expand $f$ to second order. In ordinary calculus the second-order term $\tfrac12 f_{xx}(dX)^2$ is $O(h^2)$ and disappears. Here $(dX)^2$ contains $b^2(dW)^2 = b^2\,dt$, which is first order, so half of it survives.

Proof sketch

1. Taylor-expand $f$ over one step: $\Delta f = f_t\,h + f_x\,\Delta X + \tfrac12 f_{xx}\,(\Delta X)^2 + \dots$
2. Insert $\Delta X = a h + b\,\Delta W$: $(\Delta X)^2 = a^2 h^2 + 2ab\,h\,\Delta W + b^2(\Delta W)^2$.
3. The first two terms are $O(h^2)$ and $O(h^{3/2})$ per step. Over $N = T/h$ steps they sum to $O(h)$ and $O(h^{1/2})$, which vanish as $h \to 0$.
4. The last term: $(\Delta W)^2$ has mean $h$ and variance $2h^2$. Summed over $N$ steps with a bounded weight $c_k = \tfrac12 f_{xx}b^2$, the sum $\sum c_k(\Delta W_k)^2$ has mean $\sum c_k h \to \int c\,dt$ and a variance of order $N h^2 = T h \to 0$. So in the limit $(dW)^2$ can be replaced by $dt$.
5. Third-order terms such as $(\Delta W)^3$ are $O(h^{3/2})$ per step and sum to $O(h^{1/2}) \to 0$.
6. Collect the $dt$ terms and the $dW$ term. $\blacksquare$

The rules to remember are $dt^2 = 0$, $dt\,dW = 0$, $(dW)^2 = dt$. With $b = 0$ the lemma is the ordinary chain rule. The new term $\tfrac12 b^2 f_{xx}$ is the Itô correction. For a convex $f$ it is positive: noise pushes $f(X)$ up on average (this is Jensen's inequality, $E f(X) \ge f(E X)$).

Example 1. $f = x^2$ and $X = W$ ($a = 0$, $b = 1$): $f_x = 2x$, $f_{xx} = 2$, so $d(W^2) = 2W\,dW + dt$. Taking expectations, $\frac{d}{dt}E[W_t^2] = 1$ and $E[W_t^2] = t$.

Example 2: geometric Brownian motion. $dX = \mu X\,dt + \sigma X\,dW$. Take $f = \log x$, so $f_x = 1/x$, $f_{xx} = -1/x^2$:

The right side has constant coefficients, so we can integrate it directly:

The median of $X_t$ is $x_0e^{(\mu - \sigma^2/2)t}$ (the median of $W_t$ is 0). The mean is $x_0e^{\mu t}$, because $E\,e^{\sigma W_t} = e^{\sigma^2 t/2}$ (the moment generating function of $\mathcal N(0, t)$).

Theorem 2 (the Ornstein–Uhlenbeck process) #

The OU process $dX = -\theta X\,dt + \sigma\,dW$ with mean-reversion rate $\theta > 0$, noise level $\sigma > 0$ and $X_0 = x_0$ has the solution

Proof

1. Use an integrating factor, as for the linear ODE $x' = -\theta x$: let $Y = e^{\theta t}X$, that is $f(x, t) = e^{\theta t}x$ with $f_t = \theta e^{\theta t}x$, $f_x = e^{\theta t}$, $f_{xx} = 0$.
2. Itô's lemma: $dY = \big(\theta e^{\theta t}X - \theta X e^{\theta t}\big)dt + \sigma e^{\theta t}\,dW = \sigma e^{\theta t}\,dW$. There is no correction term, because $f$ is linear in $x$.
3. Integrate from 0 to $t$: $e^{\theta t}X_t - x_0 = \sigma\int_0^t e^{\theta s}\,dW_s$. Multiply by $e^{-\theta t}$.
4. The integral of a deterministic function against $dW$ is a limit of sums of independent Gaussians, so it is Gaussian, with mean 0.
5. Its variance, by the Itô isometry, is $\sigma^2\int_0^t e^{-2\theta(t-s)}\,ds = \frac{\sigma^2}{2\theta}(1 - e^{-2\theta t})$. $\blacksquare$

As $t \to \infty$ the mean goes to 0 and the variance to $\sigma^2/(2\theta)$: the stationary distribution is $\mathcal N\big(0, \sigma^2/(2\theta)\big)$. The time scale of the memory is $1/\theta$.

The same argument over one step of length $\Delta$ gives the exact transition

So an OU process sampled every $\Delta$ is an AR(1) process (time series 05) with coefficient $e^{-\theta\Delta}$, and this sampling has no step-size error. With a mean level $\mu$, $dX = \theta(\mu - X)\,dt + \sigma\,dW$, replace $x$ by $x - \mu$ everywhere.

Theorem 3 (Fokker–Planck equation) #

Let $X$ solve $dX = a(x, t)\,dt + b(x, t)\,dW$ with smooth $a$ and $b$, and let $X_t$ have a smooth density $p(x, t)$. Then

The first term is advection (the drift carries probability), the second is diffusion (the noise spreads it). Note that $b^2$ stays inside the second derivative.

Proof (from Itô's lemma)

1. Take any smooth test function $\varphi(x)$ that is zero outside a bounded set. By Itô's lemma, $d\varphi(X) = \big(a\varphi' + \tfrac12 b^2\varphi''\big)dt + b\varphi'\,dW$.
2. Take expectations. The $dW$ term has mean zero: $\frac{d}{dt}E\varphi(X_t) = E\big[a\varphi' + \tfrac12 b^2\varphi''\big]$.
3. Write both sides as integrals against the density: $\int \varphi\,\partial_t p\,dx = \int \big(a\varphi' + \tfrac12 b^2\varphi''\big)\,p\,dx$.
4. Integrate by parts, once for the $\varphi'$ term and twice for the $\varphi''$ term. The boundary terms vanish because $\varphi$ does: $\int \varphi\,\partial_t p\,dx = \int \varphi\,\big[-\partial_x(a p) + \tfrac12\partial_{xx}(b^2 p)\big]\,dx$.
5. This holds for every test function $\varphi$, so the two functions in brackets are equal. $\blacksquare$

The operator $\mathcal A\varphi = a\varphi' + \tfrac12 b^2\varphi''$ in step 2 is the generator of the SDE. The Fokker–Planck operator is its adjoint: integration by parts moves the derivatives from $\varphi$ to $p$. This is the same "move the derivative to the other factor" idea as in the adjoint methods of lesson 03.

Fokker–Planck in $d$ dimensions

For $dX = a(X, t)\,dt + B(X, t)\,dW$ with $X \in \R^d$, $W \in \R^m$ and $B$ a $d \times m$ matrix, let $D = BB^\top$. Then

$$\partial_t p = -\nabla\cdot(a\,p) + \tfrac12\sum_{i,j}\partial_i\partial_j\big(D_{ij}\,p\big).$$

For the diffusion-model case $B = g(t)\,I$ (the same noise in every direction, independent of $x$) this is $\partial_t p = -\nabla\cdot(a\,p) + \tfrac12 g^2\Delta p$, where $\Delta$ is the Laplacian. Lesson 13 starts from this form.

Stationary densities and probability flux #

Write the 1D equation as a conservation law, $\partial_t p = -\partial_x J$, with the probability flux

$J(x)$ is the rate at which probability crosses the point $x$ from left to right. A stationary density has $\partial_x J = 0$, so $J$ is a constant. On the whole real line, with $p$ and $J$ decaying in the tails, that constant is 0. This is the zero-flux condition:

• OU ($a = -\theta x$, $b = \sigma$): $-\theta x\,p = \tfrac{\sigma^2}{2}p'$. So $p'/p = -2\theta x/\sigma^2$ and $p \propto e^{-\theta x^2/\sigma^2}$, which is $\mathcal N\big(0, \sigma^2/(2\theta)\big)$, as in Theorem 2.
• Overdamped Langevin dynamics $dX = -U'(X)\,dt + \sqrt2\,dW$ with a potential $U$: zero flux gives $-U'p = p'$, so $p \propto e^{-U}$, the Gibbs density. With $U = -\log \pi$ the SDE $dX = (\log\pi)'(X)\,dt + \sqrt2\,dW$ has $\pi$ as its stationary density. It needs only the derivative of $\log \pi$, not the normalising constant of $\pi$. This derivative is the score, the central object of lesson 13.
• With multiplicative noise ($b$ depends on $x$), solve the zero-flux equation for $q = b^2p$: $q' = (2a/b^2)\,q$, a first-order linear ODE.

Interactive figure: loading…

Itô or Stratonovich?

The Stratonovich integral $\int g \circ dW$ evaluates $g$ at the midpoint of each step instead of the left point. It obeys the ordinary chain rule, but its mean is not zero in general. Example: $\int_0^T W\,dW = \tfrac12(W_T^2 - T)$ (Itô) while $\int_0^T W\circ dW = \tfrac12W_T^2$ (Stratonovich).

A Stratonovich SDE $dX = a\,dt + b\circ dW$ is the same process as the Itô SDE with drift $a + \tfrac12 b\,b_x$. Physical noise with a short but nonzero correlation time leads to the Stratonovich form (the Wong–Zakai theorem). Euler–Maruyama converges to the Itô solution, so convert the drift first. For additive noise ($b$ does not depend on $x$, so $b_x = 0$) the two forms are the same. All diffusion models use additive noise.

Strong and weak convergence #

Let $X_N$ be a numerical solution at time $T = Nh$ and $X(T)$ the exact solution driven by the same Brownian path.

• Strong order $\gamma$: $E|X_N - X(T)| \le C h^\gamma$. This measures pathwise accuracy.
• Weak order $q$: $|E f(X_N) - E f(X(T))| \le C h^q$ for smooth test functions $f$. This measures how well the distribution is reproduced.

Scheme	Update	Strong order	Weak order
Euler–Maruyama	$X_n + a h + b\,\Delta W$	1/2 (1 if $b_x = 0$)	1
Milstein	EM $+\ \tfrac12\,b\,b_x\big((\Delta W)^2 - h\big)$	1	1

Where the Milstein term comes from (Itô–Taylor expansion)

Over one step, $X_{t_{n+1}} = X_n + \int a\,ds + \int b(X_s)\,dW_s$. EM freezes $b$ at $X_n$. Expand once more: by Itô's lemma, $b(X_s) \approx b(X_n) + b_x b\,(W_s - W_{t_n})$ to leading order. The extra term is $b\,b_x\int_{t_n}^{t_{n+1}}(W_s - W_{t_n})\,dW_s$. By Example 1 applied to the increment $W_s - W_{t_n}$, this integral equals $\tfrac12\big((\Delta W)^2 - h\big)$.

Why EM has strong order 1/2 but weak order 1. The term that EM drops, $\tfrac12 b\,b_x((\Delta W)^2 - h)$, has mean zero and standard deviation $\tfrac{1}{\sqrt2}|b\,b_x|\,h$ per step.

1. Pathwise, $N = T/h$ such errors with mean zero add like a random walk: the total is of size $\sqrt N\cdot h = \sqrt{Th} = O(h^{1/2})$.
2. In expectation they cancel. What remains is a bias of $O(h^2)$ per step, so the total is $N\cdot h^2 = O(h)$: weak order 1.
3. Milstein keeps the term. Its leftover pathwise error per step is $O(h^{3/2})$, and the strong order rises to 1. Its weak order stays 1.
4. For additive noise ($b$ constant) the term is zero. EM then already is Milstein and has strong order 1.

Interactive figure: loading…

Measuring the orders. For the strong error, every step size must see the same Brownian path: draw the finest increments once and sum blocks of them for the coarser grids. With independent paths the error would not go to zero; it would measure the spread of the process. For the weak error, the Monte Carlo noise $s/\sqrt M$ ($M$ paths, sample standard deviation $s$) must be well below the bias.

4. Worked example #

(a) One step of GBM. Let $x_0 = 1$, $\mu = 0.05$, $\sigma = 0.2$, $h = 0.25$, and suppose the Brownian increment is $\Delta W = 0.3$.

1. EM: the drift term is $\mu x_0 h = 0.05\cdot 0.25 = 0.0125$ and the noise term is $\sigma x_0\,\Delta W = 0.2\cdot 0.3 = 0.06$. So the EM value after one step is $1 + 0.0125 + 0.06 = 1.0725$.
2. Milstein term: $b = \sigma x$, so $b_x = \sigma$ and $b\,b_x = \sigma^2 x_0 = 0.04$. With $(\Delta W)^2 - h = 0.09 - 0.25 = -0.16$, the term is $\tfrac12\cdot 0.04\cdot(-0.16) = -0.0032$. Milstein: $1.0725 - 0.0032 = 1.0693$.
3. Exact (Itô's lemma): the exponent is $(0.05 - 0.02)\cdot 0.25 + 0.2\cdot 0.3 = 0.0675$, so $X(h) = e^{0.0675} = 1.06983$.
4. Errors: EM $2.7\cdot10^{-3}$, Milstein $5.3\cdot10^{-4}$. The Milstein term removed about 80 % of the EM error on this path.
5. Over $T = 1$: the mean of $X_1$ is $e^{0.05} = 1.0513$ and the median is $e^{0.03} = 1.0305$. The Itô correction $-\sigma^2/2 = -0.02$ lowers the growth rate of $\log X$.

(b) An OU process. Let $\theta = 2$, $\sigma = 1$, $x_0 = 1$.

1. At $t = 0.5$: the mean is $e^{-2\cdot 0.5} = e^{-1} = 0.368$. The variance is $\tfrac{1}{4}(1 - e^{-2})$, that is $0.25\cdot 0.8647 = 0.216$.
2. Stationary: variance $\sigma^2/(2\theta) = 1/4$, standard deviation $1/2$; the memory time is $1/\theta = 0.5$.
3. Exact transition with $\Delta = 0.1$: $e^{-0.2} = 0.819$ and $\sqrt{(1 - e^{-0.4})/4} = 0.287$, so $X_{n+1} = 0.819\,X_n + 0.287\,Z$ with $Z \sim \mathcal N(0, 1)$. This is exact for any $\Delta$.
4. Fokker–Planck: $a = -2x$, $b = 1$, so $\partial_t p = \partial_x(2x\,p) + \tfrac12\partial_{xx}p$.
5. Zero flux: $-2x\,p = \tfrac12 p'$, so $p'/p = -4x$ and $p \propto e^{-2x^2}$. Compare with $e^{-x^2/(2v)}$: $1/(2v) = 2$, so $v = 1/4$. ✓

5. In code #

The worked example, a finite-difference check that the OU density solves the Fokker–Planck equation, the quadratic variation, and the weak error of EM:

import numpy as np
from scipy.stats import chi2

# (a) One step of GBM: Euler-Maruyama, Milstein and the exact solution on the same increment
x0, mu, sigma, h, dW = 1.0, 0.05, 0.2, 0.25, 0.3
em = x0 + mu * x0 * h + sigma * x0 * dW
mil = em + 0.5 * (sigma * x0) * sigma * (dW**2 - h)        # b b_x = (sigma x) sigma
exact = x0 * np.exp((mu - sigma**2 / 2) * h + sigma * dW)  # Ito's lemma on log X
print(f"{em:.4f} {mil:.4f} {exact:.5f}")                    # 1.0725 1.0693 1.06983
print(f"{abs(em - exact):.1e} {abs(mil - exact):.1e}")      # 2.7e-03 5.3e-04

# (b) OU with theta = 2, sigma = 1, x0 = 1
theta, sig, x0 = 2.0, 1.0, 1.0
t = 0.5
print(f"{x0 * np.exp(-theta * t):.3f} {sig**2 / (2 * theta) * (1 - np.exp(-2 * theta * t)):.3f}")  # 0.368 0.216
D = 0.1
print(f"{np.exp(-theta * D):.3f} {np.sqrt(sig**2 / (2 * theta) * (1 - np.exp(-2 * theta * D))):.3f}")  # 0.819 0.287

# The OU density N(m(t), v(t)) solves p_t = d/dx(theta x p) + (sigma^2/2) p_xx: finite-difference residual
x = np.linspace(-3, 3, 601); dx = x[1] - x[0]; dt = 1e-5
def p(t):
    m, v = x0 * np.exp(-theta * t), sig**2 / (2 * theta) * (1 - np.exp(-2 * theta * t))
    return np.exp(-(x - m)**2 / (2 * v)) / np.sqrt(2 * np.pi * v)
pt = (p(t + dt) - p(t - dt)) / (2 * dt)
rhs = np.gradient(theta * x * p(t), dx) + sig**2 / 2 * np.gradient(np.gradient(p(t), dx), dx)
print(f"max |p_t| = {abs(pt).max():.2f}, max residual = {abs(pt - rhs).max():.1e}")
# max |p_t| = 0.86, max residual = 5.4e-04

# (c) Quadratic variation: the sum of N squared increments is h * chi2(N)
for N in (10, 1_000, 100_000):
    hN = 1.0 / N
    print(f"N = {N:6d}: QV mean {chi2(N, scale=hN).mean():.3f}, QV std {chi2(N, scale=hN).std():.4f}, "
          f"E[TV] {N * np.sqrt(2 * hN / np.pi):6.1f}")
# N =     10: QV mean 1.000, QV std 0.4472, E[TV]    2.5
# N =   1000: QV mean 1.000, QV std 0.0447, E[TV]   25.2
# N = 100000: QV mean 1.000, QV std 0.0045, E[TV]  252.3

# (d) Weak error of EM in the mean of GBM: E[X_{n+1}] = (1 + mu h) E[X_n] exactly
mu, T = 2.0, 1.0
errs = np.array([abs((1 + mu * T / N)**N - np.exp(mu * T)) for N in (16, 32, 64, 128)])
print(np.round(errs, 4), np.round(errs[:-1] / errs[1:], 3))
# [0.8058 0.4304 0.2228 0.1134] [1.872 1.932 1.965]

The finite-difference residual ($5\cdot10^{-4}$) is the discretisation error of the check itself; it falls when you refine dx. The weak-error ratios approach 2: weak order 1. The strong orders, measured on coupled paths:

rng = np.random.default_rng(0)
mu, sigma, T, M, NF = 1.0, 1.0, 1.0, 4000, 512
dWf = np.sqrt(T / NF) * rng.standard_normal((M, NF))        # one set of fine increments
XT = np.exp((mu - sigma**2 / 2) * T + sigma * dWf.sum(1))    # exact X(T) on every path
hs, e_em, e_mil = [], [], []
for N in (8, 16, 32, 64, 128, 256):
    dW = dWf.reshape(M, N, -1).sum(2)                        # the same paths on a coarser grid
    h = T / N
    xe = np.ones(M); xm = np.ones(M)
    for n in range(N):
        d = dW[:, n]
        xe = xe + mu * xe * h + sigma * xe * d
        xm = xm + mu * xm * h + sigma * xm * d + 0.5 * sigma**2 * xm * (d**2 - h)
    hs.append(h); e_em.append(np.mean(abs(xe - XT))); e_mil.append(np.mean(abs(xm - XT)))
slope = lambda e: np.polyfit(np.log(hs), np.log(e), 1)[0]
print(f"strong order: EM {slope(e_em):.2f}, Milstein {slope(e_mil):.2f}")
# strong order: EM 0.49, Milstein 0.96

6. Where it shows up #

• Diffusion models. The variance-preserving forward SDE $dX = -\tfrac12\beta(t)X\,dt + \sqrt{\beta(t)}\,dW$ is an OU process with a time-dependent rate. The proof of Theorem 2 gives $X_t \mid X_0 \sim \mathcal N\big(\alpha_t X_0, (1 - \alpha_t^2)I\big)$ with $\alpha_t = \exp\big(-\tfrac12\int_0^t\beta(s)\,ds\big)$. Training draws $x_t$ from this formula in one step (generative models 07).
• The probability-flow ODE. For $dX = a\,dt + g(t)\,dW$ in $d$ dimensions, the diffusion term can be written as a flux: $\tfrac12 g^2\Delta p = \nabla\cdot\big(\tfrac12 g^2\,p\,\nabla\log p\big)$. So the Fokker–Planck equation becomes a pure transport equation with velocity $a - \tfrac12 g^2\nabla\log p$, and the ODE with this velocity has the same marginals as the SDE. This is the probability-flow ODE of lesson 13.
• Langevin samplers. EM applied to $dX = \nabla\log\pi\,dt + \sqrt2\,dW$ is the unadjusted Langevin algorithm (ULA). Its stationary distribution is off by $O(h)$, a weak error. MALA removes this bias with a Metropolis–Hastings step (probability 15); score-based samplers use the same update with a learned score (generative models 06).
• Climate. Hasselmann's stochastic climate model (1976) treats the slow ocean temperature as integrating fast weather noise, $dT = -\lambda T\,dt + \sigma\,dW$: an OU process. Its spectrum is "red" (large at low frequencies), as observed for sea-surface temperature. Ensemble weather prediction perturbs the tendencies of the physics parametrisations with random patterns that evolve in time as AR(1) processes, that is, discretised OU processes (the SPPT scheme: Buizza et al. 1999; Palmer et al. 2009).
• Seismology. The Kanai–Tajimi model of earthquake ground acceleration is the output of a damped oscillator (the soil layer) driven by white noise: a linear SDE in two variables, whose stationary covariance follows from the same Itô and Fokker–Planck tools. In earthquake forecasting, the Brownian passage time model (Matthews et al. 2002) describes the time between large earthquakes on a fault as the first time a Brownian motion with drift reaches a failure threshold.
• Strong or weak? Diffusion sampling and climate statistics need only the distribution: weak order. Pathwise questions (a given noise realisation, coupled comparisons, multilevel Monte Carlo) need strong order.
• Finance. GBM is the Black–Scholes model of a stock price; the $-\sigma^2/2$ in its log-return is Itô's lemma.

7. Common mistakes #

• Forgetting the Itô term. For example, writing $d\log X = \mu\,dt + \sigma\,dW$ for GBM. The correct drift is $\mu - \sigma^2/2$.
• Increments with standard deviation $h$. A Brownian increment over a step $h$ has standard deviation $\sqrt h$: np.sqrt(h) * rng.standard_normal().
• Evaluating the coefficients at the new point or the midpoint in EM. The new point makes the scheme implicit, the midpoint changes the limit toward Stratonovich. Itô means the left point.
• Mixing Itô and Stratonovich models. Convert the drift ($a \to a + \tfrac12 b\,b_x$) before you use EM on a Stratonovich SDE.
• Measuring the strong error with independent paths. The numerical and the exact solution must use the same Brownian increments; otherwise the "error" does not go to zero.
• Weak errors hidden by Monte Carlo noise. You need $s/\sqrt M$ well below $Ch$; compare schemes on the same paths to cancel most of the noise.
• Expecting Milstein to help with additive noise. $b_x = 0$, so Milstein equals EM, which then already has strong order 1.
• OU stationary variance $\sigma^2/\theta$. The correct value is $\sigma^2/(2\theta)$.
• Writing $\tfrac12 b^2\,\partial_{xx}p$ in the Fokker–Planck equation instead of $\tfrac12\partial_{xx}(b^2p)$. The two differ when $b$ depends on $x$.

8. Check yourself #

1. Compute $d(W_t^3)$ with Itô's lemma and use it to show $E[W_t^3] = 0$.
   Answer

   $f = x^3$, $f_x = 3x^2$, $f_{xx} = 6x$: $d(W^3) = 3W^2\,dW + 3W\,dt$. Taking expectations, $\frac{d}{dt}E[W_t^3] = 3E[W_t] = 0$, and $W_0^3 = 0$, so $E[W_t^3] = 0$ for all $t$.
2. For GBM, find $E[X_t^2]$ with Itô's lemma.
   Answer

   $f = x^2$: $d(X^2) = (2\mu X^2 + \sigma^2X^2)\,dt + 2\sigma X^2\,dW$. So $\frac{d}{dt}E X^2 = (2\mu + \sigma^2)\,E X^2$ and $E X_t^2 = x_0^2e^{(2\mu + \sigma^2)t}$. The variance is $x_0^2e^{2\mu t}(e^{\sigma^2t} - 1)$.
3. Write the Fokker–Planck equation of $dX = -X^3\,dt + \sqrt2\,dW$ and its stationary density.
   Answer

   $\partial_t p = \partial_x(x^3p) + \partial_{xx}p$. Zero flux: $-x^3p - p' = 0$, so $p \propto e^{-x^4/4}$. This is Langevin dynamics with $U = x^4/4$.
4. A study on shared Brownian paths gives errors $0.040$, $0.028$, $0.020$ for $h$, $h/2$, $h/4$. Which order is this, and which scheme and noise type could produce it?
   Answer

   Each ratio is about $\sqrt2$, so the order is $1/2$. This is the strong error of Euler–Maruyama with multiplicative noise ($b_x \ne 0$), for example on GBM.
5. Why does Milstein give no gain over EM for the OU process?
   Answer

   The noise is additive: $b = \sigma$ is constant, $b_x = 0$, so the Milstein term vanishes. EM is then the same scheme as Milstein and already has strong order 1.
6. Is $M_t = W_t^2 - t$ a martingale (a process with no drift)? Is $e^{\sigma W_t}$?
   Answer

   $dM = 2W\,dW + dt - dt = 2W\,dW$: no drift, so yes. For $e^{\sigma W_t}$ Itô gives the drift $\tfrac12\sigma^2e^{\sigma W_t} \ne 0$, so no. The process $e^{\sigma W_t - \sigma^2t/2}$ is a martingale: its $f_t$ term cancels the Itô term.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. What is the distribution of a Brownian increment $W_{t+h} - W_t$?

   Show answer

2. What does the quadratic variation $\sum_k(\Delta W_k)^2$ of Brownian motion on $[0, T]$ converge to, and what rule does this give?

   Show answer

3. State Itô's lemma for $Y = f(X, t)$ with $dX = a\,dt + b\,dW$.

   Show answer

4. Write one Euler–Maruyama step for $dX = a(X, t)\,dt + b(X, t)\,dW$.

   Show answer

5. Why does the median of GBM grow at rate $\mu - \sigma^2/2$ and not $\mu$?

   Show answer

6. What is the stationary distribution of the OU process $dX = -\theta X\,dt + \sigma\,dW$?

   Show answer

7. Write the 1D Fokker–Planck equation of $dX = a\,dt + b\,dW$.

   Show answer

8. What is the stationary density of $dX = -U'(X)\,dt + \sqrt 2\,dW$, and how do you get it?

   Show answer

9. Strong vs weak order: what does each measure, and what are they for Euler–Maruyama?

   Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Why is there a $-\sigma^2/2$ in the solution of GBM?" Itô's lemma on $\log X$: $f_{xx} = -1/x^2$ gives $-\tfrac12\sigma^2$. The mean grows at rate $\mu$, the median at rate $\mu - \sigma^2/2$.
• "What is the Fokker–Planck equation, and how is it used in diffusion models?" The PDE for the density of an SDE. Rewritten as a continuity equation, it gives the probability-flow ODE with the same marginals.
• "Strong or weak convergence: which one matters for a generative sampler?" Weak: you care about the distribution of the samples, not about matching one particular noise path.
• "Itô or Stratonovich?" Itô: left point, zero-mean integrals, EM converges to it. Stratonovich: midpoint, ordinary chain rule, the limit of smooth noise. Convert with $a \to a + \tfrac12 b\,b_x$; for additive noise they agree.
• "How do you test an SDE solver?" Moments against closed forms, a strong-order study on coupled paths, a weak-order study with error bars, and stationary histograms against the Fokker–Planck solution.

10. Practice #

The notebook for this lesson runs in Colab or in the browser. It has five exercises, each with a check cell and a folded solution:

1. Brownian paths from Gaussian increments; total and quadratic variation as the grid is refined.
2. Euler–Maruyama and Milstein for a general scalar SDE; check them on the worked example.
3. A convergence study on GBM with coupled paths: measured strong orders near 0.5 (EM) and 1.0 (Milstein), and the weak order of EM in $E[X_T]$.
4. Itô's lemma as a function: drift and diffusion of $f(X, t)$; check it on $\log$ of GBM and by Monte Carlo for $e^{W_t}$.
5. Exact OU sampling with the Gaussian transition; the long-run histogram against $\mathcal N(0, \sigma^2/2\theta)$.

11. Watch and read #

18. Itō CalculusMIT OpenCourseWare · 78 min

Itô's lemma from the second-order Taylor expansion and the rule (dW)² = dt (MIT 18.S096).

21. Stochastic Differential EquationsMIT OpenCourseWare · 56 min

What it means to solve an SDE, and how Itô's lemma solves simple ones.

Ito Integration Clearly and Visually ExplainedRoman Paolucci · 32 min

A visual introduction to the Itô integral as a limit of left-point sums.

12. Brownian Motion (Wiener Process), Ornstein-Uhlenbeck Process, Einstein-Stokes RelationPhysicsOfLifeLMU · 35 min

Brownian motion and the OU process from the physics side (Langevin equation, Einstein relation).

MIT 6.S184: Flow Matching and Diffusion Models - Lecture 01 - Flow and Diffusion Models (2026)Peter Holderrieth · 78 min

ODEs, SDEs, Euler–Maruyama and the OU process as the building blocks of diffusion models.

• Särkkä and Solin, Applied Stochastic Differential Equations (free PDF, 2019): §4.1–4.3 (Itô integral and Itô formula, the OU solution), §4.6 (Stratonovich), §5.2 (Fokker–Planck–Kolmogorov equation, stationary solutions), §8.2–8.3 (Euler–Maruyama, Milstein, strong and weak order).
• D. J. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Review 43 (2001): §4 (EM), §5 (strong and weak convergence), §6 (Milstein), §8 (stochastic chain rule). Short, with MATLAB code for every figure.
• Holderrieth and Erives, MIT 6.S184 lecture notes (2026): §2.2 (SDEs, the OU process, Euler–Maruyama), §4.2 (Fokker–Planck and Langevin dynamics); lab one of the course labs simulates the same processes.