Heston–Nandi (2000) — closed-form GARCH option pricing

Replicating "A Closed-Form GARCH Option Valuation Model" (Review of Financial Studies)

The first application of the volatility arc to derivatives. Heston & Nandi (2000) specify a GARCH whose generating function is closed-form, so European options price by Fourier inversion — no simulation. Physical dynamics (log return $R_t$, per day):

$$R_t = r + \lambda h_t + \sqrt{h_t}\,z_t,\qquad h_t=\omega+\beta h_{t-1}+\alpha\big(z_{t-1}-\gamma\sqrt{h_{t-1}}\big)^2,\quad z_t\sim N(0,1),$$

with leverage $\gamma>0$ (a large $z_{t-1}$ drop raises $h_t$). Under the locally risk-neutral measure ($\lambda^{*}=-\tfrac12,\ \gamma^{*}=\gamma+\lambda+\tfrac12$) the moment generating function is exponential-affine in the current variance, $$f(\varphi)\equiv E^{*}_t\big[S_T^{\varphi}\big]=S_t^{\varphi}\,e^{A_t(\varphi)+B_t(\varphi)\,h_{t+1}},$$ where the scalar $A_t$ and the variance-loading $B_t$ accumulate the effect of every period between now and expiry. They are built by stepping backward from the terminal condition $A_T=B_T=0$, one period at a time: $$B_t=\varphi\,(\lambda^{*}+\gamma^{*})-\tfrac12(\gamma^{*})^2+\beta B_{t+1}+\frac{\tfrac12(\varphi-\gamma^{*})^2}{1-2\alpha B_{t+1}},\qquad A_t=A_{t+1}+\varphi r+\omega B_{t+1}-\tfrac12\ln\!\big(1-2\alpha B_{t+1}\big).$$ (this backward recursion, seeded at expiry and iterated back to today, is hn_gf in the module — $B_t$ is the loading on $h_{t+1}$, $A_t$ collects the rest). The European call then follows by Fourier inversion of $f$: $$C=\tfrac12 S+\frac{e^{-rT}}{\pi}\!\int_0^\infty\!\mathrm{Re}\!\Big[\tfrac{K^{-i\varphi}f(i\varphi+1)}{i\varphi}\Big]d\varphi-Ke^{-rT}\Big(\tfrac12+\frac1\pi\!\int_0^\infty\!\mathrm{Re}\!\Big[\tfrac{K^{-i\varphi}f(i\varphi)}{i\varphi}\Big]d\varphi\Big).$$

All of this is in heston_nandi.py. We fit it on S&P 500 returns and read off the implied-volatility skew the model generates.

Notation (used throughout)

symbol	meaning
$S$	spot — current price of the underlying (S&P 500 / SPY); $S_t$ at time $t$, $S_T$ at expiry
$K$	strike (exercise price) — the fixed price in the contract
$F=S\,e^{(r-q)T}$	forward — price for delivery at expiry ($\approx S$ at short horizons)
$T,\ r,\ q$	time to expiry (years), risk-free rate, dividend yield
$\sigma$	volatility (annualized); the one input Black–Scholes cannot observe

A call pays $\max(S_T-K,\,0)$ at expiry, a put $\max(K-S_T,\,0)$.

