u_analytics/smoothing/

forecast.rs

//! Forecasts of the smoothers with their prediction intervals.
//!
//! The additive smoothers are the exponential smoothing state-space models
//! with additive errors -- SES is ETS(A,N,N), Holt ETS(A,A,N), additive
//! Holt-Winters ETS(A,A,A) -- whose h-step forecast variance has a closed form
//! (Hyndman, Koehler, Ord & Snyder 2008, Table 6.1, class 1):
//!
//! ```text
//! Var(h) = σ² (1 + Σ_{j=1}^{h-1} c_j²)
//! SES:          c_j = α
//! Holt:         c_j = α (1 + j β)
//! Holt-Winters: c_j = α (1 + j β) + γ (1 − α) [j mod m = 0]
//! ```
//!
//! with the smoothing constants of the component form used here (the
//! state-space form's trend and seasonal constants are αβ and γ(1 − α)), and σ²
//! the mean squared one-step error over the observations the smoother updates
//! on (not those its initialisation reproduces). Multiplicative seasonality has
//! no such closed form (class 3), so its forecasts carry no interval.
//!
//! # Reference
//!
//! Hyndman, R.J., Koehler, A.B., Ord, J.K. & Snyder, R.D. (2008).
//! *Forecasting with Exponential Smoothing: The State Space Approach*.
//! Springer, chapter 6.

/// One forecast of a smoother: `step` periods after the last observation.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct SmoothingForecast {
    /// Steps ahead, from 1.
    pub step: usize,
    /// The point forecast.
    pub mean: f64,
    /// Lower end of the prediction interval; `None` when the smoother has no
    /// closed-form interval (multiplicative seasonality) or no one-step errors
    /// to estimate σ from.
    pub lower: Option<f64>,
    /// Upper end of the prediction interval.
    pub upper: Option<f64>,
}

/// Mean squared one-step error `data[t] − fitted[t]` for `t >= first`, or
/// `None` when there is none.
pub(crate) fn one_step_mse(data: &[f64], fitted: &[f64], first: usize) -> Option<f64> {
    let errors: Vec<f64> = data
        .iter()
        .zip(fitted)
        .skip(first)
        .map(|(x, f)| x - f)
        .collect();
    (!errors.is_empty()).then(|| errors.iter().map(|e| e * e).sum::<f64>() / errors.len() as f64)
}

/// The one-step errors from `first` on, `None` before it (where the fitted
/// value is the initialisation).
pub(crate) fn residuals(data: &[f64], fitted: &[f64], first: usize) -> Vec<Option<f64>> {
    data.iter()
        .zip(fitted)
        .enumerate()
        .map(|(t, (x, f))| (t >= first).then(|| x - f))
        .collect()
}

/// Forecasts `mean(h)` for h = 1..=horizon with the class-1 interval at
/// two-sided `level` when `mse` and `c` are given.
pub(crate) fn with_intervals(
    horizon: usize,
    level: f64,
    mse: Option<f64>,
    c: Option<&dyn Fn(usize) -> f64>,
    mean: impl Fn(usize) -> f64,
) -> Vec<SmoothingForecast> {
    let z = u_numflow::special::inverse_normal_cdf(0.5 + level / 2.0);
    let mut sum_c2 = 0.0;
    (1..=horizon)
        .map(|h| {
            let m = mean(h);
            let half = match (mse, c) {
                (Some(mse), Some(c)) => {
                    if h > 1 {
                        sum_c2 += c(h - 1).powi(2);
                    }
                    Some(z * (mse * (1.0 + sum_c2)).sqrt())
                }
                _ => None,
            };
            SmoothingForecast {
                step: h,
                mean: m,
                lower: half.map(|w| m - w),
                upper: half.map(|w| m + w),
            }
        })
        .collect()
}

/// `level` strictly inside (0, 1).
pub(crate) fn valid_level(level: f64) -> bool {
    level > 0.0 && level < 1.0
}

