Browse lessons
On this page

Learn · Lesson

Stochastic Differential Equations: Black–Scholes, Heat Kernels, and Mean Reversion

A stochastic differential equation (SDE) specifies how a state changes through a predictable drift and unpredictable shocks. Solving it gives an entire random trajectory, rather than one deterministic curve. A related partial differential equation (PDE) can describe the distribution of that trajectory or the value of a claim on it.

Learning goals

  • Derive the Black–Scholes PDE using a self-financing hedge.
  • Connect a terminal payoff to Gaussian smoothing through the heat equation.
  • State sufficient conditions for a well-defined, unique SDE solution.
  • Solve geometric Brownian motion and the Ornstein–Uhlenbeck process.
  • Distinguish expected growth, typical path growth, and stationary distributions.
  • Interpret a matrix of diffusion coefficients in a multidimensional system.

1. Reading an SDE correctly

On a filtered probability space, let BtB_t be standard Brownian motion relative to the information filtration Ft\mathcal F_t. The equation

dXt=b(t,Xt) dt+a(t,Xt) dBtdX_t=b(t,X_t)\,dt+a(t,X_t)\,dB_t

means the integral identity

Xt=X0+∫0tb(s,Xs) ds+∫0ta(s,Xs) dBs.X_t=X_0+\int_0^t b(s,X_s)\,ds+\int_0^t a(s,X_s)\,dB_s.

The last term is an Itô integral, not an ordinary derivative multiplied by a small number. Locally, conditional on today’s state, the drift contributes approximately b Δtb\,\Delta t and the shock has variance approximately a2Δta^2\Delta t. Brownian shocks scale like Δt\sqrt{\Delta t}.

The process must use information already available. Future Brownian increments cannot be used to choose today’s coefficient or position.

2. Delta hedging and the Black–Scholes PDE

Assume a non-dividend-paying asset follows

dPt=μPt dt+σPt dBt,dP_t=\mu P_t\,dt+\sigma P_t\,dB_t,

with constant volatility σ>0\sigma>0, constant risk-free rate rr, frictionless continuous trading, and borrowing and lending at that rate. For a sufficiently smooth derivative value G(t,P)G(t,P), Itô’s formula gives

dG=(Gt+μPGP+12σ2P2GPP)dt+σPGP dBt.dG=\left(G_t+\mu P G_P+\frac12\sigma^2P^2G_{PP}\right)dt+\sigma P G_P\,dB_t.

Here GtG_t denotes the time partial derivative; G(t,Pt)G(t,P_t) denotes the value process.

Use a self-financing replicating strategy

Let a portfolio hold Δt\Delta_t shares and a cash balance CtC_t, with value G=ΔP+CG=\Delta P+C. Self-financing means that rebalancing is paid for from the portfolio:

dG=Δ dP+rC dt=Δ dP+r(G−ΔP)dt.dG=\Delta\,dP+rC\,dt=\Delta\,dP+r(G-\Delta P)dt.

Matching the Brownian coefficients requires Δ=GP\Delta=G_P. Matching the drift terms then yields

Gt+μPGP+12σ2P2GPP=μPGP+r(G−PGP).G_t+\mu PG_P+\frac12\sigma^2P^2G_{PP}=\mu PG_P+r(G-PG_P).

Therefore

Gt+rPGP+12σ2P2GPP−rG=0.\boxed{G_t+rPG_P+\frac12\sigma^2P^2G_{PP}-rG=0.}

The real-world drift μ\mu cancels because the hedge removes exposure to the asset’s instantaneous random move.

Equivalently, short one derivative and hold GPG_P shares. Over an instant with the holdings fixed, the hedged gain is

−dG+GP dP=−(Gt+12σ2P2GPP)dt.-dG+G_P\,dP=-\left(G_t+\frac12\sigma^2P^2G_{PP}\right)dt.

Its value is −G+PGP-G+PG_P, so equating its gain to interest on that value gives the same PDE. Do not differentiate −G+PGP-G+PG_P as though rebalancing shares were free: the cash financing must be included.

Terminal and spatial conditions

The contract supplies the terminal payoff:

G(T,P)=(P−K)+orG(T,P)=(K−P)+.G(T,P)=(P-K)^+\quad\text{or}\quad G(T,P)=(K-P)^+.

For a European call, G(t,0)=0G(t,0)=0 and G(t,P)∼P−Ke−r(T−t)G(t,P)\sim P-Ke^{-r(T-t)} as P→∞P\to\infty. A numerical PDE solver also needs appropriate spatial boundary treatment. American exercise and path-dependent claims require additional structure; the displayed terminal-value problem describes European claims depending on the terminal asset price.

