Option Pricing under GARCH & Stochastic Volatility

Overview

The application phase of the volatility arc: the same models that fit equity returns (GARCH, GJR, stochastic volatility) are used to price options. Three seminal papers are replicated, each pricing under a volatility model that produces the empirical equity-index skew — and the punchline is that the same leverage parameter runs through all of them: the $\gamma$ that GJR-GARCH found in the returns, the $\rho$ in stochastic volatility, and Heston's continuous-time correlation are one asymmetry, seen three ways — and it is exactly what tilts a symmetric option smile into a skew.

European call by Fourier inversion of the (model-specific) characteristic function $f_j$ — the machinery shared by Heston–Nandi and Heston (1993).

Three papers, one skew

Heston & Nandi (2000) — closed-form GARCH. A GARCH whose moment generating function is exponential-affine in the current variance, so European options price by Fourier inversion — no simulation. Under the locally risk-neutral measure the coefficients $A_t, B_t$ accumulate backward from expiry one period at a time. The risk-neutral MGF satisfies the martingale identity $f(1)=Se^{rT}$ to machine precision (a hard correctness check), and setting the leverage $\gamma=0$ collapses the skew to a near-symmetric smile.

Duan (1995) — Monte-Carlo GARCH. Prices options under any GARCH and any payoff by simulation, via the locally risk-neutral valuation relationship (LRNVR): under the pricing measure the conditional mean becomes $r-\tfrac12 h_t$ and the asymmetry shifts by the unit risk premium. It reproduces the Heston–Nandi closed form to Monte-Carlo error on the same model (the CI band shrinks onto the analytic line), then goes where no closed form exists — pricing a path-dependent Asian call from the same simulated paths.

Heston (1993) — continuous-time stochastic volatility. The continuous-time counterpart and the bridge to the SV phase: variance is its own diffusion $dv_t=\kappa(\theta-v_t)dt+\sigma\sqrt{v_t}\,dW_2$, correlated with the price by $\rho$. The characteristic function of $\ln S_T$ is closed-form, so the call prices by Fourier inversion (the numerically stable "little trap" form). It reduces to Black–Scholes as the vol-of-vol vanishes and matches a direct Euler simulation — with $\rho$ again the parameter that generates the skew. Calibrated to the live SPY surface it fits all four maturities to 0.39% IV RMSE; fitting the same model to returns instead lays bare the variance risk premium — physical vol-of-vol ($\sigma\approx 0.3$) is a fifth of what options price in ($\approx 1.6$), so a returns-fit smile is far flatter (a ~2.5-point skew vs the market's ~7) even though its leverage $\rho\approx-0.6$ matches the option-implied $-0.65$.

Notebooks

Downloads

References

• Black, F. & Scholes, M. (1973). The pricing of options and corporate liabilities. Journal of Political Economy 81(3), 637–654. — the constant-volatility benchmark recovered in the limits
• Cox, J. C., Ingersoll, J. E. & Ross, S. A. (1985). A theory of the term structure of interest rates. Econometrica 53(2), 385–407. — the square-root (CIR) variance process underlying Heston, and the Feller condition
• Heston, S. L. (1993). A closed-form solution for options with stochastic volatility. Review of Financial Studies 6(2), 327–343. — the continuous-time stochastic-volatility pricer
• Duan, J.-C. (1995). The GARCH option pricing model. Mathematical Finance 5(1), 13–32. — the LRNVR risk-neutralization and Monte-Carlo GARCH pricing
• Heston, S. L. & Nandi, S. (2000). A closed-form GARCH option valuation model. Review of Financial Studies 13(3), 585–625. — the affine-GARCH closed-form (Fourier-inversion) pricer
• Carr, P. & Madan, D. B. (1999). Option valuation using the fast Fourier transform. Journal of Computational Finance 2(4), 61–73. — characteristic-function inversion for option prices

Heston–Nandi Pricer — Source Code

"""
Heston-Nandi (2000) closed-form GARCH option valuation, from scratch.

Physical dynamics (log return R_t = ln S_t/S_{t-1}, per period):
    R_t = r + lambda * h_t + sqrt(h_t) * z_t,   z_t ~ N(0,1)
    h_t = omega + beta*h_{t-1} + alpha*(z_{t-1} - gamma*sqrt(h_{t-1}))^2
Variance persistence  beta + alpha*gamma^2 (< 1);  unconditional var = (omega+alpha)/(1 - beta - alpha*gamma^2).

Risk-neutral (LRNVR): lambda* = -1/2,  gamma* = gamma + lambda + 1/2, same (omega,alpha,beta).
The generating function E*_t[S_T^phi] = S_t^phi * exp(A_t + B_t h_{t+1}) has a closed recursion, so European
options price by Fourier inversion (Heston-Nandi 2000, eqs 5-6, 14).
"""
import numpy as np
from scipy.optimize import minimize
from scipy.stats import norm
from scipy.optimize import brentq

_trapz = getattr(np, "trapezoid", getattr(np, "trapz", None))    # numpy 2.x renamed trapz -> trapezoid


# ---------------------------------------------------------------- estimation (physical measure)
def hn_filter(theta, R, r):
    omega, alpha, beta, gamma, lam = theta
    T = len(R); h = np.empty(T); z = np.empty(T)
    h[0] = (omega + alpha) / max(1 - beta - alpha * gamma**2, 1e-6)        # unconditional start
    for t in range(T):
        if h[t] <= 0: h[t] = 1e-12
        z[t] = (R[t] - r - lam * h[t]) / np.sqrt(h[t])
        if t + 1 < T:
            h[t + 1] = omega + beta * h[t] + alpha * (z[t] - gamma * np.sqrt(h[t]))**2
    return h, z


def hn_loglik(theta, R, r):
    omega, alpha, beta, gamma, lam = theta
    if omega <= 0 or alpha < 0 or beta < 0 or beta + alpha * gamma**2 >= 1:
        return -np.inf
    h, z = hn_filter(theta, R, r)
    if np.any(h <= 0) or not np.all(np.isfinite(z)):
        return -np.inf
    return -0.5 * np.sum(np.log(2 * np.pi * h) + z**2)


def hn_mle(R, r=0.0):
    v = np.var(R)
    x0 = np.array([v * 0.1, v * 0.2, 0.6, 2.0, 0.5])                       # omega, alpha, beta, gamma, lambda
    nll = lambda th: -hn_loglik(th, R, r) if np.isfinite(hn_loglik(th, R, r)) else 1e10
    res = minimize(nll, x0, method="Nelder-Mead", options=dict(xatol=1e-10, fatol=1e-6, maxiter=20000))
    return res.x


# ---------------------------------------------------------------- option pricing (risk-neutral)
def hn_gf(cphi, S, h_next, r, omega, alpha, beta, gamma, lam, T):
    """Risk-neutral generating function E*[S_T^phi] for complex phi, T periods ahead."""
    gs = gamma + lam + 0.5                                                 # risk-neutral gamma*
    ls = -0.5
    A = np.zeros_like(cphi); B = np.zeros_like(cphi)
    for _ in range(int(T)):
        Bn = cphi * (ls + gs) - 0.5 * gs**2 + beta * B + 0.5 * (cphi - gs)**2 / (1 - 2 * alpha * B)
        A = A + cphi * r + B * omega - 0.5 * np.log(1 - 2 * alpha * B)
        B = Bn
    return np.exp(cphi * np.log(S) + A + B * h_next)


def hn_call(S, K, r, T, theta, h_next, phimax=120.0, nphi=400):
    """European call price(s); K may be scalar or array. h_next = h_{t+1} (next-period variance)."""
    omega, alpha, beta, gamma, lam = theta
    phi = np.linspace(1e-8, phimax, nphi)
    f1 = hn_gf(1j * phi + 1, S, h_next, r, omega, alpha, beta, gamma, lam, T)
    f0 = hn_gf(1j * phi, S, h_next, r, omega, alpha, beta, gamma, lam, T)
    disc = np.exp(-r * T)
    K = np.atleast_1d(np.asarray(K, float)); out = np.empty(len(K))
    for i, k in enumerate(K):
        i1 = (k ** (-1j * phi) * f1 / (1j * phi)).real
        i2 = (k ** (-1j * phi) * f0 / (1j * phi)).real
        P1 = 0.5 * S + disc / np.pi * _trapz(i1, phi)
        P2 = 0.5 + 1.0 / np.pi * _trapz(i2, phi)
        out[i] = P1 - k * disc * P2
    return out if out.size > 1 else float(out[0])


# ---------------------------------------------------------------- Black-Scholes + implied vol
def bs_call(S, K, r, T, sig):
    if sig <= 0 or T <= 0: return max(S - K * np.exp(-r * T), 0.0)
    d1 = (np.log(S / K) + (r + 0.5 * sig**2) * T) / (sig * np.sqrt(T))
    d2 = d1 - sig * np.sqrt(T)
    return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)


