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.
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.