Forecasting with machine learning: lag features, gradient boosting and multi-step strategies

Before this lesson: lesson 03 (rolling-origin backtesting, MASE, leakage in evaluation) and lesson 05 (an AR model is a regression on lags). Differencing is in lesson 02. Decision trees, gradient boosting, cross-validation and leakage in pipelines are in the ml-fundamentals track. Time: about 90 minutes, plus the notebook.

---

Why this matters #

When you forecast thousands of related series with many covariates (retail sales, electricity load, web traffic), gradient-boosted trees on lag and calendar features are among the strongest methods. In the M5 competition (2020, 42,840 Walmart sales series), "most of the methods examined utilized LightGBM", a gradient-boosting library, and essentially all of the 50 best methods used cross-learning: one model trained on many series. The winning entry averaged LightGBM models trained per store, per store and category, and per store and department, in a recursive and a non-recursive variant. It was 22.4 % more accurate than the best statistical benchmark (Makridakis, Spiliotis and Assimakopoulos 2022). The same "window to table" step is the input format of the deep networks of lesson 09 and of the quantile models of lesson 10.

The risk is just as large. A forecasting table is easy to build with a hidden leak, for example a rolling mean that contains today's value. Then the backtest looks excellent and the model fails in production. This lesson shows how to build the table correctly and how to test it.

1. Predict first #

2. The picture #

From a series to a table #

Take a series $y_1, \dots, y_n$, a lookback of $L$ lags and a horizon $h$. For each target time $s$, one row of the table has

The row pretends that we stand at the forecast origin $s - h$: we know everything up to $y_{s-h}$ and we forecast $h$ steps ahead. Slide the window one step to the right and you get the next row. A lag-$k$ feature of row $s$ is the value $y_{s-k}$, $k$ steps before the target time. A rolling-mean feature is the mean of the last $w$ values that are known at the origin.

Interactive figure: loading…

Which values are known at the origin? #

At the origin $t$ we forecast $y_{t+h}$. The lag-$k$ value of the target time $t + h$ is $y_{t+h-k}$. It is known only if $t + h - k \le t$, that is, if $k \ge h$. So the set of usable lags depends on the horizon. Covariates (other variables) follow the same logic, but some of them are known in advance: the day of the week, a holiday, a price that was planned before the origin.

Interactive figure: loading…

Trees flatten trends #

A regression tree splits the inputs into boxes and predicts, in each box, the mean of the training targets in that box. On a rising series the future inputs (recent lags) are larger than every training input. They all fall into the "largest values" boxes, so the forecast stays near the top of the training range and stops rising. Boosting adds many trees, but the sum is still constant outside the training range of the inputs.

Interactive figure: loading…

3. Definitions and main results #

Reduction to supervised learning #

Let $y_1, \dots, y_n$ be the observed series. Fix a lookback $L \ge 1$ and a horizon $h \ge 1$.

• The information set $\mathcal I_t$ at the origin $t$ is everything known at time $t$: the values $y_1, \dots, y_t$, past values of the observed covariates, the known-in-advance covariates (also for future times) and the static features.
• A feature for forecasting $y_{t+h}$ is a function of $\mathcal I_t$ only. Leakage is a feature that uses anything outside $\mathcal I_t$.
• The lag table for horizon $h$ has one row per target time $s = L + h, \dots, n$: features $(y_{s-h}, \dots, y_{s-h-L+1})$ and target $y_s$. It has $n - L - h + 1$ rows: the first target needs $s - h - L + 1 \ge 1$, so $s \ge L + h$, and the last one is $s = n$.

In pandas, y.shift(k) puts $y_{s-k}$ in row $s$. So the lag-$k$ column is y.shift(k), and every feature that is computed with a window must be computed after a shift by $h$.

Feature	Value in row $s$ (target $y_s$, origin $s - h$)	pandas	Known at the origin?
lag $k$	$y_{s-k}$	y.shift(k)	only if $k \ge h$
rolling mean over $w$	$\frac1w \sum_{j=0}^{w-1} y_{s-h-j}$	y.shift(h).rolling(w).mean()	yes
expanding mean	mean of $y_1, \dots, y_{s-h}$	y.shift(h).expanding().mean()	yes
month, day of week, holiday	calendar of time $s$	index.month, ...	yes (known in advance)
cyclic month	$\sin(2\pi\,\mathrm{month}/12)$, $\cos(2\pi\,\mathrm{month}/12)$	np.sin, np.cos	yes
Fourier terms	$\sin(2\pi j s/m)$, $\cos(2\pi j s/m)$, $j = 1, \dots, K$	NumPy	yes
y.rolling(w).mean()	contains $y_s$ itself	—	no: leakage