def bs_iv(price, S, K, r, T):
    intrinsic = max(S - K * np.exp(-r * T), 0.0)
    if price <= intrinsic + 1e-10: return np.nan
    try:
        return brentq(lambda s: bs_call(S, K, r, T, s) - price, 1e-6, 5.0, maxiter=200)
    except Exception:
        return np.nan


# ---------------------------------------------------------------- Duan (1995) Monte-Carlo (risk-neutral simulation)
def hn_mc_call(S, K, r, T, theta, h_next, M=200000, seed=0, antithetic=True):
    """Monte-Carlo European call under the risk-neutral GARCH (Duan 1995 LRNVR); returns (price, std-err)."""
    omega, alpha, beta, gamma, lam = theta
    gs = gamma + lam + 0.5                                        # risk-neutral gamma*
    rng = np.random.default_rng(seed)
    if antithetic:
        Zh = rng.standard_normal((int(T), M // 2)); Z = np.concatenate([Zh, -Zh], axis=1)
    else:
        Z = rng.standard_normal((int(T), M))
    Mn = Z.shape[1]
    h = np.full(Mn, h_next); logS = np.full(Mn, np.log(S))
    for t in range(int(T)):
        z = Z[t]
        logS += r - 0.5 * h + np.sqrt(h) * z                     # R_t = r - h/2 + sqrt(h) z  (risk-neutral)
        h = omega + beta * h + alpha * (z - gs * np.sqrt(h))**2   # next-period variance
    ST = np.exp(logS); disc = np.exp(-r * T)
    K = np.atleast_1d(np.asarray(K, float))
    pay = np.maximum(ST[None, :] - K[:, None], 0.0)
    price = disc * pay.mean(1); se = disc * pay.std(1) / np.sqrt(Mn)
    return (price, se) if K.size > 1 else (float(price[0]), float(se[0]))