Overview:

Given the current SPY option chain, the pipeline performs the following steps:

1. Downloads the current SPY price and all listed call and put option chains.
2. Uses the midpoint of the bid and ask as the option's fair price and removes quotes with invalid or empty markets.
3. Checks that each option price lies inside the Black-Scholes no arbitrage bounds before attempting to solve for volatility.
4. Keeps out of the money puts below the forward and out of the money calls at or above the forward, avoiding duplicate information from the in the money side.
5. Solves for the implied volatility of every retained option using Newton-Raphson, with Brent's method as a fallback when Newton-Raphson becomes unstable.
6. Converts implied volatility into total variance and expresses strike through log-moneyness.
7. Fits a separate raw SVI smile to each expiration as a baseline model.
8. Fits one SSVI parameterization across the full set of expirations, with penalties for sufficient butterfly arbitrage conditions.
9. Interpolates the ATM total variance curve and plots the fitted surface together with the original market implied volatility points.
10. Reports the fitted SSVI parameters, optimizer status, and implied volatility RMSE.

Methodology:

Market data

Option chains are downloaded through yfinance. The current implementation uses SPY, although the data function accepts a different ticker. Time to expiry is measured in whole calendar days:

T = days_to_expiry / 365

Same day expirations are excluded because they produce T = 0. For each remaining contract, the midpoint of the displayed market is used as the fair price:

market_price = (bid + ask) / 2

Quotes are removed when the bid is nonpositive or the ask is not above the bid.

Black-Scholes pricing

The final SSVI implementation uses the dividend adjusted Black-Scholes model:

d1 = [ ln(S0 / K) + (r - q + sigma^2 / 2)T ] / (sigma * sqrt(T))
d2 = d1 - sigma * sqrt(T)

call = S0 * exp(-qT) * N(d1) - K * exp(-rT) * N(d2)
put  = K * exp(-rT) * N(-d2) - S0 * exp(-qT) * N(-d1)

where S0 is spot price, K is strike, T is time to expiry, r is the risk free rate, q is dividend yield, and sigma is volatility. The current implementation uses r = 0.04 and q = 0.012.

Price bounds

A market price must lie inside the Black-Scholes no arbitrage bounds before an implied volatility can exist. The bounds used are:

call: max(S0 * exp(-qT) - K * exp(-rT), 0) < C < S0 * exp(-qT)
put:  max(K * exp(-rT) - S0 * exp(-qT), 0) < P < K * exp(-rT)

Options outside these intervals are removed before root finding. This prevents the solver from searching for a volatility that cannot reproduce the observed price.

Implied volatility

Implied volatility is found by solving:

BlackScholes(sigma) - market_price = 0

Newton-Raphson updates volatility according to:

sigma_(n+1) = sigma_n - [ BlackScholes(sigma_n) - market_price ] / Vega(sigma_n)

with Black-Scholes vega:

Vega = S0 * exp(-qT) * normal_pdf(d1) * sqrt(T)

The final implementation starts every option at sigma = 0.30. A fixed seed avoids the near zero initial values that can result from an ATM approximation when it is applied to far from the money options. If vega becomes too small, Newton-Raphson produces a nonpositive volatility, or the iteration does not converge within 50 steps, the code falls back to Brent's method on the interval:

0.0001 <= sigma <= 5.0

The baseline SVI script uses a Brenner-Subrahmanyam starting approximation, while the SSVI version replaces it with the more reliable fixed seed and one centralized Brent fallback.

Out of the money selection

The SSVI model uses the forward price:

F(T) = S0 * exp((r - q)T)

and retains:

puts  when K < F(T)
calls when K >= F(T)

By put-call parity, the in the money call and out of the money put at the same strike contain equivalent volatility information, and likewise on the opposite side of the forward. Keeping only the out of the money contract avoids double counting and generally uses the more actively quoted side of the market.

Log-moneyness and total variance

The SSVI model is expressed using forward log-moneyness:

k = ln(K / F(T))

so k = 0 is at the money forward, negative k is below the forward, and positive k is above the forward.

Implied volatility is converted into total implied variance:

w(k, T) = implied_volatility^2 * T

Total variance is the quantity modeled directly by both SVI and SSVI.

Baseline SVI model

The first version fits the raw SVI parameterization independently at each expiration:

w(x) = a + b [ rho(x - m) + sqrt((x - m)^2 + sigma_bar^2) ]

