Forecasting: Method

Exponential Smoothing, Prediction Intervals, Accuracy Measures and Decomposition

What the tool does

A single series of equally spaced observations is pasted in, oldest first. The tool fits eight forecasting methods to it, reports the one chosen in detail (parameters, the equations of the last step, the forecast table with prediction intervals), scores every method in sample and, when a holdout is set, on the last observations it was not shown, checks the residuals of the chosen method and, for a seasonal series, decomposes it into trend, seasonal and remainder. Nothing beyond the series itself is used: no regressors, no calendar, no judgement.

The smoothing family

Exponential smoothing forecasts with weighted averages of the past in which the weights decay geometrically: the most recent observation carries weight α, the one before it α(1 − α), the next α(1 − α)², and so on. A large α tracks the series closely and forgets quickly; a small α smooths heavily and reacts slowly. The same idea applied to the level, the trend and the seasonal component gives the family (Hyndman and Athanasopoulos 2018, chapter 7), written here in the recursive form the tool computes:

• Simple exponential smoothing (SES): ℓ_t = α y_t + (1 − α) ℓ_t−1; forecast ŷ_T+h = ℓ_T, flat. For data with no trend and no season.
• Holt's linear trend: ℓ_t = α y_t + (1 − α)(ℓ_t−1 + b_t−1), b_t = β (ℓ_t − ℓ_t−1) + (1 − β) b_t−1; forecast ℓ_T + h b_T. The trend parameter is written β* in Hyndman and Athanasopoulos; the tool's β is that β*. The damped trend multiplies the carried trend by φ at every step (0.8 to 0.98), so the forecast approaches ℓ_T + φ b_T/(1 − φ) instead of growing without limit; it is the safer choice for long horizons.
• Holt-Winters additive: the seasonal component s_t is added, ℓ_t = α (y_t − s_t−m) + (1 − α)(ℓ_t−1 + b_t−1), s_t = γ (y_t − ℓ_t−1 − b_t−1) + (1 − γ) s_t−m; forecast ℓ_T + h b_T + s_T+h−m(k+1) with k the number of whole seasons in h − 1. For seasonal swings of constant size.
• Holt-Winters multiplicative: the seasonal component multiplies, ℓ_t = α y_t/s_t−m + (1 − α)(ℓ_t−1 + b_t−1), s_t = γ y_t/(ℓ_t−1 + b_t−1) + (1 − γ) s_t−m; forecast (ℓ_T + h b_T) s_T+h−m(k+1). For seasonal swings that grow with the level, which is the usual case for sales, traffic and demand. Needs a strictly positive series.

Three benchmarks are fitted alongside: the naive forecast (the last value), the seasonal naive (the value one season back) and the moving average of the last k values. A method that cannot beat them on the holdout has not earned its parameters. The linear trend regression y_t = a + b t (+ a dummy for each season but the last) is the global alternative to local smoothing: every observation counts equally, so it is right when the trend really is a straight line and wrong when the level has shifted.

The state-space view

Each smoothing method is the point forecast of an innovations state-space model, ETS(error, trend, season), in which the observation equals the forecast plus a white-noise error ε_t and the states move by the smoothing equations driven by that same error (Hyndman, Koehler, Ord and Snyder 2008). SES with additive errors is ETS(A,N,N), Holt's method ETS(A,A,N), the damped trend ETS(A,A_d,N), Holt-Winters additive ETS(A,A,A). The model is what makes the method a statistical statement: it gives the likelihood, the variance of an h-step forecast and so the prediction intervals. The tool uses the models only for the intervals; the parameters are fitted by least squares on the one-step errors, which for the additive-error models is the maximum-likelihood fit.

Initialisation and optimisation

The smoothing parameters α, β, γ (and φ for the damped trend) are bounded between 0 and 1 (φ between 0.8 and 0.98) and, when left blank, chosen to minimise the sum of squared one-step errors Σ(y_t − ŷ_t|t−1)² over the whole series: a coarse grid first, then the Nelder-Mead simplex from the best grid point, restarted once from its optimum. An entered parameter is held fixed and the rest are fitted. The initial states are treated as parameters as well, as the forecast package does for the examples of Hyndman and Athanasopoulos:

• SES: ℓ₀ starts at y₁ and is optimised with α.
• Holt: ℓ₀ = y₁ and b₀ = y₂ − y₁ as the start, both optimised.
• Holt-Winters: the initial seasonals are the first season's deviations from (additive) or ratios to (multiplicative) its average, normalised to sum zero or average one; a fitted smoothing parameter within 0.0001 of a bound is reported at the bound; the initial trend is the difference between the averages of the first two seasons divided by m; the initial level is the first-season average moved back half a season by that trend. Level and trend are then optimised with the smoothing parameters, and the seasonals too for m ≤ 12 (normalised at every evaluation); for longer seasons (weekly data, m = 52) the seasonals keep their conventional values, since fifty extra parameters would overfit.

The one-step fitted values start at t = 1 for the smoothing methods (the initial state forecasts the first observation), at t = 2 for the naive, t = m + 1 for the seasonal naive and t = k + 1 for the moving average.

Prediction intervals and their assumptions

The residual standard deviation σ is the root of the sum of squared one-step errors over n − p, p the number of fitted parameters and states. For the additive-error models the variance of the h-step forecast error is σ² c_h with the multipliers of Hyndman et al. (2008, Table 6.2; Hyndman and Athanasopoulos 2018, §7.5), written with β = αβ* and k the number of whole seasons in h − 1:

• SES: c_h = 1 + α²(h − 1).
• Holt: c_h = 1 + (h − 1)[α² + αβh + β²h(2h − 1)/6].
• Damped trend: c_h = 1 + α²(h − 1) + βφh/(1 − φ)² [2α(1 − φ) + βφ] − βφ(1 − φ^h)/[(1 − φ)²(1 − φ²)] [2α(1 − φ²) + βφ(1 + 2φ − φ^h)].
• Holt-Winters additive: c_h = 1 + (h − 1)[α² + αβh + β²h(2h − 1)/6] + γk[2α + γ + βm(k + 1)].
• Naive: c_h = h (a random walk); seasonal naive: c_h = k + 1; moving average of k: the SES multiplier with the equivalent α = 2/(k + 1), since the average of k values has the same mean age as SES at that α.

The interval is ŷ ± z σ √c_h with z the normal quantile of the chosen level. The multiplicative seasonal method has no closed form; the tool simulates 2,000 futures from the final states, adding a normal error of standard deviation σ to each simulated observation and updating the states by the method's own equations, and reads the limits off the percentiles at each step. The regression uses the usual interval ŷ ± t_n−p σ √(1 + x′(X′X)⁻¹x). All of these assume uncorrelated, homoscedastic, normally distributed errors and that the fitted parameters are the true ones; the intervals are therefore too narrow in practice, more so for long horizons and when the Ljung-Box test fails. They say nothing about a change of regime.

Accuracy measures

With errors e_t = y_t − ŷ_t: ME is the mean error (bias); MAE the mean absolute error; RMSE the root mean squared error, which the fit minimises and which penalises large misses; MAPE the mean of 100 |e_t/y_t|, undefined at y_t = 0 and inflated near it, and asymmetric (it punishes a forecast above the truth more than one below). MASE (Hyndman and Koehler 2006) divides the MAE by the in-sample mean absolute error of the naive forecast (seasonal naive when m > 1), Q = mean |y_t − y_t−m|: it is scale-free, defined for any series, and reads directly as "this method's error is x times the benchmark's", so values below 1 beat the benchmark and values can be compared across series of different scale. The scaling uses the training data even for the holdout score, as the definition requires.

Holdout evaluation against in-sample fit

In-sample accuracy is the fit to data the method has already seen; a method with more parameters will always fit better, and a one-step fit says little about an h-step forecast. The holdout keeps the last H observations out of the fit, forecasts them from the rest and scores the forecasts: that is the only measure of how the method would have done. It is one trial on one stretch of the series, so it is noisy; it still ranks the methods more honestly than the in-sample numbers. The tool ranks by the holdout RMSE when a holdout is set and by the in-sample RMSE otherwise, highlighting the winner, and draws the holdout forecasts of every method against what happened. The reported method is always refitted to the whole series for its forecast.

Classical decomposition

For a seasonal series the classical additive decomposition (Hyndman and Athanasopoulos 2018, §6.3) estimates the trend-cycle T_t by the centred moving average of one season (a 2 × m-MA when m is even, so the average is centred on an observation), the seasonal component as the mean of y_t − T_t for each season, centred to sum to zero, and the remainder R_t = y_t − T_t − S_t. The multiplicative indices, the mean of y_t/T_t per season scaled to average one, are shown beside the additive values because demand planners use them as the ratio-to-moving-average seasonal factors. The seasonal strength F_S = max(0, 1 − Var(R)/Var(S + R)) is 0 for no seasonality and approaches 1 when the remainder is small beside the seasonal swing; it tells whether a seasonal method is worth its three parameters. The decomposition loses half a season at each end and assumes the seasonal pattern is the same every year.

Residual diagnostics and the Ljung-Box test

A method that has extracted what the series can tell leaves one-step residuals without autocorrelation. The tool lists the residual autocorrelations r_k for the first ten lags against the ±1.96/√n bound of white noise and computes the Ljung-Box portmanteau statistic Q = n(n + 2) Σ_k=1..L r_k²/(n − k) (Ljung and Box 1978) over L = min(10, n/5) lags, or min(2m, n/5) for a seasonal series, with the p-value from χ² on L minus the number of smoothing parameters degrees of freedom. A small p (below 0.05) says the residuals are not white noise: the method leaves structure unused and its intervals are optimistic; the cure is usually a different method (a seasonal one when the lag-m autocorrelation is large, a trend when the low lags are) rather than different parameters.

Reading the default case

The default series is the quarterly international visitor nights in Australia from 2005 to 2015 used by Hyndman and Athanasopoulos (2018, §7.3). The seasonal swing grows with the level, so the multiplicative Holt-Winters method fits better than the additive one (RMSE 1.55 against 1.76 in sample; the book prints 1.576 and 1.763 from a slightly different optimiser), and the seasonal strength is high. The test script test-scripts/test_forecasting.php reproduces the book's SES (oil production: α = 0.834, ℓ₀ = 446.6, forecast 542.68), Holt (air passengers: forecasts 74.60 to 83.01) and Holt-Winters examples within the stated tolerances.

Limits

• No ARIMA, no regression on explanatory variables, no dynamic or hierarchical models; the eight methods above are the whole menu.
• No intermittent-demand methods (Croston and its variants): a series with many zeros should not be smoothed this way, and MAPE is undefined on it.
• No calendar effects (trading days, holidays, leap years) and a single seasonal period; daily data with weekly and yearly patterns needs a different tool.
• No transformation (log, Box-Cox); the multiplicative method covers the common case of a seasonal swing proportional to the level.
• Observations must be equally spaced with no gaps; the series is limited to 5,000 values and the season length to 366.
• The intervals assume normal, independent errors and known parameters; they are too narrow for long horizons and say nothing about structural change.