Forecasting with XGBoost: turning a time series into a supervised problem without cheating

LEARN · GRADIENT BOOSTING WITH XGBOOST

XGBoost does not know what a time series is.

It accepts a matrix of rows and columns, then learns a mapping from predictors to a target. Forecasting with gradient-boosted trees therefore begins by converting historical timestamps into supervised-learning examples.

That conversion is also where most serious forecasting mistakes happen.

A model can report excellent validation scores while quietly using information that would not exist when a real forecast is issued. Random cross-validation, rolling averages that include the target, retrospectively recorded promotions, and realized future weather can all turn an honest forecasting problem into a disguised reconstruction problem.

This tutorial builds a leakage-safe daily retail-sales forecast with:

Lagged observations.
Past-only rolling statistics.
Calendar and Fourier seasonality features.
Planned promotions, prices, and holidays as exogenous predictors.
Recursive and direct multi-step strategies.
Expanding-window backtesting.
A seasonal-naive baseline.
Seasonal differencing as an extrapolation workaround.
Quantile objectives for prediction intervals.
Automatic recording of package versions and backtest results.

The example uses a generated retail series, so it runs without downloads, credentials, or sensitive data.

The supervised-learning view of forecasting

For a one-step-ahead model, one row can represent this relationship:

today's sales = f(
    yesterday's sales,
    sales seven days ago,
    recent average sales,
    promotion today,
    price today,
    day of week
)

In compact plain-text notation:

y_t = f(y_t−1, y_t−7, rolling_mean_t−1, promo_t, price_t, calendar_t)

The critical rule is that every predictor for y_t must be available before or at time t, according to the real operational forecasting process.

Suppose the forecast is issued at the end of Sunday.

At that origin:

Monday’s date and day of week are known.
A promotion scheduled for Monday may be known.
Monday’s committed price may be known.
Sunday’s sales may be known, depending on reporting latency.
Monday’s realized sales are not known.
Monday’s realized weather is not known.
A seven-day average containing Monday’s sales is not known.

The model is usually the easy part. Reconstructing what was genuinely knowable at each historical forecast origin is the real forecasting work.

Define the forecast contract before writing features

A forecasting project needs more than a target column. It needs an explicit contract.

Record at least:

Forecast origin: When is the forecast issued?
Horizon: How many future periods are required?
Target availability: How late do final target observations arrive?
Exogenous availability: Which future predictors are known, planned, estimated, or unavailable?
Retraining policy: Is the model refitted daily, weekly, or only after a trigger?
Decision metric: Does the business care about average error, large misses, service levels, or interval coverage?

Two teams can use the same historical table and still have different valid feature sets.

A forecast issued Monday at 06:00 may not be allowed to use Sunday’s finalized sales if the accounting pipeline completes at noon. Timestamp order alone does not determine availability.

What “without cheating” means

A leakage-safe experiment obeys four broad rules.

Rolling features must stop before the target

This calculation is wrong for predicting the current row:

df["rolling_mean_7"] = df["sales"].rolling(7).mean()

For the row representing day t, that average includes sales_t, the value being predicted.

Shift first, then roll:

df["rolling_mean_7"] = (
    df["sales"]
    .shift(1)
    .rolling(7)
    .mean()
)

Now the row for day t uses only t−1 through t−7.

The order matters:

target series → shift into the past → calculate rolling feature

Not:

target series → calculate rolling feature → hope the split prevents leakage

Validation observations must occur after training observations

A random split places past and future rows in both partitions.

Even when the target is not directly leaked, adjacent rows often share almost identical:

Lags.
Rolling windows.
Prices.
Promotion states.
Calendar context.
Long-running trends.

A random split answers:

Can the model interpolate among shuffled historical states?

Forecasting requires a different question:

Can a model trained only through this cutoff predict observations that happen afterward?

Use expanding or rolling temporal windows.

Exogenous variables must be available at the forecast origin

A future variable is not automatically valid merely because it appears in the completed historical dataset.

Commonly valid future features include:

Day of week.
Month.
Public holidays.
Store opening schedules.
Promotions approved before the forecast.
Prices committed before the forecast.

Potentially invalid future features include:

Realized weather.
Final competitor prices.
Future inventory calculated after deliveries occur.
Promotions entered retrospectively.
Stockout flags derived from future sales.
Revised business plans that were not available historically.

For weather-driven demand, the correct historical predictor is usually the archived weather forecast available at each origin, not the weather that was eventually observed.

Learned transformations belong inside each fold

Anything estimated from data can leak.

That includes:

Trend coefficients.
Target normalization.
Feature selection.
Lag selection.
Hyperparameter tuning.
Missing-value imputation.
Quantile calibration.
Outlier thresholds.

A global trend fitted before splitting has already seen the validation future.

Each fold must estimate its own transformations using only that fold’s training history.

Recursive versus direct multi-step forecasting

A one-step regressor is not yet a 28-day forecasting system.

You must decide how it will produce multiple future values.

Strategy	Models trained	How day 28 is produced	Main risk
Recursive	1	Predict days 1 through 28 sequentially	Earlier errors become later features
Direct	28	A separate model predicts each horizon	Higher compute and less data per horizon
Hybrid	Several	Combines direct and recursive steps	More implementation complexity

Recursive forecasting

A recursive model predicts day 1, appends that prediction to the history, recalculates lag features, and predicts day 2.

The loop continues until the entire horizon is produced.

Advantages:

Only one model is maintained.
All one-step training rows contribute to the same estimator.
Training is comparatively inexpensive.
Arbitrary horizons can be requested.

Disadvantages:

Training uses real historical lags.
Inference eventually uses the model’s own predictions as lags.
An early error can influence several later predictions.
Long-horizon forecasts may flatten, drift, or become unstable.

The difference between training on observed lags and predicting from estimated lags is often called exposure bias.

Direct forecasting

A direct system trains one model for every required horizon:

model_1 predicts y_t+1
model_2 predicts y_t+2
model_7 predicts y_t+7
model_28 predicts y_t+28

Each model receives information available at forecast origin t, plus known future attributes for its target date.

Advantages:

Day 28 does not consume predictions for days 1 through 27.
Every model can learn horizon-specific relationships.
Direct quantile forecasts have a clear interpretation.
Recursive error propagation is avoided.

Disadvantages:

