Advanced models: Heston, SABR and local volatility.

Why other models are needed

Black-Scholes assumes a constant volatility. Lesson 3.4.5 showed that the market does not quote it that way: there is a smile and a term structure. A parametric surface such as SVI describes those prices, but it says nothing about how they will move tomorrow, nor does it serve to price an option whose payoff depends on the path (a barrier, a forward-start option) or to compute consistent hedges. That requires a model of the underlying that reproduces today's surface. There are two families:

flowchart TB
  B[Black-Scholes] --> L[Local volatility]
  B --> E[Stochastic vol.]
  L --> D[Dupire 1994]
  E --> H[Heston 1993]
  E --> S[SABR 2002]
  D --> SLV[Stochastic local]
  H --> SLV

This lesson is about “concept and code”: each model is programmed, checked with QuantLib and calibrated on the ORX surface of June 4, 2025 from the course's synthetic data, whose truth is a known SSVI. That makes it possible to separate model error from quote noise.

The starting surface

As in lesson 3.4.5: the implied volatility of each out-of-the-money option with a positive bid, with the forward and the rate of each slice taken from the truth.

import math
import numpy as np
import pandas as pd
from scipy.special import ndtr
from scipy.optimize import least_squares

sup = pd.read_parquet("opciones3_superficie.parquet")
ver = pd.read_parquet("opciones3_verdad.parquet")


def bs_norm(k, T, s):
    """Undiscounted Black-76 call divided by F."""
    sd = s * np.sqrt(T)
    d1 = -k / sd + sd / 2
    return ndtr(d1) - np.exp(k) * ndtr(d1 - sd)


def vega_norm(k, T, s):
    sd = s * np.sqrt(T)
    d1 = -k / sd + sd / 2
    return np.exp(-d1 * d1 / 2) * np.sqrt(T / (2 * math.pi))


def iv_norm(c, k, T):
    """Safeguarded Newton (lesson 3.4.4), normalized call."""
    lo, hi = np.full(len(c), 1e-4), np.full(len(c), 5.0)
    s = np.full(len(c), 0.3)
    for _ in range(60):
        d = bs_norm(k, T, s) - c
        lo, hi = np.where(d < 0, s, lo), np.where(d > 0, s, hi)
        paso = s - d / np.maximum(vega_norm(k, T, s), 1e-12)
        mal = (paso <= lo) | (paso >= hi) | ~np.isfinite(paso)
        s = np.where(mal, (lo + hi) / 2, paso)
    return s


hoy = pd.Timestamp("2025-06-04")
x = sup[sup["fecha"] == hoy].merge(ver, on=["fecha", "vencimiento"])
x["T"] = x["dias"] / 365
x["k"] = np.log(x["strike"] / x["forward"])
x = x[np.where(x["tipo"] == "C", x["k"] >= 0, x["k"] < 0)
      & (x["bid"] > 0)].copy()
D = np.exp(-x["tipo_libre"] * x["T"])
medio = (x["bid"] + x["ask"]) / 2 / (D * x["forward"])
# put -> call by parity (in units of F, undiscounted)
x["c"] = np.where(x["tipo"] == "C", medio, medio + 1 - np.exp(x["k"]))
x["iv"] = iv_norm(x["c"].to_numpy(), x["k"].to_numpy(),
                  x["T"].to_numpy())
print(len(x), "options in", x["dias"].nunique(), "expiries")

341 options in 12 expiries

These are the same 341 options of the twelve expiries of lesson 3.4.5. To work with a single formula, each price is converted to an undiscounted call in units of the forward (the puts, by the parity of lesson 3.4.1): this way Heston, SABR and Black-76 are compared on the same magnitude.

Heston: the variance is random too

Heston (1993) adds to Black-Scholes a second equation for the instantaneous variance $v_t$:

• $v_0$: today's variance. $\theta$: the long-term variance it reverts to. $\kappa$: the speed of reversion (half-life $\ln 2/\kappa$ years).
• $\xi$: volatility of variance. It gives curvature to the smile (more expensive wings).
• $\rho$: correlation between price and variance. Negative, it tilts the smile: it is the leverage effect, the skew.
• Feller condition, $2\kappa\theta \ge \xi^2$: if it holds, the variance never touches zero.

