Posterior predictive confidence intervals differ to calculation by `predict.bsts`
- Dominant language
- R
- Stars
- 1.9k
- Forks
- 264
- PR merge metrics
- No merged PRs in 30d
Description
Hi, from looking at the code I can see that `CausalImpact` generates the posterior predictive samples by getting the state value draws and sampling noise with variance equal to the `sigma.obs` draws from the bsts model object.
When I try to replicate the inference step using `predict.bsts`, I get different results. In particular, the credible intervals produced by generating the posterior predictive with `bsts.predict` are larger than `CausalImpact` suggests. The bounds of the relative difference credible interval can differ by 2-3%. Is there inherently a difference between how `CausalImpact` prediction and `predict.bsts` works? I've attached example code below.
```
library(magrittr)
library(CausalImpact)
# Dummy data
set.seed(1)
x1 <- 8000 + arima.sim(model = list(ar = 0.99), n = 100)
y <- 1.2 * x1 + rnorm(100, 0, 1)
y[71:100] <- y[71:100] + 10
data <- cbind(y, x1)
pre.period <- c(1, 70)
post.period <- c(71, 100)
# Run CausalImpact
impact <- CausalImpact(data, pre.period, post.period, model.args = list(niter = 10000, standardize.data = F))
# Predictions from bsts
bsts_model <- impact$model$bsts.model
bsts_predict <-
predict.bsts(bsts_model,
newdata = data[post.period[1]:post.period[2], "x1"])
# Actuals as a matrix
observed_post <-
data[post.period[1]:post.period[2], "y"] %>%
rep(nrow(bsts_predict$distribution)) %>%
matrix(nrow = nrow(bsts_predict$distribution), byrow = T)
# Calculate differences
ppd_diff <- observed_post - bsts_predict$distribution
# Cumulative to last day of post period
diff_cum <- rowSums(ppd_diff)
reldiff_cum <- rowSums(ppd_diff)/rowSums(observed_post)
# Calculated stats
cat("Absolute difference:\n")
c(quantile(diff_cum, 0.025), mn = mean(diff_cum), quantile(diff_cum, 0.975)) %>%
print()
# Returns 75.6, 204.3, 342.3
cat("Relative difference:\n")
c(quantile(reldiff_cum, 0.025), mn = mean(reldiff_cum), quantile(reldiff_cum, 0.975)) %>%
print()
# Returns 0.026%, 0.071%, 0.12%
# Compare to CausalImpact CIs
summary(impact)
# Returns 123.2, 203.7, 288.0
# 0.043%, 0.071%, 0.1%
```
Contributor guide
Assessment
This issue has not been assessed yet.