quant/src/generator.rs
The generator of a diffusion, and the duality that organises the rest.
//! The generator of a diffusion, and the duality that organises the rest.//!//! The generator chapter's claim is that one operator underlies most of what follows. For//! `dX = mu dt + sigma dW` it is//!//! ```text//! L f = mu f' + (1/2) sigma^2 f''//! ```//!//! and it is the drift of `f(X)`: Ito's lemma says exactly that. Pricing runs it//! one way, on functions of the state, and densities run its adjoint the other,//! and the two are related by integration by parts and nothing else.//!//! This module exists to make that checkable rather than decorative. The adjoint//! identity is verified on a grid, and the duality it implies --- that an//! expectation can be computed by evolving the payoff backwards or the density//! forwards, with the same answer --- is verified by doing both. /// A one-dimensional diffusion, as the two coefficient functions.pub struct Diffusion<M, S> { pub drift: M, pub volatility: S,} impl<M: Fn(f64) -> f64, S: Fn(f64) -> f64> Diffusion<M, S> { /// `L f` on the interior of a grid, by central differences. /// /// The coefficients sit *outside* the derivatives. That is the whole visual /// difference from the adjoint below, and it is what integration by parts /// moves. pub fn generator(&self, grid: &[f64], f: &[f64], dx: f64) -> Vec<f64> { let n = grid.len(); let mut out = vec![0.0; n]; for i in 1..n - 1 { let x = grid[i]; let first = (f[i + 1] - f[i - 1]) / (2.0 * dx); let second = (f[i + 1] - 2.0 * f[i] + f[i - 1]) / (dx * dx); let s = (self.volatility)(x); out[i] = (self.drift)(x) * first + 0.5 * s * s * second; } out } /// `L* p`, the adjoint, with the coefficients *inside* the derivatives. /// /// ```text /// L* p = -d/dx [ mu p ] + (1/2) d2/dx2 [ sigma^2 p ] /// ``` /// /// This is the right hand side of the Fokker-Planck equation of the local /// volatility chapter. Writing it beside [`Diffusion::generator`] is the /// point: the two differ only in whether the coefficients are inside, and /// that is exactly what one integration by parts does to each term. pub fn adjoint(&self, grid: &[f64], p: &[f64], dx: f64) -> Vec<f64> { let n = grid.len(); let flux: Vec<f64> = grid.iter().zip(p).map(|(&x, &q)| (self.drift)(x) * q).collect(); let spread: Vec<f64> = grid .iter() .zip(p) .map(|(&x, &q)| { let s = (self.volatility)(x); 0.5 * s * s * q }) .collect(); let mut out = vec![0.0; n]; for i in 1..n - 1 { let d_flux = (flux[i + 1] - flux[i - 1]) / (2.0 * dx); let d2_spread = (spread[i + 1] - 2.0 * spread[i] + spread[i - 1]) / (dx * dx); out[i] = -d_flux + d2_spread; } out } /// One explicit step of the backward evolution, `f <- f + dt * L f`. /// /// The Black-Scholes chapter's pricing equation run backwards from a payoff. Explicit /// stepping is used because the point is to display the operator, not to be /// a good solver; the step is kept inside the stability limit by the caller. pub fn step_backward(&self, grid: &[f64], f: &[f64], dx: f64, dt: f64) -> Vec<f64> { let lf = self.generator(grid, f, dx); let mut out: Vec<f64> = f.iter().zip(&lf).map(|(a, b)| a + dt * b).collect(); // Held fixed at the ends, which is harmless because everything here is // supported well inside the domain. out[0] = f[0]; let n = out.len(); out[n - 1] = f[n - 1]; out } /// One explicit step of the forward evolution, `p <- p + dt * L* p`. pub fn step_forward(&self, grid: &[f64], p: &[f64], dx: f64, dt: f64) -> Vec<f64> { let lp = self.adjoint(grid, p, dx); let mut out: Vec<f64> = p.iter().zip(&lp).map(|(a, b)| a + dt * b).collect(); out[0] = 0.0; let n = out.len(); out[n - 1] = 0.0; out }} /// A uniform grid and its spacing.pub fn uniform_grid(lo: f64, hi: f64, points: usize) -> (Vec<f64>, f64) { let dx = (hi - lo) / (points - 1) as f64; ((0..points).map(|i| lo + i as f64 * dx).collect(), dx)} /// The `L^2` pairing `integral of f p`, by the trapezoid rule.pub fn pairing(f: &[f64], p: &[f64], dx: f64) -> f64 { let n = f.len(); let mut total = 0.5 * (f[0] * p[0] + f[n - 1] * p[n - 1]); for i in 1..n - 1 { total += f[i] * p[i]; } total * dx} /// Probabilists' Hermite polynomial `He_n`, by the standard recurrence.////// These are the eigenfunctions of the Ornstein-Uhlenbeck generator, which is/// the one case where the solvable models chapter's spectral picture can be written down in/// full rather than computed.pub fn hermite(n: usize, y: f64) -> f64 { let (mut previous, mut current) = (1.0, y); if n == 0 { return previous; } for k in 1..n { let next = y * current - k as f64 * previous; previous = current; current = next; } current} /// The Ornstein-Uhlenbeck semigroup by eigenfunction expansion.////// ```text/// P_t f = sum over n of c_n exp(-n kappa t) He_n(x/s)/// ```////// with `s^2 = sigma^2 / (2 kappa)` the stationary variance and `c_n` the/// coefficients of `f` in the Hermite basis. Every maturity costs one/// exponential once the coefficients are known, which is the practical appeal:/// the expensive part does not depend on `t`.pub fn ou_spectral(coefficients: &[f64], kappa: f64, s: f64, t: f64, x: f64) -> f64 { coefficients .iter() .enumerate() .map(|(n, c)| c * (-(n as f64) * kappa * t).exp() * hermite(n, x / s)) .sum()} /// Hermite coefficients of `f` against the stationary Gaussian.////// `c_n = <f, He_n> / n!`, the inner product taken against the density the/// process settles into, which is the measure in which the generator is/// self-adjoint and the eigenfunctions are orthogonal.pub fn hermite_coefficients(f: impl Fn(f64) -> f64, s: f64, terms: usize) -> Vec<f64> { const LIMIT: f64 = 9.0; const STEPS: usize = 40_000; let h = 2.0 * LIMIT / STEPS as f64; (0..terms) .map(|n| { let mut total = 0.0; let mut factorial = 1.0; for k in 1..=n { factorial *= k as f64; } for i in 0..=STEPS { let y = -LIMIT + i as f64 * h; let weight = if i == 0 || i == STEPS { 0.5 } else { 1.0 }; let density = (-0.5 * y * y).exp() / (2.0 * std::f64::consts::PI).sqrt(); total += weight * f(y * s) * hermite(n, y) * density; } total * h / factorial }) .collect()} #[cfg(test)]mod tests { use super::*; /// An Ornstein-Uhlenbeck generator, which is the running example of the /// chapter and has a state-dependent drift, so the adjoint is genuinely /// different from the generator rather than accidentally equal. fn ou() -> Diffusion<impl Fn(f64) -> f64, impl Fn(f64) -> f64> { Diffusion { drift: |x: f64| -1.5 * x, volatility: |_x: f64| 0.6 } } /// A diffusion whose volatility varies with the state, so that the term the /// adjoint moves inside the derivative actually contributes. fn state_dependent() -> Diffusion<impl Fn(f64) -> f64, impl Fn(f64) -> f64> { Diffusion { drift: |x: f64| 0.4 - 0.8 * x, volatility: |x: f64| 0.3 + 0.15 * x.tanh() } } fn bump(centre: f64, width: f64) -> impl Fn(f64) -> f64 { move |x: f64| (-((x - centre) / width).powi(2)).exp() } #[test] fn the_hermite_polynomials_are_eigenfunctions_of_the_ou_generator() { // The spectral claim, checked directly rather than by finding // eigenvalues numerically: L He_n(x/s) = -n kappa He_n(x/s), so the // generator is diagonal in this basis and the eigenvalues are the // integers times the mean reversion rate. let (kappa, sigma) = (1.5f64, 0.6f64); let s = (sigma * sigma / (2.0 * kappa)).sqrt(); let d = Diffusion { drift: |x: f64| -1.5 * x, volatility: |_x: f64| 0.6 }; let (grid, dx) = uniform_grid(-4.0 * s, 4.0 * s, 20_001); for n in 1..=5 { let f: Vec<f64> = grid.iter().map(|&x| hermite(n, x / s)).collect(); let lf = d.generator(&grid, &f, dx); // Compare on the interior, away from the differencing edges. let lo = grid.len() / 5; let hi = grid.len() - lo; let mut worst: f64 = 0.0; for i in lo..hi { let expected = -(n as f64) * kappa * f[i]; worst = worst.max((lf[i] - expected).abs() / (1.0 + expected.abs())); } // The tolerance loosens with the degree because the second difference // is exact only up to cubics; beyond that it carries an O(h^2) error // scaled by the fourth derivative, which grows with n. The identity // is exact -- the grid is not. assert!(worst < 2e-5, "He_{n} was not an eigenfunction: worst {worst}"); } } #[test] fn the_eigenfunction_expansion_reproduces_the_semigroup() { // End to end. Expand a payoff in Hermites, decay each coefficient by // exp(-n kappa t), and resum -- against the exact Gaussian transition // law integrated numerically. The two share no step beyond the // definition of the process. let (kappa, sigma) = (1.2f64, 0.5f64); let s = (sigma * sigma / (2.0 * kappa)).sqrt(); let payoff = |x: f64| (x - 0.1).max(0.0); let coefficients = hermite_coefficients(payoff, s, 40); for &t in &[0.15f64, 0.5, 1.5] { for &x in &[-0.3f64, 0.0, 0.25] { let spectral = ou_spectral(&coefficients, kappa, s, t, x); // Exact transition: N(x e^{-kappa t}, s^2 (1 - e^{-2 kappa t})). let mean = x * (-kappa * t).exp(); let sd = s * (1.0 - (-2.0 * kappa * t).exp()).sqrt(); let (steps, limit) = (200_000, 9.0); let h = 2.0 * limit / steps as f64; let mut quadrature = 0.0; for i in 0..=steps { let z = -limit + i as f64 * h; let w = if i == 0 || i == steps { 0.5 } else { 1.0 }; let density = (-0.5 * z * z).exp() / (2.0 * std::f64::consts::PI).sqrt(); quadrature += w * payoff(mean + sd * z) * density; } quadrature *= h; assert!( (spectral - quadrature).abs() < 2e-4, "at t={t}, x={x}: spectral {spectral} against quadrature {quadrature}" ); } } } #[test] fn the_spectral_gap_sets_the_rate_of_forgetting() { // The eigenvalue that matters, stated correctly. The gap does not say // that any particular deviation halves in a half-life -- at short times // every mode contributes. It says the *slowest surviving* mode decays at // exp(-kappa t), so once the faster ones have died the whole deviation // decays at that rate, and that asymptotic rate is what a half-life // means. let kappa = 1.2f64; let s = (0.5f64 * 0.5 / (2.0 * kappa)).sqrt(); let coefficients = hermite_coefficients(|x: f64| (x - 0.1).max(0.0), s, 40); let deviation = |t: f64| { (ou_spectral(&coefficients, kappa, s, t, 0.6) - coefficients[0]).abs() }; // Late enough that only the first mode is left. At t = 3 the second // mode is still about three percent of the first and the ratio comes out // half a percent low, which is the transient rather than an error. let ratio = deviation(7.0) / deviation(6.0); let predicted = (-kappa * 1.0f64).exp(); assert!( (ratio - predicted).abs() < 1e-3, "asymptotic decay {ratio} against exp(-kappa) = {predicted}" ); // And at short times it is faster than the gap alone, because the higher // modes have not yet died -- which is why the naive half-life reading is // wrong. let early = deviation(0.2) / deviation(0.0); assert!( early < (-kappa * 0.2f64).exp(), "early decay {early} was not faster than the gap predicts" ); } #[test] fn the_adjoint_is_the_adjoint() { // The identity the chapter is built on, and the only content of the // phrase "integration by parts": <L f, p> = <f, L* p>, provided both // vanish at the ends. Checked with a state-dependent volatility, so the // terms that move inside the derivative are not zero. let d = state_dependent(); let (grid, dx) = uniform_grid(-8.0, 8.0, 4001); let f: Vec<f64> = grid.iter().map(|&x| bump(-0.7, 1.3)(x)).collect(); let p: Vec<f64> = grid.iter().map(|&x| bump(0.9, 1.1)(x)).collect(); let lhs = pairing(&d.generator(&grid, &f, dx), &p, dx); let rhs = pairing(&f, &d.adjoint(&grid, &p, dx), dx); assert!( (lhs - rhs).abs() < 1e-6 * lhs.abs().max(1.0), "<Lf,p> = {lhs} against <f,L*p> = {rhs}" ); } #[test] fn putting_the_coefficients_outside_breaks_it() { // The mistake the chapter warns about, made deliberately. Using L in // place of L* -- the coefficients left outside the derivatives -- fails // the same identity by a wide margin, which is why the placement is not // a matter of taste. let d = state_dependent(); let (grid, dx) = uniform_grid(-8.0, 8.0, 4001); let f: Vec<f64> = grid.iter().map(|&x| bump(-0.7, 1.3)(x)).collect(); let p: Vec<f64> = grid.iter().map(|&x| bump(0.9, 1.1)(x)).collect(); let correct = pairing(&f, &d.adjoint(&grid, &p, dx), dx); let wrong = pairing(&f, &d.generator(&grid, &p, dx), dx); assert!( (correct - wrong).abs() > 0.05 * correct.abs(), "the two agreed, so the example cannot show the difference" ); } #[test] fn the_generator_is_the_drift_of_f_of_x() { // Ito's lemma, read as a definition: L f is the rate at which E[f(X)] // moves. Checked against the definition itself, the short-time limit of // the expectation, computed by evolving a point mass forward. let d = ou(); let (grid, dx) = uniform_grid(-6.0, 6.0, 2001); let f: Vec<f64> = grid.iter().map(|&x| (0.7 * x).sin()).collect(); let analytic = d.generator(&grid, &f, dx); // Start from a narrow bump at the origin and evolve it a little. let start = 0.0; let width = 0.25; let mut p: Vec<f64> = grid.iter().map(|&x| bump(start, width)(x)).collect(); let mass = pairing(&vec![1.0; grid.len()], &p, dx); for q in p.iter_mut() { *q /= mass; } let before = pairing(&f, &p, dx); let (dt, steps) = (2e-5, 200); for _ in 0..steps { p = d.step_forward(&grid, &p, dx, dt); } let after = pairing(&f, &p, dx); let measured = (after - before) / (dt * steps as f64); // The bump is not a point mass, so the comparison is against L f // averaged over it rather than at a point. let expected = pairing(&analytic, &{ let mut q: Vec<f64> = grid.iter().map(|&x| bump(start, width)(x)).collect(); let m = pairing(&vec![1.0; grid.len()], &q, dx); for v in q.iter_mut() { *v /= m; } q }, dx); assert!( (measured - expected).abs() < 0.02 * expected.abs().max(0.1), "measured drift {measured} against L f averaged {expected}" ); } #[test] fn an_expectation_can_be_computed_from_either_end() { // The duality, which is the chapter's practical point. The same number // -- E[f(X_T)] -- comes out of evolving the payoff backwards against the // starting density, or evolving the density forwards against the payoff. // One operator, two directions, and the choice is a matter of what is // being varied rather than of what is true. let d = ou(); let (grid, dx) = uniform_grid(-7.0, 7.0, 3001); let f: Vec<f64> = grid.iter().map(|&x| bump(0.0, 2.0)(x)).collect(); let mut p: Vec<f64> = grid.iter().map(|&x| bump(1.2, 0.8)(x)).collect(); let mass = pairing(&vec![1.0; grid.len()], &p, dx); for q in p.iter_mut() { *q /= mass; } let p0 = p.clone(); let (dt, steps) = (1e-5, 4000); let mut u = f.clone(); for _ in 0..steps { u = d.step_backward(&grid, &u, dx, dt); p = d.step_forward(&grid, &p, dx, dt); } let forwards = pairing(&f, &p, dx); let backwards = pairing(&u, &p0, dx); assert!( (forwards - backwards).abs() < 2e-4, "forward route {forwards} against backward route {backwards}" ); } #[test] fn the_forward_evolution_conserves_probability() { // As the local volatility chapter argues it must: the adjoint is a // divergence, so it moves mass rather than creating it. let d = state_dependent(); let (grid, dx) = uniform_grid(-9.0, 9.0, 3001); let mut p: Vec<f64> = grid.iter().map(|&x| bump(0.0, 1.0)(x)).collect(); let ones = vec![1.0; grid.len()]; let mass = pairing(&ones, &p, dx); for q in p.iter_mut() { *q /= mass; } for _ in 0..3000 { p = d.step_forward(&grid, &p, dx, 1e-5); } let after = pairing(&ones, &p, dx); assert!((after - 1.0).abs() < 1e-4, "mass drifted to {after}"); }} /// The generator's matrix on polynomials of degree at most `degree`.////// The solvable models chapter shows that an invariant family which is a linear/// *subspace* turns the generator into a matrix, and the flow inside it into a/// linear system. This builds that matrix, for the only coefficients that can/// preserve polynomials: affine drift `mu(x) = a + b x` and quadratic variance/// `v(x) = c + d x + e x^2`.////// Applying the generator to a single power,////// ```text/// L x^n = n (a + b x) x^{n-1} + (n(n-1)/2) (c + d x + e x^2) x^{n-2}/// = [b n + e n(n-1)/2] x^n/// + [a n + d n(n-1)/2] x^{n-1}/// + [ c n(n-1)/2] x^{n-2},/// ```////// so column `n` of the matrix has at most three entries and never reaches/// above row `n`. The matrix is upper triangular in the usual orientation, its/// diagonal is `b n + e n(n-1)/2`, and those are the generator's eigenvalues on/// this subspace --- for an Ornstein-Uhlenbeck process, `-kappa n`, which is the/// spectrum the chapter meets again through Hermite polynomials.////// Returned as `m[out][in]`, so a coefficient vector evolves by `a' = m a`.pub fn polynomial_generator_matrix( drift: (f64, f64), variance: (f64, f64, f64), degree: usize,) -> Vec<Vec<f64>> { let (a, b) = drift; let (c, d, e) = variance; let mut m = vec![vec![0.0; degree + 1]; degree + 1]; for n in 0..=degree { let (n_f, pairs) = (n as f64, (n * n.saturating_sub(1)) as f64 / 2.0); m[n][n] = b * n_f + e * pairs; if n >= 1 { m[n - 1][n] = a * n_f + d * pairs; } if n >= 2 { m[n - 2][n] = c * pairs; } } m} /// The moments `E[x_t^n]` for `n` up to `degree`, by integrating the linear/// system the previous function produces.////// Moments move by the transpose: `d/dt E[x^n] = E[L x^n] = sum_m m[m][n] E[x^m]`./// Runge-Kutta, because the point is that an ODE is all that is left --- no/// simulation and no grid.pub fn polynomial_moments( drift: (f64, f64), variance: (f64, f64, f64), start: f64, horizon: f64, degree: usize, steps: usize,) -> Vec<f64> { let m = polynomial_generator_matrix(drift, variance, degree); let apply = |v: &[f64]| -> Vec<f64> { (0..=degree) .map(|n| (0..=degree).map(|k| m[k][n] * v[k]).sum()) .collect() }; let mut moments: Vec<f64> = (0..=degree).map(|n| start.powi(n as i32)).collect(); let dt = horizon / steps as f64; for _ in 0..steps { let k1 = apply(&moments); let mid1: Vec<f64> = moments.iter().zip(&k1).map(|(m, k)| m + 0.5 * dt * k).collect(); let k2 = apply(&mid1); let mid2: Vec<f64> = moments.iter().zip(&k2).map(|(m, k)| m + 0.5 * dt * k).collect(); let k3 = apply(&mid2); let end: Vec<f64> = moments.iter().zip(&k3).map(|(m, k)| m + dt * k).collect(); let k4 = apply(&end); for n in 0..=degree { moments[n] += dt / 6.0 * (k1[n] + 2.0 * k2[n] + 2.0 * k3[n] + k4[n]); } } moments} /// The moments of a normal distribution, by the standard recursion/// `E[X^n] = m E[X^{n-1}] + (n-1) v E[X^{n-2}]`.pub fn gaussian_moments(mean: f64, variance: f64, degree: usize) -> Vec<f64> { let mut moments = vec![1.0]; for n in 1..=degree { let previous = moments[n - 1]; let two_back = if n >= 2 { moments[n - 2] } else { 0.0 }; moments.push(mean * previous + (n as f64 - 1.0) * variance * two_back); } moments} #[cfg(test)]mod invariant_family_tests { use super::*; const KAPPA: f64 = 0.8; const SIGMA: f64 = 0.35; const START: f64 = 0.6; /// Ornstein-Uhlenbeck: `mu(x) = -kappa x`, `v(x) = sigma^2`. fn ou() -> ((f64, f64), (f64, f64, f64)) { ((0.0, -KAPPA), (SIGMA * SIGMA, 0.0, 0.0)) } #[test] fn the_polynomial_matrix_is_triangular_with_the_ou_spectrum_on_it() { let (drift, variance) = ou(); let m = polynomial_generator_matrix(drift, variance, 5); for n in 0..=5 { assert!( (m[n][n] + KAPPA * n as f64).abs() < 1e-12, "diagonal at {n} is {}, expected {}", m[n][n], -KAPPA * n as f64 ); // Nothing below the diagonal: the generator cannot raise a degree. for row in (n + 1)..=5 { assert_eq!(m[row][n], 0.0, "entry [{row}][{n}] should vanish"); } } // And the only other entries are two rows up, from the second // derivative: L x^n = -kappa n x^n + sigma^2 n(n-1)/2 x^{n-2}. for n in 2..=5 { let expected = SIGMA * SIGMA * (n * (n - 1)) as f64 / 2.0; assert!((m[n - 2][n] - expected).abs() < 1e-12); if n >= 1 { assert_eq!(m[n - 1][n], 0.0, "no odd coupling for a zero drift constant"); } } } #[test] fn integrating_the_linear_system_gives_the_right_moments() { // The claim in full: because polynomials are invariant, every moment is // available from an ODE rather than from a simulation. Checked against // the exact Gaussian moments, which share none of the machinery. let (drift, variance) = ou(); let degree = 6; for horizon in [0.25, 1.0, 3.0] { let solved = polynomial_moments(drift, variance, START, horizon, degree, 4000); let mean = START * (-KAPPA * horizon).exp(); let var = SIGMA * SIGMA * (1.0 - (-2.0 * KAPPA * horizon).exp()) / (2.0 * KAPPA); let exact = gaussian_moments(mean, var, degree); for n in 0..=degree { assert!( (solved[n] - exact[n]).abs() < 1e-9 * exact[n].abs().max(1.0), "t={horizon}, n={n}: solved {} against exact {}", solved[n], exact[n] ); } } } /// One step of `(u,v)' = M (u,v)` by Runge-Kutta. fn linear_step(m: [[f64; 2]; 2], p: (f64, f64), dt: f64) -> (f64, f64) { let f = |(u, v): (f64, f64)| (m[0][0] * u + m[0][1] * v, m[1][0] * u + m[1][1] * v); let k1 = f(p); let k2 = f((p.0 + 0.5 * dt * k1.0, p.1 + 0.5 * dt * k1.1)); let k3 = f((p.0 + 0.5 * dt * k2.0, p.1 + 0.5 * dt * k2.1)); let k4 = f((p.0 + dt * k3.0, p.1 + dt * k3.1)); ( p.0 + dt / 6.0 * (k1.0 + 2.0 * k2.0 + 2.0 * k3.0 + k4.0), p.1 + dt / 6.0 * (k1.1 + 2.0 * k2.1 + 2.0 * k3.1 + k4.1), ) } #[test] fn the_ratio_of_a_linear_flow_obeys_a_riccati_equation() { // The projection the solvable models chapter is built on. With // (u,v)' = M (u,v) and psi = v/u, the quotient rule gives // // psi' = v'/u - psi u'/u = c + (d - a) psi - b psi^2, // // the quadratic arriving from the second term, where the ratio // multiplies its own denominator's growth rate. let m = [[0.4, -0.9], [1.3, -0.2]]; let (a, b, c, d) = (m[0][0], m[0][1], m[1][0], m[1][1]); let (mut u, mut v) = (1.0f64, 0.3f64); let dt = 1e-6; let mut t = 0.0; while t < 0.4 { let psi = v / u; let (u_next, v_next) = linear_step(m, (u, v), dt); let measured = (v_next / u_next - psi) / dt; let predicted = c + (d - a) * psi - b * psi * psi; assert!( (measured - predicted).abs() < 1e-4 * predicted.abs().max(1.0), "t={t}: ratio moved at {measured}, Riccati says {predicted}" ); u = u_next; v = v_next; t += dt; } } #[test] fn the_riccati_equilibria_are_the_eigendirections() { // The projection is the blow-up of the plane at the origin, and this is // what that buys: the single singularity of the linear flow at 0 becomes // isolated singularities on the exceptional line, one per eigendirection. // // For psi' = c + (d-a) psi - b psi^2 the roots are psi = (lambda - a)/b, // because substituting turns the quadratic into the characteristic // polynomial (lambda - a)(lambda - d) - bc divided by b. let m = [[0.4f64, 0.9], [0.3, -0.2]]; let (a, b, c, d) = (m[0][0], m[0][1], m[1][0], m[1][1]); let (trace, det) = (a + d, a * d - b * c); let discriminant = trace * trace - 4.0 * det; assert!(discriminant > 0.0, "want two real eigendirections for this test"); let root = discriminant.sqrt(); for lambda in [0.5 * (trace + root), 0.5 * (trace - root)] { let psi = (lambda - a) / b; let drift = c + (d - a) * psi - b * psi * psi; assert!( drift.abs() < 1e-12, "eigenvalue {lambda} gives direction {psi}, where the flow moves at {drift}" ); } } #[test] fn the_riccati_flow_is_a_mobius_map_of_its_starting_point() { // The other consequence of the projection. The solution map is induced // by exp(T M) acting on the plane, and a linear map of the plane is a // Mobius map of the line of ratios. Mobius maps are exactly the ones // preserving the cross-ratio, so integrating four starting points // forward and taking their cross-ratio should return what it started as // -- an invariant no property of the differential equation alone would // suggest. let m = [[0.4f64, -0.9], [1.3, -0.2]]; let (a, b, c, d) = (m[0][0], m[0][1], m[1][0], m[1][1]); let advance = |mut psi: f64| { let (dt, steps) = (1e-6, 300_000); let f = |p: f64| c + (d - a) * p - b * p * p; for _ in 0..steps { let k1 = f(psi); let k2 = f(psi + 0.5 * dt * k1); let k3 = f(psi + 0.5 * dt * k2); let k4 = f(psi + dt * k3); psi += dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4); } psi }; let cross = |z: [f64; 4]| { ((z[0] - z[2]) * (z[1] - z[3])) / ((z[0] - z[3]) * (z[1] - z[2])) }; let start = [0.1, 0.45, -0.3, 1.2]; let finish = [advance(start[0]), advance(start[1]), advance(start[2]), advance(start[3])]; // The points genuinely moved, so this is not invariance by accident. for (s, f) in start.iter().zip(&finish) { assert!((s - f).abs() > 0.05, "start {s} barely moved, reaching {f}"); } let (before, after) = (cross(start), cross(finish)); assert!( (before - after).abs() < 1e-6 * before.abs(), "cross-ratio {before:.12} became {after:.12}" ); } /// Why the quadratic is a ceiling and not a convenience. /// /// Vector fields on a line bracket as `[f d, g d] = (f g' - g f') d`, so for /// monomials `[x^i d, x^j d] = (j - i) x^(i+j-1) d`. On `{1, x, x^2}` that /// never leaves the span --- the highest it produces is `x^2` --- so those /// three fields close into a Lie algebra, the one whose group acts on the /// projective line by Mobius maps. Admit a cubic and it escapes at once, /// since `[x^2 d, x^3 d]` is a quartic, which brackets to a quintic, and so /// on without end. /// /// The consequence is checkable rather than abstract. A Riccati flow is a /// Mobius map and so preserves the cross-ratio; a cubic flow belongs to no /// finite-dimensional group and preserves nothing. Same integrator, same /// four starting points, one extra term. #[test] fn a_cubic_term_destroys_the_cross_ratio_a_riccati_preserves() { let advance = |mut psi: f64, cubic: f64| { let (dt, steps) = (1e-5, 30_000); let f = |p: f64| 0.4 - 0.3 * p + 0.5 * p * p + cubic * p * p * p; for _ in 0..steps { let k1 = f(psi); let k2 = f(psi + 0.5 * dt * k1); let k3 = f(psi + 0.5 * dt * k2); let k4 = f(psi + dt * k3); psi += dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4); } psi }; let cross = |z: [f64; 4]| ((z[0] - z[2]) * (z[1] - z[3])) / ((z[0] - z[3]) * (z[1] - z[2])); let start = [0.1, 0.45, -0.3, 1.2]; let run = |cubic: f64| { let f: Vec<f64> = start.iter().map(|&s| advance(s, cubic)).collect(); assert!(f.iter().all(|x| x.is_finite()), "the flow escaped at cubic {cubic}"); (cross(f.as_slice().try_into().unwrap()) / cross(start) - 1.0).abs() }; // No cubic term: a Riccati, and the cross-ratio survives to machine // precision. assert!(run(0.0) < 1e-12, "the quadratic case drifted by {:.2e}", run(0.0)); // With one, it does not, and the damage grows with the coefficient. let gentle = run(-0.2); let firm = run(-0.5); assert!((gentle - 0.0305).abs() < 0.003, "cubic -0.2 moved it by {gentle:.4}"); assert!((firm - 0.0698).abs() < 0.005, "cubic -0.5 moved it by {firm:.4}"); assert!(firm > gentle); } #[test] fn the_riccati_solution_blows_up_where_the_linear_one_does_not() { // The consequence worth carrying: the ratio explodes when its // denominator passes through zero, while nothing at all happens to the // linear flow underneath. Not a numerical failure, and no step size // fixes it. let m = [[0.0, -1.0], [1.0, 0.0]]; // a rotation: u = cos t, v = sin t let (mut u, mut v) = (1.0f64, 0.0f64); let dt = 1e-5; let (mut biggest_ratio, mut t) = (0.0f64, 0.0f64); while t < 1.6 { let (nu, nv) = linear_step(m, (u, v), dt); u = nu; v = nv; t += dt; // The pair stays on the unit circle throughout. assert!(((u * u + v * v) - 1.0).abs() < 1e-8, "the linear flow is bounded"); biggest_ratio = biggest_ratio.max((v / u).abs()); } // And yet the ratio has run past any bound as u approaches zero at // t = pi/2, which is inside the interval just integrated. assert!(biggest_ratio > 1e4, "the ratio should explode, reached {biggest_ratio:.3e}"); } #[test] fn a_parametrisation_absorbs_exactly_the_span_of_its_exponent() { // What "the parametrisation can absorb" means, made checkable. For a // family u = exp(sum theta_i g_i), the multipliers the parameters can // produce are exactly the combinations of the g_i, since // d/dt log u = sum theta_i' g_i. So invariance is: L u / u lies in that // span. // // Take the quadratic exponent q = phi + psi x + gamma x^2, whose span is // {1, x, x^2}. Then // // L u / u = mu q' + v (q'^2 + q'') / 2, // // and the degrees say deg mu <= 1 and deg v = 0: a quadratic exponent // needs a *constant* variance, which is stronger than the affine family // asks for. The test checks the identity, and then that a linear // variance really does push the ratio out of the span. let (psi, gamma) = (0.4f64, -0.25f64); let q_prime = |x: f64| psi + 2.0 * gamma * x; let q_second = 2.0 * gamma; let (grid, dx) = uniform_grid(-1.5, 1.5, 4001); let u: Vec<f64> = grid.iter().map(|&x| (psi * x + gamma * x * x).exp()).collect(); // A constant variance keeps the ratio inside the span; a linear one // does not. for (name, variance, stays) in [ ("constant variance", (0.3f64, 0.0f64), true), ("linear variance", (0.3, 0.22), false), ] { let (v0, v1) = variance; let diffusion = Diffusion { drift: |x: f64| 0.1 - 0.6 * x, volatility: |x: f64| (v0 + v1 * x).max(1e-12).sqrt(), }; let lu = diffusion.generator(&grid, &u, dx); // The algebraic identity, wherever the differencing is two-sided. for i in (400..grid.len() - 400).step_by(211) { let x = grid[i]; let predicted = (0.1 - 0.6 * x) * q_prime(x) + 0.5 * (v0 + v1 * x) * (q_prime(x) * q_prime(x) + q_second); assert!( (lu[i] / u[i] - predicted).abs() < 1e-6 * predicted.abs().max(1.0), "{name} at x={x}: ratio {} against {predicted}", lu[i] / u[i] ); } // And whether it is a quadratic. A third difference annihilates // quadratics and not cubics, which is exactly the distinction // wanted: v q'^2 is cubic when v is linear. let ratio: Vec<f64> = lu.iter().zip(&u).map(|(l, u)| l / u).collect(); let (i, step) = (grid.len() / 2, 400); let third = ratio[i + 2 * step] - 3.0 * ratio[i + step] + 3.0 * ratio[i] - ratio[i - step]; if stays { assert!(third.abs() < 1e-6, "{name} should stay quadratic, got {third:.2e}"); } else { assert!(third.abs() > 1e-3, "{name} should leave the span, got {third:.2e}"); } } } #[test] fn the_exponential_family_makes_the_parameters_move_quadratically() { // The other half of the chapter's distinction. For u = exp(phi + psi x), // // L u / u = mu(x) psi + v(x) psi^2 / 2, // // which for affine mu and v is affine in x. Matching the two powers // against d/dt log u = phi' + psi' x gives // // phi' = a psi + c psi^2 / 2, psi' = b psi + d psi^2 / 2, // // and the psi^2 --- the reason the equation is Riccati and not linear // --- comes from differentiating an exponential twice. // // Verified numerically against a differenced generator that knows // nothing about the algebra. let (a, b, c, d) = (0.15, -0.5, 0.04, 0.3); let diffusion = Diffusion { drift: |x: f64| a + b * x, volatility: |x: f64| (c + d * x).max(0.0).sqrt(), }; let (grid, dx) = uniform_grid(0.2, 2.2, 2001); for psi in [-0.7, -0.2, 0.4, 1.1] { let u: Vec<f64> = grid.iter().map(|&x| (psi * x).exp()).collect(); let lu = diffusion.generator(&grid, &u, dx); // Away from the boundaries, where the differencing is one-sided. for i in (200..grid.len() - 200).step_by(97) { let x = grid[i]; let ratio = lu[i] / u[i]; let predicted = (a + b * x) * psi + 0.5 * (c + d * x) * psi * psi; assert!( (ratio - predicted).abs() < 1e-6 * predicted.abs().max(1.0), "psi={psi}, x={x}: L u / u = {ratio}, expected {predicted}" ); // And it really is affine in x, which is what lets the two // coefficients be read off separately. let constant = a * psi + 0.5 * c * psi * psi; let slope = b * psi + 0.5 * d * psi * psi; assert!((predicted - (constant + slope * x)).abs() < 1e-12); } } }} /// The stationary density of a one-dimensional diffusion, on a grid.////// A stationary density has no net flux, and in one dimension that forces the/// flux to vanish pointwise rather than merely to have zero divergence --- there/// is nowhere for a current to circulate to. Setting `mu p = (v p)' / 2` and/// solving,////// ```text/// p(x) proportional to exp( integral 2 mu / v ) / v(x),/// ```////// normalised here to integrate to one over the grid supplied.pub fn stationary_density( drift: impl Fn(f64) -> f64, variance: impl Fn(f64) -> f64, grid: &[f64], dx: f64,) -> Vec<f64> { let mut exponent = vec![0.0; grid.len()]; for i in 1..grid.len() { let mid = 0.5 * (grid[i - 1] + grid[i]); exponent[i] = exponent[i - 1] + 2.0 * drift(mid) / variance(mid) * dx; } let peak = exponent.iter().cloned().fold(f64::NEG_INFINITY, f64::max); let mut density: Vec<f64> = grid .iter() .zip(&exponent) .map(|(&x, &e)| (e - peak).exp() / variance(x)) .collect(); let mass: f64 = density.iter().sum::<f64>() * dx; for p in &mut density { *p /= mass; } density} /// The inner product of `L^2(pi)`: `integral f g pi`.////// The weighting is the whole point. A function that is enormous where the/// process never goes is not a large function as far as the process is/// concerned, and it is in this space, not in the flat one, that a reversible/// generator is self-adjoint.pub fn weighted_pairing(f: &[f64], g: &[f64], density: &[f64], dx: f64) -> f64 { f.iter() .zip(g) .zip(density) .map(|((f, g), p)| f * g * p) .sum::<f64>() * dx} #[cfg(test)]mod reversibility_tests { use super::*; /// Trim the ends, where the differenced generator is one-sided. fn interior(v: &[f64]) -> &[f64] { &v[3..v.len() - 3] } fn is_self_adjoint( drift: impl Fn(f64) -> f64 + Copy, variance: impl Fn(f64) -> f64 + Copy, lo: f64, hi: f64, ) -> (f64, f64) { let (grid, dx) = uniform_grid(lo, hi, 4001); let density = stationary_density(drift, variance, &grid, dx); let diffusion = Diffusion { drift, volatility: |x| variance(x).sqrt() }; // Two unrelated test functions. let f: Vec<f64> = grid.iter().map(|&x| (0.7 * x).sin() + 0.3 * x).collect(); let g: Vec<f64> = grid.iter().map(|&x| (-0.2 * x * x).exp() + 0.1 * x).collect(); let lf = diffusion.generator(&grid, &f, dx); let lg = diffusion.generator(&grid, &g, dx); ( weighted_pairing(interior(&lf), interior(&g), interior(&density), dx), weighted_pairing(interior(&f), interior(&lg), interior(&density), dx), ) } #[test] fn every_scalar_diffusion_is_reversible() { // The claim that the spectral method is essentially free in one // dimension. Not just for Ornstein-Uhlenbeck: for anything, because the // zero-flux condition is pointwise and turns the generator into a // divergence, which is symmetric against its own stationary density. let cases: [(&str, fn(f64) -> f64, fn(f64) -> f64, f64, f64); 3] = [ ("Ornstein-Uhlenbeck", |x| -0.8 * x, |_| 0.35 * 0.35, -3.0, 3.0), ("a double well", |x| 2.0 * x - x * x * x, |_| 0.6, -3.0, 3.0), ("state-dependent volatility", |x| -0.5 * x, |x| 0.2 + 0.05 * x * x, -3.0, 3.0), ]; for (name, drift, variance, lo, hi) in cases { let (left, right) = is_self_adjoint(drift, variance, lo, hi); assert!( (left - right).abs() < 1e-4 * left.abs().max(1.0), "{name}: <Lf,g> = {left:.9} but <f,Lg> = {right:.9}" ); } } #[test] fn a_swirl_in_two_dimensions_is_not_reversible() { // And the boundary. Take dZ = -A Z dt + dW in the plane. On linear // functions f(x) = alpha . x the generator acts as -A^T, and the // stationary covariance solves A S + S A^T = I, so // // <Lf, g>_pi = -alpha^T A S beta, <f, Lg>_pi = -alpha^T S A^T beta, // // which agree for every alpha and beta only when A S is symmetric. // // A = I + theta J with J the rotation generator leaves S = I/2 --- the // same isotropic Gaussian a reversible model would have --- while the // process circulates. Same stationary picture, distinguishable film. let theta = 0.9; let a = [[1.0, theta], [-theta, 1.0]]; let s = [[0.5, 0.0], [0.0, 0.5]]; // Confirm S really is stationary: A S + S A^T = I. for i in 0..2 { for j in 0..2 { let lhs: f64 = (0..2).map(|k| a[i][k] * s[k][j] + s[i][k] * a[j][k]).sum(); let identity = if i == j { 1.0 } else { 0.0 }; assert!((lhs - identity).abs() < 1e-12, "A S + S A^T is not the identity"); } } // Now the pairing, for a specific pair of linear test functions. let (alpha, beta) = ([1.0, 0.0], [0.0, 1.0]); let bilinear = |m: [[f64; 2]; 2]| -> f64 { (0..2).map(|i| (0..2).map(|j| alpha[i] * m[i][j] * beta[j]).sum::<f64>()).sum() }; let product = |x: [[f64; 2]; 2], y: [[f64; 2]; 2]| -> [[f64; 2]; 2] { let mut out = [[0.0; 2]; 2]; for i in 0..2 { for j in 0..2 { out[i][j] = (0..2).map(|k| x[i][k] * y[k][j]).sum(); } } out }; let transpose = |x: [[f64; 2]; 2]| [[x[0][0], x[1][0]], [x[0][1], x[1][1]]]; let left = -bilinear(product(a, s)); let right = -bilinear(product(s, transpose(a))); assert!( (left - right).abs() > 0.1, "the swirl should break self-adjointness: {left} against {right}" ); // Symmetric A, same stationary covariance, and the two agree again. let symmetric = [[1.0, 0.0], [0.0, 1.0]]; let l2 = -bilinear(product(symmetric, s)); let r2 = -bilinear(product(s, transpose(symmetric))); assert!((l2 - r2).abs() < 1e-12, "a symmetric drift is reversible"); // And the swirl puts the spectrum off the real line, so there is no // real orthogonal expansion to be had: the eigenvalues of -A are // -1 -/+ i theta. let trace = -(a[0][0] + a[1][1]); let det = a[0][0] * a[1][1] - a[0][1] * a[1][0]; let discriminant = trace * trace - 4.0 * det; assert!(discriminant < 0.0, "expected a complex pair, discriminant {discriminant}"); }} #[cfg(test)]mod domain_tests { use super::*; /// `E[f(W_t)]` by quadrature against the normal density. /// /// `half_width` has to be set from the width of the *product*, not of the /// density: `exp(x^2)` against a narrow Gaussian leaves an integrand far /// wider than the Gaussian, and truncating at a few of the density's /// standard deviations quietly loses a part of it. fn expectation(f: impl Fn(f64) -> f64, t: f64, half_width: f64) -> f64 { let sd = t.sqrt(); let (grid, dx) = uniform_grid(-half_width, half_width, 400_001); grid.iter() .map(|&x| f(x) * crate::black::norm_pdf(x / sd) / sd) .sum::<f64>() * dx } #[test] fn a_kink_is_not_in_the_domain() { // The generator chapter's first illustration of why a domain is needed. // For f(x) = |x| the defining limit at the origin is // // (E|W_t| - 0) / t = sqrt(2 / (pi t)), // // which does not converge -- it grows without bound as t shrinks. So a // function with a kink has no generator at the kink, however harmless it // looks, and the formal expression mu f' + sigma^2 f'' / 2 is not merely // hard to evaluate there but meaningless. for t in [1e-1, 1e-2, 1e-3, 1e-4] { let quadrature = expectation(f64::abs, t, 12.0 * t.sqrt()); let closed_form = (2.0 * t / std::f64::consts::PI).sqrt(); assert!( (quadrature - closed_form).abs() < 1e-6 * closed_form, "t={t}: quadrature {quadrature} against sqrt(2t/pi) {closed_form}" ); let difference_quotient = closed_form / t; let predicted = (2.0 / (std::f64::consts::PI * t)).sqrt(); assert!((difference_quotient - predicted).abs() < 1e-9 * predicted); } // And it diverges: each factor of a hundred in t multiplies the quotient // by ten. let quotient = |t: f64| (2.0 / (std::f64::consts::PI * t)).sqrt(); assert!((quotient(1e-4) / quotient(1e-2) - 10.0).abs() < 1e-9); } #[test] fn fast_growth_leaves_the_semigroup_undefined() { // The second illustration, and a different failure. Here f is perfectly // smooth -- f(x) = exp(x^2) -- so the formal generator is fine at every // point. What fails is the expectation itself: // // E[exp(W_t^2)] = 1 / sqrt(1 - 2t) for t < 1/2, // // and infinite beyond. So P_t f does not exist for t >= 1/2, and the // domain has to exclude growth as well as roughness. for t in [0.05f64, 0.15, 0.3, 0.45] { // Computed with the exponents combined. Evaluating exp(x^2) and // the density separately overflows one and underflows the other, // and their product comes back as a NaN; the integrand itself is // perfectly well behaved. let sd = t.sqrt(); let width = 14.0 * (t / (1.0 - 2.0 * t)).sqrt(); let (grid, dx) = uniform_grid(-width, width, 400_001); let quadrature: f64 = grid .iter() .map(|&x| { let exponent = x * x - 0.5 * x * x / t; exponent.exp() / (sd * (2.0 * std::f64::consts::PI).sqrt()) }) .sum::<f64>() * dx; let closed_form = 1.0 / (1.0 - 2.0 * t).sqrt(); assert!( (quadrature - closed_form).abs() < 1e-5 * closed_form, "t={t}: quadrature {quadrature} against 1/sqrt(1-2t) {closed_form}" ); } // The blow-up is at t = 1/2 exactly, and it is not gradual. let value = |t: f64| 1.0 / (1.0 - 2.0 * t).sqrt(); assert!(value(0.49) < 8.0); assert!(value(0.4999) > 70.0); assert!(!(1.0f64 - 2.0 * 0.5).sqrt().is_normal(), "at t = 1/2 there is nothing left"); }}