Gradient Boosting for Time Series: Using LightGBM to Forecast

Last time, Nora built lag and rolling-window features for the bakery’s daily sales data — and watched even a simple regression model improve once it could “see” recent history. Now she wants a second opinion from a different model family. She’s stress-testing gradient boosting on a synthetic energy-consumption benchmark (a parallel test case, not the bakery’s own data) before trusting it with the sourdough forecast.

Why Gradient Boosting Works for Time Series (When Most People Get It Wrong)

Why would a model built for classification and regression work for forecasting at all? Most people assume time series needs its own toolkit — ARIMA, exponential smoothing, some specialized architecture. But gradient boosting doesn’t care about time order. It cares about patterns in features.

Say you’re predicting tomorrow’s energy consumption. A traditional model like ARIMA looks at the pattern of past values and extrapolates. Gradient boosting takes a different approach — give it features that describe the current situation, and it learns which combinations predict the next value.

Gradient boosting works by iteratively correcting its own mistakes. The first tree makes a prediction and gets it wrong. The second tree looks at those errors and tries to fix them. The third looks at what’s left and corrects again. After dozens or hundreds of trees, the accumulated corrections produce a strong forecast. The model doesn’t assume your data follows a specific mathematical pattern — it learns from examples.

Once you engineer the right features — which we’ll do in the next section — time series becomes a standard supervised learning problem. You have inputs (lag features, rolling statistics) and outputs (the next value). Boosting learns the relationship. Done.

Here’s where it gets tricky: boosting can overfit to patterns that won’t repeat. If your training data has a weird spike, the model might learn to predict spikes even when none are coming. We’ll cover how to detect this later with residual analysis.

A simple example makes this concrete:

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

# Create a simple synthetic time series: sine wave with noise
np.random.seed(42)
t = np.arange(0, 100, 0.5)  # 200 time steps
y = np.sin(t / 10) + np.random.normal(0, 0.1, len(t))  # sine wave + noise

# Create lag features: t-1, t-2, t-3
df = pd.DataFrame({'y': y})
df['lag1'] = df['y'].shift(1)
df['lag2'] = df['y'].shift(2)
df['lag3'] = df['y'].shift(3)
df = df.dropna()  # Remove rows with NaN from shifting

print(f"Dataset shape: {df.shape}")
print(f"First few rows:\n{df.head()}")

# Split by time (not randomly!)
train_size = int(0.8 * len(df))
X_train = df[['lag1', 'lag2', 'lag3']].iloc[:train_size]
y_train = df['y'].iloc[:train_size]
X_test = df[['lag1', 'lag2', 'lag3']].iloc[train_size:]
y_test = df['y'].iloc[train_size:]

# Train LightGBM
model = LGBMRegressor(num_leaves=31, learning_rate=0.05, n_estimators=100, verbose=-1)
model.fit(X_train, y_train)

# Predict and evaluate
y_pred = model.predict(X_test)
mae = mean_absolute_error(y_test, y_pred)

print(f"\nTest MAE with lag features: {mae:.4f}")
print(f"This means predictions are off by about {mae:.4f} units on average.")

# Now let's see what happens without engineered features
# (This won't work because we have no features!)
print("\nWithout lag features, boosting has nothing to learn from.")
print("Time series becomes a supervised learning problem only when you engineer features.")

• np.random.seed(42) — fixes the random number generator so the synthetic sine-wave dataset is identical every run, making results reproducible.
• np.sin(t / 10) + np.random.normal(0, 0.1, len(t)) — generates a noisy sine wave: a smooth periodic signal plus small Gaussian perturbations, mimicking a real time series with regular cyclical behavior.
• df['y'].shift(1) / shift(2) / shift(3) — creates lag features by shifting the target column back 1, 2, and 3 time steps; the first 3 rows become NaN because there’s no prior history for them.
• df.dropna() — removes the rows where lag values don’t exist yet, so the model never trains on NaN inputs.
• int(0.8 * len(df)) — computes the 80% split point for a chronological (not random) train/test split that respects time order.
• LGBMRegressor(num_leaves=31, learning_rate=0.05, n_estimators=100, verbose=-1) — creates a LightGBM regressor with 31 leaves per tree, a 0.05 learning rate, 100 boosting rounds, and silent output (no training logs).
• mean_absolute_error(y_test, y_pred) — computes the Mean Absolute Error between actual and predicted values on the held-out test set, measuring average prediction error in the same units as the data.

Run this and you’ll see a mean absolute error (MAE) around 0.11 (this run: 0.1142) on a signal that swings roughly between -1 and 1. That’s small — the model captures the sine wave pattern using nothing but its own recent history. Without lag features, there’s nothing for it to learn. Boosting needs features to work with.

Feature Engineering: The Real Secret to Boosting Time Series

Boosting’s power comes from features, not tuning. Feature engineering is unglamorous — not as flashy as neural networks — but it’s where the real work happens, and where you’ll spend most of your time on a real project. Five categories cover most of what you need:

• Lag features are the foundation. A lag-1 feature is the previous value. Lag-2 is two steps back, lag-3 three steps back. These capture recent history. If energy consumption was high yesterday, it’s likely to be high today. Boosting learns this relationship.
• Rolling statistics capture seasonality and volatility. A 7-day rolling mean smooths out daily noise and reveals the weekly trend. A 7-day rolling standard deviation shows how volatile the data is. If volatility spikes, the model can learn to be cautious.
• Trend features help boosting see direction. A simple linear regression slope over a 30-day window tells the model whether the data is trending up, down, or flat.
• Interaction features let boosting find complex patterns. A lag-1 value multiplied by the 7-day rolling mean might reveal that the model’s behavior changes depending on the baseline level.
• Domain features inject business logic. Day of week, month, holiday flags — these tell the model about external structure.

Let’s build a feature engineering pipeline:

import numpy as np
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt

# Create a more realistic time series: energy consumption with trend and seasonality
np.random.seed(42)
t = np.arange(0, 365)  # One year of daily data
trend = t * 0.01  # Slight upward trend
seasonality = 10 * np.sin(2 * np.pi * t / 365)  # Yearly pattern
daily_seasonality = 5 * np.sin(2 * np.pi * t / 7)  # Weekly pattern
noise = np.random.normal(0, 1, len(t))
y = 50 + trend + seasonality + daily_seasonality + noise

df = pd.DataFrame({
    'date': pd.date_range('2023-01-01', periods=len(t), freq='D'),
    'y': y
})

print(f"Original time series: {len(df)} days")
print(f"Mean: {df['y'].mean():.2f}, Std: {df['y'].std():.2f}")

# Lag features (1 to 7 days back)
for i in range(1, 8):
    df[f'lag_{i}'] = df['y'].shift(i)

# Rolling statistics (7-day and 30-day windows)
for window in [7, 30]:
    df[f'rolling_mean_{window}'] = df['y'].rolling(window).mean()
    df[f'rolling_std_{window}'] = df['y'].rolling(window).std()
    df[f'rolling_min_{window}'] = df['y'].rolling(window).min()
    df[f'rolling_max_{window}'] = df['y'].rolling(window).max()

# Trend feature: linear regression slope over a 30-day window
def rolling_slope(series, window=30):
    slopes = []
    for i in range(len(series)):
        if i < window:
            slopes.append(np.nan)
        else:
            x = np.arange(window)
            y_window = series.iloc[i-window:i].values
            slope, _ = np.polyfit(x, y_window, 1)
            slopes.append(slope)
    return slopes

df['trend_slope_30'] = rolling_slope(df['y'], window=30)

# Domain features
df['day_of_week'] = df['date'].dt.dayofweek
df['month'] = df['date'].dt.month
df['day_of_year'] = df['date'].dt.dayofyear

# Interaction feature: lag1 * rolling_mean_7
df['lag1_x_rolling_mean_7'] = df['lag_1'] * df['rolling_mean_7']

# Drop rows with NaN (from shifting and rolling windows)
df_clean = df.dropna()

print(f"\nAfter feature engineering: {len(df_clean)} rows, {len(df_clean.columns)} columns")
print(f"\nFeature columns:")
feature_cols = [col for col in df_clean.columns if col not in ['date', 'y']]
for col in feature_cols:
    print(f"  - {col}")

print(f"\nFirst few rows of engineered features:")
print(df_clean[['y', 'lag_1', 'lag_7', 'rolling_mean_7', 'rolling_std_7', 'trend_slope_30']].head(10))

