XGBoost for Time Series Forecasting

XGBoost is not a time series model — it has no built-in understanding of temporal order, seasonality, or autocorrelation. But it becomes a powerful time series forecaster once you add that structure through feature engineering. The approach is to transform a forecasting problem into a supervised regression problem: create lag features, rolling statistics, and calendar features from your time series, then train XGBoost to predict the target from those engineered features. Done correctly, this approach is competitive with or better than classical methods (ARIMA, exponential smoothing) and traditional ML baselines on many real-world datasets.

The Core Idea: Supervised Regression on Lag Features

A time series y_t can be predicted from its own past values: y_t = f(y_{t-1}, y_{t-2}, ..., y_{t-k}). These past values are lag features. Add rolling means, rolling standard deviations, and calendar features (hour, day of week, month, is_weekend) and you have a rich feature set that captures trend, seasonality, and local volatility — everything XGBoost needs to learn patterns from.

import pandas as pd
import numpy as np

def create_time_features(df, target_col, lags, rolling_windows):
    """Create lag and rolling features from a time series DataFrame."""
    df = df.copy()

    # Lag features
    for lag in lags:
        df[f'lag_{lag}'] = df[target_col].shift(lag)

    # Rolling statistics
    for window in rolling_windows:
        df[f'rolling_mean_{window}'] = df[target_col].shift(1).rolling(window).mean()
        df[f'rolling_std_{window}']  = df[target_col].shift(1).rolling(window).std()
        df[f'rolling_min_{window}']  = df[target_col].shift(1).rolling(window).min()
        df[f'rolling_max_{window}']  = df[target_col].shift(1).rolling(window).max()

    # Calendar features (assumes DatetimeIndex)
    df['hour']        = df.index.hour
    df['dayofweek']   = df.index.dayofweek
    df['month']       = df.index.month
    df['quarter']     = df.index.quarter
    df['is_weekend']  = (df.index.dayofweek >= 5).astype(int)

    df.dropna(inplace=True)
    return df

# Example usage
df = pd.read_csv('sales.csv', index_col='date', parse_dates=True)
df = create_time_features(df, target_col='sales',
                          lags=[1, 2, 3, 7, 14, 28],
                          rolling_windows=[7, 14, 28])

Train/Test Split: Never Shuffle Time Series Data

The most critical rule for time series with XGBoost: split chronologically, never randomly. Random splitting leaks future information into the training set — the model sees future values through lag features and produces impossibly good validation scores that do not hold in production. Always use a cutoff date: everything before it trains the model, everything after tests it.

cutoff = '2024-01-01'
feature_cols = 

train = df[df.index < cutoff]
test  = df[df.index >= cutoff]

X_train, y_train = train[feature_cols], train['sales']
X_test,  y_test  = test[feature_cols],  test['sales']

Training and Evaluating the Model

import xgboost as xgb
from sklearn.metrics import mean_absolute_error, mean_squared_error

model = xgb.XGBRegressor(
    n_estimators=1000,
    learning_rate=0.05,
    max_depth=5,
    subsample=0.8,
    colsample_bytree=0.8,
    early_stopping_rounds=50,
    eval_metric='mae',
    random_state=42,
    verbosity=0
)
model.fit(X_train, y_train, eval_set=[(X_test, y_test)], verbose=200)

preds = model.predict(X_test)
mae  = mean_absolute_error(y_test, preds)
rmse = np.sqrt(mean_squared_error(y_test, preds))
mape = np.mean(np.abs((y_test - preds) / y_test)) * 100
print(f"MAE: {mae:.2f}  RMSE: {rmse:.2f}  MAPE: {mape:.2f}%")

Time Series Cross-Validation

Standard k-fold cross-validation shuffles data and must not be used with time series. Use TimeSeriesSplit from scikit-learn, which always trains on past data and validates on future data, simulating production deployment across multiple periods.

from sklearn.model_selection import TimeSeriesSplit, cross_val_score

tscv = TimeSeriesSplit(n_splits=5, gap=0)

model_cv = xgb.XGBRegressor(
    n_estimators=500, learning_rate=0.05, max_depth=5,
    random_state=42, verbosity=0
)

scores = cross_val_score(model_cv, df[feature_cols], df['sales'],
                         cv=tscv, scoring='neg_mean_absolute_error')
print(f"CV MAE: {-scores.mean():.2f} (+/- {scores.std():.2f})")

Figure 1 — TimeSeriesSplit: each fold trains on past, validates on future

Multi-Step Forecasting

