Calibrating Heston to an Option Chain With Weighted Pricing Errors

Price a grid of liquid European options with one shared Heston parameter vector, then minimise scaled model-minus-market residuals with a bounded nonlinear least-squares solver. Use multiple strikes and maturities, inspect residuals by quote, and run the fit from several starts. A small average error can coexist with weakly identified parameters or a poor fit at particular maturities.

Decide what a large pricing error means before optimising

For parameters $\beta=(v_0,\kappa,\theta,\xi,\rho)$ , use a characteristic-function pricer to evaluate $C_{\mathrm{H}}(K_i,T_i;\beta)$ . Let $M_i$ be the quote midpoint and $s_i$ a positive scaling factor, such as the bid–ask half-spread with a floor. Fit

\widehat\beta=\arg\min_\beta\sum_i \left[\frac{C_{\mathrm{H}}(K_i,T_i;\beta)-M_i}{s_i}\right]^2.

Optimise over positive variances, mean reversion, and volatility of variance, with

|\rho|<1

. The example evaluates the Heston call as

SP_1-Ke^{-rT}P_2

using two characteristic-function integrals.

P_2

weights the exercise event under the pricing measure.

P_1

weights that same event using the terminal stock as numeraire. The integral for

P_1

evaluates the characteristic function at a shifted complex argument and divides by the model's expected terminal stock return factor. Check whether your pricing or simulation convention requires an additional constraint such as the Feller condition. A characteristic-function fit can use parameters for which the continuous variance may touch zero. An imposed Feller restriction can change the best fit and should be reported as part of the model specification.

Separate a good fit from stable parameters

Use spot $100$ , a flat rate of $3\%$ , five strikes from $85$ to $115$ , and maturities of three months, six months, and one year. Generate a synthetic surface with $(0.045,1.60,0.040,0.35,-0.65)$ . With a constant residual scale of $0.10$ and a different starting point, the code recovers those inputs to the displayed precision. The synthetic quotes and fit use the same quadrature implementation, so this is a plumbing check. A market-data fit should also report dollar RMSE, quote errors relative to spreads, fit by maturity, and whether the result is sensitive to the start or bound choices. Run a separate pricer or vary quadrature resolution to check the fitted prices.

Inspect residuals by strike and expiry before ranking fits by a single objective value. A small average error can conceal a systematic miss on short-dated out-of-the-money options. Run the optimiser from several admissible starting points. If very different parameters produce nearly the same residual pattern, the chosen chain may identify its prices much better than it identifies each parameter. To assess whether the fit is useful, hold out strikes or expiries and check how the model reprices them. Treat that as a separate test from reproducing synthetic input prices.

The residual scale determines what the optimiser is being asked to fit. Dividing a price error by the bid–ask width gives a narrow, liquid quote more weight than a wide quote with the same dollar error. Very narrow or unreliable quotes need a floor on that scale, or one point can dominate the result. If you instead divide by vega to approximate an implied-volatility error, be careful near strikes where vega is small. Inspect both the weighted residuals and the unweighted price errors to see what the fit actually improved.

A synthetic chain and a calibration loop

import numpy as np
from scipy.optimize import least_squares
from numpy.polynomial.legendre import leggauss

spot, rate = 100.0, 0.03
maturity = np.array([0.25, 0.5, 1.0])[:, None]
strike = np.array([85, 95, 100, 105, 115])[None, :]
nodes, weights = leggauss(100)
u = 90 * (nodes + 1)
weights = 90 * weights

def log_return_cf(z, t, parameters):
    v0, kappa, theta, xi, rho = parameters
    a = kappa - rho * xi * 1j * z
    d = np.sqrt(a**2 + xi**2 * (z**2 + 1j * z))
    g = (a - d) / (a + d)
    decay = np.exp(-d * t)
    c = 1j * z * rate * t + (kappa * theta / xi**2) * (
        (a - d) * t - 2 * np.log((1 - g * decay) / (1 - g))
    )
    b = (a - d) * (1 - decay) / (xi**2 * (1 - g * decay))
    return np.exp(c + b * v0)

def heston_calls(parameters):
    t = maturity[..., None]
    strike_phase = np.exp(-1j * np.log(strike / spot)[..., None] * u)
    integrand = strike_phase / (1j * u)
    q1 = 0.5 + np.sum(weights * (integrand * log_return_cf(
        u - 1j, t, parameters
    ) / np.exp(rate * t)).real, axis=-1) / np.pi
    q2 = 0.5 + np.sum(weights * (integrand * log_return_cf(
        u, t, parameters
    )).real, axis=-1) / np.pi
    return spot * q1 - strike * np.exp(-rate * maturity) * q2

synthetic = heston_calls([0.045, 1.60, 0.040, 0.35, -0.65])
quote_scale = 0.10

def residuals(parameters):
    return ((heston_calls(parameters) - synthetic) / quote_scale).ravel()

fit = least_squares(
    residuals, x0=[0.03, 1.10, 0.05, 0.25, -0.40],
    bounds=([0.005, 0.3, 0.005, 0.08, -0.95], [0.15, 5.0, 0.15, 0.8, 0.30]),
    x_scale="jac",
)
print("parameters", fit.x.round(5))
print("RMSE", np.sqrt(np.mean((heston_calls(fit.x) - synthetic)**2)))

The volatility surface tool makes it easier to see which features of a chain a single calibration error might conceal.

Decide what a large pricing error means before optimising

Separate a good fit from stable parameters

A synthetic chain and a calibration loop

References and next steps