import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression, Lasso, Ridge
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import TimeSeriesSplit, cross_val_score
import warnings
warnings.filterwarnings('ignore')
from pathlib import Path
def data_path(name):
for base in (Path.cwd(), *Path.cwd().parents):
p = base / "data" / "processed" / name
if p.exists():
return p
raise FileNotFoundError(
f"Could not find data/processed/{name}. "
"Run `python scripts/prep_m5.py` from the repository root first."
)
RANDOM_STATE = 42
daily = pd.read_csv(data_path("m5_daily.csv"), parse_dates=["date"]).sort_values("date")
u = daily.set_index("date")["units"]
for L in list(range(1, 29)) + [35, 42, 49, 56, 364, 365, 366]:
daily[f"lag{L}"] = u.shift(L).values
for W in (7, 14, 28, 91, 364):
daily[f"roll{W}_mean"] = u.shift(1).rolling(W).mean().values
daily[f"roll{W}_std"] = u.shift(1).rolling(W).std().values
daily["dow"] = daily["date"].dt.dayofweek
daily["doy"] = daily["date"].dt.dayofyear
for k in (1, 2, 3):
daily[f"sin{k}"] = np.sin(2*np.pi*k*daily.doy/365)
daily[f"cos{k}"] = np.cos(2*np.pi*k*daily.doy/365)
daily = pd.concat([daily, pd.get_dummies(daily.dow, prefix="dow")], axis=1)
FEATURES = ([c for c in daily.columns if c.startswith(("lag", "roll", "sin", "cos", "dow_"))]
+ ["snap", "is_event", "doy", "avg_price"])
df = daily.dropna(subset=[c for c in FEATURES if c.startswith(("lag", "roll"))]).copy()
N_TEST = 365
train, test = df.iloc[:-N_TEST], df.iloc[-N_TEST:]
X_train, y_train = train[FEATURES].values.astype(float), train["units"].values
X_test, y_test = test[FEATURES].values.astype(float), test["units"].values
rmse = lambda a, b: float(np.sqrt(np.mean((np.asarray(a) - np.asarray(b)) ** 2)))
print(f"k = {len(FEATURES)} features, n_train = {len(X_train)}, n/k = {len(X_train)/len(FEATURES):.1f}")6 Lab 6: Regularization — Breaking OLS, Then Rescuing It
Run this one live: the executable lab.
6.1 What we are doing
Lecture 6 claimed that ordinary least squares breaks when you have many predictors relative to observations — coefficient variance grows with \(k\), correlated predictors produce wild offsetting estimates, and past \(k > n\) there is no unique solution at all.
Today you will break it on purpose, watch exactly how it fails, and then fix it with one line.
Same series as Labs 4 and 5: daily units for CA_1 FOODS. But we now engineer the feature set an analyst would actually build — 28 consecutive lags, weekly and annual lags, rolling statistics, Fourier terms, day-of-week dummies — which is how you end up with more predictors than you meant to have.
Suggested time: ~40 minutes · You will need: pandas, numpy, matplotlib, scikit-learn, and data/processed/m5_daily.csv.
6.2 Setup
Note how easily 62 predictors appeared. Nobody set out to build a 62-parameter model; it is what routine feature engineering produces.
Walkthrough
Look at what the setup does, because nothing in it is unreasonable.
Twenty-eight consecutive lags, because you don’t know which one matters. A few weekly lags — 35, 42, 49, 56. Last year’s lags, 364, 365, 366. Rolling means and standard deviations at five windows. Three Fourier pairs for annual seasonality. Day-of-week dummies. Then SNAP, events, day-of-year and price.
k = 62 features, n_train = 1210, n/k = 19.5Sixty-two predictors. Thirty-five lags, ten rolling statistics, six Fourier terms, seven day-of-week dummies, four others.
Nobody sat down and decided to fit a sixty-two-parameter model. Every one of those columns was a defensible small decision, and the total is what routine feature engineering produces. That is the honest version of how this problem arrives: not from ambition, from diligence.
Right now the ratio is comfortable — about twenty observations per parameter. Hold onto that number, because Step 1 is going to take it apart.
6.3 Step 1 — Break it
With 1,210 training rows and 62 features, OLS is fine. Now shrink the training window and watch.
rows = []
for n in (1210, 400, 200, 120, 80):
sub = train.iloc[-n:]
sc = StandardScaler().fit(sub[FEATURES].values.astype(float))
Xs = sc.transform(sub[FEATURES].values.astype(float))
Xt = sc.transform(X_test)
y = sub["units"].values
ols = LinearRegression().fit(Xs, y)
rid = Ridge(alpha=100).fit(Xs, y)
rows.append({"n": n, "n/k": round(n/len(FEATURES), 1),
"OLS_max_coef": f"{np.abs(ols.coef_).max():,.0f}",
"OLS_test_RMSE": f"{rmse(y_test, ols.predict(Xt)):,.0f}",
"Ridge_test_RMSE": f"{rmse(y_test, rid.predict(Xt)):,.0f}"})
print(pd.DataFrame(rows).to_string(index=False))Read that table across, then down. As the training window shrinks, the largest OLS coefficient grows from about 150 to the hundreds of thousands, and the test RMSE from 280 into the millions. Nothing about the problem changed — same features, same target, same test set. Only the ratio of observations to parameters.
Ridge, with a one-line change, degrades gracefully across the entire range.
That is Lecture 6’s opening argument, and you just measured it.
Walkthrough
We keep the features and the test set fixed and shrink only the training window. Four years, then about one, then two hundred days, then a hundred and twenty, then eighty. Same problem throughout. Only the ratio changes.
n n/k OLS_max_coef OLS_test_RMSE Ridge_test_RMSE
1210 19.5 153 280 278
400 6.5 636 2,071 323
200 3.2 1,364 5,874 363
120 1.9 1,049,739 3,578,141 621
80 1.3 355,629 2,072,167 770Read across the OLS columns. At twenty observations per parameter, the largest coefficient is 153 and the test RMSE is 280 — a perfectly respectable model. At six observations per parameter, the coefficient is 636 and the error has multiplied by seven. At three, it is 1,364 and 5,874.
And then the bottom two rows stop being a degradation and become an explosion. The largest coefficient is over a million. The test RMSE is three and a half million — on a series where a busy day is four thousand units.
Now here is the part I want you to take away, and it is not in the printed table.
At n = 120, that model’s training RMSE is about 145. It looks fine. If you fitted it and looked only at how it fits the data you gave it, you would ship it. It is the best-fitting model on this slide and it is wrong by a factor of a thousand on data it has not seen.
One technical point, because a careful student will ask. We never crossed the line the lecture warned about. Sixty-two features and eighty rows — there are still more observations than parameters, so OLS has a unique solution the whole way down this table. You do not need k greater than n to get this. Twenty-eight consecutive lags of a smooth series are nearly collinear with one another, and near-collinearity is enough. The matrix does not have to be singular, only close to it.
Now the right-hand column. Same data, same features, same rows — Ridge(alpha=100) instead of LinearRegression(). One word changed. 278, 323, 363, 621, 770.
It gets worse as data gets scarcer, which it should. But it degrades gracefully, and at the row where OLS is out by three and a half million, Ridge is at 621 — roughly the benchmark. Not good, but the model is still a model. That is Lecture 6’s opening argument, and you just measured it.
6.4 Step 2 — Two penalties, two behaviors
Back to the full training set. Fit Ridge and LASSO at a comparable penalty and count what survives.
sc = StandardScaler().fit(X_train)
Xs, Xts = sc.transform(X_train), sc.transform(X_test)
for name, model in [("OLS", LinearRegression()),
("Ridge", Ridge(alpha=100)),
("LASSO", Lasso(alpha=20, max_iter=50000))]:
m = model.fit(Xs, y_train)
nz = int(np.sum(m.coef_ != 0))
print(f" {name:<6} test RMSE {rmse(y_test, m.predict(Xts)):5,.0f} non-zero coefficients {nz:3d}/{len(FEATURES)}")Ridge keeps every coefficient and shrinks them. LASSO sets most of them to exactly zero. Neither is “better” — they encode different beliefs about how many predictors genuinely matter.
Walkthrough
Back to the full training set, where OLS was fine, and now compare the two penalties.
OLS test RMSE 280 non-zero coefficients 62/62
Ridge test RMSE 278 non-zero coefficients 62/62
LASSO test RMSE 286 non-zero coefficients 21/62I want you to notice the disappointing thing first, because if I skip it you will notice it anyway and trust me less.
On the full training set, regularization buys essentially nothing. Ridge is two units better than OLS. LASSO is six units worse. With nineteen observations per parameter there was no problem to solve, and the penalty solved it.
That is the correct lesson, not a failed demo. Regularization is insurance, and Step 1 is the fire. You do not judge insurance by what it pays out in a year when the house doesn’t burn down.
But look at the last column, because that is where the two penalties differ in kind and not just in degree.
Ridge keeps all sixty-two coefficients and shrinks them toward zero. LASSO sets forty-one of them to exactly zero and keeps twenty-one.
So for six units of RMSE — about two percent — LASSO handed you a model with a third as many moving parts. That is not an accuracy trade, it is a communicability trade. Twenty-one coefficients is a model you can print and discuss. Sixty-two is a model you can only run.
Neither penalty is better. They encode different beliefs: Ridge believes many predictors each matter a little, LASSO believes a few matter and the rest are noise. Which belief is right is a question about your data, and in Step 3 you get to see which one this series votes for.
6.5 Step 3 — The path, and the order features die in
Lecture 6 argued that the order predictors leave is a more stable finding than any single coefficient table. Produce that ordering.
alphas = np.logspace(-1, 3, 60)
paths = []
for a in alphas:
paths.append(Lasso(alpha=a, max_iter=50000).fit(Xs, y_train).coef_)
paths = np.array(paths)
fig, axes = plt.subplots(1, 2, figsize=(11, 3.6))
axes[0].plot(np.log10(alphas), paths, lw=0.8)
axes[0].set_xlabel("log10(alpha)"); axes[0].set_ylabel("coefficient"); axes[0].set_title("regularization path")
axes[1].plot(np.log10(alphas), (paths != 0).sum(1), marker=".", ms=3)
axes[1].set_xlabel("log10(alpha)"); axes[1].set_ylabel("non-zero coefficients"); axes[1].set_title("model size")
plt.tight_layout()survivors, order = set(FEATURES), []
for a, coefs in zip(alphas, paths):
nz = set(np.array(FEATURES)[coefs != 0])
for gone in survivors - nz:
order.append(gone)
survivors = nz
print("the last 8 features standing, in reverse order of elimination:")
for f in reversed(order[-8:]):
print(f" {f}")
Look at what survived
Read that list of survivors and find the pattern in the lag numbers before reading on.
Almost all of them are multiples of 7. Nobody told the model this series has a weekly cycle, and there is no explicit “week” feature among the lags. Under increasing penalty the LASSO independently kept the lags that align with day-of-week and discarded the ones that do not.
That is what makes the dropout ordering worth reporting: it recovered the structure of the series, not just a ranking of columns.
Walkthrough
Sixty LASSO fits along a grid of penalties, from 0.1 up to 1,000, and we plot every coefficient’s trajectory.
Each line is one feature. At the left, weak penalty, everything is alive. Move right and lines collapse onto zero one after another. By the right-hand edge — and I checked this — every single coefficient is zero. The model at alpha = 1,000 predicts the mean and nothing else.
The right panel is the same thing counted: model size against penalty. It is a staircase down from 62 to 0. At alpha = 100, nine features are left.
Lecture 6 argued that the order features leave is a more stable finding than any one coefficient table — coefficients wobble under resampling, the elimination order mostly does not. So let’s read the order.
the last 8 features standing, in reverse order of elimination:
lag28
lag7
lag35
lag364
lag14
lag21
lag1
lag56
Try it
Pause here. Read those eight lag numbers and find the pattern before I say it.
- Every one of them is a multiple of seven except
lag1.
- Every one of them is a multiple of seven except
Now think about what was available to be chosen. The feature set had twenty-eight consecutive lags — 1, 2, 3, all the way to 28. There is no “week” feature. Nothing in the data told the model that seven days is a meaningful interval. The day-of-week dummies were there as an alternative route to the same information, and the LASSO dropped those too.
Under increasing penalty, forced to give things up, it kept exactly the lags that align with day-of-week and discarded the ones that do not. It also kept lag364 — which is fifty-two weeks, so the same-day-last-year lag that is also a multiple of seven, rather than lag365.
That is why the dropout ordering is worth reporting to someone. It did not just rank your columns. It recovered the structure of the series — a weekly cycle plus an annual echo — from a pile of undifferentiated lags. And it agrees with what the tree found on its own in Lab 4, when day-of-week took 79% of the importance. Two completely different model families, same finding.
6.6 Step 4 — Choosing alpha, and the one-standard-error rule
tscv = TimeSeriesSplit(n_splits=5)
grid = np.logspace(-1, 2.5, 30)
means, ses = [], []
for a in grid:
s = cross_val_score(make_pipeline(StandardScaler(), Lasso(alpha=a, max_iter=50000)),
X_train, y_train, cv=tscv, scoring="neg_root_mean_squared_error")
means.append(-s.mean()); ses.append(s.std()/np.sqrt(len(s)))
means, ses = np.array(means), np.array(ses)
i_min = int(means.argmin())
threshold = means[i_min] + ses[i_min]
i_1se = int(np.where(means <= threshold)[0].max())
for label, i in [("CV minimum", i_min), ("one-SE rule", i_1se)]:
m = make_pipeline(StandardScaler(), Lasso(alpha=grid[i], max_iter=50000)).fit(X_train, y_train)
print(f" {label:<12} alpha {grid[i]:7.2f} CV {means[i]:5,.0f} "
f"test {rmse(y_test, m.predict(X_test)):5,.0f} non-zero {int(np.sum(m[-1].coef_ != 0)):3d}")Your turn: plot the CV curve with its one-standard-error band, and mark both chosen alphas.
# ---- your turn: uncomment and fill in the blank ----
# fig, ax = plt.subplots(figsize=(7, 3.4))
# ax.plot(np.log10(grid), means)
# ax.fill_between(np.log10(grid), means - ses, means + ses, alpha=0.2)
# ax.axvline(np.log10(grid[i_min]), ls="--", label="CV minimum")
# ax.axvline(np.log10(___), ls=":", label="one-SE") # <- which index?
# ax.set_xlabel("log10(alpha)"); ax.set_ylabel("CV RMSE"); ax.legend()
# plt.tight_layout()The one-SE rule roughly halves the number of surviving features. Decide whether that trade is worth it here — and say what you would need to know about how the model will be used before you could answer properly.
Walkthrough
So far I have been picking penalties by hand. Now do it properly: five-fold TimeSeriesSplit, thirty candidate alphas, and — this matters — the scaler goes inside the pipeline, so it is refitted on each fold’s training rows and never sees the validation fold.
CV minimum alpha 19.64 CV 323 test 286 non-zero 21
one-SE rule alpha 59.68 CV 346 test 309 non-zero 11Two ways to read the same CV curve.
The CV minimum takes the alpha with the best cross-validated score. Standard, and it lands on 21 features.
The one-standard-error rule takes the largest alpha — the simplest model — whose CV score is still within one standard error of the best. The argument is that the CV curve is itself an estimate with noise on it, so a model that is statistically indistinguishable from the best but much simpler is the better bet.
Here that means going from 21 features to 11, and paying 23 units of test RMSE for it — 286 up to 309, about eight percent.
Try it
Two things to do. First, uncomment the plot and fill in the blank — the one-SE index — so you can see the band and both verticals. Then answer the question the lab asks: is that trade worth it here?
I’ll give you the shape of the answer rather than the answer, because the honest reply is that the lab has not told you enough to decide.
If this model is going into a dashboard that a person reads and questions, eleven coefficients against twenty-one is a real gain, and eight percent is cheap. If it is one component inside an automated replenishment system that nobody reads, you just paid eight percent for a property nobody consumes.
It is worth knowing what the eleven survivors are: lag1, lag7, lag14, lag21, lag27, lag28, lag35, lag56, lag364, the 364-day rolling mean, and snap. Nine lags, a level, and the SNAP flag. That is a model you could describe out loud in one sentence — and notice snap survived to the very end, which is the same feature permutation importance rescued in Lab 5.
6.7 Step 5 — Why your CV number looks worse than the truth
The CV RMSE above is noticeably higher than the test RMSE. That is not a mistake, and the reason is sitting in Step 1.
for i, (tr_idx, va_idx) in enumerate(tscv.split(X_train), 1):
m = make_pipeline(StandardScaler(), Lasso(alpha=20, max_iter=50000)).fit(X_train[tr_idx], y_train[tr_idx])
print(f" fold {i}: train n={len(tr_idx):4d} (n/k={len(tr_idx)/len(FEATURES):4.1f}) "
f"val n={len(va_idx):3d} RMSE {rmse(y_train[va_idx], m.predict(X_train[va_idx])):5,.0f}")Fold 1 trains on about 200 rows — the same data-starved regime that wrecked OLS in Step 1. The final model trains on all 1,210. So walk-forward CV on an expanding window is structurally pessimistic: it averages in early folds that are trained on far less data than the model you actually ship.
That is a reason to read the CV number as a conservative estimate, not a reason to abandon it — and emphatically not a reason to switch to random folds, which fix the pessimism by letting the model see the future.
Walkthrough
One loose end from Step 4 that you should not let pass. The CV score at the minimum was 323. The test RMSE of that same model was 286. Cross-validation was pessimistic by about 13%.
If those had gone the other way you would be worried. This direction has a specific cause, and it is sitting in Step 1.
fold 1: train n= 205 (n/k= 3.3) val n=201 RMSE 433
fold 2: train n= 406 (n/k= 6.5) val n=201 RMSE 294
fold 3: train n= 607 (n/k= 9.8) val n=201 RMSE 325
fold 4: train n= 808 (n/k=13.0) val n=201 RMSE 257
fold 5: train n=1009 (n/k=16.3) val n=201 RMSE 308Look at fold 1. It trains on 205 rows — three observations per parameter. That is the row from Step 1’s table where OLS was already at 5,874, and it is scored as though it were the model you are going to ship. It comes in at 433, far worse than any other fold. The model you actually ship trains on all 1,210 rows. No fold ever saw that much data.
So expanding-window cross-validation is structurally pessimistic: it averages in early folds fitted on a fraction of your data. And you can put a number on how much that costs. The five folds average to 323. Drop fold 1 alone and the remaining four average to 296 — against a test RMSE of 286. Almost the entire gap was that one starved fold. Two conclusions, and the second is the one that gets violated in practice.
First: read a walk-forward CV number as a conservative estimate. It is a lower bound on your model’s quality, not an unbiased estimate of it, and being conservative is a fine property in a forecast you are going to stake inventory on.
Second, and please do not do this: the pessimism is not a reason to switch to random k-fold splits. Random folds would fix the number, and they would fix it by training on Thursday to predict the previous Tuesday. The optimistic number you get back would be measuring your model’s ability to interpolate inside a period it has already seen, which is not the job. A pessimistic estimate of the right quantity beats an accurate estimate of the wrong one.
6.8 Before you leave
Add today’s best linear model to the scoreboard you have been building since Lab 4.
best = make_pipeline(StandardScaler(), Lasso(alpha=grid[i_min], max_iter=50000)).fit(X_train, y_train)
board = {
"Benchmark (same day last year)": rmse(y_test, test.units_lag365.values),
"Single tree (Lab 4)": 324,
"Random Forest (Lab 5)": 248,
"XGBoost (Lab 5)": 244,
"LASSO, 62 features (today)": rmse(y_test, best.predict(X_test)),
}
print(pd.Series(board).round(0).astype(int).sort_values(ascending=False).to_string())Discuss: The regularized linear model does not win here — the tree ensembles beat it on this series. Give one property of daily store sales that favors trees over an additive linear model, however many lags you give it.
Then the counterpoint, which matters for Homework 4: on the weekly, 30-series panel the ranking may not look the same. What is different about that problem that could change the answer?
And the one to argue about: you have now fitted five model families to this one series. If you had to hand exactly one of them to a store manager who will not be retrained, and who will not retune anything for a year, which do you choose — and is it the one with the lowest RMSE?
Solutions posted to Canvas after class.
Walkthrough
Benchmark (same day last year) 625
Single tree (Lab 4) 324
LASSO, 62 features (today) 286
Random Forest (Lab 5) 248
XGBoost (Lab 5) 244Today’s model, with six times as many features as the tree ensembles had, lands third.
Say that plainly rather than skating past it. Sixty-two engineered features, a properly cross-validated penalty, and a careful pipeline produced 286 — comfortably better than one tree, and comfortably worse than a Random Forest you fitted in two lines with no feature engineering at all in Lab 5.
Try it
Three questions for the board.
Why trees win here. Name one property of daily store sales that favours a tree ensemble over an additive linear model, however many lags you feed it. It is in the tree you read out loud in Lab 4 — think about what happens on a weekend when yesterday was already busy, and whether a model that adds up separate effects can express that.
The counterpoint, and this one matters for Homework 4. In Homework 4 you run this comparison on the weekly, thirty-series panel with forty-six features. The ranking may not come out the same. What is different about that problem that could change the answer? Think about how many observations per series you get at weekly frequency, and what pooling thirty series does to the noise.
The one to argue about. You have now fitted five model families to this one series — a GAM in Homework 2, a tree, a forest, boosting, and a regularized linear model. Hand exactly one of them to a store manager who will never retune it and will not be retrained. Which, and is it the one with the lowest RMSE?
Notice the two properties this lab measured that RMSE does not see. The eleven-feature LASSO is the only model of the five you could write on a whiteboard. And Step 1 showed you which models fail loudly versus quietly — OLS at n=120 had a beautiful training fit and was wrong by a factor of a thousand. In a year, if this store’s data collection changes and the model quietly gets less data than it used to, which of your five would you rather be running unattended?
Solutions for the coded parts go up on Canvas after the deadline. Next week we leave linear models behind for good and start on neural networks.
6.9 What your run should produce
Your numbers will differ where a model is stochastic; these are the values this lab produces on the pinned environment.
Deterministic; whole lab runs in about 6 seconds.
| Quantity | Value |
|---|---|
| Feature set | 62 — 35 lags, 10 rolling stats, 6 Fourier, 7 dow dummies, 4 other · n/k = 19.5 |
| OLS as n shrinks (max|coef| / test RMSE) | 1210: 153 / 280 · 400: 636 / 2,071 · 200: 1,364 / 5,874 · 120: 1,049,739 / 3,578,141 · 80: 355,629 / 2,072,167 |
| Ridge over the same range | 278 · 323 · 363 · 621 · 770 |
| OLS at n=120 — training RMSE | 145 (not printed by the lab; worth saying out loud) |
| Full data, three models | OLS 280 (62/62) · Ridge 278 (62/62) · LASSO 286 (21/62) |
| Last 8 LASSO survivors | lag28, lag7, lag35, lag364, lag14, lag21, lag1, lag56 — all multiples of 7 except lag1 |
| Coefficients alive at alpha = 1,000 | 0 (9 at alpha = 100) |
| CV minimum | alpha 19.64 · CV 323 · test 286 · 21 features |
| One-SE rule | alpha 59.68 · CV 346 · test 309 · 11 features |
| One-SE survivors | lag1, lag7, lag14, lag21, lag27, lag28, lag35, lag56, lag364, roll364_mean, snap |
| CV folds | 433 · 294 · 325 · 257 · 308 → mean 323; without fold 1, 296 (test 286) |
| Scoreboard | benchmark 625 · tree 324 · LASSO 286 · RF 248 · XGBoost 244 |