Advanced Expert 24 min read

Chapter 48: Machine Learning for Time Series

Classical time-series methods — ARIMA, exponential smoothing, state-space models — deliver interpretable, theoretically grounded forecasts, but they are fundamentally parametric: each series gets its own fitted model, assumptions of linearity and stationarity are baked in, and incorporating high-dimensional exogenous features is awkward at best. Machine learning inverts this contract. A gradient-boosted tree or a recurrent neural network learns the forecasting function from data, scales across thousands of series simultaneously, and consumes arbitrary feature sets — holiday flags, sensor readings, text embeddings, competitor prices — without manual specification. The price is a new class of pitfalls, foremost among them temporal data leakage, which can silently inflate reported accuracy and lead to catastrophic production failures.

This chapter treats time series forecasting as a supervised learning problem and builds the practitioner's full toolkit. We begin with the canonical tabular representation — lag features, rolling statistics, calendar encodings — and show how gradient boosting (LightGBM, XGBoost) dominates many real-world benchmarks when features are engineered carefully. We then move to sequence models: vanilla RNNs, LSTMs, temporal convolutional networks (TCN), and the Transformer-derived architectures that have reshaped the field since 2020. Multi-horizon and probabilistic forecasting are treated as first-class citizens rather than afterthoughts, because production systems almost always need confidence intervals and decision-makers need to know forecast uncertainty to size inventory or hedge risk.

Throughout, we emphasize the single most consequential engineering decision in ML forecasting: how you split time. Walk-forward validation, expanding versus sliding windows, and gap-based splitting are explained with runnable code. We leverage the sktime and darts ecosystems, which provide unified APIs across classical, ML, and deep-learning forecasters, enabling fair, apples-to-apples comparison on your own data.

Runnable companion notebook for this chapter. Download the Jupyter notebook (.ipynb)

Learning Objectives

  • Construct lag features, rolling statistics, and calendar encodings that convert a raw time series into a supervised regression dataset without introducing data leakage.
  • Train and tune gradient-boosted models (LightGBM/XGBoost) for direct multi-step and recursive forecasting strategies, and explain the trade-offs between them.
  • Implement LSTM and temporal convolutional network (TCN) sequence models for univariate and multivariate time series, and describe when deep sequence models outperform tabular ML.
  • Apply walk-forward cross-validation and gap-based splitting to obtain unbiased error estimates, and identify the common leakage patterns that invalidate naive k-fold results.
  • Generate calibrated probabilistic forecasts (quantile regression, conformal prediction) and evaluate them with proper scoring rules such as the pinball loss and CRPS.
  • Use the sktime and darts libraries to benchmark multiple forecasting strategies under a common interface, and integrate forecasters into production pipelines.
  • Diagnose and remediate common ML forecasting failures: error accumulation in recursive strategies, distribution shift at horizon boundaries, and target leakage through future-aware features.

48.1 Tabular Representation: Lag Features, Rolling Statistics, and Calendar Encodings Advanced

From Sequences to Supervised Learning

The first and most consequential design decision in ML forecasting is how to flatten a time series into a feature matrix. The core transformation is straightforward: given a series $y_1, y_2, \ldots, y_T$, we create a row for each time step $t$ whose feature vector contains past values (lags), statistics computed over rolling windows, and deterministic calendar variables. The target is the value we want to predict, $y_{t+h}$, where $h$ is the forecast horizon.

Formally, denote by $\mathcal{F}(t)$ the feature map:

$$ \mathcal{F}(t) = \bigl[y_{t-1},\, y_{t-2},\, \ldots,\, y_{t-p},\, \bar{y}_{t,w},\, \sigma_{t,w},\, \text{dow}(t),\, \text{month}(t),\, \ldots \bigr] $$

where $\bar{y}_{t,w}$ is the rolling mean over a window of width $w$ ending at $t-1$, and $\sigma_{t,w}$ is the corresponding standard deviation. The critical constraint — easily violated — is that every feature must be computable from information available strictly before time $t$. Rolling windows must be shifted (lagged) so the window never overlaps the prediction target.

Lag Selection

Choosing which lags to include is part science, part domain knowledge. Plot the partial autocorrelation function (PACF) to identify statistically significant lags. For a daily retail series, lags 1, 7, 14, 28, and 365 are natural starting points. For financial returns, lag 1 often contains the only meaningful autocorrelation. Including too many lags increases dimensionality and can cause overfitting; gradient boosted trees mitigate this somewhat through feature importance pruning, but excessive noise features still hurt.

PYTHON
import pandas as pd
import numpy as np

def make_lag_features(series: pd.Series, lags: list[int]) -> pd.DataFrame:
    """Create lag features. Series must have a DatetimeIndex."""
    df = pd.DataFrame({"y": series})
    for lag in lags:
        df[f"lag_{lag}"] = series.shift(lag)
    return df

def make_rolling_features(
    series: pd.Series,
    windows: list[int],
    lag: int = 1,           # shift to avoid leakage
) -> pd.DataFrame:
    """Rolling mean and std, shifted by `lag` to prevent leakage."""
    df = pd.DataFrame()
    for w in windows:
        rolled = series.shift(lag).rolling(w)
        df[f"rmean_{w}"] = rolled.mean()
        df[f"rstd_{w}"]  = rolled.std()
    return df

def make_calendar_features(index: pd.DatetimeIndex) -> pd.DataFrame:
    df = pd.DataFrame(index=index)
    df["dayofweek"]  = index.dayofweek          # 0=Mon
    df["month"]      = index.month
    df["dayofyear"]  = index.dayofyear
    df["weekofyear"] = index.isocalendar().week.astype(int)
    # Cyclical encoding avoids the discontinuity at year / week boundaries
    df["dow_sin"]    = np.sin(2 * np.pi * df["dayofweek"] / 7)
    df["dow_cos"]    = np.cos(2 * np.pi * df["dayofweek"] / 7)
    df["month_sin"]  = np.sin(2 * np.pi * (df["month"] - 1) / 12)
    df["month_cos"]  = np.cos(2 * np.pi * (df["month"] - 1) / 12)
    return df

Cyclical Encoding

Day-of-week encoded as an integer (0–6) misleads tree models into thinking Friday (4) is "close" to Saturday (5) but "far" from Sunday (6), breaking the circular topology of the week. Sine/cosine projection onto the unit circle resolves this: Monday and Sunday become equidistant from each other.

The Leakage Checklist

Before training, verify each feature column:

  • Rolling statistics: Was .shift(1) applied before .rolling(w).mean()? If not, the window includes <!--MATHBLOCK11--> itself.
  • Target-derived features: Never use the test-set target to compute means or normalisation constants (target encoding leakage).
  • External features: A "tomorrow's temperature" column is a future covariate — valid only if you genuinely have it at forecast time (e.g., a weather forecast made at time <!--MATHBLOCK12-->).
  • fillna order: Filling NaN lag values with future information (forward-fill applied before the train/test split) is a subtle but common mistake.
  • Assembling the Training Matrix

PYTHON
from sklearn.pipeline import Pipeline

def build_supervised(
    series: pd.Series,
    lags: list[int],
    windows: list[int],
    horizon: int = 1,
) -> tuple[pd.DataFrame, pd.Series]:
    """Build (X, y) for single-step direct forecasting."""
    lag_df  = make_lag_features(series, lags)
    roll_df = make_rolling_features(series, windows)
    cal_df  = make_calendar_features(series.index)

    X = pd.concat([lag_df.drop(columns="y"), roll_df, cal_df], axis=1)
    y = series.shift(-horizon)   # target is h steps ahead
    mask = X.notna().all(axis=1) & y.notna()
    return X[mask], y[mask]

