bsts.R

source("common.R", encoding = "utf-8")

df$RP <- calc_RP(df$RP095, df$AP, .95)
z <- calc_Z(RP = df$RP, p = df$AP)
df <- mutate(df, z1=z$z1, z2=z$z2)
df <- mutate(df, z1E=z1*end, z2E=z2*end)

ss <- AddLocalLevel(list(), y = df$logPI) # c
ss <- AddAr(ss, lags=2, y = df$logPI) # AR(2)
# time-varying regression
ss <- AddDynamicRegression(ss, formula = logPI ~ z1 + z1E + z2 + z2E, data=df)

# calculate
model <- bsts(formula=df$logPI,
              state.specification = ss,
              niter = 500, seed=42)
plot(model) # it returns error if observed data contains NA value because `ylim` is specified incorrectly.
plot(model, "residuals")
plot(model, "prediction.errors")
plot(model, "components")
plot(model, "dynamic")