Michael Brenndoerfer | Data & AI, Private Equity, Technology

Michael Brenndoerfer

Calibration and Parameter Estimation in Models › Calibrating the Heston Stochastic Volatility Model

Michael BrenndoerferDecember 13, 202512 min read

Quantitative FinanceAll chapters

Part 8 of 15Econometrics and Statistical Modeling

Chapter 45Calibration and Parameter Estimation in Models › Calibrating the Heston Stochastic Volatility Model

Reading section 9 of 17 · Calibrating the Heston Stochastic Volatility Model · 12 minSection 9 of 17 · 12 min

Section 9 of 17 · 12 min read

Calibration and Parameter Estimation in Models › Calibrating the Heston Stochastic Volatility ModelOne section · Sections in this chapter9 / 17 ↓

Sections in this chapter17

Reading tools Definitions & inline links

Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.

Article links

Make inline references clickable

Heston is a standard stochastic-volatility model for equity derivatives, with joint dynamics for price and variance. In the no-dividend pricing example below, use the risk-neutral measure $Q$ , set $\mu=r$ , and interpret the variance coefficients as pricing-measure parameters. Historical coefficients require a separate measure-change model.

\begin{aligned} dS_t &= \mu S_t dt + \sqrt{v_t} S_t dW_t^S \\ dv_t &= \kappa(\theta - v_t)dt + \xi \sqrt{v_t} dW_t^v \end{aligned}

where:

$S_t$ : asset price at time $t$
$v_t$ : instantaneous variance at time $t$ (note: variance, not volatility)
$\mu$ : drift rate of the asset price
$dW_t^S, dW_t^v$ : Brownian motions with correlation $\rho$
$v_0$ : initial variance, the starting level for the variance process
$\kappa$ : speed of mean reversion for variance, determining how fast variance returns to its long-run level
$\theta$ : long-term variance level, the equilibrium to which variance is pulled
$\xi$ : volatility of variance (vol of vol), controlling how much the variance itself fluctuates
$\rho$ : correlation between spot price and variance, typically negative for equities

The variance process $v_t$ follows a Cox-Ingersoll-Ross (CIR) process and exhibits mean reversion to the long-run level $\theta$ . With positive parameters and positive initial variance, the strict Feller inequality $2\kappa\theta>\xi^2$ is sufficient for the exact continuous-time variance to stay strictly positive. This is a sufficient restriction chosen for the calibration below, not a requirement for every valid Heston specification. The square root in the variance diffusion term ensures that variance fluctuations scale appropriately with the variance level: when variance is high, it can fluctuate more; when variance is low, fluctuations are dampened.

Negative spot-variance correlation tends to associate downward price innovations with increases in variance. It contributes to the downside implied-volatility skew often observed in equity options: lower-strike options have higher implied volatility than higher-strike options at the same maturity. This does not universally rank raw premiums of puts and calls at different strikes.

The Heston model has a semi-closed-form solution for European option prices via characteristic functions. The call price is:

C = S_t P_1 - K e^{-rT} P_2

where:

$C$ : call option price
$S_t$ : asset price at time $t$
$K$ : strike price
$r$ : risk-free interest rate
$T$ : time to expiration
$P_1, P_2$ : probabilities derived from the characteristic function, representing the delta and the risk-neutral exercise probability respectively

The terms $P_1$ and $P_2$ are obtained by characteristic-function quadrature. For European options, this can cost substantially less than Monte Carlo at a comparable accuracy target, but performance depends on parameters and implementation. Numerical integration still requires convergence and price-accuracy checks.

Let's implement calibration to option prices. The numerical helper below compares successive integration cutoffs and lower endpoints and retains quadrature error checks. These are convergence diagnostics, not a rigorous universal bound on truncation error. Its supported parameter range and small-price accuracy must be checked independently before interpreting implied-volatility features:

In[21]:

Code

from scipy.integrate import quad


def heston_char_func(u, S0, K, T, r, v0, kappa, theta, xi, rho):
    """
    Heston characteristic function for pricing.
    """
    i = complex(0, 1)
    a = kappa * theta
    b = kappa

    d = np.sqrt((rho * xi * i * u - b) ** 2 + xi**2 * (i * u + u**2))
    g = (b - rho * xi * i * u - d) / (b - rho * xi * i * u + d)

    C = r * i * u * T + (a / xi**2) * (
        (b - rho * xi * i * u - d) * T - 2 * np.log((1 - g * np.exp(-d * T)) / (1 - g))
    )
    D = (
        (b - rho * xi * i * u - d)
        / xi**2
        * (1 - np.exp(-d * T))
        / (1 - g * np.exp(-d * T))
    )

    return np.exp(C + D * v0 + i * u * np.log(S0))


