quant/src/pathwise.rs
Brownian paths, and the two things the Brownian motion and no-arbitrage chapters do with them.
//! Brownian paths, and the two things the Brownian motion and//! no-arbitrage chapters do with them.//!//! Everything here is deterministic given a seed. That is not a convenience: a//! figure in the notes has to show the same picture in the PDF, on the site, and//! next year, and the prose beside it quotes numbers off it. A figure driven by//! an unseeded generator would quietly contradict its own caption. use crate::black::{black76, norm_cdf, Side}; /// A small counter-based generator.////// Not cryptographic and not trying to be. What it must do is produce the same/// stream everywhere from the same seed — including in wasm, which rules out/// anything reaching for the system — and pass well enough for drawing a/// Brownian path. This is `splitmix64`, which is four lines and does both.pub struct Rng(u64); impl Rng { pub fn new(seed: u64) -> Self { Rng(seed) } fn next_u64(&mut self) -> u64 { self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15); let mut z = self.0; z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9); z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB); z ^ (z >> 31) } /// Uniform on (0, 1), open at both ends so a log of it is finite. pub fn next_uniform(&mut self) -> f64 { // 53 bits, the most a double holds exactly. let bits = self.next_u64() >> 11; (bits as f64 + 0.5) / (1u64 << 53) as f64 } /// A standard normal, by Box-Muller. /// /// The second variate of each pair is thrown away. Wasteful, and irrelevant /// here — the cost of a figure is the drawing, not the sampling — and it /// keeps the generator stateless between calls, which keeps the stream /// reproducible under any call pattern. pub fn next_normal(&mut self) -> f64 { let u1 = self.next_uniform(); let u2 = self.next_uniform(); (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos() }} /// A Brownian path on `[0, t]`, sampled at `steps + 1` points including zero.pub fn brownian_path(t: f64, steps: usize, seed: u64) -> Vec<f64> { let dt = t / steps as f64; let sd = dt.sqrt(); let mut rng = Rng::new(seed); let mut w = Vec::with_capacity(steps + 1); let mut current = 0.0; w.push(0.0); for _ in 0..steps { current += sd * rng.next_normal(); w.push(current); } w} /// The three quantities the Brownian motion chapter compares, computed on one path at one/// resolution.#[derive(Clone, Copy, Debug)]pub struct Sums { /// Number of intervals the path was chopped into. pub steps: usize, /// The Ito sum, taking the integrand at the *start* of each interval. pub ito: f64, /// The Stratonovich sum, taking the average of the two endpoints. pub stratonovich: f64, /// The quadratic variation, the sum of squared increments. pub quadratic_variation: f64,} /// Both Riemann sums for the integral of `W` against `dW`, and the quadratic/// variation, over the same path at a given coarseness.////// `w` is a path sampled finely; `stride` says which of its points to use, so/// that every resolution is a coarsening of one path rather than a new draw./// That is the whole point of the comparison: if each resolution used a fresh/// path, the figure would be showing sampling noise rather than the fact that/// the two definitions disagree.pub fn sums_over_path(w: &[f64], t: f64, stride: usize) -> Sums { let steps = (w.len() - 1) / stride; let mut ito = 0.0; let mut strat = 0.0; let mut qv = 0.0; for i in 0..steps { let a = w[i * stride]; let b = w[(i + 1) * stride]; let dw = b - a; ito += a * dw; strat += 0.5 * (a + b) * dw; qv += dw * dw; } let _ = t; Sums { steps, ito, stratonovich: strat, quadratic_variation: qv }} /// How strongly the realised variance in one window predicts the next.////// The price process chapter uses this to settle a question the/// characterization theorems raise but do not answer: a continuous local/// martingale is always an integral against a Brownian motion, so does it/// inherit Brownian motion's independent increments? It does not, and this/// measures the failure.////// The experiment builds `M = integral of sigma dW` on `[0,1]`, splits it in/// half, and correlates the realised variance of the first half against the/// second across many paths. Increments independent of the past would force this/// to be zero, because the two windows are disjoint.////// With `vol_of_vol` at zero the volatility is constant and the correlation is/// zero, as it must be. Turn it up and the correlation is large and positive:/// the quadratic variation is now random and persistent, so a violent first half/// predicts a violent second. Under its market name this is volatility/// clustering, and it is one of the least controversial facts about returns.////// The volatility is driven by its own Brownian motion, independent of the one/// driving `M`, so `M` remains a perfectly good continuous local martingale/// throughout. Nothing has been broken to produce the effect.pub fn variance_clustering( vol_of_vol: f64, paths: usize, steps: usize, seed: u64,) -> f64 { let dt = 1.0 / steps as f64; let sqrt_dt = dt.sqrt(); let mut rng = Rng::new(seed); let (mut sx, mut sy, mut sxx, mut syy, mut sxy) = (0.0, 0.0, 0.0, 0.0, 0.0); for _ in 0..paths { let mut sigma = 0.20f64; let (mut first, mut second) = (0.0, 0.0); for i in 0..steps { // Drawn unconditionally so that the two cases consume the same // random numbers and differ only in what they do with them. let z = rng.next_normal(); sigma *= (-0.5 * vol_of_vol * vol_of_vol * dt + vol_of_vol * sqrt_dt * z).exp(); let increment = sigma * sqrt_dt * rng.next_normal(); if i < steps / 2 { first += increment * increment; } else { second += increment * increment; } } sx += first; sy += second; sxx += first * first; syy += second * second; sxy += first * second; } let n = paths as f64; let (mx, my) = (sx / n, sy / n); let (vx, vy) = (sxx / n - mx * mx, syy / n - my * my); if vx <= 0.0 || vy <= 0.0 { return 0.0; } (sxy / n - mx * my) / (vx.sqrt() * vy.sqrt())} /// Realised quadratic variation of `M = integral of sigma dW`, and of the/// Brownian motion `W` recovered from it, measured across many paths.////// The price process chapter uses this to answer the obvious objection to the/// characterization theorems. Levy's theorem needs a deterministic quadratic/// variation, and a stochastic volatility model plainly has a random one, so/// how is the theory available at all?////// Because Levy's theorem is never applied to the price. It is applied to `W`,/// and `W` is built to satisfy the hypothesis: dividing each increment by/// `sigma` divides the randomness out of the variation, whatever `sigma` was/// doing. The numbers here are that statement measured.////// Returns each variation's mean and its standard deviation across paths as a/// fraction of that mean. The fraction is the quantity to watch, and what/// matters is how it behaves as the grid is refined: shrinking means the/// variation is converging to a constant, holding steady means it is converging/// to a genuinely random limit.////// One honest limit on the experiment. It builds `M` out of `W`, so recovering/// `W` by dividing is exact rather than estimated, and the recovery itself is/// not what is being tested. What is being tested is the variance structure --/// that `<M>` stays random under refinement while `<W>` does not -- and that is/// the part the theorems turn on.pub struct VariationSpread { pub price_mean: f64, pub price_relative_sd: f64, pub brownian_mean: f64, pub brownian_relative_sd: f64,} pub fn variation_spread( vol_of_vol: f64, paths: usize, steps: usize, seed: u64,) -> VariationSpread { let dt = 1.0 / steps as f64; let sqrt_dt = dt.sqrt(); let mut rng = Rng::new(seed); let (mut sm, mut smm, mut sw, mut sww) = (0.0, 0.0, 0.0, 0.0); for _ in 0..paths { let mut sigma = 0.20f64; let (mut price, mut brownian) = (0.0, 0.0); for _ in 0..steps { let z = rng.next_normal(); sigma *= (-0.5 * vol_of_vol * vol_of_vol * dt + vol_of_vol * sqrt_dt * z).exp(); let increment = sigma * sqrt_dt * rng.next_normal(); price += increment * increment; // The increment of the recovered Brownian motion, dM / sigma. let rescaled = increment / sigma; brownian += rescaled * rescaled; } sm += price; smm += price * price; sw += brownian; sww += brownian * brownian; } let n = paths as f64; let (mm, mw) = (sm / n, sw / n); VariationSpread { price_mean: mm, price_relative_sd: (smm / n - mm * mm).max(0.0).sqrt() / mm, brownian_mean: mw, brownian_relative_sd: (sww / n - mw * mw).max(0.0).sqrt() / mw, }} #[cfg(test)]mod tests { use super::*; #[test] // The control. With sigma constant the two windows are functions of // disjoint Brownian increments, so they are genuinely independent and // the measured correlation is sampling noise about zero. let c = variance_clustering(0.0, 100_000, 400, 20_260_804); assert!(c.abs() < 0.02, "correlation was {c}, expected zero"); } #[test] fn stochastic_volatility_makes_the_increments_dependent() { // The price process chapter's point. M is still a continuous local // martingale and still an integral against a Brownian motion -- the // representation theorem does not care -- and yet its increments are // plainly not independent of the past. let c = variance_clustering(1.0, 100_000, 400, 20_260_804); assert!(c > 0.4, "correlation was only {c}"); } #[test] fn any_randomness_in_the_variation_destroys_the_independence() { // It is the randomness of the quadratic variation that does it, and it // takes very little: most of the effect is already there at a vol-of-vol // the market would call low. // // The threshold is loose at the top of the range on purpose. Once the // vol-of-vol is large the realised variances are heavy tailed enough // that the correlation *estimator* is itself high variance, so a tight // bound here would be a flaky test rather than a stronger claim. for vol_of_vol in [0.25, 0.5, 1.0, 1.5, 2.0] { let c = variance_clustering(vol_of_vol, 60_000, 300, 11_223_344); assert!(c > 0.25, "at vol_of_vol={vol_of_vol} the correlation was {c}"); } } #[test] fn the_rescaled_variation_is_deterministic_whatever_the_volatility_does() { // The answer to "Levy needs a deterministic quadratic variation, so how // is any of this available to a stochastic volatility model?". // // The recovered Brownian motion's variation converges to 1 on every // path, in both the constant and the stochastic case, and its spread is // the chi-square sampling spread sqrt(2/steps) and nothing else. That is // the hypothesis of Levy's theorem being manufactured rather than // assumed. for vol_of_vol in [0.0, 1.0] { for steps in [400usize, 1600] { let s = variation_spread(vol_of_vol, 40_000, steps, 20_260_804); assert!( (s.brownian_mean - 1.0).abs() < 0.01, "vol_of_vol={vol_of_vol} steps={steps}: mean {}", s.brownian_mean ); let sampling = (2.0 / steps as f64).sqrt(); assert!( (s.brownian_relative_sd / sampling - 1.0).abs() < 0.05, "vol_of_vol={vol_of_vol} steps={steps}: spread {} against {sampling}", s.brownian_relative_sd ); } } } #[test] fn the_price_variation_stays_random_however_finely_it_is_sampled() { // And the contrast that makes the previous test mean something. With // constant volatility the price's own variation concentrates like any // sample average; with stochastic volatility it does not concentrate at // all, because its limit is the random variable integral of sigma^2. let constant_coarse = variation_spread(0.0, 40_000, 400, 20_260_804); let constant_fine = variation_spread(0.0, 40_000, 1600, 20_260_804); assert!( constant_fine.price_relative_sd < 0.6 * constant_coarse.price_relative_sd, "constant volatility did not concentrate: {} then {}", constant_coarse.price_relative_sd, constant_fine.price_relative_sd ); let random_coarse = variation_spread(1.0, 40_000, 400, 20_260_804); let random_fine = variation_spread(1.0, 40_000, 1600, 20_260_804); assert!( random_fine.price_relative_sd > 0.8 * random_coarse.price_relative_sd, "stochastic volatility concentrated away: {} then {}", random_coarse.price_relative_sd, random_fine.price_relative_sd ); assert!(random_fine.price_relative_sd > 1.0, "{}", random_fine.price_relative_sd); } #[test] fn the_measured_correlation_peaks_and_then_falls_back() { // Recorded because it is a trap. The dependence does not weaken as the // vol-of-vol rises -- it strengthens. What falls is the ability of a // *linear* correlation to see it: the realised variances become more // extremely lognormal, and the dependence chapter shows that the attainable Pearson // correlation between lognormals collapses as their volatility grows. // // So this is the same phenomenon as the dependence chapter's Frechet bounds, met // fourteen chapters early, and it is a good reason not to read a // correlation as a measure of dependence. let low = variance_clustering(0.5, 200_000, 400, 20_260_804); let high = variance_clustering(2.0, 200_000, 400, 20_260_804); assert!(low > high, "expected the measured correlation to fall: {low}, {high}"); assert!(low > 0.6 && high > 0.4, "{low}, {high}"); } const T: f64 = 1.0; const FINE: usize = 1 << 16; #[test] fn the_generator_is_reproducible() { let a = brownian_path(T, 1000, 42); let b = brownian_path(T, 1000, 42); assert_eq!(a, b); assert_ne!(brownian_path(T, 1000, 43), a); } #[test] fn the_increments_look_standard_normal() { let n = 200_000; let w = brownian_path(1.0, n, 7); let dt = 1.0 / n as f64; let d: Vec<f64> = w.windows(2).map(|p| (p[1] - p[0]) / dt.sqrt()).collect(); let mean = d.iter().sum::<f64>() / n as f64; let var = d.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / n as f64; // Three standard errors of the mean is 3/sqrt(n). assert!(mean.abs() < 3.0 / (n as f64).sqrt(), "mean was {mean}"); assert!((var - 1.0).abs() < 0.02, "variance was {var}"); } #[test] fn stratonovich_is_exactly_half_w_squared_at_every_resolution() { // The midpoint sum telescopes: sum of (a+b)/2 * (b-a) is sum of // (b^2 - a^2)/2, which collapses to W_T^2/2 whatever the partition. So // the Stratonovich integral obeys the ordinary chain rule exactly, and // its flatness in the figure is an identity rather than a convergence. let w = brownian_path(T, FINE, 11); let target = 0.5 * w[FINE] * w[FINE]; for k in 0..12 { let s = sums_over_path(&w, T, 1 << k); assert!( (s.stratonovich - target).abs() < 1e-9, "at {} steps it was {}, expected {target}", s.steps, s.stratonovich ); } } #[test] fn the_gap_between_the_two_is_exactly_half_the_quadratic_variation() { // Algebraically: b*(b-a) - a*(b-a) = (b-a)^2, so the midpoint sum // exceeds the left sum by half the sum of squared increments. This is // the identity the whole of Ito calculus is built on, and it holds term // by term, before any limit is taken. let w = brownian_path(T, FINE, 3); for k in 0..12 { let s = sums_over_path(&w, T, 1 << k); let gap = s.stratonovich - s.ito; assert!( (gap - 0.5 * s.quadratic_variation).abs() < 1e-9, "at {} steps the gap was {gap}", s.steps ); } } #[test] fn the_ito_isometry_holds() { // E[(integral of H dW)^2] = E[integral of H^2 ds]. // // Checked at H = W, where both sides are known exactly. The integral is // (W_T^2 - T)/2, so the left side is E[(W_T^2 - T)^2]/4 = T^2/2 using // E[W_T^4] = 3 T^2; and the right side integrates E[W_s^2] = s, which is // also T^2/2. // // The tolerance is computed from the sample rather than chosen. This // estimator is badly behaved — the left side is a fourth moment of a // Gaussian, so its own variance is large and it converges slowly — and a // hand-picked tolerance here is either so tight it fails on a different // seed or so loose it would accept a wrong answer. Four standard errors // is a statement about the estimator, not a guess. const T: f64 = 1.5; const STEPS: usize = 2000; const PATHS: usize = 40_000; let dt = T / STEPS as f64; let mut left = Vec::with_capacity(PATHS); let mut right = Vec::with_capacity(PATHS); for path in 0..PATHS { let w = brownian_path(T, STEPS, 1_000 + path as u64); let ito = sums_over_path(&w, T, 1).ito; left.push(ito * ito); // The Riemann sum the chapter's calculation ends on, width included. right.push(w[..STEPS].iter().map(|x| x * x * dt).sum::<f64>()); } let stats = |v: &[f64]| { let mean = v.iter().sum::<f64>() / v.len() as f64; let var = v.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / v.len() as f64; (mean, (var / v.len() as f64).sqrt()) }; let (lhs, se_lhs) = stats(&left); let (rhs, se_rhs) = stats(&right); let exact = T * T / 2.0; assert!( (lhs - exact).abs() < 4.0 * se_lhs, "left side {lhs} against {exact}, standard error {se_lhs}" ); assert!( (rhs - exact).abs() < 4.0 * se_rhs, "right side {rhs} against {exact}, standard error {se_rhs}" ); assert!( (lhs - rhs).abs() < 4.0 * (se_lhs + se_rhs), "the two sides differ: {lhs} against {rhs}" ); } #[test] fn the_quadratic_variation_converges_to_the_elapsed_time() { // And so the Ito integral converges to W_T^2/2 - T/2, the extra term // that ordinary calculus does not have. let w = brownian_path(T, FINE, 5); let coarse = sums_over_path(&w, T, 1 << 10); let fine = sums_over_path(&w, T, 1); assert!( (fine.quadratic_variation - T).abs() < (coarse.quadratic_variation - T).abs(), "refining did not improve the quadratic variation" ); assert!( (fine.quadratic_variation - T).abs() < 0.02, "at the finest resolution it was {}", fine.quadratic_variation ); let expected = 0.5 * w[FINE] * w[FINE] - 0.5 * T; assert!((fine.ito - expected).abs() < 0.02, "Ito sum was {}", fine.ito); }} // ---------------------------------------------------------------------------// Delta hedging.// --------------------------------------------------------------------------- /// The outcome of hedging one short call to expiry.#[derive(Clone, Copy, Debug)]pub struct HedgeResult { /// Profit and loss at expiry: the premium taken in, plus everything the /// hedge made, less what the option cost to settle. pub pnl: f64, /// The same position left unhedged, for comparison. pub unhedged_pnl: f64,} /// Black-Scholes delta of a call, with zero rates.////// Zero rates throughout this module. It costs nothing in generality — the/// forward-measure version of the Black-Scholes chapter says the interesting/// content is in the forward, not the discounting — and it removes the/// financing terms that would otherwise clutter every line of the hedge/// accounting without changing what the figure shows.pub fn call_delta(s: f64, k: f64, sigma: f64, t: f64) -> f64 { if t <= 0.0 { return if s > k { 1.0 } else { 0.0 }; } let v = sigma * t.sqrt(); norm_cdf(((s / k).ln() + 0.5 * v * v) / v)} /// Black-Scholes gamma of a call, with zero rates. See [`call_delta`].pub fn call_gamma(s: f64, k: f64, sigma: f64, t: f64) -> f64 { if t <= 0.0 || s <= 0.0 { return 0.0; } let v = sigma * t.sqrt(); let d1 = ((s / k).ln() + 0.5 * v * v) / v; crate::black::norm_pdf(d1) / (s * v)} /// What selling a delta-hedged option earns, and what the earning is made of.pub struct VolTrade { /// Premium taken in, plus the hedge, less the settlement. pub pnl: f64, /// The same number predicted a completely different way: the gamma weighted /// difference between the variance sold and the variance that arrived, /// accumulated along the same path. See [`sell_variance`]. pub gamma_weighted: f64, /// Realised variance over the path, annualised, for comparison with the /// volatility the option was sold at. pub realised_variance: f64,} /// Sell an option at one volatility into a world that realises another.////// This is the volatility relative value trade of the alpha chapter: the option/// is sold at `implied` and hedged with a delta computed at `implied`, because/// that is what a desk actually does, while the path arrives at `realised`.////// The point of the function is not the profit but its decomposition. Selling/// the option and hedging it is not a bet on realised variance; it is a bet on/// gamma weighted realised variance, and the two differ whenever the moves do/// not arrive uniformly. `gamma_weighted` accumulates////// ```text/// (1/2) Gamma S^2 (implied^2 - realised^2) dt/// ```////// along the path, which is a different computation from the hedge accounting/// in every respect --- it never forms the payoff, and it uses the second/// derivative where the hedge uses the first. That the two agree is the content.pub fn sell_variance( s0: f64, k: f64, implied: f64, realised: f64, t: f64, steps: usize, seed: u64,) -> VolTrade { let premium = black76(s0, k, implied, t, Side::Call); let dt = t / steps as f64; let mut rng = Rng::new(seed); let mut s = s0; let mut hedge_pnl = 0.0; let mut gamma_weighted = 0.0; let mut sum_squared_returns = 0.0; for i in 0..steps { let remaining = t - i as f64 * dt; // Both Greeks at the volatility the position is marked at, which is the // implied one. Marking at the realised volatility would be assuming the // answer the trade is trying to earn. let delta = call_delta(s, k, implied, remaining); let gamma = call_gamma(s, k, implied, remaining); gamma_weighted += 0.5 * gamma * s * s * (implied * implied - realised * realised) * dt; let z = rng.next_normal(); let next = s * ((-0.5 * realised * realised) * dt + realised * dt.sqrt() * z).exp(); hedge_pnl += delta * (next - s); sum_squared_returns += (next / s).ln().powi(2); s = next; } VolTrade { pnl: premium + hedge_pnl - (s - k).max(0.0), gamma_weighted, realised_variance: sum_squared_returns / t, }} /// Sell a call, hedge it `steps` times, and see what is left.////// The strategy is the one the no-arbitrage chapter's replication argument describes: hold/// `delta` shares, funded from the premium, and adjust at each rebalancing date./// With zero rates the accounting is a single sum — the premium taken in, plus/// the gains on the shares held over each interval, less the payoff owed.////// `drift` is the *real world* drift of the stock, and it is a parameter on/// purpose. The Black-Scholes price does not contain it, which is the single/// most surprising claim in the no-arbitrage chapter, and the way to see that/// the claim is true rather than merely derived is to change the drift and/// watch the hedged result not move.pub fn delta_hedge( s0: f64, k: f64, sigma: f64, t: f64, drift: f64, steps: usize, seed: u64,) -> HedgeResult { let premium = black76(s0, k, sigma, t, Side::Call); let dt = t / steps as f64; let mut rng = Rng::new(seed); let mut s = s0; let mut hedge_pnl = 0.0; for i in 0..steps { let remaining = t - i as f64 * dt; let delta = call_delta(s, k, sigma, remaining); // One step of geometric Brownian motion, exactly rather than by Euler: // the discretisation error being measured is the hedging error, and it // would be contaminated by an approximation to the path itself. let z = rng.next_normal(); let next = s * ((drift - 0.5 * sigma * sigma) * dt + sigma * dt.sqrt() * z).exp(); hedge_pnl += delta * (next - s); s = next; } let payoff = (s - k).max(0.0); HedgeResult { pnl: premium + hedge_pnl - payoff, unhedged_pnl: premium - payoff, }} /// Run many hedges and report the spread of the outcomes.////// Returns the mean and the standard deviation of the profit and loss. The mean/// says whether the premium was right; the standard deviation says how well the/// hedge worked.pub fn hedge_statistics( s0: f64, k: f64, sigma: f64, t: f64, drift: f64, steps: usize, paths: usize, seed: u64,) -> (f64, f64) { let results: Vec<f64> = (0..paths) .map(|i| delta_hedge(s0, k, sigma, t, drift, steps, seed.wrapping_add(i as u64 * 7919)).pnl) .collect(); let mean = results.iter().sum::<f64>() / paths as f64; let var = results.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / paths as f64; (mean, var.sqrt())} #[cfg(test)]mod hedging_tests { use super::*; const S0: f64 = 100.0; const K: f64 = 100.0; const SIGMA: f64 = 0.2; const T: f64 = 1.0; const PATHS: usize = 4000; #[test] // The hedged position should break even in expectation. If it did not, // the Black-Scholes price would be the wrong price, and the whole of // the no-arbitrage chapter with it. let (mean, sd) = hedge_statistics(S0, K, SIGMA, T, 0.0, 250, PATHS, 1); let standard_error = sd / (PATHS as f64).sqrt(); assert!( mean.abs() < 3.0 * standard_error, "mean P&L {mean} against a standard error of {standard_error}" ); } #[test] fn hedging_error_falls_as_the_square_root_of_the_rebalancing_count() { // The classical rate. Quadrupling the number of rebalances should halve // the spread of the outcomes, and does — which is also the reason the // continuous-time argument of the no-arbitrage chapter is a limit rather than a // description of anything anybody does. let sd = |steps: usize| hedge_statistics(S0, K, SIGMA, T, 0.0, steps, PATHS, 2).1; for steps in [8, 32, 128] { let ratio = sd(steps) / sd(steps * 4); assert!( (ratio - 2.0).abs() < 0.25, "quadrupling from {steps} rebalances scaled the error by {ratio}" ); } } #[test] fn the_hedge_removes_the_drift_and_the_unhedged_position_does_not() { // The claim that makes the no-arbitrage chapter surprising: the price // contains no drift and does not need one. Change the real-world drift // over a wide range and the hedged result barely moves, while the // unhedged result moves a great deal. let hedged = |mu: f64| hedge_statistics(S0, K, SIGMA, T, mu, 500, PATHS, 3).0; let unhedged = |mu: f64| { let total: f64 = (0..PATHS) .map(|i| { delta_hedge(S0, K, SIGMA, T, mu, 500, 3u64.wrapping_add(i as u64 * 7919)) .unhedged_pnl }) .sum(); total / PATHS as f64 }; let spread_hedged = (hedged(0.20) - hedged(-0.20)).abs(); let spread_unhedged = (unhedged(0.20) - unhedged(-0.20)).abs(); assert!( spread_unhedged > 10.0 * spread_hedged, "hedged moved by {spread_hedged} and unhedged by {spread_unhedged}" ); } #[test] fn a_short_option_hedged_discretely_loses_on_the_tails() { // The hedging error is not symmetric. A short gamma position rebalances // against itself — buying after a rise, selling after a fall — so a // large move over one interval costs more than a quiet one saves, and // the distribution of outcomes is left-skewed. Anyone selling options // and hedging discretely is short that skew, which is why the error // matters beyond its standard deviation. let steps = 12; let pnl: Vec<f64> = (0..PATHS) .map(|i| delta_hedge(S0, K, SIGMA, T, 0.0, steps, 4u64.wrapping_add(i as u64 * 7919)).pnl) .collect(); let mean = pnl.iter().sum::<f64>() / PATHS as f64; let sd = (pnl.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / PATHS as f64).sqrt(); let skew = pnl.iter().map(|x| ((x - mean) / sd).powi(3)).sum::<f64>() / PATHS as f64; assert!(skew < -0.3, "the hedging error was not left-skewed, skew was {skew}"); }} // ---------------------------------------------------------------------------// Jumps.// --------------------------------------------------------------------------- /// Arrival times of a Poisson process of intensity `lambda` on `[0, t]`.////// Exactly, by drawing exponential gaps, rather than by thinning a fine grid./// The figure this feeds is about the difference between a process that moves/// continuously and one that does not, and a jump smeared across a grid cell/// would be arguing the opposite of the point.pub fn poisson_arrivals(t: f64, lambda: f64, seed: u64) -> Vec<f64> { let mut rng = Rng::new(seed); let mut arrivals = Vec::new(); let mut clock = 0.0; loop { // Inverse transform: -ln(U)/lambda is exponential with rate lambda. let u = { let bits = rng.next_u64() >> 11; (bits as f64 + 0.5) / (1u64 << 53) as f64 }; clock += -u.ln() / lambda; if clock > t { return arrivals; } arrivals.push(clock); }} /// The compensated Poisson process `M_t = N_t - lambda t`, sampled on a grid.////// The price process chapter's counterexample to Levy's characterization: a/// martingale whose variance accumulates at exactly the rate Brownian motion's/// does, and which is nothing like a Brownian motion.pub fn compensated_poisson(times: &[f64], arrivals: &[f64], lambda: f64) -> Vec<f64> { times .iter() .map(|&t| { let n = arrivals.partition_point(|&a| a <= t) as f64; n - lambda * t }) .collect()} #[cfg(test)]mod jump_tests { use super::*; const T: f64 = 5.0; const LAMBDA: f64 = 1.0; const TRIALS: usize = 20_000; #[test] fn the_count_has_the_right_mean_and_variance() { // Both equal lambda t for a Poisson process, which is the coincidence // the counterexample turns on. let counts: Vec<f64> = (0..TRIALS) .map(|i| poisson_arrivals(T, LAMBDA, 500 + i as u64).len() as f64) .collect(); let mean = counts.iter().sum::<f64>() / TRIALS as f64; let var = counts.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / TRIALS as f64; let expected = LAMBDA * T; // Standard error of the mean is sqrt(lambda t / trials). let se = (expected / TRIALS as f64).sqrt(); assert!((mean - expected).abs() < 4.0 * se, "mean {mean} against {expected}"); assert!((var - expected).abs() < 0.1 * expected, "variance {var} against {expected}"); } #[test] fn the_compensated_process_is_a_martingale() { // Zero mean at every horizon, which is what compensating buys. let grid: Vec<f64> = (0..=10).map(|i| i as f64 * T / 10.0).collect(); let mut totals = vec![0.0; grid.len()]; for i in 0..TRIALS { let a = poisson_arrivals(T, LAMBDA, 900 + i as u64); for (acc, m) in totals.iter_mut().zip(compensated_poisson(&grid, &a, LAMBDA)) { *acc += m; } } for (t, total) in grid.iter().zip(&totals) { let mean = total / TRIALS as f64; let se = (LAMBDA * t / TRIALS as f64).sqrt().max(1e-9); assert!(mean.abs() < 4.0 * se, "at t={t} the mean was {mean}"); } } #[test] fn it_has_brownian_variance_but_is_not_brownian() { // The counterexample, as two assertions. // // Its variance at T matches a Brownian motion's, so the hypothesis of // Levy's characterization about the accumulated variance is satisfied. let finals: Vec<f64> = (0..TRIALS) .map(|i| { let a = poisson_arrivals(T, LAMBDA, 1_300 + i as u64); a.len() as f64 - LAMBDA * T }) .collect(); let var = finals.iter().map(|x| x * x).sum::<f64>() / TRIALS as f64; assert!((var - T).abs() < 0.1 * T, "variance {var} against Brownian's {T}"); // And it is not normal: a Poisson count is skewed, where a Gaussian is // symmetric. Skewness of N is 1/sqrt(lambda t), about 0.45 here. let mean = finals.iter().sum::<f64>() / TRIALS as f64; let sd = var.sqrt(); let skew = finals.iter().map(|x| ((x - mean) / sd).powi(3)).sum::<f64>() / TRIALS as f64; assert!(skew > 0.25, "expected visible positive skew, got {skew}"); } #[test] fn the_realised_and_predictable_variations_differ() { // For a continuous martingale these agree, which is why the chapter can // write one symbol. Here the realised variation is the jump count --- a // random staircase --- while the predictable one is the straight line // lambda t. let a = poisson_arrivals(T, LAMBDA, 7); let realised = a.len() as f64; // every jump contributes 1^2 let predictable = LAMBDA * T; assert!( (realised - predictable).abs() > 1e-9, "this path happened to have exactly {predictable} jumps; pick another seed" ); }} /// What a hedged short option does when the world has the variance the price/// assumed and a different shape.////// The Black-Scholes chapter argues that the premium compensates the seller for/// the average of `(dS)^2` and for nothing else, so a true distribution with the/// same average and fatter tails hands over the same premium against a worse/// experience. That is a claim about two distributions rather than about one/// model, and it needs both simulated to be seen.////// The comparison holds *total variance* fixed. The diffusive case is geometric/// Brownian motion at volatility `sigma`. The jump case splits the same variance/// between a smaller diffusion and compensated lognormal jumps arriving at rate/// `lambda`, so that////// ```text/// sigma_diffusive^2 + lambda * E[(jump return)^2] = sigma^2 ,/// ```////// which leaves the option's Black-Scholes price unchanged --- the seller/// charges the same premium in both worlds, because the premium sees only the/// variance.pub struct HedgeShape { pub mean: f64, pub sd: f64, /// The worst one per cent of outcomes, averaged: the seller's tail. pub expected_shortfall: f64, /// The single worst outcome seen. pub worst: f64, /// Fraction of paths that made money, which is the other half of the story. pub win_rate: f64, /// The outcome at each percentile, from worst to best: the quantile curve, /// which is the shape the chapter's figure draws. pub quantiles: Vec<(f64, f64)>,} /// Hedge a short call under a jump diffusion whose total variance is `sigma^2`.////// `jump_share` is the fraction of the variance carried by the jumps, so zero/// recovers ordinary geometric Brownian motion and the two cases run through one/// piece of code. Jumps are lognormal with zero mean return and are compensated,/// so the underlying stays a martingale and the comparison is not confounded by/// a drift.pub fn hedge_under_jumps( s0: f64, k: f64, sigma: f64, t: f64, jump_share: f64, lambda: f64, steps: usize, paths: usize, seed: u64,) -> HedgeShape { let premium = black76(s0, k, sigma, t, Side::Call); let dt = t / steps as f64; // Split the variance. A jump has lognormal return exp(y) - 1 with y centred // so that E[exp(y)] = 1; its second moment is then exp(nu^2) - 1, and // matching lambda (exp(nu^2) - 1) to the jump share of the variance fixes // nu. let jump_variance = sigma * sigma * jump_share; let diffusive = (sigma * sigma - jump_variance).sqrt(); let nu = if jump_share > 0.0 { (1.0 + jump_variance / lambda).ln().sqrt() } else { 0.0 }; let mut rng = Rng::new(seed); let mut results = Vec::with_capacity(paths); for _ in 0..paths { let mut s = s0; let mut hedge_pnl = 0.0; for i in 0..steps { let remaining = t - i as f64 * dt; let delta = call_delta(s, k, sigma, remaining); // Diffusive step, compensated for the jump drift so the whole thing // is a martingale. let z = rng.next_normal(); let mut next = s * ((-0.5 * diffusive * diffusive) * dt + diffusive * dt.sqrt() * z).exp(); // At most one jump per step, which is the standard thinning and is // accurate when lambda * dt is small. if jump_share > 0.0 && rng.next_uniform() < lambda * dt { let y = -0.5 * nu * nu + nu * rng.next_normal(); next *= y.exp(); } hedge_pnl += delta * (next - s); s = next; } results.push(premium + hedge_pnl - (s - k).max(0.0)); } results.sort_by(|a, b| a.partial_cmp(b).unwrap()); let n = results.len(); let mean = results.iter().sum::<f64>() / n as f64; let var = results.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / n as f64; let tail = (n / 100).max(1); // One point per half per cent, which is enough to show the tail without // carrying forty thousand numbers into a chart. let quantiles = (1..200) .map(|i| { let q = i as f64 / 200.0; (q, results[((q * n as f64) as usize).min(n - 1)]) }) .collect(); HedgeShape { mean, sd: var.sqrt(), expected_shortfall: results[..tail].iter().sum::<f64>() / tail as f64, worst: results[0], win_rate: results.iter().filter(|x| **x > 0.0).count() as f64 / n as f64, quantiles, }} #[cfg(test)]mod shape_tests { use super::*; #[test] // The claim the Black-Scholes chapter makes and had no figure for. // Total variance is held fixed, so the premium is identical by // construction; what changes is the shape of what the seller lives // through. let (s0, k, sigma, t) = (100.0, 100.0, 0.20, 1.0); let smooth = hedge_under_jumps(s0, k, sigma, t, 0.0, 0.0, 250, 40_000, 20260813); let jumpy = hedge_under_jumps(s0, k, sigma, t, 0.5, 2.0, 250, 40_000, 20260813); // Both are fair: the premium compensates the average of (dS)^2 and the // average is the same. assert!(smooth.mean.abs() < 0.15, "smooth mean {:.4}", smooth.mean); assert!(jumpy.mean.abs() < 0.25, "jumpy mean {:.4}", jumpy.mean); // And the experience is not the same. assert!( jumpy.expected_shortfall < smooth.expected_shortfall * 1.8, "the jump tail should be far worse: {:.3} against {:.3}", jumpy.expected_shortfall, smooth.expected_shortfall ); assert!( jumpy.win_rate > smooth.win_rate, "and it should win more often while losing more: {:.3} against {:.3}", jumpy.win_rate, smooth.win_rate ); println!( "smooth: mean {:.3} sd {:.3} ES1% {:.3} worst {:.2} win {:.1}%", smooth.mean, smooth.sd, smooth.expected_shortfall, smooth.worst, smooth.win_rate * 100.0 ); println!( "jumpy: mean {:.3} sd {:.3} ES1% {:.3} worst {:.2} win {:.1}%", jumpy.mean, jumpy.sd, jumpy.expected_shortfall, jumpy.worst, jumpy.win_rate * 100.0 ); }} #[cfg(test)]mod discretisation_skew_tests { use super::*; /// Skewness of the hedged profit and loss at a given rebalancing frequency. fn skewness(n: usize, paths: usize) -> f64 { let v: Vec<f64> = (0..paths) .map(|i| delta_hedge(100.0, 100.0, 0.2, 1.0, 0.0, n, 202u64.wrapping_add(i as u64 * 7919)).pnl) .collect(); let m = v.iter().sum::<f64>() / v.len() as f64; let sd = (v.iter().map(|x| (x - m).powi(2)).sum::<f64>() / v.len() as f64).sqrt(); v.iter().map(|x| ((x - m) / sd).powi(3)).sum::<f64>() / v.len() as f64 } #[test] fn the_discretisation_skew_shrinks_as_the_hedge_gets_finer() { // The Black-Scholes chapter quotes both numbers, and the second is the // one that matters: this left tail is the cost of trading discretely, // and it goes away as the trading gets less discrete. That is what // separates it from the tail a jump produces, which does not. let coarse = skewness(8, 40_000); let fine = skewness(128, 40_000); assert!( (coarse - (-0.42)).abs() < 0.04, "eight rebalances: skewness {coarse:.3}" ); assert!( (fine - (-0.18)).abs() < 0.04, "a hundred and twenty-eight: skewness {fine:.3}" ); assert!(coarse < fine, "the skew must shrink with frequency"); assert!(fine < 0.0, "and remain negative"); }} #[cfg(test)]mod vol_trade_tests { use super::*; const S0: f64 = 100.0; const K: f64 = 100.0; const T: f64 = 1.0; /// The claim the alpha chapter makes about what a delta-hedged short /// volatility position actually is. The hedge accounting and the gamma /// weighted variance integral share no line of code --- one forms a payoff /// and uses the delta, the other never forms a payoff and uses the gamma --- /// so their agreement is a check on the decomposition rather than a /// restatement of it. #[test] fn the_profit_is_gamma_weighted_variance() { let mut worst = 0.0f64; for seed in 0..40u64 { let r = sell_variance(S0, K, 0.22, 0.18, T, 8_000, 3_000 + seed); // Scaled by the premium, which is the size of the thing being // decomposed; an absolute tolerance would mean nothing. let scale = black76(S0, K, 0.22, T, Side::Call); worst = worst.max((r.pnl - r.gamma_weighted).abs() / scale); } assert!(worst < 0.06, "worst relative gap between the two routes was {worst}"); } /// Selling variance above what arrives is profitable on average, and the /// average is what the premium is for. #[test] fn selling_rich_variance_pays_on_average() { let mean = |implied: f64, realised: f64| { let n = 2_000; (0..n) .map(|i| sell_variance(S0, K, implied, realised, T, 500, 700 + i as u64).pnl) .sum::<f64>() / n as f64 }; let rich = mean(0.22, 0.18); let fair = mean(0.20, 0.20); assert!(rich > 0.0, "selling rich variance averaged {rich}"); assert!(fair.abs() < 0.15, "selling fair variance averaged {fair}, expected near zero"); assert!(rich > fair, "rich {rich} should beat fair {fair}"); } /// The average profit is not a new number: it is the difference between the /// two Black-Scholes premiums, which is what "selling volatility rich" means /// once the hedge has removed everything else. #[test] let n = 4_000; let mean: f64 = (0..n) .map(|i| sell_variance(S0, K, 0.22, 0.18, T, 500, 700 + i as u64).pnl) .sum::<f64>() / n as f64; let spread = black76(S0, K, 0.22, T, Side::Call) - black76(S0, K, 0.18, T, Side::Call); assert!( (mean - spread).abs() < 0.01 * spread, "mean profit {mean} against premium difference {spread}" ); } /// And the trade is not a bet on realised variance. Restricting to paths that /// realised the same variance to four decimal places leaves a spread of /// outcomes wider than the whole of the average edge, because gamma weights /// a move by where it happened relative to the strike and to expiry. /// /// This is the number the alpha chapter quotes, and it is the reason selling /// variance and selling a hedged option are different trades. #[test] fn equal_realised_variance_does_not_mean_equal_profit() { let mut matched: Vec<f64> = Vec::new(); for seed in 0..4_000u64 { let r = sell_variance(S0, K, 0.22, 0.18, T, 500, seed); if (r.realised_variance - 0.18 * 0.18).abs() < 1e-4 { matched.push(r.pnl); } } assert!(matched.len() > 100, "only {} matched paths", matched.len()); let lo = matched.iter().cloned().fold(f64::INFINITY, f64::min); let hi = matched.iter().cloned().fold(f64::NEG_INFINITY, f64::max); let mean = matched.iter().sum::<f64>() / matched.len() as f64; assert!((lo - 0.49).abs() < 0.02, "lowest matched profit was {lo}"); assert!((hi - 2.94).abs() < 0.02, "highest matched profit was {hi}"); // The point: the dispersion left after fixing realised variance is // larger than the edge being collected. assert!( hi - lo > mean, "spread {} should exceed the average edge {mean}", hi - lo ); }}