quant/src/marketmodels.rs
The Libor market model, and the measurement that decides it against the swap market model.
//! The Libor market model, and the measurement that decides it against the swap//! market model.//!//! The market models chapter argues that forwards and swap rates are two//! coordinate systems on the same curve, and that both cannot be lognormal at//! once: a swap rate is a weighted sum of forwards, and a weighted sum of//! lognormals is not lognormal. So one of the two liquid calibration sets is//! priced by approximation, and the choice of model is the choice of which.//!//! That argument is settled by a number rather than by preference, and this//! module produces it. Simulate a lognormal Libor market model, price swaptions//! across strikes off the resulting swap rate distribution, and invert to Black//! implied volatility. If the resulting smile is flat the approximation costs//! nothing and the forward-based model's generality is free; if it is not, the//! swap-based model prices its own calibrating instruments exactly and the//! forward-based one does not.//!//! # Conventions//!//! A tenor structure `T_0 < T_1 < ... < T_N` with accruals `delta_i = T_{i+1} -//! T_i`. `F_i` is the forward rate for `[T_i, T_{i+1}]`, so there are `N` of//! them. Everything is simulated under the terminal measure, with `P(., T_N)` as//! numeraire, in which//!//! ```text//! dF_i / F_i = -sigma_i sum_{j>i} (delta_j rho_ij sigma_j F_j)/(1 + delta_j F_j) dt//! + sigma_i dW_i.//! ```//!//! The last forward has an empty sum and is therefore an exact martingale, which//! is the cheapest available check on the drift and is tested below. use crate::black::{implied_vol_black76, Side};use crate::pathwise::Rng; /// A flat-volatility, exponentially-correlated Libor market model.////// Deliberately the simplest version that still has the two features the/// chapter needs: a full set of forwards evolving jointly, and a correlation/// between them that decays with the distance between their maturities. The/// deterministic volatility is the point of the exercise rather than a/// simplification --- what is being measured is the smile a model with *no*/// smile in its forwards produces in its swap rates.pub struct Lmm { /// Accrual of each forward period. pub accrual: f64, /// Today's forward rates, one per period. pub forwards: Vec<f64>, /// Lognormal volatility of each forward. pub vols: Vec<f64>, /// Correlation decay: `rho_ij = exp(-beta |T_i - T_j|)`. pub beta: f64,} /// A single variance factor shared by every forward in the family.////// This is the stochastic volatility extension the market models chapter argues/// for: not a separate volatility process per rate, but one process scaling the/// loadings of the whole family together, so that the smile dynamics are/// specified once and jointly rather than assembled from marks.////// The variance starts and mean-reverts to one, so the deterministic `vols`/// keep their meaning as the average level and this factor supplies only the/// randomness around it.#[derive(Clone, Copy, Debug)]pub struct CommonVariance { /// Speed of mean reversion. pub theta: f64, /// Volatility of the variance. Zero recovers the deterministic model. pub eta: f64,} impl Lmm { /// A flat curve and a flat volatility, which is the configuration in which /// any smile that appears cannot have come from the inputs. pub fn flat(periods: usize, accrual: f64, rate: f64, vol: f64, beta: f64) -> Self { Lmm { accrual, forwards: vec![rate; periods], vols: vec![vol; periods], beta, } } fn correlation(&self, i: usize, j: usize) -> f64 { let gap = (i as f64 - j as f64).abs() * self.accrual; (-self.beta * gap).exp() } /// The Cholesky factor of the correlation matrix, so that correlated /// increments can be drawn from independent ones. fn cholesky(&self) -> Vec<Vec<f64>> { let n = self.forwards.len(); let mut l = vec![vec![0.0; n]; n]; for i in 0..n { for j in 0..=i { let mut sum = self.correlation(i, j); for k in 0..j { sum -= l[i][k] * l[j][k]; } if i == j { // Clamped because an exponential correlation matrix is // positive definite in exact arithmetic and can lose the // last digit of it in floating point. l[i][j] = sum.max(1e-14).sqrt(); } else { l[i][j] = sum / l[j][j]; } } } l } /// Evolve every forward to `horizon` and return the state of the curve on /// each path. /// /// Log-Euler, which keeps each forward positive by construction and is exact /// for the diffusion term; the drift is frozen across each step, which is /// the standard discretisation and the only approximation in the scheme. /// /// Paths are drawn in antithetic pairs. That matters more here than usual: /// the quantity being measured is a difference between implied volatilities /// of a few hundredths of a point, and the level of the implied volatility /// is far noisier than the difference. Antithetic sampling removes most of /// the level's noise, and drawing every strike from one set of paths removes /// the rest from the difference. pub fn simulate(&self, horizon: f64, paths: usize, steps: usize, seed: u64) -> Vec<Vec<f64>> { self.simulate_with(None, horizon, paths, steps, seed) } /// The same model with a common variance factor multiplying every forward's /// volatility, as in [`CommonVariance`]. /// /// Each forward's diffusion becomes `sqrt(V) * sigma_i`, and the /// no-arbitrage drift is unchanged in form: it is the same double sum, now /// carrying the same `V`, because the drift is quadratic in the volatilities /// it is built from. Nothing about the derivation of the drift has to be /// redone --- which is the honest answer to whether stochastic volatility /// complicates the drift. pub fn simulate_stochastic_vol( &self, variance: CommonVariance, horizon: f64, paths: usize, steps: usize, seed: u64, ) -> Vec<Vec<f64>> { self.simulate_with(Some(variance), horizon, paths, steps, seed) } fn simulate_with( &self, variance: Option<CommonVariance>, horizon: f64, paths: usize, steps: usize, seed: u64, ) -> Vec<Vec<f64>> { let n = self.forwards.len(); let l = self.cholesky(); let dt = horizon / steps as f64; let root_dt = dt.sqrt(); let mut rng = Rng::new(seed); let mut out = Vec::with_capacity(paths); let mut draws: Vec<Vec<f64>> = Vec::with_capacity(steps); for path in 0..paths { let antithetic = path % 2 == 1; if !antithetic { draws.clear(); for _ in 0..steps { // One extra draw per step for the variance factor. draws.push((0..=n).map(|_| rng.next_normal()).collect()); } } let mut f = self.forwards.clone(); let mut v = 1.0f64; for step in 0..steps { let z: Vec<f64> = if antithetic { draws[step].iter().map(|x| -x).collect() } else { draws[step].clone() }; let mut next = f.clone(); // sqrt(V) scales every forward's volatility together. let scale = v.max(0.0).sqrt(); for i in 0..n { // Terminal measure drift: everything with a later maturity // pushes this forward down. let mut drift = 0.0; for j in (i + 1)..n { drift -= self.correlation(i, j) * self.vols[j] * self.accrual * f[j] / (1.0 + self.accrual * f[j]); } drift *= self.vols[i]; // Quadratic in the volatilities, so it carries V, not sqrt(V). drift *= scale * scale; let dw: f64 = (0..=i).map(|k| l[i][k] * z[k]).sum(); let vol = self.vols[i] * scale; next[i] = f[i] * ((drift - 0.5 * vol * vol) * dt + vol * root_dt * dw).exp(); } f = next; if let Some(cv) = variance { // Full truncation, which keeps the variance from going // negative without distorting its mean the way reflection // would. let root_v = v.max(0.0).sqrt(); v += cv.theta * (1.0 - v) * dt + cv.eta * root_v * root_dt * z[n]; } } out.push(f); } out } /// Discount factor from `T_from` to `T_to`, both given as indices into the /// tenor structure, built from the forwards on one path. fn bond(&self, forwards: &[f64], from: usize, to: usize) -> f64 { (from..to).map(|j| 1.0 / (1.0 + self.accrual * forwards[j])).product() } /// Today's par swap rate for the swap over `[T_a, T_b]`. pub fn par_swap_rate(&self, a: usize, b: usize) -> f64 { let f = &self.forwards; let annuity: f64 = (a + 1..=b).map(|k| self.accrual * self.bond(f, 0, k)).sum(); (self.bond(f, 0, a) - self.bond(f, 0, b)) / annuity }} /// One point of the measured swaption smile.pub struct SmilePoint { pub strike: f64, /// Black implied volatility of the simulated swaption price. pub implied_vol: f64,} /// Price swaptions across strikes in the lognormal model and invert to implied/// volatility.////// The swaption expires at `T_a` on the swap over `[T_a, T_b]`. Simulation is/// under the terminal measure, so the payoff is divided by the numeraire/// `P(T_a, T_N)` rebuilt on each path and multiplied by `P(0, T_N)`.////// If forwards being lognormal made swap rates lognormal, every point of the/// returned curve would carry the same implied volatility. The extent to which/// it does not is the price of using a forward-based model to quote swaptions.pub fn swaption_smile( model: &Lmm, a: usize, b: usize, strikes: &[f64], paths: usize, steps: usize, seed: u64,) -> Vec<SmilePoint> { smile_with_variance(model, None, a, b, strikes, paths, steps, seed)} /// The same measurement with a common variance factor switched on.pub fn smile_with_variance( model: &Lmm, variance: Option<CommonVariance>, a: usize, b: usize, strikes: &[f64], paths: usize, steps: usize, seed: u64,) -> Vec<SmilePoint> { let n = model.forwards.len(); let expiry = a as f64 * model.accrual; let states = match variance { Some(cv) => model.simulate_stochastic_vol(cv, expiry, paths, steps, seed), None => model.simulate(expiry, paths, steps, seed), }; let numeraire_today = model.bond(&model.forwards, 0, n); let forward = model.par_swap_rate(a, b); strikes .iter() .map(|&strike| { let mut total = 0.0; for f in &states { let annuity: f64 = (a + 1..=b).map(|k| model.accrual * model.bond(f, a, k)).sum(); let swap = (1.0 - model.bond(f, a, b)) / annuity; let payoff = annuity * (swap - strike).max(0.0); // Deflate by the terminal bond, which is the numeraire. total += payoff / model.bond(f, a, n); } let price = numeraire_today * total / paths as f64; let implied_vol = implied_vol_black76(price / annuity_today(model, a, b), forward, strike, expiry, Side::Call) .unwrap_or(f64::NAN); SmilePoint { strike, implied_vol } }) .collect()} /// Today's annuity for the swap over `[T_a, T_b]`, which is what a swaption/// price is quoted against.pub fn annuity_today(model: &Lmm, a: usize, b: usize) -> f64 { (a + 1..=b) .map(|k| model.accrual * model.bond(&model.forwards, 0, k)) .sum()} /// How far the measured smile departs from flat, in volatility points.////// Returned as the spread between the highest and lowest implied volatility/// across the strikes supplied, which is the number the chapter quotes.pub fn smile_spread(smile: &[SmilePoint]) -> f64 { let hi = smile.iter().map(|p| p.implied_vol).fold(f64::MIN, f64::max); let lo = smile.iter().map(|p| p.implied_vol).fold(f64::MAX, f64::min); hi - lo} #[cfg(test)]mod tests { use super::*; use crate::black::black76; fn model() -> Lmm { // Ten semiannual periods, a flat 4% curve, 25% lognormal volatility on // every forward, correlation decaying with a five year scale. Lmm::flat(10, 0.5, 0.04, 0.25, 0.2) } #[test] fn the_last_forward_is_a_martingale_under_the_terminal_measure() { // Its drift sum is empty by construction, so this is a direct test that // the drift is indexed the way the derivation says. An error in the // limits of the sum would show up here first. let m = model(); let n = m.forwards.len(); let states = m.simulate(2.0, 40_000, 40, 20260810); let mean: f64 = states.iter().map(|f| f[n - 1]).sum::<f64>() / states.len() as f64; assert!( (mean / m.forwards[n - 1] - 1.0).abs() < 0.01, "last forward drifted: {mean:.6} against {:.6}", m.forwards[n - 1] ); } #[test] fn the_earlier_forwards_are_pushed_down() { // And the ones with maturities before the numeraire's are not // martingales in this measure: the terminal measure's drift is // one-signed, so every forward but the last has a lower mean than it // started with, and the effect grows with distance from the numeraire. let m = model(); let states = m.simulate(2.0, 40_000, 40, 20260810); let mean = |i: usize| { states.iter().map(|f| f[i]).sum::<f64>() / states.len() as f64 }; assert!(mean(0) < m.forwards[0], "the first forward should drift down"); assert!( mean(0) < mean(m.forwards.len() - 2), "the drift should grow with distance from the numeraire" ); } #[test] fn a_caplet_prices_to_black() { // The model is lognormal in each forward under that forward's own // measure, so a caplet has to come out at exactly its Black price. This // is the test that the change of numeraire is right: the drift and the // deflator have to cancel, and they only cancel if both are correct. let m = model(); let n = m.forwards.len(); let i = 4; // caplet on F_4, expiring at T_4 let expiry = i as f64 * m.accrual; let states = m.simulate(expiry, 60_000, 60, 20260811); let strike = 0.04; let numeraire_today = m.bond(&m.forwards, 0, n); let mut total = 0.0; for f in &states { let payoff = m.accrual * (f[i] - strike).max(0.0); // Pays at T_{i+1}, so it is worth that much of a bond then; deflate // by the terminal numeraire from T_i. total += payoff * m.bond(f, i, i + 1) / m.bond(f, i, n); } let simulated = numeraire_today * total / states.len() as f64; let discount = m.bond(&m.forwards, 0, i + 1); let exact = discount * m.accrual * black76(m.forwards[i], strike, m.vols[i], expiry, Side::Call); assert!( (simulated / exact - 1.0).abs() < 0.02, "caplet simulated {simulated:.8} against Black {exact:.8}" ); } #[test] fn a_lognormal_forward_model_leaves_swap_rates_nearly_lognormal() { // The measurement the chapter is built on, and the honest form of it. // // Forwards here are exactly lognormal and have no smile whatever. The // swap rate is a weighted sum of them, so it cannot be lognormal, and a // Black implied volatility read off its option prices therefore cannot // be flat across strikes. The question is how far from flat, because // that is what a forward-based model costs when it is used to quote // swaptions. // // The answer is: less than this test can resolve. Run harder -- two // million paths against the hundred thousand here, over several seeds -- // a consistent shape does appear, implied volatility rising with strike, // spanning about three hundredths of a volatility point across strikes // from seventy to a hundred and thirty per cent of the forward. At the // effort a test can afford, the shape is not stable from seed to seed // and only the bound is. So the bound is what is asserted and what the // chapter quotes, since a claim that needs eight million paths to see is // not a claim about anything a desk trades. let m = model(); let (a, b) = (4, 10); let atm = m.par_swap_rate(a, b); let strikes: Vec<f64> = [0.7, 0.85, 1.0, 1.15, 1.3].iter().map(|s| s * atm).collect(); for seed in [20260812u64, 7, 99, 12345] { let smile = swaption_smile(&m, a, b, &strikes, 120_000, 60, seed); for p in &smile { assert!(p.implied_vol.is_finite(), "no implied vol at {:.4}", p.strike); } let spread = smile_spread(&smile); assert!( spread < 0.001, "seed {seed}: induced smile spans {:.4} vol points, which is more \ than a tenth of a point and would be a real cost", spread * 100.0 ); } } #[test] fn the_swap_rate_is_less_volatile_than_the_forwards_it_averages() { // The other half of the same fact, and the one that is not small. Every // forward has 25% volatility; the swap rate averaging them has less, // because they are imperfectly correlated and an average of imperfectly // correlated things moves less than its parts. // // This is why decorrelation is the first-order modelling choice in a // market model and the smile is a second-order one: the correlation // parameter moves the at-the-money volatility by percentage points, // while the departure from lognormality moves the smile by hundredths. let m = model(); let (a, b) = (4, 10); let atm = m.par_swap_rate(a, b); let smile = swaption_smile(&m, a, b, &[atm], 120_000, 60, 20260812); let swaption_vol = smile[0].implied_vol; assert!( swaption_vol < 0.25, "the average cannot be more volatile than its parts: {swaption_vol:.4}" ); assert!( (0.20..0.245).contains(&swaption_vol), "expected a modest haircut to 25%, got {swaption_vol:.4}" ); // And it is the correlation that does it: raise the decay and the // forwards decorrelate further, so the swap rate calms down further. let decorrelated = Lmm::flat(10, 0.5, 0.04, 0.25, 0.6); let looser = swaption_smile(&decorrelated, a, b, &[atm], 120_000, 60, 20260812); assert!( looser[0].implied_vol < swaption_vol - 0.005, "more decorrelation should mean less swaption volatility: {:.4} against {:.4}", looser[0].implied_vol, swaption_vol ); }} #[cfg(test)]mod stochastic_vol_tests { use super::*; fn model() -> Lmm { Lmm::flat(10, 0.5, 0.04, 0.25, 0.2) } /// The contrast the market models chapter closes on. /// /// A lognormal forward market model induces a swaption smile of a few /// hundredths of a volatility point, which is the residue of using the wrong /// coordinates and is one to two orders of magnitude below a bid-offer. /// Attaching one common variance factor to the whole family produces a smile /// of well over a point, which is a market-sized number. So the missing /// smile is not a coordinate problem and cannot be fixed by choosing /// coordinates; it is a missing factor. #[test] fn a_common_variance_factor_produces_a_market_sized_smile() { let m = model(); let f = m.par_swap_rate(4, 10); let strikes: Vec<f64> = (0..9).map(|i| f * (0.70 + 0.075 * i as f64)).collect(); let spread = |eta: f64| { let cv = (eta > 0.0).then_some(CommonVariance { theta: 0.5, eta }); smile_spread(&smile_with_variance(&m, cv, 4, 10, &strikes, 100_000, 60, 4242)) }; let flat = spread(0.0); assert!(flat < 0.001, "the deterministic model showed {flat:.4} of smile"); let gentle = spread(0.6); let strong = spread(1.2); assert!(gentle > flat, "vol-of-vol 0.6 gave {gentle:.4} against {flat:.4}"); assert!(strong > gentle, "vol-of-vol 1.2 gave {strong:.4} against {gentle:.4}"); assert!( (0.012..0.018).contains(&strong), "vol-of-vol 1.2 gave a smile of {strong:.4}, expected about 0.014" ); // More than an order of magnitude, which is the point of the // comparison; the ratio itself is limited by the noise floor of the // deterministic case rather than by anything in the model. assert!(strong / flat > 15.0, "ratio was only {:.1}", strong / flat); } /// The drift keeps its shape. Turning the variance factor off by setting the /// vol-of-vol to zero must reproduce the deterministic model exactly, which /// checks that the factor was threaded through the drift as well as through /// the diffusion — the drift is quadratic in the volatilities, so it carries /// `V` and not its square root. #[test] fn zero_vol_of_vol_recovers_the_deterministic_model() { let m = model(); let f = m.par_swap_rate(4, 10); let strikes: Vec<f64> = (0..5).map(|i| f * (0.85 + 0.075 * i as f64)).collect(); let off = smile_with_variance(&m, None, 4, 10, &strikes, 20_000, 40, 7); let zero = smile_with_variance( &m, Some(CommonVariance { theta: 0.5, eta: 0.0 }), 4, 10, &strikes, 20_000, 40, 7, ); for (a, b) in off.iter().zip(&zero) { assert!( (a.implied_vol - b.implied_vol).abs() < 1e-12, "strike {}: {} against {}", a.strike, a.implied_vol, b.implied_vol ); } }}