quant/src/transform.rs
Prices from a characteristic function, by Fourier inversion.
//! Prices from a characteristic function, by Fourier inversion.//!//! The solvable models chapter produces a transform: two Riccati equations give//! `E[exp(i w ln S_T)]` in closed form, one `w` at a time. This turns that into//! a price, which is the step that makes the whole affine construction useful//! rather than merely elegant.//!//! The route is Carr and Madan's. A call price is not integrable in log-strike//! --- it tends to the spot as the strike goes to zero --- so it is damped by//! `exp(alpha k)` first, and the damping is what forces a condition on which//! moments of the underlying exist. That condition is the same strip the//! chapter's existence theorem needs, arriving from a different direction. /// A complex number, to the extent this module needs one.////// Only multiplication, division and the exponential of an imaginary argument/// are required, so a pair carries it without pulling in a dependency.#[derive(Clone, Copy, Debug)]pub struct C { pub re: f64, pub im: f64,} impl C { pub fn new(re: f64, im: f64) -> Self { C { re, im } } /// `exp(i theta)`. pub fn cis(theta: f64) -> Self { C { re: theta.cos(), im: theta.sin() } } /// `exp(z)`, needed because a characteristic function's exponent is complex. pub fn exp(self) -> Self { let r = self.re.exp(); C { re: r * self.im.cos(), im: r * self.im.sin() } } pub fn times(self, other: C) -> Self { C { re: self.re * other.re - self.im * other.im, im: self.re * other.im + self.im * other.re, } } pub fn over(self, other: C) -> Self { let d = other.re * other.re + other.im * other.im; C { re: (self.re * other.re + self.im * other.im) / d, im: (self.im * other.re - self.re * other.im) / d, } } pub fn scale(self, k: f64) -> Self { C { re: self.re * k, im: self.im * k } }} /// The damped call transform,////// ```text/// zeta(w) = phi(w - i(alpha+1)) / ((alpha + i w)(alpha + 1 + i w)),/// ```////// where `phi` is the characteristic function of `ln S_T`. Derived in the/// solvable models chapter by damping the call, transforming in log-strike and/// exchanging the order of integration.////// `phi` is passed as a function of a *complex* argument, because the numerator/// evaluates it off the real line --- which is exactly why the moment condition/// on `alpha` appears.pub fn damped_transform(phi: impl Fn(C) -> C, w: f64, alpha: f64) -> C { let numerator = phi(C::new(w, -(alpha + 1.0))); let a = C::new(alpha, w); let b = C::new(alpha + 1.0, w); numerator.over(a.times(b))} /// A call price by Fourier inversion,////// ```text/// C(k) = exp(-alpha k) / pi * integral_0^infinity Re[ exp(-i w k) zeta(w) ] dw,/// ```////// evaluated by the trapezoid rule out to `upper` with `points` nodes. `k` is the/// log-strike and the result carries no discounting, so multiply by the discount/// factor for a rate other than zero.////// The integrand oscillates like `exp(-i w k)` and decays like the/// characteristic function, so the truncation has to be chosen against the/// volatility rather than fixed: a low-volatility option needs a longer/// integral, since its transform decays more slowly.pub fn call_by_inversion( phi: impl Fn(C) -> C, log_strike: f64, alpha: f64, upper: f64, points: usize,) -> f64 { let dw = upper / points as f64; let mut total = 0.0; for i in 0..=points { let w = i as f64 * dw; let integrand = C::cis(-w * log_strike).times(damped_transform(&phi, w, alpha)).re; // Trapezoid: half weight at each end. The w = 0 end is finite because // the damping keeps the denominator away from zero. let weight = if i == 0 || i == points { 0.5 } else { 1.0 }; total += weight * integrand * dw; } (-alpha * log_strike).exp() * total / std::f64::consts::PI} /// The characteristic function of `ln S_T` under Black-Scholes with zero rates,/// as a function of a complex argument.////// `ln S_T` is normal with mean `ln S_0 - sigma^2 T / 2` and variance/// `sigma^2 T`, so `phi(u) = exp(i u m - u^2 v / 2)`. Written for complex `u`/// because the inversion needs it off the real line.pub fn lognormal_characteristic(spot: f64, vol: f64, expiry: f64, u: C) -> C { let mean = spot.ln() - 0.5 * vol * vol * expiry; let variance = vol * vol * expiry; // i u m - u^2 v / 2 let iu = C::new(-u.im, u.re); let u_squared = u.times(u); iu.scale(mean).times(C::new(1.0, 0.0)).exp_of_sum(u_squared.scale(-0.5 * variance))} impl C { /// `exp(self + other)`, which is how the characteristic functions here are /// assembled: a linear term plus a quadratic one. fn exp_of_sum(self, other: C) -> C { C { re: self.re + other.re, im: self.im + other.im }.exp() }} #[cfg(test)]mod tests { use super::*; use crate::black::{black_scholes, Side}; const SPOT: f64 = 100.0; #[test] fn the_inversion_reproduces_the_black_scholes_price() { // The whole point: a characteristic function and one integral give the // price a formula gives. Nothing here knows the Black-Scholes formula. for vol in [0.15, 0.25, 0.4] { for expiry in [0.25, 1.0, 3.0] { let phi = |u: C| lognormal_characteristic(SPOT, vol, expiry, u); // The transform decays like exp(-w^2 sigma^2 T / 2), so the // truncation is set from that scale rather than fixed. let upper = 20.0 / (vol * expiry.sqrt()); for strike in [70.0f64, 90.0, 100.0, 115.0, 140.0] { let inverted = call_by_inversion(&phi, strike.ln(), 1.5, upper, 20_000); let exact = black_scholes(SPOT, strike, vol, expiry, 0.0, 0.0, Side::Call); assert!( (inverted - exact).abs() < 1e-6 * exact.max(1.0), "vol={vol} T={expiry} K={strike}: inverted {inverted:.10} \ against exact {exact:.10}" ); } } } } #[test] fn the_answer_does_not_depend_on_the_damping() { // alpha is a device, not a parameter of the problem, so the price must // not move with it. A lognormal has moments of every order, so every // positive alpha is admissible here and the insensitivity is complete. let (vol, expiry, strike) = (0.25f64, 1.0f64, 110.0f64); let phi = |u: C| lognormal_characteristic(SPOT, vol, expiry, u); let exact = black_scholes(SPOT, strike, vol, expiry, 0.0, 0.0, Side::Call); for alpha in [0.5, 1.0, 1.5, 3.0, 6.0] { let price = call_by_inversion(&phi, strike.ln(), alpha, 80.0, 20_000); assert!( (price - exact).abs() < 1e-6 * exact, "alpha={alpha} gave {price:.10} against {exact:.10}" ); } } #[test] fn the_damping_is_what_makes_the_integral_exist() { // Why alpha is needed at all. At alpha = 0 the denominator // (alpha + i w)(alpha + 1 + i w) vanishes at w = 0, so the integrand has // a pole at the origin -- the undamped call price has no Fourier // transform, which is the reason for the whole construction. let phi = |u: C| lognormal_characteristic(SPOT, 0.25, 1.0, u); let size = |w: f64, alpha: f64| { let z = damped_transform(&phi, w, alpha); (z.re * z.re + z.im * z.im).sqrt() }; // The signature of a pole is how the size scales as w approaches zero, // not how large it is -- the numerator carries a factor S^(1+alpha) and // so is large whatever happens. Undamped, the denominator is i w near // the origin, so shrinking w by a thousand grows the transform by a // thousand. let undamped = size(1e-6, 0.0) / size(1e-3, 0.0); assert!( (undamped / 1000.0 - 1.0).abs() < 0.01, "expected a simple pole, size grew by {undamped:.1}" ); // Damped, the denominator tends to alpha(alpha+1) and the ratio is one: // no pole, and the integral through w = 0 is finite. for alpha in [0.5, 1.5, 3.0] { let ratio = size(1e-6, alpha) / size(1e-3, alpha); assert!( (ratio - 1.0).abs() < 1e-3, "alpha={alpha} should remove the pole, size grew by {ratio}" ); } } #[test] fn a_damping_outside_the_moment_strip_returns_a_number_the_model_lacks() { // The condition on alpha, made visible. The numerator evaluates the // characteristic function at w - i(alpha+1), which is the moment // E[S_T^{1+alpha}]. A lognormal has all of them, so to see the failure // the moment has to be capped by hand: cut the characteristic function // off outside a strip, as a model with a moment explosion does, and the // inversion returns nonsense rather than an error. let (vol, expiry, strike) = (0.25f64, 1.0f64, 110.0f64); let limit = 2.0; // moments exist only up to order 2 let capped = |u: C| { if -u.im > limit { C::new(f64::NAN, f64::NAN) } else { lognormal_characteristic(SPOT, vol, expiry, u) } }; // alpha + 1 = 1.5 is inside the strip and works. let inside = call_by_inversion(&capped, strike.ln(), 0.5, 80.0, 20_000); let exact = black_scholes(SPOT, strike, vol, expiry, 0.0, 0.0, Side::Call); assert!((inside - exact).abs() < 1e-6 * exact, "alpha=0.5 should be admissible"); // alpha + 1 = 2.5 is outside it, and the answer is not a price. let outside = call_by_inversion(&capped, strike.ln(), 1.5, 80.0, 20_000); assert!(!outside.is_finite(), "outside the strip should not return a price"); }}