12 Bayesian Statistics II — Time Series

Lecture 11, Part 1 of ECON 8310.

12.1 Lecture Outline

Why Bayesian Time Series
The Structural Decomposition
Fitting It in PyMC
Reading the Forecast
Key Takeaways

Five parts. Why you’d do this at all, given that four percent. Then the structural decomposition — writing a series as parts you can name. Then fitting it in PyMC, where our first attempt fails. Then reading the forecast, and the two very different questions you can ask about an interval. And takeaways.

12.2 Why Bayesian Time Series

Not a more accurate forecast — a different deliverable.

12.2.1 What Every Model So Far Gave You

Ten weeks of forecasting, and almost every model returned the same object: one number per period. ETS, ARIMA, random forests, XGBoost, LSTMs — a point forecast, and at best an interval derived from asymptotic theory that nobody checked.

For a planner that is often not enough. How much stock do I need so the chance of a stockout is under 5%? is a question about a distribution, and a point forecast cannot answer it however accurate it is.

A Bayesian time series model returns a full predictive distribution for every future period, plus a distribution over each component — the trend, the seasonality, and the noise, each with its own uncertainty.

Set expectations now. This will not beat XGBoost on RMSE, and on our data it barely beats the seasonal-naive benchmark. The argument for it is what it hands you, not where it ranks — and we will hold it to that standard at the end of this lecture.

Take stock for a second. Ten weeks of forecasting, and almost every model has returned the same object: one number per period.

ETS, ARIMA, random forests, XGBoost, the LSTM — a point forecast. And at best, an interval derived from asymptotic theory that, let’s be honest, nobody checked.

For a planner, that is often not enough. Consider the question: how much stock do I need so the chance of a stockout is under five percent?

That is a question about a distribution. A point forecast cannot answer it — not because the point forecast is inaccurate, but because it isn’t the right kind of object. You can’t extract a fifth percentile from a single number.

A Bayesian time series model returns a full predictive distribution for every future period, plus a distribution over each component — the trend, the seasonality, the noise, each with its own uncertainty. And let me set expectations properly, right now.

This will not beat XGBoost on RMSE. On our data it barely beats seasonal naive. The argument for it is what it hands you, not where it ranks.

I’m saying that at the start because at the end of this lecture I’m going to hold it to exactly that standard — and it’s going to partially fail. That’s a more useful hour than one where everything works.

12.2.2 Three Things That Become Possible

Priors carry knowledge	“Weekly growth is very unlikely to exceed 20%” is a sentence you can put in the model, rather than hoping the data implies it.
Components get error bars	Not just a forecast, but “how confident am I that there is a trend at all?” — separately from the seasonality.
Uncertainty compounds honestly	Forecast 52 weeks out and the interval widens as it should, because the model propagates uncertainty forward instead of assuming it away.

That third row is the one people underestimate. A prediction interval from ARIMA widens too, but according to a formula derived under assumptions you did not check. Here the widening is a consequence of the model you actually wrote down.

And the price: it is slower, it needs priors you have to defend, and — as we will see — it can fail silently in ways a least-squares fit cannot.

Three capabilities you didn’t have before.

Priors carry knowledge. “Weekly growth is very unlikely to exceed twenty percent” is a sentence you can put into the model, rather than hoping the data implies it. When you have short history, that’s the difference between a usable forecast and a wild one.

Components get error bars. Not just a forecast, but “how confident am I that there’s a trend at all?” — separately from the seasonality. You can decompose your uncertainty, not just your signal.

Uncertainty compounds honestly. Forecast fifty-two weeks out and the interval widens, because the model propagates uncertainty forward rather than assuming it away. That third row is the one people underestimate, so let me sharpen it.

An ARIMA prediction interval widens too. But it widens according to a formula, derived under assumptions about the error process that you did not check and probably could not check.

Here, the widening is a consequence of the model you wrote down. If you said the level follows a random walk, then fifty-two steps of random walk uncertainty is what you get. You can trace the interval back to a modeling choice and argue about that choice.

That traceability is worth more than it sounds. When a planner tells you the interval is too wide to use — which is exactly what will happen at the end of this lecture — you can point at the specific assumption producing it and propose changing that, rather than shrugging at a formula. An ARIMA interval you can only accept or reject. This one you can negotiate with.

And the price is real: it’s slower, it needs priors you must defend, and — as we’ll see in twenty minutes — it can fail silently in ways a least-squares fit simply cannot.

