State-space models: let the Kalman filter integrate out the latent states

Challenge C07 fits a latent time series the brute-force way: one unknown per day, all of them handed to NUTS. That is the only option there, because stochastic volatility is not linear-Gaussian. When a model is linear and Gaussian - a level that drifts, a seasonal pattern that evolves, regression effects, Normal noise - there is a much better deal on offer. The Kalman filter computes the likelihood of the data with every latent state integrated out exactly. NUTS is left with a handful of variance parameters, and the states are recovered afterwards, also exactly, by the Kalman smoother.

The question is a classic of the genre (Harvey & Durbin, 1986): did the seat-belt law reduce the number of car drivers killed or seriously injured, and by how much?

#1 · Question and data

drivers is the monthly number of car drivers killed or seriously injured (KSI) in Great Britain; rear is the same for rear-seat passengers, who were not covered by the 1983 law - keep them in mind, they will be our negative control. law is 0 until January 1983 and 1 for the 23 months from February 1983.

The file's index is a decimal year. The state-space module wants a pandas DatetimeIndex with a frequency (it uses it to label forecasts), so we build one.

Everything that makes the question hard is visible. The series has a strong seasonal pattern (winter peaks), a level that wanders (down after the 1973-74 oil crisis, up again, down again), and explanatory variables that themselves trend: petrol got dearer and people drove more. After the law the drivers series is lower - but it had been lower in the late 1970s than in the early 1970s too, and single months swing by 30%. A before/after difference in means would confound all of it.

We model log(drivers): effects on road casualties are naturally multiplicative ("10% fewer"), and on the log scale the seasonal swing has a constant size.

#2 · The idea, with nothing hidden

The simplest state-space model is the local level model. There is one latent state, the level $\mu_t$, which does a random walk; we observe it with noise:

$$ \begin{aligned} y_t &= \mu_t + \varepsilon_t, & \varepsilon_t &\sim \text{Normal}(0, \sigma_\varepsilon) && \text{(observation equation)}\\ \mu_{t+1} &= \mu_t + \eta_t, & \eta_t &\sim \text{Normal}(0, \sigma_\eta) && \text{(state equation)} \end{aligned} $$

with $\mu_1 \sim \text{Normal}(a_1, P_1)$. There are 192 unknown levels and two unknown standard deviations. But everything is Gaussian and linear, so the distribution of $\mu_t$ given the data so far is Gaussian too, and it can be carried forward one observation at a time:

1. Predict the next observation. If $\mu_t \mid y_{1:t-1} \sim \text{Normal}(a_t, P_t)$ then $y_t \mid y_{1:t-1} \sim \text{Normal}(a_t,\; F_t = P_t + \sigma_\varepsilon^2)$. The surprise $v_t = y_t - a_t$ is the innovation.
2. Update. Move the estimate towards the observation by the Kalman gain $K_t = P_t / F_t$: $a_{t|t} = a_t + K_t v_t$, $P_{t|t} = P_t (1 - K_t)$.
3. Step forward. $a_{t+1} = a_{t|t}$, $P_{t+1} = P_{t|t} + \sigma_\eta^2$.

Step 1 is the point: multiplying the one-step predictive densities gives $p(y_{1:T} \mid \sigma_\varepsilon, \sigma_\eta) = \prod_t \text{Normal}(y_t \mid a_t, F_t)$ - the likelihood of the two parameters with all 192 levels integrated out. No approximation.

A likelihood of two numbers can be maximised by any optimiser:

The filter works, and the picture shows why a local level on its own is the wrong model. With no seasonal component, the only way to follow the December peaks is to make the level very mobile: the estimated level innovation ($\sigma_\eta \approx 0.11$) is more than twice the observation noise ($\sigma_\varepsilon \approx 0.05$), the "level" chases last month's value, and the one-step predictions are simply the data shifted right by a month. The seasonality has been absorbed into the wrong component. Section 3 fixes that.

First, the promised check that nothing is magic: pymc_extras.statespace should give the same log-likelihood for the same model at the same parameter values.

They agree to three decimals (the package adds a tiny jitter to every covariance for numerical safety, which accounts for the rest). The general filter is the same three steps with matrices: a state vector $x_t$, a transition matrix $T$, a design matrix $Z_t$ that picks out what is observed, and covariances $Q$ (state shocks) and $H$ (observation noise):