Calendar encodings. A tree can split an integer month (1 to 12) directly. A linear model cannot: it would treat December (12) as "12 times January". For linear models use one-hot columns, or the cyclic encoding $(\sin(2\pi k/m), \cos(2\pi k/m))$, which puts the $m$ seasons on a circle. On the circle the distance between seasons $k$ and $k'$ is $2\,|\sin(\pi(k - k')/m)|$, so December and January are as close as January and February ($2\sin(\pi/12) \approx 0.52$). Fourier terms with $K$ pairs give a smooth seasonal shape with $2K$ columns (FPP3 §7.4).

Three kinds of covariates #

A covariate (exogenous variable) is any variable other than the target.

• Static: one value per series, for example the store id or the product category.
• Known in advance: future values are known at the origin, for example the calendar, holidays, planned prices and promotions, or a weather forecast issued before the origin.
• Observed only: values are known only up to the origin, for example the measured temperature, the sales of other products, or the number of visitors.

Proposition 1 (usable lags in the direct strategy). At the origin $t$, for the target $y_{t+h}$:

1. lag $k$ of the target, or of an observed-only covariate, is usable if and only if $k \ge h$;
2. a known-in-advance covariate is usable at every lag $k \ge 0$, including its value at the target time $t + h$;
3. a static covariate is always usable.

Proof. Lag $k$ of the target time $t + h$ is the value at time $t + h - k$. An observed-only variable is in $\mathcal I_t$ exactly at times $\le t$, and $t + h - k \le t \iff k \ge h$. A known-in-advance variable is in $\mathcal I_t$ at all times, and a static one does not depend on time. $\blacksquare$

Trees cannot extrapolate #

Proposition 2.

1. A regression tree fitted with squared loss predicts, in each leaf, the mean of the training targets in that leaf. So every prediction lies in $[\min_i y_i, \max_i y_i]$, the range of the training targets. The same holds for a random forest, which averages trees.
2. Let $F(x) = F_0 + \nu \sum_{m=1}^M T_m(x)$ be a boosted ensemble of trees. For each feature $j$, let $c_j$ be the largest split threshold on feature $j$ in all the trees. Then $F$ is constant on the region $\{x : x_j > c_j \text{ for every } j\}$.

Proof.

1. A leaf mean is an average of training targets. An average of numbers in $[a, b]$ is in $[a, b]$. An average of tree predictions (a forest) is again in $[a, b]$.
2. Each internal node of a tree asks "is $x_j \le c$?" for some $c \le c_j$. If $x_j > c_j$, the answer is always "no". So all points of the region follow the same path in every tree, reach the same leaf in every tree, and get the same sum $F(x)$. $\blacksquare$

For boosting, the predictions can leave the range of the targets a little, but part 2 is what matters: the thresholds are always between training values of the features. On a trending series the future lags are larger than all training lags, so the forecast stops following the trend. The fixes change the target so that the future looks like the past:

1. Differences. Predict $d_s = y_s - y_{s-1}$ from lags of $d$ (lesson 02). The forecast is $\hat y_{t+h\mid t} = y_t + \sum_{j=1}^{h} \hat d_{t+j\mid t}$.
2. Detrending. Fit a trend on the training data (a line, or the level of an exponential smoothing model from lesson 04), let the trees model the residual, and add the trend back.
3. Scaling by a level. Predict $y_{t+h} / \ell_t$, where $\ell_t$ is the mean of the last few known values. Use the log first if the growth is multiplicative (lesson 01).
4. A linear part. Use a linear model, or trees with linear models in the leaves (LightGBM's linear_tree option).

Multi-step strategies #

We want forecasts for $h = 1, \dots, H$ at the origin $t$. Write $x_t = (y_t, \dots, y_{t-L+1})$ for the last $L$ known values.

• Recursive (iterated) strategy: train one model $f$ for $h = 1$. Then $\hat y_{t+1\mid t} = f(y_t, \dots, y_{t-L+1})$, $\hat y_{t+2\mid t} = f(\hat y_{t+1\mid t}, y_t, \dots)$, and so on: each forecast is fed back as an input.
• Direct strategy: train $H$ models, one per horizon: $\hat y_{t+h\mid t} = f_h(x_t)$, using only values known at the origin.
• Multi-output (MIMO) strategy: one model $F: \mathbb R^L \to \mathbb R^H$ predicts the whole vector $(\hat y_{t+1\mid t}, \dots, \hat y_{t+H\mid t})$ at once (a multi-output tree, or a neural network with $H$ outputs).
• DirRec mixes the two: model $f_h$ uses $x_t$ and the forecasts of $f_1, \dots, f_{h-1}$.

Strategy	Models	Rows per model	Errors fed back?	Trained on the $h$-step loss?
recursive	1	$n - L$	yes	no (one-step loss only)
direct	$H$	$n - L - h + 1$	no	yes
multi-output	1	$n - L - H + 1$	no	yes (all $h$ together)

Proposition 3 (recursive versus direct, one lag). Let $y_t$ be stationary with mean 0, variance $\gamma(0)$ and ACF $\rho(k)$. Fit linear models with one lag by least squares on infinitely many data. Then

1. the recursive forecast is $\rho(1)^h\, y_t$ and the direct forecast is $\rho(h)\, y_t$;
2. their mean squared errors satisfy

So the direct forecast is never worse. They are equal when $\rho(h) = \rho(1)^h$ for all $h$, which is the ACF of an AR(1): then the one-lag model is correct.

Proof.

1. For any number $c$, $\mathrm E[(y_{t+h} - c\,y_t)^2] = \gamma(0) - 2c\,\gamma(h) + c^2\gamma(0)$.
2. This parabola in $c$ is smallest at $c = \gamma(h)/\gamma(0) = \rho(h)$: this is the population least-squares coefficient of the direct model for horizon $h$.
3. The one-step model has $c = \rho(1)$. Applying it $h$ times gives $\rho(1)^h y_t$, so the recursive forecast uses $c = \rho(1)^h$.
4. Insert both values of $c$ into step 1 and subtract: $\gamma(0)\bigl[\rho(1)^{2h} - 2\rho(1)^h\rho(h) + \rho(h)^2\bigr] = \gamma(0)\bigl(\rho(1)^h - \rho(h)\bigr)^2$. $\blacksquare$

The proposition is about a misspecified model. When the model is correct and the data are limited, the recursive strategy is often as good or better: it trains one model on the most rows and uses the model structure. The direct strategy needs $H$ models, and its forecasts at neighbouring horizons are not tied together, so the path can look jagged. On the 111 series of the NN5 competition, Taieb et al. (2012) found multi-output strategies the most accurate.

Interactive figure: loading…

Global and local models #

• A local model has its own parameters for each series $i$, fitted on that series only (one ARIMA or one exponential smoothing model per series).
• A global model is one function $f$ with shared parameters, trained on the rows of all series together. Cross-learning in M5 means global models.

A global model sees many more rows: with 20 series of length 24 and $L = 3$, a local AR(3) fits 4 parameters to 21 rows, while one global AR(3) fits 4 parameters to $20 \times 21 = 420$ rows. Montero-Manso and Hyndman (2021) show three things: (1) global methods are not more restrictive than local ones; both can produce the same forecasts, without any assumption that the series are similar; (2) the complexity of local methods grows with the number of series, while that of a global method stays constant, so a global model can afford to be more complex (for example, use more lags); (3) in their experiments, global linear models were competitive with local methods using two orders of magnitude fewer parameters.

Per-series scaling. Series often have very different levels (a small shop sells 10 units a day, a large one 1,000). Without scaling, a tree splits on absolute values and cannot share the shape between series, and the squared loss is dominated by the largest series. Divide each series by a scale $c_i$, for example the mean of its training values (or of the last $L$ values at each origin), fit the global model on $y_{i,s}/c_i$, and multiply the forecast by $c_i$. Compute $c_i$ from training data only.

Interactive figure: loading…

Gradient boosting in practice #

• Libraries. scikit-learn's HistGradientBoostingRegressor, LightGBM, XGBoost and CatBoost all build histogram-based boosted trees; the ideas of this lesson apply to all of them.
• Loss. Squared error estimates the conditional mean; absolute error the median. For counts with many zeros use a Poisson or Tweedie loss (the M5 winner used the Tweedie loss). The quantile (pinball) loss gives quantile forecasts: lesson 10.
• Early stopping on a time-based block. Choose the number of trees on a validation block that comes after the training rows, ideally with a gap of $H$ steps. In scikit-learn, HistGradientBoostingRegressor turns on early stopping by itself when there are more than 10,000 rows (early_stopping="auto") and then holds out a random 10 % of the rows. For time series this mixes future rows into the training data. Set early_stopping=False and choose max_iter on your own time block, or (scikit-learn ≥ 1.7) pass a later block as X_val, y_val to fit.
• Feature importance. Split-gain importance is computed on the training data and favours features with many possible split points; correlated lags share the importance between them. Use permutation importance on a time-based validation block, and do not read importance as a causal effect.
• Evaluation. Use the rolling origin of lesson 03. At each origin, recompute every feature with the data known at that origin and refit the model (or refit every few origins to save time). Choose hyperparameters inside the training part of each fold. Always compare with the benchmarks of lesson 03 (naive, seasonal naive).

A leakage checklist for features #

1. A rolling or expanding window that includes the target time (no shift, or a shift smaller than $h$).
2. Target encoding (for example the mean sales of each store) computed on all data, including the test period.
3. Scaling or normalisation (mean, standard deviation, min–max, per-series scale) fitted on all data.
4. An observed-only covariate used as if it were known in advance (measured weather; promotions decided late; data that were revised after the origin).
5. Centred filters or decompositions of the whole series used as features.
6. A random split for validation or early stopping.

A simple test catches the first five: change all values after an origin $t$ and recompute the features. No feature of a row with origin $\le t$ may change. The notebook uses this test.

Answers to section 1. (1) A flat line, near the largest training targets: trees predict means of training targets (Proposition 2). Differencing or detrending fixes it. (2) No. For $h = 7$ only lags $k \ge 7$ are known at the origin (Proposition 1). In the recursive strategy lag 1 would be a forecast, not data. (3) No. rolling(7).mean() at time $s$ contains $y_s$, the target itself; the model partly copies the answer. Use y.shift(h).rolling(7).mean().

Hierarchical forecasting and reconciliation

Sales often form a hierarchy: products add up to stores, stores to regions, regions to the total. Write the $m$ bottom-level series as a vector $b_t$ and all series (every level) as $y_t = S b_t$, where the summing matrix $S$ has one row per series. For a total $T = A + B$:

$$\begin{pmatrix} T \\ A \\ B \end{pmatrix} = \underbrace{\begin{pmatrix} 1 & 1 \\ 1 & 0 \\ 0 & 1 \end{pmatrix}}_{S} \begin{pmatrix} A \\ B \end{pmatrix}.$$

Forecasts made for each series separately ("base" forecasts $\hat y$) usually do not add up. Three ways to make them coherent (FPP3 chapter 11):

• Bottom-up: forecast only the bottom series and add them: $\tilde y = S \hat b$.
• Top-down: forecast the total and split it with historical proportions.
• MinT (minimum trace, Wickramasuriya, Athanasopoulos and Hyndman 2019): $\tilde y = S (S^\top W^{-1} S)^{-1} S^\top W^{-1} \hat y$, where $W$ is the covariance matrix of the base forecast errors. With $W = I$ this is the orthogonal projection onto the coherent subspace (OLS reconciliation).

Example. Base forecasts: total 100, $A$ 55, $B$ 50; they do not add up ($55 + 50 = 105$). With $W = I$: $S^\top S = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}$ and $S^\top \hat y = (155, 150)$, so $\tilde b = \frac13\begin{pmatrix} 2 & -1 \\ -1 & 2 \end{pmatrix}\begin{pmatrix} 155 \\ 150 \end{pmatrix} = (53.33, 48.33)$ and the total becomes 101.67. The gap of 5 is shared equally: every series moves by $5/3$. Bottom-up would give a total of 105; top-down with proportions $55/105$ and $50/105$ would give $A = 52.38$, $B = 47.62$.