# Correlation with target
correlations = df_clean[feature_cols + ['y']].corr()['y'].drop('y').sort_values(ascending=False)
print(f"\nTop 10 features by correlation with target:")
print(correlations.head(10))
print(f"\nBottom 5 features (weakest correlation):")
print(correlations.tail(5))

# Visualize
fig, axes = plt.subplots(3, 1, figsize=(12, 8))

axes[0].plot(df['date'], df['y'], label='Original', linewidth=2)
axes[0].set_ylabel('Value')
axes[0].set_title('Original Time Series')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(df['date'], df['y'], label='Original', alpha=0.5)
axes[1].plot(df['date'], df['rolling_mean_7'], label='7-day rolling mean', linewidth=2)
axes[1].plot(df['date'], df['rolling_mean_30'], label='30-day rolling mean', linewidth=2)
axes[1].set_ylabel('Value')
axes[1].set_title('Rolling Statistics Capture Trends')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

axes[2].bar(range(len(correlations.head(15))), correlations.head(15).values)
axes[2].set_xticks(range(len(correlations.head(15))))
axes[2].set_xticklabels(correlations.head(15).index, rotation=45, ha='right')
axes[2].set_ylabel('Correlation with Target')
axes[2].set_title('Top 15 Features by Correlation')
axes[2].grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig('feature_engineering.png', dpi=100, bbox_inches='tight')
print("\nVisualization saved as feature_engineering.png")

• 10 * np.sin(2 * np.pi * t / 365) — generates a yearly seasonal cycle with period 365 days and amplitude 10, modeling the slow seasonal swing across a full year.
• 5 * np.sin(2 * np.pi * t / 7) — generates a weekly seasonal cycle with period 7 days and amplitude 5, modeling the day-of-week pattern.
• for i in range(1, 8): df[f'lag_{i}'] = df['y'].shift(i) — creates 7 lag features (lag_1 through lag_7) in a loop, each shifted by a different number of days.
• rolling(window).mean() / .std() / .min() / .max() — computes four rolling statistics over 7-day and 30-day windows, giving the model both smoothed trends and volatility/range signals.
• rolling_slope() — a custom function that slides a 30-day window across the series, fits a degree-1 polynomial (a line) at each position, and extracts the slope — giving the model a local trend signal.
• np.polyfit(x, y_window, 1) — NumPy’s least-squares polynomial fit; 1 means degree-1 (a straight line), returning [slope, intercept].
• df['lag1_x_rolling_mean_7'] = df['lag_1'] * df['rolling_mean_7'] — an interaction feature that lets the model learn different behavior depending on whether recent values are high or low relative to the weekly average.
• df_clean[feature_cols + ['y']].corr()['y'] — computes the Pearson correlation matrix for all features plus the target, then extracts the column for y — a quick diagnostic for which features carry the most linear signal.

Run this and the correlations tell a clear story: feature engineering gives the model something real to work with. lag_7 tops the list at r ≈ 0.98 — a near-exact echo, because it lands exactly one period back on the 7-day daily_seasonality cycle, so it reproduces most of y’s structure by construction. The interaction feature lag1_x_rolling_mean_7 comes in second (r ≈ 0.91), followed by lag_1 (r ≈ 0.90) and lag_6 (r ≈ 0.90). The 7-day rolling mean sits at r ≈ 0.87 — still strong, not quite the “around 0.9” you might eyeball from the chart. At the bottom, rolling_std_7 and rolling_std_30 carry almost no linear correlation with the level of y (≈ 0.01 and −0.02) — expected, since volatility and level are different things — while month and day_of_year are moderately negatively correlated (≈ −0.70), an artifact of this particular year’s phase (the yearly sine cycle is heading down for most of the months sampled). That’s the case for feature engineering — you’re giving the model information it can actually use, though not every feature you add turns out to be useful in the way its name suggests.

Setting Up LightGBM for Time Series: Data Leakage and Train/Test Splits

The trap: a normal random train/test split breaks time order and leaks future information into the past. Say you’re forecasting energy consumption. You shuffle the data, put 80 percent in training and 20 percent in testing. December rows end up in training while January rows land in test. Your model learns patterns from the future and applies them backward. Those patterns won’t exist at deployment. Your metrics won’t reflect reality.

The fix — time-aware splitting: train on the past, test on the future. No shuffling. With data spanning January through December, you’d train on January through roughly October and test on the remainder. This respects time order and gives you honest metrics.

Walk-forward validation goes further. Train on months 1-3, test on month 4. Retrain on months 1-4, test on month 5. Retrain on months 1-5, test on month 6. You’re simulating real deployment: retrain as new data arrives, then test on the next batch. More compute, but the most realistic estimate of how your model will actually perform.

Why leakage is subtle: if you engineer features on the entire dataset before splitting, test-set information can leak into feature statistics computed with a centered or backward-and-forward window. A rolling mean computed with future context, then split, gives October 1st a value that includes October 2nd and beyond — data that wouldn’t exist in real time. Every rolling/lag feature in this article uses .shift() and .rolling(), which only look backward, so computing them on the full frame before splitting is safe here — but the discipline still matters: if you ever add a centered window or a target-derived aggregate, engineer it separately on training and test data.

Let’s implement this correctly:

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

# Use the energy data from the previous section
np.random.seed(42)
t = np.arange(0, 365)
trend = t * 0.01
seasonality = 10 * np.sin(2 * np.pi * t / 365)
daily_seasonality = 5 * np.sin(2 * np.pi * t / 7)
noise = np.random.normal(0, 1, len(t))
y = 50 + trend + seasonality + daily_seasonality + noise

df = pd.DataFrame({
    'date': pd.date_range('2023-01-01', periods=len(t), freq='D'),
    'y': y
})

print("Demonstrating the danger of random splits:\n")

# BAD: Random train/test split
from sklearn.model_selection import train_test_split
X_random_train, X_random_test, y_random_train, y_random_test = train_test_split(
    df[['y']].values, df['y'].values, test_size=0.2, random_state=42
)

print(f"Random split:")
print(f"  Training set size: {len(y_random_train)}")
print(f"  Test set size: {len(y_random_test)}")
print(f"  Training set mean: {y_random_train.mean():.2f}")
print(f"  Test set mean: {y_random_test.mean():.2f}")
print(f"  Difference: {abs(y_random_train.mean() - y_random_test.mean()):.2f}")
print(f"  -> The means are similar, which suggests data is mixed randomly.")

# GOOD: Time-aware split
train_size = int(0.8 * len(df))
X_time_train = df[['y']].iloc[:train_size].values
y_time_train = df['y'].iloc[:train_size].values
X_time_test = df[['y']].iloc[train_size:].values
y_time_test = df['y'].iloc[train_size:].values

print(f"\nTime-aware split:")
print(f"  Training set: {df['date'].iloc[0].date()} to {df['date'].iloc[train_size-1].date()}")
print(f"  Test set: {df['date'].iloc[train_size].date()} to {df['date'].iloc[-1].date()}")
print(f"  Training set mean: {y_time_train.mean():.2f}")
print(f"  Test set mean: {y_time_test.mean():.2f}")
print(f"  Difference: {abs(y_time_train.mean() - y_time_test.mean()):.2f}")

# Function for walk-forward validation
def walk_forward_validation(df, n_splits=5):
    """
    Perform walk-forward validation.
    Splits data into n_splits chunks.
    For each chunk, trains on all previous data and tests on that chunk.
    """
    chunk_size = len(df) // (n_splits + 1)  # Reserve last chunk for final test
    results = []
    
    for i in range(1, n_splits + 1):
        train_end = i * chunk_size
        test_end = (i + 1) * chunk_size
        
        train_data = df.iloc[:train_end]
        test_data = df.iloc[train_end:test_end]
        
        results.append({
            'fold': i,
            'train_start': train_data['date'].iloc[0].date(),
            'train_end': train_data['date'].iloc[-1].date(),
            'test_start': test_data['date'].iloc[0].date(),
            'test_end': test_data['date'].iloc[-1].date(),
            'train_size': len(train_data),
            'test_size': len(test_data)
        })
    
    return results

walk_forward_splits = walk_forward_validation(df, n_splits=4)
print(f"\nWalk-forward validation (4 folds):")
for split in walk_forward_splits:
    print(f"  Fold {split['fold']}: Train {split['train_start']} to {split['train_end']}, "
          f"Test {split['test_start']} to {split['test_end']}")
    print(f"           Train size: {split['train_size']}, Test size: {split['test_size']}")

print(f"\nKey takeaway: Each fold trains on all past data and tests on the next chunk.")
print(f"This simulates real deployment where you retrain as new data arrives.")