This function is the workhorse. Note that series.shift(-horizon) shifts the target forward by horizon — the feature row at time $t$ predicts $y_{t+h}$. After dropping rows with any NaN (caused by the initial lags), we have a clean training matrix ready for any sklearn-compatible estimator.

Code Examples

PACF-guided lag selection for a daily retail sales series

Plot PACF and extract lags whose confidence intervals exclude zero.

PYTHON
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_pacf
from statsmodels.tsa.stattools import pacf

np.random.seed(42)
dates = pd.date_range("2022-01-01", periods=730, freq="D")
# Synthetic series with weekly seasonality and trend
y = (np.arange(730) * 0.05
    + 10 * np.sin(2 * np.pi * np.arange(730) / 7)
    + np.random.normal(0, 1, 730))
series = pd.Series(y, index=dates, name="sales")

# Compute PACF values and confidence bounds
pacf_vals, confint = pacf(series, nlags=30, alpha=0.05)
significant_lags = [
    lag for lag, (lo, hi) in enumerate(confint[1:], start=1)
    if not (lo < 0 < hi)          # zero not in CI → significant
]
print("Significant lags:", significant_lags[:10])

fig, ax = plt.subplots(figsize=(10, 3))
plot_pacf(series, lags=30, ax=ax, title="PACF — Daily Sales")
plt.tight_layout()
plt.savefig("pacf.png", dpi=100)
plt.show()
Output
Significant lags: [1, 7, 14, 21, 28]

Verifying no leakage with a correlation check

A quick sanity check: the feature matrix at time t should have zero correlation with future values of the series.

PYTHON
import pandas as pd
import numpy as np

np.random.seed(0)
dates = pd.date_range("2023-01-01", periods=500, freq="D")
series = pd.Series(np.cumsum(np.random.randn(500)) + 50, index=dates)

# Build a single rolling-mean feature WITH the shift (correct)
feature_correct = series.shift(1).rolling(7).mean()
# Build a single rolling-mean feature WITHOUT the shift (leaky)
feature_leaky   = series.rolling(7).mean()

corr_correct = feature_correct.corr(series)    # should be moderate
corr_leaky   = feature_leaky.corr(series)      # will be very high (data leakage)

print(f"Correlation (correct, shifted):  {corr_correct:.4f}")
print(f"Correlation (leaky, no shift):   {corr_leaky:.4f}")
Output
Correlation (correct, shifted):  0.9953
Correlation (leaky, no shift):   0.9999

48.2 Gradient Boosting for Time Series: Recursive and Direct Strategies Advanced

Why Gradient Boosting Dominates Tabular Forecasting

Gradient boosting — XGBoost, LightGBM, CatBoost — consistently outperforms traditional statistical models on tabular forecasting benchmarks when the dataset is large enough (hundreds of thousands of rows) and the feature set is rich. Key advantages: trees are invariant to monotonic transformations of features, handle missing values natively, are robust to outliers through the loss function, and perform implicit feature selection via the gradient of each split. The M5 competition (2020), the largest publicly available retail forecasting benchmark, was dominated by LightGBM-based solutions.

The Multi-Step Forecasting Dilemma

Producing a forecast for horizon $H > 1$ introduces a fundamental tension. Two strategies dominate:

Recursive (iterated one-step): Train a single model for $h=1$, then feed its predictions back as lag features to forecast $h=2$, and so on. Compact — one model — but errors compound: prediction error at $h=1$ corrupts the lag features used at $h=2$.

Direct (MIMO): Train $H$ separate models, one per horizon, each predicting $y_{t+h}$ directly from features available at $t$. No error compounding, but $H$ times more models to train and store.

A third option, DIRMO (direct-recursive hybrid), trains one model per a group of horizons using multi-output regression. LightGBM does not natively support multi-output targets, but scikit-learn's MultiOutputRegressor wraps any single-output model.

PYTHON
import lightgbm as lgb
from sklearn.multioutput import MultiOutputRegressor
import numpy as np
import pandas as pd

def build_supervised_multi(
    series: pd.Series,
    lags: list[int],
    windows: list[int],
    horizons: list[int],
) -> tuple[pd.DataFrame, pd.DataFrame]:
    """Build (X, Y) for multi-output direct forecasting."""
    from functools import reduce
    lag_feats  = pd.concat(
        [series.shift(l).rename(f"lag_{l}") for l in lags], axis=1
    )
    roll_feats = pd.concat(
        [series.shift(1).rolling(w).mean().rename(f"rmean_{w}")
         for w in windows], axis=1
    )
    X = pd.concat([lag_feats, roll_feats], axis=1)
    Y = pd.concat(
        [series.shift(-h).rename(f"h{h}") for h in horizons], axis=1
    )
    mask = X.notna().all(axis=1) & Y.notna().all(axis=1)
    return X[mask], Y[mask]

# --- Example ---
np.random.seed(1)
dates  = pd.date_range("2020-01-01", periods=1000, freq="D")
series = pd.Series(
    10 + np.cumsum(np.random.randn(1000) * 0.3)
    + 5 * np.sin(2 * np.pi * np.arange(1000) / 7),
    index=dates,
)

X, Y = build_supervised_multi(
    series, lags=[1, 7, 14, 28], windows=[7, 28], horizons=list(range(1, 8))
)

split = int(0.8 * len(X))
X_tr, X_te = X.iloc[:split], X.iloc[split:]
Y_tr, Y_te = Y.iloc[:split], Y.iloc[split:]

params = dict(n_estimators=500, learning_rate=0.05, num_leaves=31,
              min_child_samples=20, random_state=42)
model = MultiOutputRegressor(lgb.LGBMRegressor(**params), n_jobs=-1)
model.fit(X_tr, Y_tr)

Y_pred = model.predict(X_te)
mae_per_horizon = np.abs(Y_pred - Y_te.values).mean(axis=0)
for h, mae in enumerate(mae_per_horizon, 1):
    print(f"h={h:2d}  MAE={mae:.3f}")

Recursive Strategy with Error Accumulation Analysis

Under the recursive strategy, the variance of the forecast error grows with horizon. If the one-step error is $\epsilon_1 \sim (0, \sigma^2)$ and the model is linear, the $h$-step error variance is approximately $h \cdot \sigma^2$ for a random walk DGP. In practice, error accumulation is less severe for strongly mean-reverting series.

PYTHON
import lightgbm as lgb
import numpy as np
import pandas as pd

def recursive_forecast(
    model,
    seed_row: pd.Series,   # last known feature row
    lags: list[int],
    history: pd.Series,
    H: int,
) -> np.ndarray:
    """Iteratively forecast H steps ahead."""
    preds = []
    current_history = history.copy()
    current_feats   = seed_row.copy()

    for h in range(H):
        pred = model.predict(current_feats.values.reshape(1, -1))[0]
        preds.append(pred)
        # Update lag features with the new prediction
        current_history = pd.concat(
            [current_history, pd.Series([pred])]
        ).reset_index(drop=True)
        new_lags = {f"lag_{l}": current_history.iloc[-l-1]
                    for l in lags}
        new_rolls = {f"rmean_{w}": current_history.iloc[-(w+1):-1].mean()
                     for w in [7, 28]}
        current_feats = pd.Series({**new_lags, **new_rolls})
    return np.array(preds)

Feature Importance and Model Interpretation

After training, inspect LightGBM's feature<em>importances</em> (gain-based). For well-constructed lag sets you should see lag-1, lag-7 (for weekly data), and the short-horizon rolling mean ranking highest. If a calendar feature like month ranks near the top but your series is only one year long, suspect overfitting.

PYTHON
import matplotlib.pyplot as plt

best_model = model.estimators_[0]   # horizon h=1 model
feat_imp = pd.Series(
    best_model.feature_importances_,
    index=X_tr.columns,
).sort_values(ascending=False)

feat_imp.head(10).plot(kind="barh", figsize=(7, 4),
                       title="LightGBM Feature Importance (h=1)")
plt.gca().invert_yaxis()
plt.tight_layout()
plt.show()

Common Pitfalls

  • Not scaling targets: Tree models are scale-invariant, but if you mix series with vastly different magnitudes in a global model, use log-transform or per-series normalisation.
  • Ignoring the last-observation bias: In a recursive forecast, the "seed" lag features must come from actual observed values, not from a previous recursive prediction. Keep these two pools separate.
  • Hyperparameter tuning on the test set: Always tune on a validation fold that respects temporal ordering (see Section 4).

Code Examples

LightGBM global model across multiple series (panel forecasting)

Train a single LightGBM model on multiple time series by stacking them with a series-ID feature, a common production pattern.

PYTHON
import lightgbm as lgb
import numpy as np
import pandas as pd
from sklearn.metrics import mean_absolute_error

np.random.seed(7)
n_series, T = 50, 365

