12  Lab 11 Part 1: Structural Time Series — and Two Ways to Get the Interval Wrong

Run this one live: the executable lab.

12.1 What we are doing

A structural model writes a series as parts you can name — a level that drifts, a seasonal shape, and noise — and puts a prior on each. The output is not a forecast line but a distribution for every future week.

Today you will build one, hit a failure the diagnostics catch, fix it — and then find that the correct way to forecast produces a worse-calibrated interval than the shortcut.

That last result is not a mistake in the lab. It is the model telling you something, and working out what takes the rest of the session.

Suggested time: ~40 minutes · You will need: pandas, numpy, matplotlib, pymc, arviz, and data/processed/m5_weekly.csv. Sampling takes about 20 seconds per fit.

---

12.2 Setup

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import arviz as az
import warnings, logging
warnings.filterwarnings('ignore')
logging.getLogger('pymc').setLevel(logging.ERROR)

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"])
ca1 = (weekly[(weekly.store_id == "CA_1") & (weekly.cat_id == "FOODS")]
         .sort_values("week_start").reset_index(drop=True))

N_TEST = 52
train, test = ca1.iloc[:-N_TEST], ca1.iloc[-N_TEST:]
y = train.units.values.astype(float)
n = len(train)
y_mu, y_sd = y.mean(), y.std()
y_z = (y - y_mu) / y_sd          # standardize: the priors below assume unit scale

K, P = 3, 52                     # 3 Fourier pairs, annual period on weekly data
def fourier(t):
    return np.column_stack([f(2*np.pi*k*t/P) for k in range(1, K+1) for f in (np.sin, np.cos)])
t_all = np.arange(len(ca1))
X_tr, X_te = fourier(t_all[:n]), fourier(t_all[n:])

print(f"train {n} weeks, test {N_TEST} weeks | {2*K} Fourier coefficients standing in for 52")

---

Walkthrough

train 225 weeks, test 52 weeks | 6 Fourier coefficients standing in for 52

CA_1 FOODS again — the series you have been forecasting since Homework 2 — held out for a full final year.

One line deserves attention. Six coefficients standing in for fifty-two. Three sine-cosine pairs describe the annual shape, rather than fifty-two weekly dummies. That is a modeling choice with a cost and a benefit: you cannot represent a one-week spike, but you have six parameters to estimate instead of fifty-two, on 225 observations. Remember it — it comes back in Step 3 as a suspect.

Also note the standardization. y_z is the series in standard units, because the priors below — Exponential(10), Normal(0, 0.5) — assume unit scale. Priors are not scale-free, and this is the line that makes them defensible.

12.3 Step 1 — Fit it the obvious way

The model is a drifting level plus seasonality. PyMC has a distribution for a random walk, so the natural thing is to use it.

with pm.Model() as centered:
    sigma_level = pm.Exponential("sigma_level", 10.0)
    level = pm.GaussianRandomWalk("level", sigma=sigma_level, shape=n)
    beta  = pm.Normal("beta", 0, 0.5, shape=2*K)
    sigma = pm.Exponential("sigma", 1.0)
    pm.Normal("obs", mu=level + pm.math.dot(X_tr, beta), sigma=sigma, observed=y_z)
    idata_c = pm.sample(1500, tune=2000, chains=4, progressbar=False,
                        random_seed=42, target_accept=0.95)

s_c = az.summary(idata_c, var_names=["sigma_level", "sigma", "beta"],
                 ci_kind="hdi", ci_prob=0.94)
print(f"R-hat        {s_c.r_hat.max():.3f}   (want < 1.01)")
print(f"ESS (bulk)   {s_c.ess_bulk.min():.0f}     (want > 400)")
print(f"divergences  {int(idata_c.sample_stats.diverging.sum())}")

It did not converge. No exception, a forecast is available, and the diagnostics are the only thing telling you not to use it. Note that target_accept is already at 0.95 — this is not fixed by asking the sampler to try harder.

---

Walkthrough

The model is a drifting level plus seasonality. PyMC has a distribution for a random walk, so the obvious thing is to reach for it: pm.GaussianRandomWalk. It reads exactly like the maths.

R-hat        1.022   (want < 1.01)
ESS (bulk)   274     (want > 400)
divergences  0

