quant/src/hjm.rs
The Heath-Jarrow-Morton drift condition, simulated.
//! The Heath-Jarrow-Morton drift condition, simulated.//!//! The term structure chapter derives that an arbitrage-free evolution of the//! forward curve must have//!//! ```text//! df(t,T) = sigma(t,T) (integral_t^T sigma(t,s) ds) dt + sigma(t,T) dW,//! ```//!//! with the drift determined by the volatility and no freedom left. That is a//! strong claim and it has a sharp test: the drift is exactly what makes the//! model reprice the curve it started from. Simulate with it and//! `E[exp(-integral r)]` returns today's discount factor; simulate without it//! and the model has already lost the market. use crate::pathwise::Rng; /// A Gaussian Heath-Jarrow-Morton model with an exponentially decaying/// volatility, `sigma(t,T) = sigma exp(-kappa (T-t))`.////// The maturity grid and the time grid share a spacing, so that the short rate/// `r(t) = f(t,t)` is always a point of the curve being carried rather than an/// interpolation of it.pub struct GaussianHjm { pub sigma: f64, pub kappa: f64, /// The initial instantaneous forward curve, as a function of maturity. pub initial: fn(f64) -> f64,} /// Which drift to simulate with.#[derive(Clone, Copy, PartialEq)]pub enum Drift { /// The Heath-Jarrow-Morton drift. NoArbitrage, /// No drift at all, as one might naively write down. Zero,} impl GaussianHjm { fn vol(&self, ahead: f64) -> f64 { self.sigma * (-self.kappa * ahead).exp() } /// `integral_t^T sigma(t,s) ds`, in closed form for this volatility. fn integrated_vol(&self, ahead: f64) -> f64 { self.sigma * (1.0 - (-self.kappa * ahead).exp()) / self.kappa } /// Today's discount factor, by integrating the initial curve. pub fn initial_discount(&self, maturity: f64, steps: usize) -> f64 { let dt = maturity / steps as f64; let integral: f64 = (0..steps).map(|i| (self.initial)((i as f64 + 0.5) * dt) * dt).sum(); (-integral).exp() } /// `E[exp(-integral_0^horizon r(t) dt)]` by simulating the whole curve. /// /// Returns the mean and its standard error. The curve is carried on a /// maturity grid out to `horizon`, every surviving point of it stepped every /// step, so this simulates a *surface* rather than a path --- which is the /// term structure chapter's point about what the state of the model is, and /// also why it costs what it does. /// /// Antithetic: each set of increments is used twice, once negated. The /// discount factor is close to a decreasing function of the accumulated /// noise, so the two are strongly negatively correlated and the pairing /// removes most of the variance for no extra curve stepping --- which is /// what makes this affordable in a test suite the build runs. pub fn simulated_discount( &self, horizon: f64, steps: usize, pairs: usize, drift: Drift, seed: u64, ) -> (f64, f64) { let dt = horizon / steps as f64; let root_dt = dt.sqrt(); let mut rng = Rng::new(seed); // Precompute the maturity-dependent coefficients once: they depend on // the time left, which is the same at every step for a given offset. let vol: Vec<f64> = (0..=steps).map(|k| self.vol(k as f64 * dt)).collect(); let mu: Vec<f64> = (0..=steps) .map(|k| match drift { Drift::NoArbitrage => { let ahead = k as f64 * dt; self.vol(ahead) * self.integrated_vol(ahead) * dt } Drift::Zero => 0.0, }) .collect(); let initial: Vec<f64> = (0..=steps).map(|j| (self.initial)(j as f64 * dt)).collect(); let (mut total, mut total_sq) = (0.0, 0.0); for _ in 0..pairs { let increments: Vec<f64> = (0..steps).map(|_| root_dt * rng.next_normal()).collect(); for sign in [1.0, -1.0] { let mut curve = initial.clone(); let mut accumulated = 0.0; for step in 0..steps { accumulated += curve[step] * dt; let dw = sign * increments[step]; for j in (step + 1)..=steps { let ahead = j - step; curve[j] += mu[ahead] + vol[ahead] * dw; } } let discount = (-accumulated).exp(); total += discount; total_sq += discount * discount; } } let n = 2.0 * pairs as f64; let mean = total / n; let variance = (total_sq / n - mean * mean).max(0.0); // The antithetic pairs are not independent, so the naive standard error // overstates the precision of the pair mean and understates it here. // Reported conservatively from the pooled sample. (mean, (variance / n).sqrt()) }} #[cfg(test)]mod tests { use super::*; fn model() -> GaussianHjm { GaussianHjm { sigma: 0.01, kappa: 0.15, // A gently upward sloping curve, so the test is not about a flat one. initial: |t| 0.02 + 0.01 * (1.0 - (-0.3 * t).exp()), } } /// Antithetic pairs. Twenty thousand is enough because the pairing removes /// most of the variance; the residual error below is dominated by the time /// discretisation rather than by sampling. const PAIRS: usize = 20_000; const SEED: u64 = 20260808; /// Fifty steps a year, which is what the errors below are quoted at. fn steps_for(horizon: f64) -> usize { (horizon * 20.0) as usize } #[test] fn the_drift_condition_reprices_the_initial_curve() { // The claim, tested where it bites. Nothing in the simulation is told // what the initial discount factor is; the drift condition is what makes // the simulated expectation come back to it. let m = model(); for horizon in [1.0, 5.0, 10.0] { let (mean, _) = m.simulated_discount(horizon, steps_for(horizon), PAIRS, Drift::NoArbitrage, SEED); let exact = m.initial_discount(horizon, 4000); let relative = (mean - exact).abs() / exact; assert!( relative < 4e-4, "horizon={horizon}: simulated {mean:.8} against curve {exact:.8}, relative error {relative:.2e}" ); } } #[test] fn without_the_drift_the_model_drifts_away_from_the_market() { // The converse, which is what makes the condition a condition. Drop the // drift -- the obvious thing to write down, since a forward rate looks // as though it ought to be a martingale -- and the model no longer // reprices the curve it was handed. // // The error is second order in the volatility, so it is small at short // horizons and compounds. That is the awkward case rather than a // comfortable one: a wrong model looking nearly right on the short // instruments a desk would check it against. let m = model(); let mut previous = 0.0; let mut worst = 0.0; for horizon in [1.0, 5.0, 10.0] { let (mean, _) = m.simulated_discount(horizon, steps_for(horizon), PAIRS, Drift::Zero, SEED); let exact = m.initial_discount(horizon, 4000); let error = (mean - exact).abs(); assert!(error > previous, "the error should compound with the horizon"); previous = error; worst = error / exact; // One-signed: without the convexity drift the forwards are too low, // so the discount factors are too high. assert!(mean > exact, "horizon={horizon}: the bias should be upwards"); } // By ten years it is an order of magnitude past the discretisation error // the previous test tolerates, so this is not a grid artefact. assert!(worst > 4e-3, "at ten years the relative error is only {worst:.2e}"); } #[test] fn the_drift_is_the_derivative_of_a_squared_volatility() { // The algebraic identity the proof turns on, checked directly: // // sigma(t,T) integral_t^T sigma(t,s) ds = d/dT [ (integral)^2 / 2 ]. // // This is why the drift is not an extra assumption but a consequence -- // it is what differentiating the bond volatility's square produces, and // the bond volatility is pinned by P(T,T) = 1. let m = model(); let h = 1e-6; for ahead in [0.25, 1.0, 4.0, 9.0] { let square = |a: f64| 0.5 * m.integrated_vol(a).powi(2); let differenced = (square(ahead + h) - square(ahead - h)) / (2.0 * h); let claimed = m.vol(ahead) * m.integrated_vol(ahead); assert!( (differenced - claimed).abs() < 1e-6 * claimed.abs().max(1e-12), "ahead={ahead}: differenced {differenced} against claimed {claimed}" ); } }}