Calibration of the
Heston Model.

Beyond constant volatility.

The Black-Scholes framework assumes a constant volatility, and therefore predicts a flat implied volatility surface across all strikes and maturities. Since 1987, equity options show the opposite: implied volatility decreases with the strike (the skew), with a term structure that depends on maturity. Stochastic volatility models let the variance of the underlying follow its own diffusion to capture these effects.

The Heston (1993) model remains one of the most widely used. Its main advantage is a semi-closed-form price for European options through Fourier transforms, which makes calibration to market quotes computationally efficient.

A calibration on real data also exposes a second, often underestimated, layer of complexity: the market does not quote the risk-neutral forward of the underlying. Rate and, above all, dividend assumptions must be inferred from the option prices themselves before any stochastic volatility parameter can be fitted. In practice, this issue dominated the calibration error on long-dated options far more than any deficiency of the Heston dynamics.

The study derives the characteristic function from the Riccati system, together with the Gil–Pelaez, Carr–Madan and COS pricing formulas and the Greeks. It shows that a price-space loss is, to first order, a vega-weighted implied volatility loss, and that the short-maturity smile identifies the product \(\rho\xi\) rather than \(\rho\) and \(\xi\) separately.

The Heston model equations.

Under the risk-neutral measure \(\mathbb{Q}\), the asset price \(S_t\) and its variance \(v_t\) follow the coupled stochastic differential equations:

Equivalently, \(W^v=\rho W^S+\sqrt{1-\rho^2}\,W^{\perp}\) with \(W^{\perp}\) independent of \(W^S\). This Cholesky form is the one used in simulation (section 10). Here \(r\) is the zero-coupon rate and \(q\) the cost of carry, which includes both the dividend yield and the repo/borrow spread. Both are term structures in practice, written \(r(t)\) and \(q(t)\). Applying Itô's formula to \(x_t=\ln S_t\) gives the log-price dynamics used to derive the characteristic function:

By Feynman–Kac, the value \(V(t,S,v)\) of a European derivative solves a two-dimensional backward PDE, with terminal condition \(V(T,S,v)=(S-K)^+\) for a call:

The mixed-derivative term \(\rho\xi vS\,\partial_S\partial_v V\) is what lets the model generate a skew. Unlike Black-Scholes, no closed form in elementary functions exists: the solution is obtained by Fourier methods (section 05).

Five parameters, five effects on the surface.

The model is governed by \(\Theta=(v_0,\kappa,\theta,\xi,\rho)\). The goal of the calibration is the vector \(\Theta\) that minimizes the error between market and model at every point of the volatility surface. Each parameter has a distinct effect:

For small maturities, the ATM skew of Heston has a finite limit, with \(k=\ln(K/F)\):

Two consequences matter for calibration. The skew depends on the product \(\rho\xi\), so the first-order slope cannot separate \(\rho\) from \(\xi\): only the curvature can. And for equity indices, where \(\rho\lt0\), the skew is negative, as observed.

The CIR variance process and the Feller condition.

The variance follows a Cox–Ingersoll–Ross process. Conditionally on \(v_s\), with \(\tau=t-s\), \(v_t\) is a scaled non-central chi-squared variable. This exact law is the basis of exact simulation and of the QE scheme (section 10):

As \(t\to\infty\), \(v_t\) converges to a Gamma law with shape \(\alpha=2\kappa\theta/\xi^2\) and rate \(\beta=2\kappa/\xi^2\). Its density behaves like \(v^{\alpha-1}\) near zero and blows up when \(\alpha\lt1\). This is the intuition behind the Feller condition, which keeps the variance strictly positive:

The Feller ratio \(2\kappa\theta/\xi^2\) is exactly the shape parameter of the stationary Gamma law. Market calibrations often violate it. This is not an issue for Fourier pricing, since the characteristic function remains valid, but Monte Carlo schemes then need an appropriate boundary treatment.

Integrating the mean gives the expected integrated variance, which drives the ATM term structure:

Dividing by \(T\) gives the approximate ATM implied variance: it moves from \(v_0\) at short maturity to \(\theta\) at long maturity, on a time scale \(1/\kappa\).

Characteristic function and Fourier inversion.

The European call price is written with two exercise probabilities: \(P_2=\mathbb{Q}(S_T>K)\), and \(P_1\), the same probability under the share measure. Both are obtained by Fourier inversion of the characteristic function.

