Time-series forecasting (baselines, model families, backtests)

Deep dive

How it works

Parts:

• Frame — horizon h, frequency, season length m, lead time, point vs quantile output, hierarchy.
• Baselines — mean, naive, seasonal naive, drift (naive plus average historical change per step).
• Model families — ETS, ARIMA, lag and calendar regression, global models, intermittent, hierarchical.
• Backtest — rolling-origin evaluation, same folds for every candidate.
• Metrics — MASE (points), pinball loss (quantiles).

Operations:

1. Frame from the decision: lead time and aggregation level where action is taken.
2. Run baselines first; keep them in every comparison (FPP3 section 5.2).
3. Fit families. ETS: level, trend, season states (exponential smoothing). ARIMA: difference until stationary (FPP3 chooses the number of differences with KPSS tests), then AR terms (past values) and MA terms (past errors); SARIMA adds seasonal terms; search by AICc (Hyndman and Khandakar 2008). Regression: lags, rolling means, calendar, price or promotion known at forecast time, fit by OLS or gradient-boosted trees. Global: stack all series with a series id and train once.
4. Multi-step regression: recursive (predict t+1, feed back) or direct (one model per horizon step); Ben Taieb et al. 2012. The M5 winner averaged both kinds.
5. Backtest all candidates on identical rolling origins; MASE per horizon step; pinball loss for quantiles (see Code).
6. Hierarchy: bottom-up sums bottom forecasts; MinT chooses the combination minimizing total variance of the coherent forecasts (Wickramasuriya et al. 2019).
7. Intermittent series: Croston forecasts non-zero size and inter-demand interval separately (size divided by interval); the Syntetos-Boylan approximation corrects its bias.
8. Deploy with a retrain schedule; track error by horizon; fall back to a baseline on drift.

history --> [baselines] ----------------.
   |                                     v
   +--> [ETS | ARIMA | regression | global] --> rolling-origin backtest --> MASE / pinball
                                              |                              |
                                              +--> pick or average ----> reconcile --> forecast + intervals

Code

Rolling-origin backtest on a synthetic daily series (trend, weekly pattern, seeded noise): naive, seasonal naive and OLS trend plus weekday means, scored by MASE, plus a 90 % quantile with pinball loss.

use std::time::Instant;

const M: usize = 7; // weekly season
type Fc = fn(&[f64], usize) -> Vec<f64>;

fn naive(tr: &[f64], h: usize) -> Vec<f64> { vec![tr[tr.len() - 1]; h] }
fn snaive(tr: &[f64], h: usize) -> Vec<f64> {
    (0..h).map(|k| tr[tr.len() - M + k % M]).collect() // last observed season, repeated
}

const P: usize = 1 + M; // intercept, slope, 6 weekday dummies
fn row(t: usize) -> [f64; P] {
    let mut x = [0.0; P]; x[0] = 1.0; x[1] = t as f64;
    if t % M > 0 { x[1 + t % M] = 1.0 }
    x
}

/// OLS on [1, t, d_1..d_6] via normal equations and Gauss-Jordan with partial pivoting.
fn fit(tr: &[f64]) -> [f64; P] {
    let mut a = [[0.0; P + 1]; P]; // augmented X'X | X'y
    for (t, &y) in tr.iter().enumerate() {
        let x = row(t);
        for i in 0..P { for j in 0..P { a[i][j] += x[i] * x[j]; } a[i][P] += x[i] * y; }
    }
    for c in 0..P {
        let p = (c..P).max_by(|&i, &j| a[i][c].abs().total_cmp(&a[j][c].abs())).unwrap();
        a.swap(c, p);
        for r in 0..P {
            if r == c { continue; }
            let f = a[r][c] / a[c][c];
            for k in c..=P { a[r][k] -= f * a[c][k]; }
        }
    }
    std::array::from_fn(|i| a[i][P] / a[i][i])
}
fn predict(b: &[f64; P], t: usize) -> f64 { row(t).iter().zip(b).map(|(x, w)| x * w).sum() }

