The Heston Model: Two-State SDE to Semi-Analytical Pricing

What keeps the model solvable once volatility gets its own random driver: Heston's affine structure, the full derivation of the characteristic function, the numerics of Fourier inversion, and where the five-parameter calibration turns ill-posed.

9 of 12 ~15 min

The three defects of local vol in the previous article — forward smiles too flat, spot-vol correlation pinned at $-1$ , no independent volatility risk — are really one sentence: the model lacks a second source of randomness. Adding one is easy. Adding one and still being able to compute is the hard part: bolt an arbitrary stochastic process onto $\sigma$ and the resulting two-dimensional model usually leaves you with brute-force Monte Carlo. The combination Steven Heston found in 1993 is this: let variance follow a mean-reverting square-root process, and the model lands in the affine class — the characteristic function comes out in closed form, and European pricing reduces to a one-dimensional integral. The tractability is not luck; it is structure. This article walks that structural line: how the model is assembled, how the closed form is derived, how the numerics are done, where calibration turns ill-posed, and what the model still gets wrong.

The Set-Up

The SDE system:

dS_t = r S_t \, dt + \sqrt{v_t} \, S_t \, dW^S_t

dv_t = \kappa(\theta - v_t) \, dt + \sigma \sqrt{v_t} \, dW^v_t

with $\text{Corr}(dW^S, dW^v) = \rho$ . Each of the five parameters has one job: $\kappa$ is the mean-reversion speed of variance (units $1/\text{time}$ ); $\theta$ is the long-run level variance reverts to; $\sigma$ is the vol of vol — often written $\xi$ or $\eta$ in the literature to avoid a clash with the BS $\sigma$ (it is a model parameter, not to be confused with the second-order Greek volga); $\rho$ is the correlation between stock and variance innovations; and $v_0$ is the initial variance.

The drift of the variance equation is an Ornstein-Uhlenbeck-style pullback acting on $v$ rather than on the process itself; the $\sqrt{v_t}$ in the diffusion shrinks the noise as $v$ approaches zero, and under Feller’s condition $2\kappa\theta > \sigma^2$ the process stays strictly positive. The model is usually stated directly under the risk-neutral measure $\mathbb{Q}$ ; the market price of variance risk only needs to appear explicitly if you want to go back to the physical measure $\mathbb{P}$ for scenario analysis — for pricing you can ignore it.

The parameter $\rho$ deserves its own paragraph. In equity index markets $\rho < 0$ : the index falls, variance rises. This is the leverage effect — the economic mechanism is still debated (behavioral stories and “price falls, financial leverage rises” stories both have adherents), but the empirical fact is solid, and it is exactly the dynamical source of the negatively skewed smile. FX is different: major currency pairs have roughly symmetric returns and symmetric smiles, and the same Heston model with $\rho \approx 0$ describes them well. The negative correlation is where Heston earns the equity skew and separates itself from a plain mixture-of-lognormals.

The Pricing Idea: Characteristic Functions

The model is set; how to price? Heston’s route bypasses numerical PDEs and goes through the frequency domain. Define the characteristic function of $\ln S_T$ under the risk-neutral measure:

\phi(u, t) = \mathbb{E}^{\mathbb{Q}}\left[e^{iu \ln S_T}\right]

with $t$ the time to expiry, so $\phi(u, 0) = 1$ . The characteristic function carries all the information of the distribution; if $\phi$ is available in closed form, the price of any European payoff can be recovered by a Fourier-type inversion. For a call, one convenient form (Lewis-style, along the contour $\text{Im}(u) = -1/2$ ) is:

C = S_0 - \frac{\sqrt{K}\, e^{-rT}}{\pi} \int_0^\infty \frac{\text{Re}\left[e^{-iu\ln K}\, \phi(u - i/2, T)\right]}{u^2 + 1/4}\, du

(Here $\phi$ already contains the $\ln S_0 + rT$ drift; expanding $\phi(u - i/2)$ releases a factor $\sqrt{S_0 e^{rT}}$ , which combines with $\sqrt{K}\, e^{-rT}$ into the symmetric prefactor $\sqrt{S_0 K}\, e^{-rT/2}$ — the dimensions work out.) The integrand decays like $u^{-2}$ , and standard numerical integration handles it in under a millisecond per strike.

So all the difficulty concentrates in one place: getting $\phi$ in closed form. This is where the affine structure enters.

The Backward PDE for $\phi$

As a conditional expectation, $\phi(u, t)$ satisfies the Kolmogorov backward equation:

\frac{\partial \phi}{\partial t} = \mathcal{L}\phi

where $\mathcal{L}$ is the infinitesimal generator of the joint $(S, v)$ process, under the risk-neutral measure:

\mathcal{L} = r S \frac{\partial}{\partial S} + \frac{1}{2} v S^2 \frac{\partial^2}{\partial S^2} + \kappa(\theta - v)\frac{\partial}{\partial v} + \frac{1}{2}\sigma^2 v \frac{\partial^2}{\partial v^2} + \rho \sigma v S \frac{\partial^2}{\partial S \partial v}

That is a two-dimensional PDE, no cheaper to solve directly than Monte Carlo. The decisive step is guessing the shape of the solution — try the affine ansatz $\phi = \exp(iu\ln S + A(u, t) + B(u, t) v)$ : the state variables $\ln S$ and $v$ appear only linearly in the exponent. The reason the guess is legitimate: every coefficient of the generator is linear in $v$ (note that $S^2\partial^2_S$ acting on $e^{iu\ln S}$ produces only terms linear in $v$ ) — that is the defining feature of an “affine model,” and all of Heston’s tractability flows from it.

Substituting, the PDE splits by powers of $v$ into two ODEs. The terms linear in $v$ give the equation for $B$ :

\frac{dB}{dt} = -\frac{1}{2}(u^2 + iu) - (\kappa - i\rho\sigma u) B + \frac{1}{2}\sigma^2 B^2

This is the Heston Riccati equation. Each of the three terms has a pedigree: $-\frac{1}{2}(u^2 + iu)$ comes from the log-price diffusion and the Itô drift correction; $(\kappa - i\rho\sigma u)B$ packages the mean reversion of variance together with the leverage effect (this is where $\rho$ enters); and $\frac{1}{2}\sigma^2 B^2$ is the vol-of-vol contribution. The constant terms give the equation for $A$ :

\frac{dA}{dt} = \kappa\theta B + riu

$B$ is a standard Riccati; $A$ is just an integral of $B$ . Initial conditions sit at $t = 0$ (expiry): $A(u, 0) = 0$ , $B(u, 0) = 0$ , and you integrate forward in time-to-expiry $t$ .

The Closed Form

The Riccati equation solves explicitly. Define the discriminant:

d = \sqrt{(\rho\sigma iu - \kappa)^2 + \sigma^2(iu + u^2)}

and the ratio:

g = \frac{\kappa - \rho\sigma iu - d}{\kappa - \rho\sigma iu + d}

with $|g| < 1$ on the relevant domain. The solution for $B$ (in the numerically stable $-d$ form):

B(u, t) = \frac{\kappa - \rho\sigma iu - d}{\sigma^2} \cdot \frac{1 - e^{-dt}}{1 - g e^{-dt}}

Heston’s original paper writes the algebraically equivalent form with $\tilde{g} = 1/g$ and growing exponentials:

B(u, t) = \frac{\kappa - \rho\sigma iu + d}{\sigma^2} \cdot \frac{1 - e^{dt}}{1 - \tilde{g}\, e^{dt}}

The two are equal on paper and not equivalent in floating point: $e^{dt}$ grows exponentially with maturity, and combined with the branch choice of the complex logarithm it produces discontinuities at long maturities — the most famous trap in Heston implementations. The fix is to use the $-d$ form throughout (Lord & Kahl, 2010, treat the complex-logarithm branch issue completely). The solution for $A$ (same convention):

A(u, t) = riu\, t + \frac{\kappa\theta}{\sigma^2}\left[(\kappa - \rho\sigma iu - d)t - 2\log\left(\frac{1 - ge^{-dt}}{1 - g}\right)\right]

Putting it together:

\phi(u, t) = \exp\left(A(u, t) + B(u, t) v_0 + iu \ln S_0\right)

This is the Heston characteristic function. Call pricing is valid on the strip $\text{Im}(u) \in [-1, 0]$ ; for puts, use parity or shift the contour.

Numerical Inversion and Code

Besides the Lewis form above, Heston (1993) originally used the probability decomposition:

C = S_0 P_1 - K e^{-rT} P_2, \qquad P_j = \frac{1}{2} + \frac{1}{\pi}\int_0^\infty \text{Re}\left[\frac{e^{-iu\ln K}\,\phi_j(u, T)}{iu}\right] du