where the baseline uses spot log-moneyness:

x = ln(K / S0)

Each expiration receives its own values of a, b, rho, m, and sigma_bar. The optimizer minimizes squared error in total variance, with parameter bounds and the positive minimum variance constraint:

a + b * sigma_bar * sqrt(1 - rho^2) >= 0

At least five observations are required for an expiration because the raw SVI slice has five parameters. The fit is blanked outside the observed strike range rather than extrapolated.

This produces useful individual volatility smiles, but every expiration is calibrated separately. There is no shared structure across maturity, so the slices can appear disconnected and the model does not directly control calendar consistency.

SSVI model

The final model uses the SSVI total variance parameterization:

w(k, theta) = theta / 2 * [
    1 + rho * phi(theta) * k
      + sqrt((phi(theta) * k + rho)^2 + 1 - rho^2)
]

with the power law shape function:

phi(theta) = eta * theta^(-gamma)

The model fits three global parameters:

• rho controls the direction and strength of skew.
• eta controls the overall size of the smile slope and curvature.
• gamma controls how the smile shape changes as ATM total variance changes.

theta is the ATM total variance for a given expiration:

theta(T) = ATM_implied_volatility(T)^2 * T

An exact forward ATM option is not always listed, so ATM implied volatility is linearly interpolated at k = 0. An expiration is retained only when its observed strikes bracket zero; otherwise the ATM estimate would require extrapolation.

The optimizer chooses rho, eta, and gamma to minimize squared error between market total variance and SSVI total variance across all retained options:

SSE = sum_i [ w_SSVI(k_i, theta_i) - w_market_i ]^2

Unlike the raw SVI model, the same three parameters govern every maturity. Time enters through the ATM variance curve theta(T), producing one linked surface instead of a collection of unrelated slices.

Arbitrage controls

The implementation includes the main SSVI conditions used in the model.

ATM total variance should be nondecreasing with maturity:

d theta(T) / dT >= 0

The code checks this after estimating theta(T) and prints a warning if the observed curve is not monotonic. It currently reports a violation rather than modifying the data or forcing monotonicity.

For the power law function with:

eta > 0
0 < gamma < 1
-1 < rho < 1

the second SSVI calendar condition is satisfied by construction. These restrictions are imposed through optimizer bounds.

Two sufficient butterfly arbitrage conditions are also evaluated during calibration:

theta * phi(theta) * (1 + |rho|) < 4

theta * phi(theta)^2 * (1 + |rho|) <= 4

A squared penalty is added whenever either left hand side exceeds 4:

loss = SSE + 10000 * sum(
    max(0, condition_3 - 4)^2
  + max(0, condition_4 - 4)^2
)

The penalty makes parameter sets with excessively steep or sharply curved wings unattractive to the optimizer. This makes the calibration arbitrage aware, although the conditions are implemented as a soft penalty rather than hard constraints.

Surface construction

The fitted ATM total variance values are grouped by expiration and linearly interpolated across maturity. A grid is then created over log-moneyness and time:

k_grid: 600 points
T_grid: 400 points

At each grid location, SSVI total variance is converted back to implied volatility:

implied_volatility(k, T) = sqrt(w(k, theta(T)) / T)

The Plotly figure displays the fitted surface and overlays the market-derived implied volatility observations at their original log-moneyness and expiration.

Fit evaluation

Fit quality is reported in implied volatility rather than total variance. For every market observation, the code evaluates the SSVI surface at the same k and T and computes:

RMSE = sqrt(mean((SSVI_implied_volatility - market_implied_volatility)^2))

A representative run returned:

SSVI Calibration Results:
rho:   -0.643761
eta:    0.758794
gamma:  0.579299
Optimizer success: True
vol RMSE: 0.0147

An RMSE of 0.0147 is approximately 1.47 volatility points. The exact parameters, number of observations, number of expirations, and RMSE depend on the live option chain retrieved on the day of the run.

What's left to do:

The pipeline runs end to end on live SPY option data. The remaining work is focused on market data quality and stricter surface controls:

• Quote quality. The current model treats the bid-ask midpoint as the fair price and gives each retained option equal weight. A future version could filter or weight quotes using spread, volume, open interest, or vega.
• Rates and dividends. The risk free rate and dividend yield are fixed constants. A more complete implementation would use maturity dependent rate and dividend inputs.