A clean SPX option surface dataset from paired call and put quotes, treasury rates, parity-implied forwards, and bid-ask-aware implied volatility.
A raw PCHIP benchmark surface first, then a smoother tensor-product B-spline surface fitted on log total variance.
A Dupire local-volatility surface, with real no-arbitrage diagnostics instead of only a pretty 3D plot.
A historical panel of fitted volatility surfaces, surface-shape features, local-vol stress metrics, and PCA factors of surface changes.
Surface-aware Greeks, where delta and gamma are recomputed while allowing implied volatility to move with spot.
A secondary BTC option-surface implementation using the library, where the same machinery is tested on a much rougher and more nonlinear option market.
In the earlier implied-volatility work, we inverted an option price to obtain a volatility for one contract. Here we start from a chain of options with many strikes and maturities, so we can study how those implied volatilities fit together across the market. Once we arrange them as a surface, we can ask:
How does implied volatility change across strike and maturity? Is the surface smooth enough to differentiate? Does it imply a valid risk-neutral density? What local volatility process can reproduce the market option prices? How much are Greeks wrong if we treat volatility as flat?
So the central object is no longer just \(\sigma\). The central object is a function:
\[
\sigma_{imp}(k,T)
\]
where \(k\) is log-moneyness and \(T\) is time to expiry. If we use forward log-moneyness, then:
\[
k = \log\left(\frac{K}{F_T}\right)
\]
where \(K\) is strike and \(F_T\) is the forward price for maturity \(T\). This is cleaner than using raw strike because an option with strike 4,500 means something different when SPX is at 4,800 versus 3,800. Log-moneyness normalizes the strike relative to the forward level.
We need a surface that is close enough to market quotes and smooth enough to produce meaningful derivatives. Local volatility and surface-aware Greeks both use those derivatives, so a small wiggle in fitted IV can have a much larger effect on the quantities we calculate from it. We’ll build the surface with that relationship between quote fit and smoothness in mind.
The strike and maturity directions each contribute to the shape, and we’ll examine both as we move from individual observations to a fitted surface. Every observation keeps its own strike and expiry as we compare the quotes.
Realized, implied, and local volatility
We start with three volatility quantities that appear in different parts of the analysis and describe different aspects of the underlying and its option prices.
First, realized volatility is what the underlying actually did over a time period. If SPX daily returns over the next month are volatile, realized volatility will be high.
Second, implied volatility is the volatility number that makes an option-pricing model match a market option price. It is not a direct forecast. It contains expected volatility, risk premium, supply/demand, skew demand, and liquidity effects.
Third, local volatility is a model-implied function. It answers: if the underlying followed a diffusion where volatility depended only on spot and time, what volatility function would reproduce the current option surface?
We start with market implied volatility at many strikes and maturities, smooth those observations into a surface, and differentiate option prices built from it to infer local volatility. Realized volatility describes the underlying’s historical or subsequent movement over its own observation window.
The reason this is hard is that differentiation amplifies noise. A surface can look visually acceptable, but if it contains tiny wiggles, local volatility can become nonsense. So the whole project is built around a practical tension:
Stay close enough to market quotes to be meaningful, but smooth enough that derivatives and risk measures don’t explode.
We’ll keep the raw quote information visible as we build smoother estimates. That lets us compare the fitted smiles with the market observations and then study the quantities obtained from the surface across strikes, maturities, and dates.
1) Setup and Model Objects
We’ll use numerical tools for three parts of the analysis:
Option pricing and implied volatility inversion using Black-76 and the same fast LBR-lite style solver idea from the implied-vol project.
Surface construction using PCHIP and cubic B-splines.
Automatic differentiation using JAX, which becomes important when we compute Dupire local volatility and surface-aware Greeks.
The main mathematical objects are:
Object
Meaning
Why we need it
\(C(K,T)\)
Call price as a function of strike and maturity
Dupire local volatility differentiates call prices
\(\sigma_{imp}(k,T)\)
Implied volatility surface
Market’s volatility map across moneyness and maturity
\(w(k,T)=\sigma_{imp}^2(k,T)T\)
Total implied variance
Better behaved for surface fitting than raw IV
\(\sigma_{loc}(K,T)\)
Local volatility
Instantaneous volatility function that reproduces the option surface
Needed for slope, curvature, local-vol stability, and Greeks
We use \(365.25\) days for annualization in the SPX part. That’s a practical calendar-time convention for options, especially when maturities are counted in calendar days. The BTC secondary part uses \(365\) because crypto trades every day and the code explicitly treats the BTC chain with a 365-day year.
With those objects in mind, we can start from the option quotes and put each contract on a common forward-moneyness scale. The pricing and fitting calculations below will use the rate, forward, strike, and maturity belonging to each observation.
Show code
import warningsimport timeimport numpy as npimport pandas as pdimport matplotlib.pyplot as pltfrom cycler import cyclerfrom matplotlib.colors import TwoSlopeNorm as two_slope_normfrom scipy.interpolate import PchipInterpolator as pchip_interpolatorfrom scipy.interpolate import BSpline as bsplinefrom scipy.linalg import LinAlgError as linear_errorfrom scipy.linalg import lstsq as linear_lstsqfrom scipy.linalg import solve as linear_solvefrom scipy.stats import normfrom sklearn.decomposition import PCA as pca_estimatorfrom sklearn.preprocessing import StandardScaler as standard_scalerimport jaximport jax.numpy as jnpfrom IPython.display import displayfrom quantfinlab.dataio import load_par_yield_curve, load_spx_option_pairsfrom quantfinlab.fixed_income.bootstrap import build_zero_curve_panel_from_par_yieldsfrom quantfinlab.fixed_income.discounting import discount_factor_from_rate, map_curve_rates_to_dates_and_tausfrom quantfinlab.options.iv import compute_iv_tablefrom quantfinlab.options.parity import infer_forwards_from_paired_quotesfrom quantfinlab.plotting import LAB_COLORS as lab_colors, set_plot_stylefrom quantfinlab.plotting.options import plot_forward_vs_spotjax.config.update("jax_enable_x64", True)warnings.filterwarnings("ignore")pd.set_option("display.float_format", lambda x: f"{x:,.6f}")pd.set_option("display.max_columns", 120)set_plot_style(lab_colors)plt.rcParams["axes.prop_cycle"] = cycler(color=lab_colors)plt.rcParams.update({"figure.figsize": (7.4, 4.2), "figure.dpi": 180, "savefig.dpi": 300})seed =7rng = np.random.default_rng(seed)
2) SPX Option Data and Paired Call-Put Quotes
We start from SPX option-chain data and restrict the raw chain before doing anything complicated. For reproducibility, the relevant source folder is the spx_optionsdx folder inside the repository’s data/ directory, and the rate curve comes from the us_treasury_yields folder. Each data source has a folder with a README and a script explaining how to produce the local file from its source. The option-chain and Treasury inputs can therefore be prepared using the workflows documented alongside the data.
The first cleaning layer keeps:
maturities from about 7 to 180 calendar days,
strikes inside a broad moneyness range,
option pairs with both call and put quotes,
reasonably liquid quotes based on relative spreads,
enough strike pairs per expiry to estimate a forward.
We use the call-put pairs to estimate a forward level \(F_T\) for each maturity. That forward sets the log-moneyness coordinate of every strike: if it is estimated incorrectly, the same quotes will be placed at the wrong horizontal positions in the smile.
The paired quote table has one row per call-put pair at the same strike and expiry. The important columns are:
spot: the current SPX index level,
strike: the option strike,
dte_days: days to expiry,
c_mid, p_mid: call and put mid prices,
c_rel_spread, p_rel_spread: bid-ask spread as a fraction of mid price.
The first displayed rows are all from 2022-01-03 with around 8.3 days to expiry. They show strikes slightly below spot and a clear call-put asymmetry: calls are expensive in dollar price because they are in the money, while puts are cheap because they are out of the money. That doesn’t mean calls have higher implied vol. Price level and volatility aren’t the same object. We later invert both sides to IV so quotes can be compared on a common volatility scale.
2.1 Paired quotes and forward moneyness
Before fitting IV across the option chain, we need a moneyness coordinate for each strike. We use the forward rather than the current spot to place an option relative to its maturity-specific underlying price. For strike \(K\) and forward \(F_T\), forward log-moneyness is:
\[
k=\log\left(\frac{K}{F_T}\right)
\]
At these levels, the difference between spot and the maturity-specific forward shifts the location assigned to the same strike. If \(F_T\) is too high, the same strike looks more out-of-the-money than it really is. If \(F_T\) is too low, the same strike looks more in-the-money. That shift changes the estimated skew.
For example, suppose SPX is near 4,800 and a 60-day strike is 4,500. If the correct forward is 4,790, then:
That difference looks small, but when we take derivatives of the surface, small horizontal shifts can matter. A local-volatility calculation differentiates the call price surface. If the strike coordinate is slightly misaligned across maturities, the derivative field can become contaminated by carry/forward error.
Put-call parity is the most direct way to estimate the forward from the option chain itself. We aren’t forcing an external dividend estimate into the surface. We let calls and puts reveal the maturity-specific forward, then we use treasury rates for discounting. This is more internally consistent for an options-based surface.
Before IV inversion, we need the risk-free discount factor and the forward price. We use treasury par yields, bootstrap/interpolate a zero curve, map each option maturity to a rate, and then compute:
\[
D(T)=e^{-r(T)T}
\]
where:
\(D(T)\) is the discount factor to maturity \(T\),
\(r(T)\) is the continuously compounded zero rate for that maturity,
\(T\) is time to expiry in years.
Then we infer the forward from put-call parity. In a forward-based notation, European call and put prices satisfy:
\[
C(K,T)-P(K,T)=D(T)\big(F_T-K\big)
\]
Solving for the forward gives:
\[
F_T = K + \frac{C(K,T)-P(K,T)}{D(T)}
\]
For one expiry, we have many strikes. In perfect frictionless markets, every strike would imply the same \(F_T\). Real quotes have bid-ask noise, discrete strikes, stale quotes, and microstructure errors, so we use the cross-section of paired quotes and summarize the implied forward robustly. The parity_error_iqr column is useful because it tells us how noisy the parity relationship is across strikes.
Once \(F_T\) is estimated, we can infer an implied carry rate:
where \(b(T)\) is the forward carry. For an equity index, \(q(T)\) is mostly the market-implied dividend yield plus other small carry effects. We don’t need this to be a perfect dividend model. We need it to make pricing internally consistent: the forward used in moneyness and the discount factor used in pricing must agree with the call-put data.
For example, if SPX is near 4,795 and the parity-implied 8-day forward is around 4,794, the forward is slightly below spot. That means the short-term carry is negative, which corresponds to an implied dividend yield/carry effect above the short risk-free rate. This is normal for an index with dividends.
3.1 Rate and carry curves across maturities
It is tempting to use one risk-free rate for every option, but that is too crude for a surface project. A 9-day option and a 174-day option don’t have the same discount horizon. We write the rate as \(r(T)\) because the rate is maturity-dependent:
\[
D(T)=\exp(-r(T)T)
\]
The forward is also maturity-dependent:
\[
F_T=S_0\exp(b(T)T)
\]
where \(b(T)\) is the implied carry. For an equity index, if dividends are expected before expiry, the forward can be below the spot even when the risk-free rate is positive. The carry rate packages this into one curve:
\[
b(T)=r(T)-q(T)
\]
We interpolate \(r(T)\) and \(b(T)\) over maturity \(T\) so we can evaluate them on the grid as well as at the quoted expiries. The rate controls discounting, while carry maps the spot level to the forward for each horizon. When local volatility is computed, the call price derivative with respect to maturity uses the rate and carry at each maturity node. If those curves were discontinuous, the Dupire numerator could jump for no real market reason.
Interpolating the rate and carry curves also keeps the option-pricing inputs consistent across the maturity grid. Smooth curves help prevent artificial jumps in model prices and their maturity derivatives.
From the forward/carry table, the first maturity on 2022-01-03 has spot around 4,795.57, rate around 0.05%, forward around 4,794.05, and implied carry around \(-1.39\%\) annualized. The forward is extremely close to spot in level terms because the maturity is only about 8 days, but annualizing the small difference still gives a visible carry number.
Across the sample, the parity-implied short-maturity forwards track the spot index closely. The small difference between the two levels reflects the effect of carrying the underlying until expiry, while the paired quotes give us the forward used for each maturity.
The rate/carry/dividend-yield plot for the latest date shows three curves:
\(r(T)\), the treasury rate curve,
\(b(T)\), the parity-implied carry curve,
\(q(T)=r(T)-b(T)\), the implied dividend-yield curve.
The latest rate and carry curves separate discounting from the forward implied by the paired option quotes at each expiry. We place the options in forward space using the estimated forward for their own maturity. A 30-day option and a 150-day option can then be compared at the same relative moneyness even though their forward levels and discount factors are different.
4) From Paired Quotes to an Option Market Table
After extracting forwards, we split each call-put pair into two option rows: one call row and one put row. This gives a normal option-market table where every row is one traded option quote with:
date,
expiry,
option type,
spot,
forward,
strike,
time to expiry,
rate,
implied carry,
bid, mid, and ask.
This transformation is useful because implied volatility inversion is done quote by quote. A call and a put at the same strike can both be valid, but they can also disagree slightly because of bid-ask noise. Keeping them as separate rows lets the surface fitting process see both observations and weight them by their quality.
The paired observations give us the following option-market coverage:
about 2,407,888 option rows,
505 dates,
502 expiries,
exactly half calls and half puts,
maturity range around 7.3 to 179.3 days.
We have several thousand option observations on each date and enough trading dates to study how the fitted surface changes over time. That coverage supports both the cross-section of strikes and maturities and the historical comparison across 505 dates.
We compute implied volatility with a Black-76 formula using the maturity-specific forward \(F_T\). The forward is related to spot \(S_0\) through carry, and the call price is:
\[
C = D(T)\left(F_T\Phi(d_1)-K\Phi(d_2)\right)
\]
We solve this for bid, mid, and ask prices. The mid IV is used as the main surface observation, while the bid-ask IV range becomes part of the uncertainty measure later. The LBR-lite solver is fast enough to invert the entire 2.4 million-row option table in about 24 seconds using the accelerated engine.
A useful example: two options can have very different dollar prices but similar IVs. An in-the-money call can cost hundreds of index points because most of its value is intrinsic value. A far out-of-the-money put can cost only a few points. IV removes that mechanical price-level difference and asks, “what volatility would make this price fair?” That is why all later surface work is done in IV or total variance, not raw option price.
5.1 IV bounds, existence, and quote sensitivity
An implied volatility exists only if the market price sits inside the no-arbitrage price interval. For a forward-discounted call, the basic lower bound is:
\[
C \ge D(T)\max(F_T-K,0)
\]
and the upper bound is:
\[
C \le D(T)F_T
\]
For puts:
\[
P \ge D(T)\max(K-F_T,0)
\]
and:
\[
P \le D(T)K
\]
The lower bound is approached as volatility tends toward zero, and the upper bound is approached as volatility becomes very large. If a bid, mid, or ask violates these bounds because of stale quotes, rounding, or a very wide market, the IV solver can fail or produce a nonsensical volatility. This is why a few missing IVs are not alarming. They are a data-quality diagnostic.
The sensitivity of IV to price also depends on vega. Newton-style IV solvers rely on the derivative:
When \(\mathcal{V}\) is near zero, the price-volatility curve is almost flat. A tiny price change can imply a huge volatility change. This is common for deep ITM/OTM short-dated options. We give low-vega observations less weight when fitting the surface, since their inverted IVs are especially sensitive to small price errors.
For example, if an option has vega of 2 index points per 100 vol points, a 1-point bid-ask noise can imply a 50-vol-point IV width. But if another option has vega of 500, the same 1-point price noise implies only 0.2 vol points. The surface should trust the second quote more.
5.2 Bid, mid, and ask implied volatilities
We use mid IV as the main point estimate, while the bid and ask IVs give us a range of implied volatilities consistent with the quoted market. Their difference helps us describe the uncertainty around that observation.
If the bid IV is \(\sigma_{bid}\) and the ask IV is \(\sigma_{ask}\), then the width:
\[
\Delta\sigma_{ba}=\sigma_{ask}-\sigma_{bid}
\]
is a direct market-implied uncertainty band, measured in volatility units rather than index-price points. A tight band means market makers agree fairly closely about the volatility level. A wide band means the quote is less informative.
Suppose two options have the same mid IV of 20%, while one has a bid-ask IV width of 0.5 vol points and the other has a width of 20 vol points. The tighter market gives us a more precise observation of the volatility at that strike and maturity. We can use the width to give those quotes different influence in the fitted surface.
This is also why deep wing quotes can be tricky. They sometimes carry important information about tail pricing, but their bid-ask bands can be huge. We don’t want to delete all of them because then the surface loses tail information. We also don’t want to let them dominate. The weighting system is the compromise.
The root-finding process itself also gives information. If bid IV fails but mid and ask succeed, that often means the bid is near the no-arbitrage floor. If ask IV is extremely high, it can mean the market maker is quoting a protective ask rather than a true fair value. We carry these quality fields into the surface-ready table alongside the final IV so the fitting step can distinguish a reliable quote from a wide or failed inversion.
Most quotes invert successfully. Across the full option-market sample, we have:
2,407,888 total rows,
505 dates,
502 expiries,
IV range from about 4.57% to about 100.69% before final surface filtering,
6,460 missing mid IV values.
The 6,460 missing mid IVs make up roughly 0.27% of all option rows across the 505 dates in the sample, leaving most observations available for the surface-cleaning step. They usually come from quotes near arbitrage bounds, very stale bid/ask combinations, or prices whose implied volatility is outside a stable numerical range.
We can compare quote density and median IV across log-moneyness and maturity in the two hexbin panels. The dataset is densest near the money and in the shorter-to-medium maturity region, which is expected because SPX options trade most heavily near relevant strikes. The right panel shows median implied volatility. The strongest feature is the equity-index skew: downside strikes have higher IV than ATM and upside strikes. This is the market price of crash protection.
The plot also reminds us why we later use a support mask. The surface shouldn’t be trusted equally everywhere. A node inside dense quote support is a market-informed estimate. A node outside the quote cloud is mostly extrapolation.
6) Surface-Ready Quotes and Observation Weights
Raw IV observations still need a surface-cleaning layer. We keep quotes with reasonable IV, maturity, strike, forward, rate, discount factor, and log-moneyness. The main log-moneyness filter is:
\[
-0.40 \le k \le 0.40
\]
This is broad enough to include deep wings, but it removes extreme points that are usually too illiquid for stable surface fitting.
We also assign an observation weight so that cleaner quotes can have more influence in the fit. The weight combines four pieces of information about spreads, vega, moneyness, and IV uncertainty:
6.1 Bid-ask quality
A quote with a tight relative spread is more reliable than a quote with a wide spread. If:
\[
s_i=\frac{ask_i-bid_i}{mid_i}
\]
then a basic inverse-spread weight is:
\[
w_{spread,i}\propto \frac{1}{s_i^2+c^2}
\]
where \(c\) prevents the weight from exploding when the spread is very small.
6.2 Vega information
Option prices are most sensitive to volatility near the money. The Black-76 dollar vega is roughly:
If vega is tiny, a small price error can create a huge IV error. That is why very deep ITM/OTM options can produce unstable IVs even if their price looks harmless. We give more weight to observations with meaningful vega.
6.3 Wing penalty
Even after filtering, far-wing quotes are less reliable and less dense. We still want them because they define skew and curvature, but we don’t want them to dominate the fit. So we gently penalize observations as \(|k|\) moves far from the center.
6.4 IV uncertainty
When bid, mid, and ask IVs are available, the bid-ask IV width gives a direct uncertainty measure:
\[
u_i = \sigma_{ask,i}-\sigma_{bid,i}
\]
If IV width isn’t available, price spread divided by vega gives an approximate IV uncertainty:
The final observation weight is a clipped and normalized combination of spread quality, vega information, wing penalty, and uncertainty penalty. It gives cleaner, more informative quotes more influence in the fitted surface while retaining the wing observations that pass the filters.
6.5 Observation weights and measurement uncertainty
A formal statistical surface fit would start with a measurement model:
\[
y_i=f(k_i,T_i)+\epsilon_i
\]
where \(y_i\) is observed log total variance or implied volatility, \(f(k_i,T_i)\) is the true smooth surface, and \(\epsilon_i\) is observation noise. If we knew the variance of the noise, the efficient weighted least-squares weight would be:
We don’t know \(\operatorname{Var}(\epsilon_i)\) directly, so we approximate it from market microstructure. Wide spreads, low vega, very far wings, and wide bid-ask IV intervals all imply larger observation noise.
We combine those market-based proxies into an information score for the fit. The score reflects the relative quality of observations even though we don’t directly observe a separate noise variance for every option quote.
A useful mental model is:
near-ATM, tight-spread, high-vega options are the surface anchors,
moderately OTM options shape the skew,
far-wing options inform tail behavior but get less authority,
extremely noisy quotes remain in the dataset only if they pass filters, but their weight is small.
This is also why the median weight is normalized to 1. The scale of the weights matters less than their relative values. Weighted least squares only needs the fit to know which observations are more reliable.
We can see how the weighting behaves in the surface-ready quotes. Some put quotes have very low iv_uncertainty clipped at 0.003 and low observation weights because they are far-wing or low-vega. Some call quotes have much larger IV uncertainty, including values clipped at 0.50, and get low weights.
The distribution of quote characteristics and weights is:
median relative spread is about 1.48%,
median IV uncertainty is clipped at 0.003,
median dollar vega is around 395,
median observation weight is normalized to 1.0,
95th percentile observation weight is about 2.34.
The normalized weights give more influence to clean, liquid, informative quotes while retaining the wing observations. The 95th percentile is about 2.34 times the median weight, so we can also see the scale of the upper part of the distribution.
The date-level table shows that recent dates have around 5,500 quotes per day and more than 30 expiries. The latest date, 2023-12-29, has 5,554 surface-ready quotes, 33 expiries, median IV around 11.79%, and a very wide downside log-moneyness range. This is enough coverage for a serious surface fit.
7) Main Date, High-IV Date, and Low-IV Date
We choose three important dates:
main date: 2023-12-29,
high-IV date: 2022-06-13,
low-IV date: 2023-12-14.
The main date is the most recent valid date in the SPX dataset. It has 5,554 quotes, 33 expiries, median IV around 11.79%, and near-ATM median IV around 10.82%. This is a low-volatility end-of-2023 surface.
The high-IV date, 2022-06-13, has near-ATM IV around 31.02%. This sits inside the 2022 inflation/rates shock, when equities were selling off, rates were repricing, and option demand for downside protection was much stronger. This gives us a stress reference.
The low-IV date, 2023-12-14, has near-ATM IV around 10.39%, even lower than the main date. This is the calm end-of-2023 regime after the market had largely priced a more favorable soft-landing/rate-cut path.
We mainly fit and visualize the latest date first, then fit every valid date historically. This order is useful: one date lets us understand the mechanics carefully, while the historical panel lets us ask whether the same features and risks persist through time.
Before fitting a smooth parametric surface, we build a raw PCHIP benchmark. PCHIP stands for piecewise cubic Hermite interpolating polynomial. Its important property is that it tends to preserve the shape of the data better than a standard cubic spline. A standard cubic spline can overshoot between points. PCHIP is more local and less likely to create artificial oscillations.
For one expiry, we have IV observations across log-moneyness. PCHIP gives a smooth curve:
\[
\hat{\sigma}_{PCHIP}(k \mid T_j)
\]
for each observed maturity \(T_j\). Then we interpolate across maturity to get a grid:
\[
\hat{\sigma}_{PCHIP}(k_m,T_n)
\]
This raw benchmark is useful because it stays close to the quote cloud. But it has two weaknesses:
It is built slice-by-slice, so the term structure can be less smooth than a full 2D fit.
It isn’t designed for stable derivatives, especially second derivatives with respect to strike.
That second weakness matters a lot. Dupire local volatility depends on \(\partial_T C\), \(\partial_K C\), and \(\partial_{KK}C\). A raw interpolation can match quotes nicely but explode when differentiated. So we use PCHIP as a visual/diagnostic benchmark, not as the final Dupire input.
8.1 PCHIP interpolation and derivative smoothness
PCHIP is useful because it is local. Between two adjacent strike points, it builds a cubic curve using slope information from neighboring points. Unlike a standard cubic spline, it tries to preserve monotonicity and avoid overshoot. In rough terms, for each interval \([k_i,k_{i+1}]\), PCHIP creates:
The coefficients are chosen so the curve passes through the observed points and has controlled slopes at the knots.
This is good for drawing smiles because it doesn’t create wild oscillations between strikes. But local volatility needs more than a nice curve. It needs stable derivatives of call prices across strike and maturity. PCHIP is only piecewise smooth. Its first derivative is continuous, but higher derivatives can behave less smoothly around knots. Also, when we interpolate maturity slice by slice, the two-dimensional smoothness is not globally controlled.
We use the raw PCHIP surface as a benchmark for the shape of the quote cloud. The B-spline fit will give us a separate two-dimensional estimate whose smoothness can be controlled in both moneyness and maturity before we take derivatives.
At the end of 2023, the SPX smiles have a clear downside skew across the observed maturities. Shorter maturities have a steeper left wing: as \(k\) becomes more negative, implied volatility rises quickly. Longer maturities are smoother and less steep.
This is a typical equity-index pattern. Downside crash protection is expensive, especially at shorter maturities where realized jumps matter more. Upside options have lower implied volatility because investors usually don’t pay the same premium for upside convexity as they do for downside insurance.
The raw PCHIP surface plot is useful but also shows the limitation. The surface exists only where we have quote support. There are blank or unstable regions away from the dense quote cloud. The 3D view makes the left-wing wall obvious: IV climbs sharply as strikes move below forward.
At one maturity, we can read a smile across strikes; across maturities, those smiles form the volatility surface. The shape in the \(k\) direction describes skew and the price of downside protection, while the shape in the \(T\) direction describes how implied volatility changes with the option’s horizon.
9) Fitting the Surface on Log Total Variance
Now we move from raw interpolation to a smooth fitted surface. Instead of fitting implied volatility directly, we fit log total variance:
\[
w(k,T)=\sigma_{imp}^2(k,T)T
\]
and:
\[
\ell(k,T)=\log w(k,T)
\]
There are two reasons this is a better target:
Positivity
Total variance must be positive. If we model \(\ell(k,T)\) and then exponentiate, we automatically get:
This avoids negative fitted variances, which can happen if we fit raw variance with a flexible smoother.
Better maturity behavior
For short maturities, IV can move a lot because \(T\) is small. Total variance often behaves more smoothly than volatility itself because it accumulates variance through time. A one-month implied vol and a six-month implied vol aren’t directly comparable in variance units unless we multiply by time.
For example, 20% IV for one month has total variance:
\[
0.20^2\times\frac{1}{12}=0.00333
\]
while 20% IV for six months has total variance:
\[
0.20^2\times\frac{1}{2}=0.02000
\]
Same IV, different total uncertainty. Fitting \(w(k,T)\) respects that difference.
9.1 Total variance and calendar-arbitrage intuition
Total variance is also tied to no-arbitrage logic. For the same strike/moneyness region, longer maturities should generally contain at least as much total uncertainty as shorter maturities. In a simplified setting, this means:
\[
\frac{\partial w(k,T)}{\partial T}\ge 0
\]
where \(w(k,T)=\sigma^2(k,T)T\).
This doesn’t mean implied volatility itself must always rise with maturity. IV can fall while total variance rises. For example:
1-month IV: 30%, total variance \(0.30^2\times\frac{1}{12}=0.0075\),
6-month IV: 20%, total variance \(0.20^2\times\frac{1}{2}=0.0200\).
The 6-month IV is lower, but total variance is higher. That can be perfectly reasonable. A short-term event can make front IV high while longer maturities have lower annualized IV. But the longer option still contains more total variance because it spans more time.
This is why fitting log total variance is cleaner than fitting IV directly. It aligns the fitted object with maturity accumulation. It doesn’t guarantee no calendar arbitrage by itself, but it makes the target more natural.
In local-volatility work, this point becomes critical. Dupire’s numerator is closely related to how option prices and total variance evolve with maturity. If fitted total variance wiggles downward because of noisy quotes, local variance can become negative even if IV levels look visually plausible.
10) Tensor-Product Cubic B-Splines
The smooth surface is represented with B-spline basis functions in both directions. In the log-moneyness direction, define basis functions:
The coefficients \(\beta_{ij}\) form a matrix \(B\) of size \(n_k\times n_T\). Each product \(B_i(k)A_j(T)\) combines one strike basis function with one maturity basis function. Because the individual B-splines have local support, the product and its coefficient affect a particular patch of the surface.
The local support lets us represent different shapes in different parts of the volatility surface. The left wing can be steep while the right wing is flat. The short maturities can be curved while long maturities are smoother. A global polynomial would struggle with that because changing one coefficient can distort the whole surface.
For a quote \(m\) at \((k_m,T_m)\), the fitted value is:
In matrix form, if \(X\) is the design matrix and \(\beta\) is the vectorized coefficient matrix, then:
\[
\hat{\ell}=X\beta
\]
With \(M\) quotes, \(X\) has \(M\) rows and \(n_k n_T\) columns, one for each pair of basis functions. The vector \(\beta\) contains the same coefficients as the matrix \(B\), arranged to match those columns. We can now estimate them using weighted regression with a smoothness penalty.
10.1 Local flexibility in both surface directions
A volatility surface has two axes, but the axes behave differently. The strike direction captures skew and smile curvature. The maturity direction captures term structure. A tensor-product spline lets us control both directions separately.
If \(B_i(k)\) is local in strike and \(A_j(T)\) is local in maturity, then the product:
\[
B_i(k)A_j(T)
\]
is active only in a local patch of the surface. This gives flexibility without turning the fit into a high-degree global polynomial.
For example, suppose the front-end downside wing is unusually steep on a stress date. We want the fit to bend in the short-maturity, low-strike region without forcing the entire six-month surface to bend the same way. Tensor-product B-splines can do that because the coefficient controlling that region is local.
The knot locations determine where neighboring spline pieces meet, and the number of basis functions controls flexibility. More basis functions mean the surface can follow local details. Fewer basis functions mean the surface is smoother but less adaptive. The penalty terms then control how much the coefficients are allowed to bend across neighboring basis functions.
We control the fit at two levels:
Basis size controls maximum flexibility.
Penalty strength controls how much of that flexibility we actually use.
That is more stable than simply choosing a very small basis or a very large basis with no penalty.
11) Penalized Weighted Least Squares
The fitted coefficients solve a penalized weighted least-squares problem:
\(y_m=\log(\sigma_m^2T_m)\) is observed log total variance,
\(x_m\) is the spline basis row for quote \(m\),
\(w_m\) is the observation weight from quote quality,
\(D_k\) applies second-difference penalties across log-moneyness,
\(D_T\) applies second-difference penalties across maturity,
\(\lambda_k\) controls smoothness in the strike direction,
\(\lambda_T\) controls smoothness in the maturity direction.
The second-difference penalty is the discrete version of penalizing curvature. If the coefficient surface bends too much, the penalty increases. A high \(\lambda_k\) makes the smile smoother across strike. A high \(\lambda_T\) makes the term structure smoother across maturity.
There is a real tradeoff here:
If \(\lambda\) is too small, the surface follows noise and derivatives become unstable.
If \(\lambda\) is too large, the surface is too flat and misses real skew/curvature.
For local volatility, smoothness matters more than visual closeness. A tiny wiggle in IV can become a huge spike in \(\partial_{KK}C\). So we fit two surfaces:
Visual surface: closer to market IV, useful for plots and shape summaries.
Dupire surface: smoother, slightly less close to quotes, but safer for derivatives.
We can compare the two fitted surfaces against the same quotes and use the smoother estimate for calculations that depend on first and second derivatives.
11.1 The normal equations behind the fit
We can write the penalized objective in matrix form. If \(W\) is the diagonal matrix of observation weights, \(X\) is the spline design matrix, \(y\) is the vector of observed log total variance, and \(P\) is the combined penalty matrix, then:
\[
P=\lambda_kP_k+\lambda_TP_T
\]
where \(P_k\) penalizes bending across strike coefficients and \(P_T\) penalizes bending across maturity coefficients. Then the objective is:
Taking the derivative with respect to \(\beta\) gives two contributions: the weighted residuals pull the fit toward the observations, while the penalty pulls neighboring coefficients toward a smoother pattern:
\[
\nabla_\beta = -2X^\top W(y-X\beta)+2P\beta
\]
Setting that gradient to zero gives the normal equations:
\[
(X^\top WX+P)\beta=X^\top Wy
\]
So the fitted coefficients are:
\[
\hat{\beta}=(X^\top WX+P)^{-1}X^\top Wy
\]
This equation is helpful because it shows exactly how smoothing works. The data term \(X^\top WX\) asks the surface to match quotes. The penalty term \(P\) adds curvature resistance. If \(P\) is zero, we get ordinary weighted least squares. If \(P\) is very large, the surface becomes very smooth because bending coefficients becomes expensive.
In practice, we solve the linear system rather than explicitly forming the inverse. The implementation also adds a small diagonal term, \(10^{-8}I\), to the normal matrix for numerical stability; it is much smaller than the two fitted smoothness penalties.
11.2 Strike and maturity smoothing penalties
The strike penalty and maturity penalty have different financial meanings.
The strike penalty controls how aggressively the smile can bend. A low \(\lambda_k\) allows sharp skew changes and local wiggles across strikes. That can help match far-wing quotes, but it can also create unstable densities because \(C_{KK}\) becomes noisy. A high \(\lambda_k\) makes the smile smoother, but it can underfit real crash-risk pricing.
The maturity penalty controls how quickly the term structure can change. A low \(\lambda_T\) allows the surface to jump between expiries. That can match event-specific expiries, but it can create unstable \(\partial_TC\). A high \(\lambda_T\) makes maturity evolution smoother, but it can blur real event risk.
For SPX, event risk is usually less concentrated than a single-stock earnings announcement, which gives us room to smooth across nearby expiries. At the same time, macro events and changes in market expectations can still create real differences between the front and longer maturities.
A good surface fit therefore isn’t the one with the smallest RMSE. The best derivative surface is the one that balances:
That is why the Dupire candidate selection includes invalid local-vol nodes, extreme ratios, and density problems. We choose smoothing based on the downstream task, not only the in-sample fit error.
The Dupire surface is selected from a small grid of smoothness pairs \((\lambda_k,\lambda_T)\). For each candidate, we check both fit error and local-volatility stability. The candidate table reports:
weighted RMSE,
hard invalid share,
negative density share,
extreme local-vol-to-implied-vol ratio share,
total flagged share,
median absolute deviation of local vol from implied vol,
a selection score.
The best candidate uses \(\lambda_k=90\) and \(\lambda_T=130\). It has the lowest weighted RMSE among the candidate set, about 0.00549, and the lowest total flagged share, about 4.52%. More smoothing increases the weighted RMSE and doesn’t reduce invalid nodes enough to justify the worse fit.
The visual and Dupire fit table shows the expected tradeoff:
visual weighted RMSE: about 0.00503,
Dupire weighted RMSE: about 0.00549,
both use 5,554 quotes and 33 maturities,
both have p95 absolute IV error around 0.0186.
The smoother Dupire fit has a small increase in weighted error relative to the visual fit, while the 95th-percentile absolute error remains around 0.0186 for both. We can therefore compare derivative stability without losing much accuracy in the observed quote fit.
The grid check also shows both fitted grids are finite across supported nodes, and the JAX-vs-NumPy surface evaluation check has max absolute difference 0.0. That means the automatic-differentiation version of the spline is numerically consistent with the NumPy evaluation used elsewhere.
12.1 Visual fit and Dupire fit
We fit one surface to study the market shape and another, more regularized surface to calculate the derivatives used in Dupire.
The visual fit should be close enough to the market to explain the surface shape. If a trader looks at the 30-day smile, the visual fit should represent the quote cloud clearly. Residuals should be small in the liquid center, and the wings should follow the broad market pattern.
The Dupire fit should be safe enough to differentiate. Local volatility cares about:
not only about \(\sigma(k,T)\). A small oscillation that barely changes IV RMSE can completely change \(\partial_{KK}C\). So a slightly smoother fit can be much better for Dupire even if it is a little worse visually.
This is a common pattern in quantitative finance. The estimator used for visualization isn’t always the estimator used for risk. A visually detailed surface can overfit microstructure noise. A risk surface needs to be more regularized.
The visual fit has lower weighted RMSE, while the Dupire fit is only slightly worse and is selected using the local-volatility diagnostics as well. We can read their errors alongside the stability measures for the same set of quotes.
The PCHIP-vs-spline plot compares raw market smiles with the smooth visual log-total-variance spline. The smoother curves track the broad market shape but don’t chase every local wiggle.
The most obvious feature is the steep downside wing. At short maturities, the fitted IV at \(k\approx -0.20\) is much higher than ATM IV. Across longer maturities, the smile becomes flatter, but the left wing remains above the center. This is the persistent equity-index skew.
The fitted surface follows the broad smile without passing through every raw quote. Next we can measure the difference between market IV and fitted IV across the quote cloud, including the center and the wings. Smoothing reduces the influence of small quote-to-quote fluctuations on the first and second derivatives, while preserving the downside skew that appears across the maturities.
13) Residual Diagnostics
After fitting the surface, we study residuals:
\[
e_m=\sigma_{mkt,m}-\hat{\sigma}(k_m,T_m)
\]
A residual near zero means the fitted surface matches the market quote well. A positive residual means the market quote has higher IV than the fit. A negative residual means the fit is above the market quote.
Residual diagnostics matter for two reasons:
They show whether the surface is systematically biased in a region.
They reveal where quote quality or surface smoothness is causing trouble.
The residual heatmap shows market IV minus fitted IV across log-moneyness and maturity. Most of the dense central region is close to zero, which is good. The strongest residuals are in the far downside wing and near sparse parts of the surface.
The largest residual table makes this very clear. Many of the worst residuals are around strike 3,400 on 2023-12-29, with \(k\approx -0.343\) and maturities around 26 to 44 days. Market IV is around 38% to 45%, while the fitted surface is above 50% or 60%. These are deep downside strikes. The model is smoothing the wing upward relative to some individual quotes.
The weight assigned to each observation also affects the fitted coefficients, so the largest absolute residual need not be the quote with the greatest influence on the surface. In the far-left wing, some put and call quotes disagree and bid-ask spreads can be large. The smoothness penalty pulls those observations into a common shape, so individual residuals can be large even when the fitted center is close to the more densely quoted market.
display(pd.DataFrame([fit_metrics(fit_visual, q_main, "visual"), fit_metrics(fit_dupire, q_main, "dupire")]))q_fit["maturity_bucket"] = pd.cut(q_fit["dte_days"], bins=[7, 21, 45, 90, 180], labels=["7-21d", "21-45d", "45-90d", "90-180d"], include_lowest=True)fig, axes = plt.subplots(1, 2, figsize=(12.6, 4.5))bucket_names = [x for x in q_fit["maturity_bucket"].cat.categories if (q_fit["maturity_bucket"] == x).any()]vals = [q_fit.loc[q_fit["maturity_bucket"].eq(x), "residual_visual"].dropna().to_numpy() for x in bucket_names]axes[0].boxplot(vals, labels=[str(x) for x in bucket_names], showfliers=False)axes[0].axhline(0.0, color="black", lw=0.9)axes[0].set_title("visual-fit residuals by maturity bucket")axes[0].set_xlabel("maturity bucket")axes[0].set_ylabel("iv residual")q_fit["k_bin"] = pd.cut(q_fit["k"], bins=np.linspace(-0.40, 0.40, 17), include_lowest=True)rb = q_fit.groupby("k_bin", observed=False).agg( k_mid=("k", "median"), mean_residual=("residual_visual", "mean"), q25=("residual_visual", lambda x: np.nanquantile(x, 0.25)), q75=("residual_visual", lambda x: np.nanquantile(x, 0.75)),).dropna()axes[1].plot(rb["k_mid"], rb["mean_residual"], marker="o", lw=2.0)axes[1].fill_between(rb["k_mid"], rb["q25"], rb["q75"], alpha=0.20)axes[1].axhline(0.0, color="black", lw=0.9)axes[1].set_xlim(max(float(rb["k_mid"].min()), limits_main["k_lo"]), min(float(rb["k_mid"].max()), limits_main["k_hi"]))axes[1].set_title("visual-fit wing bias by log-moneyness")axes[1].set_xlabel("log strike / forward")axes[1].set_ylabel("iv residual")plt.tight_layout()plt.show()display(q_fit.assign(abs_residual=lambda x: x["residual_visual"].abs()).sort_values("abs_residual", ascending=False)[ ["date", "expiry", "option_type", "strike", "k", "dte_days", "iv_mid", "fit_iv_visual", "fit_iv_dupire","residual_visual", "residual_dupire", "iv_uncertainty"]].head(12))
fit
date
quote_count
number_of_maturities
rmse
weighted_rmse
mae
median_absolute_error
p95_absolute_error
mean_residual
residual_std
k_min
k_max
min_dte
max_dte
0
visual
2023-12-29
5554
33
0.010989
0.005030
0.004965
0.002258
0.018665
0.001035
0.010941
-0.344795
0.117917
9.333333
174.333333
1
dupire
2023-12-29
5554
33
0.011624
0.005491
0.005490
0.002750
0.018626
0.000829
0.011594
-0.344795
0.117917
9.333333
174.333333
date
expiry
option_type
strike
k
dte_days
iv_mid
fit_iv_visual
fit_iv_dupire
residual_visual
residual_dupire
iv_uncertainty
2404433
2023-12-29
2024-01-25
put
3,400.000000
-0.342951
26.333333
0.453259
0.618881
0.627676
-0.165622
-0.174417
0.012389
2404964
2023-12-29
2024-01-30
call
3,400.000000
-0.343128
31.333333
0.419128
0.578432
0.587752
-0.159304
-0.168623
0.500000
2404965
2023-12-29
2024-01-30
put
3,400.000000
-0.343128
31.333333
0.423619
0.578432
0.587752
-0.154813
-0.164133
0.014764
2405269
2023-12-29
2024-02-01
put
3,400.000000
-0.343818
33.333333
0.418148
0.567007
0.576614
-0.148859
-0.158466
0.008424
2405852
2023-12-29
2024-02-07
call
3,400.000000
-0.344244
39.333333
0.385966
0.533230
0.543225
-0.147265
-0.157260
0.500000
2404432
2023-12-29
2024-01-25
call
3,400.000000
-0.342951
26.333333
0.475398
0.618881
0.627676
-0.143483
-0.152278
0.500000
2405753
2023-12-29
2024-02-06
put
3,400.000000
-0.344106
38.333333
0.397564
0.538062
0.547987
-0.140498
-0.150423
0.006863
2405853
2023-12-29
2024-02-07
put
3,400.000000
-0.344244
39.333333
0.394289
0.533230
0.543225
-0.138941
-0.148936
0.009881
2406131
2023-12-29
2024-02-12
put
3,400.000000
-0.344795
44.333333
0.377801
0.511296
0.521563
-0.133495
-0.143762
0.008351
2405268
2023-12-29
2024-02-01
call
3,400.000000
-0.343818
33.333333
0.466613
0.567007
0.576614
-0.100394
-0.110001
0.312205
2404072
2023-12-29
2024-01-23
call
3,700.000000
-0.257632
24.333333
0.309178
0.403774
0.401167
-0.094595
-0.091988
0.500000
2403708
2023-12-29
2024-01-18
call
5,300.000000
0.102142
19.333333
0.156919
0.068783
0.065167
0.088136
0.091751
0.005993
The residual boxplot by maturity shows that most residuals are tight around zero in all maturity buckets, but the short end has wider tails. That’s natural. Short-dated options are more sensitive to small price errors because \(T\) is small. A few cents or index points of price noise can produce a big IV difference.
The wing-bias line shows how residuals behave across log-moneyness. The central region is relatively well behaved. Bias becomes more visible in the wings, especially where quote support becomes thinner.
A useful way to think about the tradeoff is:
The middle of the surface is quote-rich and highly liquid, so the fit can be close. The wings are quote-poor and derivative-sensitive, so the fit has to be more conservative.
The support masks identify the region covered by quotes, while the separate visual and Dupire fits let us inspect the market shape and its derivatives with different smoothing levels. The sparse far wings receive less support than the central part of the surface.
14) Surface Shape Features: ATM Level, Skew, Curvature, and Term Slope
Once we have a fitted surface, we can summarize its shape with economically meaningful features. We evaluate IV at three log-moneyness points:
downside point: \(k_d=-0.20\),
ATM point: \(k_0=0\),
upside point: \(k_u=0.08\).
The ATM level is:
\[
\sigma_{ATM}(T)=\sigma(k_0,T)
\]
At \(k_0=0\), the strike equals the maturity-specific forward, so this is the central IV level for that expiry.
This is our three-point smile-shape measure. The chosen downside and upside points are at different distances from ATM, so it is a descriptive curvature feature rather than an equally spaced numerical second derivative.
Here \(T_{front}\) is the front maturity in the comparison. A positive slope means ATM IV at \(T\) is higher than at that front reference.
These features translate the surface into financial language. ATM IV measures general uncertainty. Downside skew measures crash-protection premium. Curvature measures how smile-shaped the distribution is. Term slope measures whether near-term or longer-term uncertainty is more expensive.
On 2023-12-29, the fitted smile has a particularly steep short-dated downside wing.
At 7 days, ATM IV is about 10.74%, but downside skew is enormous at about 31.68 vol points. That means the \(k=-0.20\) left-wing point is priced around 42% IV. Upside skew is negative, around -2.84 vol points, so upside calls are cheaper in IV than ATM options.
As maturity increases, the downside skew falls:
around 18.59 vol points at 28.6 days,
around 13.03 vol points at 71.9 days,
around 11.34 vol points at 180 days.
This is a classic short-dated crash-skew structure. Near-term deep downside protection is very expensive because short-dated crashes are discontinuous and hard to hedge. Longer-dated protection is still expensive, but the skew is smoother because more time allows uncertainty to spread across scenarios.
The ATM term structure is mildly upward-sloping. ATM IV rises from around 10.7% at 7 days to about 13.5% at 180 days. That tells us short-term realized-vol expectations were calm at the end of 2023, while longer maturities still priced some medium-term uncertainty.
Across the shape-point smiles, all maturities slope downward from left to right, and the short maturity has the steepest left wing. The downside points carry a much larger IV premium than the corresponding upside points in these fitted smiles.
15) Surface Derivatives from Log Total Variance
For local volatility and surface-aware Greeks, we need derivatives. Since we fit:
The \(-1/(2T)\) term appears because IV is total variance divided by time. Even if total variance rises smoothly, the volatility level can change sharply near short maturities because of the division by \(T\).
For the plots, we rescale these derivatives into units that correspond to smaller changes in moneyness and maturity:
slope for a 1% log-moneyness move: \(0.01\,\partial_k\sigma\),
curvature for a 1% squared move: \(0.01^2\,\partial_{kk}\sigma\),
term change for 30 calendar days: \((30/365.25)\,\partial_T\sigma\).
This scaling makes the derivative plots readable. Instead of showing abstract derivative units, we show approximate IV changes for practical moves.
15.1 Interpreting derivative units with a concrete move
The raw derivative \(\partial_k\sigma\) can be hard to read because \(k\) is a log variable. A one-unit move in \(k\) would be enormous. So we scale it to a 1% log-moneyness move.
If:
\[
\partial_k\sigma=-0.80
\]
then a 1% increase in log-moneyness gives approximately:
\[
\Delta\sigma\approx -0.80\times 0.01=-0.008
\]
That is a decrease of 0.8 vol points. If the current IV is 12%, it moves to about 11.2% for that local move.
For curvature, the scaled quantity is:
\[
\partial_{kk}\sigma\times(0.01)^2
\]
This is the curvature derivative scaled by the square of a 1% move. In a second-order Taylor approximation, its contribution to the IV change also has a factor of one half:
We add that contribution to the first-order slope term for a small move \(\Delta k\). The curvature can be small in IV units while still affecting gamma and the density calculations that depend on second derivatives.
For maturity, the scaled quantity:
\[
\partial_T\sigma\times\frac{30}{365.25}
\]
means “approximate IV change from extending maturity by 30 calendar days.” This is much easier to interpret than raw annualized derivative units.
Show code
k_plot_lo =max(limits_main["k_lo"], -0.22)k_plot_hi =min(limits_main["k_hi"], 0.08)k_shape_plot = np.linspace(k_plot_lo, k_plot_hi, 48)shape_marks = np.array([shape_k_down, shape_k_atm, shape_k_up], dtype=float)shape_marks = shape_marks[(shape_marks >= k_plot_lo) & (shape_marks <= k_plot_hi)]fig, ax = plt.subplots(figsize=(8.8, 4.9))for n, day inenumerate([30, 75, 150]): tau = day / ann_days y = iv_points(fit_visual, k_shape_plot, np.full_like(k_shape_plot, tau)) ax.plot(k_shape_plot, y, lw=2.0, color=lab_colors[n], label=f"{day}d")iflen(shape_marks): mark_y = iv_points(fit_visual, shape_marks, np.full_like(shape_marks, tau)) ax.scatter(shape_marks, mark_y, s=48, color=lab_colors[n], edgecolor="black", linewidth=0.6, zorder=4)for x in shape_marks: ax.axvline(x, color="black", lw=0.7, alpha=0.22)ax.set_title("smile shape points on supported log-moneyness")ax.set_xlabel("log strike / forward")ax.set_ylabel("implied vol")ax.legend(title="maturity")plt.tight_layout()plt.show()
We can now read the slope, curvature, and maturity behavior across the supported part of the fitted surface.
The slope plot is negative across most of the supported range, especially in the downside region. This means IV decreases as we move from downside strikes toward upside strikes. The front maturity has the strongest slope because the short-dated smile is steepest.
The curvature plot is positive in much of the downside/middle region and declines toward the right wing. Positive curvature means the slope increases as we move across strikes, so the smile bends upward even where the slope itself remains negative. A straight line would have zero curvature.
The contour maps add the maturity dimension. The slope map shows that the strongest negative slope sits near short-to-medium maturities and lower strikes. The curvature map shows where the surface bends most. The term-change map shows how IV changes when maturity increases by about one month. On this date, the term structure is not uniform: some strike regions gain more IV with maturity than others.
The smooth derivative fields give us the inputs needed for the local-volatility and Greek calculations. Jagged changes in slope or curvature would pass through to those calculations, particularly where the call-price second derivative becomes small.
16) Dupire Local Volatility Theory
Now we move from an implied-volatility surface to a local-volatility surface. The local-vol model assumes the underlying follows:
\(\sigma_{loc}(S_t,t)\) is volatility as a deterministic function of spot and time,
\(W_t\) is Brownian motion.
This is different from Black-Scholes. In Black-Scholes, volatility is constant. In local volatility, volatility changes depending on where the underlying is and how much time has passed.
Dupire’s formula gives us the local variance at each strike and maturity from a smooth, arbitrage-free surface of call prices \(C(K,T)\):
The derivatives are evaluated at strike \(K\) and expiry \(T\): the maturity derivative holds strike fixed, and the strike derivatives hold maturity fixed. We can read the numerator and denominator separately.
The numerator combines the change in call price with maturity and adjustments for carry. \(\partial C/\partial T\) measures the change as we extend maturity at fixed strike; \(qC\) and \((r-q)K\partial_KC\) account for the dividend and forward-drift terms at that maturity.
where \(f_Q(K,T)\) is the risk-neutral density of \(S_T\) at strike \(K\). A valid option surface must be convex in strike, so this second derivative should be nonnegative. If it becomes negative, the surface implies negative probability density, which is economically impossible.
Local variance is the ratio of the time-and-carry term to the density-related denominator. A small denominator can amplify errors in either part of the calculation, which is why smoothing and the supported region matter when we evaluate local volatility.
16.1 Deriving the intuition from the pricing PDE
The Dupire formula is easier to understand if we start from the idea that option prices must be consistent with a diffusion model. In a local-vol world, the risk-neutral dynamics are:
\[
dS_t=(r-q)S_tdt+\sigma_{loc}(S_t,t)S_tdW_t
\]
A derivative price \(V(S,t)\) satisfies the backward pricing PDE:
That is the usual option-pricing PDE written in spot-time variables. Dupire flips the problem. Instead of asking how the option price changes as current spot changes, we look at the whole surface of option prices across strikes and maturities. The call price \(C(K,T)\) contains information about the distribution of \(S_T\) under the risk-neutral measure.
The second strike derivative gives the density:
\[
\partial_{KK}C(K,T)=D(T)f_Q(K,T)
\]
The maturity derivative tells us how the distribution evolves through time. Combining those relationships gives the forward equation version of the pricing PDE, and solving it for \(\sigma_{loc}^2(K,T)\) gives Dupire.
The corresponding forward equation can be rearranged into a form for the change in the call surface with maturity:
The forward equation connects the evolution of the risk-neutral density with the curvature of call prices. Solving it for local variance gives the state-dependent volatility needed at each strike and maturity to reproduce the fitted cross-section of option prices.
A simple example makes Dupire easier to understand. Suppose downside puts are very expensive, so the left wing of implied volatility is steep. The market is saying: “if the index falls, volatility is likely higher.” A local-vol model translates that into a state-dependent volatility function where \(\sigma_{loc}(S,t)\) is higher at lower spot levels.
Local volatility can therefore sit above implied volatility in the downside region. Implied vol is the constant volatility that matches a particular option price under the pricing model, while local vol is the instantaneous diffusion coefficient assigned to a spot and time when we calibrate to the whole option surface.
The two volatility quantities have different units of interpretation:
Implied volatility is the constant volatility that reproduces a particular option price under the pricing formula.
Local volatility is a model-implied instantaneous volatility at a particular spot and time in the calibrated diffusion.
They are related, but they aren’t equal. A high local-vol-to-IV ratio means the surface requires strong state-dependent volatility in that region.
16.2 Local volatility and path dependence of the calibration
A local-vol model is calibrated so that it matches today’s vanilla option prices. That doesn’t mean it predicts future implied volatility dynamics correctly. It gives a diffusion process whose one-time marginal distributions match the current option surface.
Mathematically, vanilla prices tell us about the distribution of \(S_T\) for each maturity \(T\). The local-vol function is one way to build a Markov diffusion whose transition densities match those distributions.
This is why local vol is often called a static arbitrage calibration model. It matches today’s surface across strikes and maturities, but it doesn’t include stochastic volatility as a separate risk factor. If tomorrow’s entire implied-vol surface jumps up while spot barely moves, a pure local-vol model doesn’t naturally explain that as a new volatility state variable. It has to be recalibrated.
We use the calibrated local-volatility surface to study the state-dependent volatility implied by the observed option prices and to evaluate the behavior of the fitted surface:
What kind of state-dependent volatility is implied by today’s option surface if we force a local-vol diffusion to match it?
That answer is still informative, especially when we compare local vol to implied vol and measure where the ratio becomes extreme.
16.3 Reading Dupire through total variance intuition
We calculate Dupire from differentiated call prices, and we can also connect the formula to the total variance surface we fit. In the total-variance expression the maturity derivative is taken at fixed forward log-moneyness \(k\); in the price-space expression it is taken at fixed strike \(K\). In a simplified zero-rate, zero-dividend setting, the local variance can be written in terms of the total implied variance surface:
\[
w(k,T)=\sigma_{imp}^2(k,T)T
\]
A common form of the local variance expression is:
We don’t use this formula directly in the main implementation because rates, carry, and price differentiation are handled through the JAX call-price calculation. But the formula is very useful for intuition.
The numerator \(\partial_Tw\) says local variance needs total variance to increase with maturity. If total variance decreases with maturity, the numerator can become negative. That is one reason calendar arbitrage is dangerous for Dupire.
The denominator contains strike derivatives \(\partial_kw\) and \(\partial_{kk}w\). This tells us why skew and curvature matter. A steep skew changes the denominator. Strong curvature changes it too. If the denominator becomes very small, local variance explodes. If the denominator becomes negative, local variance becomes invalid.
This formula also explains why fitting \(w\) or \(\log w\) is better than fitting raw IV for a Dupire project. Local volatility is naturally connected to total variance dynamics. The fitted object, its maturity derivative, and its strike derivatives all enter the economic interpretation.
For example, if downside total variance rises sharply as \(k\) becomes more negative, then \(\partial_kw\) is strongly negative and \(\partial_{kk}w\) may be positive. The local-vol denominator changes, and the final local vol can become much higher in the downside region than implied vol alone suggests. The resulting local-volatility shape depends on both the slope and curvature of the fitted downside wing.
A small wiggle may barely affect the level of \(w\) while changing \(\partial_{kk}w\) substantially. Smoothing controls those high-frequency changes in the fitted surface so that the strike derivatives used in local volatility are more stable.
17) No-Arbitrage Diagnostics for Local Volatility
Before trusting local volatility, we check whether the fitted price surface behaves like a valid option surface.
Three conditions are especially important:
17.1 Total variance shouldn’t decrease with maturity
For a fixed moneyness region, total variance usually shouldn’t move downward sharply as maturity increases:
If total variance decreases with maturity, longer options can become cheaper in a way that allows calendar arbitrage.
17.2 Call prices should decrease with strike
A higher-strike call should be cheaper than a lower-strike call:
\[
\frac{\partial C}{\partial K}\le 0
\]
If this fails, one could buy the cheaper lower-strike call and sell the more expensive higher-strike call for arbitrage.
17.3 Call prices should be convex in strike
The second strike derivative should be nonnegative:
\[
\frac{\partial^2 C}{\partial K^2}\ge 0
\]
If this fails, the implied risk-neutral density is negative somewhere.
The quality output is encouraging. Inside the supported grid:
supported grid share is about 78.7%,
fitted IV finite share is 100%,
total variance maturity violation share is 0%,
call price strike monotonicity violation share is 0%,
negative density share is 0%,
total flagged share is about 4.52%.
The remaining flagged nodes mostly come from hard invalid local variance or extreme local-vol-to-IV ratios, not from negative density. That means the surface is mostly arbitrage-consistent, but Dupire still finds a few regions where derivatives become too aggressive.
17.4 Price derivatives and local-volatility diagnostics
A fitted IV surface can look smooth to the eye and still produce invalid Dupire values. We calculate local volatility from derivatives of option prices built from the fitted IV surface, so small changes in the fitted values can be amplified by those derivatives.
The conversion from IV to call price is nonlinear:
A small wave in \(\sigma(k,T)\) can become a larger wave in \(C_{KK}\), especially in low-vega or short-maturity regions. If \(C_{KK}\) becomes too small, the Dupire denominator becomes tiny. Dividing by a tiny number creates huge local vol. If \(C_{KK}\) becomes negative, the model implies negative density.
We flag grid nodes where the local-variance calculation, density, or support checks identify one of the following conditions:
local variance is negative,
density is negative,
the local-vol-to-IV ratio is extreme,
the node sits on a grid boundary,
the node is outside quote support.
The support mask is especially important. A local-vol value outside quote support is mostly an extrapolation artifact. We can plot it for visual continuity, but we shouldn’t interpret it as market information.
For the main SPX date, the supported local-volatility grid has:
median local-vol-to-IV ratio: about 1.22,
median absolute deviation from 1: about 0.231,
ATM local vol level: about 14.01%,
downside local-vol premium: about 24.61 vol points,
upside local-vol premium: about -5.15 vol points,
Dupire stress: about 0.1566.
The downside region has much higher local volatility than the ATM region, with a premium of about 24.61 vol points. That matches the fitted SPX shape, where the left-wing IV is also much higher than the central ATM level.
The local-vol surface contour shows high local volatility on the downside side of the grid and lower local volatility toward the upside. The slice plot shows that local vol and implied vol are not the same. Local vol can be above IV in the downside and deviate strongly at different maturities.
We can compare local volatility and implied volatility at each supported node through their ratio. A ratio near 1 means local vol is close to implied vol. A ratio far above 1 means the local-vol process needs stronger instantaneous volatility than the average IV suggests.
Under the calibrated local-vol model, the fitted skew corresponds to a state-dependent volatility process. The level at each node tells us the instantaneous volatility assigned to that spot and maturity by the model.
A local-volatility number by itself can be hard to interpret. Saying local vol is 25% at one node isn’t very meaningful unless we compare it with the implied vol at the same node, the ATM local vol, and nearby maturities.
The ratio tells us whether local vol is unusually high relative to the average volatility embedded in the option price. The premium tells us whether downside states require more instantaneous volatility than ATM states.
For SPX, the downside premium is the cleaner economic statistic. It says the local-vol process needs higher volatility when the index is lower. That is exactly the state-dependent volatility intuition behind equity skew.
The ratio also helps us inspect stability. Values around 1.2 or 1.3 in a well-supported region can describe a meaningful local-vol adjustment, while values of 5 or 10 in a sparse region can arise from a small or unstable derivative denominator. We flag extreme ratios when we evaluate the grid.
The flagged nodes concentrate in regions where the surface is most derivative-sensitive, especially near edges or short maturities. The rest of the supported interior gives us a clearer view of the local-volatility shape.
The invalid-by-maturity chart is also useful. Some maturities have small spikes in hard-invalid share, while most maturities are clean. This means the local-vol surface is broadly usable but shouldn’t be treated as perfectly reliable at every node.
The downside local-vol premium line stays positive across maturities. This means:
\[
\sigma_{loc}(k=-0.20,T)>\sigma_{loc}(k=0,T)
\]
for the plotted maturities. Economically, the local-vol model says future downside states have higher instantaneous volatility than ATM states. That is the local-vol translation of the market’s crash-risk pricing.
The calibrated local-vol model uses a deterministic function of spot and time to reproduce the option surface. Other market dynamics, including stochastic volatility, jumps, and volatility risk premia, can also contribute to the observed prices. We interpret the fitted map as the state-dependent diffusion implied by this calibration, evaluated from the option surface available on the date being studied.
19) Historical Surface Fitting Across All Dates
After understanding one date, we fit the visual and Dupire surfaces across all valid dates. This turns the project from a single-day calibration exercise into a time-series study of volatility-surface behavior.
The fit covers both the 2022 stress period and the lower-volatility end of 2023. We can therefore compare features through those different market conditions using the same fitting and diagnostic steps for each date.
The date-level fit summary shows average quote count around 4,749 per date, with a minimum around 3,598 and maximum around 5,714. This is dense enough to estimate a surface every trading day.
19.1 Research questions from the historical surface panel
A single-day local-vol surface is a calibration object. A historical panel turns the same machinery into a research dataset. Once we fit every date, we can ask questions like:
when is downside skew highest?
when does local-vol stress rise?
do curvature and Greek model risk move together?
is surface movement mostly level-driven or shape-driven?
how does the shape of the surface vary across high- and low-volatility dates?
With a fitted surface for each date, we can compare changes in level, skew, curvature, and maturity structure over the sample.
The historical panel also tells us whether the methodology is robust. If the fitter succeeds only on the main date, we don’t know much. If it succeeds on 505 consecutive SPX dates through 2022 and 2023, that is a stronger test. Those years include very different regimes: high inflation shock, aggressive rate hikes, equity drawdowns, and later the low-vol end-of-2023 rally.
All 505 valid dates were fitted without a skipped date, giving us a continuous set of daily surfaces under the chosen filters, basis dimensions, observation weights, and smoothing settings.
To compare surfaces through time, we need a common grid. A node is kept if it has enough support across dates. The common-node output keeps 675 nodes out of 1,225, or about 55.1% of the grid, using an 85% support threshold.
The support threshold keeps nodes that are observed across most dates, giving the historical averages and PCA a consistent set of locations. The kept log-moneyness range is roughly:
\[
-0.133 \le k \le 0.056
\]
This is narrower than the single-date surface because it requires common support across many days. That means we lose the far downside tail in the historical PCA, but we gain comparability.
The average-IV grid shows the usual SPX structure: higher vol on the downside and lower vol near the right side. The surface-volatility grid shows where the surface changes most through time. Variability is not uniform. Some regions of the smile move much more than others, especially the downside and medium maturities.
The common grid lets us compare the same surface locations on every date. Changes at different nodes can move together, so we can study both the variation at individual locations and the broader patterns shared across the fitted surfaces.
We summarize each daily surface using ATM levels, skew, curvature, term slope, and local-volatility diagnostics. Across the fitted dates, the means are:
mean short ATM IV: about 18.65%,
mean medium ATM IV: about 18.89%,
mean long ATM IV: about 19.39%,
mean term slope: about 0.74 vol points,
mean downside skew: about 7.41 vol points,
mean upside skew: about -4.52 vol points,
mean curvature: about 2.88 vol points,
mean smile asymmetry: about 11.93 vol points,
mean Dupire stress: about 0.0967.
The 5% to 95% range is more informative than the mean. Short ATM IV ranges from about 11.04% at calm points to about 27.80% in high-vol regimes. Downside skew ranges from about 4.20 to 10.63 vol points, so crash-protection premium never disappears, but it changes a lot.
The time-series plots show the regime story. Early 2022 has high and unstable short ATM IV, strong downside skew, and high curvature. Through 2023, ATM vol compresses, skew becomes smoother, and term structure changes as near-term uncertainty falls.
The term slope can turn negative during stress. A negative term slope means short-term IV is above longer-term IV. This usually happens when markets price immediate risk. In calmer periods, term slope becomes positive because longer maturities contain more accumulated uncertainty.
21.1 Economic meaning of the four main shape features
The historical features let us connect movements in the fitted surface to the pricing of different types of option exposure.
ATM IV is the general level of option-implied uncertainty. When it rises, the market is paying more for volatility across the center of the distribution. It usually rises during equity stress, macro uncertainty, liquidity shocks, or fast realized-volatility regimes.
Downside skew is the premium for left-tail protection. A high downside skew means downside options are expensive relative to ATM options. In SPX, this often reflects institutional hedging demand. Investors buy puts for crash protection, and dealers charge higher implied vol for supplying that convexity.
Curvature tells us whether the smile is strongly bowl-shaped. High curvature means both wings are far from the center, or the downside wing is extremely elevated. This is closely connected to gamma correction because curvature changes how quickly IV changes as spot moves.
Term slope tells us whether risk is front-loaded. A negative term slope means short-dated options are more expensive than longer-dated ones. That usually signals immediate event risk. A positive term slope means uncertainty is priced more gradually across time.
These four features turn a high-dimensional surface into a small set of interpretable market-state variables.
We compare medium ATM IV and downside skew across dates, with Dupire stress represented by color. Each point therefore locates one fitted surface by its overall volatility level, downside premium, and local-volatility diagnostics.
A point with high ATM IV and high downside skew means the market is pricing both broad uncertainty and strong crash asymmetry. A point with low ATM IV but still meaningful skew means the market is calm overall but still pays for left-tail insurance.
The rolling relationship between downside skew and Dupire stress changes through the sample. Sometimes they move together strongly, while at other times the relationship weakens or changes sign. Dupire stress also depends on curvature, local-vol-to-IV ratios, invalid nodes, and surface smoothness, so it can respond differently from the skew level.
For example, a surface can have high skew but be smooth and arbitrage-consistent, producing moderate Dupire stress. Another surface can have lower skew but more awkward curvature or maturity behavior, creating more local-vol stress.
Downside skew measures the relative level of the downside and ATM IV points. Dupire stress combines several local-volatility diagnostics, including invalid nodes and extreme local-vol-to-IV ratios, so it responds to more than the skew level alone.
23) PCA of Volatility Surface Changes
Now we analyze how the entire surface changes through time. If \(\sigma_t(k_i,T_j)\) is fitted IV at a common-grid node on date \(t\), the daily change at that node is:
Before fitting PCA, we center and standardize the daily changes at each retained grid node. PCA then decomposes the standardized vectors into orthogonal directions; the loading patterns can be mapped back to IV units for interpretation:
PCA gives us statistical directions of surface movement. We interpret those directions by mapping each component back onto the moneyness and maturity grid and examining which regions move together.
23.1 PCA of daily surface changes
We apply PCA to surface changes rather than raw surface levels because levels are highly persistent. If we used levels, the first component would mostly separate high-volatility dates from low-volatility dates. That is useful, but it doesn’t fully answer how the surface moves.
This is closer to a risk-management question. If we own an option book today, what kinds of surface moves might hit us tomorrow?
A broad upward move in IV is a level shock. A left-wing-only move is a skew shock. A front-end move is a term-structure shock. PCA identifies these directions statistically from the data instead of forcing us to define them manually.
The covariance matrix of flattened surface changes is:
The equation above gives the covariance of a set of surface-change vectors. In the implementation, each grid-node series is standardized before PCA, so we apply the same covariance calculation to the standardized vectors. The resulting eigenvectors describe directions in standardized space; we rescale their shocks back to IV units for the plots. A large eigenvalue means that direction explains a large share of the standardized historical variation.
The variance of the standardized surface changes is concentrated in the first few principal components:
PC1 explains about 92.55% of surface-change variance,
PC2 explains about 4.75%,
PC3 explains about 0.94%,
the first three explain about 98.25% in total.
The first component explains most of the standardized variation across the retained SPX grid nodes. Its loading pattern shows a broad level factor, with much of the surface moving in the same direction on a given day.
PC1 is the level shock. The PC1 shock map is mostly one sign across the surface, so a positive PC1 score corresponds to a broad IV increase. The top positive PC1 date is 2022-05-05, and high positive PC1 dates cluster in the volatile 2022 period. That fits the economic interpretation.
PC2 is more shape-like. It changes one region of the surface differently from another, often looking like a skew or term-structure rotation. PC3 is smaller and more localized, often closer to a term-structure or curvature adjustment.
The top-score table also gives a sanity check. PC1 positive extremes occur on dates with elevated short ATM IV, while PC1 negative extremes occur when the surface falls sharply. PCA is capturing real volatility regimes, not numerical noise.
fig, ax = plt.subplots(figsize=(10.2, 4.4))for i, arr inenumerate(pc_shocks[:3]):for day in [35, 75, 120]: idx =int(np.argmin(np.abs(tau_safe_grid * ann_days - day))) ax.plot(k_safe_grid, arr[idx], lw=1.6, label=f"pc{i +1}, {day}d")ax.axhline(0.0, color="black", lw=0.8)ax.set_title("pca shock cross-sections mapped back to iv units")set_supported_xlim(ax, k_safe_grid, common_node_mask)ax.set_xlabel("log strike / forward")ax.set_ylabel("approx iv change")ax.legend(ncol=3, fontsize=7)plt.tight_layout()plt.show()
We map the PCA shock cross-sections back into IV units. For each component, a one-standard-deviation score gives us an approximate IV change at the selected strikes and maturities.
Across the first three components:
PC1 mostly moves all strikes in the same direction.
PC2 changes the relative behavior of downside versus upside or short versus longer maturities.
PC3 is smaller and more localized.
This is similar to yield-curve PCA, where level, slope, and curvature often dominate. For volatility surfaces, the language is slightly different: level, skew/rotation, and curvature/term-structure modes.
The episode comparison uses dates selected from extreme PC1, high downside skew, high curvature, and high Dupire stress. The 60-day smile comparison shows how different episodes reshape the same surface region. The ATM term-structure comparison shows whether stress is concentrated in the front end or spread across maturities.
We can use the resulting factors to follow the broad level, skew, and term-structure movements of the fitted surfaces over time, while still inspecting their effects at individual grid nodes.
When spot \(S\) moves, the option price changes through the forward input and through the IV assigned at the new moneyness. Holding the fitted surface fixed and using \(\partial k/\partial S=-1/S\), the chain rule gives:
The first expression maps smile slope into the IV change from a spot move; the second combines smile curvature and slope for the spot second derivative.
\(\sigma_S\) and \(\sigma_{SS}\) come from the surface slope and curvature.
This formula shows why surface-aware gamma can differ a lot from flat-vol gamma. If skew is steep, \(\sigma_S\) is meaningful. If curvature is strong, \(\sigma_{SS}\) is meaningful. We use JAX automatic differentiation to calculate the full spot derivatives directly from the fitted surface and pricing function.
24.1 Connecting surface-aware Greeks to vanna and volga
The surface-aware gamma formula contains terms that option traders usually think of separately:
These terms carry the fitted smile’s slope and curvature into the spot sensitivity of the option price. The factor of two on vanna comes from the two spot-volatility cross terms in the second derivative.
If skew is steep, \(\sigma_S\) is large. If smile curvature is strong, \(\sigma_{SS}\) is large. Then vanna, volga, and vega terms become important. If the surface is flat, \(\sigma_S\approx 0\) and \(\sigma_{SS}\approx 0\), so surface-aware Greeks collapse back toward flat-vol Greeks.
We can compare the resulting corrections with the fitted slope and curvature to see where the difference between flat-vol and surface-aware Greeks is largest.
Suppose we hold a put with a fixed strike and SPX falls while the fitted surface stays fixed in forward-moneyness coordinates. The forward falls with spot, so \(k=\log(K/F)\) rises and the option moves across the fitted smile. Where the smile has negative slope in \(k\), the IV read from that fixed surface decreases as spot falls. The put also has a direct spot sensitivity, and the change in the assigned IV adds a vega-related adjustment to it. This calculation describes movement along the fitted surface; it does not assume that the entire market surface stays unchanged in the future.
Using the fitted skew, we include the change in the option’s assigned IV as it moves across the surface. The resulting delta and gamma can be compared with the corresponding flat-vol Greeks for the same starting option and spot level.
For the main SPX date, the surface-aware and flat-vol Greek comparison gives us:
delta RMS difference about 0.0612,
max absolute delta difference about 0.1061,
gamma RMS difference about 0.000242,
total Greek PnL RMS about 2.93 index points for a 1% spot shock,
relative Greek PnL RMS about 0.78% of the option price reference.
A delta difference of 0.06 corresponds to roughly six units of delta across 100 option units under the same scaling. The difference comes from the IV change along the fitted skew when spot moves, which the flat-vol calculation holds fixed.
The delta correction plot is positive across a large part of the downside region. This means surface-aware delta is higher than flat-vol delta in those regions. The gamma correction changes sign across moneyness and maturity, reflecting the interaction between skew, curvature, and option convexity.
The largest corrections are concentrated where surface slope and curvature are strongest. Their locations across moneyness and maturity let us connect the plotted Greek differences to the vega, vanna, and curvature terms in the surface-aware formulas.
We compute the same Greek-risk measures across all 505 fitted dates. The full run succeeds on all dates using JAX.
The first rows in 2022 show much higher Greek model risk than the latest main date. Total Greek PnL RMS is around 5.8 to 6.2 index points in early January 2022, compared with about 2.93 on 2023-12-29. This matches the broader regime: 2022 had steeper skew, higher IV, and more unstable surface shape.
Across the fitted dates, the Greek-risk measures have the following correlations with surface features:
total Greek PnL RMS correlation with downside skew: about 0.908,
gamma PnL RMS correlation with downside skew: about 0.925,
gamma PnL RMS correlation with curvature: about 0.979,
gamma PnL RMS correlation with short ATM IV: about -0.707.
The gamma-curvature relationship is almost direct. This makes sense mathematically. Gamma correction includes surface curvature terms, so when the smile bends more sharply, the flat-vol gamma becomes less reliable.
The high positive relationship with downside skew also makes sense. Stronger left skew means the option surface moves more aggressively as spot changes, so surface-aware delta and gamma diverge more from flat-vol Greeks.
The scatter plot between downside skew and total Greek PnL RMS shows a strong upward relationship. Higher skew means larger surface-aware Greek risk. The color dimension, Dupire stress, adds another layer: local-vol stress and Greek model risk are related but not identical.
25.1 How the Greek-risk PnL scaling works
The Greek-risk metrics translate differences in delta and gamma into approximate PnL units under a 1% spot shock. If the spot shock is:
\[
\Delta S = 0.01S
\]
then the delta correction PnL is approximately:
\[
\Delta \text{PnL}_{\Delta}=\big(\Delta_{surf}-\Delta_{flat}\big)\Delta S
\]
RMS is useful because positive and negative corrections don’t cancel. If some options have positive correction and others have negative correction, the average could look small even though the risk is large. RMS measures typical magnitude.
The relative Greek PnL RMS divides by a reference option price. That makes the metric comparable across regimes where option prices are very different. A 5-index-point correction is more serious for a cheap option than for an expensive one.
26) Secondary Application: BTC Option Surface Through the Library
We now repeat the surface workflow on BTC options using the packaged library. BTC has a different volatility level, quote structure, and economic interpretation of carry, so we can apply the same calculations to a substantially different option market.
For reproducibility, the BTC option data comes from the btc_optionsdx folder, BTC spot data from the btc_yfinance folder, and rates from the us_treasury_yields folder inside the repo’s data/ directory. Each folder explains how to rebuild the local dataset from the source/raw files.
The BTC market differs from SPX in several ways:
BTC trades every day, so calendar time is more natural than trading-day time.
Deribit-style option chains can have very wide strikes and extreme moneyness.
Volatility levels are much higher than SPX.
The surface can be rougher because crypto options have different liquidity patterns.
Carry is not a dividend yield in the equity-index sense. It comes from parity, funding, and crypto market structure.
The main BTC date is 2024-03-22. We have:
120,736 surface-ready quotes,
1,017 quotes on the main date,
6 expiries on the main date,
JAX used for Dupire and Greeks with no fallback.
The dataset is smaller and rougher per date than SPX. That matters for surface fitting because the maturity dimension is much thinner.
26.1 BTC option-market characteristics
BTC options create a very different surface problem from SPX. SPX skew is heavily shaped by institutional crash protection, index-option supply/demand, and equity risk premia. BTC options reflect crypto leverage, speculative positioning, funding conditions, event risk, and a market that trades continuously.
This changes the interpretation of the same math. For SPX, left-wing skew usually means downside crash insurance. For BTC, both upside and downside tails can be important because crypto can gap violently in either direction. A call wing can become expensive during a speculative rally, while a put wing can become expensive during liquidation risk.
BTC option demand can also produce a different shape across strikes, so the same moneyness coordinate needs to be interpreted in the context of crypto returns and option positioning. The BTC chain also has fewer expiries on the main date, only 6. That makes the maturity dimension weaker. With fewer maturities, a 2D surface fit has less information about \(\partial_T\). Since Dupire depends on \(\partial_TC\), BTC local vol can be more fragile even if strike slices look reasonable.
The BTC surface therefore gives us another setting for the full workflow, with:
higher volatility level,
wider strike range,
rougher quote structure,
fewer maturity points,
different carry interpretation,
stronger jump and event risk.
26.2 Carry and rates in BTC options
For BTC, the role of carry is different from SPX. SPX carry is tied to interest rates and dividends. BTC doesn’t pay dividends, but crypto option parity can still imply a forward different from spot. That difference can reflect funding, basis, collateral conventions, liquidity, and market demand.
We can still use the same formal relationship:
\[
F_T=S_0e^{b(T)T}
\]
but the economic interpretation of \(b(T)\) changes. In SPX, \(q(T)=r(T)-b(T)\) can be read roughly as an implied dividend yield. In BTC, reading it as a dividend yield would be wrong. It is better to call it implied carry.
We use the same forward-and-discount pricing framework for BTC, and interpret the parity-implied carry using the crypto market’s funding, basis, collateral, and liquidity conventions. Each expiry has its own relationship between the spot reference and the forward inferred from option parity. The rate and forward inputs retain their pricing roles for each expiry.
The BTC date uses far wider moneyness support than SPX. That makes the surface visually more dramatic, but it also means extrapolation risk is higher. A log-moneyness value near \(-0.79\) is an extremely deep downside region. Quotes there can be informative about disaster pricing, but they can also be sparse and noisy.
The BTC cleaning report shows how much filtering is needed:
raw rows: 258,234,
after positive price filter: 218,710,
after bid/ask/spread filter: 196,404,
after DTE filter: 136,262,
after moneyness filter: 121,361,
final rows: 121,361.
The largest removal happens at the DTE filter, which removes over 60,000 rows. That means a large part of the BTC raw chain sits outside the maturity window we want. Another 14,901 rows are removed by the moneyness filter, which is expected because crypto option chains often include very far strikes.
The BTC fit summary is much worse than SPX in absolute IV-error terms:
visual RMSE about 0.0396,
Dupire RMSE about 0.0542,
visual weighted RMSE about 0.0144,
Dupire weighted RMSE about 0.0251.
BTC implied volatilities are much higher and the fitted surface covers more extreme strikes than the SPX example. The main-date log-moneyness range runs from about -0.79 to 0.54, much wider than the SPX main surface. A few deep strikes produce residuals of 20 to 50 vol points.
The pchip-minus-spline comparison has RMSE about 0.0423, which confirms the smooth spline differs meaningfully from raw interpolation. The difference is visible in the deep wings, where the smoothing penalty reduces the influence of local quote-to-quote fluctuations on the fitted surface.
27) BTC Local Volatility, PCA, and Surface-Aware Greeks
For the main BTC date, the local-volatility diagnostics are:
hard invalid share: 0%,
extreme ratio share: 0%,
negative density share: 0%,
median local-vol-to-IV ratio: about 1.015,
median absolute local-vol-to-IV deviation: about 0.029,
ATM local vol: about 77.23%,
Dupire stress: about 0.0145.
This is very different from the SPX main date. SPX has strong downside skew and local vol sits far above IV in downside states. BTC on this date has very high volatility level, but local vol is close to implied vol across much of the supported region. High volatility level doesn’t automatically mean high Dupire stress. Stress comes from the shape and derivatives of the surface, not only from the level.
BTC PCA is also concentrated:
PC1 explains about 90.09%,
PC2 explains about 5.78%,
PC3 explains about 2.12%,
the first three explain about 97.99%.
So even in BTC, surface changes are dominated by a broad level factor, with smaller shape factors after that. The difference is that the scores are measured in a different volatility scale. BTC IV changes are large in absolute vol points, but the factor structure is still low-dimensional.
The highest BTC Greek-risk dates have total Greek PnL RMS around 11.5 to 12.9 for a 1% spot shock. The number is larger than SPX partly because BTC spot is high and BTC option vol levels are much higher. A 1% spot move scales with the underlying price, so the relative Greek PnL measure also helps us compare the correction with a reference option price. We can also compare where the corrections are concentrated: the delta and gamma maps show their largest effects in regions where the fitted smile bends or moves most aggressively with spot.
The final BTC dashboard puts the full library pipeline in one view: quote coverage, smooth surface, local vol, smile comparisons, residuals, local-vol ratio, PCA variance, PCA shocks, and Greek corrections. This confirms that the same machinery works outside SPX, but the interpretation changes. For SPX, the dominant story is downside crash skew and equity-index local-vol premium. For BTC, the dominant story is extremely high volatility with a rougher but sometimes locally cleaner surface.
27.1 BTC surface shape and risk
On the main BTC date, ATM local volatility is around 77% and Dupire stress is low. The volatility level can be high while the fitted surface remains smooth enough to produce stable derivatives in the supported region. The stress measure combines local-volatility diagnostics and deviations from implied volatility, rather than using the level alone.
The SPX main date has lower IV level but stronger local-vol-to-IV distortion because the downside skew is very steep. BTC has much higher IV level, but on this selected date the local-vol ratio is close to one in the supported region. That means the surface shape is not forcing strong state-dependent volatility beyond the already high implied-vol level.
The residual table still shows very large errors in extreme BTC strikes. Deep BTC options can be very hard to fit because their quotes can be wide and the market can price tail scenarios aggressively. The smoother surface chooses stability over exact far-tail quote matching.
The PCA result also tells us that BTC surface changes are mostly level-driven, with PC1 above 90%. That is similar to SPX, but the economic meaning is different. A BTC level shock can reflect broad crypto risk repricing, leveraged liquidation pressure, ETF/news events, or macro liquidity changes. It is not the same as SPX crash insurance, even if the PCA math is the same.
In the final dashboard, we can compare the fitted surface and its diagnostics in one place. Quote coverage tells us where data exists. The 3D smooth surface shows the fitted IV object. Local vol tells us the implied state-dependent volatility. Residuals show where the fit misses. PCA shows how the surface moves historically. Delta/gamma corrections show how much flat-vol Greeks miss once skew is allowed to move with spot.
27.2 SPX versus BTC: same machinery, different market story
We can compare the SPX and BTC applications across their volatility levels, surface shapes, support, and factor behavior.
Dimension
SPX surface
BTC surface
Main economic story
downside crash insurance
high crypto volatility and wide tails
Typical IV level
lower
much higher
Maturity support
many expiries on main date
fewer expiries on main date
Skew interpretation
institutional put demand, equity crash premium
two-sided crypto tail risk and funding/liquidity effects
Local-vol stress
driven by steep downside skew
low on selected date despite high vol level
PCA structure
PC1 level dominates
PC1 level also dominates
The shared result is that both markets have low-dimensional surface movements. PC1 dominates in both SPX and BTC. The difference is what that level shock means. In SPX, a level shock often reflects equity-risk repricing and hedging demand. In BTC, it can reflect crypto-wide leverage, liquidation risk, ETF/news events, and broader risk appetite.
We use the same PCA and Dupire concepts for both applications, then interpret their fitted factors and local-volatility shapes in the context of each market. The underlying exposures and option-market structure give the same mathematical quantities different financial meanings.