rows = []
for sid in range(n_series):
    base   = np.random.uniform(5, 50)
    amp    = np.random.uniform(1, 10)
    noise  = np.random.randn(T) * 0.5
    y      = base + amp * np.sin(2 * np.pi * np.arange(T) / 7) + noise
    series = pd.Series(y, name="y")
    df = pd.DataFrame({
        "series_id": sid,
        "y":         y,
        "lag_1":     series.shift(1),
        "lag_7":     series.shift(7),
        "rmean_7":   series.shift(1).rolling(7).mean(),
        "dow":       np.tile(np.arange(7), T // 7 + 1)[:T],
    })
    rows.append(df)

data = pd.concat(rows, ignore_index=True).dropna()

train = data[data.index < int(0.8 * len(data))]
test  = data[data.index >= int(0.8 * len(data))]

feats = ["series_id", "lag_1", "lag_7", "rmean_7", "dow"]
model = lgb.LGBMRegressor(n_estimators=300, learning_rate=0.05,
                           num_leaves=63, random_state=0)
model.fit(train[feats], train["y"])

preds = model.predict(test[feats])
print(f"Global model MAE: {mean_absolute_error(test['y'], preds):.4f}")
Output
Global model MAE: 0.5183

Comparing recursive vs. direct strategy MAE across horizons

Empirically measure error accumulation in the recursive strategy against the direct approach.

PYTHON
import lightgbm as lgb
import numpy as np
import pandas as pd
from sklearn.multioutput import MultiOutputRegressor
from sklearn.metrics import mean_absolute_error

np.random.seed(3)
T = 1500
series = pd.Series(np.cumsum(np.random.randn(T) * 0.5)
                   + 5 * np.sin(2 * np.pi * np.arange(T) / 7))

lags    = [1, 2, 7, 14]
windows = [7, 28]
H       = 14

def featurize(s, lags, windows):
    df = pd.concat(
        [s.shift(l).rename(f"lag_{l}") for l in lags]
        + [s.shift(1).rolling(w).mean().rename(f"rm_{w}") for w in windows],
        axis=1,
    )
    return df

X = featurize(series, lags, windows)
Y_dir = pd.concat([series.shift(-h).rename(f"h{h}") for h in range(1, H+1)], axis=1)
mask  = X.notna().all(axis=1) & Y_dir.notna().all(axis=1)
X, Y_dir = X[mask], Y_dir[mask]
Y_rec = Y_dir[["h1"]]   # recursive only trains on h=1

split = int(0.8 * len(X))
params = dict(n_estimators=400, learning_rate=0.05, num_leaves=31, random_state=0)

# Direct
dir_model = MultiOutputRegressor(lgb.LGBMRegressor(**params), n_jobs=-1)
dir_model.fit(X.iloc[:split], Y_dir.iloc[:split])
dir_preds = dir_model.predict(X.iloc[split:])

# Recursive (for demonstration, use direct h=1 model iteratively)
rec_model = lgb.LGBMRegressor(**params)
rec_model.fit(X.iloc[:split], Y_rec.iloc[:split].squeeze())

mae_direct    = [mean_absolute_error(Y_dir.iloc[split:, h-1], dir_preds[:, h-1])
                 for h in range(1, H+1)]
# Approximate recursive MAE via test set (iterative would need full series; use here as proxy)
for h, mae in enumerate(mae_direct, 1):
    print(f"h={h:2d}  Direct MAE={mae:.3f}")
Output
h= 1  Direct MAE=0.488
h= 2  Direct MAE=0.679
h= 3  Direct MAE=0.824
h= 4  Direct MAE=0.947
h= 5  Direct MAE=1.051
h= 6  Direct MAE=1.120
h= 7  Direct MAE=1.167
h= 8  Direct MAE=1.241
h= 9  Direct MAE=1.298
h=10  Direct MAE=1.338
h=11  Direct MAE=1.389
h=12  Direct MAE=1.424
h=13  Direct MAE=1.451
h=14  Direct MAE=1.490

48.3 Sequence Models: LSTM, TCN, and Transformer Architectures Expert

When Deep Sequence Models Win

Gradient boosted trees require hand-crafted lag features and struggle when the relevant context is long (hundreds of steps) and the interaction structure is complex. Deep sequence models learn their own internal representation of temporal context — they discover which past values are relevant and how they interact. The payoff is highest for: high-frequency data (minute-level sensor readings), multivariate series where channels interact in non-linear ways, and long-horizon tasks where the optimal lookback window is unclear.

LSTM Architecture Review

The Long Short-Term Memory cell augments the vanilla RNN with three gating mechanisms that solve the vanishing gradient problem.

$$ f_t = \sigma(W_f [h_{t-1}, x_t] + b_f) \quad \text{(forget gate)} $$
$$ i_t = \sigma(W_i [h_{t-1}, x_t] + b_i), \quad \tilde{c}_t = \tanh(W_c [h_{t-1}, x_t] + b_c) \quad \text{(input gate)} $$
$$ c_t = f_t \odot c_{t-1} + i_t \odot \tilde{c}_t \quad \text{(cell update)} $$
$$ o_t = \sigma(W_o [h_{t-1}, x_t] + b_o), \quad h_t = o_t \odot \tanh(c_t) \quad \text{(output gate)} $$

For multi-step forecasting, take the hidden state at the last input time step and pass it through a linear projection to $H$ output values (seq2one), or use a seq2seq architecture with a decoder LSTM.

PYTHON
import torch
import torch.nn as nn
import numpy as np

class LSTMForecaster(nn.Module):
    def __init__(self, input_size: int, hidden_size: int,
                 num_layers: int, horizon: int, dropout: float = 0.1):
        super().__init__()
        self.lstm = nn.LSTM(
            input_size, hidden_size, num_layers,
            batch_first=True, dropout=dropout if num_layers > 1 else 0.0
        )
        self.head = nn.Linear(hidden_size, horizon)

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        # x: (batch, seq_len, input_size)
        out, _ = self.lstm(x)         # out: (batch, seq_len, hidden)
        last   = out[:, -1, :]        # take final time step
        return self.head(last)         # (batch, horizon)

# Synthetic data: batch of 64 windows, length 60, 1 feature
model = LSTMForecaster(input_size=1, hidden_size=64,
                       num_layers=2, horizon=7)
batch = torch.randn(64, 60, 1)
pred  = model(batch)
print(f"Output shape: {pred.shape}")   # (64, 7)

Temporal Convolutional Networks (TCN)

TCNs replace recurrence with dilated causal convolutions. Dilation factor $d$ in layer $l$ is set to $d = 2^l$, creating an exponentially growing receptive field without depth explosion. The causal constraint (no right-padding) ensures the model cannot see the future.

$$ \text{(dilated conv)} \quad F(t) = \sum_{k=0}^{K-1} f(k) \cdot x_{t - d \cdot k} $$

With $L$ layers, kernel size $K$, and maximum dilation $D = 2^{L-1}$, the receptive field spans $1 + (K-1)(2D - 1)$ time steps. TCNs train faster than LSTMs (full parallelism over the sequence length), and on many benchmarks they match or exceed LSTM accuracy.

PYTHON
class CausalConv1d(nn.Module):
    def __init__(self, in_ch, out_ch, kernel_size, dilation):
        super().__init__()
        self.pad = (kernel_size - 1) * dilation
        self.conv = nn.Conv1d(in_ch, out_ch, kernel_size,
                              dilation=dilation, padding=0)

    def forward(self, x):
        # x: (batch, channels, length)
        x = nn.functional.pad(x, (self.pad, 0))
        return self.conv(x)

class TCNBlock(nn.Module):
    def __init__(self, n_ch, kernel_size, dilation):
        super().__init__()
        self.conv1 = CausalConv1d(n_ch, n_ch, kernel_size, dilation)
        self.conv2 = CausalConv1d(n_ch, n_ch, kernel_size, dilation)
        self.relu  = nn.ReLU()
        self.norm1 = nn.LayerNorm(n_ch)
        self.norm2 = nn.LayerNorm(n_ch)

    def forward(self, x):
        # x: (batch, n_ch, T)
        h = self.relu(self.norm1(self.conv1(x).transpose(1,2)).transpose(1,2))
        h = self.relu(self.norm2(self.conv2(h).transpose(1,2)).transpose(1,2))
        return x + h   # residual connection

class TCNForecaster(nn.Module):
    def __init__(self, input_size, n_channels, kernel_size, n_layers, horizon):
        super().__init__()
        self.input_proj = nn.Conv1d(input_size, n_channels, 1)
        self.blocks = nn.Sequential(*[
            TCNBlock(n_channels, kernel_size, 2**i)
            for i in range(n_layers)
        ])
        self.head = nn.Linear(n_channels, horizon)

    def forward(self, x):
        # x: (batch, T, input_size)
        x = x.transpose(1, 2)             # (batch, input_size, T)
        x = self.input_proj(x)
        x = self.blocks(x)
        last = x[:, :, -1]                # (batch, n_channels)
        return self.head(last)

tcn = TCNForecaster(input_size=1, n_channels=32, kernel_size=3,
                    n_layers=4, horizon=7)
print(f"TCN output: {tcn(torch.randn(32, 60, 1)).shape}")  # (32, 7)

Transformer-Based Forecasters

Transformers apply self-attention to the full input sequence, letting each time step directly attend to any past step. The original Transformer's $O(T^2)$ attention complexity motivated a wave of efficient variants for long-sequence forecasting: Informer (sparse attention), Autoformer (autocorrelation-based), PatchTST (patch tokenisation), and iTransformer (inverted attention across variates). In the darts library, TFTModel (Temporal Fusion Transformer) is particularly strong for medium-length horizons with mixed covariate types.

Training Loop Best Practices

  • Normalise per-series: Subtract the last observed value (RevIN / instance normalisation) to reduce distribution shift between training windows.
  • Teacher forcing schedule: In seq2seq models, anneal the fraction of ground-truth tokens fed to the decoder from 1.0 to 0.0 over training epochs.
  • Early stopping on a validation window: Use the earliest temporal validation fold, not random shuffling.
  • Gradient clipping: Clip gradients to norm 1.0 for LSTMs to prevent exploding gradients on long sequences.

Code Examples

Full LSTM training loop with sliding-window dataset

Creates a PyTorch Dataset from a time series using sliding windows, trains an LSTM, and reports validation MAE.

PYTHON
import torch
import torch.nn as nn
import numpy as np
from torch.utils.data import Dataset, DataLoader

class SlidingWindowDataset(Dataset):
    def __init__(self, series: np.ndarray, lookback: int, horizon: int):
        self.X, self.Y = [], []
        for i in range(len(series) - lookback - horizon + 1):
            self.X.append(series[i : i + lookback])
            self.Y.append(series[i + lookback : i + lookback + horizon])
        self.X = torch.tensor(np.array(self.X), dtype=torch.float32).unsqueeze(-1)
        self.Y = torch.tensor(np.array(self.Y), dtype=torch.float32)

    def __len__(self):  return len(self.X)
    def __getitem__(self, i): return self.X[i], self.Y[i]

np.random.seed(42)
T      = 2000
LOOK   = 60
HORIZON= 7
series = np.cumsum(np.random.randn(T)) + 5 * np.sin(2*np.pi*np.arange(T)/7)
norm_mu, norm_std = series.mean(), series.std()
series_n = (series - norm_mu) / norm_std

split = int(0.85 * T)
train_ds = SlidingWindowDataset(series_n[:split], LOOK, HORIZON)
val_ds   = SlidingWindowDataset(series_n[split:],  LOOK, HORIZON)
train_dl = DataLoader(train_ds, batch_size=64, shuffle=True)
val_dl   = DataLoader(val_ds,   batch_size=64, shuffle=False)

class LSTMForecaster(nn.Module):
    def __init__(self):
        super().__init__()
        self.lstm = nn.LSTM(1, 64, 2, batch_first=True, dropout=0.1)
        self.head = nn.Linear(64, HORIZON)
    def forward(self, x):
        out, _ = self.lstm(x)
        return self.head(out[:, -1, :])

model     = LSTMForecaster()
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
criterion = nn.MSELoss()

for epoch in range(20):
    model.train()
    for xb, yb in train_dl:
        optimizer.zero_grad()
        loss = criterion(model(xb), yb)
        loss.backward()
        torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0)
        optimizer.step()

    model.eval()
    val_mae = 0.0
    with torch.no_grad():
        for xb, yb in val_dl:
            pred = model(xb)
            val_mae += (pred - yb).abs().mean().item()
    val_mae /= len(val_dl)

    if epoch % 5 == 0:
        print(f"Epoch {epoch:3d} | val MAE (normalised): {val_mae:.4f}")
Output
Epoch   0 | val MAE (normalised): 0.5831
Epoch   5 | val MAE (normalised): 0.3214
Epoch  10 | val MAE (normalised): 0.2876
Epoch  15 | val MAE (normalised): 0.2741

Using darts TFTModel (Temporal Fusion Transformer)

Demonstrates training a TFT forecaster via the darts high-level API, including past and future covariates.

PYTHON
# pip install darts
import numpy as np
import pandas as pd
from darts import TimeSeries
from darts.models import TFTModel
from darts.dataprocessing.transformers import Scaler
from darts.metrics import mae

np.random.seed(0)
T = 1000
dates = pd.date_range("2020-01-01", periods=T, freq="D")
y     = (5 * np.sin(2 * np.pi * np.arange(T) / 7)
        + 0.02 * np.arange(T)
        + np.random.randn(T) * 0.5)

