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).
Numerical Implementation
Implementing the Heston formula reliably requires care. The most important consideration is to write every term in the characteristic function using e^{-\gamma T} rather than e^{\gamma T}. The complex square root \gamma has branches that can cause sign flips as \phi varies, and the formulation with e^{-\gamma T} — which decays for large T — avoids the exponential blow-up that plagues the equivalent e^{\gamma T} formulation. Failing to observe this leads to discontinuities in the integrand and wildly incorrect option prices.
A reliable way to verify the implementation is to exploit the moment-generating property of the characteristic function. Generate a large number of Monte Carlo paths for x(T) and v(T) under \operatorname{P}^*, then check that the sample average of e^{i\phi x(T) + i\varphi v(T)} matches the closed-form formula for several values of \phi and \varphi. This test should be performed across a range of parameter values to expose numerical instabilities. The same simulation can also verify the final option prices by computing the average of discounted payoffs and comparing to the Fourier-inversion formula.