3. Why a heat equation appears

The heat equation spreads initial data through a Gaussian kernel. For λ>0\lambda>0,

uτ=λuxx,u(x,0)=h(x),u_\tau=\lambda u_{xx},\qquad u(x,0)=h(x),

has the convolution solution, for suitable initial data,

u(x,τ)=∫−∞∞e−(y−x)2/(4λτ)4πλτh(y) dy.u(x,\tau)=\int_{-\infty}^{\infty}\frac{e^{-(y-x)^2/(4\lambda\tau)}}{\sqrt{4\pi\lambda\tau}}h(y)\,dy.

The fundamental solution for an initial unit point mass is the kernel itself. Its Gaussian variance is 2λτ2\lambda\tau, and its integral over space is one.

For h(x)=1{a<x<b}h(x)=\mathbf1_{\{a<x<b\}},

u(x,τ)=Φ ⁣(b−x2λτ)−Φ ⁣(a−x2λτ).u(x,\tau)=\Phi\!\left(\frac{b-x}{\sqrt{2\lambda\tau}}\right)-\Phi\!\left(\frac{a-x}{\sqrt{2\lambda\tau}}\right).

The denominator is τ\sqrt{\tau} only when λ=1/2\lambda=1/2. For example, λ=1/2\lambda=1/2, τ=1\tau=1, x=0x=0, and (a,b)=(−1,1)(a,b)=(-1,1) give u=2Φ(1)−1≈0.6827u=2\Phi(1)-1\approx0.6827.

An explicit Black–Scholes transformation

Set remaining time τ=T−t\tau=T-t, log price x=log⁡Px=\log P, and

v(τ,x)=G(T−τ,ex),c=r−12σ2.v(\tau,x)=G(T-\tau,e^x),\qquad c=r-\frac12\sigma^2.

The pricing PDE becomes

vτ=12σ2vxx+cvx−rv.v_\tau=\frac12\sigma^2v_{xx}+cv_x-rv.

Discount removal and a moving coordinate eliminate the last two terms: write

v(τ,x)=e−rτu(τ,x+cτ).v(\tau,x)=e^{-r\tau}u(\tau,x+c\tau).

Substitution gives uτ=12σ2uxxu_\tau=\tfrac12\sigma^2u_{xx} with u(0,z)=g(ez)u(0,z)=g(e^z), where gg is the terminal payoff. Thus

G(t,P)=e−rτE[g(Pe(r−σ2/2)τ+στZ)],Z∼N(0,1).G(t,P)=e^{-r\tau}E\left[g\left(Pe^{(r-\sigma^2/2)\tau+\sigma\sqrt{\tau}Z}\right)\right],\qquad Z\sim N(0,1).

This is the discounted risk-neutral expectation. The SDE describes paths, the PDE describes their valuation consequences, and the Gaussian kernel connects both views. See Robert Kohn’s NYU notes on heat equations and financial PDEs for the kernel and transformation approach.

4. Existence and uniqueness: when is the model well-defined?

A useful sufficient theorem assumes measurable deterministic coefficients, an F0\mathcal F_0-measurable initial value with E∣X0∣2<∞E|X_0|^2<\infty, and, uniformly over a finite time interval,

∣b(t,x)−b(t,y)∣2+∣a(t,x)−a(t,y)∣2≤K∣x−y∣2,|b(t,x)-b(t,y)|^2+|a(t,x)-a(t,y)|^2\le K|x-y|^2, ∣b(t,x)∣2+∣a(t,x)∣2≤K(1+∣x∣2).|b(t,x)|^2+|a(t,x)|^2\le K(1+|x|^2).

The first condition is global spatial Lipschitz continuity: nearby starting states cannot have arbitrarily different coefficients. The second is linear growth: coefficient magnitudes grow at most linearly in the state.

Under these conditions there is a continuous adapted strong solution with finite second-moment bounds on each finite horizon. Pathwise uniqueness means that solutions driven by the same Brownian motion and initial value satisfy

P(Xt=Yt for every t∈[0,T])=1.P(X_t=Y_t\text{ for every }t\in[0,T])=1.

“Strong” means the solution is constructed using the given driving noise and initial information. “Finite-horizon bounded” does not mean that its moments are bounded uniformly as T→∞T\to\infty. GBM can have rapidly growing moments while satisfying this theorem. These conditions are sufficient, not necessary: some useful models need different arguments.

5. Geometric Brownian motion by coefficient matching

Consider dXt=μXtdt+σXtdBtdX_t=\mu X_tdt+\sigma X_tdB_t, with X0=x0>0X_0=x_0>0. Try Xt=f(t,Bt)X_t=f(t,B_t). Itô’s formula requires

