solow_stats/
normality.rs

1//! Normality and serial-correlation diagnostics on a residual series.
2
3use ndarray::Array1;
4use solow_distributions::chi2_sf;
5
6/// Durbin–Watson statistic for first-order serial correlation in `resid`.
7///
8/// Defined as `sum((e_t - e_{t-1})^2) / sum(e_t^2)`; lies in `[0, 4]`, equals
9/// `2` under the null of no serial correlation.
10pub fn durbin_watson(resid: &Array1<f64>) -> f64 {
11    let n = resid.len();
12    if n < 2 {
13        return f64::NAN;
14    }
15    let mut num = 0.0;
16    let mut den = resid[0] * resid[0];
17    for i in 1..n {
18        let d = resid[i] - resid[i - 1];
19        num += d * d;
20        den += resid[i] * resid[i];
21    }
22    num / den
23}
24
25/// The `k`-th central moment of `x` about its mean (population / biased form,
26/// i.e. divided by `n`).
27fn central_moment(x: &Array1<f64>, mean: f64, k: i32) -> f64 {
28    let n = x.len() as f64;
29    x.iter().map(|&v| (v - mean).powi(k)).sum::<f64>() / n
30}
31
32/// Output of [`jarque_bera`].
33#[derive(Debug, Clone, Copy)]
34pub struct JarqueBera {
35    /// Jarque–Bera test statistic.
36    pub statistic: f64,
37    /// Two-sided p-value from the chi-squared distribution with 2 d.o.f.
38    pub pvalue: f64,
39    /// Sample skewness (biased estimator).
40    pub skew: f64,
41    /// Sample kurtosis (biased, non-excess: normal distribution gives `3`).
42    pub kurtosis: f64,
43}
44
45/// Jarque–Bera test of normality.
46///
47/// Returns the statistic `n·(S²/6 + (K−3)²/24)` where `S` is the (biased)
48/// sample skewness and `K` the (biased, non-excess) sample kurtosis, its
49/// chi-squared(2) p-value, and the underlying skewness and kurtosis.
50pub fn jarque_bera(resid: &Array1<f64>) -> JarqueBera {
51    let n = resid.len() as f64;
52    let mean = resid.sum() / n;
53    let m2 = central_moment(resid, mean, 2);
54    let m3 = central_moment(resid, mean, 3);
55    let m4 = central_moment(resid, mean, 4);
56    let skew = m3 / m2.powf(1.5);
57    let kurtosis = m4 / (m2 * m2); // non-excess kurtosis (== 3 + excess)
58    let jb = (n / 6.0) * (skew * skew + 0.25 * (kurtosis - 3.0).powi(2));
59    JarqueBera {
60        statistic: jb,
61        pvalue: chi2_sf(jb, 2.0),
62        skew,
63        kurtosis,
64    }
65}
66
67/// Z-score of D'Agostino's skewness test (transformed sample skewness).
68///
69/// Mirrors the reference `scipy.stats.skewtest` transform.
70fn skew_z(skew: f64, n: f64) -> f64 {
71    let y = skew * (((n + 1.0) * (n + 3.0)) / (6.0 * (n - 2.0))).sqrt();
72    let beta2 = 3.0 * (n * n + 27.0 * n - 70.0) * (n + 1.0) * (n + 3.0)
73        / ((n - 2.0) * (n + 5.0) * (n + 7.0) * (n + 9.0));
74    let w2 = -1.0 + (2.0 * (beta2 - 1.0)).sqrt();
75    let delta = 1.0 / (0.5 * w2.ln()).sqrt();
76    let alpha = (2.0 / (w2 - 1.0)).sqrt();
77    let y = if y == 0.0 { 1.0 } else { y };
78    let ya = y / alpha;
79    delta * (ya + (ya * ya + 1.0).sqrt()).ln()
80}
81
82/// Z-score of Anscombe–Glynn kurtosis test (transformed sample kurtosis).
83///
84/// Mirrors the reference `scipy.stats.kurtosistest` transform. `kurt` is the
85/// non-excess (Pearson) kurtosis.
86fn kurtosis_z(kurt: f64, n: f64) -> f64 {
87    let e = 3.0 * (n - 1.0) / (n + 1.0);
88    let varb2 = 24.0 * n * (n - 2.0) * (n - 3.0) / ((n + 1.0) * (n + 1.0) * (n + 3.0) * (n + 5.0));
89    let x = (kurt - e) / varb2.sqrt();
90    let sqrtbeta1 = 6.0 * (n * n - 5.0 * n + 2.0) / ((n + 7.0) * (n + 9.0))
91        * (6.0 * (n + 3.0) * (n + 5.0) / (n * (n - 2.0) * (n - 3.0))).sqrt();
92    let a =
93        6.0 + 8.0 / sqrtbeta1 * (2.0 / sqrtbeta1 + (1.0 + 4.0 / (sqrtbeta1 * sqrtbeta1)).sqrt());
94    let term1 = 1.0 - 2.0 / (9.0 * a);
95    let denom = 1.0 + x * (2.0 / (a - 4.0)).sqrt();
96    let term2 = denom.signum() * ((1.0 - 2.0 / a) / denom.abs()).powf(1.0 / 3.0);
97    (term1 - term2) / (2.0 / (9.0 * a)).sqrt()
98}
99
100/// D'Agostino–Pearson omnibus test of normality (the reference's
101/// `omni_normtest`).
102///
103/// Returns `(statistic, pvalue)` where the statistic is `Z_skew² + Z_kurt²`,
104/// asymptotically chi-squared with 2 d.o.f. Requires at least 8 observations;
105/// fewer yields `(NaN, NaN)`.
106pub fn omni_normtest(resid: &Array1<f64>) -> (f64, f64) {
107    let n = resid.len();
108    if n < 8 {
109        return (f64::NAN, f64::NAN);
110    }
111    let nf = n as f64;
112    let mean = resid.sum() / nf;
113    let m2 = central_moment(resid, mean, 2);
114    let m3 = central_moment(resid, mean, 3);
115    let m4 = central_moment(resid, mean, 4);
116    let skew = m3 / m2.powf(1.5);
117    let kurt = m4 / (m2 * m2);
118    let zs = skew_z(skew, nf);
119    let zk = kurtosis_z(kurt, nf);
120    let stat = zs * zs + zk * zk;
121    (stat, chi2_sf(stat, 2.0))
122}
123
124#[cfg(test)]
125mod tests {
126    use super::*;
127    use ndarray::array;
128
129    #[test]
130    fn dw_no_correlation_near_two() {
131        // Alternating residuals -> strong negative correlation -> dw near 4.
132        let r = array![1.0, -1.0, 1.0, -1.0, 1.0, -1.0];
133        let dw = durbin_watson(&r);
134        // Strong negative serial correlation pushes the statistic toward 4.
135        assert!(dw > 3.0, "dw = {dw}");
136    }
137
138    #[test]
139    fn jb_symmetric_low_skew() {
140        let r = array![-2.0, -1.0, 0.0, 1.0, 2.0];
141        let jb = jarque_bera(&r);
142        assert!(jb.skew.abs() < 1e-12);
143    }
144}
solow_stats/normality.rs

solow_stats/
normality.rs