It did not converge, and I want you to look carefully at how you know.

R-hat is 1.022, above the 1.01 threshold — the four chains do not agree with each other. Effective sample size is 274, below 400 — the draws are so autocorrelated that six thousand samples carry less information than a few hundred independent ones. And divergences are zero.

That last line matters, because divergences are the diagnostic everyone learns first, and here they are silent. The sampler never had to give up on a step; it just wandered badly. Different pathologies trip different instruments — which is precisely why you check all three, every time, rather than the one you remember.

Next week you will see the mirror image: a model where R-hat and ESS both pass and forty-three divergences are the only thing objecting.

And note this is not fixed by asking the sampler to try harder. target_accept is already at 0.95. The problem is not effort.

12.4 Step 2 — Reparameterize

The problem is the shape of the posterior, not the model. Write the walk as a cumulative sum of standard Normal draws scaled by sigma_level, and the sampler gets a geometry it can navigate.

with pm.Model() as model:
    sigma_level = pm.Exponential("sigma_level", 10.0)
    z = pm.Normal("z", 0, 1, shape=n)                              # standard normal
    level = pm.Deterministic("level", pm.math.cumsum(z * sigma_level))
    beta  = pm.Normal("beta", 0, 0.5, shape=2*K)
    sigma = pm.Exponential("sigma", 1.0)
    pm.Normal("obs", mu=level + pm.math.dot(X_tr, beta), sigma=sigma, observed=y_z)
    idata = pm.sample(1500, tune=2000, chains=4, progressbar=False,
                      random_seed=42, target_accept=0.95)

s = az.summary(idata, var_names=["sigma_level", "sigma", "beta"], ci_kind="hdi", ci_prob=0.94)
print(f"R-hat        {s.r_hat.max():.3f}   (want < 1.01)")
print(f"ESS (bulk)   {s.ess_bulk.min():.0f}    (want > 400)")
print(f"divergences  {int(idata.sample_stats.diverging.sum())}")

Same model, same data, same priors. cumsum(z * sigma) and GaussianRandomWalk(sigma=sigma) describe the identical distribution — only the coordinates the sampler explores have changed.

Look at the components it recovered:

lev = idata.posterior["level"].values.reshape(-1, n).mean(0) * y_sd + y_mu
sea = (idata.posterior["beta"].values.reshape(-1, 2*K) @ X_tr.T).mean(0) * y_sd

fig, ax = plt.subplots(2, 1, figsize=(10, 4.5), sharex=True)
ax[0].plot(train.week_start, train.units, lw=.8, alpha=.5, label="observed")
ax[0].plot(train.week_start, lev, lw=2, label="level")
ax[0].legend(); ax[0].set_ylabel("units")
ax[1].plot(train.week_start, sea, lw=1.5, color="darkgreen")
ax[1].set_ylabel("seasonal"); ax[1].axhline(0, color="k", lw=.5)
plt.tight_layout()

---

Walkthrough

The problem is the shape of the posterior, not the model.

A centered random walk couples the level to sigma_level: when the scale is small, the steps must be small, so the posterior narrows into a funnel that the sampler cannot traverse at a single step size. The fix is to write the walk as a cumulative sum of standard Normal draws, scaled afterwards. Now z and sigma_level are independent in the posterior geometry.

R-hat        1.002   (want < 1.01)
ESS (bulk)   2376    (want > 400)
divergences  0

R-hat 1.022 to 1.002. Effective sample size 274 to 2,376 — nearly nine times more information out of the same number of draws.

And here is the sentence to say slowly: same model, same data, same priors. cumsum(z * sigma) and GaussianRandomWalk(sigma=sigma) describe the identical distribution. Nothing about what we believe has changed. Only the coordinate system the sampler explores.

That is worth internalising, because it is unlike everything else in this course. This is not a modeling improvement or a bias-variance trade. It is a computational fix to a computational problem, and the give-away is that the posterior it targets is unchanged.

And now the payoff of a structural model — you can look at the parts. Top panel, the level: a slow drift under the noise. Bottom panel, the seasonal shape the six Fourier coefficients recovered. Neither of those is available from XGBoost at any price.

12.5 Step 3 — Forecast, two ways