fx=σf,ft+12fxx=μf.f_x=\sigma f,\qquad f_t+\frac12f_{xx}=\mu f.

The first equation gives f(t,x)=A(t)eσxf(t,x)=A(t)e^{\sigma x}. The second gives A′(t)=(μ−σ2/2)A(t)A'(t)=(\mu-\sigma^2/2)A(t). Hence

Xt=x0e(μ−σ2/2)t+σBt.\boxed{X_t=x_0e^{(\mu-\sigma^2/2)t+\sigma B_t}.}

In particular,

log⁡(Xt/x0)∼N((μ−σ2/2)t,σ2t),\log(X_t/x_0)\sim N((\mu-\sigma^2/2)t,\sigma^2t), E[Xt]=x0eμt,median⁡(Xt)=x0e(μ−σ2/2)t.E[X_t]=x_0e^{\mu t},\qquad\operatorname{median}(X_t)=x_0e^{(\mu-\sigma^2/2)t}.

Coefficient matching is an ansatz, not a universal solution algorithm: many SDE solutions depend on the Brownian history and cannot be expressed as a function of (t,Bt)(t,B_t) alone.

Expected growth versus almost-sure decay

Brownian motion satisfies Bt/t→0B_t/t\to0 almost surely. Therefore

1tlog⁡(Xt/x0)⟶μ−12σ2almost surely.\frac{1}{t}\log(X_t/x_0)\longrightarrow\mu-\frac12\sigma^2\quad\text{almost surely}.

If 0<μ<σ2/20<\mu<\sigma^2/2, then Xt→0X_t\to0 almost surely even though E[Xt]→∞E[X_t]\to\infty. Consequently, for every fixed c>0c>0, P(Xt<c)→1P(X_t<c)\to1.

There is no contradiction: an increasingly rare and extreme upper tail dominates the expectation. Interchanging limit and expectation is invalid here; the family is not uniformly integrable.

For x0=100x_0=100, μ=0.02\mu=0.02, and σ=0.30\sigma=0.30, the log-growth rate is −0.025-0.025. At ten years, the expected value is about 122.14, while the median is about 77.88. These are distribution summaries, not forecasts of a particular realized path. If μ=σ2/2\mu=\sigma^2/2, the zero limiting log-growth rate alone does not imply convergence of XtX_t to zero.

6. Ornstein–Uhlenbeck: solving mean reversion

Use mm for the long-run level to distinguish it from GBM’s drift:

dXt=−α(Xt−m)dt+σdBt,α>0.dX_t=-\alpha(X_t-m)dt+\sigma dB_t,\qquad\alpha>0.

Multiplying Xt−mX_t-m by eαte^{\alpha t} and applying the product rule yields

d(eαt(Xt−m))=σeαtdBt.d\big(e^{\alpha t}(X_t-m)\big)=\sigma e^{\alpha t}dB_t.

Integrating gives

Xt=m+(x0−m)e−αt+σ∫0te−α(t−s)dBs.\boxed{X_t=m+(x_0-m)e^{-\alpha t}+\sigma\int_0^t e^{-\alpha(t-s)}dB_s.}

The deterministic kernel makes the integral Gaussian. Its mean is zero and its variance follows from Itô’s isometry:

E[Xt]=m+(x0−m)e−αt,Var⁡(Xt)=σ22α(1−e−2αt).E[X_t]=m+(x_0-m)e^{-\alpha t},\qquad\operatorname{Var}(X_t)=\frac{\sigma^2}{2\alpha}(1-e^{-2\alpha t}).

Stationarity and the Vasicek interpretation

From a fixed initial value, the distribution approaches N(m,σ2/(2α))N(m,\sigma^2/(2\alpha)). Paths continue fluctuating; they do not settle at mm. Starting with this invariant distribution, independently of future Brownian increments, gives a stationary OU process with covariance

Cov⁡(Xt,Xt+h)=σ22αe−α∣h∣.\operatorname{Cov}(X_t,X_{t+h})=\frac{\sigma^2}{2\alpha}e^{-\alpha|h|}.

The Vasicek short-rate model uses this form for interest rates. The fitted long-run level is a model parameter, not automatically a central bank’s announced target. Physical-measure parameters describe historical dynamics; risk-neutral parameters determine bond pricing. Gaussian rates can become negative.

For x0=0.05x_0=0.05, m=0.03m=0.03, α=0.5\alpha=0.5, and σ=0.01\sigma=0.01, the two-year mean is 3.7358% and standard deviation is 0.9299 percentage points. The expected deviation halves after log⁡(2)/α≈1.3863\log(2)/\alpha\approx1.3863 years.