4. Worked example #

Use the short series

$s$	1	2	3	4	5	6	7	8
$y_s$	10	11	13	12	14	16	15	17

(a) One-step table, $L = 3$, $h = 1$. The rows are $s = L + h, \dots, n = 4, \dots, 8$, so $n - L - h + 1 = 8 - 3 - 1 + 1 = 5$ rows. Add the rolling mean of the two previous values, y.shift(1).rolling(2).mean(), which is $(y_{s-1} + y_{s-2})/2$:

$s$	$y_{s-1}$	$y_{s-2}$	$y_{s-3}$	roll2 (correct)	roll2 without shift (leaky)	target $y_s$
4	13	11	10	12.0	12.5	12
5	12	13	11	12.5	13.0	14
6	14	12	13	13.0	15.0	16
7	16	14	12	15.0	15.5	15
8	15	16	14	15.5	16.0	17

The leaky column is $(y_s + y_{s-1})/2$: in row 6 it is $(16 + 14)/2 = 15$, and it contains the target 16. To forecast $y_9$ the model gets the row $(y_8, y_7, y_6) = (17, 15, 16)$ and roll2 $= (17 + 15)/2 = 16$.

(b) Direct table for $h = 2$, $L = 2$. Row $s$ has features $(y_{s-2}, y_{s-3})$ and target $y_s$; there are $8 - 2 - 2 + 1 = 5$ rows:

$s$	$y_{s-2}$	$y_{s-3}$	target $y_s$
4	11	10	12
5	13	11	14
6	12	13	16
7	14	12	15
8	16	14	17

At the origin 8 this model forecasts $y_{10}$ from $(y_8, y_7) = (17, 15)$. It never sees $y_9$.

(c) A tree on the level versus a tree on the differences. Fit a depth-1 tree (a stump) on the one-step table with the single feature $y_{s-1}$. The 7 pairs $(y_{s-1} \to y_s)$ are $10 \to 11$, $11 \to 13$, $13 \to 12$, $12 \to 14$, $14 \to 16$, $16 \to 15$, $15 \to 17$. The best split is $y_{s-1} \le 13.5$:

• left leaf ($y_{s-1} \in \{10, 11, 13, 12\}$): targets 11, 13, 12, 14, mean $12.5$, squared error $1.5^2 + 0.5^2 + 0.5^2 + 1.5^2 = 5$;
• right leaf ($y_{s-1} \in \{14, 16, 15\}$): targets 16, 15, 17, mean $16$, squared error $0 + 1 + 1 = 2$.

The total squared error $7$ is the smallest of all splits. Forecast recursively from $y_8 = 17$: $17 > 13.5$ gives $\hat y_9 = 16$; then $16 > 13.5$ gives $\hat y_{10} = 16$; and so on. The forecast is flat at 16, below the last observed value, although the series rises by about 1 per step.

Now model the differences $d_s = y_s - y_{s-1}$: $1, 2, -1, 2, 2, -1, 2$, with mean $(17 - 10)/7 = 1$. The simplest tree (a single leaf) predicts the mean difference, $\hat d = 1$. Adding it up: $\hat y_9 = 18$, $\hat y_{10} = 19$, $\hat y_{11} = 20$. This is the drift method of lesson 03; a deeper tree on lags of $d$ would add corrections to it.

(d) Recursive versus direct with numbers. Take the AR(2) process $y_t = 1.2\,y_{t-1} - 0.5\,y_{t-2} + \varepsilon_t$, $\sigma^2 = 1$, and the (wrong) one-lag model. From the Yule–Walker equations (lesson 05):

• $\rho(1) = \phi_1/(1 - \phi_2) = 1.2/1.5 = 0.8$;
• $\rho(2) = \phi_1\rho(1) + \phi_2 = 0.96 - 0.5 = 0.46$;
• $\rho(3) = \phi_1\rho(2) + \phi_2\rho(1) = 0.552 - 0.4 = 0.152$;
• $\gamma(0) = \dfrac{(1 - \phi_2)\,\sigma^2}{(1 + \phi_2)\bigl((1 - \phi_2)^2 - \phi_1^2\bigr)} = \dfrac{1.5}{0.5 \cdot 0.81} = 3.704$.

At $h = 2$: the recursive coefficient is $0.8^2 = 0.64$; the direct coefficient is $\rho(2) = 0.46$. By step 1 of the proof,