12.3 The Structural Decomposition

Write the series as a sum of parts you can name, and put a prior on each.

12.3.1 Structural Time Series: One Equation

The whole idea is additive, and you have seen the shape before — it is the GAM decomposition from Lecture 3, with priors attached.

\[y_t \;=\; \underbrace{\mu_t}_{\text{level}} \;+\; \underbrace{s_t}_{\text{seasonality}} \;+\; \underbrace{\varepsilon_t}_{\text{noise}}, \qquad \varepsilon_t \sim \mathcal{N}(0, \sigma^2)\]

Each piece is a thing a business person can name, and each gets its own posterior. “Demand is growing” is a statement about $\mu_t$; “December is our peak” is a statement about $s_t$; “how noisy is this store?” is $\sigma$.

This is why the approach is called structural. You are not fitting a black box to a curve — you are stating what you believe the series is made of, and asking the data how much of each part there is.

Compare ARIMA, where trend and seasonality are handled by differencing and the parameters are correlations rather than quantities anyone can interpret.

The whole idea is additive, and you have seen this shape before. y-t equals mu-t, the level, plus s-t, the seasonality, plus epsilon-t, the noise.

That is the GAM decomposition from Lecture 3, with priors attached. Same additive structure, same named components — now each one is a random variable with a distribution.

And here’s why that matters practically. Each piece is a thing a business person can name, and each gets its own posterior.

“Demand is growing” is a statement about mu-t. “December is our peak” is a statement about s-t. “How noisy is this store?” is sigma. Those are three separate questions, and the model answers them separately, with uncertainty on each.

That’s why the approach is called structural. You are not fitting a flexible curve to data and hoping. You are stating what you believe the series is made of, and asking the data how much of each part there is.

Contrast ARIMA, which we spent Lecture 2 on. There, trend and seasonality are handled by differencing — you remove them rather than model them — and the parameters that come out are autocorrelation coefficients. Nobody has ever walked into a planning meeting and explained the business implications of an MA(2) coefficient. Here, every parameter is a quantity somebody outside the room can interpret. There’s a real trade being made, though, and I don’t want to sell only one side.

ARIMA’s differencing is automatic. You don’t have to know whether the trend is linear or wandering or damped — you difference, and the trend is gone, whatever shape it had.

Structural modeling makes you commit. You have to say what kind of level you believe in before you fit. Commit correctly and you get interpretable components and honest uncertainty. Commit wrongly and you have baked a false assumption into the model, and the posterior will look perfectly confident about it.

That’s the same bargain as the DAG in Lecture 12: stating your assumptions gets you more, and makes you responsible for them.

12.3.2 The Level: A Random Walk With a Prior On It

The simplest useful trend is a local level — the series has a level, and the level itself drifts:

\[\mu_t = \mu_{t-1} + \eta_t, \qquad \eta_t \sim \mathcal{N}(0, \sigma_{\text{level}}^2)\]

Read it as a sentence: this week’s level is last week’s level, plus a shock. And $\sigma_{\text{level}}$ — which you put a prior on — controls how fast the level is allowed to move.

Small $\sigma_{\text{level}}$	A nearly constant level. Rigid, but stable far into the future.
Large $\sigma_{\text{level}}$	A level that chases recent data. Fits well, and forecasts wander.

That is a familiar trade in new clothing — it is the smoothing parameter from Lecture 1’s exponential smoothing, and the penalty strength from Lecture 6, now expressed as a prior you can argue about.

The simplest useful trend is a local level: the series has a level, and the level itself drifts over time.

mu-t equals mu-t-minus-one plus a shock. Read as a sentence: this week’s level is last week’s level, plus a random disturbance.

And that disturbance has a standard deviation, sigma-level, which you put a prior on. That single parameter controls how fast the level is allowed to move. Small sigma-level: a nearly constant level. Rigid, but stable far into the future.

Large sigma-level: a level that chases recent data. Fits the training period beautifully, and the forecasts wander off. And that is a familiar trade in new clothing.

It is the smoothing parameter from Lecture 1’s exponential smoothing — alpha, deciding how fast to forget. It is the penalty strength from Lecture 6. Now it’s expressed as a prior you can argue about, and — importantly — the data gets a say in it too, because sigma-level is estimated, not fixed. That’s the sixth appearance of this idea coming into view. We’ll name it properly next week.

