## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  echo = TRUE, include = TRUE, message = FALSE, warning = FALSE, eval = FALSE
)

## -----------------------------------------------------------------------------
#  library(kimfilter)
#  library(data.table)
#  library(maxLik)
#  library(ggplot2)
#  library(gridExtra)
#  data(sw_dcf)
#  data = sw_dcf[, colnames(sw_dcf) != "dcoinc", with = F]
#  vars = colnames(data)[colnames(data) != "date"]
#  
#  #State space model for the Stock and Watson Markov Switching Dynamic Common Factor model
#  msdcf_ssm = function(par, yt, n_states = NULL){
#    #Get the number of states
#    n_states = length(unique(unlist(lapply(strsplit(names(par)[grepl("p_", names(par))], "p_"), function(x){substr(x[2], 1, 1)}))))
#  
#    #Get the parameters
#    vars = dimnames(yt)[which(unlist(lapply(dimnames(yt), function(x){!is.null(x)})))][[1]]
#    phi = par[grepl("phi", names(par))]
#    names(phi) = gsub("phi", "", names(phi))
#    gamma = par[grepl("gamma_", names(par))]
#    names(gamma) = gsub("gamma_", "", names(gamma))
#    psi = par[grepl("psi_", names(par))]
#    names(psi) = gsub("psi_", "", names(psi))
#    sig = par[grepl("sigma_", names(par))]
#    names(sig) = gsub("sigma_", "", names(sig))
#    mu = par[grepl("mu", names(par))]
#    names(mu) = gsub("mu_", "", names(mu))
#    pr = par[grepl("p_", names(par))]
#    names(pr) = gsub("p_", "", names(pr))
#    states = sort(unique(substr(names(pr), 1, 1)))
#  
#    #Steady state probabilities
#    Pm = matrix(NA, nrow = n_states, ncol = n_states)
#    rownames(Pm) = colnames(Pm) = unique(unlist(lapply(names(pr), function(x){strsplit(x, "")[[1]][2]})))
#    for(j in names(pr)){
#      Pm[strsplit(j, "")[[1]][2], strsplit(j, "")[[1]][1]] = pr[j]
#    }
#    for(j in 1:ncol(Pm)){
#      Pm[which(is.na(Pm[, j])), j] = 1 - sum(Pm[, j], na.rm = TRUE)
#    }
#  
#    #Build the transition matrix
#    phi_dim = max(c(length(phi)), length(unique(sapply(strsplit(names(gamma), "\\."), function(x){x[2]}))))
#    psi_dim = sapply(unique(sapply(strsplit(names(psi), "\\."), function(x){x[1]})), function(x){
#      max(as.numeric(sapply(strsplit(names(psi)[grepl(paste0("^", x), names(psi))], "\\."), function(x){x[2]})))
#    })
#    Fm = matrix(0, nrow = phi_dim + length(psi), ncol = phi_dim + length(psi),
#                dimnames = list(
#                  c(paste0("ct.", 0:(phi_dim - 1)),
#                    unlist(lapply(names(psi_dim), function(x){paste0("e_", x, ".", 0:(psi_dim[x] - 1))}))),
#                  c(paste0("ct.", 1:phi_dim),
#                    unlist(lapply(names(psi_dim), function(x){paste0("e_", x, ".", 1:psi_dim[x])})))
#                ))
#    Fm["ct.0", paste0("ct.", names(phi))] = phi
#    for(i in 1:length(vars)){
#      Fm[paste0("e_", i, ".0"),
#         paste0("e_", names(psi)[grepl(paste0("^", i), names(psi))])] = psi[grepl(paste0("^", i), names(psi))]
#    }
#    diag(Fm[intersect(rownames(Fm), colnames(Fm)), intersect(rownames(Fm), colnames(Fm))]) = 1
#    Fm = array(Fm, dim = c(nrow(Fm), ncol(Fm), n_states), dimnames = list(rownames(Fm), colnames(Fm), states))
#  
#    #Build the observation matrix
#    Hm = matrix(0, nrow = nrow(yt), ncol = nrow(Fm), dimnames = list(rownames(yt), rownames(Fm)))
#    for(i in 1:length(vars)){
#      Hm[i, paste0("ct.", sapply(strsplit(names(gamma)[grepl(paste0("^", i), names(gamma))], "\\."), function(x){x[2]}))] =
#        gamma[grepl(paste0("^", i), names(gamma))]
#    }
#    diag(Hm[, paste0("e_", 1:length(vars), ".0")]) = 1
#    Hm = array(Hm, dim = c(nrow(Hm), ncol(Hm), n_states), dimnames = list(rownames(Hm), colnames(Hm), states))
#  
#    #Build the state covariance matrix
#    #Set the dynamic common factor standard deviation to 1
#    Qm = matrix(0, ncol = ncol(Fm), nrow = nrow(Fm), dimnames = list(rownames(Fm), rownames(Fm)))
#    Qm["ct.0", "ct.0"] = 1
#    for(i in 1:length(vars)){
#      Qm[paste0("e_", i, ".0"), paste0("e_", i, ".0")] = sig[names(sig) == i]^2
#    }
#    Qm = array(Qm, dim = c(nrow(Qm), ncol(Qm), n_states), dimnames = list(rownames(Qm), colnames(Qm), states))
#  
#    #Build the observation equation covariance matrix
#    Rm = matrix(0, ncol = nrow(Hm), nrow = nrow(Hm), dimnames = list(vars, vars))
#    Rm = array(Rm, dim = c(nrow(Rm), ncol(Rm), n_states), dimnames = list(rownames(Rm), colnames(Rm), states))
#  
#    #State intercept matrix: the Markov switching mean matrix
#    Dm = matrix(0, nrow = nrow(Fm), ncol = 1, dimnames = list(rownames(Fm), NULL))
#    Dm = array(Dm, dim = c(nrow(Dm), 1, n_states), dimnames = list(rownames(Fm), NULL, states))
#    for(i in names(mu)){
#      Dm["ct.0", , i] = mu[i]
#    }
#  
#    #Observation equation intercept matrix
#    Am = matrix(0, nrow = nrow(Hm), ncol = 1)
#    Am = array(Am, dim = c(nrow(Am), ncol(Am), n_states), dimnames = list(vars, NULL, states))
#  
#    #Initialize the filter for each state
#    B0 = matrix(0, nrow(Fm), 1)
#    P0 = diag(nrow(Fm))
#    B0 = array(B0, dim = c(nrow(B0), ncol(B0), n_states), dimnames = list(rownames(Fm), NULL, states))
#    P0 = array(P0, dim = c(nrow(P0), ncol(P0), n_states), dimnames = list(rownames(B0), colnames(B0), states))
#    for(i in states){
#      B0[,,i] = solve(diag(ncol(Fm)) - Fm[,,i]) %*% Dm[,,i]
#      VecP_ll = solve(diag(prod(dim(Fm[,,i]))) - kronecker(Fm[,,i], Fm[,,i])) %*% matrix(as.vector(Qm[,,i]), ncol = 1)
#      P0[,,i] = matrix(VecP_ll[, 1], ncol = ncol(Fm))
#    }
#  
#    return(list(B0 = B0, P0 = P0, Am = Am, Dm = Dm, Hm = Hm, Fm = Fm, Qm = Qm, Rm = Rm, Pm = Pm))
#  }
#  
#  #Log the data
#  data.log = copy(data)
#  data.log[, c(vars) := lapply(.SD, log), .SDcols = c(vars)]
#  
#  #Difference the data
#  data.logd = copy(data.log)
#  data.logd[, c(vars) := lapply(.SD, function(x){
#    x - shift(x, type = "lag", n = 1)
#  }), .SDcols = c(vars)]
#  
#  #Center the data
#  data.logds = copy(data.logd)
#  data.logds[, c(vars) := lapply(.SD, scale, scale = FALSE), .SDcols = c(vars)]
#  
#  #Transpose the data
#  yt = t(data.logds[, c(vars), with = F])
#  
#  #Set the initial values
#  init = c(phi1 = 0.8760, phi2 = -0.2171,
#           mu_u = 0.2802, mu_d = -1.5700,
#           p_dd = 0.8406, p_uu = 0.9696,
#           psi_1.1 = 0.0364, psi_1.2 = -0.0008,
#           psi_2.1 = -0.2965, psi_2.2 = -0.0657,
#           psi_3.1 = -0.3959, psi_3.2 = -0.1903,
#           psi_4.1 = -0.2436, psi_4.2 = 0.1281,
#           gamma_1.0 = 0.0058, gamma_1.1 = -0.0033,
#           gamma_2.0 = 0.0011,
#           gamma_3.0 = 0.0051, gamma_3.1 = -0.0033 ,
#           gamma_4.0 = 0.0012, gamma_4.1 = -0.0005, gamma_4.2 = 0.0001, gamma_4.3 = 0.0002,
#           sigma_1 = 0.0048, sigma_2 = 0.0057, sigma_3 = 0.0078, sigma_4 = 0.0013)
#  
#  #Set the constraints
#  ineqA = matrix(0, nrow = 20, ncol = length(init), dimnames = list(NULL, names(init)))
#  #Stationarity constraints
#  ineqA[c(1, 2), c("phi1", "phi2")] = rbind(c(1, 1), c(-1, -1))
#  ineqA[c(3, 4), grepl("psi_1", colnames(ineqA))] = rbind(c(1, 1), c(-1, -1))
#  ineqA[c(5, 6), grepl("psi_2", colnames(ineqA))] = rbind(c(1, 1), c(-1, -1))
#  ineqA[c(7, 8), grepl("psi_3", colnames(ineqA))] = rbind(c(1, 1), c(-1, -1))
#  ineqA[c(9, 10), grepl("psi_4", colnames(ineqA))] = rbind(c(1, 1), c(-1, -1))
#  #Non-negativity constraints
#  diag(ineqA[c(11, 12, 13, 14), grepl("sigma_", colnames(ineqA))]) = 1
#  ineqA[c(15, 16), "p_dd"] = c(1, -1)
#  ineqA[c(17, 18), "p_uu"] = c(1, -1)
#  #Up/down states must be positive/negative
#  ineqA[19, "mu_u"] = 1
#  ineqA[20, "mu_d"] = -1
#  ineqB = matrix(c(rep(1, 10),
#                   rep(0, 4),
#                   c(0, 1),
#                   c(0, 1),
#                   rep(0, 2)), nrow = nrow(ineqA), ncol = 1)
#  
#  #Define the objective function
#  objective = function(par, yt){
#    ssm = msdcf_ssm(par, yt)
#    return(kim_filter(ssm, yt, smooth = FALSE)$lnl)
#  }
#  
#  #Solve the model
#  solve = maxLik(logLik = objective, start = init, method = "BFGS",
#                 finalHessian = FALSE, hess = NULL,
#                 control = list(printLevel = 2, iterlim = 10000),
#                 constraints = list(ineqA = ineqA, ineqB = ineqB),
#                 yt = yt)
#  
#  #Get the estimated state space model
#  ssm = msdcf_ssm(solve$estimate, yt)
#  
#  #Get the column means and standard deviations
#  M = matrix(unlist(data.logd[, lapply(.SD, mean, na.rm = TRUE), .SDcols = c(vars)]),
#                 ncol = 1, dimnames = list(vars, NULL))
#  
#  #Get the steady state coefficient matrices
#  Pm = matrix(ss_prob(ssm[["Pm"]]), ncol = 1, dimnames = list(rownames(ssm[["Pm"]]), NULL))
#  Hm = Reduce("+", lapply(dimnames(ssm[["Hm"]])[[3]], function(x){
#    Pm[x, ]*ssm[["Hm"]][,, x]
#  }))
#  Fm = Reduce("+", lapply(dimnames(ssm[["Fm"]])[[3]], function(x){
#    Pm[x, ]*ssm[["Fm"]][,, x]
#  }))
#  
#  #Final K_t is approximation to steady state K
#  filter = kim_filter(ssm, yt, smooth = TRUE)
#  K = filter$K_t[,, dim(filter$K_t)[3]]
#  W = solve(diag(nrow(K)) - (diag(nrow(K)) - K %*% Hm) %*% Fm) %*% K
#  d = (W %*% M)[1, 1]
#  
#  #Get the intercept terms
#  gamma = Hm[, grepl("ct", colnames(Hm))]
#  D = M - gamma %*% matrix(rep(d, ncol(gamma)))
#  
#  #Initialize first element of the dynamic common factor
#  Y1 = t(data.log[, c(vars), with = F][1, ])
#  initC = function(par){
#    return(sum((Y1 - D - gamma %*% par)^2))
#  }
#  C10 = optim(par = Y1, fn = initC, method = "BFGS", control = list(trace = FALSE))$par[1]
#  Ctt = rep(C10, ncol(yt))
#  
#  #Build the rest of the level of the dynamic common factor
#  ctt = filter$B_tt[which(rownames(Fm) == "ct.0"), ]
#  for(j in 2:length(Ctt)){
#    Ctt[j] = ctt[j] + Ctt[j - 1] + c(d)
#  }
#  Ctt = data.table(date = data$date, dcf = Ctt, d.dcf = ctt)
#  prob = data.table(date = data$date, data.table(filter$Pr_tt))
#  colnames(prob) = c("date", paste0("pr_", dimnames(ssm$Dm)[[3]]))
#  uc = merge(Ctt, prob, by = "date", all = TRUE)
#  
#  #Plot the outputs
#  g1 = ggplot(melt(data.log, id.vars = "date")[, "value" := scale(value), by = "variable"]) +
#    ggtitle("Data", subtitle = "Log Levels (Rescaled)") +
#    scale_y_continuous(name = "Value") +
#    scale_x_date(name = "") +
#    geom_line(aes(x = date, y = value, group = variable, color = variable)) +
#    theme_minimal() + theme(legend.position = "bottom") + guides(color = guide_legend(title = NULL))
#  
#  g2 = ggplot( melt(data.logds, id.vars = "date")) +
#    ggtitle("Data", subtitle = "Log Differenced & Standardized") +
#    scale_y_continuous(name = "Value") +
#    scale_x_date(name = "") +
#    geom_hline(yintercept = 0, color = "black") +
#    geom_line(aes(x = date, y = value, group = variable, color = variable)) +
#    theme_minimal() + theme(legend.position = "bottom") + guides(color = guide_legend(title = NULL))
#  
#  toplot3 = melt(uc, id.vars = "date")
#  d_range1 = range(toplot3[variable == "dcf", ]$value, na.rm = TRUE)
#  p_range1 = range(toplot3[variable %in% colnames(uc)[grepl("pr_", colnames(uc))], ]$value, na.rm = TRUE)
#  toplot3[variable %in% colnames(uc)[grepl("pr_", colnames(uc))], "value" := (value - p_range1[1])/diff(p_range1) * diff(d_range1) + d_range1[1], by = "variable"]
#  g3 = ggplot() +
#    ggtitle("Dynamic Common Factor", subtitle = "Levels") +
#    scale_x_date(name = "") +
#    geom_hline(yintercept = 0, color = "grey") +
#    geom_line(data = toplot3[variable == "dcf", ],
#              aes(x = date, y = value, group = variable, color = variable)) +
#    theme_minimal() + theme(legend.position = "bottom") +
#    guides(color = guide_legend(title = NULL), fill = guide_legend(title = NULL)) +
#    scale_color_manual(values = "black") +
#    scale_y_continuous(name = "Value", limits = range(toplot3[variable == "dcf", ]$value, na.rm = TRUE),
#                       sec.axis = sec_axis(name = "Probability", ~((. - d_range1[1])/diff(d_range1) * diff(p_range1) + p_range1[1]) * 100)) +
#    geom_ribbon(data = toplot3[variable %in% "pr_d", ],
#                aes(x = date, ymin = d_range1[1], ymax = value, group = variable, fill = variable), alpha = 0.5) +
#    scale_fill_manual(values = c("red", "green"))
#  
#  toplot4 = melt(uc, id.vars = "date")
#  d_range2 = range(toplot4[variable %in% c("d.dcf"), ]$value, na.rm = TRUE)
#  p_range2 = range(toplot4[variable %in% colnames(uc)[grepl("pr_", colnames(uc))], ]$value, na.rm = TRUE)
#  toplot4[variable %in% colnames(uc)[grepl("pr_", colnames(uc))], "value" := (value - p_range2[1])/diff(p_range2) * diff(d_range2) + d_range2[1], by = "variable"]
#  g4 = ggplot() +
#    ggtitle("Dynamic Common Factor", subtitle = "Differenced") +
#    scale_x_date(name = "") +
#    geom_hline(yintercept = 0, color = "grey") +
#    geom_line(data = toplot4[variable %in% c("d.dcf"), ],
#              aes(x = date, y = value, group = variable, color = variable)) +
#    theme_minimal() + theme(legend.position = "bottom") +
#    guides(color = guide_legend(title = NULL), fill = guide_legend(title = NULL)) +
#    scale_color_manual(values = "black") +
#    scale_y_continuous(name = "Value", limits = range(toplot4[variable %in% c("d.dcf"), ]$value, na.rm = TRUE),
#                                 sec.axis = sec_axis(name = "Probability", ~((. - d_range2[1])/diff(d_range2) * diff(p_range2) + p_range2[1]) * 100)) +
#    geom_ribbon(data = toplot4[variable %in% "pr_d", ],
#                aes(x = date, ymin = d_range2[1], ymax = value, group = variable, fill = variable), alpha = 0.5) +
#    scale_fill_manual(values = c("red", "green"))
#  
#  grid.arrange(g1, g2, g3, g4, layout_matrix = matrix(c(1, 3, 2, 4), nrow = 2))