dlm.R

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()