Since the coefficients of the log-price and variance equations are affine in \((x,v)\), the characteristic function has an exponential-affine form \(\exp\big(iux+A(\tau)+\Delta(\tau)v\big)\). Substituting into the backward Kolmogorov equation and collecting the terms in \(v^0\) and \(v^1\) gives a Riccati system with constant coefficients:

Solving it with \(A(0)=\Delta(0)=0\) gives the characteristic function of \(\ln S_T\), written in Gatheral's numerically stable form:

The two probabilities follow from the Gil–Pelaez inversion theorem. Under the share measure, the characteristic function becomes \(\phi(u-i)/\phi(-i)\), with \(\phi(-i)=\mathbb{E}^{\mathbb{Q}}[S_T]=F(T)\):

Both integrands stay bounded near \(u=0\). Puts follow from put-call parity. The integrals are computed by adaptive quadrature (scipy.integrate.quad), truncated at \(u_{\max}=100\) for every maturity: simple to implement correctly and fast enough for a few hundred quotes.

Two faster alternatives were studied. Carr–Madan prices a whole grid of \(N\) strikes with a single FFT in \(O(N\log N)\), by damping the call price with a factor \(e^{\alpha k}\):

The COS method (Fang and Oosterlee) expands the density of \(\ln(S_T/K)\) in a cosine series on a truncated interval sized from the cumulants. The series converges exponentially for smooth densities, so a few dozen terms often suffice, which makes it the natural candidate for the inner loop of a calibration.

The call is homogeneous of degree one in \((S_0,K)\), so Euler's relation gives the Greeks in closed form, the Gamma coming from the risk-neutral density:

Parameter sensitivities are obtained by differentiating under the integral sign. For instance \(\partial_{v_0}\phi=D(u,T)\,\phi\), which gives analytical gradients for the calibration.

The forward isn't quoted — it has to be rebuilt.

Index options are quoted against a forward \(F(t)=S_0e^{(r(t)-q(t))t}\) that is not directly observable. The rate \(r(t)\) is bootstrapped from EURIBOR 6M zero-coupon rates with a PCHIP interpolation. The cost of carry \(q(t)\), which combines dividends and repo, is neither flat nor quoted, and for the CAC40 it is strongly seasonal, concentrated between April and June.

The effect of a wrong forward can be quantified. Let \(\hat q\) be the carry used instead of the true \(q\), and \(\hat F\) the resulting forward. Holding the volatility dynamics fixed:

If \(\hat q=0\lt q\), the forward is too high: every call is overpriced and every put underpriced. At the longest maturity, \(T=4.84\), omitting a constant \(q=3\%\) overstates the forward by \(e^{0.1452}-1\approx15.6\%\). The bias is not only vertical: in log-moneyness, the whole smile is shifted horizontally.

Here the shift is about 0.145, a third of the total standard deviation \(\sigma\sqrt T\approx0.44\) for an illustrative 20% volatility. No Heston parameter set can reproduce a horizontal translation, so the optimizer tries to mimic it with extreme values of \(\rho\) and \(\xi\).

For each expiry, put-call parity implies that the call and put price curves cross exactly at the forward. Call and put quotes are interpolated in strike with a shape-preserving PCHIP interpolator (cubic splines would create spurious oscillations), and the crossing is found by bisection on \(K\mapsto C-P=D(T)(F(T)-K)\), which is strictly decreasing:

Parity is even linear in \(K\): \(C_i-P_i=\alpha+\beta K_i\). A least-squares fit on near-the-money strikes averages out quote noise and also returns an implied discount factor, a consistency check of the rate curve: \(\hat F(T)=-\hat\alpha/\hat\beta\) and \(\hat D(T)=-\hat\beta\approx e^{-r(T)T}\).

Assuming \(q(t)\) piecewise constant between listed expiries \(T_1\lt\cdots\lt T_n\), each segment is solved recursively:

The curve is extrapolated flat before \(T_1\) and after \(T_n\), and the forward at any maturity is \(F(T)=S_0\exp\big(r(T)T-\int_0^T q(t)\,dt\big)\). The bootstrap is a discrete derivative, so it amplifies noise: an error \(\delta_k\) on \(\ln F(T_k)\) becomes

