quant/src/quasigaussian.rs
The smile a quasi-Gaussian model produces, and where it comes from.
//! The smile a quasi-Gaussian model produces, and where it comes from.//!//! The Markovian term structure chapter claims that the shape of `sigma_r`, the//! model's volatility as a function of its own state, is where the smile lives://! constant gives none, a linear term gives a skew, a quadratic term gives//! curvature. This computes the smile so the claim can be looked at rather than//! taken on trust.//!//! # How, and why not by simulation//!//! The state `x` is a one-dimensional diffusion, so its distribution at any date//! solves the Fokker-Planck equation of the generator chapter — the same forward equation//! Dupire's formula was derived from, now being used forwards instead of//! backwards. Solving it on a grid gives the density directly, and by//! Breeden-Litzenberger in reverse the option prices are integrals against it.//!//! That is both exact, up to discretisation, and fast enough to redraw while a//! slider moves, which Monte Carlo would not be. It also keeps the figure//! honest about its subject: the smile is a property of the terminal//! distribution, and here it is computed from the terminal distribution.//!//! # What is left out//!//! The auxiliary state `y` of the Markovian term structure chapter is held on//! the deterministic path it would follow if `sigma_r` were constant. `y` has//! no diffusion term and enters only the drift, so it moves the level of the//! rate and not the shape of the smile, which is what the figure is about.//! Everything else is exact. use crate::black::{implied_vol_bachelier, Side}; /// The volatility of the state, as a polynomial in the state itself.////// Written in terms of a normalised state `u = x / scale`, so that `b` and `c`/// are dimensionless dials rather than quantities carrying powers of a rate./// `scale` is the size of the move over which the polynomial is meant to bend —/// 100 basis points is the natural unit for rates.#[derive(Clone, Copy, Debug)]pub struct QuasiGaussian { /// Overall level of volatility, in absolute rate terms. pub sigma: f64, /// Constant term. One, in the usual parametrisation. pub a: f64, /// Linear term: the skew dial. Zero recovers Hull-White. pub b: f64, /// Quadratic term: the smile dial. pub c: f64, /// Mean reversion. pub kappa: f64, /// The rate move over which `b` and `c` are measured. /// /// Set this to the typical move, `sigma * sqrt(T)`, and the normalised state /// `u = x / scale` is of order one across the range the process actually /// explores — which makes `b` and `c` dials with comparable, interpretable /// strength. Setting it to a fixed number of basis points instead makes `u` /// large at long expiries and the quadratic term swamps everything, which is /// a statement about the parametrisation rather than about the model. pub scale: f64,} impl Default for QuasiGaussian { fn default() -> Self { QuasiGaussian { sigma: 0.01, a: 1.0, b: 0.0, c: 0.0, kappa: 0.03, scale: 0.01, } }} impl QuasiGaussian { /// `sigma_r(x)`, floored. /// /// The floor is not decoration. The Markovian term structure chapter warns /// that a quadratic diffusion coefficient can go negative and can grow fast /// enough to destroy existence and uniqueness of the solution; production /// implementations cap and floor it. Doing the same here means the figure /// shows the model as it is actually used, and cannot be driven to nonsense /// by a slider. pub fn sigma_r(&self, x: f64) -> f64 { let u = x / self.scale; let raw = self.sigma * (self.a + self.b * u + self.c * u * u); // Between a twentieth and five times the level, which is the range a // desk would allow. raw.clamp(0.05 * self.sigma, 5.0 * self.sigma) } /// The density of `x_T`, by solving the Fokker-Planck equation forward. /// /// Returns the grid and the density on it. The scheme is explicit and /// conservative — it transports mass between cells rather than merely /// satisfying the equation pointwise — so the total probability stays at /// one to machine precision, which is the cheapest available check that the /// solve is behaving. pub fn density(&self, t: f64, points: usize) -> (Vec<f64>, Vec<f64>) { // Wide enough that the tails are numerically zero at the boundary. // // Sized from the volatility over the range the process plausibly // reaches — three standard deviations of the base move — rather than // from the maximum anywhere. Taking the maximum over the whole line // includes the region where the clamp has pinned `sigma_r` at its // ceiling, and sizing to that produces a domain many times too wide, // hence a grid too coarse to resolve anything. let base = self.sigma * t.max(1e-6).sqrt(); let reach = 3.0 * base; let peak = (0..64) .map(|i| -reach + 2.0 * reach * i as f64 / 63.0) .map(|x| self.sigma_r(x)) .fold(0.0, f64::max); let width = 6.0 * peak * t.max(1e-6).sqrt(); let dx = 2.0 * width / (points - 1) as f64; let grid: Vec<f64> = (0..points).map(|i| -width + i as f64 * dx).collect(); // Start from a narrow Gaussian rather than a spike on one node: a spike // is a delta of height 1/dx, which no grid resolves and no explicit // scheme handles gracefully. // // The start time is chosen so that the Gaussian is at least a few cells // wide. Any narrower and it is not resolved, its discrete mass is not // one, and — since the scheme below conserves whatever it is given — the // error never washes out. It is corrected exactly at the end regardless, // but a resolved start is the difference between a smooth density and a // staircase. let sigma0 = self.sigma_r(0.0); let t0 = (t * 0.01).max((3.0 * dx / sigma0).powi(2)).min(t * 0.5); let v0 = sigma0 * sigma0 * t0; let mut p: Vec<f64> = grid .iter() .map(|&x| (-0.5 * x * x / v0).exp() / (2.0 * std::f64::consts::PI * v0).sqrt()) .collect(); // Normalise on the grid, not analytically. The scheme conserves exactly // what it starts with, so starting from something whose discrete mass is // 0.998 gives a "density" that integrates to 0.998 forever. let mass0: f64 = p.iter().sum::<f64>() * dx; for q in p.iter_mut() { *q /= mass0; } // Explicit stability: dt <= dx^2 / (2 max sigma^2), with room to spare. let max_var = grid.iter().map(|&x| self.sigma_r(x).powi(2)).fold(0.0, f64::max); let drift_max = self.kappa * width; let dt_diff = 0.4 * dx * dx / max_var.max(1e-12); let dt_drift = 0.4 * dx / drift_max.max(1e-12); let steps = (((t - t0) / dt_diff.min(dt_drift)).ceil() as usize).max(1); let dt = (t - t0) / steps as f64; // Written as fluxes between cells rather than as derivatives at points. // // The equation is a conservation law, dp/dt = -dF/dx with // F = mu p - d(D p)/dx, and differencing the flux at the interfaces // makes the scheme conserve mass exactly: what leaves cell i enters // cell i+1, by construction, whatever the coefficients do. That matters // here because `sigma_r` is clamped and so has a kink, and the second // derivative of `D p` taken naively across that kink is unbounded — the // resulting scheme creates probability out of nothing. // // With no flux at either end the total is conserved to machine // precision, which turns the mass check in the tests into a real test of // the solve rather than a test of how wide the domain is. let d: Vec<f64> = grid.iter().map(|&x| 0.5 * self.sigma_r(x).powi(2)).collect(); let mut flux = vec![0.0; points + 1]; let mut next = p.clone(); for _ in 0..steps { for i in 1..points { // Interface between cells i-1 and i. let x_mid = 0.5 * (grid[i - 1] + grid[i]); let advect = -self.kappa * x_mid * 0.5 * (p[i - 1] + p[i]); let diffuse = (d[i] * p[i] - d[i - 1] * p[i - 1]) / dx; flux[i] = advect - diffuse; } // Closed at both ends. The domain is sized so nothing reaches them. flux[0] = 0.0; flux[points] = 0.0; for i in 0..points { next[i] = p[i] - dt * (flux[i + 1] - flux[i]) / dx; } p.copy_from_slice(&next); } (grid, p) } /// Total probability and mean of a density, for checking a solve. fn moments(grid: &[f64], p: &[f64]) -> (f64, f64) { let dx = grid[1] - grid[0]; let mass: f64 = p.iter().sum::<f64>() * dx; let mean: f64 = grid.iter().zip(p).map(|(x, q)| x * q).sum::<f64>() * dx; (mass, mean) } /// The implied volatility smile of options on the rate, in Bachelier /// (normal) terms — which is what the rates market quotes in, and what stays /// meaningful when a rate can be negative. /// /// `strikes` are offsets from today's forward rate, so zero is at the money. /// The forward is taken from the solved density rather than assumed, so the /// inversion is against the model's own forward whatever the drift does. pub fn smile(&self, t: f64, strikes: &[f64], points: usize) -> Vec<Option<f64>> { let (grid, p) = self.density(t, points); let dx = grid[1] - grid[0]; let (mass, forward) = Self::moments(&grid, &p); if !(mass > 0.5) { return vec![None; strikes.len()]; } strikes .iter() .map(|&k| { // Price the out-of-the-money side, as everywhere else in this // crate: it is where the information about volatility is. let side = if k >= forward { Side::Call } else { Side::Put }; let price: f64 = grid .iter() .zip(&p) .map(|(&x, &q)| { let payoff = match side { Side::Call => (x - k).max(0.0), Side::Put => (k - x).max(0.0), }; payoff * q }) .sum::<f64>() * dx / mass; implied_vol_bachelier(price, forward, k, t, side) }) .collect() }} #[cfg(test)]mod tests { use super::*; const T: f64 = 5.0; const N: usize = 401; /// A model with `scale` set to the typical move, which is the /// parametrisation the doc comment on `scale` argues for. fn model(b: f64, c: f64) -> QuasiGaussian { let base = QuasiGaussian::default(); QuasiGaussian { b, c, scale: base.sigma * T.sqrt(), ..base } } fn strikes() -> Vec<f64> { // Plus or minus 150 basis points around the forward. (0..13).map(|i| -0.015 + i as f64 * 0.0025).collect() } #[test] fn the_solved_density_is_a_probability_distribution() { let m = model(-0.3, 0.25); let (grid, p) = m.density(T, N); let (mass, _) = QuasiGaussian::moments(&grid, &p); // Exactly, not approximately: the scheme is conservative and the ends // are closed, so anything else is a bug rather than a discretisation. assert!((mass - 1.0).abs() < 1e-9, "mass was {mass}"); assert!(p.iter().all(|&q| q >= -1e-9), "density went negative"); } #[test] fn constant_volatility_gives_a_flat_smile() { // Hull-White is the b = c = 0 corner of the model, and it is Gaussian, // so its Bachelier smile must be flat — and at the level the model was // given, since sigma_r is then an absolute rate volatility. let m = QuasiGaussian { kappa: 0.0, ..model(0.0, 0.0) }; let vols = m.smile(T, &strikes(), N); for (k, v) in strikes().iter().zip(&vols) { let v = v.expect("a vol exists at every strike"); assert!( (v - m.sigma).abs() < 5e-5, "at k={k} the vol was {v}, expected {}", m.sigma ); } } #[test] fn the_smile_bottoms_out_where_averaging_says() { // The Markovian term structure chapter reads the vertex off the // averaging rule: sigma_r(u) = sigma(1 + b u + c u^2) averages to // sigma(1 + b u / 2 + c u^2 / 3), so the smile is least at // // u_min = -3b / 4c, // // three halves as far out as sigma_r's own minimum at -b/2c, because // averaging divides the bend by three and the tilt only by two. // // Checked against the solved density, which shares none of that // reasoning. The strike range is deliberately wider than the figure's: // at the gentler curvature the minimum lies outside what is plotted. let (sigma, expiry) = (0.01f64, 5.0f64); let scale = sigma * expiry.sqrt(); let strikes: Vec<f64> = (0..1601).map(|i| -0.06 + 0.12 * i as f64 / 1600.0).collect(); for (b, c) in [(-0.3, 0.2), (-0.3, 0.5), (-0.6, 0.5), (-0.2, 0.4)] { let m = QuasiGaussian { sigma, scale, b, c, ..Default::default() }; let vols = m.smile(expiry, &strikes, 1601); let (mut lowest, mut at) = (f64::MAX, f64::NAN); for (&k, v) in strikes.iter().zip(&vols) { if let Some(v) = v { if *v < lowest { lowest = *v; at = k; } } } let measured = at / scale; let predicted = -3.0 * b / (4.0 * c); assert!( (measured - predicted).abs() < 0.05 * predicted, "b={b} c={c}: minimum at u={measured:.3}, expected {predicted:.3}" ); // Above the forward, not at it, and not below. assert!(measured > 0.1, "b={b} c={c}: minimum at u={measured:.3} is at the money"); // And further out than sigma_r's own vertex, by about half again. let local = -b / (2.0 * c); let ratio = measured / local; assert!( (ratio - 1.5).abs() < 0.1, "b={b} c={c}: smile vertex is {ratio:.2} times sigma_r's, expected 1.5" ); } } #[test] fn the_linear_term_tilts_and_the_quadratic_term_bends() { // The claim the Markovian term structure chapter's table makes, as two // measurements. Slope and curvature are the first and second // differences of the smile. let ks = strikes(); let shape = |m: &QuasiGaussian| { let v: Vec<f64> = m.smile(T, &ks, N).into_iter().map(|x| x.unwrap()).collect(); let n = v.len(); let slope = v[n - 1] - v[0]; let curvature = v[0] - 2.0 * v[n / 2] + v[n - 1]; (slope, curvature) }; let (flat_slope, flat_curve) = shape(&model(0.0, 0.0)); assert!(flat_slope.abs() < 1e-4 && flat_curve.abs() < 1e-4, "b=c=0 should be flat"); let (skew_slope, skew_curve) = shape(&model(-0.4, 0.0)); assert!(skew_slope < -1e-4, "a negative b should tilt the smile down, got {skew_slope}"); assert!(skew_curve.abs() < skew_slope.abs(), "b alone should tilt more than it bends"); let (_, smile_curve) = shape(&model(0.0, 0.4)); assert!(smile_curve > 1e-4, "a positive c should bend the smile up, got {smile_curve}"); } #[test] fn more_curvature_in_the_state_means_more_smile() { // Monotone in the dial, which is what makes it usable as a dial. let ks = strikes(); let curvature = |c: f64| { let m = model(0.0, c); let v: Vec<f64> = m.smile(T, &ks, N).into_iter().map(|x| x.unwrap()).collect(); v[0] - 2.0 * v[v.len() / 2] + v[v.len() - 1] }; let mut previous = f64::NEG_INFINITY; for c in [0.0, 0.2, 0.4, 0.6] { let k = curvature(c); assert!(k > previous, "curvature did not increase at c={c}: {k} after {previous}"); previous = k; } } #[test] fn the_smile_is_the_shape_of_sigma_r_averaged() { // The local volatility chapter's rule of two, arriving in the rates // model: the smile's slope should be about half the slope of sigma_r, // because implied volatility averages the local volatility over the // journey. let m = QuasiGaussian { kappa: 0.0, ..model(-0.4, 0.0) }; let h = 0.005; let local_slope = (m.sigma_r(h) - m.sigma_r(-h)) / (2.0 * h); let vols = m.smile(T, &[-h, 0.0, h], N); let implied_slope = (vols[2].unwrap() - vols[0].unwrap()) / (2.0 * h); let ratio = local_slope / implied_slope; assert!( (ratio - 2.0).abs() < 0.35, "local slope {local_slope:.5} over implied slope {implied_slope:.5} is {ratio:.3}" ); }} /// The rank of a two-variable kernel, sampled on a grid.////// The Markovian term structure chapter shows that a Heath-Jarrow-Morton model/// with a deterministic volatility is Markov in `n` state variables exactly when/// the kernel `sigma(u, T)` has rank `n` --- that is, exactly when it is a sum of/// `n` products `phi_i(u) psi_i(T)`. So the number of states a volatility needs/// is a rank, and can be computed rather than argued about.////// Computed by pivoted Gram-Schmidt over the columns, which is a rank-revealing/// factorisation and needs nothing but dot products.pub fn kernel_rank(kernel: impl Fn(f64, f64) -> f64, tolerance: f64) -> usize { const N: usize = 40; const SPAN: f64 = 20.0; let grid: Vec<f64> = (0..N).map(|i| SPAN * i as f64 / (N as f64 - 1.0)).collect(); let mut columns: Vec<Vec<f64>> = grid .iter() .map(|&y| grid.iter().map(|&x| kernel(x, y)).collect()) .collect(); let norm = |v: &[f64]| v.iter().map(|x| x * x).sum::<f64>().sqrt(); let scale = columns.iter().map(|c| norm(c)).fold(0.0, f64::max); if scale == 0.0 { return 0; } let mut rank = 0; while rank < N { // Pivot on whatever is left standing furthest out of the span so far. let (best, size) = columns .iter() .enumerate() .map(|(i, c)| (i, norm(c))) .fold((0, 0.0), |acc, x| if x.1 > acc.1 { x } else { acc }); if size <= tolerance * scale { break; } let q: Vec<f64> = columns[best].iter().map(|v| v / size).collect(); for column in columns.iter_mut() { let overlap: f64 = column.iter().zip(&q).map(|(c, q)| c * q).sum(); for (c, q) in column.iter_mut().zip(&q) { *c -= overlap * q; } } rank += 1; } rank} /// How many state variables a *time-homogeneous* maturity profile needs.////// When the volatility depends on the maturity only through the time left to it,/// `sigma(u,T) = sigma_r(u) g(T-u)`, the kernel's sections are the shifts of `g`:/// writing `a = T - t` for how far ahead we look and `b = t - u` for how long ago/// it happened, the section is `b -> g(a+b)`. So this is [`kernel_rank`] of/// `(b,a) -> g(a+b)`, and it counts the dimension of the span of the shifts.////// The chapter's theorem is that this rank is finite exactly for the/// quasi-exponentials --- sums of `x^k exp(lambda x)`, with `lambda` allowed to/// be complex --- and that the rank is the order of the constant-coefficient/// differential equation the profile satisfies.////// Note how much of the work time homogeneity is doing. The same `g` used as a/// profile in *calendar* time rather than in time to maturity gives a kernel of/// rank one whatever `g` is, which the tests below show; it is the demand that/// the shape ride along with the maturity that produces the exponential.pub fn realisation_dimension(shape: impl Fn(f64) -> f64, tolerance: f64) -> usize { kernel_rank(|b, a| shape(a + b), tolerance)} /// The loading of an `n`-year par swap rate on a one-factor Gaussian state.////// A swap rate is roughly the average of the instantaneous forwards over the/// swap's life, and in a one-factor model each of those loads on the state by/// `exp(-kappa u)`, so the swap rate loads by the average of that:////// ```text/// beta(n) = (1 - exp(-kappa n)) / (kappa n)./// ```////// An approximation --- it ignores the annuity weighting --- but the point it/// is used for is a ratio of two such loadings, where the weighting largely/// cancels.pub fn swap_loading(kappa: f64, tenor: f64) -> f64 { if kappa * tenor < 1e-12 { return 1.0; } (1.0 - (-kappa * tenor).exp()) / (kappa * tenor)} /// The volatility of the spread between two rates, given their own volatilities/// and their correlation.////// `sqrt(s1^2 + s2^2 - 2 rho s1 s2)`, which is worth writing down because of/// what it does at `rho = 1`: it becomes `|s1 - s2|`, and that is the *minimum*/// over all correlations, not zero.////// A one-factor term structure model forces `rho = 1`, since every rate is an/// increasing function of the same state. So it does not make a spread/// deterministic --- the two rates load on the state by different amounts and/// the difference still moves --- but it does give the spread the least/// volatility compatible with the two rates' own, and offers no parameter with/// which to raise it.pub fn spread_volatility(sigma_one: f64, sigma_two: f64, rho: f64) -> f64 { (sigma_one * sigma_one + sigma_two * sigma_two - 2.0 * rho * sigma_one * sigma_two) .max(0.0) .sqrt()} #[cfg(test)]mod spread_tests { use super::*; const KAPPA: f64 = 0.03; #[test] fn one_factor_does_not_kill_a_spread_it_floors_it() { // The claim the Markovian term structure chapter makes. Perfect // correlation is the minimum of the spread's variance over rho, and the // minimum is not zero unless the two rates load identically. let (two, ten) = (swap_loading(KAPPA, 2.0), swap_loading(KAPPA, 10.0)); assert!(two > ten, "a shorter swap should load more heavily"); let floor = spread_volatility(two, ten, 1.0); assert!(floor > 0.0, "the spread still moves"); assert!((floor - (two - ten)).abs() < 1e-12, "and equals the difference of loadings"); // It really is the minimum: nothing below rho = 1 is smaller. for rho in [-1.0, -0.5, 0.0, 0.5, 0.9, 0.99] { assert!(spread_volatility(two, ten, rho) > floor); } // About a ninth of the two year rate's own volatility, so small but not // nothing --- which is what the chapter used to claim. let share = floor / two; assert!((share - 0.11).abs() < 0.01, "spread vol is {share:.3} of the two year's"); } #[test] fn the_shortfall_is_large_where_the_trade_actually_is() { // And how much is given up. The standard constant maturity swap spread // trade is two year against ten, which is exactly where a one-factor // model is furthest from a plausible correlation: the loadings are // close, so their difference is small, while a rho below one leaves the // two volatilities almost uncancelled. let (two, ten) = (swap_loading(KAPPA, 2.0), swap_loading(KAPPA, 10.0)); let ratio = spread_volatility(two, ten, 0.9) / spread_volatility(two, ten, 1.0); assert!((ratio - 4.0).abs() < 0.2, "expected about four times, got {ratio:.2}"); // Wider spreads suffer less, because the loadings differ enough that // the perfectly correlated case already leaves something behind. let thirty = swap_loading(KAPPA, 30.0); let wide = spread_volatility(two, thirty, 0.9) / spread_volatility(two, thirty, 1.0); assert!(wide < ratio, "the wide spread should be less distorted"); assert!((wide - 1.5).abs() < 0.2, "expected about half again, got {wide:.2}"); }} #[cfg(test)]mod realisation_tests { use super::*; const TOL: f64 = 1e-8; #[test] fn a_single_exponential_needs_one_state() { // The case the Markovian term structure chapter builds on: Hull-White, // Cheyette, and every model in the chapter. One state variable per // factor. for kappa in [0.02, 0.05, 0.1, 0.2] { assert_eq!(realisation_dimension(|x| (-kappa * x).exp(), TOL), 1, "kappa={kappa}"); } // A flat shape is the same statement at kappa = 0. assert_eq!(realisation_dimension(|_| 1.0, TOL), 1); } #[test] fn a_hump_needs_two() { // x exp(-kappa x), the shape a volatility takes when it peaks at some // tenor rather than decaying from the front. Still finite, still // Markov, but it costs a second state variable --- which is the answer // to whether the exponential is the only shape that works. let kappa = 0.05; assert_eq!(realisation_dimension(|x| x * (-kappa * x).exp(), TOL), 2); // And two decay rates cost two states for the more obvious reason. assert_eq!(realisation_dimension(|x| (-0.05 * x).exp() + (-0.3 * x).exp(), TOL), 2); // A complex pair is also two: the same theorem, with the exponential's // rate off the real axis, giving a damped oscillation. assert_eq!( realisation_dimension(|x| (-0.05 * x).exp() * (0.4 * x).cos(), TOL), 2 ); } #[test] fn the_dimension_counts_the_terms() { // Systematically: x^k exp(-kappa x) for k up to three, which spans a // space of dimension k+1 because the shifts of x^k reach every lower // power. let kappa = 0.05; for k in 0..4u32 { let d = realisation_dimension(|x: f64| x.powi(k as i32) * (-kappa * x).exp(), TOL); assert_eq!(d, k as usize + 1, "x^{k} exp(-kappa x)"); } } #[test] fn the_two_states_rebuild_the_whole_curve() { // Cheyette's theorem, checked against a direct integration of the // Heath-Jarrow-Morton equation that knows nothing about x or y. // // This is the reconstruction the Markovian term structure chapter // derives, // // f(t,T) = f(0,T) + e^{-k(T-t)} x + e^{-k(T-t)} G(t,T) y, // // and the second exponential is the part worth testing: dropping it // leaves a formula that is still exact at T = t and wrong everywhere // else, which is the easiest kind of error to keep. let kappa = 0.05f64; let sigma_r = |u: f64| 0.01 * (1.0 + 0.3 * (2.0 * u).sin()); let g = |from: f64, to: f64| (1.0 - (-kappa * (to - from)).exp()) / kappa; let (t, steps) = (5.0f64, 20_000); let dt = t / steps as f64; let mut rng = crate::pathwise::Rng::new(20260807); let increments: Vec<f64> = (0..steps).map(|_| dt.sqrt() * rng.next_normal()).collect(); // The states, built from the definitions in the chapter. `a` is the // third accumulator that the algebra cancels; it is needed to get from // the raw stochastic integral to x = r - f(0,t). let (mut stochastic, mut y, mut a) = (0.0, 0.0, 0.0); for (i, dw) in increments.iter().enumerate() { let u = i as f64 * dt; let elapsed = t - u; stochastic += sigma_r(u) * (-kappa * elapsed).exp() * dw; y += sigma_r(u).powi(2) * (-2.0 * kappa * elapsed).exp() * dt; a += sigma_r(u).powi(2) * (-kappa * elapsed).exp() * dt; } let x = stochastic + (a - y) / kappa; for maturity in [5.0, 7.0, 10.0, 20.0, 30.0] { // Straight from the HJM equation: drift plus diffusion, integrated. let direct: f64 = increments .iter() .enumerate() .map(|(i, dw)| { let u = i as f64 * dt; let vol = sigma_r(u) * (-kappa * (maturity - u)).exp(); vol * sigma_r(u) * g(u, maturity) * dt + vol * dw }) .sum(); let ahead = maturity - t; let discount = (-kappa * ahead).exp(); let rebuilt = discount * x + discount * g(t, maturity) * y; assert!( (direct - rebuilt).abs() < 1e-9, "T={maturity}: direct {direct:.12} against rebuilt {rebuilt:.12}" ); // And the version without the second exponential is wrong, by more // than a rounding error, everywhere but at the short end. if maturity > t + 1.0 { let dropped = discount * x + g(t, maturity) * y; assert!( (direct - dropped).abs() > 1e-5, "T={maturity}: dropping the discount on y should be visible" ); } } } #[test] fn markov_alone_does_not_force_an_exponential() { // The counterexample the chapter turns on, and the reason the section // cannot stop at "Markov implies separable". // // Take a volatility that depends on the maturity DATE and not on the // time left to it: sigma(u,T) = psi(T), for a psi with no finite // realisation as a time-to-maturity profile at all. As a kernel in // (u,T) it has rank one, so the model is Markov in a single state --- W // itself --- with psi entirely arbitrary. // // So finite rank is what Markov requires, and finite rank does not // require an exponential. Time homogeneity is the second, separate // demand, and it is the one the exponential comes from. let psi = |t: f64| 1.0 / (1.0 + t); let calendar = kernel_rank(|_u, t| psi(t), TOL); assert_eq!(calendar, 1, "a shape fixed in calendar time is rank one"); let time_to_maturity = realisation_dimension(psi, TOL); assert!( time_to_maturity > 8, "the same shape as a time-to-maturity profile should have no finite \ realisation, got {time_to_maturity}" ); // Not special to this psi: any shape at all is rank one when it is // pinned to the calendar. for shape in [ &(|t: f64| 1.0 / (1.0 + t * t)) as &dyn Fn(f64) -> f64, &(|t: f64| (-0.05 * t.sqrt()).exp()), &(|t: f64| (1.0 + t).ln()), ] { assert_eq!(kernel_rank(|_u, t| shape(t), TOL), 1); } } #[test] fn a_sum_of_products_has_the_rank_of_the_sum() { // The lemma in its own right, away from any shift structure: a kernel // that is a sum of n products has rank n, which is what makes "Markov" // and "separable into n terms" the same statement. let terms: [(&dyn Fn(f64) -> f64, &dyn Fn(f64) -> f64); 3] = [ (&|u: f64| (-0.1 * u).exp(), &|t: f64| 1.0 / (1.0 + t)), (&|u: f64| u, &|t: f64| (0.3 * t).sin()), (&|u: f64| (1.0 + u).ln(), &|t: f64| (-t * t / 50.0).exp()), ]; for n in 1..=3 { let rank = kernel_rank( |u, t| terms[..n].iter().map(|(phi, psi)| phi(u) * psi(t)).sum::<f64>(), TOL, ); assert_eq!(rank, n, "a sum of {n} products"); } } #[test] fn anything_else_needs_infinitely_many() { // The negative half, which is what makes the theorem worth stating. A // perfectly innocent decaying shape that is not a quasi-exponential has // no finite Markov realisation at all; the rank here is limited only by // the tolerance and the grid, not by the function. assert!(realisation_dimension(|x| 1.0 / (1.0 + x), TOL) > 8); assert!(realisation_dimension(|x: f64| 1.0 / (1.0 + x * x), TOL) > 8); assert!(realisation_dimension(|x: f64| (-0.05 * x.sqrt()).exp(), TOL) > 8); } #[test] fn the_states_of_a_quasi_exponential_reconstruct_the_curve_exactly() { // Not a rank count but the splitting identity itself, pathwise. The // Markovian term structure chapter expands g(x) = x^k exp(-kappa x) by // the binomial theorem into // // g(a+b) = sum_j C(k,j) a^{k-j} e^{-kappa a} . b^j e^{-kappa b}, // // so the stochastic integral splits into k+1 maturity-free states with // deterministic multipliers. This is algebra per increment rather than // a limit, so the two sides should agree to machine precision. let kappa = 0.05f64; let sigma_r = |u: f64| 0.01 * (1.0 + 0.3 * (2.0 * u).sin()); let (t, steps) = (5.0f64, 500); let dt = t / steps as f64; let mut rng = crate::pathwise::Rng::new(20260807); let increments: Vec<f64> = (0..steps).map(|_| dt.sqrt() * rng.next_normal()).collect(); let binomial = |k: u32, j: u32| -> f64 { (0..j).map(|i| (k - i) as f64 / (i + 1) as f64).product::<f64>() }; for k in 0..4u32 { let g = |x: f64| x.powi(k as i32) * (-kappa * x).exp(); // The k+1 states, which know nothing about any maturity: the basis // is g_j(b) = b^j exp(-kappa b). let states: Vec<f64> = (0..=k) .map(|j| { increments .iter() .enumerate() .map(|(i, dw)| { let elapsed = t - i as f64 * dt; sigma_r(i as f64 * dt) * elapsed.powi(j as i32) * (-kappa * elapsed).exp() * dw }) .sum::<f64>() }) .collect(); for maturity in [5.0, 7.0, 10.0, 20.0, 30.0] { let ahead = maturity - t; let direct: f64 = increments .iter() .enumerate() .map(|(i, dw)| sigma_r(i as f64 * dt) * g(maturity - i as f64 * dt) * dw) .sum(); // c_j(a) = C(k,j) a^{k-j} exp(-kappa a). let rebuilt: f64 = (0..=k) .map(|j| { binomial(k, j) * ahead.powi((k - j) as i32) * (-kappa * ahead).exp() * states[j as usize] }) .sum(); assert!( (direct - rebuilt).abs() < 1e-13 * direct.abs().max(1e-6), "k={k}, T={maturity}: direct {direct:.18} against rebuilt {rebuilt:.18}" ); } } }} #[cfg(test)]mod time_dependent_kappa_tests { use super::*; /// A mean reversion that varies with calendar time, staying positive. fn kappa(u: f64) -> f64 { 0.05 + 0.03 * (u / 3.0).sin() } /// `integral of kappa from 0 to x`, by Simpson on a fine grid. fn integral(x: f64) -> f64 { const N: usize = 2000; let h = x / N as f64; let mut total = kappa(0.0) + kappa(x); for i in 1..N { let w = if i % 2 == 0 { 2.0 } else { 4.0 }; total += w * kappa(i as f64 * h); } total * h / 3.0 } /// The question the quasi-Gaussian chapter's definition raises: the /// separability condition allows a time-dependent mean reversion, and the /// proposition that forces an exponential assumes a constant one. Both are /// right, and this is why. /// /// As a kernel in `(t,T)` the decay is `G(T)/G(t)`, an outer product, so it /// has rank one whatever `kappa` does. A time-dependent mean reversion /// therefore costs no state variables at all. #[test] fn a_time_dependent_mean_reversion_is_still_rank_one_in_t_and_maturity() { // In (t, T), with the two as independent coordinates rather than as a // date and a time remaining. That distinction is the whole point: the // same decay in (t, T-t) is not an outer product and does not have // rank one, which the next test shows. let sigma = |t: f64, big_t: f64| (-(integral(big_t) - integral(t))).exp(); assert_eq!(kernel_rank(sigma, 1e-9), 1); } /// What it costs instead is time homogeneity. Read as a shape in time to /// maturity, the same decay is not an exponential and its shifts span a /// space of dimension far above one, so no fixed maturity profile carries it. /// /// This is example 12.1 of the chapter in another guise: the same function /// gives one state read against the maturity date and many read against time /// remaining. #[test] fn but_read_as_a_maturity_shape_it_has_no_small_realisation() { let shape = |x: f64| (-integral(x)).exp(); let dimension = realisation_dimension(shape, 1e-9); assert!(dimension > 3, "a varying kappa gave a shape of dimension {dimension}"); // The constant case, for contrast, is exactly one. assert_eq!(realisation_dimension(|x: f64| (-0.05 * x).exp(), 1e-9), 1); }} /// Where the harmonic mean of `sigma_r(u) = 1 + b u + c u^2` over the journey/// from the forward out to `u` is least.////// The short-expiry theorem of the smile dynamics chapter gives the implied/// volatility as the harmonic mean of the local volatility along the path, not/// the arithmetic one, so this is the vertex that theorem predicts. The/// Markovian term structure chapter uses the arithmetic rule and gets/// `-3b/(4c)`; the two agree to first order in `(b, c)` and this measures the/// gap between them against a solved smile.pub fn harmonic_vertex(b: f64, c: f64) -> f64 { let mean = |u: f64| { const N: usize = 4000; let mut reciprocal = 0.0; for i in 0..N { let z = u * (i as f64 + 0.5) / N as f64; reciprocal += 1.0 / (1.0 + b * z + c * z * z); } // u divided by the integral of 1/sigma_r over [0, u]. N as f64 / reciprocal }; let (mut best, mut at) = (f64::MAX, f64::NAN); for i in 1..4000 { let u = 3.0 * i as f64 / 4000.0; let v = mean(u); if v < best { best = v; at = u; } } at} #[cfg(test)]mod averaging_rule_tests { use super::*; /// Which mean the vertex should be read off, and how much it matters. /// /// The smile dynamics chapter's short-expiry theorem gives the harmonic /// mean. Writing `sigma_r = sigma(1 + eps)`, the arithmetic mean is /// `sigma(1 + <eps>)` and the harmonic is `sigma(1 + <eps> - Var eps)`, so /// they part company only at second order in `(b, c)` — and the vertex is a /// first-order statement, being the stationary point of `<eps>`. /// /// Measured against the solved density, the harmonic rule is the better of /// the two everywhere and by very little. What is left over is the finite /// expiry: the theorem is a short-maturity limit and these are five-year /// options, which is a larger error than the choice of mean. #[test] fn the_harmonic_rule_is_better_by_less_than_the_expiry_costs() { let (sigma, expiry) = (0.01f64, 5.0f64); let scale = sigma * expiry.sqrt(); let strikes: Vec<f64> = (0..1601).map(|i| -0.06 + 0.12 * i as f64 / 1600.0).collect(); let mut worst_arithmetic: f64 = 0.0; let mut worst_harmonic: f64 = 0.0; for (b, c) in [(-0.3, 0.2), (-0.3, 0.5), (-0.6, 0.5), (-0.2, 0.4)] { let m = QuasiGaussian { sigma, scale, b, c, ..Default::default() }; let vols = m.smile(expiry, &strikes, 1601); let (mut lowest, mut at) = (f64::MAX, f64::NAN); for (&k, v) in strikes.iter().zip(&vols) { if let Some(v) = v { if *v < lowest { lowest = *v; at = k; } } } let solved = at / scale; worst_arithmetic = worst_arithmetic.max(((-3.0 * b / (4.0 * c)) / solved - 1.0).abs()); worst_harmonic = worst_harmonic.max((harmonic_vertex(b, c) / solved - 1.0).abs()); } assert!( worst_harmonic < worst_arithmetic, "harmonic {worst_harmonic:.4} should beat arithmetic {worst_arithmetic:.4}" ); // Both within a few per cent, and the gap between them smaller still. assert!((worst_arithmetic - 0.044).abs() < 0.005, "arithmetic {worst_arithmetic:.4}"); assert!((worst_harmonic - 0.034).abs() < 0.005, "harmonic {worst_harmonic:.4}"); }}