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 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))} #[cfg(test)]mod tests { use super::*; #[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 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()} #[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 ); }} // ---------------------------------------------------------------------------// 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) }} #[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) ); }} // ---------------------------------------------------------------------------// 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 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); }}