kfas.R

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

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

# specify model
model3KFAS <- SSModel(logPI ~ SSMtrend(1, Q=NA) +
                        SSMarima(ar = c(.1, .1), Q = 0) +
                        SSMregression(~z1 + z1E + z2 + z2E, Q = diag(rep(NA, 4))),
                      H=NA, data = df, distribution = "gaussian")
print(model3KFAS)
# estimate unknown parameters
update.model <- function(pars, model){
  model$H[] <- exp(pars[1])
  diag(model$Q[, , 1])[1:6] <- exp(pars[2:7])
  model$T[6:7, 6, 1] <- tanh(pars[8:9])
  return(model)
}
fit.model3 <- fitSSM(model = model3KFAS, updatefn = update.model,
                     inits = numeric(9), maxit=1000000)
# convergence check (not converged...)
fit.model3$optim.out$convergence
# log likelihood
-fit.model3$optim.out$value
# smoothed time-varying parameterss
ggplot(data.frame(coef(KFS(fit.model3$model))) %>% mutate(t=1:n()) %>% gather(key=series, value=value, -t)) +
  geom_line(aes(x=t, y=value, group=series), color="grey") + facet_wrap(~series, scales = "free") +
  labs(x="week", y="") + theme_bw()

# ---- without AR(2) ----
model2 <- SSModel(logPI ~ SSMtrend(1, Q = NA) + SSMregression(~ z1 + z1E + z2 + z2E, Q=diag(rep(NA, 4))), data=df, H = NA)
fit.model <- fitSSM(model, inits=numeric(6), method="L-BFGS-B")
# convergence check
fit.model$optim.out$convergence
# log likelihood
-fit.model$optim.out$value
# smoothed time-varying parameterss
ggplot(data.frame(coef(KFS(fit.model$model))) %>% mutate(t=1:n()) %>% gather(key=series, value=value, -t)) +
  geom_line(aes(x=t, y=value, group=series), color="grey") + facet_wrap(~series, scales = "free") +
  labs(x="week", y="") + theme_bw()