---
title: "3.4.7 Heston, SABR and local volatility"
description: "Heston with its characteristic function and QuantLib, SABR with Hagan's formula and Dupire local volatility, calibrated against a known surface."
url: "https://pythonparatrading.com/en/python-trading-course/level-3/heston-sabr-local-volatility/"
lang: "en"
image: "https://pythonparatrading.com/og/en/python-trading-course/level-3/heston-sabr-local-volatility.png"
---

# 3.4.7 · Advanced models: Heston, SABR and local volatility

> The volatility surface describes today's prices; a model says how they move. This lesson programs the three reference models, Heston, SABR and Dupire local volatility, checks them against QuantLib and calibrates them on a synthetic surface whose truth is known, to see what each one fits, what it does not and why the calibrated parameters are not those of the historical series.

- Course: [Level 3: Derivatives: futures, spreads, forex and options](https://pythonparatrading.com/en/python-trading-course/level-3/)
- Module: 3.4 · Option theory and pricing
- Duration: 75 min
- Prerequisites: [3.4.2 · Binomial trees, Black-Scholes-Merton, Black-76 and Monte Carlo](https://pythonparatrading.com/en/python-trading-course/level-3/binomial-black-scholes-black-76-monte-carlo/), [3.4.4 · Implied volatility: Newton-Raphson, Brent, py_vollib and QuantLib](https://pythonparatrading.com/en/python-trading-course/level-3/implied-volatility-newton-brent-quantlib/), [3.4.5 · Smile, skew and term structure; volatility surface with SVI](https://pythonparatrading.com/en/python-trading-course/level-3/volatility-smile-skew-svi-surface/)

## By the end of this lesson

- Explain what each model adds to Black-Scholes: stochastic variance (Heston), stochastic volatility with elasticity (SABR) or volatility as a function of price and time (Dupire).
- Price European options in Heston with its characteristic function and Lewis's formula, and check it with QuantLib's analytic engine.
- Calibrate Heston to a whole surface and SABR expiry by expiry, and measure their error against market volatility.
- Obtain the local volatility of a surface with Dupire's formula in total variance and check it with a Monte Carlo simulation.

## 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:

```mermaid
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
```

| Model | What is random | Parameters | Fits the surface | Typical use |
|---|---|---|---|---|
| Black-Scholes | Only the price | σ | No (one σ) | Quoting in volatility, reference |
| Local volatility (Dupire) | Only the price; σ a function of (S, t) | A whole surface | Exactly, by construction | Equity exotics, price consistency |
| Heston | Price and variance, correlated | v0, κ, θ, ξ, ρ | Approximately, with 5 numbers | Exotics, volatility dynamics |
| SABR | Price and volatility, correlated | α, β, ρ, ν per expiry | Very well for each expiry | Interest rates (swaptions), interpolating smiles |

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.

```python
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")
```
```output
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$:

$$
\frac{dS_t}{S_t} = (r - q)\,dt + \sqrt{v_t}\,dW^S_t, \qquad
dv_t = \kappa(\theta - v_t)\,dt + \xi\sqrt{v_t}\,dW^v_t, \qquad dW^S dW^v = \rho\,dt .
$$

- $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:

$$
\varphi(u) = \exp\big(C(u, T) + D(u, T)\,v_0\big), \quad
d = \sqrt{(\rho\xi iu - \kappa)^2 + \xi^2(iu + u^2)}, \quad
g = \frac{\kappa - \rho\xi iu - d}{\kappa - \rho\xi iu + d},
$$

$$
C = \frac{\kappa\theta}{\xi^2}\left[(\kappa - \rho\xi iu - d)\,T - 2\ln\frac{1 - g e^{-dT}}{1 - g}\right], \qquad
D = \frac{\kappa - \rho\xi iu - d}{\xi^2}\,\frac{1 - e^{-dT}}{1 - g e^{-dT}} .
$$

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

$$
C = D_T\left[F - \frac{\sqrt{FK}}{\pi}\int_0^\infty
\frac{\operatorname{Re}\left[e^{-iuk}\,\varphi(u - i/2)\right]}{u^2 + 1/4}\,du\right].
$$

```python
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}")
```
```output
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:

```python
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())
```
```output
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.

```python
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}")
```
```output
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:

$$
dF_t = \alpha_t F_t^{\beta}\,dW_t, \qquad d\alpha_t = \nu\,\alpha_t\,dZ_t, \qquad dW\,dZ = \rho\,dt .
$$

$\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]$:

$$
\sigma_{\text{imp}} \approx \frac{\alpha\,\big[1 + \big(\tfrac{(1-\beta)^2\alpha^2}{24 (FK)^{1-\beta}}
+ \tfrac{\rho\beta\nu\alpha}{4 (FK)^{(1-\beta)/2}} + \tfrac{2 - 3\rho^2}{24}\nu^2\big)T\big]}
{(FK)^{(1-\beta)/2}\big[1 + \tfrac{(1-\beta)^2}{24}\ln^2\tfrac{F}{K} + \tfrac{(1-\beta)^4}{1920}\ln^4\tfrac{F}{K}\big]}
\cdot\frac{z}{x(z)} .
$$

