tsorder_core/
levinson.rs

1//! Levinson–Durbin recursion: Yule–Walker AR fits of every order `1..=K` in
2//! `O(K²)`, yielding the PACF and the innovation variance of each order.
3//!
4//! This single recursion feeds *both* estimator families: the information
5//! criteria need `sigma2[k]` for every `k`; the Fine–Hwang estimators need the
6//! PACF (`reflection[k]`) or the long-AR coefficients (`phi_last`).
7
8use crate::{check_input, Error, Result};
9
10/// Output of [`levinson`].
11#[derive(Debug, Clone)]
12pub struct Levinson {
13    /// `reflection[k-1]` = partial autocorrelation at lag `k`, `k = 1..=K`.
14    pub pacf: Vec<f64>,
15    /// `sigma2[k]` = Yule–Walker innovation variance of the AR(`k`) fit, `k = 0..=K`.
16    pub sigma2: Vec<f64>,
17    /// Coefficients `(φ_1, …, φ_K)` of the AR(`K`) fit (the "long AR").
18    pub phi_last: Vec<f64>,
19}
20
21/// Run Levinson–Durbin on autocovariances `r[0..=K]`.
22pub fn levinson(r: &[f64]) -> Result<Levinson> {
23    if r.len() < 2 {
24        return Err(Error::TooShort {
25            n: r.len(),
26            needed: 2,
27        });
28    }
29    let k_max = r.len() - 1;
30    let mut phi: Vec<f64> = Vec::with_capacity(k_max);
31    let mut pacf = Vec::with_capacity(k_max);
32    let mut sigma2 = Vec::with_capacity(k_max + 1);
33    let mut sig = r[0];
34    if sig <= 0.0 {
35        return Err(Error::Singular);
36    }
37    sigma2.push(sig);
38    for k in 1..=k_max {
39        let mut num = r[k];
40        for j in 1..k {
41            num -= phi[j - 1] * r[k - j];
42        }
43        let kk = num / sig;
44        if !kk.is_finite() {
45            return Err(Error::Singular);
46        }
47        let prev = phi.clone();
48        for j in 1..k {
49            phi[j - 1] = prev[j - 1] - kk * prev[k - j - 1];
50        }
51        phi.push(kk);
52        sig *= 1.0 - kk * kk;
53        if sig <= 0.0 {
54            // Perfectly predictable at this order; clamp so later orders stay finite.
55            sig = f64::EPSILON * r[0];
56        }
57        pacf.push(kk);
58        sigma2.push(sig);
59    }
60    Ok(Levinson {
61        pacf,
62        sigma2,
63        phi_last: phi,
64    })
65}
66
67/// Convenience: autocovariance + Levinson in one call.
68pub fn levinson_from_series(x: &[f64], k_max: usize) -> Result<Levinson> {
69    check_input(x, k_max + 2)?;
70    let r = crate::acov::autocovariance(x, k_max)?;
71    levinson(&r)
72}
73
74/// Burg's method (Burg 1967; Kay & Marple 1981): each reflection coefficient is
75/// chosen to minimise the sum of forward and backward prediction error power,
76/// so the fitted AR is always stationary and behaves better than Yule–Walker
77/// near the unit circle. Returns the same shape as [`levinson`]; `sigma2[0]` is
78/// the biased sample variance and `sigma2[k] = sigma2[k-1] (1 - κ_k²)`.
79pub fn burg_from_series(x: &[f64], k_max: usize) -> Result<Levinson> {
80    check_input(x, k_max + 2)?;
81    let n = x.len();
82    let mean = x.iter().sum::<f64>() / n as f64;
83    let xc: Vec<f64> = x.iter().map(|v| v - mean).collect();
84    let mut sig = xc.iter().map(|v| v * v).sum::<f64>() / n as f64;
85    if sig <= 0.0 {
86        return Err(Error::Singular);
87    }
88    let r0 = sig;
89    let mut f = xc.clone();
90    let mut b = xc;
91    let mut phi: Vec<f64> = Vec::with_capacity(k_max);
92    let mut pacf = Vec::with_capacity(k_max);
93    let mut sigma2 = Vec::with_capacity(k_max + 1);
94    sigma2.push(sig);
95    for k in 1..=k_max {
96        let mut num = 0.0;
97        let mut den = 0.0;
98        for t in k..n {
99            num += f[t] * b[t - 1];
100            den += f[t] * f[t] + b[t - 1] * b[t - 1];
101        }
102        if den <= 0.0 || !den.is_finite() {
103            return Err(Error::Singular);
104        }
105        let kk = 2.0 * num / den;
106        // Descending so f[t] and b[t-1] on the right-hand side are still the
107        // order-(k-1) errors when they are read.
108        for t in (k..n).rev() {
109            let fo = f[t];
110            let bo = b[t - 1];
111            f[t] = fo - kk * bo;
112            b[t] = bo - kk * fo;
113        }
114        let prev = phi.clone();
115        for j in 1..k {
116            phi[j - 1] = prev[j - 1] - kk * prev[k - j - 1];
117        }
118        phi.push(kk);
119        sig *= 1.0 - kk * kk;
120        if sig <= 0.0 {
121            // Perfectly predictable at this order; clamp so later orders stay finite.
122            sig = f64::EPSILON * r0;
123        }
124        pacf.push(kk);
125        sigma2.push(sig);
126    }
127    Ok(Levinson {
128        pacf,
129        sigma2,
130        phi_last: phi,
131    })
132}
133
134#[cfg(test)]
135mod tests {
136    use super::*;
137
138    #[test]
139    fn ar1_exact_autocovariance_recovers_phi() {
140        // γ(k) = φ^k / (1-φ²) for AR(1) with unit innovations.
141        let phi: f64 = 0.6;
142        let r: Vec<f64> = (0..6).map(|k| phi.powi(k) / (1.0 - phi * phi)).collect();
143        let l = levinson(&r).unwrap();
144        assert!((l.pacf[0] - phi).abs() < 1e-12);
145        for k in 1..5 {
146            assert!(l.pacf[k].abs() < 1e-10, "pacf[{k}] = {}", l.pacf[k]);
147        }
148        assert!((l.sigma2[1] - 1.0).abs() < 1e-10);
149    }
150}
tsorder_core/levinson.rs

tsorder_core/
levinson.rs
