wenmin-wu/timeseries-recursive-multistep-forecasting

Overview

Direct multi-step forecasting (one model per horizon day) is expensive: 28 days = 28 models. Recursive forecasting trains one one-step model and reuses it for every horizon day by feeding its own predictions back as inputs. The trick is that lag features (sales.shift(7)) and rolling features (sales.rolling(28).mean()) at horizon day H+k depend on predictions from days H..H+k-1, so you must recompute them inside the prediction loop, not once before. Done correctly, you get 28-day forecasts with the same model that scored well on 1-step validation. Done wrong (precomputed features that don't see the predictions), you get garbage.

Quick Start

import pandas as pd
from datetime import timedelta

base_test = build_panel_with_unknown_target()   # rows for entire horizon, sales=NaN

for h in range(1, 29):                          # 28-day horizon
    day = first_forecast_day + timedelta(days=h - 1)

    # window includes max_lags days before today so we can recompute features
    window = base_test[
        (base_test.date >= day - timedelta(days=max_lags)) &
        (base_test.date <= day)
    ].copy()

    create_features(window)                     # lags + rollings using up-to-day data

    today = window.loc[window.date == day, train_cols]
    yhat = alpha * model.predict(today)         # alpha = bias-correction multiplier

    base_test.loc[base_test.date == day, 'sales'] = yhat

Workflow

1. Build a panel that already contains the future horizon rows with sales=NaN
2. Loop h over the horizon days
3. Slice a sliding window covering [day - max_lags, day] so feature creation has enough history
4. Recompute lag and rolling features on the window — this is the step everyone misses
5. Predict for the current day, optionally apply a bias-correction multiplier (alpha ≈ 1.02-1.03 for Poisson)
6. Write the prediction back into base_test so the next iteration's features see it
7. After the loop, the entire horizon column is filled

Key Decisions

• Recompute features inside the loop: precomputing lag-7 once before the loop means day 8's "lag-7" is actually a previously-predicted day-1 — but you need that prediction to exist before you compute it. The only correct order is predict → write → recompute → predict.
• Window slice for speed: don't recompute features on the entire panel — just on [day - max_lags, day]. The rest is irrelevant for today's prediction.
• Bias-correction multiplier alpha: Poisson and Tweedie LightGBM consistently underpredict by 2-3%; multiply by ~1.02-1.03 (tuned on validation) before writing back, otherwise the bias compounds across the horizon.
• Single-store batches: if you train per-store models, group by store inside the loop to amortize the model load.
• vs. direct multi-output: direct is more accurate for short horizons, recursive scales better and uses one model.

References

• M5 First Public Notebook Under 0.50
• M5 - Three shades of Dark: Darker magic