quant/src/localvol.rs
Dupire's formula, applied numerically to a model's own option prices.
//! Dupire's formula, applied numerically to a model's own option prices.//!//! The local volatility chapter derives//!//! ```text//! sigma_loc^2(T,K) = 2 (dC/dT) / (K^2 d2C/dK2)//! ```//!//! for zero rates, and this computes exactly that by differencing prices. Two//! reasons to do it the slow, obvious way rather than analytically.//!//! The first is that it is a genuine check. The derivation in the local//! volatility chapter is several integrations by parts long, and the claim it//! ends on --- that the local volatility extracted from a smile has twice the//! smile's slope --- is carried forward into the smile dynamics chapter and//! used to explain why local volatility gets the dynamics wrong. A claim doing//! that much work should be verified against something that does not share its//! algebra. The tests below do that.//!//! The second is that this is what is actually done to a real surface: the//! market gives prices at discrete strikes and expiries, and Dupire is a//! numerical operation on them. The instability the chapter warns about is//! visible here, in the need to choose the bump sizes with care. use crate::black::Side;use crate::smile::Model; /// The Dupire local volatility implied by a model's European prices.////// Note what this does *not* do: it never looks at the model's own diffusion/// coefficient. It only ever asks the model for option prices, exactly as one/// would ask the market, and then extracts the local volatility from those. So/// when the answer for a lognormal model comes back as the lognormal volatility,/// that is a statement about Dupire's formula rather than a tautology.////// `dt` and `dk` are the differencing steps. The second derivative in strike is/// the delicate one, and it sits in the denominator.pub fn dupire_local_vol(model: &Model, f: f64, t: f64, k: f64, dt: f64, dk: f64) -> Option<f64> { if t <= dt || k <= dk || f <= 0.0 { return None; } // A one-sided difference in T would be first order and, at the accuracy the // rule-of-two test needs, that is not enough. let c_up = model.price(f, k, t + dt, Side::Call); let c_dn = model.price(f, k, t - dt, Side::Call); let dc_dt = (c_up - c_dn) / (2.0 * dt); let d2c_dk2 = model.density(f, t, k, dk); if !(d2c_dk2 > 0.0) { // A non-convex price in the strike is an arbitrage, not a small number. return None; } let variance = 2.0 * dc_dt / (k * k * d2c_dk2); if variance > 0.0 { Some(variance.sqrt()) } else { None }} /// The slope of a curve at `x`, by central difference in log-moneyness.////// Used to compare the slope of the implied smile with the slope of the local/// volatility curve extracted from it, which is the comparison the rule of two/// is about.#[cfg(test)]fn slope_in_log_moneyness(f: f64, x: f64, h: f64, mut curve: impl FnMut(f64) -> Option<f64>) -> Option<f64> { let up = curve(f * (x + h).exp())?; let dn = curve(f * (x - h).exp())?; Some((up - dn) / (2.0 * h))} #[cfg(test)]mod tests { use super::*; const F: f64 = 100.0; #[test] fn dupire_returns_the_volatility_it_was_given() { // The example worked by hand in the local volatility chapter: a flat // smile must give a constant local volatility equal to it. Here the // prices go in and the volatility comes out, with the formula in // between. let m = Model::Lognormal { sigma: 0.25 }; for k in [70.0, 100.0, 140.0] { for t in [0.5, 1.0, 2.0] { let lv = dupire_local_vol(&m, F, t, k, 1e-3, 0.02).unwrap(); assert!((lv - 0.25).abs() < 2e-4, "k={k} t={t} gave {lv}"); } } } #[test] fn local_volatility_has_twice_the_slope_of_the_implied_smile() { // Equation (7.10): sigma_loc(y) ~ sigma_ATM + 2 s y, derived by hand in // the local volatility chapter by expanding Dupire in implied- // volatility coordinates to first order in the skew. // // Checked here against a numerical Dupire that shares none of that // algebra. A displaced diffusion is the right test model: its skew is // nearly a straight line, which is the assumption the expansion makes. let t = 1.0; let h = 0.02; for beta in [0.9, 0.7, 0.5] { let m = Model::DisplacedDiffusion { sigma: 0.25, beta }; let implied_slope = slope_in_log_moneyness(F, 0.0, h, |k| m.implied_vol(F, k, t)).unwrap(); let local_slope = slope_in_log_moneyness(F, 0.0, h, |k| dupire_local_vol(&m, F, t, k, 1e-3, 0.02)) .unwrap(); let ratio = local_slope / implied_slope; assert!( (ratio - 2.0).abs() < 0.06, "beta={beta}: local slope {local_slope:.5} over implied slope \ {implied_slope:.5} is {ratio:.3}, expected 2" ); } } #[test] fn the_rule_of_two_decays_away_from_the_money() { // Honesty about the approximation's range. It is a first-order result // about the slope *at* a point, and the expansion drops terms of order // y^2; so it should be sharp at the money and degrade steadily as the // strike moves away. A reader applying (7.10) far into a wing should // expect a factor nearer 1.5 than 2. let t = 1.0; let h = 0.02; let m = Model::DisplacedDiffusion { sigma: 0.25, beta: 0.5 }; let ratio_at = |y: f64| { let implied = slope_in_log_moneyness(F, y, h, |k| m.implied_vol(F, k, t)).unwrap(); let local = slope_in_log_moneyness(F, y, h, |k| dupire_local_vol(&m, F, t, k, 1e-3, 0.02)) .unwrap(); local / implied }; // Sharp at the money. assert!((ratio_at(0.0) - 2.0).abs() < 0.02, "{}", ratio_at(0.0)); // And monotonically less so away from it, rather than erratically so. let mut previous = ratio_at(0.0); for y in [0.1, 0.2, 0.3, 0.5] { let r = ratio_at(y); assert!(r < previous, "ratio rose at y={y}: {r} after {previous}"); previous = r; } assert!(previous < 1.7, "expected visible decay by y=0.5, got {previous}"); }} /// A constant elasticity of variance model,////// ```text/// dS = sigma0 * (S/S0)^(beta-1) * S dW,/// ```////// the plainest local volatility model with a skew: `beta < 1` makes volatility/// rise as the level falls, which is the equity shape.////// # What it is here to settle////// The smile dynamics chapter says a local volatility model's forward smile flattens as the/// forward start date recedes. The reason usually given --- that the model has/// one factor, so conditioning on the level leaves only the Brownian motion ---/// does not survive contact with this model. CEV is a one-factor local/// volatility model, and the tests below find that its forward skew does not/// decay at all: at four years out it is *larger* than today's.////// So one-factor-ness is not the mechanism. [`SkewSurface`] holds the one that/// is, and it is about calibration rather than about the number of factors.pub struct Cev { pub sigma0: f64, pub beta: f64, pub spot: f64,} impl Cev { fn local_vol(&self, s: f64) -> f64 { self.sigma0 * (s.max(1e-8) / self.spot).powf(self.beta - 1.0) } /// Simulate to `start`, then on to `start + tenor`, and return the implied /// volatilities of options on the *ratio* `S(start+tenor)/S(start)`. /// /// That ratio is what a forward starting option pays, and its smile is the /// forward smile. Setting `start` to zero recovers the ordinary spot smile, /// which is the comparison that makes the decay visible. pub fn forward_smile( &self, start: f64, tenor: f64, moneyness: &[f64], paths: usize, steps_per_year: usize, seed: u64, ) -> Vec<Option<f64>> { use crate::pathwise::Rng; let mut rng = Rng::new(seed); let total_steps = (((start + tenor) * steps_per_year as f64).ceil() as usize).max(2); let dt = (start + tenor) / total_steps as f64; let split = ((start / (start + tenor)) * total_steps as f64).round() as usize; let mut payoffs = vec![0.0; moneyness.len()]; for _ in 0..paths { let mut s = self.spot; let mut at_start = self.spot; for step in 0..total_steps { if step == split { at_start = s; } // Log-Euler, which keeps the level positive however coarse the // grid, and matters here because CEV with beta < 1 drives the // volatility up exactly where an arithmetic step would go // negative. let v = self.local_vol(s); s *= (-0.5 * v * v * dt + v * dt.sqrt() * rng.next_normal()).exp(); } let ratio = s / at_start; for (i, &k) in moneyness.iter().enumerate() { payoffs[i] += (ratio - k).max(0.0); } } let n = paths as f64; moneyness .iter() .zip(&payoffs) .map(|(&k, &total)| { // The ratio is a martingale under this measure, so its forward // is one and Black-76 applies with that forward. crate::black::implied_vol_black76( total / n, 1.0, k, tenor, crate::black::Side::Call, ) }) .collect() } /// The at-the-money slope of the forward smile in log-moneyness. /// /// The single number the smile dynamics chapter is about: how much skew a /// forward starting option still sees when its start date is pushed out. pub fn forward_skew( &self, start: f64, tenor: f64, paths: usize, steps_per_year: usize, seed: u64, ) -> Option<f64> { let h = 0.05f64; let moneyness = [(-h).exp(), 1.0, h.exp()]; let vols = self.forward_smile(start, tenor, &moneyness, paths, steps_per_year, seed); Some((vols[2]? - vols[0]?) / (2.0 * h)) }} /// A market-shaped implied volatility surface, and the local volatility Dupire/// extracts from it.////// The one feature of a real surface that matters for forward smiles is that/// the at-the-money skew *decays with maturity*, roughly as one over the square/// root of it. So take the simplest surface with that property,////// ```text/// sigma_BS(k, T) = sigma + psi(T) k, psi(T) = psi1 / sqrt(T + T0),/// ```////// with `k = ln(K/F)` the log-moneyness and `T0` a short-maturity floor, since/// no real skew is actually infinite at zero maturity. Nothing depends on the/// linearity in `k`; the smile's curvature plays no part in what follows.////// # The point////// Calibrating a local volatility model to this surface means running Dupire on/// it. The resulting local volatility inherits the maturity decay --- and/// inherits it *as a decay in calendar time*, which is a different thing from a/// decay across maturities. A forward starting option beginning at `t` diffuses/// through the part of the surface from `t` onwards, and that part is flat. The/// forward skew therefore dies as `t` grows, and dies as `1/sqrt(t)`.////// That is the mechanism [`Cev`] lacks. A time-homogeneous local volatility/// model has no calendar decay to inherit, so its forward skew persists --- but/// it also cannot fit the market's maturity term structure. One factor buys you/// today's surface or stationary forward dynamics, not both.pub struct SkewSurface { pub sigma: f64, pub psi1: f64, pub t0: f64, pub spot: f64,} impl SkewSurface { /// The at-the-money skew at maturity `t`, in log-moneyness. pub fn implied_skew(&self, t: f64) -> f64 { self.psi1 / (t + self.t0).sqrt() } /// The implied volatility at log-moneyness `k` and maturity `t`. pub fn implied_vol(&self, k: f64, t: f64) -> f64 { (self.sigma + self.implied_skew(t) * k).max(1e-4) } fn call(&self, k: f64, t: f64) -> f64 { let strike = self.spot * k.exp(); crate::black::black76( self.spot, strike, self.implied_vol(k, t), t, crate::black::Side::Call, ) } /// Dupire's local volatility at maturity `t` and log-moneyness `k`, by /// differencing this surface's own prices. /// /// Zero rates throughout, so the forward is the spot and /// `sigma_loc^2 = 2 (dC/dT) / (K^2 d2C/dK2)`. pub fn local_vol(&self, t: f64, k: f64) -> Option<f64> { // The maturity step has to shrink with the maturity, or the short end, // where the rule of two is sharpest, is differenced across a window as // wide as the maturity itself. let (dt, dk) = ((t / 20.0).min(1e-4), 1e-3); if t <= dt { return None; } let dc_dt = (self.call(k, t + dt) - self.call(k, t - dt)) / (2.0 * dt); // In the strike, not in the log-moneyness, because that is what Dupire's // denominator asks for. let strike = self.spot * k.exp(); let h = strike * dk; let second = |x: f64| { let kk = (x / self.spot).ln(); self.call(kk, t) }; let d2c_dk2 = (second(strike + h) - 2.0 * second(strike) + second(strike - h)) / (h * h); if !(d2c_dk2 > 0.0) { return None; } let variance = 2.0 * dc_dt / (strike * strike * d2c_dk2); (variance > 0.0).then(|| variance.sqrt()) } /// The at-the-money slope of the local volatility curve at calendar time /// `t`, in log-moneyness. pub fn local_skew(&self, t: f64) -> Option<f64> { let h = 0.02; Some((self.local_vol(t, h)? - self.local_vol(t, -h)?) / (2.0 * h)) } /// The closed form the smile dynamics chapter derives for [`Self::local_skew`], /// /// ```text /// psi1 (3T/2 + 2 T0) / (T + T0)^(3/2), /// ``` /// /// obtained by expanding Dupire in total implied variance to first order in /// the skew. Two limits are worth reading off it. At `T = 0` it is twice /// the implied skew, which is the local volatility chapter's rule of two. /// At large `T` it is `1.5 psi1 / sqrt(T)`, so the factor is 1.5 rather /// than 2: a surface whose skew is decaying gives up part of the rule of /// two, because part of Dupire's maturity derivative is now spent on the /// term structure. pub fn predicted_local_skew(&self, t: f64) -> f64 { self.psi1 * (1.5 * t + 2.0 * self.t0) / (t + self.t0).powf(1.5) } /// The skew of a short-dated forward starting option beginning at `t`. /// /// Short maturity, so Berestycki-Busca-Florent applies at the start date and /// the implied skew is half the local one --- the rule of two, read /// backwards. Hence `1/sqrt(t)` decay, against a market forward skew that is /// roughly stationary. pub fn predicted_forward_skew(&self, t: f64) -> f64 { 0.5 * self.predicted_local_skew(t) }} #[cfg(test)]mod forward_smile_tests { use super::*; /// Equity-shaped: 20% at the money, and a one-year skew of about -10 vol /// points per unit of log-moneyness. fn surface() -> SkewSurface { SkewSurface { sigma: 0.20, psi1: -0.10, t0: 1.0 / 12.0, spot: 100.0 } } #[test] fn the_cev_maturity_skew_does_not_decay_either() { // The other half of the chapter's "neither one decays" claim; the // forward-starting half is `a_time_homogeneous_local_vol_model_keeps_its_forward_skew` // below. A skew built from two five-percent moneyness bumps is a small // difference of noisy Monte Carlo quantities, so the ratio is not // reproducible to two figures even at millions of paths -- runs at // 200,000 and at 3,000,000 paths, different seeds, land anywhere from // about 0.9 to about 1.35. What is reproducible, and what the chapter's // argument only needs, is that it stays near 1 rather than decaying // toward the market's ~0.22 or the calibrated local volatility model's // ~0.11. let m = Cev { sigma0: 0.25, beta: 0.7, spot: 100.0 }; let (paths, steps, seed) = (200_000, 100, 20260807); let three_month = m.forward_skew(0.0, 0.25, paths, steps, seed).unwrap(); let five_year = m.forward_skew(0.0, 5.0, paths, steps, seed).unwrap(); let ratio = five_year / three_month; assert!( (0.7..1.4).contains(&ratio), "five-year skew is {ratio:.3} of the three-month skew, expected no strong decay" ); } #[test] fn a_time_homogeneous_local_vol_model_keeps_its_forward_skew() { // The claim the smile dynamics chapter used to make, tested and found false. If // one-factor-ness were the reason forward smiles flatten, this would // decay, and it does not. let m = Cev { sigma0: 0.25, beta: 0.7, spot: 100.0 }; let spot = m.forward_skew(0.0, 1.0, 200_000, 100, 20260807).unwrap(); assert!(spot < 0.0, "beta < 1 should skew down, got {spot}"); for start in [1.0, 2.0, 4.0] { let fwd = m.forward_skew(start, 1.0, 200_000, 100, 20260807).unwrap(); let ratio = fwd / spot; assert!( ratio > 0.9, "forward skew at {start}y is {ratio:.2} of spot, expected no decay" ); } } #[test] fn dupire_on_a_decaying_surface_gives_the_predicted_local_skew() { // The step the chapter's argument rests on, checked against a numerical // Dupire that shares none of its algebra. let s = surface(); for t in [0.25, 0.5, 1.0, 2.0, 5.0, 10.0] { let measured = s.local_skew(t).unwrap(); let predicted = s.predicted_local_skew(t); assert!( (measured - predicted).abs() < 0.05 * predicted.abs(), "t={t}: measured {measured:.4}, predicted {predicted:.4}" ); } } #[test] fn the_rule_of_two_is_a_rule_of_one_and_a_half_when_the_skew_decays() { // At zero maturity the local skew doubles the implied one, as the local // volatility chapter says. At long maturity the ratio has fallen to // 3/2, because Dupire's maturity derivative is now also picking up the // term structure. Both ends read off the same measured surface. let s = surface(); let ratio = |t: f64| s.local_skew(t).unwrap() / s.implied_skew(t); // The approach to 2 is from below and is not fast: the ratio is already // 1.95 a few days out, because it is governed by t against the floor t0 // rather than by t alone. assert!((ratio(0.001) - 2.0).abs() < 0.02, "short end gave {}", ratio(0.001)); assert!((ratio(0.01) - 1.95).abs() < 0.02, "a few days out gave {}", ratio(0.01)); assert!((ratio(20.0) - 1.5).abs() < 0.05, "long end gave {}", ratio(20.0)); // And monotonically between, rather than erratically. let mut previous = ratio(0.001); for t in [0.25, 1.0, 4.0, 20.0] { let r = ratio(t); assert!(r < previous, "ratio rose at t={t}: {r} after {previous}"); previous = r; } } #[test] fn the_forward_skew_of_a_calibrated_surface_decays_as_one_over_root_t() { // The conclusion. A forward starting option at t sees half the local // skew at t, and that has a 1/sqrt(t) tail, so doubling the start date // multiplies the forward skew by 1/sqrt(2). let s = surface(); let today = s.predicted_forward_skew(0.0); // A fifth left after a year, and under a tenth after ten. let after = |t: f64| s.predicted_forward_skew(t) / today; assert!((after(1.0) - 0.21).abs() < 0.02, "one year: {}", after(1.0)); assert!((after(10.0) - 0.07).abs() < 0.02, "ten years: {}", after(10.0)); // And the tail is the square root law, not something steeper or shallower. for t in [4.0, 8.0, 16.0] { let halving = s.predicted_forward_skew(2.0 * t) / s.predicted_forward_skew(t); let root_half = 0.5f64.sqrt(); assert!( (halving - root_half).abs() < 0.02, "t={t}: doubling the start date scaled the skew by {halving:.3}, \ expected {root_half:.3}" ); } }}