quant/src/curve.rs
Discount curves: bootstrapping one from par quotes, and the interpolation schemes that make the answer between the quotes ambiguous.
//! Discount curves: bootstrapping one from par quotes, and the interpolation//! schemes that make the answer between the quotes ambiguous.//!//! The curve construction chapter's claim is that interpolation is a modelling//! choice rather than a numerical detail, and that the evidence for it is//! invisible in the coordinate one is tempted to plot. This module is what//! draws the evidence: the schemes reprice every input instrument identically//! and agree at every node, and their instantaneous forward curves — which is//! what a product actually fixes on — do not agree anywhere else.//!//! Three of the four are one-pass and local. The fourth, Hagan and West's//! monotone convex scheme, is the standard repair for the overshoot the others//! suffer, and it pays for that with non-locality: it has to be solved globally,//! and its shape depends on the data in a way no fixed matrix captures. /// What is interpolated between the curve's nodes.////// Each is linear in *something*, and the something is the whole difference. The/// instantaneous forward rate is a derivative of the log discount factor, so/// whichever quantity is made piecewise linear, the forwards come out one degree/// rougher.#[derive(Clone, Copy, PartialEq, Eq, Debug)]pub enum Interp { /// Linear in the continuously compounded zero rate. The most commonly drawn /// curve and the worst behaved: forwards jump at every node. LinearZero, /// Linear in `ln P(0,t)`, equivalently piecewise constant forwards. LinearLogDf, /// Piecewise linear in the instantaneous forward rate itself, which is the /// only one of the three whose forwards are continuous. LinearForward, /// Hagan and West's monotone convex scheme: continuous forwards that /// reprice exactly and are not allowed to overshoot the range set by /// neighbouring quotes. /// /// The other three are local --- a node's value is determined by its own /// bucket and the one before. This one is not: the forward *at* a node is /// built from the buckets on both sides, so a quote reaches one node further /// back than it does under any of the others, and the curve has to be solved /// globally rather than swept once from the front. MonotoneConvex,} /// One bucket of the monotone convex construction, in the local coordinate/// `x = (t - start) / (end - start)`.////// Everything is expressed as `g(x) = f(t) - discrete`, the deviation of the/// instantaneous forward from the bucket's own average. The conserving condition/// is then simply that `g` integrates to zero across the bucket, which is what/// makes the scheme reprice by construction rather than by iteration.#[derive(Clone, Copy, Debug)]struct Bucket { start: f64, end: f64, discrete: f64, g0: f64, g1: f64, region: Region,} /// Which of Hagan and West's four cases a bucket falls in.////// The plain quadratic through the two node values is used where it behaves; the/// other three replace part of it with a flat section so the curve cannot leave/// the range its neighbours set. `eta` is where the pieces meet.#[derive(Clone, Copy, Debug, PartialEq)]enum Region { /// The unmodified quadratic. Quadratic, /// Flat at `g0`, then a rising quadratic. FlatThenCurve { eta: f64 }, /// A quadratic, then flat at `g1`. CurveThenFlat { eta: f64 }, /// Two quadratics meeting at a common level `a`, used when both ends sit on /// the same side of the average and something in between has to compensate. Bowl { eta: f64, a: f64 },} impl Bucket { /// `g(x)`, the deviation of the forward from the bucket average. fn g(&self, x: f64) -> f64 { let (g0, g1) = (self.g0, self.g1); match self.region { Region::Quadratic => g0 * (1.0 - 4.0 * x + 3.0 * x * x) + g1 * (-2.0 * x + 3.0 * x * x), Region::FlatThenCurve { eta } => { if x <= eta { g0 } else { let u = (x - eta) / (1.0 - eta); g0 + (g1 - g0) * u * u } } Region::CurveThenFlat { eta } => { if x >= eta { g1 } else { let u = (eta - x) / eta; g1 + (g0 - g1) * u * u } } Region::Bowl { eta, a } => { if x <= eta { let u = (eta - x) / eta; a + (g0 - a) * u * u } else { let u = (x - eta) / (1.0 - eta); a + (g1 - a) * u * u } } } } /// `integral of g from 0 to x`, in closed form. /// /// Analytic rather than numerical because the scheme's whole claim is that it /// reprices *exactly*; quadrature here would leave an error precisely where /// the selling point is. fn integral(&self, x: f64) -> f64 { let (g0, g1) = (self.g0, self.g1); match self.region { Region::Quadratic => { g0 * (x - 2.0 * x * x + x * x * x) + g1 * (-(x * x) + x * x * x) } Region::FlatThenCurve { eta } => { if x <= eta { g0 * x } else { let d = x - eta; g0 * x + (g1 - g0) * d * d * d / (3.0 * (1.0 - eta) * (1.0 - eta)) } } Region::CurveThenFlat { eta } => { let head = |x: f64| { let u = (eta - x) / eta; g1 * x + (g0 - g1) * eta / 3.0 * (1.0 - u * u * u) }; if x <= eta { head(x) } else { head(eta) + g1 * (x - eta) } } Region::Bowl { eta, a } => { let head = |x: f64| { let u = (eta - x) / eta; a * x + (g0 - a) * eta / 3.0 * (1.0 - u * u * u) }; if x <= eta { head(x) } else { let d = x - eta; head(eta) + a * d + (g1 - a) * d * d * d / (3.0 * (1.0 - eta) * (1.0 - eta)) } } } }} impl Curve { /// The buckets of the monotone convex construction, rebuilt from the yields. /// /// Derived rather than stored, because the node forwards depend on the /// buckets on *both* sides of each node: a stored copy would go stale the /// moment the bootstrap moved a later node, and the staleness would show up /// as a curve that no longer repriced. fn buckets(&self) -> Vec<Bucket> { let n = self.times.len(); if n == 0 { return Vec::new(); } // Bucket boundaries, with the origin prepended. let mut edges = Vec::with_capacity(n + 1); edges.push(0.0); edges.extend_from_slice(&self.times); // The average forward over each bucket. This is what the quotes pin -- // not the forward at any single date -- and it is the quantity the // scheme is built to conserve. let discrete: Vec<f64> = (0..n) .map(|i| { let previous = if i == 0 { 0.0 } else { self.yields[i - 1] }; (self.yields[i] - previous) / (edges[i + 1] - edges[i]) }) .collect(); // Node forwards: an interval-weighted average of the two buckets meeting // there, then collared so the curve cannot go negative. let mut node = vec![0.0; n + 1]; for j in 1..n { let (before, here, after) = (edges[j - 1], edges[j], edges[j + 1]); let span = after - before; node[j] = (here - before) / span * discrete[j] + (after - here) / span * discrete[j - 1]; } node[0] = discrete[0] - 0.5 * (node[1.min(n)] - discrete[0]); node[n] = discrete[n - 1] - 0.5 * (node[n - 1] - discrete[n - 1]); let collar = |v: f64, limit: f64| v.max(0.0).min(2.0 * limit.max(0.0)); node[0] = collar(node[0], discrete[0]); node[n] = collar(node[n], discrete[n - 1]); for j in 1..n { node[j] = collar(node[j], discrete[j - 1].min(discrete[j])); } (0..n) .map(|i| { let (g0, g1) = (node[i] - discrete[i], node[i + 1] - discrete[i]); Bucket { start: edges[i], end: edges[i + 1], discrete: discrete[i], g0, g1, region: classify(g0, g1), } }) .collect() }} /// Classify a bucket into one of Hagan and West's four cases.////// The boundaries are theirs. The cases exist so that the interpolant never/// leaves the interval set by the neighbouring inputs: wherever the plain/// quadratic would overshoot, part of it is flattened instead, and the flat part/// is sized so the area under the curve is unchanged.fn classify(g0: f64, g1: f64) -> Region { const TINY: f64 = 1e-12; // Both ends on the same side of the average: nothing can be monotone here, // so the curve dips to a common level in between and comes back. if (g0 >= 0.0 && g1 >= 0.0) || (g0 <= 0.0 && g1 <= 0.0) { let sum = g0 + g1; if sum.abs() < TINY { return Region::Quadratic; } let eta = g1 / sum; if !(eta > TINY && eta < 1.0 - TINY) { return Region::Quadratic; } return Region::Bowl { eta, a: -g0 * g1 / sum }; } // Opposite signs from here on, so g1 - g0 cannot vanish. let span = g1 - g0; if span.abs() < TINY { return Region::Quadratic; } let steep = if g0 < 0.0 { g1 > -2.0 * g0 } else { g1 < -2.0 * g0 }; if steep { let eta = (g1 + 2.0 * g0) / span; if eta > TINY && eta < 1.0 - TINY { return Region::FlatThenCurve { eta }; } return Region::Quadratic; } let shallow = if g0 > 0.0 { g1 > -0.5 * g0 && g1 < 0.0 } else { g1 < -0.5 * g0 && g1 > 0.0 }; if shallow { let eta = 3.0 * g1 / span; if eta > TINY && eta < 1.0 - TINY { return Region::CurveThenFlat { eta }; } } Region::Quadratic} /// A discount curve, held as node times and the quantity being interpolated.#[derive(Clone, Debug)]pub struct Curve { /// Node maturities in years, strictly increasing, none at zero. times: Vec<f64>, /// `-ln P(0,t_i)`, the cumulative integral of the instantaneous forwards. /// Every scheme stores this; they differ in how they read between the nodes. yields: Vec<f64>, /// Node values of the instantaneous forward rate. Only [`Interp::LinearForward`] /// uses them, and only that scheme has them as independent unknowns. forwards: Vec<f64>, pub interp: Interp,} impl Curve { /// `-ln P(0,t)`: the integral of the instantaneous forwards out to `t`. fn integrated(&self, t: f64) -> f64 { if t <= 0.0 { return 0.0; } let n = self.times.len(); if n == 0 { return 0.0; } // Index of the first node at or after t. let i = self.times.partition_point(|&x| x < t); match self.interp { Interp::LinearZero => { // z(t) linear in t, so the integral is t * z(t). let z = self.linear_zero(t); t * z } Interp::LinearLogDf => { if i == 0 { // Before the first node, extrapolate flat from the origin. return t / self.times[0] * self.yields[0]; } if i >= n { // Beyond the last node, continue at the last forward rate. let last = n - 1; let slope = if last == 0 { self.yields[0] / self.times[0] } else { (self.yields[last] - self.yields[last - 1]) / (self.times[last] - self.times[last - 1]) }; return self.yields[last] + slope * (t - self.times[last]); } let (t0, t1) = (self.times[i - 1], self.times[i]); let (y0, y1) = (self.yields[i - 1], self.yields[i]); y0 + (y1 - y0) * (t - t0) / (t1 - t0) } Interp::LinearForward => { // f piecewise linear, so its integral is piecewise quadratic. if i == 0 { // Flat at the first node's forward from the origin. return self.forwards[0] * t; } if i >= n { let last = n - 1; return self.yields[last] + self.forwards[last] * (t - self.times[last]); } let (t0, t1) = (self.times[i - 1], self.times[i]); let (f0, f1) = (self.forwards[i - 1], self.forwards[i]); let h = t - t0; let slope = (f1 - f0) / (t1 - t0); self.yields[i - 1] + f0 * h + 0.5 * slope * h * h } Interp::MonotoneConvex => { let buckets = self.buckets(); if buckets.is_empty() { return 0.0; } if i >= n { // Beyond the last node, flat at the forward the curve // actually reaches there -- not at the bucket's average, // which is a different number and would put a step in the // forwards at the very last node. let last = &buckets[buckets.len() - 1]; let terminal = last.discrete + last.g1; return self.yields[n - 1] + terminal * (t - self.times[n - 1]); } let b = &buckets[i]; let width = b.end - b.start; let x = (t - b.start) / width; let before = if i == 0 { 0.0 } else { self.yields[i - 1] }; // Conservation makes the whole-bucket integral exactly the // quoted average, so partial buckets are the only place the // shape enters at all. before + width * (b.discrete * x + b.integral(x)) } } } fn linear_zero(&self, t: f64) -> f64 { let n = self.times.len(); let i = self.times.partition_point(|&x| x < t); if i == 0 { return self.yields[0] / self.times[0]; } if i >= n { return self.yields[n - 1] / self.times[n - 1]; } let (t0, t1) = (self.times[i - 1], self.times[i]); let (z0, z1) = (self.yields[i - 1] / t0, self.yields[i] / t1); z0 + (z1 - z0) * (t - t0) / (t1 - t0) } /// The discount factor `P(0,t)`. pub fn df(&self, t: f64) -> f64 { (-self.integrated(t)).exp() } /// The continuously compounded zero rate. pub fn zero(&self, t: f64) -> f64 { if t <= 0.0 { return self.inst_forward(0.0); } self.integrated(t) / t } /// The instantaneous forward rate `f(0,t)`. /// /// Computed as the exact derivative of whatever is being interpolated rather /// than by differencing the discount factors — the point of the figure is /// that these are genuinely discontinuous under some schemes, and a finite /// difference would smear the discontinuity into a steep slope and hide it. pub fn inst_forward(&self, t: f64) -> f64 { let n = self.times.len(); let i = self.times.partition_point(|&x| x < t); match self.interp { Interp::LinearZero => { // f = z + t z', and z' jumps at each node while t does not. let z = self.linear_zero(t); let zp = if i == 0 || i >= n { 0.0 } else { let (t0, t1) = (self.times[i - 1], self.times[i]); (self.yields[i] / t1 - self.yields[i - 1] / t0) / (t1 - t0) }; z + t * zp } Interp::LinearLogDf => { if i == 0 { return self.yields[0] / self.times[0]; } if i >= n { let last = n - 1; return if last == 0 { self.yields[0] / self.times[0] } else { (self.yields[last] - self.yields[last - 1]) / (self.times[last] - self.times[last - 1]) }; } (self.yields[i] - self.yields[i - 1]) / (self.times[i] - self.times[i - 1]) } Interp::LinearForward => { if i == 0 { return self.forwards[0]; } if i >= n { return self.forwards[n - 1]; } let (t0, t1) = (self.times[i - 1], self.times[i]); let (f0, f1) = (self.forwards[i - 1], self.forwards[i]); f0 + (f1 - f0) * (t - t0) / (t1 - t0) } Interp::MonotoneConvex => { let buckets = self.buckets(); if buckets.is_empty() { return 0.0; } if i >= n { let last = &buckets[buckets.len() - 1]; return last.discrete + last.g1; } let b = &buckets[i]; let x = (t - b.start) / (b.end - b.start); b.discrete + b.g(x.clamp(0.0, 1.0)) } } }} /// Value of a par instrument's fixed leg plus redemption, given a curve.////// A bond quoted at par with coupon `rate` paid `freq` times a year prices to 1,/// which is the equation the bootstrap solves at each node.fn par_bond_price(curve: &Curve, maturity: f64, rate: f64, freq: f64) -> f64 { let n = (maturity * freq).round().max(1.0) as usize; let accrual = 1.0 / freq; let mut pv = 0.0; for k in 1..=n { let t = maturity - (n - k) as f64 * accrual; pv += rate * accrual * curve.df(t); } pv + curve.df(maturity)} /// Bootstrap a curve from par yields.////// `tenors` are maturities in years, strictly increasing; `par_rates` the quoted/// par yields as decimals; `freq` the coupon frequency of the quoted instrument./// Tenors shorter than one coupon period are treated as a single payment at/// maturity, which is what a deposit or a bill is.////// The scheme matters here as well as afterwards: the coupon dates of a ten year/// instrument fall between the nodes, so the value of a node depends on how the/// curve is read between the earlier ones. This is why the three schemes give/// three different curves rather than three readings of one curve.pub fn bootstrap_par( tenors: &[f64], par_rates: &[f64], freq: f64, interp: Interp,) -> Curve { let mut curve = single_pass(tenors, par_rates, freq, interp); if interp == Interp::MonotoneConvex { // A single sweep is not enough here, and the reason is the scheme's // defining property rather than a shortcoming of the sweep. The forward // at a node is built from the buckets on both sides of it, so while node // k was being solved the bucket beyond it did not yet exist, and the // shape assumed for it was wrong. Re-solving every node against the // finished curve and repeating converges quickly, because the // dependence on the far side is weak. // // This is the non-locality the chapter warns about, arriving as a // concrete cost: the other three schemes are done in one pass. for _ in 0..100 { let before = curve.yields.clone(); for k in 0..curve.times.len() { resolve_node(&mut curve, k, tenors[k], par_rates[k], freq); } let moved = curve .yields .iter() .zip(&before) .fold(0.0f64, |worst, (a, b)| worst.max((a - b).abs())); if moved < 1e-15 { break; } } } curve} /// One sweep of the bootstrap: solve each node in turn, left to right.fn single_pass(tenors: &[f64], par_rates: &[f64], freq: f64, interp: Interp) -> Curve { assert_eq!(tenors.len(), par_rates.len()); let mut curve = Curve { times: Vec::new(), yields: Vec::new(), forwards: Vec::new(), interp, }; for (k, (&t, &rate)) in tenors.iter().zip(par_rates).enumerate() { curve.times.push(t); curve.yields.push(rate * t); curve.forwards.push(rate); if t <= 1.0 / freq + 1e-12 { let df = 1.0 / (1.0 + rate * t); set_node(&mut curve, k, -df.ln()); continue; } let (mut lo, mut hi) = (-0.5 * t, 2.0 * t); for _ in 0..200 { let mid = 0.5 * (lo + hi); set_node(&mut curve, k, mid); if par_bond_price(&curve, t, rate, freq) > 1.0 { lo = mid; } else { hi = mid; } if hi - lo < 1e-14 { break; } } set_node(&mut curve, k, 0.5 * (lo + hi)); } curve} /// The same bootstrap, keeping a snapshot of the curve at every stage.////// The linear products chapter describes the construction as a sequence and then/// draws only its result, which leaves the reader to imagine the sequence. This/// returns it: one curve per instrument as it is added, and then one per sweep of/// the monotone convex re-solve.////// Both halves are worth seeing. Adding an instrument extends the curve to a new/// maturity and, because a par bond pays coupons between the nodes, moves the/// part already built. The re-solve afterwards is the non-locality of the scheme/// arriving as a cost: a node's forward is built from the buckets on both sides/// of it, so while it was being solved the bucket beyond it did not exist, and/// the whole curve has to be swept again until it stops moving.pub fn bootstrap_par_stages( tenors: &[f64], par_rates: &[f64], freq: f64, interp: Interp,) -> Vec<Curve> { let mut stages = Vec::new(); // One stage per instrument: bootstrap the first k quotes and keep the result. for k in 1..=tenors.len() { stages.push(bootstrap_par(&tenors[..k], &par_rates[..k], freq, interp)); } // Then, for the scheme that needs them, the sweeps. Rebuild the single pass // and re-solve it one sweep at a time so each is a stage of its own. if interp == Interp::MonotoneConvex { let mut curve = single_pass(tenors, par_rates, freq, interp); stages.push(curve.clone()); for _ in 0..8 { let before = curve.yields.clone(); for k in 0..curve.times.len() { resolve_node(&mut curve, k, tenors[k], par_rates[k], freq); } stages.push(curve.clone()); let moved = curve .yields .iter() .zip(&before) .fold(0.0f64, |worst, (a, b)| worst.max((a - b).abs())); if moved < 1e-12 { break; } } } stages} /// Re-solve one node's yield against the rest of the curve as it now stands.fn resolve_node(curve: &mut Curve, k: usize, t: f64, rate: f64, freq: f64) { if t <= 1.0 / freq + 1e-12 { let df = 1.0 / (1.0 + rate * t); set_node(curve, k, -df.ln()); return; } let (mut lo, mut hi) = (-0.5 * t, 2.0 * t); for _ in 0..200 { let mid = 0.5 * (lo + hi); set_node(curve, k, mid); if par_bond_price(curve, t, rate, freq) > 1.0 { lo = mid; } else { hi = mid; } if hi - lo < 1e-14 { break; } } set_node(curve, k, 0.5 * (lo + hi));} /// Set node `k`'s value to the integrated yield `y`, keeping the scheme's own/// parametrisation consistent.////// For the forward-interpolating scheme the unknown is the node's forward rate,/// not its yield, so the two are converted through the piecewise quadratic that/// scheme integrates to.fn set_node(curve: &mut Curve, k: usize, y: f64) { curve.yields[k] = y; if curve.interp == Interp::LinearForward { curve.forwards[k] = if k == 0 { // f linear from f_0 on [0, t_0] with the same value at both ends. y / curve.times[0] } else { // y_k - y_{k-1} = (f_{k-1} + f_k)/2 * (t_k - t_{k-1}). let dt = curve.times[k] - curve.times[k - 1]; 2.0 * (y - curve.yields[k - 1]) / dt - curve.forwards[k - 1] }; }} #[cfg(test)]mod tests { use super::*; const TENORS: [f64; 8] = [0.25, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 30.0]; const RATES: [f64; 8] = [0.0383, 0.0398, 0.0408, 0.0428, 0.0434, 0.0445, 0.0475, 0.0527]; // ----------------------------------------------------------------- // The monotone convex scheme, checked against the four properties it // is defined by rather than against a reference implementation. // ----------------------------------------------------------------- fn monotone_convex() -> Curve { bootstrap_par(&TENORS, &RATES, 2.0, Interp::MonotoneConvex) } #[test] fn monotone_convex_is_conserving() { // The property the whole construction is built on, and the one that // makes it reprice: the area under the forward curve across each bucket // is the average the quotes pinned. Checked by integrating the forwards // numerically and comparing against the closed form the curve uses -- // two routes that share only the definition of g. let c = monotone_convex(); for &t in &[0.1, 0.25, 0.7, 1.0, 1.6, 3.0, 4.4, 10.0, 19.0, 30.0] { let steps = 200_000; let h = t / steps as f64; let mut area = 0.0; for i in 0..steps { area += c.inst_forward((i as f64 + 0.5) * h) * h; } let analytic = -c.df(t).ln(); assert!( (area - analytic).abs() < 1e-8, "at {t}y: numerical {area} against analytic {analytic}" ); } } #[test] fn monotone_convex_has_continuous_forwards() { // The failing of two of the other three schemes. Approach every node // from both sides and the forward must agree, which it does because the // two buckets meeting there are built from the same node value. let c = monotone_convex(); for &t in TENORS.iter() { let (before, after) = (c.inst_forward(t - 1e-7), c.inst_forward(t + 1e-7)); assert!( (before - after).abs() < 1e-5, "forward jumped at {t}y: {before} then {after}" ); } } #[test] fn monotone_convex_keeps_forwards_positive() { // What the collar on the node values is for. A curve that dips negative // between quotes is not merely ugly, it prices a bond above par for no // reason the market gave. let c = monotone_convex(); for i in 0..3000 { let t = i as f64 * 0.01; let f = c.inst_forward(t); assert!(f >= 0.0, "forward was {f} at {t}y"); } } #[test] fn monotone_convex_does_not_overshoot_the_way_linear_forwards_do() { // The reason the scheme exists. Both interpolate forwards and both are // continuous; the difference is that linear interpolation is dragged // past the data by a distant node, and this one is not. On the long end // of this curve, where the quotes are twenty years apart, that is the // difference between a forward that reaches 7.1% and one that does not // leave the range the quotes support. let range = |c: &Curve| { let (mut lo, mut hi) = (f64::MAX, f64::MIN); for i in 0..3000 { let f = c.inst_forward(i as f64 * 0.01); lo = lo.min(f); hi = hi.max(f); } (lo, hi) }; let (_, mc_high) = range(&monotone_convex()); let (_, lf_high) = range(&bootstrap_par(&TENORS, &RATES, 2.0, Interp::LinearForward)); assert!( mc_high < lf_high - 0.005, "monotone convex peaked at {mc_high}, linear forwards at {lf_high}" ); // And it stays within reach of the longest quote rather than running // away above it. assert!(mc_high < 0.065, "monotone convex still overshot to {mc_high}"); } #[test] fn monotone_convex_is_not_local() { // The cost the chapter warns about, made concrete. Moving the thirty // year quote moves the curve at ten years, because the forward at each // node is built from the buckets on both sides of it. Under a local // scheme the same bump changes nothing before its own node. let mut bumped = RATES; bumped[7] += 0.0010; let reach = |interp: Interp| { let base = bootstrap_par(&TENORS, &RATES, 2.0, interp); let moved = bootstrap_par(&TENORS, &bumped, 2.0, interp); (moved.inst_forward(9.0) - base.inst_forward(9.0)).abs() }; assert!(reach(Interp::LinearLogDf) < 1e-12, "log-df scheme was not local"); assert!( reach(Interp::MonotoneConvex) > 1e-5, "monotone convex looked local, which would mean the node forwards \ were not seeing the bucket beyond them" ); } fn all_schemes() -> [Interp; 4] { [ Interp::LinearZero, Interp::LinearLogDf, Interp::LinearForward, Interp::MonotoneConvex, ] } #[test] fn every_scheme_reprices_every_input() { // The whole argument of the section depends on this: the schemes are not // better and worse fits, they are all exact, and they still disagree. for interp in all_schemes() { let c = bootstrap_par(&TENORS, &RATES, 2.0, interp); for (&t, &r) in TENORS.iter().zip(&RATES) { let price = if t <= 0.5 { c.df(t) * (1.0 + r * t) } else { par_bond_price(&c, t, r, 2.0) }; assert!( (price - 1.0).abs() < 1e-9, "{interp:?} misprices the {t}y at {price}" ); } } } /// Widest disagreement across the schemes, in basis points. fn spread_bp(curves: &[Curve], f: impl Fn(&Curve) -> f64) -> f64 { let v: Vec<f64> = curves.iter().map(&f).collect(); (v.iter().cloned().fold(f64::MIN, f64::max) - v.iter().cloned().fold(f64::MAX, f64::min)) * 1e4 } #[test] fn the_disagreement_grows_with_each_derivative() { let curves: Vec<Curve> = all_schemes() .iter() .map(|&i| bootstrap_par(&TENORS, &RATES, 2.0, i)) .collect(); // At a node whose instrument pays only on nodes, the schemes cannot // disagree at all: no interpolated value entered the calibration. assert!(spread_bp(&curves, |c| c.zero(1.0)) < 1e-6); // Elsewhere they do, and this is the part that is easy to get wrong: the // node values themselves differ, because a 10y instrument pays coupons // between the nodes and those payments were valued by interpolation. The // schemes are not three readings of one curve, they are three curves. assert!(spread_bp(&curves, |c| c.zero(30.0)) > 1.0); // And the forwards, being a derivative of all that, disagree by an order // of magnitude more than the zero rates that carry them. // // Compared as maxima over the curve, not point by point. Pointwise the // claim is false and instructively so: the zero rate at t accumulates // every disagreement between 0 and t, whereas the forward at t is local, // so at a t where the schemes happen to cross there is a wide zero // spread and almost no forward spread. It is the worst case along the // curve that matters, since that is where something will be priced. let mut worst_zero: f64 = 0.0; let mut worst_forward: f64 = 0.0; for i in 1..=600 { let t = i as f64 * 0.05; worst_zero = worst_zero.max(spread_bp(&curves, |c| c.zero(t))); worst_forward = worst_forward.max(spread_bp(&curves, |c| c.inst_forward(t))); } assert!( worst_forward > 3.0 * worst_zero, "forwards spread at most {worst_forward}bp against {worst_zero}bp for zeros" ); } #[test] fn linear_zero_forwards_jump_at_nodes() { // The pathology the curve construction chapter derives: z' jumps at a // node, so f = z + t z' jumps too, by t times the change in slope. let c = bootstrap_par(&TENORS, &RATES, 2.0, Interp::LinearZero); let node = 10.0; let before = c.inst_forward(node - 1e-7); let after = c.inst_forward(node + 1e-7); assert!( (after - before).abs() > 1e-3, "expected a visible jump at {node}, got {before} -> {after}" ); } #[test] fn linear_forward_scheme_is_continuous() { let c = bootstrap_par(&TENORS, &RATES, 2.0, Interp::LinearForward); for &node in &TENORS[1..] { let before = c.inst_forward(node - 1e-7); let after = c.inst_forward(node + 1e-7); assert!( (after - before).abs() < 1e-6, "forwards jumped at {node}: {before} -> {after}" ); } } #[test] fn log_df_scheme_has_piecewise_constant_forwards() { let c = bootstrap_par(&TENORS, &RATES, 2.0, Interp::LinearLogDf); // Flat strictly inside a segment. let a = c.inst_forward(5.5); let b = c.inst_forward(9.5); assert!((a - b).abs() < 1e-12, "{a} vs {b} inside one segment"); } #[test] fn zero_rate_is_the_average_of_the_forwards() { // The identity the curve construction chapter leans on to explain why // forwards travel further than zeros: z(T) = (1/T) integral of f. let c = bootstrap_par(&TENORS, &RATES, 2.0, Interp::LinearForward); let t = 7.0; let n = 20000; let mut acc = 0.0; for i in 0..n { acc += c.inst_forward((i as f64 + 0.5) / n as f64 * t); } let mean = acc / n as f64; assert!((mean - c.zero(t)).abs() < 1e-6, "{mean} vs {}", c.zero(t)); }} /// A swap's two legs and its par rate, from a discount curve.////// The curve construction chapter defines the par rate as the fixed rate making/// the swap worth zero and derives////// ```text/// S = (P(t, tau_0) - P(t, tau_N)) / sum_i delta_i P(t, tau_i),/// ```////// the numerator being the floating leg. The chapter has no code behind it, and/// the identity is worth checking rather than assumed: the floating leg/// telescoping to a difference of two discount factors is the step that does the/// work, and it holds only because each floating payment is exactly the forward/// implied by the same curve.pub struct Swap { /// Payment times of the fixed leg, in years. pub fixed_times: Vec<f64>, /// Accruals of the fixed leg. pub fixed_accruals: Vec<f64>, /// Payment times of the floating leg. pub float_times: Vec<f64>, /// Accruals of the floating leg. pub float_accruals: Vec<f64>, /// Start of both legs. pub start: f64,} impl Swap { /// The annuity: the fixed leg's value per unit of rate. pub fn annuity(&self, discount: impl Fn(f64) -> f64) -> f64 { self.fixed_times .iter() .zip(&self.fixed_accruals) .map(|(&t, &d)| d * discount(t)) .sum() } /// The floating leg, valued payment by payment from the forward rates the /// curve implies --- deliberately *not* by the telescoping shortcut, so that /// the shortcut can be tested against it. pub fn floating_leg_from_forwards(&self, discount: impl Fn(f64) -> f64) -> f64 { let mut total = 0.0; let mut previous = self.start; for (&t, &d) in self.float_times.iter().zip(&self.float_accruals) { // The simply compounded forward over (previous, t]. let forward = (discount(previous) / discount(t) - 1.0) / d; total += d * forward * discount(t); previous = t; } total } /// The same leg by the telescoping identity the chapter uses. pub fn floating_leg_telescoped(&self, discount: impl Fn(f64) -> f64) -> f64 { let last = *self.float_times.last().expect("a leg needs a payment"); discount(self.start) - discount(last) } /// The par rate. pub fn par_rate(&self, discount: impl Fn(f64) -> f64) -> f64 { self.floating_leg_telescoped(&discount) / self.annuity(&discount) }} #[cfg(test)]mod swap_tests { use super::*; /// A five year annual swap starting now. fn swap() -> Swap { let times: Vec<f64> = (1..=5).map(|i| i as f64).collect(); Swap { fixed_times: times.clone(), fixed_accruals: vec![1.0; 5], float_times: times, float_accruals: vec![1.0; 5], start: 0.0, } } /// A curve with shape, so nothing is verified only for a flat one. fn curve(t: f64) -> f64 { (-(0.02 * t + 0.004 * t * t / (1.0 + t))).exp() } #[test] fn the_floating_leg_telescopes() { // The step the chapter's formula rests on. Valuing each floating payment // from its own forward and discounting it must give the difference of two // discount factors, and it does so exactly rather than approximately -- // each term's numerator cancels the next term's denominator. let s = swap(); let explicit = s.floating_leg_from_forwards(curve); let shortcut = s.floating_leg_telescoped(curve); assert!( (explicit - shortcut).abs() < 1e-14, "explicit {explicit:.16} against telescoped {shortcut:.16}" ); } #[test] fn the_par_rate_makes_the_swap_worth_nothing() { // The definition, checked against the formula. Paying the par rate on the // fixed leg and receiving the floating leg should net to zero. let s = swap(); let rate = s.par_rate(curve); let value = s.floating_leg_from_forwards(curve) - rate * s.annuity(curve); assert!(value.abs() < 1e-14, "a par swap is worth {value:.2e}"); // And it is a sensible number for this curve rather than an artefact. assert!((0.015..0.035).contains(&rate), "par rate {rate:.5}"); } #[test] fn a_mismatched_floating_frequency_still_telescopes() { // The telescoping does not depend on the two legs agreeing, which is why // the chapter can quote one formula for a swap paying quarterly against // annual. Semi-annual floating against annual fixed: let semi: Vec<f64> = (1..=10).map(|i| i as f64 * 0.5).collect(); let s = Swap { fixed_times: (1..=5).map(|i| i as f64).collect(), fixed_accruals: vec![1.0; 5], float_times: semi, float_accruals: vec![0.5; 10], start: 0.0, }; let explicit = s.floating_leg_from_forwards(curve); let shortcut = s.floating_leg_telescoped(curve); assert!((explicit - shortcut).abs() < 1e-14, "{explicit:.16} vs {shortcut:.16}"); // The par rate differs from the annual-floating one only in the third // decimal of a per cent, since the legs discount the same cashflows -- a // reminder that the floating frequency is nearly irrelevant to the par // rate and matters for the basis rather than the level. let annual = swap().par_rate(curve); assert!((s.par_rate(curve) - annual).abs() < 1e-12); }} // ---------------------------------------------------------------------------// Conventions: the arithmetic that is not modelling and still moves the price.// --------------------------------------------------------------------------- /// The day count conventions the linear products chapter needs.////// A year fraction is not a length of time. It is what the contract says it is,/// and the conventions disagree by amounts that dwarf a bid-offer spread.#[derive(Clone, Copy, Debug, PartialEq, Eq)]pub enum DayCount { /// Actual days over 360. The money-market convention: USD Libor and SOFR /// legs, EUR floating legs. Act360, /// Actual days over 365. Sterling and several Asian markets. Act365, /// Thirty-day months over a 360-day year. US corporate and agency bonds, and /// USD swap fixed legs by convention. Thirty360,} impl DayCount { /// The year fraction between two dates, each given as `(year, month, day)`. /// /// `Act/*` count real days; `30/360` replaces the calendar with twelve equal /// months, using the US (Bond Basis) end-of-month rule. pub fn year_fraction(&self, from: (i32, u32, u32), to: (i32, u32, u32)) -> f64 { match self { DayCount::Act360 => actual_days(from, to) as f64 / 360.0, DayCount::Act365 => actual_days(from, to) as f64 / 365.0, DayCount::Thirty360 => { let (y1, m1, mut d1) = from; let (y2, m2, mut d2) = to; if d1 == 31 { d1 = 30; } if d2 == 31 && d1 == 30 { d2 = 30; } let days = 360 * (y2 - y1) + 30 * (m2 as i32 - m1 as i32) + (d2 as i32 - d1 as i32); days as f64 / 360.0 } } }} /// Days between two dates, by conversion to a day number. Proleptic Gregorian,/// which is all any contract in these notes needs.fn actual_days(from: (i32, u32, u32), to: (i32, u32, u32)) -> i64 { day_number(to) - day_number(from)} fn day_number((y, m, d): (i32, u32, u32)) -> i64 { // Fliegel and van Flandern's formula, which handles the leap rule without // a table of month lengths. let (y, m) = (y as i64, m as i64); let a = (14 - m) / 12; let y = y + 4800 - a; let m = m + 12 * a - 3; d as i64 + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100 + y / 400 - 32045} /// Accrued interest and the two prices a bond has.////// A bond quote is a *clean* price: it excludes the coupon interest earned since/// the last payment. What settles is the *dirty* price, clean plus accrued. The/// distinction exists so that a quote does not sawtooth by the coupon over the/// period, and it is a reporting convention rather than economics --- but a/// position marked on the wrong one is wrong by up to a full coupon.pub struct BondPrices { pub clean: f64, pub accrued: f64, pub dirty: f64,} /// Split a dirty price into the quote and the accrual.////// `since` and `period` are year fractions on the bond's own convention: time/// since the last coupon and the length of the current period.pub fn clean_and_dirty(dirty: f64, coupon: f64, since: f64, period: f64) -> BondPrices { // The coupon paid over the whole period is 100 * coupon * period and the // fraction earned is since / period, so the two periods cancel and the // accrual is simply the rate times the elapsed year fraction. `period` is // still taken as an argument because the cancellation only holds when both // are measured on the bond's own convention, and passing it makes that // requirement visible at every call site. let _ = period; let accrued = 100.0 * coupon * since; BondPrices { clean: dirty - accrued, accrued, dirty }} #[cfg(test)]mod convention_tests { use super::*; #[test] fn the_conventions_disagree_by_more_than_a_spread() { // A semiannual period, 15 January to 15 July, on the three conventions // in use. The chapter quotes these. let (from, to) = ((2026, 1, 15), (2026, 7, 15)); let act360 = DayCount::Act360.year_fraction(from, to); let act365 = DayCount::Act365.year_fraction(from, to); let thirty = DayCount::Thirty360.year_fraction(from, to); // 181 actual days. assert_eq!(actual_days(from, to), 181); assert!((act360 - 181.0 / 360.0).abs() < 1e-12); assert!((act365 - 181.0 / 365.0).abs() < 1e-12); assert!((thirty - 0.5).abs() < 1e-12); // On a 5% rate the interest paid over this one period differs by this // much between Act/360 and 30/360, in basis points of notional. let gap = 10_000.0 * 0.05 * (act360 - thirty); assert!( (gap - 1.39).abs() < 0.01, "Act/360 against 30/360 on one period: {gap:.2}bp" ); println!( "act360={act360:.6} act365={act365:.6} 30/360={thirty:.6}, gap {gap:.2}bp" ); } #[test] fn a_leap_day_is_worth_a_basis_point() { // The same period a year earlier and later, across a leap year, on // Act/360. February 2028 has 29 days, so an annual period spanning it is // longer and pays more on the same rate. let ordinary = DayCount::Act360.year_fraction((2026, 1, 15), (2027, 1, 15)); let leap = DayCount::Act360.year_fraction((2027, 8, 15), (2028, 8, 15)); assert_eq!( (ordinary * 360.0).round() as i64, 365, "an ordinary year has 365 days" ); assert_eq!((leap * 360.0).round() as i64, 366, "this one spans 29 February"); let gap = 10_000.0 * 0.05 * (leap - ordinary); // One day out of 360 at five per cent, which is 1.39bp and not the // 0.139 that dropping a factor of ten suggests -- the same size as the // convention gap above, and for the same reason: both are one day. assert!((gap - 1.389).abs() < 0.002, "one day at 5%: {gap:.3}bp"); } #[test] fn accrued_interest_is_the_whole_of_the_difference() { // Halfway through a semiannual period on a 4% coupon: one per cent of // face has been earned, so clean and dirty differ by exactly that. let p = clean_and_dirty(101.0, 0.04, 0.25, 0.5); assert!((p.accrued - 1.0).abs() < 1e-12, "accrued {:.6}", p.accrued); assert!((p.clean - 100.0).abs() < 1e-12); assert!((p.dirty - p.clean - p.accrued).abs() < 1e-12); }} // ---------------------------------------------------------------------------// The curve as a posterior rather than as a solution.// --------------------------------------------------------------------------- /// Curve construction done as Gaussian process regression.////// The bootstrap of the linear products chapter finds *the* curve repricing the/// quotes exactly. Two things are wrong with that as a description of the/// problem. The quotes are mids of a bid-offer, so repricing them exactly is/// fitting noise; and between and beyond the quoted maturities the curve is not/// determined at all, which an exact construction hides by producing one answer.////// The Bayesian version states both. Put a prior on the forward curve --- a/// Gaussian process, whose kernel encodes how smooth it is believed to be ---/// treat each quote as a noisy linear functional of it, and report the/// posterior. Since a zero rate is an average of forwards,////// ```text/// z(T) = (1/T) integral_0^T f(u) du,/// ```////// the observation is linear and the posterior is Gaussian in closed form.////// What comes out is not a curve but a distribution over curves: a mean, which/// is the fitted curve, and a standard deviation at every maturity, which says/// where the market has spoken and where it has not.pub struct CurvePosterior { /// Maturities of the grid the forward curve is represented on. pub grid: Vec<f64>, /// Posterior mean of the instantaneous forward rate. pub mean: Vec<f64>, /// Posterior standard deviation of it. pub sd: Vec<f64>,} /// Fit a forward curve to zero rate quotes, with a smoothness prior.////// `prior_level` is the forward rate assumed in the absence of data,/// `prior_sd` how far from it the curve is believed able to stray, and/// `length_scale` how far apart two maturities must be before their forward/// rates are treated as unrelated --- which is the quantitative form of "the/// curve is smooth".pub fn fit_curve( grid: &[f64], quote_maturities: &[f64], quote_zeros: &[f64], quote_error: f64, prior_level: f64, prior_sd: f64, length_scale: f64,) -> CurvePosterior { let n = grid.len(); let m = quote_maturities.len(); let step = grid[1] - grid[0]; // Prior covariance of the forward curve: squared exponential. let kernel = |a: f64, b: f64| { let d = (a - b) / length_scale; prior_sd * prior_sd * (-0.5 * d * d).exp() }; let k: Vec<Vec<f64>> = grid .iter() .map(|&a| grid.iter().map(|&b| kernel(a, b)).collect()) .collect(); // Observation matrix: a zero rate is the average of the forwards up to it. let mut h = vec![vec![0.0; n]; m]; for (j, &t) in quote_maturities.iter().enumerate() { let count = grid.iter().filter(|&&g| g <= t + 1e-12).count().max(1); for i in 0..count { h[j][i] = step / (count as f64 * step); } } // K H' and H K H' + R. let kh: Vec<Vec<f64>> = (0..n) .map(|i| (0..m).map(|j| (0..n).map(|l| k[i][l] * h[j][l]).sum()).collect()) .collect(); let mut s = vec![vec![0.0; m]; m]; for a in 0..m { for b in 0..m { s[a][b] = (0..n).map(|i| h[a][i] * kh[i][b]).sum::<f64>(); if a == b { s[a][b] += quote_error * quote_error; } } } // Residual of the prior against the quotes, and S^{-1} applied to it. let residual: Vec<f64> = (0..m) .map(|j| quote_zeros[j] - (0..n).map(|i| h[j][i] * prior_level).sum::<f64>()) .collect(); let alpha = solve(&s, &residual); let mean: Vec<f64> = (0..n) .map(|i| prior_level + (0..m).map(|j| kh[i][j] * alpha[j]).sum::<f64>()) .collect(); // Posterior variance: K - K H' S^{-1} H K, diagonal only. let sd: Vec<f64> = (0..n) .map(|i| { let row: Vec<f64> = (0..m).map(|j| kh[i][j]).collect(); let beta = solve(&s, &row); let reduction: f64 = (0..m).map(|j| row[j] * beta[j]).sum(); (k[i][i] - reduction).max(0.0).sqrt() }) .collect(); CurvePosterior { grid: grid.to_vec(), mean, sd }} /// Gaussian elimination with partial pivoting, for the small systems above.fn solve(matrix: &[Vec<f64>], rhs: &[f64]) -> Vec<f64> { let n = rhs.len(); let mut a: Vec<Vec<f64>> = matrix .iter() .zip(rhs) .map(|(row, &b)| { let mut r = row.clone(); r.push(b); r }) .collect(); for i in 0..n { let pivot = (i..n).max_by(|&x, &y| a[x][i].abs().partial_cmp(&a[y][i].abs()).unwrap()).unwrap(); a.swap(i, pivot); for j in (i + 1)..n { let factor = a[j][i] / a[i][i]; for c in i..=n { a[j][c] -= factor * a[i][c]; } } } let mut x = vec![0.0; n]; for i in (0..n).rev() { let sum: f64 = ((i + 1)..n).map(|j| a[i][j] * x[j]).sum(); x[i] = (a[i][n] - sum) / a[i][i]; } x} #[cfg(test)]mod posterior_tests { use super::*; fn grid() -> Vec<f64> { (1..=120).map(|i| i as f64 * 0.25).collect() } #[test] fn the_band_is_tightest_between_the_quotes_and_not_at_them() { // The finding, and it is not the obvious one. // // A zero rate is an *average* of forwards up to its maturity, not an // observation of the forward at it. So quoting the two year zero does // not pin the two year forward; it pins the average over [0, 2], which // determines the middle of that interval well and its endpoints // poorly. The posterior band is therefore narrowest in the interior of // the quoted region and widest at the quoted maturities themselves -- // the opposite of the picture a pointwise observation would give, and // a fact about what the instruments actually say. let g = grid(); let quotes = [2.0, 5.0, 10.0]; let zeros = [0.040, 0.042, 0.044]; let p = fit_curve(&g, "es, &zeros, 1e-4, 0.04, 0.02, 3.0); let at = |t: f64| { let i = g.iter().position(|&x| (x - t).abs() < 1e-9).unwrap(); p.sd[i] }; assert!( at(1.0) < at(2.0) / 5.0, "the interior should be far better determined than the quote maturity: \ {:.5} against {:.5}", at(1.0), at(2.0) ); // Past the last quote the band grows back to the prior over roughly the // length scale, which is the honest statement about extrapolation: it is // exactly as uncertain as the prior says, and a construction that // extends a flat forward past the last quote is asserting that prior // without admitting to it. assert!(at(30.0) > at(10.0) * 1.9, "extrapolation should revert to the prior"); assert!( (at(30.0) - 0.02).abs() < 1e-6, "and reach it exactly: {:.6} against the prior 0.02", at(30.0) ); assert!(p.sd.iter().all(|s| *s <= 0.02 + 1e-12), "never wider than the prior"); } #[test] fn quotes_are_reproduced_within_their_error_and_not_exactly() { // The other half of the point. A bootstrap hits the mid exactly; this // hits it to within the bid-offer, which is all the information the mid // carries. let g = grid(); let quotes = [2.0, 5.0, 10.0]; let zeros = [0.040, 0.042, 0.044]; let error = 5e-4; let p = fit_curve(&g, "es, &zeros, error, 0.04, 0.02, 3.0); let step = g[1] - g[0]; for (j, &t) in quotes.iter().enumerate() { let count = g.iter().filter(|&&x| x <= t + 1e-12).count(); let fitted: f64 = p.mean[..count].iter().sum::<f64>() / count as f64; let miss = (fitted - zeros[j]).abs(); assert!(miss < error, "quote {t}y missed by {:.6}", miss); assert!(miss > 1e-9, "and it should not be hit exactly either"); let _ = step; } } #[test] fn a_stronger_smoothness_prior_narrows_the_band_and_costs_fit() { // The tradeoff made explicit, which is the whole reason to write it this // way: the length scale is the interpolation choice, and lengthening it // buys confidence between the quotes by assuming more. let g = grid(); let quotes = [2.0, 5.0, 10.0]; let zeros = [0.040, 0.042, 0.044]; let gap = |length: f64| { let p = fit_curve(&g, "es, &zeros, 1e-4, 0.04, 0.02, length); let i = g.iter().position(|&x| (x - 7.5).abs() < 1e-9).unwrap(); p.sd[i] }; assert!(gap(8.0) < gap(1.0), "a longer length scale should narrow the gap"); }} #[cfg(test)]mod stage_tests { use super::*; fn market() -> (Vec<f64>, Vec<f64>) { ( vec![1.0, 2.0, 3.0, 5.0, 7.0, 10.0, 20.0, 30.0], vec![0.0412, 0.0408, 0.0406, 0.0411, 0.0421, 0.0435, 0.0477, 0.0488], ) } #[test] fn the_stages_end_where_the_bootstrap_does() { // The animation must show the construction the chapter describes, not a // parallel one, so its last frame has to be the curve the ordinary // bootstrap returns. let (t, r) = market(); for interp in [Interp::LinearForward, Interp::MonotoneConvex] { let stages = bootstrap_par_stages(&t, &r, 2.0, interp); let finished = bootstrap_par(&t, &r, 2.0, interp); let last = stages.last().expect("at least one stage"); for (a, b) in last.yields.iter().zip(&finished.yields) { assert!((a - b).abs() < 1e-10, "{interp:?}: {a} against {b}"); } } } #[test] fn adding_an_instrument_moves_the_curve_already_built() { // What a new quote does to the part already built, and it is worth being // exact about which part. // // The *zero* at an existing node does not move: a three year par bond // pays every coupon at or before three years, so nothing added later // enters its pricing. The *forward* does move, because under this scheme // the forward at a node is built from the buckets on both sides of it, // and until the five year quote arrives the bucket beyond three years // does not exist. That is the non-locality, and it is what the animation // shows. let (t, r) = market(); let stages = bootstrap_par_stages(&t, &r, 2.0, Interp::MonotoneConvex); // Where the conservation argument fails, measured. The scheme conserves // whole bucket areas, so a zero at a node ought to be settled once the // buckets below it are --- and it would be, if the instrument paid only // at nodes. A three year par bond pays semiannually, so its coupon at // two and a half years is discounted through half a bucket, and half a // bucket is not conserved. let full_area = |c: &Curve| 3.0 * c.zero(3.0) - 2.0 * c.zero(2.0); let half_area = |c: &Curve| 2.5 * c.zero(2.5) - 2.0 * c.zero(2.0); let full_move = (full_area(&stages[5]) - full_area(&stages[2])).abs(); let half_move = (half_area(&stages[5]) - half_area(&stages[2])).abs(); println!("bucket [2,3] area moves {full_move:.3e}, half of it moves {half_move:.3e}"); assert!( half_move > 40.0 * full_move, "the partial area must move far more than the whole one, or the \ explanation in the chapter is wrong: {half_move:.3e} against {full_move:.3e}" ); let zero_move = (stages[2].zero(3.0) - stages[5].zero(3.0)).abs() * 1e4; let fwd_move = (stages[2].inst_forward(3.0) - stages[5].inst_forward(3.0)).abs() * 1e4; println!("adding three quotes moves the 3y zero by {zero_move:.4}bp, the forward by {fwd_move:.2}bp"); // The zero barely moves --- but it does move, and that it moves at all // is the non-locality reaching further than one would guess. Every // coupon of the three year bond falls at or before three years, so the // naive expectation is that nothing added later can touch it. Under this // scheme a coupon at two and a half years is discounted through a // forward built from the bucket beyond three, so it can. assert!(zero_move > 1e-4, "the zero should move a little: {zero_move:.6}bp"); assert!(zero_move < 0.01, "but only a little: {zero_move:.6}bp"); // The forward at the same point moves by orders of magnitude more, and // that is the curve anything fixing between the nodes actually pays. assert!( fwd_move > 100.0 * zero_move, "the forward should move far more: {fwd_move:.4}bp against {zero_move:.4}bp" ); } #[test] fn only_the_non_local_scheme_needs_sweeping() { // Three of the four schemes are done in one pass; the monotone convex // one is not, and the extra stages are exactly that cost. let (t, r) = market(); let local = bootstrap_par_stages(&t, &r, 2.0, Interp::LinearForward); let nonlocal = bootstrap_par_stages(&t, &r, 2.0, Interp::MonotoneConvex); assert_eq!(local.len(), t.len(), "one stage per instrument and no more"); assert!( nonlocal.len() > t.len(), "the monotone convex build should carry sweeps beyond the instruments" ); }}