quant/src/rough.rs
A rough kernel as a sum of exponentials.
//! A rough kernel as a sum of exponentials.//!//! Rough volatility drives the variance with the Riemann-Liouville kernel//! `K(t) = t^(H - 1/2) / Gamma(H + 1/2)`, which has memory and is not Markov in//! any finite state. The smile dynamics chapter's escape is that for `H < 1/2`//! this kernel is a mixture of exponentials,//!//! ```text//! K(t) = integral of exp(-x t) mu(dx), mu(dx) = x^(-H - 1/2) dx / (Gamma(H + 1/2) Gamma(1/2 - H)),//! ```//!//! and each exponential kernel is an Ornstein-Uhlenbeck factor//! `dY = -x Y dt + dW`. Replace `mu` by a few atoms and the driven process//! `sum c_i Y_i` is Markov in `(Y_1, ..., Y_n)`.//!//! This module builds the atoms and measures what they reproduce. The finding//! that decides how the approximation is used: a kernel-shaped quantity//! converges with a modest number of nodes, but the exponent of the variance,//! which is what carries the roughness, is only recovered once the fastest node//! sits far above the reciprocal of the shortest horizon. use crate::special::ln_gamma; /// The exact kernel, `t^(H - 1/2) / Gamma(H + 1/2)`.pub fn kernel(hurst: f64, t: f64) -> f64 { t.powf(hurst - 0.5) / ln_gamma(hurst + 0.5).exp()} /// The variance of the exact process, `integral of K^2 = t^(2H) / (2H Gamma(H + 1/2)^2)`.pub fn exact_variance(hurst: f64, t: f64) -> f64 { t.powf(2.0 * hurst) / (2.0 * hurst * (2.0 * ln_gamma(hurst + 0.5)).exp())} /// `K_n(t) = sum c_i exp(-x_i t)`: weights and rates of the exponentials.pub struct ExponentialSum { pub weights: Vec<f64>, pub rates: Vec<f64>,} impl ExponentialSum { /// Atoms from cutting the spectral measure into geometric cells. /// /// The cells are `[0, x_min]` and then `n` cells with edges spaced /// geometrically up to `x_max`. Each cell becomes one exponential whose weight /// is the cell's mass and whose rate is the cell's mean, so the total mass and /// first moment of `mu` are kept exactly. The first cell is what carries the /// long end: a rate below `1 / t_max` is constant over the longest horizon. /// Mass above `x_max` is dropped, which is the only place the short end is lost. pub fn geometric(hurst: f64, n: usize, x_min: f64, x_max: f64) -> Self { assert!(hurst > 0.0 && hurst < 0.5, "the mixture exists only for H < 1/2"); assert!(n >= 1 && x_min > 0.0 && x_max > x_min); let a = 0.5 - hurst; let b = 1.5 - hurst; let norm = (ln_gamma(hurst + 0.5) + ln_gamma(a)).exp(); let mut edges = vec![0.0]; edges.extend((0..=n).map(|i| x_min * (x_max / x_min).powf(i as f64 / n as f64))); let (mut weights, mut rates) = (Vec::new(), Vec::new()); for cell in edges.windows(2) { let (lo, hi) = (cell[0], cell[1]); let mass = (hi.powf(a) - lo.powf(a)) / a / norm; let moment = (hi.powf(b) - lo.powf(b)) / b / norm; weights.push(mass); rates.push(moment / mass); } ExponentialSum { weights, rates } } /// The approximate kernel. pub fn at(&self, t: f64) -> f64 { self.weights.iter().zip(&self.rates).map(|(c, x)| c * (-x * t).exp()).sum() } /// `integral of K_n` over `[0, t]`. pub fn integrated(&self, t: f64) -> f64 { self.weights.iter().zip(&self.rates).map(|(c, x)| c * (1.0 - (-x * t).exp()) / x).sum() } /// The variance of `sum c_i Y_i` at time `t`, `integral of K_n^2`, in closed form. pub fn variance(&self, t: f64) -> f64 { let mut total = 0.0; for (ci, xi) in self.weights.iter().zip(&self.rates) { for (cj, xj) in self.weights.iter().zip(&self.rates) { total += ci * cj * (1.0 - (-(xi + xj) * t).exp()) / (xi + xj); } } total }} /// The local exponent `d log f / d log t` at `t`, by a symmetric difference in log time.pub fn local_exponent(f: impl Fn(f64) -> f64, t: f64) -> f64 { let h: f64 = 0.05; (f(t * h.exp()).ln() - f(t * (-h).exp()).ln()) / (2.0 * h)} #[cfg(test)]mod tests { use super::*; use crate::pathwise::Rng; const H: f64 = 0.1; const DAY: f64 = 1.0 / 252.0; fn log_grid(lo: f64, hi: f64, n: usize) -> Vec<f64> { (0..n).map(|i| lo * (hi / lo).powf(i as f64 / (n - 1) as f64)).collect() } /// The identity behind the whole construction, before any few-node /// approximation is trusted: with a fine cut of the spectral measure the sum /// of exponentials is the power-law kernel. #[test] fn the_kernel_is_a_mixture_of_exponentials() { let fine = ExponentialSum::geometric(H, 600, 1e-4, 1e12); for t in log_grid(1e-3, 10.0, 40) { let ratio = fine.at(t) / kernel(H, t); assert!((ratio - 1.0).abs() < 0.01, "at t = {t}: ratio {ratio}"); } } /// A couple of dozen nodes reproduce the kernel to a few per cent from a day /// to five years, the range the skew figure covers. #[test] fn a_few_nodes_reproduce_the_kernel_over_the_horizons_that_matter() { let approx = ExponentialSum::geometric(H, 24, 0.1, 1e7); let worst = log_grid(DAY, 5.0, 80) .into_iter() .map(|t| (approx.at(t) / kernel(H, t) - 1.0).abs()) .fold(0.0, f64::max); assert!(worst < 0.03, "worst relative error {worst:.4}"); // Too few nodes is visibly worse, so the check is not vacuous. let coarse = ExponentialSum::geometric(H, 4, 0.1, 1e7); let coarse_worst = log_grid(DAY, 5.0, 80) .into_iter() .map(|t| (coarse.at(t) / kernel(H, t) - 1.0).abs()) .fold(0.0, f64::max); assert!(coarse_worst > 3.0 * worst, "{coarse_worst:.4} against {worst:.4}"); } /// The roughness lives in the variance's exponent, `2H`, and that is the /// slower quantity to recover: the error is set by the share of variance /// contributed below the fastest node's time scale, which shrinks only like /// `(x_max t)^(-2H)` and `2H` is small. A ceiling near the reciprocal of the /// shortest horizon leaves the process visibly too smooth there. #[test] fn the_roughness_exponent_needs_nodes_far_above_the_shortest_horizon() { let week = 1.0 / 52.0; let target = 2.0 * H; let low = ExponentialSum::geometric(H, 16, 0.1, 2_000.0); let low_exponent = local_exponent(|t| low.variance(t), week); assert!( low_exponent > target + 0.04, "a ceiling of 2000 gave {low_exponent:.3}, too close to {target}" ); let high = ExponentialSum::geometric(H, 24, 0.1, 1e7); let high_exponent = local_exponent(|t| high.variance(t), week); assert!( (high_exponent - target).abs() < 0.02, "a ceiling of 1e7 gave {high_exponent:.3}, against {target}" ); // The exponent of the kernel's integral converges faster: at the same // low ceiling it is within a few hundredths of `H + 1/2`, while the // variance's exponent is off by more than a tenth. let integral_exponent = local_exponent(|t| low.integrated(t), week); assert!( (integral_exponent - (H + 0.5)).abs() < 0.04, "the integrated kernel's exponent was {integral_exponent:.3}" ); assert!(low_exponent - target > 0.1, "variance exponent {low_exponent:.3}"); } /// Where the missing variance comes from. The variance contributed at times /// shorter than `1 / x_max` is a share `(x_max t)^(-2H)` of the total, since /// `K^2` integrates to `tau^(2H)` up to `tau`, and a truncated kernel is smooth /// below that scale. The share the approximation actually misses is a /// constant fraction of that, because the saturated kernel gives some of it /// back, so it scales the same way. #[test] fn the_missing_variance_scales_like_the_fastest_rate_to_the_minus_two_h() { let week = 1.0 / 52.0; for x_max in [2e3, 1e5, 1e7, 1e9] { // Dense enough that the cell discretisation does not compete. let dense = ExponentialSum::geometric(H, 150, 0.1, x_max); let missing = 1.0 - dense.variance(week) / exact_variance(H, week); let ratio = missing / (x_max * week).powf(-2.0 * H); assert!( (0.7..1.0).contains(&ratio), "x_max = {x_max}: missing share {missing:.4}, ratio to the scaling law {ratio:.3}" ); } } /// The spectral measure is a pure power law, so it is self-similar, and a /// cell with a fixed ratio between its edges has the same relative error at /// every scale. So with the node ratio held fixed, covering a horizon range a /// hundred million times wider costs four times the nodes and no accuracy: /// the number of nodes grows with the logarithm of the ratio of horizons. #[test] fn the_error_is_the_same_at_every_scale_for_a_fixed_node_ratio() { let node_ratio = 1e8_f64.powf(1.0 / 24.0); let worst = |n: usize| -> f64 { let x_max = 0.1 * node_ratio.powi(n as i32); let approx = ExponentialSum::geometric(H, n, 0.1, x_max); log_grid(30.0 / x_max, 3.0, 81) .into_iter() .map(|t| (approx.at(t) / kernel(H, t) - 1.0).abs()) .fold(0.0, f64::max) }; let (small, medium, large) = (worst(12), worst(24), worst(48)); assert!(small < 0.02, "worst error {small:.4}"); for other in [medium, large] { assert!( (other / small - 1.0).abs() < 0.05, "errors differ across scales: {small:.4}, {medium:.4}, {large:.4}" ); } } /// The claim that the sum of exponentials is Markov in ordinary factors, /// checked by simulating the factors themselves. Each is an /// Ornstein-Uhlenbeck process driven by one shared Brownian motion, and the /// variance of their weighted sum has to match the closed form, which was /// derived without simulating anything. #[test] fn the_ordinary_factors_have_the_variance_the_kernel_gives() { let approx = ExponentialSum::geometric(H, 5, 0.5, 100.0); let horizon = 1.0; let steps = 1_000; let dt = horizon / steps as f64; let paths = 20_000; let mut rng = Rng::new(20260924); let (mut sum, mut sum_sq) = (0.0, 0.0); for _ in 0..paths { let mut y = vec![0.0; approx.rates.len()]; for _ in 0..steps { let dw = dt.sqrt() * rng.next_normal(); for (yi, x) in y.iter_mut().zip(&approx.rates) { *yi += -x * *yi * dt + dw; } } let value: f64 = approx.weights.iter().zip(&y).map(|(c, yi)| c * yi).sum(); sum += value; sum_sq += value * value; } let mean = sum / paths as f64; let variance = sum_sq / paths as f64 - mean * mean; let closed = approx.variance(horizon); assert!( (variance / closed - 1.0).abs() < 0.03, "simulated {variance:.5} against closed form {closed:.5}" ); }}