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

## -----------------------------------------------------------------------------
# kf = kalman_filter(ssm, yt, Xo, Xs, w, smooth)

## -----------------------------------------------------------------------------
# library(kalmanfilter)
# library(maxLik)
# library(ggplot2)
# library(data.table)
# library(gridExtra)
# data(treasuries)
# tau = unique(treasuries$maturity)
# 
# #State space model for the Nelson-Siegel dynamic factor yield curve
# yc_ssm = function(par, tau){
#   #Set factor names
#   tau = unique(tau)
#   factors = c("l", "s", "c")
# 
#   #Loading on the slope factor
#   Ls = function(par, tau){
#     return(-(1 - exp(-tau*par["lambda"]))/(tau*par["lambda"]))
#   }
# 
#   #Loading on the curvature factor
#   Lc = function(par, tau){
#     return(-Ls(par["lambda"], tau) - exp(-tau*par["lambda"]))
#   }
# 
#   ########## Define the state space model
#   #Transition matrix
#   Fm = matrix(par[grepl("phi_", names(par))], nrow = 3, ncol = 3,
#               dimnames = list(paste0(factors, "_t0"),
#                               paste0(factors, "_t1")))
#   Fm[is.na(Fm)] = 0
# 
#   #Transition intercept matrix
#   Dm = matrix(par[grepl("D_", names(par))], ncol = 1,
#               dimnames = list(rownames(Fm), NULL))
# 
#   #Transition error covariance matrix
#   Qm = matrix(par[unlist(lapply(factors, function(x){paste0("sig_", x, factors)}))], nrow = 3, ncol = 3)
#   Qm[lower.tri(Qm)] = Qm[upper.tri(Qm)]
#   dimnames(Qm) = list(rownames(Fm), rownames(Fm))
# 
#   #Observation equation matrix
#   n = length(tau)
#   Hm = matrix(NA, nrow = n, ncol = 3,
#               dimnames = list(as.character(tau), rownames(Fm)))
#   tau = as.numeric(tau)
#   if("l" %in% factors){
#     Hm[, "l_t0"] = 1
#   }
#   if("s" %in% factors){
#     Hm[, "s_t0"] = Ls(par, tau)
#   }
#   if("c" %in% factors){
#     Hm[, "c_t0"] = Lc(par, tau)
#   }
# 
#   #Observation error variance covariance matrix
#   Rm = diag(sum(grepl("sig_\\d+", names(par))))*par[grepl("sig_\\d+", names(par))]^2
#   dimnames(Rm) = list(rownames(Hm), rownames(Hm))
# 
#   #Observation intercept matrix
#   Am = matrix(0, nrow = nrow(Hm), ncol = 1,
#               dimnames = list(rownames(Hm), NULL))
# 
#   #Initial guess for the unobserved vector
#   B0 = matrix(par[names(par) %in% c("level", "slope", "curvature")], ncol = 1, nrow = nrow(Fm),
#               dimnames = list(rownames(Fm), NULL))
# 
#   #Initial guess for the variance of the unobserved vector
#   P0 = diag(par["P0"]^2, nrow = nrow(Qm), ncol = ncol(Qm))
# 
#   return(list(Fm = Fm, Dm = Dm, Qm = Qm,
#               Hm = Hm, Am = Am, Rm = Rm,
#               B0 = B0, P0 = P0))
# }
# 
# #Set the initial values
# init = c(level = 5.9030, slope = -0.7090, curvature = 0.8690, lambda = 0.0423,
#          D_l = 0.1234, D_s = -0.2285, D_c = 0.2020,
#          phi_ll = 0.9720, phi_ls = 0.1009, phi_lc = -0.1226,
#          phi_sl = -0.0209, phi_ss = 0.8189, phi_sc = 0.0192,
#          phi_cl = -0.0061, phi_cs = -0.1446, phi_cc = 0.8808,
#          sig_ll = 0.1017,
#          sig_sl = 0.0937, sig_ss = 0.2267,
#          sig_cl = 0.0303, sig_cs = 0.0351, sig_cc = 0.7964,
#          sig_1 = 0.0934, sig_3 = 0.0001, sig_6 = 0.1206, sig_12 = 0.1525,
#          sig_24 = 0.1328, sig_36 = 0.0855, sig_60 = 0.0001, sig_84 = 0.0397,
#          sig_120 = 0.0595, sig_240 = 0.1438, sig_360 = 0.1450,
#          P0 = 0.0001)
# 
# #Set up the constraints: lambda and all standard deviation parameters must be positive
# ineqA = matrix(0, nrow = 15, ncol = length(init), dimnames = list(NULL, names(init)))
# diag(ineqA[, c("lambda", "sig_ll", "sig_ss", "sig_cc", paste0("sig_", tau))]) = 1
# ineqB = matrix(0, nrow = nrow(ineqA), ncol = 1)
# 
# #Convert to an NxT matrix
# yt = dcast(treasuries, "date ~ maturity", value.var = "value")
# yt = t(yt[, 2:ncol(yt)])
# 
# #Set the objective function
# objective = function(par){
#   ssm = yc_ssm(par, tau)
#   return(kalman_filter(ssm, yt)$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))
# 
# #Get the estimated state space model
# ssm = yc_ssm(solve$estimate, tau)
# 
# #Filter the data with the model
# filter = kalman_filter(ssm, yt, smooth = TRUE)
# 
# #Get the estimated unobserved factors
# B_tt = data.table(t(filter[["B_tt"]]))[, "date" := unique(treasuries$date)]
# colnames(B_tt)[1:3] = c("level", "slope", "curvature")
# 
# #Get the estimated yields
# y_tt = data.table(t(filter[["y_tt"]]))[, "date" := unique(treasuries$date)]
# colnames(y_tt)[1:(ncol(y_tt) - 1)] = tau
# 
# #Plot the data
# g1 = ggplot(treasuries, id.vars = "date") +
#   geom_line(aes(x = date, y = value, group = factor(maturity), color = factor(maturity))) +
#   theme_minimal() + theme(legend.position = "bottom") +
#   labs(title = "Actual Treasury Yields", x = "", y = "Value") +
#   guides(color = guide_legend(title = NULL))
# g2 = ggplot(melt(y_tt, id.vars = "date")) +
#   geom_line(aes(x = date, y = value, group = variable, color = variable)) +
#   theme_minimal() + theme(legend.position = "bottom") +
#   labs(title = "Estimated Treasury Yields", x = "", y = "Value") +
#   guides(color = guide_legend(title = NULL))
# g3 = ggplot(melt(B_tt, id.vars = "date")) +
#   geom_line(aes(x = date, y = value, group = variable, color = variable)) +
#   theme_minimal() + theme(legend.position = "bottom") +
#   labs(title = "Estimated Factors", x = "", y = "Value") +
#   guides(color = guide_legend(title = NULL))
# grid.arrange(g1, g2, g3, layout_matrix = rbind(c(1, 2), c(3, 3)))

