CausalImpact: Counterfactuals from Bayesian Structural Time Series

Matthew Reda Causal inference

A campaign launches, a price changes, a competitor exits — and the measurement question arrives after the fact, with no holdout group in sight. CausalImpact answers it by forecasting the world in which the intervention never happened, then reading the effect as the gap between forecast and reality. This post walks through the machinery: the state-space model behind the forecast, the prior that picks control series automatically, and why the uncertainty bands honestly widen the longer you look.

The method comes from Brodersen, Gallusser, Koehler, Remy, and Scott (2015), and its central move is a reframing: a causal effect on a time series is the difference between the observed trajectory and the counterfactual one where the intervention never occurred. The first is data; the second is a forecast. Fit a flexible time-series model on pre-intervention data, run it forward as if nothing had changed, and difference the prediction against reality. Everything else in the paper is about making that forecast trustworthy.

The Effect Is a Forecast Error

Write \( y_t \) for the observed outcome — daily site visits, weekly sales — and \( \tilde{y}_t \) for the counterfactual outcome after the intervention at time \( n \). The causal impact at each post-intervention time point is simply \( y_t - \tilde{y}_t \). The difficulty is that \( \tilde{y}_t \) is never observed, so it must be inferred from two sources: the treated series' own pre-intervention behavior, and contemporaneous control series unaffected by the intervention.

The classical tool here, difference-in-differences, is too rigid for real outcome series: it compares pre/post means across two groups, assumes independent observations, and has no vocabulary for trends, seasonality, or serial correlation. Brodersen et al. generalize it with a diffusion-regression state-space model combining four things:

The generality is real: difference-in-differences, Abadie-style synthetic control, ARIMA models, and ordinary multiple regression all fall out as special cases under particular parameter restrictions.

A Model Assembled from Blocks

The engine is a linear-Gaussian state-space model. Two equations define it: one linking observations to a latent state, one governing how the state evolves.

Definition: BSTS state-space form

Observation equation:

$$ y_t = Z_t^\top \alpha_t + \varepsilon_t, \qquad \varepsilon_t \sim \mathcal{N}(0, \sigma_\varepsilon^2) $$

State equation:

$$ \alpha_{t+1} = T_t \alpha_t + R_t \eta_t, \qquad \eta_t \sim \mathcal{N}(0, Q_t) $$

Here \( y_t \) is the scalar outcome, \( \alpha_t \) is a \( d \)-dimensional latent state, \( Z_t \) is the output vector, \( T_t \) the transition matrix, and \( R_t \) maps a \( q \)-dimensional system error (\( q \le d \)) into the state.

The design is deliberately modular. The state vector concatenates independent components — trend, seasonality, regression — and the system matrices become block-diagonal, one block per component. Adding annual seasonality to a model that already has a trend is a matter of appending blocks to \( Z_t \), \( T_t \), and \( Q_t \).

The regression component is where the control series enter. In its static form the coefficients are fixed: the contribution \( \beta^\top \mathbf{x}_t \) folds into the observation vector with a degenerate state, so the treated series is predicted, week by week, from a weighted combination of its controls. A dynamic variant lets each coefficient evolve as its own random walk for settings where the treated–control relationship drifts over time; static coefficients are the default when the controls are the same kind of quantity as the outcome and the relationship is stable.

💡 The regression component is a synthetic control

The weighted combination of untreated series that predicts the treated one is exactly the synthetic-control idea. But where the classic construction restricts weights to be non-negative and sum to one, the BSTS regression places no sign or sum constraint on \( \beta \), and the trend and seasonal blocks absorb structure the donors alone cannot explain. Scott and Varian (2014) used the same trend-plus-regression architecture for nowcasting before it was repurposed for counterfactual inference.

Trend and Seasonality

The default trend block is the local linear trend: a level \( \mu_t \) and a slope \( \delta_t \), each evolving as a random walk,

