quant/src/special.rs
The three special functions the dependence chapter cannot be written without.
//! The three special functions the dependence chapter cannot be written without.//!//! An inverse normal, a log gamma and a regularised incomplete beta, which//! together give the Student-t distribution. Nothing here is novel and all of it//! is standard; it is in the crate because the crate has no dependencies, and it//! is in its own file because it is the one place in the crate where the code is//! numerical recipe rather than finance.//!//! Everything is tested against values that can be checked by hand or against a//! symmetry the function must obey, since an approximation that is quietly wrong//! in the tail would be wrong exactly where the dependence chapter is looking. use crate::black::norm_cdf; /// The inverse of the standard normal CDF.////// Acklam's rational approximation, with one step of Halley refinement against/// [`norm_cdf`] to clean up the last few digits. Accurate to better than 1e-15/// across the range, which matters here: the dependence chapter works in the tail, where a/// relative error in the threshold becomes a large relative error in a/// probability.pub fn norm_inv(p: f64) -> f64 { if !(p > 0.0 && p < 1.0) { return if p <= 0.0 { f64::NEG_INFINITY } else { f64::INFINITY }; } const A: [f64; 6] = [ -3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02, 1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00, ]; const B: [f64; 5] = [ -5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02, 6.680131188771972e+01, -1.328068155288572e+01, ]; const C: [f64; 6] = [ -7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00, -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00, ]; const D: [f64; 4] = [ 7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00, 3.754408661907416e+00, ]; const LOW: f64 = 0.02425; let x = if p < LOW { let q = (-2.0 * p.ln()).sqrt(); (((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5]) / ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0) } else if p <= 1.0 - LOW { let q = p - 0.5; let r = q * q; (((((A[0] * r + A[1]) * r + A[2]) * r + A[3]) * r + A[4]) * r + A[5]) * q / (((((B[0] * r + B[1]) * r + B[2]) * r + B[3]) * r + B[4]) * r + 1.0) } else { let q = (-2.0 * (1.0 - p).ln()).sqrt(); -(((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5]) / ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0) }; // One Halley step. The approximation above is good to about 1e-9; this takes // it to machine precision, and costs one normal CDF. let e = norm_cdf(x) - p; let u = e * (2.0 * std::f64::consts::PI).sqrt() * (x * x / 2.0).exp(); x - u / (1.0 + x * u / 2.0)} /// `ln(Gamma(x))` for `x > 0`, by the Lanczos approximation.pub fn ln_gamma(x: f64) -> f64 { const G: [f64; 9] = [ 0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7, ]; if x < 0.5 { // Reflection, so the series is only ever used where it converges well. (std::f64::consts::PI / (std::f64::consts::PI * x).sin()).ln() - ln_gamma(1.0 - x) } else { let x = x - 1.0; let mut a = G[0]; let t = x + 7.5; for (i, &g) in G.iter().enumerate().skip(1) { a += g / (x + i as f64); } 0.5 * (2.0 * std::f64::consts::PI).ln() + (x + 0.5) * t.ln() - t + a.ln() }} /// The regularised incomplete beta function `I_x(a,b)`.////// By the continued fraction of Lentz, with the standard symmetry applied so the/// fraction is only evaluated where it converges quickly.pub fn inc_beta(x: f64, a: f64, b: f64) -> f64 { if x <= 0.0 { return 0.0; } if x >= 1.0 { return 1.0; } let front = (ln_gamma(a + b) - ln_gamma(a) - ln_gamma(b) + a * x.ln() + b * (1.0 - x).ln()).exp(); if x < (a + 1.0) / (a + b + 2.0) { front * beta_cf(x, a, b) / a } else { 1.0 - front * beta_cf(1.0 - x, b, a) / b }} /// The continued fraction behind [`inc_beta`], evaluated by the modified Lentz/// method.////// The two halves of each iteration are the even and odd coefficients of the/// fraction, which have different forms and so cannot be folded into one loop/// body without getting the first step wrong.fn beta_cf(x: f64, a: f64, b: f64) -> f64 { const TINY: f64 = 1e-300; const EPS: f64 = 1e-15; let (qab, qap, qam) = (a + b, a + 1.0, a - 1.0); let mut c = 1.0; let mut d = 1.0 - qab * x / qap; if d.abs() < TINY { d = TINY; } d = 1.0 / d; let mut h = d; for m in 1..300 { let m = m as f64; let m2 = 2.0 * m; // Even coefficient. let aa = m * (b - m) * x / ((qam + m2) * (a + m2)); d = 1.0 + aa * d; if d.abs() < TINY { d = TINY; } c = 1.0 + aa / c; if c.abs() < TINY { c = TINY; } d = 1.0 / d; h *= d * c; // Odd coefficient. let aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2)); d = 1.0 + aa * d; if d.abs() < TINY { d = TINY; } c = 1.0 + aa / c; if c.abs() < TINY { c = TINY; } d = 1.0 / d; let delta = d * c; h *= delta; if (delta - 1.0).abs() < EPS { break; } } h} /// The CDF of a Student-t with `nu` degrees of freedom.////// Accurate to machine precision in the tails, which is where the dependence chapter uses/// it. Within about `1e-8` of the median it is not: the argument `nu/(nu + x^2)`/// rounds to exactly one once `x^2` falls below the double epsilon, so the/// function has a flat spot of that width sitting on `0.5`. Harmless for/// everything here, and stated because a caller inverting it near the median/// would otherwise be quietly disappointed.pub fn t_cdf(x: f64, nu: f64) -> f64 { let p = 0.5 * inc_beta(nu / (nu + x * x), 0.5 * nu, 0.5); if x > 0.0 { 1.0 - p } else { p }} /// The inverse of [`t_cdf`], by bisection.////// Bisection rather than a rational approximation because it is called a handful/// of times per figure and never in a loop, and because it cannot be wrong in a/// way [`t_cdf`] is not already wrong.pub fn t_inv(p: f64, nu: f64) -> f64 { if !(p > 0.0 && p < 1.0) { return if p <= 0.0 { f64::NEG_INFINITY } else { f64::INFINITY }; } let (mut lo, mut hi) = (-1e6, 1e6); for _ in 0..200 { let mid = 0.5 * (lo + hi); if t_cdf(mid, nu) < p { lo = mid; } else { hi = mid; } if hi - lo < 1e-13 * (1.0 + hi.abs()) { break; } } 0.5 * (lo + hi)} #[cfg(test)]mod tests { use super::*; #[test] fn norm_inv_inverts_norm_cdf() { for i in 1..1000 { let p = i as f64 / 1000.0; let round_trip = norm_cdf(norm_inv(p)); assert!((round_trip - p).abs() < 1e-14, "at p={p} got {round_trip}"); } } #[test] fn norm_inv_is_accurate_deep_in_the_tail() { // Where the dependence chapter actually works. A senior tranche is // priced off events at these probabilities, so an approximation that // gives up at 1e-4 is useless here. for p in [1e-3, 1e-5, 1e-8, 1e-12] { let round_trip = norm_cdf(norm_inv(p)); assert!( (round_trip / p - 1.0).abs() < 1e-9, "at p={p} got {round_trip}" ); } // Known quantiles, to the digits everybody remembers. assert!((norm_inv(0.975) - 1.959963984540054).abs() < 1e-12); assert!((norm_inv(0.99) - 2.326347874040841).abs() < 1e-12); } #[test] fn norm_inv_is_antisymmetric() { for p in [0.001, 0.05, 0.2, 0.45] { assert!((norm_inv(p) + norm_inv(1.0 - p)).abs() < 1e-12); } } #[test] fn ln_gamma_reproduces_the_factorials() { // Gamma(n) = (n-1)! let mut factorial = 1.0f64; for n in 1..15 { let got = ln_gamma(n as f64).exp(); assert!( (got / factorial - 1.0).abs() < 1e-11, "Gamma({n}) gave {got}, wanted {factorial}" ); factorial *= n as f64; } // And the half-integer value, which the reflection branch produces. assert!( (ln_gamma(0.5) - std::f64::consts::PI.sqrt().ln()).abs() < 1e-12, "Gamma(1/2) was not sqrt(pi)" ); } #[test] fn inc_beta_matches_the_cases_that_are_elementary() { // I_x(1,1) = x, and I_x(1,2) = 1-(1-x)^2, both integrable by hand. for i in 1..20 { let x = i as f64 / 20.0; assert!((inc_beta(x, 1.0, 1.0) - x).abs() < 1e-12); let expected = 1.0 - (1.0 - x) * (1.0 - x); assert!((inc_beta(x, 1.0, 2.0) - expected).abs() < 1e-12); } // And the symmetry I_x(a,b) = 1 - I_{1-x}(b,a). for &(x, a, b) in &[(0.3, 2.5, 4.0), (0.7, 0.5, 0.5), (0.1, 9.0, 1.5)] { let lhs = inc_beta(x, a, b); let rhs = 1.0 - inc_beta(1.0 - x, b, a); assert!((lhs - rhs).abs() < 1e-13, "at ({x},{a},{b}): {lhs} vs {rhs}"); } } #[test] fn the_t_distribution_is_symmetric_and_has_the_right_median() { for nu in [1.0, 2.5, 4.0, 30.0] { assert!((t_cdf(0.0, nu) - 0.5).abs() < 1e-13); for x in [0.3, 1.0, 2.5, 6.0] { assert!((t_cdf(x, nu) + t_cdf(-x, nu) - 1.0).abs() < 1e-13); } } } #[test] fn the_cauchy_case_is_the_one_with_a_closed_form() { // nu = 1 is Cauchy, whose CDF is 1/2 + arctan(x)/pi. An independent // formula for one slice of the same function. for x in [-4.0f64, -0.7, 0.0, 0.7, 4.0, 25.0] { let expected = 0.5 + x.atan() / std::f64::consts::PI; assert!( (t_cdf(x, 1.0) - expected).abs() < 1e-12, "at x={x}: {} vs {expected}", t_cdf(x, 1.0) ); } } #[test] fn the_t_approaches_the_normal_as_the_tails_thin() { // The sanity check that the whole of the dependence chapter leans on: a t with many // degrees of freedom is a normal, so a t copula with a high nu is a // Gaussian copula and its extra tail dependence must vanish. for x in [-3.0, -1.0, 0.5, 2.0] { assert!((t_cdf(x, 1e6) - norm_cdf(x)).abs() < 1e-5, "at x={x}"); } } #[test] fn t_inv_inverts_t_cdf() { for nu in [3.0, 5.0, 12.0] { for i in 1..200 { let p = i as f64 / 200.0; let round_trip = t_cdf(t_inv(p, nu), nu); // Away from the median this is good to 1e-13; at the median // itself it is limited by the flat spot documented on t_cdf. let tolerance = if (p - 0.5).abs() < 1e-9 { 1e-7 } else { 1e-11 }; assert!( (round_trip - p).abs() < tolerance, "nu={nu} p={p} gave {round_trip}" ); } } } #[test] fn the_flat_spot_at_the_median_is_where_it_is_claimed_to_be() { // Pinned so the limitation stays a known one. The cancellation in // nu/(nu+x^2) kills the last digits for x below about 1e-8, and nothing // outside that window is affected. assert_eq!(t_cdf(1e-10, 5.0), 0.5); assert!((t_cdf(1e-6, 5.0) - 0.5).abs() > 1e-7, "resolution lost too early"); }}