quant/src/estimation.rs
Estimating parameters from a history, and testing whether a model fits.
//! Estimating parameters from a history, and testing whether a model fits.//!//! Everything before the fitting and testing chapter in these notes takes a//! model's parameters as given. This module is about where they come from when//! they come from data rather than from prices, which is the real-world measure//! of the relative value chapter and a different problem from calibration in//! every respect that matters.//!//! Two results here are worth the whole module. The first is that the estimator//! of a mean reversion rate is badly biased upwards, and that the bias depends//! on the *length* of the history and not on how finely it was sampled --- so//! the usual reflex of pulling higher frequency data does nothing at all. The//! second is that whether a model fits is a question with a general answer, the//! probability integral transform, which works for any model whose conditional//! distribution can be written down and does not care where the model came from. use crate::black::norm_cdf;use crate::pathwise::Rng; // ---------------------------------------------------------------------------// Ornstein-Uhlenbeck: the process every relative value trade is a bet on.// --------------------------------------------------------------------------- /// A mean reverting path, `dX = kappa (theta - X) dt + sigma dW`.////// Sampled exactly rather than by an Euler step: the transition law of an/// Ornstein-Uhlenbeck process is known in closed form, so there is no reason to/// introduce a discretisation error that would then be confused with the/// estimation error this module is about.pub fn ou_path( kappa: f64, theta: f64, sigma: f64, start: f64, dt: f64, steps: usize, rng: &mut Rng,) -> Vec<f64> { let decay = (-kappa * dt).exp(); // Stationary-conditional variance of the exact transition. let variance = if kappa > 0.0 { sigma * sigma * (1.0 - decay * decay) / (2.0 * kappa) } else { sigma * sigma * dt }; let sd = variance.sqrt(); let mut path = Vec::with_capacity(steps + 1); let mut x = start; path.push(x); for _ in 0..steps { x = theta + (x - theta) * decay + sd * rng.next_normal(); path.push(x); } path} /// What an Ornstein-Uhlenbeck fit returns.pub struct OuFit { pub kappa: f64, pub theta: f64, pub sigma: f64,} /// Fit an Ornstein-Uhlenbeck process to an observed path by matching moments/// to their stationary values, rather than by regression.////// The stationary mean is `theta`, the stationary variance is/// `sigma^2 / (2 kappa)`, and the lag-one autocovariance at spacing `dt` is the/// stationary variance times `exp(-kappa dt)` --- three moments computed/// directly from the path, inverted for the three parameters. Classical/// (Yule-Walker, generalised to continuous time), and a genuinely different/// estimator from [`fit_ou`]: the two agree in the limit, and there is no/// reason to expect them to agree away from it, since it is not the same/// computation on the same numbers.pub fn fit_ou_moments(path: &[f64], dt: f64) -> OuFit { let n = path.len(); if n < 3 || dt <= 0.0 { return OuFit { kappa: f64::NAN, theta: f64::NAN, sigma: f64::NAN }; } let count = n as f64; let mean = path.iter().sum::<f64>() / count; let variance = path.iter().map(|&x| (x - mean).powi(2)).sum::<f64>() / count; let autocovariance = path.windows(2).map(|w| (w[0] - mean) * (w[1] - mean)).sum::<f64>() / (count - 1.0); if variance <= 0.0 { return OuFit { kappa: f64::NAN, theta: mean, sigma: f64::NAN }; } let ratio = autocovariance / variance; if ratio <= 0.0 { // No decay at all in the sample autocovariance: as with fit_ou, report // the absence of mean reversion rather than a kappa that isn't there. return OuFit { kappa: f64::INFINITY, theta: mean, sigma: f64::NAN }; } let kappa = -ratio.min(1.0).ln() / dt; let sigma = (2.0 * kappa * variance).sqrt(); OuFit { kappa, theta: mean, sigma }} /// What a Student-t fit of the Ornstein-Uhlenbeck regression returns: the same/// three parameters as [`OuFit`], plus the degrees of freedom the innovations/// were fitted with.pub struct StudentOuFit { pub kappa: f64, pub theta: f64, pub sigma: f64, pub nu: f64,} /// The AR(1) regression coefficients, scale and log-likelihood at a fixed/// `nu`, by iteratively reweighted least squares.////// Ordinary least squares is the maximum likelihood fit under a Gaussian/// innovation; this is the maximum likelihood fit under a Student-t one. The/// two differ only in the weight each residual carries: OLS gives every/// observation weight one, this gives observation `i` weight/// `(nu+1)/(nu + (e_i/s)^2)`, which is small exactly where the residual is/// large --- a fat-tailed shock is down-weighted rather than allowed to drag/// the fitted line towards it, which is the whole content of "robust".fn fit_ar1_student_t(x: &[f64], y: &[f64], nu: f64) -> (f64, f64, f64, f64) { let n = x.len() as f64; let mean_x = x.iter().sum::<f64>() / n; let mean_y = y.iter().sum::<f64>() / n; let (mut sxx, mut sxy) = (0.0, 0.0); for (&xi, &yi) in x.iter().zip(y) { let dx = xi - mean_x; sxx += dx * dx; sxy += dx * (yi - mean_y); } let mut phi = sxy / sxx; let mut c = mean_y - phi * mean_x; let mut resid: Vec<f64> = x.iter().zip(y).map(|(&xi, &yi)| yi - c - phi * xi).collect(); let mut s2 = resid.iter().map(|e| e * e).sum::<f64>() / n; for _ in 0..300 { let s = s2.sqrt().max(1e-12); let weights: Vec<f64> = resid.iter().map(|&e| (nu + 1.0) / (nu + (e / s).powi(2))).collect(); let sw: f64 = weights.iter().sum(); let (mut wx, mut wy) = (0.0, 0.0); for i in 0..x.len() { wx += weights[i] * x[i]; wy += weights[i] * y[i]; } wx /= sw; wy /= sw; let (mut swxx, mut swxy) = (0.0, 0.0); for i in 0..x.len() { let (dx, dy) = (x[i] - wx, y[i] - wy); swxx += weights[i] * dx * dx; swxy += weights[i] * dx * dy; } let new_phi = swxy / swxx; let new_c = wy - new_phi * wx; let new_resid: Vec<f64> = x.iter().zip(y).map(|(&xi, &yi)| yi - new_c - new_phi * xi).collect(); let new_s2 = weights.iter().zip(&new_resid).map(|(&w, &e)| w * e * e).sum::<f64>() / n; let moved = (new_phi - phi).abs() + (new_c - c).abs() + (new_s2 - s2).abs(); phi = new_phi; c = new_c; s2 = new_s2; resid = new_resid; if moved < 1e-13 * (1.0 + s2) { break; } } let s = s2.sqrt(); let ll = resid.iter().map(|&e| crate::dependence::ln_t_pdf(e / s, nu) - s.ln()).sum::<f64>(); (phi, c, s, ll)} /// Fit an Ornstein-Uhlenbeck process to an observed path by maximum/// likelihood, allowing the innovations to be Student-t rather than assuming/// them Gaussian.////// [`fit_ou`] is exactly this with `nu` fixed at infinity. Letting the data/// choose `nu` costs nothing when the innovations really are Gaussian, and/// buys a fit that is not dragged around by the occasional large shock when/// they are not --- which \chapref{dependence} gives every reason to expect/// of a real spread.pub fn fit_ou_student_t(path: &[f64], dt: f64) -> StudentOuFit { let n = path.len().saturating_sub(1); if n < 5 || dt <= 0.0 { return StudentOuFit { kappa: f64::NAN, theta: f64::NAN, sigma: f64::NAN, nu: f64::NAN }; } let (x, y) = (&path[..n], &path[1..]); let nu = crate::dependence::maximise_over_nu(|nu| fit_ar1_student_t(x, y, nu).3); let (phi, c, s, _) = fit_ar1_student_t(x, y, nu); if phi <= 0.0 { return StudentOuFit { kappa: f64::INFINITY, theta: c, sigma: f64::NAN, nu }; } let kappa = -phi.ln() / dt; let theta = c / (1.0 - phi); // The scale s is the t-distribution's own scale, not its standard // deviation; the variance of a standard t(nu) is nu/(nu-2), finite only // for nu > 2. Converting to the variance fit_ou's sigma is stated in // keeps the two directly comparable. let residual_variance = if nu > 2.0 { s * s * nu / (nu - 2.0) } else { f64::INFINITY }; let sigma = (residual_variance * 2.0 * kappa / (1.0 - phi * phi)).sqrt(); StudentOuFit { kappa, theta, sigma, nu }} /// Fit an Ornstein-Uhlenbeck process to an observed path, by least squares.////// The discrete transition is an AR(1),////// ```text/// X[i+1] = c + phi X[i] + noise, phi = exp(-kappa dt),/// ```////// so regressing each observation on the previous one recovers everything. This/// is also the maximum likelihood estimator, because the transition is Gaussian/// with constant variance and least squares and MLE coincide there.////// It is the estimator every desk uses, and the fitting and testing chapter is/// largely about the trouble it causes.pub fn fit_ou(path: &[f64], dt: f64) -> OuFit { let n = path.len().saturating_sub(1); if n < 2 || dt <= 0.0 { return OuFit { kappa: f64::NAN, theta: f64::NAN, sigma: f64::NAN }; } let (x, y) = (&path[..n], &path[1..]); let count = n as f64; let mean_x = x.iter().sum::<f64>() / count; let mean_y = y.iter().sum::<f64>() / count; let mut sxx = 0.0; let mut sxy = 0.0; for i in 0..n { let dx = x[i] - mean_x; sxx += dx * dx; sxy += dx * (y[i] - mean_y); } if sxx <= 0.0 { return OuFit { kappa: f64::NAN, theta: f64::NAN, sigma: f64::NAN }; } let phi = sxy / sxx; let intercept = mean_y - phi * mean_x; // Residual variance, which inverts to the diffusion coefficient. let mut residual = 0.0; for i in 0..n { let e = y[i] - intercept - phi * x[i]; residual += e * e; } let residual_variance = residual / (count - 2.0); // A path that failed to mean revert at all comes back as phi >= 1, which has // no kappa behind it. Reported as infinite rather than silently clamped. if phi <= 0.0 { return OuFit { kappa: f64::INFINITY, theta: mean_x, sigma: f64::NAN }; } let kappa = -phi.ln() / dt; let theta = intercept / (1.0 - phi); let sigma = if kappa > 0.0 && phi < 1.0 { (residual_variance * 2.0 * kappa / (1.0 - phi * phi)).sqrt() } else { (residual_variance / dt).sqrt() }; OuFit { kappa, theta, sigma }} /// The average estimated mean reversion over many independent histories.////// `window` is the length of each history in years and `per_year` how often it/// is sampled, so the number of observations is their product. The fitting and/// testing chapter's claim is that the answer depends on `window` and/// essentially not on `per_year`, which is the opposite of the usual instinct.pub fn mean_reversion_estimate( true_kappa: f64, window: f64, per_year: f64, trials: usize, seed: u64,) -> f64 { let dt = 1.0 / per_year; let steps = (window * per_year).round() as usize; if steps < 3 { return f64::NAN; } let mut rng = Rng::new(seed); let mut total = 0.0; let mut counted = 0.0; for _ in 0..trials { // Started in the stationary distribution, so nothing here is a burn-in // artefact. let stationary_sd = 1.0 / (2.0 * true_kappa).sqrt(); let start = stationary_sd * rng.next_normal(); let path = ou_path(true_kappa, 0.0, 1.0, start, dt, steps, &mut rng); let fit = fit_ou(&path, dt); if fit.kappa.is_finite() { total += fit.kappa; counted += 1.0; } } total / counted} // ---------------------------------------------------------------------------// Does the model fit?// --------------------------------------------------------------------------- /// The probability integral transform of a sample under a proposed model.////// If `cdf` really is the conditional distribution the observations were drawn/// from, then applying it to them gives numbers that are uniform on `[0,1]` and/// independent of each other. That is true whatever the model is --- a diffusion,/// a jump process, a neural network --- which makes it the one goodness of fit/// test in these notes that does not depend on what is being tested.////// The fitting and testing chapter uses it in both directions: on a series that/// really did come from the proposed model, and on one that did not.pub fn probability_integral_transform(sample: &[f64], cdf: impl Fn(f64) -> f64) -> Vec<f64> { sample.iter().map(|&x| cdf(x).clamp(0.0, 1.0)).collect()} /// The Kolmogorov-Smirnov distance of a sample from the uniform distribution.////// The largest gap between the empirical distribution of the transformed sample/// and the straight line it should lie on. Under the null it shrinks like/// `1/sqrt(n)`, so the scaled statistic below is the one to compare against a/// fixed threshold.pub fn ks_uniform(values: &[f64]) -> f64 { let n = values.len(); if n == 0 { return f64::NAN; } let mut sorted = values.to_vec(); sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)); let mut worst: f64 = 0.0; for (i, &v) in sorted.iter().enumerate() { let below = i as f64 / n as f64; let above = (i + 1) as f64 / n as f64; worst = worst.max((v - below).abs()).max((above - v).abs()); } worst} /// The Kolmogorov-Smirnov statistic scaled by `sqrt(n)`.////// Comparable across sample sizes. The asymptotic five percent critical value is/// `1.358`, so anything much beyond that is a model the data rejects.pub fn ks_statistic(values: &[f64]) -> f64 { ks_uniform(values) * (values.len() as f64).sqrt()} /// The transformed log returns of a path, under a proposed constant-volatility/// lognormal model.////// The workhorse of the fitting and testing chapter's model test. Under the/// proposal, each log return is normal with a known mean and standard/// deviation, so the normal CDF applied to it should give a uniform. Deviations/// from uniformity are the model being wrong, and their *shape* says how.pub fn lognormal_pit(path: &[f64], dt: f64, mu: f64, sigma: f64) -> Vec<f64> { if path.len() < 2 || sigma <= 0.0 || dt <= 0.0 { return Vec::new(); } let mean = (mu - 0.5 * sigma * sigma) * dt; let sd = sigma * dt.sqrt(); let returns: Vec<f64> = path.windows(2).map(|w| (w[1] / w[0]).ln()).collect(); probability_integral_transform(&returns, |r| norm_cdf((r - mean) / sd))} /// A path with genuinely mean-reverting stochastic volatility, Heston's own/// variance SDE, `dv = kappa(theta - v) dt + xi sqrt(v) dW`, full truncation/// Euler so a negative step is floored at zero rather than breaking the square/// root.////// Independent of the return's own driver, since the shape argument below is/// about a moving variance as such, not about the leverage correlation.fn cir_path(kappa: f64, theta: f64, xi: f64, mu: f64, dt: f64, n: usize, seed: u64) -> Vec<f64> { let mut rng = Rng::new(seed); let mut path = Vec::with_capacity(n + 1); path.push(100.0); let mut v = theta; for _ in 0..n { let vpos = v.max(0.0); let last = *path.last().unwrap(); let r = (mu - 0.5 * vpos) * dt + vpos.sqrt() * dt.sqrt() * rng.next_normal(); path.push(last * r.exp()); v = vpos + kappa * (theta - vpos) * dt + xi * vpos.sqrt() * dt.sqrt() * rng.next_normal(); } path} /// The shape of the probability integral transform of a mean-reverting/// stochastic-volatility path, tested against the constant volatility that/// matches its own realised variance exactly.////// Three numbers beyond the histogram matter here, because the deviation from/// uniform is not one bump but three: an excess at the centre that a mixture's/// typical volatility sits below the level fitted to its mean square (Jensen,/// unconditionally); a genuine excess right at the extreme tails, from the/// minority of days loud enough to reach them, invisible at any bin width/// wider than the tail itself; and a deficit in the shoulders between the two./// `tail_below` and `tail_above` count observations beyond the 1st and 99th/// percentile directly, rather than relying on a histogram bin to catch them.pub struct PitShape { /// `bins` densities, normalised so a perfectly uniform sample reads as 1 /// everywhere. pub density: Vec<f64>, pub ks: f64, pub tail_below: usize, pub tail_above: usize, /// What either tail count would be under the null, `n / 100`. pub tail_expected: f64,} pub fn cir_pit_shape( kappa: f64, theta: f64, xi: f64, mu: f64, dt: f64, n: usize, bins: usize, seed: u64,) -> PitShape { let path = cir_path(kappa, theta, xi, mu, dt, n, seed); let returns: Vec<f64> = path.windows(2).map(|w| (w[1] / w[0]).ln()).collect(); let realised = (returns.iter().map(|r| r * r).sum::<f64>() / returns.len() as f64 / dt).sqrt(); let pit = lognormal_pit(&path, dt, mu, realised); let mut counts = vec![0usize; bins]; for &u in &pit { let b = ((u * bins as f64) as usize).min(bins - 1); counts[b] += 1; } let scale = bins as f64 / pit.len() as f64; let density = counts.iter().map(|&c| c as f64 * scale).collect(); let tail_below = pit.iter().filter(|&&u| u < 0.01).count(); let tail_above = pit.iter().filter(|&&u| u > 0.99).count(); PitShape { density, ks: ks_statistic(&pit), tail_below, tail_above, tail_expected: pit.len() as f64 / 100.0, }} #[cfg(test)]mod tests { use super::*; /// A standard Student-t draw via a normal over the root of a scaled /// chi-square, matching the construction used in the dependence module. fn student_t_draw(rng: &mut Rng, nu: f64) -> f64 { let z = rng.next_normal(); let mut w = 0.0; // nu need not be an integer; approximate the chi-square with the // nearest integer degrees of freedom, close enough for a simulation. let dof = nu.round().max(2.0) as usize; for _ in 0..dof { w += rng.next_normal().powi(2); } z / (w / dof as f64).sqrt() } fn ou_path_student_t( kappa: f64, theta: f64, sigma: f64, nu: f64, start: f64, dt: f64, steps: usize, rng: &mut Rng, ) -> Vec<f64> { let decay = (-kappa * dt).exp(); let variance = sigma * sigma * (1.0 - decay * decay) / (2.0 * kappa); // Scaled so the innovation has the same variance as the Gaussian case, // for a fair comparison at matched sigma. let scale = (variance * (nu - 2.0) / nu).sqrt(); let mut path = Vec::with_capacity(steps + 1); let mut x = start; path.push(x); for _ in 0..steps { x = theta + (x - theta) * decay + scale * student_t_draw(rng, nu); path.push(x); } path } /// The RMSE of `fit(path).kappa` against the truth, over many independent /// paths with Student-t innovations at the given `nu` (large `nu` is /// effectively Gaussian). fn kappa_rmse( nu: f64, seed_base: u64, fit: fn(&[f64], f64) -> OuFit, ) -> f64 { let (true_kappa, true_theta, true_sigma) = (1.0f64, 0.05, 0.2); let dt = 1.0 / 252.0; let steps = (5.0 * 252.0) as usize; let trials = 800; let mut err = 0.0; for tr in 0..trials { let mut rng = Rng::new(seed_base + tr as u64); let path = ou_path_student_t( true_kappa, true_theta, true_sigma, nu, true_theta, dt, steps, &mut rng, ); err += (fit(&path, dt).kappa - true_kappa).powi(2); } (err / trials as f64).sqrt() } #[test] fn student_t_mle_is_more_accurate_than_ols_when_innovations_are_fat_tailed() { let ols = kappa_rmse(4.0, 5_000_000, fit_ou); let t_mle = kappa_rmse(4.0, 5_000_000, |p, dt| { let f = fit_ou_student_t(p, dt); OuFit { kappa: f.kappa, theta: f.theta, sigma: f.sigma } }); assert!( t_mle < 0.85 * ols, "t-MLE RMSE {t_mle} was not clearly below OLS's {ols}" ); } #[test] fn student_t_mle_recovers_the_degrees_of_freedom_and_costs_nothing_under_gaussian_innovations() { let (true_kappa, true_theta, true_sigma, true_nu) = (1.0f64, 0.05, 0.2, 4.0); let dt = 1.0 / 252.0; let steps = (5.0 * 252.0) as usize; let mut nu_sum = 0.0; for tr in 0..800u64 { let mut rng = Rng::new(5_000_000 + tr); let path = ou_path_student_t( true_kappa, true_theta, true_sigma, true_nu, true_theta, dt, steps, &mut rng, ); nu_sum += fit_ou_student_t(&path, dt).nu; } let mean_nu = nu_sum / 800.0; assert!((mean_nu - true_nu).abs() < 0.5, "mean fitted nu {mean_nu} against truth {true_nu}"); // Under near-Gaussian innovations the two estimators should be almost // indistinguishable: nothing was given up by allowing for fat tails // that were not there. let ols_gaussian = kappa_rmse(200.0, 5_000_000, fit_ou); let t_mle_gaussian = kappa_rmse(200.0, 5_000_000, |p, dt| { let f = fit_ou_student_t(p, dt); OuFit { kappa: f.kappa, theta: f.theta, sigma: f.sigma } }); assert!( (t_mle_gaussian - ols_gaussian).abs() < 0.05 * ols_gaussian, "t-MLE {t_mle_gaussian} against OLS {ols_gaussian} under Gaussian innovations" ); } #[test] fn method_of_moments_is_biased_the_same_way_ols_is() { // A genuinely different computation -- no regression anywhere in it -- // and yet it fails the same way, by about the same amount, closing at // about the same rate. That is the stronger claim: the upward bias is // a property of the estimation problem, not an artefact of least // squares in particular. let true_kappa = 1.0f64; let dt = 1.0 / 252.0; let trials = 4000; let mean_kappa = |window: f64, fit: fn(&[f64], f64) -> OuFit| -> f64 { let steps = (window * 252.0).round() as usize; let mut total = 0.0; for t in 0..trials { let mut rng = Rng::new(9_000_000 + t as u64); let stationary_sd = 0.2 / (2.0 * true_kappa).sqrt(); let start = 0.05 + stationary_sd * rng.next_normal(); let path = ou_path(true_kappa, 0.05, 0.2, start, dt, steps, &mut rng); total += fit(&path, dt).kappa; } total / trials as f64 }; for window in [2.0f64, 10.0, 40.0] { let ols = mean_kappa(window, fit_ou); let mom = mean_kappa(window, fit_ou_moments); assert!(ols > true_kappa * 1.05, "OLS not biased up at window={window}: {ols}"); assert!(mom > true_kappa * 1.05, "MoM not biased up at window={window}: {mom}"); // Neither is more than about fifteen percent away from the other: // the same problem, not a fix. assert!( (mom - ols).abs() < 0.15 * ols, "window={window}: OLS={ols} MoM={mom}, further apart than expected" ); } // And both close towards the truth at a comparable rate as the window // lengthens, rather than one flattening out while the other keeps // improving. let ols_40 = mean_kappa(40.0, fit_ou) - true_kappa; let mom_40 = mean_kappa(40.0, fit_ou_moments) - true_kappa; assert!(ols_40 < 0.15 && mom_40 < 0.15, "OLS={ols_40} MoM={mom_40} still large at 40y"); } #[test] fn the_fit_recovers_its_own_parameters_given_enough_history() { // Consistency: the estimator is not wrong, it is biased, and the two are // different complaints. Given a long enough window it converges. let mut rng = Rng::new(31_415); let path = ou_path(1.5, 0.04, 0.20, 0.04, 1.0 / 252.0, 252 * 400, &mut rng); let fit = fit_ou(&path, 1.0 / 252.0); assert!((fit.kappa - 1.5).abs() < 0.10, "kappa {}", fit.kappa); assert!((fit.theta - 0.04).abs() < 0.01, "theta {}", fit.theta); assert!((fit.sigma - 0.20).abs() < 0.01, "sigma {}", fit.sigma); } #[test] fn mean_reversion_is_overestimated_on_a_short_history() { // The result the fitting and testing chapter is built around, and the // one that costs money: the estimate says the spread snaps back faster // than it does, so the trade is sized larger and held with more // confidence than the evidence supports. let truth = 1.0; let estimated = mean_reversion_estimate(truth, 2.0, 252.0, 4_000, 20_260_804); assert!( estimated > 1.5 * truth, "a two year window gave {estimated} against a true {truth}" ); } #[test] fn sampling_more_often_does_not_fix_it() { // The part that surprises people. The bias is a function of how many // mean reversion times the window covers, and sampling the same window // more finely adds observations without adding any of those. let truth = 1.0; let daily = mean_reversion_estimate(truth, 4.0, 252.0, 4_000, 11_111); let hourly = mean_reversion_estimate(truth, 4.0, 252.0 * 8.0, 4_000, 22_222); assert!( (daily - hourly).abs() < 0.25 * (daily - truth), "daily {daily}, hourly {hourly}: sampling changed the bias materially" ); assert!(daily > 1.2, "daily {daily}"); assert!(hourly > 1.2, "hourly {hourly}"); } #[test] fn a_longer_window_does_fix_it() { // And the corresponding statement in the other direction, which is the // only lever that works. let truth = 1.0; let mut previous = f64::INFINITY; for window in [2.0, 5.0, 10.0, 20.0, 40.0] { let estimated = mean_reversion_estimate(truth, window, 252.0, 3_000, 777); assert!(estimated < previous, "bias not falling at window={window}"); previous = estimated; } assert!(previous < 1.15, "even forty years left a bias of {previous}"); } #[test] fn the_bias_scales_like_one_over_the_window() { // The sharp version, which is what makes the two tests above a law // rather than an observation: the bias in kappa is about 4/T, with T the // window in years, and the 4 does not depend on the true kappa. for truth in [0.5, 1.0, 2.0] { for window in [10.0, 20.0] { let estimated = mean_reversion_estimate(truth, window, 252.0, 6_000, 4_242); let scaled = (estimated - truth) * window; assert!( (scaled - 4.0).abs() < 1.6, "kappa={truth} window={window}: bias*T was {scaled}, not about 4" ); } } } #[test] fn the_transform_of_a_correct_model_is_uniform() { // The null case. Data really from the model, tested against the model, // passes -- so a rejection below means something. let mut rng = Rng::new(90_210); let (dt, mu, sigma) = (1.0f64 / 252.0, 0.05f64, 0.20f64); let mut path = vec![100.0]; for _ in 0..4000 { let last = *path.last().unwrap(); let r = (mu - 0.5 * sigma * sigma) * dt + sigma * dt.sqrt() * rng.next_normal(); path.push(last * r.exp()); } let pit = lognormal_pit(&path, dt, mu, sigma); let statistic = ks_statistic(&pit); assert!(statistic < 1.358, "correct model rejected, statistic {statistic}"); } #[test] fn a_realistic_mean_reverting_volatility_gives_three_bumps_not_one() { // Heston's own calibration, used elsewhere in these notes: kappa=0.3, // theta=0.02, xi=0.3. Fitted with exactly the right average variance, // so nothing about the level is wrong -- only the shape. let shape = cir_pit_shape(0.3, 0.02, 0.3, 0.05, 1.0 / 252.0, 4_000, 10, 20260925); assert!(shape.ks > 1.358, "not rejected at all, KS={}", shape.ks); // The centre: bins 4 and 5 (u in [0.4, 0.6]) well above uniform. assert!(shape.density[4] > 1.5 && shape.density[5] > 1.5, "{:?}", shape.density); // The shoulders: depleted relative to uniform, on both sides. assert!(shape.density[1] < 1.0 && shape.density[2] < 1.0, "{:?}", shape.density); assert!(shape.density[7] < 1.0 && shape.density[8] < 1.0, "{:?}", shape.density); // The tails: a genuine excess beyond the 1st and 99th percentile, // invisible in the coarse bins above but real at the resolution that // actually catches it. assert!( shape.tail_below as f64 > 1.5 * shape.tail_expected, "below 1%: {} against {} expected", shape.tail_below, shape.tail_expected ); assert!( shape.tail_above as f64 > 1.5 * shape.tail_expected, "above 99%: {} against {} expected", shape.tail_above, shape.tail_expected ); } #[test] fn faster_mean_reversion_shrinks_the_centre_and_the_tails_together() { // The question worth asking rather than assuming: does speeding up the // mean reversion trade the centre bump for a cleaner U, or does it // shrink both symmetrically towards uniform? A faster-reverting CIR // has a tighter stationary variance for the same vol-of-vol, so both // the Jensen gap at the centre and the mass reaching the tail should // fall together, not swap places. Same seed throughout, so kappa is // the only thing that differs between runs. let (theta, xi, mu, dt, n, seed) = (0.02f64, 0.3f64, 0.05f64, 1.0f64 / 252.0, 4_000usize, 20260925u64); let slow = cir_pit_shape(0.1, theta, xi, mu, dt, n, 10, seed); let fast = cir_pit_shape(10.0, theta, xi, mu, dt, n, 10, seed); let centre = |s: &PitShape| s.density[4].max(s.density[5]); let tail = |s: &PitShape| (s.tail_below.max(s.tail_above)) as f64 / s.tail_expected; assert!( centre(&fast) < centre(&slow), "centre did not shrink: slow={:.2} fast={:.2}", centre(&slow), centre(&fast) ); assert!( tail(&fast) < tail(&slow), "tail excess did not shrink: slow={:.2} fast={:.2}", tail(&slow), tail(&fast) ); // And with reversion fast enough, both are close to gone -- the shape // converges towards uniform rather than towards a clean U. assert!(centre(&fast) < 1.3, "centre still elevated at fast reversion: {}", centre(&fast)); assert!(tail(&fast) < 2.0, "tail still elevated at fast reversion: {}", tail(&fast)); } #[test] fn the_transform_detects_a_volatility_that_moves() { // The alternative. The same test on data whose volatility is stochastic, // fitted with the *correct average* volatility so that nothing is wrong // with the level -- only with the shape. let mut rng = Rng::new(13_579); let dt = 1.0f64 / 252.0; let (mu, base) = (0.05f64, 0.20f64); let mut path = vec![100.0]; let mut sigma = base; for _ in 0..4000 { sigma *= (-0.5 * 0.04 * dt + 1.2 * dt.sqrt() * rng.next_normal()).exp(); let last = *path.last().unwrap(); let r = (mu - 0.5 * sigma * sigma) * dt + sigma * dt.sqrt() * rng.next_normal(); path.push(last * r.exp()); } // Fit the constant volatility that best matches the realised variance, // so the model is being given every chance. let returns: Vec<f64> = path.windows(2).map(|w| (w[1] / w[0]).ln()).collect(); let realised = (returns.iter().map(|r| r * r).sum::<f64>() / returns.len() as f64 / dt).sqrt(); let pit = lognormal_pit(&path, dt, mu, realised); let statistic = ks_statistic(&pit); assert!( statistic > 1.358, "stochastic volatility was not detected, statistic {statistic}" ); } #[test] fn the_ks_distance_behaves_on_cases_that_can_be_checked_by_hand() { // A perfectly spread sample is nearly on the line, and a sample crammed // into half the interval is half a unit away from it. let even: Vec<f64> = (0..1000).map(|i| (i as f64 + 0.5) / 1000.0).collect(); assert!(ks_uniform(&even) < 1e-3, "{}", ks_uniform(&even)); let crammed: Vec<f64> = (0..1000).map(|i| (i as f64 + 0.5) / 2000.0).collect(); assert!( (ks_uniform(&crammed) - 0.5).abs() < 1e-2, "{}", ks_uniform(&crammed) ); }} /// A convergence trade on a mean-reverting spread, and whether it survives.////// The relative value chapter's thesis is that a mispricing is tradeable only/// with a structural reason *and* a horizon. This is the horizon half, made/// quantitative.////// Enter when an Ornstein-Uhlenbeck spread sits `entry` stationary standard/// deviations from its mean, betting it returns. Two things can happen first: the/// spread converges, or it widens to `stop` deviations and the position is closed/// at a loss. The trade is right about direction and can still lose, and how/// often is a first-passage question rather than a forecasting one.pub struct ConvergenceTrade { /// Mean reversion speed. The half-life is `ln 2 / kappa`. pub kappa: f64, pub sigma: f64, /// Entry distance, in stationary standard deviations. pub entry: f64, /// Stop-out distance, in the same units. Must exceed `entry`. pub stop: f64,} /// What a convergence trade did.pub struct TradeOutcome { /// Fraction of trades stopped out before converging. pub stopped: f64, /// Mean time to converge, over the trades that converged. pub mean_time: f64, /// Mean worst adverse excursion, in stationary standard deviations, over the /// trades that converged --- the drawdown a survivor had to sit through. pub mean_excursion: f64,} impl ConvergenceTrade { /// The stationary standard deviation of the spread, `sigma / sqrt(2 kappa)`. pub fn stationary_sd(&self) -> f64 { self.sigma / (2.0 * self.kappa).sqrt() } pub fn half_life(&self) -> f64 { std::f64::consts::LN_2 / self.kappa } /// Simulate the trade to its conclusion. pub fn run(&self, paths: usize, steps_per_year: usize, seed: u64) -> TradeOutcome { use crate::pathwise::Rng; let sd = self.stationary_sd(); let (start, barrier) = (self.entry * sd, self.stop * sd); let dt = 1.0 / steps_per_year as f64; let root_dt = dt.sqrt(); let mut rng = Rng::new(seed); let (mut stopped, mut converged) = (0usize, 0usize); let (mut total_time, mut total_excursion) = (0.0, 0.0); for _ in 0..paths { let mut x = start; let mut worst = start; let mut t = 0.0; loop { x += -self.kappa * x * dt + self.sigma * root_dt * rng.next_normal(); t += dt; worst = worst.max(x); if x >= barrier { stopped += 1; break; } if x <= 0.0 { converged += 1; total_time += t; total_excursion += worst / sd; break; } } } TradeOutcome { stopped: stopped as f64 / paths as f64, mean_time: if converged > 0 { total_time / converged as f64 } else { f64::NAN }, mean_excursion: if converged > 0 { total_excursion / converged as f64 } else { f64::NAN }, } }} #[cfg(test)]mod convergence_tests { use super::*; /// A spread with a one-year half-life, entered two standard deviations wide. fn trade(stop: f64) -> ConvergenceTrade { ConvergenceTrade { kappa: std::f64::consts::LN_2, sigma: 0.01, entry: 2.0, stop } } const PATHS: usize = 40_000; const SEED: u64 = 20260808; #[test] fn a_tight_stop_loses_a_trade_that_was_right() { // Being right about direction is not enough. A spread two deviations wide // and mean reverting with certainty still widens further a great deal of // the time, so a stop placed close to the entry is hit more often than // not before the spread converges. let close = trade(2.5).run(PATHS, 500, SEED); assert!( close.stopped > 0.3, "a stop half a deviation away should be hit on well over a quarter of \ trades, got {:.3}", close.stopped ); // Widening the stop converts those losses into survivals, monotonically. let mut previous = close.stopped; for stop in [3.0, 4.0, 6.0] { let outcome = trade(stop).run(PATHS, 500, SEED); assert!(outcome.stopped < previous, "a wider stop should be hit less often"); previous = outcome.stopped; } // Past about four deviations the stop stops mattering: the spread almost // never gets there, so the trade's risk is the holding period rather than // the loss. That is the regime the chapter argues a desk should be in, // and it is only reachable if the position is small enough to sit // through twice the entry width. assert!(trade(4.0).run(PATHS, 500, SEED).stopped < 0.02); } #[test] fn the_survivors_still_sat_through_a_drawdown() { // The number that sizes the position. Conditioning on the trades that did // converge, the spread first went further against them, and by an amount // that is not small next to the entry level. let outcome = trade(6.0).run(PATHS, 500, SEED); assert!( outcome.mean_excursion > 2.3, "survivors entered at 2 and should have seen worse, got {:.2}", outcome.mean_excursion ); } #[test] fn convergence_takes_longer_than_the_half_life() { // The horizon, and the trap. A half-life describes how an *expectation* // decays, not how long a path takes to reach the mean, and the two are // not close: first passage to the mean from two deviations takes well over // a half-life on average. // // So a trade sized on the half-life is sized on the wrong number, and it // is wrong in the dangerous direction -- the holding period is longer // than the estimate suggests. let t = trade(6.0); let outcome = t.run(PATHS, 500, SEED); assert!( outcome.mean_time > 2.0 * t.half_life(), "mean first passage {:.3} against a half-life of {:.3}", outcome.mean_time, t.half_life() ); } #[test] fn an_overstated_mean_reversion_understates_the_holding_period() { // The chain this chapter closes. The fitting chapter measures that // estimating mean reversion from a finite sample overstates it, and the // solvable models chapter explains why from the spectrum. An overstated // kappa is an understated half-life, so a desk plans for a shorter trade // than it gets. // // Here the true half-life is one year and the estimate is taken to be // thirty per cent fast, which is well inside what a decade of data // produces. The planned horizon and the realised one are compared. let truth = trade(6.0); let believed = ConvergenceTrade { kappa: truth.kappa * 1.3, ..trade(6.0) }; let realised = truth.run(PATHS, 500, SEED).mean_time; let planned = believed.run(PATHS, 500, SEED).mean_time; assert!( realised > 1.2 * planned, "realised {realised:.3} should overrun the planned {planned:.3}" ); // And the error compounds with the bias rather than washing out. let worse = ConvergenceTrade { kappa: truth.kappa * 1.6, ..trade(6.0) }; assert!(worse.run(PATHS, 500, SEED).mean_time < planned); }} /// A convergence trade's profit rate, net of the cost of getting in and out.////// The relative value chapter's horizon analysis is frictionless. A round-trip/// cost does not scale the return down; it removes trades, and this measures the/// threshold.////// Entering at `entry` deviations, the trade collects `entry` on convergence and/// pays `stop - entry` when stopped, minus `cost` either way. Since the edge grows/// with the entry width and the cost does not, there is a minimum width below/// which no amount of mean reversion pays --- and it is close to the cost itself.////// Returns expected profit per year of holding, in stationary standard deviations./// Note what is *not* in the objective: the time spent waiting for the spread to/// reach `entry` in the first place. Including it would produce an interior/// optimum in the entry level; without it the rate rises with width indefinitely,/// so this quantity is useful for locating the breakeven and not for choosing an/// entry.pub fn profit_rate( kappa: f64, entry: f64, stop: f64, cost: f64, paths: usize, seed: u64,) -> f64 { let trade = ConvergenceTrade { kappa, sigma: 0.01, entry, stop }; let outcome = trade.run(paths, 500, seed); let expected_gain = (1.0 - outcome.stopped) * entry - outcome.stopped * (stop - entry) - cost; expected_gain / outcome.mean_time} #[cfg(test)]mod cost_tests { use super::*; const KAPPA: f64 = std::f64::consts::LN_2; const PATHS: usize = 20_000; const SEED: u64 = 20260809; /// The narrowest entry, to a quarter of a deviation, at which the trade pays. fn breakeven_entry(cost: f64) -> f64 { let mut e = 0.25; while e <= 6.0 { if profit_rate(KAPPA, e, 8.0, cost, PATHS, SEED) > 0.0 { return e; } e += 0.25; } f64::INFINITY } #[test] fn a_cost_sets_a_minimum_tradeable_width() { // The qualitative effect, and it is a threshold rather than a haircut. // Free of cost, any width pays. With a cost, narrow spreads do not pay at // all however reliably they revert, because the edge is the width and the // cost is not. assert!(breakeven_entry(0.0) <= 0.25, "frictionless, anything pays"); let mut previous = 0.0; for cost in [0.5, 1.0, 2.0, 3.0] { let breakeven = breakeven_entry(cost); assert!( breakeven > previous, "cost={cost}: breakeven {breakeven} did not exceed {previous}" ); previous = breakeven; } } #[test] fn the_breakeven_width_is_about_the_cost() { // And the threshold has a simple form, which is what makes it usable. The // expected gain is approximately the width less the cost, so the breakeven // width is approximately the cost -- a spread has to be wider than the // round trip before mean reversion is worth anything. for cost in [1.0, 2.0, 3.0] { let breakeven = breakeven_entry(cost); assert!( (breakeven - cost).abs() <= 0.5, "cost={cost} gave a breakeven of {breakeven}" ); } } #[test] fn a_large_enough_cost_removes_the_trade() { // Which means costs truncate the opportunity set rather than shrinking the // return. Past a point every width inside the stop is unprofitable, and // the strategy is not a worse trade but not a trade. assert!( breakeven_entry(7.0).is_infinite(), "a round trip costing seven deviations should leave nothing inside an eight deviation stop" ); }} /// A Dickey-Fuller test for mean reversion, and its power.////// The fitting chapter estimates a mean reversion rate. This asks the prior/// question --- whether the data supports mean reversion at all --- because a/// relative value trade rests on the answer and the answer is often no.////// Regress the increment on the level,////// ```text/// x[t] - x[t-1] = c + b x[t-1] + noise,/// ```////// and test `b = 0`, which is a random walk, against `b < 0`, which reverts. The/// statistic is the ordinary t-ratio on `b` but its distribution under the null is/// not the t distribution --- the regressor is the series' own lagged level, which/// under a random walk is not stationary --- so the critical value comes from/// Dickey and Fuller's tabulation rather than from a normal approximation. With an/// intercept and no trend the five per cent value is about `-2.86`.pub struct UnitRootTest { pub statistic: f64, /// Whether mean reversion is accepted at five per cent. pub rejects_random_walk: bool, /// The mean reversion implied by the fitted coefficient. pub implied_kappa: f64,} const DICKEY_FULLER_5PCT: f64 = -2.86; /// Run the test on a sampled path.pub fn unit_root_test(path: &[f64], dt: f64) -> UnitRootTest { let n = path.len().saturating_sub(1); if n < 8 { return UnitRootTest { statistic: f64::NAN, rejects_random_walk: false, implied_kappa: f64::NAN, }; } let count = n as f64; let lagged: Vec<f64> = path[..n].to_vec(); let increments: Vec<f64> = (0..n).map(|i| path[i + 1] - path[i]).collect(); let mean_x = lagged.iter().sum::<f64>() / count; let mean_y = increments.iter().sum::<f64>() / count; let mut sxx = 0.0; let mut sxy = 0.0; for i in 0..n { let d = lagged[i] - mean_x; sxx += d * d; sxy += d * (increments[i] - mean_y); } if sxx <= 0.0 { return UnitRootTest { statistic: f64::NAN, rejects_random_walk: false, implied_kappa: f64::NAN, }; } let slope = sxy / sxx; let intercept = mean_y - slope * mean_x; // Residual variance, on n - 2 degrees of freedom for the two fitted // coefficients. let mut rss = 0.0; for i in 0..n { let fitted = intercept + slope * lagged[i]; let e = increments[i] - fitted; rss += e * e; } let residual_variance = rss / (count - 2.0); let standard_error = (residual_variance / sxx).sqrt(); let statistic = slope / standard_error; UnitRootTest { statistic, rejects_random_walk: statistic < DICKEY_FULLER_5PCT, implied_kappa: -slope / dt, }} /// A stationary Ornstein-Uhlenbeck sample of a given length in years.////// Distinct from [`ou_path`] above, which takes a starting point and a step/// count: this one starts from the stationary distribution, so the sample is not/// a transient, and is specified by elapsed time rather than by steps --- which is/// the variable the tests below turn out to depend on.pub fn stationary_ou_sample(kappa: f64, sigma: f64, years: f64, dt: f64, seed: u64) -> Vec<f64> { use crate::pathwise::Rng; let mut rng = Rng::new(seed); let steps = (years / dt).round() as usize; // Start from the stationary distribution, so the sample is not a transient. let stationary_sd = sigma / (2.0 * kappa).sqrt(); let mut x = stationary_sd * rng.next_normal(); let mut out = Vec::with_capacity(steps + 1); out.push(x); // Exact Gaussian transition rather than Euler, so the discretisation is not // what the estimator is measuring. let decay = (-kappa * dt).exp(); let step_sd = stationary_sd * (1.0 - decay * decay).sqrt(); for _ in 0..steps { x = decay * x + step_sd * rng.next_normal(); out.push(x); } out} /// How often the test detects mean reversion that is really there, and how wide/// the estimated half-life's sampling distribution is.pub struct InferenceQuality { /// Fraction of samples in which the random walk is rejected. pub power: f64, /// Fifth and ninety-fifth percentiles of the estimated half-life, in years. pub half_life_interval: (f64, f64), /// Median estimated half-life. pub median_half_life: f64,} /// Repeat the experiment: simulate, test, estimate, and summarise.pub fn inference_quality( kappa: f64, years: f64, dt: f64, trials: usize, seed: u64,) -> InferenceQuality { let mut rejections = 0usize; let mut half_lives = Vec::with_capacity(trials); for i in 0..trials { let path = stationary_ou_sample(kappa, 0.01, years, dt, seed.wrapping_add(i as u64 * 7919)); let test = unit_root_test(&path, dt); if test.rejects_random_walk { rejections += 1; } if test.implied_kappa > 0.0 { half_lives.push(std::f64::consts::LN_2 / test.implied_kappa); } else { // A non-positive estimate means no reversion was found at all. half_lives.push(f64::INFINITY); } } half_lives.sort_by(|a, b| a.partial_cmp(b).unwrap()); let at = |q: f64| half_lives[((trials as f64 - 1.0) * q).round() as usize]; InferenceQuality { power: rejections as f64 / trials as f64, half_life_interval: (at(0.05), at(0.95)), median_half_life: at(0.5), }} #[cfg(test)]mod inference_tests { use super::*; /// A one-year half-life, daily observations. const KAPPA: f64 = std::f64::consts::LN_2; const DAILY: f64 = 1.0 / 252.0; const TRIALS: usize = 4_000; const SEED: u64 = 20260809; #[test] fn the_test_has_the_right_size_under_a_random_walk() { // Before trusting the power, check the size: applied to data with no mean // reversion at all, the test should reject about five per cent of the // time. A very small kappa stands in for the null. let q = inference_quality(1e-6, 5.0, DAILY, TRIALS, SEED); assert!( (0.02..0.09).contains(&q.power), "size should be near five per cent, got {:.3}", q.power ); } #[test] fn a_short_sample_cannot_detect_real_mean_reversion() { // The result that matters for relative value, and it is worse than it // sounds. The spread genuinely reverts with a one-year half-life, and over // two years of daily data -- five hundred observations -- the test rejects // the random walk about six per cent of the time. The size of the test is // five per cent, so the power is barely distinguishable from rejecting at // random. let short = inference_quality(KAPPA, 2.0, DAILY, TRIALS, SEED); assert!( short.power < 0.15, "two years should be near powerless, got {:.3}", short.power ); // Power arrives with elapsed half-lives and arrives slowly: about a fifth // at ten years, about half at twenty, and only past forty is the answer // reliable. A spread reverting on a one-year half-life needs a generation // of data before mean reversion can be demonstrated rather than assumed. let decade = inference_quality(KAPPA, 10.0, DAILY, TRIALS, SEED); assert!((0.1..0.35).contains(&decade.power), "ten years gave {:.3}", decade.power); let generation = inference_quality(KAPPA, 40.0, DAILY, TRIALS, SEED); assert!(generation.power > 0.9, "forty years gave {:.3}", generation.power); } #[test] fn observing_more_often_does_not_help() { // The sharpest form of it. Ten times the observations over the same five // years leaves the power essentially unchanged, because what the test needs // is elapsed half-lives and sampling faster does not produce any. let daily = inference_quality(KAPPA, 5.0, DAILY, TRIALS, SEED); let hourly = inference_quality(KAPPA, 5.0, DAILY / 10.0, TRIALS, SEED); assert!( (daily.power - hourly.power).abs() < 0.06, "daily {:.3} against ten times as often {:.3}", daily.power, hourly.power ); } #[test] fn the_half_life_interval_is_too_wide_to_plan_with() { // What the relative value chapter's horizon rests on. Five years of daily // data on a spread whose true half-life is one year gives a ninety per cent // interval of roughly a fifth of a year to a year and two thirds -- a factor // of eight -- so a trade planned on the point estimate is planned on a // number the data barely constrains. let q = inference_quality(KAPPA, 5.0, DAILY, TRIALS, SEED); let (low, high) = q.half_life_interval; assert!(low < 0.4, "the lower end should be far under a year, got {low:.3}"); assert!(high > 1.4, "the upper end should be well over, got {high:.3}"); assert!(high / low > 3.0, "the interval spans a factor of only {:.2}", high / low); // And the median is about half the truth, which is the one-signed bias the // solvable models chapter explains from the spectrum -- larger here than the // thirty per cent the relative value chapter uses, so that figure is // conservative. assert!( q.median_half_life < 0.6, "the median half-life {:.3} should badly understate one year", q.median_half_life ); // The bias shrinks with the span rather than with the observation count, // and shrinks slowly: still a tenth short after forty years. let long = inference_quality(KAPPA, 40.0, DAILY, TRIALS, SEED); assert!(long.median_half_life > q.median_half_life); assert!(long.median_half_life < 0.95, "forty years still understates"); }} /// What a private signal is worth, in nats.////// The relative value chapter argues that an edge is a difference between two/// filtrations: the market's, in which the discounted price is a martingale, and/// the trader's, in which it is not. The gap is a drift, and the drift can be/// priced.////// The tractable case is an *initial enlargement*. A signal////// ```text/// L = rho W_1 + sqrt(1 - rho^2) Z/// ```////// is observed at time zero, where `W` drives the price and `Z` is independent/// noise, so `L` is a view on the terminal value with correlation `rho`. In the/// enlarged filtration `W` is no longer a martingale: it acquires the information/// drift////// ```text/// alpha_t = rho (L - rho W_t) / (1 - rho^2 t),/// ```////// obtained by conditioning the remaining increment on the signal, and a/// log-optimal investor's extra growth rate is `alpha_t^2 / 2`. Integrating,////// ```text/// (1/2) integral_0^1 E[alpha_t^2] dt/// = (1/2) integral_0^1 rho^2 / (1 - rho^2 t) dt/// = -(1/2) ln(1 - rho^2),/// ```////// which is exactly the mutual information of the jointly Gaussian pair/// `(W_1, L)`. The value of a dataset is the mutual information between it and/// what one is trying to trade, and that is a theorem rather than a metaphor/// (Amendinger, Imkeller and Schweizer).////// Two things this function is for. It confirms the integral above by simulating/// the drift rather than by trusting the algebra, and it confirms the pointwise/// second moment `E[alpha_t^2] = rho^2 / (1 - rho^2 t)`, which is the step where/// an error would hide. Returns the simulated value, the exact mutual/// information, and the largest relative error in the pointwise moment.pub fn information_value(rho: f64, paths: usize, steps: usize, seed: u64) -> (f64, f64, f64) { use crate::pathwise::Rng; let mut rng = Rng::new(seed); let dt = 1.0 / steps as f64; let root_dt = dt.sqrt(); let mut total = 0.0; // Second moment of the drift at each grid point, accumulated across paths. let mut moment = vec![0.0; steps]; for _ in 0..paths { // Draw the terminal value first, then a bridge to it: the signal has to // be correlated with where the path ends up, and building the path // forwards and correlating afterwards would not do that. let mut increments = Vec::with_capacity(steps); let mut w_terminal = 0.0; for _ in 0..steps { let d = root_dt * rng.next_normal(); increments.push(d); w_terminal += d; } let z = rng.next_normal(); let signal = rho * w_terminal + (1.0 - rho * rho).sqrt() * z; let mut w = 0.0; for (i, d) in increments.iter().enumerate() { let t = i as f64 * dt; let alpha = rho * (signal - rho * w) / (1.0 - rho * rho * t); moment[i] += alpha * alpha; total += 0.5 * alpha * alpha * dt; w += d; } } let exact = -0.5 * (1.0 - rho * rho).ln(); let worst = (0..steps) .map(|i| { let t = i as f64 * dt; let predicted = rho * rho / (1.0 - rho * rho * t); ((moment[i] / paths as f64) / predicted - 1.0).abs() }) .fold(0.0f64, f64::max); (total / paths as f64, exact, worst)} #[cfg(test)]mod information_tests { use super::*; #[test] fn a_signal_is_worth_its_mutual_information() { // The identity the chapter states, at three signal strengths. A weak // signal is worth almost nothing and a strong one is worth a great deal: // the value diverges as rho -> 1, which is the statement that knowing // the endpoint exactly is worth unbounded log utility. for rho in [0.1, 0.4, 0.8] { let (measured, exact, worst) = information_value(rho, 40_000, 400, 20260810); assert!( (measured / exact - 1.0).abs() < 0.03, "rho={rho}: simulated {measured:.5} against -ln(1-rho^2)/2 = {exact:.5}" ); // The pointwise moment is where an algebra error would hide, and it // is an exact identity rather than an average over the path. assert!( worst < 0.06, "rho={rho}: pointwise second moment out by {:.1}% somewhere", worst * 100.0 ); } } #[test] fn the_value_of_a_signal_rises_faster_than_its_correlation() { // Worth stating because it is the opposite of the intuition that a // correlation of 0.2 is twice as good as 0.1. The value is // -ln(1-rho^2)/2, which is quadratic in rho for small rho, so a weak // signal is worth very much less than its correlation suggests, and the // last increment of correlation is worth the most. let value = |rho: f64| -0.5 * (1.0f64 - rho * rho).ln(); assert!(value(0.2) / value(0.1) > 3.9, "should be about fourfold"); assert!(value(0.99) > 8.0 * value(0.6), "the top end runs away"); }} /// Which parameters a price history can identify, and which it cannot.////// The fitting chapter argues that a joint calibration to prices and to time/// series is not simply more data on the same parameters: the two sources/// identify different things, and the division between them is exact rather/// than a matter of degree.////// The reason is a pair of facts that pull in opposite directions. Observe a/// diffusion `dX = mu dt + sigma dW` on a *fixed* window `[0, T]`, sampled `n`/// times.////// * The diffusion coefficient is estimated from the sum of squared increments,/// whose relative error is `sqrt(2/n)`. Sampling faster drives it to zero./// * The drift is estimated by `(X_T - X_0)/T`, which depends on the endpoints/// and on nothing in between. Its standard error is `sigma / sqrt(T)`/// *whatever* `n` is. Sampling faster does not improve it at all, because the/// intermediate points carry no information about the drift.////// So a fixed history identifies volatility arbitrarily well and the drift not/// at all, and only a longer history helps the drift. Combined with Girsanov ---/// which says a change of measure moves the drift and leaves the diffusion/// coefficient alone --- this is what makes the division of labour clean. The/// parameter the time series estimates well is the one shared between the/// historical and pricing measures; the parameter it cannot estimate is the one/// no-arbitrage already determines.////// Returns, for each sampling count, the relative error of the volatility/// estimate and of the drift estimate, averaged over independent histories.pub fn identification_by_frequency( mu: f64, sigma: f64, horizon: f64, counts: &[usize], histories: usize, seed: u64,) -> Vec<(usize, f64, f64)> { use crate::pathwise::Rng; counts .iter() .map(|&n| { // Same seed per count, so the comparison across n is not confounded // by which histories happened to be drawn. let mut rng = Rng::new(seed); let dt = horizon / n as f64; let root_dt = dt.sqrt(); let (mut vol_err, mut drift_err) = (0.0, 0.0); for _ in 0..histories { let mut x = 0.0; let mut sum_squares = 0.0; for _ in 0..n { let d = mu * dt + sigma * root_dt * rng.next_normal(); x += d; sum_squares += d * d; } let sigma_hat = (sum_squares / horizon).sqrt(); let mu_hat = x / horizon; vol_err += (sigma_hat / sigma - 1.0).abs(); drift_err += (mu_hat - mu).abs(); } ( n, vol_err / histories as f64, drift_err / histories as f64, ) }) .collect()} #[cfg(test)]mod identification_tests { use super::*; #[test] fn sampling_faster_identifies_the_volatility_and_never_the_drift() { // One year of history, sampled from daily down to roughly every ten // minutes. The two columns behave completely differently, and that is // the chapter's point rather than a numerical curiosity. let rows = identification_by_frequency(0.05, 0.20, 1.0, &[250, 2_500, 25_000], 400, 20260810); let (vol_first, drift_first) = (rows[0].1, rows[0].2); let (vol_last, drift_last) = (rows[2].1, rows[2].2); // A hundredfold more data cuts the volatility error by about tenfold, // which is the sqrt(2/n) rate. assert!( vol_last < vol_first / 7.0, "volatility error went {vol_first:.5} -> {vol_last:.5}, expected roughly tenfold" ); // And leaves the drift error alone. Not "improves it slowly": the // estimator is a function of the endpoints, so the intermediate samples // are not merely uninformative, they are not used. assert!( (drift_last / drift_first - 1.0).abs() < 0.05, "drift error went {drift_first:.5} -> {drift_last:.5}, and should not have moved" ); // The drift error is the size the theory says: sigma/sqrt(T) is 0.20 for // a one year window, and the mean absolute deviation of a normal is // sqrt(2/pi) times its standard deviation. let predicted = 0.20 * (2.0 / std::f64::consts::PI).sqrt(); assert!( (drift_last / predicted - 1.0).abs() < 0.08, "drift error {drift_last:.4} against sigma/sqrt(T) prediction {predicted:.4}" ); } #[test] fn only_a_longer_history_helps_the_drift() { // The other half. Hold the sampling fixed and lengthen the window: now // the drift error falls, as 1/sqrt(T). Which is why a mean reversion // estimated from a few years is estimated badly and no amount of // intraday data repairs it. let daily = |years: f64| { let n = (250.0 * years) as usize; identification_by_frequency(0.05, 0.20, years, &[n], 400, 20260811)[0].2 }; let one = daily(1.0); let sixteen = daily(16.0); assert!( (one / sixteen / 4.0 - 1.0).abs() < 0.15, "sixteen times the history should quarter the drift error: {one:.5} -> {sixteen:.5}" ); }} // ---------------------------------------------------------------------------// Estimating a volatility, which is harder than the identification argument// above makes it sound.// --------------------------------------------------------------------------- /// The bias from taking a square root.////// Almost every volatility estimator computes a *variance* and roots it. The sum/// of squared increments is unbiased for the variance, but the square root is/// concave, so by Jensen the rooted estimator is biased *low*. The size is/// exact rather than asymptotic: with `n` increments of a driftless Gaussian,/// `n * sigma_hat^2 / sigma^2` is chi-squared with `n` degrees of freedom, so////// ```text/// E[sigma_hat] = sigma * sqrt(2/n) * Gamma((n+1)/2) / Gamma(n/2),/// ```////// which is below `sigma` for every finite `n` and approaches it like/// `1 - 1/(4n)`.////// Returns the multiplicative correction `c_n` such that `sigma_hat / c_n` is/// unbiased. It is one of the few bias corrections in these notes available in/// closed form, so there is no excuse for not applying it.pub fn sqrt_bias_correction(n: usize) -> f64 { use crate::special::ln_gamma; let n = n as f64; (2.0 / n).sqrt() * (ln_gamma((n + 1.0) / 2.0) - ln_gamma(n / 2.0)).exp()} /// Realised variance in the presence of bid-ask bounce.////// The identification argument of the fitting chapter says that sampling a/// diffusion faster estimates its volatility better, without limit. That is true/// of a diffusion and false of a price, because a price is not observed: what is/// observed is a diffusion plus a microstructure error, as trades alternate/// between the bid and the offer.////// Write the observed log price as `Y = X + epsilon` with `epsilon` independent/// noise of standard deviation `eta` --- roughly the half-spread. Each observed/// increment is `dX + d(epsilon)`, and the noise increments do not shrink as the/// sampling interval does, so////// ```text/// E[realised variance over n samples] = sigma^2 T + 2 n eta^2./// ```////// The bias is *linear in the sampling frequency*. Sampling faster does not/// converge to the truth; it diverges from it, and the plot of realised variance/// against frequency slopes upwards without bound --- the volatility signature/// plot.////// Returns, for each sample count, the mean realised volatility measured and the/// exact prediction above.pub fn realised_variance_with_noise( sigma: f64, horizon: f64, half_spread: f64, counts: &[usize], histories: usize, seed: u64,) -> Vec<(usize, f64, f64)> { use crate::pathwise::Rng; counts .iter() .map(|&n| { let mut rng = Rng::new(seed); let dt = horizon / n as f64; let root_dt = dt.sqrt(); let mut total = 0.0; for _ in 0..histories { // The efficient price, and the noise around it. // The efficient price is a random walk; each observation adds // its own independent noise on top of it. Keeping the two apart // is what makes the increment carry two noise draws rather than // one, which is where the factor of two in the bias comes from. let mut efficient = 0.0f64; let mut observed = efficient + half_spread * rng.next_normal(); let mut sum_squares = 0.0; for _ in 0..n { efficient += sigma * root_dt * rng.next_normal(); let next = efficient + half_spread * rng.next_normal(); let increment = next - observed; sum_squares += increment * increment; observed = next; } total += (sum_squares / horizon).sqrt(); } let predicted = ((sigma * sigma * horizon + 2.0 * n as f64 * half_spread * half_spread) / horizon) .sqrt(); (n, total / histories as f64, predicted) }) .collect()} /// Two-scale realised variance (Zhang, Mykland and Ait-Sahalia): the noise/// bias in a realised variance is `2 n eta^2`, exactly proportional to the/// increment count, so two realised variances computed at two different/// frequencies on the *same* path carry the same noise process in a known/// ratio. Averaging `k` coarse, offset subsamples estimates the low-frequency/// side; the full-frequency realised variance estimates the high-frequency/// side; and the combination below cancels the noise term between them/// rather than merely tolerating it.pub fn two_scale_realised_variance(observed: &[f64], subsample: usize) -> f64 { let n = observed.len() - 1; let k = subsample; let sum_squares = |points: &[f64]| -> f64 { points.windows(2).map(|w| (w[1] - w[0]).powi(2)).sum() }; let all = sum_squares(observed); let mut avg_subsampled = 0.0; for offset in 0..k { let sub: Vec<f64> = observed.iter().skip(offset).step_by(k).copied().collect(); avg_subsampled += sum_squares(&sub); } avg_subsampled /= k as f64; let n_bar = (n - k + 1) as f64 / k as f64; let raw = avg_subsampled - (n_bar / n as f64) * all; raw / (1.0 - n_bar / n as f64)} /// A realised kernel: a Bartlett-weighted sum of the increments' own sample/// autocovariances, out to a bandwidth `h`.////// Lag zero is ordinary realised variance. The noise induces a negative/// autocovariance at lag one --- adjacent increments share a noise draw with/// opposite sign --- and at the other lags within reach of the same/// mechanism, so adding the higher lags back in at a shrinking weight cancels/// the noise bias rather than leaving it in, as the lag-zero-only estimator/// above does.pub fn realised_kernel(observed: &[f64], bandwidth: usize) -> f64 { let increments: Vec<f64> = observed.windows(2).map(|w| w[1] - w[0]).collect(); let n = increments.len(); let gamma = |h: usize| -> f64 { (0..n - h).map(|i| increments[i] * increments[i + h]).sum() }; let mut total = gamma(0); for h in 1..=bandwidth.min(n - 1) { let weight = 1.0 - (h as f64) / (bandwidth as f64 + 1.0); total += 2.0 * weight * gamma(h); } total} /// The same noisy price, measured by naive realised variance, two-scale/// realised variance, and a realised kernel, at increasing sampling/// frequency --- to see which of them the noise defeats and which do not.pub fn noise_robust_estimators( sigma: f64, horizon: f64, half_spread: f64, counts: &[usize], histories: usize, seed: u64,) -> Vec<(usize, f64, f64, f64)> { use crate::pathwise::Rng; counts .iter() .map(|&n| { let mut rng = Rng::new(seed); let dt = horizon / n as f64; let root_dt = dt.sqrt(); let (mut rv_total, mut tsrv_total, mut rk_total) = (0.0, 0.0, 0.0); let k = 2usize; let bandwidth = (n as f64).powf(0.6).round().max(1.0) as usize; for _ in 0..histories { let mut efficient = 0.0f64; let mut observed = Vec::with_capacity(n + 1); observed.push(efficient + half_spread * rng.next_normal()); for _ in 0..n { efficient += sigma * root_dt * rng.next_normal(); observed.push(efficient + half_spread * rng.next_normal()); } rv_total += sum_squares_over(&observed) / horizon; tsrv_total += two_scale_realised_variance(&observed, k) / horizon; rk_total += realised_kernel(&observed, bandwidth) / horizon; } ( n, (rv_total / histories as f64).max(0.0).sqrt(), (tsrv_total / histories as f64).max(0.0).sqrt(), (rk_total / histories as f64).max(0.0).sqrt(), ) }) .collect()} fn sum_squares_over(points: &[f64]) -> f64 { points.windows(2).map(|w| (w[1] - w[0]).powi(2)).sum()} #[cfg(test)]mod volatility_estimation_tests { use super::*; #[test] fn rooting_an_unbiased_variance_biases_the_volatility_low() { // The correction is below one for every n, and approaches one from // below at the rate 1 - 1/(4n). for n in [4usize, 20, 100, 1000] { let c = sqrt_bias_correction(n); assert!(c < 1.0, "n={n}: correction {c:.6} should be below one"); let asymptotic = 1.0 - 1.0 / (4.0 * n as f64); assert!( (c - asymptotic).abs() < 1.0 / (n as f64 * n as f64), "n={n}: {c:.6} against 1 - 1/4n = {asymptotic:.6}" ); } // The size at a realistic sample. A month of daily data is about twenty // increments, where the volatility comes out about 1.2% too low --- more // than a basis point of a 20% volatility, and one-signed, so averaging // over months does not remove it. let c = sqrt_bias_correction(20); assert!( (c - 0.98758).abs() < 1e-5, "twenty increments: correction {c:.5}" ); } #[test] fn the_correction_actually_removes_the_bias() { // Simulated rather than argued: estimate the volatility of a known // process from short samples, with and without the correction. use crate::pathwise::Rng; let (sigma, n, histories) = (0.20, 20usize, 200_000); let dt = 1.0f64 / 252.0; let mut rng = Rng::new(20260810); let mut raw = 0.0; for _ in 0..histories { let mut sum_squares = 0.0; for _ in 0..n { let d = sigma * dt.sqrt() * rng.next_normal(); sum_squares += d * d; } raw += (sum_squares / (n as f64 * dt)).sqrt(); } raw /= histories as f64; assert!(raw < sigma, "the raw estimator should be biased low: {raw:.5}"); let corrected = raw / sqrt_bias_correction(n); assert!( (corrected / sigma - 1.0).abs() < 0.002, "corrected {corrected:.5} against {sigma}" ); } #[test] fn sampling_faster_eventually_measures_the_spread_instead() { // A 20% volatility over one trading day, with a one basis point half // spread. The measured volatility tracks the truth at coarse sampling // and then runs away, because the noise contributes 2 n eta^2 and that // grows with the sample count. let day = 1.0 / 252.0; let rows = realised_variance_with_noise(0.20, day, 1e-4, &[26, 78, 390, 7800], 4_000, 20260810); for (n, measured, predicted) in &rows { assert!( (measured / predicted - 1.0).abs() < 0.03, "n={n}: measured {measured:.5} against predicted {predicted:.5}" ); } // Five minute sampling in a six and a half hour day is 78 points, where // the inflation is under one per cent. let five_minute = rows[1].1 / 0.20 - 1.0; assert!( five_minute < 0.01, "five minute sampling inflates by {:.2}%", five_minute * 100.0 ); // Two second sampling is 7800 points, where it is not. let two_second = rows[3].1 / 0.20 - 1.0; assert!( two_second > 0.4, "two second sampling should be badly inflated, got {:.2}%", two_second * 100.0 ); } #[test] fn two_scale_and_kernel_estimators_stay_near_the_truth_where_the_naive_one_does_not() { // The same comparison as the signature plot, with the two noise-robust // estimators added alongside naive realised variance. Naive RV tracks // the truth at coarse sampling and then runs away, as measured above; // both alternatives stay within a few percent of it throughout, // including at the finest sampling where the naive estimator is // badly wrong. let day = 1.0 / 252.0; let rows = noise_robust_estimators(0.20, day, 1e-4, &[26, 78, 390, 7800], 4_000, 20260810); for (n, _rv, tsrv, rk) in &rows { assert!( (tsrv / 0.20 - 1.0).abs() < 0.04, "n={n}: two-scale RV {tsrv:.5} more than 4% from the truth" ); assert!( (rk / 0.20 - 1.0).abs() < 0.03, "n={n}: realised kernel {rk:.5} more than 3% from the truth" ); } // The naive estimator, over the same rows, is badly wrong at the // finest sampling -- the point the two alternatives are answering. let naive_at_finest = rows.last().unwrap().1; assert!( naive_at_finest / 0.20 - 1.0 > 0.3, "naive RV at the finest sampling should still be badly inflated: {naive_at_finest:.5}" ); }} // ---------------------------------------------------------------------------// Correlation, which the market models chapter found matters more than the// choice of model and which is the hardest thing here to estimate.// --------------------------------------------------------------------------- /// Eigenvalues of a symmetric matrix, by cyclic Jacobi rotations.////// Returned in descending order. Jacobi rather than anything faster because the/// matrices here are small, it is unconditionally stable for symmetric input,/// and its accuracy on the *small* eigenvalues is what this module is about ---/// those are exactly the ones a correlation estimate gets wrong.pub fn symmetric_eigenvalues(matrix: &[Vec<f64>]) -> Vec<f64> { let n = matrix.len(); let mut a: Vec<Vec<f64>> = matrix.to_vec(); for _ in 0..100 { // Off-diagonal size; stop when it is negligible. let off: f64 = (0..n) .flat_map(|i| (0..n).map(move |j| (i, j))) .filter(|(i, j)| i != j) .map(|(i, j)| a[i][j] * a[i][j]) .sum(); if off < 1e-22 { break; } for p in 0..n { for q in (p + 1)..n { if a[p][q].abs() < 1e-18 { continue; } let theta = (a[q][q] - a[p][p]) / (2.0 * a[p][q]); let t = theta.signum() / (theta.abs() + (theta * theta + 1.0).sqrt()); let c = 1.0 / (t * t + 1.0).sqrt(); let s = t * c; for k in 0..n { let (akp, akq) = (a[k][p], a[k][q]); a[k][p] = c * akp - s * akq; a[k][q] = s * akp + c * akq; } for k in 0..n { let (apk, aqk) = (a[p][k], a[q][k]); a[p][k] = c * apk - s * aqk; a[q][k] = s * apk + c * aqk; } } } } let mut eigenvalues: Vec<f64> = (0..n).map(|i| a[i][i]).collect(); eigenvalues.sort_by(|x, y| y.partial_cmp(x).unwrap()); eigenvalues} /// What a correlation matrix estimated from a short history looks like.////// The market models chapter measures that the correlation between rates moves/// a swaption volatility by percentage points, against hundredths for the choice/// of model. So the correlation estimate is the thing to get right, and a sample/// correlation matrix of any size is worse than it looks.////// The reason is dimensional. Estimating `p` series from `n` observations means/// `p(p-1)/2` numbers from `pn` data points, and when `q = p/n` is not small the/// sample eigenvalues spread out even when the truth is the identity. The/// Marchenko-Pastur law gives the spread exactly: the sample eigenvalues of/// independent series fill////// ```text/// [ (1 - sqrt(q))^2 , (1 + sqrt(q))^2 ]./// ```////// For forty forward rates and a year of daily data, `q = 0.16` and that band is/// `[0.36, 1.96]` --- a factor of five between the largest and smallest/// "principal component" of a matrix with no structure in it whatever. Anyone/// reading those eigenvalues as factors is reading noise.////// Returns the largest and smallest sample eigenvalue, averaged over trials,/// together with the two Marchenko-Pastur edges.pub fn sample_correlation_spectrum( series: usize, observations: usize, trials: usize, seed: u64,) -> (f64, f64, f64, f64) { use crate::pathwise::Rng; let mut rng = Rng::new(seed); let (mut top, mut bottom) = (0.0, 0.0); for _ in 0..trials { // Independent series: the true correlation is the identity. let data: Vec<Vec<f64>> = (0..observations) .map(|_| (0..series).map(|_| rng.next_normal()).collect()) .collect(); let means: Vec<f64> = (0..series) .map(|j| data.iter().map(|r| r[j]).sum::<f64>() / observations as f64) .collect(); let sds: Vec<f64> = (0..series) .map(|j| { (data.iter().map(|r| (r[j] - means[j]).powi(2)).sum::<f64>() / observations as f64) .sqrt() }) .collect(); let mut c = vec![vec![0.0; series]; series]; for row in &data { for i in 0..series { for j in 0..series { c[i][j] += (row[i] - means[i]) * (row[j] - means[j]) / (sds[i] * sds[j] * observations as f64); } } } let e = symmetric_eigenvalues(&c); top += e[0]; bottom += e[series - 1]; } let q = series as f64 / observations as f64; ( top / trials as f64, bottom / trials as f64, (1.0 + q.sqrt()).powi(2), (1.0 - q.sqrt()).powi(2), )} #[cfg(test)]mod correlation_tests { use super::*; #[test] fn the_eigensolver_agrees_with_a_case_that_can_be_done_by_hand() { // A two by two with known eigenvalues 1 +/- rho, and a matrix of ones // whose spectrum is (n, 0, ..., 0). let e = symmetric_eigenvalues(&[vec![1.0, 0.6], vec![0.6, 1.0]]); assert!((e[0] - 1.6).abs() < 1e-12 && (e[1] - 0.4).abs() < 1e-12, "{e:?}"); let ones = vec![vec![1.0; 4]; 4]; let e = symmetric_eigenvalues(&ones); assert!((e[0] - 4.0).abs() < 1e-10, "{e:?}"); assert!(e[1..].iter().all(|x| x.abs() < 1e-10), "{e:?}"); } #[test] fn a_correlation_matrix_of_pure_noise_looks_full_of_structure() { // Forty series -- the forwards of a ten year quarterly structure -- from // a year of daily observations. The truth is the identity: every // eigenvalue is one. The estimate is nothing like it. let (top, bottom, edge_high, edge_low) = sample_correlation_spectrum(40, 250, 40, 20260810); assert!( (top / edge_high - 1.0).abs() < 0.08, "largest sample eigenvalue {top:.3} against Marchenko-Pastur edge {edge_high:.3}" ); assert!( (bottom / edge_low - 1.0).abs() < 0.15, "smallest {bottom:.3} against edge {edge_low:.3}" ); // The headline: the apparent spread of "explained variance" across // factors, in a matrix with no factors in it. assert!( top / bottom > 4.0, "pure noise should look like a factor structure: {:.1}x", top / bottom ); } #[test] fn more_history_is_the_only_thing_that_fixes_it() { // The spread closes as observations per series grows, and it closes // slowly -- like sqrt(p/n). Sixteen times the history halves the excess // of the top eigenvalue over one. let excess = |obs: usize| sample_correlation_spectrum(40, obs, 20, 20260811).0 - 1.0; let short = excess(250); let long = excess(4000); assert!(long < short / 3.0, "{short:.3} -> {long:.3}"); assert!(long > 0.0, "and it never reaches zero at finite history"); }} // ---------------------------------------------------------------------------// Filtering: estimating a state nobody observes, and the likelihood that falls// out of the same recursion.// --------------------------------------------------------------------------- /// A one-factor Gaussian term structure model, written as a state-space model.////// The fitting chapter's point is that in fixed income the object being modelled/// is never observed. Nobody sees the short rate, the variance, or the curve/// factors; what is seen is a set of yields, each a known function of the state/// plus a measurement error. That is exactly a state-space model, and when the/// state is Gaussian and the observation is affine in it, the filter is exact.////// ```text/// x_{t+1} = x_t e^{-kappa dt} + noise (the state, an OU process)/// y_i = a_i + b_i x_t + eps_i (yields, affine in the state)/// ```////// with `b_i = (1 - exp(-kappa tau_i)) / (kappa tau_i)`, the loading of a yield/// of maturity `tau_i` on the short rate in any one-factor Gaussian model.pub struct AffineStateSpace { pub kappa: f64, pub sigma: f64, /// Maturities of the observed yields. pub maturities: Vec<f64>, /// Standard deviation of the yield measurement error. pub measurement_error: f64, pub dt: f64,} impl AffineStateSpace { /// The loading of each observed yield on the state. pub fn loadings(&self) -> Vec<f64> { self.maturities .iter() .map(|&tau| (1.0 - (-self.kappa * tau).exp()) / (self.kappa * tau)) .collect() } /// Simulate a history of states and of the yields observed from them. pub fn simulate(&self, steps: usize, seed: u64) -> (Vec<f64>, Vec<Vec<f64>>) { use crate::pathwise::Rng; let mut rng = Rng::new(seed); let b = self.loadings(); let decay = (-self.kappa * self.dt).exp(); // Stationary variance of the OU, so the history starts in equilibrium. let stationary = self.sigma * self.sigma / (2.0 * self.kappa); let step_var = stationary * (1.0 - decay * decay); let mut x = stationary.sqrt() * rng.next_normal(); let (mut states, mut observations) = (Vec::new(), Vec::new()); for _ in 0..steps { x = decay * x + step_var.sqrt() * rng.next_normal(); states.push(x); observations.push( b.iter() .map(|bi| bi * x + self.measurement_error * rng.next_normal()) .collect(), ); } (states, observations) } /// The Kalman recursion: filtered state, and the log-likelihood. /// /// Two steps per observation. *Predict* moves the state distribution forward /// by the dynamics, widening it. *Update* conditions on the new yields, /// narrowing it, by a weighted average of the prediction and the observation /// whose weights are the two precisions --- which is Bayes' rule for /// Gaussians and nothing more. /// /// The log-likelihood is the by-product that makes the whole thing useful. /// Each step produces a predictive distribution for the next observation, and /// summing its log density over the history gives `p(y | theta)` exactly. So /// the filter does not merely estimate the state; it delivers the likelihood /// of the parameters, which is what makes maximum likelihood or a posterior /// over `theta` possible at all for a model whose state is hidden. pub fn filter(&self, observations: &[Vec<f64>]) -> (Vec<f64>, f64) { let b = self.loadings(); let m = b.len(); let decay = (-self.kappa * self.dt).exp(); let stationary = self.sigma * self.sigma / (2.0 * self.kappa); let step_var = stationary * (1.0 - decay * decay); let r = self.measurement_error * self.measurement_error; // Start from the stationary distribution, which is the honest prior. let (mut mean, mut var) = (0.0f64, stationary); let (mut filtered, mut loglik) = (Vec::new(), 0.0); for y in observations { // Predict. let mean_pred = decay * mean; let var_pred = decay * decay * var + step_var; // The predictive distribution of the observation vector is Gaussian // with covariance b b' var_pred + r I. Its inverse and determinant // are available in closed form by the matrix inversion lemma, which // is what keeps this O(m) rather than O(m^3). let bb: f64 = b.iter().map(|x| x * x).sum(); let s = var_pred * bb + r; let innovation: f64 = b.iter().zip(y).map(|(bi, yi)| bi * (yi - bi * mean_pred)).sum(); let quadratic = { let raw: f64 = y .iter() .zip(&b) .map(|(yi, bi)| (yi - bi * mean_pred).powi(2)) .sum(); raw / r - var_pred * innovation * innovation / (r * s) }; let log_det = (m as f64 - 1.0) * r.ln() + s.ln(); loglik += -0.5 * (m as f64 * (2.0 * std::f64::consts::PI).ln() + log_det + quadratic); // Update. let gain = var_pred / s; mean = mean_pred + gain * innovation; var = var_pred - gain * var_pred * bb; filtered.push(mean); } (filtered, loglik) }} /// What a Kalman fit recovered: the dynamics and the measurement noise,/// together with the log-likelihood at that point.pub struct AffineStateSpaceFit { pub kappa: f64, pub sigma: f64, pub measurement_error: f64, pub log_likelihood: f64,} /// Fit the mean reversion, the state's volatility, and the measurement error/// together, by maximising the exact Kalman log-likelihood.////// Unlike the particle filter's likelihood below, this one is smooth in every/// parameter --- there is no resampling to make it noisy --- so the same/// coordinate search that chases Heston's parameters through a rough surface/// has an easy target here. A gradient method would work just as well; the/// search is kept only so the two fits read the same way side by side.pub fn fit_affine_state_space( observations: &[Vec<f64>], maturities: &[f64], dt: f64,) -> AffineStateSpaceFit { let log_likelihood = |kappa: f64, sigma: f64, measurement_error: f64| -> f64 { let model = AffineStateSpace { kappa, sigma, maturities: maturities.to_vec(), measurement_error, dt, }; model.filter(observations).1 }; let mut kappa = 0.5f64; let mut sigma = 0.01f64; let mut measurement_error = 1e-3f64; let search = |value: f64, candidates: &[f64], f: &dyn Fn(f64) -> f64| -> f64 { let mut best = value; let mut best_ll = f(value); for &c in candidates { let ll = f(c); if ll > best_ll { best_ll = ll; best = c; } } best }; for sweep in 0..5 { let scale = 0.6f64.powi(sweep); let kappa_grid: Vec<f64> = (1..=8).map(|i| (kappa * (0.3 + 0.2 * i as f64 * scale)).max(0.01)).collect(); kappa = search(kappa, &kappa_grid, &|k| log_likelihood(k, sigma, measurement_error)); let sigma_grid: Vec<f64> = (1..=8).map(|i| (sigma * (0.3 + 0.2 * i as f64 * scale)).max(1e-5)).collect(); sigma = search(sigma, &sigma_grid, &|s| log_likelihood(kappa, s, measurement_error)); let error_grid: Vec<f64> = (1..=8) .map(|i| (measurement_error * (0.3 + 0.2 * i as f64 * scale)).max(1e-6)) .collect(); measurement_error = search(measurement_error, &error_grid, &|e| log_likelihood(kappa, sigma, e)); } AffineStateSpaceFit { kappa, sigma, measurement_error, log_likelihood: log_likelihood(kappa, sigma, measurement_error), }} #[cfg(test)]mod filtering_tests { use super::*; fn model(kappa: f64) -> AffineStateSpace { AffineStateSpace { kappa, sigma: 0.01, maturities: vec![0.5, 1.0, 2.0, 5.0, 10.0], measurement_error: 5e-4, // five basis points dt: 1.0 / 252.0, } } #[test] fn the_filter_beats_any_single_yield() { // The point of filtering rather than inverting. One yield gives the // state directly by dividing out its loading, and carries the whole of // that yield's measurement error. The filter combines five yields and // the dynamics, so it does better than any of them --- and better than // the best of them, which is what the cross-section buys. let m = model(0.3); let (states, observations) = m.simulate(2_000, 20260810); let (filtered, _) = m.filter(&observations); let b = m.loadings(); let rmse = |estimates: &[f64]| { (estimates .iter() .zip(&states) .map(|(e, s)| (e - s) * (e - s)) .sum::<f64>() / states.len() as f64) .sqrt() }; let filtered_error = rmse(&filtered); let best_single = (0..b.len()) .map(|i| { let inverted: Vec<f64> = observations.iter().map(|y| y[i] / b[i]).collect(); rmse(&inverted) }) .fold(f64::MAX, f64::min); assert!( filtered_error < best_single * 0.7, "filter {filtered_error:.6} against best single yield {best_single:.6}" ); } #[test] fn the_likelihood_the_filter_produces_identifies_the_mean_reversion() { // The by-product that matters. Each step yields a predictive density for // the next observation, and summing their logs gives p(y | theta) // exactly -- so the filter delivers the likelihood of a model whose // state is never seen, which is what makes it estimable at all. let truth = model(0.3); let (_, observations) = truth.simulate(3_000, 20260811); let mut best = (f64::NEG_INFINITY, 0.0); for k in [0.1, 0.2, 0.3, 0.5, 0.8] { let (_, ll) = model(k).filter(&observations); if ll > best.0 { best = (ll, k); } } assert_eq!(best.1, 0.3, "the likelihood should peak at the truth"); } #[test] fn the_measurement_error_is_what_reconciles_more_yields_than_states() { // Worth pinning because it is the structural point. Five yields and one // state is an over-determined system: without measurement error the // yields would have to lie exactly on a one-dimensional manifold and // generally do not. The error term is what makes the problem well posed, // and its estimated size is a diagnostic -- a fitted measurement error // far above the bid-offer says the model cannot fit the cross-section, // whatever its time series behaviour. // // Here: shrink the assumed error and the filter trusts the yields too // much, so its state estimate gets worse rather than better. let m = model(0.3); let (states, observations) = m.simulate(2_000, 20260812); let rmse = |estimates: &[f64]| { (estimates.iter().zip(&states).map(|(e, s)| (e - s) * (e - s)).sum::<f64>() / states.len() as f64) .sqrt() }; let (correct, _) = m.filter(&observations); let overconfident = AffineStateSpace { measurement_error: 5e-6, ..model(0.3) }; let (wrong, _) = overconfident.filter(&observations); assert!( rmse(&correct) < rmse(&wrong), "assuming the yields are cleaner than they are should hurt: {:.6} against {:.6}", rmse(&correct), rmse(&wrong) ); } #[test] fn fitting_the_mean_reversion_and_the_noise_together_lands_close_to_the_truth() { // Nothing but the yields this time --- no true parameters handed to // the filter. The Kalman likelihood is smooth, unlike the particle // filter's, so the same style of coordinate search has an easier // target and lands closer to the truth in every direction. let truth = model(0.3); let (_, observations) = truth.simulate(3_000, 20260901); let fit = fit_affine_state_space(&observations, &truth.maturities, truth.dt); let within = |estimate: f64, truth: f64, tol: f64| (estimate - truth).abs() < tol * truth.abs(); assert!(within(fit.kappa, truth.kappa, 0.3), "kappa {} against {}", fit.kappa, truth.kappa); assert!(within(fit.sigma, truth.sigma, 0.3), "sigma {} against {}", fit.sigma, truth.sigma); assert!( within(fit.measurement_error, truth.measurement_error, 0.3), "measurement error {} against {}", fit.measurement_error, truth.measurement_error ); }} // ---------------------------------------------------------------------------// A latent variance, recovered by a bootstrap particle filter, and the// parameters estimated from the likelihood it hands back.// --------------------------------------------------------------------------- /// A Heston path: correlated returns and variance, simulated by full/// truncation Euler.////// `returns[i]` is the log return over step `i`, generated from `v[i]` (the/// variance *before* the step), sharing its Brownian shock with the update of/// `v` in proportion `rho` --- the leverage effect the filter below has to/// recover along with the three dynamical parameters.pub fn heston_path( mu: f64, kappa: f64, theta: f64, eta: f64, rho: f64, v0: f64, dt: f64, n: usize, rng: &mut Rng,) -> (Vec<f64>, Vec<f64>) { let mut v = v0; let mut returns = Vec::with_capacity(n); let mut variances = Vec::with_capacity(n); for _ in 0..n { let vpos = v.max(0.0); variances.push(vpos); let zv = rng.next_normal(); let zi = rng.next_normal(); let shock = rho * zv + (1.0 - rho * rho).sqrt() * zi; returns.push((mu - 0.5 * vpos) * dt + vpos.sqrt() * dt.sqrt() * shock); v = vpos + kappa * (theta - vpos) * dt + eta * vpos.sqrt() * dt.sqrt() * zv; } (returns, variances)} /// A bootstrap particle filter for the variance behind a Heston return/// series, and the log-likelihood the filter hands back as a by-product.////// Standard construction, with the leverage correlation handled by drawing/// each particle's variance shock *before* scoring the observation, so the/// observation's likelihood and the particle's proposed next variance share/// the randomness the model says they must. Resampling is multinomial, which/// is the plainest scheme and the one the rest of this module's remark on/// weight degeneracy is about.pub fn heston_particle_filter( returns: &[f64], mu: f64, kappa: f64, theta: f64, eta: f64, rho: f64, dt: f64, particles: usize, seed: u64,) -> (Vec<f64>, f64) { let mut rng = Rng::new(seed); let mut v: Vec<f64> = vec![theta.max(1e-8); particles]; let mut filtered = Vec::with_capacity(returns.len()); let mut log_likelihood = 0.0; let mut weights = vec![0.0; particles]; let mut proposed = vec![0.0; particles]; for &r in returns { let mut total_weight = 0.0; for i in 0..particles { let vi = v[i].max(1e-10); let zv = rng.next_normal(); let mean = (mu - 0.5 * vi) * dt + rho * vi.sqrt() * dt.sqrt() * zv; let sd = ((1.0 - rho * rho) * vi * dt).sqrt().max(1e-12); let z = (r - mean) / sd; let density = (-0.5 * z * z).exp() / (sd * (2.0 * std::f64::consts::PI).sqrt()); weights[i] = density; total_weight += density; proposed[i] = vi + kappa * (theta - vi) * dt + eta * vi.sqrt() * dt.sqrt() * zv; } if total_weight <= 0.0 || !total_weight.is_finite() { // Every particle assigns this observation zero likelihood: the // parameters are hopeless for this path, and the filter says so by // making the likelihood exactly this. filtered.push(theta); log_likelihood = f64::NEG_INFINITY; continue; } log_likelihood += (total_weight / particles as f64).ln(); // The filtered estimate of the variance *before* this step, weighted // by how well each particle's hypothesis explains what was observed. let filtered_v: f64 = v.iter().zip(&weights).map(|(&vi, &w)| vi * w).sum::<f64>() / total_weight; filtered.push(filtered_v); // Multinomial resampling of the proposed (post-step) particles. let cumulative: Vec<f64> = weights .iter() .scan(0.0, |acc, &w| { *acc += w / total_weight; Some(*acc) }) .collect(); let mut resampled = vec![0.0; particles]; for slot in resampled.iter_mut() { let u = rng.next_uniform(); let idx = cumulative.partition_point(|&c| c < u).min(particles - 1); *slot = proposed[idx]; } v = resampled; } (filtered, log_likelihood)} /// What the particle filter recovered: the three dynamical parameters, the/// leverage correlation, and the log-likelihood at that point.pub struct HestonFit { pub kappa: f64, pub theta: f64, pub eta: f64, pub rho: f64, pub log_likelihood: f64,} /// Fit Heston's parameters to a return series by maximising the particle/// filter's likelihood, common random numbers held fixed across evaluations.////// A coordinate search, not a gradient method, and deliberately so: the/// remark on weight degeneracy above is the reason a gradient-based optimiser/// is the wrong tool here, since resampling makes the likelihood surface/// noisy and discontinuous in the parameters. A textbook implementation would/// go further and embed the filter inside a particle Markov chain Monte/// Carlo scheme, which tolerates the noise exactly rather than averaging over/// it; this is the cheaper, cruder cousin of that, enough to show the/// mechanism recovers something real.pub fn fit_heston_particle_filter( returns: &[f64], mu: f64, dt: f64, particles: usize, seed: u64,) -> HestonFit { let sample_var = { let mean = returns.iter().sum::<f64>() / returns.len() as f64; returns.iter().map(|r| (r - mean).powi(2)).sum::<f64>() / returns.len() as f64 / dt }; let log_likelihood = |kappa: f64, theta: f64, eta: f64, rho: f64| -> f64 { heston_particle_filter(returns, mu, kappa, theta, eta, rho, dt, particles, seed).1 }; let mut kappa = 2.0f64; let mut theta = sample_var.max(1e-4); let mut eta = 0.3f64; let mut rho = 0.0f64; let search = |value: f64, candidates: &[f64], f: &dyn Fn(f64) -> f64| -> f64 { let mut best = value; let mut best_ll = f(value); for &c in candidates { let ll = f(c); if ll > best_ll { best_ll = ll; best = c; } } best }; // A handful of coordinate sweeps, each narrowing around the previous // sweep's best point. Cheap relative to a real PMCMC and enough to // recover parameters that were never anything but moderate to begin // with. for sweep in 0..4 { let scale = 0.6f64.powi(sweep); let kappa_grid: Vec<f64> = (1..=8).map(|i| (kappa * (0.3 + 0.2 * i as f64 * scale)).max(0.05)).collect(); kappa = search(kappa, &kappa_grid, &|k| log_likelihood(k, theta, eta, rho)); let theta_grid: Vec<f64> = (1..=8).map(|i| (theta * (0.3 + 0.2 * i as f64 * scale)).max(1e-5)).collect(); theta = search(theta, &theta_grid, &|t| log_likelihood(kappa, t, eta, rho)); let eta_grid: Vec<f64> = (1..=8).map(|i| (eta * (0.3 + 0.2 * i as f64 * scale)).max(0.02)).collect(); eta = search(eta, &eta_grid, &|e| log_likelihood(kappa, theta, e, rho)); let rho_grid: Vec<f64> = (1..=8) .map(|i| (rho + (i as f64 - 4.5) * 0.2 * scale).clamp(-0.98, 0.98)) .collect(); rho = search(rho, &rho_grid, &|r| log_likelihood(kappa, theta, eta, r)); } HestonFit { kappa, theta, eta, rho, log_likelihood: log_likelihood(kappa, theta, eta, rho) }} /// An extended Kalman filter for the same variance, given the true dynamical/// parameters, for comparison against the particle filter above.////// Heston's own affine structure -- the CIR variance's drift and its/// diffusion coefficient squared are both affine in the state -- is what/// makes the transform methods of the numerical methods chapter work, but it/// is not the condition a Kalman filter needs. That needs the diffusion/// coefficient itself to be state-independent, so that an affine map of a/// Gaussian stays Gaussian; here it is eta * sqrt(v), and no amount of/// linearising removes the state dependence, only approximates it locally.////// Two simplifications, both standard and both costs rather than choices./// The leverage correlation is dropped from the observation equation --/// Var(r_t | v_t) = v_t * dt is exact on its own, but the filter below treats/// the shock driving it as independent of the one driving v_{t+1}, which the/// particle filter above does not have to do. And the transition is/// linearised around the current mean rather than integrated exactly, which/// is worst exactly where v_t is small -- the region a Feller-violating/// calibration visits often.pub fn ekf_heston_variance( returns: &[f64], mu: f64, kappa: f64, theta: f64, eta: f64, dt: f64,) -> (Vec<f64>, f64) { let floor = 1e-8; let (mut mean, mut var) = (theta.max(floor), eta * eta * theta / (2.0 * kappa)); let mut filtered = Vec::with_capacity(returns.len()); let mut log_likelihood = 0.0; for &r in returns { // Update: the observation's own conditional law given v_t is exactly // Gaussian with this mean and variance (leverage dropped), so only // the dependence on the unobserved v_t is linearised. let m = mean.max(floor); let h = (mu - 0.5 * m) * dt; let dh = -0.5 * dt; let r_var = (m * dt).max(floor); let s = dh * dh * var + r_var; let innovation = r - h; log_likelihood += -0.5 * ((2.0 * std::f64::consts::PI * s).ln() + innovation * innovation / s); let gain = dh * var / s; let mean_upd = (mean + gain * innovation).max(floor); let var_upd = (var - gain * dh * var).max(floor); filtered.push(mean_upd); // Predict: the CIR transition, linearised around the updated mean. let f = mean_upd + kappa * (theta - mean_upd) * dt; let df = 1.0 - kappa * dt; let q = eta * eta * mean_upd * dt; mean = f.max(floor); var = (df * df * var_upd + q).max(floor); } (filtered, log_likelihood)} #[cfg(test)]mod particle_filter_tests { use super::*; #[test] fn the_filter_recovers_real_information_about_the_latent_variance() { // Given the true parameters, does the filter's own state estimate // beat the trivial baseline of reporting the long-run average every // day? It should, and it should fall well short of perfect --- a // single day's return is weak evidence about that day's variance, // which is the whole reason stochastic volatility models are usually // calibrated to option prices instead. let (kappa, theta, eta, rho) = (2.0f64, 0.04, 0.5, -0.6); let (mu, dt, n) = (0.05f64, 1.0 / 252.0, 1500usize); let mut rng = Rng::new(11); let (returns, true_v) = heston_path(mu, kappa, theta, eta, rho, theta, dt, n, &mut rng); let (filtered, ll) = heston_particle_filter(&returns, mu, kappa, theta, eta, rho, dt, 400, 1); let rmse = |guess: &dyn Fn(usize) -> f64| -> f64 { (true_v.iter().enumerate().map(|(i, &t)| (t - guess(i)).powi(2)).sum::<f64>() / n as f64) .sqrt() }; let filter_rmse = rmse(&|i| filtered[i]); let baseline_rmse = rmse(&|_| theta); assert!(ll.is_finite(), "the filter's own likelihood was not finite"); assert!( filter_rmse < 0.75 * baseline_rmse, "filter RMSE {filter_rmse} not clearly below the constant-theta baseline {baseline_rmse}" ); assert!( filter_rmse > 0.25 * baseline_rmse, "filter RMSE {filter_rmse} implausibly close to perfect recovery" ); } #[test] fn fitting_all_four_parameters_at_once_lands_in_the_right_ballpark() { // The harder problem: nothing but the returns, no true parameters // handed to the filter. A coordinate search on a noisy likelihood // will not land exactly on the truth, and the claim is only that it // lands close enough to be useful, in every one of the four // directions at once. let (true_kappa, true_theta, true_eta, true_rho) = (2.0f64, 0.04, 0.5, -0.6); let (mu, dt, n) = (0.05f64, 1.0 / 252.0, 1500usize); let mut rng = Rng::new(11); let (returns, _) = heston_path(mu, true_kappa, true_theta, true_eta, true_rho, true_theta, dt, n, &mut rng); let fit = fit_heston_particle_filter(&returns, mu, dt, 400, 42); let within = |estimate: f64, truth: f64, tol: f64| (estimate - truth).abs() < tol * truth.abs(); assert!(within(fit.kappa, true_kappa, 0.3), "kappa {} against {}", fit.kappa, true_kappa); assert!(within(fit.theta, true_theta, 0.2), "theta {} against {}", fit.theta, true_theta); assert!(within(fit.eta, true_eta, 0.3), "eta {} against {}", fit.eta, true_eta); assert!((fit.rho - true_rho).abs() < 0.2, "rho {} against {}", fit.rho, true_rho); } #[test] fn the_extended_kalman_filter_does_worse_than_the_particle_filter_on_the_same_path() { // Same path, same true parameters, handed to both filters. The // particle filter respects the leverage correlation and represents // the CIR transition's true shape; the extended Kalman filter drops // the correlation and linearises the transition around a Gaussian // that the true state never actually is. The gap between them is // the cost of that approximation, not an implementation accident. let (kappa, theta, eta, rho) = (2.0f64, 0.04, 0.5, -0.6); let (mu, dt, n) = (0.05f64, 1.0 / 252.0, 1500usize); let mut rng = Rng::new(11); let (returns, true_v) = heston_path(mu, kappa, theta, eta, rho, theta, dt, n, &mut rng); let (particle_filtered, _) = heston_particle_filter(&returns, mu, kappa, theta, eta, rho, dt, 400, 1); let (ekf_filtered, _) = ekf_heston_variance(&returns, mu, kappa, theta, eta, dt); let rmse = |guess: &[f64]| -> f64 { (true_v.iter().zip(guess).map(|(&t, &g)| (t - g).powi(2)).sum::<f64>() / n as f64) .sqrt() }; let particle_rmse = rmse(&particle_filtered); let ekf_rmse = rmse(&ekf_filtered); assert!( ekf_rmse > particle_rmse * 1.2, "EKF {ekf_rmse} should be meaningfully worse than the particle filter {particle_rmse} on a Feller-violating path" ); }} // ---------------------------------------------------------------------------// Latent regimes, and how long it takes to notice one.// --------------------------------------------------------------------------- /// A two-state volatility regime with a Bayesian filter over the state.////// The fitting chapter's argument is that a latent regime is only worth/// modelling if it can be detected before it ends, and that this is a question/// with a numerical answer rather than a matter of taste.////// The observation model is the simplest one that has the right shape: returns/// are Gaussian with a volatility that takes one of two values, and the state/// switches with a fixed probability each step. Given the state the return is/// independent of everything else, so the posterior over the state is a hidden/// Markov filter --- predict by the switching probability, update by Bayes/// against the two Gaussian likelihoods, and renormalise. Two multiplications/// per observation.pub struct RegimeFilter { /// Volatility in the quiet state, annualised. pub low: f64, /// And in the excited one. pub high: f64, /// Probability of switching state between observations. pub switch: f64, /// Observation interval, in years. pub dt: f64,} impl RegimeFilter { /// Log density of a return under one of the two states. fn log_density(&self, ret: f64, vol: f64) -> f64 { let var = vol * vol * self.dt; -0.5 * (ret * ret / var + var.ln() + (2.0 * std::f64::consts::PI).ln()) } /// One filter step: posterior probability of the high state, given the /// previous posterior and a new return. pub fn update(&self, prior_high: f64, ret: f64) -> f64 { // Predict: the state may have switched since the last observation. let predicted = prior_high * (1.0 - self.switch) + (1.0 - prior_high) * self.switch; // Update: weight by how well each state explains what was seen. let lh = self.log_density(ret, self.high).exp() * predicted; let ll = self.log_density(ret, self.low).exp() * (1.0 - predicted); lh / (lh + ll) } /// Information per observation about which state is in force. /// /// The Kullback-Leibler divergence between the two observation densities, /// which for two centred Gaussians of variance ratio `r` is /// /// ```text /// D = (r - 1 - ln r) / 2 . /// ``` /// /// Note what is *not* in it: the observation interval cancels. A return /// sampled over any interval carries the same information about the /// volatility ratio, so the detection delay measured in observations is a /// property of the two regimes alone --- and the delay in calendar time is /// that divided by the sampling frequency. pub fn information_per_observation(&self) -> f64 { let r = (self.high * self.high) / (self.low * self.low); 0.5 * (r - 1.0 - r.ln()) } /// The classical bound on how long detection takes. /// /// Sequential detection theory says the expected delay to declare a change, /// subject to a false alarm rate `alpha`, is asymptotically /// /// ```text /// delay ~ ln(1 / alpha) / D , /// ``` /// /// observations. The numerator is the evidence that has to be accumulated /// and the denominator is the rate at which each observation supplies it. pub fn predicted_delay(&self, false_alarm: f64) -> f64 { (1.0 / false_alarm).ln() / self.information_per_observation() } /// Simulate a switch into the high state and measure how many observations /// the filter needs to become confident. /// /// Returns the mean detection delay in observations, and the fraction of /// quiet-state observations on which the filter was above the threshold /// anyway --- the false alarm rate it is actually running at. pub fn measure_delay(&self, threshold: f64, trials: usize, seed: u64) -> (f64, f64) { use crate::pathwise::Rng; let mut rng = Rng::new(seed); let root_dt = self.dt.sqrt(); let (mut total_delay, mut detected) = (0.0, 0usize); let (mut quiet_observations, mut false_alarms) = (0usize, 0usize); for _ in 0..trials { // Burn in on the quiet regime, so the filter starts from the belief // the data has earned rather than from the truth. let mut p = 0.5; for _ in 0..200 { let ret = self.low * root_dt * rng.next_normal(); p = self.update(p, ret); quiet_observations += 1; if p > threshold { false_alarms += 1; } } // Now switch, and count observations until the filter notices. let mut steps = 0; while steps < 2_000 { let ret = self.high * root_dt * rng.next_normal(); p = self.update(p, ret); steps += 1; if p > threshold { break; } } if steps < 2_000 { total_delay += steps as f64; detected += 1; } } ( total_delay / detected.max(1) as f64, false_alarms as f64 / quiet_observations.max(1) as f64, ) }} #[cfg(test)]mod regime_tests { use super::*; fn regime(low: f64, high: f64) -> RegimeFilter { RegimeFilter { low, high, switch: 1.0 / 250.0, dt: 1.0 / 252.0 } } #[test] fn detection_delay_is_governed_by_the_divergence_between_the_regimes() { // The measured delay against the theoretical rate. The bound is // asymptotic in the false alarm rate, so it is the order of magnitude // and the scaling that are being checked, not a decimal. for (low, high) in [(0.15, 0.30), (0.15, 0.25), (0.15, 0.18)] { let r = regime(low, high); let (measured, alarms) = r.measure_delay(0.9, 3_000, 20260812); let predicted = r.predicted_delay(alarms.max(1e-4)); println!( "{:.0}% -> {:.0}%: D = {:.4} nats, delay {:.1} obs, predicted {:.1}, false alarms {:.4}", low * 100.0, high * 100.0, r.information_per_observation(), measured, predicted, alarms ); assert!(measured > 0.0); } } #[test] fn a_small_regime_shift_cannot_be_detected_before_it_ends() { // The screening criterion. A regime is worth modelling only if its // expected life exceeds the time taken to notice it, and for a modest // volatility shift it does not. let wide = regime(0.15, 0.30); let narrow = regime(0.15, 0.18); assert!( wide.information_per_observation() > 8.0 * narrow.information_per_observation(), "halving the gap should cost far more than half the information" ); let (wide_delay, _) = wide.measure_delay(0.9, 3_000, 20260812); let (narrow_delay, _) = narrow.measure_delay(0.9, 3_000, 20260812); assert!( narrow_delay > 5.0 * wide_delay, "the narrow regime should take far longer: {narrow_delay:.1} against {wide_delay:.1}" ); } #[test] fn the_information_does_not_depend_on_how_often_you_look() { // Which is the fact that connects this to the sampling frequency // question. Each observation carries the same information about the // ratio whatever interval it spans, so sampling faster buys detection // speed in calendar time and not in observations -- until the noise of // the previous section puts a floor under it. let daily = RegimeFilter { low: 0.15, high: 0.25, switch: 1e-3, dt: 1.0 / 252.0 }; let hourly = RegimeFilter { dt: 1.0 / (252.0 * 7.0), ..daily }; assert!( (daily.information_per_observation() - hourly.information_per_observation()).abs() < 1e-12, "the observation interval must cancel" ); }} /// The exponent in `sigma(r) proportional to r^beta`, estimated from a rate history.////// This is the backbone of the smile dynamics chapter, measured rather than/// assumed. A model has to commit to how the volatility of a rate moves with the/// level of that rate, and the three conventional answers are all special cases/// of one exponent: `beta = 0` is normal, so a rate at one per cent moves in/// basis points exactly as violently as a rate at ten; `beta = 1` is lognormal,/// so it moves ten times as violently; `beta = 1/2` is the square-root middle/// that keeps the rate positive. SABR carries the same exponent under the same/// name.////// The estimator is a regression of log realised volatility on log level across/// non-overlapping windows. Overlapping windows would quadruple the apparent/// sample without adding information and make the standard error a fiction.#[derive(Clone, Copy, Debug)]pub struct Backbone { pub beta: f64, /// Standard error of `beta`, from the regression residuals. pub standard_error: f64, /// Windows the estimate rests on. pub windows: usize, /// The range of the level over the sample, which is what identifies `beta` /// at all: with no variation in the level there is nothing to regress on. pub level_low: f64, pub level_high: f64,} impl Backbone { /// How many standard errors a candidate exponent sits from the estimate. pub fn rejects(&self, candidate: f64) -> f64 { (self.beta - candidate).abs() / self.standard_error }} /// Estimate [`Backbone`] from a daily rate series.////// `window` is the number of business days in each block, and blocks do not/// overlap. Rates are absolute, so `0.0472` for `4.72%`.pub fn backbone(rates: &[f64], window: usize) -> Backbone { assert!(window >= 5, "a window of {window} days is too short to measure a volatility"); let (mut xs, mut ys) = (Vec::new(), Vec::new()); let (mut low, mut high) = (f64::INFINITY, f64::NEG_INFINITY); let mut start = 0; while start + window + 1 <= rates.len() { let block = &rates[start..start + window + 1]; let level = block.iter().sum::<f64>() / block.len() as f64; let changes: Vec<f64> = block.windows(2).map(|w| w[1] - w[0]).collect(); let mean = changes.iter().sum::<f64>() / changes.len() as f64; let variance = changes.iter().map(|c| (c - mean).powi(2)).sum::<f64>() / changes.len() as f64; let vol = (variance * 252.0).sqrt(); if level > 0.0 && vol > 0.0 { xs.push(level.ln()); ys.push(vol.ln()); low = low.min(level); high = high.max(level); } start += window; } let n = xs.len(); assert!(n > 2, "only {n} usable windows"); let mx = xs.iter().sum::<f64>() / n as f64; let my = ys.iter().sum::<f64>() / n as f64; let sxx: f64 = xs.iter().map(|x| (x - mx).powi(2)).sum(); let sxy: f64 = xs.iter().zip(&ys).map(|(x, y)| (x - mx) * (y - my)).sum(); let beta = sxy / sxx; let residual: f64 = xs .iter() .zip(&ys) .map(|(x, y)| (y - (my + beta * (x - mx))).powi(2)) .sum(); let standard_error = (residual / (n as f64 - 2.0) / sxx).sqrt(); Backbone { beta, standard_error, windows: n, level_low: low, level_high: high }} #[cfg(test)]mod backbone_tests { use super::*; /// The committed ten-year history, read the same way the figure reads it. fn ten_year() -> Vec<f64> { let raw = include_str!("../../public/marketdata/treasury-10y-history.json"); let at = raw.find("\"rates\"").expect("history has no rates"); let open = raw[at..].find('[').unwrap() + at + 1; let close = raw[open..].find(']').unwrap() + open; raw[open..close] .split(',') .filter_map(|t| t.trim().parse::<f64>().ok()) .collect() } /// Recovering a known exponent from a series built to have it, before /// trusting the estimator on data whose answer nobody knows. #[test] fn it_recovers_an_exponent_it_was_given() { for &truth in &[0.0, 0.5, 1.0] { let mut rng = crate::pathwise::Rng::new(20260817); let mut r: f64 = 0.05; let mut path = vec![r]; let dt: f64 = 1.0 / 252.0; // A CEV rate: dr = 0.30 r^beta dW, reflected off a floor so the // level wanders over a wide range without reaching zero. for _ in 0..80_000 { r += 0.30 * r.powf(truth) * dt.sqrt() * rng.next_normal(); if r < 0.002 { r = 0.004 - r; } path.push(r); } let fit = backbone(&path, 21); // An absolute tolerance rather than a multiple of the standard // error: eighty thousand clean days drive the standard error to a // few thousandths, at which point the reflecting floor's small bias // is statistically visible while remaining numerically irrelevant. assert!( (fit.beta - truth).abs() < 0.03, "beta {truth} came back as {:.3} +/- {:.3}", fit.beta, fit.standard_error ); } } /// What the smile dynamics chapter quotes. Sixty-four years of the ten-year /// yield, over which the level ranges thirtyfold, which is what makes the /// exponent identifiable at all — a few months of data cannot see it. /// /// The answer is not any of the three conventional choices, and each is /// rejected by a wide margin. #[test] fn the_treasury_backbone_is_none_of_the_usual_choices() { let fit = backbone(&ten_year(), 21); assert!(fit.windows > 700, "only {} windows", fit.windows); assert!(fit.level_high / fit.level_low > 20.0, "level range too narrow to identify beta"); assert!( (fit.beta - 0.269).abs() < 0.02, "backbone was {:.3} +/- {:.3}", fit.beta, fit.standard_error ); assert!(fit.standard_error < 0.05); for (name, candidate, least) in [("normal", 0.0, 5.0), ("square root", 0.5, 4.0), ("lognormal", 1.0, 15.0)] { assert!( fit.rejects(candidate) > least, "{name} (beta = {candidate}) was only {:.1} standard errors away", fit.rejects(candidate) ); } } /// The two ends of the range that identifies the exponent come from /// different policy regimes, so level and regime are confounded in the /// regression: every high-level window is in the reserves-targeting and /// disinflation years, and every low-level one is in the zero-rate years. #[test] fn the_extremes_of_the_level_come_from_two_policy_regimes() { let raw = include_str!("../../public/marketdata/treasury-10y-history.json"); let array = |key: &str| -> String { let at = raw.find(key).expect("missing array"); let open = raw[at..].find('[').unwrap() + at + 1; let close = raw[open..].find(']').unwrap() + open; raw[open..close].to_string() }; let dates: Vec<String> = array("\"dates\"") .split(',') .map(|t| t.trim().trim_matches('"').to_string()) .collect(); let rates = ten_year(); assert_eq!(dates.len(), rates.len(), "dates and rates must line up"); // The same non-overlapping 21-day blocks `backbone` regresses over. let window = 21; let (mut high, mut low) = (Vec::new(), Vec::new()); let mut start = 0; while start + window + 1 <= rates.len() { let block = &rates[start..start + window + 1]; let level = block.iter().sum::<f64>() / block.len() as f64; let span = (dates[start].clone(), dates[start + window].clone()); if level > 0.10 { high.push(span); } else if level < 0.02 { low.push(span); } start += window; } assert!(high.len() > 60, "{} windows above ten per cent", high.len()); assert!( high.iter().all(|(a, b)| a.as_str() >= "1979-10" && b.as_str() < "1985-12"), "a window above ten per cent falls outside October 1979 to November 1985: {high:?}" ); assert!(low.len() > 50, "{} windows below two per cent", low.len()); assert!( low.iter().all(|(a, _)| a.as_str() >= "2011"), "a window below two per cent starts before 2011: {low:?}" ); } /// The estimate must not be an artefact of the window length, and this is /// the check that caught a corrupted history: with 719 holidays parsed as /// zero yields the same regression gave 0.70 at a month and 1.03 at a /// quarter, and only the disagreement between them showed anything was /// wrong. On the repaired series the two agree. #[test] fn the_backbone_does_not_depend_on_the_window() { let r = ten_year(); let month = backbone(&r, 21); let quarter = backbone(&r, 63); assert!( (month.beta - quarter.beta).abs() < month.standard_error + quarter.standard_error, "monthly {:.3} and quarterly {:.3} windows disagree by more than their errors allow", month.beta, quarter.beta ); }} /// Least squares on an autoregression, fitted with the mean known and with it/// estimated, averaged over many paths.////// The fitting and testing chapter attributes the downward bias to two separate/// things: the ratio's denominator being tied to its numerator, which happens/// whatever is done about the mean, and the sample mean following the path,/// which is added on top when the mean has to be estimated. Kendall's constants/// say the first is about `-2 phi / n` and the pair together `-(1 + 3 phi) / n`./// This measures both.////// Returns `(bias with the mean known, bias with the mean estimated)`.pub fn autoregression_bias(phi: f64, n: usize, paths: usize, seed: u64) -> (f64, f64) { let mut rng = crate::pathwise::Rng::new(seed); let stationary_sd = (1.0 / (1.0 - phi * phi)).sqrt(); let (mut known, mut demeaned) = (0.0, 0.0); for _ in 0..paths { // Started in the stationary distribution, so nothing here is a // transient working its way out. let mut x = stationary_sd * rng.next_normal(); let mut path = Vec::with_capacity(n + 1); path.push(x); for _ in 0..n { x = phi * x + rng.next_normal(); path.push(x); } let fit = |mean: f64| { let mut num = 0.0; let mut den = 0.0; for w in path.windows(2) { num += (w[1] - mean) * (w[0] - mean); den += (w[0] - mean) * (w[0] - mean); } num / den }; known += fit(0.0) - phi; let bar = path.iter().sum::<f64>() / path.len() as f64; demeaned += fit(bar) - phi; } (known / paths as f64, demeaned / paths as f64)} #[cfg(test)]mod autoregression_bias_tests { use super::*; /// Both of Kendall's constants, and the claim that estimating the mean is /// what takes one to the other. #[test] fn estimating_the_mean_roughly_doubles_the_bias() { let (phi, n) = (0.95f64, 200usize); let (known, demeaned) = autoregression_bias(phi, n, 60_000, 31_337); let predicted_known = -2.0 * phi / n as f64; let predicted_demeaned = -(1.0 + 3.0 * phi) / n as f64; assert!(known < 0.0 && demeaned < 0.0, "{known} {demeaned}"); assert!( (known / predicted_known - 1.0).abs() < 0.15, "known mean: measured {known:.5} against {predicted_known:.5}" ); assert!( (demeaned / predicted_demeaned - 1.0).abs() < 0.15, "estimated mean: measured {demeaned:.5} against {predicted_demeaned:.5}" ); // Which is the chapter's "roughly doubles" at phi near one. Kendall's // leading terms predict (1 + 3 phi) / (2 phi) = 2.03 here; the measured // 2.36 is that plus the corrections of order 1/n^2 the formulas drop, // which is why the chapter says roughly rather than exactly. let ratio = demeaned / known; assert!((ratio - 2.36).abs() < 0.15, "the ratio was {ratio:.3}"); }} /// How often the unit root test finds mean reversion that is really there.////// The size of the test is fixed at five per cent by construction; what decides/// whether it is any use is its *power*, and the relative value chapter needs/// that number because it is the difference between a spread that can be/// verified and one that has to be taken on structural grounds.////// Returns the fraction of paths on which a genuinely mean-reverting series is/// correctly identified. Passing `half_life = 0` gives a random walk, so the/// answer is then the size rather than the power.pub fn unit_root_power(half_life: f64, years: f64, trials: usize, seed: u64) -> f64 { let dt = 1.0 / 252.0; let kappa = if half_life > 0.0 { (2.0f64).ln() / half_life } else { 0.0 }; let steps = (years / dt) as usize; let mut rng = crate::pathwise::Rng::new(seed); let mut rejects = 0; for _ in 0..trials { let path = ou_path(kappa, 0.0, 1.0, 0.0, dt, steps, &mut rng); if unit_root_test(&path, dt).rejects_random_walk { rejects += 1; } } rejects as f64 / trials as f64} #[cfg(test)]mod unit_root_power_tests { use super::*; /// The test is correctly sized: a true random walk is called mean reverting /// about five per cent of the time, which is what was asked of it. #[test] fn the_size_is_what_it_claims() { let size = unit_root_power(0.0, 5.0, 3_000, 4_242); assert!((size - 0.05).abs() < 0.015, "size was {size:.3}"); } /// And it is nearly powerless at the horizons a desk has. A spread with a /// one year half-life, watched for five years, is identified as reverting /// eight times in a hundred; at a two year half-life it is indistinguishable /// from a random walk over any history anyone holds. #[test] fn it_is_powerless_over_a_realistic_history() { let one_year = unit_root_power(1.0, 5.0, 3_000, 77); let two_year = unit_root_power(2.0, 10.0, 3_000, 78); assert!((one_year - 0.08).abs() < 0.03, "one year half-life gave {one_year:.3}"); assert!(two_year < 0.12, "two year half-life gave {two_year:.3}"); } /// What power depends on is the number of half-lives the sample spans, not /// the number of observations in it and not the calendar length. /// /// Five half-lives is five half-lives: a two year half-life over ten years /// and a one year half-life over five give the same answer, though one /// sample is twice the length of the other. #[test] fn power_is_a_function_of_half_lives_spanned() { let pairs = [((1.0, 5.0), (2.0, 10.0)), ((0.5, 5.0), (1.0, 10.0))]; for ((h1, y1), (h2, y2)) in pairs { let a = unit_root_power(h1, y1, 3_000, 311); let b = unit_root_power(h2, y2, 3_000, 312); assert!( (a - b).abs() < 0.035, "{h1}y over {y1}y gave {a:.3} but {h2}y over {y2}y gave {b:.3}" ); } // And it does rise once enough of them are in the sample. let twenty = unit_root_power(0.5, 10.0, 3_000, 313); assert!(twenty > 0.4, "twenty half-lives gave only {twenty:.3}"); }} /// Tracking a hedge ratio that moves, three ways.////// The relative value chapter's butterfly weights come from a covariance/// estimated over a window, and the window is a choice with no good answer: too/// short is noisy, too long is stale. The alternative is to stop choosing one/// and treat the ratio as a state that drifts,////// ```text/// beta[t] = beta[t-1] + eta, y[t] = beta[t] x[t] + epsilon,/// ```////// which is the filter of this chapter's Kalman section applied to a different/// state. Its gain settles at a value determined by the ratio of the two/// variances, so it is an exponentially weighted estimate whose effective memory/// is set by how fast the ratio is believed to move rather than by a window/// picked by hand.////// Returns the root mean squared hedge residual under a short window, a long/// window, and the filter, on the same simulated path.pub fn hedge_ratio_tracking( drift_sd: f64, // What the filter is told the drift is, which need not be the truth. assumed_drift_sd: f64, noise_sd: f64, steps: usize, short_window: usize, long_window: usize, seed: u64,) -> (f64, f64, f64) { let mut rng = crate::pathwise::Rng::new(seed); // A hedge ratio that genuinely wanders, and a regression that sees it only // through noise. let mut beta = 1.0; let (mut xs, mut ys, mut betas) = (Vec::new(), Vec::new(), Vec::new()); for _ in 0..steps { beta += drift_sd * rng.next_normal(); let x = rng.next_normal(); ys.push(beta * x + noise_sd * rng.next_normal()); xs.push(x); betas.push(beta); } let rolling = |window: usize| { let mut total = 0.0; let mut counted = 0; for t in window..steps { let (mut sxx, mut sxy) = (0.0, 0.0); for i in t - window..t { sxx += xs[i] * xs[i]; sxy += xs[i] * ys[i]; } let estimate = if sxx > 0.0 { sxy / sxx } else { 1.0 }; // The residual left by hedging the next observation with it. let residual = (betas[t] - estimate) * xs[t]; total += residual * residual; counted += 1; } (total / counted as f64).sqrt() }; // The filter, with the variances it is actually given. let filtered = { let (q, r) = (assumed_drift_sd * assumed_drift_sd, noise_sd * noise_sd); let (mut mean, mut variance) = (1.0, 1.0); let mut total = 0.0; let mut counted = 0; for t in 0..steps { // Predict: the state has wandered since the last observation. variance += q; if t >= long_window { let residual = (betas[t] - mean) * xs[t]; total += residual * residual; counted += 1; } // Update on this observation. let gain = variance * xs[t] / (xs[t] * xs[t] * variance + r); mean += gain * (ys[t] - mean * xs[t]); variance -= gain * xs[t] * variance; } (total / counted as f64).sqrt() }; (rolling(short_window), rolling(long_window), filtered)} #[cfg(test)]mod hedge_tracking_tests { use super::*; const NOISE: f64 = 0.5; const STEPS: usize = 4_000; const SHORT: usize = 30; const LONG: usize = 250; const SEED: u64 = 8_888; fn run(truth: f64, told: f64) -> (f64, f64, f64) { hedge_ratio_tracking(truth, told, NOISE, STEPS, SHORT, LONG, SEED) } /// No window is the right window: which of the two wins depends on how fast /// the ratio is moving, and that is not observed. #[test] fn the_better_window_depends_on_the_drift_rate() { let (slow_short, slow_long, _) = run(0.002, 0.002); assert!(slow_long < slow_short, "slow drift: {slow_long:.4} vs {slow_short:.4}"); let (fast_short, fast_long, _) = run(0.010, 0.010); assert!(fast_short < fast_long, "fast drift: {fast_short:.4} vs {fast_long:.4}"); } /// Told the truth, the filter beats either window at every rate, because it /// is not committed to one memory. #[test] fn a_correctly_specified_filter_beats_both_windows() { for drift in [0.002, 0.005, 0.010] { let (short, long, filter) = run(drift, drift); assert!( filter < short.min(long), "drift {drift}: filter {filter:.4} against {short:.4} and {long:.4}" ); } } /// And the caveat that keeps it honest. The advantage rests on knowing the /// rate, which is the window question wearing different clothes. Within /// about a factor of two it survives; five times too slow and the filter is /// worse than having simply taken the long window. #[test] fn the_advantage_does_not_survive_a_badly_wrong_rate() { let truth = 0.005; let best_window = run(truth, truth).1; for told in [0.0025, 0.010] { let filter = run(truth, told).2; assert!( filter <= best_window * 1.05, "told {told}: filter {filter:.4} against the best window {best_window:.4}" ); } let far_too_slow = run(truth, 0.001).2; assert!( far_too_slow > best_window, "five times too slow still beat the window: {far_too_slow:.4}" ); // Being told too fast is the more forgiving error of the two. assert!(run(truth, 0.010).2 < run(truth, 0.0025).2); }}