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 }} #[cfg(test)]mod factor_tests { use super::*; /// The committed Treasury history, as a matrix of rates. fn history() -> (Vec<f64>, Vec<Vec<f64>>) { let raw = std::fs::read_to_string(concat!( env!("CARGO_MANIFEST_DIR"), "/../public/marketdata/treasury-history.json" )) .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]) ); } #[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 ); }} /// 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() ); } } #[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"); } } /// 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 }}