quant/src/heston.rs
Heston, simulated, to separate two things a smile can do.
//! Heston, simulated, to separate two things a smile can do.//!//! The smile dynamics chapter needs a fair comparison. A local volatility model//! can be made to reproduce the market's *maturity* term structure of skew, but//! only by acquiring a dependence on calendar time, which then destroys its//! forward smile. The question is whether that trade is forced or whether it is//! a consequence of having one factor.//!//! Heston answers it. Its parameters carry no calendar time at all --- the//! model is time homogeneous, so it looks the same from every start date --- and//! yet its at-the-money skew still decays with expiry. The decay comes from the//! dynamics over the option's life rather than from a dated parameter, and that//! is the distinction the chapter turns on.//!//! Everything here is by simulation, including the spot smile, so that the spot//! and forward numbers are produced by the same code and can be compared without//! a scheme difference standing between them. use crate::black::{implied_vol_black76, Side};use crate::pathwise::Rng; /// The Heston model,////// ```text/// dS = sqrt(v) S dW,/// dv = kappa (theta - v) dt + eta sqrt(v) dZ, dW dZ = rho dt./// ```////// Note what is *not* in this list: any function of `t`. Every parameter is a/// constant, which is what makes the model time homogeneous.pub struct Heston { pub v0: f64, pub kappa: f64, pub theta: f64, pub eta: f64, pub rho: f64, pub spot: f64,} impl Heston { /// The implied volatilities of options on `S(start + tenor) / S(start)`, at /// the given moneyness. /// /// `start = 0` gives the ordinary spot smile at expiry `tenor`, which is why /// one routine serves both halves of the comparison. /// /// The variance uses full truncation --- the drift sees `max(v, 0)` and so /// does the diffusion --- which is the standard way to keep an Euler scheme /// honest when the variance can be pushed negative by a discrete step. pub fn forward_smile( &self, start: f64, tenor: f64, moneyness: &[f64], paths: usize, steps_per_year: usize, seed: u64, ) -> Vec<Option<f64>> { let mut rng = Rng::new(seed); let horizon = start + tenor; let total_steps = ((horizon * steps_per_year as f64).ceil() as usize).max(2); let dt = horizon / total_steps as f64; let split = ((start / horizon) * total_steps as f64).round() as usize; let root_dt = dt.sqrt(); let root_one_minus = (1.0 - self.rho * self.rho).sqrt(); let mut payoffs = vec![0.0; moneyness.len()]; for _ in 0..paths { let mut s = self.spot; let mut v = self.v0; let mut at_start = self.spot; for step in 0..total_steps { if step == split { at_start = s; } let z1 = rng.next_normal(); let z2 = rng.next_normal(); // The spot's Brownian motion and the variance's are correlated; // build the second from the first so rho enters exactly once. let dw = z1; let dz = self.rho * z1 + root_one_minus * z2; let v_plus = v.max(0.0); let root_v = v_plus.sqrt(); s *= (-0.5 * v_plus * dt + root_v * root_dt * dw).exp(); v += self.kappa * (self.theta - v_plus) * dt + self.eta * root_v * root_dt * dz; } let ratio = s / at_start; for (i, &k) in moneyness.iter().enumerate() { payoffs[i] += (ratio - k).max(0.0); } } let n = paths as f64; moneyness .iter() .zip(&payoffs) .map(|(&k, &total)| { // The ratio is a martingale here, so its forward is one. implied_vol_black76(total / n, 1.0, k, tenor, Side::Call) }) .collect() } /// The at-the-money slope of that smile in log-moneyness. /// /// With `start = 0` this is the ordinary skew at expiry `tenor`, and /// sweeping `tenor` traces the maturity term structure. With `tenor` held /// fixed and `start` swept, it traces the forward skew. pub fn skew( &self, start: f64, tenor: f64, paths: usize, steps_per_year: usize, seed: u64, ) -> Option<f64> { let h = 0.05f64; let moneyness = [(-h).exp(), 1.0, h.exp()]; let vols = self.forward_smile(start, tenor, &moneyness, paths, steps_per_year, seed); Some((vols[2]? - vols[0]?) / (2.0 * h)) }} #[cfg(test)]mod tests { use super::*; /// Equity shaped, and deliberately started at the long-run variance so the /// variance process is already stationary. Feller is satisfied /// (`2 kappa theta = 0.16` against `eta^2 = 0.09`), so the truncation is a /// safeguard rather than something the results lean on. fn model() -> Heston { Heston { v0: 0.04, kappa: 2.0, theta: 0.04, eta: 0.3, rho: -0.7, spot: 100.0 } } const PATHS: usize = 400_000; const SEED: u64 = 20260807; #[test] fn the_skew_decays_with_expiry() { // Half of the comparison: constant parameters, and still a maturity // term structure. The decay comes from the variance mean reverting over // the option's life, not from anything dated. let m = model(); let mut previous = f64::NEG_INFINITY; let mut seen = Vec::new(); for t in [0.25, 0.5, 1.0, 2.0, 5.0] { let s = m.skew(0.0, t, PATHS, 200, SEED).unwrap(); assert!(s < 0.0, "negative correlation should skew down, got {s} at {t}"); assert!(s > previous, "skew did not shrink at {t}: {s} after {previous}"); previous = s; seen.push(s); } // And it is a real decay, not a drift: five years keeps under half of // what three months has. let ratio = seen[4] / seen[0]; assert!(ratio < 0.5, "five-year skew is {ratio:.2} of the three-month one"); } #[test] fn the_forward_skew_does_not_decay() { // The other half. Same model, same code, one-year tenor throughout, and // the start date pushed out to four years. let m = model(); let spot = m.skew(0.0, 1.0, PATHS, 200, SEED).unwrap(); for start in [1.0, 2.0, 4.0] { let fwd = m.skew(start, 1.0, PATHS, 200, SEED).unwrap(); let ratio = fwd / spot; assert!( (ratio - 1.0).abs() < 0.1, "forward skew starting at {start}y is {ratio:.3} of today's one-year skew" ); } } #[test] fn stationarity_is_what_holds_the_forward_skew_up() { // The mechanism, isolated. Time homogeneity alone does not pin the // forward smile; it pins it once the variance has reached its stationary // distribution. Start the variance away from theta and the forward skew // has to travel before it settles, which is visible as a drift in the // first year and its absence later. let hot = Heston { v0: 0.09, ..model() }; let near = hot.skew(0.0, 1.0, PATHS, 200, SEED).unwrap(); let mid = hot.skew(1.0, 1.0, PATHS, 200, SEED).unwrap(); let far = hot.skew(4.0, 1.0, PATHS, 200, SEED).unwrap(); // Starting hot, today's one-year skew is diluted by a high variance that // has not yet reverted, so it is flatter than the settled one. assert!(near > mid, "expected the hot start to flatten today's skew: {near} vs {mid}"); // Once reverted, successive start dates agree: 1/kappa is half a year, // so by one year the variance is already close to stationary. assert!( (far / mid - 1.0).abs() < 0.1, "four years should look like one: {far} against {mid}" ); }} /// The variance risk premium, as the two measures' disagreement about variance.////// The fitting chapter asks whether implied volatility exceeding realised is a/// risk premium or a flow effect. The risk-premium half of the answer is a/// statement about Girsanov, and it can be written down exactly.////// Take the standard specification in which the market price of variance risk is/// proportional to the volatility, so that the change of measure adds/// `-lambda v dt` to the variance's drift. Then////// ```text/// dv = kappa (theta - v) dt + eta sqrt(v) dZ under P/// dv = kappa* (theta* - v) dt + eta sqrt(v) dZ* under Q/// ```////// with////// ```text/// kappa* = kappa + lambda eta, theta* = kappa theta / kappa*./// ```////// Three features of that pair are the content.////// * `eta` is **the same under both measures**. It is a diffusion coefficient,/// and a change of measure moves drifts only. So vol-of-vol is not where the/// premium lives, and it is jointly identified by prices and by history --- an/// over-identifying restriction that is testable and rarely tested./// * The premium sits entirely in the drift, as a lower mean reversion and a/// higher long-run level when `lambda < 0`. Which is the same formal position/// the equity risk premium occupies: a drift that differs between measures./// * Its size is exactly computable, because the expected integrated variance of/// a square-root process has a closed form.pub struct VariancePremium { /// Expected annualised volatility over the horizon under the historical /// measure --- what a realised volatility estimate converges to. pub realised: f64, /// The same under the pricing measure --- what an at-the-money implied /// volatility of that maturity reflects. pub implied: f64, /// The gap, in volatility points. pub premium: f64, /// The risk-neutral parameters, for inspection. pub kappa_q: f64, pub theta_q: f64,} impl Heston { /// Expected integrated variance over `[0, T]`, divided by `T` and rooted, so /// the answer is in the units an implied volatility is quoted in. /// /// For `dv = kappa (theta - v) dt` the mean solves an ordinary differential /// equation and the integral is exact: /// /// ```text /// E[ integral_0^T v ] = theta T + (v0 - theta) (1 - exp(-kappa T)) / kappa. /// ``` pub fn expected_volatility(&self, horizon: f64) -> f64 { let integrated = self.theta * horizon + (self.v0 - self.theta) * (1.0 - (-self.kappa * horizon).exp()) / self.kappa; (integrated / horizon).sqrt() } /// The premium implied by a market price of variance risk `lambda`. /// /// `lambda` negative is the empirically relevant sign: long variance is a /// hedge against bad states, so bearing it earns a negative expected return /// and paying for it costs a premium. let kappa_q = self.kappa + lambda * self.eta; let theta_q = self.kappa * self.theta / kappa_q; let risk_neutral = Heston { kappa: kappa_q, theta: theta_q, ..*self }; let realised = self.expected_volatility(horizon); let implied = risk_neutral.expected_volatility(horizon); VariancePremium { realised, implied, premium: implied - realised, kappa_q, theta_q, } }} #[cfg(test)]mod premium_tests { use super::*; fn historical() -> Heston { // A variance starting at its own long-run level, so the premium below is // not confused with the transient from a state away from equilibrium. Heston { v0: 0.04, kappa: 2.0, theta: 0.04, eta: 0.3, rho: -0.7, spot: 100.0, } } #[test] // The structural claim. Changing measure moves kappa and theta and // leaves eta alone, so a vol-of-vol estimated from the history of // realised volatility and one implied by the smile's convexity are // estimates of the same number. let h = historical(); let p = h.variance_premium(-1.0, 1.0); assert!(p.kappa_q < h.kappa, "a negative lambda should slow the reversion"); assert!(p.theta_q > h.theta, "and raise the long-run level"); // eta is untouched by construction; assert it explicitly so that a // future edit adding a premium to it fails here. let risk_neutral = Heston { kappa: p.kappa_q, theta: p.theta_q, ..h }; assert_eq!(risk_neutral.eta, h.eta, "vol-of-vol must be measure-invariant"); } #[test] // With the variance at equilibrium, expected volatility under P is // sqrt(theta) = 20%. A market price of variance risk of -1 lifts the // risk-neutral long-run variance and the one year at-the-money implied // volatility with it. The chapter quotes these. let h = historical(); let p = h.variance_premium(-1.0, 1.0); assert!( (p.realised - 0.20).abs() < 1e-12, "at equilibrium the realised expectation is sqrt(theta): {:.6}", p.realised ); // 0.90 volatility points, on a 20% realised level. Which is the // right order for an equity index, where the observed gap is one to // three points -- so a market price of variance risk of a few units is // what the market is telling us, not a fitted fudge. assert!( (p.premium - 0.008962).abs() < 1e-5, "premium {:.4} volatility points", p.premium * 100.0 ); // And it scales with the price of risk, so the observed gap measures // lambda rather than merely being consistent with one. let bigger = h.variance_premium(-2.0, 1.0); assert!(bigger.premium > p.premium * 1.8, "the premium should roughly double"); } #[test] fn a_longer_horizon_shows_more_of_it() { // The transient matters: the premium is a change to the long-run level, // so it shows up more fully the longer the horizon over which variance // is averaged. Which is a testable signature -- a risk premium has a // term structure fixed by the mean reversion, and a flow effect has // whatever term structure the flow has. let h = historical(); let short = h.variance_premium(-1.0, 0.25).premium; let long = h.variance_premium(-1.0, 5.0).premium; assert!(long > short, "{long:.5} against {short:.5}"); }}