Exponential smoothing (ETS, Holt-Winters)

Deep dive

How it works

In: new observation y_t, smoothing constants alpha, beta, gamma, phi · Out: one-step forecast and the updated states · State: level, trend and m seasonal values (m+2 floats) · Cost: O(1) per update

State-space ("innovations") form of additive damped-trend, additive-season ETS(A,Ad,A), the same recursion as additive Holt-Winters with reparametrized smoothing constants. Notation: l level, b trend, s_i seasonal, m season length, phi damping.

1. Predict: yhat_t = l + phi*b + s[t-m]; error e_t = y_t - yhat_t.
2. Update all states from the same error: l += phi*b + alpha*e, b = phi*b + beta*e, s[t] = s[t-m] + gamma*e.
3. Forecast h steps: l + (phi + phi^2 + ... + phi^h)*b + s[t+h-m].
4. Initialize states from the first two seasons (season means give level and trend, detrended first season gives seasonal), or estimate them jointly with the parameters.
5. Fit (alpha, beta, gamma, phi) by minimizing SSE or the likelihood; typical bounds are 0<alpha<1, 0<beta<alpha, 0<gamma<1-alpha, 0.8<=phi<=0.98.
6. Multiplicative variants replace + by * for error and/or season (needs strictly positive data); choose among models by AICc.

y: ...   y_t   -> e_t = y_t - (l + phi b + s_{t-m}) ->  l, b, s_i updated (O(1))
state = [l, b, s_0..s_{m-1}]   (m+2 floats, the entire history)

Code

/// Additive ETS(A,Ad,A): damped-trend Holt-Winters in innovations form.
/// e = y - (l + phi*b + s[t-m]); l += phi*b + a*e; b = phi*b + be*e; s[t] = s[t-m] + g*e
#[derive(Clone, Copy, Debug)]
pub struct Params { pub alpha: f64, pub beta: f64, pub gamma: f64, pub phi: f64 }

pub struct State { pub level: f64, pub trend: f64, pub season: Vec<f64>, pub pos: usize } // season[pos] = oldest

pub fn init(y: &[f64], m: usize) -> State {
    let mean = |s: &[f64]| s.iter().sum::<f64>() / s.len() as f64;
    let (m1, m2) = (mean(&y[..m]), mean(&y[m..2 * m]));
    let b = (m2 - m1) / m as f64;                 // trend: change of season means per step
    let c = (m as f64 - 1.0) / 2.0;                // centre of the first season
    State { level: m1 - b * (c + 1.0), trend: b,   // level just before y[0]
            season: y[..m].iter().enumerate().map(|(i, v)| v - (m1 + b * (i as f64 - c))).collect(), pos: 0 }
}

/// One step: returns the one-step-ahead forecast error and advances the state.
pub fn step(st: &mut State, p: &Params, y: f64) -> f64 {
    let s = st.season[st.pos];
    let pred = st.level + p.phi * st.trend + s;
    let e = y - pred;
    st.level += p.phi * st.trend + p.alpha * e;
    st.trend = p.phi * st.trend + p.beta * e;
    st.season[st.pos] = s + p.gamma * e;
    st.pos = (st.pos + 1) % st.season.len();
    e
}

/// Sum of squared one-step errors over the series (what the optimizer minimizes).
pub fn sse(y: &[f64], m: usize, p: &Params) -> f64 {
    let mut st = init(y, m);
    y.iter().map(|&v| { let e = step(&mut st, p, v); e * e }).sum()
}

/// h-step forecast: l + (phi + ... + phi^h) b + s[t+h-m]
pub fn forecast(st: &State, p: &Params, h: usize) -> Vec<f64> {
    let m = st.season.len();
    let (mut damp, mut acc) = (1.0, 0.0);
    (1..=h).map(|k| { damp *= p.phi; acc += damp;
        st.level + acc * st.trend + st.season[(st.pos + k - 1) % m] }).collect()
}

