Python + C++● Live

Implied Vol Surface Calibrator

SVI parametrisation calibrated to market quotes via Levenberg-Marquardt. Newton-Raphson implied vol inversion with Brent fallback. Butterfly and calendar spread arbitrage detection. Full pytest suite.

PythonSVICalibrationArbitrage-freePlotlyscipy

1. Market Context and Assumptions

Under Black-Scholes (1973), implied volatility is a constant σ. Real equity option markets quote different implied vols for different strikes and maturities — the volatility smile. This arises because real returns exhibit fat tails, negative skew (crash risk), and stochastic volatility — all features absent from log-normal dynamics.

An options desk cannot use a single constant vol to price a book of vanilla options. It needs a calibrated implied vol surface — a smooth, arbitrage-free function σ_IV(K, T) that reproduces market prices across all observable strikes and maturities. That surface is the input to exotic option pricing, delta hedging, and risk management.

Model Assumptions

  • Log-moneyness: k = log(K/F), where F = S · exp((r − q) · T) is the forward price. This normalisation makes the smile comparable across spot levels and maturities.
  • Total variance: w(k) = σ_IV²(k) · T. Working in total variance space linearises the no-arbitrage conditions and improves numerical conditioning of the calibration.
  • Rates: continuously compounded, annualised, decimal (e.g. 0.05 = 5%).
  • Day count: ACT/365 for T throughout.
  • Scope: European exercise. No jump or stochastic vol model is implied — SVI is a parametrisation of the surface, not a model.

2. SVI Parametrisation

Gatheral (2004) proposed the Stochastic Volatility Inspired (SVI) parametrisation as a five-parameter family that is consistent with the large-strike behaviour of implied vols in the Heston model (Roger Lee 2004) and easy to calibrate to market data.

The SVI raw parametrisation for a single maturity slice is:

w(k;a,b,ρ,m,σ)=a+b(ρ(km)+(km)2+σ2)w(k;\,a,b,\rho,m,\sigma) = a + b\Bigl(\rho(k-m) + \sqrt{(k-m)^2+\sigma^2}\Bigr)

where w(k)=σIV2(k)Tw(k) = \sigma_{\text{IV}}^2(k)\cdot T is the total implied variance, and:

  • aRa \in \mathbb{R} — overall variance level (controls the minimum of the smile)
  • b0b \geq 0 — slope and curvature (controls both wings simultaneously)
  • ρ(1,1)\rho \in (-1,1) — skew (negative for equity markets: puts more expensive than calls)
  • mRm \in \mathbb{R} — horizontal shift of the ATM point
  • σ>0\sigma > 0 — ATM curvature / smoothness (prevents a kink at the minimum)

The implied volatility for a given maturity TT is then recovered as:

σIV(k,T)=w(k)T\sigma_{\text{IV}}(k, T) = \sqrt{\frac{w(k)}{T}}

Parameter constraints for no-arbitrage (necessary, not sufficient):

  • b0b \geq 0
  • ρ<1|\rho| < 1
  • a+bσ1ρ20a + b\sigma\sqrt{1-\rho^2} \geq 0 (ensures w(k)0w(k) \geq 0 for all kk)

