tsorder_core/
integration.rs

1//! Order of integration `d` by sequential unit-root testing with a vanishing
2//! level: difference until the test stops indicating a unit root, with
3//! `α_n → 0` so the number of differences is estimated consistently.
4//!
5//! Two tests are offered. KPSS (default; null = stationarity): under a true
6//! stationary level the rejection probability tends to `α_n → 0`, under a
7//! unit root the statistic diverges, so `d̂ → d` either way. ADF (null =
8//! unit root): the type-I error vanishes and the power tends to one, same
9//! conclusion from the other side. Fractional `d` is out of scope.
10
11use crate::criteria::{ic_order, Criterion};
12use crate::dist::{df_tau_c_critical, kpss_critical};
13use crate::linalg::ols;
14use crate::{check_input, Error, Result};
15
16/// Which unit-root test drives the sequence.
17#[derive(Debug, Clone, Copy, PartialEq, Eq)]
18pub enum UnitRootTest {
19    /// Kwiatkowski–Phillips–Schmidt–Shin level statistic with Bartlett/
20    /// Newey–West long-run variance, bandwidth `⌊4 (T/100)^{1/4}⌋`.
21    Kpss,
22    /// Augmented Dickey–Fuller with constant; augmentation lag chosen by BIC
23    /// on the differenced series (the "lag by tsorder" rule).
24    Adf,
25}
26
27/// Configuration of [`integration_order`].
28#[derive(Debug, Clone, Copy)]
29pub struct IntegrationConfig {
30    /// Largest `d` considered (default 2).
31    pub max_d: usize,
32    /// Level; `None` → `0.05·log 50 / log n` (the same vanishing-level rule
33    /// as the rank tests; satisfies Hosoya's consistency condition).
34    pub alpha: Option<f64>,
35    pub test: UnitRootTest,
36}
37
38impl Default for IntegrationConfig {
39    fn default() -> Self {
40        IntegrationConfig {
41            max_d: 2,
42            alpha: None,
43            test: UnitRootTest::Kpss,
44        }
45    }
46}
47
48/// Output of [`integration_order`].
49#[derive(Debug, Clone)]
50pub struct IntegrationResult {
51    /// Estimated order of integration (`max_d` if every stage still indicates
52    /// a unit root — inspect `stats` in that case).
53    pub d: usize,
54    /// Test statistic at each stage `0..=d` actually run.
55    pub stats: Vec<f64>,
56    /// Critical value used at each stage.
57    pub crits: Vec<f64>,
58    pub alpha: f64,
59    pub test: UnitRootTest,
60}
61
62/// Estimate the order of integration by sequential testing with a vanishing
63/// level. The ARMA order selectors in this crate assume stationarity; run
64/// this first and difference `d` times.
65pub fn integration_order(x: &[f64], config: IntegrationConfig) -> Result<IntegrationResult> {
66    let n = x.len();
67    check_input(x, 20 + 4 * config.max_d)?;
68    let alpha = config
69        .alpha
70        .unwrap_or(0.05 * 50f64.ln() / (n as f64).ln().max(1.0));
71    let mut series: Vec<f64> = x.to_vec();
72    let mut stats = Vec::new();
73    let mut crits = Vec::new();
74    let mut d = config.max_d;
75    for stage in 0..=config.max_d {
76        let (stat, crit, unit_root) = match config.test {
77            UnitRootTest::Kpss => {
78                let s = kpss_stat(&series)?;
79                let c = kpss_critical(alpha);
80                (s, c, s > c) // rejecting stationarity = unit root indicated
81            }
82            UnitRootTest::Adf => {
83                let s = adf_stat(&series)?;
84                let c = df_tau_c_critical(alpha);
85                (s, c, s >= c) // failing to reject the unit root
86            }
87        };
88        stats.push(stat);
89        crits.push(crit);
90        if !unit_root {
91            d = stage;
92            break;
93        }
94        if stage < config.max_d {
95            series = series.windows(2).map(|w| w[1] - w[0]).collect();
96        }
97    }
98    Ok(IntegrationResult {
99        d,
100        stats,
101        crits,
102        alpha,
103        test: config.test,
104    })
105}
106
107/// KPSS level statistic: `T^{-2} Σ S_t² / σ̂²_lr`, Bartlett long-run variance.
108fn kpss_stat(x: &[f64]) -> Result<f64> {
109    let t = x.len();
110    let mean = x.iter().sum::<f64>() / t as f64;
111    let e: Vec<f64> = x.iter().map(|v| v - mean).collect();
112    let mut s = 0.0;
113    let mut sum_s2 = 0.0;
114    for v in &e {
115        s += v;
116        sum_s2 += s * s;
117    }
118    let bandwidth = (4.0 * (t as f64 / 100.0).powf(0.25)).floor() as usize;
119    let mut lr = e.iter().map(|v| v * v).sum::<f64>() / t as f64;
120    for j in 1..=bandwidth.min(t - 1) {
121        let w = 1.0 - j as f64 / (bandwidth as f64 + 1.0);
122        let g = (0..t - j).map(|i| e[i] * e[i + j]).sum::<f64>() / t as f64;
123        lr += 2.0 * w * g;
124    }
125    if lr <= 0.0 {
126        return Err(Error::Singular);
127    }
128    Ok(sum_s2 / (t as f64 * t as f64 * lr))
129}
130
131/// ADF t-statistic (with constant): `Δx_t = α + ρ x_{t−1} + Σ δ_j Δx_{t−j}`,
132/// `k` lags by BIC on the differenced series.
133fn adf_stat(x: &[f64]) -> Result<f64> {
134    let t = x.len();
135    let dx: Vec<f64> = x.windows(2).map(|w| w[1] - w[0]).collect();
136    let k = ic_order(&dx, Criterion::Bic, None)?.order;
137    let start = k + 1; // Δx index offset
138    let rows = dx.len() - start;
139    let cols = 2 + k;
140    if rows <= cols + 2 {
141        return Err(Error::TooShort {
142            n: t,
143            needed: cols + start + 3,
144        });
145    }
146    let mut design = vec![0.0; rows * cols];
147    let mut y = vec![0.0; rows];
148    for r in 0..rows {
149        let ti = start + r; // index into dx; level index = ti (x[ti] is the lagged level)
150        y[r] = dx[ti];
151        design[r * cols] = 1.0;
152        design[r * cols + 1] = x[ti];
153        for j in 0..k {
154            design[r * cols + 2 + j] = dx[ti - j - 1];
155        }
156    }
157    let (beta, _s2, se) = ols(&design, rows, cols, &y)?;
158    Ok(beta[1] / se[1])
159}
tsorder_core/integration.rs

tsorder_core/
integration.rs