• train_test_split(..., random_state=42) — the BAD approach: train_test_split with default shuffle=True randomly permutes the rows, mixing past and future data into both the training and test sets.
• df[['y']].iloc[:train_size] — the GOOD approach: slices the first 80% of rows chronologically for training and the last 20% for testing, preserving time order.
• walk_forward_validation(df, n_splits=4) — a custom function that creates expanding-window splits: each fold trains on all data up to a growing cutoff and tests on the next non-overlapping chunk.
• chunk_size = len(df) // (n_splits + 1) — divides the dataset into equal-sized segments so each fold’s test window is a distinct, non-overlapping slice of the timeline.
• train_data = df.iloc[:train_end] — the training set for each fold includes all data from the beginning up to train_end, simulating the expanding-window approach where you retrain on all available history.
• test_data = df.iloc[train_end:test_end] — the test set for each fold is the next chunk after the training cutoff, never overlapping with training data.

Run this and you’ll see the difference — and one result that isn’t what you’d naively expect. The random split’s means are close (51.55 train vs. 52.93 test, a 1.38 difference on 292/73 rows) simply because shuffling scatters both halves of the yearly cycle into both sets. The time-aware split (train: 2023-01-01 to 2023-10-19; test: 2023-10-20 to 2023-12-31) shows a much bigger gap — 52.86 vs. 47.70, a 5.16 difference — but notice the direction: the test set mean is lower, not higher. The small upward trend (t * 0.01, worth about +3.6 over the year) is real, but the Oct–Dec test window falls right in the trough of the 365-day yearly seasonality term, and that swing (±10) dwarfs the trend. This is exactly the point of a time-aware split: it exposes the regime your model will actually face at deployment — including regimes, like a seasonal trough, that a random split would have averaged away. Walk-forward validation (4 folds, each testing on the next ~73-day chunk) is the gold standard for honest evaluation precisely because it forces the model through several such regimes instead of just one.

Building Your First LightGBM Forecaster

Here’s the end-to-end code. We load data, engineer features, split by time, train LightGBM, and make predictions.

LightGBM is fast and handles categorical features natively — unlike XGBoost. You use it like any sklearn model: fit on training data, predict on test. The hyperparameters below are sensible defaults for time series:

• num_leaves: Controls tree complexity. More leaves means more complex trees and higher overfitting risk. We start with 31.
• learning_rate: Step size for each boosting round. Smaller values (0.01–0.05) are safer; larger ones (0.1+) train faster but risk overfitting.
• n_estimators: Number of boosting rounds (trees). More trees usually improves performance, with diminishing returns.

LightGBM gives you point forecasts — a single predicted value per time step. We’ll add uncertainty bands later.

import numpy as np
import pandas as pd
from lightgbm import LGBMRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error
import matplotlib.pyplot as plt

# Create realistic time series data
np.random.seed(42)
t = np.arange(0, 365)
trend = t * 0.01
seasonality = 10 * np.sin(2 * np.pi * t / 365)
daily_seasonality = 5 * np.sin(2 * np.pi * t / 7)
noise = np.random.normal(0, 1, len(t))
y = 50 + trend + seasonality + daily_seasonality + noise

df = pd.DataFrame({
    'date': pd.date_range('2023-01-01', periods=len(t), freq='D'),
    'y': y
})

# Feature engineering
for i in range(1, 8):
    df[f'lag_{i}'] = df['y'].shift(i)

for window in [7, 30]:
    df[f'rolling_mean_{window}'] = df['y'].rolling(window).mean()
    df[f'rolling_std_{window}'] = df['y'].rolling(window).std()

df['day_of_week'] = df['date'].dt.dayofweek
df['month'] = df['date'].dt.month

df_clean = df.dropna()

# Time-aware split: 80% train, 20% test
train_size = int(0.8 * len(df_clean))

feature_cols = [col for col in df_clean.columns if col not in ['date', 'y']]

X_train = df_clean[feature_cols].iloc[:train_size]
y_train = df_clean['y'].iloc[:train_size]

X_test = df_clean[feature_cols].iloc[train_size:]
y_test = df_clean['y'].iloc[train_size:]

print(f"Training set: {len(X_train)} samples")
print(f"Test set: {len(X_test)} samples")
print(f"Features: {len(feature_cols)}")

# Train LightGBM with sensible defaults
model = LGBMRegressor(
    num_leaves=31,
    learning_rate=0.05,
    n_estimators=100,
    random_state=42,
    verbose=-1
)

model.fit(X_train, y_train)

print(f"\nModel trained successfully.")

# Make predictions
y_pred_train = model.predict(X_train)
y_pred_test = model.predict(X_test)

# Evaluate
mae_train = mean_absolute_error(y_train, y_pred_train)
mae_test = mean_absolute_error(y_test, y_pred_test)
rmse_test = np.sqrt(mean_squared_error(y_test, y_pred_test))

print(f"\nPerformance:")
print(f"  Training MAE: {mae_train:.4f}")
print(f"  Test MAE: {mae_test:.4f}")
print(f"  Test RMSE: {rmse_test:.4f}")

# Feature importance
feature_importance = pd.DataFrame({
    'feature': feature_cols,
    'importance': model.feature_importances_
}).sort_values('importance', ascending=False)

print(f"\nTop 10 most important features:")
print(feature_importance.head(10))

# Plot predictions vs actuals
fig, axes = plt.subplots(2, 1, figsize=(14, 8))

# Full time series
axes[0].plot(df_clean['date'].iloc[:train_size], y_train, label='Training data', linewidth=1.5)
axes[0].plot(df_clean['date'].iloc[train_size:], y_test, label='Test data (actual)', linewidth=1.5)
axes[0].plot(df_clean['date'].iloc[train_size:], y_pred_test, label='Test data (predicted)', 
             linewidth=1.5, linestyle='--')