Its advantage is that the characteristic function of $x_T = \ln(S_T/F)$ has a closed form. In the form of Albrecher et al. (2007), which avoids branch jumps of the complex logarithm:

With it, Lewis's (2001) formula gives the call as a real integral, with $k = \ln(K/F)$ and the discount $D_T$:

from scipy.integrate import quad
import QuantLib as ql


def heston_cf(u, T, v0, kappa, theta, xi, rho):
    """E[exp(i·u·ln(S_T/F))], form of Albrecher et al. (2007)."""
    iu = 1j * u
    b = kappa - rho * xi * iu
    d = np.sqrt(b ** 2 + xi ** 2 * (iu + u * u))
    g = (b - d) / (b + d)
    e = np.exp(-d * T)
    C = kappa * theta / xi ** 2 * (
        (b - d) * T - 2 * np.log((1 - g * e) / (1 - g)))
    Dv = (b - d) / xi ** 2 * (1 - e) / (1 - g * e)
    return np.exp(C + Dv * v0)


def heston_call(F, K, T, DT, p):
    k = math.log(K / F)
    f = lambda u: (np.exp(-1j * u * k)
                   * heston_cf(u - 0.5j, T, *p)).real / (u * u + 0.25)
    integral = quad(f, 0, 200, limit=400)[0]
    return DT * (F - math.sqrt(F * K) / math.pi * integral)


p = (0.04, 2.0, 0.04, 0.5, -0.7)        # v0, κ, θ, ξ, ρ
F, r, dias = 100.0, 0.03, 183
T, DT = dias / 365, math.exp(-0.03 * dias / 365)

d0 = ql.Date(4, 6, 2025)
ql.Settings.instance().evaluationDate = d0
dc = ql.Actual365Fixed()
plana = lambda t: ql.YieldTermStructureHandle(
    ql.FlatForward(d0, t, dc, ql.Continuous))
proc = ql.HestonProcess(plana(r), plana(r),        # q = r: F = S
                        ql.QuoteHandle(ql.SimpleQuote(F)), *p)
motor = ql.AnalyticHestonEngine(ql.HestonModel(proc))
for K in (80.0, 100.0, 120.0):
    op = ql.VanillaOption(ql.PlainVanillaPayoff(ql.Option.Call, K),
                          ql.EuropeanExercise(d0 + dias))
    op.setPricingEngine(motor)
    propio = heston_call(F, K, T, DT, p)
    print(f"K = {K:.0f}: own {propio:.6f}, QuantLib "
          f"{op.NPV():.6f}, difference {propio - op.NPV():.1e}")

K = 80: own 20.458189, QuantLib 20.458189, difference -9.4e-13
K = 100: own 5.109865, QuantLib 5.109865, difference -2.0e-12
K = 120: own 0.158180, QuantLib 0.158180, difference -2.5e-12

The two implementations agree to the twelfth digit: the difference is the numerical integration error. A practical warning: QuantLib's time to expiry comes from the dates (183 calendar days with Actual365Fixed), so our own code has to use exactly the same $T$; with $T = 0.5$ the difference would be in the thousandths.

Calibrating Heston to the surface

Calibrating requires pricing hundreds of options at each iteration. quad is exact but slow; with a fixed Gauss-Legendre quadrature (160 nodes between 0 and 150) the integral is vectorized for all the options at once. The objective is the volatility error, approximated by the price error divided by the vega, which is much cheaper than inverting the volatility at each step:

nodos, pesos = np.polynomial.legendre.leggauss(160)
U, W = (nodos + 1) * 75, pesos * 75          # from [-1, 1] to [0, 150]


def heston_norm(k, T, p):
    """Undiscounted Heston call, in units of F."""
    f = (np.exp(-1j * U * k[:, None])
         * heston_cf(U - 0.5j, T[:, None], *p)).real / (U * U + 0.25)
    return 1 - np.exp(k / 2) / math.pi * (f * W).sum(axis=1)


