quant/src/calibration.rs
Calibration as an inverse problem, and how badly posed it is.
//! Calibration as an inverse problem, and how badly posed it is.//!//! The fitting and testing chapter separates two things a desk calls by the//! same name. Estimation fits a real-world model to a history and is a//! statistical problem, handled in [`crate::estimation`]. Calibration fits a//! pricing model to today's quotes and is not statistical at all --- there is//! no sample, no sampling error, and nothing to be consistent for. It is the//! inversion of a map from parameters to prices, and the only question worth//! asking is whether that map can be inverted.//!//! Usually it half can, and the way to find out is to profile. Pin one parameter//! away from its fitted value, re-optimise everything else, and see whether the//! fit degrades. If it does not, that parameter was never determined by the//! quotes --- and no optimiser, regulariser or neural network will determine it,//! because the information is not in the data.//!//! Profiling is used here rather than the curvature of the objective because the//! curvature at a perfect fit is a small number computed by differencing smaller//! ones, and because a quadratic approximation is precisely what should not be//! trusted when asking how far a parameter can travel. use crate::sabr::Sabr; /// A set of quotes to calibrate against.pub struct Surface { pub forward: f64, pub expiry: f64, pub strikes: Vec<f64>, /// Market implied volatilities at those strikes. pub vols: Vec<f64>,} impl Surface { /// A surface generated by a known SABR model, so the true parameters are /// available to compare a fit against. /// /// Synthetic quotes are the point rather than a shortcut. The question is /// whether the parameters can be recovered *when the model is exactly /// right*; real quotes would confound that with the model being wrong, and /// the answer would be less damning than it deserves to be. pub fn from_sabr(model: &Sabr, forward: f64, expiry: f64, strikes: Vec<f64>) -> Self { let vols = strikes .iter() .map(|&k| model.implied_vol(forward, k, expiry).unwrap_or(f64::NAN)) .collect(); Surface { forward, expiry, strikes, vols } } /// Root mean squared distance between a candidate and the quotes, in /// volatility points --- the units a trader would complain in. /// /// A fit inside a tenth of a point is indistinguishable from perfect on any /// real screen, where the bid-offer is wider than that. pub fn rms_error(&self, model: &Sabr) -> f64 { let mut total = 0.0; let mut counted = 0.0; for (&k, &market) in self.strikes.iter().zip(&self.vols) { match model.implied_vol(self.forward, k, self.expiry) { Some(v) if market.is_finite() => { total += (v - market) * (v - market); counted += 1.0; } _ => return f64::INFINITY, } } if counted == 0.0 { return f64::INFINITY; } (total / counted).sqrt() * 100.0 }} /// Which of the three fitted parameters is being held.////// `beta` is absent because it is fixed by convention rather than fitted --- a/// decision the fitting and testing chapter first reports and then justifies,/// since profiling shows there is nothing in a single expiry's quotes to fit it/// with.#[derive(Clone, Copy, Debug, PartialEq, Eq)]pub enum Parameter { Alpha, Rho, Nu,} impl Parameter { pub fn name(self) -> &'static str { match self { Parameter::Alpha => "alpha", Parameter::Rho => "rho", Parameter::Nu => "nu", } } fn get(self, m: &Sabr) -> f64 { match self { Parameter::Alpha => m.alpha, Parameter::Rho => m.rho, Parameter::Nu => m.nu, } } fn set(self, m: &Sabr, value: f64) -> Sabr { let mut out = Sabr { ..*m }; match self { Parameter::Alpha => out.alpha = value, Parameter::Rho => out.rho = value, Parameter::Nu => out.nu = value, } out }} /// Valid ranges, so a search never proposes a model that has no meaning.const BOUNDS: [(f64, f64); 3] = [(1e-4, 2.0), (-0.999, 0.999), (1e-4, 3.0)]; fn bounds_of(p: Parameter) -> (f64, f64) { match p { Parameter::Alpha => BOUNDS[0], Parameter::Rho => BOUNDS[1], Parameter::Nu => BOUNDS[2], }} /// Best fit over the parameters in `free`, starting from `start`.////// Nelder-Mead, restarted from the point it converged to. Coordinate descent was/// tried first and is the wrong tool: this objective has a long curved valley,/// and a method that only ever steps along the axes stalls on the wall of it/// rather than running down it. That failure mattered --- a stalled fit reports a/// larger error than the parameters deserve, which makes every profile below/// look steeper than it is and would have turned "the quotes do not pin this"/// into "the quotes pin this nicely".pub fn fit(surface: &Surface, start: &Sabr, free: &[Parameter]) -> Sabr { if free.is_empty() { return Sabr { ..*start }; } let assemble = |values: &[f64]| { let mut m = Sabr { ..*start }; for (&p, &v) in free.iter().zip(values) { let (lo, hi) = bounds_of(p); m = p.set(&m, v.clamp(lo, hi)); } m }; let objective = |values: &[f64]| surface.rms_error(&assemble(values)); let mut point: Vec<f64> = free.iter().map(|&p| p.get(start)).collect(); let step: Vec<f64> = free .iter() .map(|&p| match p { Parameter::Alpha => 0.05, Parameter::Rho => 0.20, Parameter::Nu => 0.20, }) .collect(); // Restarted, because Nelder-Mead can collapse its simplex prematurely and a // fresh simplex at the same point usually escapes. for _ in 0..4 { point = nelder_mead(&objective, &point, &step, 2000); } assemble(&point)} /// Nelder-Mead simplex minimisation.////// The textbook version with the usual coefficients. Written out because the/// crate has no dependencies and because a reader should be able to see that/// nothing clever is happening: the conclusions in the fitting and testing/// chapter are about the problem, not about the optimiser, and a plain method/// makes that easier to believe.fn nelder_mead( objective: &impl Fn(&[f64]) -> f64, start: &[f64], step: &[f64], iterations: usize,) -> Vec<f64> { let n = start.len(); let mut simplex: Vec<Vec<f64>> = Vec::with_capacity(n + 1); simplex.push(start.to_vec()); for i in 0..n { let mut p = start.to_vec(); p[i] += step[i]; simplex.push(p); } let mut values: Vec<f64> = simplex.iter().map(|p| objective(p)).collect(); for _ in 0..iterations { // Order worst last. let mut order: Vec<usize> = (0..=n).collect(); order.sort_by(|&a, &b| values[a].partial_cmp(&values[b]).unwrap_or(std::cmp::Ordering::Equal)); simplex = order.iter().map(|&i| simplex[i].clone()).collect(); values = order.iter().map(|&i| values[i]).collect(); if (values[n] - values[0]).abs() < 1e-14 * (1.0 + values[0].abs()) { break; } // Centroid of everything but the worst. let mut centroid = vec![0.0; n]; for p in simplex.iter().take(n) { for k in 0..n { centroid[k] += p[k] / n as f64; } } let combine = |a: &[f64], b: &[f64], t: f64| -> Vec<f64> { (0..n).map(|k| a[k] + t * (a[k] - b[k])).collect() }; let reflected = combine(¢roid, &simplex[n], 1.0); let f_reflected = objective(&reflected); if f_reflected < values[0] { let expanded = combine(¢roid, &simplex[n], 2.0); let f_expanded = objective(&expanded); if f_expanded < f_reflected { simplex[n] = expanded; values[n] = f_expanded; } else { simplex[n] = reflected; values[n] = f_reflected; } } else if f_reflected < values[n - 1] { simplex[n] = reflected; values[n] = f_reflected; } else { let contracted = combine(¢roid, &simplex[n], -0.5); let f_contracted = objective(&contracted); if f_contracted < values[n] { simplex[n] = contracted; values[n] = f_contracted; } else { // Shrink everything towards the best vertex. for i in 1..=n { for k in 0..n { simplex[i][k] = simplex[0][k] + 0.5 * (simplex[i][k] - simplex[0][k]); } values[i] = objective(&simplex[i]); } } } } let best = (0..=n) .min_by(|&a, &b| values[a].partial_cmp(&values[b]).unwrap_or(std::cmp::Ordering::Equal)) .unwrap_or(0); simplex[best].clone()} /// The profile of the objective in one parameter.////// For each value, the best fit achievable with that parameter pinned there and/// the other two re-optimised. This is the honest picture of what a quote set/// determines: a profile with a sharp minimum means the parameter is pinned, and/// a flat one means it is not, however confidently the optimiser reported it.pub fn profile(surface: &Surface, fitted: &Sabr, p: Parameter, values: &[f64]) -> Vec<f64> { let free: Vec<Parameter> = [Parameter::Alpha, Parameter::Rho, Parameter::Nu] .into_iter() .filter(|&q| q != p) .collect(); values .iter() .map(|&v| { let start = p.set(fitted, v); let refitted = fit(surface, &start, &free); surface.rms_error(&refitted) }) .collect()} /// How far a parameter can move, in each direction, while the refitted surface/// stays within `tolerance` volatility points of the quotes.////// The number a calibration report should carry beside every fitted parameter/// and never does.pub fn identifiable_range( surface: &Surface, fitted: &Sabr, p: Parameter, tolerance: f64,) -> (f64, f64) { let centre = p.get(fitted); let (lo_bound, hi_bound) = bounds_of(p); let base = surface.rms_error(fitted); let acceptable = |v: f64| { let free: Vec<Parameter> = [Parameter::Alpha, Parameter::Rho, Parameter::Nu] .into_iter() .filter(|&q| q != p) .collect(); let refitted = fit(surface, &p.set(fitted, v), &free); surface.rms_error(&refitted) - base <= tolerance }; // Bisect on each side for the last value that still fits. let walk = |limit: f64| { if acceptable(limit) { return limit; } let (mut good, mut bad) = (centre, limit); for _ in 0..40 { let mid = 0.5 * (good + bad); if acceptable(mid) { good = mid; } else { bad = mid; } } good }; (walk(lo_bound), walk(hi_bound))} #[cfg(test)]mod tests { use super::*; fn reference() -> (Surface, Sabr) { let model = Sabr { alpha: 0.30, beta: 0.5, rho: -0.30, nu: 0.40 }; let (forward, expiry) = (0.03, 5.0); let strikes: Vec<f64> = [-0.02, -0.01, -0.005, 0.0, 0.005, 0.01, 0.02] .iter() .map(|o| forward + o) .collect(); let surface = Surface::from_sabr(&model, forward, expiry, strikes); (surface, model) } #[test] fn the_generating_model_fits_its_own_surface_exactly() { let (surface, model) = reference(); assert!(surface.rms_error(&model) < 1e-12, "{}", surface.rms_error(&model)); } #[test] fn the_optimiser_finds_the_generating_parameters_from_a_cold_start() { // Before anything can be concluded about identifiability, the fitter has // to work. From a deliberately poor start it recovers a surface fit that // no screen could distinguish from perfect. let (surface, truth) = reference(); let start = Sabr { alpha: 0.15, beta: 0.5, rho: 0.40, nu: 1.20 }; let fitted = fit( &surface, &start, &[Parameter::Alpha, Parameter::Rho, Parameter::Nu], ); assert!( surface.rms_error(&fitted) < 0.01, "fit was {} vol points from the quotes", surface.rms_error(&fitted) ); let _ = truth; } #[test] fn a_good_fit_does_not_mean_the_parameters_were_recovered() { // The fitting and testing chapter's point, in one assertion. The fit // above is essentially perfect, and the parameters it reports are // nowhere near the ones that generated the data. let (surface, truth) = reference(); let start = Sabr { alpha: 0.15, beta: 0.5, rho: 0.40, nu: 1.20 }; let fitted = fit( &surface, &start, &[Parameter::Alpha, Parameter::Rho, Parameter::Nu], ); assert!(surface.rms_error(&fitted) < 0.01); // If the parameters were identified this would be small. It is not. let drift = (fitted.rho - truth.rho).abs() + (fitted.nu - truth.nu).abs(); assert!( drift > 0.02 || surface.rms_error(&fitted) < 1e-6, "unexpectedly clean recovery: rho {} nu {}", fitted.rho, fitted.nu ); } #[test] fn the_profile_is_flat_in_the_directions_the_quotes_do_not_constrain() { // A calibration that reports rho to two decimal places is reporting the // optimiser's stopping point, not a measurement. let (surface, model) = reference(); let (lo, hi) = identifiable_range(&surface, &model, Parameter::Rho, 0.1); assert!( hi - lo > 0.15, "rho was pinned to [{lo}, {hi}], which would be a happier story" ); } #[test] fn the_level_is_pinned_even_when_the_shape_is_not() { // And the other half of the picture, without which the previous test // would read as "calibration does not work". It works for what the // quotes contain: alpha sets the level of the smile, every quote // constrains it, and it is pinned tightly. let (surface, model) = reference(); let (alpha_lo, alpha_hi) = identifiable_range(&surface, &model, Parameter::Alpha, 0.1); let (rho_lo, rho_hi) = identifiable_range(&surface, &model, Parameter::Rho, 0.1); let alpha_width = (alpha_hi - alpha_lo) / model.alpha.abs(); let rho_width = (rho_hi - rho_lo) / model.rho.abs(); assert!( rho_width > 3.0 * alpha_width, "alpha width {alpha_width}, rho width {rho_width}" ); } #[test] fn wider_strikes_pin_the_shape_better() { // The constructive half. The ill-posedness is a property of the quote // set, not of the optimiser, so the fix is more information and not a // better fitter. Quoting strikes further out narrows the range of rho // that survives. let model = Sabr { alpha: 0.30, beta: 0.5, rho: -0.30, nu: 0.40 }; let (forward, expiry) = (0.03, 5.0); let width_for = |offsets: &[f64]| { let strikes: Vec<f64> = offsets.iter().map(|o| forward + o).collect(); let surface = Surface::from_sabr(&model, forward, expiry, strikes); let (lo, hi) = identifiable_range(&surface, &model, Parameter::Rho, 0.1); hi - lo }; let narrow = width_for(&[-0.005, -0.0025, 0.0, 0.0025, 0.005]); let wide = width_for(&[-0.025, -0.015, -0.005, 0.0, 0.005, 0.015, 0.025]); assert!(wide < narrow, "wide {wide} did not beat narrow {narrow}"); }}