The Heston Model
The Model
The Heston (1993) model is one of the most cited papers in finance, and for good reason: it provides a tractable framework for pricing options under stochastic volatility. In the Black-Scholes world, the volatility of returns is constant, which is at odds with the empirical observation that implied volatilities vary across strikes and maturities—the so-called volatility smile or skew. The Heston model addresses this by allowing the instantaneous variance of returns to evolve as a mean-reverting stochastic process.
Specifically, the stock price S and its instantaneous variance v follow the joint diffusion \begin{aligned} \frac{dS}{S} &= \mu \, dt + \sqrt{v} \, dB_1, \\ dv &= \kappa(\theta - v) \, dt + \sigma \sqrt{v} \, dB_2, \end{aligned} where B_1 and B_2 are correlated Brownian motions satisfying dB_1 \, dB_2 = \rho \, dt. The variance process follows a Cox-Ingersoll-Ross (CIR) process, which keeps the variance non-negative and strictly positive whenever the Feller condition 2\kappa\theta > \sigma^2 holds. The parameter \kappa > 0 governs the speed of mean-reversion of v toward its long-run mean \theta, and \sigma > 0 is the volatility of volatility. The correlation \rho between the two Brownian motions captures the empirically documented leverage effect: in equity markets \rho < 0, so that price declines tend to coincide with increases in variance.
The Stochastic Discount Factor
A natural SDF specification in this two-factor setting is \frac{d\Lambda}{\Lambda} = -r \, dt - \frac{\lambda_1}{\sqrt{v}} \, dB_1 - \frac{\lambda_2}{\sqrt{v}} \, dB_2. The market prices of risk are \lambda_1 / \sqrt{v} for the stock-price risk factor and \lambda_2 / \sqrt{v} for the variance risk factor. Scaling by 1/\sqrt{v} keeps the risk prices bounded as v \to 0 and produces an affine structure under the risk-neutral measure.
By Girsanov’s theorem, define the risk-neutral measure \operatorname{P}^* via d\operatorname{P}^*/d\operatorname{P}= \Lambda_T/\Lambda_0, so that dB_1^* = dB_1 + \frac{\lambda_1}{\sqrt{v}} \, dt, \qquad dB_2^* = dB_2 + \frac{\lambda_2}{\sqrt{v}} \, dt define Brownian motions under \operatorname{P}^* with dB_1^* \, dB_2^* = \rho \, dt. Substituting dB_i = dB_i^* - \lambda_i / \sqrt{v} \, dt into the physical dynamics gives the risk-neutral dynamics of the stock and variance: \begin{aligned} \frac{dS}{S} &= (r - q) \, dt + \sqrt{v} \, dB_1^*, \\ dv &= \kappa(\theta^* - v) \, dt + \sigma \sqrt{v} \, dB_2^*, \end{aligned} where \mu - \lambda_1 = r - q is the stock’s equilibrium risk-premium condition, and \theta^* \equiv \theta - \sigma\lambda_2/\kappa is the risk-neutral long-run variance mean. Under \operatorname{P}^*, every derivative on S earns exactly the risk-free rate in expectation, \operatorname{E}^*(dF) = rF \, dt.
The Partial Differential Equation
Consider a European-style derivative with price F(S, v, T) written on the stock S and expiring at time T. Under the risk-neutral measure \operatorname{P}^* defined by the SDF, every derivative earns exactly the risk-free rate in expectation, \operatorname{E}^*(dF) = rF \, dt. The risk-neutral dynamics of (S, v) are \begin{aligned} \frac{dS}{S} &= (r - q) \, dt + \sqrt{v} \, dB_1^*, \\ dv &= \kappa(\theta^* - v) \, dt + \sigma \sqrt{v} \, dB_2^*, \end{aligned} with dB_1^* \, dB_2^* = \rho \, dt. Applying Ito’s lemma to F(S, v, T) under \operatorname{P}^*: \begin{aligned} dF &= F_S \, dS + F_v \, dv + \tfrac{1}{2} F_{SS} (dS)^2 + \tfrac{1}{2} F_{vv} (dv)^2 + F_{Sv} (dS)(dv) - F_T \, dt \\ &= \left((r-q) S F_S + \kappa(\theta^* - v) F_v + \tfrac{1}{2} v S^2 F_{SS} + \tfrac{1}{2} \sigma^2 v F_{vv} + \sigma v \rho S F_{Sv} - F_T\right) dt \\ &\qquad + F_S S \sqrt{v} \, dB_1^* + F_v \sigma \sqrt{v} \, dB_2^*. \end{aligned} Setting the drift equal to rF yields the Heston PDE: (r - q) S F_S + \kappa(\theta^* - v) F_v + \tfrac{1}{2} v S^2 F_{SS} + \tfrac{1}{2} \sigma^2 v F_{vv} + \sigma v \rho S F_{Sv} - F_T - r F = 0. \tag{1} This PDE, subject to the terminal condition F(S, v, 0) = f(S) for the specific payoff function f, fully characterizes the price of any European-style derivative written on S in the Heston model.
For most payoffs—including vanilla calls and puts—equation (1) cannot be solved in closed form directly. This was Heston’s key insight: rather than solving the PDE for the price itself, one can solve it for the characteristic function of the log stock price, which is available in closed form and can then be used to recover option prices via Fourier inversion.
Option Pricing
The characteristic function of a random variable X is defined as f_X(\phi) = \operatorname{E}(e^{i \phi X}), where i = \sqrt{-1}. Knowing the characteristic function is equivalent to knowing the full distribution of X, and the CDF can be recovered via the Gil-Pelaez inversion formula. For option pricing, if we know the characteristic function of \ln S_T, the in-the-money probabilities can be computed by fast numerical integration.
The price of a European call with strike K and maturity T is C = S e^{-q T} P_1 - K e^{-r T} P_2, where P_1 and P_2 are risk-adjusted probabilities that the option expires in the money, given by \begin{aligned} P_1 &= \frac{1}{2} + \frac{1}{\pi F} \int_{0}^{\infty} \operatorname{Re}\left(\frac{e^{i \phi \ln K} f(\phi - i, 0)}{i \phi}\right) d\phi, \\ P_2 &= \frac{1}{2} + \frac{1}{\pi} \int_{0}^{\infty} \operatorname{Re}\left(\frac{e^{i \phi \ln K} f(\phi, 0)}{i \phi}\right) d\phi, \end{aligned} where F = S e^{(r - q) T} is the forward price and f(\phi, \varphi) is the joint characteristic function of \ln S(T) and v(T). The probability P_2 is the risk-neutral probability that S_T > K, while P_1 is the corresponding probability under the stock-numeraire measure. Unlike Heston (1993), who used two distinct characteristic functions, the formulas above use a single characteristic function evaluated at different arguments by switching to the forward measure.
Define x = \ln S so that the risk-neutral dynamics become \begin{aligned} dx &= (r - q) \, dt + \sqrt{v} \, dB_1^*, \\ dv &= \kappa(\theta^* - v) \, dt + \sigma \sqrt{v} \, dB_2^*, \end{aligned} where B_1^* and B_2^* are correlated Brownian motions under \operatorname{P}^* with dB_1^* \, dB_2^* = \rho \, dt. The joint characteristic function of x(T) and v(T), f(\phi, \varphi) = \operatorname{E}^*\left(e^{i\phi x(T) + i\varphi v(T)}\right), satisfies the same PDE (1) with terminal condition f(\phi, \varphi) = e^{i\phi x + i\varphi v} at T = 0. Because the model is affine, the characteristic function has a closed-form solution: \begin{gathered} f(\phi, \varphi) = \exp\left(i(r - q)\phi T + \left(\delta - \frac{2\gamma}{\sigma^2}\right)\kappa\theta^* T + i\phi x + \delta v\right) \\ \times \left(e^{-\gamma T} - \frac{\sigma^2}{2\gamma}(1 - e^{-\gamma T})(i\varphi - \delta)\right)^{-2\kappa\theta^*/\sigma^2} \\ \times \exp\left(\frac{v(i\varphi - \delta)}{e^{-\gamma T} - \dfrac{\sigma^2}{2\gamma}(1 - e^{-\gamma T})(i\varphi - \delta)}\right), \end{gathered} where \begin{aligned} \gamma &= \sqrt{\kappa^2 + (1 - \rho^2)\sigma^2 \phi^2 + i(\sigma - 2\kappa\rho)\sigma\phi}, \\ \delta &= \frac{\kappa + \gamma - i\rho\sigma\phi}{\sigma^2}. \end{aligned} The expectation is taken under the risk-neutral measure given current values of x = \ln S and v. To price a call option, one only needs f(\phi, 0).
Variance Swaps
A variance swap is a forward contract on realized variance. The holder pays a fixed strike K_{var} and receives the continuously sampled realized variance over [0, T], defined as V_T = \frac{1}{T}\int_0^T v_t \, dt. The payoff at maturity is V_T - K_{var}, so the fair strike is K_{var} = \operatorname{E}^*\!\left[\frac{1}{T}\int_0^T v_t \, dt\right].
The Heston model delivers a closed-form expression for K_{var}. Since v_t follows a CIR process under \operatorname{P}^*, its conditional expectation is \operatorname{E}^*[v_t] = \theta^* + (v_0 - \theta^*)e^{-\kappa t}. Applying Fubini’s theorem gives K_{var} = \theta^* + \frac{v_0 - \theta^*}{\kappa T}(1 - e^{-\kappa T}). \tag{2}
This result also follows directly from the characteristic function. Setting \phi = 0 reduces f(0, \varphi) to the moment-generating function of v(T): the expressions for \gamma and \delta simplify to \gamma|_{\phi=0} = \kappa and \delta|_{\phi=0} = 2\kappa / \sigma^2, so f(0, \varphi) = \operatorname{E}^*[e^{i\varphi v(T)}]. Differentiating with respect to i\varphi and evaluating at \varphi = 0 recovers \operatorname{E}^*[v(T)]; integrating over t \in [0, T] reproduces (2).
When v_0 = \theta^*, the fair strike collapses to the long-run mean, K_{var} = \theta^*. More generally, K_{var} is a weighted average of v_0 and \theta^*, with the weight on v_0 given by (1 - e^{-\kappa T})/(\kappa T). This weight lies in (0, 1], equals one as T \to 0, and decays to zero as T \to \infty, so K_{var} \to \theta^* for long maturities.
Numerical Implementation
The Calibrating the Heston Model notebook covers the full numerical implementation of the option pricing formula, including reliable evaluation of the characteristic function and calibration to observed option prices.