quant/src/sabr.rs
SABR, and Hagan's implied volatility expansion.
//! SABR, and Hagan's implied volatility expansion.//!//! The model of the smile dynamics chapter://!//! ```text//! dF = alpha F^beta dW, dalpha = nu alpha dZ, dW dZ = rho dt//! ```//!//! Four parameters, and each does one recognisable thing to the smile: `alpha`//! sets its level, `beta` and `rho` between them set its slope, and `nu` sets//! its curvature. That separation is the reason the model is used — a trader can//! mark a surface by moving parameters whose effects they can predict.//!//! What makes it usable at all is that Hagan and coauthors found an asymptotic//! formula for the implied volatility, so no simulation or PDE is needed to//! price a vanilla. The formula is an expansion in the time to expiry, which is//! also its limitation: it is excellent for short and moderate expiries and it//! degrades for long ones, where it can eventually imply a negative density. /// SABR parameters.#[derive(Clone, Copy, Debug)]pub struct Sabr { /// Initial volatility level. Not an implied volatility — for `beta < 1` it /// carries different units, which is why it is fitted rather than quoted. pub alpha: f64, /// Backbone exponent, in `[0, 1]`. Fixes how volatility scales with the /// level of the forward, and hence how the smile moves when the forward /// moves. The smile dynamics chapter is largely about this parameter. pub beta: f64, /// Correlation between the forward and its volatility, in `(-1, 1)`. pub rho: f64, /// Volatility of volatility. Zero recovers a deterministic-volatility model /// and a smile with no curvature. pub nu: f64,} impl Sabr { /// Hagan's lognormal (Black) implied volatility. /// /// The standard 2002 expansion. Written in the usual two pieces: a /// leading-order term that carries the smile's shape, and a bracket of /// first-order corrections in the expiry. /// /// The at-the-money case is handled separately because the general formula /// has a removable singularity there — `z / x(z)` is `0/0` when the strike /// equals the forward. Taking the limit rather than nudging the strike is /// what keeps the at-the-money volatility exact, and the at-the-money point /// is the one the whole surface is anchored on. pub fn implied_vol(&self, f: f64, k: f64, t: f64) -> Option<f64> { if !(f > 0.0) || !(k > 0.0) || !(t > 0.0) { return None; } let Sabr { alpha, beta, rho, nu } = *self; if !(alpha > 0.0) || !(0.0..=1.0).contains(&beta) || rho.abs() >= 1.0 || nu < 0.0 { return None; } let one_m_b = 1.0 - beta; let fk = (f * k).powf(one_m_b); let log_fk = (f / k).ln(); // The correction bracket, common to both branches. let term1 = one_m_b * one_m_b / 24.0 * alpha * alpha / fk; let term2 = 0.25 * rho * beta * nu * alpha / fk.sqrt(); let term3 = (2.0 - 3.0 * rho * rho) / 24.0 * nu * nu; let correction = 1.0 + (term1 + term2 + term3) * t; // At the money, in the limit k -> f. if log_fk.abs() < 1e-9 { return Some(alpha / f.powf(one_m_b) * correction); } let z = nu / alpha * fk.sqrt() * log_fk; // x(z), the integral that converts the log-moneyness into the // volatility's own coordinate. let x = (((1.0 - 2.0 * rho * z + z * z).sqrt() + z - rho) / (1.0 - rho)).ln(); if !x.is_finite() || x.abs() < 1e-12 { return Some(alpha / f.powf(one_m_b) * correction); } let denom = fk.sqrt() * (1.0 + one_m_b * one_m_b / 24.0 * log_fk * log_fk + one_m_b.powi(4) / 1920.0 * log_fk.powi(4)); let vol = alpha / denom * (z / x) * correction; if vol.is_finite() && vol > 0.0 { Some(vol) } else { None } } /// The at-the-money volatility, which is the formula worth carrying in the /// head: everything else is a correction to it. pub fn atm_vol(&self, f: f64, t: f64) -> Option<f64> { self.implied_vol(f, f, t) } /// The smile over a strike grid. pub fn smile(&self, f: f64, t: f64, strikes: &[f64]) -> Vec<Option<f64>> { strikes.iter().map(|&k| self.implied_vol(f, k, t)).collect() } /// The at-the-money volatility's sensitivity to the forward: the backbone. /// /// This is the quantity the smile dynamics chapter is about. It is not the /// smile's slope in the strike — that is the skew, measured at one instant /// with the forward held still. The backbone is what happens to the at-the- /// money volatility when the market actually moves, and it is what /// determines whether the model's delta is right. /// /// Computed by moving the forward and asking the model again, rather than by /// differentiating the formula, so that it measures the model's behaviour /// and not our algebra. pub fn backbone(&self, f: f64, t: f64, h: f64) -> Option<f64> { let up = self.atm_vol(f * (1.0 + h), t)?; let dn = self.atm_vol(f * (1.0 - h), t)?; Some((up - dn) / (2.0 * h * f)) } /// The smile's slope in the strike, at the money: the skew. /// /// Measured with the forward held fixed, which is the whole distinction from /// [`Sabr::backbone`]. pub fn skew(&self, f: f64, t: f64, h: f64) -> Option<f64> { let up = self.implied_vol(f, f * (1.0 + h), t)?; let dn = self.implied_vol(f, f * (1.0 - h), t)?; Some((up - dn) / (2.0 * h * f)) }} /// Fit `alpha` and `rho` so that the model reproduces a given at-the-money/// volatility and at-the-money skew, with `beta` and `nu` held where they were/// put.////// This is the operation the smile dynamics chapter turns into its main point./// `beta` is not determined by the smile: for any `beta` one likes, there is a/// `rho` that reproduces today's skew, and the resulting models are near enough/// indistinguishable from today's quotes. They then disagree completely about/// what happens when the forward moves, because `beta` is the backbone and/// `rho` is not. So calibrating perfectly to the market leaves the dynamics/// undetermined, and the choice has to be made on other grounds.////// Alternates two one-dimensional solves rather than searching in two/// dimensions at once. It converges quickly because the coupling is weak in the/// right direction: `alpha` mostly sets the level and `rho` mostly sets the/// slope.pub fn calibrate_alpha_rho( target_atm: f64, target_skew: f64, beta: f64, nu: f64, f: f64, t: f64,) -> Option<Sabr> { let h = 0.01; let mut s = Sabr { alpha: target_atm * f.powf(1.0 - beta), beta, rho: 0.0, nu }; for _ in 0..60 { // alpha for the level: the at-the-money volatility is increasing in it. let (mut lo, mut hi) = (1e-6, 10.0 * target_atm.max(1.0) * f.powf(1.0 - beta)); for _ in 0..100 { let mid = 0.5 * (lo + hi); let trial = Sabr { alpha: mid, ..s }; match trial.atm_vol(f, t) { Some(v) if v < target_atm => lo = mid, Some(_) => hi = mid, None => lo = mid, } } s.alpha = 0.5 * (lo + hi); // rho for the slope: the skew is increasing in it. let (mut lo, mut hi) = (-0.999, 0.999); for _ in 0..100 { let mid = 0.5 * (lo + hi); let trial = Sabr { rho: mid, ..s }; match trial.skew(f, t, h) { Some(v) if v < target_skew => lo = mid, Some(_) => hi = mid, None => break, } } s.rho = 0.5 * (lo + hi); } let ok = s.atm_vol(f, t).is_some_and(|v| (v - target_atm).abs() < 1e-6) && s.skew(f, t, h).is_some_and(|v| (v - target_skew).abs() < 1e-7); ok.then_some(s)} #[cfg(test)]mod tests { use super::*; #[test] fn the_smile_dynamics_backbone_table_holds_at_every_beta() { // The "Smile Does Not Determine the Dynamics" table: fix a market (25% // at-the-money, a skew of six volatility points per unit of log-moneyness // -- which at F=100 is -0.0006 per unit of strike, since `Sabr::skew` is // measured per strike -- and a vol-of-vol of 0.35, then fit beta at four // values and read off rho, the backbone, the backbone-to-skew ratio, and // the at-the-money volatility after a ten percent rise in the forward. // Every column of the table, not only the beta=1 and beta=0 rows the // interactive figure computes. let (atm, skew, nu, f, t) = (0.25, -0.0006, 0.35, 100.0, 1.0); let expected = [ // beta, rho, backbone, b/s, atm after +10% (1.0, -0.343, 0.0, 0.0, 0.2500), (0.7, -0.128, -0.000749, 1.248, 0.2430), (0.5, 0.015, -0.001252, 2.086, 0.2384), (0.0, 0.371, -0.002513, 4.188, 0.2272), ]; for (beta, rho, b, b_over_s, atm_after) in expected { let fitted = calibrate_alpha_rho(atm, skew, beta, nu, f, t) .unwrap_or_else(|| panic!("beta={beta}: fit failed")); let backbone = fitted.backbone(f, t, 0.01).unwrap(); let s = fitted.skew(f, t, 0.01).unwrap(); let bumped = fitted.atm_vol(f * 1.10, t).unwrap(); assert!((fitted.rho - rho).abs() < 0.001, "beta={beta}: rho was {}", fitted.rho); assert!((backbone - b).abs() < 5e-6, "beta={beta}: backbone was {backbone}"); if beta < 1.0 { assert!( (backbone / s - b_over_s).abs() < 0.001, "beta={beta}: b/s was {}", backbone / s ); } assert!( (bumped - atm_after).abs() < 5e-5, "beta={beta}: atm after the bump was {bumped}, not {atm_after}" ); } } #[test] fn the_backbone_splits_into_the_skew_and_the_fixed_strike_move() { // The smile dynamics chapter's identity: b = s + d(sigma)/d ln F at // fixed strike. It is only the chain rule, but it is the thing that // keeps the backbone and the fixed-strike move from being confused with // each other -- they are both "how volatility responds when the market // moves" and they differ, under a local volatility model, by a factor // of two and a sign. // // Checked here because `skew` and `backbone` are computed separately and // nothing else forces them to stay consistent. let m = Sabr { alpha: 0.05, beta: 0.5, rho: -0.30, nu: 0.40 }; let (f, t, h) = (0.03f64, 5.0, 1e-4); let skew = f * m.skew(f, t, h).expect("a skew"); let backbone = f * m.backbone(f, t, h).expect("a backbone"); // The third derivative: the forward moves, the strike does not. let up = m.implied_vol(f * (1.0 + h), f, t).expect("a vol"); let down = m.implied_vol(f * (1.0 - h), f, t).expect("a vol"); let fixed_strike = (up - down) / (2.0 * h); assert!( (backbone - skew - fixed_strike).abs() < 1e-7, "b={backbone} against s={skew} plus {fixed_strike}" ); } const F: f64 = 100.0; const T: f64 = 1.0; fn base() -> Sabr { Sabr { alpha: 0.25 * F.powf(0.4), beta: 0.6, rho: -0.3, nu: 0.4 } } #[test] fn zero_vol_of_vol_with_beta_one_is_lognormal() { // With no randomness in the volatility and beta = 1, SABR *is* Black-76, // so the smile must be flat and equal to alpha. let s = Sabr { alpha: 0.3, beta: 1.0, rho: 0.0, nu: 0.0 }; for k in [60.0, 80.0, 100.0, 130.0, 170.0] { let v = s.implied_vol(F, k, T).unwrap(); assert!((v - 0.3).abs() < 1e-9, "k={k} gave {v}"); } } #[test] fn the_at_the_money_limit_is_continuous() { // The general formula divides by x(z), which vanishes at the money. The // separate branch must agree with the limit of the general one, or the // smile has a notch exactly where the market is most liquid. let s = base(); let at = s.atm_vol(F, T).unwrap(); for eps in [1e-3, 1e-4, 1e-5, 1e-6] { let near = s.implied_vol(F, F * (1.0 + eps), T).unwrap(); assert!((near - at).abs() < 1e-3 * (1.0 + eps / 1e-3), "eps={eps}: {near} vs {at}"); } } #[test] fn vol_of_vol_adds_curvature_and_correlation_adds_slope() { // The parameter separation the model is used for. let curvature = |s: &Sabr| { let v: Vec<f64> = [80.0, 100.0, 125.0] .iter() .map(|&k| s.implied_vol(F, k, T).unwrap()) .collect(); v[0] - 2.0 * v[1] + v[2] }; let slope = |s: &Sabr| s.skew(F, T, 0.01).unwrap(); let flat = Sabr { alpha: 0.25, beta: 1.0, rho: 0.0, nu: 0.0 }; let bendy = Sabr { nu: 0.6, ..flat }; let tilted = Sabr { rho: -0.6, nu: 0.3, ..flat }; assert!(curvature(&bendy) > curvature(&flat) + 1e-3, "nu should bend the smile"); assert!(slope(&tilted) < slope(&flat) - 1e-4, "negative rho should tilt it down"); } #[test] fn beta_one_has_no_backbone_and_beta_zero_has_a_strong_one() { // The claim the smile dynamics chapter turns into its central argument. // With beta = 1 the volatility is proportional to the forward, so the // *relative* volatility does not move when the forward does: the at- // the-money implied volatility is nearly flat in F. With beta = 0 the // volatility is absolute, so a rise in the forward divides it by a // bigger number and the at-the-money implied volatility falls steeply. let lognormal = Sabr { alpha: 0.25, beta: 1.0, rho: 0.0, nu: 0.3 }; let normal = Sabr { alpha: 0.25 * F, beta: 0.0, rho: 0.0, nu: 0.3 }; let b_lognormal = lognormal.backbone(F, T, 0.01).unwrap(); let b_normal = normal.backbone(F, T, 0.01).unwrap(); assert!(b_lognormal.abs() < 1e-4, "beta=1 backbone was {b_lognormal}"); assert!(b_normal < -1e-3, "beta=0 backbone was {b_normal}, expected clearly negative"); } #[test] fn the_smile_does_not_determine_the_dynamics() { // The smile dynamics chapter's central claim, as an executable statement. // // Fit the same at-the-money volatility and the same at-the-money skew // with three different backbone exponents. All three reproduce today's // quotes. Their backbones differ by an order of magnitude, so they // disagree about tomorrow's at-the-money volatility, and therefore about // the delta to hedge with today. let (atm, skew, nu) = (0.25, -0.0006, 0.35); let fits: Vec<Sabr> = [1.0, 0.5, 0.0] .iter() .map(|&beta| { calibrate_alpha_rho(atm, skew, beta, nu, F, T) .unwrap_or_else(|| panic!("no fit at beta={beta}")) }) .collect(); // Same smile: they agree away from the money too, not only where they // were fitted. for k in [85.0, 95.0, 105.0, 118.0] { let v: Vec<f64> = fits.iter().map(|s| s.implied_vol(F, k, T).unwrap()).collect(); let spread = v.iter().cloned().fold(f64::MIN, f64::max) - v.iter().cloned().fold(f64::MAX, f64::min); assert!(spread < 0.004, "at k={k} the fits differ by {spread}, too much to call equal"); } // Different dynamics. let backbones: Vec<f64> = fits.iter().map(|s| s.backbone(F, T, 0.01).unwrap()).collect(); assert!(backbones[0].abs() < 1e-4, "beta=1 should be nearly flat, got {}", backbones[0]); assert!( backbones[2] < backbones[0] - 1e-3, "beta=0 backbone {} should be far below beta=1 backbone {}", backbones[2], backbones[0] ); // And the middle one sits between them, as the parameter suggests. assert!(backbones[2] < backbones[1] && backbones[1] < backbones[0]); } #[test] fn backbone_and_skew_are_different_numbers() { // The distinction the chapter exists to make. Two models can be given // the same skew today and still disagree about what happens tomorrow, // because the skew is a snapshot and the backbone is a dynamic. let a = Sabr { alpha: 0.25, beta: 1.0, rho: -0.5, nu: 0.35 }; let b = Sabr { alpha: 0.25 * F.powf(0.5), beta: 0.5, rho: -0.24, nu: 0.35 }; let (sa, sb) = (a.skew(F, T, 0.01).unwrap(), b.skew(F, T, 0.01).unwrap()); let (ba, bb) = (a.backbone(F, T, 0.01).unwrap(), b.backbone(F, T, 0.01).unwrap()); // Close on today's skew... assert!((sa - sb).abs() < 0.2 * sa.abs().max(sb.abs()), "skews {sa} vs {sb}"); // ...and far apart on how the smile will move. assert!((ba - bb).abs() > 1e-4, "backbones {ba} vs {bb} are too close to make the point"); }}