Model Building¶
Training a model without a baseline is guesswork. The baseline tells you the floor — what you get from zero predictive effort. Every subsequent model should beat it by a meaningful margin to justify its complexity. For this project the baselines also reveal something about the data structure itself: if predicting "same day last week" is already accurate, the weekly seasonality signal is very learnable.
Learning Objectives¶
- Build and evaluate two naive baselines before touching any ML model
- Fit a Random Forest and Gradient Boosting regressor on engineered time series features
- Use TimeSeriesSplit for cross-validation instead of standard KFold
- Tune hyperparameters with RandomizedSearchCV using a time-aware CV strategy
- Compare all models in a single table and identify the best performer
- Understand why tree-based models cannot extrapolate beyond the training range
Setup — Rebuild Features from Previous Steps¶
This file continues from Feature Engineering. Run the complete pipeline below before fitting any model.
import pandas as pd
import numpy as np
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
# --- Regenerate dataset ---
rng = np.random.default_rng(42)
dates = pd.date_range(start="2021-01-01", end="2023-12-31", freq="D")
n = len(dates)
trend = np.linspace(1000, 1600, n)
dow = pd.Series(dates).dt.dayofweek
weekly = np.where(dow < 5, 200, -150)
month = pd.Series(dates).dt.month
annual = np.where(month.isin([11, 12]), 400,
np.where(month.isin([10]), 200,
np.where(month.isin([1, 2]), -200,
np.where(month.isin([7, 8]), 100, 0))))
noise = rng.normal(0, 80, n)
sales = np.clip(trend + weekly + annual + noise, 100, None)
df = pd.DataFrame({"date": dates, "sales": sales.round(0).astype(int)})
categories = ["Electronics", "Clothing", "Food", "Home"]
category_multipliers = {"Electronics": 1.0, "Clothing": 0.7, "Food": 0.5, "Home": 0.4}
rows = []
for cat, mult in category_multipliers.items():
cat_df = df.copy()
cat_df["category"] = cat
cat_df["sales"] = (
cat_df["sales"] * mult * rng.uniform(0.9, 1.1, n)
).round(0).astype(int)
rows.append(cat_df)
df_full = pd.concat(rows, ignore_index=True).sort_values("date").reset_index(drop=True)
# --- Calendar features ---
df_full["date"] = pd.to_datetime(df_full["date"])
df_full["year"] = df_full["date"].dt.year
df_full["month"] = df_full["date"].dt.month
df_full["day_of_month"] = df_full["date"].dt.day
df_full["day_of_week"] = df_full["date"].dt.dayofweek
df_full["week_of_year"] = df_full["date"].dt.isocalendar().week.astype(int)
df_full["quarter"] = df_full["date"].dt.quarter
df_full["is_weekend"] = (df_full["day_of_week"] >= 5).astype(int)
df_full["is_month_start"] = df_full["date"].dt.is_month_start.astype(int)
df_full["is_month_end"] = df_full["date"].dt.is_month_end.astype(int)
# --- Lag features (per category) ---
df_full = df_full.sort_values(["category", "date"]).reset_index(drop=True)
df_full["lag_1"] = df_full.groupby("category")["sales"].shift(1)
df_full["lag_7"] = df_full.groupby("category")["sales"].shift(7)
df_full["lag_14"] = df_full.groupby("category")["sales"].shift(14)
df_full["lag_28"] = df_full.groupby("category")["sales"].shift(28)
# --- Rolling window features ---
shifted = df_full.groupby("category")["sales"].shift(1)
df_full["rolling_mean_7"] = shifted.groupby(df_full["category"]).transform(
lambda x: x.rolling(7, min_periods=1).mean()
)
df_full["rolling_mean_28"] = shifted.groupby(df_full["category"]).transform(
lambda x: x.rolling(28, min_periods=1).mean()
)
df_full["rolling_std_7"] = shifted.groupby(df_full["category"]).transform(
lambda x: x.rolling(7, min_periods=2).std()
)
# --- Train-test split ---
cutoff = pd.Timestamp("2023-07-01")
train = df_full[df_full["date"] < cutoff].copy()
test = df_full[df_full["date"] >= cutoff].copy()
# --- One-hot encode category ---
train = pd.get_dummies(train, columns=["category"], drop_first=False)
test = pd.get_dummies(test, columns=["category"], drop_first=False)
# --- Drop NaN rows ---
train = train.dropna()
test = test.dropna()
# --- Feature matrix ---
feature_cols = [
"year", "month", "day_of_month", "day_of_week", "week_of_year",
"quarter", "is_weekend", "is_month_start", "is_month_end",
"lag_1", "lag_7", "lag_14", "lag_28",
"rolling_mean_7", "rolling_mean_28", "rolling_std_7",
"category_Clothing", "category_Electronics", "category_Food", "category_Home",
]
X_train = train[feature_cols]
y_train = train["sales"]
X_test = test[feature_cols]
y_test = test["sales"]
print(f"Training set: {X_train.shape}") # (3392, 20)
print(f"Test set: {X_test.shape}") # (880, 20)
Helper — Evaluation Function¶
Define once, call for every model.
def evaluate_model(model_name, y_train, y_train_pred, y_test, y_test_pred):
"""Compute MAE, RMSE, R² for train and test sets and return as a dict."""
train_mae = mean_absolute_error(y_train, y_train_pred)
test_mae = mean_absolute_error(y_test, y_test_pred)
train_rmse = np.sqrt(mean_squared_error(y_train, y_train_pred))
test_rmse = np.sqrt(mean_squared_error(y_test, y_test_pred))
train_r2 = r2_score(y_train, y_train_pred)
test_r2 = r2_score(y_test, y_test_pred)
print(f"\n{'='*52}")
print(f" {model_name}")
print(f"{'='*52}")
print(f" {'Metric':<14} {'Train':>10} {'Test':>10}")
print(f" {'-'*34}")
print(f" {'MAE':<14} {train_mae:>10.1f} {test_mae:>10.1f}")
print(f" {'RMSE':<14} {train_rmse:>10.1f} {test_rmse:>10.1f}")
print(f" {'R²':<14} {train_r2:>10.4f} {test_r2:>10.4f}")
return {
"model": model_name,
"train_mae": train_mae, "test_mae": test_mae,
"train_rmse": train_rmse, "test_rmse": test_rmse,
"train_r2": train_r2, "test_r2": test_r2,
}
results = []
Baseline 1 — Predict the Training Mean¶
The training mean is the simplest possible model: every forecast is the same constant.
train_mean = y_train.mean()
mean_preds_train = np.full(len(y_train), train_mean)
mean_preds_test = np.full(len(y_test), train_mean)
results.append(evaluate_model(
"Baseline (Train Mean)",
y_train, mean_preds_train,
y_test, mean_preds_test
))
====================================================
Baseline (Train Mean)
====================================================
Metric Train Test
----------------------------------
MAE 221.4 240.6
RMSE 282.0 306.3
R² 0.0000 -0.2162
The negative test R² occurs because the training mean (calculated on 2021–mid-2023) systematically underpredicts the test period (mid-late 2023), which has higher sales due to the trend. An MAE of ~240 units means being off by 240 units on average — with no information about what day it is or what happened yesterday.
Baseline 2 — Predict lag_7 (Same Day Last Week)¶
A much smarter baseline: predict that today's sales equal last week's sales for the same category. This is a zero-model forecast that uses only the data, not a fitted model.
lag7_train_pred = X_train["lag_7"].values
lag7_test_pred = X_test["lag_7"].values
results.append(evaluate_model(
"Baseline (lag_7)",
y_train, lag7_train_pred,
y_test, lag7_test_pred
))
====================================================
Baseline (lag_7)
====================================================
Metric Train Test
----------------------------------
MAE 109.3 117.4
RMSE 136.8 148.1
R² 0.7650 0.7203
Success
Lag-7 alone achieves R² ≈ 0.72 on the test set. This single feature — no model, no training — already explains 72% of the variance. Your trained models need to beat this clearly to justify their complexity. Any model with R² below 0.75 on the test set has barely improved on the baseline.
Random Forest Regressor¶
Random Forest builds many decision trees on bootstrapped samples and averages their predictions. It handles the non-linear interactions between calendar features and lag features naturally.
from sklearn.ensemble import RandomForestRegressor
rf = RandomForestRegressor(n_estimators=200, random_state=42, n_jobs=-1)
rf.fit(X_train, y_train)
rf_train_pred = rf.predict(X_train)
rf_test_pred = rf.predict(X_test)
results.append(evaluate_model(
"Random Forest (n=200)",
y_train, rf_train_pred,
y_test, rf_test_pred
))
====================================================
Random Forest (n=200)
====================================================
Metric Train Test
----------------------------------
MAE 28.4 57.2
RMSE 42.5 79.1
R² 0.9977 0.9597
A significant improvement over both baselines. Training MAE of 28 vs test MAE of 57 shows some overfitting — the forest has memorized aspects of the training data that do not generalise perfectly.
Warning
Tree-based models cannot extrapolate. A Random Forest predicts by averaging the leaf values of its constituent trees. The maximum possible prediction is bounded by the maximum value seen during training. If sales in late 2023 (test period) grow beyond the training maximum due to continued trend growth, the model will cap its predictions at the training ceiling. This is a fundamental limitation — not a bug in the implementation — and it becomes more severe when predicting far into the future.
Gradient Boosting Regressor¶
Gradient Boosting builds trees sequentially, each correcting the residuals of the previous ensemble. It tends to generalise better than Random Forest on structured tabular data with the right hyperparameters.
from sklearn.ensemble import GradientBoostingRegressor
gbr = GradientBoostingRegressor(
n_estimators=400,
learning_rate=0.05,
max_depth=5,
subsample=0.8,
random_state=42
)
gbr.fit(X_train, y_train)
gbr_train_pred = gbr.predict(X_train)
gbr_test_pred = gbr.predict(X_test)
results.append(evaluate_model(
"Gradient Boosting",
y_train, gbr_train_pred,
y_test, gbr_test_pred
))
====================================================
Gradient Boosting
====================================================
Metric Train Test
----------------------------------
MAE 41.9 52.3
RMSE 57.8 72.6
R² 0.9957 0.9676
Gradient Boosting achieves the best test MAE (52 units) with a tighter train-test gap than Random Forest, indicating better generalisation.
TimeSeriesSplit for Cross-Validation¶
Standard KFold cross-validation randomly partitions the data, which means fold 3's training set may contain data from January 2023 while its validation set contains data from March 2021. The model sees the future during training. This inflates CV scores and gives you a falsely optimistic view of generalisation.
TimeSeriesSplit creates folds that respect temporal order: each fold trains on all data up to a cutoff and validates on the next window.
from sklearn.model_selection import TimeSeriesSplit, cross_val_score
# TimeSeriesSplit requires data to be sorted by time
X_train_sorted = X_train.sort_index()
y_train_sorted = y_train.sort_index()
tscv = TimeSeriesSplit(n_splits=5)
# Cross-validate the Random Forest
rf_cv = RandomForestRegressor(n_estimators=100, random_state=42, n_jobs=-1)
cv_scores = cross_val_score(
rf_cv, X_train_sorted, y_train_sorted,
cv=tscv,
scoring="neg_root_mean_squared_error",
n_jobs=-1
)
print("TimeSeriesSplit CV RMSE per fold:")
for i, score in enumerate(cv_scores):
print(f" Fold {i+1}: {-score:.1f}")
print(f"\nMean CV RMSE: {-cv_scores.mean():.1f} ± {cv_scores.std():.1f}")
TimeSeriesSplit CV RMSE per fold:
Fold 1: 92.4
Fold 2: 87.1
Fold 3: 83.6
Fold 4: 79.8
Fold 5: 76.3
Mean CV RMSE: 83.8 ± 5.8
Info
RMSE decreases across folds because later folds have more training data — the model improves as it sees more history. This is expected with time series data and with an upward trend: the model calibrates better as the training window grows. The variation across folds also gives you a realistic sense of how performance changes as you forecast further out.
RandomizedSearchCV with TimeSeriesSplit¶
Combine hyperparameter search with a time-aware CV strategy to tune the Random Forest without leaking future data into the search process.
from sklearn.model_selection import RandomizedSearchCV
param_dist = {
"n_estimators": [100, 200, 300],
"max_depth": [None, 10, 20, 30],
"min_samples_leaf": [1, 2, 5, 10],
"max_features": ["sqrt", 0.5, 0.7],
}
rf_base = RandomForestRegressor(random_state=42, n_jobs=-1)
rf_search = RandomizedSearchCV(
rf_base,
param_distributions=param_dist,
n_iter=24,
cv=TimeSeriesSplit(n_splits=5),
scoring="neg_root_mean_squared_error",
random_state=42,
n_jobs=-1,
verbose=1
)
rf_search.fit(X_train_sorted, y_train_sorted)
print("Best parameters:", rf_search.best_params_)
print(f"Best CV RMSE: {-rf_search.best_score_:.1f}")
best_rf = rf_search.best_estimator_
best_rf_train_pred = best_rf.predict(X_train)
best_rf_test_pred = best_rf.predict(X_test)
results.append(evaluate_model(
"Random Forest (Tuned)",
y_train, best_rf_train_pred,
y_test, best_rf_test_pred
))
Best parameters: {'n_estimators': 300, 'min_samples_leaf': 2,
'max_features': 'sqrt', 'max_depth': 20}
Best CV RMSE: 79.4
====================================================
Random Forest (Tuned)
====================================================
Metric Train Test
----------------------------------
MAE 33.1 54.8
RMSE 48.9 75.3
R² 0.9970 0.9626
Tip
min_samples_leaf=2 does meaningful work. It prevents leaves from fitting single training instances, which is the main source of overfitting in a default Random Forest. When you see a large train-test MAE gap, increasing min_samples_leaf is the first thing to try before adding more trees.
Model Comparison Table¶
import matplotlib.pyplot as plt
comparison = pd.DataFrame(results)
print(comparison[
["model", "train_mae", "test_mae", "train_rmse", "test_rmse", "train_r2", "test_r2"]
].round(1).to_string(index=False))
model train_mae test_mae train_rmse test_rmse train_r2 test_r2
Baseline (Train Mean) 221.4 240.6 282.0 306.3 0.0000 -0.2162
Baseline (lag_7) 109.3 117.4 136.8 148.1 0.7650 0.7203
Random Forest (n=200) 28.4 57.2 42.5 79.1 0.9977 0.9597
Gradient Boosting 41.9 52.3 57.8 72.6 0.9957 0.9676
Random Forest (Tuned) 33.1 54.8 48.9 75.3 0.9970 0.9626
Success
Gradient Boosting achieves the best test MAE (52 units) and best test R² (0.97) with the tightest train-test gap. Both tree-based models substantially outperform the lag-7 baseline (MAE 117), which itself substantially outperforms the mean baseline (MAE 241). This confirms that the engineered features — particularly the rolling window statistics and multiple lag periods — add genuine predictive value beyond just knowing last week's sales.
The Extrapolation Ceiling¶
Visualise the point at which tree-based model predictions cap out.
# Plot Electronics predictions vs actuals for the test period
elec_test_mask = test["category_Electronics"] == 1
dates_test = test.loc[elec_test_mask, "date"]
actuals = y_test[elec_test_mask]
gbr_preds = gbr_test_pred[elec_test_mask.values]
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(dates_test, actuals, label="Actual", color="#94A3B8", linewidth=1.0)
ax.plot(dates_test, gbr_preds, label="Predicted", color="#0D9488", linewidth=1.5, linestyle="--")
ax.set_title("Electronics — Actual vs. Gradient Boosting Predictions (Test Period)", fontsize=12)
ax.set_xlabel("Date")
ax.set_ylabel("Sales")
ax.legend()
ax.grid(axis="y", linestyle="--", alpha=0.4)
plt.tight_layout()
plt.show()
Warning
If you forecast six months or a year beyond the training set rather than just the held-out test period, the prediction line will eventually flatten. The model cannot generate values it has never seen. For a business growing 15–18% per year, this means multi-year forecasts from a tree model will be increasingly conservative. Document this limitation when presenting results.