Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment