Implied Volatility Surface Calibration & Cross-Model Repricing

📊 Description

From a real SPY option-chain snapshot, this project builds an arbitrage-aware implied-volatility surface (SVI), derives the Dupire local-volatility surface from it, and reprices the whole chain under three models (Black-Scholes, Dupire local volatility via forward PDE, and Heston stochastic volatility via characteristic function), comparing them by implied-vol RMSE across moneyness and maturity buckets.

All quantities live in log-forward-moneyness $k = \ln(K / F)$ with $F = S_0 e^{r\tau}$, which centres every smile and is the natural coordinate for both SVI and Dupire.

🎯 Pipeline

Stage	Module	Role
Data fetch	fetch_data.py	Pull the SPY option chain from yfinance, save a dated raw CSV snapshot.
Market data	market_data.py	Load, clean and filter quotes; compute mid price, log-moneyness, maturity.
Implied vol	black_scholes.py	Closed-form price, vega, and implied-vol inversion (Newton + Brent fallback).
Surface	vol_surface.py	Slice-by-slice SVI calibration; linear-in-tau interpolation of total variance.
Local vol	dupire.py	Gatheral's local-variance formula → pre-computed sigma_loc(k, tau) grid.
PDE pricer	pde_pricer.py	Forward Crank-Nicolson scheme; prices the whole call surface in one sweep.
Heston	heston.py, heston_calibration.py	Lewis characteristic-function pricer and calibration.
Repricing	repricing.py	Cross-model repricing and bucketed RMSE comparison.
Plots	plotting.py	Surface, smiles, term structure.
Orchestration	run_pipeline.py	End-to-end run: data → IV → calibration → repricing.

📐 Mathematical Model

1. Implied-volatility surface (SVI)

Each maturity slice is fitted in total variance $w(k) = \sigma_{\text{iv}}^2(k),\tau$ with Gatheral's raw SVI parametrisation:

$$w(k) = a + b\left[\rho(k - m) + \sqrt{(k - m)^2 + \sigma^2}\right]$$

Five parameters per slice, calibrated by least squares with bounds ($a \ge 0$, $b \ge 0$, $\rho \in (-1, 1)$, $\sigma > 0$). Between slices, total variance is interpolated linearly in $\tau$, which preserves calendar monotonicity $\partial_\tau w \ge 0$ and hence rules out calendar-spread arbitrage. SVI is chosen over splines because it is analytic in $k$ (essential for the Dupire derivatives below) and smooth by construction.

2. Dupire local volatility

The local variance is computed from Gatheral's formula in $(k, w)$ coordinates:

$$\sigma_{\text{loc}}^2(k, \tau) = \frac{\partial_\tau w}{1 - \frac{k}{w}\partial_k w + \frac{1}{4}\left(-\frac{1}{4} - \frac{1}{w} + \frac{k^2}{w^2}\right)(\partial_k w)^2 + \frac{1}{2}\partial_{kk} w}$$

Because the surface is SVI, $\partial_k w$ and $\partial_{kk} w$ are closed-form:

$$\partial_k w = b\left[\rho + \frac{k - m}{\sqrt{(k-m)^2 + \sigma^2}}\right], \qquad \partial_{kk} w = \frac{b,\sigma^2}{\left[(k - m)^2 + \sigma^2\right]^{3/2}}$$

The $\tau$-derivative uses a central finite difference of the interpolated surface. The result is a 2D grid sigma_loc[i, j], floored at a small positive value where the denominator degenerates (deep wings).

3. Forward PDE pricer (Crank-Nicolson)

The call surface $C(k, \tau)$ is propagated forward in maturity from the payoff, with a Crank-Nicolson discretisation whose diffusion coefficient is the state- and time-dependent local variance $\sigma_{\text{loc}}^2(k, \tau)$. At each step a tridiagonal system is solved (scipy.linalg.solve_banded); unlike the constant-volatility scheme, the coefficient vectors are rebuilt at every step from the local-vol column. Dirichlet boundaries: discounted intrinsic $S_0 - K e^{-r\tau}$ deep ITM, $0$ deep OTM.

4. Heston stochastic volatility

$$dS_t = r S_t,dt + \sqrt{v_t},S_t,dW_t^S, \qquad dv_t = \kappa(\theta - v_t),dt + \sigma_v \sqrt{v_t},dW_t^v, \qquad d\langle W^S, W^v\rangle = \rho,dt$$

European options are priced via the Lewis (2001) characteristic-function integral (Little Heston Trap formulation for branch-cut stability). The five parameters $(\kappa, \theta, \sigma_v, \rho, v_0)$ are calibrated to the SVI surface by vega-weighted least squares.

📈 Results

Snapshot: SPY, 12 June 2026, spot ≈ 740, 16 calibrated maturities (0.13–1.26y). Implied-vol RMSE of each model against the market:

Model	Global RMSE
Black-Scholes (ATM vol)	15.4%
Heston (stochastic vol, CF)	8.0%
Dupire (local vol, PDE)	5.8%

Bucketed RMSE (maturity × moneyness):

Maturity	Moneyness	BS	Dupire	Heston
Short (<3M)	ITM Put / OTM Call (k<0)	22.3%	6.2%	13.4%
Short (<3M)	ATM	2.1%	1.4%	1.7%
Short (<3M)	OTM Put / ITM Call (k>0)	7.9%	8.4%	8.6%
Medium (3M-1Y)	ITM Put / OTM Call (k<0)	17.2%	4.8%	4.6%
Medium (3M-1Y)	ATM	1.4%	1.2%	1.2%
Medium (3M-1Y)	OTM Put / ITM Call (k>0)	10.8%	11.1%	11.2%
Long (>1Y)	ITM Put / OTM Call (k<0)	19.5%	7.5%	3.9%
Long (>1Y)	ATM	2.8%	1.4%	2.0%
Long (>1Y)	OTM Put / ITM Call (k>0)	7.6%	5.5%	7.8%

Reading the table. Black-Scholes here is a deliberate naive baseline: a single implied vol (the 6-month ATM level) applied to every quote, ignoring both the term structure and the skew. That is the point - it sets the floor the smile-aware models improve on. Black-Scholes collapses on the low-strike side (up to 22%, the k<0 column), where the put skew lives - a single flat vol cannot bend. Both smile models recover it. Dupire wins in-sample by construction: the local-vol surface is built to reproduce the very quotes it is scored against, so its edge is expected, not surprising - it is a consistency check that the calibration → Dupire → PDE chain is sound. Heston, with only five parameters for the whole surface, cannot fit every slice but captures the skew through $\rho \approx -0.73$ and is strongest on the long-dated wing.

A note on Heston: the calibrated parameters violate the Feller condition $2\kappa\theta &gt; \sigma_v^2$, which is common on equity-index surfaces - the high vol-of-vol needed to match the skew pushes the variance process toward zero. The CF pricer remains stable thanks to the Little Heston Trap formulation.

Benchmark prices on a few reference quotes (market mid vs models) make the wing effect concrete:

Benchmark	Market	BS	Dupire	Heston
3M ATM	23.63	25.77	24.12	24.91
1Y ATM	51.88	52.05	55.19	53.70
3M 10% OTM put	8.62	3.85	8.72	9.14
1Y 10% OTM put	33.72	22.27	33.66	32.60
3M 10% OTM call	2.00	4.90	2.02	1.86

BS underprices the 1Y OTM put by more than 11 (22.27 vs 33.72, it cannot see the put skew) while Dupire and Heston track the market closely.

🧮 Data & Numerical Boundary

The snapshot covers all listed SPY expirations (calls and puts), fetched via yfinance and committed under data/ for reproducibility. Prices are the bid-ask mid; spot is the last close at snapshot time.

SPY (the ETF) lists American-style options (exercisable any time before expiry), while the whole pipeline assumes European-style exercise (payoff at expiry only). The early-exercise premium is small for index options with moderate dividends, most visible on deep in-the-money puts; switching to SPX (European-style, cash-settled) would remove the approximation (see Perspectives).

pandas is confined to the I/O and data-cleaning boundary (fetch_data, market_data, the repricing result tables). Every numerical core (black_scholes, vol_surface, dupire, pde_pricer, heston) operates on NumPy arrays exposed by the MarketSnapshot dataclass.

🚀 Usage

python fetch_data.py     # one-off: refresh the raw snapshot
python run_pipeline.py   # full pipeline

✅ Tests

pytest -q

🔜 Perspectives

• Butterfly (static) arbitrage control on the SVI fit (Gatheral–Jacquier conditions).
• Out-of-sample repricing (calibrate on a subset, score on held-out quotes).
• SPX (European-style, cash-settled) snapshot to remove the early-exercise approximation.
• Heston Monte Carlo cross-check of the CF pricer (already implemented in heston.py).

🔗 Related Repositories

Finite-difference schemes for the Black-Scholes PDE (explicit, implicit with Thomas vs SciPy solve_banded, and Crank-Nicolson), whose Crank-Nicolson solver is generalised here to the Dupire local-volatility PDE.

👨‍💻 Author

Alexandre R. - Université Paris Cité