gam_math/
serial_dependence.rs

1//! Dependence-corrected summaries of a serially correlated sample.
2//!
3//! Both statistics read a sequence `x_1..x_n` whose terms may be autocorrelated
4//! (held-out per-row losses in row order, chain draws) and correct the naive
5//! i.i.d. summary for that dependence with the lag window `L = ⌊√n⌋`, which
6//! grows with `n` while its share `L/n` vanishes — the standard consistent
7//! bandwidth for a sample of unknown correlation length.
8
9/// Effective sample size `n / (1 + 2 Σ_{k≥1} ρ_k)` from the initial positive
10/// sequence of sample autocorrelations `ρ_k`, truncated at the first
11/// non-positive lag (Geyer's rule) and at the lag window `⌊√n⌋`. Returns `n`
12/// itself for a degenerate or constant sample and never less than one.
13pub fn autocorr_ess(x: &[f64]) -> f64 {
14    let n = x.len();
15    if n <= 1 {
16        return n as f64;
17    }
18    let mean = x.iter().sum::<f64>() / n as f64;
19    let var = x.iter().map(|v| (v - mean) * (v - mean)).sum::<f64>() / n as f64;
20    if var <= 0.0 {
21        return n as f64;
22    }
23    let lag_cap = (n as f64).sqrt() as usize;
24    let mut rho_sum = 0.0;
25    for lag in 1..=lag_cap.max(1).min(n - 1) {
26        let mut cov = 0.0;
27        for i in lag..n {
28            cov += (x[i] - mean) * (x[i - lag] - mean);
29        }
30        cov /= (n - lag) as f64;
31        let rho = cov / var;
32        if rho <= 0.0 || !rho.is_finite() {
33            break;
34        }
35        rho_sum += rho;
36    }
37    (n as f64 / (1.0 + 2.0 * rho_sum)).max(1.0)
38}
39
40/// Newey–West (Bartlett-kernel) standard error of the sample mean with lag
41/// window `⌊√n⌋`: `√(γ_0 + 2 Σ_k w_k γ_k) / √n`, `w_k = 1 − k/(L+1)`. Infinite
42/// for a sample of fewer than two terms, where no dispersion is measurable.
43pub fn newey_west_se(x: &[f64]) -> f64 {
44    let n = x.len();
45    if n <= 1 {
46        return f64::INFINITY;
47    }
48    let mean = x.iter().sum::<f64>() / n as f64;
49    let lag_cap = (n as f64).sqrt() as usize;
50    let mut gamma0 = 0.0;
51    for v in x {
52        gamma0 += (v - mean) * (v - mean);
53    }
54    gamma0 /= n as f64;
55    let mut var = gamma0;
56    for lag in 1..=lag_cap.max(1).min(n - 1) {
57        let mut gamma = 0.0;
58        for i in lag..n {
59            gamma += (x[i] - mean) * (x[i - lag] - mean);
60        }
61        gamma /= n as f64;
62        let w = 1.0 - lag as f64 / (lag_cap as f64 + 1.0);
63        var += 2.0 * w * gamma;
64    }
65    (var.max(0.0) / n as f64).sqrt()
66}
gam_math/serial_dependence.rs

gam_math/
serial_dependence.rs