To forecast you need the level at each future week. There are two candidates, and one of them is theoretically right.

Held constant: carry the last estimated level forward unchanged. Simple, and it asserts the level stops drifting the moment your data ends — which is not what the model says.

Propagated: let the level keep random-walking forward, exactly as it did in-sample. This is what the fitted model actually implies.

rng = np.random.default_rng(42)
post = idata.posterior
lev_last = post["level"].values[:, :, -1].ravel()
sl       = post["sigma_level"].values.ravel()
b        = post["beta"].values.reshape(-1, 2*K)
sg       = post["sigma"].values.ravel()
D = len(lev_last)
actual = test.units.values

# The POINT forecast is identical for both methods, and not approximately so:
# propagating the level adds zero-mean steps, so E[level] is unchanged. Only the
# spread differs. Compute the point forecast once, from the deterministic part.
point = (lev_last[:, None] + b @ X_te.T).mean(0) * y_sd + y_mu
rmse_point = float(np.sqrt(np.mean((actual - point) ** 2)))

def score(draws, label):
    lo, hi = np.percentile(draws, [3, 97], axis=0) * y_sd + y_mu
    return {"method": label,
            "RMSE": round(rmse_point),
            "coverage_%": round(100*float(np.mean((actual >= lo) & (actual <= hi)))),
            "mean_width": round(float(np.mean(hi - lo)))}

# (a) level held constant at its last value  -- the tempting shortcut
mu_a = lev_last[:, None] + b @ X_te.T
draws_a = mu_a + rng.normal(0, 1, size=mu_a.shape) * sg[:, None]

# (b) level keeps random-walking into the future -- what the model actually says
steps = rng.normal(0, 1, size=(D, N_TEST)) * sl[:, None]
mu_b = lev_last[:, None] + np.cumsum(steps, axis=1) + b @ X_te.T
draws_b = mu_b + rng.normal(0, 1, size=mu_b.shape) * sg[:, None]

print(pd.DataFrame([score(draws_a, "level held constant"),
                    score(draws_b, "level propagated")]).to_string(index=False))

The theoretically correct choice calibrated worse. Work out why.

The two methods have the same point forecast — not nearly the same, exactly the same, because propagating the level adds steps with mean zero. So RMSE cannot distinguish them, and everything that separates them lives in the interval — which differs enormously, and not in the direction you would predict.

Propagating the walk is the defensible choice: the level genuinely does keep drifting, and freezing it asserts a certainty the model does not have. Yet the propagated interval over-covers, while the frozen one lands close to its nominal 94%.

Over-covering means too wide. Over 52 weeks the propagated level accumulates roughly \(\sqrt{52}\,\sigma_{\text{level}}\) of drift, so if \(\sigma_{\text{level}}\) is too large the interval inflates fast.

And there is good reason to think it is too large. The random walk is the most flexible part of this model. Anything the Fourier terms fail to capture — a promotion, a holiday shift, an outlier — gets absorbed into the level, inflating \(\sigma_{\text{level}}\) beyond genuine drift.

So the frozen version is not right. It is wrong in a way that happens to cancel an overestimated drift on this series and this horizon. Rely on that and it will fail elsewhere.

Your turn: check the size of the estimated drift against the noise it is competing with.

# ---- your turn: uncomment and fill in ----
# print(az.summary(idata, var_names=[___, ___],        # <- the two scale parameters
#                  ci_kind="hdi", ci_prob=0.94).round(4).to_string())
# print("\ndrift accumulated over 52 weeks, in standardized units:")
# print(f"  {np.sqrt(52) * float(post['sigma_level'].values.mean()):.2f}")

---

Now split coverage by horizon, which is where the difference actually shows.

def by_horizon(draws, label):
    lo, hi = np.percentile(draws, [3, 97], axis=0) * y_sd + y_mu
    inside = (actual >= lo) & (actual <= hi)
    segs = {"wk 1-13": slice(0,13), "wk 14-39": slice(13,39), "wk 40-52": slice(39,52)}
    row = {"method": label}
    for name, sl_ in segs.items():
        row[name] = f"{100*inside[sl_].mean():.0f}% / {int((hi-lo)[sl_].mean()):,}"
    return row

print("coverage / mean width, by horizon (nominal 94%)")
print(pd.DataFrame([by_horizon(draws_a, "level held constant"),
                    by_horizon(draws_b, "level propagated")]).to_string(index=False))

Read the widths across each row. The frozen version’s interval is essentially the same width at week 50 as at week 1 — it has no mechanism to widen. The propagated version’s roughly doubles.

That is the real distinction, and it is structural rather than a matter of which scored better here: one of these models knows that forecasting further ahead is harder, and the other does not. On this series, at this horizon, that ignorance happened not to cost anything. On a series with a genuinely drifting level it would.

---

Walkthrough

To forecast you need the level at each future week, and there are two candidates.

Held constant. Carry the last estimated level forward unchanged. Simple — and it asserts that the level stops drifting the moment your data ends, which is not what the model says.

Propagated. Let the level keep random-walking forward exactly as it did in-sample. This is what the fitted model actually implies. It is the defensible choice.

             method  RMSE  coverage_%  mean_width
level held constant  1589          96        6340
   level propagated  1589         100       11667

The point forecasts are identical — 1,589 on both rows, and not by coincidence. Propagating the level adds steps whose mean is zero, so the expected level is unchanged and RMSE simply cannot tell these two apart. Everything separating them is in the interval, which differs by a factor of nearly two, and not in the direction you would predict.

The theoretically correct method over-covers: 100% of weeks fell inside a 94% interval. The shortcut lands at 96%, close to nominal. The defensible choice calibrated worse.

Try it

Pause and work out why before I tell you. Then do the first “your turn” — print the two scale parameters and the drift accumulated over 52 weeks.

Over 52 weeks a random walk accumulates about \(\sqrt{52}\) — call it seven — times sigma_level of drift. So if sigma_level is even somewhat too large, the interval inflates fast, and it inflates with the square root of the horizon.

And there is good reason to think it is too large, which brings back the Fourier choice from the setup. The random walk is the most flexible component in this model. Anything the six coefficients cannot represent — a promotion, a holiday that moves, a single bad week — has nowhere else to go, so it gets absorbed into the level. That inflates sigma_level beyond genuine drift, and the forecast then propagates that inflated drift for a year.

So — and this is the sentence that matters — the frozen version is not right. It is wrong in a way that happens to cancel an overestimated drift, on this series, at this horizon. Rely on that and it will fail somewhere else, silently, and you will not know which of the two errors moved.

coverage / mean width, by horizon (nominal 94%)
             method      wk 1-13      wk 14-39      wk 40-52
level held constant 100% / 6,255   92% / 6,537  100% / 6,027
   level propagated 100% / 8,081 100% / 11,974 100% / 14,637

Read the widths across each row, not the coverage.

The frozen version: 6,255 at weeks 1–13, and 6,027 at weeks 40–52. It is narrower a year out than a month out. It has no mechanism to widen, because nothing in it knows that time has passed. The propagated version: 8,081, then 11,974, then 14,637. It nearly doubles.

One of these models knows that forecasting further ahead is harder, and the other does not. That is a structural difference, not a scoring difference, and it is why the aggregate table on its own would have misled you. On this series, at this horizon, the ignorance happened to cost nothing. On a series with a genuinely drifting level it would cost a great deal.

And notice the only sub-nominal cell in the table: the frozen version drops to 92% in weeks 14–39. Its respectable 96% average is hiding a stretch where it is overconfident. Aggregates hide things; that is the third time this course has said so.

12.6 Step 4 — Calibrated is not the same as useful

Take the propagated version — the one the model implies — and ask two separate questions of it.

res = score(draws_b, "level propagated")
print(f"coverage      {res['coverage_%']}%   against a nominal 94%")
print(f"mean width    {res['mean_width']:,} units")
print(f"              = {100*res['mean_width']/ca1.units.mean():.0f}% of average weekly units")
print(f"RMSE          {res['RMSE']:,}")
print(f"seasonal naive{np.sqrt(np.mean((actual - train.units.values[-52:])**2)):>7,.0f}")

Question one: is it calibrated? Compare coverage against 94%.

Question two: is it useful? Compare the width against the size of the thing you are forecasting, and ask whether you would place an inventory order against it.

Your turn: state, in one sentence each, your answer to both — and whether passing the first excuses failing the second.