The difference $3.704 \cdot (0.64 - 0.46)^2 = 3.704 \cdot 0.0324 = 0.12$ agrees with Proposition 3. At $h = 3$: $3.704\,(1 - 2 \cdot 0.512 \cdot 0.152 + 0.512^2) = 4.10$ versus $3.704\,(1 - 0.152^2) = 3.62$. The recursive error at $h = 3$ is even larger than $\gamma(0) = 3.70$, the error of forecasting the mean 0: the wrong model applied three times is worse than no model.

5. In code #

pandas builds the table with shift; scikit-learn's tree confirms part (c); a long simulation confirms part (d).

import numpy as np
import pandas as pd
from sklearn.tree import DecisionTreeRegressor

# (a) series -> table: lags 1..3 and a rolling mean of the 2 previous values (h = 1)
y = pd.Series([10, 11, 13, 12, 14, 16, 15, 17.], index=range(1, 9))
X = pd.DataFrame({f"lag{k}": y.shift(k) for k in (1, 2, 3)})
X["roll2"] = y.shift(1).rolling(2).mean()       # shift first: uses y_{s-1}, y_{s-2}
X["roll2_leak"] = y.rolling(2).mean()           # WRONG: contains the target y_s
print(pd.concat([X, y.rename("target")], axis=1).dropna().to_string())
#    lag1  lag2  lag3  roll2  roll2_leak  target
# 4  13.0  11.0  10.0   12.0        12.5    12.0
# 5  12.0  13.0  11.0   12.5        13.0    14.0
# 6  14.0  12.0  13.0   13.0        15.0    16.0
# 7  16.0  14.0  12.0   15.0        15.5    15.0
# 8  15.0  16.0  14.0   15.5        16.0    17.0

# (c) a depth-1 tree (stump) on lag 1: the forecast stays flat
stump = DecisionTreeRegressor(max_depth=1).fit(X[["lag1"]].iloc[1:], y.iloc[1:])
print(stump.tree_.threshold[0], stump.tree_.value[1:, 0, 0])   # 13.5 [12.5 16. ]
x, path = 17.0, []
for h in range(3):                                # recursive: feed the forecast back
    x = stump.predict(pd.DataFrame({"lag1": [x]}))[0]
    path.append(float(x))
print("raw target:", path)                                       # [16.0, 16.0, 16.0]
print("differences:", y.iloc[-1] + np.cumsum([np.diff(y).mean()] * 3))   # [18. 19. 20.]

# (d) recursive vs direct for a misspecified AR(1) model on AR(2) data
rng = np.random.default_rng(0)
z = np.zeros(200_000)
e = rng.normal(size=z.size)
for t in range(2, z.size):
    z[t] = 1.2 * z[t - 1] - 0.5 * z[t - 2] + e[t]
z = z[1000:]
def ls_coef(h):                                   # least squares of z_{t+h} on z_t
    return np.dot(z[:-h], z[h:]) / np.dot(z[:-h], z[:-h])
a = ls_coef(1)
for h in (2, 3):
    b = ls_coef(h)
    mse_rec = np.mean((z[h:] - a**h * z[:-h])**2)
    mse_dir = np.mean((z[h:] - b * z[:-h])**2)
    print(h, round(a**h, 3), round(b, 3), round(mse_rec, 2), round(mse_dir, 2))
# 2 0.639 0.458 3.04 2.92
# 3 0.51 0.15 4.09 3.61

The simulated values match the theory: $0.64$ and $0.46$ at $h = 2$, $0.512$ and $0.152$ at $h = 3$, and MSEs $3.04, 2.92, 4.10, 3.62$ up to simulation noise. The notebook builds the lag table with a leak test, fits HistGradientBoostingRegressor on raw and differenced targets, compares the strategies in a backtest, trains a global model, and repairs a leaky pipeline.

6. Where it shows up #

• Retail demand. M5 (Walmart, 2020): most methods used LightGBM on lags, rolling statistics, prices and calendar events, trained across many series. The winners of the earlier Kaggle competitions "Corporación Favorita Grocery Sales Forecasting" and "Recruit Restaurant Visitor Forecasting" also built on LightGBM (as reported in the M5 paper).
• Energy. Electricity load forecasts use lags of the load, calendar features (hour, day of week, holidays) and temperature. In production the temperature of the forecast period is a weather forecast: a known-in-advance covariate only if it was issued before the origin.
• Software. Libraries implement the reduction for you: sktime's make_reduction (strategies "recursive", "direct", "multioutput"), skforecast, and Nixtla's mlforecast. They still need you to decide which features are known at the origin.
• ML for science. In hydrology, Kratzert et al. (2019) trained one LSTM on 531 river basins of the CAMELS data set, with meteorological series and static basin attributes as inputs. This global model was more accurate than hydrological models calibrated basin by basin. The same global-model idea is used for air-quality stations and other sensor networks.
• Deep learning. The window-to-table step, per-series scaling and the multi-output strategy are the standard set-up of the neural forecasters in lesson 09.