def heston_call_price(
    S0, K, T, r, v0, kappa, theta, xi, rho, price_tol=1e-7, max_cutoff=1600
):
    """Approximate no-dividend European-call price with convergence checks."""
    inputs = (S0, K, T, r, v0, kappa, theta, xi, rho, price_tol)
    if not all(np.isfinite(x) for x in inputs):
        raise ValueError("Heston inputs must be finite")
    if min(S0, K, T, v0, kappa, theta, xi, price_tol) <= 0 or abs(rho) >= 1:
        raise ValueError("Invalid Heston parameter domain")
    discounted_strike = K * np.exp(-r * T)
    intrinsic = max(S0 - discounted_strike, 0.0)

    def integrand1(u):
        cf = heston_char_func(u - 1j, S0, K, T, r, v0, kappa, theta, xi, rho)
        return np.real(np.exp(-1j * u * np.log(K)) * cf / (1j * u * S0 * np.exp(r * T)))

    def integrand2(u):
        cf = heston_char_func(u, S0, K, T, r, v0, kappa, theta, xi, rho)
        return np.real(np.exp(-1j * u * np.log(K)) * cf / (1j * u))

    def integrate(lower, cutoff):
        scale = max(S0 + abs(discounted_strike), 1.0)
        epsabs = price_tol * np.pi / (32 * scale)
        values = [
            quad(f, lower, cutoff, epsabs=epsabs, epsrel=1e-10, limit=500)
            for f in (integrand1, integrand2)
        ]
        p1, p2 = (0.5 + item[0] / np.pi for item in values)
        price = S0 * p1 - discounted_strike * p2
        error = (S0 * values[0][1] + abs(discounted_strike) * values[1][1]) / np.pi
        if not np.isfinite(price) or not np.isfinite(error):
            raise ValueError("Nonfinite Heston quadrature")
        return price, error

    # Check both the omitted near-zero interval and increasing tail cutoffs.
    lower = 1e-7
    previous, _ = integrate(lower, 100)
    cutoff = 200
    while cutoff <= max_cutoff:
        coarse, _ = integrate(lower, cutoff)
        current, error = integrate(lower / 100, cutoff)
        if (
            abs(current - previous) <= price_tol
            and abs(current - coarse) <= price_tol / 4
            and error <= price_tol / 4
        ):
            if current < intrinsic or current >= S0:
                raise ValueError("Numerical Heston price violates call bounds")
            return float(current)
        previous = current
        lower /= 100
        cutoff *= 2
    raise RuntimeError("Heston quadrature did not converge to the requested tolerance")

We generate synthetic option quotes with 1% multiplicative noise and then calibrate. The generating parameters have $2\kappa\theta=0.20$ and $\xi^2=0.25$ , so the strict Feller restriction imposed below excludes the generating parameter vector. This example is a constrained fit to noisy synthetic inputs, not a test of recovering the generating parameters inside the feasible set:

In[22]:

Code

# Market setup
S0 = 100
r = 0.03

# True Heston parameters
v0_true = 0.04  # 20% initial vol
kappa_true = 2.0
theta_true = 0.05  # 22.4% long-term vol
xi_true = 0.5  # vol of vol
rho_true_h = -0.7  # Heston correlation, separate from the SABR example

# Generate market option prices for different strikes and maturities
option_data = []
maturities_opt = [0.25, 0.5, 1.0]
strikes_by_moneyness = [0.85, 0.90, 0.95, 1.00, 1.05, 1.10, 1.15]

np.random.seed(123)
for T in maturities_opt:
    for m in strikes_by_moneyness:
        K = S0 * m
        true_price = heston_call_price(
            S0, K, T, r, v0_true, kappa_true, theta_true, xi_true, rho_true_h
        )
        # Add small noise
        noisy_price = true_price * (1 + np.random.normal(0, 0.01))
        option_data.append({"T": T, "K": K, "price": noisy_price})

In[23]:

Code

from scipy.optimize import differential_evolution


def heston_calibration_objective(params, S0, r, option_data):
    """Relative squared price errors under the chosen strict Feller restriction."""
    if not option_data or not np.isfinite(S0) or S0 <= 0 or not np.isfinite(r):
        raise ValueError("Require finite setup and nonempty quotes")
    for opt in option_data:
        if (
            not all(np.isfinite(opt[k]) for k in ("K", "T", "price"))
            or min(opt["K"], opt["T"], opt["price"]) <= 0
        ):
            raise ValueError("Require finite positive strikes, maturities and prices")
    params = np.asarray(params, dtype=float)
    if params.shape != (5,) or not np.all(np.isfinite(params)):
        return np.inf
    v0, kappa, theta, xi, rho = params
    if min(v0, kappa, theta, xi) <= 0 or abs(rho) >= 1:
        return np.inf
    if 2 * kappa * theta <= xi**2:
        return np.inf
    try:
        model_prices = np.array(
            [
                heston_call_price(S0, opt["K"], opt["T"], r, v0, kappa, theta, xi, rho)
                for opt in option_data
            ]
        )
    except (ValueError, RuntimeError, FloatingPointError, OverflowError):
        return np.inf
    targets = np.array([opt["price"] for opt in option_data])
    if not np.all(np.isfinite(model_prices)):
        return np.inf
    return float(np.sum(((model_prices - targets) / targets) ** 2))


