Lab 2: Stationarity, ARIMA, and Granger Causality

ECON 8310: Business Forecasting — University of Nebraska at Omaha

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).

---

Setup

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)}")

277 weeks — train 225, test 52

---

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")

levels                       ADF   -1.33   p  0.614   lags 12   UNIT ROOT

Now difference it and re-test.

dy = y.diff().dropna()
report_adf(dy, "first difference")

first difference             ADF   -5.00   p  0.000   lags 11   stationary

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.

---

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.

---

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}")

                               SARIMAX Results                                
==============================================================================
Dep. Variable:                  units   No. Observations:                  225
Model:                 ARIMA(1, 1, 1)   Log Likelihood               -1616.638
Date:                Mon, 19 Oct 2026   AIC                           3239.275
Time:                        16:24:56   BIC                           3249.510
Sample:                    01-31-2011   HQIC                          3243.407
                         - 05-18-2015                                         
Covariance Type:                  opg                                         
==============================================================================

52-week RMSE: 569

---

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}")

walk-forward 1-step RMSE : 352   (172 origins)
single-split 52-step RMSE: 569

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.

---

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")

      food -> household   F =  3.83   p = 0.0048   significant

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(___, ___)

---

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)

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   significant

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.

Dropping the rows first and then lagging would splice January 2020 straight onto July 2021, as if they were neighboring months. So build every lag on the unbroken series, then drop each row whose target or any of its lags falls inside the window.

import statsmodels.api as sm
from scipy import stats

def granger_excluding(effect, cause, data, start, end, lag=4):
    d = pd.DataFrame({"y": data[effect]})
    for k in range(1, lag + 1):                       # lags on the unbroken series
        d[f"y{k}"], d[f"x{k}"] = data[effect].shift(k), data[cause].shift(k)
    d = d.dropna()
    d = d[(d.index < start) | (d.index > pd.Timestamp(end) + pd.DateOffset(months=lag))]
    own = sm.OLS(d["y"], sm.add_constant(d[[f"y{k}" for k in range(1, lag + 1)]])).fit()
    both = sm.OLS(d["y"], sm.add_constant(d.drop(columns="y"))).fit()
    F = ((own.ssr - both.ssr) / lag) / (both.ssr / both.df_resid)
    p = stats.f.sf(F, lag, both.df_resid)
    print(f"  {cause:>12} -> {effect:<14} F = {F:6.2f}   p = {p:.4f}"
          f"   {'significant' if p < 0.05 else 'NOT significant'}")

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_excluding("retail_growth", "unrate_diff", macro, "2020-02-01", "2021-06-01")
granger_excluding("unrate_diff", "retail_growth", macro, "2020-02-01", "2021-06-01")

Excluding Feb 2020 - Jun 2021 (395 months)
   unrate_diff -> retail_growth  F =   2.42   p = 0.0480   significant
  retail_growth -> unrate_diff    F =   9.96   p = 0.0000   significant

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)

---

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.