Heston (1993) — stochastic-volatility closed-form option pricing

Replicating "A Closed-Form Solution for Options with Stochastic Volatility" (Review of Financial Studies)

The continuous-time counterpart of the GARCH option models — and the bridge to the SV phase of the arc. Volatility is its own diffusion, correlated with the price: $$dS_t=rS_t\,dt+\sqrt{v_t}\,S_t\,dW_1,\qquad dv_t=\kappa(\theta-v_t)\,dt+\sigma\sqrt{v_t}\,dW_2,\qquad \mathrm{corr}(dW_1,dW_2)=\rho.$$

$\kappa$ is the mean-reversion speed, $\theta$ the long-run variance, $\sigma$ the vol-of-vol, and $\rho$ the leverage correlation (negative for equities). Heston's contribution is that the characteristic function of $\ln S_T$ is closed-form, so the European call is $$C=S\,P_1-Ke^{-rT}P_2,\qquad P_j=\tfrac12+\frac1\pi\int_0^\infty\mathrm{Re}\Big[\frac{e^{-i\varphi\ln K}f_j(\varphi)}{i\varphi}\Big]d\varphi,$$ priced by Fourier inversion (heston1993.py, using the numerically stable "little trap" form). This is the same characteristic-function machinery as Heston–Nandi, now in continuous time.

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 heston1993 as H
S, r = 100.0, 0.02
v0, kappa, theta, sigma, rho = 0.04, 2.0, 0.04, 0.4, -0.7        # ~20% long-run vol, strong leverage
# validation 1: sigma -> 0 collapses Heston to Black-Scholes with vol sqrt(theta)
cf0 = H.heston_call(S, 100.0, r, 1.0, theta, kappa, theta, 1e-6, 0.0)
print("sigma->0 limit:  Heston %.4f  vs  Black-Scholes(vol=sqrt(theta)) %.4f  (must match)"
      % (cf0, H.bs_call(S, 100.0, r, 1.0, np.sqrt(theta))))
# validation 2: closed form vs Euler Monte-Carlo of the SDE
cf = H.heston_call(S, 100.0, r, 1.0, v0, kappa, theta, sigma, rho)
mc, se = H.heston_mc(S, 100.0, r, 1.0, v0, kappa, theta, sigma, rho, M=300000, steps=300, seed=1)
print("closed form vs Monte-Carlo (ATM, 1y):  %.4f  vs  %.4f +/- %.4f  (agree within MC error)" % (cf, mc, se))

sigma->0 limit:  Heston 8.9161  vs  Black-Scholes(vol=sqrt(theta)) 8.9160  (must match)

closed form vs Monte-Carlo (ATM, 1y):  8.5230  vs  8.5399 +/- 0.0194  (agree within MC error)

mony = np.linspace(0.80, 1.20, 25); Kg = S * mony
def iv_curve(T, rr=rho, ss=sigma, v=v0):
    C = H.heston_call(S, Kg, r, T, v, kappa, theta, ss, rr)
    return np.array([H.bs_iv(C[i], S, Kg[i], r, T) for i in range(len(Kg))]) * 100
fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
for T, col, lb in [(1/12, "steelblue", "1-month"), (0.25, "seagreen", "3-month"), (0.5, "goldenrod", "6-month"), (1.0, "firebrick", "1-year")]:
    ax[0].plot(mony, iv_curve(T), color=col, lw=1.8, label=lb)
ax[0].axvline(1, color="0.7", lw=.7); ax[0].set_xlabel("moneyness  K / S"); ax[0].set_ylabel("implied vol (%)")
ax[0].set_title("Heston implied-vol SKEW by maturity\n(steep and short-dated, flattening with T)"); ax[0].legend(fontsize=8)
# the skew is the leverage rho (3-month)
for rr, col, ls in [(0.0, "steelblue", "--"), (-0.4, "seagreen", "-"), (-0.7, "firebrick", "-")]:
    ax[1].plot(mony, iv_curve(0.25, rr=rr), color=col, ls=ls, lw=2, label="$\\rho = %.1f$" % rr)
ax[1].axvline(1, color="0.7", lw=.7); ax[1].set_xlabel("moneyness  K / S"); ax[1].set_ylabel("implied vol (%)")
ax[1].set_title("The skew IS the leverage correlation $\\rho$ (3-month)\n($\\rho=0$: symmetric smile;  $\\rho<0$: skew)"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig("heston1993.png", dpi=120, bbox_inches="tight"); plt.show()

Reading the figure

Parameters are illustrative. These are representative equity-index values — long-run vol $\sqrt{\theta}=20\%$, mean-reversion $\kappa=2$, vol-of-vol $\sigma=0.4$, leverage $\rho=-0.7$ — not a calibration to a live option surface. The aim is to show what shapes the model can produce, not to fit a particular day's prices. As before, each option price is inverted to its Black–Scholes implied volatility and plotted against moneyness $K/S$.

Left — skew and term structure. Each curve is a maturity (1 month → 1 year). Heston reproduces the same qualitative surface as the GARCH models: a downward skew, steepest at 1 month and flattening toward 1 year as the variance mean-reverts to $\theta$ — the continuous-time analogue of GARCH leverage.

Right — two shapes, two sources. Three 3-month curves at different $\rho$. With $\rho=0$ (blue dashed) the price and its variance move independently, so the surface is a symmetric smile — and that U-shape comes purely from $\sigma>0$: random volatility fattens both tails and lifts the wings. Making $\rho$ more negative ($-0.4$, then $-0.7$) tilts that smile into a skew by concentrating probability on joint (price-down, vol-up) moves. So the two features have two distinct causes: $\sigma$ (vol-of-vol) sets the curvature/smile; $\rho$ (leverage) sets the tilt/skew. This $\rho$ is the same leverage seen as $\gamma$ in the GARCH options above and as the return/vol-shock correlation in the discrete SV-leverage model.

Does it fit reality? Calibrating to a live option surface

Everything above used illustrative parameters. The real test of an option model is whether it can match a market surface — the grid of implied vols across strikes and maturities that actually trades. Here Heston is calibrated to a real SPY option chain (a Yahoo Finance snapshot cached in spy_options.csv; the mechanics live in spx_options.py).

Three steps turn raw quotes into something calibratable:

1. Forward, not spot. Index options live on the forward $F$. For each maturity we recover $F$ as the strike where the call and put mid-prices cross ($C-P=0$) — a discount-free way to locate it that sidesteps the noise in free quote data.
2. Prices → implied vols. Each out-of-the-money option's mid price is inverted (forward Black formula) to a market implied vol, giving the observed smile at each maturity.
3. One model, whole surface. A single parameter set $(v_0,\kappa,\theta,\sigma,\rho)$ is fit to all four maturities at once by least squares. Fitting the surface jointly (not one smile in isolation) is what identifies the mean-reversion structure — the term structure of the skew is what pins down $\kappa$ and $\theta$.

import spx_options as X                                    # helpers: market smiles + Heston calibration
df, spot, asof = X.load("spy_options.csv")                 # SPY option-chain snapshot (pulled from yfinance)
DTES = [15, 29, 60, 90]
smiles = {d: X.smile(df, d) for d in DTES}                 # forward via put-call parity; market IV from mid prices
print("SPY spot %.2f  as of %s" % (spot, asof))
for d in DTES:
    s = smiles[d]
    print("  %2d-day: F %.2f  ATM %.1f%%  skew(0.90-1.05) %.1f pts  n=%d"
          % (d, s["F"], 100*np.interp(1, s["mny"], s["iv"]),
             100*(np.interp(0.90, s["mny"], s["iv"]) - np.interp(1.05, s["mny"], s["iv"])), len(s["K"])))

cal = X.calibrate_surface([smiles[d] for d in DTES])       # ONE (v0,kappa,theta,sigma,rho) for the whole surface
print("\nHeston calibrated JOINTLY to the full SPY surface (15/29/60/90-day):")
print("  v0 %.4f (%.1f%%)  kappa %.2f  theta %.4f (%.1f%%)  sigma(vol-of-vol) %.2f  rho %+.2f"
      % (cal["v0"], 100*np.sqrt(cal["v0"]), cal["kappa"], cal["theta"], 100*np.sqrt(cal["theta"]), cal["sigma"], cal["rho"]))
print("  overall IV RMSE %.2f%%   per-maturity %s" % (cal["rmse"], ["%.2f" % r for r in cal["rmse_per"]]))

SPY spot 744.78  as of 2026-07-02
  15-day: F 745.53  ATM 12.0%  skew(0.90-1.05) 14.8 pts  n=136
  29-day: F 746.68  ATM 13.4%  skew(0.90-1.05) 12.3 pts  n=168
  60-day: F 749.05  ATM 14.1%  skew(0.90-1.05) 9.6 pts  n=173
  90-day: F 750.06  ATM 14.7%  skew(0.90-1.05) 8.5 pts  n=176

Heston calibrated JOINTLY to the full SPY surface (15/29/60/90-day):
  v0 0.0132 (11.5%)  kappa 13.06  theta 0.0359 (19.0%)  sigma(vol-of-vol) 1.59  rho -0.65
  overall IV RMSE 0.39%   per-maturity ['0.60', '0.41', '0.12', '0.28']

p = (cal["v0"], cal["kappa"], cal["theta"], cal["sigma"], cal["rho"])
cols = {15: "steelblue", 29: "seagreen", 60: "goldenrod", 90: "firebrick"}; LO, HI = 0.92, 1.08
fig, ax = plt.subplots(1, 2, figsize=(13.5, 4.8))
for d in DTES:                                             # (A) the surface fit: market points vs one Heston
    s = smiles[d]; F, D, T = s["F"], s["D"], s["T"]
    m = (s["mny"] >= LO) & (s["mny"] <= HI)
    ax[0].plot(s["mny"][m][::2], 100*s["iv"][m][::2], "o", ms=3.4, color=cols[d], alpha=.45)
    mg = np.linspace(LO, HI, 55)
    ivm = [X.heston_iv(F, mm*F, D, T, *p, nphi=1600) for mm in mg]
    ax[0].plot(mg, 100*np.array(ivm), "-", color=cols[d], lw=2.0, label="%d-day" % d)
ax[0].axvline(1, color="0.7", lw=.7); ax[0].set_xlim(LO, HI)
ax[0].set_xlabel("moneyness  K / F"); ax[0].set_ylabel("implied vol (%)")
ax[0].set_title("Heston calibrated to the REAL SPY surface\n(points = market, lines = one Heston fit;  "
                r"$\rho=%+.2f$, RMSE %.2f%%)" % (cal["rho"], cal["rmse"])); ax[0].legend(fontsize=8, title="maturity")
dd = np.array(DTES)                                        # (B) ATM term structure: market vs model
atm_mkt = [100*np.interp(1, smiles[d]["mny"], smiles[d]["iv"]) for d in DTES]
atm_mod = [100*X.heston_iv(smiles[d]["F"], smiles[d]["F"], smiles[d]["D"], smiles[d]["T"], *p, nphi=1600) for d in DTES]
ax[1].plot(dd, atm_mkt, "o-", color="steelblue", lw=1.6, label="market ATM vol")
ax[1].plot(dd, atm_mod, "s--", color="firebrick", lw=1.6, label="Heston ATM vol")
ax[1].axhline(100*np.sqrt(cal["theta"]), color="0.6", ls=":", lw=1, label=r"$\sqrt{\theta}$ (long-run)")
ax[1].set_xlabel("maturity (days)"); ax[1].set_ylabel("at-the-money implied vol (%)")
ax[1].set_title("ATM term structure: model tracks the market\n(rises toward the long-run level, calm-market contango)")
ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()

Reading the calibration

• Heston fits the real surface. One parameter set reproduces all four SPY smiles to an average 0.39% implied-vol error (left panel) and tracks the at-the-money term structure rising from ~12% to ~15% toward the long-run level $\sqrt\theta=19\%$ (right panel) — the gentle upward slope of a calm market.
• The leverage is real, and it is the same number. The fitted $\rho=-0.65$ sits squarely in the range we have seen all along: the illustrative $\rho=-0.7$ above, the Heston–Nandi / GJR $\gamma>0$ estimated from returns, and the SV-leverage $\rho\approx-0.6$. The asymmetry the option market prices and the asymmetry in the return series are the same economic effect, recovered from completely different data.
• Sensible, with the usual caveats. Short-run variance $v_0$ ($11.5\%$) sits below the long-run $\sqrt\theta$ ($19\%$), consistent with the upward vol term structure; the vol-of-vol $\sigma=1.59$ is large enough to break the Feller condition ($2\kappa\theta<\sigma^2$), a well-known feature of equity-index Heston fits (the variance can graze zero — it does not break the pricing). Heston also cannot bend the very short-dated skew as sharply as the 15-day market shows — the limitation that motivates rough-volatility models.
• What "real data" costs. This is a single snapshot of SPY (an ETF, ~1/10th of SPX) from a free feed with unsynchronized quotes; it demonstrates the method, not a production calibration. A trading desk would use synchronized SPX quotes, more maturities, and liquidity weighting. Re-running spx_options on a fresh snapshot repeats the whole exercise.

The same model, two calibration targets: returns vs the surface

The fit above targeted what options cost. What if we point the same Heston model at how returns actually moved instead — estimating its parameters from the SPY return series (2020–2026) rather than the option chain? Level $\theta$, current variance $v_0$, mean-reversion $\kappa$ and vol-of-vol $\sigma$ come from moment-matching an EWMA variance proxy. The leverage $\rho$ cannot be read off a magnitude-based variance proxy (it needs the return sign), so we take $\rho=-0.6$ — the value the SV-leverage notebook (sv_ext_numpyro) recovers from these same returns and, as the reading above noted, essentially the option-implied $-0.65$. Pricing the smile with those physical parameters as-is gives the smile the return dynamics alone imply.

# ── The SAME Heston, two calibration targets: the option surface vs the return series ──
r = pd.read_csv("spy_returns.csv")["ret"].values / 100.0                    # decimal daily, same era as the chain
dt, lam = 1/252.0, 0.94
h = np.empty(len(r)); h[0] = np.var(r)
for t in range(1, len(r)): h[t] = lam*h[t-1] + (1-lam)*r[t-1]**2            # RiskMetrics EWMA variance
vp = h * 252.0                                                              # annualized variance proxy
theta_p, v0_p = vp.mean(), vp[-1]                                           # long-run & today's variance
kappa_p = -np.log(np.polyfit(vp[:-1], vp[1:], 1)[0]) * 252.0                # mean reversion from AR(1) persistence
eps = (vp[1:]-vp[:-1]) - kappa_p*(theta_p-vp[:-1])*dt                       # variance innovations
sigma_p = np.std(eps / (np.sqrt(vp[:-1])*np.sqrt(dt)))                      # vol-of-vol
rho_p = -0.60                                                              # physical leverage (SV-leverage nb; = option-implied)

pr  = (v0_p, kappa_p, theta_p, sigma_p,     rho_p)                          # returns-fit  (physical)
ps  = (cal["v0"], cal["kappa"], cal["theta"], cal["sigma"], cal["rho"])     # surface-cal   (risk-neutral)
prX = (v0_p, kappa_p, theta_p, cal["sigma"], rho_p)                        # returns-fit but with the OPTION vol-of-vol

def sk(p, s): F,D,T=s["F"],s["D"],s["T"]; return 100*X.heston_iv(F,.95*F,D,T,*p) - 100*X.heston_iv(F,1.05*F,D,T,*p)
def mkt_sk(s): return 100*(np.interp(.95,s["mny"],s["iv"]) - np.interp(1.05,s["mny"],s["iv"]))

print("Same Heston, two targets   (v0 / kappa / theta / vol-of-vol / rho):")
print("  option SURFACE :  %.1f%% / %5.2f / %.1f%% / %.2f / %+.2f   (IV RMSE %.2f%%)"
      % (100*ps[0]**.5, ps[1], 100*ps[2]**.5, ps[3], ps[4], cal["rmse"]))
print("  2020-26 RETURNS:  %.1f%% / %5.2f / %.1f%% / %.2f / %+.2f"
      % (100*pr[0]**.5, pr[1], 100*pr[2]**.5, pr[3], pr[4]))
print("\n29-day skew (IV.95-IV1.05):  market %.1f   surface %.1f   returns %.1f   returns+option-vol-of-vol %.1f  pts"
      % (mkt_sk(smiles[29]), sk(ps,smiles[29]), sk(pr,smiles[29]), sk(prX,smiles[29])))

fig, ax = plt.subplots(1, 2, figsize=(13.5, 4.8))
for axi, d in zip(ax, [29, 60]):
    s = smiles[d]; F, D, T = s["F"], s["D"], s["T"]
    m = (s["mny"] >= 0.90) & (s["mny"] <= 1.10); mg = np.linspace(0.90, 1.10, 45)
    axi.plot(s["mny"][m], 100*s["iv"][m], "o", ms=3.5, color="0.55", alpha=.55, label="market (SPY)")
    axi.plot(mg, [100*X.heston_iv(F, mm*F, D, T, *ps) for mm in mg], "-",  color="firebrick", lw=2, label="calibrated to surface")
    axi.plot(mg, [100*X.heston_iv(F, mm*F, D, T, *pr) for mm in mg], "--", color="steelblue", lw=2, label="fit to returns")
    axi.axvline(1, color="0.8", lw=.7)
    axi.set_title("%d-day skew:  market %.1f  ·  surface %.1f  ·  returns %.1f pts" % (d, mkt_sk(s), sk(ps,s), sk(pr,s)))
    axi.set_xlabel("moneyness  K / F"); axi.set_ylabel("implied vol (%)"); axi.legend(fontsize=8)
plt.tight_layout(); plt.show()

Same Heston, two targets   (v0 / kappa / theta / vol-of-vol / rho):
  option SURFACE :  11.5% / 13.06 / 19.0% / 1.59 / -0.65   (IV RMSE 0.39%)
  2020-26 RETURNS:  15.3% /  4.21 / 16.8% / 0.33 / -0.60

29-day skew (IV.95-IV1.05):  market 7.2   surface 7.2   returns 2.8   returns+option-vol-of-vol 7.0  pts

Reading the comparison — the skew premium is a vol-of-vol premium

Point the same Heston at returns instead of options and it lands somewhere quite different:

• Level and leverage roughly agree. The returns-fit ATM vol (~15%) is in the market's neighbourhood (~13–14%, a touch high because 2020–26 realized vol was lifted by the COVID crash and the 2022 bear), and the physical leverage $\rho\approx-0.6$ essentially equals the option-implied $-0.65$ — the return series and the option market agree on the direction of the asymmetry.
• But the vol-of-vol does not. The realized return path implies $\sigma\approx0.33$; the option surface prices in $\sigma\approx1.6$ — nearly five times more randomness-of-volatility. So the returns-fit smile is far flatter: a ~2.5–2.8-point skew against the market's ~6–7.
• And that one parameter is essentially the whole gap. Swap only the vol-of-vol from its physical value to the option-implied one — holding level, mean-reversion and leverage fixed — and the returns-fit skew jumps from ~2.8 to ~7.0, right onto the market (the returns + option vol-of-vol figure in the printout).

That difference is the variance risk premium: options embed compensation for the volatility- and jump-risk that the realized path — even one spanning 2020 and 2022 — does not itself exhibit. Fitting to returns recovers the level and the leverage direction but not the price of volatility risk, which is exactly why one calibrates to the option surface, not to returns. It is the pricing-side counterpart to the sv_vs_garch horse race: the same model, judged on what it is asked to explain.

Results

• The closed form is correct. It reduces to Black–Scholes when the vol-of-vol vanishes ($\sigma\to0$, matching to 4 dp) and agrees with a direct Euler simulation of the SDE within Monte-Carlo error — two independent checks on the characteristic-function inversion.
• Stochastic volatility generates the skew. Implied volatility slopes down in strike, steepest at short maturities and flattening as $T$ grows — the same qualitative surface as Heston–Nandi's GARCH and as the SPX market. With $\sigma>0$ the model also curves the far wings into a smile (excess kurtosis from random volatility).
• The skew is the correlation $\rho$. Setting $\rho=0$ gives a symmetric smile; making $\rho$ more negative tilts it into a skew. This $\rho$ is the continuous-time leverage — the same effect that appears as $\gamma$ in GJR/Heston–Nandi GARCH and as the return/vol-shock correlation in the discrete SV-leverage model (sv_ext_numpyro.ipynb). One leverage parameter, three representations: it prices the option skew, sharpens the GARCH vol response to bad news, and correlates the SV shocks.

This closes the option-pricing set — GARCH closed form (Heston–Nandi), GARCH Monte-Carlo (Duan), and SV closed form (Heston) — all pricing the same skew from the same leverage, and all tying back to the volatility models fit earlier in the arc.