where $P_2$ is the risk-neutral exercise probability ( $\phi_2$ is the $\phi$ above) and $P_1$ the exercise probability under the stock numéraire ( $\phi_1(u) = \phi(u - i)/\phi(-i)$ ). Authors differ in contours and normalization; the content is equivalent. The integrand decays like $u^{-2}$ , so a few hundred quadrature points give six digits; truncating at $u \leq 200$ is usually ample.

A runnable implementation (stable $-d$ form, checked against Black-Scholes in the $\sigma \to 0$ limit):

import numpy as np
from scipy.integrate import quad

def heston_phi(u, S0, r, T, v0, kappa, theta, sigma, rho):
    d = np.sqrt((rho * sigma * 1j * u - kappa)**2
                + sigma**2 * (1j * u + u**2))
    g = (kappa - rho * sigma * 1j * u - d) / (kappa - rho * sigma * 1j * u + d)
    e = np.exp(-d * T)
    B = (kappa - rho * sigma * 1j * u - d) / sigma**2 * (1 - e) / (1 - g * e)
    A = (1j * u * (np.log(S0) + r * T)
         + kappa * theta / sigma**2
         * ((kappa - rho * sigma * 1j * u - d) * T
            - 2 * np.log((1 - g * e) / (1 - g))))
    return np.exp(A + B * v0)

def heston_call(S0, K, r, T, v0, kappa, theta, sigma, rho):
    """One-dimensional integral on the Lewis contour Im(u) = -1/2.
    phi(u - i/2) carries a sqrt(S0 e^{rT}) factor, so the prefactor
    reduces to sqrt(K) e^{-rT}."""
    def integrand(u):
        return np.real(np.exp(-1j * u * np.log(K))
                       * heston_phi(u - 0.5j, S0, r, T, v0, kappa, theta, sigma, rho)
                       ) / (u**2 + 0.25)
    val, _ = quad(integrand, 0, 200, limit=400)
    return S0 - np.sqrt(K) * np.exp(-r * T) / np.pi * val

Under fifty lines, and European options under Heston are done. Calling quad strike by strike is fine for research; to price a whole chain, use FFT-based methods (Carr & Madan, 1999) or the COS method (Fang & Oosterlee, 2008) — both below.

Calibration

With a pricer in hand, the next step is inverting the market: optimize over $(v_0, \kappa, \theta, \sigma, \rho)$ to fit market IVs. Five parameters against hundreds of quotes sounds safely overdetermined; in practice the problem is ill-conditioned in characteristic ways. The procedure itself is standard: pick a calibration set (say weekly expiries, five strikes each), define the loss $\text{RMSE} = \sqrt{\frac{1}{N}\sum_i (\sigma_{\text{mkt}}^i - \sigma_{\text{Heston}}^i)^2}$ , and optimize with Levenberg-Marquardt locally or differential evolution globally.

Where does it hurt? First, identifiability: different parameter combinations grow nearly identical surfaces. The ATM level is set mainly by $\theta$ and $v_0$ , the skew mainly by $\rho$ , the wing steepness by the product of $\sigma$ and $\rho$ — near the optimum the parameters are highly collinear. Second, multiple local minima; without a global search you will settle in the wrong basin. Third, $\kappa$ is the weakest-identified of the five: the data cannot pin it down, the realistic equity range is roughly $0.5$ to $5$ , and short-dated options barely feel it — which is the calibration-side reflection of the fact that $\kappa$ drops out of expected variance entirely when $v_0 = \theta$ .

A typical SPX calibration starts from $\kappa=2, \theta=0.04, \sigma=0.3, \rho=-0.7, v_0=0.04$ . After convergence, the residuals usually concentrate in deep OTM puts — exactly where the model is least flexible: a continuous-path diffusion cannot produce a fat enough short-dated left tail. Fixing that requires jumps (Merton, 1976), which is the next article’s subject.

What Heston Itself Gets Wrong

Heston cures half of local vol’s disease and leaves three of its own, in order of severity.

First, skew dynamics. $\rho$ is a constant, so the relationship between skew steepness and the vol level is welded shut inside the model. Empirically, skew (suitably measured) tends to be flatter in low-vol periods and steeper in high-vol periods, and the co-movement itself changes. A constant $\rho$ cannot express that.

Second, the forward smile. Heston’s forward smile is closer to reality than local vol’s (it does not flatten indefinitely), but parameters calibrated to today’s vanilla surface still produce forward-starting option prices that miss the market quotes for those products — the forward smile dynamics remain too tame.