12.3.3 Seasonality as Fourier Terms

Fifty-two weekly seasonal parameters would be a lot to estimate. Instead, build seasonality out of a few smooth waves:

\[s_t = \sum_{k=1}^{K}\left[a_k \sin\!\left(\tfrac{2\pi k t}{P}\right) + b_k \cos\!\left(\tfrac{2\pi k t}{P}\right)\right]\]

$P$ is the period — for our weekly data, $P = 52$. $K$ is how many waves you allow: $K=1$ is a single smooth annual cycle, larger $K$ lets the shape get more detailed. Each $a_k, b_k$ is just a coefficient with a Normal prior.

You have met this exact device twice. It is how Prophet encodes seasonality in Lecture 3, and it is the Fourier feature block we engineered by hand for Homework 4. The difference now is that $a_k$ and $b_k$ come back as distributions.

$K$ trades flexibility against overfitting, exactly like spline knots. We use $K = 3$ — six coefficients standing in for fifty-two.

Now seasonality. Our data is weekly with an annual cycle, so a naive approach would estimate fifty-two separate seasonal parameters. That’s a lot to estimate, and most of it would be noise. Instead, build the seasonal shape out of a few smooth waves.

s-t is a sum over k of sine and cosine terms at increasing frequencies. P is the period — fifty-two for us. K is how many waves you allow. K equals one is a single smooth annual cycle; larger K lets the shape get more detailed, with sharper peaks. And each a-k and b-k is just a coefficient with a Normal prior. Nothing exotic. You have met this exact construction twice already.

It is how Prophet encodes seasonality, in Lecture 3. And it is the Fourier feature block you engineered by hand for Homework 4 — those sin_1, cos_1, sin_2 columns.

Same device, third appearance. The difference now is that a-k and b-k come back as distributions, so you can ask how confident you are that the seasonal peak is where you think it is.

K trades flexibility against overfitting, exactly like spline knots in Lecture 3. We use K equals three — six coefficients standing in for fifty-two parameters.

Notice what that construction assumes, because it’s an assumption and not a free lunch: it assumes the seasonal shape is smooth.

Three sine-cosine pairs can describe a broad annual hump with a couple of gentle bumps. What they cannot describe is a single enormous one-week spike — Black Friday, say — sitting on an otherwise flat year. Fourier terms will smear that spike across neighboring weeks, because smooth curves are all they can make.

If your business has sharp calendar events, the honest answer is to model them as separate indicator variables rather than raising K until the waves can approximate a spike. That’s what Prophet does with its holiday regressors, and it’s what you did by hand with the SNAP and event flags in Homework 2.

12.4 Fitting It in PyMC

The model in code, and the check you run before believing any of it.

12.4.1 The Model in PyMC

Local level $+$ Fourier seasonality

with pm.Model() as model:
    sigma_level = pm.Exponential("sigma_level", 10.0)
    z           = pm.Normal("z", 0, 1, shape=n)     # non-centered
    level = pm.Deterministic("level", pm.math.cumsum(z * sigma_level))
    beta  = pm.Normal("beta", 0, 0.5, shape=2*K)    # Fourier
    sigma = pm.Exponential("sigma", 1.0)

    pm.Normal("obs", mu=level + pm.math.dot(X, beta),
              sigma=sigma, observed=y_scaled)
    idata = pm.sample(1500, tune=2000, chains=4, target_accept=0.95)

Line for line this is the equations above — cumsum of scaled Normal draws is a random walk. Standardize $y$ first: these priors assume unit scale, and on raw units of 20,000 they would be absurdly tight.

Here is everything from the last three slides, in about eight lines.

sigma_level with an Exponential prior — it must be positive, so Exponential, not Normal. Then z, a vector of standard Normals. Then the level is the cumulative sum of z times sigma_level.

Stop on that line for a second, because it’s the clever bit: cumsum of scaled Normal draws is a random walk. That’s the definition. Each step adds a scaled Normal shock to the running total.

Then beta, the Fourier coefficients. Then sigma, the observation noise. Then the likelihood: the observed data is Normal around level plus the Fourier terms. One thing you must do, and it’s bolded for a reason. Standardize y first.

These priors assume the data is on roughly unit scale. Our weekly units are in the tens of thousands. An Exponential prior with rate ten on sigma_level is a sensible prior for a standardized series and an absurdly tight one for raw units of twenty thousand — you’d be telling the model the level essentially cannot move.

