ECON 8310: Business Forecasting — University of Nebraska at Omaha
What we are doing
At the end of Lab 4 you built a bagged ensemble by accident. You drew twenty bootstrap samples, fitted a tree to each, averaged the predictions, and the RMSE fell from about 324 to 266 — with nothing tuned and no new information added. What you did not do is the one thing that turns bagging into a Random Forest.
Today we do that properly, then find out what else the real algorithm does, then replace averaging with something better.
Same series throughout: daily units for CA_1 FOODS. You have now fitted a GAM to it (Homework 2) and a single tree (Lab 4), so every number today lands on a scale you already know.
Suggested time: ~40 minutes · You will need:pandas, numpy, matplotlib, scikit-learn, xgboost, and data/processed/m5_daily.csv.
Setup
import pandas as pdimport numpy as npimport matplotlib.pyplot as pltfrom sklearn.ensemble import RandomForestRegressorfrom sklearn.inspection import permutation_importanceimport xgboost as xgbimport warningswarnings.filterwarnings('ignore')from pathlib import Pathdef data_path(name):for base in (Path.cwd(), *Path.cwd().parents): p = base /"data"/"processed"/ nameif p.exists():return praiseFileNotFoundError(f"Could not find data/processed/{name}. ""Run `python scripts/prep_m5.py` from the repository root first." )RANDOM_STATE =42daily = pd.read_csv(data_path("m5_daily.csv"), parse_dates=["date"])daily["dow"] = daily["date"].dt.dayofweekdaily["doy"] = daily["date"].dt.dayofyearFEATURES = ["units_lag1", "units_lag7", "units_lag14", "units_lag28", "units_lag365","snap", "is_event", "dow", "doy", "avg_price"]df = daily.dropna(subset=FEATURES).copy()N_TEST =365train, test = df.iloc[:-N_TEST], df.iloc[-N_TEST:]X_train, y_train = train[FEATURES].values, train["units"].valuesX_test, y_test = test[FEATURES].values, test["units"].valuesrmse =lambda a, b: float(np.sqrt(np.mean((np.asarray(a) - np.asarray(b)) **2)))scoreboard = {"Benchmark (same day last year)": rmse(y_test, test.units_lag365.values),"Single tree (Lab 4)": 324,"20 trees averaged by hand (Lab 4)": 266,}for k, v in scoreboard.items():print(f"{k:<36}{v:>6,.0f}")
Benchmark (same day last year) 625
Single tree (Lab 4) 324
20 trees averaged by hand (Lab 4) 266
Better than your hand-built version, and much better than one tree. Two differences account for it: more trees (300 against 20), and deeper ones (Lab 4’s trees stopped at depth 5). Note what is not different yet: scikit-learn’s regression forest defaults to max_features=1.0, every feature at every split, so strictly this is still bagging. The second source of randomness, the one that makes it a random forest, comes in Step 3.
Step 2 — How many trees is enough?
The usual worry about a big model is overfitting. Adding trees to a forest is the rare knob where that worry does not apply — but it does stop paying.
counts = [5, 20, 50, 100, 300, 600]curve = []for n in counts: r = RandomForestRegressor(n_estimators=n, random_state=RANDOM_STATE, n_jobs=-1).fit(X_train, y_train) curve.append(rmse(y_test, r.predict(X_test)))for n, v inzip(counts, curve):print(f" n_estimators={n:<4} test RMSE {v:,.0f}")fig, ax = plt.subplots(figsize=(6, 3.2))ax.plot(counts, curve, marker="o")ax.set_xscale("log"); ax.set_xlabel("n_estimators (log scale)"); ax.set_ylabel("test RMSE")plt.tight_layout()
n_estimators=5 test RMSE 288
n_estimators=20 test RMSE 257
n_estimators=50 test RMSE 251
n_estimators=100 test RMSE 249
n_estimators=300 test RMSE 248
n_estimators=600 test RMSE 248
Say what the curve does after about fifty trees, and what that means for how you would spend a compute budget. More trees never hurt accuracy — but every one costs time to fit and to predict.
Step 3 — The “random” in Random Forest
Bagging alone gives every tree a different sample of rows. A Random Forest also gives each split a different random subset of columns — max_features. Lecture 5 argued this decorrelates the trees, so their errors cancel more completely when averaged.
Test it.
for mf in [None, 0.6, "sqrt", 0.2]: r = RandomForestRegressor(n_estimators=300, max_features=mf, random_state=RANDOM_STATE, n_jobs=-1).fit(X_train, y_train) label ="all features"if mf isNoneelsestr(mf)print(f" max_features={label:<13} test RMSE {rmse(y_test, r.predict(X_test)):,.0f}")
max_features=all features test RMSE 248
max_features=0.6 test RMSE 247
max_features=sqrt test RMSE 254
max_features=0.2 test RMSE 257
The textbook rule does not win here, and you should know why
max_features="sqrt" is the standard recommendation, and on this dataset it is worse than using every feature. That is not a bug — it is the rule meeting its own boundary condition.
You have 10 features, and one of them (dow) carries most of the signal. sqrt(10) is about 3, so at any given split most trees are not even allowed to look at day-of-week. Decorrelation is only a bargain when you have enough redundant, correlated features that losing a few at each split costs you little.
Do not expect Homework 3 to restore the textbook answer: its forest has 27 features, but units_lag1 dominates there the way dow dominates here, and sqrt loses there too. The lesson is the condition, not the default.
Step 4 — Validation for free
Each bootstrap sample leaves out about a third of the rows. Those out-of-bag rows never trained that tree, so they can score it — a validation estimate with no holdout set.
Compare the two numbers and account for the gap. The OOB estimate is computed on training-period rows scattered through 2012–2015; the test RMSE is a solid block of the final year. Which one would you quote to someone asking how the model will do next month, and why?
Permutation importance is in units of RMSE: shuffle that column on the test set, and this is how much worse the forecast gets. That is directly interpretable in a way MDI is not.
Find the disagreements. Three are worth naming before you read on:
one lag that MDI ranks aboveunits_lag1, while permutation ranks it far below
snap, which MDI places near the bottom and permutation places near the top
avg_price against is_event: MDI ranks the price column above the event flag, while permutation scores price at exactly zero and the event flag at 22 units of RMSE
Your turn: plot the two rankings side by side as horizontal bar charts so the disagreements are visible at a glance.
# ---- your turn: uncomment and fill in the blank ----# fig, axes = plt.subplots(1, 2, figsize=(11, 3.5))# mdi.sort_values().plot.barh(ax=axes[0], title="MDI (impurity)")# ___.sort_values().plot.barh(ax=axes[1], title="Permutation (ΔRMSE)") # <- which Series?# plt.tight_layout()
Step 6 — Stop averaging. Start correcting.
A forest builds many trees in parallel and averages them. Boosting builds them in sequence, each one fitted to what the previous ones got wrong.
The first four rows hold the learning rate fixed and add rounds — the same experiment you ran on the forest in Step 2. The last two change the rate instead.
for lr, n in [(0.1, 100), (0.1, 300), (0.1, 1000), (0.1, 3000), (0.3, 100), (0.05, 1000)]: m = xgb.XGBRegressor(n_estimators=n, learning_rate=lr, max_depth=4, subsample=0.8, colsample_bytree=0.8, random_state=RANDOM_STATE).fit(X_train, y_train)print(f" lr={lr:<5} n={n:<5} train {rmse(y_train, m.predict(X_train)):5,.0f}"f" test {rmse(y_test, m.predict(X_test)):5,.0f}")
lr=0.1 n=100 train 171 test 244
lr=0.1 n=300 train 95 test 243
lr=0.1 n=1000 train 19 test 250
lr=0.1 n=3000 train 0 test 253
lr=0.3 n=100 train 95 test 238
lr=0.05 n=1000 train 54 test 239
Benchmark (same day last year) 625
Single tree (Lab 4) 324
20 trees averaged by hand (Lab 4) 266
Random Forest (300 trees) 248
XGBoost (lr=0.1, 100 trees) 244
Read the first four rows as a curve, then compare it with Step 2. Training error falls to zero. Test error bottoms out around 300 rounds and then climbs back up — 3,000 rounds is worse than 100. In Step 2 you added trees to a forest all the way to 600 and it never once hurt.
That is the practical difference between the two ensembles. A forest’s trees are fitted independently and averaged, so extra ones can only sharpen the average. Boosting’s trees are fitted to the previous ensemble’s residuals, so extra ones keep going — eventually correcting noise. n_estimators is a free choice in a forest and a tuning parameter in boosting.
The last two rows say the obvious complement: a big rate with few rounds and a small rate with many land in much the same place, because roughly it is the product that sets the budget. What you cannot do is set a rate and then let the rounds run.
Note which configuration we put on the scoreboard below — lr=0.1, n=100, the conventional default, not the row with the best test RMSE. Picking the winner out of a column computed on the test set is the mistake Lecture 1 warned about, and it does not stop being a mistake because the grid is small.
Before you leave
Your scoreboard should run from a benchmark near 625 down to roughly 245, across five models fitted to one series.
Discuss (the ranking): Boosting beat the forest here. Name one property of this dataset that favors boosting, and one situation where you would still reach for the forest first — thinking about tuning effort and about what happens when nobody retunes the model for a year.
Discuss (the importances): You have two importance tables that disagree about snap. A category manager asks whether SNAP timing matters for stocking. Which table do you show them, what do you say, and what would you want to check before committing to an answer?
Discuss (the whole arc): Across Labs 4 and 5 the RMSE fell from 324 to about 245 on the same series with the same features. What was actually purchased with that improvement, and what was given up along the way? Be specific about the thing you can no longer do that you could do in Step 1 of Lab 4.
