Lab 9 Part 1: Recurrence — and When Gates Are Not Worth It

ECON 8310: Business Forecasting — University of Nebraska at Omaha

What we are doing

Labs 7 and 8 gave you two models that could not properly use a window. The feedforward network flattened it and lost the ordering; the CNN kept ordering but could only see seven weeks at a time.

A recurrent network has neither limit. It reads the sequence in order and carries a hidden state, so in principle nothing stops it using the whole window.

Today you build two of them — a plain RNN and an LSTM — and find that at the training budget anyone would have used, the more sophisticated one loses badly. Then you vary the budget and the result changes. Working out which claim the evidence actually supports, and why, is the point of the lab, and it is the most transferable thing in it.

Suggested time: ~40 minutes · You will need: pandas, numpy, matplotlib, torch, and data/processed/m5_weekly.csv. Most fits are quick — about six seconds each. The exception is Step 4, which trains four models to 100 epochs and takes two to three minutes; start it and read the text underneath while it runs.


Setup

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import torch
import torch.nn as nn
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.")

weekly = (pd.read_csv(data_path("m5_weekly.csv"), parse_dates=["week_start"])
            .sort_values(["store_id", "cat_id", "week_start"]))
CHANNELS = ["units", "snap_days", "event_days", "avg_price"]
CUTOFF   = np.datetime64("2015-05-17")
rmse = lambda a, b: float(np.sqrt(np.mean((np.asarray(a) - np.asarray(b)) ** 2)))

def build(W):
    X, Y, D = [], [], []
    for _, g in weekly.groupby(["store_id", "cat_id"]):
        arr, dates = g[CHANNELS].values.astype("float32"), g.week_start.values
        for i in range(len(g) - W):
            X.append(arr[i:i+W]); Y.append(arr[i+W, 0]); D.append(dates[i+W])
    X = np.stack(X); Y = np.array(Y, dtype="float32"); D = np.array(D)
    tr, te = D < CUTOFF, D >= CUTOFF
    mu, sd = X[tr].mean((0,1), keepdims=True), X[tr].std((0,1), keepdims=True) + 1e-8
    ym, ys = Y[tr].mean(), Y[tr].std()
    Xn, Yn = (X - mu)/sd, (Y - ym)/ys
    # recurrent layers want (batch, time, channels) -- NOT the Conv1d layout
    return (torch.tensor(Xn[tr]), torch.tensor(Yn[tr]), torch.tensor(Xn[te]), Y[te], ym, ys)

DATA = {W: build(W) for W in (13, 26, 52)}
print("shape at W=26:", DATA[26][0].shape, "(batch, time, channels)")
shape at W=26: torch.Size([5940, 26, 4]) (batch, time, channels)
The layout flips back

Conv1d wanted (batch, channels, time). nn.RNN and nn.LSTM with batch_first=True want (batch, time, channels) — the transpose you added in Lab 8 must come out again here.

Neither library will complain. Both will train on nonsense.


Step 1 — A plain RNN

One recurrent layer, then a linear layer on the final hidden state — the summary after reading the whole window.

class Recurrent(nn.Module):
    def __init__(self, kind, n_channels=4, hidden=64):
        super().__init__()
        layer = nn.RNN if kind == "rnn" else nn.LSTM
        self.rnn = layer(n_channels, hidden, batch_first=True)
        self.fc  = nn.Linear(hidden, 1)
    def forward(self, x):
        out, _ = self.rnn(x)          # nn.LSTM returns out, (h, c) -- the _ absorbs both
        return self.fc(out[:, -1, :]).squeeze(-1)

def train(kind, W, seed, epochs=30):
    Xtr, Ytr, Xte, Yte, ym, ys = DATA[W]
    torch.manual_seed(seed)
    model = Recurrent(kind)
    opt = torch.optim.Adam(model.parameters(), lr=1e-3, weight_decay=1e-4)
    lossfn = nn.MSELoss()
    loader = torch.utils.data.DataLoader(
        torch.utils.data.TensorDataset(Xtr, Ytr), batch_size=128, shuffle=False)
    for _ in range(epochs):
        model.train()
        for xb, yb in loader:
            opt.zero_grad()
            lossfn(model(xb), yb).backward()
            torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0)   # see below
            opt.step()
    model.eval()
    with torch.no_grad():
        score = rmse(Yte, model(Xte).numpy() * ys + ym)
    return score, sum(p.numel() for p in model.parameters())

