quant/src/risk.rs
Risk measures, and the arithmetic of the risk management chapter's counterexample.
//! Risk measures, and the arithmetic of the risk management chapter's counterexample.//!//! Two measures, defined on a discrete loss distribution because that is what a//! real risk system has: a set of scenarios with weights, not a formula.//!//! The discreteness is not incidental. Value at Risk fails to be subadditive,//! and the cleanest demonstration of the failure needs a distribution with//! atoms; on smooth distributions the failure is real but easy to miss.//! Expected shortfall on such a distribution also has to be defined by its//! integral rather than as a conditional expectation, and the difference is not//! cosmetic — see [`expected_shortfall`]. /// A loss distribution, as pairs of loss and probability.////// Losses are positive numbers: a loss of 100 is worse than a loss of 0. Signs/// in risk are a perennial source of error, so this module never uses a P&L/// convention anywhere.pub struct Losses(pub Vec<(f64, f64)>); impl Losses { /// Build from scenarios, sorting and merging as it goes. pub fn new(mut scenarios: Vec<(f64, f64)>) -> Self { scenarios.retain(|&(_, p)| p > 0.0); scenarios.sort_by(|a, b| a.0.partial_cmp(&b.0).expect("no NaN losses")); Losses(scenarios) } /// Total probability, which should be one. pub fn mass(&self) -> f64 { self.0.iter().map(|&(_, p)| p).sum() } pub fn mean(&self) -> f64 { self.0.iter().map(|&(l, p)| l * p).sum() } /// Value at Risk at confidence `alpha`: the `alpha` quantile of the loss. /// /// The smallest loss level that the loss fails to exceed with probability at /// least `alpha`. With atoms this is a step function of `alpha`, which is /// where its bad behaviour comes from. pub fn value_at_risk(&self, alpha: f64) -> f64 { let mut cumulative = 0.0; for &(loss, p) in &self.0 { cumulative += p; if cumulative >= alpha { return loss; } } self.0.last().map(|&(l, _)| l).unwrap_or(0.0) } /// Expected shortfall at confidence `alpha`. /// /// Defined as the average of the Value at Risk over the worst `1 - alpha` of /// the distribution: /// /// ```text /// ES_alpha = 1/(1-alpha) * integral over u in [alpha, 1] of VaR_u du /// ``` /// /// Not as `E[L | L >= VaR_alpha]`. The two agree when the loss has a /// continuous distribution and disagree when it has atoms, which is exactly /// the case the counterexample is built from. Getting this wrong makes /// expected shortfall look as badly behaved as Value at Risk: on a single /// defaultable bond the conditional-expectation version returns the mean /// loss of 4, where the correct answer is 80. /// /// The integral form is also the one that is coherent, which is the whole /// reason for preferring the measure. pub fn expected_shortfall(&self, alpha: f64) -> f64 { if !(0.0..1.0).contains(&alpha) { return f64::NAN; } let mut cumulative = 0.0; let mut integral = 0.0; for &(loss, p) in &self.0 { let (lo, hi) = (cumulative, cumulative + p); cumulative = hi; // The part of this atom's slab of quantiles that lies beyond alpha. let (a, b) = (lo.max(alpha), hi.min(1.0)); if b > a { integral += loss * (b - a); } } integral / (1.0 - alpha) }} /// The loss distribution of `n` independent defaultable bonds, equally/// weighted, with total notional `notional`.////// The risk management chapter's example. Each bond defaults with probability/// `p` and is then worth nothing; the portfolio loses `notional * k / n` when/// `k` of them default, and `k` is binomial.pub fn defaultable_portfolio(n: usize, p: f64, notional: f64) -> Losses { let mut scenarios = Vec::with_capacity(n + 1); for k in 0..=n { scenarios.push((notional * k as f64 / n as f64, binomial_pmf(n, k, p))); } Losses::new(scenarios)} /// `P(k successes in n trials)`, computed multiplicatively.////// Not by way of a factorial: at the sizes a figure uses, `n!` overflows a/// double long before the answer does.fn binomial_pmf(n: usize, k: usize, p: f64) -> f64 { let mut out = (1.0 - p).powi((n - k) as i32) * p.powi(k as i32); for i in 0..k { out *= (n - i) as f64 / (i + 1) as f64; } out} #[cfg(test)]mod tests { use super::*; const ALPHA: f64 = 0.95; const P: f64 = 0.04; #[test] fn the_portfolio_is_a_probability_distribution() { for n in 1..=12 { let l = defaultable_portfolio(n, P, 100.0); assert!((l.mass() - 1.0).abs() < 1e-12, "mass at n={n} was {}", l.mass()); // Expected loss does not depend on how finely it is split. assert!((l.mean() - 100.0 * P).abs() < 1e-9, "mean at n={n} was {}", l.mean()); } } #[test] fn value_at_risk_is_not_subadditive() { // The risk management chapter's counterexample, as a computation. // // One bond: it survives with probability 0.96, so the 95% quantile of // the loss is zero. Two half-sized bonds: both survive with probability // 0.9216, which is below 95%, so the quantile jumps to the one-default // loss of 50. // // Splitting one position into two independent halves therefore takes the // measured risk from nothing to half the notional. A measure that // punishes diversification will push a firm holding it towards // concentration. let one = defaultable_portfolio(1, P, 100.0); let two = defaultable_portfolio(2, P, 100.0); assert_eq!(one.value_at_risk(ALPHA), 0.0); assert_eq!(two.value_at_risk(ALPHA), 50.0); assert!( two.value_at_risk(ALPHA) > 2.0 * one.value_at_risk(ALPHA) / 2.0, "the counterexample has stopped working" ); } #[test] fn expected_shortfall_rewards_diversification() { // The same sequence under the coherent measure. It falls at every step, // which is what one wants a risk number to do as a position is spread // over more independent names. let mut previous = f64::INFINITY; for n in 1..=10 { let es = defaultable_portfolio(n, P, 100.0).expected_shortfall(ALPHA); assert!(es < previous, "ES rose at n={n}: {es} after {previous}"); previous = es; } } #[test] fn expected_shortfall_is_not_the_conditional_expectation_here() { // The trap the doc comment warns about. On a single bond the correct // answer is 80 — the average of the worst 5% of outcomes, which are // dominated by the 4% chance of losing everything — and the // conditional-expectation formula gives the unconditional mean of 4, // because the whole distribution lies at or above the zero quantile. let one = defaultable_portfolio(1, P, 100.0); assert!((one.expected_shortfall(ALPHA) - 80.0).abs() < 1e-9); assert!((one.mean() - 4.0).abs() < 1e-9); } #[test] fn expected_shortfall_dominates_value_at_risk() { // It is an average of quantiles at least as extreme, so it can never be // smaller. A risk system reporting otherwise has a sign or an index off. for n in 1..=12 { let l = defaultable_portfolio(n, P, 100.0); assert!( l.expected_shortfall(ALPHA) >= l.value_at_risk(ALPHA) - 1e-12, "at n={n}" ); } } #[test] fn both_measures_are_translation_invariant() { // Adding a certain loss of c must add exactly c. This is the axiom that // makes a risk number readable as an amount of capital. let base = defaultable_portfolio(4, P, 100.0); let shifted = Losses::new(base.0.iter().map(|&(l, p)| (l + 25.0, p)).collect()); assert!((shifted.value_at_risk(ALPHA) - base.value_at_risk(ALPHA) - 25.0).abs() < 1e-9); assert!( (shifted.expected_shortfall(ALPHA) - base.expected_shortfall(ALPHA) - 25.0).abs() < 1e-9 ); }} /// Where the risk in a hedged book actually comes from.////// The risk management chapter argues that a risk number without a decomposition/// is not a risk number, and that the three sources have to be separated because/// only one of them is a live decision. This measures the split on the simplest/// book that has all three: a short call, delta hedged.////// * *Chosen* --- the position the desk meant to hold. For a hedged option this/// is nearly nothing, which is the point of hedging./// * *Discretisation* --- the hedging error of the no-arbitrage chapter, which/// arises even with a perfectly correct model because rebalancing is not/// continuous. Unavoidable, and shrinking as the square root of the frequency./// * *Model* --- the hedge computed with a wrong volatility. Not a sampling error/// and not reducible by trading more often.////// The three are separated by running the same paths three ways, so the/// comparison is not contaminated by different noise.pub struct RiskSources { /// Standard deviation of profit and loss with the correct hedge. pub with_correct_model: f64, /// The same, hedging with a mis-specified volatility. pub with_wrong_model: f64, /// The part attributable to the mis-specification, in quadrature. pub model_component: f64,} /// Hedge a short call over `steps` rebalances, once with the true volatility and/// once with `hedging_vol`, on identical paths.////// `true_vol` drives the paths; `hedging_vol` is what the desk believes. The gap/// stands in for every reason a delta can be wrong --- the wrong backbone of the/// smile dynamics chapter, a stale surface, a mis-estimated mean reversion --- and/// the point is the size of its contribution rather than its cause.pub fn hedge_risk_sources( spot: f64, strike: f64, true_vol: f64, hedging_vol: f64, expiry: f64, steps: usize, paths: usize, seed: u64,) -> RiskSources { use crate::black::{black76, Side}; use crate::pathwise::call_delta; use crate::pathwise::Rng; let dt = expiry / steps as f64; let premium = black76(spot, strike, hedging_vol, expiry, Side::Call); let mut correct = Vec::with_capacity(paths); let mut wrong = Vec::with_capacity(paths); for path in 0..paths { // The same increments for both hedges, so the difference between them is // the mis-specification and nothing else. let mut rng = Rng::new(seed.wrapping_add(path as u64 * 7919)); let increments: Vec<f64> = (0..steps).map(|_| rng.next_normal()).collect(); for (which, vol) in [(0usize, true_vol), (1, hedging_vol)] { let mut s = spot; let mut hedge = 0.0; for (i, z) in increments.iter().enumerate() { let remaining = expiry - i as f64 * dt; let delta = call_delta(s, strike, vol, remaining); let next = s * ((-0.5 * true_vol * true_vol) * dt + true_vol * dt.sqrt() * z).exp(); hedge += delta * (next - s); s = next; } let pnl = premium + hedge - (s - strike).max(0.0); if which == 0 { correct.push(pnl); } else { wrong.push(pnl); } } } let sd = |v: &[f64]| { let n = v.len() as f64; let mean = v.iter().sum::<f64>() / n; (v.iter().map(|x| (x - mean) * (x - mean)).sum::<f64>() / n).sqrt() }; let (a, b) = (sd(&correct), sd(&wrong)); RiskSources { with_correct_model: a, with_wrong_model: b, model_component: (b * b - a * a).max(0.0).sqrt(), }} #[cfg(test)]mod decomposition_tests { use super::*; const SPOT: f64 = 100.0; const STRIKE: f64 = 100.0; const TRUE_VOL: f64 = 0.20; const EXPIRY: f64 = 1.0; const PATHS: usize = 20_000; const SEED: u64 = 20260809; #[test] fn trading_more_often_removes_one_source_and_not_the_other() { // The distinction the chapter turns on. Discretisation error falls as the // square root of the rebalancing frequency, as the no-arbitrage chapter // measures. Model error does not fall at all: it is a bias in the hedge // ratio, and rebalancing more often applies the wrong ratio more often. let mut correct = Vec::new(); let mut model = Vec::new(); for steps in [50usize, 200, 800] { let r = hedge_risk_sources( SPOT, STRIKE, TRUE_VOL, 0.24, EXPIRY, steps, PATHS, SEED, ); correct.push(r.with_correct_model); model.push(r.model_component); } // Quadrupling the frequency halves the discretisation error. for i in 0..2 { let ratio = correct[i] / correct[i + 1]; assert!( (ratio - 2.0).abs() < 0.25, "correct-model risk fell by {ratio:.2}, expected about 2" ); } // And the model component is essentially unchanged across a sixteenfold // increase in trading. let drift = model[2] / model[0]; assert!( (drift - 1.0).abs() < 0.15, "model risk moved by {drift:.3} over sixteen times the trading" ); } #[test] fn in_a_well_hedged_book_the_model_is_most_of_the_risk() { // The chapter's headline. At a realistic rebalancing frequency and a // volatility error well inside what a surface is marked to, the risk of a // "fully hedged" book is dominated by the hedge being wrong rather than by // anything the position was chosen to hold. let r = hedge_risk_sources( SPOT, STRIKE, TRUE_VOL, 0.24, EXPIRY, 250, PATHS, SEED, ); assert!( r.model_component > r.with_correct_model, "model component {:.4} against discretisation {:.4}", r.model_component, r.with_correct_model ); // And a risk measure applied to the position sees only the smaller part: // the book's stated risk is the correct-model number, which understates // the truth by a factor worth naming. let understatement = r.with_wrong_model / r.with_correct_model; assert!( understatement > 1.5, "the hedged book's risk is understated by {understatement:.2}" ); } #[test] fn the_model_component_scales_with_the_mis_specification() { // A sanity check that the component is measuring what it claims. A larger // volatility error should produce a larger component, and no error should // produce none. let none = hedge_risk_sources( SPOT, STRIKE, TRUE_VOL, TRUE_VOL, EXPIRY, 250, PATHS, SEED, ); assert!(none.model_component < 1e-9, "no error should leave no component"); let mut previous = 0.0; for hedging_vol in [0.21, 0.24, 0.30] { let r = hedge_risk_sources( SPOT, STRIKE, TRUE_VOL, hedging_vol, EXPIRY, 250, PATHS, SEED, ); assert!(r.model_component > previous, "vol={hedging_vol} did not increase it"); previous = r.model_component; } }} /// How much of a book's profit and loss its own risk factors can explain.////// The risk management chapter's attribution is computed before the fact. This is/// the check after it, and it is what a desk and its regulator actually run: take/// the realised profit and loss, predict it from the sensitivities the risk system/// carries, and look at what is left over. Unexplained profit and loss is not/// noise. It is the part of the book's behaviour that its risk factors do not/// span, which is a statement about the model rather than about the market.////// The setting is the one the term structure chapters warn about. Curve moves have/// at least two factors --- a level and a slope --- and a one-factor risk model/// carries a sensitivity to the first only. Whatever the second does is/// unexplained by construction.pub struct ExplainTest { /// Standard deviation of the level factor over the horizon. pub level_sd: f64, /// Standard deviation of the slope factor. pub slope_sd: f64, /// The book's exposure to the level factor. pub level_exposure: f64, /// The book's exposure to the slope factor. pub slope_exposure: f64,} impl ExplainTest { /// The fraction of profit and loss variance a level-only risk model leaves /// unexplained. /// /// With independent factors the answer is available in closed form, and is /// worth writing down because it says the interesting thing directly: the /// unexplained share depends on the book's exposures and not only on how much /// of the curve's variance each factor carries. let level = (self.level_exposure * self.level_sd).powi(2); let slope = (self.slope_exposure * self.slope_sd).powi(2); slope / (level + slope) } /// The same quantity by simulation and regression, which is how a desk /// measures it: regress realised profit and loss on the predicted profit and /// loss and report one minus the coefficient of determination. use crate::pathwise::Rng; let mut rng = Rng::new(seed); let (mut actual, mut predicted) = (Vec::new(), Vec::new()); for _ in 0..paths { let level = self.level_sd * rng.next_normal(); let slope = self.slope_sd * rng.next_normal(); actual.push(self.level_exposure * level + self.slope_exposure * slope); // The risk system knows about the level factor only. predicted.push(self.level_exposure * level); } let n = paths as f64; let mean = |v: &[f64]| v.iter().sum::<f64>() / n; let (ma, mp) = (mean(&actual), mean(&predicted)); let var_a: f64 = actual.iter().map(|x| (x - ma) * (x - ma)).sum::<f64>() / n; let residual: f64 = actual .iter() .zip(&predicted) .map(|(a, p)| { let r = (a - ma) - (p - mp); r * r }) .sum::<f64>() / n; residual / var_a }} #[cfg(test)]mod explain_tests { use super::*; /// Curve factor sizes measured from the committed Treasury history by /// [`curve_factors`], annualised: the level moves about 184 basis points a /// year and carries 79 per cent of the variance, the slope about 69 and /// carries 11. const LEVEL_SD: f64 = 0.0184; const SLOPE_SD: f64 = 0.0069; #[test] fn the_closed_form_matches_the_regression() { // The two routes to the same number, so the formula can be used in the // chapter and the regression is what a desk would run. for (level_exposure, slope_exposure) in [(1.0, 0.0), (1.0, 1.0), (0.2, 1.0), (0.0, 1.0)] { let t = ExplainTest { level_sd: LEVEL_SD, slope_sd: SLOPE_SD, level_exposure, slope_exposure, }; let measured = t.measured_unexplained_share(200_000, 20260809); let exact = t.unexplained_share(); assert!( (measured - exact).abs() < 0.01, "exposures ({level_exposure}, {slope_exposure}): measured {measured:.4} \ against {exact:.4}" ); } } #[test] fn an_outright_position_explains_well_and_a_hedged_one_does_not() { // The chapter's point, and it is the same shape as the attribution // result: hedging removes the factor the model understands and leaves the // ones it does not, so the explain gets worse as the book gets better // hedged. let outright = ExplainTest { level_sd: LEVEL_SD, slope_sd: SLOPE_SD, level_exposure: 1.0, slope_exposure: 0.3, }; assert!( outright.unexplained_share() < 0.02, "an outright book should explain well, got {:.4}", outright.unexplained_share() ); // Now hedge the level exposure down to a tenth without touching the // slope exposure, which is what a duration hedge does to a curve trade. let hedged = ExplainTest { level_exposure: 0.1, slope_exposure: 1.0, ..outright }; assert!( hedged.unexplained_share() > 0.9, "a level-hedged book should explain very badly, got {:.4}", hedged.unexplained_share() ); // Fully level-neutral and it is unexplained entirely, however small the // slope factor's share of curve variance happens to be. let neutral = ExplainTest { level_exposure: 0.0, ..outright }; assert!((neutral.unexplained_share() - 1.0).abs() < 1e-12); } #[test] fn the_slope_factor_being_small_does_not_save_it() { // Why the usual defence fails. The slope factor carries a small share of // the curve's variance -- here about a tenth of the level's -- and it is // tempting to conclude it can be neglected. It cannot, because what // matters is the product of the factor's size with the book's exposure to // it, and a hedged book has arranged for the other product to be small. let ratio = SLOPE_SD / LEVEL_SD; assert!(ratio < 0.45, "the slope factor really is the smaller one"); let hedged = ExplainTest { level_sd: LEVEL_SD, slope_sd: SLOPE_SD, level_exposure: 0.1, slope_exposure: 1.0, }; assert!( hedged.unexplained_share() > 0.9, "and it still accounts for nearly all of a hedged book's variance: {:.3}", hedged.unexplained_share() ); }} /// The factors of a curve, from a history of curves.////// The risk chapter's profit-and-loss explain needs the sizes of the level and/// slope factors, and quoted them as "the kind a principal component analysis/// returns". This computes them from the Treasury history committed in the/// repository instead, so the numbers in the chapter are measurements.////// The analysis is of daily *changes* rather than levels, because a risk model's/// factors describe how a curve moves and a covariance of levels would be/// dominated by where rates happen to have been.pub struct CurveFactors { /// Standard deviations of the factors, largest first, in rate units per day. pub sizes: Vec<f64>, /// The loading of each tenor on each factor, `loadings[factor][tenor]`. pub loadings: Vec<Vec<f64>>, /// Share of total variance explained by each factor. pub shares: Vec<f64>,} /// Principal components of the daily changes, by the power method on the/// covariance matrix with deflation.////// Three factors is all the chapters use and all a curve of this length supports/// distinguishing; the fourth is already at the level of the data's rounding.pub fn curve_factors(observations: &[Vec<f64>], factors: usize) -> CurveFactors { let n = observations.len(); assert!(n > 2, "a covariance needs observations"); let m = observations[0].len(); // Daily changes. let changes: Vec<Vec<f64>> = (1..n) .map(|i| (0..m).map(|j| observations[i][j] - observations[i - 1][j]).collect()) .collect(); let count = changes.len() as f64; let means: Vec<f64> = (0..m).map(|j| changes.iter().map(|c| c[j]).sum::<f64>() / count).collect(); let mut cov = vec![vec![0.0; m]; m]; for c in &changes { for i in 0..m { for j in 0..m { cov[i][j] += (c[i] - means[i]) * (c[j] - means[j]) / count; } } } let total: f64 = (0..m).map(|i| cov[i][i]).sum(); let (mut sizes, mut loadings, mut shares) = (Vec::new(), Vec::new(), Vec::new()); for f in 0..factors.min(m) { // Power iteration. A curve's leading eigenvalue is well separated, so // this converges quickly and needs no more machinery. let mut v: Vec<f64> = (0..m).map(|i| 1.0 / ((i + 1) as f64).sqrt()).collect(); let mut eigenvalue = 0.0; for _ in 0..2000 { let w: Vec<f64> = (0..m).map(|i| (0..m).map(|j| cov[i][j] * v[j]).sum()).collect(); let norm = w.iter().map(|x| x * x).sum::<f64>().sqrt(); if norm <= 0.0 { break; } v = w.iter().map(|x| x / norm).collect(); eigenvalue = norm; } // Sign convention: the first factor should load positively, so that // "level" means what it sounds like. if v.iter().sum::<f64>() < 0.0 && f == 0 { for x in &mut v { *x = -*x; } } sizes.push(eigenvalue.max(0.0).sqrt()); shares.push(eigenvalue / total); loadings.push(v.clone()); // Deflate and find the next. for i in 0..m { for j in 0..m { cov[i][j] -= eigenvalue * v[i] * v[j]; } } } CurveFactors { sizes, loadings, shares }} /// Whether the curve's factors are real, tested against the same panel rather/// than against simulated noise.////// [`sample_correlation_spectrum`](crate::estimation::sample_correlation_spectrum)/// makes the Marchenko-Pastur point on synthetic series with no structure in/// them at all. This applies the same edge to the real thing: the sample/// correlation matrix of daily changes in the observed curve, exactly as/// [`curve_factors`] builds it, and asks how many of its eigenvalues clear the/// edge a structureless panel of this shape would produce.////// Returns the sample eigenvalues (largest first) and the Marchenko-Pastur/// upper edge for this panel's own ratio of tenors to observations.pub fn curve_spectrum_against_noise(observations: &[Vec<f64>]) -> (Vec<f64>, f64) { let n = observations.len(); let p = observations[0].len(); let corr = correlation_of_changes(observations); let eigenvalues = crate::estimation::symmetric_eigenvalues(&corr); let q = p as f64 / (n - 1) as f64; let edge = (1.0 + q.sqrt()).powi(2); (eigenvalues, edge)} /// The sample correlation matrix of daily changes across tenors, exactly as/// [`curve_spectrum_against_noise`] builds it internally, exposed so a/// structural form can be fitted against the same matrix the noise edge is/// tested against.pub fn correlation_of_changes(observations: &[Vec<f64>]) -> Vec<Vec<f64>> { let n = observations.len(); assert!(n > 2, "a covariance needs observations"); let p = observations[0].len(); let changes: Vec<Vec<f64>> = (1..n) .map(|i| (0..p).map(|j| observations[i][j] - observations[i - 1][j]).collect()) .collect(); let count = changes.len() as f64; let means: Vec<f64> = (0..p).map(|j| changes.iter().map(|c| c[j]).sum::<f64>() / count).collect(); let sds: Vec<f64> = (0..p) .map(|j| (changes.iter().map(|c| (c[j] - means[j]).powi(2)).sum::<f64>() / count).sqrt()) .collect(); let mut corr = vec![vec![0.0; p]; p]; for c in &changes { for i in 0..p { for j in 0..p { corr[i][j] += (c[i] - means[i]) * (c[j] - means[j]) / (sds[i] * sds[j] * count); } } } corr} /// Fit the two-parameter correlation form `rho_inf + (1 - rho_inf) * exp(-beta/// |T_i - T_j|)` to a sample correlation matrix's off-diagonal entries, by/// coordinate search on the sum of squared residuals.////// This is the estimation chapter's question turned into code: does the/// market-model chapter's structural form actually match the shape a/// principal component analysis finds, or only its rank?pub fn fit_maturity_correlation(corr: &[Vec<f64>], maturities: &[f64]) -> (f64, f64) { let p = corr.len(); let sse = |rho_inf: f64, beta: f64| -> f64 { let mut total = 0.0; for i in 0..p { for j in 0..p { if i == j { continue; } let model = rho_inf + (1.0 - rho_inf) * (-beta * (maturities[i] - maturities[j]).abs()).exp(); total += (corr[i][j] - model).powi(2); } } total }; let search = |value: f64, candidates: &[f64], f: &dyn Fn(f64) -> f64| -> f64 { let mut best = value; let mut best_sse = f(value); for &c in candidates { let s = f(c); if s < best_sse { best_sse = s; best = c; } } best }; let (mut rho_inf, mut beta) = (0.5f64, 0.2f64); for sweep in 0..6 { let scale = 0.5f64.powi(sweep); let rho_grid: Vec<f64> = (1..=10).map(|i| (rho_inf + (i as f64 - 5.5) * 0.15 * scale).clamp(0.0, 0.98)).collect(); rho_inf = search(rho_inf, &rho_grid, &|r| sse(r, beta)); let beta_grid: Vec<f64> = (1..=10).map(|i| (beta * (0.3 + 0.15 * i as f64 * scale)).max(0.01)).collect(); beta = search(beta, &beta_grid, &|b| sse(rho_inf, b)); } (rho_inf, beta)} /// The kernel a rank-truncated PCA gives, built from exactly the pieces/// [`curve_factors`] already computes.////// This is the curve-construction chapter's "infinitely strong prior on a/// subspace" made concrete on real data: keep `rank` factors --- chosen by/// [`curve_spectrum_against_noise`], not guessed --- and the returned/// covariance is exactly that rank, with zero variance in every direction the/// kept factors do not span. In the Karhunen-Loeve form,/// `K = sum_r sizes[r]^2 * loadings[r] loadings[r]^T`, truncated to `rank`/// terms; taking `rank` equal to the full tenor count reproduces the sample/// covariance exactly, which/// `factor_tests::the_full_rank_kernel_reproduces_the_sample_covariance`/// checks.pub fn empirical_kernel(observations: &[Vec<f64>], rank: usize) -> Vec<Vec<f64>> { let f = curve_factors(observations, rank); let m = f.loadings[0].len(); let mut k = vec![vec![0.0; m]; m]; for r in 0..f.sizes.len() { let var = f.sizes[r] * f.sizes[r]; for i in 0..m { for j in 0..m { k[i][j] += var * f.loadings[r][i] * f.loadings[r][j]; } } } k} #[cfg(test)]mod factor_tests { use super::*; /// The committed Treasury history, as a matrix of rates. fn history() -> (Vec<f64>, Vec<Vec<f64>>) { history_from(concat!(env!("CARGO_MANIFEST_DIR"), "/../public/marketdata/treasury-history.json")) } /// The long multi-tenor panel: eight tenors common to all of them back to /// 1981, spanning several rate regimes rather than one trading year. fn history_long() -> (Vec<f64>, Vec<Vec<f64>>) { history_from(concat!( env!("CARGO_MANIFEST_DIR"), "/../public/marketdata/treasury-history-long.json" )) } fn history_from(path: &str) -> (Vec<f64>, Vec<Vec<f64>>) { let raw = std::fs::read_to_string(path).expect("run `npm run marketdata`"); // Tenor labels, in the order the file lists them. let tenor_block = &raw[raw.find("\"tenors\"").unwrap()..raw.find("\"observations\"").unwrap()]; let labels: Vec<String> = tenor_block .match_indices("\"label\":") .map(|(i, _)| { let rest = &tenor_block[i + 8..]; let a = rest.find('"').unwrap() + 1; let b = rest[a..].find('"').unwrap(); rest[a..a + b].to_string() }) .collect(); let years: Vec<f64> = tenor_block .match_indices("\"years\":") .map(|(i, _)| { let rest = &tenor_block[i + 8..]; let end = rest.find(|c| c == ',' || c == '}').unwrap(); rest[..end].trim().parse().unwrap() }) .collect(); let obs_block = &raw[raw.find("\"observations\"").unwrap()..]; let mut rows = Vec::new(); for (i, _) in obs_block.match_indices("\"rates\"") { let rest = &obs_block[i..]; let end = rest.find('}').unwrap(); let seg = &rest[..end]; let mut row = Vec::new(); for label in &labels { let key = format!("\"{label}\":"); match seg.find(&key) { Some(k) => { let tail = &seg[k + key.len()..]; let stop = tail.find(|c| c == ',' || c == '\n').unwrap_or(tail.len()); row.push(tail[..stop].trim().parse::<f64>().unwrap_or(f64::NAN)); } None => row.push(f64::NAN), } } if row.iter().all(|x| x.is_finite()) { rows.push(row); } } (years, rows) } #[test] fn the_curve_has_a_level_a_slope_and_a_curvature() { // The standard finding, on this year's data rather than on assertion. let (years, rows) = history(); assert!(rows.len() > 100, "expected a year of curves, got {}", rows.len()); let f = curve_factors(&rows, 3); // The first factor explains most of the variance and loads with one sign // everywhere: every point of the curve moves together. That is the level. assert!(f.shares[0] > 0.7, "level explains {:.3}", f.shares[0]); assert!( f.loadings[0].iter().all(|x| *x > 0.0), "the level factor should load positively at every tenor" ); // The second changes sign exactly once across the curve: the short end and // the long end move oppositely. That is the slope. let crossings = f.loadings[1] .windows(2) .filter(|w| w[0] * w[1] < 0.0) .count(); assert_eq!(crossings, 1, "the slope factor should change sign once"); // The third changes sign twice: the wings against the middle. let curvature = f.loadings[2].windows(2).filter(|w| w[0] * w[1] < 0.0).count(); assert_eq!(curvature, 2, "the curvature factor should change sign twice"); assert_eq!(years.len(), f.loadings[0].len()); } #[test] fn the_slope_is_the_smaller_factor_but_not_negligible() { // The quantity the risk chapter's explain rests on. The slope carries far // less variance than the level and is nowhere near zero, which is exactly // the combination that makes a level-hedged book's profit and loss // unexplainable by a one-factor report. let (_, rows) = history(); let f = curve_factors(&rows, 3); let ratio = f.sizes[1] / f.sizes[0]; assert!( (0.1..0.6).contains(&ratio), "the slope should be a minority of the level, at {ratio:.3}" ); // Annualised, both are basis-point quantities a reader can hold on to. let annual = |daily: f64| daily * (252.0f64).sqrt() * 10_000.0; assert!( (30.0..200.0).contains(&annual(f.sizes[0])), "level {:.0}bp a year", annual(f.sizes[0]) ); } /// The bond-curve tenors \chapref{curve-construction}'s own example uses, /// one year and beyond, dropping the money-market front end. fn bond_curve_panel() -> Vec<Vec<f64>> { let (years, rows) = history(); let keep: Vec<usize> = years.iter().enumerate().filter(|(_, &y)| y >= 1.0).map(|(i, _)| i).collect(); rows.iter().map(|r| keep.iter().map(|&i| r[i]).collect()).collect() } #[test] fn only_level_and_slope_clear_the_noise_edge_on_the_full_panel() { // The estimation chapter's Marchenko-Pastur point, run on the real // panel rather than on simulated noise, on the same thirteen tenors // curve_factors' own tests use. let (_, rows) = history(); let (eig, edge) = curve_spectrum_against_noise(&rows); let real = eig.iter().filter(|&&x| x > edge).count(); assert_eq!(real, 2, "eigenvalues {eig:?} against edge {edge:.3}"); } #[test] fn only_level_and_slope_clear_the_noise_edge_across_four_decades() { // The regime question. curve_spectrum_against_noise found the same // count on eight months of one year; here it is asked again on // forty-four years and several rate cycles, so the finding is not an // artefact of the one window that happens to be committed. let (_, rows) = history_long(); assert!(rows.len() > 10_000, "expected decades of daily observations, got {}", rows.len()); let (eig, edge) = curve_spectrum_against_noise(&rows); let real = eig.iter().filter(|&&x| x > edge).count(); assert_eq!(real, 2, "eigenvalues {eig:?} against edge {edge:.3}"); // The margin matters as much as the count: slope clears by a little, // curvature misses by a lot, so the result is not one noisy day from // flipping. assert!(eig[1] / edge < 1.2, "slope should clear narrowly: {:.3} against {:.3}", eig[1], edge); assert!(eig[2] / edge < 0.5, "curvature should miss by a wide margin: {:.3} against {:.3}", eig[2], edge); } #[test] fn the_two_parameter_correlation_form_gets_the_rank_right_but_understates_the_slope() { // Does the market-model chapter's structural form match the shape a // principal component analysis finds on this book's own curve, or // only its rank? Fit it to the same correlation matrix the noise-edge // test above uses, then compare the fitted matrix's own eigenvalues // against the real ones. let (years, rows) = history_long(); let corr = correlation_of_changes(&rows); let (eig, edge) = curve_spectrum_against_noise(&rows); let (rho_inf, beta) = fit_maturity_correlation(&corr, &years); assert!(rho_inf > 0.3 && rho_inf < 0.7, "rho_inf {rho_inf:.3} should be a genuine positive floor"); assert!(beta > 0.05 && beta < 0.5, "beta {beta:.3} should be a genuine decay rate"); let p = corr.len(); let mut fitted = vec![vec![0.0; p]; p]; for i in 0..p { for j in 0..p { fitted[i][j] = if i == j { 1.0 } else { rho_inf + (1.0 - rho_inf) * (-beta * (years[i] - years[j]).abs()).exp() }; } } let fitted_eig = crate::estimation::symmetric_eigenvalues(&fitted); // The level eigenvalue is matched closely. assert!( (fitted_eig[0] / eig[0] - 1.0).abs() < 0.05, "fitted level {:.3} against real {:.3}", fitted_eig[0], eig[0] ); // The slope eigenvalue is not: the fitted form understates it, because // a single exponential decay is too rigid to match both eigenvalues' // sizes at once. Pinned rather than asserted, so the size does not // drift silently. let slope_ratio = fitted_eig[1] / eig[1]; assert!( slope_ratio > 0.6 && slope_ratio < 0.8, "fitted slope {:.3} against real {:.3}, ratio {:.3}", fitted_eig[1], eig[1], slope_ratio ); // Neither matrix manufactures a third real factor: the fitted form's // own remaining eigenvalues stay below the noise edge, matching the // rank the real data supports. assert!( fitted_eig[2] < edge, "fitted third eigenvalue {:.3} should stay below the edge {:.3}", fitted_eig[2], edge ); } #[test] fn only_the_level_clears_the_noise_edge_on_the_bond_curve_panel() { // The narrower tenor set \chapref{curve-construction} actually fits a // curve to. Fewer, more homogeneous series concentrate more of the // variance into the first factor and leave less for a second to clear // the edge with. let narrow = bond_curve_panel(); let (eig, edge) = curve_spectrum_against_noise(&narrow); let real = eig.iter().filter(|&&x| x > edge).count(); assert_eq!(real, 1, "eigenvalues {eig:?} against edge {edge:.3}"); } #[test] fn the_full_rank_kernel_reproduces_the_sample_covariance() { // Keeping every factor is not a truncation at all, so the // Karhunen-Loeve sum has to add back up to the covariance it was // decomposed from. let (_, rows) = history(); let m = rows[0].len(); let n = rows.len(); let k = empirical_kernel(&rows, m); let changes: Vec<Vec<f64>> = (1..n).map(|i| (0..m).map(|j| rows[i][j] - rows[i - 1][j]).collect()).collect(); let count = changes.len() as f64; let means: Vec<f64> = (0..m).map(|j| changes.iter().map(|c| c[j]).sum::<f64>() / count).collect(); let mut cov = vec![vec![0.0; m]; m]; for c in &changes { for i in 0..m { for j in 0..m { cov[i][j] += (c[i] - means[i]) * (c[j] - means[j]) / count; } } } for i in 0..m { for j in 0..m { assert!( (k[i][j] - cov[i][j]).abs() < 1e-8, "({i},{j}): kernel {} against covariance {}", k[i][j], cov[i][j] ); } } } #[test] fn a_rank_truncated_kernel_has_exactly_that_rank() { // The point of truncating at all: the discarded directions are not // merely small, they are gone. let narrow = bond_curve_panel(); let k = empirical_kernel(&narrow, 1); let eig = crate::estimation::symmetric_eigenvalues(&k); let nonzero = eig.iter().filter(|&&x| x > 1e-10).count(); assert_eq!(nonzero, 1, "eigenvalues of the rank-1 kernel: {eig:?}"); assert!(eig.iter().all(|&x| x >= -1e-10), "kernel should be positive semidefinite: {eig:?}"); } #[test] fn the_butterfly_is_where_the_third_factor_lives() { // The chapter's claim about curve trades, on this year's curve. The // 1:-2:1 weighting of 2s, 5s and 10s is meant to be neutral to the level // and the slope, leaving curvature; the comparison that gives the number // meaning is the outright, which is almost all level. let (years, rows) = history(); let find = |t: f64| years.iter().position(|y| (y - t).abs() < 1e-9).expect("tenor"); let n = factor_neutrality(&rows, (find(2.0), find(5.0), find(10.0)), 3); assert!( n.outright_explained > 0.9, "an outright should be almost entirely the first two factors, got {:.3}", n.outright_explained ); // Most of it, but not all: the equal-weight fly is only approximately // factor-neutral, and how far from neutral is the measurement the // chapter's remark about the weights being a modelling choice needs. assert!( n.butterfly_explained < 0.4, "the fly should have weighted most of them out, got {:.3}", n.butterfly_explained ); assert!( n.butterfly_explained < n.outright_explained / 2.0, "the fly should be far less factor-driven than the outright" ); // And the residual is not zero: a fly that moved not at all would be no // trade either. assert!( n.butterfly_volatility > 0.0, "the fly has to move to be tradeable" ); println!( "fly explained {:.3}, outright {:.3}, fly vol {:.4} per day", n.butterfly_explained, n.outright_explained, n.butterfly_volatility ); } #[test] fn the_condor_is_genuinely_neutral_and_its_own_return_is_unmeasurable() { // The butterfly's third factor is one factor short of neutral, per the // test above. A fourth point removes curvature too, and what is left // is exactly what the chapter's own theme about drift says it should // be: too little history to tell whether it has a return at all. let (years, rows) = history(); let find = |t: f64| years.iter().position(|y| (y - t).abs() < 1e-9).expect("tenor"); let columns = [find(2.0), find(5.0), find(10.0), find(30.0)]; let last = rows.last().unwrap(); let durations = [ par_duration(last[columns[0]], 2.0), par_duration(last[columns[1]], 5.0), par_duration(last[columns[2]], 10.0), par_duration(last[columns[3]], 30.0), ]; let cd = condor_neutrality(&rows, columns, durations); // Genuinely neutral, not merely approximately: the fourth constraint // removes what the butterfly's 26.4 per cent left behind. assert!( cd.explained < 0.05, "the condor should be far more neutral than the fly, got {:.3}", cd.explained ); // The point worth measuring: does this residual carry a return? A // t-statistic on one year of daily data, computed rather than assumed. let n = (rows.len() - 1) as f64; let standard_error = cd.daily_volatility / n.sqrt(); let t_stat = cd.daily_mean / standard_error; assert!( t_stat.abs() < 2.0, "a year of daily data should not be able to tell this mean from zero, got t={:.2}", t_stat ); println!( "condor explained {:.4}, mean {:.6}, vol {:.4}, t-stat {:.3}, weights {:?}, normalised {:?}", cd.explained, cd.daily_mean, cd.daily_volatility, t_stat, cd.weights, cd.weights.map(|w| w / cd.weights[0].abs()) ); }} /// Carry and roll-down against the risk of holding, from a real curve.////// The relative value chapter decomposes a bond's return into carry, roll-down/// and a yield change, and observes that the first two are known today while the/// third is not. It does not say how they compare, and the comparison is what/// decides whether a carry trade is a harvest or a bet.////// `curve` is a set of (maturity, par yield) points, `funding` the repo rate, and/// `yield_vol` the standard deviation of that maturity's yield change over the/// horizon --- measured from history rather than assumed.pub struct CarryTrade { /// Carry plus roll-down over the horizon, as a return. pub known: f64, /// One standard deviation of the yield-change term, as a return. pub uncertain: f64,} impl CarryTrade { /// The ratio of what is known to what is not. Below one, the trade's edge is /// smaller than a single standard deviation of its risk over the same /// horizon. pub fn ratio(&self) -> f64 { self.known / self.uncertain }} /// Duration of a par bond at yield `y` and maturity `t`, annual coupons.fn par_duration(y: f64, t: f64) -> f64 { if y.abs() < 1e-9 { return t; } (1.0 - (1.0 + y).powf(-t)) / y} /// Evaluate carry and roll for one maturity.////// `yield_at` interpolates the curve, so the roll-down term is the actual slope/// between `t` and `t - horizon` rather than a local derivative.pub fn carry_and_roll( yield_at: impl Fn(f64) -> f64, maturity: f64, funding: f64, horizon: f64, yield_vol: f64,) -> CarryTrade { let y = yield_at(maturity); let rolled = yield_at((maturity - horizon).max(0.01)); let duration = par_duration(y, maturity - horizon); let carry = (y - funding) * horizon; let roll = -duration * (rolled - y); CarryTrade { known: carry + roll, uncertain: duration * yield_vol }} #[cfg(test)]mod carry_tests { use super::*; /// An upward sloping curve, so roll-down is positive and the test is not /// about a degenerate shape. fn curve(t: f64) -> f64 { 0.037 + 0.013 * (1.0 - (-0.25 * t).exp()) } #[test] fn roll_down_is_positive_on_an_upward_sloping_curve() { // The sign the chapter claims, and the reason: a bond ages into a lower // yield, so its price rises even if nothing moves. let flat = carry_and_roll(|_| 0.04, 10.0, 0.037, 0.25, 0.008); let sloped = carry_and_roll(curve, 10.0, 0.037, 0.25, 0.008); assert!(sloped.known > flat.known, "slope should add roll-down"); } #[test] fn the_known_part_is_small_against_the_risk() { // The comparison the chapter leaves out. Over a quarter, with a yield // volatility of the size the Treasury panel shows, carry and roll on any // maturity are a fraction of one standard deviation of the yield-change // term. A carry trade is not a harvest with a little noise on it; it is a // bet with a small tilt. for maturity in [2.0, 5.0, 10.0, 30.0] { let t = carry_and_roll(curve, maturity, 0.037, 0.25, 0.008); assert!(t.known > 0.0, "maturity {maturity} should have positive carry and roll"); assert!( t.ratio() < 0.6, "maturity {maturity}: known {:.5} against risk {:.5}, ratio {:.2}", t.known, t.uncertain, t.ratio() ); } } /// The committed Treasury par curve `latest.json` points at, tenors and /// rates in quoted order. fn treasury_snapshot() -> (Vec<f64>, Vec<f64>) { let dir = concat!(env!("CARGO_MANIFEST_DIR"), "/../public/marketdata"); let latest = std::fs::read_to_string(format!("{dir}/latest.json")) .expect("run `npm run marketdata`"); let key = "\"treasuryParYields\":"; let at = latest.find(key).unwrap() + key.len(); let rest = &latest[at..]; let open = rest.find('"').unwrap() + 1; let close = rest[open..].find('"').unwrap() + open; let file = &rest[open..close]; let body = std::fs::read_to_string(format!("{dir}/{file}")).unwrap(); let numbers = |needle: &str| -> Vec<f64> { let pat = format!("\"{needle}\":"); let mut out = Vec::new(); let mut rest = body.as_str(); while let Some(i) = rest.find(&pat) { let tail = &rest[i + pat.len()..]; let end = tail.find(|c: char| c == ',' || c == '}').unwrap(); out.push(tail[..end].trim().parse::<f64>().unwrap()); rest = &tail[end..]; } out }; (numbers("years"), numbers("rate")) } #[test] fn the_real_curve_has_a_local_inversion_past_twenty_years() { // The committed panel is not just upward sloping and flattening: the // twenty year point sits above both ten and thirty, so the curve // genuinely inverts between twenty and thirty years. Pinned because // the next test's surprising result rests entirely on it being real // rather than an artefact of interpolation. let (tenors, rates) = treasury_snapshot(); let at = |t: f64| -> f64 { let i = tenors.iter().position(|&x| (x - t).abs() < 1e-9).unwrap(); rates[i] }; assert!(at(20.0) > at(10.0), "20y {} should exceed 10y {}", at(20.0), at(10.0)); assert!(at(20.0) > at(30.0), "20y {} should exceed 30y {}", at(20.0), at(30.0)); } #[test] fn carry_and_roll_on_the_real_curve_is_not_monotone_and_can_go_negative() { // The figure and the chapter's table were computed by two different // routes and disagreed; this settles it by construction, using the // figure's own exact discount-factor arithmetic rather than the // duration-linearised approximation `carry_and_roll` above uses, // against the same committed panel and the same three month horizon. use crate::curve::{bootstrap_par, Interp}; let (tenors, rates) = treasury_snapshot(); let curve = bootstrap_par(&tenors, &rates, 2.0, Interp::LinearForward); let horizon = 0.25; let excess = |maturity: f64| -> f64 { let hold = curve.df(maturity - horizon) / curve.df(maturity) - 1.0; let fund = 1.0 / curve.df(horizon) - 1.0; (hold - fund) / horizon }; let thirty = excess(30.0); let twenty = excess(20.0); let ten = excess(10.0); assert!(twenty > ten, "twenty year excess {twenty:.5} should exceed ten year {ten:.5}"); assert!( thirty < 0.0, "thirty year excess should be negative on this inverted panel, got {thirty:.5}" ); assert!( thirty < ten, "thirty year excess {thirty:.5} should be worse than ten year {ten:.5}, not merely smaller than twenty's" ); } #[test] fn the_ratio_improves_with_a_shorter_horizon_and_worsens_with_duration() { // Two dependencies worth having. Carry accrues linearly in the horizon // while the risk grows as its square root, so a shorter hold is a worse // ratio; and duration multiplies the risk while adding only roll, so the // long end is the wrong place to look for a carry edge. let quarter = carry_and_roll(curve, 10.0, 0.037, 0.25, 0.008); let year = carry_and_roll(curve, 10.0, 0.037, 1.0, 0.016); assert!(year.ratio() > quarter.ratio(), "a longer hold should improve the ratio"); let short = carry_and_roll(curve, 2.0, 0.037, 0.25, 0.008); let long = carry_and_roll(curve, 30.0, 0.037, 0.25, 0.008); assert!(short.ratio() > long.ratio(), "duration should hurt the ratio"); } /// Daily multi-tenor Treasury history from the named committed file, /// tenors and rates per observation in quoted order. fn treasury_history_from(path: &str) -> (Vec<f64>, Vec<Vec<f64>>) { let raw = std::fs::read_to_string(path).expect("run `npm run marketdata`"); let tenor_block = &raw[raw.find("\"tenors\"").unwrap()..raw.find("\"observations\"").unwrap()]; let labels: Vec<String> = tenor_block .match_indices("\"label\":") .map(|(i, _)| { let rest = &tenor_block[i + 8..]; let a = rest.find('"').unwrap() + 1; let b = rest[a..].find('"').unwrap(); rest[a..a + b].to_string() }) .collect(); let years: Vec<f64> = tenor_block .match_indices("\"years\":") .map(|(i, _)| { let rest = &tenor_block[i + 8..]; let end = rest.find(|c| c == ',' || c == '}').unwrap(); rest[..end].trim().parse().unwrap() }) .collect(); let obs_block = &raw[raw.find("\"observations\"").unwrap()..]; let mut rows = Vec::new(); for (i, _) in obs_block.match_indices("\"rates\"") { let rest = &obs_block[i..]; let end = rest.find('}').unwrap(); let seg = &rest[..end]; let mut row = Vec::new(); for label in &labels { let key = format!("\"{label}\":"); let k = seg.find(&key).unwrap(); let tail = &seg[k + key.len()..]; let stop = tail.find(|c| c == ',' || c == '\n').unwrap_or(tail.len()); row.push(tail[..stop].trim().parse::<f64>().unwrap()); } rows.push(row); } (years, rows) } /// Quarterly yield volatility at a tenor, from daily changes in the /// named history file: sample standard deviation of the change (not the /// raw move, so a sustained trend is not counted as risk), scaled to a /// quarter by the usual square-root-of-time in trading days. fn quarterly_vol_from(path: &str, maturity: f64) -> f64 { let (years, history) = treasury_history_from(path); let i = years.iter().position(|&y| (y - maturity).abs() < 1e-9).unwrap(); let changes: Vec<f64> = history.windows(2).map(|w| w[1][i] - w[0][i]).collect(); let mean = changes.iter().sum::<f64>() / changes.len() as f64; let daily_sd = (changes.iter().map(|c| (c - mean).powi(2)).sum::<f64>() / changes.len() as f64) .sqrt(); daily_sd * (63.0f64).sqrt() } #[test] fn what_the_curve_currently_pays_to_hold_by_maturity() { // Calculation "How the known part compares with the unknown" as // measured today, replacing whatever an earlier snapshot produced. // The known part is the figure's own exact discount-factor // arithmetic, not the duration-linearised approximation // `carry_and_roll` above uses, since that approximation is exactly // what breaks down against the real curve's local inversion. The // risk side keeps the duration approximation as a multiplier, but on // volatility measured from four decades of history where that // history exists (one, two, five and ten years); twenty and thirty // year Treasuries were not issued continuously over that span, so // those two fall back to the year the committed multi-tenor panel // actually has. use crate::curve::{bootstrap_par, Interp}; let (tenors, rates) = treasury_snapshot(); let curve = bootstrap_par(&tenors, &rates, 2.0, Interp::LinearForward); let horizon = 0.25; let long = concat!( env!("CARGO_MANIFEST_DIR"), "/../public/marketdata/treasury-history-long.json" ); let short = concat!( env!("CARGO_MANIFEST_DIR"), "/../public/marketdata/treasury-history.json" ); let known = |maturity: f64| -> f64 { let hold = curve.df(maturity - horizon) / curve.df(maturity) - 1.0; let fund = 1.0 / curve.df(horizon) - 1.0; (hold - fund) / horizon }; let mut rows = Vec::new(); for &maturity in &[1.0, 2.0, 5.0, 10.0, 30.0] { let par_yield = { let i = tenors.iter().position(|&t| (t - maturity).abs() < 1e-9).unwrap(); rates[i] }; let history_used = if maturity <= 10.0 { long } else { short }; let duration = par_duration(par_yield, maturity - horizon); let uncertain = duration * quarterly_vol_from(history_used, maturity); rows.push((maturity, par_yield, known(maturity), uncertain)); } // One year is a genuine, currently measured exception to "every // ratio is well below one": the front end of this snapshot's curve // is unusually steep, and even four decades of yield volatility does // not dominate the carry and roll it offers. let one_year_ratio = rows[0].2 / rows[0].3; assert!( one_year_ratio > 1.5, "one year ratio {one_year_ratio:.3} should still be a clear exception" ); // From two years out the original picture holds: known is a small // fraction of the risk. for row in &rows[1..4] { let ratio = (row.2 / row.3).abs(); assert!(ratio < 0.6, "maturity {}: ratio {ratio:.3} should be well below one", row.0); } // The finding worth pinning: thirty years is not the best place on // the curve to hold, and its known return is comfortably negative, // on a curve that inverts between twenty and thirty years. let thirty = rows.last().unwrap(); let best_known = rows.iter().map(|r| r.2).fold(f64::NEG_INFINITY, f64::max); assert!( thirty.2 < best_known, "thirty year known return {:.5} should not be the best on the curve", thirty.2 ); assert!(thirty.2 < -0.01, "thirty year known return {:.5} should be clearly negative", thirty.2); } } /// What a butterfly buys you, measured rather than asserted.////// The relative value chapter claims that a butterfly is a trade in the third/// factor because the first two are weighted out of it. That is a claim about/// this year's curve and it can be checked on it.////// Take the `1 : -2 : 1` weighting of three tenors and regress its daily changes/// on the leading principal components of the whole curve. The share of variance/// the factors explain is the exposure the weighting failed to remove. The same/// regression on an outright position in the belly tenor is the comparison that/// makes the number mean something: an outright is almost entirely the level/// factor, and if the butterfly were not different there would be no trade.pub struct FactorNeutrality { /// Share of the butterfly's daily variance explained by the first two factors. pub butterfly_explained: f64, /// The same for an outright position in the middle tenor. pub outright_explained: f64, /// Daily standard deviation of the butterfly, in the units of the input. pub butterfly_volatility: f64,} /// Regress the `1 : -2 : 1` fly on the leading factors. `columns` are indices/// into each observation, short to long.pub fn factor_neutrality( observations: &[Vec<f64>], columns: (usize, usize, usize), factors: usize,) -> FactorNeutrality { let f = curve_factors(observations, factors); let (a, b, c) = columns; // Daily changes of the two positions, and of every tenor, on the same days. let changes: Vec<Vec<f64>> = (1..observations.len()) .map(|i| { (0..observations[0].len()) .map(|j| observations[i][j] - observations[i - 1][j]) .collect() }) .collect(); // A position's exposure to a factor is the loading contracted with its // weights, so the explained share is that contracted variance over the total. let explained = |weights: &[(usize, f64)]| { let series: Vec<f64> = changes .iter() .map(|d| weights.iter().map(|(j, w)| w * d[*j]).sum::<f64>()) .collect(); let mean = series.iter().sum::<f64>() / series.len() as f64; let total: f64 = series.iter().map(|x| (x - mean) * (x - mean)).sum::<f64>() / series.len() as f64; // Project each day's change onto the factors, rebuild the position from // the projection, and see how much of its variance survives. let mut from_factors = 0.0; for k in 0..factors.min(2) { let load: f64 = weights.iter().map(|(j, w)| w * f.loadings[k][*j]).sum(); from_factors += load * load * f.sizes[k] * f.sizes[k]; } (from_factors / total, total.sqrt()) }; let (butterfly_explained, butterfly_volatility) = explained(&[(a, 1.0), (b, -2.0), (c, 1.0)]); let (outright_explained, _) = explained(&[(b, 1.0)]); FactorNeutrality { butterfly_explained, outright_explained, butterfly_volatility }} /// What is left over once a curve position is null to level, slope \emph{and}/// curvature: four points instead of three, one more constraint than the/// butterfly imposes.pub struct Condor { /// Notional weights on the four points, in the order given. pub weights: [f64; 4], /// Share of daily variance the first three factors still explain --- the /// three-factor analogue of `FactorNeutrality::butterfly_explained`. pub explained: f64, /// Daily mean and standard deviation of the position's own P\&L, in the /// input's units. pub daily_mean: f64, pub daily_volatility: f64,} /// The unique (up to scale) weighting of four curve points with zero net/// loading on level, slope and curvature at once, found as the null vector of/// the $3\times4$ matrix of duration-weighted loadings --- the generalised/// cross product, the same construction a normal to three vectors in four/// dimensions always has.pub fn condor_neutrality(observations: &[Vec<f64>], columns: [usize; 4], durations: [f64; 4]) -> Condor { let f = curve_factors(observations, 3); let row = |k: usize| -> [f64; 4] { let mut r = [0.0; 4]; for i in 0..4 { r[i] = durations[i] * f.loadings[k][columns[i]]; } r }; let (level, slope, curvature) = (row(0), row(1), row(2)); let det3 = |m: [[f64; 3]; 3]| -> f64 { m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]) }; let minor = |skip: usize| -> f64 { let cols: Vec<usize> = (0..4).filter(|&j| j != skip).collect(); det3([ [level[cols[0]], level[cols[1]], level[cols[2]]], [slope[cols[0]], slope[cols[1]], slope[cols[2]]], [curvature[cols[0]], curvature[cols[1]], curvature[cols[2]]], ]) }; let weights = [minor(0), -minor(1), minor(2), -minor(3)]; // Daily P&L of the weighted position, in the same duration-weighted units // the null vector was solved in --- a notional $w_i$ at duration $D_i$ // moves by $w_i D_i$ per unit yield change --- and that same series // projected onto the three factors, exactly as `factor_neutrality` // compares them. let changes: Vec<Vec<f64>> = (1..observations.len()) .map(|i| { (0..observations[0].len()) .map(|j| observations[i][j] - observations[i - 1][j]) .collect() }) .collect(); let series: Vec<f64> = changes .iter() .map(|d| (0..4).map(|i| weights[i] * durations[i] * d[columns[i]]).sum::<f64>()) .collect(); let n = series.len() as f64; let mean = series.iter().sum::<f64>() / n; let total: f64 = series.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / n; let mut from_factors = 0.0; for k in 0..3 { let load: f64 = (0..4).map(|i| weights[i] * durations[i] * f.loadings[k][columns[i]]).sum(); from_factors += load * load * f.sizes[k] * f.sizes[k]; } Condor { weights, explained: from_factors / total, daily_mean: mean, daily_volatility: total.sqrt(), }}