u-analytics 0.33.0

//! Estimation of smoothing constants by least squares.
//!
//! The constants left unset are those minimising the sum of squared one-step
//! errors over the observations the smoother updates on (Hyndman et al. 2008,
//! §5.1 — the least-squares estimate, as statsmodels' `ExponentialSmoothing.fit`
//! computes by default). Each free constant is searched on the logistic scale,
//! so it stays strictly inside (0, 1): a 5-point grid per constant gives the
//! start, then the Nelder–Mead simplex method refines it. The search is
//! deterministic.

use u_numflow::optimize::{nelder_mead, NelderMeadOptions};

/// A logistic-scale value far enough out that the constant is 2e-9 from the
/// end of (0, 1) without reaching it.
const EDGE: f64 = 20.0;

fn constant(u: f64) -> f64 {
    1.0 / (1.0 + (-u.clamp(-EDGE, EDGE)).exp())
}

fn logit(p: f64) -> f64 {
    (p / (1.0 - p)).ln()
}

/// The constants `fixed` with every `None` replaced by the value that
/// minimises `sse`; `sse` is called with all of them.
pub(crate) fn least_squares(fixed: &[Option<f64>], sse: impl Fn(&[f64]) -> f64) -> Vec<f64> {
    let free: Vec<usize> = (0..fixed.len()).filter(|&i| fixed[i].is_none()).collect();
    let full = |u: &[f64]| -> Vec<f64> {
        let mut c: Vec<f64> = fixed.iter().map(|v| v.unwrap_or(0.0)).collect();
        for (k, &i) in free.iter().enumerate() {
            c[i] = constant(u[k]);
        }
        c
    };
    if free.is_empty() {
        return full(&[]);
    }
    let grid: Vec<f64> = [0.1, 0.3, 0.5, 0.7, 0.9]
        .iter()
        .map(|&p| logit(p))
        .collect();
    let d = free.len();
    let mut start = vec![grid[0]; d];
    let mut best = f64::INFINITY;
    let mut point = vec![0usize; d];
    loop {
        let u: Vec<f64> = point.iter().map(|&g| grid[g]).collect();
        let v = sse(&full(&u));
        if v < best {
            best = v;
            start = u;
        }
        // Next grid point, odometer order.
        let mut k = 0;
        while k < d {
            point[k] += 1;
            if point[k] < grid.len() {
                break;
            }
            point[k] = 0;
            k += 1;
        }
        if k == d {
            break;
        }
    }
    let m = nelder_mead(
        |u| sse(&full(u)),
        &start,
        NelderMeadOptions::for_dimension(d),
    );
    if m.value <= best {
        full(&m.x)
    } else {
        full(&start)
    }
}

/// The sum of squared one-step errors `data[t] − fitted[t]` for `t >= first`.
pub(crate) fn one_step_sse(data: &[f64], fitted: &[f64], first: usize) -> f64 {
    data.iter()
        .zip(fitted)
        .skip(first)
        .map(|(x, f)| (x - f).powi(2))
        .sum()
}

/// The least-squares constants against statsmodels 0.15.0 with the same
/// initial states, on seeded series rounded to 3 places: a noisy level
/// (simple) and a noisy line (Holt) against `ExponentialSmoothing(...,
/// initialization_method="known").fit()`; a period-6 seasonal line (additive
/// Holt-Winters) against `ETSModel(..., error="add").fit()` -- whose Gaussian
/// likelihood peaks where the sum of squares is least -- with its state-space
/// constants converted (beta = beta*/alpha, gamma = gamma*/(1 - alpha)), since
/// `ExponentialSmoothing` updates the season on y - l(t-1) - b(t-1), not on
/// Winters' y - l(t) as here.
#[cfg(test)]
#[allow(clippy::approx_constant)] // seeded data, not mathematical constants
mod statsmodels_reference {
    use super::*;
    use crate::smoothing::{
        HoltLinear, HoltResult, HoltWinters, Seasonality, SesResult, SimpleExponentialSmoothing,
    };