$$ y_t = Z_t x_t + \varepsilon_t,\quad \varepsilon_t \sim \text{N}(0, H), \qquad x_{t+1} = T x_t + R\,\eta_t,\quad \eta_t \sim \text{N}(0, Q), \qquad x_1 \sim \text{N}(x_0, P_0). $$

#3 · A structural time-series model with pymc_extras.statespace

A structural model writes the series as a sum of interpretable pieces, each of which is a small state-space model; stacking their state vectors gives one big one:

$$ \log(\text{drivers}_t) = \underbrace{\mu_t}_{\text{level}} + \underbrace{\gamma_t}_{\text{seasonal}} + \underbrace{\beta_p \log(\text{petrol}_t) + \beta_k \log(\text{kms}_t)}_{\text{explanatory}} + \underbrace{\lambda \cdot \text{law}_t}_{\text{intervention}} + \varepsilon_t $$

• Level: the random walk of section 2. A stochastic slope could be added (LevelTrend(order=2)); a series with no persistent direction does not need one, and it is the first exercise at the end.
• Seasonal, "dummy" form: twelve monthly effects constrained to sum to zero over any year, up to a shock: $\gamma_t = -\sum_{j=1}^{11} \gamma_{t-j} + \omega_t$. With $\sigma_\omega = 0$ it is a fixed monthly pattern; with $\sigma_\omega > 0$ the pattern evolves. The state holds the last 11 effects. The alternative trigonometric form (FrequencySeasonality) writes the pattern as a sum of up to six harmonics; its advantage is that you may keep only the first few when the seasonal shape is smooth or the period is long (weekly data with a yearly cycle). For monthly data the full dummy form is cheap.
• Regression and intervention: a coefficient is a state that never moves (zero innovation variance) and a design matrix $Z_t$ that contains the regressor. The law enters as a step dummy, i.e. a permanent level shift $\lambda$ from February 1983. That is the right shape for a law with immediate, near-universal compliance (belt wearing went from about 40% to over 90% within the month). A gradual effect would use a ramp or a decaying pulse as the regressor instead; the shape is an assumption you choose, not something the model discovers.
• Measurement error $\varepsilon_t$: the irregular month-to-month noise.

Components are combined with + and compiled into one model with .build(), which prints what it expects you to supply:

How the package works, in one paragraph. build() creates the symbolic matrices $x_0, P_0, T, Z, R, H, Q$ with placeholders for the parameters in the table. Inside a pm.Model you create a variable with exactly that name (and, ideally, the listed dims; ss_mod.coords supplies the coordinates) for every row, and register the regressors as pm.Data("data_exog", ...) with dims ("time", "state_exog"). Then ss_mod.build_statespace_graph(y) substitutes your variables into the matrices, runs the Kalman filter as a pytensor.scan, and registers the resulting log-likelihood as the observed variable obs. A missing or extra name is an error, which is a good thing.

#Priors, on the log scale

All three standard deviations are in log units, so 0.01 means "about 1%":

• sigma_obs ~ HalfNormal(0.1): irregular noise of a few percent up to 20% a month.
• sigma_level ~ HalfNormal(0.05): the underlying level moves by at most a few percent a month. A level that moved 10% a month would be doing the seasonal component's job, as in section 2.
• sigma_seasonal ~ HalfNormal(0.02): the seasonal pattern changes slowly if at all.

#Where do the coefficients go? Two ways to declare the initial state

The state at the first time point is $x_1 \sim \text{N}(x_0, P_0)$, and $x_0$ is made of the parameters initial_level, params_seasonal and beta_exog. The package documentation's pattern is to give each of those a prior and set P0 to something small: NUTS then samples 1 + 11 + 3 = 15 location parameters plus the three standard deviations.

