rustyqlib/core/montecarlo/
variance_reduction.rs

1//! Variance-reduction building blocks: moment matching and the generic
2//! regression-based control-variate estimator. (Antithetic pairing lives
3//! where draws are generated — see
4//! [`pseudo_normals`](super::rng::pseudo_normals) — and low-discrepancy
5//! sampling in [`sobol`](super::sobol) / [`halton`](super::halton) is
6//! itself the strongest variance reduction for smooth payoffs.)
7
8/// Rescale draws in place to sample mean 0 and variance 1 exactly —
9/// removes the O(1/sqrt(n)) noise in the first two sample moments.
10pub fn moment_match(draws: &mut [f64]) {
11    let n = draws.len() as f64;
12    let mean = draws.iter().sum::<f64>() / n;
13    let var = draws.iter().map(|z| (z - mean) * (z - mean)).sum::<f64>() / n;
14    let std = var.sqrt();
15    if std > 0.0 {
16        for z in draws.iter_mut() {
17            *z = (*z - mean) / std;
18        }
19    }
20}
21
22/// Control-variate estimate of `E[payoff]` given per-path values of a
23/// control with known expectation `control_mean`.
24///
25/// The optimal coefficient `beta = cov(payoff, control) / var(control)`
26/// is estimated from the same sample, and the adjusted estimator is
27/// `mean(payoff - beta (control - control_mean))`. Returns
28/// `(estimate, standard_error, beta)`; the standard error is that of the
29/// adjusted per-path values, so it shows the variance actually removed.
30pub fn control_variate_estimate(
31    payoffs: &[f64],
32    controls: &[f64],
33    control_mean: f64,
34) -> (f64, f64, f64) {
35    assert_eq!(payoffs.len(), controls.len());
36    let n = payoffs.len() as f64;
37    assert!(n > 1.0, "need at least two paths");
38    let mean_y = payoffs.iter().sum::<f64>() / n;
39    let mean_c = controls.iter().sum::<f64>() / n;
40    let mut cov = 0.0;
41    let mut var_c = 0.0;
42    for (y, c) in payoffs.iter().zip(controls) {
43        cov += (y - mean_y) * (c - mean_c);
44        var_c += (c - mean_c) * (c - mean_c);
45    }
46    let beta = if var_c > 0.0 { cov / var_c } else { 0.0 };
47
48    let mut sum = 0.0;
49    let mut sum_sq = 0.0;
50    for (y, c) in payoffs.iter().zip(controls) {
51        let adjusted = y - beta * (c - control_mean);
52        sum += adjusted;
53        sum_sq += adjusted * adjusted;
54    }
55    let mean = sum / n;
56    let var = (sum_sq / n - mean * mean).max(0.0);
57    (mean, (var / n).sqrt(), beta)
58}
59
60#[cfg(test)]
61mod tests {
62    use super::super::rng::path_normals;
63    use super::*;
64
65    #[test]
66    fn moment_matching_standardizes_exactly() {
67        let mut z = vec![0.0; 1001];
68        path_normals(3, 0, &mut z);
69        moment_match(&mut z);
70        let n = z.len() as f64;
71        let mean: f64 = z.iter().sum::<f64>() / n;
72        let var: f64 = z.iter().map(|x| x * x).sum::<f64>() / n;
73        assert!(mean.abs() < 1e-12 && (var - 1.0).abs() < 1e-12);
74    }
75
76    #[test]
77    fn control_variate_removes_correlated_noise() {
78        // payoff = 3 + 2 C + independent noise, control C with known mean 0:
79        // the fit must find beta ~ 2 and cut the standard error sharply
80        let n = 4000;
81        let mut c = vec![0.0; n];
82        let mut eps = vec![0.0; n];
83        path_normals(11, 0, &mut c);
84        path_normals(11, 1, &mut eps);
85        let payoffs: Vec<f64> =
86            c.iter().zip(&eps).map(|(ci, ei)| 3.0 + 2.0 * ci + 0.1 * ei).collect();
87
88        let (est, se, beta) = control_variate_estimate(&payoffs, &c, 0.0);
89        assert!((beta - 2.0).abs() < 0.05, "beta {beta}");
90        assert!((est - 3.0).abs() < 0.01, "estimate {est}");
91        // raw standard error of the payoffs for comparison
92        let mean_y: f64 = payoffs.iter().sum::<f64>() / n as f64;
93        let var_y: f64 =
94            payoffs.iter().map(|y| (y - mean_y) * (y - mean_y)).sum::<f64>() / n as f64;
95        let raw_se = (var_y / n as f64).sqrt();
96        assert!(se < 0.1 * raw_se, "cv se {se} vs raw {raw_se}");
97    }
98
99    #[test]
100    fn zero_variance_control_degrades_gracefully() {
101        let payoffs = [1.0, 2.0, 3.0];
102        let controls = [5.0, 5.0, 5.0];
103        let (est, _, beta) = control_variate_estimate(&payoffs, &controls, 5.0);
104        assert_eq!(beta, 0.0);
105        assert!((est - 2.0).abs() < 1e-12);
106    }
107}
rustyqlib/core/montecarlo/variance_reduction.rs

rustyqlib/core/montecarlo/
variance_reduction.rs