score, n_par = train("rnn", 26, seed=1)
print(f"vanilla RNN — {n_par:,} parameters, test RMSE {score:,.0f}")
vanilla RNN — 4,545 parameters, test RMSE 860

clip_grad_norm_ is not decoration. Backpropagating through 26 steps multiplies gradients repeatedly, and the product can explode. Clipping caps the update size. It is one line, it costs nothing, and recurrent models are the place you need it.


Step 2 — An LSTM

Identical code, one word changed.

score_lstm, n_par_lstm = train("lstm", 26, seed=1)
print(f"LSTM        — {n_par_lstm:,} parameters, test RMSE {score_lstm:,.0f}")
print(f"\nparameter ratio: {n_par_lstm / n_par:.1f}x the RNN")
LSTM        — 17,985 parameters, test RMSE 1,004

parameter ratio: 4.0x the RNN

Account for that ratio before reading the RMSE. An LSTM has a forget gate, an input gate, an output gate and a candidate — four sets of weights where the RNN has one. Same hidden size, four times the parameters.


Step 3 — Which wins, and does it hold?

One comparison at one window is not evidence. Vary the window and use two seeds.

rows = []
for W in (13, 26, 52):
    r = [train("rnn",  W, s)[0] for s in (1, 2)]
    l = [train("lstm", W, s)[0] for s in (1, 2)]
    rows.append({"window": W,
                 "RNN": round(np.mean(r)),  "RNN_sd": round(np.std(r)),
                 "LSTM": round(np.mean(l)), "LSTM_sd": round(np.std(l)),
                 "winner": "RNN" if np.mean(r) < np.mean(l) else "LSTM"})
print(pd.DataFrame(rows).to_string(index=False))
 window  RNN  RNN_sd  LSTM  LSTM_sd winner
     13  878      33   990       21    RNN
     26  842      19   987       17    RNN
     52  805       5   998        3    RNN
Read the columns in this order

The winner column first. At this training budget, on this panel, the plain RNN wins at every window length. That part is not a seed artifact — the run-to-run spreads (3 to 33) are tiny next to the gaps (112 to 193).

Then the parameter count. The model that lost has four times as many parameters as the model that won.

Then be suspicious. Both models trained for exactly 30 epochs, because that is what train() defaults to. Nothing above establishes that 30 epochs is enough for both of them — and a model with four times the parameters is not obviously finished learning at the same point as the smaller one. Step 4 checks, and the answer is not the one this table implies.

Your turn: compute the gap at each window and say whether it is closing or widening.

# ---- your turn: uncomment and fill in ----
# df = pd.DataFrame(rows)
# df["gap"] = ___ - ___          # <- which two columns?
# print(df[["window", "gap"]].to_string(index=False))

Step 4 — Was that a fair fight?

Step 3 compared two models at one training budget. Lab 7 taught you not to trust a comparison run at one seed; this is the same objection one level up. An LSTM has four times the parameters, and bigger models generally take longer to converge.

So train both to 100 epochs, and record the test RMSE along the way instead of only at the end.

CKPT = [10, 20, 30, 45, 60, 80, 100]

def curve(kind, W, seed):
    """Train once to 100 epochs, scoring the test set at each checkpoint."""
    Xtr, Ytr, Xte, Yte, ym, ys = DATA[W]
    torch.manual_seed(seed)
    model = Recurrent(kind)
    opt = torch.optim.Adam(model.parameters(), lr=1e-3, weight_decay=1e-4)
    lossfn = nn.MSELoss()
    loader = torch.utils.data.DataLoader(
        torch.utils.data.TensorDataset(Xtr, Ytr), batch_size=128, shuffle=False)
    out = {}
    for ep in range(1, max(CKPT) + 1):
        model.train()
        for xb, yb in loader:
            opt.zero_grad(); lossfn(model(xb), yb).backward()
            torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0); opt.step()
        if ep in CKPT:
            model.eval()
            with torch.no_grad():
                out[ep] = rmse(Yte, model(Xte).numpy() * ys + ym)
    return out