/// Linear trend + weekday means (direct multi-step: one fit, any horizon).
fn trend_season(tr: &[f64], h: usize) -> Vec<f64> {
    let b = fit(tr);
    (0..h).map(|k| predict(&b, tr.len() + k)).collect()
}

/// Pinball (quantile) loss for quantile level q; q = 0.5 gives half the absolute error.
fn pinball(y: f64, f: f64, q: f64) -> f64 { if y >= f { q * (y - f) } else { (1.0 - q) * (f - y) } }

/// Rolling-origin backtest (expanding window); returns (MASE, fits). Scale = in-sample
/// seasonal-naive MAE of each training slice (Hyndman & Athanasopoulos, FPP3 5.8).
fn backtest(y: &[f64], min_train: usize, h: usize, step: usize, f: Fc) -> (f64, usize) {
    let (mut sum, mut cnt, mut folds) = (0.0, 0, 0);
    let mut o = min_train;
    while o + h <= y.len() {
        let scale = (M..o).map(|t| (y[t] - y[t - M]).abs()).sum::<f64>() / (o - M) as f64;
        let fc = f(&y[..o], h);
        for k in 0..h { sum += (y[o + k] - fc[k]).abs() / scale; cnt += 1; }
        folds += 1; o += step;
    }
    (sum / cnt as f64, folds)
}

fn main() {
    let mut s: u64 = 7; // xorshift64 seeded RNG, noise ~ U(-1.5, 1.5)
    let mut rnd = || { s ^= s << 13; s ^= s >> 7; s ^= s << 17; (s >> 11) as f64 / (1u64 << 53) as f64 - 0.5 };
    let wk = [0.0, 1.0, 2.0, 1.5, 3.0, 8.0, 6.0]; // weekend peak
    let y: Vec<f64> = (0..420).map(|t| 100.0 + 0.08 * t as f64 + wk[t % M] + 3.0 * rnd()).collect();

    let mut res = Vec::new();
    for (name, f) in [("naive", naive as Fc), ("seasonal-naive", snaive), ("trend+weekday OLS", trend_season)] {
        let t0 = Instant::now();
        let reps = 200;
        let mut out = (0.0, 0);
        for _ in 0..reps { out = backtest(&y, 140, 14, 7, f); }
        let dt = t0.elapsed().as_secs_f64();
        println!("{name:18} MASE {:.3}  folds {}  {:.0} fits/s", out.0, out.1, (reps * out.1) as f64 / dt);
        res.push(out.0);
    }
    assert!(res[1] < res[0], "seasonal-naive must beat naive");
    assert!(res[2] < res[1], "fitted model must beat seasonal-naive");

    // q90 = point forecast + 0.9-quantile of in-sample residuals
    let tr = &y[..350];
    let b = fit(tr);
    let mut r: Vec<f64> = tr.iter().enumerate().map(|(t, v)| v - predict(&b, t)).collect();
    r.sort_by(f64::total_cmp);
    let q90 = r[(0.9 * r.len() as f64) as usize];
    let fc = trend_season(tr, 70);
    let loss: f64 = (0..70).map(|k| pinball(y[350 + k], fc[k] + q90, 0.9)).sum::<f64>() / 70.0;
    let cover = (0..70).filter(|&k| y[350 + k] <= fc[k] + q90).count() as f64 / 70.0;
    println!("q90 offset {q90:.2}  pinball {loss:.3}  coverage {cover:.2}");
}

The OLS model is direct (one fit serves any horizon). MASE divides by the in-sample seasonal-naive MAE of the fold's training slice.

Complexity and constants

• Baselines O(1) per step; ETS and ARIMA O(n) per likelihood evaluation; OLS with p features O(n p^2) to build, O(p^3) to solve.
• Backtest = folds x (fit + predict + scale); 39 folds here (n=420, first origin 140, h=14, step 7).
• Measured on this Mac (Apple M4 Pro, release, one thread, p=8): naive about 9 M folds/s, seasonal naive about 6 M folds/s, OLS about 0.48 M fits/s (about 2 µs per fit). Example results: MASE naive 2.86, seasonal naive 1.11, OLS 0.69; q90 pinball 0.157, coverage 0.90 on 70 held-out days.