Standardize, fit, then transform the forecasts back. Getting this wrong produces a model that fits terribly for reasons that look mysterious.

This is the same lesson as Lecture 6’s warning about standardizing before Ridge, and it will come back one more time in Lecture 12 with the opposite sign — where standardizing a variable you then interpret silently destroys the interpretation.

So the rule isn’t “always standardize” or “never.” It’s: know what scale each variable is on, and know what your priors assume about that scale. A prior is a statement about magnitudes. If you don’t know the magnitudes, you don’t know what you’ve assumed.

12.4.2 The First Fit Failed — Here Is What That Looks Like

Written the obvious way, with pm.GaussianRandomWalk, this model samples badly — even at target_accept=0.95. Lecture 10’s checks caught it:

	First attempt	After the fix	Want
$\hat{R}$	1.022	1.002	$< 1.01$
ESS (bulk)	252	2,376	$> 400$
Divergences	0	0	$0$

Nothing crashed. The model returned a forecast, and the forecast looked reasonable. Only the diagnostics said the chains had not agreed with each other.

The fix was reparameterization, not more samples. Writing the walk as cumsum(z * sigma) with $z \sim \mathcal{N}(0,1)$ — the non-centered form — gives the sampler a geometry it can navigate. BMCP §4.6.1 covers exactly this, and you will meet it again next week.

Now something worth your attention.

The obvious way to write this model uses pm.GaussianRandomWalk, which is exactly what it sounds like and exactly what the equation says. We wrote it that way first. And it sampled badly — even with target_accept raised to point-nine-five.

R-hat: one-point-oh-two-two, against a threshold of one-point-oh-one. Failed. ESS: two hundred fifty-two, against four hundred. Failed. Divergences: zero — that one passed. After the fix: R-hat one-point-oh-oh-two, ESS two thousand six hundred seventeen. Here is what I want you to take from this.

Nothing crashed. The model returned a forecast. The forecast, plotted, looked entirely reasonable — a sensible level, a sensible seasonal shape, plausible intervals. If we had not run the diagnostics, we would have shipped it and never known.

The only thing that told us was R-hat saying the four chains had not agreed with each other.

That is precisely why Lecture 10 made those three numbers non-negotiable. This is the failure mode they exist for. And the fix is worth understanding, because you’ll meet it again next week.

The fix was reparameterization, not more samples. Not a longer run, not a smaller step size — a different way of writing the same model.

Writing the walk as cumsum(z * sigma), with z a standard Normal, is called the non-centered form. It gives the sampler a geometry it can actually navigate. Identical model, identical posterior in principle, completely different sampling behavior.

BMCP section four-point-six-point-one covers exactly this, and next week we’ll see the same fix applied to a hierarchy — where it will matter even more.

12.5 Reading the Forecast

The point estimate, the interval, and whether either is any use.

12.5.1 Does It Work? The Point Forecast

CA_1 FOODS, 52 weeks held out, forecast from the fitted model.

Model	Test RMSE
Bayesian structural TS	1,589
Seasonal naive (same week last year)	1,660

It beats the benchmark by about 4%. That is a real improvement and a small one, and after nine weeks of this course you should be suspicious of anyone who reports it as a triumph.

It is also worth remembering what this model was given: the series and a calendar. No lag features, no price, no SNAP flags, no store identity. The XGBoost model in Homework 4 had 46 engineered predictors and thirty series to learn from.

So the point forecast is not the reason to do this. The next slide is.

CA_1 FOODS, fifty-two weeks held out.

Bayesian structural model: RMSE one thousand five hundred eighty-seven. Seasonal naive — same week last year, one line of code: one thousand six hundred sixty. So we beat the benchmark by about four percent.

That is a real improvement and a small one. And I’ll repeat what I said at the start: after ten weeks of this course, you should be suspicious of anyone who reports four percent as a triumph. But let me also be fair to the model, because the comparison isn’t like-for-like.

Look at what this model was given: the series, and a calendar. That’s it. No lag features. No price. No SNAP flags. No store identity. No cross-series information.

The XGBoost model in Homework 4 had forty-six engineered predictors and thirty series to learn from, and beat naive by a much wider margin. It also had a great deal more information.