$$ \mu_{t+1} = \mu_t + \delta_t + \eta_{\mu,t}, \qquad \delta_{t+1} = \delta_t + \eta_{\delta,t}, $$

with \( \eta_{\mu,t} \sim \mathcal{N}(0, \sigma_\mu^2) \) and \( \eta_{\delta,t} \sim \mathcal{N}(0, \sigma_\delta^2) \). The two variances have clean interpretations: \( \sigma_\mu^2 \) sets how quickly the level jumps; \( \sigma_\delta^2 \) sets how quickly the slope reorients. This component adapts fast to local variation, but it pays for the flexibility at long horizons: a slope that is itself a random walk drifts without bound, and prediction intervals inflate accordingly. For longer-range work the paper offers a stabler variant in which the slope reverts to a long-term value \( D \) at rate \( \rho \), with \( |\rho| \lt 1 \):

$$ \delta_{t+1} = D + \rho\,(\delta_t - D) + \eta_{\delta,t}. $$

Seasonality is handled by a sum-to-zero component. With \( S \) seasons per cycle, each new seasonal effect is minus the sum of the previous \( S-1 \), plus noise:

$$ \gamma_{t+1} = -\sum_{s=0}^{S-2} \gamma_{t-s} + \eta_{\gamma,t}, $$

so the effects average to zero over a full cycle and the level component keeps custody of the series' overall height. Daily data commonly carries two such blocks: day-of-week (\( S = 7 \)) and annual. Each component's diffusion variance gets an inverse-Gamma prior; the default is weakly informative, centering the prior diffusion variance at roughly 10% of the sample variance with the weight of a single observation.

Spike-and-Slab Control Selection

A practical counterfactual analysis rarely has one obvious control series. It has dozens — every untreated market, every industry index, every related search-volume series. Including all of them invites overfitting; hand-picking a subset injects the researcher degrees of freedom the analysis is supposed to avoid. The BSTS answer makes selection part of the inference, via a spike-and-slab prior on the regression (George & McCulloch, 1993).

Attach to each candidate a binary inclusion indicator \( \varrho_j \), with \( \varrho_j = 1 \) meaning \( \beta_j \neq 0 \). The prior factorizes as

$$ p(\varrho, \beta, 1/\sigma_\varepsilon^2) \;=\; p(\varrho)\; p(\sigma_\varepsilon^2 \mid \varrho)\; p(\beta_\varrho \mid \varrho, \sigma_\varepsilon^2), $$

where the spike is an independent Bernoulli on each indicator,

$$ p(\varrho) \;=\; \prod_{j=1}^{J} \pi_j^{\varrho_j} (1 - \pi_j)^{1 - \varrho_j}, $$

and the slab is a weakly informative Gaussian on the coefficients that survive, \( \beta_\varrho \mid \varrho, \sigma_\varepsilon^2 \sim \mathcal{N}(\mathbf{b}_\varrho,\, \sigma_\varepsilon^2 (\Sigma_\varrho^{-1})^{-1}) \). Rather than eliciting each \( \pi_j \), the recommended shortcut is to name an expected model size \( M \) and set \( \pi_j = M/J \); setting an individual \( \pi_j \) to 1 or 0 forces a control in or out. The slab's precision follows Zellner's g-prior, \( \Sigma^{-1} = \tfrac{g}{n}\{ w\, X^\top X + (1-w)\, \mathrm{diag}(X^\top X) \} \), worth \( g \) observations of prior information (default \( g = 1 \)), with the diagonal average (default \( w = 0.5 \)) keeping the prior proper when candidates are collinear.

Two properties make this machinery earn its keep. Computationally, the Gibbs sampler flips one indicator at a time, and every matrix it touches has dimension equal to the number of included controls, so a sparse model over hundreds of candidates stays fast. Statistically, the counterfactual never commits to any single subset: the posterior predictive averages over control sets weighted by posterior probability — Bayesian model averaging doing the robustness work that manual sensitivity analyses approximate.

