Time-Series Forecasting

Forecasting is the job of saying what a quantity will be in the future and how sure you are. Retailers forecast demand per product and store to set inventory, grid operators forecast load every few minutes, cloud teams forecast capacity, and finance teams forecast cash. Forecasting is different from ordinary supervised learning in one way that matters a lot: the rows are ordered in time, so the past can be used to predict the future but never the other way round. Most forecasting failures in industry come from ignoring that fact.

Intuition: weather, habits and luck

Think of a coffee shop's daily sales. There is a slow drift up as the neighbourhood grows (trend). There is a weekly rhythm, busy on weekdays and quiet on Sunday (seasonality). And there is day-to-day noise that nobody could have predicted: a rainstorm, a bus strike (noise or the irregular component). A good forecaster separates these three. The trend and seasonality can be projected forward; the noise cannot, but its size tells you how wide your uncertainty band should be.

Components and decomposition

A time series is a sequence $y_1, y_2, \dots, y_T$ observed at regular intervals. A classical additive decomposition writes

where $T_t$ is the trend-cycle, $S_t$ the seasonal component with period $m$ (for example $m = 12$ for monthly data with a yearly pattern, $m = 7$ for daily data with a weekly pattern), and $R_t$ the remainder. When seasonal swings grow with the level of the series, a multiplicative decomposition $y_t = T_t \times S_t \times R_t$ fits better; taking logs turns it into an additive one. Robust tools such as STL (Seasonal and Trend decomposition using Loess) estimate these components and are the standard first look at a new series.

Stationarity and autocorrelation

A series is (weakly) stationary if its mean, variance and autocovariance do not change over time. Many classical models assume stationarity, and trends and seasonality break it. The usual fix is differencing: the first difference $y'_t = y_t - y_{t-1}$ removes a linear trend, and the seasonal difference $y_t - y_{t-m}$ removes a stable seasonal pattern. Unit-root tests (ADF, KPSS) help decide how much differencing is needed, but plots usually tell you first.

Autocorrelation measures how a series relates to its own past. The sample autocorrelation at lag $k$ is

where $\bar{y}$ is the sample mean. Plotting $r_k$ against $k$ gives the ACF. A slowly decaying ACF signals a trend; spikes at $k = m, 2m, \dots$ signal seasonality. The partial autocorrelation function (PACF) measures the correlation at lag $k$ after removing the effect of shorter lags, and together the ACF and PACF guide the choice of ARIMA orders. For a white-noise series, about 95% of $r_k$ values should fall inside $\pm 1.96/\sqrt{T}$.

Baselines first

Every forecasting project starts with baselines, because they are surprisingly hard to beat and they give your metrics a meaning.

Seasonal naive simply repeats the last observed season: next Monday equals last Monday. If a sophisticated model cannot beat it in a fair backtest, ship the baseline.

Exponential smoothing

Simple exponential smoothing (SES) forecasts with a weighted average of past observations whose weights decay geometrically. With smoothing parameter $0 < \alpha \le 1$ and level $\ell_t$:

A large $\alpha$ reacts quickly to new data; a small $\alpha$ smooths heavily. Holt's method adds a trend term $b_t$, and Holt–Winters adds a seasonal term $s_t$. The modern formulation wraps these in state-space ETS models (Error, Trend, Seasonal), which give likelihoods, automatic model selection by AIC, and proper prediction intervals.

Worked example: SES by hand

Take $y = (10, 12, 11, 14)$, $\alpha = 0.3$, and initialise $\ell_1 = y_1 = 10$.

• $\ell_2 = 0.3 \times 12 + 0.7 \times 10 = 3.6 + 7.0 = 10.6$
• $\ell_3 = 0.3 \times 11 + 0.7 \times 10.6 = 3.3 + 7.42 = 10.72$
• $\ell_4 = 0.3 \times 14 + 0.7 \times 10.72 = 4.2 + 7.504 = 11.704$

So the forecast for every future step is $\hat{y}_{5}, \hat{y}_{6}, \dots = 11.704$. SES produces flat forecasts; that is why Holt and Holt–Winters exist.

ARIMA in one page of notation

ARIMA combines three ideas. An autoregressive term of order $p$ regresses $y_t$ on its own past values; a moving-average term of order $q$ regresses on past forecast errors $\varepsilon_t$; and integration of order $d$ means the model is fitted to the $d$-times differenced series. Using the backshift operator $B y_t = y_{t-1}$, an ARIMA($p$, $d$, $q$) model is

where $\phi_i$ are AR coefficients, $\theta_j$ are MA coefficients, $c$ is a constant and $\varepsilon_t$ is white noise. Seasonal ARIMA, written ARIMA$(p,d,q)(P,D,Q)_m$, adds the same three pieces at seasonal lags $m, 2m, \dots$. Two famous special cases: ARIMA(0,1,0) is the random walk (the naive forecast), and ARIMA(0,1,1) is equivalent to simple exponential smoothing. In practice, automatic order search (for example auto.arima in R or AutoARIMA in the Python statsforecast package) picks orders by information criteria. ARIMA shines on single, well-behaved series with a few hundred points; it scales poorly to millions of series and ignores external drivers unless you add them as regressors (ARIMAX).

