Forecasting models, from simple to learned

You know how to validate a forecast honestly. This is what to put inside it. The classical methods, exponential smoothing and ARIMA, are fast, interpretable and genuinely hard to beat on many business series, so they are a baseline rather than a formality. The alternative is to reshape the problem into ordinary supervised learning with lag features and fit a tree, which handles many related series and external variables easily. The section also covers the horizon, since accuracy one day ahead says little about thirty, and giving an interval rather than a single number nobody should trust.

What you have                                  Start with
one series, a few hundred points               exponential smoothing or ARIMA
one series, external drivers you know ahead    gradient boosting on lag features
thousands of related series                    one global gradient boosting model
thousands of series, years of history, a team  a deep global model, and prove it
demand that is zero most days                  an intermittent-demand method
a series you cannot explain at all             the seasonal naive baseline, honestly

Classical models#

Exponential smoothing weights recent observations more heavily and handles trend and seasonality. ARIMA models a series through its own past values and past errors. Both are strong, fast and interpretable baselines, and on many business series they are hard to beat.

Simple exponential smoothing is one line, and writing it out shows exactly what the single parameter does. The level is a running blend of what just happened and what you already believed.

import numpy as np

def ses(series, alpha):
    """level_t = alpha * observation + (1 - alpha) * level_{t-1}"""
    level = series[0]
    out = [level]
    for x in series[1:]:
        level = alpha * x + (1 - alpha) * level
        out.append(level)
    return np.array(out)

s = np.array([10.0, 12, 11, 40, 13, 12, 11])     # a series with one spike
print(ses(s, 0.3).round(2))
print(ses(s, 0.8).round(2))

# [10.   10.6  10.72 19.5  17.55 15.89 14.42]   alpha 0.3: slow, ignores the spike
# [10.   11.6  11.12 34.22 17.24 13.05 11.41]   alpha 0.8: fast, chases the spike

The fuller versions add a trend term and a seasonal term with their own smoothing parameters, which is Holt-Winters. ARIMA arrives at similar forecasts from a different direction: differencing the series until it is stationary, then modelling what is left as a combination of its own past values and its own past errors. The practical point is not which is theoretically nicer but that both fit in milliseconds, both produce prediction intervals for free from their error model, and both are hard to beat on a single business series with a few hundred points.

Treat it as supervised learning#

Build lag and rolling features, add calendar features, and fit a gradient-boosted tree. This handles many related series at once and incorporates external variables easily, which is why it is common in practice.

date        sales     ->   lag_1  lag_7  roll_7  dow  holiday  promo   y
2026-03-08   120            118    115    117.4    6     0        0     120
2026-03-09   118            120    119    117.9    0     0        1     118
2026-03-10   131            118    122    118.6    1     0        1     131

One row per date, every feature known before that date. Now it is a table,
and everything you know about tabular modelling applies.

Gradient boosting on that table is the honest default for applied work, for four reasons. It takes external variables — price, promotion, weather, a competitor's opening — which classical models handle awkwardly. It fits thousands of series in one model, so a product with six weeks of history borrows the shape of everything similar. It handles non-linear effects and interactions without being told they exist. And it is fast enough to backtest properly, which matters more than any single point of accuracy, because a model you can only evaluate once is a model you cannot trust.

import numpy as np
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(0)
n = 300
t = np.arange(n).astype(float)
y = 50 + 0.4 * t + 6 * np.sin(2 * np.pi * t / 7) + rng.normal(0, 3, n)

X = np.column_stack([t % 7, np.sin(2 * np.pi * t / 7), np.cos(2 * np.pi * t / 7)])
train, test = slice(0, 250), slice(250, 300)

def rmse(a, b):
    return round(float(np.sqrt(np.mean((a - b) ** 2))), 2)

gbm = GradientBoostingRegressor(random_state=0).fit(X[train], y[train])
print("GBM on the level    ", rmse(y[test], gbm.predict(X[test])))

trend = LinearRegression().fit(t[train].reshape(-1, 1), y[train])
resid = y - trend.predict(t.reshape(-1, 1))
gbm2 = GradientBoostingRegressor(random_state=0).fit(X[train], resid[train])
pred = trend.predict(t[test].reshape(-1, 1)) + gbm2.predict(X[test])
print("GBM on the residual ", rmse(y[test], pred))
print("the trend alone     ", rmse(y[test], trend.predict(t[test].reshape(-1, 1))))

# GBM on the level     59.81
# GBM on the residual   3.11
# the trend alone       4.93

When a deep model earns its place#

Not on one series with three years of monthly data. The forecasting competitions that made global deep models credible were run on tens of thousands of related series with years of history each, and that is the setting where the extra capacity has something to learn from. Below roughly that scale, a deep model is a slower, less interpretable way of arriving at a worse number.

• Many related series, so the model can learn shape from thousands of examples rather than one. This is the main condition, and it is the one people skip.
• Inputs that do not fit in a table: raw text about a product, an image of a store shelf, a long history where the relevant pattern is a hundred steps back.
• A need for the whole horizon at once. Sequence models output all thirty days in one pass, which avoids the compounding errors of feeding forecasts back in.
• Enough traffic that a point of accuracy is worth a retraining pipeline, a GPU and someone who understands it at three in the morning.

The bar to clear is a tuned gradient boosting model on the same features, backtested over the same rolling origins, at the same horizons. That comparison is easy to arrange and frequently not arranged, which is why published deep forecasting results and applied experience disagree as often as they do.

Horizon matters#

Forecasting tomorrow is easier than forecasting next month, and errors grow with the horizon. Evaluate at every horizon you will actually use rather than reporting a single averaged number that hides how quickly accuracy decays.

There are two ways to get a multi-step forecast and the choice has consequences. The recursive way forecasts one step, appends the prediction, and forecasts again; it needs one model and it compounds its own errors, so by step thirty it may be forecasting from a history it invented. The direct way trains a separate model per horizon, each predicting h steps ahead from features available now; it costs h models and does not compound, and it is usually the better answer once the horizon is more than a few steps.

                 Recursive                   Direct
models            1                          one per horizon
error behaviour   compounds with each step   independent per step
features          can use lag_1              only lags at least h old
consistency       forecasts flow smoothly    steps can disagree with each other
use when          short horizons             the horizon is the thing you care about

Intervals, not points#

Coverage is a measurement, not an assumption. Take every backtest fold, count how often the actual value fell inside the stated interval, and compare that with the number on the label. An eighty percent interval that contains sixty percent of outcomes is not a cautious forecast, it is a wrong one, and the planner who sized their safety stock from it will find out in a way you do not want to be present for.

import numpy as np

def coverage(actual, low, high):
    return float(np.mean((actual >= low) & (actual <= high)))

rng = np.random.default_rng(1)
actual = rng.normal(100, 10, 500)              # the truth: standard deviation 10

honest = coverage(actual, 100 - 1.2816 * 10, 100 + 1.2816 * 10)
optimistic = coverage(actual, 100 - 1.2816 * 6, 100 + 1.2816 * 6)
print(round(honest, 3), round(optimistic, 3))

# 0.852 0.616
# The second interval is labelled 80% and covers 62%. Nothing in the point
# forecast changed; only the claimed uncertainty was wrong.

Intervals from a model's own error assumptions are usually too narrow, because they account for the noise and not for the model being wrong about the level. Two fixes work in practice. Fit quantile models directly, one per quantile, so the interval is learned rather than derived. Or take the residuals from the backtest, use their empirical quantiles as the interval width, and widen with the horizon — the crude version of conformal prediction, and considerably better than a formula that assumes the model is correct.