/// Derivative-free fit: coordinate search with shrinking step, inside the usual bounds
/// (0<alpha<1, 0<beta<alpha, 0<gamma<1-alpha, 0.8<=phi<=0.98).
pub fn fit(y: &[f64], m: usize) -> Params {
    let ok = |p: &Params| p.alpha > 0.0 && p.alpha < 1.0 && p.beta >= 0.0 && p.beta <= p.alpha
        && p.gamma >= 0.0 && p.gamma <= 1.0 - p.alpha && p.phi >= 0.8 && p.phi <= 0.98;
    let starts = [(0.1, 0.01, 0.1), (0.3, 0.03, 0.3), (0.6, 0.05, 0.3)]; // multi-start: SSE is not convex
    let mut overall: Option<(f64, Params)> = None;
    for (alpha, beta, gamma) in starts {
        let mut best = Params { alpha, beta, gamma, phi: 0.95 };
        let mut best_sse = sse(y, m, &best);
        let mut step_size = 0.2;
        while step_size > 1e-4 {
            let mut improved = false;
            for i in 0..8 {
                let mut c = best;
                let d = if i % 2 == 0 { step_size } else { -step_size };
                match i / 2 { 0 => c.alpha += d, 1 => c.beta += d * 0.2, 2 => c.gamma += d, _ => c.phi += d * 0.1 }
                if !ok(&c) { continue; }
                let s = sse(y, m, &c);
                if s < best_sse { best = c; best_sse = s; improved = true; }
            }
            if !improved { step_size *= 0.5; }
        }
        if overall.map_or(true, |(s, _)| best_sse < s) { overall = Some((best_sse, best)); }
    }
    overall.unwrap().1
}

fn main() {
    let m = 12;
    // Synthetic monthly series: trend + seasonal + deterministic pseudo-noise (LCG).
    let mut x: u64 = 12345;
    let mut noise = || { x = x.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
                         ((x >> 33) as f64 / (1u64 << 31) as f64 - 0.5) * 2.0 };
    let f = |t: usize| 100.0 + 0.8 * t as f64 + 10.0 * (2.0 * std::f64::consts::PI * t as f64 / m as f64).sin();
    let y: Vec<f64> = (0..96).map(|t| f(t) + noise()).collect();
    let (train, test) = (&y[..84], &y[84..]);

    let p = fit(train, m);
    let mut st = init(train, m);
    for &v in train { step(&mut st, &p, v); }
    let fc = forecast(&st, &p, 12);
    let mae = fc.iter().zip(test).map(|(a, b)| (a - b).abs()).sum::<f64>() / 12.0;
    let default = Params { alpha: 0.3, beta: 0.03, gamma: 0.1, phi: 0.95 };
    println!("{:?}\nsse fit {:.1} vs default {:.1}\nMAE 12 ahead {:.2}", p, sse(train, m, &p), sse(train, m, &default), mae);
    assert!(sse(train, m, &p) <= sse(train, m, &default));
    assert!(mae < 3.0);
    // damped forecast of a pure level series is flat
    let flat = vec![5.0; 30];
    let (pp, mut s2) = (Params { alpha: 0.5, beta: 0.1, gamma: 0.1, phi: 0.9 }, init(&flat, 4));
    for &v in &flat { step(&mut s2, &pp, v); }
    assert!(forecast(&s2, &pp, 8).iter().all(|v| (v - 5.0).abs() < 1e-9));
}

The optimizer is a toy coordinate search with three starts; production code uses Nelder-Mead or L-BFGS on the likelihood. SSE is not convex in (alpha, beta, gamma), hence the multi-start.

Complexity and constants

• Update: about 10 flops and one seasonal slot read/write per observation; forecasting h steps is O(h).
• State: (m+2) f64, e.g. 14 doubles = 112 bytes for monthly data; 8760+2 doubles (about 70 KB) for hourly with a yearly period, which is where the model stops being sensible.
• Measured by me (Rust, Apple silicon, m=12, scalar, release build): about 4 ns per step for the SSE pass; a full fit on 120 points took about 0.3 ms. Order of magnitude only, one machine.
• Fit cost = (number of optimizer evaluations, order of hundreds to thousands) x n x per-step cost.

Efficiency

• The recursion is a serial dependency chain over time (level depends on the last level), so one series runs latency-bound, not throughput-bound. Parallelism comes from running thousands of independent series side by side, one per core or per SIMD lane (structure-of-arrays: all levels in one array, all trends in another).
• Seasonal states live in a ring buffer of m doubles: a single cache line pair for monthly data, no allocation in the loop.
• Online serving: persist only the state vector and parameters, update in O(1) per new point instead of refitting; refit on a schedule (daily or weekly), not per point.
• Cheap training: fit on the last few seasons only (old data has exponentially small weight anyway), warm-start the optimizer from yesterday's parameters, and fit one model family per series group rather than searching all ETS variants for every series.
• Energy: no GPU and no training set to hold; a CPU forecasts millions of series in a batch window. Statistical baselines are also the reference for judging whether a neural forecaster's much higher compute is worth it (see Deep forecasting models).
• Numbers above are measured on one machine; no benchmark claims are made for other libraries.

Implementation notes

• Use f64; f32 accumulates visible drift in the trend over long series.
• Seasonal init with fewer than two full seasons is unreliable: fall back to non-seasonal Holt or estimate the initial states in the optimizer.
• Constrain parameters (the "usual" region above, or the admissible region) so the recursion stays stable; unconstrained fits give alpha near 1 (naive forecast) on noisy data.
• Prediction intervals for additive models are closed form: for ETS(A,N,N) the h-step variance is sigma^2 * (1 + (h-1)*alpha^2); for multiplicative models simulate sample paths.
• Estimate on one-step errors but evaluate on the horizon you care about, with rolling-origin validation (time-series cross-validation).
• Reconcile forecasts across a hierarchy (store, region, total) afterwards: hierarchical reconciliation.
• The one-step errors e_t are a ready anomaly signal: large standardized errors flag outliers (anomaly detection); clip or skip them in the update so one outlier does not shift the level.

Pitfalls

• Multiplicative seasonality or error on data with zeros or negatives breaks (division by level near zero).
• Undamped trend extrapolated over long horizons overshoots; damp it or drop the trend.
• alpha fitted at the boundary (0 or 1) signals a wrong model (deterministic level, or random walk), not a precise one.
• Non-integer or very long seasonal periods are not representable; m must be an integer number of steps.
• Missing values: skip the update (keep states, advance the seasonal index) rather than inserting zeros.
• Comparing SSE or AIC across series transformed differently (Box-Cox vs raw) is meaningless.

Variants and alternatives

In production

Demand and capacity forecasting baselines (R forecast::ets, statsmodels, Nixtla statsforecast), sales and inventory planning (Winters' 1960 method was built for sales forecasting), metric monitoring (Prometheus double_exponential_smoothing), and a component of the winning hybrid method in the M4 competition (Smyl, ES-RNN).

Implementations

Sources

• Rob J. Hyndman, George Athanasopoulos, Forecasting: Principles and Practice, 3rd ed., OTexts 2021, ch. 8 (ETS). book
• Rob J. Hyndman, Anne B. Koehler, Ralph D. Snyder, Simone Grose, A state space framework for automatic forecasting using exponential smoothing methods, Int. J. Forecasting 18(3), 2002. paper
• Rob J. Hyndman, Anne B. Koehler, J. Keith Ord, Ralph D. Snyder, Forecasting with Exponential Smoothing: The State Space Approach, Springer 2008. book
• Charles C. Holt, Forecasting seasonals and trends by exponentially weighted moving averages, Int. J. Forecasting 20(1), 2004 (reprint of the 1957 report). paper
• Peter R. Winters, Forecasting sales by exponentially weighted moving averages, Management Science 6(3), 1960. paper
• Everette S. Gardner Jr., Ed McKenzie, Forecasting trends in time series, Management Science 31(10), 1985. paper
• Slawek Smyl, A hybrid method of exponential smoothing and recurrent neural networks for time series forecasting, Int. J. Forecasting 36(1), 2020. paper