steve-the-bayesian / BOOM

A C++ library for Bayesian modeling, mainly through Markov chain Monte Carlo, but with a few other methods supported. BOOM = "Bayesian Object Oriented Modeling". It is also the sound your computer makes when it crashes.
GNU Lesser General Public License v2.1
35 stars 14 forks source link

Wouldn't the predict function from the bsts package return the same values as the state contributions for missing y values? #75

Open maferre1ra opened 1 year ago

maferre1ra commented 1 year ago

I am comparing the estimates obtained from running a bsts model on a variable with the last 4 observations missing (held out) to what we can obtain from the predict function using the same model (horizon = 4), and the two set of estimates do not match.

Could anyone help shed some light?

Here is a reproducible example.

############################ library(bsts) library(CausalImpact) library(tidyverse) library(lubridate)

n = 33 # observations

GENERATE SOME DATA

y2 = rnorm(n, 1542429, 922877.8) post.period <- c(30, 33) # hold out 4 observations post.period.response <- y2[post.period[1] : post.period[2]] y2[post.period[1] : post.period[2]] <- NA

################

RUN A BSTS MODEL

niter = 15000 prior.level.sd = 0.1 sample.size = nr/10

sdy <- sd(y2, na.rm = TRUE) ss <- list() sd.prior <- SdPrior(sigma.guess = prior.level.sd * sdy, upper.limit = sdy, sample.size = sample.size) ss <- AddSeasonal(ss, y, nseasons = 12, season.duration = 1)

p = 1 ss2 <- bsts::AddAr(ss, y2, lags = p, sigma.prior = sd.prior)

local2 <- bsts(y2, state.specification = ss2, niter = niter, seed = 1, ping = 0, model.options = BstsOptions(save.prediction.errors = TRUE))

############################

GET STATE CONTRIBUTIONS

bsts.model = local2 burn <- SuggestBurn(0.1, bsts.model) state.contributions <- bsts.model$state.contributions[-seq_len(burn), , , drop = FALSE] state.samples <- rowSums(aperm(state.contributions, c(1, 3, 2)), dims = 2) point.pred.mean <- colMeans(state.samples)

##############

GET PREDITION

bsts_predict <- predict.bsts(local2, horizon = 4)

From the example above, the point predictions from state contributions do not match the 4 point predictions from the predict function.