## -----------------------------------------------------------------------------
# #Set the objective function
# objective = function(par){
#   ssm = yc_ssm(par, tau)
#   for(x in names(ssm)){
#     if(!x %in% c("B0", "P0")){
#       ssm[[x]] = array(ssm[[x]], dim = c(dim(ssm[[x]]), ncol(yt)))
#     }
#   }
#   return(kalman_filter(ssm, yt)$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))

## -----------------------------------------------------------------------------
# #Note the outcome does not work as expected anymore
# library(kalmanfilter)
# 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 Dynamic Common Factor model
# dcf_ssm = function(par, yt){
#   #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))
# 
#   #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
# 
#   #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
# 
#   #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
#   }
# 
#   #Build the observation equation covariance matrix
#   Rm = matrix(0, ncol = nrow(Hm), nrow = nrow(Hm), dimnames = list(vars, vars))
# 
#   #Transition equation intercept matrix
#   Dm = matrix(0, nrow = nrow(Fm), ncol = 1, dimnames = list(rownames(Fm), NULL))
# 
#   #Observation equation intercept matrix
#   Am = matrix(0, nrow = nrow(Hm), ncol = 1)
# 
#   #Initialize the filter for each state
#   B0 = matrix(0, nrow(Fm), 1)
#   P0 = diag(nrow(Fm))
#   dimnames(B0) = list(rownames(Fm), NULL)
#   dimnames(P0) = list(rownames(Fm), rownames(Fm))
# 
#   B0 = solve(diag(ncol(Fm)) - Fm) %*% Dm
#   VecP_ll = solve(diag(prod(dim(Fm))) - kronecker(Fm, Fm)) %*% matrix(as.vector(Qm), ncol = 1)
#   P0 = 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))
# }
# 
# #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)]
# 
# #Standardize 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 = FALSE])
# 
# #Set the initial values
# init = c(phi1 = 0.8588, phi2 = -0.1526,
#          psi_1.1 = -0.1079, psi_1.2 = -0.1415,
#          psi_2.1 = -0.3166, psi_2.2 = -0.0756,
#          psi_3.1 = -0.3994, psi_3.2 = -0.2028,
#          psi_4.1 = -0.0370, psi_4.2 = 0.3646,
#          gamma_1.0 = 0.0064, gamma_1.1 = -0.0020,
#          gamma_2.0 = 0.0018,
#          gamma_3.0 = 0.0059, gamma_3.1 = -0.0027,
#          gamma_4.0 = 0.0014, gamma_4.1 = -0.0004, gamma_4.2 = 0.0001, gamma_4.3 = 0.0004,
#          sigma_1 = 0.0049, sigma_2 = 0.0056, sigma_3 = 0.0077, sigma_4 = 0.0013)
# 
# #Set the constraints
# ineqA = matrix(0, nrow = 14,
#                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
# ineqB = matrix(c(rep(1, 10),
#                  rep(0, 4)), nrow = nrow(ineqA), ncol = 1)
# 
# #Define the objective function
# objective = function(par, yt){
#   ssm = dcf_ssm(par, yt)
#   return(kalman_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 = dcf_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 coefficient matrices
# Hm = ssm[["Hm"]]
# Fm = ssm[["Fm"]]
# 
# #Final K_t is approximation to steady state K
# filter = kalman_filter(ssm, yt, smooth = T)
# 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)
# 
# #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))
# 
# g3 = ggplot(melt(Ctt, id.vars = "date")[variable == "dcf", ]) +
#   ggtitle("Dynamic Common Factor", subtitle = "Levels") +
#   scale_x_date(name = "") +
#   geom_hline(yintercept = 0, color = "grey") +
#   geom_line(aes(x = date, y = value, group = variable, color = variable)) +
#   theme_minimal() + theme(legend.position = "bottom") + scale_color_manual(values = "black") +
#   guides(color = guide_legend(title = NULL), fill = guide_legend(title = NULL)) +
#   scale_y_continuous(name = "Value", limits = range(Ctt$dcf, na.rm = TRUE))
# 
# g4 = ggplot(melt(Ctt, id.vars = "date")[variable == "d.dcf", ]) +
#   ggtitle("Dynamic Common Factor", subtitle = "Differenced") +
#   scale_x_date(name = "") +
#   geom_hline(yintercept = 0, color = "grey") +
#   geom_line(aes(x = date, y = value, group = variable, color = variable)) +
#   theme_minimal() + theme(legend.position = "bottom") + scale_color_manual(values = "black") +
#   guides(color = guide_legend(title = NULL), fill = guide_legend(title = NULL)) +
#   scale_y_continuous(name = "Value", limits = range(Ctt$d.dcf, na.rm = TRUE))
# 
# grid.arrange(g1, g2, g3, g4, layout_matrix = matrix(c(1, 3, 2, 4), nrow = 2))