XGBoost predicts one step at a time by default. For multi-step forecasts (predict 7 days ahead), use one of two strategies. The direct strategy trains a separate model for each horizon — one model for h=1, another for h=2, etc. The recursive strategy uses a single model and feeds predictions back as lag features for subsequent steps. Direct is more accurate but requires training N models; recursive is simpler but accumulates errors.

# Direct strategy — separate model per horizon
horizons = [1, 3, 7, 14]
models = {}

for h in horizons:
    df_h = create_time_features(df.copy(), 'sales',
                                lags=[h, h+1, h+6, h+13],
                                rolling_windows=[7, 28])
    X, y = df_h[feature_cols], df_h['sales'].shift(-h)
    df_h = df_h[:-h]  # drop last h rows (no target)
    X, y = df_h[feature_cols], df_h['sales']

    model_h = xgb.XGBRegressor(n_estimators=500, verbosity=0)
    model_h.fit(X, y)
    models[h] = model_h
    print(f"Trained h={h} model")

Feature Importance for Time Series

import shap
import matplotlib.pyplot as plt

explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X_test)

# Which features drive the forecast?
shap.summary_plot(shap_values, X_test, plot_type='bar', max_display=15)

# Built-in importance
fi = pd.Series(model.feature_importances_, index=feature_cols).sort_values(ascending=False)
print(fi.head(10))

In most time series applications, lag features dominate SHAP importance — recent lags (lag_1, lag_2) usually outrank all other features. If rolling statistics outrank recent lags, your series likely has strong seasonal patterns that the short-term lags miss. If calendar features (dayofweek, month) rank highly, the series has strong calendar-driven seasonality worth modelling explicitly with interaction terms between calendar features and lags.

Handling Multiple Time Series

When you have many related time series — sales across 100 stores, energy usage across 50 buildings — train a single global XGBoost model with a series identifier as a categorical feature. This often outperforms individual per-series models because XGBoost can learn patterns shared across series while using the identifier to distinguish series-specific behaviour.

# Stack all series into one DataFrame
dfs = []
for series_id in df_all['store_id'].unique():
    df_s = df_all[df_all['store_id']==series_id].copy()
    df_s = create_time_features(df_s, 'sales', lags=[1,7,28], rolling_windows=[7,28])
    df_s['series_id'] = series_id
    dfs.append(df_s)

df_global = pd.concat(dfs).sort_index()
# Train one XGBoost on df_global — series_id acts as a categorical identifier

XGBoost is a practical first choice for tabular time series forecasting — it requires no distributional assumptions, handles mixed feature types natively, generalises well across series types, and can be tuned with Optuna like any other XGBoost model. The key discipline is temporal integrity: chronological splits everywhere, lag features computed without peeking at future values, and validation on future periods that genuinely test generalisation rather than memorisation. Get those right and XGBoost is a strong baseline that dedicated time series methods often struggle to beat.

Figure 1 — TimeSeriesSplit: each fold trains on past, validates on future

When XGBoost Beats Classical Methods

XGBoost outperforms ARIMA and exponential smoothing most reliably in three situations. First, when external regressors matter — price, promotions, weather, economic indicators — XGBoost incorporates them naturally as features while ARIMA handles them awkwardly. Second, when the relationship between predictors and target is non-linear or involves interactions that linear models cannot capture. Third, when you have many related series and can train a global model that shares patterns across series. ARIMA and ETS retain advantages for short, clean, well-behaved univariate series with strong linear autocorrelation structure and no external regressors — classic inventory or financial series. For everything messier, XGBoost is usually worth trying first.

Common Mistakes

The most costly mistakes when using XGBoost for time series. Using random train/test splits — this leaks future information through lag features and produces validation scores that do not reflect production performance at all. Creating lag features before splitting — always create lags after splitting or use .shift() which respects the DataFrame index ordering. Using StandardScaler.fit() on the full dataset — fit the scaler on training data only and transform both train and test. Forgetting to retrain on the full dataset before deployment — models tuned on a train/val split should be retrained on all available data before serving predictions in production. And using the same lag depth for all horizons in multi-step forecasting — a model predicting 14 days ahead should not use lag_1 as a feature since that value will not be available 14 days in advance; minimum lag should be at least equal to the forecast horizon.

Feature Engineering Deep Dive