And there’s nothing stopping you from combining them, which is what a serious production system usually does. Add price and promotion regressors to this structural model and the level has less work to do, so it drifts less, so the interval narrows. You keep the predictive distribution and you buy back some accuracy.

We’re not doing that today because I want the mechanism visible rather than buried under covariates. But it is the obvious next step, and it’s a very reasonable direction for a final project.

So: the point forecast is not the reason to do this. Which brings us to the slide that is.

12.5.2 The Interval Is the Product

Every forecast comes with a 94% predictive interval. Two questions about it, and they are different questions.

	Measured	Reading
Calibrated?	100% of weeks fell inside	Not overconfident — but it should be about 94%, so it is conservative
Useful?	Width $\approx$ 11,700 units, on a mean of 19,700	About $$30%. Honest, and too wide to order against

Passing the first does not excuse failing the second. A model saying “somewhere between 14,000 and 26,000” is telling the truth and helping nobody. The cause is structural — a level free to walk for 52 steps drifts, and the model says so. The fix is a model with a stronger opinion: a damped level, or covariates that pin it down.

But notice what just happened. We could ask whether the uncertainty was honest. For every point-forecast model in this course, that question had no answer.

Every forecast comes with a ninety-four percent predictive interval. And there are two questions you can ask about it, which are different questions, and this is the heart of the lecture.

Is it calibrated? Measured: one hundred percent of held-out weeks fell inside the interval. Reading: not overconfident. Though it should be about ninety-four percent, so it’s actually conservative — a bit too wide.

Is it useful? Measured: the width is roughly eleven thousand seven hundred units, on a mean of nineteen thousand seven hundred. That’s about plus or minus thirty percent. Honest, and far too wide to order against. Passing the first does not excuse failing the second.

A model that says “demand next week will be somewhere between fourteen thousand and twenty-six thousand” is telling you the truth and helping nobody. You cannot write a purchase order against that.

And the cause is structural, not a bug. A level that is free to random-walk for fifty-two steps drifts, and the accumulated uncertainty is genuinely that large. The model is correctly reporting the consequence of what we told it to believe.

So the fix is a model with a stronger opinion: a damped level that reverts rather than wandering, or covariates that pin it down — price, promotions, the things we withheld. But notice what just happened, because this is the argument for the entire lecture.

We could ask whether the uncertainty was honest. We computed coverage against held-out data, and got an answer.

For every point-forecast model in this course — every single one, from ETS through the Transformer — that question had no answer at all. There was nothing to check. You got a number, and the number was right or wrong, and that was the end of the conversation.

A calibrated-but-too-wide interval is a worse forecast and a better deliverable, because it tells you something true about how much you actually know. That is worth paying four percent for.

12.6 Key Takeaways

What structural models are for.

12.6.1 Lecture 11 Part 1: Key Takeaways

A structural model writes the series as named parts — level, seasonality, noise — and gives each one a prior and a posterior.
The local level is a random walk whose step size $\sigma_{\text{level}}$ you put a prior on. Small means rigid, large means it chases the data — the same trade as a smoothing parameter or a penalty.
Fourier terms encode seasonality in a few coefficients instead of 52. Same device as Prophet in Lecture 3 and the features you built in Homework 4.
Check the sampler before the forecast. Our first fit had $\hat{R}$ = 1.02 and ESS = 252 while returning a plausible-looking answer. Reparameterizing fixed it.
The point forecast beat seasonal naive by 4%. That is not the argument.
The interval was calibrated but too wide to act on — and being able to establish that is the actual product.

Next: hierarchical models — what to do when you have thirty series and some of them barely have data.

One. A structural model writes the series as named parts — level, seasonality, noise — and gives each one a prior and a posterior.

Two. The local level is a random walk whose step size you put a prior on. Small means rigid, large means it chases the data. Same trade as a smoothing parameter or a penalty.

Three. Fourier terms encode seasonality in a few coefficients instead of fifty-two. Same device as Prophet in Lecture 3 and the features you built in Homework 4 — third appearance.

Four. Check the sampler before the forecast. Our first fit had R-hat one-point-oh-two and ESS two fifty-two, while returning a perfectly plausible-looking answer. Reparameterizing fixed it. Five. The point forecast beat seasonal naive by four percent. That is not the argument.

And six. The interval was calibrated but too wide to act on — and being able to establish that is the actual product.

Next time: hierarchical models. What to do when you have thirty series and some of them barely have any data — and a question you have to answer before you fit anything.