Machine learning: gradient boosting with lag features

When you have many related series (thousands of products) and external drivers (price, promotions, holidays, weather), global machine-learning models trained across all series usually win. You turn forecasting into tabular regression by building features from the past:

• Lags: $y_{t-1}, y_{t-7}, y_{t-364}$.
• Rolling statistics: mean and standard deviation over the last 7 or 28 steps, computed only from data before $t$.
• Calendar features: day of week, month, holiday flags.
• Known-future covariates: planned price or promotion.
• Series identifiers or static attributes (store, category).

Gradient-boosted trees (LightGBM, XGBoost, scikit-learn's HistGradientBoostingRegressor) dominated the top of the M5 retail forecasting competition this way. Two design decisions matter. First, recursive forecasting predicts one step, feeds the prediction back as a lag, and repeats, which compounds errors; direct forecasting trains a separate model (or a horizon feature) for each step ahead. Second, trees cannot extrapolate beyond the target values seen in training, so a trending series should be modelled as differences or ratios to a recent level.

import numpy as np
from sklearn.ensemble import HistGradientBoostingRegressor

rng = np.random.default_rng(0)
n, m = 156, 12
t = np.arange(n)
y = 100 + 0.5 * t + 10 * np.sin(2 * np.pi * t / m) + rng.normal(0, 3, n)

LAGS = [1, 2, 3, 12]
def make_xy(series, start):
    X, Y = [], []
    for i in range(start, len(series)):
        X.append([series[i - L] for L in LAGS] + [i % m])   # lags + month-of-year
        Y.append(series[i])
    return np.array(X), np.array(Y)

d = np.diff(y)                  # model first differences: trees cannot extrapolate a trend
split = n - 12                  # hold out the last 12 points of y
X, Y = make_xy(d[:split - 1], max(LAGS))
model = HistGradientBoostingRegressor(max_iter=300, learning_rate=0.05).fit(X, Y)

hist, level, preds = list(d[:split - 1]), y[split - 1], []
for k in range(12):             # recursive forecast: feed predictions back as lags
    i = len(hist)
    step = model.predict([[hist[i - L] for L in LAGS] + [i % m]])[0]
    hist.append(step); level += step; preds.append(level)

print("GBM MAE:", np.mean(np.abs(y[split:] - np.array(preds))).round(2))
print("Seasonal naive MAE:", np.mean(np.abs(y[split:] - y[split - 12:split])).round(2))

On this single synthetic series the boosted model (MAE about 3.6) does not beat seasonal naive (about 3.5). That is a realistic result, not a bug: a pure sine wave plus noise is exactly what seasonal naive is built for, and one series of 156 points is too little data for trees to shine. Gradient boosting earns its keep across many series with covariates.

Deep learning and pretrained time-series models

Neural forecasters learn across many series at once. Notable designs include DeepAR (Salinas et al., Amazon), an autoregressive recurrent network that outputs a probability distribution at each step; N-BEATS, a deep stack of fully connected blocks; the Temporal Fusion Transformer, which handles static, known-future and observed covariates with attention; and PatchTST, which splits a series into patches that act like tokens in a Transformer.

Since 2023 a new family has appeared: pretrained time-series foundation models, trained on very large collections of real and synthetic series and used zero-shot on new series, much as a language model answers questions it was not fine-tuned for. Examples with published papers include Amazon's Chronos (Ansari et al., 2024), which scales and quantises values into a vocabulary of tokens and trains a T5-style language model on them; Google's TimesFM (Das et al., ICML 2024), a decoder-only Transformer over patches; and Salesforce's Moirai (Woo et al., 2024). These models give strong baselines with no training, which is valuable for cold-start series and quick prototypes. They do not make backtesting optional: on a series with strong domain-specific drivers (promotions, outages), a tuned global gradient-boosting model with covariates is often still better. Treat a foundation model as one more candidate in your backtest. Model families and versions change fast; check the current model cards before you adopt one.

Backtesting with rolling-origin evaluation

A random train/test split leaks the future into training. The correct procedure is rolling-origin evaluation (also called time-series cross-validation or backtesting): choose several forecast origins $T_1 < T_2 < \dots$; for each origin, train only on data up to $T_i$, forecast $h$ steps ahead, and score against the actuals. The training window can expand (keep all history) or slide (fixed length, useful when old data is no longer relevant). scikit-learn's TimeSeriesSplit implements the expanding version. Use the same horizon your business needs: if purchasing decisions lock in two weeks ahead, the one-day-ahead error is irrelevant. Leave a gap between training data and the forecast window when features are not available in real time.

Metrics: MAE, MASE, sMAPE

Let $y_t$ be actuals and $\hat{y}_t$ forecasts over a test window of $H$ points.

• MAE $= \frac{1}{H}\sum |y_t - \hat{y}_t|$, in the units of the data. Optimising MAE targets the median.
• RMSE $= \sqrt{\frac{1}{H}\sum (y_t - \hat{y}_t)^2}$, which targets the mean and punishes large errors.
• sMAPE $= \frac{100\%}{H}\sum \frac{2|y_t - \hat{y}_t|}{|y_t| + |\hat{y}_t|}$, a percentage error that is bounded but unstable near zero and asymmetric despite its name.
• MASE (Hyndman & Koehler, 2006) scales the MAE by the in-sample MAE of the seasonal naive method on the training data of length $T$:

MASE below 1 means you beat the in-sample seasonal naive forecast. It is scale-free, so you can average it across thousands of series, and it is well defined when values are zero.

Worked example: computing the metrics

Non-seasonal training data ($m = 1$): $20, 22, 21, 25, 24$. Test actuals:$26, 27$. Forecasts:$25, 25$.

1. Absolute errors: $|26-25| = 1$, $|27-25| = 2$, so MAE $= 1.5$.
2. MASE scale: one-step naive in-sample errors are $|22-20| = 2$, $|21-22| = 1$, $|25-21| = 4$, $|24-25| = 1$; mean $= 8/4 = 2$. MASE $= 1.5/2 = 0.75$.
3. sMAPE: $\frac{2 \times 1}{26+25} = 0.0392$ and $\frac{2 \times 2}{27+25} = 0.0769$; mean $= 0.0581$, so sMAPE $\approx 5.81\%$.

Prediction intervals

A point forecast without uncertainty is hard to act on: a warehouse needs the 95th percentile of demand, not the mean. For the naive method with one-step residual standard deviation $\sigma$, the $h$-step forecast error variance grows as $h\sigma^2$, so an approximate 95% interval is

With $\sigma = 2$, the half-width is $3.92$ at $h = 1$ and $7.84$ at $h = 4$. Intervals widen with horizon; a model whose intervals do not is suspicious. Other routes to intervals: ETS and ARIMA give analytic intervals under Gaussian errors; quantile regression trains models on the pinball loss for chosen quantiles (gradient boosting supports loss="quantile"); probabilistic deep models output distributions; and conformal prediction calibrates intervals from backtest residuals with coverage guarantees under exchangeability (which time series only approximately satisfy). Always check coverage in the backtest: a nominal 90% interval should contain about 90% of actuals.

import numpy as np

def mae(a, f):   return np.mean(np.abs(a - f))
def smape(a, f): return 100 * np.mean(2 * np.abs(a - f) / (np.abs(a) + np.abs(f)))
def mase(a, f, train, m):
    return mae(a, f) / np.mean(np.abs(train[m:] - train[:-m]))

def seasonal_naive(train, h, m):
    return np.array([train[-m:][i % m] for i in range(h)])

rng = np.random.default_rng(0)
n, m, h = 156, 12, 12
t = np.arange(n)
y = 100 + 0.5 * t + 10 * np.sin(2 * np.pi * t / m) + rng.normal(0, 3, n)

scores = []
for origin in range(n - 4 * h, n, h):            # 4 rolling origins, expanding window
    train, test = y[:origin], y[origin:origin + h]
    f = seasonal_naive(train, h, m)
    scores.append((mae(test, f), smape(test, f), mase(test, f, train, m)))
print("MAE=%.2f sMAPE=%.2f%% MASE=%.2f" % tuple(np.mean(scores, axis=0)))

This backtest prints roughly MAE 6.0, sMAPE 3.7% and MASE 0.95 for seasonal naive, and the plain naive method scores MASE about 1.16 on the same origins. Any new model has to beat those numbers on the same origins to be worth deploying.

Common pitfalls

• Leakage through features. A rolling mean that includes the current value, or a feature computed on the whole dataset (scaling, target encoding) before splitting, quietly uses the future.
• Random cross-validation. Shuffled K-fold on time series gives optimistic scores. Use rolling origins.
• No baseline. Without seasonal naive in the table, a MAPE of 12% means nothing.
• MAPE on intermittent demand. Percentage errors explode or are undefined when actuals are zero; use MASE or a scaled error.
• Ignoring the horizon. Reporting one-step-ahead accuracy for a model used four weeks ahead.
• Trees and trend. Raw-level gradient boosting flattens out at the highest value seen in training.
• Overconfident intervals. Residual-based intervals ignore parameter and model uncertainty; check empirical coverage.
• Structural breaks. COVID-era data, a price change or a new competitor can make history misleading; flag or down-weight those periods.

In practice

Start with plots, a decomposition and the ACF. Build seasonal naive, ETS and (for few series) ARIMA baselines with a library such as statsforecast or R's fable. For many series with covariates, build a global gradient-boosting model with lag, rolling and calendar features using a framework that handles the bookkeeping safely (for example mlforecast, skforecast or darts). Add a pretrained foundation model as a zero-shot candidate. Compare everything in one rolling-origin backtest at the business horizon, with MASE for accuracy and coverage for intervals, and pick per segment if needed. In production, forecast quantiles rather than points, retrain on a schedule, monitor error drift against the baseline, and keep a human override for known events the model cannot see. Simple ensembles (the average of a few good models) are a reliable last step.