k, T, iv = (x[c].to_numpy() for c in ("k", "T", "iv"))
objetivo, vg = bs_norm(k, T, iv), vega_norm(k, T, iv)
res = least_squares(
    lambda p: (heston_norm(k, T, p) - objetivo) / vg,
    [0.05, 3.0, 0.05, 0.8, -0.7],
    bounds=([1e-4, 0.1, 1e-4, 0.05, -0.99], [1, 30, 1, 5, 0.99]))
v0, kappa, theta, xi, rho_h = ph = res.x
print(f"v0 = {v0:.4f} (vol {math.sqrt(v0):.1%}), κ = {kappa:.2f}, "
      f"θ = {theta:.4f} (vol {math.sqrt(theta):.1%})")
print(f"ξ = {xi:.3f}, ρ = {rho_h:.3f}, "
      f"Feller 2κθ − ξ² = {2 * kappa * theta - xi ** 2:.2f}")
x["iv_heston"] = iv_norm(heston_norm(k, T, ph), k, T)
err = 100 * (x["iv_heston"] - x["iv"])
print(err.groupby(x["dias"]).agg(
    lambda z: np.sqrt(np.mean(z ** 2))).round(2).to_dict())

v0 = 0.0691 (vol 26.3%), κ = 7.55, θ = 0.0669 (vol 25.9%)
ξ = 1.995, ρ = -0.750, Feller 2κθ − ξ² = -2.97
{16: 2.51, 44: 0.52, 72: 0.76, 107: 0.87, 135: 0.86, 170: 0.78, 198: 0.71, 289: 0.49, 379: 0.6, 471: 0.88, 562: 1.28, 653: 1.61}

Calibrated Heston starts with a volatility of 26.3% ($v_0$) that reverts with $\kappa = 7.55$ (a half-life of about five weeks) to a long-term level of 25.9%, with $\xi = 1.995$ and $\rho = -0.750$. The Feller condition does not hold ($-2.97$), as is usual when calibrating indices: to generate the steep skew of the short tenors, the model needs a very high volatility of variance. The error per expiry draws a U: 2.51 points at 16 days, between 0.49 and 0.88 from 44 to 471 days and 1.61 at 653. Five parameters are not enough for the 16-day skew and the almost-two-year one at the same time.

S = pd.read_csv("opciones3_indice.csv", parse_dates=["fecha"]) \
    .set_index("fecha").loc[hoy, "cierre"]
c44 = x[x["dias"] == 44].iloc[0]
r44 = c44["tipo_libre"]
q44 = r44 - math.log(c44["forward"] / S) / c44["T"]
proc = ql.HestonProcess(plana(r44), plana(q44),
                        ql.QuoteHandle(ql.SimpleQuote(S)),
                        0.05, 3.0, 0.05, 0.8, -0.7)
modelo = ql.HestonModel(proc)
motor = ql.AnalyticHestonEngine(modelo)
ayudas = []
for f in x.itertuples():
    K = S * math.exp((r44 - q44) * f.T) * math.exp(f.k)
    a = ql.HestonModelHelper(
        ql.Period(int(f.dias), ql.Days), ql.NullCalendar(), S, K,
        ql.QuoteHandle(ql.SimpleQuote(f.iv)), plana(r44), plana(q44),
        ql.BlackCalibrationHelper.ImpliedVolError)
    a.setPricingEngine(motor)
    ayudas.append(a)
modelo.calibrate(ayudas, ql.LevenbergMarquardt(),
                 ql.EndCriteria(500, 50, 1e-8, 1e-8, 1e-8))
th_q, ka_q, xi_q, rho_q, v0_q = modelo.params()
print(f"QuantLib: v0 = {v0_q:.4f}, κ = {ka_q:.2f}, θ = {th_q:.4f}, "
      f"ξ = {xi_q:.3f}, ρ = {rho_q:.3f}")
iv_q = iv_norm(heston_norm(k, T, (v0_q, ka_q, th_q, xi_q, rho_q)),
               k, T)