Third, vol-of-vol moves too. The parameter $\sigma$ is constant, while in reality the volatility of volatility is low in calm periods and spikes under stress. Heston has none of that layer; rough volatility and multi-factor Bergomi-type models exist to fix exactly this.

Beyond Vanilla: Simulation

European options go through the characteristic function; path-dependent products — Asians, barriers, lookbacks — force you back onto the paths themselves: either a two-dimensional PDE in $(S, v)$ (feasible for low-dimensional payoffs, but slow) or Monte Carlo. The first lesson of simulating Heston: naive Euler-Maruyama dies. Discretized, the square-root diffusion samples negative $v$ , and $\sqrt{v}$ produces NaNs. Production schemes are designed around this. Andersen’s (2008) Quadratic-Exponential (QE) scheme branches on the size of $v$ — for large $v$ , a moment-matched quadratic (squared Gaussian) approximation; for small $v$ , a moment-matched exponential distribution, non-negative by construction; the same paper gives a martingale-corrected variant that removes the discretization bias in the forward. Alfonsi’s (2005) implicit scheme is another positivity-preserving route.

def simulate_heston_qe(S0, r, T, v0, kappa, theta, sigma, rho, n_paths, n_steps):
    """QE scheme skeleton (Andersen, 2008). psi <= psi_c takes the quadratic
    branch, otherwise the exponential branch."""
    dt = T / n_steps
    psi_c = 1.5
    S = np.full(n_paths, S0)
    v = np.full(n_paths, v0)
    for step in range(n_steps):
        # conditional moments of v (closed form for the CIR process)
        m  = theta + (v - theta) * np.exp(-kappa * dt)
        s2 = (v * sigma**2 * np.exp(-kappa * dt) / kappa * (1 - np.exp(-kappa * dt))
              + theta * sigma**2 / (2 * kappa) * (1 - np.exp(-kappa * dt))**2)
        psi = s2 / m**2
        # quadratic branch: v' = a (b + Z)^2 with a, b moment-matched to (m, s2)
        # exponential branch: v' = 0 with prob p, else exponential sample, moment-matched
        # log-price step: Andersen's K0..K4 coefficients wire in the integrated
        # variance and the correlation
        ...
    return S

A full implementation is about a hundred lines of careful numerical work, each step straightforward; QuantLib and similar libraries ship ready versions. The branch constants and the $K_0 \ldots K_4$ coefficients of the log-price step are all tabulated in Andersen’s paper.

Engineering the Inversion: Damping and COS

Beyond throwing quadrature points at the pricing integral, two pieces of engineering are worth knowing.

First, damping. The call payoff itself is not integrable (it grows linearly in $S_T$ ); Carr and Madan (1999) handle this by multiplying the price by $e^{\alpha k}$ ( $k = \ln K$ ) before transforming, and the resulting transform has a closed form:

\psi(u) = \frac{e^{-rT}\,\phi\big(u - (\alpha + 1)i,\; T\big)}{\alpha^2 + \alpha - u^2 + i(2\alpha + 1)u}

Invert $\psi$ and strip off $e^{\alpha k}$ to recover the price. The choice of $\alpha$ controls the decay and oscillation of the integrand; chosen well, the integration range shrinks dramatically; chosen badly (especially at long maturity and extreme strikes) you hit complex-log branch trouble and underflow — one more reason the characteristic function must be implemented in the stable $-d$ form.

Second, the COS method (Fang & Oosterlee, 2008): expand the density in a cosine basis on a truncated interval $[a, b]$ , and the price becomes a fast-converging series:

V(S_0, K, T) = e^{-rT} \sum_{n=0}^{N-1}{}' \; \text{Re}\left[\phi\!\left(\frac{n\pi}{b-a}\right) e^{i n\pi \frac{x - a}{b - a}}\right] V_n, \qquad x = \ln(S_0 / K)

where $V_n$ are the cosine coefficients of the payoff (closed form for vanillas), and the prime halves the first term. The characteristic-function evaluations do not depend on the strike, so one pass feeds an entire chain; the coefficient structure also vectorizes with FFTs. Typical settings: $N = 1024$ for 4 digits, $N = 4096$ for 6, milliseconds for a full chain.

import numpy as np

