A clean SPX option modeling panel, reusing the option-chain, forward, IV, and Greek machinery from the earlier option projects without re-teaching those pieces.
A sequence of smile and stochastic-volatility models: SVI, SSVI, SABR, Merton jump-diffusion, Heston, and Bates.
A model comparison layer that separates fast surface fit, global surface discipline, stochastic dynamics, jump-tail behavior, runtime, Monte Carlo noise, and parameter warnings.
A model-risk-aware fair-value engine that turns fitted surfaces into residual scores only when pricing error, bid-ask uncertainty, and model disagreement leave enough room.
A next-day validation test that checks whether residuals have any forward information after hedging and costs.
A scheduled delta-hedged trading experiment comparing model-selected residual trades against matched ATM-style baseline trades.
A secondary BTC implementation using the library, where the same model stack faces a market with much wider smiles, rougher liquidity, and 24/7 trading.
We start with different ways of explaining the same implied-volatility surface.
A surface model says, “give me today’s market quotes and I’ll draw a smooth smile.”
A stochastic volatility model says, “give me a stochastic process for the underlying and its variance, and I’ll generate option prices from the distribution of future outcomes.”
A jump model adds the idea that part of the smile comes from discontinuous moves that a continuous diffusion can’t reproduce efficiently.
We already spent earlier option projects on Black-Scholes/Black-76, implied-volatility inversion, Greeks, and surface/local-vol construction. Here we treat those as tools. We now build the model layer above the surface and study how each model represents skew, convexity, tails, uncertainty, and tradability.
The goal isn’t only to find the smallest calibration error. We want to know which model errors are economically meaningful, which are just quote noise, which are too expensive to trade, and which survive a hedged next-day test.
The comparison brings the fitted surface, model dynamics, and hedged outcomes into the same workflow.
Surface Shape, Volatility Dynamics, and Their Use in Trading
We can organize volatility modeling into three connected layers: the quotes we observe, the values implied by a model, and the validation of a possible trade.
Layer A: quoted implied volatility
The market gives us option prices. After cleaning and IV inversion, we can write the quote as a point on a surface:
\[
(k,T,\sigma_{imp})
\]
where \(k=\log(K/F)\) is log-forward moneyness and \(T\) is time to expiry. This is the descriptive layer. It tells us what the market is charging today.
and we compare that against the quote. This is the calibration layer. SVI, SSVI, and SABR mostly live here as fast surface/smile explanations. Merton, Heston, and Bates also live here, but they generate prices through a distributional story.
Layer C: trade validation
A residual becomes useful only if it has forward information after costs:
\[
e_i=V_i^{model}-V_i^{mkt}
\]
If \(e_i\) is positive, the model says the option is cheap. If \(e_i\) is negative, it says the option is rich. But that residual still has to survive bid-ask, model disagreement, hedge cost, and next-day validation.
After fitting the models, we can also test their residuals in a scheduled hedged experiment. A close fit tells us about today’s surface; the subsequent payoff tells us whether the residual ranking contains useful information after trading costs.
Each layer has a different output: observed quotes, fitted prices, and subsequent trade outcomes.
Model Families in One Map
The six models can be organized by what they believe the source of the smile is.
Family
Models here
Main belief about smile
Static smile geometry
SVI, SSVI
The surface shape can be parameterized directly in \(k\) and \(T\)
Approximate stochastic smile
SABR
Volatility is random and correlated with the forward, but we use a fast IV approximation
Jump distribution
Merton
Tail prices come from rare discontinuous moves
Continuous stochastic variance
Heston
Tail prices and skew come from random mean-reverting variance
Hybrid jump + variance
Bates
Both stochastic variance and discontinuous jumps are needed
With these model families in mind, we can compare fit and dynamics on their own terms. SVI can beat Heston on RMSE because it is designed to fit the observed smile directly. That doesn’t mean SVI gives a better scenario generator. Heston can have worse fit but a richer path story. Bates can improve the stochastic fit but become harder to identify.
So the practical question becomes: which role do we need?
For daily fair-value residuals, we need tight, stable cross-sectional pricing.
For understanding tail mechanisms, we need jump and stochastic-volatility diagnostics.
For avoiding false trades, we need model disagreement and uncertainty filters.
For backtesting, we need a hedged P&L workflow that doesn’t assume calibration error is automatically tradable.
We will use each family for the role its assumptions and fitted behavior support.
Surface Fit Versus Path Consistency
We can model the same option quotes with a cross-sectional surface or with a stochastic process for the underlying.
A surface model is cross-sectional. It asks:
\[
(k,T)\mapsto \sigma_{imp}(k,T)
\]
The fitted object is a map across strikes and maturities at one date. It can be extremely accurate because it is built to match the observed surface.
A path model is distributional. It asks for a stochastic law:
\[
F_0 \rightarrow F_T
\]
or, with stochastic variance:
\[
(F_0,v_0)\rightarrow (F_T,v_T)
\]
The option price then comes from the implied terminal payoff distribution:
\[
V=D\,E^Q[g(F_T)]
\]
Path consistency is valuable because the same process can price many payoffs, simulate scenarios, compute path-dependent risks, and connect smile shape to forward dynamics. But a low-dimensional path model may fail to match every quote on one surface.
We therefore compare calibration RMSE alongside the assumptions behind each model. A static model can fit today’s surface beautifully but say little about tomorrow’s smile dynamics. A path model can teach us about skew, jumps, and variance clustering while still losing the fit contest.
1) Reused Data and the New Modeling Target
We use the SPX option-chain setting from the earlier option projects. The data sources and their preparation scripts are documented in these folders:
The option-chain cleaning, put-call parity forward inference, Black-76 implied-volatility inversion, and Greek computation have already been covered before. We only remind the minimum structure needed for the models in this project.
At every quote we need:
option type: call or put,
strike \(K\),
expiry time \(T\),
discount factor \(D(T)\),
forward price \(F_t(T)\),
mid, bid, ask, and bid-ask spread,
implied volatility \(\sigma_{imp}\),
log-forward moneyness \(k=\log(K/F_t)\),
Greek-based sensitivity, especially vega.
The pricing object underneath most of the project is the forward-discounted risk-neutral expectation:
\[
C_t(K,T)=D(t,T)\,E_t^Q[(F_T-K)^+]
\]
and for puts:
\[
P_t(K,T)=D(t,T)\,E_t^Q[(K-F_T)^+]
\]
Here \(E_t^Q[\cdot]\) means expectation under the risk-neutral pricing measure, not under the real-world return distribution. The models are calibrated to option prices, while historical SPX returns are observed under the real-world measure. If the market charges a volatility risk premium or jump risk premium, the calibrated parameters reflect the pricing distribution the market implies.
The discount factor converts a payoff at expiry into its value today. For example, \(D(t,T)\) is the price at \(t\) of one unit paid at \(T\), while \(F_t(T)\) is the forward for the same expiry. Once those inputs are fixed, all six models can be compared on the same strike and maturity rather than each model choosing its own carry assumption.
The expectation is over possible terminal outcomes under \(Q\). A put gives a positive payoff only when the terminal forward falls below \(K\); a call pays only above \(K\). A change in the distribution’s tails therefore changes option values even when the expected forward remains the same.
Risk-Neutral Model Parameters Versus Real-World Volatility
A common beginner mistake is to read calibrated option-model parameters as direct historical forecasts. We’ll avoid that. In this project, parameters are fitted to option prices, so they live under the risk-neutral measure \(Q\).
The real-world return process is usually written under \(P\):
\[
dS_t = \mu_P S_t dt + \sigma_t S_t dW_t^P
\]
Under risk-neutral pricing, the drift changes so discounted prices become martingales:
\[
dS_t = rS_tdt+\sigma_t S_t dW_t^Q
\]
For forwards, the drift is usually removed and we price with:
\[
\frac{dF_t}{F_t} = \sigma_t dW_t^Q
\]
The volatility process can also differ under \(P\) and \(Q\). For Heston, the variance mean-reversion level and speed can contain volatility risk premia. A fitted \(\theta\) is therefore the long-run variance level implied by option prices, not necessarily the long-run realized variance we’d estimate from SPX returns.
We use the fitted models for pricing and relative value, keeping option-implied and realized future volatility separate. When a fitted Bates jump intensity is high, the clean interpretation is: the option surface is pricing jump risk strongly. It doesn’t automatically mean the real-world probability of a jump is exactly the fitted \(\lambda\).
The change of measure changes probabilities used to value payoffs. It doesn’t mean that the realized index has risk-free drift. In a simple setting without dividends, discounting the spot by the money-market account gives the martingale condition \(E_t^Q[e^{-r(T-t)}S_T]=S_t\); the notebook works with forwards and discount factors so that carry can be handled separately.
For stochastic variance, the parameters of the variance equation under \(Q\) describe the variance risk priced by options. We can compare these parameters across dates as model-implied quantities, while historical estimates under \(P\) answer a different statistical question.
Risk Premia and Calibrated Parameters
Calibrated volatility models often produce parameters that look strange if we read them as historical statistics. That happens because option prices include risk premia.
For example, equity-index options often price downside risk more heavily than the realized frequency of crashes alone might suggest. Investors are willing to pay for crash protection, so the risk-neutral left tail can be heavier than the physical left tail. In model terms, that can show up as:
more negative \(\rho\) in Heston,
higher \(\xi\) in Heston,
negative jump means in Merton/Bates,
high jump intensity under \(Q\),
steep SVI left-wing slope.
The pricing measure doesn’t need to match historical frequencies exactly. It needs to reproduce the market prices of payoffs. This is why a model can imply a large crash premium even in a calm realized-volatility period.
We can also use historical outcomes to check whether the resulting residuals contain useful information. If the model’s residual mostly reflects a risk premium that the market intentionally charges, fading it can lose money. A residual becomes useful only if it predicts subsequent repricing or hedged P&L, not simply because it exists.
We start by loading the option pairs, Treasury rates, pricing utilities, and plotting style. The models need a consistent quote panel: duplicated, illiquid, or unusually wide quotes can make fitted parameters respond to microstructure noise instead of the volatility patterns in the market.
We attach a forward and discount factor to each quote, giving the later model prices a consistent set of pricing inputs.
2) Paired Calls, Rates, Forwards, and the Pricing Coordinate System
We start with paired calls and puts. Put-call parity gives us the forward needed to put the different expiries on a common pricing scale.
For European-style index options, the parity relationship in discounted form is:
\[
C_t(K,T)-P_t(K,T)=D(t,T)\,[F_t(T)-K]
\]
Solving for the forward gives:
\[
F_t(T)=K+\frac{C_t(K,T)-P_t(K,T)}{D(t,T)}
\]
We do this across many paired strikes for the same date and expiry, then use a robust estimate rather than trusting one strike. The reason is simple: a single quote can be stale or wide. If 150 paired quotes imply almost the same forward but a few bad pairs are outliers, the robust estimate keeps the forward anchored to the liquid part of the chain.
Once we have \(F_t(T)\), the natural moneyness coordinate is:
\[
k=\log\!\left(\frac{K}{F_t(T)}\right)
\]
This coordinate is better than raw strike because it puts every expiry on a comparable scale:
\(k=0\) means at-the-money forward,
\(k<0\) means downside puts / low strikes,
\(k>0\) means upside calls / high strikes.
For example, a 5% OTM put and a 5% OTM call are roughly symmetric around \(k=0\), while raw strikes aren’t comparable when the index level changes through time.
Forward Moneyness and Model Stability
We express the smiles in log-forward moneyness \(k\) so that the strike coordinate can be compared across expiries and index levels. The forward price already contains the impact of rates, dividends/carry, and maturity. Without that transformation, the same economic strike location would shift mechanically when rates or carry change.
The transformation is:
\[
k=\log(K/F)
\]
A few reference points:
\(k=0\) means \(K=F\).
\(k=-0.10\) means \(K\approx 0.905F\), a roughly 9.5% downside strike.
\(k=0.10\) means \(K\approx 1.105F\), a roughly 10.5% upside strike.
Because the logarithm is symmetric in return space, a downside move and upside move become more comparable. A 10% downside strike has:
\[
\log(0.90)\approx -0.105
\]
and a 10% upside strike has:
\[
\log(1.10)\approx 0.095
\]
That is why surfaces are plotted and fitted in \(k\). The shape near \(k=-0.20\) always means “deep downside protection,” regardless of whether SPX is at 3,800 or 5,000.
Strike \(K\) and the forward \(F\) have the same price units, making their ratio dimensionless.
The first paired quotes contain deep ITM calls and very cheap OTM puts. That’s normal for low strikes when spot is around \(4{,}795\). The relative spread is already telling us something useful: the call side is tight because the option is large intrinsic value, while the matching tiny put can have a much larger relative spread. Later, this is exactly why we don’t treat all quotes equally.
Loading the rate curve takes about 138 seconds in the executed run. That timing isn’t the economic result, but it tells us this project is heavy from the start: large SPX option panels plus curve mapping plus IV solving create millions of pricing rows.
The quote size and relative spread can therefore tell different stories for the two legs.
On the first SPX date, the inferred forwards are just below spot for the short expiries. The implied carry values are slightly negative. That can come from index dividends, rate assumptions, and parity noise. Several maturities have more than 150 paired quotes per expiry. That gives us enough cross-sectional support to infer forwards rather than imposing them mechanically.
The parity error IQR is below a few index points for many expiries. For SPX, where the index level is near \(4{,}800\), that is small enough for surface modeling. If parity errors were huge, the model layer would be contaminated because the same strike would be assigned the wrong \(k\) coordinate.
Forward estimates are grouped by quote date and expiry, so each option is assigned the forward for its own maturity.
3) From Paired Quotes to a Long Option Panel
After the forward is inferred, we transform paired call-put rows into a long option table with one row per option quote. This is the format needed for IV solving, Greeks, and model fitting.
The implied volatility is the value of \(\sigma\) that solves:
\[
V^{B76}(F,K,T,D,\sigma)=V^{mkt}
\]
where \(V^{B76}\) is the Black-76 price and \(V^{mkt}\) is the option mid, bid, or ask. We don’t re-derive the Black-76 solver here because that was already done in the earlier option project. What matters here is how the solved IV becomes the modeling target.
We will move between three spaces throughout the model comparison:
Space
Object
Why we use it
Price space
\(V^{mkt}\) and \(V^{model}\)
Needed for real P&L, bid-ask costs, and Monte Carlo pricing models
Implied-vol space
\(\sigma_{imp}\)
Makes model errors comparable across strikes and maturities
Total-variance space
\(w(k,T)=\sigma_{imp}^2(k,T)T\)
Usually smoother and safer for surface parameterizations
The fitted target depends on the model. SVI and SSVI fit total variance. SABR produces implied volatility. Merton, Heston, and Bates naturally produce prices, then we invert those prices back into IV for diagnostics.
Price Error, IV Error, and Total-Variance Error
We use three error scales because each one answers a different question.
A price error is:
\[
e_i^V=V_i^{model}-V_i^{mkt}
\]
This is closest to actual tradable money. If the model is $2 above the market mid, that can look like a cheap option. But price errors are hard to compare across strikes. A $2 error on a $300 option is small. A $2 error on a $1 option is huge.
An implied-vol error is:
\[
e_i^\sigma=\sigma_i^{model}-\sigma_i^{mkt}
\]
This makes errors more comparable across strikes and maturities. A 50 vol-bp error means something interpretable. But IV error can exaggerate very low-vega regions because tiny price changes can imply large IV changes.
A total-variance error is:
\[
e_i^w=w_i^{model}-w_i^{mkt}
\]
with:
\[
w_i=\sigma_i^2T_i
\]
Total variance is useful for surface parameterizations because variance accumulates with time. If the one-month vol is \(20\%\), the total variance is roughly:
\[
0.20^2 \times \frac{1}{12}=0.0033
\]
If the one-year vol is \(20\%\), the total variance is:
\[
0.20^2 \times 1=0.04
\]
The volatility number is the same, but the variance budget is very different. SVI and SSVI work in this accumulated-variance space because it behaves more smoothly across maturities.
The transformations can also be read locally. At a fixed maturity, \(dw/d\sigma=2\sigma T\), so a small volatility change \(\Delta\sigma\) gives \(\Delta w\approx2\sigma T\Delta\sigma\). The same IV error therefore corresponds to different total-variance errors for different maturities. Dividing by an uncertainty scale in the same units as the target lets the calibration compare those observations consistently.
Error Translation Through Vega
The link between price error and IV error is approximately vega:
Comparing models in both price and IV space keeps those different sensitivities visible. A model can have a small price error on a low-vega option but a large IV error. It can also have a large dollar error on a high-vega option but a small IV error.
For example, a $1 price error on an option with vega \(500\) is:
\[
\frac{1}{500}=0.002
\]
or about 20 vol basis points. The same $1 error on an option with vega \(20\) is:
\[
\frac{1}{20}=0.05
\]
or 5 vol points. Those two errors have completely different meaning.
This is also why the project never treats deep wings casually. Deep OTM options can have low vega, wide spreads, and jump-sensitive prices. They are exactly where stochastic models are interesting, but also where diagnostics have to be strict.
Vega here is the derivative with respect to volatility expressed as a decimal: an IV move from \(0.20\) to \(0.21\) is \(0.01\), or one volatility point. The approximation comes from the first-order Taylor expansion \(V(\sigma+\Delta\sigma)\approx V(\sigma)+\mathcal{V}\Delta\sigma\). It is local; a large move in IV also depends on the curvature of price with respect to volatility.
Calls and puts now appear as separate rows at the same strike and expiry. Very deep OTM options can have tiny prices and wide relative spreads, so a valid row in the table isn’t automatically a good modeling quote. We still need filtering, OTM selection, vega checks, and bid-ask-aware weights.
The IV solving step is large: more than 4.0 million rows are processed, with about 3.99 million valid mid IV quotes. The engine is numba, and the solve takes about 50.8 seconds in the executed run. Processing this panel efficiently makes the later model comparisons practical across many dates.
Bid and ask IVs are solved alongside the mid IV and will be used to estimate quote uncertainty.
4) Greeks, Data Quality, and the Model-Ready Quote Set
After IV solving, we compute Greeks because they enter two later tasks:
uncertainty scaling, because a price error translates into an IV error through vega,
delta-hedged backtesting, because the trading experiment needs hedge ratios.
The basic local approximation is:
\[
dV \approx \mathcal{V}\,d\sigma
\]
where \(\mathcal{V}\) is vega. Rearranging gives the rough IV uncertainty implied by a price uncertainty:
\[
d\sigma \approx \frac{dV}{\mathcal{V}}
\]
This is why a $0.50 bid-ask half-spread isn’t equally meaningful everywhere. If vega is very high, that half-spread corresponds to a small IV uncertainty. If vega is almost zero, the same half-spread can mean the implied volatility is poorly pinned down.
A simple example makes this clear. If an option has vega \(200\) and the half-spread is \(0.50\), the spread-implied IV uncertainty is about \(0.0025\), or 25 volatility basis points:
\[
\frac{0.50}{200}=0.0025
\]
If another option has vega \(5\) with the same half-spread, the implied IV uncertainty becomes \(0.10\), which is 10 vol points. The second quote is much less useful for calibration even if the dollar spread looks the same.
That is a strong option panel: wide enough in date and expiry dimensions, dense enough for daily surface fitting, and liquid enough for calibration work. The valid-quote plot rises through the sample, which means later dates have more quote support. We select the main calibration date using quote support across strikes and maturities, including the wings.
The spread and vega summaries describe typical observations across the panel, rather than a single maturity.
5) OTM Selection, Quote Support, and Observation Uncertainty
We now build the surface-ready panel. The main filters are:
DTE between about 14 and 150 days for calibration,
log-moneyness window around \(|k|\le 0.32\),
minimum vega,
maximum relative spread,
one OTM option leg per strike.
We select one out-of-the-money leg per strike. Calls and puts with the same strike and expiry contain almost the same volatility information once the forward and discount factor are known. Keeping both sides can double-count strikes and overweight some parts of the surface. We keep the out-of-the-money side because it’s usually more liquid in volatility terms and more directly reflects the part of the smile traders care about.
Mathematically, after parity adjustment, a call and put at the same \(K,T\) imply the same volatility if markets are frictionless. In real data they can disagree because of bid-ask, stale quotes, and exercise/microstructure details. We therefore use the OTM side as one observation of the smile:
The model panel is slightly looser than the main calibration panel because stochastic models like Merton, Heston, and Bates need tail points to identify jumps and stochastic-volatility skew. If we make the quote set too strict, the tails disappear and those models become underidentified.
OTM Quotes and Smile Construction
For a given strike and expiry, call and put prices are tied by parity:
\[
C-P=D(F-K)
\]
If parity holds exactly, the implied volatility from the call and the put must match. If they don’t match, the difference usually comes from market friction rather than a new volatility signal. It is often bid-ask noise, stale quotes, or parity friction.
OTM options are usually preferred for implied-volatility surfaces because:
OTM options are more actively quoted for volatility exposure,
ITM options have large intrinsic value that can dominate the option premium,
small errors in parity or forwards can distort ITM implied vol more visibly,
the OTM side lines up with how traders quote downside puts and upside calls.
So for \(K<F\), the put is OTM and is usually kept. For \(K>F\), the call is OTM and is usually kept. Around ATM, either side can be used, but the implementation keeps one clean quote per strike.
For stochastic models, coverage across the wings also affects which parameters can be identified. If we double-count call and put at every strike, the optimization can overweight the center and underweight the actual tail information. A model like Bates then calibrates to duplicated structure rather than to a balanced view of the surface.
main_date: 2023-12-20; clean quotes: 3,873; model quotes: 3,885; expiries: 29
The filters retain a large share of the available option data. From about 3.99 million surface-ready quotes, the DTE window keeps about 3.06 million, the moneyness window keeps about 3.02 million, and the final OTM panel has about 1.52 million quotes. The model panel is almost the same size but slightly looser on vega/spread.
The selected main date is 2023-12-20, with 3,873 clean calibration quotes, 3,885 model quotes, and 29 expiries. We have a dense surface snapshot for the selected date. It gives us enough cross-section to fit slice models, global models, and price-space stochastic models on the same day.
Show code
expiry_support = q0.groupby("expiry").agg(n=("quote_id", "size"), dte=("dte_days", "median")).sort_values("dte")fig, ax = plt.subplots(figsize=(9.0, 4.2))ax.bar(expiry_support["dte"], expiry_support["n"], width=3.5, alpha=0.82)ax.set_title(f"main-date quote support by expiry ({main_date.date()})")ax.set_xlabel("days to expiry")ax.set_ylabel("otm quotes")plt.tight_layout()plt.show()fig, ax = plt.subplots(figsize=(8.4, 4.8))sc = ax.scatter(q0["k"], q0["iv_mid"], c=q0["dte_days"], s=20, alpha=0.82, cmap="viridis")ax.set_title("main-date implied volatility by k and maturity")ax.set_xlabel("k = log strike / forward")ax.set_ylabel("implied volatility")fig.colorbar(sc, ax=ax, pad=0.01, label="days to expiry")plt.tight_layout()plt.show()
Quote support is dense across several expiries. The scatter of IV by \(k\) and maturity has the familiar SPX structure: downside puts have much higher implied vol than ATM and upside calls. The front part of the surface is also more elevated and more curved. This is the empirical shape the models have to explain.
By maturity and \(k\) bucket, we can also see which parts of the surface have the most observations. The middle-left region is the most populated: near-ATM to moderately downside strikes and short-to-medium maturities. The far right wing and far left wing are thinner. That means a model can look impressive visually while still being mostly determined by the center of the surface. Tail diagnostics later are needed because the tails don’t have the same quote density as the center.
The heatmap counts also help us distinguish dense interpolation regions from sparsely observed wings.
Support Density Across the Surface
We can use the quote counts to see which parts of the fitted surface are well supported by observations.
A dense area of the surface means many liquid quotes are telling the model the same local story. A sparse area means the model is extrapolating more. In a dense region, a residual is more likely to be meaningful. In a sparse wing, a residual can be an artifact of interpolation, quote noise, or model assumptions.
This is why the later fair-value logic includes model disagreement and quote uncertainty. The surface can always return a number, even where data support is weak. In a thin region, the fitted value depends more heavily on the chosen model and the available nearby quotes.
On the main SPX date, the dense region is short-to-medium maturity and near-ATM to moderate downside. That is the strongest calibration region. The extreme right wing and very far downside wing are thinner, so tail interpretations need more caution.
Show code
z = q0.copy()z["k_bin"] = pd.cut(z["k"], np.linspace(-max_abs_k, max_abs_k, 15))z["dte_bin"] = pd.cut(z["dte_days"], np.linspace(min_dte, max_dte, 11))h = z.pivot_table(index="dte_bin", columns="k_bin", values="quote_id", aggfunc="count", observed=False)fig, ax = plt.subplots(figsize=(9.8, 4.8))im = ax.imshow(h.fillna(0).to_numpy(), aspect="auto", origin="lower", cmap="magma")ax.set_title("quote count by maturity and k bucket")ax.set_xlabel("k bucket")ax.set_ylabel("maturity bucket")ax.set_xticks(np.arange(h.shape[1]))ax.set_xticklabels([f"{b.mid:.2f}"for b in h.columns], rotation=45, ha="right")ax.set_yticks(np.arange(h.shape[0]))ax.set_yticklabels([f"{b.mid:.0f}"for b in h.index])fig.colorbar(im, ax=ax, pad=0.01, label="quotes")plt.tight_layout()plt.show()
Bid-Ask-Aware Uncertainty
We assign each quote an uncertainty scale based on its spread and vega. This lets the calibration account for differences in the precision of observed option prices.
For implied volatility we combine two uncertainty estimates:
This way, a model isn’t punished heavily for missing a quote by less than the market’s own uncertainty, and it is punished more when it misses a tight, liquid quote.
We clip \(\epsilon_\sigma\) between 0.0025 and 0.40, and apply a small floor to \(\epsilon_w\). These floors keep a nearly zero spread or vega from giving a single quote unlimited influence.
For the price-model calibration, we use a separate scale, \(s_{V,i}=\max(\text{half-spread}_i,0.25\,\mathcal{V}_i\epsilon_{\sigma,i},0.02)\). The full \(\epsilon_V\) above is retained as a quote-uncertainty measure. Keeping the two names separate makes it clear which scale enters the optimizer and which describes the wider quote uncertainty.
Weighted Least Squares as a Statistical Model
The calibration objective can be read as a weighted Gaussian measurement model. Suppose the observed target \(y_i\) is the true model value plus quote noise:
\[
y_i=f(x_i;\theta)+\varepsilon_i
\]
with:
\[
\varepsilon_i\sim N(0,\epsilon_i^2)
\]
Then the negative log-likelihood is proportional to:
This is the weighted residual logic behind the standardized calibration residuals. A quote with large \(\epsilon_i\) has a wider noise distribution, so a larger miss is statistically less surprising. A tight quote has a small \(\epsilon_i\), so even a small miss is meaningful.
When we add quote weights \(\omega_i\), the objective becomes:
The weight can also rebalance maturity buckets or de-emphasize poor quote regions. So the calibration objective is doing two things at once:
measuring residuals in uncertainty units,
controlling which regions of the surface have influence.
The standardized residual and the observation weight therefore determine how each noisy quote contributes to the fitted model.
To see where the weight enters, define \(r_i=\sqrt{\omega_i}(f(x_i;\theta)-y_i)/\epsilon_i\). Squaring and adding the \(r_i\) gives the weighted least-squares expression above. We build \(\omega_i\) from relative spread, vega, IV uncertainty, distance into the wings, and the number of quotes in an expiry. The components are normalized and the final weight is clipped, so dense expiries and extremely tight quotes do not dominate without limit.
The Gaussian least-squares calculation is a useful starting point. The actual fits use the same standardized residuals with a soft-\(L_1\) loss, which behaves like least squares close to zero and gives less influence to large residuals.
Robust Loss and Parameter Bounds
Several calibrations use bounded least squares and robust losses. The bounded part means parameters are constrained to plausible ranges:
\[
\theta_{min}\le \theta \le \theta_{max}
\]
The bounds keep the search within the parameter ranges used for these models. For example:
\(\rho\) outside \([-1,1]\) is impossible as a correlation.
negative vol-of-vol has no meaning.
negative jump volatility has no meaning.
extremely large jump intensity can compensate for a tiny jump size in unstable ways.
The robust loss, such as a soft-\(L_1\) style loss, reduces the influence of extreme residuals. In ordinary least squares, one extreme quote can dominate because squared errors grow very fast:
\[
L(r)=r^2
\]
A robust loss behaves quadratically near zero but grows more slowly for large residuals. That is useful for option chains because some quotes are simply bad or stale even after cleaning.
The bound-warning diagnostics later tell us whether the optimizer is pushing against the edges of the parameter box. If a parameter is at a bound, the fit may be saying: the model wants something more extreme than we allowed. That is a warning about misspecification or weak identification.
For the standardized residual \(r\), the soft-\(L_1\) loss used with unit scale is proportional to
\[
\ell(r)=2\left(\sqrt{1+r^2}-1\right).
\]
For small \(r\), \(\ell(r)\approx r^2\); for large \(|r|\), it grows approximately like \(2|r|\). The optimizer can still respond to an outlying quote, but its influence no longer grows as quickly as an ordinary squared error. The reported weighted RMSE is calculated separately from the fitted errors and remains a familiar accuracy measure.
Across the maturity buckets, the uncertainty measures have a clear pattern. Median IV uncertainty is clipped at 0.0025 across buckets, while median price uncertainty rises with maturity: about 0.44 for the shortest bucket, then 0.79, 1.31, 1.74, and 1.92 for longer maturities. That is economically sensible because longer-dated options have more dollar value and higher vega, so price-scale uncertainty naturally grows even if IV uncertainty stays tight.
The observation weights vary with spread and vega. Very wide relative spreads receive lower weights, while liquid high-vega quotes receive higher weights. This is exactly what we want: the models should listen most to quotes where the market price is informative.
6) Model Targets: What Each Model Is Trying to Explain
Before fitting anything, we define the target space for each model:
Model
Main fitted object
Why this target fits the model
SVI
total variance \(w(k,T)\)
SVI is a smile parameterization in total variance
SSVI
total variance \(w(k,T)\)
SSVI imposes a global surface structure around ATM total variance
SABR
implied volatility \(\sigma_{imp}\)
Hagan’s approximation directly gives implied vol
Merton
option price \(V\)
jump-diffusion prices are generated from a distribution
Heston
simulated option price \(V\)
stochastic variance defines a risk-neutral future distribution
Bates
simulated option price \(V\)
stochastic variance plus jumps defines a richer distribution
The fitted target changes how we interpret the errors. A model that minimizes IV error and a model that minimizes price error aren’t solving exactly the same problem. IV error treats a one-vol-point miss similarly across the surface. Price error matters more where the option has high dollar sensitivity. Price-space models are closer to tradable P&L, but they are slower and can be noisier because the pricing engine has to simulate or integrate a distribution.
We use both spaces later. Fit quality in IV terms tells us how well a model explains the surface shape. Fit quality in price terms tells us how much money the miss might matter in actual option prices.
Calibration as an Inverse Problem
Every model calibration here is an inverse problem. The forward problem is:
\[
\theta \rightarrow V^{model}(\theta)
\]
Given parameters \(\theta\), we can compute model prices or implied vols. Calibration reverses this:
\[
V^{mkt} \rightarrow \hat{\theta}
\]
This inverse problem is often ill-conditioned. Small quote changes can move parameters a lot, especially when parameters affect similar parts of the surface. For example:
in Heston, more negative \(\rho\) and higher \(\xi\) can both steepen the left skew,
in Bates, negative jumps and stochastic variance can both raise OTM put prices,
in SABR, \(\rho\) and \(\nu\) can trade off between skew and curvature,
in SVI, \(b\), \(\rho\), \(m\), and \(\sigma\) can combine to produce similar-looking smiles.
So the diagnostic tables are part of the modeling, not decoration. Success flags, bound warnings, runtime, residual heatmaps, Monte Carlo standard error, and holdout error are all ways of checking whether the inverse problem produced a stable result.
If two nearby parameter vectors produce almost the same fitted prices, the local price sensitivities with respect to those parameters point in similar directions. The optimizer can move along that nearly flat direction with little change in loss. Adding strikes and expiries where those sensitivities differ helps identify the parameters, while multiple starting values help us see whether a nonlinear fit settles on a better solution.
Parameter Identification and Economic Plausibility
Identification describes whether the observed quotes can distinguish one parameter effect from another. Several parameters can produce similar changes in a fitted volatility smile.
For SVI, identification is usually acceptable because each parameter has a geometric role and we have many strike observations per expiry. But even there, parameters can trade off. A shift in \(m\) and a change in \(\rho\) can both affect where the smile is steepest.
For Heston, identification is harder. The option surface may say “I need a steep left wing,” and the model can create that through several combinations:
These changes don’t have identical effects, but they overlap enough that a small quote panel can’t always separate them cleanly. Bates adds even more overlap because jumps can also steepen tails:
A parameter set is more believable when it fits well, avoids bounds, has reasonable standard errors/noise, and remains stable across nearby dates. This project doesn’t estimate full parameter uncertainty distributions, but it does record warning flags and MC standard errors so we don’t overstate precision.
7) SVI: Slice-Level Smile Modeling in Total Variance
The first model family is SVI, which stands for stochastic volatility inspired. Despite the name, in this implementation SVI is a static smile parameterization. For each expiry, it fits a curve in log-moneyness \(k\):
where \(w(k)=\sigma_{imp}^2(k,T)T\) is total implied variance for a fixed maturity. Each parameter has a visible shape role:
\(a\) lifts the whole curve up or down.
\(b\) controls the overall slope/wing steepness.
\(\rho\) controls skew direction. For SPX, \(\rho\) is usually negative because downside puts are expensive.
\(m\) shifts the center of the smile left or right in moneyness.
\(\sigma\) controls the smoothness/curvature around the smile center.
The square-root term is the main reason SVI is useful. It behaves smoothly near the center but becomes close to linear in the wings. That lets the curve fit a U-shaped smile or skew without exploding as badly as a high-order polynomial.
For a fixed expiry, the market gives us observations:
\[
(k_i,w_i), \quad i=1,\ldots,n
\]
and we can express the weighted least-squares objective for parameters \(\theta=(a,b,\rho,m,\sigma)\) as:
The term \(\epsilon_{w,i}\) is the total-variance uncertainty from bid-ask and vega. So the residual is measured in “number of quote-uncertainty units,” not raw variance units.
If a quote is very precise, missing it matters. If a quote is noisy, missing it slightly doesn’t mean the model is bad.
Once a slice is fitted, its implied volatility at maturity \(T\) is recovered as \(\sigma_{SVI}(k,T)=\sqrt{w(k)/T}\). The optimization is carried out in total variance, but the resulting IV and Black-76 price can be compared directly with the observed IV and mid price.
The displayed squared objective explains the standardized residual. In the actual slice fit we pass \(\sqrt{\omega_i}(w_{SVI}-w_i)/\epsilon_{w,i}\) to a bounded soft-\(L_1\) optimizer and add a penalty when fitted variance becomes nonpositive on a grid over the observed slice.
SVI Wing Slopes and Curvature
SVI is especially useful because the left and right wings become almost linear for large \(|k-m|\). For large positive \(x=k-m\):
\[
\sqrt{x^2+\sigma^2}\approx |x|
\]
so on the right wing:
\[
w(k)\approx a+b(\rho x+x)=a+b(1+\rho)x
\]
and on the left wing, where \(x<0\):
\[
w(k)\approx a+b(\rho x-x)=a+b(\rho-1)x
\]
Because \(x\) is negative on the left, a negative value of \((\rho-1)x\) raises total variance. If \(\rho<0\), the left wing becomes steeper than the right wing. That is the equity-index skew.
The wing slope intuition is:
\[
\text{right wing slope}\approx b(1+\rho)
\]
\[
\text{left wing slope}\approx b(\rho-1)
\]
A more negative \(\rho\) lowers the right wing slope and steepens the left wing. That is why \(\rho\) is one of the most economically meaningful SVI parameters for SPX.
The parameter \(\sigma\) controls how sharp the transition is around \(m\). If \(\sigma\) is small, the curve can bend more sharply near the center. If it is large, the curve is smoother and flatter around the minimum.
So SVI has a clear geometry: level, wing steepness, skew direction, center shift, and curvature width.
We can make the shape more explicit by differentiating with respect to \(k\). With \(x=k-m\),
For \(b>0\) and \(\sigma>0\) the total-variance curve is convex in \(k\). If \(|\rho|<1\), its minimum is at \(k_*=m-\rho\sigma/\sqrt{1-\rho^2}\), so \(m\) locates the center of the square-root term but need not be the smile minimum. As \(x\) moves far from zero, the first derivative approaches the two wing slopes above.
SVI as a Fair-Value Core
SVI is a useful candidate for the fair-value core because it combines flexible slice fits, speed, and interpretable local residuals.
First, it is flexible enough to fit the smile shape in each expiry without needing too many parameters. Five parameters per expiry is a good compromise: more flexible than a quadratic curve, less fragile than a high-order polynomial.
Second, it is fast. Even fitting 29 expiries and thousands of quotes is feasible. That matters because the validation section fits many daily surfaces.
Third, it has interpretable residuals. If SVI misses a point, that miss is measured relative to a local smooth smile. This is close to how option traders think: they build a smooth volatility curve and look for quotes that sit above or below it.
The limitation is that SVI residuals can reflect the model’s own flexibility. A flexible model can bend around noise. That is why SVI is combined with holdout scoring, model disagreement, and validation rather than being treated as a final arbitrage signal by itself.
Show code
def svi_w_demo(k, a, b, rho, m, sigma):return a + b * (rho * (k - m) + np.sqrt((k - m) **2+ sigma **2))k_line = np.linspace(-0.35, 0.35, 240)tau_demo =60.0/ ann_daysp_demo = [ (0.010, 0.030, -0.60, 0.000, 0.120, "balanced skew"), (0.008, 0.045, -0.85, -0.030, 0.080, "strong downside wing"), (0.012, 0.025, 0.20, 0.020, 0.160, "gentler curvature"),]fig, ax = plt.subplots(figsize=(8.4, 4.4))for a, b, rho, m, sigma, label in p_demo: ax.plot(k_line, np.sqrt(np.maximum(svi_w_demo(k_line, a, b, rho, m, sigma), 1e-8) / tau_demo), lw=2.0, label=label)ax.set_title("raw svi total-variance smiles")ax.set_xlabel("k")ax.set_ylabel("implied volatility")ax.legend()plt.tight_layout()plt.show()
Changing the SVI parameters changes the shape of each smile. The strong downside-wing curve has a much steeper left wing, which is the SPX crash-protection shape. The balanced version has a smoother skew, and the gentler curvature version is flatter around ATM.
The five SVI parameters summarize the shape of one expiry’s smile. They describe the current cross-section without specifying how volatility evolves through time. That is already valuable because a raw expiry can have 100 to 200 quotes. Compressing those quotes into five interpretable parameters gives us a much more manageable object for diagnostics.
The plot converts fitted total variance back to IV at the fixed demo maturity.
SVI Fitting Details
For every expiry we:
sort quotes by \(k\),
fit total variance \(w(k)\),
use multiple starting points to avoid a weak local minimum,
include previous expiry parameters as a warm start when possible,
record weighted IV RMSE, median absolute IV error, P90 absolute IV error, and parameter-bound warnings.
Adjacent expiries often have related smile shapes, so the previous fit is a useful additional starting point. If the 45-day smile and the 50-day smile are wildly different, that’s often a data issue or optimization issue rather than true market structure. Reusing the previous slice as a candidate start helps stabilize the term structure of parameters.
SVI and Static Arbitrage
Option surfaces have basic no-arbitrage requirements. We don’t re-build the full no-arbitrage machinery here, but it helps to understand what can go wrong.
For a fixed maturity, call prices must be decreasing and convex in strike:
\[
\frac{\partial C}{\partial K}\le 0
\]
\[
\frac{\partial^2 C}{\partial K^2}\ge 0
\]
The second condition is butterfly arbitrage: if the call price curve isn’t convex, a butterfly spread can have negative cost and nonnegative payoff.
Across maturities, total variance should usually be nondecreasing in maturity for the same moneyness coordinate:
\[
\frac{\partial w(k,T)}{\partial T}\ge 0
\]
This is linked to calendar-spread arbitrage. Longer maturity options shouldn’t be cheaper in a way that creates free time value.
Raw SVI slice fitting doesn’t automatically enforce all cross-maturity conditions. It can fit each expiry beautifully and still create an irregular term structure. We next fit SSVI to connect the expiry smiles through a shared surface structure. We keep both models: SVI gives a flexible local fit, while SSVI uses common shape parameters across maturities.
The convexity requirement concerns call price as a function of strike; convexity of the SVI total-variance formula in \(k\) alone doesn’t establish it. Likewise, the calendar condition is considered at fixed forward moneyness under the usual carry assumptions. We fit individual slices and check variance positivity on an observed-strike grid; this fit doesn’t explicitly impose every butterfly or calendar-arbitrage condition.
Raw SVI for Local Fair Values
Raw slice SVI can create cross-maturity irregularities, but it remains extremely useful because this project’s trading signal is local and short-horizon. We aren’t using raw SVI to price long-dated exotic books. We are using it to estimate a fair value for listed options on a dense surface date.
The model’s intended use determines which properties we need to examine most closely. For exotic pricing, an arbitrage-free global surface is needed. For daily residual ranking, local fit and holdout validation can be more important.
We assess the raw SVI fit alongside three additional checks:
SSVI is fitted as a disciplined global cross-check.
Residuals are uncertainty-adjusted with model disagreement.
Historical validation checks whether residuals work out of sample.
So raw SVI is used where it is strongest, and its weaknesses are controlled by diagnostics rather than ignored.
Local Flexibility and Overfitting Risk
SVI’s flexibility is useful, but it creates overfitting risk. If a slice has a few noisy quotes, SVI might bend toward them. That can create a false residual nearby or hide a true residual by absorbing it into the fitted curve.
We assess the fit using:
quote weights,
multiple starting values,
parameter-bound checks,
residual plots,
SSVI/SABR comparison,
holdout scoring,
historical validation.
A good in-sample curve is a starting point. The additional comparisons and holdout results help us assess whether the fitted residuals remain useful beyond the quotes used to fit the slice.
This is also why raw calibration RMSE isn’t enough. A model with very low in-sample error can be worse for trading if it overfits noise. The holdout and validation sections are the protection against that.
The diagnostics also help separate a failed numerical fit from a fitted curve that misses specific regions.
Show code
@njitdef svi_w(k, a, b, rho, m, sigma): out = np.empty(k.size, dtype=np.float64)for i inrange(k.size): x = k[i] - m out[i] = a + b * (rho * x + math.sqrt(x * x + sigma * sigma))return out@njitdef svi_iv(k, tau, a, b, rho, m, sigma): w = svi_w(k, a, b, rho, m, sigma) out = np.empty(k.size, dtype=np.float64)for i inrange(k.size): out[i] = math.sqrt(max(w[i], 1e-12) /max(tau[i], 1e-12))return outt0 = tic()_ = svi_w(np.array([0.0, 0.1]), 0.01, 0.04, -0.4, 0.0, 0.12)svi_compile_time = toc("compile svi", t0, 2)
Show code
def wrmse(x, wt): x = np.asarray(x, dtype=float) wt = np.asarray(wt, dtype=float) ok = np.isfinite(x) & np.isfinite(wt) & (wt >0)if ok.sum() ==0:return np.nanreturnfloat(np.sqrt(np.average(x[ok] **2, weights=wt[ok])))def near_bound(p, lo, hi, tol=2e-3): p = np.asarray(p, dtype=float) lo = np.asarray(lo, dtype=float) hi = np.asarray(hi, dtype=float) scale = np.maximum(hi - lo, 1e-12)returnbool(np.any((p - lo) / scale < tol) or np.any((hi - p) / scale < tol))
Show code
def svi_starts(q, p_prev=None): k = q["k"].to_numpy(dtype=float) w = q["total_var"].to_numpy(dtype=float) iv = q["iv_mid"].to_numpy(dtype=float) i0 =int(np.nanargmin(np.abs(k))) w0 =float(np.clip(w[i0], 1e-5, 2.0)) slope =float(np.polyfit(k, w, 1)[0]) iflen(q) >=5else0.0 starts = [ np.array([0.35* w0, max(0.02, 1.2*abs(slope)), -0.50, 0.0, 0.10]), np.array([0.20* w0, max(0.03, 0.8* np.nanmedian(iv) **2), -0.75, float(k[i0]), 0.08]), np.array([0.55* w0, max(0.01, 0.5*abs(slope) +0.015), -0.25, float(np.nanmedian(k)), 0.18]), ]if p_prev isnotNoneand np.all(np.isfinite(p_prev)): starts.insert(0, np.asarray(p_prev, dtype=float))return startsdef fit_svi_slice(q, p_prev=None): z = q.sort_values("k").copy() k = z["k"].to_numpy(dtype=float) y = z["total_var"].to_numpy(dtype=float) eps = z["eps_w"].to_numpy(dtype=float) eps_med =float(np.nanmedian(eps[np.isfinite(eps) & (eps >0)])) if np.isfinite(eps).any() else1e-4 eps = np.clip(eps, max(1e-5, 0.25* eps_med), 0.25) root_wt = np.sqrt(z["weight"].to_numpy(dtype=float)) lo = np.array([-0.20, 1e-5, -0.999, -0.80, 1e-4]) hi = np.array([2.00, 5.00, 0.999, 0.80, 2.00])def err(p): w = np.asarray(svi_w(k, p[0], p[1], p[2], p[3], p[4]), dtype=float) r = root_wt * (w - y) / eps r = np.where((~np.isfinite(w)) | (w <=1e-8), 1e3, r) k_grid = np.linspace(np.nanmin(k), np.nanmax(k), 25) w_grid = np.asarray(svi_w(k_grid, p[0], p[1], p[2], p[3], p[4]), dtype=float) pen =max(0.0, 1e-7-float(np.nanmin(w_grid))) *1e5return np.r_[r, pen] best =Nonefor p0 in svi_starts(z, p_prev): t0 = tic() res = optimize.least_squares(err, np.clip(p0, lo +1e-8, hi -1e-8), bounds=(lo, hi), loss="soft_l1", f_scale=1.0, max_nfev=180, xtol=1e-8, ftol=1e-8) item = {"res": res, "cost": float(res.cost), "seconds": time.perf_counter() - t0}if best isNoneor item["cost"] < best["cost"]: best = itemreturn best["res"].x, best["res"].success, best["seconds"], best["res"].nfev, lo, hi
The SVI fits are close to the observed smiles. The first displayed expiries all succeed, and many weighted IV RMSE values are tiny. For example, several short expiries have weighted IV RMSE around 0.0005 to 0.0011, and even the weaker short slice shown has about 0.0039. Median absolute IV errors are mostly below one vol point and often only a few basis points of vol.
Across the selected high-support expiries, SVI follows the observed market smiles closely. The residual scatter is centered around zero and mostly tight. The few larger residuals appear in thinner parts of the wing, which is expected because tails are harder to fit and quotes are less dense there.
The fitted parameters move unevenly across maturity, even when the fitted smiles themselves remain close to their observed quotes. That is the classic SVI tradeoff: slice-level flexibility gives excellent fit, but it doesn’t automatically guarantee a smooth, arbitrage-aware surface across maturities.
The comparison is made slice by slice, so each maturity has its own fitted parameters.
Reading SVI Residuals Like a Surface Diagnostic
We can use the SVI residual pattern to locate the surface features that are harder to fit.
If residuals are randomly scattered around zero, the model is flexible enough and remaining errors are mostly quote noise. If residuals form a pattern by \(k\), the smile shape is wrong. If residuals form a pattern by maturity, the term structure is unstable. If residuals are concentrated in one wing, the model is underfitting tail curvature.
In the displayed residual scatter, most observations are close to zero, with a few larger wing misses. The close fit supports using SVI as the fair-value core, while the larger wing misses give us a reason to compare the tail estimates with other models.
The parameter-term plots are also useful. Smooth parameter paths across maturity would be ideal, but the output shows some movement and noise. We shouldn’t overinterpret individual SVI parameters at one expiry. The fitted surface and residuals are more reliable than a narrative about one noisy parameter series.
The residuals use model IV minus market IV; positive values correspond to a model curve above the quote.
Show code
selected_expiries = svi_diag.sort_values("n_quotes", ascending=False).head(5)["expiry"].tolist()fig, ax = plt.subplots(figsize=(8.8, 5.0))for expiry in selected_expiries: g = svi_fit[svi_fit["expiry"].eq(expiry)].sort_values("k")if g.empty:continue ax.scatter(g["k"], g["iv_mid"], s=18, alpha=0.55) ax.plot(g["k"], g["model_iv"], lw=2.0, label=f"{g['dte_days'].median():.0f}d")ax.set_title("raw svi fitted smiles versus market iv")ax.set_xlabel("k")ax.set_ylabel("implied volatility")ax.legend(title="maturity")plt.tight_layout()plt.show()fig, ax = plt.subplots(figsize=(8.4, 4.4))ax.scatter(svi_fit["k"], svi_fit["iv_residual"], c=svi_fit["dte_days"], cmap="viridis", s=18, alpha=0.80)ax.axhline(0.0, lw=0.9, alpha=0.70)ax.set_title("raw svi iv residuals by k")ax.set_xlabel("k")ax.set_ylabel("model iv - market iv")plt.tight_layout()plt.show()
Show code
fig, axes = plt.subplots(2, 3, figsize=(13.4, 6.2), sharex=True)for ax, col inzip(axes.ravel(), ["a", "b", "rho", "m", "sigma"]): ax.plot(svi_params["dte_days"], svi_params[col], marker="o", lw=1.6) ax.set_title(col) ax.set_xlabel("days to expiry")axes.ravel()[-1].axis("off")plt.tight_layout()plt.show()fig, ax = plt.subplots(figsize=(8.6, 4.2))ax.bar(svi_diag["dte_days"], svi_diag["weighted_rmse_iv"], width=3.0, alpha=0.82)ax.set_title("raw svi weighted iv rmse by expiry")ax.set_xlabel("days to expiry")ax.set_ylabel("weighted iv rmse")plt.tight_layout()plt.show()
8) SSVI: Global Surface Discipline
SVI fits each expiry separately. SSVI adds a shared global structure across maturities. Instead of giving each expiry its own five parameters, SSVI expresses the total variance smile through the ATM total variance \(\theta\) and a smaller set of global parameters.
\(\theta\) says how much variance the market prices at the center for each maturity.
\(\rho\) says whether downside or upside strikes are more expensive.
\(\eta\) controls how steep the smile is.
\(\gamma\) controls how that steepness decays or changes as total variance grows with maturity.
This is more restrictive than SVI. That restriction is useful when we want a smooth surface, but it can be costly when the actual market has maturity-specific irregularities.
Substituting \(k=0\) into the surface formula gives \(w(0,\theta)=\theta\), since the term under the square root becomes one. This is why \(\theta\) is the ATM total variance. Differentiating at the center gives \(\partial w/\partial k\vert_{k=0}=\rho\theta\phi(\theta)\): \(\rho\) determines the sign of the ATM skew, while \(\theta\phi(\theta)\) sets its magnitude.
The Role of \(\theta\) in SSVI
The key SSVI idea is to separate the surface into:
ATM total variance by maturity,
a global rule that shapes the smile around it.
The ATM total variance is:
\[
\theta(T)=\sigma_{ATM}^2(T)T
\]
If \(\theta\) rises with maturity, the total variance budget rises. SSVI then asks how skew changes as that total variance budget changes. That is much more structured than fitting every expiry independently.
makes skew level depend on maturity. When \(\theta\) is small, short maturity smiles can be steep. When \(\theta\) is large, longer maturity smiles usually become smoother. The parameter \(\gamma\) controls how quickly that transition happens.
This helps with calendar stability because the model ties smiles together through \(\theta(T)\). If the 21-day and 30-day smiles are very different, SSVI doesn’t freely fit both unless the ATM total variance curve and global parameters justify it. That can reduce overfitting, especially for sparse or noisy surfaces.
The cost is bias. Real SPX smiles can have expiration-specific effects, event risk, dealer positioning, and liquidity shocks. SSVI smooths through some of those details.
We estimate \(\theta\) for each expiry from the weighted average of quotes within \(|k|\le0.045\) when there are at least three; otherwise we use up to the five nearest quotes. Those expiry values are inputs to the joint fit of \(\rho,\eta,\gamma\). Tying the smile parameters together provides a common surface structure, but the fitted ATM values are not forced to be monotone and the optimizer doesn’t explicitly enforce all no-arbitrage conditions.
Show code
def theta_from_atm(q): rows = []for expiry, g in q.groupby("expiry", sort=True): near = g[np.abs(g["k"]) <=0.045].copy() use = near iflen(near) >=3else g.iloc[np.argsort(np.abs(g["k"].to_numpy(dtype=float)))[: min(5, len(g))]].copy() wt = np.maximum(use["weight"].to_numpy(dtype=float), 1e-8) theta =float(np.average(use["total_var"].to_numpy(dtype=float), weights=wt)) rows.append({"expiry": expiry, "tau": float(g["tau"].median()), "dte_days": float(g["dte_days"].median()),"theta_ssvi": max(theta, 1e-8)})return pd.DataFrame(rows).sort_values("tau").reset_index(drop=True)@njitdef ssvi_w(k, theta, rho, eta, gamma): out = np.empty(k.size, dtype=np.float64)for i inrange(k.size): th =max(theta[i], 1e-12) phi = eta / (th ** gamma * (1.0+ th) ** (1.0- gamma)) x = phi * k[i] inside = (x + rho) * (x + rho) +1.0- rho * rho out[i] =0.5* th * (1.0+ rho * x + math.sqrt(max(inside, 1e-12)))return outt0 = tic()_ = ssvi_w(np.array([0.0, 0.1]), np.array([0.02, 0.03]), -0.5, 1.0, 0.5)ssvi_compile_time = toc("compile ssvi", t0, 2)
Show code
def fit_ssvi(q): theta_curve = theta_from_atm(q) z = q.merge(theta_curve[["expiry", "theta_ssvi"]], on="expiry", how="left").copy() k = z["k"].to_numpy(dtype=float) theta = z["theta_ssvi"].to_numpy(dtype=float) y = z["total_var"].to_numpy(dtype=float) eps = z["eps_w"].to_numpy(dtype=float) eps_med =float(np.nanmedian(eps[np.isfinite(eps) & (eps >0)])) if np.isfinite(eps).any() else1e-4 eps = np.clip(eps, max(1e-5, 0.25* eps_med), 0.25) root_wt = np.sqrt(z["weight"].to_numpy(dtype=float)) lo = np.array([-0.999, 0.001, 0.00]) hi = np.array([0.999, 25.000, 1.00])def err(p): w = np.asarray(ssvi_w(k, theta, p[0], p[1], p[2]), dtype=float) r = root_wt * (w - y) / eps r = np.where((~np.isfinite(w)) | (w <=1e-8), 1e3, r)return np.r_[r, max(0.0, 1e-8-float(np.nanmin(w))) *1e5] starts = [np.array([-0.55, 1.2, 0.50]), np.array([-0.75, 2.0, 0.35]), np.array([-0.35, 0.8, 0.70]), np.array([0.05, 1.0, 0.50])] best =Nonefor p0 in starts: t0 = tic() res = optimize.least_squares(err, np.clip(p0, lo +1e-8, hi -1e-8), bounds=(lo, hi), loss="soft_l1", f_scale=1.0, max_nfev=160, xtol=1e-8, ftol=1e-8) item = {"res": res, "cost": float(res.cost), "seconds": time.perf_counter() - t0}if best isNoneor item["cost"] < best["cost"]: best = itemreturn best["res"].x, theta_curve, best["res"].success, best["seconds"], best["res"].nfev, lo, hi
For SSVI, we obtain \(\rho\approx -0.497\), \(\eta\approx 2.98\), and \(\gamma\approx 0.31\), with success and no parameter-bound warning. The negative \(\rho\) is exactly what we expect for SPX: downside skew is structurally expensive.
Fit quality is weaker than raw SVI: weighted IV RMSE is around 0.0083, median absolute IV error around 0.0038, and P90 error around 0.0165. The residual heatmap contains systematic patches across maturity and moneyness buckets. That’s the price of global discipline.
SSVI is useful as a surface regularizer and a comparison with the more flexible slice fits. Its shared parameters provide a structured surface representation, while the fitted term structure still needs the no-arbitrage checks discussed above. If we need quote-level relative value signals, raw SVI is much sharper in this SPX snapshot.
SSVI as a Global Cross-Check
SSVI’s larger error comes with a shared surface structure that gives us another way to assess local residuals. If SVI says an option is cheap but SSVI says that quote is only cheap because SVI is bending too much locally, then the trade is less convincing.
A rigid model can act as a regularization benchmark. We can think of SSVI as asking:
\[
\text{Can a globally disciplined surface explain this quote?}
\]
SVI asks:
\[
\text{Can a flexible local smile explain this quote?}
\]
If both agree that the quote is rich or cheap, the residual is more robust. If only the flexible model sees it, the residual may just be an artifact of local overfitting.
That is why the later ensemble gives SSVI a tiny weight rather than zero. It is too biased for primary pricing, but it can still contribute model disagreement and discipline.
Show code
fig, ax = plt.subplots(figsize=(8.8, 5.0))for expiry in selected_expiries: g = ssvi_fit[ssvi_fit["expiry"].eq(expiry)].sort_values("k")if g.empty:continue ax.scatter(g["k"], g["iv_mid"], s=18, alpha=0.50) ax.plot(g["k"], g["model_iv"], lw=2.0, label=f"{g['dte_days'].median():.0f}d")ax.set_title("ssvi fitted smiles across expiries")ax.set_xlabel("k")ax.set_ylabel("implied volatility")ax.legend(title="maturity")plt.tight_layout()plt.show()z = ssvi_fit.copy()z["k_bin"] = pd.cut(z["k"], np.linspace(-max_abs_k, max_abs_k, 15))z["dte_bin"] = pd.cut(z["dte_days"], np.linspace(min_dte, max_dte, 11))h = z.pivot_table(index="dte_bin", columns="k_bin", values="iv_residual", aggfunc="median", observed=False)lim =float(np.nanquantile(np.abs(h.to_numpy(dtype=float)), 0.95)) if np.isfinite(h.to_numpy(dtype=float)).any() else0.02fig, ax = plt.subplots(figsize=(9.8, 4.8))im = ax.imshow(h.to_numpy(dtype=float), aspect="auto", origin="lower", cmap="coolwarm", vmin=-lim, vmax=lim)ax.set_title("ssvi median iv residual heatmap")ax.set_xlabel("k bucket")ax.set_ylabel("maturity bucket")ax.set_xticks(np.arange(h.shape[1]))ax.set_xticklabels([f"{b.mid:.2f}"for b in h.columns], rotation=45, ha="right")ax.set_yticks(np.arange(h.shape[0]))ax.set_yticklabels([f"{b.mid:.0f}"for b in h.index])fig.colorbar(im, ax=ax, pad=0.01, label="model iv - market iv")plt.tight_layout()plt.show()
Show code
svi_ssvi = svi_diag[["expiry", "dte_days","weighted_rmse_iv"]].rename(columns={"weighted_rmse_iv": "svi"}).merge(ssvi_by_expiry[["expiry","weighted_rmse_iv"]].rename(columns={"weighted_rmse_iv": "ssvi"}), on="expiry", how="inner")fig, ax = plt.subplots(figsize=(8.4, 4.4))ax.plot(svi_ssvi["dte_days"], svi_ssvi["svi"], marker="o", lw=1.8, label="raw svi")ax.plot(svi_ssvi["dte_days"], svi_ssvi["ssvi"], marker="o", lw=1.8, label="ssvi")ax.set_title("raw svi versus ssvi error by expiry")ax.set_xlabel("days to expiry")ax.set_ylabel("weighted iv rmse")ax.legend()plt.tight_layout()plt.show()
Across maturities, SVI errors stay much closer to zero than SSVI errors. SSVI errors are especially large at the short end and remain above SVI through most of the term structure. The shape is economically informative: short-maturity SPX smiles often have event risk, crash protection demand, and discrete supply-demand imbalances that a global three-parameter structure can’t fully absorb.
A clean reading is:
SVI is better for local quote-level pricing on this day.
SSVI imposes a coherent shared shape across maturities.
Using SSVI alone for trading residuals would create many false signals because some residuals are just model rigidity.
The comparison uses weighted IV RMSE for each expiry, on a common volatility scale.
9) SABR: Stochastic Smile Dynamics Through an Approximation
SABR starts from stochastic dynamics for both the forward and its volatility scale. The standard SABR dynamics for a forward \(F_t\) are:
\[
dF_t=\alpha_t F_t^\beta dW_t^F
\]
\[
d\alpha_t=\nu\alpha_t dW_t^\alpha
\]
with correlation:
\[
dW_t^F dW_t^\alpha=\rho dt
\]
The parameters mean:
\(F_t\) is the forward price.
\(\alpha_t\) is the instantaneous volatility scale.
\(\beta\) controls how volatility changes with the level of the forward.
\(\nu\) is volatility-of-volatility.
\(\rho\) controls the skew through correlation between forward shocks and volatility shocks.
If \(\rho<0\), a negative forward move tends to come with a positive volatility shock. That creates left skew, which is exactly the equity-index behavior we usually see.
SABR is usually used through an approximate implied-volatility formula. We don’t simulate the SABR process here. We fit the implied-volatility approximation to the market smile. The approximation gives us a fast way to compare a dynamically motivated smile with the observed quotes.
Hagan-Style SABR Implied Volatility
The full SABR process doesn’t have a simple exact closed-form implied-volatility formula for every case, so practitioners often use the Hagan approximation. The formula is built around:
The first fraction scales the level and adjusts for moneyness; \(z/x(z)\) shapes the smile away from ATM; the final bracket adds the maturity correction. The expression is an approximation, so the fitted parameters are interpreted through the accuracy of this formula as well as the observed quotes.
The implied volatility is roughly an ATM scaling term multiplied by a correction involving \(z/x(z)\) and maturity adjustments. The approximation combines these terms with the following parameter effects:
\(\alpha\) sets the ATM volatility level,
\(\beta\) controls how the level changes with the forward,
\(\rho\) controls skew through spot-vol correlation,
\(\nu\) controls curvature through volatility-of-volatility.
Near ATM, \(\log(F/K)\) is close to zero, so \(z\) is close to zero and the formula behaves smoothly. Away from ATM, \(z/x(z)\) adjusts the vol to create skew and curvature.
This is why SABR is popular: it gives a fast implied-volatility formula with parameters that still have dynamic meaning. It sits between a pure surface fit and a full Monte Carlo stochastic-volatility model.
At the money, \(z\to0\) and \(x(z)\to0\), but their ratio tends to one. We use the ATM limit directly to avoid a numerically unstable \(0/0\) calculation. Its leading term is \(\alpha/F^{1-\beta}\), followed by the maturity correction in the approximation. Thus \(\alpha\) itself equals the leading ATM volatility only when \(\beta=1\); for other beta values it also depends on the forward’s scale.
For strikes away from ATM, the implementation uses \(z/x(z)\) together with the log-moneyness denominator and the first-order maturity correction. These terms produce a fast approximate Black implied volatility, which is then converted to a price for comparison.
SABR and Smile Movement
SABR became popular partly because it gives more realistic smile movement than simple local-volatility intuition. The core reason is that volatility has its own random driver and is correlated with the forward.
If the forward falls and \(\rho<0\), volatility tends to rise. That means the smile can move with the underlying in a way that looks closer to observed equity and rate markets. The correlation connects changes in the forward to changes in its volatility scale.
In formula form, the covariance between the forward shock and volatility shock is:
When \(\rho\) is negative, negative forward moves are associated with positive volatility moves. That creates skew and affects hedging. A delta hedge derived from a model with stochastic smile movement can behave differently from one derived from a sticky-strike or sticky-delta rule.
We don’t build the full dynamic SABR hedging system here. We use SABR as a structural smile model, but understanding its dynamic origin helps explain why it is more meaningful than a generic curve fit.
SABR Parameter Intuition
The parameter \(\beta\) deserves careful explanation. It controls the elasticity of volatility to the forward level.
If \(\beta=0\), the forward has approximately normal volatility. Absolute price moves are more stable.
If \(\beta=1\), the forward has lognormal volatility. Percentage moves are more stable.
If \(0<\beta<1\), the model sits between normal and lognormal behavior.
For equity index options, \(\beta=1\) can be natural because prices are often modeled in log-return space. But equity smiles can be better fit by lower beta values because the model then gets a different balance between skew and curvature.
The fitted SABR formula is sensitive to \(\rho\) and \(\nu\):
More negative \(\rho\) steepens downside skew.
Larger \(\nu\) increases smile curvature because volatility itself can move more violently.
For example, two SABR smiles can have the same ATM volatility but different wing behavior. A higher \(\nu\) makes both wings more expensive because future volatility is more uncertain. A more negative \(\rho\) makes the left wing more expensive than the right wing because bad spot moves tend to coincide with higher volatility.
SABR Calibration by Expiry
SABR is usually fitted slice by slice, just like SVI. For a fixed expiry, the target is implied volatility:
In this project, \(\beta\) is chosen from a small grid and then held fixed inside each slice fit. This is common because \(\beta\) can be weakly identified. If we try to estimate \(\alpha\), \(\beta\), \(\rho\), and \(\nu\) all at once for every expiry, the optimizer can trade off parameters in unstable ways.
The beta grid is:
\[
\beta\in \{0.5,\ 0.7,\ 1.0\}
\]
We compare the median errors and bound warnings across the beta grid, and use the fitted surface to choose a practical setting.
For each fixed \(\beta\), we start the nonlinear fit from several parameter combinations. The starting \(\alpha\) is scaled from the near-ATM IV as \(\alpha_0\approx\sigma_{ATM}F^{1-\beta}\). Once the three beta fits are compared, we retain \(\beta=1\) when its median error is within 8% of the best median error; otherwise we use the best beta. The quoted squared objective describes the residuals, and the implemented fit applies soft-\(L_1\) to those standardized residuals.
SABR and SVI in Signal Design
SABR and SVI can both fit smiles, but they give different types of residual information.
SVI residuals are relative to a very flexible local curve. If a quote is far away from SVI, it is genuinely off the local smooth smile. SABR residuals are relative to a more structural stochastic-volatility approximation. If SABR misses a quote, the miss may reflect a real market irregularity or simply SABR’s lower flexibility.
Combining them is useful. If both models agree a quote is cheap, then the signal is more robust. If SVI says cheap but SABR says rich or fair, model disagreement rises and the total-error filter weakens the trade. The comparison lets us reflect that disagreement in the residual score.
This is also why the ensemble weight is capped. SVI is strong enough to lead, but allowing a second model to contribute prevents the signal from being purely one-model residual chasing.
Show code
def sabr_iv_demo(f, k_strike, tau, alpha, beta, rho, nu): k_strike = np.asarray(k_strike, dtype=float) log_fk = np.log(np.clip(f / k_strike, 1e-12, None)) fk_beta = (f * k_strike) ** ((1.0- beta) /2.0) z = (nu /max(alpha, 1e-12)) * fk_beta * log_fk xz = np.log((np.sqrt(1.0-2.0* rho * z + z * z) + z - rho) / (1.0- rho)) z_xz = np.where(np.abs(z) <1e-8, 1.0, z / xz) log_adj =1.0+ ((1.0- beta) **2/24.0) * log_fk **2+ ((1.0- beta) **4/1920.0) * log_fk **4 term = ((1.0- beta) **2/24.0) * alpha **2/ np.maximum(fk_beta **2, 1e-12) term +=0.25* rho * beta * nu * alpha / np.maximum(fk_beta, 1e-12) term += (2.0-3.0* rho * rho) * nu * nu /24.0return np.maximum(alpha / np.maximum(fk_beta * log_adj, 1e-12) * z_xz * (1.0+ term * tau), 1e-5)f_demo =100.0strike_demo = f_demo * np.exp(k_line)fig, axes = plt.subplots(1, 3, figsize=(14.2, 4.2), sharey=True)for rho_value in [-0.85, -0.40, 0.10]: axes[0].plot(k_line, sabr_iv_demo(f_demo, strike_demo, tau_demo, 0.24, 1.0, rho_value, 0.75), lw=2.0, label=f"rho={rho_value:+.2f}")for nu_value in [0.25, 0.75, 1.35]: axes[1].plot(k_line, sabr_iv_demo(f_demo, strike_demo, tau_demo, 0.24, 1.0, -0.55, nu_value), lw=2.0, label=f"nu={nu_value:.2f}")for beta_value in [1.0, 0.7, 0.5]: axes[2].plot(k_line, sabr_iv_demo(f_demo, strike_demo, tau_demo,0.24* f_demo ** (1.0- beta_value), beta_value, -0.55,0.75), lw=2.0, label=f"beta={beta_value:.1f}")axes[0].set_title("sabr rho and skew")axes[1].set_title("sabr nu and curvature")axes[2].set_title("fixed beta scenarios")for ax in axes: ax.set_xlabel("k") ax.legend(fontsize=8)axes[0].set_ylabel("implied volatility")plt.tight_layout()plt.show()
The SABR demonstrations make the parameter effects visible. The first panel varies \(\rho\) to change the skew, the second varies \(\nu\) to change curvature, and the third compares \(\beta\) values. The scale \(\alpha\) is adjusted in the beta comparison so that the starting ATM levels remain comparable.
For the skew panel, we compare \(\rho=-0.85\), \(\rho=-0.40\), and \(\rho=0.10\). More negative correlation makes the downside side of the smile steeper. In the curvature panel, we vary \(\nu\) from 0.25 to 1.35 while keeping the other parameters fixed, so we can see the effect of additional volatility-of-volatility.
The last panel compares \(\beta=0.5\), \(\beta=0.7\), and \(\beta=1.0\). We set \(\alpha\) using the same leading ATM volatility for each beta, then compare how the curves differ as strikes move away from the forward. The full surface fit below will estimate \(\rho\) and \(\nu\) for each expiry and compare those beta choices across the observed SPX smiles.
For \(\beta=0.5\), the demo scales \(\alpha\) by \(F^{0.5}\) to preserve the same leading ATM volatility.
Show code
@njitdef sabr_iv_hagan(fwd, strike, tau, alpha, beta, rho, nu): n = strike.size out = np.empty(n, dtype=np.float64)for i inrange(n): f =max(fwd[i], 1e-12) k =max(strike[i], 1e-12) t =max(tau[i], 1e-12) a =max(alpha, 1e-12) b = beta r =min(max(rho, -0.999999), 0.999999) v =max(nu, 1e-12) log_fk = math.log(f / k)ifabs(log_fk) <1e-7: f_beta = f ** (1.0- b) term = ((1.0- b) **2/24.0) * a * a /max(f_beta * f_beta, 1e-12) term +=0.25* r * b * v * a /max(f_beta, 1e-12) term += (2.0-3.0* r * r) * v * v /24.0 out[i] =max(a /max(f_beta, 1e-12) * (1.0+ term * t), 1e-8)else: fk_beta = (f * k) ** ((1.0- b) /2.0) z = (v / a) * fk_beta * log_fk root = math.sqrt(max(1.0-2.0* r * z + z * z, 1e-12)) xz = math.log(max((root + z - r) / (1.0- r), 1e-12)) z_xz =1.0ifabs(z) <1e-8orabs(xz) <1e-10else z / xz log_adj =1.0+ ((1.0- b) **2/24.0) * log_fk * log_fk + ((1.0- b) **4/1920.0) * log_fk **4 term = ((1.0- b) **2/24.0) * a * a /max(fk_beta * fk_beta, 1e-12) term +=0.25* r * b * v * a /max(fk_beta, 1e-12) term += (2.0-3.0* r * r) * v * v /24.0 out[i] =max((a /max(fk_beta * log_adj, 1e-12)) * z_xz * (1.0+ term * t), 1e-8)return outt0 = tic()_ = sabr_iv_hagan(np.array([100.0, 100.0]), np.array([95.0, 105.0]), np.array([0.2, 0.2]), 0.22, 1.0, -0.5, 0.8)sabr_compile_time = toc("compile sabr", t0, 2)
t0 = tic()sabr_runs = []for beta in [1.0, 0.7, 0.5]: sabr_runs.append(fit_sabr(q0, beta))sabr_params = pd.concat([x[0] for x in sabr_runs], ignore_index=True)sabr_fit = pd.concat([x[1] for x in sabr_runs], ignore_index=True)sabr_diag = pd.concat([x[2] for x in sabr_runs], ignore_index=True)toc("calibrate sabr beta scenarios", t0, len(sabr_fit))sabr_beta_compare = sabr_diag.groupby("beta").agg(success_rate=("success", "mean"), median_weighted_iv_rmse=("weighted_rmse_iv", "median"), median_abs_iv_error=("median_abs_iv_error", "median"), rho_bound_share=("rho_bound_flag", "mean"), nu_bound_share=("nu_bound_flag", "mean"), calibration_time=("calibration_time", "sum")).reset_index()best_beta =float(sabr_beta_compare.sort_values(["median_weighted_iv_rmse", "rho_bound_share", "nu_bound_share"]).iloc[0]["beta"])rmse_best =float(sabr_beta_compare["median_weighted_iv_rmse"].min())rmse_beta_one =float(sabr_beta_compare.loc[sabr_beta_compare["beta"].eq(1.0), "median_weighted_iv_rmse"].iloc[0])primary_beta =1.0if rmse_beta_one <=1.08* rmse_best else best_betasabr_fit_main = sabr_fit[sabr_fit["beta"].eq(primary_beta)].copy()sabr_params_main = sabr_params[sabr_params["beta"].eq(primary_beta)].copy()display(sabr_beta_compare)print(f"primary_beta: {primary_beta:.1f}")
beta
success_rate
median_weighted_iv_rmse
median_abs_iv_error
rho_bound_share
nu_bound_share
calibration_time
0
0.500000
1.000000
0.002032
0.001344
0.000000
0.000000
0.468387
1
0.700000
1.000000
0.002160
0.001381
0.000000
0.000000
0.450547
2
1.000000
1.000000
0.002356
0.001437
0.000000
0.000000
0.873524
primary_beta: 0.5
Show code
fig, ax = plt.subplots(figsize=(8.8, 5.0))for expiry in selected_expiries: g = sabr_fit_main[sabr_fit_main["expiry"].eq(expiry)].sort_values("k")if g.empty:continue ax.scatter(g["k"], g["iv_mid"], s=18, alpha=0.50) ax.plot(g["k"], g["model_iv"], lw=2.0, label=f"{g['dte_days'].median():.0f}d")ax.set_title(f"sabr fitted smiles, beta={primary_beta:.1f}")ax.set_xlabel("k")ax.set_ylabel("implied volatility")ax.legend(title="maturity")plt.tight_layout()plt.show()rep_expiry = selected_expiries[0]fig, ax = plt.subplots(figsize=(8.4, 4.4))for beta in [1.0, 0.7, 0.5]: g = sabr_fit[(sabr_fit["expiry"].eq(rep_expiry)) & (sabr_fit["beta"].eq(beta))].sort_values("k") ax.plot(g["k"], g["model_iv"], lw=2.0, label=f"beta={beta:.1f}")market = q0[q0["expiry"].eq(rep_expiry)].sort_values("k")ax.scatter(market["k"], market["iv_mid"], s=18, alpha=0.55, label="market")ax.set_title("sabr beta sensitivity for one expiry")ax.set_xlabel("k")ax.set_ylabel("implied volatility")ax.legend()plt.tight_layout()plt.show()
SABR tracks the main skew well around the liquid center and moderate downside region. The beta-sensitivity plot shows that the selected \(\beta=0.5\) curve is a meaningful choice rather than a random setting.
The residual comparison with SVI shows that the two models disagree most where the smile is hardest to fit: short maturity and wing regions. Disagreement between the flexible static fit and the stochastic-volatility approximation becomes a model-risk signal in the fair-value calculation. If SVI says a quote is cheap but SABR says it is fair, we shouldn’t jump directly into a trade. The disagreement is itself information.
Both residuals use model IV minus market IV on the same quoted contracts.
The beta grid chooses \(\beta=0.5\) as the primary setting. All three beta values succeed, with median weighted IV RMSE around 0.00203 for \(\beta=0.5\), 0.00216 for \(\beta=0.7\), and 0.00236 for \(\beta=1.0\). The displayed results have no bound-share warnings for \(\rho\) or \(\nu\). SABR fits this SPX snapshot well, with more error than the flexible raw SVI slice fits.
Show code
res_xy = svi_fit[["quote_id", "iv_residual"]].rename(columns={"iv_residual": "svi"}).merge(sabr_fit_main[["quote_id","iv_residual"]].rename(columns={"iv_residual": "sabr"}), on="quote_id", how="inner")fig, ax = plt.subplots(figsize=(7.2, 5.0))ax.scatter(res_xy["svi"], res_xy["sabr"], s=18, alpha=0.72)lim =float(np.nanquantile(np.abs(res_xy[["svi", "sabr"]].to_numpy()), 0.98))ax.plot([-lim, lim], [-lim, lim], lw=1.0, alpha=0.70)ax.axhline(0.0, lw=0.8, alpha=0.65)ax.axvline(0.0, lw=0.8, alpha=0.65)ax.set_title("sabr versus raw svi iv residuals")ax.set_xlabel("svi residual")ax.set_ylabel("sabr residual")plt.tight_layout()plt.show()fig, axes = plt.subplots(1, 3, figsize=(13.2, 4.1), sharex=True)for ax, col inzip(axes, ["alpha", "rho", "nu"]): ax.plot(sabr_params_main["dte_days"], sabr_params_main[col], marker="o", lw=1.7) ax.set_title(col) ax.set_xlabel("days to expiry")axes[0].set_ylabel("parameter value")plt.tight_layout()plt.show()
SABR Parameter Term Structures
Each expiry gets its own SABR fit, so we can compare how \(\alpha\), \(\rho\), and \(\nu\) vary with maturity.
A smooth \(\alpha(T)\) curve means the ATM volatility term structure is coherent. A strongly negative \(\rho(T)\) at short maturities means short-dated downside skew is intense. A high \(\nu(T)\) means volatility-of-volatility is needed to explain curvature.
In the output, the term structures are not flat, especially at the short end. That matches SPX reality: short-dated options can reprice sharply around macro events, Fed days, CPI releases, and dealer-flow imbalances. Longer maturities usually smooth some of that out.
This is also why SABR works as a structural cross-check rather than a final truth. If \(\rho\) or \(\nu\) jump around too much between adjacent expiries, the model is fitting local surface features rather than estimating stable dynamics.
Each parameter path comes from independent fits to the quotes available for that expiry.
10) Balanced Price-Space Benchmark Panel
Before fitting Merton, Heston, and Bates, we create a much smaller balanced model panel. The surface models can fit thousands of quotes. The stochastic models are slower, especially Monte Carlo Heston and Bates, so we fit them to a representative grid of maturities and moneyness points.
The target maturities are around:
\[
21,\ 30,\ 45,\ 60,\ 90,\ 120 \text{ days}
\]
and the target \(k\) values include left wing, ATM, and right wing points. The idea is to preserve shape information while keeping the calibration computationally possible.
A random sample would mostly pick quotes from the dense near-ATM region, so we balance the selected panel across moneyness. Then a stochastic-volatility model could fit the center and ignore the tails. Here we intentionally include:
left-wing points,
near-ATM points,
right-wing points,
multiple maturities.
That gives jump and stochastic-vol models a fair chance to identify skew and tail behavior.
Price-Space Calibration for Stochastic Models
Merton, Heston, and Bates generate terminal distributions. Their natural output is an option price:
\[
V^{model}=D\,E^Q[\text{payoff}]
\]
We can invert that price into implied volatility for diagnostics, but the model itself doesn’t start as an implied-volatility formula. That makes price-space calibration more natural.
The tradeoff is that price errors are scale-dependent. A deep ITM option and a tiny OTM option can have very different dollar prices. We can first express the error in units of quote uncertainty:
This standardized residual is the starting point; the actual fit also uses observation weights and the calibration scale described below.
The balanced panel also limits computational cost. A Monte Carlo Heston price has to simulate many paths for every quote and every parameter trial. If we tried to fit Heston to all 3,873 quotes on the main date, the calibration would become too slow for an exploratory project. A representative 54-quote panel gives us enough structure to diagnose the model.
For the actual price calibrations, the denominator is the calibration scale \(s_V=\max(\text{half-spread},0.25\,\mathcal{V}\epsilon_\sigma,0.02)\) defined earlier. The diagnostic quote uncertainty \(\epsilon_V\) uses the full vega-scaled IV uncertainty. Both are in price units, but they serve different roles. We multiply each standardized price residual by \(\sqrt{\omega_i}\) and apply the same soft-\(L_1\) loss used in the smile fits.
The illustration uses \(r_i\) to express a price residual in units of quote uncertainty.
The balanced panel contains six expiries and 54 benchmark quotes. Each expiry has a small number of quotes spread across wings and near-ATM. The selected strikes give the diagnostic stochastic models enough cross-sectional structure to calibrate, while keeping the number of Monte Carlo prices in each optimization trial manageable.
The synthetic jump-distribution plot and sample path visualization introduce the jump idea before we fit Merton. The histogram shows that jumps create a distribution with fatter tails than a pure normal diffusion. The sample path shows discontinuous drops that can’t be generated by a continuous Brownian path. This is exactly the motivation for jump-diffusion option models: short-dated OTM options are often too expensive to be explained by smooth diffusion alone.
Heston and Bates use a smaller Monte Carlo panel capped at 34 selected quotes for their fits.
11) Merton Jump-Diffusion
The Merton jump-diffusion model extends a lognormal diffusion by adding random jumps. In forward form, a simplified risk-neutral log process can be written as:
\(\sigma\) is the continuous diffusion volatility,
\(N_t\) is a Poisson jump counter,
\(\lambda\) is the expected number of jumps per year,
\(J_t\) is the log jump size when a jump occurs,
\(J_t\sim N(\mu_J,\sigma_J^2)\),
\(\kappa_J=E[e^{J_t}-1]\) is the average percentage jump effect.
The compensator \(-\lambda\kappa_J\) is essential. Without it, adding jumps would change the risk-neutral drift of the forward. Under risk-neutral pricing, the discounted forward must be a martingale, so the drift has to offset the average jump effect.
If jumps are usually negative, then \(\mu_J<0\). If jump sizes are uncertain, then \(\sigma_J\) is positive. Both raise the value of OTM puts because the left tail becomes fatter.
For example, if \(\lambda=0.4\), the model expects about 0.4 jumps per year, or one jump every 2.5 years on average. If \(\mu_J=-0.15\), the typical jump is strongly negative in log terms. Even if such jumps are rare, short-dated OTM puts can become much more valuable because they pay exactly in those rare states.
Merton as a Poisson Mixture of Black Prices
Merton’s jump-diffusion model has a useful pricing interpretation. Conditional on the number of jumps \(n\) over the option life, log returns are still normally distributed, but with shifted mean and extra variance from jumps. The probability of \(n\) jumps is Poisson:
So the option price can be interpreted as a weighted sum of Black-like prices over possible jump counts:
\[
V^{Merton}=\sum_{n=0}^{\infty} P(N_T=n)\,V^{Black}(\text{adjusted parameters for } n)
\]
The exact adjustment depends on convention, but the intuition is simple. The no-jump scenario gets the largest probability when \(\lambda T\) is small. One-jump and two-jump scenarios get smaller probabilities, but they can contribute a lot to OTM option values because their payoffs are convex.
This is why rare jumps matter so much for options. A 5% probability of a large downside move can have a large effect on a deep OTM put, even if it barely changes the average forward.
Conditional on \(N_T=n\), the diffusion and \(n\) independent normal log jumps give
A normal log price with mean \(m_n\) and variance \(v_n\) has conditional forward mean \(F_n=e^{m_n+v_n/2}\). Here that becomes \(F_n=F_0\exp[-\lambda\kappa_JT+n\mu_J+\tfrac12 n\sigma_J^2]\), with \(v_n=\sigma^2T+n\sigma_J^2\). We pass \(F_n\) and \(v_n\) into each Black price and weight the result by the Poisson probability. Numerically the infinite sum is truncated at a chosen jump count.
Merton and Short-Maturity Smile Convexity
Jump models are especially important at short maturities. Diffusion moves scale with \(\sqrt{T}\), while jump arrival probability scales with \(T\).
For a diffusion:
\[
\text{typical move} \sim \sigma\sqrt{T}
\]
For a Poisson jump:
\[
P(\text{one jump})\approx \lambda T
\]
At very short maturities, a continuous diffusion needs very high volatility to create deep tail probability. A jump model can put a small amount of probability mass far away from the center. That can make deep OTM options valuable without raising the whole ATM volatility too much.
We use Merton as a tail diagnostic for this reason. If the market’s short-dated left wing is much richer than ATM Black can explain, a negative jump component can bridge that gap. It won’t fit every detail, but it tells us whether discontinuity is an important part of the surface.
The jump probability and diffusion scale both depend on maturity, but affect different parts of the return distribution.
Jump Parameters and Tail Shape
The three jump parameters affect the tail differently.
The intensity \(\lambda\) controls how often jumps arrive. Higher \(\lambda\) raises the probability of one or more jump events:
\[
P(N_T\ge1)=1-e^{-\lambda T}
\]
The mean jump \(\mu_J\) controls the direction and average size. A negative \(\mu_J\) creates downside skew. A positive \(\mu_J\) would make upside calls more expensive.
The jump volatility \(\sigma_J\) controls the dispersion of jump sizes. Larger \(\sigma_J\) makes both tails fatter because jumps can be more extreme in either direction.
For SPX, a negative fitted \(\mu_J\) would be consistent with downside jump risk being priced into the surface. We can compare the estimated mean with the jump volatility to see whether the distribution is mainly one-sided or has substantial dispersion.
The resulting mixture can change both the frequency and the size of the tail outcomes relevant to an option payoff.
Merton Calibration Objective
Merton is calibrated in price space. For selected quotes \(i=1,\ldots,n\), we minimize weighted pricing residuals:
Each residual is divided by the quote’s price uncertainty \(\epsilon_{V,i}\). This matters because a $2 pricing error on a wide-tail option is less serious than a $2 error on a tight, liquid quote.
In this implementation, Merton acts as a jump-tail diagnostic bridge. It is much simpler than Heston and Bates, and it calibrates fast. If it improves tail pricing compared with ATM Black, that tells us jumps are useful for explaining the observed short-maturity wings.
The price scale actually passed to the fit is \(s_V\) from the earlier uncertainty section, so the residual is \(\sqrt{\omega_i}(V_i^{Merton}-V_i^{mkt})/s_{V,i}\). We apply soft-\(L_1\) to that residual and compare several starting vectors within the parameter bounds. The squared expression above shows how the weighting works near zero; the robust loss reduces the influence of larger standardized misses.
Jump Compensation and the Forward
The compensator \(\lambda\kappa_J\) can feel abstract, so here is the financial logic.
If jumps have average multiplicative effect \(E[e^J]\), then:
If jumps are negative on average, \(\kappa_J<0\). If we simply added those jumps without adjusting the drift, the forward would drift downward under the risk-neutral measure. That would violate the martingale condition for a forward.
So the drift includes:
\[
-\lambda\kappa_J
\]
If \(\kappa_J\) is negative, this term is positive. The continuous part drifts slightly upward to offset the average negative jumps. The total risk-neutral process remains correctly centered, while the distribution becomes more negatively skewed and fat-tailed.
The drift adjustment keeps the expected forward anchored while the jump component changes the shape of the payoff distribution. For example, if \(\mu_J=-0.10\) and \(\sigma_J=0.10\), then \(\kappa_J=e^{-0.095}-1\approx-0.091\): the average percentage jump is about \(-9.1\%\), so the compensating drift is positive.
The compensator depends on the mean percentage jump, not just the mean log jump.
Reading the parameters together, continuous volatility is low compared with the surface wings, while the model adds rare negative jumps to explain downside protection. A jump intensity around 0.39 means less than one jump per year on average, but the average log jump is materially negative. That is a natural SPX shape: ordinary daily variation is moderate, but the market prices crash states.
Fit quality is reasonable for this compact price-space panel: weighted price RMSE around 3.33, weighted IV RMSE around 0.0092, median absolute IV error around 0.0044, and P90 around 0.0323. The tail diagnostic is valid and tail error is about 0.0133. The Merton curve improves the tail fit compared with a simple expiry-ATM Black benchmark, where the tail weighted price RMSE is about 16.58 versus 3.25 for Merton.
Both errors are measured on the fitted price-panel observations.
Merton Parameter Interpretation in Annual Units
The fitted jump intensity \(\lambda\approx 0.386\) is annualized. For a 30-day option, the approximate probability of one or more jumps is:
\[
1-e^{-0.386(30/365)}\approx 3.1\%
\]
For a 90-day option:
\[
1-e^{-0.386(90/365)}\approx 9.1\%
\]
These probabilities are small but not negligible. For OTM puts, a few percent probability of a large negative jump can dominate value.
The fitted \(\mu_J\approx -0.146\) corresponds to a typical multiplicative jump around:
\[
e^{-0.146}-1\approx -13.6\%
\]
The jump volatility \(\sigma_J\approx 0.126\) makes jump sizes dispersed around the mean log jump. The \(-13.6\%\) figure is the median percentage jump implied by the mean log jump; the average percentage jump differs because exponentiation is nonlinear. The fitted risk-neutral distribution contains a range of possible jump sizes.
Merton captures the broad left-wing richness much better than the flat ATM-vol comparison. The jump component gives the distribution additional tail mass. It provides a compact way to explain why OTM puts have higher implied volatility.
The residual panels show that Merton is still limited. A single jump distribution across the selected panel can’t fully match every maturity and wing point. It can thicken tails, but it doesn’t create a full stochastic term structure of variance. That limitation is the bridge into Heston: instead of adding only discontinuous moves, we allow variance itself to be random and mean reverting.
The comparison uses the same selected expiries and observed quotes.
Merton Tail Fit and Surface Fit
Merton has only four main parameters here, so it can’t match every maturity-specific smile feature. Its strength is narrower: it can add fat tails.
The model changes the terminal distribution from a single lognormal into a mixture. The left tail gets heavier when jumps are negative or volatile. OTM puts then rise in price.
But Merton lacks stochastic variance. That means it doesn’t naturally create a realistic variance term structure. If the market has different short- and medium-term skew dynamics, Merton has to approximate them through one jump intensity and one jump-size distribution.
That is why Merton’s tail improvement is meaningful but limited. It tells us jump risk is a useful ingredient. It doesn’t mean a pure jump-diffusion is enough for the whole surface.
The same four parameters must account for the selected strikes and maturities simultaneously.
12) Heston: Stochastic Variance as a State Variable
Heston is the central stochastic-volatility model in this project. The risk-neutral dynamics are:
how quickly variance moves back toward long-run variance
\(\theta\)
long-run variance
variance level the process reverts toward
\(\xi\)
vol-of-vol
how violently variance itself moves
\(\rho\)
spot-vol correlation
creates leverage-effect skew
This is a much richer model than Black-Scholes. Volatility isn’t a constant parameter. It is a random process. A low \(v_t\) today doesn’t guarantee low volatility tomorrow because the variance process can receive shocks. The square-root term \(\sqrt{v_t}\) makes variance shocks state-dependent: when variance is high, variance itself is more volatile.
The mean-reversion term \(\kappa(\theta-v_t)\) pulls variance toward \(\theta\):
if \(v_t<\theta\), the drift is positive,
if \(v_t>\theta\), the drift is negative.
This produces volatility clustering without letting variance wander forever.
Here \(v_t\) is variance, so \(\sqrt{v_t}\) is instantaneous volatility. Over a short interval \(dt\), the forward’s conditional return variance is approximately \(v_tdt\). The two Brownian shocks have covariance \(\rho\,dt\), so they can move together even though variance has a separate source of randomness.
We use one parameter vector for the selected strikes and maturities. Each option sees the same variance dynamics over a different remaining lifetime, which links the shape of the smile to the maturity dimension.
The conditional mean gives us a direct way to read the effect of mean reversion. If \(v_t\) is above \(\theta\), the expected future variance decays toward \(\theta\). If \(v_t\) is below \(\theta\), it rises toward \(\theta\). The speed is controlled by \(\kappa\).
For example, if \(\kappa=2\), the half-life of a variance shock is approximately:
So variance shocks are expected to decay over a few months, not instantly. That is why Heston can create term-structure effects: short options see current variance strongly, while longer options average over future mean reversion.
To obtain the conditional mean, take expectations in the variance equation. The stochastic integral has zero expectation under the usual integrability conditions, leaving \(m'(u)=\kappa(\theta-m(u))\) with \(m(0)=v_t\). Solving that ordinary differential equation gives the exponential expression above. The half-life is the horizon at which \(e^{-\kappa u}=1/2\).
Heston Term Structure Through Expected Variance
The conditional mean equation can be integrated to approximate expected average variance over the option life:
This expression explains how Heston creates an implied-volatility term structure. Short maturities are heavily influenced by \(v_0\). Long maturities move closer to \(\theta\) because the variance process has more time to mean revert.
If \(v_0<\theta\), longer options can have higher average variance than very short options. If \(v_0>\theta\), the opposite can happen. When \(v_0\) is much lower than \(\theta\), the expected variance component rises with maturity, while \(\rho\) and \(\xi\) influence the distribution’s skew and curvature.
The term structure comes from the current variance state and the mean-reversion dynamics shared across maturities.
Integrating the conditional mean gives \(E[\int_0^T v_sds]=\theta T+(v_0-\theta)(1-e^{-\kappa T})/\kappa\). Dividing by \(T\) produces the average above. As \(T\to0\), the average tends to \(v_0\); as \(T\) grows, the mean-reversion contribution tends toward \(\theta\). This is an expected variance calculation, not a direct formula for the Black implied volatility of every strike.
Heston and the Leverage Effect in Covariance Form
The spot-vol correlation can also be read through covariance. From the Heston equations:
\[
d\log F_t \approx \sqrt{v_t}dW_t^F
\]
and:
\[
dv_t \approx \xi\sqrt{v_t}dW_t^v
\]
The instantaneous covariance is:
\[
\text{Cov}(d\log F_t,dv_t)=\xi v_t \rho\,dt
\]
If \(\rho<0\), this covariance is negative. That means negative returns tend to come with positive variance shocks. The strength of the effect increases with \(\xi\) and with current variance \(v_t\).
The parameters \(\rho\) and \(\xi\) therefore work together in the covariance. A mildly negative \(\rho\) with small \(\xi\) produces weak skew. A strongly negative \(\rho\) with high \(\xi\) produces strong skew. A fit with both strongly negative \(\rho\) and high \(\xi\) uses a stronger leverage channel to explain the observed downside skew.
How Heston Creates Skew
The correlation \(\rho\) is the skew engine. If \(\rho<0\), negative forward shocks tend to occur with positive variance shocks. In equity index markets, that matches the leverage effect: when equities fall, volatility often rises.
If \(\rho=-0.7\) and the forward shock \(dW_t^F\) is negative, then \(dW_t^v\) tends to be positive. That means the price path falls while future variance rises. A put option benefits from both effects: the underlying moves toward the strike and volatility increases. This is why negative \(\rho\) produces downside skew.
We can examine both effects with simulated paths. Variance paths mean-revert but fluctuate because variance has its own Brownian shock. A correlated forward/variance simulation lets us see the forward and instantaneous volatility moving together through time.
These paths use correlated shocks at each time step, so the relationship can be seen along the same simulated trajectory.
Integrated Variance and Option Prices
An option payoff depends on the terminal forward, but under stochastic volatility the terminal distribution is driven by the whole path of variance. The key path quantity is integrated variance:
\[
I_T=\int_0^T v_s ds
\]
The path of \(v_t\) determines the accumulated quadratic variation \(I_T\). Options therefore depend on how variance evolves through the full holding period, as well as on the starting variance \(v_0\).
This is why Heston can create smiles. Random integrated variance changes the shape of the final return distribution. When the return and variance shocks are independent, we can view that distribution as a mixture of conditional normal distributions with different variances. Negative correlation \(\rho\) then makes the left tail even heavier because bad forward moves tend to occur with higher future variance.
A simple mental picture:
Black-76 uses one variance bucket.
Heston randomizes the variance bucket.
Negative \(\rho\) links the bucket size to the direction of the move.
That combination produces skew and smile while keeping the path continuous.
The integrated variance controls the quadratic variation of the continuous log-forward, \(\langle\log F\rangle_T=I_T\). When \(\rho\ne0\), the forward shock is correlated with the variance path, so conditioning on the full variance path does not generally leave a centered normal return with variance exactly \(I_T\). The correlation also shifts the conditional return distribution. The pathwise quantity is still useful for understanding how changing variance accumulates through the option’s life.
The Feller Condition
The Heston variance process is a square-root process. A classic condition for strict positivity is:
If this ratio is above 1, the continuous-time process stays strictly positive under the classical condition. If it is below 1, variance can approach zero in theory, and numerical schemes need care.
The Feller ratio is one diagnostic to consider alongside the calibration fit. Real market calibrations often violate it because the market wants high vol-of-vol \(\xi\) to fit steep smiles. A low Feller ratio says: the parameter set is using very aggressive variance randomness relative to mean reversion and long-run variance. We can still use it diagnostically, but we have to read it as a warning.
Monte Carlo Pricing Under Heston
Although Heston has a semi-closed-form characteristic-function solution, here we price it by Monte Carlo. For every path:
start at the current forward \(F_0\) and variance \(v_0\),
The use of \(\max(v_t,0)\) is a practical full-truncation-style protection. It prevents the square root from breaking when the Euler step makes variance slightly negative.
This standard error is central to reading the residuals. If the model price is $0.40 away from the market but Monte Carlo standard error is $0.45, the residual is barely meaningful.
The correlated normals can be constructed from independent standard normals \(Z_1,Z_2\) as \(Z_F=Z_1\) and \(Z_v=\rho Z_1+\sqrt{1-\rho^2}Z_2\). Both have unit variance and their correlation is \(\rho\). The fit uses 1,024 paths per trial; final prices use 8,192 paths, with 52 time steps per year before rounding each maturity to an integer number of steps.
The displayed standard-error formula is the payoff-variance estimate used by the pricing routine. Since the paths include antithetic pairs, it is a useful diagnostic of sampling scale rather than a separate estimate of the covariance within each pair.
Common Random Numbers and Calibration Noise
Monte Carlo calibration can become unstable if every parameter trial uses new random numbers. The objective function would be noisy: a parameter vector might look better simply because it got luckier simulated paths.
To reduce this, we use common random numbers. The same random shock arrays are reused across parameter evaluations. Then the difference between two parameter vectors is mostly due to the parameters, not due to fresh simulation noise.
If the objective is:
\[
L(\theta)=\sum_i r_i(\theta)^2
\]
and \(r_i(\theta)\) is estimated by Monte Carlo, common random numbers make the estimated loss surface smoother. That helps least-squares optimization converge.
We also use antithetic shocks. If \(Z\) is a normal shock, we also simulate \(-Z\). Averaging the pair reduces simulation variance because positive and negative random fluctuations partially cancel. This doesn’t remove Monte Carlo error, but it makes the same number of paths more efficient.
If prices under two parameter vectors are \(\widehat V(\theta_1)\) and \(\widehat V(\theta_2)\), their difference has variance \(\operatorname{Var}(\widehat V_1)+\operatorname{Var}(\widehat V_2)-2\operatorname{Cov}(\widehat V_1,\widehat V_2)\). Reusing shocks generally makes that covariance positive for nearby parameter values, so the difference can be much less noisy than two independently simulated estimates. We use a separate shock set for final pricing after fitting.
Optimization Noise and False Parameter Precision
Monte Carlo calibrated parameters can look more precise than they really are. The optimizer returns a parameter vector, but that vector depends on:
selected quotes,
parameter bounds,
random path set,
number of paths,
time discretization,
objective scaling,
starting values.
If the MC standard error is large relative to pricing residuals, then small differences in objective value are not meaningful. Two parameter sets can produce prices that differ by less than simulation noise.
We use Heston and Bates as diagnostic models in this comparison. The parameter values are useful for economic interpretation, but we shouldn’t treat the last decimal as meaningful. For example, a value such as \(\rho=-0.6286\) describes a strongly negative spot-vol correlation; the final decimal should be read in the context of the simulation and calibration precision.
We report Monte Carlo standard errors alongside the fitted prices so the size of sampling noise is visible next to the pricing residuals.
The initial volatility implied by \(v_0\) is roughly:
\[
\sqrt{0.0117}\approx 10.8\%
\]
The long-run volatility implied by \(\theta\) is roughly:
\[
\sqrt{0.0836}\approx 28.9\%
\]
So the model says current instantaneous variance is much lower than the long-run variance level needed to explain the option panel. That sounds extreme, but in a short-dated option surface it can happen: the model may use high long-run variance and high vol-of-vol to create skew and convexity even if the current variance starts low.
The Feller ratio is only about 0.212, so the Feller warning is active. This means the fitted process has very high vol-of-vol relative to mean reversion and long-run variance. It is trying hard to fit the smile dynamics with variance randomness.
The variance parameters are annualized; their square roots are volatility levels.
Heston’s Economic Reading on This SPX Date
The Heston fit implies a market that is calm in current instantaneous variance but still prices substantial future variance risk. That combination often appears when spot markets are calm but option markets still price a nontrivial probability of volatility expansion.
The negative \(\rho\) tells us the model expects bad equity states and high variance states to arrive together. The large \(\xi\) tells us the variance state can move sharply. The low Feller ratio tells us those variance shocks are aggressive enough to push the square-root process near zero in the numerical approximation.
The fitted surface uses a strong leverage channel and volatile variance dynamics. Its residual diagnostics also show how much cross-sectional error remains compared with the SVI and SABR fits.
The parameter values and the residuals should be read together when comparing the continuous-variance model with the surface fits.
Reading the Heston Parameters Together
We can read the parameters together to understand how the fitted variance level and leverage effect combine.
means the model starts with low instantaneous variance but expects much higher long-run variance. The combination of:
\[
\xi\approx 1.25,\quad \rho\approx -0.63
\]
means variance shocks are large and strongly negatively correlated with the forward. That creates a steep left skew.
The Feller ratio:
\[
\frac{2\kappa\theta}{\xi^2}\approx 0.21
\]
says the variance process is operating in a high-vol-of-vol regime. In calibration language, the model is stretching the variance process to reproduce option smile convexity. That is a warning about parameter stability, especially if we wanted to use these parameters for long-horizon simulation or risk scenarios.
Discretization Bias and Full Truncation
Euler simulation of square-root variance can produce negative variance if the random shock is large. The continuous-time model is designed for nonnegative variance, but the discrete approximation can step below zero.
We apply a practical numerical protection:
\[
v^+=\max(v,0)
\]
inside square roots and drift terms. Then the simulation update uses \(v^+\) rather than raw \(v\).
This keeps the square root numerically valid while leaving a discretization approximation in the simulated paths. With more time steps per year, the approximation improves, but computation becomes slower. This is another reason Monte Carlo Heston is diagnostic here. The model has both parameter uncertainty and numerical uncertainty.
The MC convergence plot helps quantify one side of this. It shows sampling error as paths increase. It doesn’t fully remove discretization bias, but it tells us whether path count is large enough for price noise to be acceptable.
Fit quality is weaker than the best surface models: weighted price RMSE around 3.86, weighted IV RMSE around 0.0149, median absolute IV error around 0.0096, and P90 around 0.0392. The median Monte Carlo standard error is about 0.455 in price units.
The comparison reflects the different restrictions placed on the surface. Heston is theoretically richer than SVI, but it doesn’t automatically fit a single-day quote surface better. SVI is a flexible interpolating smile form. Heston is a low-dimensional stochastic process forced to explain many strikes and maturities through a single variance dynamic. That restriction is valuable for dynamics, but it costs cross-sectional precision.
As the number of Monte Carlo paths increases, the estimated price and standard error change. The standard error falls as expected with approximately:
\[
SE \propto \frac{1}{\sqrt{M}}
\]
where \(M\) is the number of paths. Doubling paths doesn’t halve the error; it only reduces it by about \(\sqrt{2}\). That is why Monte Carlo stochastic-volatility calibration is expensive.
The Heston smile plot captures the broad negative skew but misses some detailed curvature. The residual heatmap shows systematic regions of underfit/overfit rather than pure random noise. The parameter diagnostic highlights the Feller warning, which matches the intuition that the fitted model is relying heavily on aggressive variance shocks.
Show code
fig, ax = plt.subplots(figsize=(8.8, 4.8))for expiry in heston_fit["expiry"].drop_duplicates().head(5): g = heston_fit[heston_fit["expiry"].eq(expiry)].sort_values("k")iflen(g) <3:continue ax.scatter(g["k"], g["iv_mid"], s=22, alpha=0.52) ax.plot(g["k"], g["model_iv"], marker="o", lw=1.6, label=f"{g['dte_days'].median():.0f}d")ax.set_title("heston simulation model iv")ax.set_xlabel("k")ax.set_ylabel("implied volatility")ax.legend(title="maturity", fontsize=8)plt.tight_layout()plt.show()z = heston_fit.copy()z["k_bin"] = pd.cut(z["k"], np.linspace(-0.26, 0.26, 11))z["dte_bin"] = pd.cut(z["dte_days"], np.linspace(21.0, 135.0, 8))h = z.pivot_table(index="dte_bin", columns="k_bin", values="iv_residual", aggfunc="median", observed=False)lim =float(np.nanquantile(np.abs(h.to_numpy(dtype=float)), 0.90)) if np.isfinite(h.to_numpy(dtype=float)).any() else0.03fig, ax = plt.subplots(figsize=(8.8, 4.6))im = ax.imshow(h.to_numpy(dtype=float), origin="lower", aspect="auto", cmap="coolwarm", vmin=-lim, vmax=lim)ax.set_title("heston residual heatmap")ax.set_xlabel("k bucket")ax.set_ylabel("maturity bucket")ax.set_xticks(np.arange(h.shape[1]))ax.set_xticklabels([f"{b.mid:.2f}"for b in h.columns], rotation=45, ha="right")ax.set_yticks(np.arange(h.shape[0]))ax.set_yticklabels([f"{b.mid:.0f}"for b in h.index])fig.colorbar(im, ax=ax, pad=0.01, label="model iv - market iv")plt.tight_layout()plt.show()
Show code
fig, axes = plt.subplots(1, 2, figsize=(11.6, 4.3))axes[0].scatter(np.abs(heston_fit["price_residual"]), heston_fit["mc_standard_error"], s=24, alpha=0.75)axes[0].set_title("heston mc error versus residual size")axes[0].set_xlabel("absolute price residual")axes[0].set_ylabel("mc standard error")p_plot = heston_params[["v0", "kappa", "theta", "xi", "rho"]].iloc[0]axes[1].bar(p_plot.index, p_plot.values, alpha=0.82)axes[1].set_title("heston parameter diagnostic")axes[1].tick_params(axis="x", rotation=45)plt.tight_layout()plt.show()
13) Bates: Heston Plus Jumps
Bates combines Heston stochastic variance with Merton-style jumps. The log-forward update adds a jump term:
continuous stochastic variance, where volatility rises and creates fat tails gradually,
jumps, where the underlying can move discontinuously.
The model can use Heston variance to explain persistent skew and term-structure effects, while using jumps to handle short-dated tail prices. That is exactly why Bates is often a useful extension: short-dated OTM options can be too expensive for pure continuous stochastic volatility, while long-dated smile dynamics need more than a one-shot jump distribution.
The jump compensator uses the same \(\kappa_J=e^{\mu_J+\sigma_J^2/2}-1\) as Merton. In each numerical step we simulate the correlated forward and variance shocks and allow a jump with probability approximately \(\lambda\Delta t\), capped at 0.35 by the implementation. This one-jump-per-step approximation is part of the simulation; the continuous-time model above has a Poisson jump count.
This nesting is useful because the calibration can decide whether the surface needs jumps, stochastic variance, or both. But nesting also creates identification risk. If a left-wing price can be explained by either high jump intensity or very negative spot-vol correlation, the optimizer may find several parameter sets with similar prices.
We start Bates from the Heston and Merton fits and add regularization to the calibration. We want Bates to improve the fit by adding an economically meaningful jump component, not by randomly redistributing error across too many parameters.
The fitted starts put both channels on the same price-panel scale.
Bates as a Diagnostic Model
Bates improves the compact stochastic-panel fit because it has two mechanisms. But better compact-panel fit doesn’t automatically make it better for daily signals.
There are three reasons:
Monte Carlo noise: each Bates price has simulation error.
Parameter identification: jumps and stochastic variance can substitute for each other.
Daily runtime: fitting Bates repeatedly is expensive.
So Bates is best used here as a diagnostic model. If Bates improves tail fit relative to Heston, we learn that jumps are useful. If Bates and Merton agree that jump risk matters, we learn the tail is more discontinuity-driven. But the actual fair-value engine still leans on SVI/SABR because those models are faster, cleaner, and easier to validate daily.
We use the richer process to study the surface and the faster smile models to construct and validate repeated daily signals.
The daily fair-value fit therefore uses the faster surface models while the jump and variance processes remain available for comparison.
More parameters mean more flexibility, but also more calibration risk. The implementation starts from the fitted Heston parameters and fitted Merton jump parameters. It also includes small regularization terms that discourage Bates from moving too far from those starts unless the pricing residuals justify it.
Heston variance and jumps can explain similar parts of the smile. For example, left-tail put richness can come from:
very negative \(\rho\) and high \(\xi\) in Heston,
negative jumps with positive intensity in Merton/Bates,
or both.
When two mechanisms explain the same quotes, the calibration becomes weakly identified. Regularization prevents the optimizer from producing extreme parameter combinations that fit the panel but have poor economic stability.
The penalty enters as additional residuals in the same soft-\(L_1\) fit. We anchor \(v_0\) and \(\xi\) to the Heston fit and \(\lambda,\mu_J,\sigma_J\) to the Merton fit, scaling each difference before applying its small coefficient. The optimizer can move away from those starts when the improvement in the price residuals is large enough; \(\kappa\), \(\theta\), and \(\rho\) are still fitted within their bounds.
Jump and Volatility Channels in Bates
Bates gives us two different channels for tail pricing.
The continuous stochastic-volatility channel affects many maturities:
The jump channel affects short-dated wings especially strongly:
\[
N_T \sim \text{Poisson}(\lambda T)
\]
For small \(T\), the probability of one jump is approximately:
\[
P(N_T=1)\approx \lambda T
\]
So a one-month option with \(\lambda=0.46\) has a one-jump probability around:
\[
0.46\times\frac{1}{12}\approx 3.8\%
\]
That may sound small, but for an OTM put the payoff is concentrated exactly in those rare bad states. A small jump probability can materially raise the option price.
These two channels give Bates additional flexibility in short-dated tail pricing. Heston can generate fat tails through random variance, but jumps can put probability mass farther into the tail more directly.
Compared with Heston, Bates has lower long-run variance and slightly less negative \(\rho\), while the jump component takes part of the tail burden. The jump mean is less negative than in pure Merton, and jump volatility is smaller. This makes sense: once stochastic variance is present, jumps don’t need to explain the entire left wing by themselves.
Fit quality improves materially versus Heston: weighted price RMSE falls to about 2.09, weighted IV RMSE falls to about 0.0090, median absolute IV error is about 0.0053, and P90 IV error is about 0.0217. It is still less accurate than SVI on the full surface, but for the compact stochastic panel Bates is clearly stronger than Heston.
Bates Parameter Reading on This SPX Date
The Bates jump intensity \(\lambda\approx 0.459\) is slightly higher than Merton’s fitted jump intensity, but the average jump size is much less negative. The fitted jump component changes alongside the variance parameters.
Pure Merton had to use larger negative jumps to explain the left tail because it had no stochastic variance. Bates can let stochastic variance handle part of the skew, so jumps become smaller and more frequent/moderate. That is a more balanced decomposition:
The fitted \(\rho\approx -0.571\) remains strongly negative, so the leverage effect is still active. The jump mean \(\mu_J\approx -0.0396\) corresponds to a typical jump around:
\[
e^{-0.0396}-1\approx -3.9\%
\]
with jump volatility about 5.9%. That is a smaller jump distribution than pure Merton, because Bates no longer needs jumps to explain everything.
For Bates, as for Heston, more paths reduce standard error, but slowly. The Heston-vs-Bates residual plot shows Bates reducing the miss in several regions, especially where jump tails help.
The Bates residual and parameter diagnostics add to the Heston comparison. However, the model still carries Monte Carlo noise, and the tail diagnostic is not fully valid in the displayed table because only a few tail quotes are available in the balanced stochastic panel. That means Bates is better as a diagnostic model here than as the main signal generator.
The tail statistic needs sufficient tail observations before it is reported as valid.
Bates Residuals Versus Heston Residuals
Comparing the Heston and Bates residuals helps us see where jumps improve the fit. If Bates improves the same quotes where Heston struggles, then jumps are adding useful explanatory power. If Bates improves only through unstable parameters or huge jump intensity, then the improvement may be overfit.
Here Bates improves weighted price and IV errors without hitting parameter bounds. The jump intensity is positive but not extreme, and jump mean/volatility are moderate. That suggests the jump component is being used as a controlled tail adjustment rather than as a wild parameter patch.
Still, Bates remains a compact-panel diagnostic. It isn’t used as the main daily fair-value engine because MC noise and identification risk are too high compared with SVI/SABR.
The comparison holds the compact pricing panel fixed.
14) Comparing Models Across Accuracy, Runtime, and Reliability
We now compare all six models across their fitted errors, runtime, and parameter diagnostics.
The fit objective isn’t one-dimensional. A model can be:
accurate but too flexible,
stable but biased,
structural but too slow,
theoretically attractive but weak on the current surface,
useful for diagnostics but risky for trading signals.
The main comparison metrics are:
\[
\text{weighted IV RMSE}=\sqrt{\frac{\sum_i \omega_i(\sigma_i^{model}-\sigma_i^{mkt})^2}{\sum_i\omega_i}}
\]
The goal is to choose which models are safe enough for signals instead of relying only on the lowest table error.
The RMSE is in the units of the quantity being compared: an IV RMSE of 0.01 is one volatility point, while price RMSE is in option price units. The own-panel results use each model’s own calibration quotes; the common benchmark evaluates the fitted models on the same selected quotes. Reading both prevents a model’s score from being confused with a change in the number or location of quotes.
For the common benchmark, the current Heston and Bates simulations need a numerical correction before their longer-maturity errors can be relied on. The shock arrays are sized from the longest quote in the smaller Monte Carlo calibration panel, while the common panel includes longer expiries. Those longer quotes can require more simulation steps than the arrays contain. The compact-panel calibration results use maturities within the allocated range.
A Model Ranking Needs Several Axes
A useful model ranking needs at least four axes:
Axis
Question
Examples from this project
Accuracy
Does it fit market quotes?
SVI has the best IV RMSE
Structure
Does it encode sensible dynamics?
Heston/Bates have variance processes
Stability
Are parameters away from warnings and bounds?
Heston has a Feller warning
Speed
Can we fit it repeatedly?
SSVI/Merton are fast, MC models are slow
Each of these axes describes a different part of model performance. SVI is extremely accurate but doesn’t describe a forward variance process. Heston has a true stochastic variance story but fits worse and has Monte Carlo noise. Bates improves stochastic-panel fit but adds identification risk. SABR is fast and structural enough to act as a cross-check.
This multi-axis view is what allows model-risk-aware signal design later. A residual confirmed by SVI and SABR is stronger than a residual generated by a single fragile model.
The common benchmark keeps the quote set fixed when comparing accuracy.
If several credible models produce similar fair values, the residual is more trustworthy. If models disagree widely, then the quote may be in a region where assumptions matter more than data.
In the wings, different model assumptions can lead to materially different fair values. A deep OTM put can look cheap under SVI but fair under Merton if jump risk is priced differently. Or it can look rich under SSVI because the global model is too rigid. Disagreement tells us to reduce confidence.
We can compare the median cross-model IV disagreement with the median quote IV uncertainty on the same volatility scale. If the two are similar in size, model choice and quote precision both contribute meaningfully to the error budget.
Show code
def svi_iv_at(q, p_table): vals = []for _, row in q.iterrows(): p = p_table[p_table["expiry"].eq(row["expiry"])]if p.empty: vals.append(np.nan)else: pp = p.iloc[0] vals.append(float(svi_w(np.array([row["k"]], dtype=float), pp["a"], pp["b"], pp["rho"], pp["m"], pp["sigma"])[0]))return np.sqrt(np.maximum(np.asarray(vals, dtype=float), 1e-12) / q["tau"].to_numpy(dtype=float))def ssvi_iv_at(q, p_frame, theta_table): x = q.merge(theta_table[["expiry", "theta_ssvi"]], on="expiry", how="left") p = p_frame.iloc[0] w = ssvi_w(x["k"].to_numpy(dtype=float), x["theta_ssvi"].to_numpy(dtype=float), p["rho"], p["eta"], p["gamma"])return np.sqrt(np.maximum(w, 1e-12) / x["tau"].to_numpy(dtype=float))def sabr_iv_at(q, p_table, beta): vals = np.full(len(q), np.nan, dtype=float)for expiry, g in q.groupby("expiry"): p = p_table[p_table["expiry"].eq(expiry)]if p.empty:continue pp = p.iloc[0] loc = q.index.get_indexer(g.index) vals[loc] = sabr_iv_hagan(g["forward"].to_numpy(dtype=float), g["strike"].to_numpy(dtype=float), g["tau"].to_numpy(dtype=float), pp["alpha"], beta, pp["rho"], pp["nu"])return vals
The fitted model errors give us the following comparison:
SVI dominates raw fit: weighted IV RMSE around 0.00076, median absolute IV error about 0.00034, and P90 about 0.00102.
SABR is second among fast smile models: weighted IV RMSE around 0.00250.
SSVI is much weaker on quote-level fit: weighted IV RMSE around 0.00832.
Bates is the best among the stochastic price-space models on the compact panel, with weighted IV RMSE around 0.00905.
Merton is close to Bates in IV RMSE on the compact panel, around 0.00921, but it lacks stochastic variance.
Heston is weaker here, around 0.01487, with MC and Feller warnings.
On the common benchmark of 54 quotes, SVI and SABR remain strong relative to SSVI. The longer-maturity Heston and Bates benchmark values need the shock-array sizing correction described above before we use them as a reliable comparison; the compact-panel results are still useful for reading their fitted dynamics.
The comparison also gives each model a role in the rest of the analysis:
SVI is the fair-value core.
SABR is a structural cross-check.
SSVI is a surface discipline check.
Merton, Heston, and Bates are diagnostic engines for jumps and stochastic volatility.
These roles follow the diagnostics for the selected SPX surface and the common benchmark.
Calibration speed gives us another dimension for the comparison:
SSVI calibrates in about 0.056 seconds after compilation.
Merton calibrates in about 0.051 seconds on the compact panel.
SABR takes about 0.47 seconds.
SVI takes about 7.64 seconds across 3,873 quotes and 29 slices.
Heston and Bates take around 3.4–3.7 seconds each on only 34 quotes, and only about 9–10 quotes per second because of Monte Carlo.
This matters for real workflows. If we want to fit models daily across hundreds of dates, SVI and SABR are feasible. Heston and Bates can be used, but not carelessly across millions of quotes unless we move to faster Fourier pricing, better vectorization, or fewer calibration targets.
Accuracy Per Second
A practical way to read the speed table is accuracy per unit of compute.
SVI is slower than SSVI and SABR on this date because it fits many expiry slices, but it processes the full 3,873-quote surface and produces the best accuracy. SSVI is extremely fast, but its error is much larger. SABR is a good middle ground: fast and reasonably accurate.
Heston and Bates are expensive because every objective evaluation runs thousands of simulated paths across the compact quote set. The calibration time looks acceptable for one date, but daily historical use would scale quickly. For 500 dates, even 3.5 seconds per calibration becomes roughly half an hour per model before any tuning, plotting, validation, or reruns.
The choice of pricing method changes how often we can run each model. We use MC stochastic models where they add diagnostic value, and we use fast surface models where repeated daily fair values are needed.
svi_err =float(model_comparison.loc[model_comparison["model"].eq("svi"), "weighted_iv_rmse"].iloc[0])ssvi_err =float(model_comparison.loc[model_comparison["model"].eq("ssvi"), "weighted_iv_rmse"].iloc[0])sabr_err =float(model_comparison.loc[model_comparison["model"].eq("sabr"), "weighted_iv_rmse"].iloc[0])heston_err =float(model_comparison.loc[model_comparison["model"].eq("heston"), "weighted_iv_rmse"].iloc[0])sabr_diag_main = sabr_diag[sabr_diag["beta"].eq(primary_beta)].copy()sabr_success_rate =float(sabr_diag_main["success"].mean())rho_bound_share =float(sabr_diag_main["rho_bound_flag"].mean())nu_bound_share =float(sabr_diag_main["nu_bound_flag"].mean())ssvi_ok_for_trade =bool((ssvi_err <=0.004) or (ssvi_err <=3.0* sabr_err))sabr_ok_for_trade =bool((sabr_err <=0.0025) and (sabr_success_rate >0.85) and (rho_bound_share <0.20) and (nu_bound_share <0.20))sabr_in_signal =bool((sabr_success_rate >0.85) and (rho_bound_share <0.20) and (nu_bound_share <0.20))heston_ok_for_trade =bool( heston_diag.loc[0, "success"]and heston_err <=0.0060andnotbool(heston_params.loc[0, "parameter_at_bound"])andfloat(heston_fit["mc_standard_error"].median()) <0.50*float(np.nanmedian(np.abs(heston_fit["price_residual"]))))model_decision = pd.DataFrame([ {"model": "svi", "signal_decision": "candidate", "reason": "holdout signal core"}, {"model": "ssvi", "signal_decision": "candidate"if ssvi_ok_for_trade else"diagnostic", "reason": "surface discipline model"}, {"model": "sabr","signal_decision": "candidate"if sabr_ok_for_trade else"diagnostic"ifnot sabr_in_signal else"candidate with fit warning","reason": f"abs err={sabr_err:.4f}, success={sabr_success_rate:.2f}"}, {"model": "heston", "signal_decision": "candidate"if heston_ok_for_trade else"diagnostic","reason": "mc se, bounds, and residual checks"}, {"model": "merton", "signal_decision": "no signal", "reason": "diagnostic jump bridge only"}, {"model": "bates", "signal_decision": "no signal", "reason": "diagnostic stochastic-vol jump model only"},])display(model_decision)
model
signal_decision
reason
0
svi
candidate
holdout signal core
1
ssvi
candidate
surface discipline model
2
sabr
candidate with fit warning
abs err=0.0027, success=1.00
3
heston
diagnostic
mc se, bounds, and residual checks
4
merton
no signal
diagnostic jump bridge only
5
bates
no signal
diagnostic stochastic-vol jump model only
We use the fit and warning diagnostics to decide which models can contribute to signals. SVI and SSVI are candidates, SABR is a candidate with a fit warning, Heston is diagnostic only, and Merton/Bates are no-signal diagnostic models. Each label reflects the model’s accuracy and the warnings recorded for this fit.
The uncertainty summary shows:
median quote IV uncertainty around 0.0025,
median cross-model IV disagreement around 0.00262,
Heston MC standard error around 0.455,
Bates MC standard error around 0.474,
bootstrap price uncertainty later around 0.283.
Those numbers tell us that many small residuals aren’t tradable. If the model says an option is mispriced by $0.25, but model uncertainty and bid-ask costs are larger than that, the signal is weak.
Candidate, Diagnostic, and No-Signal Labels
We assign models practical roles based on their fit and diagnostics.
A candidate model is accurate and stable enough to contribute to fair values.
A candidate with warning can contribute, but with lower trust.
A diagnostic model teaches us about structure but shouldn’t directly drive signals.
A no-signal model is useful for comparison but too noisy, too narrow, or too risky to trade from.
This is a better design than forcing all models into one ensemble. If Heston has a Monte Carlo warning and a Feller warning, giving it a fair-value weight could inject noise into the residual score. If Merton is only a jump bridge, it can help interpret tails without becoming a daily trading model.
The resulting ensemble draws on the models selected for repeated fair-value estimates.
Show code
plot_q = model_comparison.copy()if"calibration_time"notin plot_q.columns: plot_q = plot_q.merge(speed_summary[["model", "calibration_time"]], on="model", how="left")else: time_map = speed_summary.set_index("model")["calibration_time"] plot_q["calibration_time"] = plot_q["calibration_time"].fillna(plot_q["model"].map(time_map))fig, axes = plt.subplots(1, 3, figsize=(14.2, 4.2))axes[0].bar(plot_q["model"], plot_q["weighted_iv_rmse"], alpha=0.82)axes[0].set_title("weighted iv error by model")axes[0].set_ylabel("weighted iv rmse")axes[0].tick_params(axis="x", rotation=45)axes[1].bar(plot_q["model"], plot_q["calibration_time"], alpha=0.82)axes[1].set_title("calibration time by model")axes[1].set_ylabel("seconds")axes[1].tick_params(axis="x", rotation=45)axes[2].scatter(plot_q["calibration_time"], plot_q["weighted_iv_rmse"], s=54, alpha=0.82)for _, row in plot_q.iterrows(): axes[2].annotate(row["model"], (row["calibration_time"], row["weighted_iv_rmse"]), fontsize=8)axes[2].set_title("error versus runtime")axes[2].set_xlabel("calibration seconds")axes[2].set_ylabel("weighted iv rmse")plt.tight_layout()plt.show()
Show code
z = fit_iv.merge(wide_iv[["model_iv_std"]], left_on="quote_id", right_index=True, how="left")z["k_bin"] = pd.cut(z["k"], np.linspace(-max_abs_k, max_abs_k, 15))z["dte_bin"] = pd.cut(z["dte_days"], np.linspace(min_dte, max_dte, 11))h = z.drop_duplicates("quote_id").pivot_table(index="dte_bin", columns="k_bin", values="model_iv_std", aggfunc="median", observed=False)fig, ax = plt.subplots(figsize=(9.8, 4.8))im = ax.imshow(h.to_numpy(dtype=float), origin="lower", aspect="auto", cmap="magma")ax.set_title("cross-model iv disagreement heatmap")ax.set_xlabel("k bucket")ax.set_ylabel("maturity bucket")ax.set_xticks(np.arange(h.shape[1]))ax.set_xticklabels([f"{b.mid:.2f}"for b in h.columns], rotation=45, ha="right")ax.set_yticks(np.arange(h.shape[0]))ax.set_yticklabels([f"{b.mid:.0f}"for b in h.index])fig.colorbar(im, ax=ax, pad=0.01, label="median iv std")plt.tight_layout()plt.show()
15) From Model Fit to Fair Value
We now turn fitted models into fair values. This is the step where a pricing model becomes a possible trading signal.
For each quote, the residual is:
\[
e_i=V_i^{model}-V_i^{mkt}
\]
The sign tells us the direction:
\(e_i>0\): the model fair value is above the market mid, so the option looks cheap.
\(e_i<0\): the model fair value is below the market mid, so the option looks rich.
But a raw residual is not enough. We scale it by total uncertainty:
\[
z_i=\frac{e_i}{\epsilon_i}
\]
where \(\epsilon_i\) combines bid-ask spread, model disagreement, and fit/hedge buffers. A larger \(|z_i|\) means the residual is large relative to uncertainty.
The ensemble uses SVI as the main core, with weight rules for the eligible models:
SVI has a maximum weight of 0.70,
SABR has a minimum weight of 0.20 when it is included,
SSVI has a maximum weight of 0.20 when it is included.
The rules retain a contribution from SABR and limit how much SSVI can influence the ensemble. The fitted errors determine the starting weights before those limits are applied.
Before the caps, the ensemble starts with inverse squared benchmark IV errors, \(\tilde w_m\propto1/\mathrm{RMSE}_m^2\), and normalizes the weights to sum to one. We then apply the SVI and SSVI caps and the SABR minimum when the corresponding models are eligible. The fair price is the weighted sum \(\widehat V_i=\sum_m w_mV_i^{(m)}\).
Total Error Used for Residual Signals
A raw residual becomes a standardized residual only after we divide by total error:
\(\epsilon_i^{quote}\) comes from bid-ask and vega uncertainty,
\(\epsilon_i^{model}\) comes from disagreement between models,
\(\epsilon_i^{fit}\) comes from calibration fit uncertainty,
\(\epsilon_i^{hedge}\) reflects the fact that the trade will be hedged discretely and with costs.
The residual is compared with the combined price-unit error budget. A model residual has to beat several sources of uncertainty before it becomes a candidate. That is why the candidate table later filters out most quotes.
For example, suppose an option has model residual \(0.60\), entry half-spread \(0.25\), expected exit half-spread \(0.20\), model disagreement \(0.30\), fit error \(0.10\), and hedge buffer \(0.15\). With no MC error, the weighted components add to \(0.25+0.5(0.20)+0.5(0.30)+0.5(0.10)+0.15\), giving total error \(0.70\), so:
\[
z=\frac{0.60}{0.70}\approx 0.86
\]
That is a weak signal. Even though the raw residual looks nontrivial, the uncertainty-adjusted residual is less than one error unit.
The displayed sum groups the uncertainty sources conceptually. On the main-date holdout panel, the actual price-unit denominator is
with a floor of 0.02. Here \(h_i\) is the entry half-spread, the exit term uses the expiry’s median half-spread, and \(b_i\) is the hedge buffer. The fit component comes from two perturbed-anchor refits in this run. Historical daily scoring omits that bootstrap and MC term; it retains the spread, model-disagreement, and hedge components.
Cheap and Rich Signals in Position Direction
The sign of the residual maps into a trade direction.
If:
\[
e_i=V_i^{model}-V_i^{mkt}>0
\]
then the model price is above the market mid. The option looks cheap, so the natural direction is to buy the option:
\[
q_i=+1
\]
If:
\[
e_i<0
\]
then the model price is below the market mid. The option looks rich, so the natural direction is to sell the option:
\[
q_i=-1
\]
But direction is only taken after the residual beats the no-trade filters. A cheap option with a tiny residual isn’t useful. A rich option with high model disagreement isn’t useful. The clean signal is a residual with the right sign, enough magnitude, acceptable spread, and enough historical validation support.
The watchlist uses the magnitude of the standardized residual to identify quotes for a closer look, while the strict candidate rule uses a higher threshold.
We record the residual sign and its uncertainty-adjusted magnitude separately.
No-Trade Zones
We compare each residual with its combined error estimate to define the strict candidate threshold.
A quote is inside the strict candidate band when:
\[
|e_i|\le \epsilon_i^{total}
\]
or equivalently:
\[
|z_i|\le 1
\]
When this holds, the residual is inside the combined uncertainty band. Entering a trade there means paying spread and hedge costs for a signal that may be indistinguishable from noise.
We also keep a watchlist for \(|z|\) above 0.5. A quote can therefore enter the watchlist while still remaining inside the stricter no-trade band of one total-error unit.
The watchlist is a separate, lower threshold. In the historical schedule we consider watchlist entries that also satisfy the direction, spread, maturity, and spacing rules.
The threshold is expressed in the same dimensionless residual units as the watchlist score.
models_for_trade = ["svi"]if sabr_in_signal: models_for_trade.append("sabr")if ssvi_ok_for_trade: models_for_trade.append("ssvi")if heston_ok_for_trade: models_for_trade.append("heston")err_map = model_comparison.set_index("model")["weighted_iv_rmse"].to_dict()model_wt = {name: 1.0/max(err_map.get(name, svi_err), 1e-4) **2for name in models_for_trade}wt_sum =sum(model_wt.values())for name in model_wt: model_wt[name] /=max(wt_sum, 1e-12)if"svi"in model_wt and model_wt["svi"] >0.70: extra = model_wt["svi"] -0.70 model_wt["svi"] =0.70 receivers = [m for m in model_wt if m !="svi"] recv_sum =sum(model_wt[m] for m in receivers)for m in receivers: model_wt[m] += extra * model_wt[m] /max(recv_sum, 1e-12)if"ssvi"in model_wt and model_wt["ssvi"] >0.20: extra = model_wt["ssvi"] -0.20 model_wt["ssvi"] =0.20 receivers = [m for m in model_wt if m !="ssvi"] recv_sum =sum(model_wt[m] for m in receivers)for m in receivers: model_wt[m] += extra * model_wt[m] /max(recv_sum, 1e-12)if"sabr"in model_wt and model_wt["sabr"] <0.20andlen(model_wt) >1: need =0.20- model_wt["sabr"] donors =sorted([m for m in model_wt if m !="sabr"], key=lambda m: model_wt[m], reverse=True)for m in donors: take =min(need, max(0.0, model_wt[m] -0.02)) model_wt[m] -= take model_wt["sabr"] += take need -= takeif need <=1e-12:breakwt_sum =sum(model_wt.values())for name in model_wt: model_wt[name] /=max(wt_sum, 1e-12)px_cols = [f"{name}_price"for name in models_for_trade iff"{name}_price"in fair_values.columns]iv_cols = [f"{name}_iv"for name in models_for_trade iff"{name}_iv"in fair_values.columns]fair_values["ensemble_price"] =0.0fair_values["ensemble_iv"] =0.0for name in models_for_trade:iff"{name}_price"in fair_values.columns: fair_values["ensemble_price"] += model_wt[name] * fair_values[f"{name}_price"]iff"{name}_iv"in fair_values.columns: fair_values["ensemble_iv"] += model_wt[name] * fair_values[f"{name}_iv"]fair_values["model_spread"] = fair_values[px_cols].std(axis=1).fillna(0.0)fair_values["mc_error"] = fair_values["heston_mc_standard_error"] if"heston_mc_standard_error"in fair_values.columns else0.0display(pd.DataFrame({"model": list(model_wt.keys()), "ensemble_weight": list(model_wt.values())}))
We fit 54 anchor quotes and score 60 holdout quotes across six expiries. The split is useful because it prevents us from calling a quote cheap just because the model was fitted directly to it. We fit on representative anchors, then score nearby holdout quotes.
After the uncertainty filters, the holdout quotes fall into the following groups:
41 cheap and 19 rich strict-side observations exist before filters,
36 are inside entry spread,
20 are inside model disagreement,
1 is inside fit or hedge buffer,
only 3 watchlist candidates survive.
After accounting for spread and model disagreement, only three quotes exceed the watchlist threshold on the selected date; none exceeds the stricter \(|z|>1\) threshold.
The top watchlist candidates are all modest. The largest displayed cheap call has a price residual around 0.34 with total error around 0.53, so \(z\approx 0.65\). That is a watchlist signal rather than an aggressive trade.
The stricter and watchlist thresholds describe different levels of the same uncertainty-scaled residual.
The fitted ensemble assigns 0.70 to SVI, about 0.27 to SABR, and about 0.03 to SSVI. SVI leads the fair-value estimate, with SABR providing a structural cross-check and SSVI adding a smaller global-surface contribution.
Holdout Quotes and Relative Value
If a model is fitted directly to a quote, the residual on that quote is partly in-sample noise. A flexible model can reduce residuals on fitted points even when it has no predictive value. Holdout scoring helps avoid that.
We fit on anchor quotes and score holdout quotes. The holdout residual is closer to a true fair-value estimate because the model had to infer that quote from the surrounding surface. If a holdout quote is far from the model and outside uncertainty bands, it is more interesting than an in-sample residual.
The split also creates a practical trading logic: market makers and relative-value traders often build a smooth surface from liquid anchor points, then look for off-surface quotes away from the anchors. We use a simplified version of that idea for the holdout comparison.
The holdout quotes are separate observations from the same expiry and date as the fitted anchors.
Show code
z = fair_values.copy()z["k_bin"] = pd.cut(z["k"], np.linspace(-0.26, 0.26, 13))z["dte_bin"] = pd.cut(z["dte_days"], np.linspace(21.0, 135.0, 9))h = z.pivot_table(index="dte_bin", columns="k_bin", values="z_residual", aggfunc="median", observed=False)lim =float(np.nanquantile(np.abs(h.to_numpy(dtype=float)), 0.90)) if np.isfinite(h.to_numpy(dtype=float)).any() else2.0fig, axes = plt.subplots(1, 3, figsize=(15.0, 4.5))im = axes[0].imshow(h.to_numpy(dtype=float), origin="lower", aspect="auto", cmap="coolwarm", vmin=-lim, vmax=lim)axes[0].set_title("candidate residual heatmap")axes[0].set_xlabel("k bucket")axes[0].set_ylabel("maturity bucket")axes[0].set_xticks(np.arange(h.shape[1]))axes[0].set_xticklabels([f"{b.mid:.2f}"for b in h.columns], rotation=45, ha="right")axes[0].set_yticks(np.arange(h.shape[0]))axes[0].set_yticklabels([f"{b.mid:.0f}"for b in h.index])fig.colorbar(im, ax=axes[0], pad=0.01)axes[1].scatter(fair_values["total_error"], fair_values["price_residual"], c=fair_values["candidate"].astype(int), cmap="viridis", s=26, alpha=0.78)axes[1].axhline(0.0, lw=0.9, alpha=0.65)axes[1].set_title("residual versus total uncertainty")axes[1].set_xlabel("total uncertainty")axes[1].set_ylabel("model fair - market mid")cheap = fair_values.nlargest(8, "z_residual")rich = fair_values.nsmallest(8, "z_residual")bars = pd.concat([cheap, rich]).sort_values("z_residual")labels = [f"{r.option_type[0]}{r.dte_days:.0f}d {r.k:+.2f}"for r in bars.itertuples()]axes[2].barh(labels, bars["z_residual"], alpha=0.82)axes[2].axvline(0.0, lw=0.9, alpha=0.65)axes[2].set_title("top rich and cheap options")axes[2].set_xlabel("uncertainty-adjusted residual")plt.tight_layout()plt.show()
Most holdout residuals fall within the uncertainty-adjusted no-trade region. Raw residuals exist across the surface, but once we account for entry spread and model disagreement, most points sit inside the no-trade zone. We can compare their magnitude with both entry spread and model disagreement. A surface model will always produce residuals because real quotes are noisy. The useful question is whether those residuals exceed the uncertainty required to trade them.
The Main-Date Watchlist
The small watchlist is a sign that the filters are strict. We start from a fitted surface with thousands of quotes, but after holdout scoring, spread adjustment, model disagreement, and error buffers, only a few quotes survive as watchlist candidates.
This is realistic for liquid SPX options. In a competitive market, obvious model-free mispricings shouldn’t appear everywhere. Most deviations from a fitted model are simply the market’s bid-ask, uncertainty, or a different model assumption.
The small watchlist also protects the later backtest. If the strategy selected hundreds of trades per day, the result would be dominated by correlated exposures and transaction costs. A sparse schedule gives a cleaner test of whether residuals have any useful information.
The watchlist threshold is lower than the strict candidate threshold.
16) Historical Residual Validation
A model-selected residual only matters if it predicts something out of sample. We therefore fit trade models across historical signal dates and test whether residuals are related to next-day hedged P&L.
The validation object uses a next-day hedged payoff. The idea is:
score a quote at date \(t\),
hold it for one day,
hedge delta exposure,
subtract scheduled costs,
check whether the residual direction predicted the next outcome.
The delta hedge removes first-order spot exposure, so the remaining P&L is more about volatility, curvature, residual repricing, and model error. It is never perfectly pure because delta changes, volatility moves, spreads matter, and discrete hedging creates tracking error.
The first two terms separate the option mark change and the underlying hedge. The sign \(q\) changes both exposures together, while costs are paid on entry and exit regardless of position direction. The one-day validation uses the next available matching contract mark; the scheduled experiment later follows its own fixed-hold hedging rules.
For the historical panel, we refit the eligible smile models on each signal date. The primary beta, eligible-model decisions, and ensemble weights are chosen on the selected main date and then held fixed across those historical fits.
Spearman IC and Signal Ranks
The validation uses rank correlation because residual magnitudes can be noisy. Spearman correlation asks whether larger signals tend to rank with better future outcomes. If \(R(z_i)\) is the rank of the residual score and \(R(\Pi_i)\) is the rank of future hedged payoff, then Spearman correlation is the Pearson correlation of those ranks:
\[
\rho_s=\text{Corr}(R(z_i),R(\Pi_i))
\]
This is more robust than raw correlation when the payoff distribution has outliers. Options have nonlinear payoffs, spreads, and jump moves, so rank-based diagnostics are often more stable.
The daily IC is the same idea applied date by date:
\[
IC_t=\rho_s(z_{i,t},\Pi_{i,t+1})
\]
We then average earlier daily ICs to estimate whether the residual direction has worked over the available history. If training IC is negative, the residual sign is flipped. This is why the strategy can adapt when the model residual behaves more like a contrarian signal than a direct fair-value signal.
Ranks describe ordering rather than the size of each price change. An IC of \(+1\) corresponds to the same ordering, \(-1\) to the opposite ordering, and a value near zero to little monotonic association in that cross-section. Dates with fewer than five matched observations do not receive a daily IC. Training uses an expanding mean of available earlier daily ICs, shifted by one date so that the current day’s outcome is not included.
Training IC and Direction Learning
The signal direction is learned because residuals can behave in either direction. In theory, if a model says an option is cheap, buying it should make money as price converges upward. In practice, a residual can also indicate a risk premium or model bias. Then the profitable direction may be the opposite.
We use the mean of earlier daily IC values to estimate this:
If recent IC is negative, the strategy treats raw residuals as contrarian in the opposite direction. This is a very practical idea. It lets the data decide whether model residuals are converging or whether they are systematically pointing the wrong way under the current market microstructure.
In the implemented expanding estimate, \(M\) is the number of available earlier daily IC observations, rather than a fixed rolling-window length. At least two past ICs are required. The direction is \(+1\) above 0.02, \(-1\) below \(-0.02\), and zero in between; the sign expression above describes the direction outside that threshold band.
For the validation score we first set \(q=\operatorname{sign}(e)\) and calculate \(\Pi^{raw}=q(V_{t+1}-V_t)-q\Delta_t(S_{t+1}-S_t)-\text{costs}\). The displayed direction-adjusted quantity is \(d_t\Pi^{raw}\). Thus a negative \(d_t\) reverses the entire raw score, including its cost term; the scheduled hedging backtest calculates its own trade costs from the selected positions.
Validation Sample Size and Dependence
The validation panel has many rows, but option rows are not independent in the same way stock returns might be. Quotes from the same date and same expiry share the same underlying move, volatility move, and liquidity environment. So 28,415 validation rows don’t mean 28,415 independent experiments.
The daily IC aggregates the quote-level association into one cross-sectional measure per date. It aggregates cross-sectional predictive behavior by date. The model can have many quote-level rows but still only about 504 daily environments. If the IC flips sign over time, the strategy has to adapt.
The decile table also has dependence. Options in the same decile can cluster by maturity or option type. So the validation results are evidence, but not final proof. The scheduled trade experiment is a stricter test because it samples a limited number of actual entries with spacing.
We score 504 signal dates and obtain 28,415 matched validation rows. The Spearman correlations are negative for all residual scores:
raw residual: about -0.155,
spread-adjusted: about -0.144,
uncertainty-adjusted: about -0.147.
The negative sign means the raw residual direction needs to be interpreted carefully. The training IC column later flips signal direction based on recent performance. The daily IC table shows large variation: some dates have strongly positive IC, others strongly negative. By late 2023 the training IC is around -0.15, so the system learns that the residual should be traded in the opposite direction.
The decile results add detail to the overall correlation. Mean next hedged P&L is positive in every decile, but the middle-positive residual deciles have the highest values, with hit rates above 70% in several bins. The highest \(z\) decile has lower mean P&L than mid-high deciles, which warns us that extreme residuals may include more model error or liquidity noise rather than cleaner opportunity.
Daily IC is calculated across the matched quotes available on each date.
Interpreting the Non-Monotonic Deciles
The decile table doesn’t show a perfect monotonic edge. A steadily ordered set of decile outcomes would be easier to interpret as a consistent ranking signal. Option residual trading is messier.
A mid-level residual can be tradable because it is large enough to matter but still inside a liquid, stable region. An extreme residual can be dangerous because it may come from:
stale quote,
low vega,
jump risk,
wide spread,
model extrapolation,
or a true market premium that doesn’t mean-revert.
So the fact that the middle-high deciles show strong hit rates while the highest decile weakens is economically believable. It tells us to use residuals with cost and model-risk gates rather than just buying the largest absolute residual.
The decile averages are calculated from the matched next-day observations.
The validation results have three features worth examining together:
the distribution of daily ICs is wide,
training IC evolves through time,
decile performance isn’t monotonic in a clean textbook way.
That is a realistic result. Option relative value is noisy. A residual can mean mispricing, but it can also mean quote staleness, model misspecification, or temporary supply-demand pressure that doesn’t revert quickly. The validation layer stops the project from assuming every model residual is alpha.
17) Scheduled Delta-Hedged Trading Experiment
After validation, we run a scheduled trading experiment. The goal is to compare two entry selectors:
surface residual selector, which chooses trades based on the model residual logic,
matched ATM baseline, which chooses similar option type, date, sign, and maturity but targets a more ATM-like contract.
The residual strategy schedules up to 80 entries, spaced at least five calendar days apart. This prevents the backtest from taking dozens of overlapping signals from the same date and pretending they are independent.
Each trade is held for three days. The backtest evaluates delta-hedged P&L with bid-ask and hedge costs. The delta hedge is meant to reduce directional exposure:
\[
\text{hedge shares}\approx -q\Delta
\]
If we are long one call with \(\Delta=0.45\), we short about \(0.45\) units of the underlying per option multiplier. If the index rises, the call gains but the short hedge loses. The remaining P&L comes from volatility movement, gamma, theta, residual convergence, and execution costs.
For this SPX application, the hedge uses SPY rather than direct index units. A rolling 60-day SPX-to-SPY return beta and the relative index/ETF levels translate the index delta into the hedge instrument; the backtest also applies a delta band, share lots, and transaction costs. Entry spacing is at least five calendar days in the scheduling function, with entries drawn from available signal dates.
Delta-Hedged P&L Through a Taylor Expansion
The delta-hedged P&L can be understood through the option’s local expansion:
\[
dV \approx \Delta dS+\frac{1}{2}\Gamma(dS)^2+\mathcal{V}d\sigma+\Theta dt+\text{higher-order terms}
\]
If we hedge with \(-\Delta\) units of the underlying, the first-order spot term is reduced:
After removing first-order direction, the hedged position still has:
gamma exposure to realized movement,
vega exposure to implied-volatility changes,
theta decay,
residual convergence or divergence,
re-hedging costs,
bid-ask costs,
discrete hedging error.
A relative-value option signal is usually trying to earn from residual convergence after controlling delta. It is not guaranteed to be positive because gamma/theta/vega can dominate over a short holding period.
The expansion is a local approximation for one interval. Gamma and vega are evaluated at the starting state, and larger moves also bring cross-effects and changes in the Greeks. Rebalancing changes the hedge for the next interval, while the realized backtest P&L is computed from option and hedge marks rather than by summing the Taylor terms.
Gamma, Theta, and Vega Interactions
A three-day option trade can lose money even when the fair-value residual is correct because the Greeks interact.
If we buy a cheap option, we are usually long gamma and long vega:
\[
\Gamma>0,\quad \mathcal{V}>0
\]
Long gamma benefits from realized movement:
\[
\frac{1}{2}\Gamma(dS)^2
\]
Long vega benefits from implied volatility rising:
\[
\mathcal{V}d\sigma
\]
But long options usually pay theta:
\[
\Theta dt<0
\]
and the bid-ask spread is paid at entry and exit. So the residual has to converge fast enough to beat theta and costs. Short options reverse these exposures: they collect theta but suffer from gamma and vega shocks.
This is why residual trading is difficult. A model can identify a cheap option, but the position still has Greek exposures over the holding window. The hedged backtest is needed because static residual ranking alone can’t answer that.
P&L Components
The attribution is split into:
\[
\Pi_{\text{net}}=\Pi_{\text{option mid}}+\Pi_{\text{delta hedge}}-\text{bid-ask and hedge costs}
\]
where:
\(\Pi_{\text{option mid}}\) is the change in option mid prices,
\(\Pi_{\text{delta hedge}}\) is the P&L from rebalancing the underlying hedge,
costs include option entry/exit spread and hedge trading costs.
Breaking down P&L helps us see the different sources of a loss:
the option signal can be wrong,
the delta hedge can lose money,
costs can overwhelm a small edge,
the selected options can carry too much vega or gamma risk.
A model-based strategy that produces lower gross loss and lower exposure than the matched baseline can still lose net money if spreads and hedging costs are too large.
The individual components use the same option and hedge positions as the reported net P&L.
Matched ATM Baseline
A model-selected option can have a different delta, maturity, spread, and vega than a generic option. If we compare it only against zero, we can’t tell whether the result came from the residual signal or from general option exposure.
We choose a matched option with similar broad characteristics, without using the residual-selection rule. The comparison asks:
The comparison can therefore reveal differences in both P&L and exposure. We will examine the net outcomes alongside premium, notional, vega, delta, and cost measures for the fixed three-day schedule.
The baseline retains the selected trade date and position sign while choosing an ATM-like option of the same type.
The matched-entry ranking also includes the number of available future quote marks, which helps select markable contracts in this historical comparison.
We schedule 80 residual entries and 80 matched baseline entries, with a same-type match rate of 100%. That gives a fair comparison: the matched strategy isn’t a random portfolio; it is designed to share broad trade characteristics with the residual schedule.
The hedge-book summary shows 160 scheduled entries and 6,658 hedge-book rows across 159 contracts. The residual strategy and matched strategy both have enough marks to run.
After costs, the performance comparison is:
matched ATM delta strategy loses about $26,621,
surface residual delta strategy loses about $9,032,
both are negative after costs,
the residual strategy has lower drawdown and much lower traded notional.
So the model-selected residual strategy improves materially over the matched ATM baseline, but it still doesn’t produce a positive net strategy in this configuration. That is a much stronger conclusion than simply saying “the model worked.” The signal has some selection value, but trading frictions and short-horizon option noise are still too strong.
The comparison uses matched entry dates, option types, signs, and broadly similar maturities.
Exposure Differences Between Residual and Matched Trades
The matched baseline trades much larger notional than the residual strategy: about $37.37 million traded notional for matched ATM versus about $4.80 million for surface residual. That difference is material. It means the residual selector naturally chooses a different exposure profile.
Mean absolute delta after hedging is also much lower for the residual strategy, around 8.0 versus 28.7. Mean absolute vega is lower too, around 17,677 versus 45,980.
The residual selector changes the exposure profile of the chosen options as well as the entry direction. The strategy loses less partly because it takes smaller and more controlled exposures. That is a valid risk-management benefit, even though the net P&L remains negative.
The notional and Greek measures describe different aspects of the positions held through the backtest.
Reading Negative P&L Honestly
Both scheduled strategies lose money after costs, even though the residual-selected trades have smaller losses and a different exposure profile.
A model can be useful and still fail as a naive trading strategy. The SPX residual selector:
loses much less than the matched baseline,
has lower notional,
has lower absolute delta after hedging,
has lower vega,
has lower total cost,
still loses net money.
That combination says the model is doing some useful selection, but the execution rule is not yet strong enough. Possible reasons include:
three-day horizon is too short or too long,
thresholds are too loose,
spreads are too expensive,
signals are too small,
delta hedging is too crude,
residual direction is unstable.
These observations give us specific directions to investigate while keeping the reported three-day experiment and its results intact.
We can use the separate exposure and cost measures to understand how the two selections differ.
The normalized performance table shows the residual strategy has worse P&L per option premium than the matched baseline, but it uses much less notional and has smaller absolute losses. The cost percentage of premium is higher for the residual strategy, around 3.4% versus 1.8% for the matched baseline, which means it trades options where spreads are more painful relative to premium.
The P&L attribution gives us another view of the difference:
matched ATM option mid P&L is about -$17,135,
residual option mid P&L is about -$13,475,
residual delta hedge P&L is much better, around +$8,446 versus +$3,594,
residual costs are much lower, around $3,910 versus $14,372.
The option, hedge, and cost components explain much of the difference in losses. It isn’t because every option residual converges. It is because the selected trades have better hedge behavior and lower total cost exposure than the matched baseline.
The fixed-hold equity curves follow the same pattern as the performance table. Both strategies draw down, but the matched ATM baseline loses more and has a much larger drawdown. The residual selector is less damaging and more controlled. That still matters for research: a relative-value model can first prove that it selects better bad trades before it becomes good enough to generate positive trades.
For a production-style version, the next improvements would be tighter entry thresholds, fewer trades, more liquidity constraints, and perhaps longer or adaptive holding windows. We don’t add those here because the point of this project is model comparison and signal validation, not overfitting a trading rule until it looks good.
BTC is a much harder market than SPX for this kind of workflow. The surface has wider wings, fewer institutional-style quote regularities, 24/7 trading, different settlement conventions, and a much more jump-sensitive underlying. The same model stack therefore faces a different economic problem:
SPX smiles are dominated by downside equity crash protection and structured product flows.
BTC smiles often reflect jump risk in both directions, leverage liquidations, weekend gap risk, and unstable liquidity.
Delta hedging is harder because the underlying trades continuously but option marks can be missing or irregular.
These differences give us a useful second market for applying the same library workflow.
Model Hierarchy in BTC Options
BTC options force the same models into a different volatility ecosystem.
SPX has a structurally persistent downside skew because investors buy crash protection and dealers manage equity-index risk. BTC can have much more two-sided jump risk: upside squeezes, downside liquidations, weekend gaps, exchange-specific liquidity effects, and sharp volatility repricings.
That changes which model families look useful:
Merton-style jumps can become more competitive because discontinuous moves are central to BTC.
Heston can struggle if the surface wants abrupt jumps rather than continuous variance shocks.
SSVI can struggle because one global smooth surface discipline may be too restrictive.
SVI can still fit well, but its static fit doesn’t guarantee tradability.
SABR can help as a structural fast check, especially when \(\beta=1\) is used for lognormal-style crypto moves.
The secondary implementation is therefore a real market stress test: it checks whether the same architecture adapts when the market’s volatility economics change.
Each model still uses the corresponding BTC quote panel and the same error metrics for comparison.
BTC Surface Economics
BTC options tend to price several risks that are less central in SPX:
weekend movement while traditional markets are closed,
liquidation cascades from leverage,
exchange-specific liquidity changes,
jumps after ETF, regulatory, or macro headlines,
volatility demand from both upside and downside speculation,
unstable term structure when spot trends strongly.
This can make BTC smiles less smooth and less one-sided than SPX. In SPX, the left wing is often structurally rich. In BTC, both tails can matter depending on the regime. A huge upside squeeze can be just as important as a crash in some periods.
A jump distribution gives us a way to represent discontinuities and gap risk when we fit the BTC option panel.
We will compare how well the jump and continuous-variance models fit the selected BTC quotes.
With the library, we repeat the main components of the workflow:
clean BTC option quotes,
attach BTC spot and U.S. rates,
infer forwards from parity,
compute IVs and Greeks,
fit SVI, SSVI, SABR, Merton, Heston, and Bates,
compare model errors and runtime,
validate residual scores,
schedule residual and matched baseline trades,
run unhedged and delta-hedged backtests under 24/7 crypto calendar assumptions.
We infer forwards from available call-put pairs and track the share of observations that need a fallback. Those forwards determine the log-moneyness used by the surface models.
The subsequent model fits and trading experiments use the resulting BTC panel and its observed quote support.
On the selected BTC benchmark panel, the model comparison differs from SPX:
Merton has the lowest weighted IV RMSE, about 0.0048,
SVI is around 0.0156,
SABR is around 0.0182,
Bates is around 0.0213,
Heston is around 0.0369,
SSVI is weakest at about 0.0492.
The jump-diffusion fit is especially strong on this benchmark. For BTC, the pure jump-diffusion diagnostic is very strong on the common benchmark. That makes sense because BTC option smiles often price abrupt jumps and gap-like moves very heavily. A jump component can explain a lot of the surface shape that SPX explained through skewed equity crash demand and smoother stochastic variance.
SSVI struggles because BTC smiles are too irregular for a globally disciplined surface to fit tightly. Heston also struggles because continuous stochastic variance alone isn’t enough to match the jump-like shape in this panel.
On this selected date, the reported forward fallback share is approximately 0.053%, so parity inference supplies almost all of the forward inputs.
BTC Residual Validation and Hedging Friction
Mean next-period hedged P&L is positive in most BTC residual deciles, although the pattern is irregular. The highest \(z\) decile is negative, which fits the earlier warning that extreme crypto residuals can be contaminated by microstructure and jump risk.
The hedging comparison also shows that delta hedging in BTC can reduce drawdown but also damage net P&L through trading costs and hedge losses. In the matched ATM case, unhedged P&L is positive while delta-hedged P&L is negative. That means the underlying directional exposure helped the unhedged result, and hedging removed that benefit while adding costs.
For residual trades, delta hedging also worsens total P&L, but it reduces daily volatility and drawdown. This creates a real strategy design choice. If the goal is pure option relative value, delta hedging is necessary. If the goal is total return, the unhedged exposure may dominate, but then the result is less about model residuals.
The table compares the same scheduled entries with and without the hedge.
The BTC residual deciles have an uneven relationship with next-period hedged P&L. Mean next hedged P&L is positive in most deciles, but the highest residual decile has negative mean P&L. That again warns us against blindly trading the most extreme residuals. In BTC, extreme residuals may be exactly where liquidity is weakest, marks are less reliable, or jump risk is most mismeasured.
The schedule has 41 residual entries and 41 matched entries. The BTC hedging comparison separates unhedged and delta-hedged versions:
matched ATM unhedged gains about $1,705,
matched ATM delta-hedged loses about $6,395,
residual selector unhedged loses about $560,
residual selector delta-hedged loses about $1,277.
Delta hedging reduces volatility and drawdown for BTC, but it also turns the matched unhedged profit into a loss because hedge costs and hedge P&L dominate. For residual trades, both unhedged and delta versions lose, but the delta version has smaller daily volatility and drawdown.
The BTC backtest diagnostics record the following details:
the residual run has missing option marks on some days,
the matched run has no missing marks in the displayed diagnostics,
the calendar is 24/7,
the contract size is 1,
valuation is USD equivalent,
both unhedged and delta strategies are run.
The BTC calendar, contract conventions, and mark availability affect the backtest mechanics. The trading calendar, liquidity, surface shape, and jump behavior change the meaning of the same model residual.
The final multi-panel figure gives a compact visual summary: BTC quote map, smile term structure, SVI/SSVI/SABR fits, Merton tail fit, Heston convergence, Heston/Bates comparison, model runtime/error, residual deciles, and scheduled hedge comparison. The strongest interpretation is that jump-style models are more competitive in BTC than in SPX, but model residuals still need strict cost and validation gates before they become tradable.
BTC Model Ranking and Surface Shape
The BTC repeat is useful because it puts the same volatility-modeling stack into a market with a very different risk structure. SPX options are dominated by institutional index hedging, crash insurance demand, and a long history of left-skew behavior. BTC options also price downside risk, but they add another layer: large jump moves, weekend gaps, liquidity differences across expiries, and a market where the underlying itself can move like a high-volatility growth asset and a macro liquidity proxy at the same time.
The model comparison reflects that difference. In the SPX run, SVI is the clean surface-marking core. It fits the dense SPX smile with very low IV error and gives a stable fair-value engine. In the BTC run, the Merton jump-diffusion becomes much more competitive on the benchmark error table. That makes economic sense: BTC return distributions often have abrupt large moves, and a jump component can explain tail option prices with fewer contortions than a smooth diffusion-only smile.
For this sample, BTC rewards jump-aware explanations more strongly than SPX does in this sample. SVI and SABR still matter because they give fast and flexible surface marks, but the jump model carries more information about the option cross-section than it did for SPX.
The cleaned BTC panel also shows which quotes enter the secondary application. The raw chain loses many rows through price, spread, DTE, and moneyness filters before we reach the final panel. That filtering is good research hygiene, but it also means the final comparison is based on the liquid and usable part of the BTC surface. We are modeling the tradable surface, not every quote shown by the venue.
The common benchmark gives the six models the same selected quotes for the BTC comparison.
BTC Residual Validation
We use the BTC residual deciles to check the same relationship as in SPX: after a model marks a quote as cheap or rich, does that residual have any relationship with the next-period hedged payoff?
For each option observation, the residual is built as:
where \(V_i^{\text{model}}\) is the model fair value, \(V_i^{\text{mid}}\) is the market mid price, and the denominator combines quote noise, model disagreement, and fit uncertainty. A positive \(z_i\) means the model value is above the market mid after scaling by uncertainty, so the option looks cheap. A negative \(z_i\) means the model value is below the market mid, so the option looks rich.
The decile table doesn’t show a clean monotonic pattern across all buckets. Some middle deciles have strong average next-period payoffs, while the most extreme positive bucket is weaker. The uneven ordering is useful to keep in mind when interpreting a large residual. A residual can be mathematically large because the model is confident, because the quote is wide, because the tail is difficult, or because the option is genuinely mispriced. The validation step separates these possibilities by asking what happened afterward.
For BTC, the evidence is mixed rather than clean. There is useful information in the residuals, but the relationship is less stable than the surface error table alone would suggest. That’s why the secondary run uses scheduled trade tests instead of stopping at model fit. A model that fits the smile well still needs the residual to survive spreads, hedge noise, and the large realized volatility of the underlying.
The denominator has price units, making \(z_i\) dimensionless. A value of \(z_i=1\) means that the positive price residual equals one unit of the chosen uncertainty scale; it doesn’t by itself imply a statistical one-standard-deviation probability because the denominator combines several different buffers.
BTC Scheduled Trades and Delta Hedging
The BTC scheduled-trade section compares two ways of entering option positions. The matched ATM strategy gives us a simple reference trade around the same date and maturity region, while the residual-selector strategy uses the model-risk-aware residual filter. Both are then evaluated with and without delta hedging.
A delta-hedged option position tries to remove the first-order exposure to BTC price moves. If the option position has delta \(\Delta_t\), the hedge holds approximately \(-\Delta_t\) units of BTC exposure:
For a perfectly delta-neutral world with stable implied volatility and continuous rebalancing, the first-order BTC move would mostly disappear. In reality, BTC jumps, the hedge is discrete, spreads are meaningful, and the option surface itself moves. So the delta hedge reduces directional exposure, but it also adds hedge turnover and extra transaction costs.
The BTC results show exactly that tradeoff. The unhedged matched ATM run has positive total P&L in this sample, while the delta-hedged version turns negative after costs. That means the unhedged benchmark likely benefited from directional BTC moves or favorable option exposure, and the hedge removed part of that benefit while adding cost. The residual-selector trades are smaller and more selective, but they also struggle after costs.
The model fit and the hedged trading result answer different questions. A model residual can look interesting in the surface, and a hedge can look theoretically cleaner, but the actual trading layer can still lose money once execution and discrete rebalancing are included.
The unhedged and hedged runs use the same scheduled option entries.
BTC Diagnostic Counters
We can use the diagnostic counters to examine missing marks and early exits in the backtest. BTC options can disappear from the markable panel, shift liquidity across strikes, or reach the DTE exit before the fixed holding period completes. These events change how a clean theoretical trade becomes an actual backtest trade.
The residual-selector run has 30 missing option marks, while the matched ATM run has no missing marks in the displayed diagnostics. That difference is meaningful. A missing mark is a serious accounting and execution issue. The backtest couldn’t observe a reliable exit value for the contract on a scheduled day, so the realized trade path becomes less clean. In option strategies, missing marks often happen in the places where the surface is thinner, spreads are wider, or strikes move away from the active trading region.
The DTE exits also deserve attention. An option position close to expiry becomes extremely sensitive to small underlying moves because gamma increases as time shrinks. If the trade reaches the DTE stop, the backtest exits before the option becomes too short-dated and noisy. This is a realistic control because most systematic volatility strategies avoid holding tiny-dated options unless they are explicitly designed for short-gamma or event-volatility trading.
The final BTC hedge comparison therefore has to be read through three layers:
model layer: which model explains the surface best,
signal layer: whether residuals predict future option payoffs,
execution layer: whether marks, spreads, and hedge costs allow the signal to survive.
BTC performs well as a secondary stress test because it pressures all three layers at once. It has enough option data to run the workflow, but the market is still noisy enough that weak signals are punished quickly.
Secondary BTC Reading
We have applied the full modeling and trading workflow to a second market with different surface and execution characteristics. The same cleaning logic, model comparison, residual scaling, validation, and scheduled-trade engine can run on BTC options, but the model ranking changes and the trading results become more sensitive to execution.
The most important BTC result is that jump-aware modeling becomes more valuable, while hedged residual trading remains difficult. Merton’s strong benchmark fit says the BTC smile contains jump-like information. The residual validation says that model mispricing signals exist but aren’t uniformly clean across deciles. The hedge comparison says that removing directional BTC exposure can reduce drawdown in some places, but the extra costs can dominate the benefit.
The secondary results connect the surface fit with the signal and trading outcomes. BTC options are rich enough for advanced volatility modeling, but the trading edge is harder than the model fit. The surface can be explained; the residual can be ranked; the trade can still fail after costs. That final separation between model fit, signal value, and tradable P&L is the main lesson from the BTC implementation.
