quant/src/numerics.rs
The two ways a derivative actually gets priced, and what each costs.
//! The two ways a derivative actually gets priced, and what each costs.//!//! The numerical methods chapter's argument is that the models in these notes//! were selected by a computational constraint --- affine, Gaussian, tractable//! --- and that the constraint has largely lifted. Making that argument//! honestly means being precise about what the numerical methods cost, so this//! module measures it rather than asserting it.//!//! Three things are built here. A Monte Carlo, whose error falls as one over the//! square root of the path count and no faster. A binomial tree, which handles//! early exercise exactly by backward induction and is the reference the others//! are judged against. And Longstaff-Schwartz, which is how early exercise is//! done when the state is too large for a tree, and which is biased --- in a//! direction that depends on a detail people get wrong. use crate::black::{black_scholes, Side};use crate::pathwise::Rng; /// A Black-Scholes market, in the parameters every routine here shares.#[derive(Clone, Copy, Debug)]pub struct Market { pub spot: f64, pub strike: f64, pub rate: f64, pub vol: f64, pub expiry: f64,} impl Market { /// The European put price, in closed form. The thing a numerical method has /// to reproduce before it can be trusted on anything harder. pub fn european_put(&self) -> f64 { black_scholes(self.spot, self.strike, self.vol, self.expiry, self.rate, 0.0, Side::Put) }} /// A Monte Carlo price and its own estimate of how wrong it is.pub struct MonteCarlo { pub price: f64, /// Standard error: the sample standard deviation over the square root of the /// path count. Reported alongside the price because a Monte Carlo number /// without one is not a number, it is a draw. pub standard_error: f64,} /// Price a European put by simulation.////// Deliberately the plainest possible estimator --- no control variate, no/// antithetic sampling, no low-discrepancy sequence. The point of the figure it/// feeds is the convergence rate, and every one of those techniques improves the/// constant while leaving the rate alone.pub fn monte_carlo_put(m: &Market, paths: usize, seed: u64) -> MonteCarlo { let mut rng = Rng::new(seed); let drift = (m.rate - 0.5 * m.vol * m.vol) * m.expiry; let diffusion = m.vol * m.expiry.sqrt(); let discount = (-m.rate * m.expiry).exp(); let (mut total, mut total_sq) = (0.0, 0.0); for _ in 0..paths { let terminal = m.spot * (drift + diffusion * rng.next_normal()).exp(); let payoff = discount * (m.strike - terminal).max(0.0); total += payoff; total_sq += payoff * payoff; } let n = paths as f64; let mean = total / n; let variance = (total_sq / n - mean * mean).max(0.0); MonteCarlo { price: mean, standard_error: (variance / n).sqrt() }} /// A Bermudan put on a Cox-Ross-Rubinstein tree.////// Backward induction handles early exercise exactly: at each node the value is/// the larger of exercising and continuing, and continuing is known because the/// step after has already been solved. There is no estimation anywhere in it,/// which is what makes this the reference.////// `exercises` is the number of equally spaced exercise dates, the last at/// expiry. The step count is rounded up so every exercise date lands on a step.pub fn binomial_bermudan(m: &Market, exercises: usize, steps_per_exercise: usize) -> f64 { let steps = exercises * steps_per_exercise; let dt = m.expiry / steps as f64; let up = (m.vol * dt.sqrt()).exp(); let down = 1.0 / up; let discount = (-m.rate * dt).exp(); let p = ((m.rate * dt).exp() - down) / (up - down); // Terminal layer. let mut value: Vec<f64> = (0..=steps) .map(|j| { let s = m.spot * up.powi(j as i32) * down.powi((steps - j) as i32); (m.strike - s).max(0.0) }) .collect(); for step in (0..steps).rev() { let exercisable = (step + 1) % steps_per_exercise == 0 && step + 1 != steps; for j in 0..=step { let continuation = discount * (p * value[j + 1] + (1.0 - p) * value[j]); value[j] = if exercisable { let s = m.spot * up.powi(j as i32) * down.powi((step - j) as i32); continuation.max((m.strike - s).max(0.0)) } else { continuation }; } value.truncate(step + 1); } value[0]} /// Which Longstaff-Schwartz estimator to run.#[derive(Clone, Copy, Debug, PartialEq, Eq)]pub enum Lsm { /// Regress and exercise on the same paths. This is the estimator as it is /// usually written down, and it is biased *high*: the regression has seen /// the payoffs it is deciding about, so the exercise rule is fitted to the /// noise as well as the signal, and a rule fitted to noise looks better on /// the data it was fitted to than it is. InSample, /// Fit the rule on one set of paths and apply it to a fresh set. Any fixed /// exercise rule is suboptimal, and a suboptimal rule undervalues the /// option, so this is biased *low*. OutOfSample,} /// Price a Bermudan put by Longstaff-Schwartz.////// The method every desk uses when the state is too large for a tree: simulate/// forward, then step backwards estimating the continuation value by regressing/// discounted future cashflows on functions of the current state, and exercise/// where intrinsic beats the estimate.////// The regression is on `1, S, S^2` over the in-the-money paths only, which is/// the original recipe. Out-of-the-money paths carry no information about the/// exercise boundary and including them spends the fit on the wrong region.pub fn longstaff_schwartz( m: &Market, exercises: usize, paths: usize, mode: Lsm, seed: u64,) -> f64 { let dt = m.expiry / exercises as f64; let discount = (-m.rate * dt).exp(); let simulate = |seed: u64| -> Vec<Vec<f64>> { let mut rng = Rng::new(seed); let drift = (m.rate - 0.5 * m.vol * m.vol) * dt; let diffusion = m.vol * dt.sqrt(); (0..paths) .map(|_| { let mut s = m.spot; (0..exercises) .map(|_| { s *= (drift + diffusion * rng.next_normal()).exp(); s }) .collect() }) .collect() }; let fitting = simulate(seed); // The rule is fitted on `fitting` and, out of sample, applied to `pricing`. let pricing = match mode { Lsm::InSample => fitting.clone(), Lsm::OutOfSample => simulate(seed ^ 0x5DEE_CE66_D5B1_1F17), }; let intrinsic = |s: f64| (m.strike - s).max(0.0); // Cashflow carried backwards along each path, on both sets. let mut fit_value: Vec<f64> = fitting.iter().map(|p| intrinsic(p[exercises - 1])).collect(); let mut price_value: Vec<f64> = pricing.iter().map(|p| intrinsic(p[exercises - 1])).collect(); for step in (0..exercises - 1).rev() { for v in fit_value.iter_mut() { *v *= discount; } for v in price_value.iter_mut() { *v *= discount; } // Fit the continuation value on the in-the-money paths of the fitting set. let mut rows: Vec<[f64; 3]> = Vec::new(); let mut targets: Vec<f64> = Vec::new(); for (path, &carried) in fitting.iter().zip(&fit_value) { let s = path[step]; if intrinsic(s) > 0.0 { rows.push([1.0, s, s * s]); targets.push(carried); } } let beta = match least_squares(&rows, &targets) { Some(b) => b, // Nothing in the money at this date: nobody exercises, carry on. None => continue, }; let continuation = |s: f64| beta[0] + beta[1] * s + beta[2] * s * s; // Apply the rule. On the fitting set this is the in-sample estimator; on // a fresh set it is the out-of-sample one. let apply = |paths: &Vec<Vec<f64>>, values: &mut Vec<f64>| { for (path, v) in paths.iter().zip(values.iter_mut()) { let s = path[step]; let exercise = intrinsic(s); if exercise > 0.0 && exercise > continuation(s) { *v = exercise; } } }; apply(&fitting, &mut fit_value); if mode == Lsm::OutOfSample { apply(&pricing, &mut price_value); } } let values = match mode { Lsm::InSample => &fit_value, Lsm::OutOfSample => &price_value, }; discount * values.iter().sum::<f64>() / paths as f64} /// Least squares for three basis functions, by the normal equations.////// Three by three and symmetric, so Gaussian elimination is the whole solver./// Returns `None` when the system is singular, which happens when too few paths/// are in the money to determine a quadratic.fn least_squares(rows: &[[f64; 3]], targets: &[f64]) -> Option<[f64; 3]> { if rows.len() < 3 { return None; } let mut a = [[0.0f64; 4]; 3]; for (row, &t) in rows.iter().zip(targets) { for i in 0..3 { for j in 0..3 { a[i][j] += row[i] * row[j]; } a[i][3] += row[i] * t; } } for i in 0..3 { // Partial pivoting, since the columns differ by orders of magnitude when // the spot is far from one. let pivot = (i..3).max_by(|&x, &y| { a[x][i].abs().partial_cmp(&a[y][i].abs()).unwrap_or(std::cmp::Ordering::Equal) })?; a.swap(i, pivot); if a[i][i].abs() < 1e-12 { return None; } for k in (i + 1)..3 { let factor = a[k][i] / a[i][i]; for j in i..4 { a[k][j] -= factor * a[i][j]; } } } let mut x = [0.0; 3]; for i in (0..3).rev() { let mut sum = a[i][3]; for j in (i + 1)..3 { sum -= a[i][j] * x[j]; } x[i] = sum / a[i][i]; } Some(x)} #[cfg(test)]mod tests { use super::*; fn market() -> Market { Market { spot: 100.0, strike: 100.0, rate: 0.05, vol: 0.25, expiry: 1.0 } } #[test] fn monte_carlo_finds_the_closed_form() { let m = market(); let mc = monte_carlo_put(&m, 400_000, 20_260_806); let exact = m.european_put(); assert!( (mc.price - exact).abs() < 4.0 * mc.standard_error, "{} against {exact}, standard error {}", mc.price, mc.standard_error ); } #[test] fn monte_carlo_error_falls_as_one_over_root_n() { // The rate that decides everything about when Monte Carlo is affordable: // four times the work to halve the error. Measured rather than asserted, // because it is the basis of the comparison the chapter draws. let m = market(); let coarse = monte_carlo_put(&m, 25_000, 11_111).standard_error; let fine = monte_carlo_put(&m, 400_000, 11_111).standard_error; // Sixteen times the paths should be four times as accurate. let ratio = coarse / fine; assert!( (ratio - 4.0).abs() < 0.4, "sixteen times the paths gave {ratio} times the accuracy, not four" ); } #[test] fn the_tree_reproduces_the_closed_form_when_exercise_is_european() { // One exercise date, at expiry, is a European option -- so the tree can // be checked against Black-Scholes before being trusted as a reference // for anything it cannot be checked against. let m = market(); let tree = binomial_bermudan(&m, 1, 4000); assert!( (tree - m.european_put()).abs() < 0.005, "tree {tree} against closed form {}", m.european_put() ); } #[test] fn early_exercise_is_worth_something() { let m = market(); let european = binomial_bermudan(&m, 1, 2000); let bermudan = binomial_bermudan(&m, 50, 40); assert!(bermudan > european, "bermudan {bermudan} not above european {european}"); } /// Average the estimator over independent runs, since a single run cannot /// see its own bias: at these path counts the standard deviation of one /// estimate is several times the bias being looked for, and its sign /// changes with the seed. fn mean_bias(m: &Market, exercises: usize, paths: usize, mode: Lsm, reps: usize) -> (f64, f64) { let reference = binomial_bermudan(m, exercises, 400); let (mut sum, mut sum_sq) = (0.0, 0.0); for r in 0..reps { let seed = 1_000_003u64.wrapping_mul(r as u64 + 1); let e = longstaff_schwartz(m, exercises, paths, mode, seed) - reference; sum += e; sum_sq += e * e; } let n = reps as f64; let mean = sum / n; (mean, ((sum_sq / n - mean * mean) / n).sqrt()) } #[test] fn one_run_cannot_see_its_own_bias() { // Worth pinning first, because it is the practical point and it defeats // the obvious experiment. The spread of a single estimate is far wider // than the bias underneath it, so comparing one Longstaff-Schwartz run // against a tree says nothing about which way the method leans. let m = market(); let mut lo = f64::MAX; let mut hi = f64::MIN; for r in 0..12 { let v = longstaff_schwartz(&m, 10, 1_000, Lsm::InSample, 7919 * (r + 1)); lo = lo.min(v); hi = hi.max(v); } let (bias, _) = mean_bias(&m, 10, 1_000, Lsm::InSample, 60); assert!( hi - lo > 3.0 * bias.abs(), "run-to-run spread {} was not far wider than the bias {bias}", hi - lo ); } #[test] fn reusing_the_paths_biases_the_answer_high() { // The foresight bias. The regression has seen the payoffs it is deciding // about, so the exercise rule is fitted to the noise as well as to the // signal, and a rule fitted to noise flatters itself on the data it was // fitted to. let m = market(); let (bias, err) = mean_bias(&m, 10, 1_000, Lsm::InSample, 100); assert!( bias > 2.5 * err, "in-sample bias {bias} +/- {err} was not clearly positive" ); } #[test] fn a_fixed_rule_applied_to_fresh_paths_biases_it_low() { // The other bias, and the one with a different cause: any exercise rule // that is not the optimal stopping rule undervalues the option, and a // quadratic in the spot is not the optimal rule. let m = market(); for paths in [1_000usize, 4_000] { let (bias, err) = mean_bias(&m, 10, paths, Lsm::OutOfSample, 100); assert!( bias < -2.0 * err, "at {paths} paths the out-of-sample bias {bias} +/- {err} was not clearly negative" ); } } #[test] fn the_two_biases_behave_differently_as_paths_are_added() { // The distinction that matters when deciding what to spend on. Foresight // is a finite sample effect and dies as paths are added; the loss from a // poor exercise rule is the basis's fault and does not, so buying paths // fixes one and not the other. let m = market(); let (foresight_few, _) = mean_bias(&m, 10, 1_000, Lsm::InSample, 100); let (foresight_many, _) = mean_bias(&m, 10, 8_000, Lsm::InSample, 100); assert!( foresight_many < foresight_few, "the high bias did not decay: {foresight_few} then {foresight_many}" ); let (policy_few, _) = mean_bias(&m, 10, 1_000, Lsm::OutOfSample, 100); let (policy_many, err) = mean_bias(&m, 10, 8_000, Lsm::OutOfSample, 100); assert!( policy_many < -1.5 * err, "the low bias vanished with paths, which would mean the basis was rich \ enough after all: {policy_few} then {policy_many}" ); }} /// A controlled experiment in dimension.////// The numerical methods chapter claims that a grid pays for accuracy in dimensions and a/// simulation pays in variance. The `1/sqrt(N)` half of that is measured above./// This is the other half, and it needs an integral whose *dimension* can be/// varied while everything else about it is held still.////// The construction: drive a single lognormal with `d` independent normals/// combined into one,////// ```text/// Z = (z_1 + ... + z_d) / sqrt(d),/// F_T = F exp(-sigma^2 T / 2 + sigma sqrt(T) Z),/// ```////// so `Z` is standard normal for every `d`. The call on `F_T` therefore has the/// *same* Black-76 price at every dimension, and the payoff has the same/// distribution, so a Monte Carlo estimator has the same variance too. Only the/// dimension of the integral being done changes. Anything that then varies with/// `d` is a property of the method rather than of the problem.pub struct DimensionTest { pub forward: f64, pub strike: f64, pub vol: f64, pub maturity: f64, pub dimension: usize,} impl DimensionTest { /// The exact price, which does not depend on the dimension. pub fn exact(&self) -> f64 { crate::black::black76(self.forward, self.strike, self.vol, self.maturity, Side::Call) } fn payoff(&self, sum_of_normals: f64) -> f64 { let z = sum_of_normals / (self.dimension as f64).sqrt(); let terminal = self.forward * (-0.5 * self.vol * self.vol * self.maturity + self.vol * self.maturity.sqrt() * z) .exp(); (terminal - self.strike).max(0.0) } /// Monte Carlo, drawing `dimension` normals per path and using every one. /// /// Costs `paths * dimension` evaluations of the random number generator but /// only `paths` evaluations of the payoff, which is the count that matters /// for the comparison against quadrature. pub fn monte_carlo(&self, paths: usize, seed: u64) -> MonteCarlo { let mut rng = Rng::new(seed); let (mut total, mut total_sq) = (0.0, 0.0); for _ in 0..paths { let mut sum = 0.0; for _ in 0..self.dimension { sum += rng.next_normal(); } let payoff = self.payoff(sum); total += payoff; total_sq += payoff * payoff; } let n = paths as f64; let mean = total / n; let variance = (total_sq / n - mean * mean).max(0.0); MonteCarlo { price: mean, standard_error: (variance / n).sqrt() } } /// The same integral by a tensor product trapezoid rule against the Gaussian /// weight, `nodes_per_axis` points on each of the `dimension` axes. /// /// This is the deterministic method in its plainest form: lay a grid over /// the state, evaluate everywhere, weight and sum. It costs /// `nodes_per_axis^dimension` payoff evaluations, and that is the whole /// point of the exercise. /// /// The axes are truncated at eight standard deviations, where the Gaussian /// weight has less than `1e-15` of its mass left, so truncation is not what /// limits the answer. pub fn tensor_trapezoid(&self, nodes_per_axis: usize) -> f64 { const LIMIT: f64 = 8.0; let n = nodes_per_axis.max(2); let d = self.dimension; let h = 2.0 * LIMIT / (n as f64 - 1.0); let node: Vec<f64> = (0..n).map(|j| -LIMIT + h * j as f64).collect(); let weight: Vec<f64> = (0..n) .map(|j| { let edge = if j == 0 || j == n - 1 { 0.5 } else { 1.0 }; edge * h * crate::black::norm_pdf(node[j]) }) .collect(); let mut index = vec![0usize; d]; let mut total = 0.0; 'sweep: loop { let mut w = 1.0; let mut sum = 0.0; for &j in &index { w *= weight[j]; sum += node[j]; } total += w * self.payoff(sum); // Odometer over the multi-index: increment, carry, stop on overflow // of the last axis. let mut axis = 0; loop { if axis == d { break 'sweep; } index[axis] += 1; if index[axis] < n { break; } index[axis] = 0; axis += 1; } } total } /// The most nodes per axis a tensor grid can afford within `budget` payoff /// evaluations. Two is the floor, since a trapezoid needs both ends. pub fn nodes_for_budget(&self, budget: usize) -> usize { let mut n = 2usize; while ((n + 1) as f64).powi(self.dimension as i32) <= budget as f64 { n += 1; } n }} #[cfg(test)]mod dimension_tests { use super::*; fn test(dimension: usize) -> DimensionTest { DimensionTest { forward: 100.0, strike: 100.0, vol: 0.2, maturity: 1.0, dimension } } const BUDGET: usize = 100_000; const SEED: u64 = 20260807; #[test] fn the_quadrature_is_right_where_it_can_be_checked() { // Before comparing methods, confirm the deterministic one is correct at // all. In one dimension it has the whole budget on a single axis and // should reproduce Black-76 to several decimals. let t = test(1); assert!((t.tensor_trapezoid(BUDGET) - t.exact()).abs() < 1e-5); } #[test] fn monte_carlo_does_not_notice_the_dimension() { // The claim, isolated. The problem is built so that the payoff has the // same distribution at every dimension, so anything the estimator does // differently is the estimator's doing --- and it does nothing // differently. Same standard error to three figures across four // doublings. let base = test(1).monte_carlo(BUDGET, SEED).standard_error; for d in [2, 4, 8, 16] { let t = test(d); let mc = t.monte_carlo(BUDGET, SEED); assert!( (mc.standard_error / base - 1.0).abs() < 0.02, "d={d}: standard error {} against {base} in one dimension", mc.standard_error ); assert!( (mc.price - t.exact()).abs() < 4.0 * mc.standard_error, "d={d}: priced {} against exact {}", mc.price, t.exact() ); } } #[test] fn a_tensor_grid_falls_off_a_cliff_around_five_dimensions() { // The other half. Same budget of payoff evaluations, spread over a // tensor grid instead of drawn at random. // // Low dimensions: the grid is not merely better, it is in a different // league, because a deterministic rule converges in the mesh rather than // in the square root of a sample. for d in [1, 2] { let t = test(d); let error = (t.tensor_trapezoid(t.nodes_for_budget(BUDGET)) - t.exact()).abs(); assert!(error < 1e-3, "d={d} gave {error}"); assert!(error < t.monte_carlo(BUDGET, SEED).standard_error); } // High dimensions: the same budget buys five or six nodes an axis, and // the answer is not wrong by a rounding error, it is wrong by most of // the price. for d in [6, 7] { let t = test(d); let error = (t.tensor_trapezoid(t.nodes_for_budget(BUDGET)) - t.exact()).abs(); assert!(error > 1.0, "d={d} gave {error}, expected a collapse"); assert!(error > 0.5 * t.exact(), "d={d}: {error} against a price of {}", t.exact()); } } #[test] fn the_crossover_is_where_the_nodes_per_axis_run_out() { // Why the cliff is where it is, rather than anywhere else. It is set by // the budget's d-th root: a fixed number of evaluations spread over more // axes leaves too few points on each to resolve anything. let expected = [100_000, 316, 46, 17, 10, 6, 5]; for (d, &nodes) in (1..=7).zip(expected.iter()) { assert_eq!(test(d).nodes_for_budget(BUDGET), nodes, "at d={d}"); } }} /// Whether a set of simulated prices is free of static arbitrage.////// The numerical methods chapter asks whether arbitrage-free pricing survives/// without a solvable model. The answer turns on where the guarantee comes from,/// and this isolates it.////// A call payoff is convex in the strike, so for any single path////// ```text/// (S-K1)^+ - 2 (S-K2)^+ + (S-K3)^+ >= 0, K2 = (K1+K3)/2,/// ```////// and averaging over a *common* set of paths preserves the inequality term by/// term. Price the three strikes on *different* paths and nothing preserves it:/// each price carries its own sampling error, the butterfly is a difference of/// nearly equal numbers, and the error is of the same order as the quantity.////// So no-arbitrage under simulation is not automatic and not a property of the/// model alone. It is a property of the estimator, and common random numbers are/// what supply it.////// Returns the fraction of strike triples whose butterfly comes out negative.pub fn butterfly_violations( m: &Market, width: f64, triples: usize, paths: usize, common_paths: bool, seed: u64,) -> f64 { let mut rng = Rng::new(seed); let drift = (m.rate - 0.5 * m.vol * m.vol) * m.expiry; let diffusion = m.vol * m.expiry.sqrt(); let discount = (-m.rate * m.expiry).exp(); let shared: Vec<f64> = (0..paths) .map(|_| m.spot * (drift + diffusion * rng.next_normal()).exp()) .collect(); let mut violations = 0usize; for i in 0..triples { let centre = m.strike * (0.7 + 0.6 * i as f64 / triples.max(2) as f64); let strikes = [centre - width, centre, centre + width]; let mut price = [0.0f64; 3]; for (j, &k) in strikes.iter().enumerate() { if common_paths { price[j] = discount * shared.iter().map(|s| (s - k).max(0.0)).sum::<f64>() / paths as f64; } else { let mut total = 0.0; for _ in 0..paths { let s = m.spot * (drift + diffusion * rng.next_normal()).exp(); total += (s - k).max(0.0); } price[j] = discount * total / paths as f64; } } if price[0] - 2.0 * price[1] + price[2] < 0.0 { violations += 1; } } violations as f64 / triples as f64} #[cfg(test)]mod arbitrage_under_simulation_tests { use super::*; fn market() -> Market { Market { spot: 100.0, strike: 100.0, rate: 0.02, vol: 0.2, expiry: 1.0 } } #[test] fn common_paths_make_simulated_prices_arbitrage_free() { // Priced on shared paths the butterfly is a sum of non-negative per-path // terms, so it cannot be negative however few paths are used. Ten // thousand is far too few for accurate prices and is still enough for // exact convexity: the guarantee is structural, not statistical. let m = market(); for width in [0.5, 2.0, 5.0] { let rate = butterfly_violations(&m, width, 40, 10_000, true, 20260808); assert_eq!(rate, 0.0, "width={width}: common paths should never violate"); } } #[test] fn independent_paths_do_not() { // The same model and the same path count, priced one option at a time. A // narrow butterfly is a difference of nearly equal numbers, so sampling // error swamps it and the price set admits a static arbitrage. let m = market(); let narrow = butterfly_violations(&m, 0.5, 200, 10_000, false, 20260808); assert!(narrow > 0.25, "narrow butterflies should invert often, got {narrow:.3}"); // Widening makes the true value large against the noise, so violations // recede -- the failure is worst exactly where a desk cares, at adjacent // strikes. let wide = butterfly_violations(&m, 8.0, 200, 10_000, false, 20260808); assert!(wide < narrow, "wide {wide:.3} should beat narrow {narrow:.3}"); } #[test] fn more_paths_do_not_repair_it() { // And the fix is not to work harder. Sampling error falls as one over the // square root of the paths while a butterfly's value falls as the square // of the width, so for any budget there is a width at which the // violations return. The discipline has to be the estimator. let m = market(); for (paths, width) in [(10_000usize, 0.5f64), (160_000, 0.125)] { let rate = butterfly_violations(&m, width, 200, paths, false, 20260808); assert!( rate > 0.25, "paths={paths} width={width}: expected frequent violations, got {rate:.3}" ); } }}