quant/src/dependence.rs
Dependence between defaults, and what it does to a tranche.
//! Dependence between defaults, and what it does to a tranche.//!//! The credit chapter ended on the observation that single-name curves pin every//! marginal default distribution and say nothing whatever about the joint one.//! This module is the arithmetic of that gap: what the attainable correlations//! actually are, what tail dependence is and which models have any, and what the//! difference is worth on a tranche of a portfolio.//!//! The portfolio results are all in the large-portfolio limit, where the loss//! fraction converges to the conditional default probability given the common//! factors. That is not a simplification made for convenience — it is the regime//! the instruments live in, an index has a hundred and twenty five names, and it//! removes idiosyncratic noise so that what is left is the dependence, which is//! the only thing under discussion. use crate::black::norm_cdf;use crate::special::{ln_gamma, norm_inv, t_cdf, t_inv}; // ---------------------------------------------------------------------------// What linear correlation can and cannot say.// --------------------------------------------------------------------------- /// The attainable range of linear correlation between two lognormals.////// Given `exp(X)` and `exp(Y)` with `X, Y` normal of volatilities `s1, s2`, the/// correlation of the *exponentials* cannot reach one even when `X` and `Y` are/// the same random variable. The Hoeffding-Frechet bounds put it at////// ```text/// rho_max = (exp( s1*s2) - 1) / sqrt((exp(s1^2)-1)(exp(s2^2)-1))/// rho_min = (exp(-s1*s2) - 1) / sqrt((exp(s1^2)-1)(exp(s2^2)-1))/// ```////// The dependence chapter makes this the first argument against thinking in/// correlations: a Pearson correlation of 0.9 between two volatile lognormals is/// a number that no joint distribution can produce.pub fn lognormal_correlation_bounds(s1: f64, s2: f64) -> (f64, f64) { let scale = lognormal_scale(s1, s2); if scale <= 0.0 { return (-1.0, 1.0); } ( ((-s1 * s2).exp() - 1.0) / scale, ((s1 * s2).exp() - 1.0) / scale, )} /// `sqrt((exp(s1^2) - 1)(exp(s2^2) - 1))`, the denominator of every lognormal/// correlation in this module.fn lognormal_scale(s1: f64, s2: f64) -> f64 { (((s1 * s1).exp() - 1.0) * ((s2 * s2).exp() - 1.0)).sqrt()} /// The Pearson correlation of two lognormals coupled by a Gaussian copula with/// parameter `rho`.////// ```text/// Corr = (exp(rho s1 s2) - 1) / sqrt((exp(s1^2)-1)(exp(s2^2)-1))/// ```////// Increasing in `rho`, and equal to [`lognormal_correlation_bounds`] at/// `rho = -1` and `rho = 1`. The copula parameter is always attainable; the/// Pearson correlation it produces is not `rho`.pub fn gaussian_copula_pearson(rho: f64, s1: f64, s2: f64) -> f64 { let scale = lognormal_scale(s1, s2); if scale <= 0.0 { return rho; } ((rho * s1 * s2).exp() - 1.0) / scale} /// The Gaussian copula parameter that produces a target Pearson correlation/// between two lognormals, or `None` when no joint distribution can.////// This is [`gaussian_copula_pearson`] inverted, and it fails exactly where the/// attainable range ends: the logarithm has no argument once the target is below/// the lower bound, and the parameter leaves `[-1, 1]` once it is above the upper/// one.pub fn gaussian_copula_parameter(target: f64, s1: f64, s2: f64) -> Option<f64> { let scale = lognormal_scale(s1, s2); if scale <= 0.0 { return (-1.0..=1.0).contains(&target).then_some(target); } let argument = 1.0 + target * scale; if argument <= 0.0 { return None; } let rho = argument.ln() / (s1 * s2); (-1.0..=1.0).contains(&rho).then_some(rho)} /// Kendall's tau of a Gaussian copula with correlation `rho`.////// `tau = (2/pi) arcsin(rho)`. A rank statistic, so unlike linear correlation it/// is unchanged by any increasing transformation of either margin — which is/// what makes it a property of the copula rather than of the marginals that/// happen to be attached to it.pub fn gaussian_kendall_tau(rho: f64) -> f64 { 2.0 / std::f64::consts::PI * rho.asin()} /// The coefficient of upper tail dependence.////// ```text/// lambda = lim as u -> 1 of P(U > u | V > u)/// ```////// The probability that one variable is extreme *given* the other is, in the/// limit of extreme. This is the number a senior tranche is a bet on, and the/// central fact of the dependence chapter is that the Gaussian copula sets it to zero for/// every correlation below one, while the Student-t sets it to////// ```text/// lambda = 2 * t_{nu+1}( -sqrt( (nu+1)(1-rho) / (1+rho) ) )./// ```////// Returns the Gaussian value, which is zero. Kept as a function rather than a/// constant because being able to write the two side by side is the point.pub fn gaussian_tail_dependence(_rho: f64) -> f64 { 0.0} /// The coefficient of upper tail dependence of a Student-t copula.////// Positive for every `rho > -1` and every finite `nu`, and it approaches the/// Gaussian zero only as `nu` grows. Two names with the same correlation and the/// same marginals can therefore have wildly different probabilities of failing/// together, which is precisely the freedom the credit chapter said was left open.pub fn t_tail_dependence(rho: f64, nu: f64) -> f64 { if rho >= 1.0 { return 1.0; } let argument = -((nu + 1.0) * (1.0 - rho) / (1.0 + rho)).sqrt(); 2.0 * t_cdf(argument, nu + 1.0)} // ---------------------------------------------------------------------------// Assembling pairwise correlations.// --------------------------------------------------------------------------- /// A Cholesky factorisation of a symmetric matrix, or `None` if it is not/// positive semi-definite.////// Plain Cholesky, no pivoting. For a matrix that really is PSD this succeeds;/// for one that is not, some diagonal entry under the square root goes/// negative, which is the algebraic signature of the failure and exactly what/// "no real random vector has this covariance" comes down to.pub fn cholesky(matrix: &[Vec<f64>]) -> Option<Vec<Vec<f64>>> { let n = matrix.len(); let mut l = vec![vec![0.0; n]; n]; for i in 0..n { for j in 0..=i { let mut sum = matrix[i][j]; for k in 0..j { sum -= l[i][k] * l[j][k]; } if i == j { if sum < -1e-9 { return None; } l[i][j] = sum.max(0.0).sqrt(); } else { if l[j][j] < 1e-12 { return None; } l[i][j] = sum / l[j][j]; } } } Some(l)} /// Whether three pairwise correlations can belong to any real joint/// distribution at all.////// Standardise the three variables and read each as a unit vector, so a/// correlation is the cosine of an angle between two of them (the inner/// product of the Brownian motion chapter's `L^2(Omega)`). Three vectors always/// have a positive semi-definite Gram matrix of their pairwise inner products,/// so the converse is the test: `(rho12, rho13, rho23)` are consistent exactly/// when the matrix with these off-diagonal entries and a unit diagonal is PSD,/// which for three variables is one inequality beyond `|rho_ij| <= 1`,////// ```text/// 1 - rho12^2 - rho13^2 - rho23^2 + 2 rho12 rho13 rho23 >= 0 ./// ```////// Assembling three pairwise correlations from three separate calibrations,/// as a quanto CMS spread option needs, has no reason to satisfy this.pub fn pairwise_correlations_are_consistent(rho12: f64, rho13: f64, rho23: f64) -> bool { let m = vec![ vec![1.0, rho12, rho13], vec![rho12, 1.0, rho23], vec![rho13, rho23, 1.0], ]; cholesky(&m).is_some()} // ---------------------------------------------------------------------------// A common intensity factor with jumps.// --------------------------------------------------------------------------- /// The common factor `Y` of `lambda_i = a_i Y + Z_i`, taken as an/// Ornstein-Uhlenbeck process driven by a compound Poisson process:////// ```text/// dY = -kappa Y dt + dJ, J = compound Poisson, rate `jump_rate`,/// jump sizes exponential with mean `jump_mean`./// ```////// Only the integral `Lambda_T = int_0^T Y dt` matters for survival and for/// tranches, and////// ```text/// Lambda_T = y0 B(T) + int_0^T B(T - s) dJ_s, B(t) = (1 - exp(-kappa t)) / kappa,/// ```////// a sum over the jumps, so it can be sampled exactly and its Laplace/// transform is in closed form.pub struct JumpFactor { pub kappa: f64, pub jump_rate: f64, pub jump_mean: f64, pub y0: f64,} impl JumpFactor { fn b(&self, t: f64) -> f64 { (1.0 - (-self.kappa * t).exp()) / self.kappa } /// `E[exp(-a Lambda_T)]`, the survival probability of a name that loads on /// the factor with weight `a` and has no other source of default. /// /// By the exponential formula for a Poisson process the jumps contribute /// `exp(rate int_0^T (E[exp(-a B(s) X)] - 1) ds)` with `X` exponential, and /// `E[exp(-u X)] = 1 / (1 + mean u)`. With `c = mean a / kappa` the integrand /// is `1 / ((1 + c) - c exp(-kappa s)) - 1`, which integrates in closed form. pub fn laplace(&self, a: f64, t: f64) -> f64 { let c = self.jump_mean * a / self.kappa; let integral = ((1.0 + c) * (self.kappa * t).exp() - c).ln() / (self.kappa * (1.0 + c)) - t; (-a * self.y0 * self.b(t) + self.jump_rate * integral).exp() } /// The mean and variance of `Lambda_T`, from the jump compensator. pub fn mean_and_variance(&self, t: f64) -> (f64, f64) { // int_0^T B(s) ds and int_0^T B(s)^2 ds, by the trapezoid rule; the // integrands are smooth. let n = 20_000; let h = t / n as f64; let (mut m1, mut m2) = (0.0, 0.0); for i in 0..=n { let w = if i == 0 || i == n { 0.5 } else { 1.0 }; let b = self.b(i as f64 * h); m1 += w * b * h; m2 += w * b * b * h; } let mean = self.y0 * self.b(t) + self.jump_rate * self.jump_mean * m1; let variance = self.jump_rate * 2.0 * self.jump_mean * self.jump_mean * m2; (mean, variance) } /// One exact draw of `Lambda_T`: a Poisson number of jumps at uniform times, /// each contributing `B(T - s)` times an exponential size. pub fn sample_integral(&self, rng: &mut crate::pathwise::Rng, t: f64) -> f64 { let mut total = self.y0 * self.b(t); let mut time = 0.0; loop { time += -rng.next_uniform().ln() / self.jump_rate; if time >= t { break; } let size = -rng.next_uniform().ln() * self.jump_mean; total += self.b(t - time) * size; } total }} /// The expected loss of a tranche `[attach, detach]`, as a fraction of its/// width, in the large-portfolio limit where every name loads `a` on the factor/// and defaults with probability `1 - exp(-a Lambda - z T)` given it.////// The loss fraction is then a deterministic function of `Lambda`, so the/// tranche is an expectation over the law of `Lambda` and nothing else.pub fn factor_tranche_loss( lambdas: &[f64], a: f64, z_t: f64, recovery: f64, attach: f64, detach: f64,) -> f64 { let total: f64 = lambdas .iter() .map(|&l| { let loss = (1.0 - recovery) * (1.0 - (-a * l - z_t).exp()); ((loss - attach).max(0.0) - (loss - detach).max(0.0)) / (detach - attach) }) .sum(); total / lambdas.len() as f64} // ---------------------------------------------------------------------------// Fitting the t.// --------------------------------------------------------------------------- /// The largest number of degrees of freedom the fits will report. A sample that/// looks Gaussian has a likelihood in `nu` that is still creeping upwards at/// any finite value, so the search stops here and the fit says so by returning/// this number.pub const NU_CAP: f64 = 200.0; /// A Student-t fitted to a sample: location, scale and degrees of freedom.pub struct StudentFit { pub location: f64, pub scale: f64, pub nu: f64, pub log_likelihood: f64,} /// The log density of a standard Student-t with `nu` degrees of freedom.pub(crate) fn ln_t_pdf(x: f64, nu: f64) -> f64 { ln_gamma((nu + 1.0) / 2.0) - ln_gamma(nu / 2.0) - 0.5 * (nu * std::f64::consts::PI).ln() - (nu + 1.0) / 2.0 * (1.0 + x * x / nu).ln()} /// Location, scale and log-likelihood of a t sample at fixed `nu`.////// The t is a normal with variance `nu / W`, `W ~ chi-square(nu)`, and given the/// latent `W` the location and scale are weighted normal estimates. EM/// alternates the two: the E step sets each observation's weight to the/// posterior mean of `W / nu`, `(nu + 1) / (nu + z^2)`, so an observation far/// out is judged to have come from a wide draw and counts for less, and the M/// step is the weighted mean and variance.fn t_fixed_nu(data: &[f64], nu: f64) -> (f64, f64, f64) { let n = data.len() as f64; let mut location = data.iter().sum::<f64>() / n; let mut scale2 = data.iter().map(|x| (x - location).powi(2)).sum::<f64>() / n; for _ in 0..500 { let (mut sw, mut swx) = (0.0, 0.0); let weights: Vec<f64> = data .iter() .map(|x| (nu + 1.0) / (nu + (x - location).powi(2) / scale2)) .collect(); for (x, w) in data.iter().zip(&weights) { sw += w; swx += w * x; } let new_location = swx / sw; let new_scale2 = data .iter() .zip(&weights) .map(|(x, w)| w * (x - new_location).powi(2)) .sum::<f64>() / n; let moved = (new_location - location).abs() + (new_scale2 - scale2).abs(); location = new_location; scale2 = new_scale2; if moved < 1e-12 * (1.0 + scale2) { break; } } let scale = scale2.sqrt(); let ll = data .iter() .map(|x| ln_t_pdf((x - location) / scale, nu) - scale.ln()) .sum::<f64>(); (location, scale, ll)} /// Maximise a function of `nu` on `[2, NU_CAP]`: a coarse search over a/// logarithmic grid, then golden section inside the bracket round the best/// point. The likelihood in `nu` is flat for large `nu`, which is why the grid/// is logarithmic.pub(crate) fn maximise_over_nu(mut f: impl FnMut(f64) -> f64) -> f64 { let (lo, hi) = (2.0f64.ln(), NU_CAP.ln()); let points = 24; let grid: Vec<f64> = (0..points) .map(|i| lo + (hi - lo) * i as f64 / (points - 1) as f64) .collect(); let values: Vec<f64> = grid.iter().map(|&l| f(l.exp())).collect(); let best = (0..points) .max_by(|&a, &b| values[a].partial_cmp(&values[b]).unwrap()) .unwrap(); let (mut a, mut b) = (grid[best.saturating_sub(1)], grid[(best + 1).min(points - 1)]); let phi = (5.0f64.sqrt() - 1.0) / 2.0; let (mut c, mut d) = (b - phi * (b - a), a + phi * (b - a)); let (mut fc, mut fd) = (f(c.exp()), f(d.exp())); for _ in 0..40 { if fc > fd { b = d; d = c; fd = fc; c = b - phi * (b - a); fc = f(c.exp()); } else { a = c; c = d; fc = fd; d = a + phi * (b - a); fd = f(d.exp()); } } ((a + b) / 2.0).exp()} /// Fit a Student-t to a sample by maximum likelihood.////// For each `nu` the location and scale come from EM ([`t_fixed_nu`]), which/// leaves a one-dimensional profile likelihood in `nu` to maximise. This/// measures the tails of a *marginal*. It is not the `nu` that governs tail/// dependence, which belongs to the copula: see [`fit_t_copula`].pub fn fit_student_t(data: &[f64]) -> StudentFit { let nu = maximise_over_nu(|nu| t_fixed_nu(data, nu).2); let (location, scale, log_likelihood) = t_fixed_nu(data, nu); StudentFit { location, scale, nu, log_likelihood }} /// A Student-t copula fitted to a pair of samples.pub struct TCopulaFit { pub rho: f64, pub nu: f64, pub log_likelihood: f64,} /// Midranks scaled into `(0, 1)`: `rank / (n + 1)`, ties sharing the average.////// Discards the marginals entirely, which is the point: the copula is what is/// left when only the order is kept.fn pseudo_observations(data: &[f64]) -> Vec<f64> { let n = data.len(); let mut order: Vec<usize> = (0..n).collect(); order.sort_by(|&a, &b| data[a].partial_cmp(&data[b]).unwrap()); let mut ranks = vec![0.0; n]; let mut i = 0; while i < n { let mut j = i; while j + 1 < n && data[order[j + 1]] == data[order[i]] { j += 1; } let mid = (i + j) as f64 / 2.0 + 1.0; for &k in &order[i..=j] { ranks[k] = mid / (n as f64 + 1.0); } i = j + 1; } ranks} /// Log-likelihood of the t copula at `(rho, nu)`, given the t quantiles `x, y`/// of the two pseudo-observation series at that `nu`.////// The copula density is the bivariate t density divided by its two marginal/// t densities, so the marginals cancel and only the dependence is scored.fn t_copula_log_likelihood(x: &[f64], y: &[f64], rho: f64, nu: f64) -> f64 { let constant = ln_gamma((nu + 2.0) / 2.0) - ln_gamma(nu / 2.0) - (nu * std::f64::consts::PI).ln() - 0.5 * (1.0 - rho * rho).ln(); x.iter() .zip(y) .map(|(&a, &b)| { let q = (a * a - 2.0 * rho * a * b + b * b) / (1.0 - rho * rho); constant - (nu + 2.0) / 2.0 * (1.0 + q / nu).ln() - ln_t_pdf(a, nu) - ln_t_pdf(b, nu) }) .sum()} /// The best correlation at a given `nu`, by golden section.fn best_rho(x: &[f64], y: &[f64], nu: f64) -> (f64, f64) { let phi = (5.0f64.sqrt() - 1.0) / 2.0; let (mut a, mut b) = (-0.995f64, 0.995f64); let (mut c, mut d) = (b - phi * (b - a), a + phi * (b - a)); let (mut fc, mut fd) = ( t_copula_log_likelihood(x, y, c, nu), t_copula_log_likelihood(x, y, d, nu), ); for _ in 0..45 { if fc > fd { b = d; d = c; fd = fc; c = b - phi * (b - a); fc = t_copula_log_likelihood(x, y, c, nu); } else { a = c; c = d; fc = fd; d = a + phi * (b - a); fd = t_copula_log_likelihood(x, y, d, nu); } } let rho = (a + b) / 2.0; (rho, t_copula_log_likelihood(x, y, rho, nu))} /// Fit a Student-t copula to two paired samples by maximum pseudo-likelihood.////// Each series is replaced by its ranks, so nothing about either marginal/// survives, and the copula's correlation and degrees of freedom are chosen to/// maximise the copula likelihood of the result. This is the `nu` of the tail/// dependence coefficient. It is not the `nu` from [`fit_student_t`] on each/// series, which describes how thick each marginal's tails are and says nothing/// about how the two move together in them.pub fn fit_t_copula(a: &[f64], b: &[f64]) -> TCopulaFit { assert_eq!(a.len(), b.len(), "paired samples must have the same length"); let (u, v) = (pseudo_observations(a), pseudo_observations(b)); let profile = |nu: f64| { let x: Vec<f64> = u.iter().map(|&p| t_inv(p, nu)).collect(); let y: Vec<f64> = v.iter().map(|&p| t_inv(p, nu)).collect(); best_rho(&x, &y, nu) }; let nu = maximise_over_nu(|nu| profile(nu).1); let (rho, log_likelihood) = profile(nu); TCopulaFit { rho, nu, log_likelihood }} // ---------------------------------------------------------------------------// Portfolio loss.// --------------------------------------------------------------------------- /// Which copula links the defaults.#[derive(Clone, Copy, Debug, PartialEq)]pub enum Copula { /// The one-factor Gaussian copula. Zero tail dependence by construction. Gaussian, /// A one-factor Student-t, whose common shock is scaled by an independent /// chi-square. Same correlation, same marginals, tails that do not vanish. StudentT { nu: f64 },} /// The distribution of the loss on a large homogeneous portfolio.////// Every name defaults by the horizon with probability `p`, loses `1 - recovery`/// when it does, and is linked to the others by `copula` with correlation `rho`.#[derive(Clone, Copy, Debug)]pub struct Portfolio { pub default_probability: f64, pub recovery: f64, pub correlation: f64, pub copula: Copula,} impl Portfolio { /// The default rate the portfolio settles at, given the common factor `m` /// and, for the Student-t, the mixing variable `w` drawn from a chi-square /// with `nu` degrees of freedom. /// /// In the large-portfolio limit the idiosyncratic risk averages away, so this /// *is* the realised default fraction rather than its expectation. The /// Gaussian case is Vasicek's formula. pub fn conditional_default_rate(&self, m: f64, w: f64) -> f64 { let rho = self.correlation.clamp(0.0, 0.999_999); let sqrt_rho = rho.sqrt(); let sqrt_one_minus = (1.0 - rho).sqrt(); match self.copula { Copula::Gaussian => { let threshold = norm_inv(self.default_probability); norm_cdf((threshold - sqrt_rho * m) / sqrt_one_minus) } Copula::StudentT { nu } => { // The t variate is sqrt(nu/w) * (sqrt(rho) m + sqrt(1-rho) z), // so conditioning on w turns the threshold into a scaled one and // the rest is Gaussian again. let threshold = t_inv(self.default_probability, nu); norm_cdf((threshold * (w / nu).sqrt() - sqrt_rho * m) / sqrt_one_minus) } } } /// The probability that the loss fraction exceeds `level`. /// /// For the Gaussian copula this is Vasicek's closed form, inverted. For the /// Student-t the mixing variable has to be integrated out, which is done on a /// grid over the chi-square density; the integrand is smooth and one /// dimensional, so a few hundred points is far more accuracy than the model /// deserves. pub fn exceedance(&self, level: f64) -> f64 { let loss_given_default = 1.0 - self.recovery; if loss_given_default <= 0.0 { return 0.0; } // Convert a loss level into the default rate that produces it. let rate = level / loss_given_default; if rate <= 0.0 { return 1.0; } if rate >= 1.0 { return 0.0; } let rho = self.correlation.clamp(1e-9, 0.999_999); let sqrt_rho = rho.sqrt(); let sqrt_one_minus = (1.0 - rho).sqrt(); match self.copula { Copula::Gaussian => { // Loss exceeds the level exactly when the factor is low enough. let threshold = norm_inv(self.default_probability); let critical = (threshold - sqrt_one_minus * norm_inv(rate)) / sqrt_rho; norm_cdf(critical) } Copula::StudentT { nu } => { let threshold = t_inv(self.default_probability, nu); chi_square_average(nu, |w| { // Given w, the same inversion as the Gaussian case. let critical = (threshold * (w / nu).sqrt() - sqrt_one_minus * norm_inv(rate)) / sqrt_rho; norm_cdf(critical) }) } } } /// The expected loss on the tranche covering `[attach, detach]`, as a /// fraction of the tranche's own notional. /// /// A tranche is a call spread on the portfolio loss, so its expected loss is /// /// ```text /// ( E[(L - a)+] - E[(L - d)+] ) / (d - a) /// ``` /// /// and each call is the integral of the exceedance probability above its /// strike. This is Breeden-Litzenberger from the local volatility chapter /// read backwards, and the dependence chapter makes something of that: the /// capital structure of a portfolio is a strip of call spreads on one /// variable, so a complete set of tranche quotes implies a loss /// distribution the same way a complete set of option quotes implies a /// density. pub fn tranche_expected_loss(&self, attach: f64, detach: f64) -> f64 { if detach <= attach { return 0.0; } (self.call_on_loss(attach) - self.call_on_loss(detach)) / (detach - attach) } /// The common factor above which a tranche attaching at `strike` takes no /// loss at all, given the mixing variable `w`. /// /// The payoff `(loss(m) - strike)+` has a kink exactly here, and it is worth /// knowing where: Simpson's rule integrated straight through a kink converges /// at first order instead of fourth, which showed up as the equity tranche /// disagreeing with itself in the third significant figure. Splitting the /// integration at this point removes the problem rather than out-resolving /// it. /// /// `None` when the tranche is never touched. fn critical_factor(&self, strike: f64, w: f64) -> Option<f64> { let loss_given_default = 1.0 - self.recovery; let rate = strike / loss_given_default; if rate >= 1.0 { return None; } if rate <= 0.0 { // A tranche attaching at zero always takes some loss, so there is no // factor above which the payoff switches off. return Some(f64::INFINITY); } let rho = self.correlation.clamp(1e-9, 0.999_999); let threshold = match self.copula { Copula::Gaussian => norm_inv(self.default_probability), Copula::StudentT { nu } => t_inv(self.default_probability, nu) * (w / nu).sqrt(), }; // The payoff is positive for m below this. Some((threshold - (1.0 - rho).sqrt() * norm_inv(rate)) / rho.sqrt()) } /// `E[(L - strike)+]`. /// /// Integrated over the common factors rather than over the loss level. Both /// routes are correct and the factor one is far better behaved: as a function /// of the loss level the integrand has a steep shoulder near zero that a /// fixed grid resolves badly, while as a function of the factor it is smooth /// and monotone apart from the single kink located by [`Portfolio::critical_factor`]. fn call_on_loss(&self, strike: f64) -> f64 { let loss_given_default = 1.0 - self.recovery; if strike >= loss_given_default { return 0.0; } let integrate = |w: f64| match self.critical_factor(strike, w) { None => 0.0, Some(upper) => normal_average_below(upper, |m| { self.conditional_default_rate(m, w) * loss_given_default - strike }), }; match self.copula { Copula::Gaussian => integrate(1.0), Copula::StudentT { nu } => chi_square_average(nu, integrate), } }} /// `E[ f(M) 1{M < upper} ]` for a standard normal `M`, by Simpson.////// The caller passes the point where its integrand stops contributing, so the/// grid ends exactly on the kink rather than straddling it. Eight standard/// deviations is the lower limit, past which the weight is smaller than anything/// it multiplies.fn normal_average_below(upper: f64, f: impl Fn(f64) -> f64) -> f64 { const STEPS: usize = 100; const LIMIT: f64 = 8.0; let upper = upper.min(LIMIT); if upper <= -LIMIT { return 0.0; } let h = (upper + LIMIT) / STEPS as f64; let normalisation = (2.0 * std::f64::consts::PI).sqrt(); let mut total = 0.0; for i in 0..=STEPS { let m = -LIMIT + i as f64 * h; let weight = if i == 0 || i == STEPS { 1.0 } else if i % 2 == 1 { 4.0 } else { 2.0 }; total += weight * (-0.5 * m * m).exp() * f(m); } total * h / 3.0 / normalisation} /// The average of `f(w)` over a chi-square with `nu` degrees of freedom.////// Simpson on the log of the variable, which keeps the grid tight where the/// density is and still reaches far enough into the right tail to matter: it is/// a small `w` that produces the fat tail of a Student-t, so the left end is the/// end that has to be resolved.fn chi_square_average(nu: f64, f: impl Fn(f64) -> f64) -> f64 { const STEPS: usize = 150; // Six e-folds either side of the mode covers the density to far below the // precision of anything it is multiplied by. let (lo, hi) = ((nu * 1e-4).ln(), (nu * 40.0).ln()); let h = (hi - lo) / STEPS as f64; let mut total = 0.0; let mut mass = 0.0; for i in 0..=STEPS { let log_w = lo + i as f64 * h; let w = log_w.exp(); let weight = if i == 0 || i == STEPS { 1.0 } else if i % 2 == 1 { 4.0 } else { 2.0 }; // The chi-square density times w, since the variable of integration is // log w. Normalising constants cancel in the ratio below. let density = w.powf(0.5 * nu) * (-0.5 * w).exp(); total += weight * density * f(w); mass += weight * density; } total / mass} #[cfg(test)]mod tests { use super::*; use crate::pathwise::Rng; #[test] // Worth getting right rather than assuming. Two lognormals built from // the *same* normal are the same variable up to a monotone map, so the // upper bound is genuinely one. It is the lower bound that collapses: // exp(X) and exp(-X) are the closest to opposed that two lognormals can // be, and that is not very close at all. let (lo, hi) = lognormal_correlation_bounds(1.0, 1.0); assert!((hi - 1.0).abs() < 1e-12, "upper bound was {hi}"); assert!((lo + 0.3679).abs() < 1e-3, "lower bound was {lo}"); } fn example_factor() -> JumpFactor { JumpFactor { kappa: 0.5, jump_rate: 0.4, jump_mean: 0.05, y0: 0.01 } } #[test] fn the_factor_transform_matches_exact_simulation() { let (factor, t) = (example_factor(), 5.0f64); let mut rng = Rng::new(2024); let n = 400_000; let draws: Vec<f64> = (0..n).map(|_| factor.sample_integral(&mut rng, t)).collect(); for a in [1.0f64, 5.0, 20.0] { let sampled = draws.iter().map(|l| (-a * l).exp()).sum::<f64>() / n as f64; let closed = factor.laplace(a, t); assert!((sampled - closed).abs() < 3e-3, "a={a}: sampled {sampled}, closed form {closed}"); } let (mean, variance) = factor.mean_and_variance(t); let sample_mean = draws.iter().sum::<f64>() / n as f64; let sample_var = draws.iter().map(|l| (l - sample_mean).powi(2)).sum::<f64>() / n as f64; assert!((sample_mean - mean).abs() < 2e-3, "mean {sample_mean} against {mean}"); assert!((sample_var / variance - 1.0).abs() < 0.05, "variance {sample_var} against {variance}"); } #[test] fn matching_mean_and_variance_a_jump_factor_loads_the_senior_tranche_more_than_a_diffusive_one() { // The diffusive comparison has the same mean and variance of the integral // Lambda and a Gaussian law, which is what the integral of a diffusive // Ornstein-Uhlenbeck intensity has. The jump factor's law has a heavy // right tail, and the senior tranche is a bet on it. let (factor, t) = (example_factor(), 5.0f64); let (mean, variance) = factor.mean_and_variance(t); let mut rng = Rng::new(77); let n = 400_000; let jump: Vec<f64> = (0..n).map(|_| factor.sample_integral(&mut rng, t)).collect(); let gauss: Vec<f64> = (0..n).map(|_| (mean + variance.sqrt() * rng.next_normal()).max(0.0)).collect(); // Same mean and variance of Lambda, so the same first two moments of // everything upstream, and the middle of the capital structure loses a // little less under jumps: the jump factor is quiet most of the time and // makes up for it in the far tail, where the senior tranche loses // several times what it does under the diffusive factor. let ratio = |attach: f64, detach: f64| { let j = factor_tranche_loss(&jump, 1.0, 0.02, 0.4, attach, detach); let g = factor_tranche_loss(&gauss, 1.0, 0.02, 0.4, attach, detach); j / g }; assert!(ratio(0.03, 0.10) < 1.0, "mezzanine ratio {}", ratio(0.03, 0.10)); assert!(ratio(0.10, 0.20) < 1.0, "upper mezzanine ratio {}", ratio(0.10, 0.20)); assert!(ratio(0.20, 0.60) > 2.5, "senior ratio {}", ratio(0.20, 0.60)); } /// A standard t draw: a normal over the root of an independent chi-square /// divided by its degrees of freedom. fn draw_t(rng: &mut Rng, nu: usize) -> (f64, f64) { let w: f64 = (0..nu).map(|_| rng.next_normal().powi(2)).sum(); (rng.next_normal(), (nu as f64 / w).sqrt()) } #[test] fn the_marginal_fit_recovers_the_degrees_of_freedom_of_a_simulated_t() { let mut rng = Rng::new(555); let sample: Vec<f64> = (0..20_000) .map(|_| { let (z, s) = draw_t(&mut rng, 4); 2.0 + 3.0 * z * s }) .collect(); let fit = fit_student_t(&sample); assert!((fit.nu - 4.0).abs() < 0.4, "nu {}", fit.nu); assert!((fit.location - 2.0).abs() < 0.1, "location {}", fit.location); assert!((fit.scale - 3.0).abs() < 0.1, "scale {}", fit.scale); // A Gaussian sample has no finite best nu: the fit runs to large values. let normal: Vec<f64> = (0..20_000).map(|_| rng.next_normal()).collect(); assert!(fit_student_t(&normal).nu > 40.0, "a Gaussian sample fitted nu = {}", fit_student_t(&normal).nu); } #[test] fn the_copula_fit_recovers_nu_and_ignores_the_marginals() { // A bivariate t with nu = 4 and rho = 0.5, then a monotone map applied to // each coordinate. The ranks are unchanged, so the fit must not move, and // the marginal fit of the mapped data must be nothing like nu = 4. let (nu, rho) = (4usize, 0.5f64); let mut rng = Rng::new(808); let (mut a, mut b) = (Vec::new(), Vec::new()); for _ in 0..4_000 { let (z1, s) = draw_t(&mut rng, nu); let z2 = rho * z1 + (1.0 - rho * rho).sqrt() * rng.next_normal(); a.push((0.7 * z1 * s).exp()); b.push((z2 * s).powi(3)); } let fit = fit_t_copula(&a, &b); assert!((fit.nu - 4.0).abs() < 1.0, "copula nu {}", fit.nu); assert!((fit.rho - 0.5).abs() < 0.05, "copula rho {}", fit.rho); let marginal = fit_student_t(&a).nu; assert!((marginal - 4.0).abs() > 1.0, "the marginal fit of a lognormal-ish series gave {marginal}"); } /// Two columns of a committed Treasury panel, as daily changes in basis /// points. fn treasury_changes(labels: [&str; 2]) -> (Vec<f64>, Vec<f64>) { let path = concat!(env!("CARGO_MANIFEST_DIR"), "/../public/marketdata/treasury-history-long.json"); let raw = std::fs::read_to_string(path).expect("run `npm run marketdata`"); let obs = &raw[raw.find("\"observations\"").unwrap()..]; let mut levels = [Vec::new(), Vec::new()]; for (i, _) in obs.match_indices("\"rates\"") { let rest = &obs[i..]; let seg = &rest[..rest.find('}').unwrap()]; let read = |label: &str| -> Option<f64> { let key = format!("\"{label}\":"); let at = seg.find(&key)? + key.len(); let tail = &seg[at..]; let stop = tail.find(|c| c == ',' || c == '\n').unwrap_or(tail.len()); tail[..stop].trim().parse().ok() }; if let (Some(x), Some(y)) = (read(labels[0]), read(labels[1])) { levels[0].push(x); levels[1].push(y); } } let diff = |v: &Vec<f64>| v.windows(2).map(|w| (w[1] - w[0]) * 1e4).collect::<Vec<f64>>(); (diff(&levels[0]), diff(&levels[1])) } #[test] fn fitting_a_t_to_treasury_yield_changes() { let (two, ten) = treasury_changes(["2 Yr", "10 Yr"]); let (m2, m10) = (fit_student_t(&two), fit_student_t(&ten)); let copula = fit_t_copula(&two, &ten); // The panel is a committed snapshot, so these are fixed numbers to the // precision quoted in the chapter. The point is the last comparison: // the copula's nu is not either marginal's. assert_eq!(two.len(), 11_241); assert!((m2.nu - 2.5).abs() < 0.15, "2y marginal nu {}", m2.nu); assert!((m10.nu - 4.1).abs() < 0.15, "10y marginal nu {}", m10.nu); assert!((copula.nu - 5.1).abs() < 0.3, "copula nu {}", copula.nu); assert!((copula.rho - 0.82).abs() < 0.02, "copula rho {}", copula.rho); assert!((copula.nu - m2.nu).abs() > 1.5, "copula nu {} against 2y marginal {}", copula.nu, m2.nu); } /// `P(Y < x | X = s)` for the bivariate t with correlation `rho` and `nu` /// degrees of freedom, from the conditional law derived in the chapter: a /// t with `nu + 1` degrees of freedom, mean `rho s` and scale /// `sqrt((1 - rho^2)(nu + s^2)/(nu + 1))`. fn conditional_below(x: f64, s: f64, rho: f64, nu: f64) -> f64 { let scale = ((1.0 - rho * rho) * (nu + s * s) / (nu + 1.0)).sqrt(); t_cdf((x - rho * s) / scale, nu + 1.0) } #[test] fn the_conditional_law_of_a_bivariate_t_is_a_t_with_one_more_degree_of_freedom() { // Simulate the construction (X, Y) = sqrt(nu / W) (Z1, Z2) with // W ~ chi-square(nu), keep the draws whose X lies in a thin window round // x0, and compare what Y does there with the conditional law. let (nu, rho, x0, half) = (4usize, 0.3f64, -2.0f64, 0.05f64); let mut rng = Rng::new(31337); let (mut kept, mut below, mut sum) = (0usize, 0usize, 0.0f64); for _ in 0..3_000_000 { let w: f64 = (0..nu).map(|_| rng.next_normal().powi(2)).sum(); let (z1, e) = (rng.next_normal(), rng.next_normal()); let z2 = rho * z1 + (1.0 - rho * rho).sqrt() * e; let scale = (nu as f64 / w).sqrt(); let (x, y) = (scale * z1, scale * z2); if (x - x0).abs() < half { kept += 1; sum += y; if y < x0 { below += 1; } } } assert!(kept > 15_000, "only {kept} draws in the window"); let empirical = below as f64 / kept as f64; let formula = conditional_below(x0, x0, rho, nu as f64); assert!((empirical - formula).abs() < 0.01, "sample {empirical}, formula {formula}"); let mean = sum / kept as f64; assert!((mean - rho * x0).abs() < 0.05, "conditional mean {mean}, expected {}", rho * x0); } #[test] fn the_tail_dependence_coefficient_is_twice_the_limit_of_the_conditional_probability() { // lambda = lim P(Y < x | X < x), and by l'Hopital that is twice the // limit of P(Y < x | X = x), because the joint tail P(X < x, Y < x) moves // with both coordinates. Compute the joint tail by quadrature from the // conditional law, at an extreme x, and compare. let (nu, rho, x) = (4.0f64, 0.3f64, -1000.0f64); let density = |s: f64| { (ln_gamma((nu + 1.0) / 2.0) - ln_gamma(nu / 2.0)).exp() / (nu * std::f64::consts::PI).sqrt() * (1.0 + s * s / nu).powf(-(nu + 1.0) / 2.0) }; // s = x (1 + v) runs from x down the tail as v runs from zero. let (top, n) = (60.0f64, 60_000usize); let h = top / n as f64; let joint: f64 = (0..=n) .map(|i| { let v = i as f64 * h; let s = x * (1.0 + v); let weight = if i == 0 || i == n { 0.5 } else { 1.0 }; weight * density(s) * conditional_below(x, s, rho, nu) * x.abs() * h }) .sum(); let given_below = joint / t_cdf(x, nu); let given_equal = conditional_below(x, x, rho, nu); let lambda = t_tail_dependence(rho, nu); assert!((given_below - lambda).abs() < 2e-3, "P(Y<x|X<x) = {given_below}, lambda = {lambda}"); assert!( (given_below / given_equal - 2.0).abs() < 0.02, "ratio {}", given_below / given_equal ); } #[test] fn three_pairwise_correlations_calibrated_separately_need_not_be_consistent() { // rho12: the two rates, as calibrated to the CMS spread option. rho13, // rho23: each rate's correlation with the exchange rate, as calibrated // to a quanto adjustment on that rate alone. Nothing ties the three // together, and this triple fails the determinant test. let (r12, r13, r23) = (0.8f64, 0.7f64, -0.1f64); let det = 1.0 - r12 * r12 - r13 * r13 - r23 * r23 + 2.0 * r12 * r13 * r23; assert!(det < 0.0, "expected an inconsistent triple, det = {det}"); assert!(!pairwise_correlations_are_consistent(r12, r13, r23)); // A triple that does satisfy the inequality passes, and Cholesky // succeeds on it. let (r12, r13, r23) = (0.8f64, 0.7f64, 0.5f64); let det = 1.0 - r12 * r12 - r13 * r13 - r23 * r23 + 2.0 * r12 * r13 * r23; assert!(det > 0.0, "expected a consistent triple, det = {det}"); assert!(pairwise_correlations_are_consistent(r12, r13, r23)); } #[test] fn opposed_lognormals_at_equal_volatility_correlate_at_minus_exp_of_minus_variance() { // exp(sZ) and exp(-sZ) multiply to one, so E[XY] = 1 and the covariance // is 1 - exp(s^2) against a variance of exp(s^2)(exp(s^2) - 1). The // ratio is -exp(-s^2), and the table in the chapter is this formula. for sigma in [0.25f64, 1.0, 2.0] { let (lo, _) = lognormal_correlation_bounds(sigma, sigma); let closed = -(-sigma * sigma).exp(); assert!((lo - closed).abs() < 1e-12, "sigma={sigma}: {lo} against {closed}"); } // The mean and standard deviation of exp(2Z), quoted in the chapter. let mean = 2.0f64.exp(); let sd = (2.0f64.exp()) * ((4.0f64).exp() - 1.0).sqrt(); assert!((mean - 7.39).abs() < 0.01, "mean {mean}"); assert!((sd - 54.1).abs() < 0.1, "standard deviation {sd}"); } #[test] fn the_lower_bound_collapses_towards_zero_as_volatility_rises() { // At a 200% volatility two lognormals cannot be more than two percent // negatively correlated, and at 300% they cannot be negatively // correlated in any meaningful sense at all. A risk system that accepts // -0.5 here has accepted a number no joint distribution can produce. let mut previous = -1.0; for sigma in [0.25, 0.5, 1.0, 2.0, 3.0] { let (lo, _) = lognormal_correlation_bounds(sigma, sigma); assert!(lo > previous, "lower bound not rising at sigma={sigma}"); assert!(lo < 0.0); previous = lo; } let (lo, _) = lognormal_correlation_bounds(3.0, 3.0); assert!(lo > -0.001, "lower bound at sigma=3 was {lo}"); } #[test] fn the_gaussian_copula_parameter_is_not_the_pearson_correlation() { // Simulate the coupling directly, Z2 = rho Z1 + sqrt(1 - rho^2) W, and // measure the Pearson correlation of the exponentials. The formula was // derived without a sample, so agreement is not by construction. let (s1, s2) = (0.5f64, 0.5f64); for rho in [-0.5f64, 0.6] { let mut rng = Rng::new(9001); let n = 400_000; let (mut sx, mut sy, mut sxx, mut syy, mut sxy) = (0.0, 0.0, 0.0, 0.0, 0.0); for _ in 0..n { let (z1, w) = (rng.next_normal(), rng.next_normal()); let z2 = rho * z1 + (1.0 - rho * rho).sqrt() * w; let (x, y) = ((s1 * z1).exp(), (s2 * z2).exp()); sx += x; sy += y; sxx += x * x; syy += y * y; sxy += x * y; } let n = n as f64; let (mx, my) = (sx / n, sy / n); let corr = (sxy / n - mx * my) / ((sxx / n - mx * mx).sqrt() * (syy / n - my * my).sqrt()); let formula = gaussian_copula_pearson(rho, s1, s2); assert!((corr - formula).abs() < 0.01, "rho={rho}: sample {corr}, formula {formula}"); // Already visible at 50% volatility; the 200% case below is the one // where it is large. assert!((formula - rho).abs() > 0.02, "rho={rho}: the two coincided at {formula}"); } } #[test] fn the_copula_parameter_sweeps_exactly_the_attainable_range() { // rho = -1 and 1 are the two ends of the bounds theorem, and the map is // increasing between them. for (s1, s2) in [(0.25, 0.25), (1.0, 1.0), (0.5, 2.0)] { let (lo, hi) = lognormal_correlation_bounds(s1, s2); assert!((gaussian_copula_pearson(-1.0, s1, s2) - lo).abs() < 1e-12); assert!((gaussian_copula_pearson(1.0, s1, s2) - hi).abs() < 1e-12); let mut previous = lo - 1.0; for i in 0..=20 { let rho = -1.0 + 0.1 * i as f64; let value = gaussian_copula_pearson(rho, s1, s2); assert!(value > previous, "not increasing at rho={rho}"); previous = value; } } } #[test] fn calibrating_to_an_unattainable_pearson_correlation_fails() { // A target inside the range inverts and round-trips; a target outside it // has no copula parameter, which is the bounds theorem from the other // side. At 25% volatility -0.5 is attainable, at 200% it is not. let rho = gaussian_copula_parameter(-0.5, 0.25, 0.25).expect("attainable"); assert!((gaussian_copula_pearson(rho, 0.25, 0.25) + 0.5).abs() < 1e-12); assert!(gaussian_copula_parameter(-0.5, 2.0, 2.0).is_none()); assert!(gaussian_copula_parameter(0.9, 0.5, 2.0).is_none()); } #[test] fn a_gaussian_copula_parameter_of_minus_a_half_at_200_percent_volatility_is_almost_zero() { // The number quoted in the chapter's example: a perfectly valid copula // whose Pearson correlation is nowhere near its parameter. let pearson = gaussian_copula_pearson(-0.5, 2.0, 2.0); assert!((pearson + 0.0161).abs() < 5e-4, "Pearson correlation was {pearson}"); } #[test] fn unequal_volatilities_cannot_reach_one_either() { // The upper bound only survives because the two volatilities were equal. // Make them differ and the ceiling drops fast: a 50% name and a 200% // name cannot be more than 44% correlated however they are coupled. let (lo, hi) = lognormal_correlation_bounds(0.5, 2.0); assert!((hi - 0.4404).abs() < 1e-3, "upper bound {hi}"); assert!((lo + 0.1620).abs() < 1e-3, "lower bound {lo}"); // And the ceiling is monotone in how far apart the two volatilities are. let mut previous = 1.0; for spread in [1.0, 2.0, 4.0, 8.0] { let (_, hi) = lognormal_correlation_bounds(0.5, 0.5 * spread); assert!(hi <= previous + 1e-12, "ceiling rose at spread={spread}"); previous = hi; } } #[test] fn the_bounds_are_reached_by_the_comonotone_pair() { // The upper bound is the correlation of exp(X) with exp(X) rescaled, so // simulating the comonotone pair has to reproduce it. let (s1, s2) = (0.8f64, 1.2f64); let (_, analytic) = lognormal_correlation_bounds(s1, s2); let mut rng = Rng::new(4242); let n = 150_000; let (mut sx, mut sy, mut sxx, mut syy, mut sxy) = (0.0, 0.0, 0.0, 0.0, 0.0); for _ in 0..n { let z = rng.next_normal(); let (x, y) = ((s1 * z).exp(), (s2 * z).exp()); sx += x; sy += y; sxx += x * x; syy += y * y; sxy += x * y; } let n = n as f64; let (mx, my) = (sx / n, sy / n); let cov = sxy / n - mx * my; let simulated = cov / ((sxx / n - mx * mx).sqrt() * (syy / n - my * my).sqrt()); assert!( (simulated - analytic).abs() < 0.02, "simulated {simulated} against analytic {analytic}" ); } #[test] fn the_gaussian_copula_has_no_tail_dependence_and_the_t_does() { // The single most consequential fact in the chapter. for rho in [0.1, 0.3, 0.5, 0.9, 0.99] { assert_eq!(gaussian_tail_dependence(rho), 0.0); let t = t_tail_dependence(rho, 4.0); assert!(t > 0.0, "t copula had no tail dependence at rho={rho}"); } // At a correlation of 0.3, two names under a t(4) copula fail together // in the limit about a sixth of the time, against never. let lambda = t_tail_dependence(0.3, 4.0); assert!(lambda > 0.10 && lambda < 0.25, "lambda was {lambda}"); } #[test] fn tail_dependence_vanishes_as_the_t_becomes_normal() { // The continuity check: a t copula with many degrees of freedom is a // Gaussian copula, so its tail dependence must go to the Gaussian zero. let mut previous = 1.0; for nu in [2.0, 4.0, 10.0, 40.0, 200.0] { let lambda = t_tail_dependence(0.5, nu); assert!(lambda < previous, "not decreasing at nu={nu}"); previous = lambda; } assert!(t_tail_dependence(0.5, 5000.0) < 1e-3); } #[test] fn the_gaussian_loss_distribution_is_vasicek() { // Checked against a simulation of the factor, which shares only the // conditional default rate with the closed form. let p = Portfolio { default_probability: 0.05, recovery: 0.4, correlation: 0.3, copula: Copula::Gaussian, }; let mut rng = Rng::new(918_273); let n = 150_000; for level in [0.005, 0.02, 0.05, 0.10] { let mut count = 0.0; for _ in 0..n { let m = rng.next_normal(); if p.conditional_default_rate(m, 1.0) * 0.6 > level { count += 1.0; } } let simulated = count / n as f64; let analytic = p.exceedance(level); assert!( (simulated - analytic).abs() < 0.006, "at {level}: simulated {simulated} against analytic {analytic}" ); } } #[test] fn the_t_loss_distribution_matches_a_simulation() { // Same check for the branch that integrates out the mixing variable, // which is the one with room to be wrong. The chi-square is simulated as // a sum of squared normals, so the quadrature and the simulation share // no code at all. let nu = 6.0; let p = Portfolio { default_probability: 0.05, recovery: 0.4, correlation: 0.3, copula: Copula::StudentT { nu }, }; let mut rng = Rng::new(555_111); let n = 150_000; for level in [0.005, 0.02, 0.05, 0.10] { let mut count = 0.0; for _ in 0..n { let m = rng.next_normal(); let w: f64 = (0..nu as usize).map(|_| rng.next_normal().powi(2)).sum(); if p.conditional_default_rate(m, w) * 0.6 > level { count += 1.0; } } let simulated = count / n as f64; let analytic = p.exceedance(level); assert!( (simulated - analytic).abs() < 0.007, "at {level}: simulated {simulated} against analytic {analytic}" ); } } #[test] fn the_whole_capital_structure_sums_to_the_portfolio_loss() { // The identity that makes tranches a strip of call spreads: the // notional-weighted expected losses of a partition of [0, 1-R] have to // add back up to the expected loss of the portfolio itself. It is also // the strongest available check on the integration, since it has to hold // for either copula and any correlation. for copula in [Copula::Gaussian, Copula::StudentT { nu: 5.0 }] { for correlation in [0.05, 0.3, 0.7] { let p = Portfolio { default_probability: 0.05, recovery: 0.4, correlation, copula, }; let edges = [0.0, 0.03, 0.07, 0.15, 0.30, 0.60]; let total: f64 = edges .windows(2) .map(|w| p.tranche_expected_loss(w[0], w[1]) * (w[1] - w[0])) .sum(); // The portfolio's own expected loss, which needs no model. let expected = p.default_probability * (1.0 - p.recovery); assert!( (total - expected).abs() < 1e-6, "{copula:?} at rho={correlation}: tranches gave {total}, portfolio {expected}" ); } } } #[test] fn the_senior_tranche_is_where_the_copula_choice_shows_up() { // The dependence chapter's headline number. Same marginals, same correlation, same // everything a risk report records -- and a senior tranche that is worth // a different order of magnitude. let gaussian = Portfolio { default_probability: 0.05, recovery: 0.4, correlation: 0.3, copula: Copula::Gaussian, }; let student = Portfolio { copula: Copula::StudentT { nu: 4.0 }, ..gaussian }; // The effect is not a uniform shift: it is a transfer up the capital // structure. Equity gets *safer* under the fat-tailed copula, because // more of the probability sits at low losses; everything senior gets // worse, and the further up the worse it gets. let ratio = |a: f64, d: f64| { student.tranche_expected_loss(a, d) / gaussian.tranche_expected_loss(a, d) }; let equity = ratio(0.0, 0.03); assert!(equity < 0.85, "equity ratio was {equity}, expected a fall"); let mut previous = equity; for &(a, d) in &[(0.03, 0.07), (0.07, 0.15), (0.15, 0.30), (0.30, 0.60)] { let r = ratio(a, d); assert!(r > previous, "ratio not rising at {a}-{d}: {r} after {previous}"); previous = r; } // And the super senior, the tranche the whole argument is about, is out // by an order of magnitude rather than a few percent. assert!(previous > 5.0, "super senior ratio was only {previous}"); }} /// A spread option on two rates, priced from their marginals and a copula.////// This is the rates counterpart of the tranche above and it is set up the same/// way: both rates keep exactly the same marginal distribution whatever the/// copula, and the correlation is held fixed too, so anything that moves between/// two runs is the dependence structure and nothing else.////// The marginals are normal, which is the rates convention and also convenient/// here — with a Gaussian copula the spread is then normal and the option has a/// closed form, so one of the two columns can be checked against Bachelier.#[derive(Clone, Copy, Debug)]pub struct SpreadOption { pub forward1: f64, pub forward2: f64, /// Absolute volatilities, in rate units per root year. pub vol1: f64, pub vol2: f64, pub expiry: f64, pub correlation: f64, pub copula: Copula,} impl SpreadOption { /// `E[(S1 - S2 - K)^+]`, undiscounted, by simulation. pub fn price(&self, strike: f64, paths: usize, seed: u64) -> f64 { let mut rng = crate::pathwise::Rng::new(seed); let rho = self.correlation.clamp(-0.999_999, 0.999_999); let (s1, s2) = (self.vol1 * self.expiry.sqrt(), self.vol2 * self.expiry.sqrt()); let mut total = 0.0; for _ in 0..paths { let z1 = rng.next_normal(); let z2 = rho * z1 + (1.0 - rho * rho).sqrt() * rng.next_normal(); // The copula enters only here, as the map from the correlated pair // to uniforms. The marginals are applied afterwards and are the same // in both branches. let (u1, u2) = match self.copula { Copula::Gaussian => (norm_cdf(z1), norm_cdf(z2)), Copula::StudentT { nu } => { // Chi-square with nu degrees of freedom, nu a whole number. let k = nu.round().max(1.0) as usize; let w: f64 = (0..k).map(|_| rng.next_normal().powi(2)).sum::<f64>() / k as f64; let scale = w.max(1e-12).sqrt(); (t_cdf(z1 / scale, nu), t_cdf(z2 / scale, nu)) } }; let r1 = self.forward1 + s1 * norm_inv(u1); let r2 = self.forward2 + s2 * norm_inv(u2); total += (r1 - r2 - strike).max(0.0); } total / paths as f64 } /// The Bachelier price the Gaussian case must reproduce, since normal /// marginals joined by a Gaussian copula make the spread normal. pub fn gaussian_closed_form(&self, strike: f64) -> f64 { let (s1, s2) = (self.vol1 * self.expiry.sqrt(), self.vol2 * self.expiry.sqrt()); let sd = (s1 * s1 + s2 * s2 - 2.0 * self.correlation * s1 * s2).sqrt(); let forward = self.forward1 - self.forward2; if sd <= 0.0 { return (forward - strike).max(0.0); } let d = (forward - strike) / sd; (forward - strike) * norm_cdf(d) + sd * crate::black::norm_pdf(d) }} #[cfg(test)]mod spread_option_tests { use super::*; fn base(copula: Copula) -> SpreadOption { SpreadOption { forward1: 0.040, forward2: 0.025, vol1: 0.010, vol2: 0.010, expiry: 5.0, correlation: 0.8, copula, } } /// Normal marginals joined by a Gaussian copula give a normal spread, so the /// simulation has a closed form to answer to. #[test] fn the_gaussian_case_reproduces_bachelier() { let m = base(Copula::Gaussian); for &k in &[0.0, 0.015, 0.03] { let simulated = m.price(k, 400_000, 11); let exact = m.gaussian_closed_form(k); assert!( (simulated - exact).abs() < 0.02 * exact, "strike {k}: simulated {simulated:.6} against {exact:.6}" ); } } /// The dependence chapter's claim for rates. Both rates keep the same /// marginal and the same correlation; only the copula changes, and the far /// strike spread option moves by a multiple of its value. /// /// The direction is the same as the credit case, though the mechanism reads /// differently at first. A Student-t copula is a bivariate normal divided by /// a common random scale, and a small draw of that scale makes both rates /// extreme without making them equal — so the *spread* inherits the heavy /// tail too. The Gaussian copula therefore understates a far strike spread /// option for the same reason it understated a senior tranche. /// /// Near the money it goes the other way and by much less: the extra mass in /// the tails has come from somewhere, and it is taken out of the middle. #[test] fn the_copula_moves_the_far_strike_spread_option() { let gaussian = base(Copula::Gaussian); let student = base(Copula::StudentT { nu: 4.0 }); let ratio = |k: f64| student.price(k, 600_000, 7) / gaussian.price(k, 600_000, 7); // Slightly cheaper in the middle, where the mass was taken from. let middle = ratio(0.015); assert!( (middle - 0.949).abs() < 0.02, "at the money the ratio was {middle:.3}, expected about 0.95" ); // And a large multiple far out, which is where the instrument is a bet // on the two rates coming apart. let far = ratio(0.045); assert!((far - 2.60).abs() < 0.25, "far strike ratio was {far:.3}, expected about 2.6"); assert!(ratio(0.035) > 1.1, "the crossover should be well inside the far strike"); }}