    // Generated by fit_ref.py - statsmodels 0.15.0, seed 20261029.
    const SES_DATA: &[f64] = &[
        29.475, 31.322, 31.938, 29.368, 31.634, 28.627, 28.92, 27.876, 29.177, 28.265, 29.61,
        29.944, 31.558, 30.053, 28.515, 30.12, 30.059, 29.895, 31.831, 29.918, 32.733, 31.246,
        30.705, 30.139, 28.105, 29.329, 32.088, 30.443, 30.486, 31.99, 32.107, 30.395, 29.677,
        32.174, 32.153, 32.518, 31.03, 33.966, 30.395, 31.598,
    ];
    const SES_PARAMS: &[f64] = &[0.18426008603813843];
    const SES_SSE: f64 = 72.19522040043387;
    const HOLT_DATA: &[f64] = &[
        10.49, 8.627, 11.0, 11.983, 10.598, 10.914, 12.878, 13.463, 14.365, 13.203, 14.253, 14.004,
        17.141, 17.345, 16.5, 18.067, 19.242, 18.887, 20.942, 20.084, 19.726, 21.863, 22.02,
        22.071, 24.73, 25.327, 26.286, 23.974, 26.131, 27.991, 28.502, 28.554, 30.023, 31.005,
        33.029, 30.639, 31.791, 33.536, 33.809, 34.901,
    ];
    const HOLT_PARAMS: &[f64] = &[0.6250044815126141, 0.37229625369293723];
    const HOLT_SSE: f64 = 90.06578406159154;
    const HW_DATA: &[f64] = &[
        54.041, 52.526, 49.687, 46.041, 48.653, 53.357, 59.842, 54.564, 52.027, 51.279, 50.912,
        56.422, 58.583, 57.37, 54.356, 51.113, 54.242, 57.343, 62.011, 58.622, 57.0, 53.567,
        54.037, 59.86, 64.447, 60.463, 59.352, 56.267, 59.016, 63.27, 68.032, 65.084, 61.567,
        58.568, 61.01, 64.063, 68.956, 66.112, 63.484, 61.251, 63.204, 67.139, 71.318, 70.241,
        66.962, 65.065, 65.899, 70.269, 74.911, 70.691, 69.215, 66.879, 67.432, 72.993, 76.296,
        74.423, 71.018, 69.041, 69.612, 72.759,
    ];
    const HW_PARAMS: &[f64] = &[
        0.31474194439015335,
        0.05809934357601471,
        0.41314819780867096,
    ];
    const HW_SSE: f64 = 82.94643429427231;

    /// The least-squares constants reach the reference's sum of squares (or
    /// below it), and at the reference's own constants the smoother gives the
    /// same sum -- the same model, the same objective.
    fn check(sse_at: impl Fn(&[f64]) -> f64, ours: &[f64], params: &[f64], sse: f64) {
        let at_reference = sse_at(params);
        assert!(
            (at_reference - sse).abs() <= 1e-9 * sse,
            "objective differs from the reference's: {at_reference} vs {sse}"
        );
        let at_ours = sse_at(ours);
        assert!(
            at_ours <= sse * (1.0 + 1e-6),
            "{ours:?}: {at_ours} vs {params:?}: {sse}"
        );
        for (a, b) in ours.iter().zip(params) {
            assert!((a - b).abs() < 1e-3, "{ours:?} vs {params:?}");
        }
    }

    #[test]
    fn ses_least_squares_matches_statsmodels() {
        let s = SimpleExponentialSmoothing::fit(SES_DATA).unwrap();
        let sse = |c: &[f64]| {
            let r = SimpleExponentialSmoothing::new(c[0])
                .unwrap()
                .smooth(SES_DATA)
                .unwrap();
            one_step_sse(SES_DATA, &r.fitted, SesResult::FIRST_ERROR)
        };
        check(sse, &[s.alpha()], SES_PARAMS, SES_SSE);
    }

    #[test]
    fn holt_least_squares_matches_statsmodels() {
        let h = HoltLinear::fit(HOLT_DATA, None, None).unwrap();
        let sse = |c: &[f64]| {
            let r = HoltLinear::new(c[0], c[1])
                .unwrap()
                .smooth(HOLT_DATA)
                .unwrap();
            one_step_sse(HOLT_DATA, &r.fitted, HoltResult::FIRST_ERROR)
        };
        check(sse, &[h.alpha(), h.beta()], HOLT_PARAMS, HOLT_SSE);
    }

    #[test]
    fn holt_winters_least_squares_matches_statsmodels() {
        let w = HoltWinters::fit(HW_DATA, None, None, None, 6, Seasonality::Additive).unwrap();
        let sse = |c: &[f64]| {
            let r = HoltWinters::new(c[0], c[1], c[2], 6, Seasonality::Additive)
                .unwrap()
                .smooth(HW_DATA)
                .unwrap();
            one_step_sse(HW_DATA, &r.fitted, 6)
        };
        let r = w.smooth(HW_DATA).unwrap();
        check(sse, &[r.alpha, r.beta, r.gamma], HW_PARAMS, HW_SSE);
    }

    #[test]
    fn a_given_constant_is_held() {
        let h = HoltLinear::fit(HOLT_DATA, Some(0.3), None).unwrap();
        assert_eq!(h.alpha(), 0.3);
        let free = HoltLinear::fit(HOLT_DATA, None, None).unwrap();
        assert_ne!(h.beta(), free.beta());
    }
}