Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
source("common.R", encoding = "utf-8")
# ---- dlm ------
# model 3
res <- data.frame()
for( a in seq(from=.1, to=.95, by=.05) ){
RP <- calc_RP(RP = df$RP095, p = df$AP, a = a)
z <- calc_Z(RP = RP, p = df$AP)
z1 <- z$z1
z2 <- z$z2
z1E <- z$z1 * df$end
z2E <- z$z2 * df$end
# create the design matrix
regressors <- cbind(z1, z1E, z2, z2E)
# model settings
# params[1] = obs. noise variance
# params[2] = local trend noise variance
# params[3:6] = noises variance of alpha_1, beta_1, alpha_2, beta_2
# params[7:8] = AR coefs
# params[9] = AR(2) noise variance
buildModel3 <- function(params){
model <- dlmModReg(regressors, addInt = F,
dV = exp(params[1]), dW = exp(params[3:6])) + # time-varying regression
dlmModPoly(order=1, dW = exp(params[2])) + # time-varying intercept
dlmModARMA(ar=c(tanh(params[7]), tanh(params[8])), sigma=exp(params[9])) # AR(2) component
return(model)
}
# estimate structural params
e <- try({res.mle <- dlmMLE(df$logPI, rep(0, 9), buildModel3, method = "L-BFGS-B", hessian=T,
control = list(maxit = 1000000, trace=1))
})
if (class(e) == "try-error") {
}
else {
res <- bind_rows(res,
data.frame(
a=a,
conv=res.mle$convergence, loglike=-res.mle$value,
AIC=2*res.mle$value+2*10,
matrix(res.mle$par, nrow=1))
)
}
}
res[,5:10] <- exp(res[,5:10])
res[,11:12] <- tanh(res[,11:12])
res[,13] <- exp(res[,13])
(result <- filter(res, conv == 0) %>% arrange(desc(loglike)))
a <- result$a[1]
RP <- calc_RP(RP = df$RP095, p = df$AP, a = a)
z <- calc_Z(RP = RP, p = df$AP)
# create the design matrix
x<-cbind(z$z1, z$z2, z$z1 * df$end, z$z2 * df$end)
parameters <- result[1,5:ncol(res)] %>% as.numeric
print(paste("AIC=", round(2*result$loglike[1]+2*(ncol(res) - 4), 2)))
# filtering
mod3 <- buildModel3(parameters)
mod.filt3 <- dlmFilter(df$logPI, mod3)
# smoothing
mod.smooth3 <- dlmSmooth(mod.filt3)
# plot smoothed parameters
df.params <- data.frame(dropFirst(mod.smooth3$s[,1:7]))
colnames(df.params) <- c("alpha1", "beta1", "alpha2", "beta2", "c", "phi1", "phi2")
df.params$t <- 1:nrow(df.params)
df.params <- gather(df.params, key=series, value=value, -t)
ggplot(df.params) + geom_line(aes(x=t, y=value, group=series), color="grey") +
facet_wrap(~series, scales = "free") + scale_color_discrete(guide=F) +
labs(x="week", y="") + theme_bw()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment