bsts.R

require(bsts) # 0.7.1

data(iclaims)     # bring the initial.claims data into scope

# --- model 1 ----
ss <- AddLocalLinearTrend(list(), initial.claims$iclaimsNSA)
ss <- AddSeasonal(ss, initial.claims$iclaimsNSA, nseasons = 52)
model1 <- bsts(initial.claims$iclaimsNSA,
               state.specification = ss,
               niter = 1000)
plot(model1, "comp")
plot(model1, "resid")

# forecasting
pred1 <- predict(model1, horizon = 12)
plot(pred1, plot.original = 156)


# ---- model 2 -----
# Fit a bsts model with expected model size 1, the default.
model2 <- bsts(iclaimsNSA ~ .,
               state.specification = ss,
               niter = 1000,
               data = initial.claims)
plot(model2, "comp")
plot(model2, "coef")
plot(model2, "resid")


# ----- model 3 -----
# Fit a bsts model with expected model size 5, to include more coefficients.
model3 <- bsts(iclaimsNSA ~ .,
               state.specification = ss,
               niter = 1000,
               data = initial.claims,
               expected.model.size = 5)  # Passed to SpikeSlabPrior.
plot(model3, "comp")
plot(model3, "coef")
plot(model3, "resid")

# ----- model 4 -----
ss4 <- AddAutoAr(ss, y = initial.claims$iclaimsNSA, lags=2)
model4 <- bsts(formula = iclaimsNSA ~ .,
               state.specification = ss4, niter=1000, data=initial.claims)
plot(model4, "comp")
plot(model4, "coef")
plot(model4, "resid")

# ----- model 5 -----
ss5 <- AddDynamicRegression(ss4, formula=iclaimsNSA ~ ., data = initial.claims)
model5 <- bsts(formula = initial.claims$iclaimsNSA,
               state.specification = ss5, niter=1000)
plot(model5, "comp")
plot(model5, "dyn")
plot(model5, "resid")

CompareBstsModels(list("Model 1" = model1,
                       "Model 2" = model2,
                       "Model 3" = model3,
                       "Model 4" = model4,
                       "Model 5" = model5),
                  colors = c("black", "red", "blue", "green", "purple"))