import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller, grangercausalitytests
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima.model import ARIMA
import warnings
warnings.filterwarnings('ignore')
import io
from contextlib import redirect_stdout
from pathlib import Path
def data_path(name):
"""Find the processed data file whether you render from the repo root or from
this file's own folder."""
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 the matching prep script from the repository root first: "
"`prep_m5.py` for m5_*.csv, `prep_fred.py` for fred_monthly.csv."
)
weekly = pd.read_csv(data_path("m5_weekly.csv"), parse_dates=["week_start"])
def get_series(store, cat):
return (weekly[(weekly.store_id == store) & (weekly.cat_id == cat)]
.set_index("week_start")["units"].sort_index().astype(float))
y = get_series("CA_1", "HOUSEHOLD")
N_TEST = 52
train, test = y.iloc[:-N_TEST], y.iloc[-N_TEST:]
print(f"{len(y)} weeks — train {len(train)}, test {len(test)}")2 Lab 2: Stationarity, ARIMA, and Granger Causality
Run this one live: the executable lab.
2.1 What we are doing
Lecture 2 gave you a rule — test for stationarity before you read an ACF — and a caveat: Granger causality is not structural causality. In the next 35 minutes you will do both, on data where each one actually bites.
The business question: you are forecasting weekly HOUSEHOLD demand at store CA_1, and you want to know whether watching the food aisle helps you predict it. Then, in Step 6, you run the same test on the macro question from the lecture slide — does unemployment predict retail sales? — and find that the answer depends on which years you include.
Suggested time: ~40 minutes · You will need: pandas, numpy, matplotlib, statsmodels, and two prepared files: data/processed/m5_weekly.csv (prep_m5.py) and data/processed/fred_monthly.csv (prep_fred.py).
2.2 Setup
Walkthrough
Same pattern as Lab 1 — imports, a path helper, a series accessor, a split.
277 weeks — train 225, test 52Same store dataset, different store and same category as Lab 1: CA_1 household. Fifty-two weeks held out.
The one new import worth noticing is redirect_stdout from contextlib, and io. That’s not statistics — that’s damage control, and I’ll show you why in Step 5.
2.3 Step 1 — Look, then test
fig, ax = plt.subplots(figsize=(11, 3.2))
ax.plot(y.index, y.values, lw=1.3)
ax.set_title("CA_1 HOUSEHOLD — weekly units"); plt.tight_layout()Notice: the level wanders. Your eye says non-stationary. Now make the test agree.
def report_adf(series, label):
stat, p, lags, nobs, crit, _ = adfuller(series.dropna(), autolag="AIC")
print(f"{label:<28} ADF {stat:>7.2f} p {p:>6.3f} lags {lags:>2} "
f"{'stationary' if p < 0.05 else 'UNIT ROOT'}")
report_adf(y, "levels")Now difference it and re-test.
dy = y.diff().dropna()
report_adf(dy, "first difference")You should see the statistic move decisively. Write down both numbers — you will be asked to justify a choice of \(d\) in the homework, and this is what a clean case looks like.
Walkthrough
Look first. Always look first.
The level wanders. It drifts up, it drifts back, there’s no fixed value it returns to. Your eye says non-stationary. Now make the test agree — because your eye is not evidence you can put in a report.
levels ADF -1.33 p 0.614 lags 12 UNIT ROOT
first difference ADF -5.00 p 0.000 lags 11 stationaryRead those two lines carefully, because this is what a clean case looks like and you want it in your memory for comparison later.
In levels: ADF statistic minus one point three three, p-value point six one. Nowhere near significance. We fail to reject the null of a unit root — the series is non-stationary, exactly as your eye said. First difference: minus five point zero zero, p-value effectively zero. Decisive. One difference was enough.
That tells you d equals one, and you now have an evidence trail for that choice rather than a hunch. You’ll be asked to justify a d in the homework. This is what the justification looks like.
2.4 Step 2 — Now you are allowed to read the ACF
Lecture 2’s rule: the ACF of a non-stationary series decays slowly no matter what the underlying process is. Run both and see it.
fig, axes = plt.subplots(2, 2, figsize=(11, 5))
plot_acf(y, lags=30, ax=axes[0,0]); axes[0,0].set_title("ACF — levels (uninformative)")
plot_pacf(y, lags=30, ax=axes[0,1]); axes[0,1].set_title("PACF — levels")
plot_acf(dy, lags=30, ax=axes[1,0]); axes[1,0].set_title("ACF — differenced")
plot_pacf(dy, lags=30, ax=axes[1,1]); axes[1,1].set_title("PACF — differenced")
plt.tight_layout()Compare the top row against the bottom. The top ACF crawls down — that is the signature of an undifferenced series, and it tells you nothing about \(p\) or \(q\). The bottom row is interpretable.
Your turn: using the bottom row and the four reference patterns from Lecture 2, propose a \((p, d, q)\).
In the room: say your proposed order out loud to a classmate and defend it in one sentence.
On your own: type it into a comment before you fit anything, in this form:
# (p, d, q) = (_, _, _) — PACF cuts off after lag ____, ACF cuts off after lag ____Naming the lag that decided each term is what turns a guess into a reading you can be wrong about — and being wrong about a specific lag is something you can learn from.
Either way the point is committing to an order before you see how it scores.
Walkthrough
Four panels. Top row is the raw series, bottom row is differenced.
Look at the top ACF. It crawls down — a long slow decay, bar after bar, staying well above the significance band for twenty-plus lags.
Here’s the trap. That pattern looks like enormously strong autocorrelation, and a reasonable person reads it as “there’s a lot of structure here to model.” It isn’t. That slow decay is the signature of an undifferenced series, and it appears no matter what the underlying process is. A random walk produces it. An AR(1) with a unit root produces it. White noise plus a trend produces it. It tells you nothing about p or q. It tells you only that you skipped a step.
Now the bottom row, after differencing. That is interpretable. The bars drop into the band quickly and you can actually read structure off it.
Try it
Pause the video. Using the bottom row and the four reference patterns from Lecture 2, propose a .
Write it down before you continue, and write down why — which panel, which lags. Then post it in the Week 2 thread. I want to see the range of proposals, because there genuinely is one, and that’s part of the lesson.
2.5 Step 3 — Fit it
order = (1, 1, 1) # <- REPLACE with your proposal from Step 2
model = ARIMA(train, order=order).fit()
print(model.summary().tables[0])
fc = model.forecast(N_TEST)
rmse = float(np.sqrt(np.mean((test.values - np.asarray(fc)) ** 2)))
print(f"\n52-week RMSE: {rmse:,.0f}")
Walkthrough
Line one of this cell is order = (1, 1, 1). Replace it with your proposal. The lab ships with (1,1,1) as a placeholder so the file renders for everyone — it is not the answer.
(Note for the narrator: with the shipped (1,1,1) the 52-week RMSE is 569. Say your own number aloud, then say the following.)
Your number will differ from mine if you chose a different order, and that is fine — that is the exercise. What you should be able to do is say why you chose what you chose.
2.6 Step 4 — Break it: the horizon trap
Now the same model, evaluated a different way.
errs = []
vals = y.values
for i in range(104, len(vals) - 1):
m = ARIMA(vals[:i+1], order=order).fit()
errs.append(vals[i+1] - float(m.forecast(1)[0]))
wf = float(np.sqrt(np.mean(np.array(errs) ** 2)))
print(f"walk-forward 1-step RMSE : {wf:,.0f} ({len(errs)} origins)")
print(f"single-split 52-step RMSE: {rmse:,.0f}")Notice how much smaller the walk-forward number is. It is tempting to quote it — it makes your model look better.
Do not. They answer different questions. One-step-ahead, knowing everything up to last week, is a fundamentally easier problem than forecasting a year out from a fixed origin. If a planner is ordering inventory twelve months ahead, the 52-step number is the honest one.
Walkthrough
Same model. Different evaluation. Watch what happens to the number.
walk-forward 1-step RMSE : 352 (172 origins)
single-split 52-step RMSE: 569Three fifty-two against five sixty-nine. The walk-forward number is nearly forty percent smaller, and it is computed from a hundred and seventy-two separate forecasts instead of one — so it even looks more rigorous. More origins, more evaluation, smaller error. It is enormously tempting to quote it. Do not.
They are not two estimates of the same thing. They answer different questions. One-step-ahead — knowing everything up to last week, forecasting next week — is a fundamentally easier problem than standing at a fixed origin and forecasting fifty-two weeks out. Of course the error is smaller. It would be alarming if it weren’t.
The question is which one matches the decision. If a planner is ordering inventory twelve months ahead, the 52-step number is the honest one, and quoting 352 to that person is not optimistic — it’s wrong.
This is the most common way I see forecast accuracy overstated in practice, and almost nobody does it dishonestly. They compute the number that’s easy to compute, and it happens to be the flattering one.
2.7 Step 5 — Granger causality, and the trap in it
Does the food aisle help predict the household aisle?
grangercausalitytests prints four test variants for every lag up to maxlag — a wall of output you did not ask for. The helper below silences it and pulls out the single F-test.
food = get_series("CA_1", "FOODS")
pair = pd.concat([y.rename("household"), food.rename("food")], axis=1).diff().dropna()
def granger(effect, cause, lag=4):
"""F-test for: do lags of `cause` help predict `effect`?"""
with redirect_stdout(io.StringIO()): # statsmodels prints 4 variants per lag
res = grangercausalitytests(pair[[effect, cause]], maxlag=lag)
F, p = res[lag][0]["ssr_ftest"][0], res[lag][0]["ssr_ftest"][1]
verdict = "significant" if p < 0.05 else "not significant"
print(f"{cause:>10} -> {effect:<10} F = {F:5.2f} p = {p:.4f} {verdict}")
granger("household", "food")Note the argument order: the second argument is the candidate cause. Getting this backwards is the single most common mistake with this function.
Your turn: run the other direction and report both F-statistics.
# ---- your turn: uncomment, then swap the two arguments ----
# granger(___, ___)
Walkthrough
Does the food aisle help predict the household aisle?
First, the practical note. grangercausalitytests prints four test variants for every lag up to maxlag — for lag 4, that’s a screen and a half of output nobody asked for. The helper swallows it with redirect_stdout and pulls out the single F-test we want. That’s what those imports were for.
food -> household F = 3.83 p = 0.0048 significantNow the argument order, and please write this down: the second argument is the candidate cause. granger("household", "food") asks whether food helps predict household. Getting this backwards is the single most common mistake with this function, and it is silent — you get a perfectly reasonable-looking result to the wrong question.
Try it
Your turn. Run the other direction and report both F-statistics. It’s one line — uncomment and swap the arguments.
household -> food F = 2.80 p = 0.0263 significantBoth directions are significant. Food predicts household, and household predicts food.
Sit with how strange that is. Household goods do not cause groceries. Groceries do not cause household goods. Nobody buys paper towels because they bought bananas. And yet the test says each helps predict the other, at conventional significance, in both directions. We’ll come back to what’s generating that. Hold the question.
2.8 Step 6 — The same test on macro data, and why the sample matters
Step 5 used two aisles of one store. Now the question from the lecture slide: does unemployment Granger-cause retail sales?
Run python scripts/prep_fred.py once before class — it caches three FRED series to data/processed/fred_monthly.csv, so this step needs no network and no API key.
fred = pd.read_csv(data_path("fred_monthly.csv"), parse_dates=["date"]).set_index("date")
# Both series are already stationary transforms:
# retail_growth = monthly % change in retail sales (log-differenced)
# unrate_diff = monthly change in the unemployment rate
macro = fred[["retail_growth", "unrate_diff"]].dropna()
def granger_macro(effect, cause, data, lag=4):
with redirect_stdout(io.StringIO()):
res = grangercausalitytests(data[[effect, cause]], maxlag=lag)
F, p = res[lag][0]["ssr_ftest"][0], res[lag][0]["ssr_ftest"][1]
print(f" {cause:>12} -> {effect:<14} F = {F:6.2f} p = {p:.4f}"
f" {'significant' if p < 0.05 else 'NOT significant'}")
print(f"Full sample: {macro.index.min():%Y-%m} to {macro.index.max():%Y-%m} ({len(macro)} months)")
granger_macro("retail_growth", "unrate_diff", macro)
granger_macro("unrate_diff", "retail_growth", macro)Both directions are significant, and unemployment looks like the stronger driver. That matches the story everyone expects: people lose jobs, people stop spending.
Now drop the COVID window and run it again.
ex_covid = macro[(macro.index < "2020-02-01") | (macro.index > "2021-06-01")]
print(f"Excluding Feb 2020 - Jun 2021 ({len(ex_covid)} months)")
granger_macro("retail_growth", "unrate_diff", ex_covid)
granger_macro("unrate_diff", "retail_growth", ex_covid)Your turn: restrict the sample to everything before February 2020 and run both directions once more.
# ---- your turn: uncomment and fill in the filter ----
# pre_covid = macro[___] # <- rows with index < "2020-02-01"
#
# granger_macro("retail_growth", "unrate_diff", pre_covid)
# granger_macro("unrate_diff", "retail_growth", pre_covid)
Walkthrough
Now the question from the lecture slide, on real macro data. Does unemployment Granger-cause retail sales?
Both series arrive already stationary — retail growth is log-differenced monthly percent change, unemployment is the monthly change in the rate. So we can test directly.
Full sample: 1992-02 to 2026-06 (412 months)
unrate_diff -> retail_growth F = 18.22 p = 0.0000 significant
retail_growth -> unrate_diff F = 12.80 p = 0.0000 significantBoth directions significant, and unemployment is the stronger of the two — F of eighteen against thirteen. That matches the story everybody expects. People lose jobs, people stop spending. You could write that paragraph without running anything, which is exactly why you should be suspicious of how satisfying it is.
Now drop the COVID window.
Excluding Feb 2020 - Jun 2021 (395 months)
unrate_diff -> retail_growth F = 2.91 p = 0.0214 significant
retail_growth -> unrate_diff F = 6.74 p = 0.0000 significantSeventeen months removed out of four hundred and twelve — four percent of the sample. The unemployment-to-retail F-statistic falls from 18.22 to 2.91. It’s still significant, barely, but it has lost about six-sevenths of its strength. The reverse direction went from 12.80 to 6.74. Weaker, but still overwhelming.
Try it
Your turn. Restrict to everything before February 2020 — no COVID period at all, no recovery either — and run both directions.
Pre-COVID, before Feb 2020 (336 months)
unrate_diff -> retail_growth F = 1.85 p = 0.1196 NOT significant
retail_growth -> unrate_diff F = 9.67 p = 0.0000 significantEighteen point two two. Two point nine one. One point eight five.
On twenty-eight years of pre-pandemic data, unemployment does not Granger-cause retail sales. The intuitive direction — the one everybody expects, the one that has a mechanism, the one you’d put in a slide without checking — is gone.
And the direction nobody leads with, retail predicting unemployment, survives every single cut: 12.80, then 6.74, then 9.67. Never above a p-value of point zero zero zero one.
2.9 Before you leave
You should have found, in Step 5, that the relationship runs both ways, at least in part.
Discuss (Step 5): Household goods do not cause groceries, and groceries do not cause household goods. So what is generating a two-way predictive relationship between two aisles of the same store?
And the question that matters commercially: given what you found, is there anything a planner is entitled to do with it? Name one legitimate use and one thing they must not conclude.
In Step 6 the intuitive direction — unemployment driving retail — is strong on the full sample, weak once COVID is removed, and gone entirely before 2020. The reverse direction survives every cut.
Discuss (Step 6): One analyst runs the full sample and reports that unemployment predicts retail sales. Another runs 1992–2019 and reports it does not. Neither has made an error.
What actually happened here — and what does it tell you about a Granger result computed on a sample containing one enormous shared shock?
Retail sales are a coincident-to-leading indicator and unemployment is a lagging one. Does that help explain which direction survived?
Solutions posted to Canvas after class.
Walkthrough
Two questions, and I want both in the Week 2 thread before Lecture 3.
Try it
First, the aisles. You found a two-way relationship between household goods and groceries. Neither causes the other. So what is generating it?
And the commercial half, which is the one that matters: given what you found, is there anything a planner is entitled to do with it? Name one legitimate use and one thing they must not conclude.
Try it
Second, the macro result. One analyst runs the full sample and reports that unemployment predicts retail sales. Another runs 1992 to 2019 and reports that it does not.
Neither of them has made an error. Both ran the correct test, correctly, on a defensible sample. That’s what makes this worth twenty minutes of your attention rather than being a cautionary tale about sloppy work.
What actually happened — and what does it tell you about a Granger result computed on a sample containing one enormous shared shock?
One hint worth using: retail sales are a coincident-to-leading indicator, and unemployment is a lagging one. Firms cut hours and stop hiring after demand falls, not before. Does that help explain which direction survived every cut?
Solutions for the coded parts go up on Canvas after the deadline. Lecture 3 next, where we leave retail data entirely — because it can’t show you what we need to see.
2.10 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.
Reproducible; verified across repeated runs. Step 3 excepted — that number depends on the order the student chooses.
| Quantity | Value |
|---|---|
| Series | CA_1 HOUSEHOLD, 277 weeks (train 225, test 52) |
| ADF levels | −1.33, p = 0.614, 12 lags → unit root |
| ADF first difference | −5.00, p < 0.001, 11 lags → stationary |
| ARIMA(1,1,1) 52-step RMSE | 569 (shipped placeholder order only) |
| Walk-forward 1-step RMSE | 352, over 172 origins |
| food → household | F = 3.83, p = 0.0048 |
| household → food | F = 2.80, p = 0.0263 |
| Macro, full sample (412 mo) | unrate→retail F = 18.22 · retail→unrate F = 12.80 |
| Ex-COVID (395 mo) | unrate→retail F = 2.91 (p = 0.0214) · retail→unrate F = 6.74 |
| Pre-COVID (336 mo) | unrate→retail F = 1.85 (p = 0.1196, not significant) · retail→unrate F = 9.67 |
Note for the instructor. The FRED full sample ends 2026-06, so it moves as prep_fred.py is re-run. The three F-statistics will drift slightly between terms; the pattern — collapse in one direction, survival in the other — is what’s stable, and it’s what the discussion questions rest on. Re-run this appendix before recording if the cache is more than a term old.