Bayesian structural time-series intervention analysis

Native Bayesian structural time-series counterfactual in the manner of CausalImpact (Brodersen et al. 2015). The response is modelled over the pre-period as a local level plus a static regression on the concurrent covariates, $$y_t = \mu_t + x_t'\beta + \varepsilon_t,\quad \mu_{t+1} = \mu_t + \eta_t,$$ with $\varepsilon_t \sim N(0, \sigma^2)$ and $\eta_t \sim N(0, \sigma_\mu^2)$, after standardising every column on the pre-period (CausalImpact's standardize.data = TRUE). A Gibbs sampler alternates (i) forward-filter backward-sampling of the level path, with the post-period response treated as missing, (ii) the conjugate normal draw of $\beta$, and (iii) inverse-gamma draws of $\sigma^2$ and $\sigma_\mu^2$. Priors follow the bsts defaults CausalImpact uses: level sd prior with guess prior.level.sd (0.01) times the sd of y, 32 prior observations and upper limit sd(y); observation sd prior with guess $\sqrt{1 - 0.8}$ sd(y) and 50 prior observations (expected $R^2 = 0.8$) when there are covariates, upper limit 1.2 sd(y); initial level $N(y_1, sd(y)^2)$. Simplification: bsts' spike-and-slab coefficient prior is replaced by its slab alone, $\beta \mid \sigma^2 \sim N(0, \sigma^2 \Omega^{-1})$ with $\Omega = 0.01\,(0.5 X'X/n + 0.5\,\mathrm{diag}(X'X/n))$, so every covariate is kept (no variable selection). The first 10 percent of the draws are discarded as burn-in. The posterior predictive of the counterfactual (level path + regression + observation noise, back on the original scale) gives the summary table: Actual, Pred, Pred.lower/upper/sd, AbsEffect and RelEffect with interval and sd, alpha, and p, the posterior tail-area probability of the observed post-period sum (smaller tail, CausalImpact's (+1)-corrected rule). Results are reproducible under set.seed() or model_args$seed. CausalImpact itself (with its Boom dependency) is not used and was not available to cross-validate against; the tests check coverage of a known effect, the no-effect case and agreement with a Kalman filter.

Usage

morie_causal_impact(
  data,
  pre_period,
  post_period,
  model_args = NULL,
  alpha = 0.05
)

Arguments

data: A data frame or matrix (or a zoo series) whose first column is the outcome and remaining columns are concurrent covariate predictors.
pre_period: Length-2 vector giving the first and last row (or, for a zoo series, time index) of the pre-intervention window.
post_period: Length-2 vector giving the first and last row (or time index) of the post-intervention window.
model_args: Optional named list of sampler settings, as in CausalImpact's model.args: niter (default 1000), prior.level.sd (0.01), standardize.data (TRUE), seed (NULL: use the session RNG). nseasons > 1 and dynamic.regression = TRUE are not supported and error.
alpha: Posterior credible-interval coverage (default 0.05, meaning 95 percent intervals).

Value

Named list with elements average_effect, cumulative_effect, ci_lower, ci_upper (interval of the average effect), posterior_prob_causal (the tail-area probability p, as before), summary (the Average / Cumulative summary data frame), and impact: a native "morie_causal_impact" list with summary, series (pointwise predictions and effects with intervals over the analysed window), y_samples (posterior predictive draws for the post-period), model (posterior draws of the coefficients and standard deviations, burn-in, settings), pre_period, post_period and alpha. It replaces the former CausalImpact object.

References

Brodersen KH, Gallusser F, Koehler J, Remy N, Scott SL (2015). Inferring causal impact using Bayesian structural time-series models. Annals of Applied Statistics, 9(1):247-274.

Scott SL, Varian HR (2014). Predicting the present with Bayesian structural time series. International Journal of Mathematical Modelling and Numerical Optimisation, 5(1-2):4-23.

Examples

set.seed(1)
x <- cumsum(rnorm(60))
y <- 2 + 0.8 * x + rnorm(60, sd = 0.3)
y[41:60] <- y[41:60] + 2
res <- morie_causal_impact(data.frame(y = y, x = x),
  pre_period = c(1, 40), post_period = c(41, 60),
  model_args = list(niter = 300)
)
res$summary
#>              Actual       Pred Pred.lower Pred.upper   Pred.sd AbsEffect
#> Average      7.9354   5.921068   5.622831   6.289495 0.1881097  2.014332
#> Cumulative 158.7080 118.421356 112.456624 125.789897 3.7621948 40.286635
#>            AbsEffect.lower AbsEffect.upper AbsEffect.sd RelEffect
#> Average           1.645905        2.312568    0.1881097 0.3401974
#> Cumulative       32.918094       46.251368    3.7621948 0.3401974
#>            RelEffect.lower RelEffect.upper RelEffect.sd alpha           p
#> Average          0.2616911       0.4112821   0.04245159  0.05 0.003690037
#> Cumulative       0.2616911       0.4112821   0.04245159  0.05 0.003690037