Efficiency

• Cost: a baseline copies a few floats; the OLS fit above is about 2 µs (measured, M4 Pro). At scale the bound is series x folds x fits; memory O(n) per series, or one n x features table.
• Levers:
  • fewer, wider folds (step = h) — cost falls with folds (cross-validation) — stable error over several seasons.
  • one global model — the M5 paper notes one model lowers overall compute and eases short histories — aligned related series.
  • state update instead of refit — O(1) per point for ETS-type models — between scheduled retrains.
  • drop candidates that lose to seasonal naive on a few folds — large grids.
• Naive → standard → state of the art:
  • seasonal naive · O(1) per step, no fit
  • per-series ETS/ARIMA selected by AICc · one fit per series per retrain
  • global LightGBM on lag and calendar features, recursive and direct averaged · one training run over all series (M5 winner family)

Implementation notes

• Features must be known at forecast time (forecast weather, planned promotions); in direct mode lags shorter than h are unavailable.
• One transform (log, Box-Cox) for all candidates; compare errors on the original scale.
• Quantiles: residual quantiles (as above), simulation, or pinball-loss training (quantile regression); check coverage per horizon, intervals tend to be too narrow far out (calibration).
• Holidays: calendar features, or decompose first (STL); outliers distort fits (robust statistics).

Pitfalls

• Leakage from features not known at forecast time (realized prices, actual weather): the backtest looks great, production fails.
• MAPE and sMAPE with zeros or tiny actuals are infinite or extreme (FPP3 section 5.8); use MASE.
• Forecasting through a structural break (price change, new store, policy): the old trend is wrong afterwards; shorten the window, add a level-shift feature or refit.
• Over-tuning on one backtest: the best of many configurations looks better than it is; keep an untouched final tail.
• Averaging errors across series in different units: use scaled errors.

Variants and alternatives

In production

Retail replenishment, grid load forecasting, cloud capacity planning, cash-flow forecasts, hospital and call-centre staffing, web traffic forecasting.

M4 (100,000 series, 61 methods): the winner was a hybrid of exponential smoothing and a recurrent network; the runner-up used XGBoost to weight statistical forecasts; both were similar in accuracy to a simple median of four statistical methods (per the M5 paper). M5 (42,840 retail series, 5,507 teams): 7.5 % of teams beat the best benchmark; the winner was 22.4 % better in WRMSSE; most top methods used LightGBM and all top 50 used cross-learning.

Implementations

Sources

• Rob J. Hyndman, George Athanasopoulos, Forecasting: Principles and Practice, 3rd ed., OTexts 2021: baselines (5.2), accuracy measures (5.8), distributional accuracy (5.9), ARIMA (9), hierarchies and MinT (11). book
• Rob J. Hyndman, Anne B. Koehler, Another look at measures of forecast accuracy, Int. J. Forecasting 22(4), 2006. paper
• Rob J. Hyndman, Yeasmin Khandakar, Automatic time series forecasting: the forecast package for R, J. Statistical Software 27(3), 2008. paper
• Spyros Makridakis, Evangelos Spiliotis, Vassilios Assimakopoulos, The M4 Competition: 100,000 time series and 61 forecasting methods, Int. J. Forecasting 36(1), 2020. paper
• Spyros Makridakis, Evangelos Spiliotis, Vassilios Assimakopoulos, M5 accuracy competition: Results, findings, and conclusions, Int. J. Forecasting 38(4), 2022. paper
• Shanika L. Wickramasuriya, George Athanasopoulos, Rob J. Hyndman, Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization, JASA 114(526), 2019. paper
• J. D. Croston, Forecasting and stock control for intermittent demands, Operational Research Quarterly 23(3), 1972. paper
• Souhaib Ben Taieb, Gianluca Bontempi, Amir F. Atiya, Antti Sorjamaa, A review and comparison of strategies for multi-step ahead time series forecasting based on the NN5 forecasting competition, Expert Systems with Applications 39(8), 2012. paper