quant/src/black.rs
Black-Scholes, Black-76 and Bachelier, plus the inversion from price back to implied volatility.
//! Black-Scholes, Black-76 and Bachelier, plus the inversion from price back to//! implied volatility.//!//! The Black-Scholes chapter derives the Black-Scholes formula; this is that//! formula, the forward-measure version of it the rates market quotes in, and//! the normal model that replaced the lognormal one when rates went negative.//! The local volatility chapter needs the inversion, because the implied//! volatility surface is *defined* by it: the smile is what you get when you//! insist on quoting every option price through a model that cannot produce//! those prices at a single volatility. /// Standard normal density.pub fn norm_pdf(x: f64) -> f64 { const INV_SQRT_2PI: f64 = 0.398_942_280_401_432_7; INV_SQRT_2PI * (-0.5 * x * x).exp()} /// Standard normal cumulative distribution function.////// Hart's rational approximation, in the arrangement West (2005) gives. Good to/// roughly machine precision across the whole line, which matters here: the/// implied volatility solve differentiates this, and a cheap approximation with/// a kink in the tail shows up as a visible wobble in the wings of a smile.pub fn norm_cdf(x: f64) -> f64 { let a = x.abs(); if a > 37.0 { return if x > 0.0 { 1.0 } else { 0.0 }; } let e = (-0.5 * a * a).exp(); let tail = if a < 7.071_067_811_865_47 { // Horner, written out a term at a time rather than as one nested // expression — the coefficients are impossible to check otherwise. let mut n = 3.526_249_659_989_11e-2; for c in [ 0.700_383_064_443_688, 6.373_962_203_531_65, 33.912_866_078_383, 112.079_291_497_871, 221.213_596_169_931, 220.206_867_912_376, ] { n = n * a + c; } let mut d = 8.838_834_764_831_84e-2; for c in [ 1.755_667_163_182_64, 16.064_177_579_207, 86.780_732_202_946_1, 296.564_248_779_674, 637.333_633_378_831, 793.826_512_519_948, 440.413_735_824_752, ] { d = d * a + c; } e * n / d } else { // Continued fraction; the rational form loses its accuracy out here. let f = a + 1.0 / (a + 2.0 / (a + 3.0 / (a + 4.0 / (a + 0.65)))); e / (f * 2.506_628_274_631) }; if x > 0.0 { 1.0 - tail } else { tail }} /// Which way round an option pays.#[derive(Clone, Copy, PartialEq, Eq, Debug)]pub enum Side { Call, Put,} impl Side { /// +1 for a call, -1 for a put. Both formulas below are the same expression /// with this factor threaded through, which is why they are written once. fn sign(self) -> f64 { match self { Side::Call => 1.0, Side::Put => -1.0, } }} /// Black-76: an option on a forward `f` under the measure in which the forward/// is a martingale, discounted separately by the caller.////// This is the undiscounted value — multiply by the numeraire (a bond price for/// a caplet, the annuity for a swaption) to get a price. Keeping the/// discounting outside is not tidiness: the curve construction chapter prices/// swaptions under the annuity measure precisely so that the payoff is a/// Black-76 payoff on the swap rate, and baking in a discount factor would hide/// that.pub fn black76(f: f64, k: f64, sigma: f64, t: f64, side: Side) -> f64 { let w = side.sign(); let v = sigma * t.sqrt(); // Zero variance: the option is worth its intrinsic value. if v <= 0.0 || f <= 0.0 || k <= 0.0 { return (w * (f - k)).max(0.0); } let d1 = ((f / k).ln() + 0.5 * v * v) / v; let d2 = d1 - v; w * (f * norm_cdf(w * d1) - k * norm_cdf(w * d2))} /// Vega of [`black76`]: the derivative of the price with respect to `sigma`.////// The implied volatility solve needs it, and it is also the reason the solve is/// well posed — vega is strictly positive, so price is strictly increasing in/// volatility and the inverse exists and is unique.pub fn black76_vega(f: f64, k: f64, sigma: f64, t: f64) -> f64 { let v = sigma * t.sqrt(); if v <= 0.0 || f <= 0.0 || k <= 0.0 { return 0.0; } let d1 = ((f / k).ln() + 0.5 * v * v) / v; f * norm_pdf(d1) * t.sqrt()} /// Black-Scholes on a spot underlying with continuous rate `r` and dividend/// yield `q`. Expressed through [`black76`] on the forward, which is what the/// change of numeraire in the Black-Scholes chapter says it is.pub fn black_scholes(s: f64, k: f64, sigma: f64, t: f64, r: f64, q: f64, side: Side) -> f64 { let df = (-r * t).exp(); let f = s * ((r - q) * t).exp(); df * black76(f, k, sigma, t, side)} /// Bachelier (normal) model: `dF = sigma_n dW`, so the forward is Gaussian.////// The rates market quotes in this because it stays well defined when `f` or `k`/// is negative, which the lognormal model does not. `sigma_n` is an absolute/// rate volatility — basis points per root year, not a percentage.pub fn bachelier(f: f64, k: f64, sigma_n: f64, t: f64, side: Side) -> f64 { let w = side.sign(); let v = sigma_n * t.sqrt(); if v <= 0.0 { return (w * (f - k)).max(0.0); } let d = (f - k) / v; w * (f - k) * norm_cdf(w * d) + v * norm_pdf(d)} /// Vega of [`bachelier`].pub fn bachelier_vega(f: f64, k: f64, sigma_n: f64, t: f64) -> f64 { let v = sigma_n * t.sqrt(); if v <= 0.0 { return 0.0; } norm_pdf((f - k) / v) * t.sqrt()} /// The volatility that reproduces `price` in the Black-76 model.////// Newton from a Brenner-Subrahmanyam start, kept inside a bracket that is/// halved whenever Newton tries to leave it. Newton alone is not safe here: vega/// collapses in the deep wings, and a step divided by a vega of 1e-12 leaves the/// positive half-line entirely.////// The solve is always done on the out-of-the-money option, moving to it by/// put-call parity when necessary. An in-the-money price is mostly intrinsic/// value, and the part that carries the volatility is the small remainder;/// inverting it directly throws away most of the significant digits before the/// solver starts.////// Returns `None` when the price is outside the no-arbitrage bounds, which is/// not a numerical failure but the useful answer — the local volatility chapter/// uses exactly this to show that a surface a reader might have written down by/// hand admits arbitrage. It also returns `None` when the out-of-the-money/// price has underflowed to zero, where the option carries no recoverable/// information about volatility at all.pub fn implied_vol_black76( price: f64, f: f64, k: f64, t: f64, side: Side,) -> Option<f64> { let w = side.sign(); let intrinsic = (w * (f - k)).max(0.0); // Upper bound: a call is worth at most the forward, a put at most the strike. let cap = match side { Side::Call => f, Side::Put => k, }; if !(price > intrinsic && price < cap) || t <= 0.0 || f <= 0.0 || k <= 0.0 { return None; } // Move to the out-of-the-money option: c - p = f - k. let in_the_money = w * (f - k) > 0.0; let (price, side) = if in_the_money { (price - w * (f - k), if side == Side::Call { Side::Put } else { Side::Call }) } else { (price, side) }; if !(price > 0.0) { return None; } let (mut lo, mut hi) = (1e-9_f64, 10.0_f64); // Brenner-Subrahmanyam: exact at the money, and a decent start elsewhere. let mut sigma = (2.0 * core::f64::consts::PI / t).sqrt() * price / f; sigma = sigma.clamp(lo, hi); for _ in 0..100 { let diff = black76(f, k, sigma, t, side) - price; if diff.abs() < 1e-12 { return Some(sigma); } // Price is increasing in sigma, so the sign of the error says which side // of the root we are on and tightens the bracket for free. if diff > 0.0 { hi = sigma; } else { lo = sigma; } let vega = black76_vega(f, k, sigma, t); let next = if vega > 1e-12 { sigma - diff / vega } else { f64::NAN }; sigma = if next.is_finite() && next > lo && next < hi { next } else { 0.5 * (lo + hi) }; } Some(sigma)} /// The volatility that reproduces `price` in the Bachelier model. Same scheme as/// [`implied_vol_black76`].pub fn implied_vol_bachelier(price: f64, f: f64, k: f64, t: f64, side: Side) -> Option<f64> { let w = side.sign(); let intrinsic = (w * (f - k)).max(0.0); if !(price > intrinsic) || t <= 0.0 { return None; } // Same move to the out-of-the-money option, by the same parity relation — // it holds in the normal model too, being a statement about payoffs. let in_the_money = w * (f - k) > 0.0; let (price, side) = if in_the_money { (price - w * (f - k), if side == Side::Call { Side::Put } else { Side::Call }) } else { (price, side) }; if !(price > 0.0) { return None; } let (mut lo, mut hi) = (1e-12_f64, 1.0_f64); // Normal vols are absolute rates, so a wide bracket in the same units as the // forward is safe; grow it rather than assume a range. while bachelier(f, k, hi, t, side) < price && hi < 1e6 { hi *= 2.0; } let mut sigma = 0.5 * (lo + hi); for _ in 0..100 { let diff = bachelier(f, k, sigma, t, side) - price; if diff.abs() < 1e-14 { return Some(sigma); } if diff > 0.0 { hi = sigma; } else { lo = sigma; } let vega = bachelier_vega(f, k, sigma, t); let next = if vega > 1e-14 { sigma - diff / vega } else { f64::NAN }; sigma = if next.is_finite() && next > lo && next < hi { next } else { 0.5 * (lo + hi) }; } Some(sigma)} /// The first-order risk report of the risk management chapter, for one/// Black-Scholes option.////// `theta` is the derivative with respect to calendar time, not with respect to/// time remaining, so it is negative for a long option and the theta-gamma/// identity of the Black-Scholes chapter can be read off directly.#[derive(Clone, Copy, Debug)]pub struct Greeks { pub price: f64, pub delta: f64, pub gamma: f64, pub vega: f64, pub theta: f64,} /// Price and sensitivities together, because they share `d1` and `d2` and/// because a risk report wants all of them at once.////// At zero time or zero volatility the option is worth its intrinsic value and/// the sensitivities are the derivatives of that: delta is a step, gamma is/// unbounded at the strike and zero elsewhere, and vega and theta vanish. The/// step's value exactly at the strike is a convention --- the derivative does/// not exist there --- and `0.5` is chosen because it is the limit approached/// from either side of a symmetric perturbation.pub fn greeks(s: f64, k: f64, sigma: f64, t: f64, r: f64, q: f64, side: Side) -> Greeks { let w = side.sign(); let v = sigma * t.sqrt(); let price = black_scholes(s, k, sigma, t, r, q, side); if v <= 0.0 || s <= 0.0 || k <= 0.0 { let in_the_money = w * (s - k); let delta = if in_the_money > 0.0 { w } else if in_the_money < 0.0 { 0.0 } else { 0.5 * w }; let gamma = if s == k { f64::INFINITY } else { 0.0 }; return Greeks { price, delta, gamma, vega: 0.0, theta: 0.0 }; } let df_r = (-r * t).exp(); let df_q = (-q * t).exp(); let f = s * ((r - q) * t).exp(); let d1 = ((f / k).ln() + 0.5 * v * v) / v; let d2 = d1 - v; let n_d1 = norm_pdf(d1); Greeks { price, delta: w * df_q * norm_cdf(w * d1), gamma: df_q * n_d1 / (s * v), vega: s * df_q * n_d1 * t.sqrt(), // d/dt, hence the leading minus on the decay term. theta: -s * df_q * n_d1 * sigma / (2.0 * t.sqrt()) - w * r * k * df_r * norm_cdf(w * d2) + w * q * s * df_q * norm_cdf(w * d1), }} #[cfg(test)]mod greek_tests { use super::*; /// The Black-Scholes chapter reduces its PDE to `Theta = -sigma^2 S^2 Gamma / 2` /// at zero rates. That is an identity rather than an approximation, so it is /// testable to machine precision, and it is the one relation between the /// Greeks the chapter asks the reader to carry. #[test] fn theta_is_gamma_at_zero_rates() { for &s in &[60.0, 90.0, 100.0, 110.0, 160.0] { for &t in &[0.05, 0.5, 2.0] { for side in [Side::Call, Side::Put] { let g = greeks(s, 100.0, 0.2, t, 0.0, 0.0, side); let want = -0.5 * 0.2 * 0.2 * s * s * g.gamma; assert!( (g.theta - want).abs() < 1e-9 * want.abs().max(1.0), "S={s} T={t} {side:?}: theta {} but -sigma^2 S^2 gamma/2 {want}", g.theta ); } } } } /// Every closed form against a differenced price, which shares none of the /// algebra above. #[test] fn the_closed_forms_agree_with_a_differenced_price() { const S: f64 = 100.0; const K: f64 = 95.0; const SIG: f64 = 0.25; const T: f64 = 0.75; const R: f64 = 0.03; const Q: f64 = 0.01; for side in [Side::Call, Side::Put] { let g = greeks(S, K, SIG, T, R, Q, side); let price = |s: f64, sig: f64, t: f64| black_scholes(s, K, sig, t, R, Q, side); let h = 1e-4 * S; let delta = (price(S + h, SIG, T) - price(S - h, SIG, T)) / (2.0 * h); let gamma = (price(S + h, SIG, T) - 2.0 * price(S, SIG, T) + price(S - h, SIG, T)) / (h * h); let hv = 1e-6; let vega = (price(S, SIG + hv, T) - price(S, SIG - hv, T)) / (2.0 * hv); let ht = 1e-6; // Calendar time runs against time to maturity. let theta = -(price(S, SIG, T + ht) - price(S, SIG, T - ht)) / (2.0 * ht); for (name, exact, differenced, tol) in [ ("delta", g.delta, delta, 1e-6), ("gamma", g.gamma, gamma, 1e-4), ("vega", g.vega, vega, 1e-4), ("theta", g.theta, theta, 1e-4), ] { assert!( (exact - differenced).abs() < tol * differenced.abs().max(1.0), "{side:?} {name}: closed form {exact}, differenced {differenced}" ); } } } /// What the figure in the Black-Scholes chapter shows: as expiry approaches /// the price falls onto the payoff, delta onto a step and gamma onto a spike. #[test] fn the_short_dated_limit_is_the_payoff() { const K: f64 = 100.0; for &s in &[80.0, 120.0] { let intrinsic = (s - K).max(0.0); let near = greeks(s, K, 0.2, 1e-6, 0.0, 0.0, Side::Call); assert!((near.price - intrinsic).abs() < 1e-6); assert!((near.delta - if s > K { 1.0 } else { 0.0 }).abs() < 1e-6); assert!(near.gamma < 1e-6, "gamma away from the strike was {}", near.gamma); } // At the strike the peak grows without bound as the maturity shortens. let mut previous = 0.0; for &t in &[0.25, 0.05, 0.01, 0.001] { let g = greeks(K, K, 0.2, t, 0.0, 0.0, Side::Call); assert!(g.gamma > previous, "gamma at the money fell at T={t}"); previous = g.gamma; } }} #[cfg(test)]mod tests { use super::*; #[test] fn norm_cdf_matches_known_values() { assert!((norm_cdf(0.0) - 0.5).abs() < 1e-15); assert!((norm_cdf(1.0) - 0.841_344_746_068_543).abs() < 1e-13); assert!((norm_cdf(-1.0) - 0.158_655_253_931_457).abs() < 1e-13); // Symmetry has to hold to full precision or the wings of a smile bend. for i in 0..80 { let x = i as f64 * 0.25; assert!((norm_cdf(x) + norm_cdf(-x) - 1.0).abs() < 1e-14); } } #[test] fn put_call_parity_holds() { let (f, t) = (100.0, 1.5); for k in [50.0, 90.0, 100.0, 110.0, 200.0] { let c = black76(f, k, 0.22, t, Side::Call); let p = black76(f, k, 0.22, t, Side::Put); assert!((c - p - (f - k)).abs() < 1e-10, "parity failed at k={k}"); } } #[test] fn implied_vol_inverts_black76() { // Deliberately includes in-the-money strikes: the parity step inside the // solver is what makes those recoverable at all. let (f, t) = (0.03, 5.0); let mut checked = 0; for k in [0.01, 0.02, 0.03, 0.045, 0.08] { for sigma in [0.2, 0.6, 1.2] { let price = black76(f, k, sigma, t, Side::Call); let back = implied_vol_black76(price, f, k, t, Side::Call) .unwrap_or_else(|| panic!("no inverse at k={k} sigma={sigma}")); assert!((back - sigma).abs() < 1e-8, "k={k} sigma={sigma} back={back}"); checked += 1; } } assert_eq!(checked, 15); } #[test] fn implied_vol_inverts_bachelier() { let (f, t) = (0.02, 3.0); for k in [-0.01, 0.0, 0.02, 0.05] { for sigma_n in [0.008, 0.03] { let price = bachelier(f, k, sigma_n, t, Side::Call); let back = implied_vol_bachelier(price, f, k, t, Side::Call) .unwrap_or_else(|| panic!("no inverse at k={k} sigma_n={sigma_n}")); assert!((back - sigma_n).abs() < 1e-10, "k={k} back={back}"); } } } #[test] fn implied_vol_is_none_where_the_price_carries_no_information() { // Five years, but a volatility so low that the option is worth its // intrinsic value to the last bit of a double. There is no inverse to // find, and reporting one would be reporting noise. let (f, k, t, sigma) = (0.03, 0.01, 5.0, 0.05); let price = black76(f, k, sigma, t, Side::Call); assert_eq!(price, f - k); assert!(implied_vol_black76(price, f, k, t, Side::Call).is_none()); } #[test] fn implied_vol_agrees_across_the_two_sides() { // A call and the put at the same strike must imply the same volatility, // since parity says they are the same price. let (f, t, sigma) = (100.0, 2.0, 0.28); for k in [70.0, 100.0, 140.0] { let c = implied_vol_black76(black76(f, k, sigma, t, Side::Call), f, k, t, Side::Call); let p = implied_vol_black76(black76(f, k, sigma, t, Side::Put), f, k, t, Side::Put); assert!((c.unwrap() - p.unwrap()).abs() < 1e-10, "k={k}"); } } #[test] fn arbitrageable_prices_have_no_implied_vol() { let (f, k, t) = (100.0, 90.0, 1.0); // Below intrinsic, and above the forward: both outside the bounds. assert!(implied_vol_black76(9.0, f, k, t, Side::Call).is_none()); assert!(implied_vol_black76(101.0, f, k, t, Side::Call).is_none()); } #[test] fn bachelier_handles_negative_rates() { // The case the lognormal model cannot express at all. let v = bachelier(-0.004, -0.002, 0.005, 2.0, Side::Call); assert!(v > 0.0 && v.is_finite()); }}