quant/src/measure.rs
What a change of measure actually does to a distribution.
//! What a change of measure actually does to a distribution.//!//! The numeraires chapter is the most abstract chapter in the notes and the one//! whose content is easiest to state in a picture. Girsanov's theorem does not//! move any path. It reweights them. The set of futures the world can take is//! fixed; changing measure changes only how much each one counts.//!//! So the arithmetic here is deliberately arranged to make that checkable rather//! than assertable. A sample of terminal prices is drawn under the real-world//! measure and never touched again; each draw is multiplied by the Radon-Nikodym//! derivative; and the reweighted histogram is compared against the risk-neutral//! density computed independently in closed form. If the two agree, the claim in//! the previous paragraph is true, and if they did not, no amount of prose would//! make it true. use crate::black::norm_cdf;use crate::pathwise::Rng; /// The lognormal density of `S_T` when the drift is `drift`.////// The same function serves both measures: real-world and risk-neutral differ/// only in what is passed here, which is itself the point of the chapter.pub fn terminal_density(x: f64, spot: f64, drift: f64, sigma: f64, t: f64) -> f64 { if x <= 0.0 || sigma <= 0.0 || t <= 0.0 { return 0.0; } let variance = sigma * sigma * t; let z = (x / spot).ln() - (drift - 0.5 * sigma * sigma) * t; (-z * z / (2.0 * variance)).exp() / (x * variance.sqrt() * (2.0 * std::f64::consts::PI).sqrt())} /// The market price of risk, `(mu - r) / sigma`.////// How much excess return the asset offers per unit of the risk it carries./// Girsanov charges exactly this to move between the two measures, which is the/// formal statement of the idea that the risk-neutral measure is the real one/// with the risk premium removed.pub fn market_price_of_risk(mu: f64, r: f64, sigma: f64) -> f64 { (mu - r) / sigma} /// What a market of several assets over several factors says about risk prices.////// With `n` assets driven by `d` Brownian motions, writing `sigma` for the/// matrix of loadings and `mu - r` for the vector of excess returns, requiring/// every asset to drift at `r` under one measure is the linear system////// ```text/// mu - r = sigma theta/// ```////// for a single vector `theta` of factor prices. The numeraires chapter's point/// is that this is the whole of the fundamental theorem written in coordinates:/// a risk neutral measure exists exactly when the system is solvable, and is/// unique exactly when the solution is.pub struct RiskPrices { /// One price per Brownian factor. A property of the factor, not of any /// asset: every asset's excess return is its own loadings against this. pub theta: Vec<f64>, /// What the factor prices cannot explain. Zero when a risk neutral measure /// exists; otherwise it is itself the arbitrage, as portfolio weights. pub residual: Vec<f64>, /// Whether `theta` is pinned. False when there are more factors than the /// assets can span, which is an incomplete market and a family of measures. pub unique: bool,} impl RiskPrices { /// Whether the excess returns are consistent with a single set of factor /// prices, and so with no arbitrage. pub fn consistent(&self) -> bool { self.residual.iter().all(|r| r.abs() < 1e-9) } /// The excess return of the residual portfolio, which is `|residual|^2`. /// /// Positive exactly when the market is inconsistent, and it is riskless --- /// see the test, which checks the portfolio has no factor exposure at all. pub fn arbitrage_profit(&self) -> f64 { self.residual.iter().map(|r| r * r).sum() }} /// Solve `sigma theta = excess` in the least squares sense, and report what is/// left over.////// The leftover is the interesting part. Least squares makes the residual/// orthogonal to the column space of `sigma`, which is to say a portfolio with/// those weights has *no exposure to any factor* --- and its excess return is/// the squared length of the residual, which is positive whenever the residual/// is. So when no set of factor prices explains the market, the failure hands/// over the arbitrage directly rather than merely asserting one exists.////// `loadings[i]` is asset `i`'s row of `sigma`.pub fn factor_risk_prices(loadings: &[Vec<f64>], excess: &[f64]) -> RiskPrices { let n = loadings.len(); let d = if n == 0 { 0 } else { loadings[0].len() }; assert_eq!(excess.len(), n); // Normal equations: (sigma^T sigma) theta = sigma^T excess. let mut a = vec![vec![0.0f64; d + 1]; d]; for i in 0..d { for j in 0..d { a[i][j] = (0..n).map(|k| loadings[k][i] * loadings[k][j]).sum(); } a[i][d] = (0..n).map(|k| loadings[k][i] * excess[k]).sum(); } // Gaussian elimination with pivoting, tracking rank so a flat direction is // reported rather than silently resolved. let mut theta = vec![0.0; d]; let mut rank = 0; let mut pivot_of = vec![usize::MAX; d]; let mut row = 0; for col in 0..d { let best = (row..d).max_by(|&x, &y| { a[x][col].abs().partial_cmp(&a[y][col].abs()).unwrap_or(std::cmp::Ordering::Equal) }); let Some(best) = best else { break }; if a[best][col].abs() < 1e-12 { continue; } a.swap(row, best); pivot_of[col] = row; for k in 0..d { if k != row { let factor = a[k][col] / a[row][col]; for j in col..=d { a[k][j] -= factor * a[row][j]; } } } rank += 1; row += 1; if row == d { break; } } // Free variables are left at zero: one solution out of the family. for col in 0..d { if pivot_of[col] != usize::MAX { let r = pivot_of[col]; theta[col] = a[r][d] / a[r][col]; } } let residual: Vec<f64> = (0..n) .map(|k| excess[k] - (0..d).map(|j| loadings[k][j] * theta[j]).sum::<f64>()) .collect(); RiskPrices { theta, residual, unique: rank == d }} /// The Radon-Nikodym derivative `dQ/dP`, as a function of the terminal price.////// Every path ending at `x` carries this weight. It is a decreasing function of/// `x` whenever the asset earns a risk premium: the risk-neutral measure counts/// the good outcomes for less, which is the whole of what "removing the risk/// premium" means once it is written down.pub fn radon_nikodym(x: f64, spot: f64, mu: f64, r: f64, sigma: f64, t: f64) -> f64 { if x <= 0.0 || sigma <= 0.0 || t <= 0.0 { return 1.0; } let theta = market_price_of_risk(mu, r, sigma); // The standard normal draw that produced this terminal price under P. let z = ((x / spot).ln() - (mu - 0.5 * sigma * sigma) * t) / (sigma * t.sqrt()); (-theta * t.sqrt() * z - 0.5 * theta * theta * t).exp()} /// The reweighted histogram of a sample drawn under the real-world measure.////// `grid` supplies the bin centres. Returns a density: the total weight landing/// in each bin, divided by the bin width and by the total weight of the sample.////// Nothing here knows the risk-neutral density. The sample is drawn once under/// `mu`, and the only risk-neutral quantity used is the weight. The numeraires/// chapter's claim is that this reproduces the risk-neutral density anyway, and/// [`tests::reweighting_a_real_world_sample_gives_the_risk_neutral_density`]/// checks it.pub fn reweighted_histogram( grid: &[f64], spot: f64, mu: f64, r: f64, sigma: f64, t: f64, paths: usize, seed: u64,) -> Vec<f64> { let n = grid.len(); if n < 2 { return vec![0.0; n]; } let width = (grid[n - 1] - grid[0]) / (n - 1) as f64; let (lo, hi) = (grid[0] - 0.5 * width, grid[n - 1] + 0.5 * width); let mut bins = vec![0.0; n]; let mut total = 0.0; let mut rng = Rng::new(seed); for _ in 0..paths { // One draw under the real-world measure. This is the sample, and it is // never redrawn: the risk-neutral column below comes from reweighting // these same numbers. let z = rng.next_normal(); let x = spot * ((mu - 0.5 * sigma * sigma) * t + sigma * t.sqrt() * z).exp(); let weight = radon_nikodym(x, spot, mu, r, sigma, t); total += weight; if x >= lo && x < hi { let bin = (((x - lo) / width) as usize).min(n - 1); bins[bin] += weight; } } if total <= 0.0 { return bins; } for b in &mut bins { *b /= total * width; } bins} /// The probability of finishing above `strike`, under two different measures.////// Returns `(N(d1), N(d2))`. These are the two terms of the Black formula, and/// the numeraires chapter derives them as the probability of the *same event*/// under the share measure and under the `T`-forward measure. They are/// different numbers because the measures are different, not because the event/// is.pub fn exercise_probabilities(forward: f64, strike: f64, sigma: f64, t: f64) -> (f64, f64) { if strike <= 0.0 { return (1.0, 1.0); } if sigma <= 0.0 || t <= 0.0 { let exercised = if forward > strike { 1.0 } else { 0.0 }; return (exercised, exercised); } let vol = sigma * t.sqrt(); let d1 = ((forward / strike).ln() + 0.5 * vol * vol) / vol; (norm_cdf(d1), norm_cdf(d1 - vol))} #[cfg(test)]mod tests { use super::*; use crate::black::{black76, Side}; const SPOT: f64 = 100.0; const MU: f64 = 0.12; const R: f64 = 0.03; const SIGMA: f64 = 0.25; const T: f64 = 1.0; #[test] fn one_factor_two_assets_agreeing_on_its_price() { // The simplest market where the constraint has teeth. Two assets, one // Brownian motion, and the same Sharpe ratio: the factor has one price // and both assets are quoted consistently with it. let loadings = vec![vec![0.20], vec![0.35]]; let theta = 0.4; let excess = vec![0.20 * theta, 0.35 * theta]; let prices = factor_risk_prices(&loadings, &excess); assert!(prices.consistent(), "residual {:?}", prices.residual); assert!(prices.unique); assert!((prices.theta[0] - theta).abs() < 1e-12, "theta {:?}", prices.theta); } #[test] fn disagreeing_on_the_price_of_a_factor_is_an_arbitrage() { // The claim behind the whole construction. Two assets on one Brownian // motion with different Sharpe ratios cannot both be right, and the // failure is not abstract -- the residual is a portfolio, and it has no // exposure to the factor while earning a positive excess return. let loadings = vec![vec![0.20], vec![0.35]]; let excess = vec![0.20 * 0.4, 0.35 * 0.7]; let prices = factor_risk_prices(&loadings, &excess); assert!(!prices.consistent(), "the market was consistent after all"); // Riskless: the portfolio's loading on the factor is zero. let exposure: f64 = (0..2).map(|i| prices.residual[i] * loadings[i][0]).sum(); assert!(exposure.abs() < 1e-12, "portfolio kept exposure {exposure}"); // And strictly profitable, by exactly the squared length of the residual. let profit: f64 = (0..2).map(|i| prices.residual[i] * excess[i]).sum(); assert!(profit > 1e-6, "profit was {profit}"); assert!( (profit - prices.arbitrage_profit()).abs() < 1e-12, "profit {profit} against |residual|^2 {}", prices.arbitrage_profit() ); } #[test] fn a_complete_market_pins_every_factor_price() { // As many independent assets as factors: the loadings are invertible, // the system has one solution, and there is one risk neutral measure. // That is the second fundamental theorem in coordinates. let loadings = vec![vec![0.20, 0.05], vec![0.10, 0.30]]; let theta = [0.4, -0.2]; let excess: Vec<f64> = loadings .iter() .map(|row| row[0] * theta[0] + row[1] * theta[1]) .collect(); let prices = factor_risk_prices(&loadings, &excess); assert!(prices.consistent()); assert!(prices.unique, "a complete market should pin the prices"); for j in 0..2 { assert!((prices.theta[j] - theta[j]).abs() < 1e-10, "{:?}", prices.theta); } } #[test] fn fewer_assets_than_factors_leaves_the_prices_free() { // One asset, two factors. There is no arbitrage -- the single excess // return is easily explained -- but a whole line of factor prices // explains it, so there is a whole family of risk neutral measures. // That is an incomplete market, and it is why the unhedgeable part of a // payoff has no price rather than a wrong one. let loadings = vec![vec![0.20, 0.10]]; let excess = vec![0.05]; let prices = factor_risk_prices(&loadings, &excess); assert!(prices.consistent(), "a single asset cannot be arbitraged against itself"); assert!(!prices.unique, "the prices should not have been pinned"); } #[test] fn the_one_asset_case_is_the_ratio_everyone_writes() { // Consistency with the scalar definition used elsewhere in this module: // with one asset and one factor the factor price is (mu - r)/sigma. let (mu, r, sigma) = (0.12, 0.03, 0.25); let prices = factor_risk_prices(&vec![vec![sigma]], &[mu - r]); assert!( (prices.theta[0] - market_price_of_risk(mu, r, sigma)).abs() < 1e-12, "{:?}", prices.theta ); } #[test] fn the_weight_turns_one_density_into_the_other() { // The identity the whole chapter rests on: // // f_P(x) * (dQ/dP)(x) = f_Q(x) // // pointwise, for every terminal price. Not approximately, and not in // distribution -- at every single point. for i in 1..400 { let x = i as f64 * 1.0; let p = terminal_density(x, SPOT, MU, SIGMA, T); let q = terminal_density(x, SPOT, R, SIGMA, T); let w = radon_nikodym(x, SPOT, MU, R, SIGMA, T); assert!( (p * w - q).abs() < 1e-15 + 1e-12 * q, "at x={x}: {} against {q}", p * w ); } } #[test] fn the_weight_averages_to_one_under_the_real_world_measure() { // A Radon-Nikodym derivative is a probability density ratio, so it has // to integrate to one against the measure it is defined over. If it did // not, the change of measure would be creating or destroying // probability. let (lo, hi, steps) = (1e-6f64, 2000.0f64, 400_000); let h = (hi - lo) / steps as f64; let mut total = 0.0; for i in 0..steps { let x = lo + (i as f64 + 0.5) * h; total += terminal_density(x, SPOT, MU, SIGMA, T) * radon_nikodym(x, SPOT, MU, R, SIGMA, T) * h; } assert!((total - 1.0).abs() < 1e-6, "the weight integrated to {total}"); } #[test] fn each_density_integrates_to_one_and_has_the_right_mean() { // The sanity check on terminal_density itself, so that the identity // above is not two errors cancelling. for (drift, expected_mean) in [(MU, SPOT * (MU * T).exp()), (R, SPOT * (R * T).exp())] { let (lo, hi, steps) = (1e-6f64, 3000.0f64, 600_000); let h = (hi - lo) / steps as f64; let (mut mass, mut mean) = (0.0, 0.0); for i in 0..steps { let x = lo + (i as f64 + 0.5) * h; let d = terminal_density(x, SPOT, drift, SIGMA, T) * h; mass += d; mean += x * d; } assert!((mass - 1.0).abs() < 1e-6, "mass {mass} at drift {drift}"); assert!( (mean / expected_mean - 1.0).abs() < 1e-5, "mean {mean} against {expected_mean} at drift {drift}" ); } } #[test] fn reweighting_a_real_world_sample_gives_the_risk_neutral_density() { // The figure's claim, tested. Draw under P, never redraw, multiply by // the weight, and the histogram lands on a Q density that the sampling // never saw. let grid: Vec<f64> = (0..80).map(|i| 20.0 + i as f64 * 2.5).collect(); let histogram = reweighted_histogram(&grid, SPOT, MU, R, SIGMA, T, 2_000_000, 20_260_804); // Compared where there is enough sample to compare: the far tail of a // lognormal has too few draws in it to say anything about. The cutoff is // on the density rather than on the count, so it does not depend on how // many paths were run. let mut compared = 0; for (&x, &h) in grid.iter().zip(&histogram) { let analytic = terminal_density(x, SPOT, R, SIGMA, T); if analytic < 5e-4 { continue; } compared += 1; assert!( (h - analytic).abs() < 0.06 * analytic, "at x={x}: histogram {h} against analytic {analytic}" ); } assert!(compared > 40, "only {compared} bins had enough sample"); // And the test has to be capable of failing. The same histogram against // the density it was *drawn* from is out by a factor of three, so the // agreement above is the reweighting working rather than the two // lognormals being hard to tell apart. let worst = grid .iter() .zip(&histogram) .map(|(&x, &h)| { let drawn_from = terminal_density(x, SPOT, MU, SIGMA, T); if drawn_from < 5e-4 { 0.0 } else { (h - drawn_from).abs() / drawn_from } }) .fold(0.0f64, f64::max); assert!(worst > 1.0, "the two densities are indistinguishable: {worst}"); } #[test] fn reweighting_moves_the_mean_from_the_real_world_one_to_the_forward() { // The same statement at the level of the first moment, which is the one // a trader would check: under P the asset is expected to earn mu, and // after reweighting the same sample it earns r. let mut rng = Rng::new(77_777); let n = 500_000; let (mut plain, mut weighted, mut mass) = (0.0, 0.0, 0.0); for _ in 0..n { let z = rng.next_normal(); let x = SPOT * ((MU - 0.5 * SIGMA * SIGMA) * T + SIGMA * T.sqrt() * z).exp(); let w = radon_nikodym(x, SPOT, MU, R, SIGMA, T); plain += x; weighted += w * x; mass += w; } let real_world = plain / n as f64; let risk_neutral = weighted / mass; assert!( (real_world / (SPOT * (MU * T).exp()) - 1.0).abs() < 0.01, "real-world mean {real_world}" ); assert!( (risk_neutral / (SPOT * (R * T).exp()) - 1.0).abs() < 0.01, "risk-neutral mean {risk_neutral}" ); } #[test] fn the_two_exercise_probabilities_reproduce_the_black_formula() { // The cross-check that ties this module to the rest of the crate: if // N(d1) and N(d2) are what the numeraires chapter says they are, then // the Black price must fall out of them, and black76 was written // independently. let forward = SPOT * (R * T).exp(); for strike in [60.0, 90.0, 100.0, 115.0, 160.0] { let (n1, n2) = exercise_probabilities(forward, strike, SIGMA, T); let assembled = forward * n1 - strike * n2; let direct = black76(forward, strike, SIGMA, T, Side::Call); assert!( (assembled - direct).abs() < 1e-12, "at K={strike}: {assembled} against {direct}" ); } } #[test] // N(d1) > N(d2) at every strike, and the gap is the measure change. The // share measure weights by the terminal price itself, so it counts the // paths that finish high for more -- and those are exactly the paths on // which the option is exercised. let forward = SPOT * (R * T).exp(); for strike in [50.0, 80.0, 100.0, 130.0, 200.0] { let (n1, n2) = exercise_probabilities(forward, strike, SIGMA, T); assert!(n1 > n2, "at K={strike}: {n1} against {n2}"); } // And the two collapse together as the volatility vanishes, because with // no randomness there is nothing for the reweighting to bite on. let (n1, n2) = exercise_probabilities(forward, 100.0, 1e-6, T); assert!((n1 - n2).abs() < 1e-5, "{n1} against {n2}"); }}