Duan (1995) — Monte-Carlo GARCH option pricing¶

Replicating "The GARCH Option Pricing Model" (Mathematical Finance)¶

Duan's seminal paper priced options under GARCH by simulation, using the locally risk-neutral valuation relationship (LRNVR): under the pricing measure the conditional mean becomes $r-\tfrac12 h_t$ and the volatility asymmetry shifts by the unit risk premium ($\gamma^{*}=\gamma+\lambda+\tfrac12$ in Heston–Nandi notation). One then simulates risk-neutral price paths and averages the discounted payoff, $$C = e^{-rT}\,\mathbb E^{*}\big[(S_T-K)^+\big]\approx e^{-rT}\frac1M\sum_{m=1}^M (S_T^{(m)}-K)^+ .$$

Unlike a closed form, this works for any GARCH and any payoff. We (1) validate it against the Heston–Nandi closed form on the same model, (2) show convergence, (3) reproduce the smile, and (4) price a path-dependent option that has no closed form at all.

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.
In [1]:
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
rf = 0.05 / 252
theta = HN.hn_mle(R, r=rf); o, a, b, g, lam = theta
h1 = (o + a) / (1 - b - a * g**2); S = 100.0; T = 30
Ks = np.array([90., 95., 100., 105., 110.])
Ccf = HN.hn_call(S, Ks, rf, T, theta, h1)                        # closed form (Heston-Nandi)
Cmc, se = HN.hn_mc_call(S, Ks, rf, T, theta, h1, M=400000, seed=1)   # Duan Monte-Carlo (risk-neutral sim)
print("Duan Monte-Carlo vs Heston-Nandi closed form (30-day, S=100):")
print("%6s  %11s  %-18s %8s" % ("strike", "closed-form", "MC (+/- SE)", "diff"))
for i, K in enumerate(Ks):
    print("%6.0f  %11.4f  %8.4f +/- %.4f   %+.4f" % (K, Ccf[i], Cmc[i], se[i], Cmc[i] - Ccf[i]))
print("\n-> every MC price agrees with the closed form to within Monte-Carlo error.")

Duan Monte-Carlo vs Heston-Nandi closed form (30-day, S=100):
strike  closed-form  MC (+/- SE)            diff
    90      10.6231   10.6231 +/- 0.0090   +0.0000
    95       6.0921    6.0927 +/- 0.0080   +0.0006
   100       2.6241    2.6239 +/- 0.0058   -0.0002
   105       0.7805    0.7814 +/- 0.0033   +0.0009
   110       0.1625    0.1630 +/- 0.0015   +0.0005

-> every MC price agrees with the closed form to within Monte-Carlo error.
In [2]:
fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
# (A) convergence of the ATM MC price to the closed form as M grows
cf_atm = HN.hn_call(S, 100.0, rf, T, theta, h1)
Ms = [1000, 3000, 10000, 30000, 100000, 300000, 1000000]
mc = []; err = []
for m in Ms:
    p, s = HN.hn_mc_call(S, 100.0, rf, T, theta, h1, M=m, seed=7); mc.append(p); err.append(1.96 * s)