• Moneyness $K/S$ (or $K/F$) places the strike relative to spot: $\;=1$ is at-the-money (ATM); $\;<1$ = low strikes (out-of-the-money (OTM) puts / in-the-money calls); $\;>1$ = high strikes (OTM calls / ITM puts). It is scale-free, so one axis compares any underlying or date.
• Implied volatility (IV) is the single $\sigma$ that, put into Black–Scholes, reproduces an option's price. If returns were exactly log-normal it would be flat across strikes; its downward slope (the skew) or U-shape (the smile) is precisely the model's / market's departure from log-normality — which is what every figure below plots.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import heston_nandi as HN
R = pd.read_csv("sp500ret.csv")["ret"].values / 100.0            # decimal daily log-returns
rf = 0.05 / 252                                                  # 5% annual risk-free, per day
theta = HN.hn_mle(R, r=rf); o, a, b, g, lam = theta
persist = b + a * g**2; h_unc = (o + a) / (1 - persist)
S = 100.0
mart = HN.hn_gf(np.array([1.0 + 0j]), S, h_unc, rf, o, a, b, g, lam, 30)[0].real
print("Heston-Nandi GARCH fit:")
print("  omega %.2e  alpha %.2e  beta %.4f  gamma %.2f  lambda %.2f" % (o, a, b, g, lam))
print("  variance persistence beta+alpha*gamma^2 = %.4f   unconditional annual vol %.1f%%" % (persist, 100*np.sqrt(252*h_unc)))
print("  risk-neutral martingale check: f(1) = %.4f  vs  S e^{rT} = %.4f  (must match)" % (mart, S*np.exp(rf*30)))

Heston-Nandi GARCH fit:
  omega 2.98e-19  alpha 8.44e-06  beta 0.9231  gamma 24.71  lambda -6.37
  variance persistence beta+alpha*gamma^2 = 0.9283   unconditional annual vol 17.2%
  risk-neutral martingale check: f(1) = 100.5970  vs  S e^{rT} = 100.5970  (must match)

# implied-volatility skew the model generates, across maturities; and the role of leverage (gamma)
mony = np.linspace(0.85, 1.15, 31); Kg = S * mony
def iv_curve(th, T, h1=h_unc):
    C = HN.hn_call(S, Kg, rf, T, th, h1)
    return np.array([HN.bs_iv(C[i], S, Kg[i], rf, T) for i in range(len(Kg))]) * np.sqrt(252) * 100
fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
for T, col in [(30, "steelblue"), (90, "seagreen"), (180, "firebrick"), (360, "purple")]:
    ax[0].plot(mony, iv_curve(theta, T), color=col, lw=1.8, label="%d-day" % T)