for nombre, ivm in [("own", x["iv_heston"]), ("QuantLib", iv_q)]:
    e = 100 * np.sqrt(np.mean((ivm - x["iv"]) ** 2))
    print(f"Root mean squared error, {nombre} calibration: {e:.3f}")

QuantLib: v0 = 0.0710, κ = 9.16, θ = 0.0661, ξ = 2.261, ρ = -0.737
Root mean squared error, own calibration: 1.126
Root mean squared error, QuantLib calibration: 1.075

QuantLib arrives at different parameters ($\kappa = 9.16$, $\xi = 2.261$) with a slightly smaller error (1.075 against 1.126 points). Both things have an explanation. QuantLib minimizes the exact volatility error, and our own calibration approximates it with price divided by vega. And the objective function is almost flat along a ridge on which $\kappa$ and $\xi$ rise together: faster reversion offsets a more volatile variance. Two very different parameter sets give almost the same prices, and that is why Heston's calibrated parameters jump from one day to the next even when the surface barely changes.

SABR: one model per expiry

SABR (Hagan, Kumar, Lesniewski and Woodward, 2002) models the forward of an expiry and its volatility:

$\beta$ sets how the volatility depends on the level (1 lognormal, 0 normal), $\alpha$ is the level, $\rho$ the tilt and $\nu$ the volatility of volatility (the curvature). It has no mean reversion and does not link one expiry with another: it is calibrated slice by slice. Its appeal is that Hagan et al. gave a closed approximation of the implied volatility; with $z = (\nu/\alpha)(FK)^{(1-\beta)/2}\ln(F/K)$ and $x(z) = \ln\big[(\sqrt{1 - 2\rho z + z^2} + z - \rho)/(1 - \rho)\big]$:

def hagan(K, F, T, a, b, rho, nu):
    K = np.asarray(K, float)
    lfk = np.log(F / K)
    fk = (F * K) ** ((1 - b) / 2)
    z = nu / a * fk * lfk
    xz = np.log((np.sqrt(1 - 2 * rho * z + z * z) + z - rho)
                / (1 - rho))
    zx = np.where(np.abs(z) < 1e-8, 1.0, z / np.where(xz == 0, 1, xz))
    den = fk * (1 + (1 - b) ** 2 / 24 * lfk ** 2
                + (1 - b) ** 4 / 1920 * lfk ** 4)
    corr = 1 + ((1 - b) ** 2 / 24 * a * a / fk ** 2
                + rho * b * nu * a / (4 * fk)
                + (2 - 3 * rho * rho) / 24 * nu * nu) * T
    return a / den * zx * corr


for K, b in [(90.0, 1.0), (105.0, 1.0), (105.0, 0.5)]:
    a = 0.2 if b == 1 else 2.0
    propio = float(hagan(K, 100.0, 0.5, a, b, -0.5, 0.6))
    # QuantLib arguments: strike, forward, T, α, β, ν, ρ
    qlv = ql.sabrVolatility(K, 100.0, 0.5, a, b, 0.6, -0.5)
    print(f"K = {K:.0f}, β = {b}: own {propio:.8f}, "
          f"QuantLib {qlv:.8f}")

K = 90, β = 1.0: own 0.21787385, QuantLib 0.21787385
K = 105, β = 1.0: own 0.19353697, QuantLib 0.19353697
K = 105, β = 0.5: own 0.19187243, QuantLib 0.19187243

Our own implementation agrees with ql.sabrVolatility in the eight digits shown, also with $\beta = 0.5$. Watch the order of the arguments in QuantLib: strike, forward, tenor, $\alpha$, $\beta$, $\nu$ and $\rho$, with $\nu$ before $\rho$, the reverse of the usual notation.

filas = []
for dias, g in x.groupby("dias"):
    F, T = g["forward"].iloc[0], g["T"].iloc[0]
    K, ivm = g["strike"].to_numpy(), g["iv"].to_numpy()
    s = least_squares(
        lambda q: hagan(K, F, T, q[0], 1.0, q[1], q[2]) - ivm,
        [0.2, -0.5, 1.0], bounds=([1e-3, -0.999, 1e-3], [3, 0.999, 20]))
    a, rho_s, nu = s.x
    x.loc[g.index, "iv_sabr"] = hagan(K, F, T, a, 1.0, rho_s, nu)
    filas.append({"dias": dias, "alfa": a, "rho": rho_s, "nu": nu,
                  "nu_raizT": nu * math.sqrt(T),
                  "rmse_sabr": 100 * np.sqrt(np.mean(s.fun ** 2)),
                  "rmse_heston": 100 * np.sqrt(np.mean(
                      (g["iv"] - x.loc[g.index, "iv_heston"]) ** 2))})
sabr = pd.DataFrame(filas).set_index("dias")
sabr.round(3)

       alfa    rho     nu  nu_raizT  rmse_sabr  rmse_heston
dias
16    0.228 -0.797  3.045     0.638      0.945        2.506
44    0.230 -0.786  1.834     0.637      1.001        0.522
72    0.231 -0.778  1.501     0.667      0.746        0.759
107   0.229 -0.760  1.309     0.709      0.499        0.875
135   0.229 -0.758  1.189     0.723      0.425        0.855
170   0.228 -0.756  1.077     0.735      0.365        0.780
198   0.228 -0.755  1.014     0.747      0.302        0.707
289   0.227 -0.753  0.856     0.762      0.231        0.489
379   0.225 -0.750  0.758     0.773      0.183        0.596
471   0.225 -0.752  0.682     0.775      0.166        0.882
562   0.223 -0.752  0.628     0.779      0.144        1.281
653   0.222 -0.754  0.583     0.780      0.128        1.609

$\alpha$ stays between 0.222 and 0.231 at all tenors: with $\beta = 1$ it is approximately the ATM volatility. $\rho$ ranges from −0.80 to −0.75 and $\nu$ falls from 3.0 to 0.58, while $\nu\sqrt{T}$ moves much less, from 0.64 to 0.78: the surface scales with the standard deviation (lesson 3.4.5). The SABR error falls with tenor, from 0.94 and 1.00 points in the first two slices to 0.13 in the last one. In the short slices, the strikes reach 10 standard deviations below the forward, and Hagan's form with $\beta = 1$ cannot follow the wing and the center at the same time: it leaves errors that change sign along the smile, of ±1.4 points in the center and up to 2.9 in the farthest wing. Heston, with a single parameter set for everything, beats SABR only at 44 days (0.52 against 1.00).

The figure compares the two errors by tenor: Heston's U, highlighted, against the steady decline of SABR. A model with per-expiry parameters fits each slice better; one with global parameters is consistent across tenors and serves to price products that depend on several. The header figure shows the detail: at 16 days Heston falls 6.9 points short in the put wing and 2 points high at the money; at 289 days, both models almost touch the points.

Local volatility: Dupire's surface

The opposite path is to add no source of randomness and let the volatility be a deterministic function of price and time, $dS_t/S_t = (r - q)\,dt + \sigma_L(S_t, t)\,dW_t$. Dupire (1994) showed that there is a unique function $\sigma_L$ that reproduces all of today's European prices, and that it is read from the surface itself. Written in total variance $w(k, T)$ and moneyness relative to the forward of each tenor (Gatheral, 2006):

The denominator is the $g(k)$ of lesson 3.4.5 and the numerator, the calendar slope: with no arbitrage, both are positive and the local volatility exists. With arbitrage, it comes out negative or infinite, and that is why it is computed on a clean surface. Here, the true SSVI of the day (its ρ, η = 0.95 and γ = 0.5, which lesson 3.4.5 recovered with the fit) with $\theta_T$ linearly interpolated between expiries:

v = ver[ver["fecha"] == hoy]
RHO, ETA, GAM = float(v["svi_rho"].iloc[0]), 0.95, 0.5
T_n = np.r_[0, v["dias"] / 365]
TH_n = np.r_[0, v["theta"]]