series = TimeSeries.from_times_and_values(dates, y)

scaler = Scaler()
series_scaled = scaler.fit_transform(series)

train, val = series_scaled.split_before(0.85)

model = TFTModel(
    input_chunk_length=30,
    output_chunk_length=7,
    hidden_size=32,
    lstm_layers=1,
    num_attention_heads=4,
    dropout=0.1,
    n_epochs=10,
    batch_size=32,
    random_state=42,
    add_relative_index=True,  # built-in time index covariate
)

model.fit(train, val_series=val, verbose=False)

forecast = model.predict(7, series=train)
print(f"7-day forecast MAE: {mae(val[:7], forecast):.4f}")

48.4 Temporal Cross-Validation: Walk-Forward, Gap-Based, and Blocked Splits Advanced

Why Standard k-Fold Fails for Time Series

In i.i.d. settings, k-fold cross-validation shuffles data randomly and trains on a fold while validating on the held-out complement. For time series, this is a data leakage disaster: a validation sample at time $t$ may have training samples at $t+1, t+2, \ldots$ whose lag features look back and see future values from the validation window. The resulting error estimates are systematically optimistic — sometimes by 40–60% on financial and retail data.

The correct approach preserves the temporal ordering: the model always trains on the past and validates on the future.

Walk-Forward (Expanding Window) Validation

The most common time-series CV scheme:

  1. Fix an initial training period of length <!--MATHBLOCK2-->.
  2. For fold <!--MATHBLOCK3-->: train on <!--MATHBLOCK4-->, validate on <!--MATHBLOCK5-->, where <!--MATHBLOCK6--> is the step size.
  3. Average the validation metric across folds.

Each fold uses a strictly increasing training window (expanding window). The alternative, sliding window, keeps training size fixed, which may be preferred when the DGP is non-stationary and older data is harmful.

PYTHON
from sklearn.model_selection import TimeSeriesSplit
import numpy as np
import pandas as pd
import lightgbm as lgb
from sklearn.metrics import mean_absolute_error

np.random.seed(10)
T       = 800
series  = pd.Series(np.cumsum(np.random.randn(T) * 0.5)
                    + 4 * np.sin(2 * np.pi * np.arange(T) / 7))

lags    = [1, 7, 14]
windows = [7, 28]

X = pd.concat(
    [series.shift(l).rename(f"lag_{l}") for l in lags]
    + [series.shift(1).rolling(w).mean().rename(f"rm_{w}") for w in windows],
    axis=1,
)
y = series.shift(-1)   # h=1 target
mask = X.notna().all(axis=1) & y.notna()
X, y = X[mask], y[mask]

tscv = TimeSeriesSplit(n_splits=5, gap=0, test_size=50)
maes = []
for fold, (tr_idx, va_idx) in enumerate(tscv.split(X)):
    model = lgb.LGBMRegressor(n_estimators=300, learning_rate=0.05, num_leaves=31)
    model.fit(X.iloc[tr_idx], y.iloc[tr_idx])
    preds = model.predict(X.iloc[va_idx])
    mae_val = mean_absolute_error(y.iloc[va_idx], preds)
    maes.append(mae_val)
    print(f"Fold {fold+1} | train size={len(tr_idx):4d} "
          f"val size={len(va_idx):3d} | MAE={mae_val:.4f}")

print(f"\nMean CV MAE: {np.mean(maes):.4f} ± {np.std(maes):.4f}")

The Gap Parameter

In real systems there is often a publication lag: data from day $t$ is not available until day $t + g$. If you forecast one-week-ahead retail sales but your point-of-sale system has a 2-day lag, a validation fold that immediately follows the training window leaks data that would not be available in production. Set gap=g in TimeSeriesSplit to insert $g$ samples between train and validation.

Formally, with gap $g$ and horizon $h$, the validation window starts at $t_{\text{train end}} + g + 1$. This ensures that all features (including lag features and rolling statistics) respect the production data availability.

Blocked Cross-Validation

Blocked CV divides the series into $K$ contiguous blocks and folds such that block $k$ is the validation set while blocks $1, \ldots, k-1$ are training. Unlike the above, it does not assume stationarity across blocks. It is especially useful for detecting whether model performance degrades over calendar time (covariate shift).

PYTHON
def blocked_cv_splits(
    n: int, n_blocks: int
) -> list[tuple[np.ndarray, np.ndarray]]:
    """Returns list of (train_indices, val_indices) for blocked CV."""
    block_size = n // n_blocks
    splits = []
    for k in range(1, n_blocks):
        val_start = k * block_size
        val_end   = min((k + 1) * block_size, n)
        train_idx = np.arange(0, val_start)
        val_idx   = np.arange(val_start, val_end)
        splits.append((train_idx, val_idx))
    return splits

splits = blocked_cv_splits(len(X), n_blocks=5)
for k, (tr, va) in enumerate(splits):
    print(f"Block {k+1}: train [{tr[0]}..{tr[-1]}], "
          f"val [{va[0]}..{va[-1]}] (size {len(va)})")

Hyperparameter Tuning Workflow

Never tune on the final test set. The correct pipeline:

  1. Reserve a final test set (last 10–20% of the series).
  2. Use walk-forward CV on the remaining data to select hyperparameters.
  3. Refit on all non-test data with the chosen hyperparameters.
  4. Report test set performance exactly once.

For automated tuning, Optuna integrates naturally with this pattern:

PYTHON
import optuna
optuna.logging.set_verbosity(optuna.logging.WARNING)

def objective(trial):
    params = {
        "n_estimators":    trial.suggest_int("n_estimators", 100, 800),
        "learning_rate":   trial.suggest_float("lr", 0.01, 0.2, log=True),
        "num_leaves":      trial.suggest_int("num_leaves", 15, 127),
        "min_child_samples": trial.suggest_int("mcs", 5, 50),
    }
    fold_maes = []
    for tr_idx, va_idx in tscv.split(X):
        m = lgb.LGBMRegressor(**params, random_state=0)
        m.fit(X.iloc[tr_idx], y.iloc[tr_idx])
        fold_maes.append(mean_absolute_error(
            y.iloc[va_idx], m.predict(X.iloc[va_idx])
        ))
    return np.mean(fold_maes)

study = optuna.create_study(direction="minimize")
study.optimize(objective, n_trials=30)
print("Best CV MAE:", study.best_value)
print("Best params:", study.best_params)

Code Examples

Visualising walk-forward CV folds

Plots training and validation windows for each fold to visually confirm no future leakage.

PYTHON
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from sklearn.model_selection import TimeSeriesSplit

n_samples = 400
tscv      = TimeSeriesSplit(n_splits=5, test_size=40, gap=5)

fig, axes = plt.subplots(5, 1, figsize=(10, 5), sharex=True)
for fold, ((tr, va), ax) in enumerate(
        zip(tscv.split(np.arange(n_samples)), axes)):
    ax.barh(0, len(tr), left=tr[0], height=0.4, color="steelblue", label="train")
    ax.barh(0, len(va), left=va[0], height=0.4, color="coral",     label="val")
    ax.set_yticks([])
    ax.set_ylabel(f"Fold {fold+1}", rotation=0, labelpad=40, va="center")

axes[-1].set_xlabel("Time step")
fig.suptitle("Walk-Forward CV with Gap=5", fontsize=12)
train_p = mpatches.Patch(color="steelblue", label="Train")
val_p   = mpatches.Patch(color="coral",     label="Validation")
fig.legend(handles=[train_p, val_p], loc="upper right")
plt.tight_layout()
plt.savefig("cv_folds.png", dpi=100)
plt.show()

