vol_surface.py

"""Implied-volatility surface via slice-by-slice SVI calibration (Gatheral)

The surface stores one calibrated SVI slice per maturity. Total variance is
SVI in k within a slice and linearly interpolated in tau between slices, which
preserves calendar monotonicity and hence no calendar arbitrage
"""

from dataclasses import astuple, dataclass

import numpy as np
from scipy.optimize import least_squares

from market_data import MarketSnapshot


@dataclass(frozen=True)
class SVIParams:
    """Raw SVI parameters for a single maturity slice

    Total variance: w(k) = a + b * (rho * (k - m) + sqrt((k - m)^2 + sigma^2))
    """

    a: float
    b: float
    rho: float
    m: float
    sigma: float


def svi_total_variance(k: np.ndarray, params: SVIParams) -> np.ndarray:
    """Evaluate the SVI total variance w(k) for a slice"""
    a, b, rho, m, sigma = astuple(params)
    centred = k - m
    return a + b * (rho * centred + np.sqrt(centred**2 + sigma**2))


def fit_svi_slice(k: np.ndarray, w: np.ndarray) -> SVIParams:
    """Calibrate SVI parameters on one maturity slice by least squares"""
    eps = 1e-6

    # data-driven initial guess
    a0 = 0.5 * float(np.min(w))
    m0 = float(k[np.argmin(w)])
    sigma0 = 0.1
    b0 = 0.1
    rho0 = -0.5
    p0 = [a0, b0, rho0, m0, sigma0]

    lower = [0.0, 0.0, -1 + eps, -np.inf, eps]
    upper = [np.inf, np.inf, 1 - eps, np.inf, np.inf]

    def residual(p: np.ndarray) -> np.ndarray:
        return svi_total_variance(k, SVIParams(*p)) - w

    result = least_squares(residual, p0, bounds=(lower, upper))
    return SVIParams(*result.x)


class VolSurface:
    """Calibrated implied-volatility surface (one SVI slice per maturity)"""

    def __init__(self, taus: np.ndarray, slices: list[SVIParams]) -> None:
        """Store maturities and their SVI parameters, sorted by increasing tau"""
        order = np.argsort(taus)
        self.taus = np.asarray(taus)[order]
        self.slices = [slices[i] for i in order]

    def total_variance(self, k: float, tau: float) -> float:
        """Total variance w(k, tau) with linear-in-tau interpolation

        Three regimes:
        - tau <= tau_min : interpolate between the origin (0, 0) and the first
          slice (vol extrapolation towards expiry)
        - tau_min < tau < tau_max : linear interpolation in tau of the total
          variance evaluated at k on the two bracketing slices
        - tau >= tau_max : vol extrapolation, w grows proportionally to tau
        """
        taus = self.taus

        # short end: interpolate from (0, 0) to the first slice
        if tau <= taus[0]:
            w1 = svi_total_variance(k, self.slices[0])
            return float((tau / taus[0]) * w1)

        # long end:implied vol -> w proportional to tau
        if tau >= taus[-1]:
            wn = svi_total_variance(k, self.slices[-1])
            return float((tau / taus[-1]) * wn)

        # interior: bracket tau and interpolate the total variance linearly
        j = int(np.searchsorted(taus, tau)) - 1
        t_lo, t_hi = taus[j], taus[j + 1]
        w_lo = svi_total_variance(k, self.slices[j])
        w_hi = svi_total_variance(k, self.slices[j + 1])
        weight = (tau - t_lo) / (t_hi - t_lo)
        return float(w_lo + weight * (w_hi - w_lo))

    def implied_vol(self, k: float, tau: float) -> float:
        """Implied volatility sigma = sqrt(w / tau)"""
        return float(np.sqrt(self.total_variance(k, tau) / tau))


def otm_mask(snapshot: MarketSnapshot) -> np.ndarray:
    """Keep out-of-the-money quotes: puts for k < 0, calls for k > 0"""
    is_put = snapshot.option_type == "put"
    is_call = snapshot.option_type == "call"
    k = snapshot.log_moneyness
    return (is_put & (k < 0)) | (is_call & (k > 0))


def build_surface(
    snapshot: MarketSnapshot,
    iv: np.ndarray,
    k_max: float = 0.25,
    tau_min: float = 0.1,
    tau_max: float = 1.5,
    min_points: int = 10,
    rmse_max: float = 0.06,
) -> VolSurface:
    """Calibrate one SVI slice per maturity and assemble the surface

    Parameters -->
    snapshot : MarketSnapshot
        Clean market data
    iv : np.ndarray
        Implied vols aligned with ``snapshot`` (nan where inversion failed)
    k_max : float
        Keep only quotes with |log-moneyness| < k_max for the fit
    tau_min, tau_max : float
        Maturity window where SVI is trusted. Very short maturities are
        microstructure-dominated and too peaked for raw SVI; very long ones
        are sparsely quoted. Outside this window the surface extrapolates
    min_points : int
        Skip a maturity with fewer than this many usable quotes
    rmse_max : float
        Skip a slice whose relative fit RMSE exceeds this (pathological fit)
    """
    base = (
        ~np.isnan(iv)
        & (np.abs(snapshot.log_moneyness) < k_max)
        & (snapshot.tau >= tau_min)
        & (snapshot.tau <= tau_max)
        & otm_mask(snapshot)
    )

    taus, slices = [], []
    for tau in np.unique(snapshot.tau[base]):
        sl = base & (snapshot.tau == tau)
        if int(sl.sum()) < min_points:
            print(f"skip tau={tau:.3f}: only {int(sl.sum())} points")
            continue

        k = snapshot.log_moneyness[sl]
        w = iv[sl] ** 2 * tau
        order = np.argsort(k)
        k, w = k[order], w[order]
        params = fit_svi_slice(k, w)

        # quality gate: reject pathological fits
        rmse_rel = np.sqrt(np.mean((svi_total_variance(k, params) - w) ** 2)) / np.mean(w)
        if rmse_rel > rmse_max:
            print(f"skip tau={tau:.3f}: fit RMSE {rmse_rel:.1%} > {rmse_max:.0%}")
            continue

        taus.append(float(tau))
        slices.append(params)

    return VolSurface(np.array(taus), slices)