```python
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}")
```
```output
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.

```python
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)
```
```output
       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).

*Fig. — Root mean square error of Heston and SABR (β = 1) against market volatility, by expiry, ORX, 4 June 2025, in volatility points.* Two curves against maturity in days. Heston (highlighted) is U-shaped: 2.51 points at 16 days, below 0.9 between 44 and 471 days and 1.61 at 653. SABR falls steadily, from 0.95 at 16 days to 0.13 at 653.

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):

$$
\sigma_L^2(k, T) = \frac{\partial w / \partial T}
{1 - \frac{k}{w}\frac{\partial w}{\partial k} + \frac{1}{4}\left(-\frac{1}{4} - \frac{1}{w}
+ \frac{k^2}{w^2}\right)\left(\frac{\partial w}{\partial k}\right)^2 + \frac{1}{2}\frac{\partial^2 w}{\partial k^2}} .
$$

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:

```python
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}")
```
```output
   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$).

*Fig. — Dupire local volatility against implied volatility at 30 days, true SSVI surface of ORX on 4 June 2025: local volatility has almost twice the slope at the money.* Two curves against k, from −0.4 to 0.25. They cross at the money, at 22%. Towards the puts, local volatility (highlighted) rises much faster: 65% against 47% at k = −0.3. Towards the calls, local volatility falls less and turns up earlier: 23% against 19% at k = 0.2.

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

```python
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)
```
```output
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.

<aside class="callout" data-kind="aviso">

For a system that trades listed vanilla options, these models are above all tools for interpolation,
consistency and Greeks. They give no edge on their own: a calibrated model reproduces market prices, it does not
correct them. A strategy that “buys what the model says is cheap” is betting that the model is more
correct than the market, and that has to be tested with a backtest (lesson 3.6.4).

</aside>

## 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?

**Solution.**

```python
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}")
```
```output
calibrated: [38.54 22.03 16.66], Feller -2.97
ξ / 2: [34.2  24.87 19.65], Feller 0.02
```

With half of $\xi$, the smile flattens: the volatility of the put wing falls from 38.5% to 34.2%, that of the
money rises from 22.0% to 24.9% and that of the call wing, from 16.7% to 19.7%. The volatility of variance gives
curvature and, with a negative $\rho$, slope; in addition, in Heston more $\xi$ cheapens the money (ATM implied
volatility is a concave function of variance). Feller goes from −2.97 to 0.02: it just holds.

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?

**Solution.**

```python
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})")
```
```output
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)
```

With $\beta = 0.5$ the error rises a little (from 0.945 to 1.088 at 16 days and from 0.231 to 0.375 at 289) and $\rho$ adjusts
to compensate (−0.79 and −0.69). $\beta$ and $\rho$ both generate slope: $\beta < 1$ makes volatility rise
when the price falls, and a negative $\rho$ does too. With a single expiry, the two effects almost substitute for each other and the
fit does not separate them well; that is why $\beta$ is fixed in advance (1 in equities, often 0.5 or less in rates) and
$\alpha$, $\rho$ and $\nu$ are calibrated.

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*).

**Solution.**

A forward-start option sets its strike in three months at the price of that day, so its value depends on
the smile that will exist then, conditional on the path: on the future skew. Today's European prices only fix the
distribution of the price at each expiry separately, not how they chain together. Local volatility and stochastic
volatility calibrated to the same surface give equal marginal distributions and different future smiles: the
local one, flatter (the function $\sigma_L$ flattens with tenor); the stochastic one, regenerated. Pricing the product
requires choosing the dynamics, and checking it against prices of products that reveal it (the *forward start* itself, cliquets
or options on the VIX) or using a stochastic local model.

## Self-assessment

1. Heston calibrated to the June 4, 2025 surface gives a volatility of variance ξ ≈ 2.0, while the variance of the index in the data was simulated with ξ = 0.40. What does it indicate?

   - A. That the calibration has failed
   - B. That the index has more volatility than was simulated
   - C. That the calibrated parameters are those that reproduce prices (pricing measure) with the Heston structure, not statistical estimates from the series: to generate the short-term skew, Heston needs a much higher volatility of variance
   - D. That more expiries need to be used

   *Answer C: That the calibrated parameters are those that reproduce prices (pricing measure) with the Heston structure, not statistical estimates from the series: to generate the short-term skew, Heston needs a much higher volatility of variance.* The surface in the data is not Heston: it is an SSVI with a risk premium. Heston approximates it with the parameters that best reproduce the prices, and it only achieves the steep skew of the short tenors with a large ξ. Calibrating and estimating are different questions.