ax[0].axhline(cf_atm, color="firebrick", lw=1.5, label="closed form")
ax[0].errorbar(Ms, mc, yerr=err, fmt="o-", color="steelblue", capsize=3, label="Monte-Carlo $\\pm$95% CI")
ax[0].set_xscale("log"); ax[0].set_xlabel("simulated paths M"); ax[0].set_ylabel("ATM call price")
ax[0].set_title("Monte-Carlo converges to the closed form"); ax[0].legend(fontsize=8)
# (B) the smile: MC vs closed form
mony = np.linspace(0.90, 1.10, 15); Kg = S * mony
Ccf_g = HN.hn_call(S, Kg, rf, T, theta, h1); Cmc_g, _ = HN.hn_mc_call(S, Kg, rf, T, theta, h1, M=800000, seed=3)
iv_cf = np.array([HN.bs_iv(Ccf_g[i], S, Kg[i], rf, T) for i in range(len(Kg))]) * np.sqrt(252) * 100
iv_mc = np.array([HN.bs_iv(Cmc_g[i], S, Kg[i], rf, T) for i in range(len(Kg))]) * np.sqrt(252) * 100
ax[1].plot(mony, iv_cf, color="firebrick", lw=2, label="closed form")
ax[1].plot(mony, iv_mc, "o", color="steelblue", ms=5, label="Monte-Carlo")
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("Same implied-vol skew, two ways"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig("duan_mc.png", dpi=120, bbox_inches="tight"); plt.show()
No description has been provided for this image

Reading the figure¶

Left — Monte-Carlo convergence. The horizontal red line is the exact closed-form at-the-money price; the blue points are the Duan Monte-Carlo estimate as the number of simulated paths $M$ grows (log scale), with 95% confidence intervals. Two things to notice: the estimate is unbiased — points scatter around the true line at every $M$, not creeping up from one side — and its error shrinks like $1/\sqrt{M}$ (the interval halves for each 4× increase in paths, the signature of Monte-Carlo). By $M=10^6$ the band is a fraction of a cent. Since simulation and closed form price the same model, they must agree, and they do.

Right — the same skew, two ways. The red line is the closed-form implied-vol skew (identical to heston_nandi); the blue dots are the Monte-Carlo prices inverted to implied vol. They sit on top of each other. The point is not a new fact about the skew — it is the same by construction — but that the simulation engine is trustworthy, which licenses using it where no closed form exists.

That payoff is the next cell (no figure needed): it prints a European call — which has a closed form to check against — and an arithmetic-average Asian call, which does not. Monte-Carlo prices both from the same paths; the Asian is cheaper because averaging the price over the month dampens the terminal dispersion the payoff sees.

Data note: same as heston_nandi — parameters come from S&P 500 returns, not option quotes; this notebook validates the simulation engine against the closed form, then extends it to a payoff the closed form cannot handle.

European vs Asian: a path-dependent payoff¶

The options priced so far are European — the payoff looks at the underlying at a single instant, expiry: $\max(S_T-K,\,0)$. Nothing about the path in between matters, which is exactly why closed forms exist.

An Asian option instead pays on the average price over its life: $\max(\bar S-K,\,0)$ with $\bar S=\frac1n\sum_{t}S_t$. That payoff is path-dependent — the whole trajectory matters — so the average of (log-normal) prices is no longer log-normal and no closed form exists. Monte-Carlo prices it for nothing extra: it already simulates every path, so we just average along each path before applying the payoff.

An Asian is cheaper than the otherwise-identical European because averaging smooths the terminal dispersion the payoff sees (a lucky spike on the last day is diluted by every earlier day). They are common where the natural exposure is an average rather than an endpoint — commodities, FX, energy — and they blunt expiry-day manipulation. The next cell prices both from the same risk-neutral GARCH paths.

In [3]:
# the payoff of Monte-Carlo: price a PATH-DEPENDENT option with no closed form (arithmetic Asian call)
def hn_mc_asian(S, K, r, T, theta, h1, M=400000, seed=0):
    o, a, b, g, lam = theta; gs = g + lam + 0.5
    rng = np.random.default_rng(seed)
    Zh = rng.standard_normal((int(T), M // 2)); Z = np.concatenate([Zh, -Zh], axis=1); Mn = Z.shape[1]
    h = np.full(Mn, h1); logS = np.full(Mn, np.log(S)); acc = np.zeros(Mn)
    for t in range(int(T)):
        z = Z[t]; logS += r - 0.5 * h + np.sqrt(h) * z; h = o + b * h + a * (z - gs * np.sqrt(h))**2
        acc += np.exp(logS)
    avg = acc / T; disc = np.exp(-r * T); pay = np.maximum(avg - K, 0.0)
    return disc * pay.mean(), disc * pay.std() / np.sqrt(Mn)
euro, se_e = HN.hn_mc_call(S, 100.0, rf, T, theta, h1, M=400000, seed=5)
asian, se_a = hn_mc_asian(S, 100.0, rf, T, theta, h1, M=400000, seed=5)
print("30-day ATM, GARCH risk-neutral:")
print("  European call  %.4f +/- %.4f   (also has a closed form)" % (euro, se_e))
print("  Asian call     %.4f +/- %.4f   (average-price; NO closed form -- only simulation)" % (asian, se_a))
print("  Asian is cheaper (averaging over the path dampens the terminal dispersion) -- ratio %.2f" % (asian / euro))

30-day ATM, GARCH risk-neutral:
  European call  2.6215 +/- 0.0058   (also has a closed form)
  Asian call     1.5306 +/- 0.0034   (average-price; NO closed form -- only simulation)
  Asian is cheaper (averaging over the path dampens the terminal dispersion) -- ratio 0.58

Results¶

  • Duan's Monte-Carlo reproduces the closed form exactly. On the Heston–Nandi model every simulated price matches the analytic value to within the Monte-Carlo standard error, and the ATM estimate converges to it as $M$ grows (the CI band shrinks onto the closed-form line). The implied-vol skew comes out identically the two ways.
  • The point of simulation is generality. The closed form needs the affine GARCH structure; Duan's LRNVR + Monte-Carlo prices any GARCH (GJR, EGARCH, Student-t innovations) and any payoff. We price a path-dependent Asian call — which has no closed form under GARCH — from the same simulated paths; it is cheaper than the European because averaging over the path dampens the terminal dispersion.
  • Together with heston_nandi.ipynb this is the two-sided view of GARCH option pricing: the closed form (fast, exact, restricted to affine models) and Monte-Carlo (universal, with a controllable error). Next: Heston (1993) — the continuous-time stochastic-volatility closed form (the SV bridge).