def w_sup(k, T):
    th = np.interp(T, T_n, TH_n)
    phi = ETA / (th ** GAM * (1 + th) ** (1 - GAM))
    return th / 2 * (1 + RHO * phi * k
                     + np.sqrt((phi * k + RHO) ** 2 + 1 - RHO ** 2))


def var_local(k, T, h=1e-4):
    w = w_sup(k, T)
    w1 = (w_sup(k + h, T) - w_sup(k - h, T)) / (2 * h)
    w2 = (w_sup(k + h, T) - 2 * w + w_sup(k - h, T)) / h ** 2
    wt = (w_sup(k, T + h) - w_sup(k, T - h)) / (2 * h)
    g = (1 - k * w1 / (2 * w)) ** 2 - w1 ** 2 / 4 * (1 / w + 0.25) \
        + w2 / 2
    return wt / g


kk = np.array([-0.3, -0.2, -0.1, 0.0, 0.1, 0.2])
T30 = 30 / 365
print(pd.DataFrame({
    "k": kk, "implícita": 100 * np.sqrt(w_sup(kk, T30) / T30),
    "local": 100 * np.sqrt(var_local(kk, T30))}).round(2)
    .to_string(index=False))
for dias in (8, 30, 180):
    T, h = dias / 365, 0.005
    p_loc = (np.sqrt(var_local(h, T))
             - np.sqrt(var_local(-h, T))) / (2 * h)
    p_imp = (np.sqrt(w_sup(h, T) / T) - np.sqrt(w_sup(-h, T) / T)) \
        / (2 * h)
    print(f"{dias} days: local ATM slope {p_loc:.2f}, implied "
          f"{p_imp:.2f}, ratio {p_loc / p_imp:.2f}")

   k  implícita  local
-0.3      47.37  65.08
-0.2      40.52  54.41
-0.1      32.34  40.91
 0.0      22.13  22.09
 0.1      16.65  16.52
 0.2      18.96  22.79
8 days: local ATM slope -3.89, implied -2.11, ratio 1.84
30 days: local ATM slope -2.04, implied -1.09, ratio 1.86
180 days: local ATM slope -0.80, implied -0.44, ratio 1.80

At the money the two agree (22.1%), but toward the puts the local rises much faster: 65.1% against 47.4% at $k = -0.3$. The at-the-money slope of the local is 1.84, 1.86 and 1.80 times the implied at 8, 30 and 180 days: it is Derman's rule, which says that the implied volatility of a strike is approximately the average of the local along the path between the current price and that strike, and that is why it has half the slope. Toward the calls, the local drops below the implied and turns up again earlier (22.8% against 19.0% at $k = 0.2$).

The figure draws the two curves at 30 days: the local one, highlighted, is the implied one “unrolled”.

rng = np.random.default_rng(347)
N, T44, pasos = 100_000, 44 / 365, 176      # 4 steps per day
dt = T44 / pasos
rejilla = np.linspace(-2.0, 1.5, 701)
X = np.zeros(N)                               # ln(S_t / F_t)
for i in range(pasos):
    s2 = np.interp(X, rejilla, var_local(rejilla, (i + 0.5) * dt))
    X += -0.5 * s2 * dt + np.sqrt(s2 * dt) * rng.standard_normal(N)
ST = np.exp(X) / np.exp(X).mean()          # exact martingale
filas = []
for k_ in (-0.2, -0.1, -0.05, 0.0, 0.05):
    c = np.maximum(ST - math.exp(k_), 0).mean()
    filas.append({"k": k_,
                  "MC": 100 * iv_norm(np.array([c]), np.array([k_]),
                                      np.array([T44]))[0],
                  "superficie": 100 * math.sqrt(w_sup(k_, T44) / T44)})
mc = pd.DataFrame(filas)
mc["dif"] = mc["MC"] - mc["superficie"]
print(f"E[S_T/F] before normalizing: {np.exp(X).mean():.5f}")
mc.round(2)

E[S_T/F] before normalizing: 0.99983
      k     MC  superficie   dif