# Use differential evolution for global optimization
bounds_heston = [
    (0.001, 0.5),  # v0
    (0.1, 5.0),  # kappa
    (0.001, 0.5),  # theta
    (0.1, 2.0),  # xi
    (-0.99, 0.99),  # rho
]

result_heston = differential_evolution(
    heston_calibration_objective,
    bounds_heston,
    args=(S0, r, option_data),
    maxiter=200,
    seed=42,
    polish=True,  # Finish with L-BFGS-B
)

require_successful_fit(result_heston, "Heston calibration")
v0_cal, kappa_cal_h, theta_cal_h, xi_cal, rho_cal_h = result_heston.x
if 2 * kappa_cal_h * theta_cal_h <= xi_cal**2:
    raise RuntimeError("Heston fit failed the full-precision strict Feller check")

Out[24]:

Console

Heston Model Calibration Results:
--------------------------------------------------
Parameter   |   True    | Calibrated |  Error
--------------------------------------------------
v0          |  0.0400   |   0.0378   | 0.0022
kappa       |  2.0000   |   1.4323   | 0.5677
theta       |  0.0500   |   0.0527   | 0.0027
xi          |  0.5000   |   0.3884   | 0.1116
rho         | -0.7000   |  -0.7300   | 0.0300

Full-precision Feller check: 2κθ = 0.1508535690, ξ² = 0.1508321712, margin = 2.140e-05: ✓

The Heston calibration gives a constrained fit to the noisy synthetic quotes. The fitted parameters differ from the generating values in a calibration that includes both added quote noise and the strict Feller restriction. The generating vector is excluded; this example does not isolate how those factors affect the individual parameter estimates. The printed check uses full-precision calibrated parameters to test $2\kappa\theta>\xi^2$ . The rounded displayed values alone do not establish the margin of that inequality. With the positive parameters and initial variance used here, this restriction is sufficient for strictly positive exact continuous-time variance.

To visualize the model's behavior, we compute the implied volatility surface corresponding to the calibrated parameters. This involves generating Heston prices across a grid of strikes and maturities, then inverting the Black-Scholes formula to find the implied volatility for each.

In[25]:

Code

from scipy.optimize import brentq
from scipy.stats import norm


def black_scholes_call(S, K, T, r, sigma):
    if sigma == 0:
        return max(S - K * np.exp(-r * T), 0.0)
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)


def implied_vol_from_price(price, S, K, T, r, tol=1e-8):
    """Invert an unmodified valid call price with an expanding bracket."""
    if not all(np.isfinite(x) for x in (price, S, K, T, r, tol)):
        raise ValueError("Option inputs must be finite")
    if min(S, K, T, tol) <= 0:
        raise ValueError("Require S, K, T and volatility tolerance > 0")
    intrinsic = max(S - K * np.exp(-r * T), 0.0)
    if price < intrinsic or price >= S:
        raise ValueError("Target is outside the finite-IV call price domain")
    if price == intrinsic:
        return 0.0

    def price_error(sigma):
        return black_scholes_call(S, K, T, r, sigma) - price

    upper_vol = 1.0
    for _ in range(20):
        upper_error = price_error(upper_vol)
        if not np.isfinite(upper_error):
            raise ValueError("Nonfinite Black-Scholes bracket evaluation")
        if upper_error >= 0:
            return brentq(price_error, 0.0, upper_vol, xtol=tol)
        upper_vol *= 2
    raise RuntimeError(
        "Valid target was not bracketed within the numerical search limit"
    )


# Create volatility surface
# Selected display domain; each price must pass the quadrature checks.
# A positive price or a successful IV inversion alone does not establish accuracy.
strikes_surface = np.linspace(80, 118, 20)
maturities_surface = np.linspace(0.1, 1.5, 15)

vol_surface = np.zeros((len(maturities_surface), len(strikes_surface)))

for i, T in enumerate(maturities_surface):
    for j, K in enumerate(strikes_surface):
        price = heston_call_price(
            S0, K, T, r, v0_cal, kappa_cal_h, theta_cal_h, xi_cal, rho_cal_h
        )
        vol_surface[i, j] = implied_vol_from_price(price, S0, K, T, r)

Out[27]:

Visualization

3D surface plot of implied volatility varying across strikes and maturities from the calibrated Heston model. — Heston implied volatilities for the fitted parameters on the selected grid, using numerically checked prices. The shape describes this example; numerical accuracy must be sufficient to interpret small tail features.

The surface and cross-sections describe the joint effect of this fitted parameter vector over the selected grid. Negative correlation contributes to downside skew, while mean reversion and vol of vol affect the maturity profile. One fitted surface does not isolate those effects, and small numerical tail features require separate accuracy checks.

Out[29]:

Visualization

Line chart showing Heston implied volatility smiles at four maturities from 0.1 to 1.5 years, with the smile becoming progressively flatter at longer maturities. — Selected maturity cross-sections of the fitted Heston example. Compare their strike ranges and numerical accuracy before interpreting skew and term-structure differences.