# src/pysvi/surface.py
"""Fitted volatility surface: evaluation, interpolation, and pricing.
`VolSurface` turns per-slice calibration results into the object quant
work actually consumes: model -> calibration -> fitted surface. It owns
calibrated slices across maturities and exposes vectorized evaluation
(IVs, total variance, ATM level/skew/curvature), maturity interpolation
between slices, arbitrage verification, and a Black-76 pricing and
Greeks layer on the slice forwards.
`calibrate_surface` is the calendar-aware fitter: it orders expiries,
derives per-slice inputs, chains the prior slice's total variance into
each NO_CALENDAR penalty automatically, and — for eSSVI — fits the
global term structure jointly across all slices.
Conventions
-----------
* Pricing is Black-76 on the slice forward with a flat continuously
compounded rate ``r`` (default 0): ``C = e^{-rT}[F N(d1) - K N(d2)]``.
* Greeks hold the implied volatility fixed (sticky-strike): delta and
gamma are with respect to the forward, vega is per unit volatility,
theta is per year of calendar time.
* Between fitted maturities the surface interpolates (linearly in total
variance at fixed log-moneyness by default; see ``interp_method``);
forwards interpolate log-linearly. Extrapolation beyond the fitted
maturity range raises.
"""
from bisect import bisect_left
from dataclasses import dataclass
from typing import Dict, Iterable, Mapping, Optional, Tuple, Union
import numpy as np
from loguru import logger
from scipy.special import ndtr
from . import _kernels
from .models import (
ArbitrageFreedom, DirectSVI, ESSVI, JumpWings, NaturalSVI, Parametrization,
SABR, SSVI, SVI,
_initialization, _mad_scale, _minimize_with_starts, _multistart_variants,
_penalty_grid, _prepare_loss_inputs, essvi_total_variance,
)
from .calibration import calibrate_slice, get_model, prepare_slice
from .diagnostics import ArbitrageReport, check_arbitrage
from .report import (
SLICE_FAILED, SLICE_INSUFFICIENT, SLICE_OK,
SurfaceDiagnostics, SurfaceFitReport,
build_slice_report, build_surface_report, validate_mode,
)
_SQRT_2PI = np.sqrt(2.0 * np.pi)
_INTERP_METHODS = ("total_variance", "theta", "monotone_cubic")
#: Smoothness in maturity guaranteed by each interpolation method.
_REGULARITY = {"total_variance": "C0", "theta": "C0", "monotone_cubic": "C1"}
#: Serialization schema version written by VolSurface.save.
_SCHEMA_VERSION = 1
_MODEL_CLASSES = {
cls.__name__: cls
for cls in (SVI, NaturalSVI, SSVI, ESSVI, JumpWings, DirectSVI, SABR)
}
def _npdf(x):
return np.exp(-0.5 * x * x) / _SQRT_2PI
def _is_call(cp: str) -> bool:
flag = str(cp).lower()
if flag in ("call", "c"):
return True
if flag in ("put", "p"):
return False
raise ValueError(f"cp must be 'call' or 'put', got {cp!r}")
def _shape_like(values, original):
"""Return a float for scalar input, the array otherwise."""
return float(values[0]) if np.ndim(original) == 0 else values
[docs]
@dataclass(frozen=True)
class VarianceEvent:
"""A discrete variance event (earnings, FOMC, CPI, election).
Real term structures contain known jumps: generic maturity
interpolation smooths across them, silently asserting all
term-structure curvature is continuous variance. With events, total
variance decomposes as::
w(k, T) = w_cont(k, T) + sum of event variances with t_e <= T
Fitting subtracts each expiry's cumulative event variance before
calibration (slices store the CONTINUOUS component), interpolation
acts on the continuous component, and evaluation adds the events
back on the correct side of each event time -- so the surface
reproduces the jump exactly instead of smearing it, and the
calendar diagnostics (which see the continuous slices) raise no
false violation across an event.
Attributes
----------
time : float
Event time as a year fraction (same clock as maturities).
variance : float
Total variance added by the event (ATM units, e.g.
sigma_event^2 * dt), non-negative.
label : str
Optional tag ("AAPL earnings", "FOMC").
"""
time: float
variance: float
label: str = ""
def __post_init__(self):
if not (self.time > 0 and np.isfinite(self.time)):
raise ValueError(f"event time must be positive, got {self.time}")
if not (self.variance >= 0 and np.isfinite(self.variance)):
raise ValueError(
f"event variance must be non-negative, got {self.variance}"
)
def _cum_event_variance(events, T) -> float:
"""Total event variance realized by expiry T (events at t_e <= T)."""
return float(sum(e.variance for e in events if e.time <= T + 1e-12))
def _atm_theta(g, event_var: float = 0.0) -> float:
"""ATM total variance of a slice: w interpolated at k = 0.
SSVI fixes w(0) = theta exactly, so theta must be the ATM level --
the smile MINIMUM sits away from k = 0 on any skewed smile and
systematically understates the ATM vol, which the (rho, eta) fit
can never correct. Reads the cleaned inputs (prepare_slice) so a
junk quote cannot drive theta toward zero; falls back to the raw
minimum only for slices too thin to clean.
"""
k, w, _ = prepare_slice(g)
if k is None:
return float(np.nanmin(g["iv"] ** 2 * g["maturity"])) - event_var
order = np.argsort(k)
return float(np.interp(0.0, k[order], w[order])) - event_var
def _subtract_events(w, events, T, mode):
"""Continuous total variance: quoted w minus event variance by T.
Returns None (slice unusable) when the event variance exceeds the
quoted total variance anywhere -- the specified events are
inconsistent with the market; strict mode raises instead.
"""
ev = _cum_event_variance(events, T)
if ev == 0.0:
return w
w_cont = w - ev
if np.any(w_cont <= 0):
msg = (
f"slice T={T:g}: cumulative event variance {ev:g} exceeds "
"the quoted total variance at some strikes -- the events "
"are inconsistent with the quoted term structure"
)
if mode == "strict":
raise ValueError(f"(mode='strict') {msg}")
if mode == "warn":
logger.warning(msg)
return None
return w_cont
def _shift_band_events(kwargs, events, T) -> None:
"""Shift bid/ask bands into continuous-variance space alongside w."""
ev = _cum_event_variance(events, T)
if ev and "w_bid" in kwargs:
kwargs["w_bid"] = kwargs["w_bid"] - ev
kwargs["w_ask"] = kwargs["w_ask"] - ev
[docs]
def implied_event_variances(df, events) -> list:
"""What the quoted term structure says about each event's variance.
For every event, finds the quoted expiries straddling it and
reports the jump in ATM total variance across the event
(``atm_w(T_after) - atm_w(T_before)``) next to the variance the
event specifies. The quoted jump also contains the continuous
variance accrued between the two expiries, so it is an UPPER bound
on the event variance; a specified variance far above it is
inconsistent with the market.
Returns a list of dicts: {event, T_before, T_after, quoted_jump,
specified}. Events without straddling quotes report None bounds.
"""
events = tuple(
(e if isinstance(e, VarianceEvent) else VarianceEvent(*e))
for e in events
)
atm = {}
for T, g in df.groupby("maturity"):
atm[float(T)] = _atm_theta(g)
Ts = sorted(atm)
out = []
for e in events:
before = [T for T in Ts if T < e.time]
after = [T for T in Ts if T >= e.time]
if not before or not after:
out.append({"event": e, "T_before": None, "T_after": None,
"quoted_jump": None, "specified": e.variance})
continue
T1, T2 = before[-1], after[0]
out.append({
"event": e, "T_before": T1, "T_after": T2,
"quoted_jump": atm[T2] - atm[T1], "specified": e.variance,
})
return out
def _prior_slice_kwargs(prior, anchor, T, kwargs) -> None:
"""Attach the matching prior slice (by maturity) for temporal anchoring.
``prior`` is a previously fitted VolSurface (or mapping T -> params);
slices are matched by np.isclose on maturity, unmatched slices fit
unanchored. See the calibrate() docs for 'prior'/'anchor'.
"""
if prior is None:
return
items = prior._slices if isinstance(prior, VolSurface) else list(
prior.items() if hasattr(prior, "items") else prior
)
for T_p, params_p in items:
if np.isclose(float(T_p), T, rtol=1e-9, atol=1e-12):
kwargs.setdefault("prior", params_p)
if anchor:
kwargs.setdefault("anchor", float(anchor))
return
def _band_kwargs(g, T, k, w, sel, kwargs) -> None:
"""Derive per-slice w_bid/w_ask for the bid_ask objective from the
panel's iv_bid/iv_ask columns (as OptionChain produces), filtered
exactly like the quotes. Rows whose band is missing or crossed
degenerate to a zero-width band at the mid (fit-to-mid there).
Explicit w_bid/w_ask kwargs win; absent columns leave the kwargs
untouched (the model then raises its usual requirement error).
"""
if kwargs.get("objective") != "bid_ask":
return
if "w_bid" in kwargs or "w_ask" in kwargs:
return
if not ("iv_bid" in g.columns and "iv_ask" in g.columns):
return
iv_bid = g["iv_bid"].to_numpy(dtype=float)[sel]
iv_ask = g["iv_ask"].to_numpy(dtype=float)[sel]
w_bid = iv_bid ** 2 * T
w_ask = iv_ask ** 2 * T
bad = ~np.isfinite(w_bid) | ~np.isfinite(w_ask) | (w_bid > w_ask)
w_bid[bad] = w[bad]
w_ask[bad] = w[bad]
kwargs["w_bid"] = w_bid
kwargs["w_ask"] = w_ask
def _auto_slice_kwargs(instance, T, df_slice, model_kwargs, theta_by_T, theta_ref):
"""Derive the per-slice calibrate kwargs for a model instance."""
kwargs = dict(model_kwargs)
if isinstance(instance, ESSVI):
kwargs["theta"] = theta_by_T[T]
kwargs.setdefault("theta_ref", theta_ref)
elif isinstance(instance, SSVI):
kwargs["theta"] = theta_by_T[T]
elif isinstance(instance, JumpWings):
kwargs["T"] = T
elif isinstance(instance, SABR):
kwargs["T"] = T
kwargs["F"] = float(df_slice["implied_forward"].iloc[0])
kwargs.setdefault("beta", 0.5)
return kwargs
def _data_thetas(instance, groups, events=()):
"""Per-slice CONTINUOUS ATM total variance for SSVI/eSSVI (event
variance subtracted), plus the median ref."""
theta_by_T: Dict[float, float] = {}
theta_ref = None
if isinstance(instance, (SSVI, ESSVI)):
for T, g in groups:
theta_by_T[T] = _atm_theta(g, _cum_event_variance(events, T))
# One junk expiry (all-NaN ivs, flat garbage) must not poison
# every other slice: the reference is the median of the FINITE,
# positive per-slice thetas only. Slices whose own theta is bad
# fail individually and are recorded.
good = [v for v in theta_by_T.values() if np.isfinite(v) and v > 0]
theta_ref = float(np.median(good)) if good else None
return theta_by_T, theta_ref
[docs]
class VolSurface:
"""A fitted implied-volatility surface across maturities.
Construct via :meth:`fit` (independent per-slice calibration),
:func:`calibrate_surface` (calendar-aware), or directly from
calibrated slices::
surface = VolSurface.fit(df, model="svi")
surface = VolSurface(model, {0.25: params_1, 0.5: params_2})
Direct construction requires each params dict to carry ``'forward'``
(as returned by ``calibrate_slice``).
Parameters
----------
model : Parametrization
The model instance all slices were calibrated with.
slices : mapping or iterable of (maturity, params)
Calibrated slices; sorted by maturity internally.
r : float, default 0.0
Flat continuously compounded discount rate used by the pricing
layer.
interp_method : str, default "total_variance"
Maturity interpolation between fitted slices. "total_variance"
blends w(k) linearly in T at fixed log-moneyness (model-agnostic;
calendar-free whenever the bracketing slices are). "theta"
(SSVI/eSSVI only) interpolates the ATM total variance and shape
parameters, yielding a genuine parametric slice at any maturity.
"monotone_cubic" is a shape-preserving cubic (PCHIP,
Fritsch-Carlson) in maturity at fixed log-moneyness across ALL
fitted slices: continuously differentiable in T (C1, see
:attr:`regularity` and :meth:`dw_dT`), exact at fitted
maturities, and monotone in T wherever the fitted slices are --
so calendar-free slices stay calendar-free between expiries.
Requires at least two fitted slices.
fit_report : SurfaceFitReport, optional
Calibration provenance and per-slice evidence; populated by
:meth:`fit` and :func:`calibrate_surface`, None on direct
construction.
"""
def __init__(
self,
model: Parametrization,
slices: Union[
Mapping[float, Dict[str, float]],
Iterable[Tuple[float, Dict[str, float]]],
],
r: float = 0.0,
interp_method: str = "total_variance",
fit_report: Optional["SurfaceFitReport"] = None,
events=(),
) -> None:
if isinstance(slices, Mapping):
slices = slices.items()
ordered = sorted(
((float(T), dict(params)) for T, params in slices),
key=lambda item: item[0],
)
if not ordered:
raise ValueError("VolSurface requires at least one calibrated slice")
maturities = [T for T, _ in ordered]
if len(set(maturities)) != len(maturities):
raise ValueError(f"duplicate maturities in slices: {maturities}")
for T, params in ordered:
if T <= 0:
raise ValueError(f"maturities must be positive, got {T}")
forward = params.get("forward")
if forward is None or not np.isfinite(forward) or forward <= 0:
raise ValueError(
f"slice T={T:g} lacks a positive 'forward' entry; "
"calibrate_slice adds it, direct construction must supply it"
)
if interp_method not in _INTERP_METHODS:
raise ValueError(
f"unknown interp_method {interp_method!r}; choose from {_INTERP_METHODS}"
)
if interp_method == "theta" and not isinstance(model, (SSVI, ESSVI)):
raise ValueError(
"interp_method='theta' requires an SSVI or eSSVI model"
)
if interp_method == "monotone_cubic" and len(ordered) < 2:
raise ValueError(
"interp_method='monotone_cubic' requires at least two "
"fitted slices"
)
self.events = tuple(sorted(
((e if isinstance(e, VarianceEvent) else VarianceEvent(*e))
for e in events),
key=lambda e: e.time,
))
self.model = model
self.r = float(r)
self.interp_method = interp_method
self.fit_report = fit_report
self._slices = ordered
# ── Construction ─────────────────────────────────────────────────
[docs]
@classmethod
def fit(
cls,
df,
model: Union[str, Parametrization] = "svi",
arbitrage_condition: ArbitrageFreedom = ArbitrageFreedom.QUASI,
r: float = 0.0,
interp_method: str = "total_variance",
mode: str = "warn",
prior: "VolSurface" = None,
anchor: float = 0.0,
events=(),
**model_kwargs,
) -> "VolSurface":
"""Calibrate every maturity slice of an option panel independently.
Expects the ``calibrate_slice`` schema: columns ``strike``,
``iv``, ``maturity``, ``implied_forward``, with multiple
maturities in one DataFrame. Model-specific per-slice arguments
are derived automatically: ``theta`` (ATM total variance) for
SSVI and eSSVI -- plus, for eSSVI only, ``theta_ref``
defaulting to the median across slices -- ``T`` for
jump-wings, and ``T``/``F`` for SABR (``beta``
defaults to 0.5 — override via ``model_kwargs``).
Calibration controls (``objective``, ``loss``, ``f_scale``,
``initialization``) pass through to every slice via
``model_kwargs``. Slices that fail to calibrate are skipped with
a warning; fitting fails only if no slice succeeds. For
cross-slice calendar enforcement use :func:`calibrate_surface`.
Parameters
----------
df : pd.DataFrame
Multi-expiry option panel.
model : str or Parametrization, default "svi"
Factory name or a model instance.
arbitrage_condition : ArbitrageFreedom, default QUASI
Used when ``model`` is a factory name.
r : float, default 0.0
Flat discount rate for the pricing layer.
interp_method : str, default "total_variance"
Maturity interpolation method (see the class docstring).
mode : str, default "warn"
Failure handling for slices. ``"strict"``: the first slice
that is too thin to calibrate or fails to converge raises
with its maturity. ``"warn"``: skipped with a logged
warning. ``"lenient"``: skipped silently. Every mode
records the failure on the fit report.
**model_kwargs
Forwarded to every per-slice calibration.
Returns
-------
VolSurface
"""
validate_mode(mode)
instance = (
get_model(model, arbitrage_condition)
if isinstance(model, str) else model
)
groups = sorted(
((float(T), g) for T, g in df.groupby("maturity")),
key=lambda item: item[0],
)
if not groups:
raise ValueError("VolSurface.fit: empty input panel")
events = tuple(
(e if isinstance(e, VarianceEvent) else VarianceEvent(*e))
for e in events
)
theta_by_T, theta_ref = _data_thetas(instance, groups, events)
slices = []
slice_reports = []
for T, g in groups:
kwargs = _auto_slice_kwargs(
instance, T, g, model_kwargs, theta_by_T, theta_ref
)
n_quotes = len(g)
k, w, F, sel = prepare_slice(g, return_index=True)
if k is not None:
w = _subtract_events(w, events, T, mode)
if k is None or w is None:
if mode == "strict":
raise ValueError(
f"VolSurface.fit(mode='strict'): slice T={T:g} has "
f"insufficient data ({n_quotes} quotes before cleaning)"
)
if mode == "warn":
logger.warning(
f"VolSurface.fit: slice T={T:g} has insufficient data; skipping"
)
slice_reports.append(
build_slice_report(T, SLICE_INSUFFICIENT, n_quotes)
)
continue
_band_kwargs(g, T, k, w, sel, kwargs)
_shift_band_events(kwargs, events, T)
_prior_slice_kwargs(prior, anchor, T, kwargs)
try:
params = instance.calibrate(k, w, **kwargs)
except np.linalg.LinAlgError as exc:
# numerical failure (LinAlgError subclasses ValueError!):
# honor the skip-and-record contract
logger.warning(
f"slice T={T:g} raised LinAlgError: {exc}; skipping"
)
params = None
except ValueError:
# caller bugs (unknown loss/objective, missing kwargs)
# must surface, not dissolve into a skipped slice
raise
except Exception as exc: # noqa: BLE001 -- contract: skip + record
# 'Slices that fail to calibrate are skipped with a
# warning' must hold for exceptions too (e.g.
# DirectSVI's closed form hits a singular matrix on a
# flat junk slice), not only for a None return -- and
# strict mode surfaces it instead of swallowing it.
if mode == "strict":
raise
if mode == "warn":
logger.warning(
f"VolSurface.fit: slice T={T:g} raised "
f"{type(exc).__name__}: {exc}; skipping"
)
params = None
if params is None:
if mode == "strict":
raise ValueError(
f"VolSurface.fit(mode='strict'): slice T={T:g} "
"failed to calibrate"
)
if mode == "warn":
logger.warning(
f"VolSurface.fit: slice T={T:g} failed to calibrate; skipping"
)
slice_reports.append(
build_slice_report(T, SLICE_FAILED, n_quotes, k=k)
)
continue
params["forward"] = F
w_fit = instance.total_variance(k, params)
slice_reports.append(build_slice_report(
T, SLICE_OK, n_quotes, k=k,
iv_mkt=np.sqrt(np.maximum(w, 0.0) / T),
iv_fit=np.sqrt(np.maximum(w_fit, 0.0) / T),
))
slices.append((T, params))
if not slices:
raise ValueError("VolSurface.fit: no slice calibrated successfully")
report = build_surface_report(
slice_reports, type(instance).__name__, model_kwargs,
calendar_enforced=False,
)
return cls(instance, slices, r=r, interp_method=interp_method,
fit_report=report, events=events)
# ── Slice access and maturity location ───────────────────────────
@property
def maturities(self) -> np.ndarray:
"""Fitted maturities, ascending."""
return np.array([T for T, _ in self._slices])
def _locate(self, maturity):
"""Locate a maturity: ("exact", T, params) or ("interp", lo, hi, lam)."""
T = float(maturity)
for Ti, params in self._slices:
if abs(Ti - T) <= 1e-12 * max(1.0, abs(T)):
return ("exact", Ti, params)
Ts = [Ti for Ti, _ in self._slices]
if T < Ts[0] or T > Ts[-1]:
fitted = ", ".join(f"{Ti:g}" for Ti in Ts)
raise ValueError(
f"maturity {T:g} is outside the fitted range [{Ts[0]:g}, {Ts[-1]:g}]"
f" (fitted: {fitted}); extrapolation is not supported"
)
j = bisect_left(Ts, T)
T_lo, p_lo = self._slices[j - 1]
T_hi, p_hi = self._slices[j]
lam = (T - T_lo) / (T_hi - T_lo)
return ("interp", (T_lo, p_lo), (T_hi, p_hi), lam)
def _slice(self, maturity) -> Tuple[float, Dict[str, float]]:
loc = self._locate(maturity)
if loc[0] != "exact":
fitted = ", ".join(f"{Ti:g}" for Ti, _ in self._slices)
raise ValueError(
f"maturity {float(maturity):g} is not a fitted slice"
f" (fitted: {fitted})"
)
return loc[1], loc[2]
def _theta_params(self, lo, hi, lam) -> Dict[str, float]:
"""Synthetic SSVI/eSSVI params at an interpolated maturity."""
(T_lo, p_lo), (T_hi, p_hi) = lo, hi
blend = {
key: (1.0 - lam) * p_lo[key] + lam * p_hi[key]
for key in p_lo
if key != "forward" and key in p_hi
}
blend["forward"] = float(np.exp(
(1.0 - lam) * np.log(p_lo["forward"])
+ lam * np.log(p_hi["forward"])
))
if (
"rho_theta" in blend and isinstance(self.model, ESSVI)
and all(key in blend for key in ("theta_ref", "rho0", "rho1", "alpha"))
):
# keep the stored rho(theta) consistent with the blended theta
blend["rho_theta"] = self.model._rho_of(
blend["theta"], blend["theta_ref"],
blend["rho0"], blend["rho1"], blend["alpha"],
)
# slices carrying only rho_theta (a shape _rho_phi documents as
# sufficient on its own) keep the linearly blended rho_theta
return blend
[docs]
def slice_at(self, maturity) -> Dict[str, float]:
"""Parameter dict of the (possibly synthetic) slice at ``maturity``.
Exact for fitted maturities. Between slices this requires a
parametric interpolation method (``interp_method="theta"``);
the default total-variance blend has no parameter representation
and raises here (evaluation methods still work at any maturity).
"""
loc = self._locate(maturity)
if loc[0] == "exact":
return dict(loc[2])
if self.interp_method != "theta":
raise ValueError(
"slice_at between fitted maturities requires a parametric "
"interpolation (interp_method='theta'); the total-variance "
"and monotone_cubic blends have no parameter "
"representation. Evaluation methods such as "
"iv/total_variance/price work at any maturity in range"
)
return self._theta_params(loc[1], loc[2], loc[3])
[docs]
def params(self, maturity) -> Dict[str, float]:
"""Calibrated parameter dict of the fitted slice at ``maturity`` (a copy)."""
return dict(self._slice(maturity)[1])
[docs]
def forward(self, maturity) -> float:
"""Forward price at ``maturity`` (log-linear between fitted slices)."""
loc = self._locate(maturity)
if loc[0] == "exact":
return float(loc[2]["forward"])
(_, p_lo), (_, p_hi), lam = loc[1], loc[2], loc[3]
return float(np.exp(
(1.0 - lam) * np.log(p_lo["forward"]) + lam * np.log(p_hi["forward"])
))
# ── Evaluation ───────────────────────────────────────────────────
def _pchip(self, k: np.ndarray):
"""Shape-preserving cubic in T at fixed k, over all fitted slices.
Surfaces are immutable after construction, so interpolators are
cached per k-grid (small LRU): pricing a book at repeated
strikes no longer re-evaluates every slice and rebuilds the
PCHIP setup on each call.
"""
from scipy.interpolate import PchipInterpolator
cache = getattr(self, "_pchip_cache", None)
if cache is None:
cache = {}
object.__setattr__(self, "_pchip_cache", cache)
key = (k.shape, k.tobytes())
hit = cache.get(key)
if hit is not None:
return hit
T_knots = np.array([T for T, _ in self._slices])
W = np.vstack([
self.model.total_variance(k, params) for _, params in self._slices
])
interp = PchipInterpolator(T_knots, W, axis=0, extrapolate=False)
if len(cache) >= 16: # bound memory for churning strike grids
cache.pop(next(iter(cache)))
cache[key] = interp
return interp
def _w_at(self, k: np.ndarray, maturity) -> np.ndarray:
"""Total variance at maturity: the interpolated CONTINUOUS
component of the stored slices, plus the cumulative variance of
events realized by ``maturity`` (see :class:`VarianceEvent`)."""
ev = _cum_event_variance(self.events, float(maturity))
return self._w_cont_at(k, maturity) + ev
def _w_cont_at(self, k: np.ndarray, maturity) -> np.ndarray:
loc = self._locate(maturity)
if loc[0] == "exact":
return self.model.total_variance(k, loc[2])
lo, hi, lam = loc[1], loc[2], loc[3]
if self.interp_method == "theta":
return self.model.total_variance(k, self._theta_params(lo, hi, lam))
if self.interp_method == "monotone_cubic":
return self._pchip(k)(float(maturity))
w_lo = self.model.total_variance(k, lo[1])
w_hi = self.model.total_variance(k, hi[1])
return (1.0 - lam) * w_lo + lam * w_hi
def _dw_at(self, k: np.ndarray, maturity, second: bool = False) -> np.ndarray:
deriv = self.model.d2w_dk2 if second else self.model.dw_dk
loc = self._locate(maturity)
if loc[0] == "exact":
return deriv(k, loc[2])
lo, hi, lam = loc[1], loc[2], loc[3]
if self.interp_method == "theta":
return deriv(k, self._theta_params(lo, hi, lam))
if self.interp_method == "monotone_cubic":
# PCHIP slopes are nonlinear in the knot values, so the
# k-derivative of the interpolant is not the interpolant of
# the k-derivatives; central differences on the surface.
h = self.model.fd_step
fn = self._pchip(np.concatenate([k - h, k, k + h]))
row = fn(float(maturity))
n = k.shape[0]
w_m, w_0, w_p = row[:n], row[n:2 * n], row[2 * n:]
if second:
return (w_p - 2.0 * w_0 + w_m) / (h * h)
return (w_p - w_m) / (2.0 * h)
# derivative of the linear blend is the blend of derivatives
return (1.0 - lam) * deriv(k, lo[1]) + lam * deriv(k, hi[1])
@property
def regularity(self) -> str:
"""Smoothness guarantee of the surface in maturity: "C0" or "C1".
"C0" (the default total-variance blend and the theta method):
w(k, T) is continuous in T but its maturity derivative jumps at
every fitted slice. Ready for implied vols, prices, and
sticky-strike Greeks; NOT ready for quantities that consume
dw/dT -- Dupire local volatility, forward variance, PDE
coefficients -- whose inputs would be discontinuous.
"C1" (interp_method="monotone_cubic"): dw/dT exists and is
continuous everywhere in the fitted maturity range (exposed via
:meth:`dw_dT`), making the surface Dupire-ready in maturity.
Smoothness in strike comes from the model itself and is
analytic (C-infinity) for the SVI family either way.
"""
return _REGULARITY[self.interp_method]
[docs]
def dw_dT(self, k, maturity):
"""Maturity derivative of total variance, dw/dT at fixed k.
The Dupire numerator. Only available when the interpolation
method is C1 in maturity (interp_method="monotone_cubic");
the C0 methods have jump discontinuities at the fitted slices,
and a one-sided number there would be silently wrong.
"""
if self.regularity != "C1":
raise ValueError(
f"dw_dT requires a C1 maturity interpolation; this surface "
f"uses interp_method={self.interp_method!r} (regularity "
f"{self.regularity}). Refit or construct with "
"interp_method='monotone_cubic'."
)
T = float(maturity)
self._locate(T) # range check (raises outside the fitted range)
k_arr = np.atleast_1d(np.asarray(k, dtype=np.float64))
values = self._pchip(k_arr).derivative()(T)
return _shape_like(values, k)
[docs]
def total_variance(self, k, maturity):
"""Total variance w(k) at any maturity in the fitted range."""
values = self._w_at(np.atleast_1d(np.asarray(k, dtype=np.float64)), maturity)
return _shape_like(values, k)
[docs]
def iv(self, strike, maturity, return_status: bool = False):
"""Implied volatility at absolute strike(s), any maturity in range.
With ``return_status=True`` also returns a domain-of-validity
label per point -- ``"observed"`` (a fitted maturity, inside
that slice's quoted strike range), ``"interpolated"`` (between
fitted maturities, inside the bracketing slices' joint quoted
range), or ``"extrapolated"`` (outside the quoted domain: the
model's wings, not market information). Requires a fit report
(surfaces constructed directly from parameter dicts carry no
quoted ranges and raise).
"""
T = float(maturity)
F = self.forward(T)
K = np.atleast_1d(np.asarray(strike, dtype=np.float64))
if np.any(K <= 0):
raise ValueError("strikes must be positive")
k = np.log(K / F)
w = self._w_at(k, T)
sigma = np.sqrt(np.maximum(w, 0.0) / T)
if not return_status:
return _shape_like(sigma, strike)
status = self._eval_status(k, T)
if np.ndim(strike) == 0:
return float(sigma[0]), str(status[0])
return sigma, status
def _quoted_range_of(self, T: float):
"""Quoted (k_min, k_max) of the fitted slice at T, from the report."""
for s in self.fit_report.slices:
if np.isclose(s.maturity, T, rtol=1e-9, atol=1e-12) and s.k_min is not None:
return s.k_min, s.k_max
return None
def _eval_status(self, k: np.ndarray, T: float) -> np.ndarray:
if self.fit_report is None:
raise ValueError(
"evaluation status requires a fit report; surfaces "
"constructed directly from parameter dicts carry no "
"quoted ranges"
)
loc = self._locate(T)
if loc[0] == "exact":
rng = self._quoted_range_of(float(T))
in_T = "observed"
else:
(T_lo, _), (T_hi, _) = loc[1], loc[2]
r_lo = self._quoted_range_of(float(T_lo))
r_hi = self._quoted_range_of(float(T_hi))
if r_lo is None or r_hi is None:
rng = r_lo or r_hi
else:
rng = (max(r_lo[0], r_hi[0]), min(r_lo[1], r_hi[1]))
in_T = "interpolated"
status = np.full(k.shape, "extrapolated", dtype="<U12")
if rng is not None:
inside = (k >= rng[0]) & (k <= rng[1])
status[inside] = in_T
return status
[docs]
def atm_vol(self, maturity) -> float:
"""At-the-money (k = 0) implied volatility."""
T = float(maturity)
w0 = float(self._w_at(np.array([0.0]), T)[0])
return float(np.sqrt(max(w0, 0.0) / T))
[docs]
def skew(self, maturity) -> float:
"""ATM total-variance skew dw/dk at k = 0."""
return float(self._dw_at(np.array([0.0]), maturity)[0])
[docs]
def curvature(self, maturity) -> float:
"""ATM total-variance curvature d2w/dk2 at k = 0."""
return float(self._dw_at(np.array([0.0]), maturity, second=True)[0])
# ── Diagnostics ──────────────────────────────────────────────────
[docs]
def check_arbitrage(self, **kwargs) -> ArbitrageReport:
"""Run the arbitrage diagnostics over all fitted slices.
Forwards to :func:`pysvi.diagnostics.check_arbitrage` (butterfly,
Lee wing bounds, calendar); keyword arguments (``k_min``,
``k_max``, ``n_grid``, ``tol``, ``k_data``) pass through.
"""
return check_arbitrage(self.model, self._slices, **kwargs)
[docs]
def diagnose(self, **kwargs) -> SurfaceDiagnostics:
"""Fit report and arbitrage diagnostics in one result block.
Combines :attr:`fit_report` (calibration status, quote
accounting, residuals, settings, provenance) with
:meth:`check_arbitrage`. When no explicit grid is given and a
fit report is available, the arbitrage checks run on the quoted
log-moneyness range (the surface's domain of validity) rather
than the wide default grid. ``print(surface.diagnose())`` renders
the formatted block; all fields are individually accessible.
"""
if (
self.fit_report is not None
and "k_data" not in kwargs
and "k_min" not in kwargs
and "k_max" not in kwargs
):
quoted = self.fit_report.quoted_range()
if quoted is not None:
kwargs["k_data"] = np.array(quoted)
return SurfaceDiagnostics(
fit=self.fit_report,
arbitrage=self.check_arbitrage(**kwargs),
)
# ── Serialization ────────────────────────────────────────────────
[docs]
def save(self, path) -> None:
"""Write the surface to ``path`` as versioned JSON.
The schema captures everything evaluation needs — model name and
arbitrage condition, per-slice maturities and parameters
(forwards included), the flat rate, the interpolation method —
plus the fit report and provenance when present. Calibrating is
expensive and evaluating is cheap: save once, distribute, and
:meth:`load` reproduces evaluation exactly.
"""
import json
from dataclasses import asdict
from .report import _pysvi_version
payload = {
"schema_version": _SCHEMA_VERSION,
"pysvi_version": _pysvi_version(),
"model": type(self.model).__name__,
"arbitrage_condition": int(self.model.arbitrage_condition.value),
"r": self.r,
"interp_method": self.interp_method,
"events": [
{"time": e.time, "variance": e.variance, "label": e.label}
for e in self.events
],
"slices": [
{"maturity": T, "params": params} for T, params in self._slices
],
"fit_report": asdict(self.fit_report) if self.fit_report else None,
}
with open(path, "w", encoding="utf-8") as fh:
json.dump(payload, fh, indent=2)
[docs]
@classmethod
def load(cls, path) -> "VolSurface":
"""Reconstruct a surface saved by :meth:`save`.
Validates the schema version and model name; raises ValueError
on an unknown schema or model.
"""
import json
from .report import SliceFitReport, SurfaceFitReport
with open(path, "r", encoding="utf-8") as fh:
payload = json.load(fh)
version = payload.get("schema_version")
if version != _SCHEMA_VERSION:
raise ValueError(
f"unsupported surface schema_version {version!r} "
f"(this pysvi reads version {_SCHEMA_VERSION})"
)
model_cls = _MODEL_CLASSES.get(payload.get("model"))
if model_cls is None:
raise ValueError(f"unknown model {payload.get('model')!r} in surface file")
model = model_cls(
arbitrage_condition=ArbitrageFreedom(payload.get("arbitrage_condition", 0))
)
report = None
if payload.get("fit_report"):
raw = dict(payload["fit_report"])
raw["slices"] = tuple(SliceFitReport(**item) for item in raw["slices"])
report = SurfaceFitReport(**raw)
return cls(
model,
[(item["maturity"], item["params"]) for item in payload["slices"]],
r=payload.get("r", 0.0),
interp_method=payload.get("interp_method", "total_variance"),
fit_report=report,
events=[
VarianceEvent(e["time"], e["variance"], e.get("label", ""))
for e in payload.get("events", [])
],
)
# ── Black-76 pricing and Greeks ──────────────────────────────────
def _black_inputs(self, strike, maturity):
T = float(maturity)
F = self.forward(T)
K = np.atleast_1d(np.asarray(strike, dtype=np.float64))
if np.any(K <= 0):
raise ValueError("strikes must be positive")
k = np.log(K / F)
w = np.maximum(self._w_at(k, T), 1e-16)
sqrt_w = np.sqrt(w)
d1 = (-k + 0.5 * w) / sqrt_w
d2 = d1 - sqrt_w
return T, F, K, sqrt_w, d1, d2
[docs]
def price(self, strike, maturity, cp: str = "call"):
"""Black-76 option price at absolute strike(s).
::
call = e^{-rT} [F N(d1) - K N(d2)]
put = e^{-rT} [K N(-d2) - F N(-d1)]
with d1 = (log(F/K) + w/2)/sqrt(w), d2 = d1 - sqrt(w), and w the
surface total variance at the strike.
"""
is_call = _is_call(cp)
T, F, K, sqrt_w, d1, d2 = self._black_inputs(strike, maturity)
disc = np.exp(-self.r * T)
if is_call:
value = disc * (F * ndtr(d1) - K * ndtr(d2))
else:
value = disc * (K * ndtr(-d2) - F * ndtr(-d1))
return _shape_like(value, strike)
[docs]
def delta(self, strike, maturity, cp: str = "call"):
"""Black-76 forward delta, e^{-rT} N(d1) (call) or -e^{-rT} N(-d1) (put).
Holds the implied volatility fixed (sticky-strike).
"""
is_call = _is_call(cp)
T, F, K, sqrt_w, d1, _ = self._black_inputs(strike, maturity)
disc = np.exp(-self.r * T)
value = disc * ndtr(d1) if is_call else -disc * ndtr(-d1)
return _shape_like(value, strike)
[docs]
def gamma(self, strike, maturity):
"""Black-76 gamma, e^{-rT} phi(d1) / (F sqrt(w)). Same for calls and puts."""
T, F, K, sqrt_w, d1, _ = self._black_inputs(strike, maturity)
value = np.exp(-self.r * T) * _npdf(d1) / (F * sqrt_w)
return _shape_like(value, strike)
[docs]
def vega(self, strike, maturity):
"""Black-76 vega per unit volatility, e^{-rT} F phi(d1) sqrt(T)."""
T, F, K, sqrt_w, d1, _ = self._black_inputs(strike, maturity)
value = np.exp(-self.r * T) * F * _npdf(d1) * np.sqrt(T)
return _shape_like(value, strike)
[docs]
def theta(self, strike, maturity, cp: str = "call"):
"""Black-76 theta per year of calendar time, holding sigma fixed.
::
theta = r V - e^{-rT} F phi(d1) sigma / (2 sqrt(T))
"""
T, F, K, sqrt_w, d1, _ = self._black_inputs(strike, maturity)
decay = np.exp(-self.r * T) * F * _npdf(d1) * sqrt_w / (2.0 * T)
value = self.r * np.atleast_1d(self.price(strike, maturity, cp)) - decay
return _shape_like(value, strike)
# ── Calendar-aware surface calibration ───────────────────────────────
def _warn_ssvi_admissibility(slices) -> None:
"""Gatheral-Jacquier sufficient no-butterfly bounds for SSVI-form slices.
One aggregated warning when theta*phi*(1+abs(rho)) > 4 or
theta*phi^2*(1+abs(rho)) > 4 (phi = eta/sqrt(theta)) on any slice.
The bounds are sufficient, not necessary — violating them does not
prove arbitrage, it only leaves the slice outside the proven-safe
region, so verify with check_arbitrage().
"""
outside = []
for T, params in slices:
theta = params["theta"]
rho = params.get("rho_theta", params.get("rho"))
phi = params["eta"] / np.sqrt(theta)
lvl = theta * phi * (1.0 + abs(rho))
crv = theta * phi * phi * (1.0 + abs(rho))
if lvl > 4.0 or crv > 4.0:
outside.append(f"T={T:g}")
if outside:
logger.warning(
"SSVI admissibility: slices outside the Gatheral-Jacquier "
f"sufficient no-butterfly bounds: {', '.join(outside)} "
"(the bounds are conservative; verify with check_arbitrage())"
)
[docs]
def calibrate_surface(
df,
model: Union[str, Parametrization] = "ssvi",
enforce_calendar: bool = True,
arbitrage_condition: ArbitrageFreedom = ArbitrageFreedom.QUASI,
r: float = 0.0,
interp_method: str = "total_variance",
mode: str = "warn",
prior: "VolSurface" = None,
anchor: float = 0.0,
events=(),
**model_kwargs,
) -> VolSurface:
"""Calendar-aware multi-expiry calibration returning a VolSurface.
Orders the expiries and calibrates them oldest-first, automatically
evaluating each fitted slice on the next slice's penalty grid and
passing it as ``w_prev`` — the manual chaining the per-slice API
requires. With ``enforce_calendar`` the NO_CALENDAR flag is added to
the model's arbitrage condition, and for SSVI/eSSVI the per-slice
ATM total variances are made non-decreasing before fitting.
For eSSVI the global term structure (rho0, rho1, alpha, eta) is
fitted **jointly across all slices** against the shared theta_ref
(median ATM total variance), rather than independently per slice —
every returned slice carries identical shape parameters.
After fitting, SSVI-form slices are checked against the
Gatheral-Jacquier sufficient no-butterfly bounds and a warning is
logged for slices outside the proven-safe region.
Parameters
----------
df : pd.DataFrame
Multi-expiry option panel (``calibrate_slice`` schema).
model : str or Parametrization, default "ssvi"
Factory name or model instance. DirectSVI is rejected when
``enforce_calendar`` is set (no penalty support).
enforce_calendar : bool, default True
Add NO_CALENDAR to the arbitrage condition and chain w_prev.
arbitrage_condition : ArbitrageFreedom, default QUASI
Base condition when ``model`` is a factory name (NO_BUTTERFLY
may be OR-ed in; NO_CALENDAR is added by ``enforce_calendar``).
r : float, default 0.0
Flat discount rate for the pricing layer.
interp_method : str, default "total_variance"
Maturity interpolation method for the returned surface.
**model_kwargs
Calibration controls and model extras, forwarded to every slice
(``beta`` for SABR, etc.). The bid_ask objective is not
supported in the joint eSSVI path.
Returns
-------
VolSurface
"""
validate_mode(mode)
if isinstance(model, str):
condition = arbitrage_condition
if enforce_calendar:
condition |= ArbitrageFreedom.NO_CALENDAR
instance = get_model(model, condition)
# the converse must hold too: asking for NO_CALENDAR via the
# arbitrage condition means calendar enforcement, with chaining
enforce_calendar = (
enforce_calendar
or ArbitrageFreedom.NO_CALENDAR in instance.arbitrage_condition
)
else:
instance = model
if enforce_calendar and ArbitrageFreedom.NO_CALENDAR not in instance.arbitrage_condition:
instance = type(instance)(
arbitrage_condition=instance.arbitrage_condition
| ArbitrageFreedom.NO_CALENDAR
)
enforce_calendar = (
enforce_calendar
or ArbitrageFreedom.NO_CALENDAR in instance.arbitrage_condition
)
if enforce_calendar and isinstance(instance, DirectSVI):
raise ValueError(
"DirectSVI does not support penalty-based calendar enforcement; "
"use enforce_calendar=False or an iterative model"
)
if interp_method not in _INTERP_METHODS:
raise ValueError(
f"unknown interp_method {interp_method!r}; choose from {_INTERP_METHODS}"
)
if interp_method == "theta" and not isinstance(instance, (SSVI, ESSVI)):
raise ValueError("interp_method='theta' requires an SSVI or eSSVI model")
groups = sorted(
((float(T), g) for T, g in df.groupby("maturity")),
key=lambda item: item[0],
)
if not groups:
raise ValueError("calibrate_surface: empty input panel")
events = tuple(
(e if isinstance(e, VarianceEvent) else VarianceEvent(*e))
for e in events
)
prepared = []
slice_reports = []
for T, g in groups:
k_i, w_i, F_i, sel_i = prepare_slice(g, return_index=True)
if k_i is not None:
w_i = _subtract_events(w_i, events, T, mode)
if k_i is None or w_i is None:
if mode == "strict":
raise ValueError(
f"calibrate_surface(mode='strict'): slice T={T:g} has "
f"insufficient data ({len(g)} quotes before cleaning)"
)
if mode == "warn":
logger.warning(
f"calibrate_surface: slice T={T:g} has insufficient data; skipping"
)
slice_reports.append(
build_slice_report(T, SLICE_INSUFFICIENT, len(g))
)
continue
prepared.append((T, g, k_i, w_i, F_i, sel_i))
if not prepared:
raise ValueError("calibrate_surface: no usable slice in the panel")
theta_by_T, theta_ref = _data_thetas(
instance, (item[:2] for item in prepared), events
)
if enforce_calendar and theta_by_T:
ordered_T = [item[0] for item in prepared]
raw = np.array([theta_by_T[T] for T in ordered_T])
monotone = np.maximum.accumulate(raw)
if np.any(monotone > raw):
if mode == "strict":
bad_T = [f"{T:g}" for T, r_i, m_i in
zip(ordered_T, raw, monotone) if m_i > r_i]
raise ValueError(
"calibrate_surface(mode='strict'): per-slice ATM total "
"variances are not non-decreasing (calendar-arbitrageable "
f"ATM term structure) at T = {', '.join(bad_T)}"
)
if mode == "warn":
logger.warning(
"calibrate_surface: per-slice ATM total variances were not "
"non-decreasing; clipped upward to enforce a monotone theta(T)"
)
theta_by_T = dict(zip(ordered_T, (float(x) for x in monotone)))
if isinstance(instance, ESSVI):
if prior is not None and mode == "warn":
logger.warning(
"calibrate_surface: prior/anchor are not supported in the "
"joint eSSVI fit yet; fitting unanchored"
)
slices = _calibrate_essvi_global(
instance, prepared, theta_by_T, theta_ref, enforce_calendar, model_kwargs
)
else:
slices = []
prev_params = None
for T, g, k_i, w_i, F_i, sel_i in prepared:
kwargs = _auto_slice_kwargs(
instance, T, g, model_kwargs, theta_by_T, theta_ref
)
_band_kwargs(g, T, k_i, w_i, sel_i, kwargs)
_shift_band_events(kwargs, events, T)
_prior_slice_kwargs(prior, anchor, T, kwargs)
if enforce_calendar and prev_params is not None:
grid = _penalty_grid(k_i)
kwargs["w_prev"] = instance.total_variance(grid, prev_params)
try:
params = instance.calibrate(k_i, w_i, **kwargs)
except np.linalg.LinAlgError as exc:
# numerical failure (LinAlgError subclasses ValueError!):
# honor the skip-and-record contract
logger.warning(
f"slice T={T:g} raised LinAlgError: {exc}; skipping"
)
params = None
except ValueError:
# caller bugs (unknown loss/objective, missing band or
# prior kwargs) must surface, not dissolve into
# 'no slice calibrated successfully'
raise
except Exception as exc: # noqa: BLE001 -- contract: skip + record
if mode == "strict":
raise
if mode == "warn":
logger.warning(
f"calibrate_surface: slice T={T:g} raised "
f"{type(exc).__name__}: {exc}; skipping"
)
params = None
if params is None:
if mode == "strict":
raise ValueError(
f"calibrate_surface(mode='strict'): slice T={T:g} "
"failed to calibrate"
)
if mode == "warn":
logger.warning(
f"calibrate_surface: slice T={T:g} failed to calibrate; skipping"
)
continue
params["forward"] = F_i
slices.append((T, params))
prev_params = params
if not slices:
raise ValueError("calibrate_surface: no slice calibrated successfully")
if isinstance(instance, (SSVI, ESSVI)) and mode != "lenient":
_warn_ssvi_admissibility(slices)
params_by_T = dict(slices)
for T, g, k_i, w_i, F_i, _sel in prepared:
if T in params_by_T:
w_fit = instance.total_variance(k_i, params_by_T[T])
slice_reports.append(build_slice_report(
T, SLICE_OK, len(g), k=k_i,
iv_mkt=np.sqrt(np.maximum(w_i, 0.0) / T),
iv_fit=np.sqrt(np.maximum(w_fit, 0.0) / T),
))
else:
slice_reports.append(build_slice_report(T, SLICE_FAILED, len(g), k=k_i))
slice_reports.sort(key=lambda item: item.maturity)
report = build_surface_report(
slice_reports, type(instance).__name__, model_kwargs,
calendar_enforced=enforce_calendar,
)
return VolSurface(instance, slices, r=r, interp_method=interp_method,
fit_report=report, events=events)
def _calibrate_essvi_global(
instance, prepared, theta_by_T, theta_ref, enforce_calendar, model_kwargs
):
"""Joint eSSVI fit: one (rho0, rho1, alpha, eta) across all slices."""
if model_kwargs.get("objective") == "bid_ask":
raise ValueError(
"objective='bid_ask' is not supported in the joint eSSVI fit"
)
check_bf = ArbitrageFreedom.NO_BUTTERFLY in instance.arbitrage_condition
init = _initialization(model_kwargs)
ctx = []
for T, g, k_i, w_i, F_i, _sel in prepared:
mode, loss_code, weights, w_lo, w_hi = _prepare_loss_inputs(
k_i, w_i, model_kwargs
)
ctx.append((T, k_i, w_i, F_i, theta_by_T[T], mode, loss_code, weights))
# one common grid spanning all slices, for butterfly and calendar
# (the same grid policy as every per-slice penalty: _penalty_grid)
common_grid = _penalty_grid(np.concatenate([item[1] for item in ctx]))
empty = np.empty(0)
core = _kernels.resolve("essvi_obj")
theta_ref = float(theta_ref)
def make_objective(kernel, loss_override=None, f_scale_val=1.0,
w_fn=None):
def objective(p):
p = np.asarray(p, dtype=np.float64)
total = 0.0
w_prev = empty
has_prev = False
for T, k_i, w_i, F_i, theta_i, mode, loss_code, weights in ctx:
lc = loss_code if loss_override is None else loss_override
total += kernel(
p, k_i, w_i, theta_i, theta_ref, common_grid, w_prev,
check_bf, enforce_calendar, has_prev,
mode, weights, empty, empty, lc, f_scale_val,
)
if enforce_calendar:
rho_t = ESSVI._rho_of(theta_i, theta_ref, p[0], p[1], p[2])
phi = p[3] / np.sqrt(theta_i)
# the w kernel must match the objective kernel:
# the plain-rank path stays fully plain, or the
# chained w_prev reintroduces fastmath platform
# noise into the candidate ranking
w_eval = w_fn if w_fn is not None else essvi_total_variance
w_prev = w_eval(common_grid, theta_i, rho_t, phi)
has_prev = True
return total
return objective
x0 = np.array([0.0, -0.5, 0.5, 1.0])
bounds = [(-0.999, 0.999), (-2.0, 2.0), (-2.0, 2.0), (1e-8, None)]
validity = lambda x: x[3] > 0
tight = {"ftol": 1e-15, "gtol": 1e-12, "maxiter": 1000}
# Robust-loss scale: same policy as every per-slice path
# (_resolve_f_scale) -- explicit kwarg wins; l2 needs none; else
# 1.4826 * MAD of the joint mode-space residuals at a pilot l2 fit.
f_scale = model_kwargs.get("f_scale")
is_l2 = all(item[6] == 0 for item in ctx)
if f_scale is not None:
f_scale = float(f_scale)
elif is_l2:
f_scale = 1.0
else:
pilot = _minimize_with_starts(
make_objective(core, loss_override=0, f_scale_val=1.0),
[x0], bounds, lbfgs_options=tight, nm_options={"maxiter": 2000},
validity=validity,
)
if pilot is not None:
rho0_p, rho1_p, alpha_p, eta_p = pilot.x
resid_k, resid_wm, resid_wt, resid_wgt = [], [], [], []
for T, k_i, w_i, F_i, theta_i, mode, loss_code, weights in ctx:
rho_t = ESSVI._rho_of(theta_i, theta_ref, rho0_p, rho1_p, alpha_p)
resid_k.append(k_i)
resid_wm.append(essvi_total_variance(
k_i, theta_i, rho_t, eta_p / np.sqrt(theta_i)))
resid_wt.append(w_i)
resid_wgt.append(weights)
all_wgt = np.concatenate(resid_wgt) if resid_wgt[0].size else empty
f_scale = _mad_scale(
np.concatenate(resid_k), np.concatenate(resid_wm),
np.concatenate(resid_wt), ctx[0][5], all_wgt, empty, empty,
)
else:
f_scale = 1.0
objective = make_objective(core, f_scale_val=f_scale)
starts = _multistart_variants(x0, 0, 3) if init == "multi_start" else [x0]
res = _minimize_with_starts(
objective, starts, bounds,
lbfgs_options=tight,
nm_options={"maxiter": 2000},
validity=validity,
rank_objective=make_objective(
_kernels._PLAIN["essvi_obj"], f_scale_val=f_scale,
w_fn=_kernels._PLAIN["essvi_w"],
),
)
if res is None:
raise ValueError("calibrate_surface: joint eSSVI fit failed to converge")
rho0, rho1, alpha, eta = (float(x) for x in res.x)
if eta <= 0:
raise ValueError("calibrate_surface: joint eSSVI fit returned eta <= 0")
slices = []
for T, k_i, w_i, F_i, theta_i, mode, loss_code, weights in ctx:
slices.append((T, {
"rho0": rho0, "rho1": rho1, "alpha": alpha, "eta": eta,
"theta": float(theta_i), "theta_ref": theta_ref,
"rho_theta": ESSVI._rho_of(theta_i, theta_ref, rho0, rho1, alpha),
"forward": float(F_i),
}))
return slices