quant/src/credit.rs
Credit default swaps, and the hazard curve behind them.
//! Credit default swaps, and the hazard curve behind them.//!//! The price process chapter built the default time as the first jump of a Cox process and//! observed that the survival probability//!//! ```text//! Q(t) = E[ exp(-integral of lambda) ]//! ```//!//! is formally the bond price of the curve construction chapter with the//! intensity in place of the short rate. The credit chapter cashes that//! observation in: this module prices the instrument that trades on it, and//! bootstraps a hazard curve from quotes by the same procedure the curve//! construction chapter uses on par swaps.//!//! Rates are taken flat here. Nothing depends on it --- the derivation in the//! chapter shows the discounting cancels between the two legs of a swap --- and//! it keeps the code about credit rather than about a second curve. /// A piecewise constant hazard rate curve.////// Piecewise constant in the *forward* hazard, which is the credit analogue of/// the curve construction chapter's piecewise constant instantaneous forward/// rate, and is what makes the bootstrap triangular: each new quote fixes the/// hazard on one new segment and disturbs nothing before it.#[derive(Clone, Debug)]pub struct CreditCurve { /// Segment end times, increasing. times: Vec<f64>, /// Hazard rate on each segment. hazards: Vec<f64>,} impl CreditCurve { pub fn new(times: Vec<f64>, hazards: Vec<f64>) -> Self { assert_eq!(times.len(), hazards.len()); CreditCurve { times, hazards } } /// Flat curve, which is the case the credit triangle is exact for. pub fn flat(hazard: f64, horizon: f64) -> Self { CreditCurve::new(vec![horizon], vec![hazard]) } /// The hazard rate in force at `t`. pub fn hazard(&self, t: f64) -> f64 { let i = self.times.partition_point(|&x| x < t).min(self.hazards.len() - 1); self.hazards[i] } /// The survival probability `Q(t) = exp(-integral of the hazard)`. pub fn survival(&self, t: f64) -> f64 { if t <= 0.0 { return 1.0; } let mut accumulated = 0.0; let mut previous = 0.0; for (&end, &h) in self.times.iter().zip(&self.hazards) { let upper = end.min(t); if upper > previous { accumulated += h * (upper - previous); previous = upper; } if end >= t { break; } } // Beyond the last node, continue at the final hazard. if t > previous { accumulated += self.hazards[self.hazards.len() - 1] * (t - previous); } (-accumulated).exp() }} /// The two legs of a credit default swap, per unit notional and per unit of/// spread on the premium leg.////// `freq` is the number of premium payments a year; the market convention is/// four. `steps` sets the grid the protection leg is integrated on, which needs/// to be finer than the payment schedule because default can happen any day.pub struct CdsLegs { /// Present value of receiving one unit of loss on default. pub protection: f64, /// Present value of paying one unit of spread until default or maturity, /// including the fraction of a coupon accrued when default falls mid-period. /// The market calls this the risky annuity, and it is the credit analogue of /// the curve construction chapter's annuity. pub risky_annuity: f64,} pub fn cds_legs( curve: &CreditCurve, maturity: f64, recovery: f64, rate: f64, freq: f64, steps: usize,) -> CdsLegs { let df = |t: f64| (-rate * t).exp(); // Protection: the loss, times the probability of defaulting in each small // interval, discounted. Integrated on a fine grid rather than at the payment // dates, since default does not wait for a coupon. let dt = maturity / steps as f64; let mut protection = 0.0; for i in 0..steps { let (a, b) = (i as f64 * dt, (i + 1) as f64 * dt); let default_prob = curve.survival(a) - curve.survival(b); protection += (1.0 - recovery) * df(0.5 * (a + b)) * default_prob; } // Premium: a coupon on each payment date if still alive, plus the accrued // fraction if default happens partway through a period. Half a period is the // standard approximation for where in the period it lands. let accrual = 1.0 / freq; let payments = (maturity * freq).round() as usize; let mut risky_annuity = 0.0; for i in 1..=payments { let t = i as f64 * accrual; let previous = t - accrual; risky_annuity += accrual * df(t) * curve.survival(t); risky_annuity += 0.5 * accrual * df(t) * (curve.survival(previous) - curve.survival(t)); } CdsLegs { protection, risky_annuity }} /// The spread at which a credit default swap is worth zero.pub fn par_spread( curve: &CreditCurve, maturity: f64, recovery: f64, rate: f64, freq: f64, steps: usize,) -> f64 { let legs = cds_legs(curve, maturity, recovery, rate, freq, steps); legs.protection / legs.risky_annuity} /// The credit triangle: the hazard rate implied by a spread, `lambda = s/(1-R)`.////// The credit chapter shows this is not a rule of thumb but exact, for a flat hazard and/// continuously paid premium, whatever the interest rate does --- the/// discounting divides out between the two legs. Against a real contract paying/// every `1/freq` of a year the error is about `rate / (2 * freq)`, the interest/// on half a payment period, and does not depend on the spread. See the tests.pub fn triangle_hazard(spread: f64, recovery: f64) -> f64 { spread / (1.0 - recovery)} /// Bootstrap a piecewise constant hazard curve from par spreads.////// The same procedure as the curve construction chapter's curve bootstrap and/// for the same reason: each quote, given everything shorter, involves one new/// unknown. Solved by bisection on the new segment's hazard, since the par/// spread is increasing in it and no derivative is needed.pub fn bootstrap_hazards( tenors: &[f64], spreads: &[f64], recovery: f64, rate: f64, freq: f64,) -> CreditCurve { assert_eq!(tenors.len(), spreads.len()); let mut curve = CreditCurve::new(Vec::new(), Vec::new()); for (k, (&t, &s)) in tenors.iter().zip(spreads).enumerate() { curve.times.push(t); // Start from the triangle, which is the answer when the curve is flat // and a good starting bracket when it is not. curve.hazards.push(triangle_hazard(s, recovery)); let (mut lo, mut hi) = (1e-8, 5.0); for _ in 0..200 { let mid = 0.5 * (lo + hi); curve.hazards[k] = mid; // More steps at longer maturities, so the protection integral stays // accurate as the horizon grows. let steps = (t * 200.0).ceil() as usize; if par_spread(&curve, t, recovery, rate, freq, steps) < s { lo = mid; } else { hi = mid; } if hi - lo < 1e-14 { break; } } curve.hazards[k] = 0.5 * (lo + hi); } curve} #[cfg(test)]mod tests { use super::*; const R: f64 = 0.4; const TENORS: [f64; 5] = [1.0, 3.0, 5.0, 7.0, 10.0]; const SPREADS: [f64; 5] = [0.0080, 0.0105, 0.0125, 0.0138, 0.0150]; #[test] fn survival_is_a_decreasing_probability() { let c = bootstrap_hazards(&TENORS, &SPREADS, R, 0.03, 4.0); assert!((c.survival(0.0) - 1.0).abs() < 1e-12); let mut previous = 1.0; for i in 1..=100 { let q = c.survival(i as f64 * 0.1); assert!(q > 0.0 && q <= 1.0, "survival was {q}"); assert!(q < previous, "survival rose"); previous = q; } } #[test] fn the_credit_triangle_is_exact_in_the_continuous_limit() { // The credit chapter's derivation: with premium paid continuously and a flat // hazard, the discounting cancels between the two legs and the par // spread is exactly lambda times the loss given default -- for any // interest rate, which is the part that surprises. for rate in [0.0, 0.03, 0.08] { for lambda in [0.01, 0.05, 0.20] { let c = CreditCurve::flat(lambda, 30.0); // 2000 payments a year stands in for continuous. let s = par_spread(&c, 5.0, R, rate, 2000.0, 200_000); let triangle = lambda * (1.0 - R); assert!( (s - triangle).abs() / triangle < 2e-3, "rate={rate} lambda={lambda}: spread {s} against {triangle}" ); } } } #[test] fn discrete_payment_costs_half_a_period_of_interest_and_nothing_else() { // The credit chapter's equation for the convention error. Paying at period end // rather than continuously delays the premium by half a period on // average; the accrued-interest convention compensates the survival part // of that delay, and the interest on it is what is left over. // // Checked at the market's quarterly convention, where the first order // term is the whole story to within a twentieth of a percent. for rate in [0.0, 0.03, 0.08] { for lambda in [0.005, 0.02, 0.20] { let c = CreditCurve::flat(lambda, 30.0); let quoted = par_spread(&c, 5.0, R, rate, 4.0, 50_000); let error = (quoted - lambda * (1.0 - R)) / (lambda * (1.0 - R)); let predicted = rate / (2.0 * 4.0); assert!( (error - predicted).abs() < 5e-4, "rate={rate} lambda={lambda}: error {error} against {predicted}" ); } } } #[test] fn the_convention_error_does_not_depend_on_the_spread() { // The part that is usually got backwards. The triangle is not something // that "breaks down for distressed names": the hazard rate cancels // between the two legs whatever it is, so nothing that survives the // cancellation can be a function of it. Thirty basis points and three // thousand carry the same relative error. let relative_error = |lambda: f64| { let c = CreditCurve::flat(lambda, 30.0); let quoted = par_spread(&c, 5.0, R, 0.03, 4.0, 50_000); (quoted - lambda * (1.0 - R)) / (lambda * (1.0 - R)) }; let tight = relative_error(0.005); for lambda in [0.02, 0.0667, 0.20, 0.50] { let wide = relative_error(lambda); assert!( (wide - tight).abs() < 2e-3, "at lambda={lambda} the error was {wide}, against {tight} at 30bp" ); } } #[test] fn annual_payment_of_a_distressed_name_picks_up_a_second_order_term() { // Where the first order law does run out: with a whole year between // payments and a hazard of 20%, half a period is no longer a good // stand-in for where in the period default lands, and the leftover is // the (lambda*delta)^2/12 that a midpoint rule drops. Recorded here so // that the tolerance in the test above is a known quantity rather than a // tuned one. let (rate, lambda, delta) = (0.0, 0.20, 1.0); let c = CreditCurve::flat(lambda, 30.0); let quoted = par_spread(&c, 5.0, R, rate, 1.0 / delta, 50_000); let error = (quoted - lambda * (1.0 - R)) / (lambda * (1.0 - R)); let second_order = -(lambda * delta).powi(2) / 12.0; assert!( (error - second_order).abs() < 1e-4, "error {error} against second order term {second_order}" ); } #[test] fn the_bootstrap_reprices_every_quote() { // The defining property, exactly as for the curve construction // chapter's rate curve. let c = bootstrap_hazards(&TENORS, &SPREADS, R, 0.03, 4.0); for (&t, &s) in TENORS.iter().zip(&SPREADS) { let repriced = par_spread(&c, t, R, 0.03, 4.0, (t * 400.0) as usize); assert!( (repriced - s).abs() < 1e-7, "at {t}y: repriced {repriced} against quoted {s}" ); } } #[test] fn the_worked_example_in_the_chapter_still_holds() { // The credit chapter prints this table and marks a position off it, so the // numbers are pinned here. If the model changes, this fails and the // chapter gets corrected rather than quietly disagreeing with the code. let c = bootstrap_hazards(&TENORS, &SPREADS, R, 0.03, 4.0); let expected_hazards = [133.0, 197.0, 264.0, 295.0, 311.0]; for (&t, &bp) in TENORS.iter().zip(&expected_hazards) { let hazard = c.hazard(t - 0.01) * 1e4; assert!((hazard - bp).abs() < 0.5, "at {t}y the hazard was {hazard}bp, not {bp}bp"); } // The five year risky annuity, and the mark to market of a contract // struck at 100bp against a market at 125bp. let legs = cds_legs(&c, 5.0, R, 0.03, 4.0, 5000); assert!( (legs.risky_annuity - 4.43).abs() < 0.005, "risky annuity {}", legs.risky_annuity ); let mark_to_market = (0.0125 - 0.0100) * legs.risky_annuity; assert!( (mark_to_market - 0.0111).abs() < 5e-5, "mark to market {mark_to_market}" ); // And the claim that credit shortens the annuity by about four percent. let riskless: f64 = (1..=20).map(|i| 0.25 * (-0.03 * i as f64 * 0.25).exp()).sum(); assert!((riskless - 4.63).abs() < 0.005, "riskless annuity {riskless}"); } #[test] fn a_rising_spread_curve_implies_rising_forward_hazards() { // And by more than the spreads rise, for the same reason the curve // construction chapter's forward rates move further than the zero rates // that carry them: a spread is an average of hazards up to its // maturity, and an average moves less than the thing averaged. let c = bootstrap_hazards(&TENORS, &SPREADS, R, 0.03, 4.0); for w in c.hazards.windows(2) { assert!(w[1] > w[0], "forward hazards not increasing: {:?}", c.hazards); } let hazard_range = c.hazards[4] / c.hazards[0]; let spread_range = SPREADS[4] / SPREADS[0]; assert!( hazard_range > spread_range, "hazards spread by {hazard_range}, quotes by only {spread_range}" ); } #[test] fn recovery_and_spread_trade_off_against_each_other() { // A spread pins down lambda times (1-R), not lambda. Assuming a // different recovery moves the implied hazard and leaves the price // unchanged, which is why recovery is conventionally fixed rather than // fitted -- there is nothing in a single spread to fit it with. let s = 0.0125; for r in [0.2, 0.4, 0.6] { let lambda = triangle_hazard(s, r); let c = CreditCurve::flat(lambda, 30.0); let repriced = par_spread(&c, 5.0, r, 0.03, 2000.0, 100_000); assert!((repriced - s).abs() / s < 2e-3, "R={r} gave {repriced}"); } }} /// A square-root intensity, which is the credit chapter's stochastic hazard rate.////// ```text/// d lambda = kappa (theta - lambda) dt + eta sqrt(lambda) dZ/// ```////// The survival probability is `E[exp(-int lambda)]`, which is the same/// expectation as a discount bond with the intensity in place of the short rate,/// so the whole of the term structure machinery transfers. It is affine, so by/// the solvable models chapter it comes out as an exponential-affine function of/// the current intensity with coefficients solving a Riccati pair — and here the/// pair has a closed form.////// Square root rather than Gaussian because an intensity is an arrival rate and/// cannot be negative, which a Hull-White intensity would not respect.#[derive(Clone, Copy, Debug)]pub struct CirIntensity { pub lambda0: f64, pub kappa: f64, pub theta: f64, pub eta: f64,} impl CirIntensity { /// `E[exp(-int_0^t lambda_s ds)]`, in closed form. pub fn survival(&self, t: f64) -> f64 { if t <= 0.0 { return 1.0; } let (k, th, e2) = (self.kappa, self.theta, self.eta * self.eta); if e2 <= 0.0 { // Deterministic limit: the intensity relaxes from lambda0 to theta. let mean = if k.abs() < 1e-12 { self.lambda0 * t } else { th * t + (self.lambda0 - th) * (1.0 - (-k * t).exp()) / k }; return (-mean).exp(); } let h = (k * k + 2.0 * e2).sqrt(); let expht = (h * t).exp(); let denom = 2.0 * h + (k + h) * (expht - 1.0); let b = 2.0 * (expht - 1.0) / denom; let a = (2.0 * h * ((k + h) * t / 2.0).exp() / denom).powf(2.0 * k * th / e2); a * (-b * self.lambda0).exp() } /// The flat par spread this intensity implies over `[0, t]`, by the credit /// triangle applied to the average hazard the survival curve corresponds to. pub fn implied_flat_spread(&self, t: f64, recovery: f64) -> f64 { let average_hazard = -self.survival(t).ln() / t; average_hazard * (1.0 - recovery) }} #[cfg(test)]mod intensity_tests { use super::*; use crate::pathwise::Rng; fn model(eta: f64) -> CirIntensity { CirIntensity { lambda0: 0.02, kappa: 0.5, theta: 0.02, eta } } /// The closed form against a simulation, which shares none of its algebra. #[test] fn the_closed_form_matches_a_simulation() { let m = model(0.12); let t = 5.0; let steps = 2_000; let dt = t / steps as f64; let mut rng = Rng::new(20260820); let paths = 40_000; let mut total = 0.0; for _ in 0..paths { let mut l = m.lambda0; let mut integral = 0.0; for _ in 0..steps { integral += l * dt; // Full truncation, so the intensity cannot go negative. let root = l.max(0.0).sqrt(); l += m.kappa * (m.theta - l) * dt + m.eta * root * dt.sqrt() * rng.next_normal(); } total += (-integral).exp(); } let simulated = total / paths as f64; let exact = m.survival(t); assert!( (simulated - exact).abs() < 0.002 * exact, "simulation {simulated:.6} against closed form {exact:.6}" ); } /// Zero vol-of-intensity returns the deterministic curve. #[test] fn no_randomness_is_the_deterministic_curve() { let m = model(0.0); for &t in &[0.5, 2.0, 10.0] { assert!((m.survival(t) - (-0.02 * t).exp()).abs() < 1e-12); } } /// The credit chapter's convexity claim. Starting at the long-run level so /// that the intensity's mean is flat, randomness still moves the curve: /// survival is `E[exp(-X)]`, convex in `X`, so it exceeds `exp(-E[X])` and /// the implied spread comes out *below* the deterministic one. /// /// The effect is a Jensen gap and grows with the vol of the intensity, which /// is the same shape of argument as every other convexity in these notes. #[test] fn a_random_intensity_lowers_the_spread() { let t = 5.0; let flat = model(0.0).implied_flat_spread(t, 0.4); let mut previous = flat; for &eta in &[0.05, 0.10, 0.20] { let s = model(eta).implied_flat_spread(t, 0.4); assert!(s < previous, "eta {eta}: spread {s:.6} did not fall below {previous:.6}"); previous = s; } // Sized, so the prose can quote it. let widest = model(0.20).implied_flat_spread(t, 0.4); let drop = (flat - widest) * 1e4; assert!( (drop - 4.19).abs() < 0.1, "the spread fell by {drop:.2} basis points, expected about 4.2" ); // Against a base spread of 120 basis points, so a few per cent. assert!((flat * 1e4 - 120.0).abs() < 0.1, "base spread was {:.2}", flat * 1e4); }}