With two expiries 4 days apart, a 1 basis-point error on \(\ln F\) turns into about 0.9% on \(q_k\). The bootstrapped segments at short maturities must therefore be checked, or the cumulative carry smoothed, before use.

The machinery is split into two single-responsibility classes, usable by any pricing engine:

The pipeline: (1) solve parity for each expiry to get \(F(T_k)\); (2) bootstrap \(q_k\); (3) evaluate \(F(T)\), \(q(T)\) and \(D(T)\) at each quote's maturity; (4) invert prices into implied volatilities with these inputs; (5) minimize the loss.

Why calibrate on volatility, not price.

The first implementation minimized the mean squared error on option prices:

It did not give good results. A 5-year CAC40 option costs about €1,000, a 9-day option a dozen euros: minimizing the price error makes the calibration concentrate on long maturities, which weigh far more. The loss retained is vega-weighted on implied volatilities, so that every maturity carries the same weight:

It has two limits: options with a small vega (deep out of the money) barely count, and the vega used comes from market data, not from the model.

The initial point is \(v_0=0.04\), \(\theta=0.04\), \(\kappa=2\), \(\xi=0.3\), \(\rho=-0.7\): a 20% volatility, usual for major indices, and a return to equilibrium in about 4 months (\(\ln2/\kappa=0.35\) years). L-BFGS-B, however, is a local optimizer, and several runs showed that it does not reliably reach the global minimum.

Starting points \(\Theta\) are therefore drawn in a Latin hypercube within the bounds below, each start is refined by L-BFGS-B, and a Tikhonov penalty is added:

Seven starts out of ten fall into the same minimum, with the same RMSE of 0.4211 vol points. This point is used to rebuild the implied volatility surface. The other runs stop in worse local minima, with \(\kappa\) up to 14.9 and \(\xi\) up to 4.0: a sign of the ill-posedness discussed in section 11.

Interpolating the surface without arbitrage.

The market implied volatility surface is obtained by inverting Black-Scholes at every quoted point, then interpolated linearly across strikes and in total variance (V2T) across expiries. A quoted surface is free of static arbitrage if and only if:

The second condition says that the risk-neutral density is non-negative (Breeden–Litzenberger). The calendar condition is written in forward-normalized form, one more reason to get \(F(t)\) right. Heston prices are arbitrage-free by construction, but noisy quotes are not, and must be filtered before calibration.

Impact of correcting the forward.

Price RMSE on the 142-quote CAC40 set across the successive stages of debugging:

The drop from €135 to €32 with a crude constant dividend confirms that the forward, not the stochastic volatility dynamics, was the dominant source of error at long maturity. The fully bootstrapped, seasonal \(q(T)\) is expected to close most of the remaining gap and to bring \(\rho\) and \(\xi\) back from their bounds toward economically plausible values.

On the smiles, the model reproduces the market well from the second expiry onward. The shortest expiry (\(T=0.02\) years) is completely mispriced: the structural limit of a pure diffusion (section 11). As a sanity check, semi-analytical prices are compared with Monte Carlo prices from HestonModel.generate_paths.

The exercise suggests a checklist linking what the calibration shows to its likely cause, in the order in which to investigate:

Simulating the variance process.

Calibration relies on Fourier pricing, but path-dependent products need Monte Carlo under the calibrated parameters. A naive Euler scheme can produce negative variances when the Feller condition is violated. Full truncation replaces \(v\) with \(v^+=\max(v,0)\) in the drift and diffusion, with time step \(\Delta\):

with \(Z_v=\rho Z_S+\sqrt{1-\rho^2}\,Z_{\perp}\). The scheme is simple and robust, but biased for large \(\Delta\) or when \(2\kappa\theta/\xi^2\ll1\).

Andersen's Quadratic-Exponential (QE) scheme uses the exact transition law by moment matching. Given \(v_t\), it computes the exact conditional mean and variance, then \(\psi=s^2/m^2\):

Below a switching level \(\psi_c\) (typically 1.5), the next variance is a squared shifted Gaussian, \(a(b+Z)^2\). Above it, it is drawn from a mixture of a point mass at zero and an exponential law. The scheme matches the first two conditional moments exactly and behaves much better than Euler when Feller is violated, which is precisely the regime met in practice.

What the model still gets wrong.

From dynamics to robust calibration.

Paper and scripts.