A 28-day horizon requires 28 models.
The horizon-28 model loses the final 28 possible target rows.
Neighboring horizon predictions are not guaranteed to be smooth.
Training, tuning, storage, and governance cost more.

Neither strategy wins by definition. Backtest both.

Create a current, reproducible environment

The following top-level pins reflect current stable releases on August 4, 2026:

numpy==2.5.1
pandas==3.0.5
scikit-learn==1.9.0
xgboost==3.4.0
matplotlib==3.11.1

NumPy 2.5.1, pandas 3.0.5, scikit-learn 1.9.0, XGBoost 3.4.0, and Matplotlib 3.11.1 were the current published versions at the time of writing. XGBoost 3.4.0 requires Python 3.12 or newer.

Save the package list as requirements.lock, then create the environment:

python --version
python -m venv .venv
source .venv/bin/activate
python -m pip install --upgrade pip
python -m pip install -r requirements.lock

On Windows PowerShell, activate the environment with:

.venv\Scripts\Activate.ps1

Release-day documentation can briefly lag package publication. An earlier review snapshot saw XGBoost 3.4.0 documentation while PyPI still listed 3.3.0. As of August 4, 2026, PyPI publishes XGBoost 3.4.0; Read the Docs may temporarily show the stable branch as 3.3.0 and the latest branch as 3.4.0 development documentation while its channels update. Use PyPI metadata as the source of truth for the installable release.

Every experiment should print its exact runtime versions rather than relying only on the lockfile:

from __future__ import annotations

import json
import platform
from importlib.metadata import version


def version_report() -> dict[str, str]:
    packages = [
        "numpy",
        "pandas",
        "scikit-learn",
        "xgboost",
        "matplotlib",
    ]

    report = {
        "python": platform.python_version(),
        "platform": platform.platform(),
    }

    report.update(
        {
            package: version(package)
            for package in packages
        }
    )

    return report


print(
    json.dumps(
        version_report(),
        indent=2,
        sort_keys=True,
    )
)

Keep this output with the backtest scores. A table of metrics without its package versions, cutoff dates, seed, and feature definition is not a fully reproducible result.

Import the libraries and define the feature contract

Run the remaining blocks in order in a Python script or notebook.

from __future__ import annotations

import json
from pathlib import Path
from typing import Any

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn.metrics import (
    mean_absolute_error,
    mean_squared_error,
)
from xgboost import XGBRegressor


LAGS = (1, 7, 14, 28)
ROLL_WINDOWS = (7, 28)

CALENDAR_COLS = [
    "dow",
    "month",
    "is_weekend",
    "year_sin_1",
    "year_cos_1",
    "year_sin_2",
    "year_cos_2",
]

EXOG_COLS = [
    "promo",
    "price",
    "holiday",
]

RECURSIVE_FEATURES = [
    *EXOG_COLS,
    *CALENDAR_COLS,
    *(f"lag_{lag}" for lag in LAGS),
    *(
        f"roll_mean_{window}"
        for window in ROLL_WINDOWS
    ),
    *(
        f"roll_std_{window}"
        for window in ROLL_WINDOWS
    ),
]

DIRECT_HISTORY_FEATURES = [
    "lag_0",
    "lag_1",
    "lag_7",
    "lag_14",
    "lag_28",
    "roll_mean_7",
    "roll_mean_28",
    "roll_std_7",
    "roll_std_28",
]

DIRECT_FEATURES = [
    *DIRECT_HISTORY_FEATURES,
    *EXOG_COLS,
    *CALENDAR_COLS,
]

Keeping feature names in explicit lists reduces alignment bugs. Training and prediction must use the same columns in the same order.

Generate a neutral retail-sales series

The generated process contains:

A gradual upward trend.
Weekly seasonality.
Annual seasonality.
Planned promotional campaigns.
Occasional extra promotions.
Price variation.
Holiday effects.
Random noise.

def make_retail_series() -> pd.DataFrame:
    rng = np.random.default_rng(42)

    dates = pd.date_range(
        "2021-01-01",
        "2026-07-31",
        freq="D",
    )

    n_rows = len(dates)
    time_index = np.arange(n_rows)

    day_of_week = dates.dayofweek.to_numpy()
    day_of_year = dates.dayofyear.to_numpy()

    scheduled_campaign = (
        (
            (dates.day >= 3)
            & (dates.day <= 6)
            & dates.month.isin([2, 5, 8, 11])
        )
        | (
            (dates.month == 12)
            & (dates.day >= 15)
        )
    ).astype(int)

    occasional_promo = rng.binomial(
        1,
        0.05,
        n_rows,
    )

    promo = np.maximum(
        scheduled_campaign,
        occasional_promo,
    ).astype(int)

    holiday = (
        (
            (dates.month == 12)
            & (dates.day >= 24)
        )
        | (
            (dates.month == 1)
            & (dates.day <= 2)
        )
        | (
            (dates.month == 11)
            & (dates.day >= 25)
            & (dates.day <= 30)
        )
    ).astype(int)

    price = (
        10.5
        + 0.25
        * np.sin(
            2
            * np.pi
            * day_of_year
            / 365.25
        )
        + 0.0007 * time_index
        - 0.8 * promo
        + rng.normal(
            0,
            0.08,
            n_rows,
        )
    )

    weekly_pattern = (
        7.5
        * np.sin(
            2
            * np.pi
            * day_of_week
            / 7
        )
        + 4.0
        * (day_of_week >= 5)
    )

    annual_pattern = (
        6.0
        * np.sin(
            2
            * np.pi
            * day_of_year
            / 365.25
        )
        + 3.0
        * np.cos(
            4
            * np.pi
            * day_of_year
            / 365.25
        )
    )

    trend = 0.018 * time_index

    noise = rng.normal(
        0,
        3.2,
        n_rows,
    )

    sales = (
        75
        + trend
        + weekly_pattern
        + annual_pattern
        + 17 * promo
        + 9 * holiday
        - 3.8 * (price - 10.5)
        + noise
    )

    sales = np.maximum(
        0,
        np.round(sales, 1),
    )

    return pd.DataFrame(
        {
            "date": dates,
            "sales": sales,
            "promo": promo,
            "price": price.round(3),
            "holiday": holiday,
        }
    )


data = make_retail_series()

print(data.head().to_string(index=False))
print(data.tail().to_string(index=False))

The generated future price is treated as a committed plan. That assumption is valid for this exercise, but it must be audited in a real deployment.

Add calendar and Fourier seasonality features

A raw day-of-week integer is useful for tree splits, but it does not express cyclic distance.

Sunday, encoded as 6, is chronologically adjacent to Monday, encoded as 0. Numerically, however, those values look far apart.

Fourier features map a repeating cycle onto sine and cosine coordinates.

For annual seasonality:

year_sin_1 and year_cos_1 represent one cycle per year
year_sin_2 and year_cos_2 represent two cycles per year

Higher harmonics allow more complicated seasonal shapes, but they also increase the model’s capacity to fit noise. Choose the number of harmonics through temporal backtesting.

def add_calendar_features(
    frame: pd.DataFrame,
) -> pd.DataFrame:
    out = frame.copy()

    dates = pd.DatetimeIndex(out["date"])

    day_of_year = dates.dayofyear.to_numpy()

    out["dow"] = (
        dates.dayofweek
        .to_numpy(dtype=np.int8)
    )

    out["month"] = (
        dates.month
        .to_numpy(dtype=np.int8)
    )

    out["is_weekend"] = (
        dates.dayofweek
        .to_numpy()
        >= 5
    ).astype(np.int8)

    for harmonic in (1, 2):
        angle = (
            2
            * np.pi
            * harmonic
            * day_of_year
            / 365.25
        )

        out[
            f"year_sin_{harmonic}"
        ] = np.sin(angle)

        out[
            f"year_cos_{harmonic}"
        ] = np.cos(angle)

    return out

Calendar features describe when an observation occurs. Lagged features describe the series’ recent state.

You usually need both.

Create leakage-safe recursive features

The recursive design uses:

Lags 1, 7, 14, and 28.
Seven- and 28-day rolling means.
Seven- and 28-day rolling standard deviations.
Known future exogenous variables.
Calendar and Fourier features.

The rolling calculations begin with sales.shift(1).

def make_recursive_features(
    frame: pd.DataFrame,
) -> pd.DataFrame:
    out = add_calendar_features(frame)

    for lag in LAGS:
        out[f"lag_{lag}"] = (
            out["sales"]
            .shift(lag)
        )

    past_only = out["sales"].shift(1)

    for window in ROLL_WINDOWS:
        out[f"roll_mean_{window}"] = (
            past_only
            .rolling(window)
            .mean()
        )

        out[f"roll_std_{window}"] = (
            past_only
            .rolling(window)
            .std(ddof=0)
        )

    return out

A useful audit question for every training row is:

Could this value have been calculated immediately before the target occurred?

When the answer is no, shift the feature, replace it with an archived forecast, or remove it.

Configure the XGBoost regressor

The following settings are conservative tutorial defaults, not universal optimum values.

def new_model(
    **overrides: Any,
) -> XGBRegressor:
    params: dict[str, Any] = {
        "n_estimators": 180,
        "learning_rate": 0.05,
        "max_depth": 5,
        "min_child_weight": 4,
        "subsample": 0.9,
        "colsample_bytree": 0.9,
        "reg_lambda": 4.0,
        "objective": "reg:squarederror",
        "tree_method": "hist",
        "random_state": 42,
        "n_jobs": 4,
    }

    params.update(overrides)

    return XGBRegressor(**params)

Important interactions include:

Smaller learning_rate generally requires more trees.
Deeper trees capture more interactions but can fit temporal noise.
Larger min_child_weight discourages highly specific leaves.
Row and column subsampling provide regularization.
reg_lambda shrinks leaf weights.
Lag selection can matter more than a fine-grained parameter sweep.

Do not tune these parameters with shuffled folds.

A practical tuning design uses either:

An inner trailing validation window inside each training fold.
A set of historical forecast origins reserved for tuning.
Nested temporal backtesting when the compute budget permits it.

Train and run a recursive forecast

Training is ordinary supervised regression after feature construction.

Prediction is sequential. Every new estimate is appended before the next feature row is calculated.

The clip_lower parameter is important. Raw sales forecasts can be clipped at zero, but forecasts of a differenced target must be allowed to become negative.

def fit_recursive(
    train: pd.DataFrame,
) -> XGBRegressor:
    supervised = (
        make_recursive_features(train)
        .dropna(
            subset=[
                *RECURSIVE_FEATURES,
                "sales",
            ]
        )
    )

    model = new_model()

    model.fit(
        supervised[RECURSIVE_FEATURES],
        supervised["sales"],
    )

    return model


def recursive_forecast(
    model: XGBRegressor,
    train: pd.DataFrame,
    future: pd.DataFrame,
    *,
    clip_lower: float | None = 0.0,
) -> np.ndarray:
    history = train[
        [
            "date",
            "sales",
            *EXOG_COLS,
        ]
    ].copy()

    predictions: list[float] = []

    future_columns = [
        "date",
        *EXOG_COLS,
    ]

    for row in future[
        future_columns
    ].itertuples(index=False):
        next_row = {
            "date": row.date,
            "sales": np.nan,
            "promo": row.promo,
            "price": row.price,
            "holiday": row.holiday,
        }

        history = pd.concat(
            [
                history,
                pd.DataFrame([next_row]),
            ],
            ignore_index=True,
        )

        feature_row = (
            make_recursive_features(
                history
            )
            .iloc[[-1]]
        )

        prediction = float(
            model.predict(
                feature_row[
                    RECURSIVE_FEATURES
                ]
            )[0]
        )

        if clip_lower is not None:
            prediction = max(
                clip_lower,
                prediction,
            )

        history.loc[
            history.index[-1],
            "sales",
        ] = prediction

        predictions.append(prediction)

    return np.asarray(
        predictions,
        dtype=float,
    )

The earlier draft always applied max(0.0, prediction) inside this function. That is correct for sales levels but wrong when the same function predicts weekly changes.

A weekly difference can legitimately be negative:

sales this Tuesday − sales last Tuesday = −12

Clipping that change to zero would systematically prevent the differenced model from forecasting declines.

Construct the direct strategy correctly

For horizon h, a direct training row at origin t contains:

Sales observed through t.
Rolling features ending at t.
Calendar fields for t+h.
Known exogenous fields for t+h.
Target sales at t+h.

The origin’s current sales is valid because the forecast is issued after time t has been observed.

def make_direct_origin_features(
    frame: pd.DataFrame,
) -> pd.DataFrame:
    out = frame.copy()

    out["lag_0"] = out["sales"]

    for lag in (1, 7, 14, 28):
        out[f"lag_{lag}"] = (
            out["sales"]
            .shift(lag)
        )

    for window in ROLL_WINDOWS:
        out[f"roll_mean_{window}"] = (
            out["sales"]
            .rolling(window)
            .mean()
        )

        out[f"roll_std_{window}"] = (
            out["sales"]
            .rolling(window)
            .std(ddof=0)
        )

    return out


def make_direct_design(
    train: pd.DataFrame,
    step: int,
) -> pd.DataFrame:
    if step < 1:
        raise ValueError(
            "step must be at least 1"
        )

    origin = make_direct_origin_features(
        train
    )

    future_features = (
        add_calendar_features(train)[
            [
                *EXOG_COLS,
                *CALENDAR_COLS,
            ]
        ]
    )

    design = origin[
        DIRECT_HISTORY_FEATURES
    ].copy()

    for column in [
        *EXOG_COLS,
        *CALENDAR_COLS,
    ]:
        design[column] = (
            future_features[column]
            .shift(-step)
        )

    design["target"] = (
        train["sales"]
        .shift(-step)
    )

    return design.dropna(
        subset=[
            *DIRECT_FEATURES,
            "target",
        ]
    )


def fit_direct_models(
    train: pd.DataFrame,
    horizon: int,
) -> list[XGBRegressor]:
    models: list[XGBRegressor] = []

    for step in range(
        1,
        horizon + 1,
    ):
        design = make_direct_design(
            train,
            step,
        )

        model = new_model()

        model.fit(
            design[DIRECT_FEATURES],
            design["target"],
        )

        models.append(model)

    return models


def direct_forecast(
    models: list[XGBRegressor],
    train: pd.DataFrame,
    future: pd.DataFrame,
) -> np.ndarray:
    if len(models) != len(future):
        raise ValueError(
            "The number of direct models "
            "must equal the forecast horizon."
        )

    origin_row = (
        make_direct_origin_features(
            train
        )
        .iloc[-1]
    )

    future_features = (
        add_calendar_features(future)
    )

    predictions: list[float] = []

    for step, model in enumerate(
        models,
        start=1,
    ):
        row = {
            name: float(
                origin_row[name]
            )
            for name
            in DIRECT_HISTORY_FEATURES
        }

        future_row = (
            future_features
            .iloc[step - 1]
        )

        for column in [
            *EXOG_COLS,
            *CALENDAR_COLS,
        ]:
            row[column] = float(
                future_row[column]
            )

        design_row = pd.DataFrame(
            [row]
        )[DIRECT_FEATURES]

        prediction = float(
            model.predict(
                design_row
            )[0]
        )

        predictions.append(
            max(0.0, prediction)
        )

    return np.asarray(
        predictions,
        dtype=float,
    )

A frequent direct-forecasting bug is to construct the row at the target date instead of at the forecast origin.

That can accidentally introduce values such as y_t+h−1 even though they would not be known when the original h-step forecast was issued.

Always build the seasonal-naive baseline first

For daily retail data with weekly seasonality, a serious baseline predicts each future day from the corresponding day in the latest observed week.

next Monday = most recently observed Monday
next Tuesday = most recently observed Tuesday

def seasonal_naive(
    train: pd.DataFrame,
    horizon: int,
    season: int = 7,
) -> np.ndarray:
    if len(train) < season:
        raise ValueError(
            "Training data is shorter "
            "than the seasonal period."
        )

    last_season = (
        train["sales"]
        .iloc[-season:]
        .to_numpy(dtype=float)
    )

    return np.resize(
        last_season,
        horizon,
    )

This model:

Has no learned parameters.
Trains instantly.
Is easy to explain.
Often performs surprisingly well.
Exposes whether the complex model adds genuine predictive value.

Never declare victory because a forecast looks smooth.

Report the learned model and seasonal-naive baseline under the same:

Cutoffs.
Horizon.
Metrics.
Target availability assumptions.
Exogenous scenarios.

A useful skill score based on MAE is:

MAE skill versus naive = 1 − model MAE / seasonal-naive MAE

Interpretation:

Positive skill means the model beats seasonal naive.
Zero means it ties.
Negative skill means the model is worse.

Why trees struggle to extrapolate trend

A decision tree divides feature space into regions and assigns a learned value to each terminal leaf.

Suppose it splits an integer time feature at:

time_index < 500
time_index < 1,000
time_index < 1,500

A future value such as time_index = 2,000 follows the same final branch as all sufficiently large historical values.

The leaf returns a learned constant. It does not extend a slope.

Adding a raw time index can help a tree distinguish historical regimes, but it does not transform piecewise-constant leaves into a linear extrapolator.

Practical workarounds include:

Fit and remove a trend, model residuals, then add the future trend back.
Difference the target.
Predict change relative to a seasonal lag.
Supply an externally extrapolated trend as a feature.
Combine a trend model with boosted-tree residual correction.

Every fitted trend or transformation must be estimated inside the training fold.

Seasonal differencing as an extrapolation workaround

For weekly seasonal differencing:

difference_t = sales_t − sales_t−7

The model predicts change relative to the same weekday one week earlier.

Reconstruction uses:

forecast_t = predicted_difference_t + sales_t−7

For horizons beyond seven days, the seasonal reference eventually becomes an earlier reconstructed forecast. Error propagation therefore still exists.

def fit_seasonally_differenced_recursive(
    train: pd.DataFrame,
    season: int = 7,
) -> tuple[
    XGBRegressor,
    pd.DataFrame,
]:
    if season < 1:
        raise ValueError(
            "season must be positive"
        )

    differenced_train = train.copy()

    differenced_train["sales"] = (
        train["sales"]
        .diff(season)
    )

    model = fit_recursive(
        differenced_train
    )

    return model, differenced_train


def seasonally_differenced_forecast(
    model: XGBRegressor,
    differenced_train: pd.DataFrame,
    raw_train: pd.DataFrame,
    future: pd.DataFrame,
    season: int = 7,
) -> np.ndarray:
    difference_predictions = (
        recursive_forecast(
            model,
            differenced_train,
            future,
            clip_lower=None,
        )
    )

    raw_history = (
        raw_train["sales"]
        .astype(float)
        .tolist()
    )

    forecasts: list[float] = []

    for predicted_change in (
        difference_predictions
    ):
        seasonal_reference = (
            raw_history[-season]
        )

        prediction = (
            seasonal_reference
            + float(predicted_change)
        )

        prediction = max(
            0.0,
            prediction,
        )

        raw_history.append(
            prediction
        )

        forecasts.append(
            prediction
        )

    return np.asarray(
        forecasts,
        dtype=float,
    )

Notice where clipping occurs:

predicted weekly change: not clipped
reconstructed sales level: clipped at zero

That distinction fixes the correctness problem in the earlier draft.

Differencing is not automatically beneficial. It can:

Reduce trend extrapolation risk.
Remove useful level information.
Amplify short-term noise.
Complicate reconstruction.
Introduce its own recursive dependency.

Treat it as one candidate and let the backtest decide.

Implement expanding-window backtesting

The following experiment uses four non-overlapping 28-day test windows.

For every fold, it:

Trains only on observations before the cutoff.
Fits the raw-target recursive model.
Fits 28 direct models.
Fits the seasonally differenced recursive model.
Forecasts the same 28-day period.
Scores all methods, including seasonal naive.
Expands the training history for the next fold.

def score_forecast(
    y_true: np.ndarray,
    y_pred: np.ndarray,
) -> dict[str, float]:
    if y_true.shape != y_pred.shape:
        raise ValueError(
            "Actual and predicted arrays "
            "must have equal shapes."
        )

    mae = mean_absolute_error(
        y_true,
        y_pred,
    )

    rmse = mean_squared_error(
        y_true,
        y_pred,
    ) ** 0.5

    return {
        "mae": float(mae),
        "rmse": float(rmse),
    }


def expanding_window_backtest(
    data: pd.DataFrame,
    horizon: int = 28,
    folds: int = 4,
) -> pd.DataFrame:
    if horizon < 1:
        raise ValueError(
            "horizon must be positive"
        )

    if folds < 1:
        raise ValueError(
            "folds must be positive"
        )

    required_test_rows = (
        horizon * folds
    )

    if len(data) <= required_test_rows:
        raise ValueError(
            "Not enough data for the "
            "requested backtest."
        )

    first_cutoff = (
        len(data)
        - required_test_rows
    )

    records: list[
        dict[str, object]
    ] = []

    cutoffs = range(
        first_cutoff,
        len(data),
        horizon,
    )

    for fold, cutoff in enumerate(
        cutoffs,
        start=1,
    ):
        train = (
            data.iloc[:cutoff]
            .copy()
        )

        test = (
            data.iloc[
                cutoff:
                cutoff + horizon
            ]
            .copy()
        )

        if len(test) != horizon:
            raise ValueError(
                "Every test fold must "
                f"contain {horizon} rows."
            )

        y_true = (
            test["sales"]
            .to_numpy(dtype=float)
        )

        recursive_model = (
            fit_recursive(train)
        )

        direct_models = (
            fit_direct_models(
                train,
                horizon,
            )
        )

        (
            diff_model,
            differenced_train,
        ) = (
            fit_seasonally_differenced_recursive(
                train
            )
        )

        forecasts = {
            "seasonal_naive": (
                seasonal_naive(
                    train,
                    horizon,
                )
            ),
            "recursive": (
                recursive_forecast(
                    recursive_model,
                    train,
                    test,
                )
            ),
            "direct": (
                direct_forecast(
                    direct_models,
                    train,
                    test,
                )
            ),
            "recursive_diff7": (
                seasonally_differenced_forecast(
                    diff_model,
                    differenced_train,
                    train,
                    test,
                )
            ),
        }

        for method, prediction in (
            forecasts.items()
        ):
            metrics = score_forecast(
                y_true,
                prediction,
            )

            records.append(
                {
                    "fold": fold,
                    "cutoff": (
                        train["date"]
                        .iloc[-1]
                    ),
                    "test_start": (
                        test["date"]
                        .iloc[0]
                    ),
                    "test_end": (
                        test["date"]
                        .iloc[-1]
                    ),
                    "method": method,
                    **metrics,
                }
            )

    return pd.DataFrame(
        records
    )

Run the experiment and calculate skill against seasonal naive:

results = expanding_window_backtest(
    data,
    horizon=28,
    folds=4,
)

summary = (
    results
    .groupby(
        "method",
        as_index=False,
    )[
        [
            "mae",
            "rmse",
        ]
    ]
    .mean()
    .sort_values("mae")
)

baseline_mae = float(
    summary.loc[
        summary["method"]
        == "seasonal_naive",
        "mae",
    ]
    .iloc[0]
)

summary[
    "mae_skill_vs_naive"
] = (
    1
    - summary["mae"]
    / baseline_mae
)

print(
    summary
    .round(3)
    .to_string(index=False)
)

The exact decimals can vary slightly across processor architectures, thread scheduling, and package builds. Do not separate the scores from the environment that produced them.

Record both in a machine-readable report:

def records_with_dates_as_text(
    frame: pd.DataFrame,
) -> list[dict[str, object]]:
    serializable = frame.copy()

    for column in [
        "cutoff",
        "test_start",
        "test_end",
    ]:
        if column in serializable:
            serializable[column] = (
                pd.to_datetime(
                    serializable[column]
                )
                .dt.strftime("%Y-%m-%d")
            )

    return serializable.to_dict(
        orient="records"
    )


report = {
    "environment": version_report(),
    "configuration": {
        "seed": 42,
        "horizon": 28,
        "folds": 4,
        "lags": list(LAGS),
        "rolling_windows": list(
            ROLL_WINDOWS
        ),
        "exogenous_columns": EXOG_COLS,
    },
    "fold_results": (
        records_with_dates_as_text(
            results
        )
    ),
    "summary": (
        summary
        .round(6)
        .to_dict(orient="records")
    ),
}

output_dir = Path("artifacts")
output_dir.mkdir(
    parents=True,
    exist_ok=True,
)

report_path = (
    output_dir
    / "forecast_backtest_report.json"
)

report_path.write_text(
    json.dumps(
        report,
        indent=2,
        sort_keys=True,
    ),
    encoding="utf-8",
)

print(
    f"Wrote {report_path}"
)

On the seeded synthetic process, a healthy run should normally show both raw-target tree strategies beating seasonal naive. Direct and recursive forecasting are often close, while seasonal differencing may help or hurt.

Do not treat that expected pattern as a universal benchmark. The important result is what happens under your own historical cutoffs.

How to interpret the comparison

Several outcomes are possible.

Direct beats recursive

That suggests avoiding estimated intermediate lags is valuable enough to justify separate horizon models.

It is especially plausible when:

The horizon is long.
Short-term errors strongly affect later lags.
Horizon-specific exogenous effects matter.
Sufficient training data exists for every direct model.

Recursive is close to direct

A single recursive model may offer a better production trade-off when:

Retraining is frequent.
There are thousands of series.
Forecast latency matters.
The requested horizon changes.
Model governance makes large bundles expensive.

A tiny accuracy gain from 28 models may not justify 28 times the maintenance surface.

Differencing wins

The raw target may contain a trend that tree leaves cannot extend reliably.

Differencing can make the transformed series more stationary and easier to model.

Confirm that the gain survives:

Several folds.
Several horizons.
Different seasons.
Promotion and non-promotion periods.

Differencing loses

That is also a valid result.

The transformation may have:

Removed useful level information.
Increased target noise.
Made exogenous relationships harder to learn.
Added reconstruction error.

A theoretically motivated transformation still needs empirical evidence.

Seasonal naive wins

Stop and investigate before tuning harder.

Possible explanations include:

Leakage was removed and the original result was unrealistic.
The series is already highly stable by weekday.
Promotions or prices are not genuinely known.
The model is overfitting old regimes.
The chosen horizon is too long for autoregressive features.
Training data is insufficient.
The feature set does not capture important seasonal structure.

The baseline is not ceremonial. It is a deployment gate.

Inspect error by forecast horizon

Average error can hide recursive degradation.

A model may perform well for the first week and deteriorate sharply afterward.

Collect every forecast step:

def collect_backtest_predictions(
    data: pd.DataFrame,
    horizon: int = 28,
    folds: int = 4,
) -> pd.DataFrame:
    first_cutoff = (
        len(data)
        - folds * horizon
    )

    records: list[
        dict[str, object]
    ] = []

    for fold, cutoff in enumerate(
        range(
            first_cutoff,
            len(data),
            horizon,
        ),
        start=1,
    ):
        train = (
            data.iloc[:cutoff]
            .copy()
        )

        test = (
            data.iloc[
                cutoff:
                cutoff + horizon
            ]
            .copy()
        )

        recursive_model = (
            fit_recursive(train)
        )

        direct_models = (
            fit_direct_models(
                train,
                horizon,
            )
        )

        forecasts = {
            "seasonal_naive": (
                seasonal_naive(
                    train,
                    horizon,
                )
            ),
            "recursive": (
                recursive_forecast(
                    recursive_model,
                    train,
                    test,
                )
            ),
            "direct": (
                direct_forecast(
                    direct_models,
                    train,
                    test,
                )
            ),
        }

        actual = (
            test["sales"]
            .to_numpy(dtype=float)
        )

        for method, prediction in (
            forecasts.items()
        ):
            for step, (
                date,
                observed,
                predicted,
            ) in enumerate(
                zip(
                    test["date"],
                    actual,
                    prediction,
                    strict=True,
                ),
                start=1,
            ):
                records.append(
                    {
                        "fold": fold,
                        "date": date,
                        "horizon": step,
                        "method": method,
                        "actual": observed,
                        "prediction": (
                            predicted
                        ),
                        "absolute_error": abs(
                            observed
                            - predicted
                        ),
                    }
                )

    return pd.DataFrame(
        records
    )


prediction_log = (
    collect_backtest_predictions(
        data,
        horizon=28,
        folds=4,
    )
)

horizon_errors = (
    prediction_log
    .groupby(
        [
            "method",
            "horizon",
        ],
        as_index=False,
    )["absolute_error"]
    .mean()
)

print(
    horizon_errors
    .head(12)
    .round(3)
    .to_string(index=False)
)

Plot the horizon profile:

plt.figure(
    figsize=(10, 5)
)

for method, group in (
    horizon_errors
    .groupby("method")
):
    plt.plot(
        group["horizon"],
        group["absolute_error"],
        marker="o",
        markersize=3,
        label=method,
    )

plt.xlabel(
    "Forecast horizon in days"
)

plt.ylabel(
    "Mean absolute error"
)

plt.title(
    "Backtest error by forecast horizon"
)

plt.legend()
plt.tight_layout()
plt.show()

A recursive curve often rises because later rows contain more estimated lags.

Direct error can rise too. Farther targets are intrinsically harder, and current observations may have weaker relationships with distant outcomes.

Diagnose error compounding from the horizon profile rather than assuming it from the strategy name.

Use exogenous regressors without creating an oracle

The example passes promotion, price, and holiday fields for each target date.

That is valid only because the exercise assumes they are available when the forecast is issued.

Variable	Availability	Treatment
Day of week	Deterministic	Safe future feature
Public holiday	Known calendar	Safe future feature
Planned promotion	Usually scheduled	Safe only when plans are frozen
Planned price	Business-dependent	Safe only when committed
Weather forecast	Available but uncertain	Use archived forecasts
Realized weather	Unknown at origin	Leakage
Future stockout	Usually unknown	Model separately or use scenarios

When an exogenous variable is uncertain, produce conditional scenario forecasts.

For example:

Base case: The promotion runs as approved.
Low case: The promotion is canceled.
High case: The promotion receives extra advertising.

A forecast using a planned promotion is not unconditional demand. It is demand conditional on that plan.

Also audit data latency.

Suppose Sunday’s final sales arrive on Tuesday. A forecast issued Monday morning cannot use Sunday as lag_1, even though Sunday is chronologically in the past.

Your historical feature pipeline must reproduce actual reporting delays.

Quantile objectives for prediction intervals

A point forecast answers only one part of an operational question.

Inventory, staffing, and capacity planning also need uncertainty estimates:

How much stock covers a high-demand outcome?
How uncertain is the promotion period?
How wide should a safety margin be?
Is uncertainty increasing with the horizon?

XGBoost provides the reg:quantileerror objective and the quantile_alpha parameter. The alpha value may be a scalar or an ascending list of quantiles. Current documentation describes the objective as a smooth approximation to pinball loss.

For a nominal central 80% interval, estimate:

0.10 for the lower bound.
0.50 for the median.
0.90 for the upper bound.

Direct quantile models are convenient because each horizon is estimated from the same observed origin rather than from a recursively generated quantile path.

def fit_direct_quantile_models(
    train: pd.DataFrame,
    horizon: int,
    alphas: tuple[
        float,
        ...,
    ] = (
        0.10,
        0.50,
        0.90,
    ),
) -> list[XGBRegressor]:
    alpha_values = np.asarray(
        alphas,
        dtype=float,
    )

    if np.any(
        np.diff(alpha_values) <= 0
    ):
        raise ValueError(
            "alphas must be strictly "
            "increasing"
        )

    if np.any(
        (alpha_values <= 0)
        | (alpha_values >= 1)
    ):
        raise ValueError(
            "alphas must be between "
            "0 and 1"
        )

    models: list[
        XGBRegressor
    ] = []

    for step in range(
        1,
        horizon + 1,
    ):
        design = make_direct_design(
            train,
            step,
        )

        model = new_model(
            objective=(
                "reg:quantileerror"
            ),
            quantile_alpha=(
                alpha_values
            ),
            n_estimators=220,
            tree_method="hist",
        )

        model.fit(
            design[DIRECT_FEATURES],
            design["target"],
        )

        models.append(model)

    return models


def direct_quantile_forecast(
    models: list[XGBRegressor],
    train: pd.DataFrame,
    future: pd.DataFrame,
) -> np.ndarray:
    if len(models) != len(future):
        raise ValueError(
            "The number of models "
            "must equal the horizon."
        )

    origin_row = (
        make_direct_origin_features(
            train
        )
        .iloc[-1]
    )

    future_features = (
        add_calendar_features(future)
    )

    predictions: list[
        np.ndarray
    ] = []

    for step, model in enumerate(
        models,
        start=1,
    ):
        row = {
            name: float(
                origin_row[name]
            )
            for name
            in DIRECT_HISTORY_FEATURES
        }

        future_row = (
            future_features
            .iloc[step - 1]
        )

        for column in [
            *EXOG_COLS,
            *CALENDAR_COLS,
        ]:
            row[column] = float(
                future_row[column]
            )

        design_row = pd.DataFrame(
            [row]
        )[DIRECT_FEATURES]

        quantiles = np.asarray(
            model.predict(
                design_row
            ),
            dtype=float,
        ).reshape(-1)

        predictions.append(
            quantiles
        )

    return np.vstack(
        predictions
    )

Test the final 28 days:

interval_horizon = 28

interval_train = (
    data.iloc[
        :-interval_horizon
    ]
    .copy()
)

interval_test = (
    data.iloc[
        -interval_horizon:
    ]
    .copy()
)

quantile_models = (
    fit_direct_quantile_models(
        interval_train,
        horizon=interval_horizon,
        alphas=(
            0.10,
            0.50,
            0.90,
        ),
    )
)

quantile_predictions = (
    direct_quantile_forecast(
        quantile_models,
        interval_train,
        interval_test,
    )
)

lower = (
    quantile_predictions[:, 0]
)

median = (
    quantile_predictions[:, 1]
)

upper = (
    quantile_predictions[:, 2]
)

actual = (
    interval_test["sales"]
    .to_numpy(dtype=float)
)

crossings = int(
    np.sum(
        (lower > median)
        | (median > upper)
    )
)

coverage = float(
    np.mean(
        (actual >= lower)
        & (actual <= upper)
    )
)

mean_width = float(
    np.mean(
        upper - lower
    )
)

print(
    f"Quantile crossings: "
    f"{crossings}"
)

print(
    "Empirical interval coverage: "
    f"{coverage:.3f}"
)

print(
    "Mean interval width: "
    f"{mean_width:.3f}"
)

Evaluate each quantile with pinball loss:

def pinball_loss(
    y_true: np.ndarray,
    quantile_prediction: (
        np.ndarray
    ),
    alpha: float,
) -> float:
    error = (
        y_true
        - quantile_prediction
    )

    loss = np.maximum(
        alpha * error,
        (alpha - 1) * error,
    )

    return float(
        np.mean(loss)
    )


for alpha, prediction in zip(
    (
        0.10,
        0.50,
        0.90,
    ),
    quantile_predictions.T,
    strict=True,
):
    loss = pinball_loss(
        actual,
        prediction,
        alpha,
    )

    print(
        f"alpha={alpha:.2f}, "
        f"pinball_loss={loss:.3f}"
    )

Plot the interval:

plt.figure(
    figsize=(11, 5)
)

plt.plot(
    interval_test["date"],
    actual,
    marker="o",
    label="actual",
)

plt.plot(
    interval_test["date"],
    median,
    label="predicted median",
)

plt.fill_between(
    interval_test["date"],
    lower,
    upper,
    alpha=0.25,
    label=(
        "10th–90th quantile "
        "interval"
    ),
)

plt.xlabel("Date")
plt.ylabel("Sales")

plt.title(
    "Direct quantile forecast"
)

plt.legend()
plt.tight_layout()
plt.show()

The plotting block is executable Python and is correctly fenced as python; it is not plain-text output.

Do not call the interval calibrated merely because it uses the 10th and 90th quantiles.

Backtest:

Empirical coverage.
Mean interval width.
Pinball loss for every quantile.
Coverage by horizon.
Coverage during promotions.
Coverage during holidays.
Quantile-crossing frequency.

Quantile crossing occurs when estimates are not ordered:

predicted 10th percentile > predicted median

or:

predicted median > predicted 90th percentile

Silently sorting crossed values makes a chart prettier, but it changes the model output without fixing calibration.

Investigate:

Training sample size.
Model capacity.
Horizon-specific instability.
Alpha spacing.
Feature leakage.
Post-hoc calibration.
Separate models for different quantiles.

Also distinguish marginal intervals from a simultaneous path interval.

An 80% interval at each individual day does not imply an 80% probability that all 28 observations will remain inside their respective bands.

Why random cross-validation can look absurdly good

Consider a slowly changing series with lag_1.

A shuffled split may place Tuesday in training and Wednesday in validation. Wednesday’s row contains Tuesday’s target as lag_1.

Using Tuesday’s target is legitimate at deployment, so this is not direct leakage. The problem is that the validation row may be almost identical to nearby training contexts.

Random splitting becomes especially misleading with:

Long rolling windows.
Strong autocorrelation.
Slowly changing prices.
Multi-day promotions.
Encoded store or product identities.
Global models trained across related series.
Features generated before splitting.
Regime changes hidden by averaging.

Temporal backtesting also reveals instability across eras.

A model may work before a pricing-policy change and fail afterward. One shuffled average conceals that.

Feature-engineering ideas worth testing

The tutorial deliberately uses a compact feature set.

Real retail systems may benefit from additional past-only predictors.

More autoregressive features

Candidates include:

Lags 2 through 6.
Lag 364 for approximately annual daily seasonality.
Rolling minimum and maximum.
Rolling median.
Exponentially weighted means.
Short-minus-long rolling averages.
Same-weekday rolling statistics.
Days since the last nonzero sale.
Recent promotion-adjusted sales.

Do not blindly create hundreds of correlated lags.

Extra features increase:

Memory use.
Training time.
Tuning cost.
Opportunities to overfit.
The difficulty of leakage audits.

Promotion structure

A single binary flag may be too crude.

Useful planned features can include:

Promotion type.
Discount percentage.
Days since campaign start.
Days until campaign end.
Advertising channel.
Promotion depth relative to normal price.
Whether related products are promoted.

These fields remain valid only when they are known at the forecast origin.

Relative calendar events

Demand often responds to proximity rather than a single-day flag.

Candidates include:

Days before a holiday.
Days after a holiday.
Payday proximity.
School vacation periods.
Month-end.
Quarter-end.
Event lead and lag indicators.

Cross-series information

For many products or locations, a global model can pool observations.

Possible static or identifying features include:

Store.
Product.
Category.
Department.
Region.
Package size.
Normal price band.

The split must still occur by time.

Randomly holding out rows from a panel creates temporal leakage at a larger scale.

Scaling beyond manual feature code

The manual implementation is valuable because every alignment decision is visible.

In production, the feature-generation boilerplate becomes repetitive.

MLForecast

The current MLForecast API uses recursive forecasting by default. Direct one-model-per-horizon training is available through max_horizon, while horizons can train only selected forecast steps. Its API also supports lag transformations, cross-validation, exogenous variables, and prediction-interval workflows.

That is useful when you need:

Many related series.
Efficient lag generation.
Reusable cross-validation windows.
Global models.
Selected direct horizons.
Distributed or alternative dataframe backends.

Skforecast

Current Skforecast forecasters wrap estimators compatible with the scikit-learn API and automate recursive or direct forecasting, lag construction, exogenous variables, window features, and backtesting.

Its differentiation parameter means the order of ordinary consecutive differencing:

first-order difference_t = y_t − y_t−1

The transformation is reversed when forecasts are returned.

It should not be described as automatically reproducing this tutorial’s lag-7 seasonal difference:

seasonal difference_t = y_t − y_t−7

Those are different transformations. Use a separately supported seasonal transformer or implement the lag-7 reconstruction explicitly, as shown earlier.

Forecasting libraries remove boilerplate, not responsibility.

You still need to define:

The forecast origin.
Target latency.
Future-feature availability.
Horizon.
Backtest windows.
Baseline.
Metrics.
Retraining policy.

Cherry on the cake: the M5 result that surprised the field

The M5 Accuracy competition asked participants to forecast 28 days of hierarchical Walmart unit sales.

Its winning submission was not a giant end-to-end deep neural network.

The winner used an equal-weighted ensemble of LightGBM models trained at multiple retail aggregation levels. It included recursive and non-recursive variants, built 220 models in total, and averaged six contributing models for each series. The models used sales history, identifiers, calendar information, events, promotions, and prices.

The third-place system was an equal-weighted ensemble of 43 deep neural networks containing multiple LSTM layers.

The lesson is not that deep learning is ineffective.

The lesson is that large-scale retail forecasting rewards system design:

Careful supervised reformulation.
Strong lag and calendar features.
Cross-learning across related series.
Historical forecast cutoffs.
Appropriate objectives.
Recursive and non-recursive variants.
Ensembling.
Comparison against strong baselines.

Gradient-boosted trees did not win because a tree inherently understands time. They won because the surrounding forecasting system presented time correctly.

A production checklist

Before shipping a boosted-tree forecast, verify each item.

Forecast definition

Record the exact issue time.
Record the operational horizon.
Document target publication latency.
State whether forecasts are conditional on future plans.

Feature causality

Shift rolling target features.
Audit every exogenous field.
Use archived forecasts for uncertain future variables.
Reproduce reporting latency.
Fit transformations only inside each fold.

Evaluation

Use expanding or rolling temporal windows.
Report seasonal naive beside every model.
Inspect fold-level scores, not only averages.
Inspect error by horizon.
Evaluate important business segments.
Compare recursive and direct strategies empirically.

Trend handling

Test the raw target.
Test ordinary or seasonal differencing.
Test detrending.
Test a trend-plus-residual hybrid.
Do not assume a time-index feature solves extrapolation.

Uncertainty

Evaluate pinball loss.
Measure empirical interval coverage.
Measure interval width.
Check quantile crossing.
Evaluate coverage by horizon and event type.
Distinguish pointwise intervals from path uncertainty.

Reproducibility

Save package versions.
Save cutoff dates.
Save feature definitions.
Save model parameters.
Save random seeds.
Save fold-level predictions.
Save the baseline forecasts.
Save the exact exogenous scenario.

Monitoring

Track live error against the trained model.
Track live error against seasonal naive.
Detect changes in target latency.
Detect missing future regressors.
Detect drift in prices, promotions, and demand.
Re-run historical backtests after major pipeline changes.

Take the next step

Run the complete versioned backtest before changing a single hyperparameter.

Then replace the generated series with one non-sensitive business series while preserving the same contract:

date
sales
known future promotion
known future price
known holiday calendar

Document exactly when each predictor becomes available. Refuse to tune the model until the seasonal-naive score is visible beside it.

After the point forecast beats that baseline across several historical cutoffs, inspect horizon-level degradation, test direct against recursive forecasting, and add backtested quantile intervals.

That sequence turns an attractive XGBoost demonstration into an evidence-based forecasting system.