Coverage for src / qsmile / models / sabr.py: 94%
71 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-01 22:47 +0000
1"""SABR stochastic volatility model — Hagan et al. (2002) lognormal approximation."""
3from __future__ import annotations
5from dataclasses import dataclass
6from typing import ClassVar
8import numpy as np
9from numpy.typing import ArrayLike, NDArray
11from qsmile.core.coords import XCoord, YCoord
12from qsmile.models.base import SmileModel
15@dataclass
16class SABRModel(SmileModel):
17 """SABR model with Hagan (2002) lognormal implied volatility approximation.
19 The SABR model describes the dynamics of a forward rate F and its
20 stochastic volatility alpha via:
22 dF = alpha * F^beta * dW_1
23 dalpha = nu * alpha * dW_2
24 <dW_1, dW_2> = rho * dt
26 The Hagan et al. (2002) formula maps these parameters to a
27 closed-form lognormal implied volatility approximation.
29 Fitted parameters (included in the parameter vector):
30 alpha, beta, rho, nu
32 Context (provided via ``metadata``):
33 expiry and forward are read from ``metadata.texpiry``
34 and ``metadata.forward``.
36 Parameters
37 ----------
38 alpha : float
39 Initial volatility. Must be > 0.
40 beta : float
41 CEV exponent. Must be in [0, 1].
42 rho : float
43 Correlation between forward and vol. Must be in (-1, 1).
44 nu : float
45 Vol-of-vol. Must be >= 0.
46 metadata : SmileMetadata
47 Market context containing expiry, forward, etc.
48 """
50 alpha: float
51 beta: float
52 rho: float
53 nu: float
55 # -- Class-level model metadata --
57 native_x_coord: ClassVar[XCoord] = XCoord.LogMoneynessStrike
58 native_y_coord: ClassVar[YCoord] = YCoord.Volatility
59 param_names: ClassVar[tuple[str, ...]] = ("alpha", "beta", "rho", "nu")
60 bounds: ClassVar[tuple[list[float], list[float]]] = (
61 [1e-8, 0.0, -0.999, 0.0],
62 [np.inf, 1.0, 0.999, np.inf],
63 )
65 def __post_init__(self) -> None:
66 """Validate SABR parameter constraints."""
67 super().__post_init__()
68 if self.alpha <= 0:
69 msg = f"alpha must be positive, got {self.alpha}"
70 raise ValueError(msg)
71 if not (0 <= self.beta <= 1):
72 msg = f"beta must be in [0, 1], got {self.beta}"
73 raise ValueError(msg)
74 if not (-1 < self.rho < 1):
75 msg = f"rho must be in (-1, 1), got {self.rho}"
76 raise ValueError(msg)
77 if self.nu < 0:
78 msg = f"nu must be non-negative, got {self.nu}"
79 raise ValueError(msg)
81 def _evaluate(self, x: ArrayLike) -> NDArray[np.float64] | np.float64:
82 """Compute Hagan (2002) lognormal implied volatility at log-moneyness values.
84 Parameters
85 ----------
86 x : ArrayLike
87 Log-moneyness k = ln(K/F).
89 Returns:
90 -------
91 NDArray[np.float64] | np.float64
92 Implied volatility (lognormal).
93 """
94 k = np.asarray(x, dtype=np.float64)
95 forward = self.metadata.forward
96 if forward is None:
97 msg = "forward must be set in metadata before evaluating SABR"
98 raise ValueError(msg)
99 expiry = self.metadata.texpiry
100 strikes = forward * np.exp(k)
101 return self._hagan_implied_vol(forward, strikes, expiry, self.alpha, self.beta, self.rho, self.nu)
103 @staticmethod
104 def _hagan_implied_vol(
105 fwd: float,
106 strikes: NDArray[np.float64] | float,
107 expiry: float,
108 alpha: float,
109 beta: float,
110 rho: float,
111 nu: float,
112 ) -> NDArray[np.float64] | np.float64:
113 """Hagan et al. (2002) lognormal implied vol approximation.
115 Handles ATM (K ≈ F) and OTM/ITM cases separately for numerical
116 stability.
117 """
118 strikes = np.asarray(strikes, dtype=np.float64)
119 eps = 1e-12
121 # ATM mask
122 is_atm = np.abs(strikes - fwd) < eps * fwd
124 # --- ATM formula ---
125 fb = fwd ** (1 - beta)
126 atm_vol = (alpha / fb) * (
127 1
128 + (
129 ((1 - beta) ** 2 / 24) * alpha**2 / fwd ** (2 * (1 - beta))
130 + 0.25 * rho * beta * nu * alpha / fb
131 + (2 - 3 * rho**2) / 24 * nu**2
132 )
133 * expiry
134 )
136 # --- OTM/ITM formula ---
137 fk_mid = np.where(is_atm, fwd, np.sqrt(fwd * strikes))
138 fk_mid_b = fk_mid ** (1 - beta)
139 log_fk = np.where(is_atm, 0.0, np.log(fwd / strikes))
141 z = np.where(is_atm, 0.0, (nu / alpha) * fk_mid_b * log_fk)
142 x_z = np.where(
143 is_atm,
144 1.0,
145 np.where(
146 np.abs(z) < eps,
147 1.0,
148 z / np.log((np.sqrt(1 - 2 * rho * z + z**2) + z - rho) / (1 - rho + eps)),
149 ),
150 )
152 term1 = alpha / (fk_mid_b * (1 + (1 - beta) ** 2 / 24 * log_fk**2 + (1 - beta) ** 4 / 1920 * log_fk**4))
153 correction = (
154 1
155 + (
156 (1 - beta) ** 2 / 24 * alpha**2 / fk_mid ** (2 * (1 - beta))
157 + 0.25 * rho * beta * nu * alpha / fk_mid_b
158 + (2 - 3 * rho**2) / 24 * nu**2
159 )
160 * expiry
161 )
163 otm_vol = term1 * x_z * correction
165 result = np.where(is_atm, atm_vol, otm_vol)
166 # Clamp to avoid negative implied vol from numerical issues
167 return np.maximum(result, eps)
169 @staticmethod
170 def initial_guess(x: NDArray[np.float64], y: NDArray[np.float64]) -> NDArray[np.float64]:
171 """Compute a heuristic initial guess for SABR parameters from market data.
173 Parameters
174 ----------
175 x : NDArray[np.float64]
176 Log-moneyness values.
177 y : NDArray[np.float64]
178 Observed implied volatility values.
179 """
180 # alpha: ATM implied vol is a reasonable starting point
181 atm_idx = int(np.argmin(np.abs(x)))
182 alpha0 = max(float(y[atm_idx]), 0.01)
184 # beta: start at 0.5 (between normal and lognormal)
185 beta0 = 0.5
187 # rho: estimate skew direction from slope of iv vs moneyness
188 if len(x) >= 3:
189 coeffs = np.polyfit(x, y, 1)
190 rho0 = float(np.clip(coeffs[0] / (alpha0 + 1e-8), -0.9, 0.9))
191 else:
192 rho0 = 0.0
194 # nu: estimate from curvature (smile convexity)
195 if len(x) >= 3:
196 coeffs2 = np.polyfit(x, y, 2)
197 nu0 = max(float(np.abs(coeffs2[0])) * 2, 0.1)
198 else:
199 nu0 = 0.3
201 return np.array([alpha0, beta0, rho0, nu0])