The quality of your lag and rolling features is the largest driver of XGBoost time series performance — more than any hyperparameter. A few principles that consistently help. Use multiple lag depths: short lags (1–3) capture recent momentum, medium lags (7, 14) capture weekly patterns, long lags (28, 90, 365) capture monthly and seasonal patterns. For hourly data, lags of 24, 48, 168 (one week in hours) are natural. Rolling statistics should use the shifted series (shift(1) before rolling) to avoid including the current value in the window. Beyond mean and std, consider rolling median (robust to spikes), rolling skewness (detects distributional shifts), and rolling first differences (captures acceleration/deceleration of trend).

For seasonal series, Fourier features often outperform simple calendar dummies. Instead of a binary is_monday feature, encode day-of-week as sine and cosine pairs — this lets the model learn smooth, continuous seasonality patterns rather than treating Monday and Sunday as completely unrelated despite being adjacent. The same applies to hour-of-day and month-of-year.

import numpy as np

def add_fourier_features(df, period, n_terms):
    """Add Fourier terms for a given seasonal period."""
    t = np.arange(len(df))
    for k in range(1, n_terms + 1):
        df[f'sin_{period}_{k}'] = np.sin(2 * np.pi * k * t / period)
        df[f'cos_{period}_{k}'] = np.cos(2 * np.pi * k * t / period)
    return df

# Weekly seasonality (7 days), 3 Fourier pairs
df = add_fourier_features(df, period=7, n_terms=3)
# Annual seasonality (365 days), 5 Fourier pairs
df = add_fourier_features(df, period=365, n_terms=5)

Hyperparameter Tuning for Time Series

Tuning XGBoost for time series uses the same Optuna pattern as tabular data, but with TimeSeriesSplit for cross-validation. One additional consideration: max_depth for time series should typically be shallower (3–5) than for tabular data because time series patterns are often smoother and less categorical. Deep trees can overfit noise in recent lags while missing the longer-term seasonal signal.

import optuna
from sklearn.model_selection import TimeSeriesSplit, cross_val_score

tscv = TimeSeriesSplit(n_splits=5)

def objective(trial):
    params = {
        'n_estimators': trial.suggest_int('n_estimators', 100, 1000),
        'learning_rate': trial.suggest_float('lr', 0.01, 0.3, log=True),
        'max_depth': trial.suggest_int('max_depth', 2, 6),
        'subsample': trial.suggest_float('subsample', 0.5, 1.0),
        'colsample_bytree': trial.suggest_float('colsample_bytree', 0.5, 1.0),
        'reg_lambda': trial.suggest_float('reg_lambda', 1e-2, 10, log=True),
        'random_state': 42, 'verbosity': 0
    }
    model = xgb.XGBRegressor(**params)
    scores = cross_val_score(model, df[feature_cols], df['sales'],
                             cv=tscv, scoring='neg_mean_absolute_error')
    return -scores.mean()

study = optuna.create_study(direction='minimize')
study.optimize(objective, n_trials=50)
print(study.best_params)

Residual Analysis and Model Validation

After training, plot the residuals (actual minus predicted) over time. Well-behaved residuals should look like white noise — no trend, no seasonality, no autocorrelation. If residuals show clear patterns, the model has systematic blind spots that more features or a different architecture could address. A systematic under-prediction at Christmas and over-prediction in January suggests the model needs stronger seasonal features for year-end effects. Gradually widening residual variance suggests the model handles high-volume periods poorly — consider log-transforming the target before training and exponentiating predictions afterward. The Ljung-Box test checks whether residual autocorrelation is statistically significant: if it is, the model is leaving predictable signal on the table.

import matplotlib.pyplot as plt
from statsmodels.stats.diagnostic import acorr_ljungbox

residuals = y_test.values - preds

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].plot(y_test.index, residuals)
axes[0].axhline(0, color='red', linestyle='--')
axes[0].set_title('Residuals over time')

axes[1].stem(range(1, 21),
             [pd.Series(residuals).autocorr(lag=i) for i in range(1, 21)])
axes[1].set_title('Residual autocorrelation')
plt.tight_layout(); plt.show()

lb = acorr_ljungbox(residuals, lags=10)
print(lb['lb_pvalue'])  # p < 0.05 = significant autocorrelation remains

Running this analysis after every model iteration keeps you honest about what patterns remain unmodelled. It is faster than looking at individual prediction errors and more actionable than a single MAE score — because it tells you not just how wrong the model is but where and in what pattern, pointing directly to what to fix next. Pair residual analysis with SHAP feature importance and you have a complete diagnostic loop: SHAP tells you what the model learned, residual patterns tell you what it failed to learn, and together they drive the next round of feature engineering.