chronos_ts/
stat_tests.rs

1#![allow(non_snake_case)]
2use crate::linalg;
3use ndarray::{s, Array1, Array2};
4
5/// Augmented Dickey-Fuller (ADF) Test
6/// Null Hypothesis (H0): The series has a unit root (is non-stationary).
7/// Alternative (H1): The series is stationary.
8pub struct AdfTestResult {
9    pub stat: f64,
10    pub p_value: f64,
11    pub used_lags: usize,
12}
13
14pub fn adf_test(series: &Array1<f64>, max_lags: Option<usize>) -> AdfTestResult {
15    let n = series.len();
16
17    // 1. Minimum length guard: need at least 4-5 points to fit even a zero-lag regression
18    if n < 4 {
19        return AdfTestResult {
20            stat: 0.0,
21            p_value: 1.0,
22            used_lags: 0,
23        };
24    }
25
26    // 2. Zero-variance / constant series guard
27    let std_dev = series.std(0.0);
28    if std_dev.abs() < f64::EPSILON || std_dev.is_nan() {
29        return AdfTestResult {
30            stat: 0.0,
31            p_value: 0.0, // Constant series has no unit root (stationary)
32            used_lags: 0,
33        };
34    }
35
36    let lags = max_lags.unwrap_or_else(|| ((n as f64 - 1.0).powf(1.0 / 3.0)) as usize);
37
38    // Delta Y_t = Y_t - Y_{t-1}
39    let dy = &series.slice(s![1..]) - &series.slice(s![..-1]);
40    let dy_len = dy.len();
41
42    if dy_len <= lags {
43        return AdfTestResult {
44            stat: 0.0,
45            p_value: 1.0,
46            used_lags: lags,
47        };
48    }
49
50    let effective_n = dy_len - lags;
51    let cols = 2 + lags;
52
53    // Ensure degrees of freedom are positive
54    if effective_n <= cols {
55        return AdfTestResult {
56            stat: 0.0,
57            p_value: 1.0,
58            used_lags: lags,
59        };
60    }
61
62    // Dependent variable: Delta Y_t from t = lags to dy_len
63    let y_dep = dy.slice(s![lags..dy_len]).to_owned();
64
65    // Design Matrix X: [Y_{t-1}, 1 (constant), Delta Y_{t-1}, ..., Delta Y_{t-lags}]
66    let mut x = Array2::<f64>::zeros((effective_n, cols));
67
68    for i in 0..effective_n {
69        let idx = i + lags;
70        x[[i, 0]] = series[idx]; // Level term Y_{t-1}
71        x[[i, 1]] = 1.0; // Intercept term
72
73        for j in 0..lags {
74            x[[i, 2 + j]] = dy[idx - 1 - j]; // Lagged differences
75        }
76    }
77
78    // OLS via the normal equations (pure Rust, no LAPACK backend required).
79    let beta = match linalg::lstsq(&x, &y_dep) {
80        Ok(b) => b,
81        Err(_) => {
82            return AdfTestResult {
83                stat: 0.0,
84                p_value: 1.0,
85                used_lags: lags,
86            };
87        }
88    };
89
90    // Compute standard error of gamma (beta[0])
91    let residuals = &y_dep - &x.dot(&beta);
92    let sse = residuals.iter().map(|r| r.powi(2)).sum::<f64>();
93    let df = effective_n - cols;
94    let mse = sse / (df as f64);
95
96    let xtx_inv = match linalg::inv(&x.t().dot(&x)) {
97        Ok(m) => m,
98        Err(_) => {
99            return AdfTestResult {
100                stat: 0.0,
101                p_value: 1.0,
102                used_lags: lags,
103            };
104        }
105    };
106
107    let se_gamma = (mse * xtx_inv[[0, 0]]).sqrt();
108
109    if se_gamma <= 0.0 || se_gamma.is_nan() {
110        return AdfTestResult {
111            stat: 0.0,
112            p_value: 0.0,
113            used_lags: lags,
114        };
115    }
116
117    let t_stat = beta[0] / se_gamma;
118
119    // The ADF statistic does NOT follow a Student-t distribution under the null;
120    // it follows the (non-standard) Dickey-Fuller distribution. Map the statistic
121    // to a p-value using tabulated DF quantiles for the constant-only case.
122    let p_value = dickey_fuller_pvalue(t_stat);
123
124    AdfTestResult {
125        stat: t_stat,
126        p_value,
127        used_lags: lags,
128    }
129}
130
131/// Approximate p-value of the Augmented Dickey-Fuller statistic for the
132/// "constant, no trend" regression case.
133///
134/// The ADF `tau` statistic follows the Dickey-Fuller distribution rather than a
135/// Student-t. This uses the well-established asymptotic quantiles of that
136/// distribution (Fuller 1976 / MacKinnon) and performs monotone linear
137/// interpolation of the CDF, which is accurate around the decision region
138/// (~1%-10%) that matters for differencing decisions. Values outside the table
139/// are clamped to `[0, 1]`.
140fn dickey_fuller_pvalue(tau: f64) -> f64 {
141    // (tau quantile, cumulative probability) pairs, ascending in tau.
142    // Left tail => small p (reject unit root / stationary).
143    const TABLE: [(f64, f64); 8] = [
144        (-3.43, 0.01),
145        (-3.12, 0.025),
146        (-2.86, 0.05),
147        (-2.57, 0.10),
148        (-0.44, 0.90),
149        (-0.07, 0.95),
150        (0.23, 0.975),
151        (0.60, 0.99),
152    ];
153
154    if tau <= TABLE[0].0 {
155        return 0.01;
156    }
157    let last = TABLE[TABLE.len() - 1];
158    if tau >= last.0 {
159        return 0.99;
160    }
161
162    for w in TABLE.windows(2) {
163        let (t0, p0) = w[0];
164        let (t1, p1) = w[1];
165        if tau >= t0 && tau <= t1 {
166            let frac = (tau - t0) / (t1 - t0);
167            return (p0 + frac * (p1 - p0)).clamp(0.0, 1.0);
168        }
169    }
170
171    // Unreachable given the bounds checks above, but stay safe.
172    1.0
173}
174
175/// Automatically determines required non-seasonal differencing d
176pub fn estimate_d(series: &Array1<f64>, max_d: usize, alpha: f64) -> usize {
177    let mut current = series.clone();
178    let mut d = 0;
179
180    while d < max_d {
181        let res = adf_test(&current, None);
182        if res.p_value < alpha {
183            // Reject H0 -> Series is stationary
184            break;
185        }
186        // Fail to reject H0 -> Need differencing
187        if current.len() <= 2 {
188            break;
189        }
190        current = &current.slice(s![1..]) - &current.slice(s![..-1]);
191        d += 1;
192    }
193    d
194}
195
196/// Estimates the required seasonal differencing order `D` using a seasonal-strength
197/// heuristic (not a formal OCSB / Canova-Hansen unit-root test).
198///
199/// At each step it seasonally differences the series and measures how much variance
200/// that removes: `F_s = max(0, 1 - Var(seasonally differenced) / Var(current))`. If
201/// `F_s` exceeds a fixed threshold (0.64) the seasonal component is deemed strong
202/// enough to warrant another seasonal difference. This mirrors the strength-based
203/// rule popularised by Wang, Smith & Hyndman and is cheaper than a full unit-root
204/// test, at the cost of some statistical rigor.
205pub fn estimate_D(series: &Array1<f64>, m: usize, max_D: usize) -> usize {
206    if m <= 1 || series.len() < 2 * m {
207        return 0;
208    }
209
210    let mut current = series.clone();
211    let mut D = 0;
212
213    while D < max_D {
214        let n = current.len();
215        if n <= 2 * m {
216            break;
217        }
218
219        // Calculate seasonal strength index
220        let m_neg = -(m as isize);
221        let seasonal_diff = &current.slice(s![m..]) - &current.slice(s![..m_neg]);
222        let var_orig = crate::utils::variance(&current);
223        let var_sdiff = crate::utils::variance(&seasonal_diff);
224
225        // Seasonal strength F_s = max(0, 1 - Var(res) / Var(res + seasonal))
226        let seasonal_strength = (1.0 - (var_sdiff / var_orig)).max(0.0);
227
228        if seasonal_strength < 0.64 {
229            break;
230        }
231
232        current = seasonal_diff;
233        D += 1;
234    }
235    D
236}
chronos_ts/stat_tests.rs

chronos_ts/
stat_tests.rs