7. Multidimensional systems and correlated shocks

For dd state variables and kk independent Brownian drivers,

dXt=b(t,Xt)dt+A(t,Xt)dWt,d\mathbf X_t=\mathbf b(t,\mathbf X_t)dt+A(t,\mathbf X_t)d\mathbf W_t,

where AA is a d×kd\times k matrix. The dimensions need not be equal. Instantaneous covariance is

d[Xi,Xj]t=(AAT)ijdt.d[X_i,X_j]_t=(AA^T)_{ij}dt.

Independent drivers can therefore produce correlated states through shared matrix loadings. For example,

dX1=b1dt+σ1dW1,dX2=b2dt+σ2(ρdW1+1−ρ2dW2)dX_1=b_1dt+\sigma_1dW_1,\qquad dX_2=b_2dt+\sigma_2\left(\rho dW_1+\sqrt{1-\rho^2}dW_2\right)

has instantaneous covariance ρσ1σ2dt\rho\sigma_1\sigma_2dt. This construction supports cross-asset and multi-factor interest-rate models.

A position–velocity system illustrates unequal state and noise dimensions:

dQt=Vtdt,dVt=−αVtdt+σdBt.dQ_t=V_tdt,\qquad dV_t=-\alpha V_tdt+\sigma dB_t.

Two state variables share one noise driver. Position depends on the integrated velocity history.

For a smooth function, the generator is

Lf=b⋅∇f+12tr⁡(AAT∇2f).\mathcal L f=\mathbf b\cdot\nabla f+\frac12\operatorname{tr}(AA^T\nabla^2f).

Under suitable regularity and integrability, a discounted terminal expectation solves the backward PDE ft+Lf−rf=0f_t+\mathcal L f-rf=0. In pricing, use the risk-neutral generator; a physical-measure generator generally gives a different expectation.

8. Exact transitions before numerical approximations

For GBM, simulate exactly over a step hh using

Xt+h=Xtexp⁡((μ−σ2/2)h+σhZ).X_{t+h}=X_t\exp((\mu-\sigma^2/2)h+\sigma\sqrt h Z).

For OU, use

Xt+h=m+(Xt−m)e−αh+σ1−e−2αh2αZ,X_{t+h}=m+(X_t-m)e^{-\alpha h}+\sigma\sqrt{\frac{1-e^{-2\alpha h}}{2\alpha}}Z,

with fresh independent standard normals across successive steps. These transitions avoid discretization error at sampled times. They do not reveal every intervening barrier crossing.

For a general SDE, Euler–Maruyama uses Xn+1=Xn+b(tn,Xn)h+a(tn,Xn)hZnX_{n+1}=X_n+b(t_n,X_n)h+a(t_n,X_n)\sqrt h Z_n. Check convergence and stability: a coarse Euler step can make a positive GBM approximation negative or distort OU mean reversion.

Practice questions

  1. Why does μ\mu disappear from the Black–Scholes PDE?
  2. What variance corresponds to ut=3uxxu_t=3u_{xx} after two time units?
  3. What extra scaling is needed in the step-payoff CDF formula for general λ\lambda?
  4. For GBM with μ=0.03\mu=0.03 and σ=0.40\sigma=0.40, what is the almost-sure long-run log-growth rate?
  5. If OU has α=2\alpha=2 and σ=3\sigma=3, what is its invariant variance?
  6. Does convergence to an OU stationary distribution imply paths converge to a constant?
  7. If AA is 5×25\times2, how many states and independent Brownian drivers are present?
  8. Why is differentiating −G+PGP-G+PG_P alone insufficient to justify a dynamic hedge?

Answers

  1. The replicating stock position cancels the asset drift on both sides of the self-financing equation, together with matching its Brownian exposure.
  2. 2λt=122\lambda t=12.
  3. The denominator is 2λt\sqrt{2\lambda t}.
  4. 0.03−0.402/2=−0.050.03-0.40^2/2=-0.05 per year; paths tend toward zero almost surely in this model.
  5. σ2/(2α)=9/4=2.25\sigma^2/(2\alpha)=9/4=2.25.
  6. No. The limiting distribution retains positive variance when σ>0\sigma>0.
  7. Five states and two independent drivers; covariance has rank at most two at each instant.
  8. Changing the share holding requires cash financing. The self-financing gains equation accounts for that rebalancing.

Further study

The MIT lecture on stochastic differential equations provides a companion treatment of hedging, diffusion equations, and OU dynamics. For the underlying calculus and SDE framework, consult NYU’s stochastic calculus course resources. Continue with Volatility Modeling to examine how constant-volatility diffusion assumptions are extended.

Back to Learn