quant/src/roughheston.rs
The rough Heston model's Riccati equation, and its Markov lift.
//! The rough Heston model's Riccati equation, and its Markov lift.//!//! The rough Heston variance is//!//! ```text//! V_t = V_0 + integral_0^t K(t - s) [ lambda (theta - V_s) ds + nu sqrt(V_s) dB_s ],//! K(t) = t^(alpha - 1) / Gamma(alpha), alpha = H + 1/2,//! ```//!//! and `d log S = -V/2 dt + sqrt(V) dW` with `dW dB = rho dt`. The exponent of//! `E[S_T^a]` is affine in the state, as in the classical model, but the function//! that enters it solves a *Volterra integral equation* rather than an ordinary//! differential equation://!//! ```text//! h(t) = integral_0^t K(t - s) F(a, h(s)) ds,//! F(a, x) = (a^2 - a)/2 + (a rho nu - lambda) x + nu^2 x^2 / 2.//! ```//!//! For `H = 1/2` the kernel is one and this is `h' = F(a, h)`, Heston's own//! Riccati equation. For `H < 1/2` it has no closed form and is solved//! numerically. Replacing `K` by the sum of exponentials in [`crate::rough`]//! turns it back into a finite system of ordinary equations, one per factor,//! which is the sense in which the lifted model is affine. use crate::rough::ExponentialSum;use crate::special::ln_gamma; pub struct RoughHeston { pub hurst: f64, /// Mean reversion, `lambda`. pub reversion: f64, /// Long-run variance, `theta`. pub level: f64, /// Volatility of variance, `nu`. pub vol_of_vol: f64, pub correlation: f64, pub initial_variance: f64,} /// The trapezoid rule for samples on a uniform grid of step `dt`.fn trapezoid(samples: &[f64], dt: f64) -> f64 { let n = samples.len(); let inner: f64 = samples[1..n - 1].iter().sum(); dt * (inner + 0.5 * (samples[0] + samples[n - 1]))} impl RoughHeston { fn alpha(&self) -> f64 { self.hurst + 0.5 } /// `F(a, x)`, the right-hand side of the Riccati equation. pub fn riccati_rhs(&self, a: f64, x: f64) -> f64 { 0.5 * (a * a - a) + (a * self.correlation * self.vol_of_vol - self.reversion) * x + 0.5 * self.vol_of_vol * self.vol_of_vol * x * x } /// The solution `h` of the fractional Riccati equation on a uniform grid of /// `steps` intervals over `[0, horizon]`, by the Adams-Bashforth-Moulton /// predictor-corrector for Volterra equations (Diethelm, Ford and Freed). /// /// Quadratic cost in `steps`, which is the price of a kernel with memory: the /// whole history of `F(a, h)` enters every step. pub fn fractional_riccati(&self, a: f64, horizon: f64, steps: usize) -> Vec<f64> { let alpha = self.alpha(); let dt = horizon / steps as f64; let predictor_scale = dt.powf(alpha) / ln_gamma(alpha + 1.0).exp(); let corrector_scale = dt.powf(alpha) / ln_gamma(alpha + 2.0).exp(); let mut h = vec![0.0; steps + 1]; let mut f = vec![0.0; steps + 1]; f[0] = self.riccati_rhs(a, 0.0); for n in 0..steps { let nf = n as f64; let mut predictor = 0.0; for j in 0..=n { let jf = j as f64; predictor += ((nf + 1.0 - jf).powf(alpha) - (nf - jf).powf(alpha)) * f[j]; } let predicted = predictor_scale * predictor; let mut corrector = (nf.powf(alpha + 1.0) - (nf - alpha) * (nf + 1.0).powf(alpha)) * f[0]; for j in 1..=n { let d = nf - j as f64; corrector += ((d + 2.0).powf(alpha + 1.0) + d.powf(alpha + 1.0) - 2.0 * (d + 1.0).powf(alpha + 1.0)) * f[j]; } h[n + 1] = corrector_scale * (self.riccati_rhs(a, predicted) + corrector); f[n + 1] = self.riccati_rhs(a, h[n + 1]); } h } /// `log E[S_T^a]` for the rough model, from `h`: /// `integral_0^T F(a, h(tau)) g(T - tau) dtau` with /// `g(t) = V_0 + lambda theta t^alpha / Gamma(alpha + 1)`. pub fn log_moment(&self, a: f64, horizon: f64, steps: usize) -> f64 { let alpha = self.alpha(); let dt = horizon / steps as f64; let h = self.fractional_riccati(a, horizon, steps); let scale = self.reversion * self.level / ln_gamma(alpha + 1.0).exp(); let integrand: Vec<f64> = (0..=steps) .map(|i| { let remaining = horizon - i as f64 * dt; self.riccati_rhs(a, h[i]) * (self.initial_variance + scale * remaining.powf(alpha)) }) .collect(); trapezoid(&integrand, dt) } /// The same quantities for the lifted model, whose kernel is the sum of /// exponentials `kernel`: the vector `psi` of one component per factor solves /// `psi_i' = -x_i psi_i + c_i F(a, sum psi)`, and `sum psi` is the lifted /// counterpart of `h`. Stepped by a second-order exponential integrator, since /// the fastest rates make the system stiff and an explicit scheme would be /// unstable at any step worth using. /// /// Returns `(h_lifted on the grid, log E[S_T^a])`. pub fn lifted(&self, a: f64, horizon: f64, steps: usize, kernel: &ExponentialSum) -> (Vec<f64>, f64) { let dt = horizon / steps as f64; let n = kernel.rates.len(); // phi1 = (1 - e^-z)/x and phi2 = (e^-z - 1 + z)/(x^2 dt), with z = x dt, // by series where the subtraction would lose everything. let decay: Vec<f64> = kernel.rates.iter().map(|x| (-x * dt).exp()).collect(); let phi1: Vec<f64> = kernel .rates .iter() .map(|&x| { let z = x * dt; if z < 1e-3 { dt * (1.0 - z / 2.0 + z * z / 6.0) } else { (1.0 - (-z).exp()) / x } }) .collect(); let phi2: Vec<f64> = kernel .rates .iter() .map(|&x| { let z = x * dt; if z < 1e-3 { dt * (0.5 - z / 6.0 + z * z / 24.0) } else { ((-z).exp() - 1.0 + z) / (x * x * dt) } }) .collect(); let mut psi = vec![0.0; n]; let mut h = vec![0.0; steps + 1]; let mut f = vec![0.0; steps + 1]; f[0] = self.riccati_rhs(a, 0.0); for step in 0..steps { let total: f64 = psi.iter().sum(); let f_now = self.riccati_rhs(a, total); let predicted: f64 = (0..n) .map(|i| decay[i] * psi[i] + kernel.weights[i] * f_now * phi1[i]) .sum(); let f_next = self.riccati_rhs(a, predicted); for i in 0..n { psi[i] = decay[i] * psi[i] + kernel.weights[i] * (f_now * phi1[i] + (f_next - f_now) * phi2[i]); } h[step + 1] = psi.iter().sum(); f[step + 1] = self.riccati_rhs(a, h[step + 1]); } let scale = self.reversion * self.level; let integrand: Vec<f64> = (0..=steps) .map(|i| { let remaining = horizon - i as f64 * dt; f[i] * (self.initial_variance + scale * kernel.integrated(remaining)) }) .collect(); (h, trapezoid(&integrand, dt)) }} #[cfg(test)]mod tests { use super::*; use crate::pathwise::Rng; fn model(hurst: f64) -> RoughHeston { RoughHeston { hurst, reversion: 0.3, level: 0.02, vol_of_vol: 0.3, correlation: -0.7, initial_variance: 0.02, } } /// Heston's Riccati solution in closed form, for `a` with `a^2 - a <= 0`. fn classical_h(m: &RoughHeston, a: f64, t: f64) -> f64 { let nu2 = m.vol_of_vol * m.vol_of_vol; let beta = m.reversion - a * m.correlation * m.vol_of_vol; let d = (beta * beta - nu2 * (a * a - a)).sqrt(); let ratio = (beta - d) / (beta + d); (beta - d) / nu2 * (1.0 - (-d * t).exp()) / (1.0 - ratio * (-d * t).exp()) } /// At `H = 1/2` the kernel is one and the Volterra equation is Heston's own /// Riccati equation, whose solution is known in closed form. The exponent is /// then `V_0 h(T) + lambda theta integral h`, which is Heston's `C + D V_0`. #[test] fn at_one_half_the_fractional_equation_is_hestons_riccati() { let m = model(0.5); let (a, horizon, steps) = (0.5, 1.0, 400); let h = m.fractional_riccati(a, horizon, steps); for i in [50, 100, 200, 400] { let t = horizon * i as f64 / steps as f64; assert!( (h[i] - classical_h(&m, a, t)).abs() < 1e-6, "at t = {t}: {} against {}", h[i], classical_h(&m, a, t) ); } let fine = 4_000; let integral = trapezoid( &(0..=fine).map(|i| classical_h(&m, a, horizon * i as f64 / fine as f64)).collect::<Vec<_>>(), horizon / fine as f64, ); let expected = m.initial_variance * classical_h(&m, a, horizon) + m.reversion * m.level * integral; assert!( (m.log_moment(a, horizon, steps) - expected).abs() < 1e-7, "{} against {expected}", m.log_moment(a, horizon, steps) ); } /// The price is a martingale, so `E[S_T] = 1` for every `H`: `a = 1` makes /// `F(1, 0) = 0`, hence `h = 0`, and the exponent vanishes. Likewise `a = 0`. #[test] fn the_price_is_a_martingale_at_every_hurst_parameter() { for hurst in [0.1, 0.3, 0.5] { let m = model(hurst); for a in [0.0, 1.0] { assert!(m.log_moment(a, 1.0, 200).abs() < 1e-12, "H = {hurst}, a = {a}"); } } } /// The solver's own convergence, at a Hurst parameter with no closed form to /// compare against. #[test] fn the_fractional_solution_settles_as_the_grid_is_refined() { let m = model(0.1); let a = 0.5; let coarse = m.fractional_riccati(a, 1.0, 250); let fine = m.fractional_riccati(a, 1.0, 2000); assert!((coarse[250] - fine[2000]).abs() < 1e-6, "{} against {}", coarse[250], fine[2000]); assert!((m.log_moment(a, 1.0, 250) - m.log_moment(a, 1.0, 2000)).abs() < 1e-6); // And the answer is not the classical one: roughness moves it. assert!((fine[2000] - classical_h(&m, a, 1.0)).abs() > 5e-4, "H = 0.1 gave the classical solution"); } /// The lift is the fractional Riccati equation with the kernel replaced by /// its exponential sum, so refining the sum has to converge to it. Each /// doubling of the nodes divides the error by roughly four. #[test] fn the_factor_lift_converges_to_the_fractional_riccati() { let m = model(0.1); let a = 0.5; let exact = m.fractional_riccati(a, 1.0, 1000)[1000]; let errors: Vec<f64> = [12usize, 24, 48, 96] .iter() .map(|&n| { let kernel = ExponentialSum::geometric(0.1, n, 0.1, 1e7); (m.lifted(a, 1.0, 4000, &kernel).0[4000] - exact).abs() }) .collect(); for pair in errors.windows(2) { assert!(pair[1] < pair[0] / 2.5, "errors did not shrink fast enough: {errors:?}"); } assert!(errors[3] / exact.abs() < 0.002, "final relative error {:.5}", errors[3] / exact.abs()); } /// The claim that the lifted model is affine with those coefficients, checked /// by simulating the factors themselves. Nothing here uses the Riccati /// system, so agreement means the coefficients were derived correctly. #[test] fn the_lifted_moment_matches_a_simulation_of_the_factors() { let m = RoughHeston { hurst: 0.3, reversion: 0.5, level: 0.09, vol_of_vol: 0.5, correlation: -0.7, initial_variance: 0.09, }; let kernel = ExponentialSum::geometric(0.3, 4, 0.5, 50.0); let (a, horizon) = (2.0, 1.0); let (_, log_moment) = m.lifted(a, horizon, 2000, &kernel); let steps = 500; let dt = horizon / steps as f64; let paths = 100_000; let independent = (1.0 - m.correlation * m.correlation).sqrt(); let mut rng = Rng::new(20260925); let mut total = 0.0; for _ in 0..paths { let mut u = vec![0.0; kernel.rates.len()]; let mut x = 0.0; for step in 0..steps { let t = step as f64 * dt; let v = m.initial_variance + m.reversion * m.level * kernel.integrated(t) + kernel.weights.iter().zip(&u).map(|(c, ui)| c * ui).sum::<f64>(); let root = v.max(0.0).sqrt(); let z1 = rng.next_normal(); let z2 = rng.next_normal(); let dw = dt.sqrt() * z1; let db = dt.sqrt() * (m.correlation * z1 + independent * z2); x += -0.5 * v.max(0.0) * dt + root * dw; for (ui, r) in u.iter_mut().zip(&kernel.rates) { *ui += (-r * *ui - m.reversion * v.max(0.0)) * dt + m.vol_of_vol * root * db; } } total += (a * x).exp(); } let simulated = total / paths as f64; let closed = log_moment.exp(); assert!( (simulated / closed - 1.0).abs() < 0.01, "simulated {simulated:.5} against the Riccati system's {closed:.5}" ); }}