From Forecast to Impact

Fitting the model on pre-intervention data yields the object everything else is computed from: the posterior predictive distribution over counterfactual trajectories,

$$ p\big(\tilde{\mathbf{y}}_{n+1:m} \mid \mathbf{y}_{1:n}, \mathbf{x}_{1:m}\big), $$

conditioned on the pre-period outcomes and the full run of control series, with all parameters — coefficients, inclusion indicators, variances — integrated out. It is a joint distribution over the whole post-period, not a set of per-week marginals, so serial correlation in the counterfactual is preserved. Each MCMC iteration draws one complete alternate history \( \tilde{\mathbf{y}}_{n+1:m}^{(\tau)} \).

Deep diveThe Gibbs sampler behind the counterfactual

The spike-and-slab prior rules out closed-form posteriors, so inference alternates two blocks. The state step draws the full latent trajectory \( \alpha_{1:m} \) given parameters, using the simulation smoother of Durbin and Koopman (2002) — a refinement of the forward-filtering, backward-sampling algorithms of Carter and Kohn (1994) and Frühwirth-Schnatter (1994). The trick is a Gaussian identity: the conditional variance of the states does not depend on the observed data, so one can simulate a fake dataset and state sequence from the model, subtract the fake sequence's smoothed mean, and add back the real data's smoothed mean from the Kalman filter — an exact draw from the state posterior.

The parameter step conditions on the states: with \( \alpha \) known the component errors are observable, so each diffusion variance is a conjugate Gamma draw, while the regression block updates each \( \varrho_j \) given the others (a two-point distribution) and draws \( \beta_\varrho \) and \( \sigma_\varepsilon^2 \) from conjugate formulae involving only the included columns. Each iteration then runs the model forward through the post-period to emit one counterfactual trajectory. The algorithm is linear in the number of time points and quadratic in the state dimension — under 30 seconds for 500 time points, 10 covariates, and 10,000 iterations.

With posterior draws in hand, the causal quantities are arithmetic. The pointwise impact at each post-period time \( t \) is, per draw,

$$ \phi_t^{(\tau)} \;=\; y_t - \tilde{y}_t^{(\tau)}, $$

and the collection over draws is a full posterior for the effect at each time point — the time-series generalization of the difference-in-differences contrast. The cumulative impact \( \sum_{t'=n+1}^{t} \phi_{t'}^{(\tau)} \) totals the effect from launch; it is the headline number for flow outcomes like weekly sales or daily clicks, but meaningless for stock outcomes like subscriber counts, where summing levels double-counts. For stocks, the running average \( \tfrac{1}{t-n}\sum_{t'=n+1}^{t} \phi_{t'}^{(\tau)} \) is the interpretable summary. Each quantity is reported as a posterior mean with a central 95% credible interval, and an effect is called significant when that interval excludes zero.

Observed minus counterfactual, pointwise and cumulative

The top panel shows the observed series (green) against the counterfactual forecast (dashed) with its 95% credible band; the band fans out with horizon, because a random-walk trend admits more alternate histories the further from the last pre-period point. The bottom panel accumulates the pointwise gap into the cumulative effect. Slide the true intervention effect toward zero and the cumulative interval swallows the zero line — the honest read of “no detectable effect.”

18
Pointwise effect (final week)—
Cumulative effect—
Interval excludes zero?—

Report cumulative impact for a flow like weekly sales — and treat any tool that shows constant-width uncertainty at week 18 as overconfident, not more precise.

Why the bands widen — and why that is correct

Every CausalImpact plot shares a visual signature: the interval around the counterfactual fans out as the post-period lengthens. This is not a defect. The local linear trend's slope is a random walk, so uncertainty about where the series would have gone compounds with horizon — the further from the last pre-period observation, the more alternate histories are consistent with the data. The paper's simulations confirm the consequence: absolute estimation error grows with campaign length, and structural breaks accelerate the decay. The design implication: measure over the shortest window the business question allows, and treat any tool reporting constant-width uncertainty at month six as overconfident, not more precise.