/// The fitted values, forecasts and prediction intervals of the additive
/// smoothers against statsmodels 0.15.0 `ETSModel(..., error="add",
/// initialization_method="known").smooth(params)` with the same initial states
/// (`get_prediction(...).summary_frame(alpha=0.1)`), on seeded series rounded to
/// 3 places: a random walk (SES), a noisy line (Holt) and a period-4 seasonal
/// line (Holt-Winters), six steps ahead.
#[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 ets_ref.py - statsmodels 0.15.0, seed 20261028.
    const SES_DATA: &[f64] = &[
        49.695, 48.755, 48.385, 48.678, 47.146, 48.813, 49.806, 49.983, 50.957, 49.731, 48.969,
        49.761, 48.5, 49.557, 50.796, 51.74, 52.288, 55.501, 55.645, 55.541, 55.655, 55.734,
        56.497, 55.544, 53.348, 52.352, 53.479, 52.5, 53.481, 55.59,
    ];
    const SES_FITTED: &[f64] = &[
        49.695,
        49.366,
        49.02265,
        48.9020225,
        48.287414625000004,
        48.47136950625,
        48.9384901790625,
        49.30406861639062,
        49.8825946006539,
        49.82953649042504,
        49.52834871877628,
        49.609776667204585,
        49.22135483368298,
        49.338830641893935,
        49.848839917231054,
        50.51074594620019,
        51.13278486503012,
        52.66166016226957,
        53.70582910547522,
        54.34813891855889,
        54.80554029706328,
        55.130501193091135,
        55.60877577550924,
        55.586104254081,
        54.80276776515265,
        53.94499904734923,
        53.781899380777,
        53.33323459750505,
        53.38495248837828,
    ];
    const SES_MEAN: &[f64] = &[
        54.15671911744588,
        54.15671911744588,
        54.15671911744588,
        54.15671911744588,
        54.15671911744588,
        54.15671911744588,
    ];
    const SES_LOWER: &[f64] = &[
        51.50271152380141,
        51.344848484798455,
        51.195388909183464,
        51.053118506866284,
        50.917089987805795,
        50.786547453385566,
    ];
    const SES_UPPER: &[f64] = &[
        56.81072671109035,
        56.96858975009331,
        57.1180493257083,
        57.26031972802548,
        57.39634824708597,
        57.5268907815062,
    ];
    const HOLT_DATA: &[f64] = &[
        18.624, 22.028, 22.022, 22.77, 23.168, 25.438, 23.856, 26.135, 25.839, 25.885, 27.652,
        29.208, 30.853, 32.14, 31.311, 31.097, 32.108, 31.913, 36.539, 34.939, 34.931, 36.448,
        37.208, 36.854, 38.229, 38.696, 41.158, 43.842, 42.405, 41.971,
    ];
    const HOLT_FITTED: &[f64] = &[
        25.432,
        27.199199999999998,
        28.204383999999997,
        28.563783679999997,
        29.4373607936,
        28.882198198272,
        29.24092478521344,
        29.065606754561227,
        28.724367395805,
        29.140634388886866,
        30.018183833625034,
        31.269498793777935,
        32.60472786636744,
        32.97076708061175,
        32.954889242709385,
        33.28201140055122,
        33.29076378321223,
        35.406274110151855,
        35.99819837750347,
        36.26477706771417,
        37.04618211642346,
        37.83197057633515,
        38.08360600617535,
        38.796218783585445,
        39.402568947344676,
        40.89161352981263,
        43.094671196908386,
        43.78653210141317,
    ];
    const HOLT_MEAN: &[f64] = &[
        43.88280607600299,
        44.70529289115808,
        45.52777970631317,
        46.35026652146826,
        47.172753336623344,
        47.99524015177843,
    ];
    const HOLT_LOWER: &[f64] = &[
        39.881805953904774,
        40.26724764211682,
        40.55622511518835,
        40.758018971303315,
        40.88214280288195,
        40.93713118743697,
    ];
    const HOLT_UPPER: &[f64] = &[
        47.88380619810121,
        49.14333814019934,
        50.49933429743799,
        51.9425140716332,
        53.46336387036474,
        55.0533491161199,
    ];
    const HW_DATA: &[f64] = &[
        107.195, 99.797, 99.111, 104.953, 111.536, 103.86, 105.869, 111.842, 117.337, 112.456,
        110.451, 115.115, 124.234, 116.183, 117.719, 123.434, 127.431, 122.666, 123.802, 130.548,
        135.536, 127.843, 129.797, 133.84, 142.645, 133.886, 134.436, 141.34, 148.881, 140.176,
        141.191, 145.134,
    ];
    const HW_FITTED: &[f64] = &[
        108.57318749999999,
        103.53110312499999,
        104.42071096874999,
        112.21758513031247,
        116.37452167699682,
        110.342324187491,
        112.08772508093094,
        118.67183234829378,
        122.04503373510259,
        116.55551472827399,
        116.797425126349,
        123.78843943015733,
        129.19938105911905,
        122.02718587096355,
        122.78208766197088,
        129.56537788125988,
        135.15618721562316,
        129.1403441795491,
        129.4136407564024,
        136.01289275800133,
        140.4703440485413,
        134.67208484162995,
        135.38085344677242,
        141.08242228813216,
        147.0508746050361,
        140.6418295369387,
        141.43920083128222,
        147.6014010922094,
    ];
    const HW_MEAN: &[f64] = &[
        152.98701643046797,
        145.53169868698552,
        146.42543580819583,
        152.19971424007915,
        158.7573450971372,
        151.30202735365475,
    ];
    const HW_LOWER: &[f64] = &[
        150.4614216897758,
        142.8721382966387,
        143.61475408277198,
        149.22143937008624,
        155.42143213962436,
        147.77782707686427,
    ];
    const HW_UPPER: &[f64] = &[
        155.51261117116016,
        148.19125907733235,
        149.23611753361968,
        155.17798911007205,
        162.09325805465005,
        154.82622763044523,
    ];

    fn close(got: f64, want: f64, what: &str) {
        let tol = 1e-9 * want.abs().max(1.0);
        assert!((got - want).abs() <= tol, "{what}: {got} vs {want}");
    }

    fn compare(
        fitted: &[f64],
        first: usize,
        forecasts: &[SmoothingForecast],
        want_fitted: &[f64],
        want: [&[f64]; 3],
    ) {
        assert_eq!(fitted.len() - first, want_fitted.len());
        for (got, want) in fitted[first..].iter().zip(want_fitted) {
            close(*got, *want, "fitted");
        }
        for (i, f) in forecasts.iter().enumerate() {
            close(f.mean, want[0][i], "mean");
            close(f.lower.expect("interval"), want[1][i], "lower");
            close(f.upper.expect("interval"), want[2][i], "upper");
        }
    }

    #[test]
    fn ses_matches_ets_ann() {
        let r = SimpleExponentialSmoothing::new(0.35)
            .unwrap()
            .smooth(SES_DATA)
            .unwrap();
        let f = r.forecasts(SES_DATA, 6, 0.9).unwrap();
        compare(
            &r.fitted,
            SesResult::FIRST_ERROR,
            &f,
            SES_FITTED,
            [SES_MEAN, SES_LOWER, SES_UPPER],
        );
    }

    #[test]
    fn holt_matches_ets_aan() {
        let r = HoltLinear::new(0.4, 0.2)
            .unwrap()
            .smooth(HOLT_DATA)
            .unwrap();
        let f = r.forecasts(HOLT_DATA, 6, 0.9).unwrap();
        compare(
            &r.fitted,
            HoltResult::FIRST_ERROR,
            &f,
            HOLT_FITTED,
            [HOLT_MEAN, HOLT_LOWER, HOLT_UPPER],
        );
    }

    #[test]
    fn additive_holt_winters_matches_ets_aaa() {
        let r = HoltWinters::new(0.3, 0.1, 0.25, 4, Seasonality::Additive)
            .unwrap()
            .smooth(HW_DATA)
            .unwrap();
        let f = r.forecasts(HW_DATA, 6, 0.9).unwrap();
        compare(
            &r.fitted,
            r.first_error(),
            &f,
            HW_FITTED,
            [HW_MEAN, HW_LOWER, HW_UPPER],
        );
    }
}