paths = {}
for kind in ("rnn", "lstm"):
    runs = [curve(kind, 26, s) for s in (1, 2)]        # fit each seed ONCE
    paths[kind] = [np.mean([r[e] for r in runs]) for e in CKPT]

table = pd.DataFrame({"epochs": CKPT,
                      "RNN":  [round(v) for v in paths["rnn"]],
                      "LSTM": [round(v) for v in paths["lstm"]]})
table["gap"] = table.LSTM - table.RNN
print(table.to_string(index=False))
 epochs  RNN  LSTM  gap
     10 1099  1002  -97
     20  922  1033  111
     30  842   987  145
     45  834   953  119
     60  828   926   98
     80  833   885   52
    100  830   821   -9
fig, ax = plt.subplots(figsize=(6.5, 3.4))
ax.plot(CKPT, paths["rnn"],  marker="o", label="vanilla RNN")
ax.plot(CKPT, paths["lstm"], marker="o", label="LSTM")
ax.axvline(30, ls="--", lw=1, color="grey")
ax.set_xlabel("epochs"); ax.set_ylabel("test RMSE"); ax.legend()
ax.set_title("Step 3 measured at the dashed line", fontsize=9)
plt.tight_layout()

The gap is largest at exactly the budget we used

Read the gap column downwards. It is negative at 10 epochs — the LSTM is ahead early — then peaks at 30, the setting Step 3 used, and shrinks at every checkpoint after until the two models are level.

Look at the shapes, not just the endpoint. The RNN is finished by about epoch 45 and flat thereafter. The LSTM is still improving at 100 epochs; its curve has not turned over. We ran out of patience, not the model.

So the honest verdict is not “the LSTM loses.” It is:

  • at a fixed 30-epoch budget the RNN is far ahead — and 30 epochs is the point of its maximum advantage anywhere in this range;
  • given enough epochs the two are tied;
  • the RNN reaches its best score about three times faster, on a quarter of the parameters.

That last line is a real advantage and worth having. It is simply a different claim from the one the Step 3 table appears to make — and you could not have told the two apart without varying the budget.

Your turn: the run-to-run spread in Step 3 was roughly 20 RMSE. From which checkpoint on is the gap smaller than that?

# ---- your turn: uncomment and fill in ----
# noise = ___                       # <- the seed spread from Step 3
# print(table[table.gap.abs() < noise].epochs.tolist())

Step 5 — Why the gates buy nothing here

The LSTM’s gates exist to solve one problem: in a plain RNN, gradients shrink as they propagate backward, so signal from far in the past never reaches the weights. Put a number on it.

for factor in (0.7, 0.9):
    print(f"gradient shrink factor {factor} per step:")
    for steps in (1, 5, 10, 26, 52):
        print(f"   {steps:>3} steps back -> {factor**steps:,.6f}")
    print()
gradient shrink factor 0.7 per step:
     1 steps back -> 0.700000
     5 steps back -> 0.168070
    10 steps back -> 0.028248
    26 steps back -> 0.000094
    52 steps back -> 0.000000

gradient shrink factor 0.9 per step:
     1 steps back -> 0.900000
     5 steps back -> 0.590490
    10 steps back -> 0.348678
    26 steps back -> 0.064611
    52 steps back -> 0.004175

At 0.7 per step the signal is essentially gone by 26 steps. That is the problem gating fixes.

Two things about that table before we go on. At a shrink factor of 0.7 the signal is gone by 26 steps. At 0.9 it is still 0.065 at 26 steps and 0.004 at 52 — small, but not zero. Vanishing gradients are a matter of degree and of sequence length, and 26 weeks is short.

Now stop asserting and go measure whether this dataset has the problem gating solves. Three questions, three numbers.

w = weekly.assign(series=weekly.store_id + "_" + weekly.cat_id)
w["level"] = w.groupby("series").units.transform("mean")
w["dev"]   = w.units - w.level

