March 16, 2017

Fit Model

library(plotly)
library(dlm)

buildFun <- function(x) {
     m <- dlmModPoly(1, dV = exp(x[1]))
     m$JW <- matrix(1)
     m$X <- matrix(exp(x[2]), nc = 1, nr = length(Nile))
     j <- which(time(Nile) == 1899)
     m$X[j,1] <- m$X[j,1] * (1 + exp(x[3]))
     return(m)
}
fit <- dlmMLE(Nile, parm = c(0,0,0), build = buildFun)
dlmNileJump <- buildFun(fit$par)
nileJumpFilt <- dlmFilter(Nile, dlmNileJump)

Calculate probability intervals

attach(nileJumpFilt)
v <- unlist(dlmSvd2var(U.C, D.C))

pl <- dropFirst(m) + qnorm(0.05, sd = sqrt(v[-1]))
pu <- dropFirst(m) + qnorm(0.95, sd = sqrt(v[-1]))
detach()


p <- plot_ly() %>%
  add_lines(x = time(Nile), y = Nile,
            color = I("black"), name = "observed") %>%
  add_ribbons(x = time(dropFirst(nileJumpFilt$m)), ymin = pl, ymax = pu,
              color = I("gray95"), name = "95% probability") %>%
  add_lines(x = time(dropFirst(nileJumpFilt$m)), y = dropFirst(nileJumpFilt$m), color = I("blue"), name = "filtered observation")

Go Bayesians!