7. Common mistakes #

• A rolling feature without a shift. y.rolling(7).mean() contains $y_s$. Use y.shift(h).rolling(7).mean(), with the shift equal to the horizon.
• The same lags for every horizon. In the direct strategy the model for horizon $h$ can use only lags $k \ge h$. A table built for $h = 1$ leaks when it is reused for $h = 7$.
• Trees on a trending level. The forecast goes flat. Model differences, ratios or residuals from a trend instead, and check the forecast against the drift benchmark.
• Statistics computed on all data. Target encodings, per-series scales and normalisation constants must be computed from the training part of each fold only.
• Random validation split. A shuffled hold-out (also inside early stopping) puts future rows into training. Validate on a later time block.
• Measured covariates as if they were known. Train with what was known at the origin, for example a weather forecast, not the measured weather.
• One global model on raw levels. Large series dominate the loss and patterns are not shared. Scale each series first and undo the scaling after the forecast.

8. Check yourself #

1. A series has $n = 100$ values. How many rows has the lag table with $L = 12$ lags for the horizon $h = 3$?
   Answer

   $n - L - h + 1 = 100 - 12 - 3 + 1 = 86$. The first target time is $s = L + h = 15$, the last is $100$.
2. Write the pandas code for a 28-day rolling mean feature for a direct model at horizon $h = 7$. Which values does it use for the target $y_s$?
   Answer

   `y.shift(7).rolling(28).mean()`. It is the mean of $y_{s-7}, y_{s-8}, \dots, y_{s-34}$: 28 values, the most recent of which is at the origin $s - 7$.
3. The ACF of a series is $\rho(1) = 0.5$, $\rho(2) = 0.6$. With one lag, what are the recursive and the direct coefficients at $h = 2$, and how large is the gap in MSE?
   Answer

   Recursive: $\rho(1)^2 = 0.25$. Direct: $\rho(2) = 0.6$. Gap: $\gamma(0)(0.25 - 0.6)^2 = 0.1225\,\gamma(0)$. The one-lag model misses the strong lag-2 dependence, so iterating it is a poor 2-step forecast.
4. Which of these features are leaks for a one-day-ahead forecast of today's sales made last night: (i) yesterday's sales; (ii) today's planned price; (iii) today's measured temperature; (iv) the store's mean sales over the whole data set, test period included?
   Answer

   (iii) and (iv). Yesterday's sales are known (lag 1 with $h = 1$). A planned price is known in advance. Today's measured temperature is observed only after the origin (use the weather forecast). The store mean over all data includes the test targets (compute it on the training part only).
5. Why does a global model need per-series scaling, while a local model does not?
   Answer

   A local linear model with an intercept adapts to its own series' level, and a local tree splits on that series' own values. A global model shares one function across series: on raw values the same split or intercept means different things for different series, and the loss is dominated by the series with the largest scale. Scaling puts all series on one scale so the shared patterns can be learned.
6. You fit HistGradientBoostingRegressor() with default settings on 50,000 rows of a lag table. What is wrong with the early stopping, and how do you fix it?
   Answer

   With more than 10,000 rows, `early_stopping="auto"` holds out a random 10 % of rows, so the validation rows are mixed in time with the training rows. Set `early_stopping=False` and choose `max_iter` on a later time block, or pass a later block as `X_val, y_val` (scikit-learn ≥ 1.7).

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. What is "reduction to supervised learning" in forecasting?

   Show answer

2. How many rows has a lag table with $n$ values, $L$ lags and horizon $h$?

   Show answer

3. Write a leak-free 7-step rolling mean feature for horizon $h$ in pandas.

   Show answer

4. In the direct strategy, which lags $k$ of the target are usable at horizon $h$?

   Show answer

