solow_stats/
diagnostic.rs

1//! Heteroscedasticity and serial-correlation diagnostic tests.
2
3use ndarray::{Array1, Array2};
4use solow_core::error::Result;
5use solow_distributions::chi2_sf;
6use solow_regression::LinearModel;
7
8/// Breusch–Pagan Lagrange-multiplier test for heteroscedasticity.
9///
10/// Regresses the squared residuals on `exog_het` (which must contain a
11/// constant) and reports the Koenker (studentized, robust) LM statistic
12/// `nobs · R²` with its chi-squared p-value (`k − 1` d.o.f.), together with the
13/// auxiliary regression's F statistic and its p-value. Returns
14/// `(lm, lm_pvalue, fvalue, f_pvalue)`. Mirrors the reference
15/// `het_breuschpagan` with `robust=True`.
16pub fn het_breuschpagan(
17    resid: &Array1<f64>,
18    exog_het: &Array2<f64>,
19) -> Result<(f64, f64, f64, f64)> {
20    let (nobs, nvars) = exog_het.dim();
21    let y: Array1<f64> = resid.mapv(|r| r * r);
22    let res = LinearModel::ols(y, exog_het.clone())?.fit()?;
23    let lm = nobs as f64 * res.rsquared;
24    let lm_pvalue = chi2_sf(lm, (nvars - 1) as f64);
25    Ok((lm, lm_pvalue, res.fvalue, res.f_pvalue))
26}
27
28/// White's Lagrange-multiplier test for heteroscedasticity.
29///
30/// Builds the auxiliary design from all squares and pairwise cross-products of
31/// the columns of `exog` (which must contain a constant), regresses the squared
32/// residuals on it, and reports the LM statistic `nobs · R²` with a chi-squared
33/// p-value using the auxiliary model's degrees of freedom (`rank − 1`),
34/// together with that model's F statistic and p-value. Returns
35/// `(lm, lm_pvalue, fvalue, f_pvalue)`. Mirrors the reference `het_white`.
36pub fn het_white(resid: &Array1<f64>, exog: &Array2<f64>) -> Result<(f64, f64, f64, f64)> {
37    let (nobs, k) = exog.dim();
38    // Upper-triangular index pairs (i <= j), matching numpy's triu_indices.
39    let mut pairs: Vec<(usize, usize)> = Vec::new();
40    for i in 0..k {
41        for j in i..k {
42            pairs.push((i, j));
43        }
44    }
45    let mut aux = Array2::<f64>::zeros((nobs, pairs.len()));
46    for (col, &(i, j)) in pairs.iter().enumerate() {
47        for row in 0..nobs {
48            aux[[row, col]] = exog[[row, i]] * exog[[row, j]];
49        }
50    }
51    let y: Array1<f64> = resid.mapv(|r| r * r);
52    let res = LinearModel::ols(y, aux)?.fit()?;
53    let lm = nobs as f64 * res.rsquared;
54    // Degrees of freedom take a possibly reduced rank into account: rank - 1.
55    let df = res.df_model;
56    let lm_pvalue = chi2_sf(lm, df);
57    Ok((lm, lm_pvalue, res.fvalue, res.f_pvalue))
58}
59
60/// Per-lag output row of [`acorr_ljungbox`].
61#[derive(Debug, Clone, Copy)]
62pub struct LjungBox {
63    /// The lag (1-based) this row reports.
64    pub lag: usize,
65    /// Ljung–Box cumulative test statistic up to this lag.
66    pub lb_stat: f64,
67    /// Chi-squared p-value of `lb_stat` with `lag` degrees of freedom.
68    pub lb_pvalue: f64,
69}
70
71/// Sample autocorrelation function of `x` for lags `0..=maxlag`.
72///
73/// Uses the biased estimator: the data are demeaned and each autocovariance is
74/// divided by `n`, matching the reference `acf` (via `acovf`).
75fn acf(x: &Array1<f64>, maxlag: usize) -> Vec<f64> {
76    let n = x.len();
77    let mean = x.sum() / n as f64;
78    let xo: Vec<f64> = x.iter().map(|&v| v - mean).collect();
79    let mut acov = vec![0.0; maxlag + 1];
80    for (lag, a) in acov.iter_mut().enumerate() {
81        let mut s = 0.0;
82        for t in lag..n {
83            s += xo[t] * xo[t - lag];
84        }
85        *a = s / n as f64;
86    }
87    let a0 = acov[0];
88    acov.iter().map(|&c| c / a0).collect()
89}
90
91/// Ljung–Box test of autocorrelation in the series `x`.
92///
93/// Computes, for every lag `1..=lags`, the cumulative Ljung–Box statistic
94/// `n(n+2) Σ_{k=1}^{lag} ρ_k² / (n−k)` and its chi-squared p-value with `lag`
95/// degrees of freedom (`model_df = 0`). Returns one [`LjungBox`] row per lag.
96/// Mirrors the reference `acorr_ljungbox` with default `model_df` and
97/// `boxpierce=False`.
98pub fn acorr_ljungbox(x: &Array1<f64>, lags: usize) -> Vec<LjungBox> {
99    let nobs = x.len();
100    let n = nobs as f64;
101    let sacf = acf(x, lags);
102    let mut cum = 0.0;
103    let mut out = Vec::with_capacity(lags);
104    for (lag, &rho) in sacf.iter().enumerate().take(lags + 1).skip(1) {
105        cum += rho * rho / (n - lag as f64);
106        let lb = n * (n + 2.0) * cum;
107        out.push(LjungBox {
108            lag,
109            lb_stat: lb,
110            lb_pvalue: chi2_sf(lb, lag as f64),
111        });
112    }
113    out
114}
115
116#[cfg(test)]
117mod tests {
118    use super::*;
119    use ndarray::array;
120
121    #[test]
122    fn ljungbox_lag0_count() {
123        let x = array![0.5, -0.2, 0.1, 0.3, -0.4, 0.2, 0.0, -0.1, 0.25, -0.15];
124        let res = acorr_ljungbox(&x, 3);
125        assert_eq!(res.len(), 3);
126        assert_eq!(res[0].lag, 1);
127        for r in &res {
128            assert!(r.lb_stat >= 0.0);
129            assert!((0.0..=1.0).contains(&r.lb_pvalue));
130        }
131    }
132
133    #[test]
134    fn acf_lag0_is_one() {
135        let x = array![1.0, 2.0, 3.0, 2.0, 1.0, 0.5];
136        let a = acf(&x, 2);
137        assert!((a[0] - 1.0).abs() < 1e-12);
138    }
139}
solow_stats/diagnostic.rs

solow_stats/
diagnostic.rs