def cos_heston_chain(S0, K_array, r, T, v0, kappa, theta, sigma, rho, N=512, L=10):
    """COS method skeleton (Fang & Oosterlee, 2008)."""
    # truncation interval [a, b]: set from the first cumulants of ln(S_T/K);
    # L controls the coverage width
    a, b = -L * np.sqrt(T) * 0.3, L * np.sqrt(T) * 0.3   # schematic; use the cumulant formula
    x = np.log(S0 / K_array)
    n = np.arange(N)
    u_n = n * np.pi / (b - a)
    phi_n = np.array([heston_phi(u, 1.0, r, T, v0, kappa, theta, sigma, rho)
                      for u in u_n])                      # strike-independent, computed once
    V_n = cos_payoff_coeffs_call(a, b, N)                 # vanilla payoff coefficients, closed form
    w = np.ones(N); w[0] = 0.5
    series = (w * np.real(phi_n[None, :] * np.exp(1j * np.outer(x - a, u_n))) * V_n).sum(axis=1)
    return np.exp(-r * T) * K_array * series

COS is the production standard for whole-chain pricing under affine models; its only requirement is a known characteristic function, so it covers Heston, Bates, and Lévy models alike.

Smile Dynamics

A calibrated Heston does not just give today’s prices; it implicitly answers “how does the smile move when spot moves?” Two reference frames: sticky-strike — spot moves, IV at each fixed strike stays put; sticky-delta / sticky-moneyness — spot moves, IV at each moneyness stays put, the whole smile translating with spot. Heston, being spot-homogeneous, has sticky-moneyness as its baseline: only $S/K$ means anything inside the model, so the smile naturally floats with spot. Local vol is the opposite — its smile responds to spot in a sticky-strike or even contrarian way. The market itself sits between the two and shifts with the regime — trending markets look more sticky-strike, panic repricings more sticky-delta. Which model you hedge with changes the delta by a “smile correction term”; this is the direct source of the hedge-ratio disagreement between the two model families. Extensions with stochastic rates (Heston–Hull-White and relatives) improve the long-maturity dynamics further, but the basic picture stands.

An Easy Trap: Backing Out Local Vol

A common operation: take the calibrated parametric surface (Heston or otherwise) as input, back out local vol via the Dupire formula, and feed it to Monte Carlo. A distinction matters here. The surface generated by the Heston model itself is arbitrage-free by construction, so the implied $\sigma^2_{\text{loc}}$ is non-negative in theory. But in a real pipeline what you invert is usually not the model surface — it is a market fit that has been interpolated and extrapolated, and when the interpolation scheme carries no no-arbitrage constraint, $\sigma^2_{\text{loc}} < 0$ pops up in grid corners. A negativity check after calibration should be a fixed station in the pipeline. What to do when it fires depends on the use: refit with a more flexible parametric form, cap the local vol at a floor, switch to Bates which accommodates steep smiles natively — or, if you only price Europeans, simply trust the Heston characteristic function, which is immune to these interpolation artifacts.

Why It Is Still Alive After 25 Years

Alternatives keep arriving: SABR, local-stochastic vol hybrids, Bergomi, rough vol. Heston remains the default calibration on most desks for practical reasons. Semi-analytical pricing is fast — orders of magnitude faster than Monte Carlo or finite differences, which is decisive inside a calibration loop. Five parameters are reasonably well-conditioned most of the time, and the toolchain is mature. For short-dated vanilla hedging books it performs no worse than many “better” models. There is also a teaching dividend: junior quants learn the affine framework through Heston, and every later model is incremental on top. The next generation of systems may run on rough Bergomi or proxy models, but Heston will live in legacy code for decades.

Nor is the modern desk’s reality “pick one model” — it runs a portfolio of them in parallel: Heston for vanilla speed, Bates (Heston plus jumps) for skewed exotics, SABR for rate vol surfaces, local vol as the deterministic fit feeding Monte Carlo, Bergomi or rough vol for the dynamics of volatility risk. Each model holds one section of the front; the combination covers the product line. The single-model era ended long ago.

The next two articles pick up the two loose ends Heston leaves: jump-diffusion (to fix the short-dated left tail) and volatility surface construction (to give the smile an arbitrage-free fitting infrastructure).

Reading list. Hull, Options, Futures, and Other Derivatives, 10th ed. (2018). Wilmott, Paul Wilmott on Quantitative Finance, 2nd ed. (2007), on stochastic volatility in general. Heston, “A Closed-Form Solution for Options with Stochastic Volatility,” Review of Financial Studies (1993) — the original is very readable. Gatheral, The Volatility Surface (2006), for the smile dynamics. Andersen, “Simple and Efficient Simulation of the Heston Model” (2008), for variance simulation schemes. Carr & Madan (1999) and Fang & Oosterlee (2008) for transform methods; Lord & Kahl (2010) for the complex-logarithm branch treatment.