2. What does the Feller condition, 2κθ ≥ ξ², do in the Heston model?

   - A. It guarantees there is no arbitrage in the surface
   - B. It makes the correlation negative
   - C. It ensures that the calibration converges
   - D. It guarantees that the variance does not reach zero; if it does not hold, it can touch zero (and rise again), which complicates simulation but does not prevent pricing with the characteristic function

   *Answer D: It guarantees that the variance does not reach zero; if it does not hold, it can touch zero (and rise again), which complicates simulation but does not prevent pricing with the characteristic function.* With 2κθ < ξ² the variance process reaches zero. The semi-analytic formula remains valid; Monte Carlo schemes need to handle zero (truncated or reflected). The lesson's calibration violates it, as is usual for indices.

3. Near the money and at short tenors, how does the skew of local volatility compare with that of implied volatility?

   - A. Local is approximately twice as steep: implied volatility is, to a first approximation, an average of local volatility along the path between spot and strike
   - B. They are equal
   - C. Local is half
   - D. Local has no skew

   *Answer A: Local is approximately twice as steep: implied volatility is, to a first approximation, an average of local volatility along the path between spot and strike.* It is Derman's rule: the implied volatility of a strike averages local volatility between the current price and that strike. In the data, the ratio of at-the-money slopes comes out between 1.80 and 1.86 at 8, 30 and 180 days.

4. In SABR with β = 1, the calibrated ν falls with tenor from 3.0 (16 days) to 0.58 (653 days). What pattern does it hide?

   - A. None: it is calibration noise
   - B. That ν·√T changes much less (from 0.64 to 0.78): the surface in the data scales the smile with the standard deviation, and a per-expiry SABR reproduces it with a ν that falls approximately like 1/√T
   - C. That the volatility of volatility rises with tenor
   - D. That β = 0 should be used

   *Answer B: That ν·√T changes much less (from 0.64 to 0.78): the surface in the data scales the smile with the standard deviation, and a per-expiry SABR reproduces it with a ν that falls approximately like 1/√T.* With β = 1 and a similar α at all tenors, the curvature and slope of the smile in k depend on ν√T. The SSVI in the data has the same shape in standard deviations at all tenors (lesson 3.4.5), so ν√T varies little; what deviates at the short tenors is the weight of the very far wings in the fit.

## Key points

- Heston makes the variance stochastic (reversion κ to θ, volatility of variance ξ, correlation ρ with the price); its characteristic function has a closed form and Lewis's formula prices a European option with one integral. It agrees with QuantLib to the twelfth digit.
- Calibrated to the whole surface, Heston leaves root mean squared errors of 0.5 to 2.5 volatility points per expiry: it cannot have both the steep 16-day skew and the gentle one at almost two years. Its calibrated parameters (ξ ≈ 2.0) are not those of the series (ξ = 0.40).
- SABR is calibrated expiry by expiry with Hagan's formula, which agrees with QuantLib; with β = 1 it fits each slice with errors of 0.13 to 1.0 points, smaller the longer the tenor.
- Dupire local volatility comes from the implied surface: its variance is the derivative of w with respect to tenor divided by g(k). It reproduces all of today's prices (a simulation checks it) and its skew near the money is 1.8 times the implied one.
- Fitting today's prices is not enough: local and stochastic models can give the same surface and different dynamics (future skew, hedging), and that is where they differ for exotic options and for the delta.

## Resources

### Research

- [Stochastic Processes - Random Walks](https://pythonparatrading.com/en/stochastic-processes/) · Research-0005 · 11 min

### Documentation

- [QuantLib 1.44: Heston model (hestonmodel.hpp)](https://github.com/lballabio/QuantLib/blob/v1.44/ql/models/equity/hestonmodel.hpp) · github.com
- [QuantLib-Python: pricing engines (AnalyticHestonEngine)](https://quantlib-python-docs.readthedocs.io/en/latest/pricing_engines.html) · quantlib-python-docs.readthedocs.io

### Books

- [Heston, S. L. (1993): A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options. Review of Financial Studies 6(2), 327-343](https://doi.org/10.1093/rfs/6.2.327) · doi.org
- [Hagan, P. S., Kumar, D., Lesniewski, A. S. and Woodward, D. E. (2002): Managing Smile Risk (Wilmott Magazine article, PDF copy)](https://www.next-finance.net/IMG/pdf/pdf_SABR.pdf) · next-finance.net
- [Dupire, B. (1994): Pricing with a Smile. Risk, January 1994](https://www.risk.net/derivatives/equity-derivatives/1500211/pricing-with-a-smile) · risk.net

### Data

- [Level 3 data: weekly option surface on ORX (Parquet)](https://pythonparatrading.com/curso/datos/nivel-3/opciones3_superficie.parquet) · Download · opciones3_superficie.parquet

Previous lesson: [3.4.6 · Realized vs implied volatility; volatility risk premium; VIX and VSTOXX](https://pythonparatrading.com/en/python-trading-course/level-3/realized-vs-implied-volatility-vix/)

Next lesson: [3.5.1 · Directional and vertical structures: call, put, covered call, collar and spreads](https://pythonparatrading.com/en/python-trading-course/level-3/directional-and-vertical-option-spreads/)