The hyperbolic form ρ(k−m) + √((k−m)²+σ²) is key: it is linear in k for large |k| (giving the correct large-strike slope from Roger Lee's moment formula), and smooth at the ATM point (no kink, controlled by σ).

3. Interactive Demo

Adjust the SVI parameters and maturity with the sliders. The Vol Smile tab shows σ_IV(k) for the chosen maturity — orange dots are simulated market quotes (SVI + small fixed noise), the purple curve is the current fit. Red bands mark butterfly arbitrage violations.

The Surface Heatmap tab shows σ_IV(k, T) across the full surface using the SSVI-style scaling w(k,T) = a·T + b·√T · (ρ(k−m) + √((k−m)²+σ²)). Red cells indicate calendar spread violations (∂w/∂T < 0).

SVI Surface Calibratorinteractive
Python / NumPy

SVI Parameters

0.040
0.150
-0.40
0.00
0.20

Maturity

1.00y
ATM vol σ_IV(0,T)
26.46%
-0.4-0.20.00.20.40%11%22%33%43%k = log(K/F)σ_IVmarket quotes

SVI fit   market quotes (SVI + noise)   butterfly violation zone (g < 0)

4. Python Implementation

4.1 SVI Kernel

The kernel module implements total variance, implied vol, and the first and second derivatives used by the arbitrage checks. All functions are vectorised over NumPy arrays of log-moneyness.

implied_vol_surface.py — SVI kernelPython 3.10+
import numpy as np
from typing import NamedTuple


class SVIParams(NamedTuple):
    """
    SVI raw parametrisation (Gatheral 2004).

    Assumptions:
        - Log-moneyness k = log(K/F), where F = S·exp((r−q)·T) is the forward.
        - Total implied variance w(k) = σ_IV²(k) · T.
        - b ≥ 0, |ρ| < 1, σ > 0, and a + b·σ·√(1−ρ²) ≥ 0.

    Reference:
        Gatheral, J. (2004). A parsimonious arbitrage-free implied volatility
        parametrisation with application to the valuation of volatility
        derivatives. Merrill Lynch, Quantitative Finance 2006.
    """
    a:   float   # level parameter  (overall variance intercept)
    b:   float   # slope/curvature  (b ≥ 0)
    rho: float   # skew correlation (|ρ| < 1)
    m:   float   # ATM location shift
    sig: float   # ATM curvature    (σ > 0)


def svi_total_variance(k: np.ndarray, p: SVIParams) -> np.ndarray:
    """
    w(k; a,b,ρ,m,σ) = a + b·(ρ·(k−m) + √((k−m)²+σ²))

    Returns total implied variance (σ_IV² · T), not implied vol itself.
    Requires division by T to obtain σ_IV².
    """
    dk = k - p.m
    return p.a + p.b * (p.rho * dk + np.sqrt(dk**2 + p.sig**2))


def svi_implied_vol(k: np.ndarray, T: float, p: SVIParams) -> np.ndarray:
    """
    Annualised implied volatility from SVI.

    Convention: continuous compounding, ACT/365, annualised decimal (0.20 = 20%).
    """
    w = svi_total_variance(k, p)
    return np.sqrt(np.maximum(w, 0.0) / T)


def svi_dw_dk(k: np.ndarray, p: SVIParams) -> np.ndarray:
    """First derivative ∂w/∂k (per-slice; T cancels in the SVI slice form)."""
    dk = k - p.m
    sq = np.sqrt(dk**2 + p.sig**2)
    return p.b * (p.rho + dk / sq)


def svi_d2w_dk2(k: np.ndarray, p: SVIParams) -> np.ndarray:
    """Second derivative ∂²w/∂k²."""
    dk = k - p.m
    sq = np.sqrt(dk**2 + p.sig**2)
    return p.b * p.sig**2 / sq**3


def butterfly_density(k: np.ndarray, T: float, p: SVIParams) -> np.ndarray:
    """
    Normalised risk-neutral density g(k) (Lee 2004).

    No butterfly arbitrage iff g(k) ≥ 0 for all k ∈ ℝ.

    g(k) = (1 − k·w′/(2w))² + w″/2 − (w′)²/4·(1/w + 1/4)

    where w, w′, w″ are evaluated at the per-slice total variance.
    """
    w  = svi_total_variance(k, p)
    wp = svi_dw_dk(k, p)
    wpp = svi_d2w_dk2(k, p)
    with np.errstate(divide='ignore', invalid='ignore'):
        term1 = (1 - k * wp / (2 * w)) ** 2
        term2 = wpp / 2
        term3 = wp**2 / 4 * (1.0 / w + 0.25)
    return np.where(w > 0, term1 + term2 - term3, -np.inf)

4.2 Implied Vol Inversion

Given a market option price, we must invert the Black-Scholes formula numerically. Newton-Raphson converges quadratically near the solution but fails when vega is near zero (deep OTM options, very short maturities). Brent's method provides a robust fallback with guaranteed convergence on [1e-6, 1000%].

implied_vol_surface.py — BS inversionPython 3.10+
from scipy import stats
from scipy.optimize import brentq


def bs_price(F: float, K: float, r: float, T: float,
             sigma: float, option_type: str = "call") -> float:
    """
    Black-Scholes price under risk-neutral measure Q.

    Assumptions:
        - dF/F = σ dW^Q  (log-normal forward dynamics)
        - r: continuously compounded risk-free rate (annualised, decimal)
        - T: time to expiry (ACT/365, years)
        - European exercise only

    Limitations:
        - Constant σ: cannot capture smile or skew.
        - Log-normal: no jump risk, no fat tails.
    """
    if T <= 0 or sigma <= 0:
        intrinsic = max(0.0, (F - K) if option_type == "call" else (K - F))
        return np.exp(-r * T) * intrinsic

    sqrt_T = np.sqrt(T)
    d1 = (np.log(F / K) + 0.5 * sigma**2 * T) / (sigma * sqrt_T)
    d2 = d1 - sigma * sqrt_T
    N = stats.norm.cdf
    disc = np.exp(-r * T)

    if option_type == "call":
        return disc * (F * N(d1) - K * N(d2))
    return disc * (K * N(-d2) - F * N(-d1))


def bs_vega(F: float, K: float, r: float, T: float, sigma: float) -> float:
    """∂Price/∂σ — used as the Newton-Raphson Jacobian."""
    if T <= 0 or sigma <= 0:
        return 0.0
    d1 = (np.log(F / K) + 0.5 * sigma**2 * T) / (sigma * np.sqrt(T))
    return np.exp(-r * T) * F * stats.norm.pdf(d1) * np.sqrt(T)


def implied_vol(
    market_price: float,
    F: float,
    K: float,
    r: float,
    T: float,
    option_type: str = "call",
    tol: float = 1e-8,
    max_iter: int = 100,
) -> float:
    """
    Newton-Raphson implied volatility inversion with Brent fallback.

    Update rule: σ_{n+1} = σ_n − (BS(σ_n) − C_mkt) / vega(σ_n)

    Falls back to Brent's method when vega < 1e-12 (deep OTM,
    very short expiry, or ATM for very small T).

    Raises:
        ValueError: price outside no-arbitrage bounds.
    """
    disc = np.exp(-r * T)
    intrinsic = max(0.0, disc * (F - K if option_type == "call" else K - F))
    max_price = disc * (F if option_type == "call" else K)

    if market_price < intrinsic - 1e-10:
        raise ValueError(f"Price {market_price:.6f} below intrinsic {intrinsic:.6f}")
    if market_price > max_price + 1e-10:
        raise ValueError(f"Price {market_price:.6f} above no-arbitrage bound {max_price:.6f}")

    sigma = 0.20  # Initial guess: 20% vol
    for _ in range(max_iter):
        price = bs_price(F, K, r, T, sigma, option_type)
        vega  = bs_vega(F, K, r, T, sigma)
        error = price - market_price
        if abs(error) < tol:
            return sigma
        if abs(vega) < 1e-12:
            break
        sigma -= error / vega
        sigma = np.clip(sigma, 1e-6, 10.0)

    # Brent fallback
    obj = lambda s: bs_price(F, K, r, T, s, option_type) - market_price
    return brentq(obj, 1e-6, 10.0, xtol=tol, maxiter=200)

4.3 SVI Calibration (Levenberg-Marquardt)

Calibration minimises the weighted squared IV residuals over the five SVI parameters. Levenberg-Marquardt (via scipy.optimize.least_squares) handles the non-linear, bounded problem robustly. Multiple random restarts guard against local minima — a real risk when the smile data is sparse or noisy.

implied_vol_surface.py — calibrationPython 3.10+
from scipy.optimize import least_squares


def calibrate_svi(
    k_mkt: np.ndarray,
    iv_mkt: np.ndarray,
    T: float,
    weights: np.ndarray | None = None,
    n_restarts: int = 5,
) -> tuple[SVIParams, float]:
    """
    Calibrate SVI parameters to a single maturity slice.

    Method: Levenberg-Marquardt least squares on weighted IV residuals.

    Minimises:  Σ_i  w_i · (σ_SVI(k_i, T; θ) − σ_mkt_i)²

    Constraints (enforced via bounded optimisation):
        b ≥ 0,  |ρ| < 1,  σ > 0,  a + b·σ·√(1−ρ²) ≥ 0

    Args:
        k_mkt:     Log-moneyness array  k = log(K/F)
        iv_mkt:    Market implied vols   (annualised, decimal)
        T:         Time to maturity       (ACT/365, years)
        weights:   Per-strike weights     (default: uniform)
        n_restarts: Number of random restarts to avoid local minima

    Returns:
        (calibrated_params, root_mean_square_error_in_vol_points)
    """
    if weights is None:
        weights = np.ones_like(k_mkt)

    # Estimate ATM vol for initial guess scaling
    atm_var = float(np.interp(0.0, k_mkt, iv_mkt**2)) * T

    def residuals(x: np.ndarray) -> np.ndarray:
        p = SVIParams(a=x[0], b=x[1], rho=x[2], m=x[3], sig=x[4])
        w = svi_total_variance(k_mkt, p)
        iv_svi = np.sqrt(np.maximum(w, 1e-12) / T)
        return weights * (iv_svi - iv_mkt)

    bounds = (
        [1e-6,  1e-4, -0.999, -0.5,  1e-4],   # lower
        [0.5,   2.0,   0.999,  0.5,  1.0 ],   # upper
    )

    best_result = None
    rng = np.random.default_rng(seed=42)  # Fixed seed: reproducibility

    for i in range(n_restarts):
        if i == 0:
            x0 = [atm_var * 0.8, 0.15, -0.3, 0.0, 0.25]
        else:
            x0 = [
                rng.uniform(1e-4, 0.3),
                rng.uniform(0.05, 0.5),
                rng.uniform(-0.9, 0.9),
                rng.uniform(-0.3, 0.3),
                rng.uniform(0.05, 0.5),
            ]
        try:
            result = least_squares(residuals, x0, bounds=bounds,
                                   method="trf", ftol=1e-10, xtol=1e-10)
            if best_result is None or result.cost < best_result.cost:
                best_result = result
        except Exception:
            continue

    if best_result is None:
        raise RuntimeError("SVI calibration failed on all restarts")

    p = SVIParams(*best_result.x)
    rmse = float(np.sqrt(np.mean(residuals(best_result.x)**2)))
    return p, rmse

4.4 Plotly Visualisation

The visualisation module renders the calibrated surface as an interactive Plotly 3-D plot. Run locally or export to surface.html for standalone sharing.

visualisation.py — Plotly 3D surface and arbitrage plotPython 3.10+
import plotly.graph_objects as go
from plotly.subplots import make_subplots


def plot_svi_surface(
    T_grid: np.ndarray,
    k_grid: np.ndarray,
    params_per_slice: list[SVIParams],
) -> go.Figure:
    """
    Render the calibrated implied vol surface as an interactive 3-D plot.

    Args:
        T_grid:           Maturity grid (years)
        k_grid:           Log-moneyness grid
        params_per_slice: Calibrated SVIParams for each maturity

    Returns:
        Plotly Figure — call .show() or .write_html("surface.html")
    """
    K, T = np.meshgrid(k_grid, T_grid)
    Z = np.zeros_like(K)
    for i, (t, p) in enumerate(zip(T_grid, params_per_slice)):
        Z[i, :] = svi_implied_vol(k_grid, t, p) * 100  # percent

    fig = make_subplots(
        rows=1, cols=2,
        specs=[[{"type": "surface"}, {"type": "scatter"}]],
        subplot_titles=["Implied Vol Surface", "Per-Slice Smiles"],
        column_widths=[0.6, 0.4],
    )

    # 3-D surface
    fig.add_trace(
        go.Surface(
            x=K, y=T, z=Z,
            colorscale="RdYlBu_r",
            colorbar=dict(title="σ_IV (%)", x=0.55),
            name="SVI Surface",
        ),
        row=1, col=1,
    )

    # 2-D smile overlays
    colors = [f"hsl({int(240 - i * 240 / len(T_grid))},70%,50%)" for i in range(len(T_grid))]
    for i, (t, p, c) in enumerate(zip(T_grid, params_per_slice, colors)):
        iv = svi_implied_vol(k_grid, t, p) * 100
        fig.add_trace(
            go.Scatter(
                x=k_grid, y=iv,
                mode="lines",
                line=dict(color=c, width=2),
                name=f"T={t:.2f}y",
            ),
            row=1, col=2,
        )

    fig.update_layout(
        title="SVI Calibrated Implied Volatility Surface",
        scene=dict(
            xaxis_title="k = log(K/F)",
            yaxis_title="Maturity T (years)",
            zaxis_title="σ_IV (%)",
            camera=dict(eye=dict(x=1.5, y=-1.8, z=0.9)),
        ),
        xaxis2_title="k = log(K/F)",
        yaxis2_title="σ_IV (%)",
        height=500,
        template="plotly_dark",
    )
    return fig


def plot_arbitrage_check(k_grid: np.ndarray, T: float, p: SVIParams) -> go.Figure:
    """
    Plot butterfly density g(k) and calendar spread condition ∂w/∂T for a given slice.
    Violations (g < 0 or ∂w/∂T < 0) are highlighted in red.
    """
    g   = butterfly_density(k_grid, T, p)
    dk  = k_grid - p.m
    sq  = np.sqrt(dk**2 + p.sig**2)
    inner = p.rho * dk + sq
    cal = p.a + p.b / (2 * np.sqrt(T)) * inner   # ∂w/∂T

    fig = make_subplots(rows=1, cols=2,
                        subplot_titles=["Butterfly density g(k)", "Calendar spread ∂w/∂T"])

    for col, (y, label) in enumerate([(g, "g(k)"), (cal, "∂w/∂T")], start=1):
        color = np.where(y < 0, "red", "royalblue")
        fig.add_trace(
            go.Scatter(x=k_grid, y=y, mode="lines",
                       line=dict(width=2, color="royalblue"), name=label),
            row=1, col=col,
        )
        fig.add_hline(y=0, line_dash="dash", line_color="white",
                      line_width=1, row=1, col=col)

    fig.update_layout(height=300, template="plotly_dark",
                      title=f"Arbitrage conditions at T={T:.2f}y")
    return fig

5. Arbitrage Conditions

5.1 Butterfly Arbitrage

A surface is free of butterfly arbitrage iff the implied risk-neutral density is non-negative everywhere. Lee (2004) shows this is equivalent to:

g(k)=(1kw2w)2+w2(w)24(1w+14)0kRg(k) = \left(1 - \frac{k w'}{2w}\right)^2 + \frac{w''}{2} - \frac{(w')^2}{4}\left(\frac{1}{w} + \frac{1}{4}\right) \geq 0 \quad \forall\, k \in \mathbb{R}

For the SVI parametrisation (single slice):

w(k)=b(ρ+km(km)2+σ2),w(k)=bσ2((km)2+σ2)3/2w'(k) = b\left(\rho + \frac{k-m}{\sqrt{(k-m)^2+\sigma^2}}\right), \qquad w''(k) = \frac{b\sigma^2}{\left((k-m)^2+\sigma^2\right)^{3/2}}

Note that w>0w'' > 0 always (SVI is always convex in kk), so violations of g0g \geq 0 arise from the ww' term dominating — typically with extreme skew (ρ|\rho| near 1) combined with small σ\sigma.

5.2 Calendar Spread Arbitrage

For a surface parametrised across maturities, the total implied variance must be non-decreasing in TT at every log-moneyness level:

w(k,T)T0kR,  T>0\frac{\partial w(k,T)}{\partial T} \geq 0 \quad \forall\, k \in \mathbb{R},\; T > 0

Violation means a longer-dated option is cheaper than a shorter-dated one — a static arbitrage exploitable by buying the near expiry and selling the far expiry.

For the SSVI-style surface w(k,T)=aT+bTϕ(k)w(k,T) = aT + b\sqrt{T}\cdot\phi(k) where ϕ(k)=ρ(km)+(km)2+σ2\phi(k) = \rho(k-m)+\sqrt{(k-m)^2+\sigma^2}:

wT=a+b2Tϕ(k)\frac{\partial w}{\partial T} = a + \frac{b}{2\sqrt{T}}\,\phi(k)

Since ϕ(k)σ>0\phi(k) \geq \sigma > 0 for ρ=0\rho = 0, but ϕ(k)\phi(k) can be negative for ρ<0\rho < 0 and large negative kk (deep OTM puts). Violations occur when aa is small and b/Tb/\sqrt{T} is large with a negatively skewed left wing.

Joint Arbitrage Across the Full Surface

The per-slice SVI parameters (a₁,b₁,ρ₁,m₁,σ₁) for T₁ and (a₂,b₂,ρ₂,m₂,σ₂) for T₂ can each be individually butterfly-free while still admitting a calendar spread arbitrage. Checking and enforcing joint no-arbitrage requires the SSVI (Gatheral & Jacquier 2014) or eSSVI (Corbetta et al. 2019) framework, which imposes parametric constraints on how SVI parameters can vary with T.

6. Test Suite

Tests are organised by concern: SVI kernel correctness, Black-Scholes inversion round-trips, arbitrage condition checks, and calibration integration tests. Run with pytest tests/test_svi.py -v.

tests/test_svi.py — full pytest suitePython 3.10+
# tests/test_svi.py
# Run: pytest tests/test_svi.py -v
import numpy as np
import pytest
from implied_vol_surface import (
    SVIParams, svi_total_variance, svi_implied_vol,
    butterfly_density, calibrate_svi, implied_vol, bs_price,
)

# ── Fixtures ─────────────────────────────────────────────────────────────────

REFERENCE_PARAMS = SVIParams(a=0.04, b=0.15, rho=-0.4, m=0.0, sig=0.20)
K_GRID = np.linspace(-0.4, 0.4, 50)
T_SLICE = 1.0


# ── Unit tests: SVI kernel ────────────────────────────────────────────────────

class TestSVIKernel:
    def test_atm_variance(self):
        """At k=0 and m=0, w(0) = a + b·σ (since ρ·0 + √σ² = σ)."""
        p = REFERENCE_PARAMS
        expected = p.a + p.b * p.sig
        assert abs(svi_total_variance(np.array([0.0]), p)[0] - expected) < 1e-12

    def test_total_variance_positive(self):
        """No negative total variance for standard parameters."""
        w = svi_total_variance(K_GRID, REFERENCE_PARAMS)
        assert np.all(w > 0), "Total variance must be positive everywhere"

    def test_skew_direction(self):
        """Negative ρ ⟹ left wing has higher vol than right wing (equity skew)."""
        p = REFERENCE_PARAMS  # rho = -0.4
        iv = svi_implied_vol(K_GRID, T_SLICE, p)
        assert iv[0] > iv[-1], "Negative skew: put wing > call wing"

    def test_symmetric_for_zero_rho_and_m(self):
        """Zero skew and shift ⟹ symmetric smile."""
        p = SVIParams(a=0.04, b=0.15, rho=0.0, m=0.0, sig=0.20)
        iv = svi_implied_vol(K_GRID, T_SLICE, p)
        # iv(k) should equal iv(-k) up to floating-point tolerance
        iv_rev = svi_implied_vol(-K_GRID, T_SLICE, p)
        np.testing.assert_allclose(iv, iv_rev, atol=1e-12)

    def test_vol_decreases_with_maturity_at_fixed_params(self):
        """
        For fixed SVI params, short-dated slices have higher ATM vol
        because σ_IV = √(w/T) and w is the same across maturities.
        """
        p = REFERENCE_PARAMS
        iv_short = svi_implied_vol(np.array([0.0]), 0.25, p)[0]
        iv_long  = svi_implied_vol(np.array([0.0]), 2.00, p)[0]
        assert iv_short > iv_long


# ── Unit tests: Black-Scholes inversion ──────────────────────────────────────

class TestImpliedVolInversion:
    @pytest.mark.parametrize("sigma_true", [0.10, 0.20, 0.35, 0.60])
    def test_roundtrip_call(self, sigma_true: float):
        """Price at σ_true → invert → should recover σ_true."""
        F, K, r, T = 100.0, 100.0, 0.05, 1.0
        price = bs_price(F, K, r, T, sigma_true, "call")
        sigma_inv = implied_vol(price, F, K, r, T, "call")
        assert abs(sigma_inv - sigma_true) < 1e-7, (
            f"Round-trip failed: true={sigma_true:.4f}, inverted={sigma_inv:.4f}"
        )

    @pytest.mark.parametrize("moneyness", [-0.3, -0.1, 0.0, 0.1, 0.3])
    def test_roundtrip_put_otm(self, moneyness: float):
        """Test OTM puts across a range of moneyness levels."""
        F, r, T, sigma_true = 100.0, 0.05, 0.5, 0.25
        K = F * np.exp(moneyness)
        price = bs_price(F, K, r, T, sigma_true, "put")
        sigma_inv = implied_vol(price, F, K, r, T, "put")
        assert abs(sigma_inv - sigma_true) < 1e-6

    def test_put_call_parity(self):
        """Call − Put = (F − K)·e^{−rT}: must hold to machine precision."""
        F, K, r, T, sigma = 100.0, 95.0, 0.05, 1.0, 0.20
        call = bs_price(F, K, r, T, sigma, "call")
        put  = bs_price(F, K, r, T, sigma, "put")
        parity = (F - K) * np.exp(-r * T)
        assert abs(call - put - parity) < 1e-10


# ── Unit tests: Arbitrage detection ──────────────────────────────────────────

class TestArbitrage:
    def test_reference_params_butterfly_free(self):
        """Reference parameters should produce no butterfly arbitrage."""
        g = butterfly_density(K_GRID, T_SLICE, REFERENCE_PARAMS)
        assert np.all(g >= -1e-9), f"Butterfly violation: min g = {g.min():.6f}"

    def test_extreme_b_causes_butterfly_arbitrage(self):
        """Very large b with extreme skew will violate the butterfly condition."""
        p = SVIParams(a=0.001, b=0.8, rho=-0.9, m=0.0, sig=0.05)
        g = butterfly_density(K_GRID, T_SLICE, p)
        assert np.any(g < 0), "Expected butterfly violation for extreme parameters"

    def test_high_b_no_arb_at_larger_sig(self):
        """Larger σ (ATM curvature) regularises the smile — reduces violations."""
        p_bad  = SVIParams(a=0.01, b=0.5, rho=-0.7, m=0.0, sig=0.05)
        p_good = SVIParams(a=0.01, b=0.5, rho=-0.7, m=0.0, sig=0.40)
        g_bad  = butterfly_density(K_GRID, T_SLICE, p_bad)
        g_good = butterfly_density(K_GRID, T_SLICE, p_good)
        # Good params should have fewer (or no) violations
        assert g_good.min() > g_bad.min()


# ── Integration test: calibration roundtrip ──────────────────────────────────

class TestCalibration:
    def test_calibration_recovers_params(self):
        """
        Generate market IVs from known SVI params, calibrate, and verify
        the calibrated params reproduce the original surface within tolerance.
        """
        true_params = SVIParams(a=0.06, b=0.20, rho=-0.3, m=0.05, sig=0.15)
        k_mkt = np.linspace(-0.35, 0.35, 15)
        iv_mkt = svi_implied_vol(k_mkt, T_SLICE, true_params)

        fitted_params, rmse = calibrate_svi(k_mkt, iv_mkt, T_SLICE)

        assert rmse < 1e-6, f"Calibration RMSE too large: {rmse:.2e}"

        # Verify surface (not just params — SVI params can be non-unique)
        iv_fit = svi_implied_vol(k_mkt, T_SLICE, fitted_params)
        np.testing.assert_allclose(iv_fit, iv_mkt, atol=1e-5,
            err_msg="Calibrated surface does not match market IVs")

    def test_calibration_rmse_units(self):
        """RMSE should be in vol units (decimal), not basis points."""
        true_params = REFERENCE_PARAMS
        k_mkt = np.linspace(-0.3, 0.3, 10)
        iv_mkt = svi_implied_vol(k_mkt, T_SLICE, true_params)
        # Add small noise (5bp)
        rng = np.random.default_rng(0)
        iv_noisy = iv_mkt + rng.normal(0, 0.0005, size=iv_mkt.shape)

        _, rmse = calibrate_svi(k_mkt, iv_noisy, T_SLICE)
        # RMSE should be roughly 5bp = 0.0005 in vol decimal
        assert rmse < 0.005, "RMSE unexpectedly large (check units)"

7. Limitations and Failure Modes

Where SVI Breaks Down

  • Non-uniqueness: Different parameter sets can produce identical smiles. This is known as the SVI parameter redundancy problem and can cause the calibration to return unstable parameters even when the fit is perfect.
  • Per-slice model: SVI has no built-in constraint linking parameters across maturities. Calendar spread arbitrage must be checked and enforced separately — it is not guaranteed by per-slice calibration.
  • Extreme strike extrapolation: SVI's linear large-strike asymptotics match Roger Lee's moment formula in form, but the specific slope depends on b and ρ. Extrapolating far beyond observed strikes should be treated with extreme caution.
  • Short-dated smiles: For very short maturities (T < 0.05 years), the smile can be very steep and SVI may require a very large b to fit. This often triggers butterfly violations. Alternative parametrisations (e.g., SABR, jump-diffusion) may be more appropriate.
  • Local minima: The LM objective is non-convex. Multiple restarts are necessary, and even then the global optimum is not guaranteed. Regularisation or Bayesian priors on parameters can help.
  • No model dynamics: SVI is a calibration tool, not a pricing model. It cannot be used to price cliquets, forward-starting options, or other path-dependent payoffs without an accompanying stochastic vol model consistent with the SVI surface.

8. Interview Angle

Vol surface calibration appears in quant interviews at every level. Junior candidates are expected to understand the smile and basic SVI geometry. Senior candidates must derive the no-arbitrage conditions and discuss numerical implementation. Researcher-level candidates are expected to critique SVI, discuss its extensions, and connect it to local vol and stochastic vol models.

L1 — Junior

  • What is the volatility smile and why does it exist under real market dynamics?
  • What is implied volatility and how do you extract it from an option price?
  • What are the five SVI parameters and what does each one control geometrically?
  • What is log-moneyness and why do we prefer it over strike/spot ratios?

L2 — Senior

  • Derive the butterfly arbitrage condition g(k) ≥ 0 from the no-negative-density requirement.
  • Explain why Newton-Raphson for implied vol inversion can fail and how you handle it.
  • What is the calendar spread no-arbitrage condition on total variance? How does it constrain joint calibration across maturities?
  • How does the Levenberg-Marquardt algorithm handle the ill-conditioned Jacobian that arises when calibrating SVI to a sparse grid of strikes?
  • Why is calibrating SVI in total variance space (w = σ²T) numerically preferable to calibrating in implied vol space?

L3 — Researcher

  • SVI has a known non-uniqueness problem: different parameter sets can produce the same smile. Describe it and explain how SSVI (Surface SVI) partially resolves it.
  • Gatheral & Jacquier's eSSVI guarantees joint no-arbitrage across the full surface. Derive the constraints it imposes on the ρ and θ functions.
  • The Dupire local vol formula involves ∂²C/∂K² and ∂C/∂T. How would you numerically stabilise both derivatives when computed from a calibrated SVI surface?
  • When would you prefer a non-parametric approach (e.g., SSVI with cubic spline θ) over raw SVI? What are the regularisation requirements?