In davidtedfordholt/fable.bsts: Bayesian Structural Time Series Modelling Interface for 'fable'

knitr::opts_chunk$set(echo = TRUE)
library(bsts)     # load the bsts package
library(tidyverse)
library(tsibble)
library(tsbox)
library(fable)
library(fable.bsts)

This is a mock-up of how fable-style formulas could be used to express state specifications for bsts models. It includes code from Steven Scott's article "Fitting Bayesian structural time series with the bsts R package."

This Rmd does not knit!

I'm using the following functions as formula "specials": intercept ar(lags = NULL) level(type = c("local", "shared")) trend(type = c("locallinear", "semilocallinear", "studentlocallinear"), ...) season(type = c("regression", "trig", "monthlyannual"), ...) xreg(data, ...) holiday(holidays, ...)

This allows the specification of the state through a formula such as response ~ trend("level") + season(24) + xreg(.) + holiday(holidays)

Example 1: Nowcasting

# DATA
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)

# FORECAST MODEL 1
pred1 <- predict(model1, horizon = 12)

# MODEL 2
model2 <- bsts(iclaimsNSA ~ .,
               state.specification = ss,
               niter = 1000,
               data = initial.claims)

# MODEL 3
model3 <- bsts(iclaimsNSA ~ .,
               state.specification = ss,
               niter = 1000,
               data = initial.claims,
               expected.model.size = 5)  # Passed to SpikeSlabPrior.

# DATA - fable
# code that would transform iclaims into a wide tsibble
iclaims <- ts_tsibble(initial.claims) %>% 
  spread(key = id, value = value)

# MODELS 1, 2, 3 - fable
mbl <- model(iclaims,
             bsts_1 = BSTS(iclaimsNSA ~ trend("local") + season(52), niter = 1000),
             bsts_2 = BSTS(iclaimsNSA ~ trend("local") + season(52) + xreg(.)),
             bsts_3 = BSTS(iclaimsNSA ~ trend("local") + season(52) + xreg(., expected_size = 5)))

# FORECAST MODEL 1 - fable
fbl <- mbl %>%
  select(-bsts_2, -bsts_3) %>%     # only because models 2 and 3 aren't modeled in the article
  forecast(mbl, h = 12)

Example 2: Long term forecasting

We assume data is set up properly for each. The mbl (model table) could be passed automatically to the forecast() function in the fable version, if the only object needed at the end was the table of forecasts.

# MODEL 1 
ss1 <- AddLocalLinearTrend(list(), sp500)
model1 <- bsts(sp500, state.specification = ss1, niter = 1000)
pred1 <- predict(model1, horizon = 360)

# MODEL 2
ss2 <- AddSemilocalLinearTrend(list(), sp500)
model2 <- bsts(sp500, state.specification = ss2, niter = 1000)
pred2 <- predict(model2, horizon = 360)

# MODELS 1, 2 - fable
mbl <- model(sp500,
             bsts_1 = BSTS(iclaimsNSA ~ trend("local"), niter = 1000),
             bsts_2 = BSTS(iclaimsNSA ~ trend("semilocal"), niter = 1000))
fbl <- forecast(mbl, h = 360)

Example 3: Recession modeling using non-Gaussian data

# MODEL
ss <- AddLocalLevel(list(),
                    sigma.prior = SdPrior(sigma.guess = .1,
                                          sample.size = 1,
                                          upper.limit = 1),
                    initial.state.prior = NormalPrior(0, 5))
ts.model <- bsts(nber ~ ., ss, data = gdp, niter = 20000,
                 family = "logit", expected.model.size = 10)

# MODEL - fable
mbl <- model(gdp,
  bsts = BSTS(
    nber ~ trend("level",
                 sigma_prior = SdPrior(sigma_guess = .1, sample_size = 1, upper_limmit = 1),
                 initial_state_prior = NormalPrior(0, 5)) + 
      xreg(., expected_size = 10), 
    niter = 20000, family = "logit"))

Conclusion

# bsts
ss <- AddSeasonal(ss, y, nseasons = 24)
# fable 1
mbl <- model(data, bsts = BSTS(value ~ season(24)))
# fable 2
mbl <- model(data, bsts = BSTS(value ~ season("1 day")))

# bsts
ss <- AddSeasonal(ss, y, nseasons = 7, season.duration = 24)
# fable 1
mbl <- model(data, bsts = BSTS(value ~ season(7, duration = 24)))
# fable 2
mbl <- model(data, bsts = BSTS(value ~ season("1 week", duration = "1 day")))

davidtedfordholt/fable.bsts documentation built on Sept. 15, 2020, 11:38 a.m.