0 -0.20  37.78       37.93 -0.15
1 -0.10  30.81       30.81  0.00
2 -0.05  26.72       26.69  0.03
3  0.00  22.29       22.22  0.07
4  0.05  18.29       18.22  0.07

Simulating the index with the local volatility of each instant and pricing the 44-day calls reproduces the smile of the surface with differences of −0.15 to 0.07 points, the Monte Carlo error with 100,000 paths and the discretization error with 176 steps. The mean of $S_T/F$ comes out at 0.99983 before normalizing: the $-\sigma^2/2$ term of the drift keeps the martingale except for sampling error. That is the property that defines Dupire's model: by construction, it returns all of today's European prices.

What today's surface does not decide

Local volatility and Heston calibrated to the same surface give the same European prices today (local exactly, Heston approximately) and yet they are very different models:

• Future skew. In local volatility the future smile flattens: the function $\sigma_L(S, t)$ is fixed today and, at long tenors, its dependence on price is smooth. In stochastic volatility models, the smile regenerates. A forward-start option or a cliquet are worth different amounts in the two.
• Delta. With local volatility, when the index falls, volatility rises along the function itself; the delta includes that effect. With Heston, the variance moves on its own, correlated. The hedge of an option (lesson 3.6.2) changes with the model.
• Calibration vs estimation. No calibrated parameter is a statistical measure. The index in the data was simulated with a Heston-type variance with κ = 8, ξ = 0.40 and ρ = −0.72 under the real-world measure (data README); calibrated to prices, Heston asks for ξ ≈ 2.0 and ρ = −0.75, because prices carry a risk premium and the surface does not have a Heston shape.

In bank practice they are combined: stochastic local volatility (SLV) models use Heston dynamics multiplied by a local function that is calibrated so that today's surface comes out exactly.

Exercises

1. With the calibrated Heston parameters, halve ξ and compute the implied volatility of the 44-day slice at $k = -0.2$, $0$ and $0.1$. What changes in the smile? Does the Feller condition hold now?

kq = np.array([-0.2, 0.0, 0.1])
Tq = np.full(3, 44 / 365)
for nombre, pp in [("calibrated", ph),
                   ("ξ / 2", (v0, kappa, theta, xi / 2, rho_h))]:
    vol = 100 * iv_norm(heston_norm(kq, Tq, pp), kq, Tq)
    feller = 2 * pp[1] * pp[2] - pp[3] ** 2
    print(f"{nombre}: {np.round(vol, 2)}, Feller {feller:.2f}")

calibrated: [38.54 22.03 16.66], Feller -2.97
ξ / 2: [34.2  24.87 19.65], Feller 0.02

2. Calibrate SABR with β = 0.5 instead of 1 on the 16- and 289-day slices and compare the error with that of β = 1. Why does the value of β matter little in equities for fitting a single expiry?

for dias in (16, 289):
    g = x[x["dias"] == dias]
    F, T = g["forward"].iloc[0], g["T"].iloc[0]
    K, ivm = g["strike"].to_numpy(), g["iv"].to_numpy()
    s = least_squares(
        lambda q: hagan(K, F, T, q[0], 0.5, q[1], q[2]) - ivm,
        [0.2 * F ** 0.5, -0.5, 1.0],
        bounds=([1e-3, -0.999, 1e-3], [500, 0.999, 20]))
    print(f"{dias} days, β = 0.5: ρ = {s.x[1]:.3f}, ν = {s.x[2]:.3f}, "
          f"error {100 * np.sqrt(np.mean(s.fun ** 2)):.3f} "
          f"(β = 1: {sabr.loc[dias, 'rmse_sabr']:.3f})")

16 days, β = 0.5: ρ = -0.792, ν = 2.818, error 1.088 (β = 1: 0.945)
289 days, β = 0.5: ρ = -0.693, ν = 0.703, error 0.375 (β = 1: 0.231)

3. A local volatility model calibrated to June 4, 2025 exactly reproduces all the European prices of that day. Explain why that is not enough to price an option that starts in three months with the strike equal to the price on that day (forward start).