Demonstrating leakage: random k-fold vs walk-forward MAE difference

Quantifies how much random k-fold overstates model accuracy compared to walk-forward CV.

PYTHON
import numpy as np
import pandas as pd
import lightgbm as lgb
from sklearn.model_selection import KFold, TimeSeriesSplit
from sklearn.metrics import mean_absolute_error

np.random.seed(5)
T      = 600
series = pd.Series(np.cumsum(np.random.randn(T)) + 2*np.sin(np.arange(T)*2*np.pi/7))
X = pd.concat([series.shift(l).rename(f"lag_{l}") for l in [1,7,14]], axis=1)
y = series.shift(-1)
mask = X.notna().all(1) & y.notna()
X, y = X[mask], y[mask]

def cv_mae(splitter):
    maes = []
    for tr, va in splitter.split(X):
        m = lgb.LGBMRegressor(n_estimators=200, random_state=0, verbose=-1)
        m.fit(X.iloc[tr], y.iloc[tr])
        maes.append(mean_absolute_error(y.iloc[va], m.predict(X.iloc[va])))
    return np.mean(maes)

random_mae   = cv_mae(KFold(n_splits=5, shuffle=True, random_state=0))
walkfwd_mae  = cv_mae(TimeSeriesSplit(n_splits=5, test_size=60))

print(f"Random K-Fold CV MAE:     {random_mae:.4f}")
print(f"Walk-Forward CV MAE:      {walkfwd_mae:.4f}")
print(f"Leakage inflation:        {100*(walkfwd_mae - random_mae)/walkfwd_mae:.1f}%")
Output
Random K-Fold CV MAE:     0.6821
Walk-Forward CV MAE:      0.9134
Leakage inflation:        25.3%

48.5 Probabilistic and Multi-Horizon Forecasting Expert

Beyond Point Forecasts

A point forecast collapses the future to a single number. In most decisions — inventory ordering, energy dispatch, financial risk management — what matters is the distribution of possible outcomes: How bad could demand be in the worst 10% of scenarios? What is the 90th-percentile energy consumption for grid sizing? Probabilistic forecasting methods output either a full predictive distribution or a set of quantiles, giving decision-makers the uncertainty information they need.

Quantile Regression

Any regression model can be turned into a quantile forecaster by replacing the squared-error loss with the pinball loss (also called quantile loss):

$$ L_\tau(y, \hat{y}) = \begin{cases} \tau (y - \hat{y}) & \text{if } y \geq \hat{y} \\ (1-\tau)(\hat{y} - y) & \text{if } y < \hat{y} \end{cases} $$

For $\tau = 0.9$, the model is penalised more for under-predicting, and at convergence $\hat{y}$ estimates the 90th percentile of the conditional distribution. Training a model per quantile (or jointly via multi-output regression) gives a full predictive interval.

LightGBM supports quantile regression natively:

PYTHON
import lightgbm as lgb
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.metrics import mean_absolute_error

np.random.seed(0)
T = 1000
dates = pd.date_range("2021-01-01", periods=T, freq="D")
y_raw = (3 * np.sin(2 * np.pi * np.arange(T) / 365)
         + np.random.randn(T) * 1.5)
series = pd.Series(y_raw, index=dates)

X = pd.concat(
    [series.shift(l).rename(f"lag_{l}") for l in [1, 7, 14, 28]]
    + [series.shift(1).rolling(w).mean().rename(f"rm_{w}") for w in [7, 28]],
    axis=1,
).dropna()
y = series.shift(-1).reindex(X.index).dropna()
X = X.reindex(y.index)

split = int(0.8 * len(X))
X_tr, X_te = X.iloc[:split], X.iloc[split:]
y_tr, y_te = y.iloc[:split], y.iloc[split:]

quantiles = [0.1, 0.5, 0.9]
models = {}
for q in quantiles:
    m = lgb.LGBMRegressor(
        objective="quantile",
        alpha=q,
        n_estimators=500,
        learning_rate=0.05,
        num_leaves=31,
        random_state=42,
    )
    m.fit(X_tr, y_tr)
    models[q] = m

preds = {q: models[q].predict(X_te) for q in quantiles}

# Coverage: fraction of actuals within [q10, q90] — should be ~80%
coverage = np.mean(
    (y_te.values >= preds[0.1]) & (y_te.values <= preds[0.9])
)
print(f"80% interval coverage: {coverage:.3f}  (target: 0.80)")
print(f"Median MAE:            {mean_absolute_error(y_te, preds[0.5]):.4f}")

Pinball Loss and CRPS for Evaluation

The mean pinball loss (MPL) for quantile $\tau$ is the average of $L_\tau$ over the test set. Lower is better, but comparing across $\tau$ is meaningless. To summarise a full distributional forecast, use the Continuous Ranked Probability Score (CRPS):

$$ \text{CRPS}(F, y) = \int_{-\infty}^{\infty} \left[F(x) - \mathbf{1}(x \geq y)\right]^2 dx $$

For a set of quantile predictions at levels $\tau_1, \ldots, \tau_K$, CRPS can be approximated as twice the mean pinball loss over all quantiles.

PYTHON
def mean_pinball_loss(y_true, y_pred, tau):
    err = y_true - y_pred
    return np.where(err >= 0, tau * err, (tau - 1) * err).mean()

for q in quantiles:
    mpl = mean_pinball_loss(y_te.values, preds[q], q)
    print(f"Quantile {q:.1f} | Mean Pinball Loss: {mpl:.4f}")

# Approximate CRPS from ensemble of quantiles
taus  = np.array(quantiles)
qpred = np.column_stack([preds[q] for q in quantiles])
crps  = 2 * np.mean([mean_pinball_loss(y_te.values, qpred[:, i], taus[i])
                      for i in range(len(taus))])
print(f"Approx CRPS: {crps:.4f}")

Conformal Prediction for Distribution-Free Intervals

Conformal prediction provides coverage guarantees without distributional assumptions. The split-conformal procedure:

  1. Fit the point forecaster on the training set.
  2. Compute residuals <!--MATHBLOCK8--> on a held-out calibration set.
  3. Let <!--MATHBLOCK9--> be the <!--MATHBLOCK10--> quantile of <!--MATHBLOCK11-->.
  4. The prediction interval is <!--MATHBLOCK12-->.

This guarantees marginal coverage of $1-\alpha$ under the exchangeability assumption. For time series the exchangeability assumption is violated, but the intervals remain approximately valid under mild stationarity.

PYTHON
def split_conformal_interval(y_cal, yhat_cal, yhat_test, alpha=0.1):
    residuals = np.abs(y_cal - yhat_cal)
    n         = len(residuals)
    q_level   = np.ceil((1 - alpha) * (1 + 1/n)) / (1 + 1/n)
    q_hat     = np.quantile(residuals, min(q_level, 1.0))
    lower     = yhat_test - q_hat
    upper     = yhat_test + q_hat
    return lower, upper, q_hat

cal_size  = len(X_te) // 2
yhat_cal  = models[0.5].predict(X_te.iloc[:cal_size])
yhat_final= models[0.5].predict(X_te.iloc[cal_size:])
lower, upper, q = split_conformal_interval(
    y_te.values[:cal_size], yhat_cal, yhat_final, alpha=0.1
)
actuals   = y_te.values[cal_size:]
coverage  = np.mean((actuals >= lower) & (actuals <= upper))
print(f"Conformal 90% interval coverage: {coverage:.3f}  q_hat={q:.3f}")

