quant/src/hullwhite.rs
The future-forward basis in the Hull-White model.
//! The future-forward basis in the Hull-White model.//!//! The term structure chapter derives this from scratch and then never draws//! it, which is a pity: the formula is opaque and the picture is not. The basis//! is small at the front end, grows roughly with the square of the maturity,//! and is then pulled back by mean reversion — three statements that take a//! paragraph to make and a glance to see.//!//! This implements the chapter's own expressions, so the figure is the algebra//! evaluated rather than a separate model that happens to agree. /// The exact basis, from the term structure chapter:////// ```text/// Fut(0) - FRA(0) = (1/delta + L) (exp(M + V/2) - 1)/// ```////// with `M` the drift of the yield away from today's forward and `V` its/// variance, both written in terms of `u = 1 - exp(-alpha T)` and/// `v = 1 - exp(-alpha delta)`.////// `forward_rate` is today's simply compounded forward over the period, which/// enters only through the leading factor and barely matters — the basis is/// essentially `M + V/2` divided by `delta`.pub fn basis(t: f64, delta: f64, sigma: f64, alpha: f64, forward_rate: f64) -> f64 { if t <= 0.0 || delta <= 0.0 || sigma < 0.0 { return 0.0; } // The alpha -> 0 limit is a genuine removable singularity: every term below // divides by a power of alpha. The term structure chapter takes the limit // by hand and gets the Ho-Lee expressions, which is what this returns. if alpha.abs() < 1e-8 { let m = 0.5 * sigma * sigma * t * delta * (t + delta); let v = sigma * sigma * t * delta * delta; return (1.0 / delta + forward_rate) * ((m + 0.5 * v).exp() - 1.0); } let u = 1.0 - (-alpha * t).exp(); let v = 1.0 - (-alpha * delta).exp(); let a3 = alpha * alpha * alpha; let m = sigma * sigma / (4.0 * a3) * u * v * (2.0 * u + 2.0 * v - u * v); let var = sigma * sigma / (2.0 * a3) * u * (2.0 - u) * v * v; (1.0 / delta + forward_rate) * ((m + 0.5 * var).exp() - 1.0)} /// The Ho-Lee approximation the term structure chapter ends on,////// ```text/// Fut - FRA ~ sigma^2 T (T + 2 delta) / 2,/// ```////// which is the no-mean-reversion limit with the exponential expanded away. It/// is drawn beside the exact curve so that the effect of mean reversion is the/// visible gap between them rather than an assertion.pub fn ho_lee_basis(t: f64, delta: f64, sigma: f64) -> f64 { 0.5 * sigma * sigma * t * (t + 2.0 * delta)} #[cfg(test)]mod tests { use super::*; const DELTA: f64 = 0.25; const SIGMA: f64 = 0.01; const F: f64 = 0.03; #[test] fn the_basis_is_positive_and_grows_with_maturity() { // The term structure chapter's sign argument: the covariance between // the fixing and the discount factor is negative, so the futures rate // sits above the forward. Nothing about the parameters should be able // to flip that. let mut previous = 0.0; for t in [0.5, 1.0, 2.0, 5.0, 10.0] { let b = basis(t, DELTA, SIGMA, 0.03, F); assert!(b > 0.0, "basis at {t}y was {b}"); assert!(b > previous, "basis fell from {previous} to {b} at {t}y"); previous = b; } } #[test] fn it_grows_like_t_times_t_plus_two_delta() { // Not like T squared, which is the loose thing one says about it. The // exact leading behaviour is T(T + 2 delta), and the difference is not // academic at the front end: doubling the maturity from one year to two // multiplies the basis by 3.33, not by 4, because the accrual period is // still a quarter of the maturity. Only well beyond the accrual does the // growth become the quadratic it is usually described as. let scaling = |t: f64| { basis(2.0 * t, DELTA, SIGMA, 1e-9, F) / basis(t, DELTA, SIGMA, 1e-9, F) }; let predicted = |t: f64| 2.0 * t * (2.0 * t + 2.0 * DELTA) / (t * (t + 2.0 * DELTA)); for t in [1.0, 2.0, 5.0, 20.0] { let (got, want) = (scaling(t), predicted(t)); assert!((got - want).abs() < 0.05, "at {t}y: scaled by {got}, expected {want}"); } // And it does approach four, from below. assert!(scaling(1.0) < 3.5 && scaling(50.0) > 3.9); } #[test] fn mean_reversion_pulls_the_basis_down() { // And increasingly so with maturity, since it is the accumulated decay // that matters. for t in [2.0, 5.0, 10.0] { let none = basis(t, DELTA, SIGMA, 0.0, F); let some = basis(t, DELTA, SIGMA, 0.05, F); assert!(some < none, "at {t}y mean reversion did not reduce the basis"); } let short = basis(2.0, DELTA, SIGMA, 0.0, F) / basis(2.0, DELTA, SIGMA, 0.05, F); let long = basis(10.0, DELTA, SIGMA, 0.0, F) / basis(10.0, DELTA, SIGMA, 0.05, F); assert!(long > short, "mean reversion should bite harder further out"); } #[test] fn the_zero_reversion_limit_matches_ho_lee() { // The term structure chapter takes alpha to zero by hand and gets the // textbook formula. The exact expression evaluated at a tiny alpha must // agree with it, or one of the two derivations is wrong. for t in [1.0, 3.0, 7.0] { let exact = basis(t, DELTA, SIGMA, 1e-9, F); let approx = ho_lee_basis(t, DELTA, SIGMA); let relative = (exact - approx).abs() / approx; assert!(relative < 0.02, "at {t}y: exact {exact} against Ho-Lee {approx}"); } } #[test] fn the_alpha_limit_is_continuous() { // The formula divides by alpha cubed, so the limit is taken separately. // The two branches have to meet, or the figure has a step in it wherever // the slider passes zero. let near = basis(5.0, DELTA, SIGMA, 1e-7, F); let at = basis(5.0, DELTA, SIGMA, 0.0, F); assert!((near - at).abs() / at < 1e-4, "{near} against {at}"); }} /// Jamshidian's decomposition of an option on a coupon bond.////// In a one-factor model every discount bond is a decreasing function of the one/// state, `P_i(x) = exp(-a_i - b_i x)` with `b_i > 0`, so a coupon bond is/// decreasing too and meets any strike at a single critical state. This finds/// that state by bisection.////// The term structure chapter uses it to price a swaption in closed form: at the/// critical state each bond has a value `P_i(x*)`, and the option on the sum/// becomes a sum of options on the parts struck there.pub fn critical_state(coupons: &[f64], a: &[f64], b: &[f64], strike: f64) -> f64 { let bond = |x: f64| -> f64 { coupons .iter() .zip(a) .zip(b) .map(|((c, ai), bi)| c * (-ai - bi * x).exp()) .sum() }; let (mut lo, mut hi) = (-5.0f64, 5.0f64); for _ in 0..200 { let mid = 0.5 * (lo + hi); // Decreasing in x, so the root moves the opposite way to a bisection on // an increasing function. if bond(mid) > strike { lo = mid; } else { hi = mid; } } 0.5 * (lo + hi)} #[cfg(test)]mod jamshidian_tests { use super::*; /// A five-year annual coupon bond, in exponential-affine form. fn bond_data() -> (Vec<f64>, Vec<f64>, Vec<f64>) { let coupons = vec![0.04, 0.04, 0.04, 0.04, 1.04]; let a: Vec<f64> = (1..=5).map(|i| 0.03 * i as f64).collect(); // Loadings rise with maturity, as they do in any one-factor model. let b: Vec<f64> = (1..=5).map(|i| i as f64 * 0.9).collect(); (coupons, a, b) } /// The identity the decomposition rests on. Because every bond moves the /// same way with the one state, the option on the coupon bond and the sum of /// options on the zeros agree at *every* state, not merely in expectation. #[test] fn the_option_on_the_sum_is_the_sum_of_the_options() { let (c, a, b) = bond_data(); let strike = 1.0; let x_star = critical_state(&c, &a, &b, strike); for i in 0..=200 { let x = -2.0 + 4.0 * i as f64 / 200.0; let coupon_bond: f64 = c.iter().zip(&a).zip(&b).map(|((c, ai), bi)| c * (-ai - bi * x).exp()).sum(); let direct = (coupon_bond - strike).max(0.0); let decomposed: f64 = c .iter() .zip(&a) .zip(&b) .map(|((c, ai), bi)| { let strike_i = (-ai - bi * x_star).exp(); c * ((-ai - bi * x).exp() - strike_i).max(0.0) }) .sum(); assert!( (direct - decomposed).abs() < 1e-12, "state {x:.3}: direct {direct:.12} against decomposed {decomposed:.12}" ); } } /// And the hypothesis is one factor, not tractability. /// /// With two states the exercise boundary is a curve rather than a point, so /// there is no single critical state at which to strike the pieces. Two /// states on that boundary give the same coupon bond value and different /// individual bond prices, which is exactly the obstruction: whichever of /// them supplied the strikes, the other would be decomposed wrongly. #[test] fn two_factors_have_no_single_critical_state() { let (c, a, b) = bond_data(); // A second factor loading differently across maturities: short end more, // long end less, which is the slope shape of the curve chapters. let d: Vec<f64> = (1..=5).map(|i| 1.2 - 0.2 * i as f64).collect(); let coupon_bond = |x: f64, y: f64| -> f64 { c.iter() .zip(&a) .zip(b.iter().zip(&d)) .map(|((c, ai), (bi, di))| c * (-ai - bi * x - di * y).exp()) .sum() }; // Two points on the same exercise boundary, found by solving in x for // two different y. let solve_x = |y: f64| { let (mut lo, mut hi) = (-5.0f64, 5.0f64); for _ in 0..200 { let mid = 0.5 * (lo + hi); if coupon_bond(mid, y) > 1.0 { lo = mid; } else { hi = mid; } } 0.5 * (lo + hi) }; let (y1, y2) = (-0.3, 0.3); let (x1, x2) = (solve_x(y1), solve_x(y2)); assert!((coupon_bond(x1, y1) - 1.0).abs() < 1e-9); assert!((coupon_bond(x2, y2) - 1.0).abs() < 1e-9); // Same coupon bond value, different constituent bonds. let piece = |x: f64, y: f64, i: usize| (-a[i] - b[i] * x - d[i] * y).exp(); let gap: f64 = (0..5) .map(|i| (piece(x1, y1, i) - piece(x2, y2, i)).abs()) .fold(0.0f64, f64::max); assert!( gap > 0.05, "the two boundary states gave nearly identical bonds, gap {gap:.4}" ); }}