But a Normal prior on $x_0$ with $P_0 = 0$ is the same model as a fixed $x_0$ at the prior mean with the prior variance in P0 - and in that form the Kalman filter integrates the initial level, the seasonal pattern and the regression coefficients out analytically along with everything else. They are Gaussian and enter linearly, so there is no reason to make NUTS do it. We use that form: Deterministic constants for the means, prior variances on the diagonal of P0, and three sampled parameters. (When this notebook was developed the documentation's form gave the same posterior and took roughly three times as long, because 18 parameters need longer NUTS trajectories and every gradient is a full pass of the filter.) The coefficients come back from the smoother in section 5.

The implied priors: initial level Normal(log 1500, 0.5), monthly effects Normal(0, 0.2), Normal(0, 1) for the two elasticities, and for the law effect a first guess of Normal(0, 1) - "wide and centred on no effect" - which the prior predictive check below will make us reconsider.

(Deterministic rather than pm.Data for the constants: the post-estimation methods look every parameter up in the posterior group, and fail with "Cannot sample from flat variable" if it is not there.)

#Prior predictive check

pm.sample_prior_predictive draws the parameters; sample_unconditional_prior then simulates whole trajectories from the state-space dynamics, using the real regressors.

The left panel is what the check is for. Everything looks like a plausible casualty series until February 1983, where the fan bursts open: with Normal(0, 1) on $\lambda$ about half of the simulated series change by more than a factor of two in the year after the law. No road-safety measure does that; "wide and centred on zero" was wide on a scale nobody had looked at. Normal(0, 0.3) puts two prior standard deviations at -45% / +82%, still generous for a seat belt, and brings the factor-two share down to a few percent (right panel).

What remains wide is the long-run spread: a random walk's variance accumulates over 192 months, and a unit elasticity on a distance index that doubled is allowed. That is acceptable - the month-to-month behaviour, the seasonal amplitudes and the noise level are what these priors need to get roughly right, and they do; about 1% of simulated months or fewer fall outside 100-25,000 drivers.

#4 · Fit, and diagnose

Three parameters adapt quickly, so 500 warm-up steps are plenty (under nutpie the default is only 400 anyway); target_accept=0.9 because sigma_seasonal has posterior mass close to zero, where the log-transformed parameter has a long flat tail and the default step size produced the odd divergence.

Clean: no divergences, r_hat of 1.00, ESS in the thousands from 4 x 1000 draws - what you expect from a three-dimensional posterior. Reading the numbers:

• sigma_obs of about 0.06: the irregular component is about 6% of the monthly total. That is close to what Poisson counting noise alone would give for counts of this size only if they were around 250; at 1,700 a month Poisson noise is 2.4%, so most of the irregular is real month-to-month variation (weather, holidays falling differently, ...).
• sigma_level of about 0.02, a sixth of the value the seasonal-free model of section 2 needed: the level is now a slowly moving thing.
• sigma_seasonal is small, with a posterior piled up against zero: the data cannot distinguish a slowly evolving seasonal pattern from a fixed one. (Harvey and Durbin used a fixed one.) The model is free to choose, and it costs nothing to let it.

#5 · Recover the states: the Kalman smoother

The posterior so far contains three standard deviations and no level, no seasonal pattern and no law effect. They were integrated out, not lost: given the parameters, the distribution of all states given all the data is Gaussian and the Kalman smoother computes it (a backward pass after the forward filter). sample_conditional_posterior does that for every posterior draw and samples one whole state trajectory from each - a joint draw, so cross-time correlations are right (that is the simulation smoother).

Memory matters here: the outputs are (draws x time x states) arrays for the predicted, filtered and smoothed states, and if you ask sample_filter_outputs for smoothed covariances you get (draws x time x states x states) - a third of a gigabyte for 1000 draws of this small model. So: thin first (every 4th draw: 1000 trajectories), keep what you plot, delete the rest.

The coefficients are back, with full posteriors, although NUTS never saw them. The petrol-price elasticity is negative as expected (dearer petrol, less and slower driving, fewer casualties). The distance elasticity is poorly determined - kms is smooth and trending, so it competes with the level for the same slow movements, and the level wins. The law effect is about -0.23 with a standard deviation of 0.05: section 9 turns it into the answer. First, is the sampler comparison as one-sided as the introduction implied?

#6 · The brute-force contrast, measured

The same model - same priors, same data - written the C07 way: every latent shock is a parameter. The level is a cumulative sum of 191 scaled standard-Normal steps (a non-centred random walk). The seasonal recursion $\gamma_t = -\sum_j \gamma_{t-j} + \sigma_\omega\omega_t$ is linear, so it can be unrolled once in NumPy into two fixed matrices, $\gamma = A\,\gamma_{\text{init}} + \sigma_\omega B\,\omega$, which keeps scan out of the brute-force model altogether. NUTS gets 3 standard deviations + 3 coefficients + 1 initial level + 11 initial seasonal effects + 191 + 191 shocks = 400 parameters. The latent shocks are not needed for the comparison, so they are dropped as soon as the summary is made.

Read this honestly, because it is not the rout the introduction may have led you to expect.

• The posteriors agree to within Monte Carlo error, for the variances and for the coefficients - which one method sampled and the other integrated out and recovered with the smoother. Same model, two ways of doing the same integral.
• On raw speed, brute force wins here: several times less wall-clock and three to four times the effective samples per second (the exact ratio moves with the machine's load). At $T = 192$ the brute-force log-density is a cumulative sum and two matrix-vector products, microseconds of vectorised arithmetic. Each Kalman gradient is a 192-step sequential scan through 15 x 15 matrix algebra, differentiated, and the graph takes a while to compile. The Kalman model needs about 6 gradient evaluations per draw against about 120 and gets twice the ESS from the same number of draws - the geometry is far better - but each of its gradients costs on the order of a hundred times more.
• On reliability, the Kalman version wins: zero divergences against dozens at the same target_accept. The brute-force posterior is a 400-dimensional hierarchy in which the shocks and their scales are entangled - the geometry C07 is about. The divergent draws are not confined to the neck of a funnel near sigma_seasonal = 0 (their median is about the same as the overall one), so there is no single cheap fix; the answers happen to agree here, but with brute force you have to earn that trust by raising target_accept, reparameterising and checking where the divergences sit, and every fix eats into the speed advantage.

So for a series this short the honest summary is: the Kalman filter buys clean geometry, exact state inference and the machinery of sections 8-10 - not speed. The balance shifts as series get longer (the filter's cost is linear in $T$ and it still has three parameters; brute force adds two parameters per time step to an already awkward posterior), but this notebook has not measured where the lines cross, so it will not tell you.

#7 · The decomposition

The reason to fit a structural model: every piece has a meaning and a credible band. The bands below are 94% intervals over the 1000 joint smoother draws, so they include parameter uncertainty as well as state uncertainty.

• Level. Rises to about 7.55 in 1973, falls by some 15% through 1974-75 (the oil crisis and its aftermath - over and above what the petrol-price regressor explains), then sits near 7.38 for eight years. Its 94% band is about 0.13 wide: the level is known to roughly ±6%.
• Seasonal. November is about 21% and December 28% above the level, April 13% below; the pattern is essentially the same in every year, which is what sigma_seasonal near zero means.
• Regression effects. Petrol price and distance together move the series by less than ±0.1; the visible step down at the start of 1974 is the petrol-price shock. The law is a step of about -0.23 with a band from -0.32 to -0.15.
• Irregular. Mostly within ±0.1, no visible pattern - but see the next section for why this panel is not a model check.

One feature deserves a second look: after the law the level band fans out and its median drifts up by about 0.09. Level and law effect are competing explanations for the last 23 months - "a big drop from the law, partly offset by a rising level" fits about as well as "a smaller drop and a flat level" - and the posterior correlation between $\lambda$ and the final level is about -0.7. The uncertainty in $\lambda$ already includes this. What separates them at all is the smoothness of the level: with monthly steps of about 0.018, a fall of 0.23 in a single month is a thirteen-sigma event for the random walk, so the model assigns the February 1983 drop to the law. That smoothness is the identifying assumption of the whole analysis, exactly as "parallel trends" is for a difference-in-differences.

#8 · The check that matters: one-step-ahead innovations

The smoothed signal above follows the data closely - of course it does, the smoother has seen the future, and a flexible enough level will "fit" anything. Residuals from a smoothed fit are the time-series version of checking a model on its training data.

The honest check uses the quantity from step 1 of the filter: the innovation $v_t = y_t - \text{E}[y_t \mid y_{1:t-1}]$, the error of a genuine one-month-ahead forecast made without seeing $y_t$. If the model is right, the standardised innovations $v_t / \sqrt{F_t}$ are independent standard Normals - no autocorrelation (especially not at lags 1 and 12), no fat tails. sample_filter_outputs returns the filter's predictive means and variances; we ask for exactly those two outputs and nothing with a states-squared dimension.

The first 15 months are dropped: the filter starts from the vague prior in P0 and needs as many observations as there are states before its predictions mean anything.

Nothing to fix. The standardised innovations have mean 0.03 and standard deviation 0.98; the autocorrelations at lags 1 and 12 - where a missing dynamic or a wrong seasonal would show up - are 0.06 and 0.04; three of 24 lags (8, 18 and 22) touch the ±2/√n band with no pattern, and the Ljung-Box test over 24 lags gives p = 0.10. The Q-Q plot is straight, Shapiro-Wilk does not reject, and fewer months than expected lie beyond ±2. The largest surprise is December 1981 at -2.8.

The first month under the law has an innovation of only about -1, although casualties fell by a fifth: the filter knew a law was arriving (the dummy switches on) and its prior on the effect was wide, so a large move was expected in that month and $F_t$ was large. From March 1983 onwards the effect has been learned and the innovations are as unremarkable as before - which is itself evidence for the step shape: a gradual effect modelled as a step would leave a run of same-signed innovations after the law, and there is none.

#9 · The answer: what did the law do?

$\lambda$ is a shift in log-casualties, so $100\,(e^{\lambda} - 1)$ is the percentage change in drivers killed or seriously injured, holding the level, the season, petrol price and distance driven fixed.

The counterfactual follows from the same additivity. In each post-law month the model says $\log y_t = (\text{everything else})_t + \lambda$; without the law, the same month - same weather, same traffic, same noise - would have produced $y_t\,e^{-\lambda}$. The casualties avoided are the difference, summed over the 23 post-law months. All of the uncertainty comes from $\lambda$, which is as it should be: we are asking about months that really happened.

The answer. The law cut the number of car drivers killed or seriously injured by about 21% (94% interval roughly 14% to 28%). A reduction is certain as far as this model is concerned, a reduction of more than 10% has probability 0.99, and "more than 20%" is close to a coin flip. Over the 23 months to December 1984 that is about 8,000 fewer drivers killed or seriously injured (94% interval about 5,000 to 11,600), some 350 a month, against the 30,399 who were.

This agrees with the classical maximum-likelihood analyses of the same series: Commandeur & Koopman's textbook treatment (2007, a level + seasonal + petrol price + intervention model) reports a coefficient of about -0.24, a fall of about 21%.

What the number does and does not mean:

• It is the effect at fixed distance driven and petrol price. If belts made people drive more, that pathway is excluded.
• It assumes a step. The innovations after the law gave no reason to doubt that.
• It counts drivers only. Front-seat passengers were covered by the law too (front).
• It inherits the identifying assumption of section 7: nothing else shifted the level abruptly in February 1983. That one can be probed.

#A negative-control outcome: rear-seat passengers

An intervention analysis on one series cannot distinguish "the law worked" from "something else changed in February 1983" - a mild winter, a recession, a change in how injuries were recorded. A negative control is an outcome exposed to all of those but not to the intervention. Rear-seat passengers ride in the same cars on the same roads in the same weather, and the 1983 law did not require them to wear a belt. Same model, same priors (the initial level re-centred on a series four to five times smaller).

The same model that finds a 21% fall for drivers finds nothing for rear-seat passengers: a change centred on zero (+0.4%), with a 94% interval of roughly -9% to +10% and a probability of about 0.02 that they fell by more than 10%. Whatever happened in February 1983 happened to the people the law applied to and not to the people sitting behind them. A mild winter, a recession or a change in injury recording would have hit both.

Two honest limits. The rear series is four to five times smaller and noisier (sigma_obs of 0.09 against 0.06), so the interval is wide: this rules out a common shock of -21%, not a modest effect in either direction - including the small increase that the risk-compensation argument predicts (belted drivers taking more risks, unbelted passengers paying for it). And a negative control can only fail to falsify; it does not prove the mechanism.

#10 · Two things the Kalman filter gives you for free

#Missing observations

If $y_t$ is missing, the filter simply skips the update step for that month: the prediction becomes the filtered estimate, and the uncertainty $P_t$ keeps growing until data return. No imputation model, no extra parameters. In pymc_extras you mark missing values as NaN in the data frame; the package warns that it will treat them as hidden states, replaces them by a sentinel internally and masks them in the filter.

To see it work, blank out 18 months (January 1977 - June 1978), refit, and compare the smoother's reconstruction with what really happened.

The reconstruction is not a straight line across the gap. The seasonal states keep rotating while no data arrive, so the December 1977 peak is where it should be, and because this is the smoother the level is pinned from both sides - by 1976 and by late 1978. The price of 18 months without data is visible in the dark band: its width inside the gap is about 0.20 against 0.11 for the same months when they are observed (and 0.12 for the observed months of this fit). The outer band adds the irregular noise and is the one to compare with data: all 18 held-out months fall inside it. The refit itself was unremarkable (no divergences, r_hat 1.00), and needed no change to the model - only NaNs in the data.

#Forecasting

A forecast is the filter's predict step applied repeatedly with no updates: the last state estimate is pushed through $T$, and every step adds another $Q$ to the state covariance, so the uncertainty grows with the horizon in exactly the way the model implies. ss_mod.forecast does this for each posterior draw. Because the model has regressors, it insists on a scenario for them: here petrol stays at its December 1984 price, distance driven repeats its last twelve months, and the law stays in force.

The forecast repeats the seasonal pattern around the last estimated level, with the law still switched on. The predictive standard deviation on the log scale grows from about 0.08 one month ahead to 0.10 at twelve months and 0.12 at twenty-four: one month ahead it is mostly irregular noise plus uncertainty about the current level, and each further month adds another sigma_level-sized step of the random walk. In drivers, the 94% interval for December 1986 runs from about 1,400 to 2,200 around a median of 1,750.

Two caveats. The interval is conditional on the scenario: uncertainty about future petrol prices and traffic is not in it. And a local level model forecasts a flat level by construction - if you believe in a drift, model a slope (exercise 1) and pay for it in wider intervals.

#11 · Limits, stated plainly

• Linear and Gaussian, or nothing. The exact marginal likelihood exists because every conditional distribution stays Gaussian. Monthly counts in the thousands are fine on the log scale; VanKilled (around ten a month) is not, and a Poisson observation equation needs an approximate filter, a particle filter, or brute force. C07's stochastic volatility has the latent state inside the variance of the observation - not linear-Gaussian - which is why it was sampled state by state. Heavy-tailed noise, regime switches and saturating effects break the assumption in the same way.
• It is not fast. Section 6: 40-50 s per fit for a three-parameter posterior, because the filter is a sequential scan and its gradient is expensive, and every new pm.Model (the rear-seat series, the series with a gap) is compiled again. The brute-force model was several times quicker on 192 observations. The Kalman filter's selling points here are geometry and exactness.
• The API, as experienced with pymc-extras 0.15.1 on PyMC 6.3.2. Everything used in this notebook worked. The friction:
  • Parameters are matched by name, and you must supply P0 - there is no default and no diffuse initialisation, so the initial-state prior is always your decision.
  • Every parameter must end up in the posterior group: constants must be pm.Deterministic, not pm.Data ("Cannot sample from flat variable").
  • Regression coefficients are states. Integrating them out (our P0 trick) means their posterior comes from the smoother; extract_components_from_idata returns the coefficient itself, not coefficient x regressor, so contributions are computed by hand.
  • Post-estimation is where memory goes: outputs are (draws x time x states) and covariances (draws x time x states x states). Thin first and name the outputs you want. The whole notebook peaks a little above 1 GB with 1000 thinned draws.
  • mvn_method="cholesky" fails for this model ("Matrix is not positive definite": constant coefficient states have singular smoothed covariances); the slow but safe default "svd" works.
  • Forecasting needs a DatetimeIndex with a frequency and an explicit scenario for the regressors.
  • pm.compute_log_likelihood works and returns one term per month, but each term is the one-step-ahead density $\log p(y_t \mid y_{1:t-1}, \theta)$. az.loo will run on it without complaint, yet it is not leave-one-out: $y_t$ also sits in the conditioning set of every later term. Treat the sum as a prequential score, or compare forecasts on held-out months.

#Try it yourself

1. A stochastic slope. Replace the level with st.LevelTrend(order=2, innovations_order=[1, 1]) (a local linear trend). The model then asks for a two-element initial_level and sigma_level, and P0 grows by one state. Does the law effect move? What happens to the 24-month forecast interval, and why is that the honest price of admitting the series might have a drift?
2. A different model class, compared by forecast score. Fit pymc_extras.statespace.BayesianSARIMAX with order=(0, 1, 1), seasonal_order=(0, 1, 1, 12) (the "airline model") and the same regressors to the data up to December 1982, forecast 1983-84 from both models with the law dummy in the scenario, and compare the log predictive density of the 24 held-out months. (Not az.loo: see the last bullet of the limits above. And compare only months after both filters have burnt in - a differenced model treats its first 13 observations differently.)
3. Front and rear jointly. Give every component observed_state_names=["front", "rear"] and model the two series together, with a full covariance for the level shocks. The law coefficient becomes a 2-vector. Is the front-seat effect sharper once the shared month-to-month shocks are modelled? And is the rear-seat effect still indistinguishable from zero - or slightly positive, as the "risk compensation" argument would predict?

Next: C07 for a latent time series that is not linear-Gaussian, or E12 for what to try when NUTS itself is the bottleneck.