axes[0].axvline(df_clean['date'].iloc[train_size], color='red', linestyle=':', alpha=0.5, label='Train/test split')
axes[0].set_ylabel('Value')
axes[0].set_title('LightGBM Time Series Forecast')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Zoomed in on test set
axes[1].plot(df_clean['date'].iloc[train_size:], y_test, label='Actual', linewidth=2, marker='o', markersize=3)
axes[1].plot(df_clean['date'].iloc[train_size:], y_pred_test, label='Predicted', linewidth=2, marker='s', markersize=3)
axes[1].fill_between(df_clean['date'].iloc[train_size:], y_test, y_pred_test, alpha=0.2, label='Error')
axes[1].set_ylabel('Value')
axes[1].set_xlabel('Date')
axes[1].set_title('Test Set: Actual vs Predicted (Zoomed)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('lightgbm_forecast.png', dpi=100, bbox_inches='tight')
print(f"\nVisualization saved as lightgbm_forecast.png")

# Feature importance plot
fig, ax = plt.subplots(figsize=(10, 6))
ax.barh(range(len(feature_importance.head(15))), feature_importance.head(15)['importance'])
ax.set_yticks(range(len(feature_importance.head(15))))
ax.set_yticklabels(feature_importance.head(15)['feature'])
ax.set_xlabel('Importance')
ax.set_title('Top 15 Feature Importances')
ax.grid(True, alpha=0.3, axis='x')
plt.tight_layout()
plt.savefig('feature_importance.png', dpi=100, bbox_inches='tight')
print(f"Feature importance plot saved as feature_importance.png")

• feature_cols = [col for col in df_clean.columns if col not in ['date', 'y']] — builds the feature list by excluding the date column (not a predictor) and the target y (what we’re predicting), so the model sees only the engineered features.
• LGBMRegressor(num_leaves=31, learning_rate=0.05, n_estimators=100, random_state=42, verbose=-1) — sensible defaults: 31 leaves per tree (moderate complexity), 0.05 learning rate (slow but safe), 100 boosting rounds, a fixed random seed for reproducibility, and silent output.
• model.feature_importances_ — LightGBM’s built-in importance scores, counting how often each feature was used in tree splits and how much it reduced error — a quick way to see which features the model actually relied on.
• fill_between(...) — shades the area between actual and predicted values on the zoomed plot, visually representing the prediction error at each time step.
• ax.axvline(df_clean['date'].iloc[train_size], ...) — draws a vertical dotted line at the train/test boundary so you can visually distinguish where training data ends and the forecast begins.
• feature_importance.head(15)['importance'] — extracts the top 15 features by importance score for the horizontal bar chart, sorted from most to least important.

Run this and you’ll get a working forecast — but a more honest read of it than “training MAE and test MAE were close.” With 268 training rows and 68 test rows across 13 features, training MAE lands at 0.6676 while test MAE comes in at 1.5900 (test RMSE 1.8863) — the test error is roughly 2.4× the training error, not “slightly higher.” That gap is the same story as the split section above: the test window (Oct–Dec) sits in a seasonal trough the model saw comparatively little of during training, and tree-based models don’t extrapolate smoothly the way a linear trend term would — they predict based on which historical leaf a new row resembles most. Feature importance also breaks from the naive intuition: lag_7 is the strongest feature by a wide margin (importance 206), not lag_1 (importance 83, ranked 4th, behind lag_6 at 109 and rolling_mean_7 at 91 too). That tracks with the correlation table from the previous section — lag_7 lands exactly one weekly cycle back, so it’s the single most informative single number the model can look up.

Tuning LightGBM: The Hyperparameters That Matter

Now that you have a working model, let’s tune it. Which knobs actually matter?

• Learning rate controls step size. Smaller values (0.01-0.05) train slower but often generalize better. Larger values (0.1+) train faster but risk overfitting. Think of it like a student: one who learns slowly and carefully tends to grasp more than one who rushes.
• Num_leaves controls tree complexity. More leaves means more capacity to fit training data — and more risk of memorizing noise. Typical range is 31 to 127. Start with 31 and increase if underfitting. (As you’ll see below, num_leaves only matters if your data is large enough to reach it — min_data_in_leaf, which defaults to 20, can cap effective complexity first.)
• Regularization (lambda_l1, lambda_l2) penalizes large tree weights. Higher values keep overfitting in check. The model pays a cost for each complex tree it builds.
• Early stopping does real work here. You monitor validation error across boosting rounds. When it stops improving, you stop training. That saves time and keeps overfitting down.

We’ll look at these as two separate exercises — a hyperparameter grid search, then early stopping — since they use different parts of the LightGBM API and answer different questions.

Grid Search: Learning Rate and Num_Leaves

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

# Recreate the dataset
np.random.seed(42)
t = np.arange(0, 365)
trend = t * 0.01
seasonality = 10 * np.sin(2 * np.pi * t / 365)
daily_seasonality = 5 * np.sin(2 * np.pi * t / 7)
noise = np.random.normal(0, 1, len(t))
y = 50 + trend + seasonality + daily_seasonality + noise

df = pd.DataFrame({
    'date': pd.date_range('2023-01-01', periods=len(t), freq='D'),
    'y': y
})

# Feature engineering
for i in range(1, 8):
    df[f'lag_{i}'] = df['y'].shift(i)
for window in [7, 30]:
    df[f'rolling_mean_{window}'] = df['y'].rolling(window).mean()
    df[f'rolling_std_{window}'] = df['y'].rolling(window).std()
df['day_of_week'] = df['date'].dt.dayofweek
df['month'] = df['date'].dt.month
df_clean = df.dropna()

# Split
train_size = int(0.8 * len(df_clean))
feature_cols = [col for col in df_clean.columns if col not in ['date', 'y']]
X_train = df_clean[feature_cols].iloc[:train_size]
y_train = df_clean['y'].iloc[:train_size]
X_test = df_clean[feature_cols].iloc[train_size:]
y_test = df_clean['y'].iloc[train_size:]

# Hyperparameter grid search
learning_rates = [0.01, 0.05, 0.1]
num_leaves_list = [31, 63, 127]

results = []

for lr in learning_rates:
    for num_leaves in num_leaves_list:
        model = LGBMRegressor(
            learning_rate=lr,
            num_leaves=num_leaves,
            n_estimators=200,
            random_state=42,
            verbose=-1
        )
        model.fit(X_train, y_train)
        y_pred = model.predict(X_test)
        mae = mean_absolute_error(y_test, y_pred)
        results.append({
            'learning_rate': lr,
            'num_leaves': num_leaves,
            'mae': mae
        })
        print(f"lr={lr}, leaves={num_leaves}: MAE={mae:.4f}")

# Find best combination
results_df = pd.DataFrame(results)
best = results_df.loc[results_df['mae'].idxmin()]
print(f"\nBest combination:")
print(f"  Learning rate: {best['learning_rate']}")
print(f"  Num leaves: {int(best['num_leaves'])}")
print(f"  Test MAE: {best['mae']:.4f}")

# Visualize grid search results
fig, ax = plt.subplots(figsize=(10, 6))
pivot = results_df.pivot(index='num_leaves', columns='learning_rate', values='mae')
im = ax.imshow(pivot.values, cmap='RdYlGn_r', aspect='auto')
ax.set_xticks(range(len(pivot.columns)))
ax.set_yticks(range(len(pivot.index)))
ax.set_xticklabels(pivot.columns)
ax.set_yticklabels(pivot.index)
ax.set_xlabel('Learning Rate')
ax.set_ylabel('Num Leaves')
ax.set_title('Hyperparameter Grid Search: Test MAE')

# Add values to cells
for i in range(len(pivot.index)):
    for j in range(len(pivot.columns)):
        text = ax.text(j, i, f'{pivot.values[i, j]:.3f}',
                      ha="center", va="center", color="black", fontsize=10)

plt.colorbar(im, ax=ax, label='MAE')
plt.tight_layout()
plt.savefig('hyperparameter_grid.png', dpi=100, bbox_inches='tight')
print(f"\nGrid search visualization saved as hyperparameter_grid.png")

• learning_rates = [0.01, 0.05, 0.1] and num_leaves_list = [31, 63, 127] — a 3×3 grid exploring how step size and tree complexity interact to affect test MAE.
• results_df.loc[results_df['mae'].idxmin()] — finds the row with the lowest MAE value, identifying the best hyperparameter combination from the grid.
• results_df.pivot(index='num_leaves', columns='learning_rate', values='mae') — reshapes the results DataFrame into a matrix (rows = num_leaves, columns = learning_rate) suitable for a heatmap.

Run this and the grid tells a sharper story than “one combination wins by a small margin.” All nine MAE values collapse into exactly three, one per learning rate — 1.1707 for every num_leaves at lr=0.01, 1.5072 at lr=0.05, and 1.5426 at lr=0.1. num_leaves had zero measurable effect here: with only 268 training rows and the default min_data_in_leaf=20, a tree can’t grow past roughly 268 ÷ 20 ≈ 13 leaves regardless of the num_leaves ceiling you set, so 31, 63, and 127 all hit the same practical wall. The learning rate, on the other hand, mattered a lot: the gap between the best (lr=0.01, MAE 1.1707) and worst (lr=0.1, MAE 1.5426) combination is about 24%, not the 5-10% you might expect from casual tuning. The lesson isn’t “tuning barely helps” — it’s that which knob helps depends on your data size, and testing that assumption (rather than grid-searching blind) is worth doing before you spend compute on a wide num_leaves sweep.

Early Stopping: Halting Before You Overfit

Grid search answers “which fixed configuration works best?” Early stopping answers a different question: for one configuration, how many boosting rounds should you actually run? Instead of comparing several fully-trained models, you watch validation error across a single long training run and stop once it stalls. That needs LightGBM’s native train() API (rather than the sklearn-style LGBMRegressor) so you can capture the per-round validation history. Picking back up with the same X_train/y_train/X_test/y_test split from the grid search above:

import numpy as np
import pandas as pd
from lightgbm import train, Dataset, record_evaluation
import matplotlib.pyplot as plt

# Recreate the same dataset and split used for the grid search
np.random.seed(42)
t = np.arange(0, 365)
trend = t * 0.01
seasonality = 10 * np.sin(2 * np.pi * t / 365)
daily_seasonality = 5 * np.sin(2 * np.pi * t / 7)
noise = np.random.normal(0, 1, len(t))
y = 50 + trend + seasonality + daily_seasonality + noise

df = pd.DataFrame({
    'date': pd.date_range('2023-01-01', periods=len(t), freq='D'),
    'y': y
})
for i in range(1, 8):
    df[f'lag_{i}'] = df['y'].shift(i)
for window in [7, 30]:
    df[f'rolling_mean_{window}'] = df['y'].rolling(window).mean()
    df[f'rolling_std_{window}'] = df['y'].rolling(window).std()
df['day_of_week'] = df['date'].dt.dayofweek
df['month'] = df['date'].dt.month
df_clean = df.dropna()

train_size = int(0.8 * len(df_clean))
feature_cols = [col for col in df_clean.columns if col not in ['date', 'y']]
X_train = df_clean[feature_cols].iloc[:train_size]
y_train = df_clean['y'].iloc[:train_size]
X_test = df_clean[feature_cols].iloc[train_size:]
y_test = df_clean['y'].iloc[train_size:]

train_data = Dataset(X_train, label=y_train)
valid_data = Dataset(X_test, label=y_test, reference=train_data)

params = {
    'objective': 'regression',
    'metric': 'mae',
    'learning_rate': 0.05,
    'num_leaves': 31,
    'verbose': -1
}

evals_result = {}
model_early = train(
    params,
    train_data,
    num_boost_round=300,
    valid_sets=[train_data, valid_data],
    valid_names=['train', 'valid'],
    callbacks=[record_evaluation(evals_result)]
)

print(f"\nTraining completed.")
print(f"Number of boosting rounds: {model_early.num_trees()}")

# Plot training progress
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(evals_result['train']['l1'], label='Training MAE', linewidth=2)
ax.plot(evals_result['valid']['l1'], label='Validation MAE', linewidth=2)
ax.set_xlabel('Boosting Round')
ax.set_ylabel('MAE')
ax.set_title('Training Progress: Early Stopping Concept')
ax.legend()
ax.grid(True, alpha=0.3)

# Mark where validation stopped improving
valid_mae = evals_result['valid']['l1']
best_round = int(np.argmin(valid_mae))
ax.axvline(best_round, color='red', linestyle='--', alpha=0.5, label=f'Best round: {best_round}')
ax.legend()

plt.tight_layout()
plt.savefig('early_stopping.png', dpi=100, bbox_inches='tight')
print(f"Early stopping plot saved as early_stopping.png")
print(f"Best round: {best_round}")
print(f"Final training MAE (round 300): {evals_result['train']['l1'][-1]:.4f}")
print(f"Final validation MAE (round 300): {evals_result['valid']['l1'][-1]:.4f}")
print(f"Best validation MAE (round {best_round}): {min(valid_mae):.4f}")

• from lightgbm import train, Dataset, record_evaluation — switches from the sklearn-compatible LGBMRegressor API to LightGBM’s native API, which provides access to evaluation logging via a callback. (LightGBM 4.x dropped the old evals_result= keyword argument on train() — you now pass record_evaluation(evals_result) as a callback instead, or train() raises TypeError: train() got an unexpected keyword argument 'evals_result'.)
• Dataset(X_train, label=y_train) — wraps the training data in LightGBM’s internal data structure, which is more memory-efficient than passing raw arrays.
• valid_sets=[train_data, valid_data] — passes both training and validation datasets so LightGBM logs MAE on both at each boosting round, enabling you to spot overfitting.
• callbacks=[record_evaluation(evals_result)] — registers a callback that fills the initially empty evals_result dict with per-round training and validation metrics as training runs, later used to plot the learning curve.
• np.argmin(valid_mae) — finds the boosting round where validation MAE was lowest, marking the “best” stopping point on the training-progress plot.

Run this and you’ll see exactly why early stopping earns its keep. Training runs the full 300 rounds regardless (this demonstration doesn’t halt training — it just records where the would-be stopping point was), and training MAE keeps falling the whole way, ending at 0.3375. Validation MAE doesn’t cooperate: it bottoms out at 1.1704 around round 36, then climbs back up to 1.4891 by round 300 — training all the way to 300 rounds leaves you about 27% worse on held-out data than stopping at the right moment would have. That’s the real argument for early stopping: not “it trains faster,” but “training longer actively makes the held-out forecast worse” once the model starts fitting training-set noise instead of signal.

Handling Seasonality and Trends: When Simple Lags Aren’t Enough

Real time series have repeating patterns (seasonality) and long-term direction (trend). Lag features alone can miss these.

Seasonality is a pattern that repeats at fixed intervals. Energy consumption shows strong daily seasonality — peaks in the evening, lows at night — and weekly seasonality, where weekends differ from weekdays. Stock prices tend toward weak seasonality.

Trend is the long-term direction, separate from noise. Energy consumption might trend up in winter (more heating) and down in summer.

Lag features capture recent history but can miss seasonality if the lag window is too short. Use only lag-1, lag-2, and lag-3, and the model sees recent values without any sense of “what happened at this time last week.” That’s why seasonal lags matter in principle — lag-7 (for daily data) captures weekly patterns; lag-365 would capture yearly patterns. Here’s how to engineer them, and — just as important — how to check whether adding them actually helped:

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

# Create data with strong seasonality
np.random.seed(42)
t = np.arange(0, 365)
trend = t * 0.01
seasonality = 10 * np.sin(2 * np.pi * t / 365)  # Yearly
daily_seasonality = 5 * np.sin(2 * np.pi * t / 7)  # Weekly
noise = np.random.normal(0, 1, len(t))
y = 50 + trend + seasonality + daily_seasonality + noise

df = pd.DataFrame({
    'date': pd.date_range('2023-01-01', periods=len(t), freq='D'),
    'y': y
})

# Basic lag features
for i in range(1, 8):
    df[f'lag_{i}'] = df['y'].shift(i)

# Rolling statistics
for window in [7, 30]:
    df[f'rolling_mean_{window}'] = df['y'].rolling(window).mean()
    df[f'rolling_std_{window}'] = df['y'].rolling(window).std()

# Seasonal lags (lag_7 duplicates the column the loop above already created;
# lag_14 and lag_21 are the genuinely new seasonal lags)
df['lag_7'] = df['y'].shift(7)    # One week ago
df['lag_14'] = df['y'].shift(14)  # Two weeks ago
df['lag_21'] = df['y'].shift(21)  # Three weeks ago

# Trend feature
def rolling_slope(series, window=30):
    slopes = []
    for i in range(len(series)):
        if i < window:
            slopes.append(np.nan)
        else:
            x = np.arange(window)
            y_window = series.iloc[i-window:i].values
            slope, _ = np.polyfit(x, y_window, 1)
            slopes.append(slope)
    return slopes

df['trend_slope_30'] = rolling_slope(df['y'], window=30)

# Domain features
df['day_of_week'] = df['date'].dt.dayofweek
df['month'] = df['date'].dt.month
df['day_of_year'] = df['date'].dt.dayofyear

df_clean = df.dropna()

# Split
train_size = int(0.8 * len(df_clean))
feature_cols = [col for col in df_clean.columns if col not in ['date', 'y']]

X_train = df_clean[feature_cols].iloc[:train_size]
y_train = df_clean['y'].iloc[:train_size]
X_test = df_clean[feature_cols].iloc[train_size:]
y_test = df_clean['y'].iloc[train_size:]

# Train WITHOUT the extra seasonal/trend features
feature_cols_basic = [col for col in feature_cols if col not in ['lag_7', 'lag_14', 'lag_21', 'trend_slope_30']]
model_basic = LGBMRegressor(learning_rate=0.05, num_leaves=31, n_estimators=100, random_state=42, verbose=-1)
model_basic.fit(X_train[feature_cols_basic], y_train)
y_pred_basic = model_basic.predict(X_test[feature_cols_basic])
mae_basic = mean_absolute_error(y_test, y_pred_basic)

print(f"Model WITHOUT seasonal/trend features:")
print(f"  Features: {len(feature_cols_basic)}")
print(f"  Test MAE: {mae_basic:.4f}")

# Train WITH the extra seasonal/trend features
model_seasonal = LGBMRegressor(learning_rate=0.05, num_leaves=31, n_estimators=100, random_state=42, verbose=-1)
model_seasonal.fit(X_train[feature_cols], y_train)
y_pred_seasonal = model_seasonal.predict(X_test[feature_cols])
mae_seasonal = mean_absolute_error(y_test, y_pred_seasonal)

print(f"\nModel WITH seasonal/trend features:")
print(f"  Features: {len(feature_cols)}")
print(f"  Test MAE: {mae_seasonal:.4f}")

change = (mae_basic - mae_seasonal) / mae_basic * 100
print(f"\nChange: {change:.1f}%")

# Feature importance comparison
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Basic model
fi_basic = pd.DataFrame({
    'feature': feature_cols_basic,
    'importance': model_basic.feature_importances_
}).sort_values('importance', ascending=False).head(10)
axes[0].barh(range(len(fi_basic)), fi_basic['importance'])
axes[0].set_yticks(range(len(fi_basic)))
axes[0].set_yticklabels(fi_basic['feature'])
axes[0].set_xlabel('Importance')
axes[0].set_title(f'Basic Model (MAE: {mae_basic:.4f})')
axes[0].grid(True, alpha=0.3, axis='x')

# Seasonal model
fi_seasonal = pd.DataFrame({
    'feature': feature_cols,
    'importance': model_seasonal.feature_importances_
}).sort_values('importance', ascending=False).head(10)
axes[1].barh(range(len(fi_seasonal)), fi_seasonal['importance'])
axes[1].set_yticks(range(len(fi_seasonal)))
axes[1].set_yticklabels(fi_seasonal['feature'])
axes[1].set_xlabel('Importance')
axes[1].set_title(f'Seasonal Model (MAE: {mae_seasonal:.4f})')
axes[1].grid(True, alpha=0.3, axis='x')

plt.tight_layout()
plt.savefig('seasonal_features_importance.png', dpi=100, bbox_inches='tight')
print(f"\nVisualization saved as seasonal_features_importance.png")

# Plot predictions
fig, ax = plt.subplots(figsize=(14, 6))
ax.plot(df_clean['date'].iloc[train_size:], y_test, label='Actual', linewidth=2, marker='o', markersize=3)
ax.plot(df_clean['date'].iloc[train_size:], y_pred_basic, label='Basic model', linewidth=1.5, linestyle='--', alpha=0.7)
ax.plot(df_clean['date'].iloc[train_size:], y_pred_seasonal, label='Seasonal model', linewidth=1.5, linestyle='--', alpha=0.7)
ax.set_ylabel('Value')
ax.set_xlabel('Date')
ax.set_title('Predictions: Basic vs Seasonal Features')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('seasonal_predictions.png', dpi=100, bbox_inches='tight')
print(f"Predictions plot saved as seasonal_predictions.png")

• df['lag_7'] = df['y'].shift(7) — re-creates the same lag_7 column the earlier loop already built (it’s listed here for clarity, not because it’s new); lag_14 and lag_21 are the two genuinely new seasonal lags at two- and three-week offsets.
• feature_cols_basic = [col for col in feature_cols if col not in ['lag_7', 'lag_14', 'lag_21', 'trend_slope_30']] — removes lag_7 along with the newer seasonal/trend features to create a “basic” set, so the comparison below is really “6 recent lags + rolling stats + calendar” vs. “the same, plus lag_7, lag_14, lag_21, and the 30-day trend slope.”
• (mae_basic - mae_seasonal) / mae_basic * 100 — computes the percentage change in MAE from adding the extra features; a negative value means the extra features made held-out error worse, not better.
• The code trains two identical, unregularized LightGBM models (default num_leaves=31, no lambda_l1/lambda_l2) — one on the basic feature set, one on the full set — so the only difference between them is those four extra, partially-correlated features.

Run this and the result cuts against the intuition the earlier correlation table built. The basic model (13 features, no lag_7/lag_14/lag_21/trend_slope_30) scores MAE 1.1586. Adding the seasonal lags and trend slope (17 features) makes it worse — MAE 1.7846, a 54% increase in error, not the improvement you’d expect from adding lag_7 (the single most-correlated feature in the dataset, r ≈ 0.98). What’s going on: these are the same unregularized settings from the “Building Your First LightGBM Forecaster” section (num_leaves=31, no lambda_l1/lambda_l2), trained on the same 268-row window and evaluated on the same seasonal-trough test period. Four extra, partly redundant features (three lags that are highly correlated with each other and with existing rolling stats, plus a noisy 30-day slope) give an unregularized tree more ways to fit patterns specific to the training window that don’t hold in the test window — the model overfits harder, not less. This is the overfitting problem the next two sections tackle directly: more features raise the ceiling on what a model can learn, but without regularization or a validation-driven stopping point, that ceiling includes memorizing noise. Feature engineering and regularization aren’t separate concerns — they’re a package deal.

Multi-Step Forecasting: Predicting Multiple Steps Ahead

So far we’ve predicted one step ahead. Given today, you predict tomorrow. Often, though, you need to forecast a week or month ahead.

There are two main strategies: direct and recursive.

Direct multi-step trains separate models for each horizon. One model predicts t+1, another predicts t+2, another t+3. Each model uses the same input features—lags from t, t-1, and so on. Since each model is independent, errors don’t accumulate step-to-step. The catch is you have to train multiple models.

Recursive multi-step trains one model for t+1. To predict t+2, you pass the predicted t+1 back in as a feature. For t+3, you use both predicted t+1 and t+2. Errors can accumulate here: if the t+1 prediction is off, it affects t+2 and beyond. But you only train one model.

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

# Create dataset
np.random.seed(42)
t = np.arange(0, 365)
trend = t * 0.01
seasonality = 10 * np.sin(2 * np.pi * t / 365)
daily_seasonality = 5 * np.sin(2 * np.pi * t / 7)
noise = np.random.normal(0, 1, len(t))
y = 50 + trend + seasonality + daily_seasonality + noise

df = pd.DataFrame({
    'date': pd.date_range('2023-01-01', periods=len(t), freq='D'),
    'y': y
})

# Feature engineering
for i in range(1, 8):
    df[f'lag_{i}'] = df['y'].shift(i)
for window in [7, 30]:
    df[f'rolling_mean_{window}'] = df['y'].rolling(window).mean()
    df[f'rolling_std_{window}'] = df['y'].rolling(window).std()
df['day_of_week'] = df['date'].dt.dayofweek
df['month'] = df['date'].dt.month
df_clean = df.dropna()

# Split
train_size = int(0.8 * len(df_clean))
feature_cols = [col for col in df_clean.columns if col not in ['date', 'y']]
X_train = df_clean[feature_cols].iloc[:train_size]
y_train = df_clean['y'].iloc[:train_size]
X_test = df_clean[feature_cols].iloc[train_size:]
y_test = df_clean['y'].iloc[train_size:]

print("DIRECT MULTI-STEP FORECASTING\n")

# Direct multi-step: train separate models for each horizon
horizon = 7  # Forecast 7 days ahead
models_direct = {}
y_pred_direct = np.zeros((len(X_test), horizon))

for h in range(1, horizon + 1):
    # Create target for h steps ahead
    y_train_h = df_clean['y'].iloc[:train_size + h].shift(-h).iloc[:train_size]
    
    # Train model for this horizon
    model = LGBMRegressor(learning_rate=0.05, num_leaves=31, n_estimators=100, random_state=42, verbose=-1)
    model.fit(X_train, y_train_h)
    models_direct[h] = model
    
    # Predict
    y_pred_direct[:, h-1] = model.predict(X_test)
    
    # Evaluate
    y_test_h = df_clean['y'].iloc[train_size:].shift(-h).iloc[:len(X_test)]
    mae_h = mean_absolute_error(y_test_h.dropna(), y_pred_direct[:len(y_test_h.dropna()), h-1])
    print(f"Horizon {h}: MAE = {mae_h:.4f}")

print(f"\nRECURSIVE MULTI-STEP FORECASTING\n")

# Recursive multi-step: use predictions as features
y_pred_recursive = np.zeros((len(X_test), horizon))

# Train one model for t+1
model_recursive = LGBMRegressor(learning_rate=0.05, num_leaves=31, n_estimators=100, random_state=42, verbose=-1)
model_recursive.fit(X_train, y_train)

# Predict step by step
for i in range(len(X_test)):
    X_current = X_test.iloc[i:i+1].copy()
    
    for h in range(horizon):
        # Predict next step
        pred = model_recursive.predict(X_current)[0]
        y_pred_recursive[i, h] = pred
        
        # Update features for next step: shift lags and add new prediction
        if h < horizon - 1:  # Don't update on last step
            X_current['lag_1'] = pred
            X_current['lag_2'] = X_current['lag_1'].values[0]
            X_current['lag_3'] = X_current['lag_2'].values[0]
            # Rolling stats would need more history, so we skip updating them

# Evaluate recursive
for h in range(1, horizon + 1):
    y_test_h = df_clean['y'].iloc[train_size:].shift(-h).iloc[:len(X_test)]
    mae_h = mean_absolute_error(y_test_h.dropna(), y_pred_recursive[:len(y_test_h.dropna()), h-1])
    print(f"Horizon {h}: MAE = {mae_h:.4f}")

print(f"\nComparison:")
print(f"Direct multi-step: Train {horizon} models, no step-to-step error accumulation")
print(f"Recursive multi-step: Train 1 model, but only lag_1-lag_3 get refreshed each step")

# Plot
fig, axes = plt.subplots(2, 1, figsize=(14, 8))

# Direct
for h in range(horizon):
    axes[0].plot(range(len(X_test)), y_pred_direct[:, h], marker='o', markersize=2, label=f'Horizon {h+1}')
axes[0].set_ylabel('Predicted Value')
axes[0].set_title('Direct Multi-Step: Separate Models for Each Horizon')
axes[0].legend(ncol=7)
axes[0].grid(True, alpha=0.3)

# Recursive
for h in range(horizon):
    axes[1].plot(range(len(X_test)), y_pred_recursive[:, h], marker='s', markersize=2, label=f'Horizon {h+1}')
axes[1].set_ylabel('Predicted Value')
axes[1].set_xlabel('Test Sample')
axes[1].set_title('Recursive Multi-Step: One Model, Predictions Feed Into Next Step')
axes[1].legend(ncol=7)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('multistep_forecasting.png', dpi=100, bbox_inches='tight')
print(f"\nVisualization saved as multistep_forecasting.png")

• for h in range(1, horizon + 1) — the direct approach: trains a separate LightGBM model for each forecast horizon (h=1 through h=7), each independently predicting h steps ahead.
• df_clean['y'].iloc[:train_size + h].shift(-h).iloc[:train_size] — shifts the target forward by h steps to create the h-step-ahead target, so each model learns to predict a different horizon using the same lag features from time t.
• model_recursive.predict(X_current)[0] — the recursive approach: predicts one step ahead, then feeds that prediction back into the feature vector for the next step.
• X_current['lag_1'] = pred — overwrites the lag_1 feature with the model’s own one-step-ahead prediction, creating the recursive chain where each prediction depends on prior predictions.
• X_current['lag_2'] = X_current['lag_1'].values[0] — shifts the old lag_1 into lag_2 and lag_2 into lag_3, maintaining the lag structure as the recursive window moves forward.
• The comment “Rolling stats would need more history, so we skip updating them” — acknowledges a simplification: rolling-window features (e.g., rolling_mean_7) and lags beyond lag_3 stay frozen at each test row’s real historical values throughout the whole 7-step chain — only lag_1, lag_2, and lag_3 actually get overwritten with predictions.

Run this and the two curves don’t tell quite the story you’d expect from “recursive errors compound over the horizon.” Direct’s MAE does climb fairly steadily as the horizon grows — 1.8817 (h=1) up to 3.5105 (h=7) — because farther-out targets are inherently harder to predict from today’s features, independent model or not. Recursive’s MAE, though, is not monotonic: it starts at 2.8367 (h=1), climbs sharply to a peak of 6.0829 at h=3, then falls back down to 2.4823 by h=7 — actually beating direct’s h=7 result (3.5105). The reason is visible in the code: only lag_1 through lag_3 get overwritten with predictions on each recursive step; lag_4 through lag_7, both rolling windows, and the calendar features all stay pinned to each test row’s real historical values for the entire 7-step chain. So past h=3, the majority of the recursive model’s inputs are still genuine data, not compounded guesses — which caps how badly it can drift, but also means it isn’t really forecasting 5-7 steps ahead using its own trajectory the way a “true” recursive forecaster would. Direct wins clearly at horizons 1 through 6; only at horizon 7 does recursive pull ahead, and for a reason (frozen features, not superior long-range modeling) that argues for caution rather than confidence in that result.

Avoiding Overfitting: Validation, Regularization, and Reality Checks

Boosting is greedy. It will memorize noise if you let it. A model that performs perfectly on training data but poorly on test data has overfit.

Overfitting happens when a model learns noise instead of patterns. Consider a model that learns “on Tuesdays in October, the value is 47.3.” It memorized the training data but won’t generalize to new Tuesdays.

How do you spot it? Plot training error and validation error over boosting rounds. If they diverge — validation error keeps rising while training error falls — you’re overfitting.

The usual fixes:

• lambda_l1 / lambda_l2 penalize large tree weights.
• max_depth limits tree depth.
• min_data_in_leaf requires a minimum number of samples in each leaf (as the grid search above showed, this can silently cap num_leaves too).
• Early stopping halts training when validation error plateaus.

import numpy as np
import pandas as pd
from lightgbm import train, Dataset, record_evaluation
from sklearn.metrics import mean_absolute_error
import matplotlib.pyplot as plt

# Create dataset
np.random.seed(42)
t = np.arange(0, 365)
trend = t * 0.01
seasonality = 10 * np.sin(2 * np.pi * t / 365)
daily_seasonality = 5 * np.sin(2 * np.pi * t / 7)
noise = np.random.normal(0, 1, len(t))
y = 50 + trend + seasonality + daily_seasonality + noise

df = pd.DataFrame({
    'date': pd.date_range('2023-01-01', periods=len(t), freq='D'),
    'y': y
})

# Feature engineering
for i in range(1, 8):
    df[f'lag_{i}'] = df['y'].shift(i)
for window in [7, 30]:
    df[f'rolling_mean_{window}'] = df['y'].rolling(window).mean()
    df[f'rolling_std_{window}'] = df['y'].rolling(window).std()
df['day_of_week'] = df['date'].dt.dayofweek
df['month'] = df['date'].dt.month
df_clean = df.dropna()

# Split
train_size = int(0.8 * len(df_clean))
feature_cols = [col for col in df_clean.columns if col not in ['date', 'y']]
X_train = df_clean[feature_cols].iloc[:train_size]
y_train = df_clean['y'].iloc[:train_size]
X_test = df_clean[feature_cols].iloc[train_size:]
y_test = df_clean['y'].iloc[train_size:]

print("Training with and without regularization:\n")

# Train WITHOUT regularization
train_data = Dataset(X_train, label=y_train)
valid_data = Dataset(X_test, label=y_test, reference=train_data)

params_no_reg = {
    'objective': 'regression',
    'metric': 'mae',
    'learning_rate': 0.05,
    'num_leaves': 127,  # High complexity
    'lambda_l1': 0,     # No L1 regularization
    'lambda_l2': 0,     # No L2 regularization
    'verbose': -1
}

evals_result_no_reg = {}
model_no_reg = train(
    params_no_reg,
    train_data,
    num_boost_round=200,
    valid_sets=[train_data, valid_data],
    valid_names=['train', 'valid'],
    callbacks=[record_evaluation(evals_result_no_reg)]
)

# Train WITH regularization
params_reg = {
    'objective': 'regression',
    'metric': 'mae',
    'learning_rate': 0.05,
    'num_leaves': 31,   # Lower complexity
    'lambda_l1': 1.0,   # L1 regularization
    'lambda_l2': 1.0,   # L2 regularization
    'min_data_in_leaf': 20,
    'verbose': -1
}

evals_result_reg = {}
model_reg = train(
    params_reg,
    train_data,
    num_boost_round=200,
    valid_sets=[train_data, valid_data],
    valid_names=['train', 'valid'],
    callbacks=[record_evaluation(evals_result_reg)]
)

print(f"Model WITHOUT regularization:")
print(f"  Final training MAE: {evals_result_no_reg['train']['l1'][-1]:.4f}")
print(f"  Final validation MAE: {evals_result_no_reg['valid']['l1'][-1]:.4f}")
print(f"  Divergence: {evals_result_no_reg['valid']['l1'][-1] - evals_result_no_reg['train']['l1'][-1]:.4f}")

print(f"\nModel WITH regularization:")
print(f"  Final training MAE: {evals_result_reg['train']['l1'][-1]:.4f}")
print(f"  Final validation MAE: {evals_result_reg['valid']['l1'][-1]:.4f}")
print(f"  Divergence: {evals_result_reg['valid']['l1'][-1] - evals_result_reg['train']['l1'][-1]:.4f}")

# Plot
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# No regularization
axes[0].plot(evals_result_no_reg['train']['l1'], label='Training', linewidth=2)
axes[0].plot(evals_result_no_reg['valid']['l1'], label='Validation', linewidth=2)
axes[0].set_xlabel('Boosting Round')
axes[0].set_ylabel('MAE')
axes[0].set_title('WITHOUT Regularization (Overfitting)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# With regularization
axes[1].plot(evals_result_reg['train']['l1'], label='Training', linewidth=2)
axes[1].plot(evals_result_reg['valid']['l1'], label='Validation', linewidth=2)
axes[1].set_xlabel('Boosting Round')
axes[1].set_ylabel('MAE')
axes[1].set_title('WITH Regularization')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('regularization_effect.png', dpi=100, bbox_inches='tight')
print(f"\nVisualization saved as regularization_effect.png")

# Sanity check: ensure predictions are reasonable
y_pred_no_reg = model_no_reg.predict(X_test)
y_pred_reg = model_reg.predict(X_test)

print(f"\nSanity checks:")
print(f"Training data range: [{y_train.min():.2f}, {y_train.max():.2f}]")
print(f"Test data range: [{y_test.min():.2f}, {y_test.max():.2f}]")
print(f"\nPredictions (no reg) range: [{y_pred_no_reg.min():.2f}, {y_pred_no_reg.max():.2f}]")
print(f"Predictions (with reg) range: [{y_pred_reg.min():.2f}, {y_pred_reg.max():.2f}]")

• params_no_reg — a parameter dict with num_leaves=127 (high complexity), lambda_l1=0, and lambda_l2=0 — intentionally no regularization, so the model is free to memorize training data.
• params_reg — a parameter dict with num_leaves=31 (lower complexity), lambda_l1=1.0, lambda_l2=1.0, and min_data_in_leaf=20 — a regularized configuration that penalizes complexity and requires a minimum number of samples per leaf.
• lambda_l1 and lambda_l2 — L1 (Lasso-style) and L2 (Ridge-style) regularization penalties applied to leaf weights; higher values discourage the model from producing extreme predictions driven by a few noisy samples.
• min_data_in_leaf=20 — requires at least 20 data points in each leaf node, preventing the tree from creating overly specific rules that split on tiny subsets of data — a direct defense against memorizing individual data points.
• callbacks=[record_evaluation(evals_result_no_reg)] / callbacks=[record_evaluation(evals_result_reg)] — the LightGBM 4.x callback that fills each dict with per-round training and validation MAE, used to plot whether the two curves diverge (overfitting) or stay close (good generalization).
• model_no_reg.predict(X_test) vs model_reg.predict(X_test) — sanity checks: both models’ predicted ranges should fall within the training data’s range; predictions far outside that range would signal a problem.

Run this and the effect is real, but modest — worth reporting exactly rather than rounding up to “regularization fixes overfitting.” Without regularization: final training MAE 0.4461, final validation MAE 1.5072, a divergence of 1.0611. With regularization: final training MAE 0.5049 (the model fits the training data slightly worse, as intended), final validation MAE 1.5273, divergence 1.0224. The gap narrows by about 4% — but validation MAE itself doesn’t improve; it’s actually a touch higher with regularization (1.5273 vs. 1.5072) in this single train/test split. Both prediction ranges stay sane (no-reg: [38.77, 56.67]; reg: [38.59, 56.09], both inside the training range of [36.21, 66.81]), so neither model is doing anything pathological — but this one fixed pair of lambda_l1=1.0/lambda_l2=1.0 values isn’t automatically the right amount of regularization for this data. The honest takeaway: regularization is a knob to search over against validation performance (the way the grid search section searched over learning_rate), not a single value you set once and trust to help.

Interpreting Results: Feature Importance and Residual Analysis

A model is only useful if you understand what it learned. Feature importance shows which inputs the model leaned on. Residual analysis shows whether the errors are random or systematic.

LightGBM reports feature importance directly: how many times each feature was used in splits, and how much it reduced error. High importance means the model leaned on that feature.

Residuals are just actual minus predicted. For a good model, they should look like random noise — no pattern, centered at zero, roughly normal.

If residuals show patterns (say, they’re high when the actual value is high), the model missed something.

Gradient boosting vs. ARIMA vs. Prophet: which should you reach for?

Approach	Strengths	Weaknesses	Best When
Gradient boosting (LightGBM)	Handles non-linear patterns and feature interactions natively; no distributional assumptions; fast training; handles external regressors (weather, promotions) easily	Requires manual feature engineering (lags, rolling stats, calendar); can overfit to noise in small datasets; point forecasts only unless you add quantile regression	You have rich features, stable seasonal patterns, and a medium-to-large dataset; you want to combine time-based features with external drivers
ARIMA / SARIMA	Discovers lag structure automatically — no manual feature engineering; interpretable parameters; well-understood confidence intervals; works well on small datasets	Assumes linear relationships; struggles with multiple overlapping seasonalities; breaks on structural changes; limited support for external regressors	You have limited data, a single seasonal cycle, and want interpretable statistical guarantees
Prophet	Handles multiple seasonalities, holidays, and structural breaks out of the box; tunable trend flexibility; produces uncertainty intervals automatically	Less control over feature interactions; can over-smooth sharp changes; slower on very large datasets; treats external regressors additively rather than through tree interactions	Your data has multiple overlapping seasonalities, holidays, or regime shifts (e.g., a new store opening) that linear models can’t handle

The key decision: gradient boosting wins when you can engineer good features and the underlying pattern is stable — it excels at learning non-linear interactions between lags, rolling statistics, and external drivers. ARIMA wins when data is limited and the pattern is simple and linear. Prophet wins when the pattern itself shifts over time — holidays, structural breaks, or multiple overlapping seasonalities that make both ARIMA’s linear assumptions and boosting’s reliance on stable patterns a liability.

Gradient boosting handles stable seasonal patterns well — though, as this article’s own numbers show, “stable” is doing real work in that sentence: even a clean synthetic signal produced a test window (a seasonal trough) where errors were notably higher than training, and extra features that looked correlated on paper made things worse without regularization to back them up. Nora’s real headache is worse still — a store whose sales shifted overnight when a new branch opened nearby — and that kind of structural break is where Prophet earns its keep.

Check Your Understanding

Questions below move from simple recall up to open-ended design, roughly following Bloom’s Taxonomy.

Remember Why does gradient boosting need engineered features (like lags and rolling statistics) to work on time series, when a model like ARIMA doesn’t need that feature engineering step?

Understand In your own words, explain why a random train/test split “leaks future information into the past” for time series data, using the article’s December/January example.

Apply Using the article’s Direct vs. Recursive multi-step framework, if you needed to forecast exactly 3 days ahead for a dashboard that refreshes hourly, which approach would you choose and why, based on the article’s stated tradeoffs (accuracy vs. number of models to maintain)?

Analyze The article’s regularization comparison shows the no-regularization model’s training and validation MAE diverging (1.0611) while the regularized model’s stay marginally closer (1.0224) — a modest effect, not a dramatic one. Walk through why lambda_l1, lambda_l2, and min_data_in_leaf specifically target the mechanism that lets boosting “memorize noise,” and why a single fixed regularization strength isn’t guaranteed to lower validation error the way tuning learning_rate did in the grid search.

Evaluate The article’s recursive multi-step code updates lag_1, lag_2, lag_3 with predictions but explicitly skips updating the rolling-window statistics and lag_4-lag_7 (“Rolling stats would need more history, so we skip updating them”). Critique this simplification: what happens to forecast quality at longer horizons when most of the model’s features are still real historical data rather than the model’s own compounding predictions — and why might that make the recursive approach look artificially better at horizon 6-7 than a true multi-step recursive forecast would be?

Create Design a walk-forward validation plan (following the article’s Section on time-aware splitting) for a new forecasting task: predicting weekly warehouse inventory needs, where you retrain monthly using all data available so far. Specify the fold structure (how many folds, what each trains on and tests on) and explain what business risk this validation approach specifically protects against compared to a single 80/20 time-aware split.

---

Related articles

• Feature Engineering for Time Series: Lags, Rolling Windows, and Seasonality) — the previous article: Nora built lag, rolling-window, and seasonal features for the bakery’s daily sales data, setting up the feature matrix that gradient boosting consumes here.
• Prophet vs. Statistical Models: A Practical Forecasting Showdown) — the next article: when ARIMA and boosting both break on a structural break (a new store opening), Nora reaches for Prophet.

References & Further reading

• Ke, G., Meng, Q., Finley, T., Wang, T., Chen, W., Ma, W., Ye, Q., & Liu, T.-Y. (2017). “LightGBM: A Highly Efficient Gradient Boosting Decision Tree.” Advances in Neural Information Processing Systems (NeurIPS). papers.nips.cc/book/lightgbm — the foundational LightGBM paper introducing leaf-wise (rather than level-wise) tree growth and the histogram-based gradient boosting algorithm.
• Kaggle “Store Sales — Time Series Forecasting” competition: kaggle.com/competitions/store-sales-time-series-forecasting — a real-world retail forecasting benchmark where gradient boosting with engineered lag/rolling features consistently tops the leaderboard.
• LightGBM official documentation: lightgbm.readthedocs.io — the library’s API reference and parameter tuning guide, covering num_leaves, lambda_l1/l2, min_data_in_leaf, and early stopping callbacks used throughout this article.

---