The Google Campaign

The paper's empirical test wraps an observational method around a randomized experiment and checks whether the two agree. A Google advertiser ran search ads for six consecutive weeks, geo-targeted to 95 of 190 US designated market areas with randomized assignment; the outcome was total search-related visits to the advertiser's site, with roughly a year of pre-campaign history. The model used a local level plus a static regression under a spike-and-slab prior with an expected model size of three, run for 10,000 MCMC iterations. Three analyses:

The simulation study backs up the calibration story: a true effect of 25% was detected in over 90% of runs, effects below 1% were correctly left unflagged about 90% of the time, and the central 95% interval covered the true effect at close to its nominal rate across campaign lengths.

The Assumptions That Carry the Weight

The counterfactual is only as good as the bridge from pre-period to post-period, and two assumptions form that bridge.

⚠️ Controls must be untouched, and relationships must hold

First: the control series must not themselves be affected by the intervention. A national campaign that spills into “control” markets, or steals share from the competitor series used as a predictor, contaminates the counterfactual and biases the estimate. Second: the relationship between controls and the treated series, learned entirely from pre-intervention data, must remain stable through the post-period. Nothing in the model can detect a post-period break in that relationship — it will silently misattribute the break to the intervention.

Neither assumption is checkable from the treated series alone, which is why the workflow leans on design: choose controls that are structurally insulated from the treatment (industry aggregates, geographically separated markets), inspect the pre-period fit — a model that cannot track the treated series before the intervention has no business forecasting it after — and run placebo analyses on units that should show nothing, the falsification check that makes an observational estimate credible.

Versus Synthetic Control

CausalImpact and Abadie-style synthetic control aim at the same estimand — the effect of an intervention on a single treated unit, inferred from untreated peers — and a practitioner should think of them as siblings rather than rivals. The differences are in the machinery:

Synthetic controlCausalImpact (BSTS)
CounterfactualConvex combination of donor units matched on the pre-periodPosterior predictive forecast from a state-space model
Donor weightsNon-negative, sum to oneUnconstrained regression coefficients; may be negative
Trend & seasonalityMust be inherited from the donorsExplicit local-trend and seasonal state components
Control selectionAnalyst curates the donor poolSpike-and-slab prior averages over control subsets
UncertaintyPlacebo permutation inferenceCredible intervals from the joint posterior

Formally, the synthetic-control construction is a special case: fix the trend variance at zero, drop seasonality, constrain the weights, and the BSTS regression collapses to it. The Bayesian version buys three practical upgrades — structure the donors cannot supply, automatic selection when the donor pool is large, and uncertainty that propagates parameter and forecast error jointly — at the price of the modeling assumptions above.

Within a measurement program, this machinery is not only for retrospective rescue jobs. It is the natural readout for experiments you did design: the counterfactual for treated geos in a geo-holdout chosen by expected information gain is exactly a BSTS forecast from the untreated markets. And the lift it produces — a posterior, with honest width — is the form of evidence that calibrating a media mix model against experiments requires: the interval, not just the point, determines how hard the MMM's parameters should update.

Takeaways

  • CausalImpact turns causal measurement into forecasting: fit a structural time-series model on pre-intervention data, forecast the no-intervention world, and read the effect as observed minus predicted.
  • The model is modular — local linear trend, sum-to-zero seasonality, and a regression on control series snap together as independent state-space blocks.
  • The spike-and-slab prior selects control series automatically and averages over subsets, replacing hand-picked donors with Bayesian model averaging.
  • Report cumulative impact for flows, running-average impact for stocks — and expect intervals to widen with horizon, because they should.
  • The estimate rests on two assumptions: controls unaffected by the treatment, and a stable pre-period relationship. Placebo tests and pre-period fit checks are the defense.

References