5. Name the three kinds of covariates.

   Show answer

6. Why can a tree ensemble not extrapolate a trend?

   Show answer

7. Give two fixes for trees on a trending series.

   Show answer

8. Recursive versus direct: what is the main trade-off?

   Show answer

9. For a one-lag linear model, what is $\mathrm{MSE}_{\text{rec}}(h) - \mathrm{MSE}_{\text{dir}}(h)$?

   Show answer

10. Why does a global model need per-series scaling?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "How do you forecast with XGBoost or LightGBM?" — Build a lag table: lags $\ge h$, rolling statistics after shift(h), calendar and known-in-advance covariates. Choose a strategy (recursive, direct, multi-output). Make the target stationary (differences or scaling) because trees do not extrapolate. Backtest with a rolling origin against naive benchmarks.
• "Your gradient-boosting forecast is flat on a growing series. Why?" — Trees predict means of training targets; inputs beyond the training range all land in the same leaves. Model differences or detrend.
• "Recursive or direct?" — Recursive: one model, errors fed back, good when the model is right. Direct: one model per $h$, robust to misspecification, more models and fewer rows. Multi-output predicts the whole path at once. Know the formula $\gamma(0)(\rho(1)^h - \rho(h))^2$.
• "Your backtest is much better than production. What do you check?" — Leaks: rolling windows without a shift, encodings or scalers fitted on all data, covariates that are only observed, random validation splits; then the perturbation test (change the future, features must not move).
• "One model for 10,000 products or 10,000 models?" — Usually one global model with per-series scaling (cross-learning, as in M5): more data per parameter and shared patterns; add series-level features (category, store) for differences between series.

10. Practice #

The notebook for this lesson runs in Colab. It has five exercises, each with a check cell and a folded solution:

1. Write make_lag_table(y, lags, horizon, rolling) with pandas; it is checked on the series of section 4 and with a leak test that changes the future.
2. Fit HistGradientBoostingRegressor on lags of the "monthly demand" series, recursively, on the raw and on the differenced target; compare the test RMSE.
3. Write least-squares recursive and direct forecasters and compare their RMSE by horizon on AR(2) data with one lag and with two lags.
4. Train one global linear model on 30 related series with per-series scaling and compare its MASE with 30 local AR models.
5. Find and fix the three leaks in a feature pipeline (unshifted rolling mean, target encoding on all data, a same-day measured covariate).

11. Watch and read #

Kishan Manani - Feature Engineering for Time Series Forecasting | PyData London 2022PyData · 43 min

Lag and window features, recursive and direct forecasting, and why trees need a detrended target.

Gradient Boost Part 1 (of 4): Regression Main IdeasStatQuest with Josh Starmer · 16 min

How boosting adds small trees fitted to residuals; useful before Proposition 2.

HSMA: Introduction to autoregressive models + iterative and direct forecastingHSMA · 39 min

A lecture on sliding windows and on the iterative (recursive) and direct strategies.

Kaggle Winning Solution Series: Walmart M5 Forecasting – AccuracyYanAITalk · 18 min

How the M5 winners built features and trained LightGBM models across series.

• Makridakis, Spiliotis and Assimakopoulos (2022), M5 accuracy competition: Results, findings, and conclusions, International Journal of Forecasting 38(4): LightGBM, cross-learning and the winning methods.
• Ben Taieb, Bontempi, Atiya and Sorjamaa (2012), A review and comparison of strategies for multi-step ahead time series forecasting based on the NN5 forecasting competition, Expert Systems with Applications.
• Montero-Manso and Hyndman (2021), Principles and algorithms for forecasting groups of time series: Locality and globality, International Journal of Forecasting.
• Bergmeir, Hyndman and Koo (2018), A note on the validity of cross-validation for evaluating autoregressive time series prediction, Computational Statistics and Data Analysis.
• Hyndman and Athanasopoulos, Forecasting: Principles and Practice (3rd ed.), free online: time series cross-validation and chapter 11, hierarchical forecasting.
• scikit-learn examples: Lagged features for time series forecasting and Time-related feature engineering (cyclic encodings, splines and HistGradientBoostingRegressor).
• Kratzert et al. (2019), Towards learning universal, regional, and local hydrological behaviors via machine learning applied to large-sample datasets, Hydrology and Earth System Sciences: one global LSTM for 531 basins.