quant/src/lsv.rs
The conditional expectation at the heart of the local-stochastic volatility chapter's leverage function.
//! The conditional expectation at the heart of the local-stochastic volatility//! chapter's leverage function.//!//! The calibration condition is//!//! ```text//! L(t,K) = sigma_loc(t,K) / sqrt( E[ v_t | S_t = K ] ),//! ```//!//! and the word doing the work is *conditional*. This computes that conditional//! expectation in the one case where it can be done in closed form — the two-//! regime model the local-stochastic volatility chapter works through by hand —//! so that the shape can be drawn rather than described.//!//! The model: a coin is tossed at time zero and the forward is lognormal for the//! rest of its life, with one of two volatilities. Crude, deliberately. It is the//! smallest thing that is genuinely stochastic volatility, and everything the//! full calibration has to cope with is already visible in it. use crate::smile::lognormal_density; /// A volatility drawn once, from two possibilities.#[derive(Clone, Copy, Debug)]pub struct TwoRegime { /// Probability of the low-volatility state. pub weight: f64, pub vol_low: f64, pub vol_high: f64,} impl TwoRegime { /// The unconditional mean variance — what one would divide by if the word /// "conditional" in the calibration condition were ignored. pub fn mean_variance(&self) -> f64 { self.weight * self.vol_low.powi(2) + (1.0 - self.weight) * self.vol_high.powi(2) } /// `E[v | F_T = k]`, by Bayes' rule. /// /// Observing where the forward ended is evidence about which state was /// drawn, and the posterior weights the two variances. A move far from the /// start is strong evidence for the noisy state, which is why the answer /// rises in both wings. pub fn conditional_variance(&self, f: f64, t: f64, k: f64) -> f64 { let lo = self.weight * lognormal_density(f, t, self.vol_low, k); let hi = (1.0 - self.weight) * lognormal_density(f, t, self.vol_high, k); let total = lo + hi; if total <= 0.0 { // Beyond the reach of both densities; the limit is the noisier one, // since it dominates arbitrarily far out. return self.vol_high.powi(2); } (lo * self.vol_low.powi(2) + hi * self.vol_high.powi(2)) / total } /// The posterior probability of the high-volatility state given where the /// forward ended. The quantity the conditional variance is an average over. pub fn posterior_high(&self, f: f64, t: f64, k: f64) -> f64 { let lo = self.weight * lognormal_density(f, t, self.vol_low, k); let hi = (1.0 - self.weight) * lognormal_density(f, t, self.vol_high, k); let total = lo + hi; if total <= 0.0 { 1.0 } else { hi / total } } /// The leverage function this regime model would need, against a target /// local volatility surface. /// /// With a flat target the shape is entirely the conditioning: `L` dips where /// the stochastic part is already supplying enough volatility and rises /// where it is not. That is the general shape, and the reason a well-chosen /// stochastic part leaves `L` close to one. pub fn leverage(&self, f: f64, t: f64, k: f64, target_local_vol: f64) -> f64 { target_local_vol / self.conditional_variance(f, t, k).sqrt() } /// The same model with its randomness scaled by the mixing weight `lambda`. /// /// The chapter's dial, in the one model where it can be turned by hand. Here /// the vol-of-vol is the gap between the two states, so `lambda` closes that /// gap towards the common mean while leaving the mean variance exactly where /// it was. At `lambda = 0` the volatility is deterministic and the model is /// local volatility; at `lambda = 1` it is the original. /// /// Holding the mean variance fixed is what makes the comparison mean /// anything: it isolates the randomness of the variance from its level, so /// that what moves across the dial is the split between `L` and `v` rather /// than the total amount of volatility in the model. pub fn with_mixing(&self, lambda: f64) -> TwoRegime { let mean = self.mean_variance(); let scale = |v: f64| { let shifted = mean + lambda * lambda * (v * v - mean); shifted.max(0.0).sqrt() }; TwoRegime { weight: self.weight, vol_low: scale(self.vol_low), vol_high: scale(self.vol_high), } }} #[cfg(test)]mod tests { use super::*; const F: f64 = 100.0; const T: f64 = 1.0; fn regime() -> TwoRegime { TwoRegime { weight: 0.5, vol_low: 0.15, vol_high: 0.35 } } #[test] fn the_conditional_variance_is_humped() { // The local-stochastic volatility chapter's table: lowest at the money, // rising in both wings, and the shape is what makes the leverage // function non-trivial. let r = regime(); let at = r.conditional_variance(F, T, F); for k in [55.0, 70.0, 140.0, 180.0] { let away = r.conditional_variance(F, T, k); assert!(away > at, "at k={k} the conditional variance {away} did not exceed {at}"); } } #[test] fn it_is_bracketed_by_the_two_regimes() { // It is an average of two numbers, so it can never leave their range — // a weak check, and the one that would catch a sign or a normalisation // error immediately. let r = regime(); for i in 0..200 { let k = 30.0 + i as f64 * 1.5; let v = r.conditional_variance(F, T, k).sqrt(); assert!( v >= r.vol_low - 1e-9 && v <= r.vol_high + 1e-9, "at k={k} the conditional vol was {v}" ); } } #[test] fn conditioning_matters_by_several_volatility_points() { // The reason the local-stochastic volatility chapter dwells on the // word. If the unconditional average were good enough, the calibration // would be a division and the particle method would not exist. let r = regime(); let unconditional = r.mean_variance().sqrt(); let at_the_money = r.conditional_variance(F, T, F).sqrt(); assert!( unconditional - at_the_money > 0.03, "unconditional {unconditional} against conditional {at_the_money}" ); } #[test] fn the_leverage_function_dips_where_the_stochastic_part_is_loudest() { // With a flat target, L is the inverse of the hump: below one where the // conditional volatility exceeds the target, above it where it falls // short. let r = regime(); let target = r.mean_variance().sqrt(); let at = r.leverage(F, T, F, target); let wing = r.leverage(F, T, 60.0, target); assert!(at > 1.0, "at the money L should exceed one, got {at}"); assert!(wing < at, "L should fall towards the wings, {wing} against {at}"); } #[test] fn the_posterior_moves_the_way_evidence_should() { let r = regime(); // A big move is evidence of the noisy state; a small one is evidence // against it. assert!(r.posterior_high(F, T, F) < 0.4, "no move should favour the quiet state"); assert!(r.posterior_high(F, T, 55.0) > 0.9, "a large fall should favour the noisy one"); assert!(r.posterior_high(F, T, 190.0) > 0.9, "as should a large rise"); }} /// Gyongi's theorem, tested by construction.////// The local volatility chapter's central result is that a process with random/// volatility has the same one-dimensional marginals --- and so the same European/// prices --- as the local volatility model whose squared volatility is the/// conditional expectation of the true one,////// ```text/// sigma_loc^2(t, x) = E[ sigma_t^2 | X_t = x ]./// ```////// The whole of the local-stochastic volatility chapter rests on it, and it had/// been stated and used and never checked. Here it is checked on the one case/// where the conditional expectation is available in closed form: a forward whose/// volatility is drawn once, at time zero, from two values. That is a genuine/// stochastic volatility --- the volatility is random and unknown --- and/// [`TwoRegime::conditional_variance`] is exactly the right-hand side above.////// Simulates the mimicking local volatility model and returns the implied/// volatilities of calls on it, to be compared against the mixture's own.pub fn mimicking_smile( regime: &TwoRegime, forward: f64, expiry: f64, strikes: &[f64], paths: usize, steps: usize, seed: u64,) -> Vec<Option<f64>> { use crate::black::{implied_vol_black76, Side}; use crate::pathwise::Rng; let dt = expiry / steps as f64; let root_dt = dt.sqrt(); let mut rng = Rng::new(seed); let mut payoffs = vec![0.0; strikes.len()]; for _ in 0..paths { let mut f = forward; for i in 0..steps { let t = i as f64 * dt; // The local volatility at this time and level, which is the // conditional expectation Gyongi's theorem prescribes. At t = 0 the // conditioning is vacuous, so fall back to the unconditional mean. let variance = if t <= 0.0 { regime.mean_variance() } else { regime.conditional_variance(forward, t, f) }; let vol = variance.max(0.0).sqrt(); f *= (-0.5 * vol * vol * dt + vol * root_dt * rng.next_normal()).exp(); } for (j, &k) in strikes.iter().enumerate() { payoffs[j] += (f - k).max(0.0); } } let n = paths as f64; strikes .iter() .zip(&payoffs) .map(|(&k, &total)| implied_vol_black76(total / n, forward, k, expiry, Side::Call)) .collect()} #[cfg(test)]mod gyongi_tests { use super::*; use crate::black::{black76, implied_vol_black76, Side}; const FORWARD: f64 = 100.0; const EXPIRY: f64 = 1.0; fn regime() -> TwoRegime { TwoRegime { weight: 0.6, vol_low: 0.15, vol_high: 0.35 } } /// The mixture's own smile: price each call under both regimes and weight. fn mixture_smile(r: &TwoRegime, strikes: &[f64]) -> Vec<f64> { strikes .iter() .map(|&k| { let price = r.weight * black76(FORWARD, k, r.vol_low, EXPIRY, Side::Call) + (1.0 - r.weight) * black76(FORWARD, k, r.vol_high, EXPIRY, Side::Call); implied_vol_black76(price, FORWARD, k, EXPIRY, Side::Call).unwrap() }) .collect() } #[test] fn the_mimicking_local_volatility_reproduces_the_smile() { // The theorem. A model with one random volatility drawn at the start, and // a local volatility model with no randomness in its volatility at all, // agree on every European price -- which is a strong statement, since the // two processes are nothing like each other path by path. let r = regime(); let strikes = [80.0, 90.0, 100.0, 110.0, 125.0]; let target = mixture_smile(&r, &strikes); let mimicked = mimicking_smile(&r, FORWARD, EXPIRY, &strikes, 400_000, 250, 20260809); for (i, &k) in strikes.iter().enumerate() { let got = mimicked[i].expect("a price should invert"); assert!( (got - target[i]).abs() < 0.006, "K={k}: mimicking gave {got:.4} against the mixture's {:.4}", target[i] ); } } #[test] fn the_smile_being_reproduced_is_a_real_smile() { // Guarding against a vacuous pass. If the mixture's smile were flat, a // constant volatility would reproduce it and the test above would prove // nothing. It is not flat: the mixture is convex in log-strike by a wide // margin, and the conditional variance the local model uses genuinely // varies with the level. let r = regime(); let strikes = [80.0, 100.0, 125.0]; let smile = mixture_smile(&r, &strikes); assert!(smile[0] > smile[1] + 0.01, "the wings should be well above the money"); assert!(smile[2] > smile[1] + 0.01); // And the local volatility it implies is not constant either. let low = r.conditional_variance(FORWARD, EXPIRY, 100.0).sqrt(); let wing = r.conditional_variance(FORWARD, EXPIRY, 140.0).sqrt(); assert!(wing > low + 0.05, "conditional vol {low:.3} at the money, {wing:.3} in the wing"); }} #[cfg(test)]mod mixing_tests { use super::*; const F: f64 = 100.0; const T: f64 = 1.0; const TARGET: f64 = 0.25; fn regime() -> TwoRegime { TwoRegime { weight: 0.5, vol_low: 0.15, vol_high: 0.35 } } /// The dial turns the randomness of the variance and nothing else. #[test] fn mixing_holds_the_mean_variance_and_moves_only_the_spread() { let base = regime(); for &lambda in &[0.0, 0.25, 0.5, 1.0] { let m = base.with_mixing(lambda); assert!( (m.mean_variance() - base.mean_variance()).abs() < 1e-12, "lambda {lambda} moved the mean variance" ); let gap = (m.vol_high.powi(2) - m.vol_low.powi(2)).abs(); let full = (base.vol_high.powi(2) - base.vol_low.powi(2)).abs(); assert!( (gap - lambda * lambda * full).abs() < 1e-12, "lambda {lambda} scaled the variance gap wrongly" ); } // At zero the two states coincide: the variance is deterministic. let off = base.with_mixing(0.0); assert!((off.vol_high - off.vol_low).abs() < 1e-12); } /// What the chapter's mixing section claims, as an identity rather than a /// description: the product of the leverage function and the conditional /// variance is the market's local variance at every setting of the dial, so /// turning it redistributes a fixed product between its two factors and /// changes no European price. #[test] fn the_dial_redistributes_a_product_it_cannot_change() { let base = regime(); for &lambda in &[0.0, 0.25, 0.5, 0.75, 1.0] { let m = base.with_mixing(lambda); for &k in &[80.0, 100.0, 130.0] { let l = m.leverage(F, T, k, TARGET); let product = l * l * m.conditional_variance(F, T, k); assert!( (product - TARGET * TARGET).abs() < 1e-12, "lambda {lambda}, K {k}: L^2 E[v|K] was {product}, not {}", TARGET * TARGET ); } } } /// And the redistribution is not cosmetic. What the dial changes is the /// *shape* of the leverage function, not merely its size. /// /// Against a flat local volatility surface, a deterministic variance leaves /// `L` with nothing strike-dependent to do and it comes out flat. Turn the /// randomness on and the stochastic part generates a smile of its own, so /// `L` has to remove it: it rises at the money and dips in both wings, and a /// function that was constant acquires a range of half its own value. #[test] fn the_dial_changes_the_shape_of_the_leverage_function() { let base = regime(); let strikes = [60.0, 70.0, 80.0, 100.0, 130.0, 160.0, 200.0]; let spread = |lambda: f64| { let m = base.with_mixing(lambda); let ls: Vec<f64> = strikes.iter().map(|&k| m.leverage(F, T, k, TARGET)).collect(); let lo = ls.iter().cloned().fold(f64::INFINITY, f64::min); let hi = ls.iter().cloned().fold(f64::NEG_INFINITY, f64::max); (lo, hi, hi / lo) }; let (lo, hi, ratio) = spread(0.0); assert!(ratio - 1.0 < 1e-9, "with no randomness L should be flat, got {lo}..{hi}"); assert!((lo - 0.9285).abs() < 1e-3, "flat leverage was {lo}"); let (lo, hi, ratio) = spread(1.0); assert!((lo - 0.714).abs() < 0.01, "wing leverage was {lo}"); assert!((hi - 1.094).abs() < 0.01, "at-the-money leverage was {hi}"); assert!((ratio - 1.53).abs() < 0.02, "leverage spread was {ratio}"); }}