ax[0].axvline(1, color="0.7", lw=.7)
ax[0].set_xlabel("moneyness  K / S"); ax[0].set_ylabel("Black–Scholes implied vol (annualized %)")
ax[0].set_title("Heston–Nandi implied-vol SKEW by maturity\n(down-sloping — crash protection is dearer)"); ax[0].legend(fontsize=8)
# leverage drives the skew: fitted gamma vs gamma = 0
th0 = np.array([o, a, b, 0.0, lam])
ax[1].plot(mony, iv_curve(theta, 30), color="firebrick", lw=2, label="fitted $\\gamma=%.0f$ (skew)" % g)
ax[1].plot(mony, iv_curve(th0, 30), color="steelblue", lw=2, ls="--", label="$\\gamma=0$ (near-symmetric)")
ax[1].axhline(100*np.sqrt(252*h_unc), color="0.6", lw=1, ls=":", label="Black–Scholes (flat)")
ax[1].axvline(1, color="0.7", lw=.7); ax[1].set_xlabel("moneyness  K / S"); ax[1].set_ylabel("implied vol (annualized %)")
ax[1].set_title("The skew IS the leverage $\\gamma$ (30-day)"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig("heston_nandi.png", dpi=120, bbox_inches="tight"); plt.show()

Reading the figure

What implied volatility is. Any option price can be quoted as the single Black–Scholes volatility that reproduces it — its implied volatility (IV). If markets were truly Black–Scholes (one constant volatility), every strike and maturity would return the same IV and these plots would be flat lines. Any slope or curvature is the model saying returns are not log-normal. Both panels invert each Heston–Nandi option price back into its BS-implied vol and plot it against moneyness $K/S$ (strike over spot: $K/S=1$ is at-the-money (ATM, the center); $K/S<1$ = low strikes / out-of-the-money (OTM) puts; $K/S>1$ = OTM calls).

Left — the skew and its term structure. Each curve is one maturity (30 → 360 days).

• Downward slope (the "equity skew"). Low-strike options carry higher IV than high-strike ones. A GARCH with leverage puts extra probability on large down moves, so out-of-the-money puts (crash insurance) are worth more than Black–Scholes says — and "worth more" reads as higher implied vol. This down-sloping shape is the single most robust feature of index-option markets since 1987.
• Flattening with maturity. The 30-day curve is steepest; by 360 days it is much flatter. Short-dated options are dominated by the current variance and its asymmetric reaction to the next few shocks, so leverage bites hardest; over long horizons the variance mean-reverts to $\omega/(1-\text{persistence})$ and the per-period asymmetry averages out, pulling the smile toward flat.

Right — the skew IS the leverage $\gamma$. Same 30-day maturity, three curves: the fitted model (red) shows a pronounced skew; setting $\gamma=0$ (blue dashed) removes the leverage and the curve collapses to a nearly symmetric smile; the dotted line is the flat Black–Scholes IV (the model's unconditional vol). The entire tilt of the surface is carried by one number — the same $\gamma$ that GJR-GARCH estimated from the returns in Phase 2. The faint residual U-shape at $\gamma=0$ is the fat tails a GARCH generates even without leverage.

Data note: the parameters are fit to S&P 500 daily returns (sp500ret.csv), not to option quotes; the notebook then reads off the option surface the fitted returns-model implies. Calibrating directly to a market option chain is a separate exercise.

Model vs market: overlaying the real skew

The skew above came from S&P returns, risk-neutralized — it is what the return dynamics imply an option surface should look like. Does the market agree? Here we overlay the real SPY option skew (spy_options.csv) on the GARCH model's skew.

To make the comparison fair the model is re-fit to contemporaneous SPY returns (2020–2026, spy_returns.csv) — the same era as the option snapshot — and started at today's conditional variance. Both curves then describe the same market at the same moment: one built from how returns actually moved, the other from what options actually cost.

Terminology: a strike at today's (forward) price is at-the-money (ATM), $K/F=1$ — the center of the plot; strikes below are OTM puts / ITM calls, above are OTM calls / ITM puts. "ATM vol" is simply the implied vol at $K/F=1$.

import heston_nandi as HN, spx_options as X
Rspy = pd.read_csv("spy_returns.csv")["ret"].values / 100.0        # contemporaneous SPY returns (2020-2026)
rf2 = 0.043 / 252
th2 = HN.hn_mle(Rspy, r=rf2); o2, a2, b2, g2, lam2 = th2
hp, _ = HN.hn_filter(th2, Rspy, rf2)
z_last = (Rspy[-1] - rf2 - lam2*hp[-1]) / np.sqrt(hp[-1])
h_now = o2 + b2*hp[-1] + a2*(z_last - g2*np.sqrt(hp[-1]))**2        # today's conditional variance
persist2 = b2 + a2*g2**2
print("HN fit to recent SPY returns: persist %.4f  gamma %.1f  lambda %.2f" % (persist2, g2, lam2))
print("  current conditional vol %.1f%%   unconditional %.1f%%"
      % (100*np.sqrt(252*h_now), 100*np.sqrt(252*(o2+a2)/(1-persist2))))

df2, spot2, asof2 = X.load("spy_options.csv")
def hn_skew(T):
    Kg = 100.0*np.linspace(0.90, 1.10, 25)
    C = HN.hn_call(100.0, Kg, rf2, T, th2, h_now)
    return Kg/100.0, np.array([HN.bs_iv(C[i], 100.0, Kg[i], rf2, T) for i in range(len(Kg))]) * np.sqrt(252) * 100

fig, ax = plt.subplots(1, 2, figsize=(13.5, 4.8))
for T, dte, axi in [(30, 29, ax[0]), (60, 60, ax[1])]:
    mm, ivmod = hn_skew(T); s = X.smile(df2, dte)
    axi.plot(s["mny"], 100*s["iv"], "o", ms=3.5, color="steelblue", alpha=.6, label="market (SPY %d-day)" % dte)
    axi.plot(mm, ivmod, "-", color="firebrick", lw=2, label="GARCH model (risk-neutral)")
    axi.axvline(1, color="0.7", lw=.7); axi.set_xlim(0.90, 1.10)
    axi.set_title("%d-day   (market ATM %.1f%%  vs  GARCH ATM %.1f%%)"
                  % (dte, 100*np.interp(1, s["mny"], s["iv"]), np.interp(1, mm, ivmod)))
    axi.set_xlabel("moneyness  K / F"); axi.set_ylabel("annualized implied vol (%)"); axi.legend(fontsize=8)
    sk_mod = np.interp(0.95, mm, ivmod) - np.interp(1.05, mm, ivmod)
    sk_mkt = 100*(np.interp(0.95, s["mny"], s["iv"]) - np.interp(1.05, s["mny"], s["iv"]))
    print("  %d-day skew (IV.95 - IV1.05):  model %.1f pts   market %.1f pts" % (dte, sk_mod, sk_mkt))
plt.tight_layout(); plt.show()

HN fit to recent SPY returns: persist 0.9210  gamma 57.8  lambda -9.69
  current conditional vol 15.0%   unconditional 17.8%
  29-day skew (IV.95 - IV1.05):  model 1.8 pts   market 7.2 pts
  60-day skew (IV.95 - IV1.05):  model 1.1 pts   market 5.7 pts

Reading the overlay — the volatility risk premium

• The model reproduces the sign of the skew, not its steepness. The GARCH skew slopes down — leverage is working — but only gently: a 29-day 0.95→1.05 slope of ~2 vol points, versus ~7 in the market. The market skew is roughly four to five times steeper at both maturities. Out-of-the-money puts (crash protection) trade far richer than the historical return asymmetry, risk-neutralized, can justify.
• That gap is a risk premium. A returns-fitted diffusive GARCH captures the physical asymmetry of returns; the market skew additionally prices the premium investors pay to hedge tail risk — and the jumps a pure diffusion omits. The distance between the red line and the blue points is, in essence, the variance / jump risk premium in the option surface — the same effect that historically makes index puts profitable to sell.
• A regime caveat on the level. Here the GARCH ATM vol (~15–16%) sits slightly above the market's (~13–14%): the model is fit to 2020–2026 returns — which include the 2020 crash and the 2022 bear market — so it remembers more volatility than the currently-calm market is pricing. The level reflects realized history; the market prices a calm forward. The robust, level-independent message is the shape: the market skew is much steeper.

This is exactly why heston1993 calibrates directly to the option surface rather than to returns — to reproduce the market's skew you must fit the price of risk, not only the dynamics of returns.

Results

• The pricing is exact. The risk-neutral generating function satisfies the martingale identity $f(1)=S e^{rT}$ to machine precision, so the closed-form call is correct — no Monte-Carlo error.
• The model produces the equity-index skew. Implied volatility slopes downward in strike (out-of-the-money puts / low strikes carry higher IV), steepest at short maturities and flattening as $T$ grows — exactly the shape Heston–Nandi obtained and that the SPX market shows. A flat Black–Scholes line cannot.
• The skew is the leverage parameter $\gamma$. Setting $\gamma=0$ collapses the curve to a near-symmetric smile; the fitted $\gamma$ is what tilts it into a skew. The same asymmetry that GJR-GARCH found in the returns (Phase 2) is what prices the options — the volatility arc and the option surface are the same $\gamma$ seen two ways.

This is the option-pricing counterpart of the risk arc: closed-form GARCH pricing from the same returns model. Natural follow-ups: Duan (1995) Monte-Carlo GARCH pricing (the LRNVR by simulation), and Heston (1993) stochastic-volatility pricing (the SV bridge) — and calibrating to a live SPX option chain.