Volatility

Breeden–Litzenberger. Since $C(K) = e^{-rT}\mathbb E^Q[(S_T - K)^+]$, differentiating twice in $K$ gives $\partial_{KK}C = e^{-rT}q(K)$. Call prices across all strikes therefore determine the entire risk-neutral distribution at that maturity.

Local volatility (Dupire). Given call prices across strikes and maturities, there is exactly one diffusion $dS = rS\,dt + \sigma_{\text{loc}}(S, t)S\,dW$ that reproduces all of them: $\sigma_{\text{loc}}^2(K, T) = \dfrac{2\,(\partial_T C + rK\,\partial_K C)}{K^2\,\partial_{KK}C}$. It fits today's smile exactly, but it tends to predict future smiles that flatten, unlike real markets.

Arbitrage bounds. A smile must keep call prices decreasing and convex in $K$ (no butterfly arbitrage) and increasing in $T$ (no calendar arbitrage). Otherwise $q$ would be negative somewhere.

References: D. Breeden & R. Litzenberger, "Prices of state-contingent claims implicit in option prices", J. Business 51 (1978), 621–651. B. Dupire, "Pricing with a smile", Risk 7 (1994), 18–20. J. Gatheral, The Volatility Surface (Wiley, 2006).

import numpy as np
from math import erf, exp, log, sqrt, factorial
N = lambda x: 0.5 * (1 + erf(x / sqrt(2)))
def bs_call(S, K, T, r, s):
    d1 = (log(S / K) + (r + s * s / 2) * T) / (s * sqrt(T))
    return S * N(d1) - K * exp(-r * T) * N(d1 - s * sqrt(T))
def merton_call(S, K, T, r, s, lam, mu, dj):                 # average of BS prices over the number of jumps
    k = exp(mu + dj * dj / 2) - 1
    return sum(exp(-lam * (1 + k) * T) * (lam * (1 + k) * T) ** n / factorial(n)
               * bs_call(S, K, T, r - lam * k + n * log(1 + k) / T, sqrt(s * s + n * dj * dj / T)) for n in range(60))
def implied_vol(price, S, K, T, r, lo=1e-4, hi=3.0):
    for _ in range(80):                                       # bisection: call prices rise with volatility
        mid = (lo + hi) / 2
        lo, hi = (mid, hi) if bs_call(S, K, T, r, mid) < price else (lo, mid)
    return (lo + hi) / 2
S, T, r, s, lam, mu, dj = 100.0, 0.5, 0.02, 0.15, 1.0, -0.10, 0.10   # one jump a year, 10% down on average
for K in (80, 90, 100, 110, 120):
    print(f"K = {K}: implied volatility {implied_vol(merton_call(S, K, T, r, s, lam, mu, dj), S, K, T, r):.4f}")
Ks = np.arange(40, 201, 0.5)                                  # Breeden-Litzenberger: density = e^{rT} C''(K)
C = np.array([merton_call(S, K, T, r, s, lam, mu, dj) for K in Ks])
density = exp(r * T) * np.diff(C, 2) / 0.5 ** 2
print(f"density integrates to {density.sum() * 0.5:.4f}")

Output when we ran it:

K = 80: implied volatility 0.2333
K = 90: implied volatility 0.2134
K = 100: implied volatility 0.1968
K = 110: implied volatility 0.1858
K = 120: implied volatility 0.1796
density integrates to 1.0000

Pricing PDE. Itô's lemma for two variables (section 4.4) and a hedge in the stock and one other option give $\partial_tV + \tfrac12vS^2V_{SS} + \rho\xi vSV_{Sv} + \tfrac12\xi^2vV_{vv} + rSV_S + \kappa(\theta - v)V_v - rV = 0$ (with zero market price of volatility risk). Volatility itself can't be traded, so the model is incomplete: the stock alone can't hedge an option.

Characteristic function. $\mathbb E[e^{iu\ln S_T}] = \exp\big(C(u,T) + D(u,T)\,v_0 + iu\ln S_0\big)$ with explicit $C$ and $D$. The price is $C = S_0P_1 - Ke^{-rT}P_2$ with $P_j = \tfrac12 + \tfrac1\pi\int_0^\infty \operatorname{Re}[\cdots]\,du$. The "little Heston trap" form, used here, avoids branch-cut errors in the complex logarithm.

Feller condition. If $2\kappa\theta \ge \xi^2$, the variance never reaches 0 (section 5.2's CIR result). Calibrated equity parameters often violate it, which is why simulations need careful schemes (full truncation here, or Andersen's QE scheme).

References: S. Heston, "A closed-form solution for options with stochastic volatility with applications to bond and currency options", Review of Financial Studies 6 (1993), 327–343. H. Albrecher et al., "The little Heston trap", Wilmott Magazine (2007). Gatheral, The Volatility Surface, ch. 2.

import numpy as np
from math import erf, exp, log, sqrt, pi
T, r, v0, kap, th, xi, rho, S0, K = 1.0, 0.02, 0.04, 2.0, 0.04, 0.5, -0.7, 100.0, 90.0
def cf(u):                                                    # characteristic function of ln S_T
    b = kap - rho * xi * 1j * u
    d = np.sqrt(b ** 2 + xi ** 2 * (1j * u + u ** 2)); g = (b - d) / (b + d); e = np.exp(-d * T)
    C = r * 1j * u * T + kap * th / xi ** 2 * ((b - d) * T - 2 * np.log((1 - g * e) / (1 - g)))
    return np.exp(C + (b - d) / xi ** 2 * (1 - e) / (1 - g * e) * v0 + 1j * u * log(S0))
u = np.linspace(1e-6, 200, 4001)
P1 = 0.5 + np.trapz(np.real(np.exp(-1j * u * log(K)) * cf(u - 1j) / (1j * u * cf(-1j))), u) / pi
P2 = 0.5 + np.trapz(np.real(np.exp(-1j * u * log(K)) * cf(u) / (1j * u)), u) / pi
exact = S0 * P1 - K * exp(-r * T) * P2
rng, n, steps = np.random.default_rng(30), 100_000, 200      # Monte Carlo, full-truncation Euler
dt, lnS, v = T / steps, np.full(n, log(S0)), np.full(n, v0)
for _ in range(steps):
    z1 = rng.normal(size=n); z2 = rho * z1 + sqrt(1 - rho ** 2) * rng.normal(size=n); vp = np.maximum(v, 0)
    lnS += (r - vp / 2) * dt + np.sqrt(vp * dt) * z1
    v += kap * (th - vp) * dt + xi * np.sqrt(vp * dt) * z2
pay = exp(-r * T) * np.maximum(np.exp(lnS) - K, 0)
print(f"Heston call K = {K}: Fourier {exact:.4f}   Monte Carlo {pay.mean():.4f} +- {pay.std() / sqrt(n):.4f}")

Output when we ran it:

Heston call K = 90.0: Fourier 15.0921   Monte Carlo 15.0726 +- 0.0414

Spanning formula. Any twice-differentiable payoff $f(S_T)$ can be written, for any $F$, as $f(F) + f'(F)(S_T - F) + \int_0^F f''(K)(K - S_T)^+dK + \int_F^\infty f''(K)(S_T - K)^+dK$ (Carr & Madan). With $f = -\ln$, $f''(K) = 1/K^2$.

Assumptions. The replication holds exactly for continuous price paths with any stochastic volatility. With jumps, realized variance differs from the log-contract value by a term of order $\sum_{\text{jumps}}\big(J^2 - 2(e^J - 1 - J)\big)$, which is positive for crashes: realized variance then exceeds what the strip pays, so the strip underprices the swap. Discretely sampled realized variance and interest rates add small corrections.

Volatility swaps. A swap on volatility $\sqrt{\text{variance}}$ is not replicable by a static strip. By Jensen's inequality its fair strike is below $\sqrt{K_{\text{var}}}$: a convexity adjustment.

References: K. Demeterfi, E. Derman, M. Kamal & J. Zou, "A guide to volatility and variance swaps", J. Derivatives 6 (1999), 9–32. P. Carr & D. Madan, "Towards a theory of volatility trading", in Volatility (Risk Books, 1998). Gatheral, The Volatility Surface, ch. 11.

import numpy as np
from math import exp, sqrt
T, v0, kap, th, xi, rho = 1.0, 0.04, 2.0, 0.04, 0.5, -0.7     # Heston, r = 0
fair = th + (v0 - th) * (1 - exp(-kap * T)) / (kap * T)      # E[(1/T) * integral of v dt]
rng, n, days = np.random.default_rng(31), 20_000, 252
dt, v, realized = T / days, np.full(n, v0), np.zeros(n)
for _ in range(days):
    z1 = rng.normal(size=n); z2 = rho * z1 + sqrt(1 - rho ** 2) * rng.normal(size=n); vp = np.maximum(v, 0)
    ret = -vp / 2 * dt + np.sqrt(vp * dt) * z1                # daily log return
    realized += ret ** 2
    v += kap * (th - vp) * dt + xi * np.sqrt(vp * dt) * z2
realized /= T
print(f"fair variance strike {fair:.5f}   average realized {realized.mean():.5f} +- {realized.std() / sqrt(n):.5f}")
print(f"realized variance ranges from {np.percentile(realized, 5):.4f} to {np.percentile(realized, 95):.4f} (5th-95th percentile)")

Output when we ran it:

fair variance strike 0.04000   average realized 0.03973 +- 0.00022
realized variance ranges from 0.0094 to 0.1019 (5th-95th percentile)