# ---- your turn: uncomment and answer in the strings ----
# calibrated = "..."
# useful     = "..."
# print(calibrated); print(useful)

---

Walkthrough

Take the propagated version — the one the model implies — and ask it two separate questions.

coverage      100%   against a nominal 94%
mean width    11,667 units
              = 59% of average weekly units
RMSE          1,589
seasonal naive  1,660

Question one: is it calibrated? Yes, and then some. 100% against a nominal 94% — not overconfident. Conservative. Question two: is it useful?

The interval is 11,667 units wide on a series averaging about 19,700 — call it plus or minus thirty percent. Now say that out loud as a business sentence: “I am 94% confident that next week’s demand is somewhere between roughly thirty percent below and thirty percent above average.”

Nobody places an inventory order against that. It is honest and it is useless, and those are not the same axis.

Try it

Write your one-sentence answer to each question, and then the third question the lab actually asks: does passing the first excuse failing the second?

And the point forecast: 1,589 against seasonal naive’s 1,660. About 4% better, for a model that took twenty seconds to sample, a reparameterization to converge, and a page of code.

If accuracy were the argument for this model, that number would be an embarrassment. Accuracy is not the argument. What you bought is a level and a seasonal component you can plot and defend, a full distribution for every future week, and — the thing the whole lab has been about — an interval that admits when it does not know. Whether that trade is worth 5% is a business question, and it depends entirely on whether anyone downstream will use the uncertainty you went to the trouble of producing.

12.7 Before you leave

Discuss (the two problems): The convergence failure was caught by a diagnostic. The over-wide interval was not — it needed a held-out set and a coverage calculation. What is the structural difference between the two, and what habit catches the second kind?

Discuss (the payoff): The point forecast beat seasonal naive by only a few percent. Given that, what did fitting this model actually buy — and name a business question it can answer that XGBoost from Homework 4 cannot.

Discuss (the fix): The interval is honest and too wide to act on. Propose one change to the model — not the forecasting code — that would narrow it, and say what you would be assuming in exchange.

---

Solutions posted to Canvas after class.

Walkthrough

Try it

Three questions for the board.

The two problems. The convergence failure was caught by a diagnostic — R-hat and ESS, computed from the fit alone, before you looked at any data you were predicting. The over-wide interval was caught by nothing of the sort: it needed a held-out year and a coverage calculation. What is the structural difference between those two kinds of error, and what habit catches the second kind?

held-out check asks “is the model any good?” — and no amount of the first substitutes for the second. A perfectly converged sampler will happily give you a useless interval.*

The payoff. The point forecast beat seasonal naive by about four percent. Given that, what did fitting this model actually buy? Name one business question it answers that XGBoost from Homework 4 cannot.

The fix. The interval is honest and too wide to act on. Propose one change to the model — not the forecasting code — that would narrow it, and say what you are assuming in exchange.

assumes the level really does drift slowly; or more seasonal flexibility so the Fourier terms stop dumping their residual into the level. The second is the more interesting answer because it addresses the cause rather than the symptom.*

Solutions for the coded parts go up on Canvas after the deadline. Next week: ten series instead of one, and what happens when some of them barely have any data.

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

random_seed=42 on both fits, default_rng(42) on the forecast draws. Two fits, ~20 s each.

Quantity	Value
Split	train 225 weeks, test 52 · 6 Fourier coefficients (3 pairs) for a 52-week period
Centered model	R-hat 1.022 · ESS 274 · divergences 0
Non-centered model	R-hat 1.002 · ESS 2,376 · divergences 0
ESS improvement	274 → 2,376, about 8.7×, from a reparameterization that changes no belief
Level held constant	RMSE 1,589 · coverage 96% · mean width 6,340
Level propagated	RMSE 1,589 (identical by construction) · coverage 100% · mean width 11,667
By horizon — frozen	100% / 6,255 · 92% / 6,537 · 100% / 6,027 (narrower a year out)
By horizon — propagated	100% / 8,081 · 100% / 11,974 · 100% / 14,637 (1.8× widening)
Interval width in context	11,667 on a mean of ~19,700 → about ±30%
Seasonal naive	1,660 → the model is about 4% better (4.28%)