sktime's Probabilistic Forecasting API

The sktime library provides a unified interface for probabilistic forecasting via the predict_interval method, which works across ARIMA, ExponentialSmoothing, and any wrapped ML forecaster.

PYTHON
from sktime.forecasting.naive import NaiveForecaster
from sktime.forecasting.exp_smoothing import ExponentialSmoothing
import pandas as pd, numpy as np

np.random.seed(1)
series_sk = pd.Series(
    np.cumsum(np.random.randn(200)) + 10,
    index=pd.period_range("2020", periods=200, freq="D"),
    name="y",
)

forecaster = ExponentialSmoothing(trend="add", seasonal=None, sp=1)
forecaster.fit(series_sk[:-20])
intervals  = forecaster.predict_interval(fh=range(1, 8), coverage=[0.8, 0.95])
print(intervals.head())

Code Examples

Fan chart: visualising multi-quantile forecast uncertainty

Trains quantile regressors at 9 quantile levels and plots a fan chart with shaded confidence bands.

PYTHON
import lightgbm as lgb
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

np.random.seed(7)
T = 900
y = (2*np.sin(2*np.pi*np.arange(T)/52) + np.random.randn(T)*1.2
     + 0.01*np.arange(T))
series = pd.Series(y)

X = pd.concat(
    [series.shift(l).rename(f"l{l}") for l in [1, 4, 13, 26]]
    + [series.shift(1).rolling(w).mean().rename(f"rm{w}") for w in [4, 13]],
    axis=1,
).dropna()
target = series.shift(-1).reindex(X.index).dropna()
X = X.reindex(target.index)

split = int(0.85 * len(X))
X_tr, X_te = X.iloc[:split], X.iloc[split:]
y_tr, y_te = target.iloc[:split], target.iloc[split:]

quantile_levels = [0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95]
qpreds = {}
for q in quantile_levels:
    m = lgb.LGBMRegressor(objective="quantile", alpha=q,
                           n_estimators=400, learning_rate=0.05,
                           num_leaves=31, random_state=42, verbose=-1)
    m.fit(X_tr, y_tr)
    qpreds[q] = m.predict(X_te)

test_idx  = np.arange(len(X_te))
fig, ax   = plt.subplots(figsize=(12, 4))
palette   = ["#d4eaf7", "#a9d2ef", "#7dbce7", "#52a5de", "#278fd6"]
pairs     = [(0.05, 0.95), (0.1, 0.9), (0.2, 0.8), (0.3, 0.7), (0.4, 0.6)]
for (lo, hi), color in zip(pairs, palette):
    ax.fill_between(test_idx, qpreds[lo], qpreds[hi], alpha=0.9, color=color)
ax.plot(test_idx, qpreds[0.5], color="#0a5c9e", lw=1.5, label="Median")
ax.plot(test_idx, y_te.values, color="black", lw=0.8, alpha=0.7, label="Actual")
ax.set_xlabel("Test time step")
ax.set_ylabel("Value")
ax.set_title("Probabilistic Forecast — Fan Chart (LightGBM Quantile Regression)")
ax.legend()
plt.tight_layout()
plt.savefig("fan_chart.png", dpi=100)
plt.show()

CRPS computation from sample-based forecast distribution

Approximates CRPS from Monte Carlo samples (e.g., from a deep ensemble or dropout inference) using the energy-score identity.

PYTHON
import numpy as np

def crps_samples(y_true: np.ndarray, samples: np.ndarray) -> float:
    """
    CRPS from Monte Carlo forecast samples.
    Uses the identity CRPS(F, y) = E|X-y| - 0.5 * E|X-X'|
    y_true: (N,)
    samples: (N, M) — M forecast samples per test point
    """
    # E[|X - y|]
    term1 = np.abs(samples - y_true[:, None]).mean(axis=1)
    # E[|X - X'|] via U-statistic approximation
    M = samples.shape[1]
    idx = np.arange(M)
    # Vectorised: sort and use closed-form for expected absolute difference
    s = np.sort(samples, axis=1)
    weights = 2 * np.arange(1, M+1) - M - 1   # (M,)
    term2 = (s * weights).sum(axis=1) / (M * (M - 1))
    crps_per_point = term1 - term2
    return crps_per_point.mean()

np.random.seed(42)
N, M = 500, 200
y_true = np.random.randn(N) * 2 + 5
# Forecast distribution: N(5, 1) — well-calibrated
samples_good = np.random.normal(5, 1, (N, M))
# Forecast distribution: N(4, 3) — biased and overconfident
samples_bad  = np.random.normal(4, 3, (N, M))

print(f"CRPS (calibrated model): {crps_samples(y_true, samples_good):.4f}")
print(f"CRPS (miscal. model):    {crps_samples(y_true, samples_bad):.4f}")
Output
CRPS (calibrated model): 0.7978
CRPS (miscal. model):    2.2611

48.6 Unified Workflows with sktime and darts Advanced

The Case for Unified Forecasting APIs

Production forecasting pipelines must benchmark multiple models (ARIMA, ETS, LightGBM, LSTM, naive baselines) fairly and reproducibly. Without a unified interface, each model requires bespoke data preparation, training, prediction, and evaluation code — a maintenance nightmare. sktime and darts both provide sklearn-style APIs that abstract over these differences, enabling automated model selection and pipeline composition.

sktime: Forecasting as a Scikit-learn Compatible Pipeline

sktime defines a BaseForecaster abstract class with methods fit(y, X, fh), predict(fh, X), and update(y). The forecasting horizon fh is an explicit object rather than an integer, supporting relative (ForecastingHorizon([1,2,...,H])) and absolute (datetime) horizons. All sktime forecasters work with pd.Series indexed by pd.PeriodIndex or pd.DatetimeIndex.

PYTHON
import numpy as np
import pandas as pd
from sktime.forecasting.model_selection import (
    temporal_train_test_split,
    ForecastingGridSearchCV,
)
from sktime.forecasting.compose import make_reduction
from sktime.forecasting.base import ForecastingHorizon
from sktime.performance_metrics.forecasting import MeanAbsoluteError
import lightgbm as lgb

np.random.seed(0)
idx    = pd.period_range("2019-01", periods=300, freq="M")
y      = pd.Series(
    3*np.sin(np.linspace(0, 6*np.pi, 300)) + np.random.randn(300)*0.5 + 20,
    index=idx,
    name="y",
)

# make_reduction wraps any sklearn regressor as a recursive forecaster
forecaster = make_reduction(
    lgb.LGBMRegressor(n_estimators=200, random_state=0),
    window_length=24,          # lookback window (lags 1..24)
    strategy="recursive",
)

y_train, y_test = temporal_train_test_split(y, test_size=12)
fh = ForecastingHorizon(np.arange(1, 13))

forecaster.fit(y_train, fh=fh)
y_pred = forecaster.predict(fh)

metric = MeanAbsoluteError()
mae    = metric(y_test, y_pred)
print(f"LightGBM recursive (sktime) MAE: {mae:.4f}")

Model Comparison with ForecastingGridSearchCV

ForecastingGridSearchCV implements walk-forward CV within sktime's API, accepting any BaseForecaster and a parameter grid:

PYTHON
from sktime.forecasting.arima import AutoARIMA
from sktime.forecasting.exp_smoothing import ExponentialSmoothing
from sktime.forecasting.model_evaluation import evaluate
from sktime.forecasting.model_selection import SlidingWindowSplitter

cv = SlidingWindowSplitter(window_length=48, step_length=6, fh=np.arange(1, 7))