grand = w.units.mean()
btw = ((w.level - grand) ** 2).sum() / ((w.units - grand) ** 2).sum()
print(f"1. between-series share of total variance : {btw:.1%}")

print("\n2. within-series autocorrelation of units")
for lag in (1, 4, 13, 26, 52):
    a = w.groupby("series").dev.apply(lambda x: x.autocorr(lag)).mean()
    print(f"      lag {lag:>2}: {a:+.3f}")

print("\n3. does a flag at t predict units at t+k?")
for ch in ("snap_days", "event_days"):
    cells = []
    for k in (0, 1, 2, 4):
        c = w.groupby("series").apply(lambda g: g[ch].corr(g.dev.shift(-k))).mean()
        cells.append(f"t+{k}: {c:+.3f}")
    print(f"      {ch:<11} " + "   ".join(cells))
1. between-series share of total variance : 91.0%

2. within-series autocorrelation of units
      lag  1: +0.828
      lag  4: +0.824
      lag 13: +0.695
      lag 26: +0.575
      lag 52: +0.432

3. does a flag at t predict units at t+k?
      snap_days   t+0: +0.143   t+1: +0.030   t+2: -0.128   t+4: +0.116
      event_days  t+0: +0.004   t+1: +0.023   t+2: -0.009   t+4: -0.009
Three reasons the gates have nothing to do

1. Most of the variance is not sequential at all. About 91% of the variation in weekly units is between series — which store, which category — and only 9% is within a series over time. The single biggest thing a model must get right is the level of the series in front of it, and that is readable from the most recent few weeks. No memory mechanism is needed for it.

2. The series are smooth, which is the opposite of needing memory. Within-series autocorrelation is 0.83 at lag 1 and still 0.43 at lag 52. That looks like long memory and is in fact the reason long memory is unnecessary: because neighboring weeks are nearly the same number, the recent past is a sufficient statistic for the distant past. Nothing has to be carried across the gap, because the gap has been filled in continuously along the way. Gating pays off when a value must survive many steps that would otherwise overwrite it — and here nothing overwrites anything.

3. Nothing in the inputs arrives early. snap_days correlates +0.14 with units in the same week. The −0.13 at t+2 and +0.12 at t+4 are not advance notice: SNAP days fall on the same dates every month, so this week’s snap_days predicts the count two and four weeks out (correlations of −0.84 and +0.82), and those weeks carry their own same-week effect. The lead correlations are the same-week +0.14 echoed by the calendar. event_days is flat at every horizon including zero. There is no announcement structure in this data — nothing shows up at time \(t\) that must be stored until \(t+k\). That is the canonical use for a forget gate and an input gate, and this panel does not contain one instance of it.

So the gates solve something real, on data that does not have it. They charge four times the parameters and — as Step 4 showed — roughly three times the epochs for the service, and by 100 epochs they return a tie.

The transferable part is the three questions, not the answer. Before reaching for a gated or attention-based model on your own series, ask: how much of my variance is cross-sectional rather than sequential? Is my series smooth enough that recent values already summarize old ones? Does anything in my inputs lead the target? If the answers look like this panel’s, the extra machinery has nothing to do.


Before you leave

Discuss (the result): Using the three diagnostics from Step 5, describe a forecasting problem where you would expect the LSTM to beat the RNN outright. Be specific: what would the between-series share look like, what would the autocorrelation look like, and what would have to be true of the exogenous inputs?

Discuss (the budget): Step 3 and Step 4 support different claims from the same two models on the same data. State the strongest claim you would be willing to publish, in one sentence — and say what a reader would need to be told alongside it for the sentence to be honest.

Discuss (the practice): You have now been caught twice: once by seeds in Lab 7, once by the training budget here. Both times a comparison looked settled and was not. What is the general form of that mistake, and what would you check the next time you see two models compared in a table?

Discuss (for Homework 5): Part 2 asks you to fit an RNN, an LSTM and a Transformer, and to explain where the extra LSTM parameters go and whether they bought anything. You now have the RNN and LSTM half. What is your prediction for the Transformer, and what would make you wrong?


Solutions posted to Canvas after class.