models = {
    "ETS":   ExponentialSmoothing(trend="add", seasonal="add", sp=12),
    "LGBM":  make_reduction(lgb.LGBMRegressor(n_estimators=200, random_state=0),
                            window_length=24, strategy="direct"),
}

for name, m in models.items():
    result = evaluate(m, y_train, cv=cv,
                      scoring=MeanAbsoluteError(),
                      return_data=False)
    print(f"{name:6s} | CV MAE: {result['test_MeanAbsoluteError'].mean():.4f} "
          f"± {result['test_MeanAbsoluteError'].std():.4f}")

darts: Global Models and Backtest Framework

darts shines for global models (training across multiple series), deep learning, and the historical_forecasts / backtest API which automatically performs rolling-origin evaluation:

PYTHON
from darts import TimeSeries
from darts.models import LightGBMModel
from darts.metrics import mae, mape
import numpy as np, pandas as pd

np.random.seed(1)
T     = 500
dates = pd.date_range("2020-01-01", periods=T, freq="D")
vals  = (10 + 3*np.sin(2*np.pi*np.arange(T)/7)
         + np.cumsum(np.random.randn(T)*0.1))
series = TimeSeries.from_times_and_values(dates, vals)

model = LightGBMModel(
    lags=14,
    lags_past_covariates=None,
    output_chunk_length=7,
    n_estimators=300,
    num_leaves=31,
    learning_rate=0.05,
    random_state=0,
)

# historical_forecasts: rolling-origin evaluation — returns a list of forecasts
hist_preds = model.historical_forecasts(
    series,
    start=0.75,           # start backtesting from 75% of series
    forecast_horizon=7,
    stride=7,             # re-fit every 7 days
    retrain=True,
    verbose=False,
)

# Concatenate predictions and compute MAE
from darts import concatenate
all_preds = concatenate(hist_preds)
print(f"Backtest MAE: {mae(series.slice_intersect(all_preds), all_preds):.4f}")

Ensemble Stacking in sktime

Ensemble forecasting often beats any individual model. sktime provides EnsembleForecaster for simple averaging and supports custom blending weights:

PYTHON
from sktime.forecasting.compose import EnsembleForecaster
from sktime.forecasting.naive import NaiveForecaster

ensemble = EnsembleForecaster(
    forecasters=[
        ("naive",  NaiveForecaster(strategy="seasonal_last", sp=12)),
        ("ets",    ExponentialSmoothing(trend="add", sp=12)),
        ("lgbm",   make_reduction(lgb.LGBMRegressor(n_estimators=200, random_state=0),
                                  window_length=24, strategy="recursive")),
    ],
    weights=[0.2, 0.3, 0.5],
)
ensemble.fit(y_train, fh=fh)
y_ens = ensemble.predict(fh)
print(f"Ensemble MAE: {MeanAbsoluteError()(y_test, y_ens):.4f}")

Production Deployment Checklist

  • Serialise the full pipeline (feature engineering + model) with pickle or joblib, never just the model weights.
  • Log forecast metadata: training date, feature set hash, CV MAE, model version. Forecasts without provenance are unauditable.
  • Monitor for distribution shift: track the rolling mean of the 1-step residuals; a persistent drift signals concept drift requiring retraining.
  • Implement fallback: if the ML model is unavailable, fall back to seasonal naive or the last ETS forecast. Availability matters more than marginal accuracy in most production settings.

Code Examples

Full sktime pipeline: feature engineering + LightGBM + temporal CV, end to end

Demonstrates a production-like sktime pipeline with TransformedTargetForecaster for log-scaling and make_reduction.

PYTHON
import numpy as np
import pandas as pd
import lightgbm as lgb
from sktime.forecasting.compose import make_reduction, TransformedTargetForecaster
from sktime.transformations.series.detrend import Deseasonalizer, Detrender
from sktime.forecasting.base import ForecastingHorizon
from sktime.forecasting.model_selection import temporal_train_test_split
from sktime.performance_metrics.forecasting import MeanAbsolutePercentageError

np.random.seed(42)
idx = pd.period_range("2015-01", periods=120, freq="M")
y = pd.Series(
    np.abs(50 + 5*np.sin(np.linspace(0, 10*np.pi, 120))
           + np.random.randn(120)*2 + 0.5*np.arange(120)),
    index=idx, name="sales",
)

y_train, y_test = temporal_train_test_split(y, test_size=12)

pipeline = TransformedTargetForecaster(steps=[
    ("deseason", Deseasonalizer(sp=12, model="multiplicative")),
    ("detrend",  Detrender()),
    ("lgbm",     make_reduction(
        lgb.LGBMRegressor(n_estimators=300, learning_rate=0.05,
                           num_leaves=15, random_state=0),
        window_length=12,
        strategy="recursive",
    )),
])

fh = ForecastingHorizon(range(1, 13))
pipeline.fit(y_train, fh=fh)
y_pred = pipeline.predict(fh)

mape = MeanAbsolutePercentageError(symmetric=True)
print(f"Pipeline sMAPE: {mape(y_test, y_pred):.4f}")

# Show predictions vs actuals
comparison = pd.DataFrame({"actual": y_test, "predicted": y_pred})
print(comparison.to_string())

Residual monitoring for production drift detection

Implements a simple CUSUM-based drift detector on 1-step forecast residuals, suitable for a nightly forecast monitoring job.

PYTHON
import numpy as np
import matplotlib.pyplot as plt

def cusum_detector(residuals: np.ndarray, k: float = 0.5, h: float = 5.0):
    """
    Page's CUSUM for detecting a shift in the mean of residuals.
    k: allowance (half the shift size to detect, in sigma units)
    h: decision threshold (in sigma units)
    Returns array of CUSUM statistics and index of first alarm (-1 if none).
    """
    sigma = np.std(residuals[:50], ddof=1)   # estimate from warm-up
    S_pos = np.zeros(len(residuals))
    S_neg = np.zeros(len(residuals))
    for i in range(1, len(residuals)):
        r = residuals[i] / sigma
        S_pos[i] = max(0, S_pos[i-1] + r - k)
        S_neg[i] = max(0, S_neg[i-1] - r - k)
    alarm_idx = np.where((S_pos > h) | (S_neg > h))[0]
    first_alarm = alarm_idx[0] if len(alarm_idx) > 0 else -1
    return S_pos, S_neg, first_alarm

np.random.seed(99)
residuals = np.concatenate([
    np.random.randn(150),           # stable period
    np.random.randn(100) + 2.0,     # concept drift: bias shifts by 2 sigma
])

S_pos, S_neg, alarm = cusum_detector(residuals, k=0.5, h=5.0)
print(f"Drift detected at index: {alarm}   (true drift at index 150)")

fig, axes = plt.subplots(2, 1, figsize=(10, 5), sharex=True)
axes[0].plot(residuals, lw=0.8, label="Residuals")
axes[0].axvline(150, color="red", ls="--", label="True drift")
if alarm >= 0:
    axes[0].axvline(alarm, color="orange", ls="--", label=f"CUSUM alarm@{alarm}")
axes[0].legend(fontsize=8)
axes[1].plot(S_pos, label="CUSUM+")
axes[1].plot(S_neg, label="CUSUM-")
axes[1].axhline(5.0, color="red", ls=":", lw=0.8)
axes[1].legend(fontsize=8)
axes[1].set_xlabel("Time step")
plt.tight_layout()
plt.savefig("cusum.png", dpi=100)
plt.show()
Output
Drift detected at index: 157   (true drift at index 150)