library("eDMA") ## This code reproduces the results of ## ## Catania and Nonejad: ## 'Dynamic Model Averaging for Practitioners in Economics and Finance: ## The eDMA Package'. ## ## The structure follows the related article. It is divided in ## sections. The sections can be extecuded independently, however, the ## initial part of the code 'Starting' has to be always executed. The ## script reproduces the table as well as the figures of the paper. ############################################# ## Starting ############################################# ## Path where all the results are saved sDir <- "./" ## Number of cores ## The number of cores affects the output of Section 5. Computational challenges ## The number of cores is bounded for some part of the analysis where the required ## computational power is low hence controlling for parallel efficiency. iNcores <- 8 ################## Create Folders ################## ## general folder sDir_Tables <- paste(sDir, "Tables/", sep = "") sDir_Figures <- paste(sDir, "Figures/", sep = "") dir.create(sDir_Tables, showWarnings = FALSE) dir.create(sDir_Figures, showWarnings = FALSE) ## simulation folder: Section 4.1 sDir_Sim <- paste(sDir, "simulation/", sep = "") dir.create(sDir_Sim, showWarnings = FALSE) ## Computational challenges : Section 5 sDir_Comp <- paste(sDir, "computational/", sep = "") dir.create(sDir_Comp, showWarnings = FALSE) ## folder empirical dir.create(paste(sDir, "empirical/", sep = ""), showWarnings = FALSE) sDir_h1 <- paste(sDir, "empirical/h1/", sep = "") sDir_h5 <- paste(sDir, "empirical/h5/", sep = "") dir.create(sDir_h1, showWarnings = FALSE) dir.create(sDir_h5, showWarnings = FALSE) ############################################# ## The Model : Section 2 ############################################# ## FIGURES data("USRecessions", package = "eDMA") iBurnPeriod <- 32 vBurnIn <- 1:iBurnPeriod rgb <- c("#D33F6A4D", "#E2E2E24D", "#4A6FE34D", "#6969694D") data("USData", package = "eDMA") dim(USData) head(round(USData[, 1:6], 2)) AR2_formula <- GDPDEF ~ Lag(GDPDEF, 1) + Lag(GDPDEF, 2) AR2_Reg_formula <- GDPDEF ~ Lag(GDPDEF, 1) + Lag(GDPDEF, 2) + Lag(UNEMP, 1) ## Do recursive ols iT <- nrow(USData) ## we write iBurnPeriod + 2 due to the fact that we use two lag of GDPDEF lFit_AR2 <- lapply((iBurnPeriod + 2):iT, function(i, formula, data) lm(formula, data = data[1:i]), data = USData, formula = AR2_formula) lFit_AR2_Reg <- lapply((iBurnPeriod + 2):iT, function(i, formula, data) lm(formula, data = data[1:i]), data = USData, formula = AR2_Reg_formula) ## a) inflation rate vGDPDEF <- USData[(iBurnPeriod + 2):iT, "GDPDEF"] ## b) sum of ar2 coefficients vSumCoef <- sapply(lFit_AR2, function(x) sum(coef(x)[2:3]) ) ## c) recurive beta estimate of unemployment for an AR(2)-unemployment model vBetaUnemp <- sapply(lFit_AR2_Reg, function(x) coef(x)[4] ) ## d) recursive pval t-test H0 : beta_unemployment = 0 vPval <- sapply(lFit_AR2_Reg, function(x) summary(x)$coefficients[4, "Pr(>|t|)"] ) ## Figure 1 vDates <- as.Date(index(vGDPDEF)) pdf(paste(sDir_Figures, "Figure0.pdf", sep = ""), width = 15, height = 10) layout(matrix(1:4, 2, 2), heights = c(2.0, 2.5, 2.0, 2.5)) par(mar = c(0, 4, 0.1, 2)) plot(vDates, vGDPDEF, type = "n", xaxt = "n", xlab = "", yaxt = "n", ylab = "", ylim = c(-1.5, 4.0)) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, vGDPDEF, lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) axis(2, at = c(-1, 0, 1, 2, 3, 4) , labels = c("-1.00", "0.00", "1.00", "2.00", "3.00", "4.00"), las = 1, cex.axis = 1.5) axis(4, at = sum(c(-1.5, 4.0))/2 , labels = "(a)", tick = FALSE, las = 1, hadj = 0.5, cex.axis = 1.5) par(mar = c(3, 4, 0, 2)) plot(vDates, vSumCoef, type = "n", xaxt = "n", xlab = "", ylab = "", yaxt = "n", ylim = c(0.8, 1.1)) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, vSumCoef, lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) axis(2, at = seq(0.8, 1.1, length.out = 6), labels = c("0.80", "0.86", "0.92", "0.98", "1.04", "1.10"), las = 1, cex.axis = 1.5) axis(4, at = (max(vSumCoef) + min(vSumCoef))/2 , labels = "(b)", tick = FALSE, las = 1, hadj = 0.5, cex.axis = 1.5) axis.Date(1, at = seq(min(vDates), max(vDates), "year"), cex.axis = 1.5) axis.Date(1, at = seq(min(vDates), max(vDates), "quarter"), labels = FALSE, tcl = -0.2) par(mar = c(0, 4, 0.1, 2)) plot(vDates, vBetaUnemp, type = "n", yaxt = "n", xaxt = "n", xlab = "", ylab = "", ylim = c(-0.5, 0.0)) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, vBetaUnemp, lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) axis(2, at = seq(-0.5, 0.0, length.out = 6), labels = c("-0.50", "-0.40", "-0.30", "-0.20", "-0.10", "0.00"), las = 1, cex.axis = 1.5) axis(4, at = (max(vBetaUnemp) + min(vBetaUnemp))/2 , labels = "(c)", tick = FALSE, las = 1, hadj = 0.5, cex.axis = 1.5) par(mar = c(3, 4, 0, 2)) plot(vDates, vPval, type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "") grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, vPval, lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) axis(2, at = seq(0.0, 0.8, length.out = 6), labels = c("0.00", "0.16", "0.32", "0.48", "0.64", "0.80"), las = 1, cex.axis = 1.5) axis(4, at = (max(vPval) + min(vPval))/2 , labels = "(d)", tick = FALSE, las = 1, hadj = 0.5, cex.axis = 1.5) axis.Date(1, at = seq(min(vDates), max(vDates), "year"), cex.axis = 1.5) axis.Date(1, at = seq(min(vDates), max(vDates), "quarter"), labels = FALSE, tcl = -0.2) dev.off() ############################################# ## Simulations : Section 4.1 ############################################# set.seed(7892) iT <- 500 iK <- 3 dV <- 0.1 mW <- diag(iK + 1) * 0.01 dPhi <- 1 vBeta0 <- rep(0, iK + 1) mX <- cbind(1, matrix(rnorm(iT * (iK)), iT, iK)) lOut <- SimulateDLM(iT, mX, vBeta0, mW, dV, dPhi) vY <- lOut$vY mX <- mX[, -1] iK_Add <- 2 mX_add <- matrix(rnorm(iT * iK_Add), iT, iK_Add) SimData <- cbind(y = vY, mX, mX_add) colnames(SimData) <- c("y", paste("x", 2:(iK + iK_Add + 1), sep = "")) ## Alternatively: ## data("SimData", package = "eDMA") Fit <- DMA(y ~ x2 + x3 + x4 + x5 + x6, data = SimData, vDelta = seq(0.9, 1, 0.01)) Fit save(Fit, file = paste(sDir_Sim, "/SimulatedFit.Rdata", sep = "")) PostProb <- as.data.frame(Fit, which = "mincpmt", iBurnPeriod = 50) round(tail(PostProb), 2) if (interactive()) { plot(Fit) } plot(Fit, which = "mincpmt", iBurnPeriod = 50) ## Figure 2 iBurnPeriod <- 50 mincpmt <- as.data.frame(Fit, which = "mincpmt", iBurnPeriod = iBurnPeriod) vLim <- c(0, 1) pdf(paste(sDir_Figures, "Figure1.pdf", sep = ""), width = 15, height = 10) layout(matrix(1:6, 3, 2), heights = c(rep(2, 2), 2.5, rep(2, 2), 2.5)) par(mar = c(0, 4, 0.1, 2)) plot(mincpmt[, 1], type = "n", xaxt = "n", xlab = "", ylab = "", las = 1, ylim = vLim, cex.axis = 1.5) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(mincpmt[, 1], col = "black", lwd = 2) axis(4, at = 0.5, labels = "(a)", tick = FALSE, las = 1, hadj = 0.5, cex.axis = 1.5) legend("bottomright", legend = expression(x[1]), bg = "white", cex = 2) par(mar = c(0, 4, 0, 2)) plot(mincpmt[, 2], type = "n", xaxt = "n", xlab = "", ylab = "", las = 1, ylim = vLim, cex.axis = 1.5) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(mincpmt[, 2], col = "black", lwd = 2) axis(4, at = 0.5, labels = "(b)", tick = FALSE, las = 1, hadj = 0.5, cex.axis = 1.5) legend("bottomright", legend = expression(x[2]), bg = "white", cex = 2) par(mar = c(3, 4, 0, 2)) plot(mincpmt[, 3], type = "n", xaxt = "n", xlab = "", ylab = "", las = 1, ylim = vLim, cex.axis = 1.5) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(mincpmt[, 3], col = "black", lwd = 2) axis(4, at = 0.5, labels = "(c)", tick = FALSE, las = 1, hadj = 0.5, cex.axis = 1.5) legend("bottomright", legend = expression(x[3]), bg = "white", cex = 2) axis(1, at = seq(0, iT - iBurnPeriod, 50), labels = seq(iBurnPeriod, iT, 50), cex.axis = 1.5) par(mar = c(0, 4, 0.1, 2)) plot(mincpmt[, 4], type = "n", xaxt = "n", xlab = "", ylab = "", las = 1, ylim = vLim, cex.axis = 1.5) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(mincpmt[, 4], col = "black", lwd = 2) axis(4, at = 0.5, labels = "(d)", tick = FALSE, las = 1, hadj = 0.5, cex.axis = 1.5) legend("bottomright", legend = expression(x[4]), bg = "white", cex = 2) par(mar = c(0, 4, 0, 2)) plot(mincpmt[, 5], type = "n", xaxt = "n", xlab = "", ylab = "", las = 1, ylim = vLim, cex.axis = 1.5) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(mincpmt[, 5], col = "black", lwd = 2) axis(4, at = 0.5, labels = "(e)", tick = FALSE, las = 1, hadj = 0.5, cex.axis = 1.5) legend("topright", legend = expression(x[5]), bg = "white", cex = 2) par(mar = c(3, 4, 0, 2)) plot(mincpmt[, 6], type = "n", xaxt = "n", xlab = "", ylab = "", las = 1, ylim = vLim, cex.axis = 1.5) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(mincpmt[, 6], col = "black", lwd = 2) axis(4, at = 0.5, labels = "(f)", tick = FALSE, las = 1, hadj = 0.5, cex.axis = 1.5) legend("topright", legend = expression(x[6]), bg = "white", cex = 2) axis(1, at = seq(0, iT - iBurnPeriod, 50), labels = seq(iBurnPeriod, iT, 50), cex.axis = 1.5) dev.off() ############################################## ## Section 4.2 ############################################## summary(Fit, iBurnPeriod = 50) ############################################## ## Computational challenges : Section 5 ############################################## ## ## Comparison with dma of McCormick et al. (2016). library("dma") library("microbenchmark") set.seed(72463) vK <- seq(4, 16, 2) vT <- c(100, 500, 1000) dV <- 0.1 mW <- diag(max(vK)) * 0.03 dPhi <- 0.98 vBeta0 <- rep(0, max(vK)) mX <- matrix(rnorm(max(vT) * max(vK)), max(vT), max(vK)) lOut <- SimulateDLM(max(vT), mX, vBeta0, mW, dV, dPhi) vY <- lOut$vY ## DMA options vDelta <- 0.95 dAlpha <- 0.99 ## build the powerset for DMA-Raftery RaftyTime <- lModel <- lPowerSet <- list() mTime_gdma <- matrix(NA, length(vT), length(vK), dimnames = list(vT, vK)) ## This is computational expensive for (iT in vT) { RaftyTime[[paste(iT)]] <- list() for (iK in vK) { RaftyTime[[paste(iT)]][[paste(iK)]] <- mean(microbenchmark({ lPowerSet[[paste(iK)]] <- PowerSet(iK)[-1] lModel[[paste(iK)]] <- matrix(0, length(lPowerSet[[paste(iK)]]), iK) for (i in 1:length(lPowerSet[[paste(iK)]])) { lModel[[paste(iK)]][i, lPowerSet[[paste(iK)]][[i]] + 1] <- 1 } dma(mX[1:iT, 1:iK], vY[1:iT], lModel[[paste(iK)]], vDelta, dAlpha, initialperiod = 1) }, times = 10)$time) } } save(RaftyTime, file = paste(sDir_Comp, "RaftyTime.Rdata", sep = "")) ## This is computational expensive for (iK in vK) { for (iT in vT) { mTime_gdma[paste(iT), paste(iK)] <- mean(microbenchmark({ if (iK <= 6) { iCores = 2 } else { iCores = iNcores } DMA(vY[1:iT] ~ mX[1:iT, 1:iK], vDelta = vDelta, dAlpha = dAlpha, iCores = iCores) }, times = 10)$time) } } save(mTime_gdma, file = paste(sDir_Comp, "mTime_gdma.Rdata", sep = "")) ## TABLE 1 load(file = paste(sDir_Comp, "RaftyTime.Rdata", sep = "")) load(file = paste(sDir_Comp, "mTime_gdma.Rdata", sep = "")) ## Comparison Raftery code mTime_Raftery <- matrix(NA, nrow(mTime_gdma), ncol(mTime_gdma), dimnames = dimnames(mTime_gdma)) mTime_Raftery["100", ] <- unlist(RaftyTime[["100"]]) mTime_Raftery["500", ] <- unlist(RaftyTime[["500"]]) mTime_Raftery["1000", ] <- unlist(RaftyTime[["1000"]]) Ratio <- format(round(mTime_Raftery/mTime_gdma, 1), scientific = FALSE) Ratio[, ncol(Ratio)] <- paste(Ratio[, ncol(Ratio)], "\\\\") write.table(Ratio, file = paste(sDir_Tables, "/Table1.txt", sep = ""), quote = FALSE, sep = " & ") rm(list = grep("(^sDir|^iNcores)", ls(), value = TRUE, invert = TRUE)) ## Computational Time eDMA ## This is computational expensive set.seed(7892) vT <- seq(100, 1000, 100) vK <- seq(2, 18, 1) lDelta <- list() for (i in 2:10) { lDelta[[paste(i)]] <- seq(0.9, 1, length.out = i) } cTime <- array(NA, dim = c(length(vT), length(vK), length(lDelta)), dimnames = list(vT, vK, names(lDelta))) dV <- 0.1 mW <- diag(max(vK)) * 0.03 dPhi <- 0.98 vBeta0 <- rep(0, max(vK)) mX <- matrix(rnorm(max(vT) * max(vK)), max(vT), max(vK)) lOut <- SimulateDLM(max(vT), mX, vBeta0, mW, dV, dPhi) vY <- lOut$vY for (iK in vK) { if (iK <= 6) { iCores = 2 } else { iCores = iNcores } for (vDelta in lDelta) { for (iT in vT) { Est <- DMA(vY[1:iT] ~ . - 1, data = mX[1:iT, 1:iK], dAlpha = 0.99, vKeep = NULL, bZellnerPrior = FALSE, dG = 100, bParallelize = TRUE, iCores = iCores, vDelta = vDelta) cTime[paste(iT), paste(iK), paste(length(vDelta))] <- as.double(Est@model$elapsedTime, units = "mins") save(cTime, file = paste(sDir_Comp, "cTime.Rdata", sep = "")) print(paste("iK =", iK, ", iT =", iT, ", iD =", length(vDelta), " - Time=", Sys.time(), "\n")) } } } ## Figure 3 load(paste(sDir_Comp, "cTime.Rdata", sep = "")) vT <- dimnames(cTime)[[1]] vK <- dimnames(cTime)[[2]] vD <- dimnames(cTime)[[3]] vLim <- c(0, max(cTime)) pdf(paste(sDir_Figures, "Figure2.pdf", sep = ""), width = 12, height = 15) iC <- 1 layout(matrix(1:10, 10, 1), heights = c(rep(2, 9), 2.4)) par(mar = c(0, 5, 0.1, 4)) plot(vK, cTime["100", , vD[1]], type = "n", ylim = vLim, xaxt = "n", yaxt = "n", xlab = "", ylab = "") for (iD in vD) { lines(vK, cTime["100", , iD], type = "l") points(vK, cTime["100", , iD], pch = 1) } axis(2, at = c(1, 5, 15, 30), labels = c("1 min", "5 min", "15 min", "30 min"), las = 1, cex.axis = 1.4) lines(0:19, rep(1, 20), type = "l", lwd = 2, lty = 2) lines(0:19, rep(5, 20), type = "l", lwd = 2, lty = 2) lines(0:19, rep(15, 20), type = "l", lwd = 2, lty = 2) lines(0:19, rep(30, 20), type = "l", lwd = 2, lty = 2) legend("topleft", paste("T = 100"), cex = 1.5) axis(4, at = mean(vLim), paste("(", letters[iC], ")", sep = ""), las = 1, tick = FALSE, hadj = 0.3, cex.axis = 2.5) iC <- iC + 1 par(mar = c(0, 5, 0, 4)) for (iT in vT[2:9]) { plot(vK, cTime[iT, , vD[1]], type = "n", ylim = vLim, xaxt = "n", yaxt = "n", xlab = "", ylab = "") for (iD in vD) { lines(vK, cTime[iT, , iD], type = "l") points(vK, cTime[iT, , iD], pch = 1) } axis(2, at = c(1, 5, 15, 30), labels = c("1 min", "5 min", "15 min", "30 min"), las = 1, cex.axis = 1.4) lines(0:19, rep(1, 20), type = "l", lwd = 2, lty = 2) lines(0:19, rep(5, 20), type = "l", lwd = 2, lty = 2) lines(0:19, rep(15, 20), type = "l", lwd = 2, lty = 2) lines(0:19, rep(30, 20), type = "l", lwd = 2, lty = 2) legend("topleft", paste("T =", iT), cex = 1.5) axis(4, at = mean(vLim), paste("(", letters[iC], ")", sep = ""), las = 1, tick = FALSE, hadj = 0.3, cex.axis = 2.5) iC <- iC + 1 } par(mar = c(2, 5, 0, 4)) plot(vK, cTime[vT[10], , vD[1]], type = "n", ylim = vLim, xaxt = "n", yaxt = "n", xlab = "", ylab = "") for (iD in vD) { lines(vK, cTime[vT[10], , iD], type = "l") points(vK, cTime[vT[10], , iD], pch = 1) } axis(2, at = c(1, 5, 15, 30), labels = c("1 min", "5 min", "15 min", "30 min"), las = 1, cex.axis = 1.4) lines(0:19, rep(1, 20), type = "l", lwd = 2, lty = 2) lines(0:19, rep(5, 20), type = "l", lwd = 2, lty = 2) lines(0:19, rep(15, 20), type = "l", lwd = 2, lty = 2) lines(0:19, rep(30, 20), type = "l", lwd = 2, lty = 2) axis(1, at = 2:(length(vK) + 1), labels = vK, padj = 0, tick = -1, cex.axis = 1.5) axis(1, at = 2, labels = "n", tick = FALSE, padj = -0.5, hadj = 5.5, font = 2, cex.axis = 1.5) legend("topleft", paste("T =", vT[10]), cex = 1.5) axis(4, at = mean(vLim), paste("(", letters[iC], ")", sep = ""), las = 1, tick = FALSE, hadj = 0.3, cex.axis = 2.5) dev.off() rm(list = grep("(^sDir)", ls(), value = TRUE, invert = TRUE)) ############################################# ## Empirical Application : Section 6 ############################################# library("eDMA") data("USData", package = "eDMA") head(round(USData[, 1:6], 2)) Lag(USData[, "GDPDEF"], 1) ############################################## ## Section 6.2 ############################################## Fit <- DMA(GDPDEF ~ Lag(GDPDEF, 1) + Lag(GDPDEF, 2) + Lag(GDPDEF, 3) + Lag(GDPDEF, 4) + Lag(ROUTP, 1) + Lag(RCONS, 1) + Lag(RINVR, 1) + Lag(PIMP, 1) + Lag(UNEMP, 1) + Lag(NFPR, 1) + Lag(HSTS, 1) + Lag(M2, 1) + Lag(OIL, 1) + Lag(RAW, 1) + Lag(FOOD, 1) + Lag(YL, 1) + Lag(TS, 1) + Lag(CS, 1) + Lag(MS, 1), data = USData, vDelta = seq(0.90, 1.00, 0.01), vKeep = 1, dBeta = 0.96, dAlpha = 0.99) Fit summary(Fit, iBurnPeriod = 32) ############################################## ## Section 6.3 ############################################## InclusionProb <- inclusion.prob(Fit, iBurnPeriod = 32) tail(round(InclusionProb[, 1:4], 2)) mTheta <- coef(Fit, iBurnPeriod = 32) plot(Fit, which = "vdeltahat", iBurnPeriod = 32) InclusionProbDelta <- as.data.frame(Fit, which = "mpmt", iBurnPeriod = 32) round(tail(InclusionProbDelta), 2) plot(Fit, which = "vsize_DMS", iBurnPeriod = 32) ############################################## ## Section 6.4 ############################################## Fit_M0 <- DMA(GDPDEF ~ Lag(GDPDEF, 1) + Lag(GDPDEF, 2) + Lag(GDPDEF, 3) + Lag(GDPDEF, 4), data = USData, vDelta = 1.00, dAlpha = 1.00, vKeep = c(1, 2, 3, 4, 5), dBeta = 1.0) ############################################## ## Section 6.5 ############################################## vPL_M0 <- pred.like(Fit_M0, iBurnPeriod = 32) vPL_M3 <- pred.like(Fit, iBurnPeriod = 32) vPLD_M3.M0 <- cumsum(vPL_M3 - vPL_M0) ############################################## ## Section 6: Additional results ############################################## ## h = 1 AR4_formula <- GDPDEF ~ Lag(GDPDEF, 1) + Lag(GDPDEF, 2) + Lag(GDPDEF, 3) + Lag(GDPDEF, 4) AR4_Reg_formula <- GDPDEF ~ Lag(GDPDEF, 1) + Lag(GDPDEF, 2) + Lag(GDPDEF, 3) + Lag(GDPDEF, 4) + Lag(ROUTP, 1) + Lag(RCONS, 1) + Lag(RINVR, 1) + Lag(PIMP, 1) + Lag(UNEMP, 1) + Lag(NFPR, 1) + Lag(HSTS, 1) + Lag(M2, 1) + Lag(OIL, 1) + Lag(RAW, 1) + Lag(FOOD, 1) + Lag(YL, 1) + Lag(TS, 1) + Lag(CS, 1) + Lag(MS, 1) vDelta <- seq(0.9, 1, 0.01) dBeta <- 0.96 Model0 <- DMA(AR4_formula, data = USData, vDelta = 1, dAlpha = 1, vKeep = c(1, 2, 3, 4, 5), bZellnerPrior = FALSE, dG = 100, dBeta = 1.0) save(Model0, file = paste(sDir_h1, "Model0.RData", sep = "")) Model1 <- DMA(AR4_formula, data = USData, vDelta = vDelta, dAlpha = 0.99, vKeep = c(1, 2, 3, 4, 5), bZellnerPrior = FALSE, dG = 100, dBeta = dBeta) save(Model1, file = paste(sDir_h1, "Model1.RData", sep = "")) Model2 <- DMA(AR4_formula, data = USData, vDelta = vDelta, dAlpha = 0.99, vKeep = c(1), bZellnerPrior = FALSE, dG = 100, dBeta = dBeta) save(Model2, file = paste(sDir_h1, "Model2.RData", sep = "")) Model3_4 <- DMA(AR4_Reg_formula, data = USData, vDelta = vDelta, dAlpha = 0.99, vKeep = c(1), bZellnerPrior = FALSE, dG = 100, dBeta = dBeta) save(Model3_4, file = paste(sDir_h1, "Model3_4.RData", sep = "")) Model5_6 <- DMA(AR4_Reg_formula, data = USData, vDelta = 1, dAlpha = 1, vKeep = c(1), bZellnerPrior = FALSE, dG = 100, dBeta = 1.0) save(Model5_6, file = paste(sDir_h1, "Model5_6.RData", sep = "")) Model7 <- DMA(AR4_Reg_formula, data = USData, vDelta = vDelta, dAlpha = 0.99, vKeep = "KS", bZellnerPrior = FALSE, dG = 100, dBeta = dBeta) save(Model7, file = paste(sDir_h1, "Model7.RData", sep = "")) ## PRIOR SENTITIVITY ANALYS 6.6 Model_3_Z01 <- DMA(AR4_Reg_formula, data = USData, vDelta = vDelta, dAlpha = 0.99, vKeep = c(1), bZellnerPrior = FALSE, dG = 0.1, dBeta = dBeta) save(Model_3_Z01, file = paste(sDir_h1, "Model_3_Z01.RData", sep = "")) Model_3_Z1 <- DMA(AR4_Reg_formula, data = USData, vDelta = vDelta, dAlpha = 0.99, vKeep = c(1), bZellnerPrior = FALSE, dG = 1, dBeta = dBeta) save(Model_3_Z1, file = paste(sDir_h1, "Model_3_Z1.RData", sep = "")) Model_3_Z20 <- DMA(AR4_Reg_formula, data = USData, vDelta = vDelta, dAlpha = 0.99, vKeep = c(1), bZellnerPrior = FALSE, dG = 20, dBeta = dBeta) save(Model_3_Z20, file = paste(sDir_h1, "Model_3_Z20.RData", sep = "")) Model_3_ZTo2 <- DMA(AR4_Reg_formula, data = USData, vDelta = vDelta, dAlpha = 0.99, vKeep = c(1), bZellnerPrior = FALSE, dG = 206/2, dBeta = dBeta) save(Model_3_ZTo2, file = paste(sDir_h1, "Model_3_ZTo2.RData", sep = "")) Model_3_Z206 <- DMA(AR4_Reg_formula, data = USData, vDelta = vDelta, dAlpha = 0.99, vKeep = c(1), bZellnerPrior = FALSE, dG = 206, dBeta = dBeta) save(Model_3_Z206, file = paste(sDir_h1, "Model_3_Z206.RData", sep = "")) ## h = 5 AR4_formula <- GDPDEF ~ Lag(GDPDEF, 5) + Lag(GDPDEF, 6) + Lag(GDPDEF, 7) + Lag(GDPDEF, 8) AR4_Reg_formula <- GDPDEF ~ Lag(GDPDEF, 5) + Lag(GDPDEF, 6) + Lag(GDPDEF, 7) + Lag(GDPDEF, 8) + Lag(ROUTP, 5) + Lag(RCONS, 5) + Lag(RINVR, 5) + Lag(PIMP, 5) + Lag(UNEMP, 5) + Lag(NFPR, 1) + Lag(HSTS, 5) + Lag(M2, 5) + Lag(OIL, 5) + Lag(RAW, 5) + Lag(FOOD, 5) + Lag(YL, 5) + Lag(TS, 5) + Lag(CS, 5) + Lag(MS, 5) vDelta <- seq(0.9, 1, 0.01) Model0 <- DMA(AR4_formula, data = USData, vDelta = 1, dAlpha = 1, vKeep = c(1, 2, 3, 4, 5), bZellnerPrior = FALSE, dG = 100, dBeta = 1.0) save(Model0, file = paste(sDir_h5, "Model0.RData", sep = "")) Model1 <- DMA(AR4_formula, data = USData, vDelta = vDelta, dAlpha = 0.99, vKeep = c(1, 2, 3, 4, 5), bZellnerPrior = FALSE, dG = 100, dBeta = dBeta) save(Model1, file = paste(sDir_h5, "Model1.RData", sep = "")) Model2 <- DMA(AR4_formula, data = USData, vDelta = vDelta, dAlpha = 0.99, vKeep = c(1), bZellnerPrior = FALSE, dG = 100, dBeta = dBeta) save(Model2, file = paste(sDir_h5, "Model2.RData", sep = "")) Model3_4 <- DMA(AR4_Reg_formula, data = USData, vDelta = vDelta, dAlpha = 0.99, vKeep = c(1), bZellnerPrior = FALSE, dG = 100, dBeta = dBeta) save(Model3_4, file = paste(sDir_h5, "Model3_4.RData", sep = "")) Model5_6 <- DMA(AR4_Reg_formula, data = USData, vDelta = 1, dAlpha = 1, vKeep = c(1), bZellnerPrior = FALSE, dG = 100, dBeta = 1.0) save(Model5_6, file = paste(sDir_h5, "Model5_6.RData", sep = "")) Model7 <- DMA(AR4_Reg_formula, data = USData, vDelta = vDelta, dAlpha = 0.99, vKeep = "KS", bZellnerPrior = FALSE, dG = 100, dBeta = dBeta) save(Model7, file = paste(sDir_h5, "Model7.RData", sep = "")) ## RMSE and PredLike tables iBurnPeriod <- 32 vFilesh1 <- list.files(sDir_h1, full.names = TRUE) vFilesh5 <- list.files(sDir_h1, full.names = TRUE) vModels <- paste("Model", 0:7, sep = "") mTab_h1 <- matrix(0, length(vModels), 2, dimnames = list(vModels, c("RMSE", "LS"))) mTab_h5 <- matrix(0, length(vModels), 2, dimnames = list(vModels, c("RMSE", "LS"))) for (i in 0:2) { load(paste(sDir_h1, "Model", i, ".RData", sep = "")) object <- get(paste("Model", i, sep = "")) Back <- BacktestDMA(object, iBurnPeriod) mTab_h1[i + 1, "RMSE"] <- Back["MSE", "DMA"] mTab_h1[i + 1, "LS"] <- Back["Log-predictive Likelihood", "DMA"] load(paste(sDir_h5, "Model", i, ".RData", sep = "")) object <- get(paste("Model", i, sep = "")) Back <- BacktestDMA(object, iBurnPeriod) mTab_h5[i + 1, "RMSE"] <- Back["MSE", "DMA"] mTab_h5[i + 1, "LS"] <- Back["Log-predictive Likelihood", "DMA"] } ## Model 3 h1 load(paste(sDir_h1, "Model3_4.RData", sep = "")) object <- get("Model3_4") Back <- BacktestDMA(object, iBurnPeriod) mTab_h1[4, "RMSE"] <- Back["MSE", "DMA"] mTab_h1[4, "LS"] <- Back["Log-predictive Likelihood", "DMA"] ## Model 4 mTab_h1[5, "RMSE"] <- Back["MSE", "DMS"] mTab_h1[5, "LS"] <- Back["Log-predictive Likelihood", "DMS"] ## h5 load(paste(sDir_h5, "Model3_4.RData", sep = "")) object <- get("Model3_4") Back <- BacktestDMA(object, iBurnPeriod) mTab_h5[4, "RMSE"] <- Back["MSE", "DMA"] mTab_h5[4, "LS"] <- Back["Log-predictive Likelihood", "DMA"] ## Model 4 mTab_h5[5, "RMSE"] <- Back["MSE", "DMS"] mTab_h5[5, "LS"] <- Back["Log-predictive Likelihood", "DMS"] ## Model 5 h1 load(paste(sDir_h1, "Model5_6.RData", sep = "")) object <- get("Model5_6") Back <- BacktestDMA(object, iBurnPeriod) mTab_h1[6, "RMSE"] <- Back["MSE", "DMA"] mTab_h1[6, "LS"] <- Back["Log-predictive Likelihood", "DMA"] ## Model 6 mTab_h1[7, "RMSE"] <- Back["MSE", "DMS"] mTab_h1[7, "LS"] <- Back["Log-predictive Likelihood", "DMS"] ## Model 5 h5 load(paste(sDir_h5, "Model5_6.RData", sep = "")) object <- get("Model5_6") Back <- BacktestDMA(object, iBurnPeriod) mTab_h5[6, "RMSE"] <- Back["MSE", "DMA"] mTab_h5[6, "LS"] <- Back["Log-predictive Likelihood", "DMA"] ## Model 6 mTab_h5[7, "RMSE"] <- Back["MSE", "DMS"] mTab_h5[7, "LS"] <- Back["Log-predictive Likelihood", "DMS"] ## ## Model 7 h1 load(paste(sDir_h1, "Model7.RData", sep = "")) object <- get("Model7") Back <- BacktestDMA(object, iBurnPeriod) mTab_h1[8, "RMSE"] <- Back["MSE", "DMA"] mTab_h1[8, "LS"] <- Back["Log-predictive Likelihood", "DMA"] ## Model 7 h5 load(paste(sDir_h5, "Model7.RData", sep = "")) object <- get("Model7") Back <- BacktestDMA(object, iBurnPeriod) mTab_h5[8, "RMSE"] <- Back["MSE", "DMS"] mTab_h5[8, "LS"] <- Back["Log-predictive Likelihood", "DMS"] ## for (i in 2:length(vModels)) { mTab_h1[i, "RMSE"] <- mTab_h1[i, "RMSE"]/mTab_h1[1, "RMSE"] mTab_h1[i, "LS"] <- mTab_h1[i, "LS"] - mTab_h1[1, "LS"] mTab_h5[i, "RMSE"] <- mTab_h5[i, "RMSE"]/mTab_h5[1, "RMSE"] mTab_h5[i, "LS"] <- mTab_h5[i, "LS"] - mTab_h5[1, "LS"] } mTab_h1[1, "RMSE"] <- 1 mTab_h1[1, "LS"] <- 0 mTab_h5[1, "RMSE"] <- 1 mTab_h5[1, "LS"] <- 0 mTab <- cbind(mTab_h1, mTab_h5) mTab <- format(round(mTab, 3), scientific = FALSE) mTab[, ncol(mTab)] <- paste(mTab[, ncol(mTab)], "\\\\") mTab <- mTab[-1, ] rownames(mTab) <- paste("$\\mathcal{M}_{", 1:7, "}$", sep = "") write.table(mTab, paste(sDir_Tables, "Table3.txt", sep = ""), sep = " & ", quote = FALSE, col.names = TRUE, row.names = TRUE) ## PRIOR SENTITIVITY load(file = paste(sDir_h1, "Model0.RData", sep = "")) load(file = paste(sDir_h1, "Model_3_Z01.RData", sep = "")) load(file = paste(sDir_h1, "Model_3_Z20.RData", sep = "")) load(file = paste(sDir_h1, "Model_3_ZTo2.RData", sep = "")) load(file = paste(sDir_h1, "Model_3_Z206.RData", sep = "")) M0 <- BacktestDMA(Model0, iBurnPeriod)[c("MSE", "Log-predictive Likelihood"), "DMA"] mTab <- matrix(NA, 4, 2, dimnames = list(c("Z01", "Z20", "ZTo2", "Z206"), c("RMSE", "PLD"))) mTab["Z01", ] <- BacktestDMA(Model_3_Z1, iBurnPeriod)[c("MSE", "Log-predictive Likelihood"), "DMA"] mTab["Z20", ] <- BacktestDMA(Model_3_Z20, iBurnPeriod)[c("MSE", "Log-predictive Likelihood"), "DMA"] mTab["ZTo2", ] <- BacktestDMA(Model_3_ZTo2, iBurnPeriod)[c("MSE", "Log-predictive Likelihood"), "DMA"] mTab["Z206", ] <- BacktestDMA(Model_3_Z206, iBurnPeriod)[c("MSE", "Log-predictive Likelihood"), "DMA"] mTab[, "RMSE"] <- mTab[, "RMSE"]/M0["MSE"] mTab[, "PLD"] <- mTab[, "PLD"] - M0["Log-predictive Likelihood"] mTab <- format(round(mTab, 3), scientific = FALSE) rownames(mTab) <- paste("$g = ", c(0.1, 20, "T/2", "T"), "$", sep = "") mTab[, ncol(mTab)] <- paste(mTab[, ncol(mTab)], "\\\\") write.table(mTab, paste(sDir_Tables, "Table4.txt", sep = ""), sep = " & ", quote = FALSE, col.names = TRUE, row.names = TRUE) # FIGURES iBurnPeriod <- 32 data("USRecessions", package = "eDMA") vBurnIn <- 1:iBurnPeriod rgb <- c("#D33F6A4D", "#E2E2E24D", "#4A6FE34D", "#6969694D") ## Model 3 load(paste(sDir_h1, "Model3_4.RData", sep = "")) Est <- Model3_4 ## Figure 4 v2Plot <- c(2, 4, 5, 9, 20, 13, 14, 17) mincpmt <- as.data.frame(Est, which = "mincpmt", iBurnPeriod = iBurnPeriod) vNames <- c("const", "lag1", "lag2", "lag3", "lag4", colnames(USData)[-1]) colnames(mincpmt) <- vNames vDates <- as.Date(index(mincpmt)) mincpmt <- mincpmt[, v2Plot] plotSeq <- seq(1, 9, 4) vLim <- c(0, 1) pdf(paste(sDir_Figures, "Figure3.pdf", sep = ""), width = 15, height = 10) layout(matrix(1:8, 4, 2), heights = c(rep(2, 3), 2.5, rep(2, 3), 2.5)) for (i in 1:length(v2Plot)) { if (any(i == plotSeq)) par(mar = c(0, 4, 0.1, 2)) if (all(i != plotSeq) & all(i != plotSeq - 1)) par(mar = c(0, 4, 0, 2)) if (any(i == plotSeq - 1)) par(mar = c(3, 4, 0, 2)) plot(vDates, mincpmt[, i], type = "n", xaxt = "n", xlab = "", ylab = "", las = 1, ylim = vLim, cex.axis = 1.5) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, mincpmt[, i], col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) if (any(i == c(1, 2, 3))) { if (i == 1) legend("bottomright", legend = expression(y[t - 1]), bg = "white") if (i == 2) legend("topright", legend = expression(y[t - 3]), bg = "white") if (i == 3) legend("topright", legend = expression(y[t - 4]), bg = "white") } else { if (i == 4) { legend("bottomright", legend = colnames(mincpmt)[i], bg = "white") } else { legend("topright", legend = colnames(mincpmt)[i], bg = "white") } } axis(4, at = mean(vLim), labels = paste("(", letters[i], ")", sep = ""), tick = FALSE, hadj = 0.5, las = 1, cex.axis = 1.5) if (any(i == plotSeq - 1)) { axis.Date(1, at = seq(min(vDates), max(vDates), "year"), cex.axis = 1.5) axis.Date(1, at = seq(min(vDates), max(vDates), "quarter"), labels = FALSE, tcl = -0.2) } } dev.off() ## Figure 5 v2Plot <- c(2, 4, 5, 9, 20, 13, 14, 17) mmhat <- coef(Est, iBurnPeriod = iBurnPeriod) colnames(mmhat) <- vNames mmhat <- mmhat[, v2Plot] vDates <- as.Date(index(mmhat)) pdf(paste(sDir_Figures, "/Figure4.pdf", sep = ""), width = 15, height = 10) layout(matrix(1:8, 4, 2), heights = c(rep(2, 3), 2.2, rep(2, 3), 2.2)) par(mar = c(0, 4, 0.1, 2)) vLim <- c(0, 0.75) plot(vDates, mmhat[, 1], type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", las = 1, ylim = vLim, cex.axis = 1.5) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, mmhat[, 1], col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) legend("topright", legend = expression(y[t - 1]), bg = "white") axis(2, at = seq(0, 0.75, 0.15), labels = c("0.00", "0.15", "0.30", "0.45", "0.60", "0.75"), las = 1, cex.axis = 1.5) axis(4, at = mean(vLim), labels = paste("(", letters[1], ")", sep = ""), tick = FALSE, hadj = 0.5, las = 1, cex.axis = 1.5) par(mar = c(0, 4, 0, 2)) vLim <- c(0, 0.35) plot(vDates, mmhat[, 2], type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", las = 1, ylim = vLim, cex.axis = 1.5) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, mmhat[, 2], col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) legend("topright", legend = expression(y[t - 3]), bg = "white") axis(2, at = seq(0, 0.35, 0.07), labels = c("0.00", "0.07", "0.14", "0.21", "0.28", "0.35"), las = 1, cex.axis = 1.5) axis(4, at = mean(vLim), labels = paste("(", letters[2], ")", sep = ""), tick = FALSE, hadj = 0.5, las = 1, cex.axis = 1.5) par(mar = c(0, 4, 0, 2)) vLim <- c(0, 0.4) plot(vDates, mmhat[, 3], type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", las = 1, ylim = vLim, cex.axis = 1.5) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, mmhat[, 3], col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) legend("topright", legend = expression(y[t - 4]), bg = "white") axis(2, at = seq(0, 0.4, 0.08), labels = c("0.00", "0.08", "0.16", "0.24", "0.32", "0.40"), las = 1, cex.axis = 1.5) axis(4, at = mean(vLim), labels = paste("(", letters[3], ")", sep = ""), tick = FALSE, hadj = 0.5, las = 1, cex.axis = 1.5) par(mar = c(3, 4, 0, 2)) vLim <- c(-0.1, 0.35) plot(vDates, mmhat[, 4], type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", las = 1, ylim = vLim, cex.axis = 1.5) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, mmhat[, 4], col = "black", lwd = 2) legend("topright", legend = colnames(mmhat)[4], bg = "white") xblocks(USRecessions == 1, col = rgb[4]) axis(2, at = seq(-0.1, 0.35, 0.09), labels = c("-0.10", "-0.01", "0.08", "0.17", "0.26", "0.35"), las = 1, cex.axis = 1.5) axis(4, at = mean(vLim), labels = paste("(", letters[4], ")", sep = ""), tick = FALSE, hadj = 0.5, las = 1, cex.axis = 1.5) axis.Date(1, at = seq(min(vDates), max(vDates), "year"), cex.axis = 1.5) axis.Date(1, at = seq(min(vDates), max(vDates), "quarter"), labels = FALSE, tcl = -0.2) par(mar = c(0, 4, 0.1, 2)) vLim <- c(-0.03, 0.4) plot(vDates, mmhat[, 5], type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", las = 1, ylim = vLim, cex.axis = 1.5) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, mmhat[, 5], col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) legend("topright", legend = colnames(mmhat)[5], bg = "white") axis(2, at = seq(-0.03, 0.4, 0.08), labels = c("-0.03", "0.05", "0.13", "0.21", "0.29", "0.37"), las = 1, cex.axis = 1.5) axis(4, at = mean(vLim), labels = paste("(", letters[5], ")", sep = ""), tick = FALSE, hadj = 0.5, las = 1, cex.axis = 1.5) par(mar = c(0, 4, 0, 2)) vLim <- c(-0.02, 0.08) plot(vDates, mmhat[, 6], type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", las = 1, ylim = vLim, cex.axis = 1.5) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, mmhat[, 6], col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) legend("topright", legend = colnames(mmhat)[6], bg = "white") axis(2, at = seq(-0.02, 0.08, 0.02), labels = c("-0.02", "0.00", "0.02", "0.04", "0.06", "0.08"), las = 1, cex.axis = 1.5) axis(4, at = mean(vLim), labels = paste("(", letters[6], ")", sep = ""), tick = FALSE, hadj = 0.5, las = 1, cex.axis = 1.5) par(mar = c(0, 4, 0, 2)) vLim <- c(-0.2, 0.4) plot(vDates, mmhat[, 7], type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", las = 1, ylim = vLim, cex.axis = 1.5) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, mmhat[, 7], col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) legend("topright", legend = colnames(mmhat)[7], bg = "white") axis(2, at = seq(-0.2, 0.4, 0.12), labels = c("-0.20", "-0.08", "0.04", "0.16", "0.28", "0.40"), las = 1, cex.axis = 1.5) axis(4, at = mean(vLim), labels = paste("(", letters[7], ")", sep = ""), tick = FALSE, hadj = 0.5, las = 1, cex.axis = 1.5) par(mar = c(3, 4, 0, 2)) vLim <- c(0, 1.5) plot(vDates, mmhat[, 8], type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", las = 1, ylim = vLim, cex.axis = 1.5) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, mmhat[, 8], col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) legend("topright", legend = colnames(mmhat)[8], bg = "white") axis(2, at = seq(0, 1.5, 0.3), labels = c("0.00", "0.30", "0.60", "0.90", "1.20", "1.50"), las = 1, cex.axis = 1.5) axis(4, at = mean(vLim), labels = paste("(", letters[8], ")", sep = ""), tick = FALSE, hadj = 0.5, las = 1, cex.axis = 1.5) axis.Date(1, at = seq(min(vDates), max(vDates), "year"), cex.axis = 1.5) axis.Date(1, at = seq(min(vDates), max(vDates), "quarter"), labels = FALSE, tcl = -0.2) dev.off() ## Figure 5 ## This is a 3x2 plot. (a): Posterior weighted average (b): vsize DMA, ## vsize DMS estimate of delta. In panel (c), we report ceil ## (0.1*prob) vhighmpTop01_DMS. In panel (d), we plot vobs, in panel ## (e), we plot vcoeff and in panel (f), we plot mod and TVP from ## panel (d) of Figure 4 in our paper ## a) vDeltaHat vDeltaHat <- as.data.frame(Est, which = "vdeltahat", iBurnPeriod = iBurnPeriod) ## b) vsize DMS vsizeDMS <- as.data.frame(Est, which = "vsize_DMS", iBurnPeriod = iBurnPeriod) ## c) vhighmpTop01 vhighmpTop01 <- as.data.frame(Est, which = "vhighmpTop01_DMS", iBurnPeriod = iBurnPeriod) ## VarianceDecomposition : vobs/vtotal, vcoeff/vtotal, vmod/vtotal, vtvp/vtotal iT <- length(vDates) VarianceDecomposition <- matrix(NA, iT, 4) colnames(VarianceDecomposition) <- c("vobs", "vcoef", "vmod", "vtvp") VarianceDecomposition[, 1] <- (Est@Est$vobs/Est@Est$vtotal)[-vBurnIn] VarianceDecomposition[, 2] <- (Est@Est$vcoeff/Est@Est$vtotal)[-vBurnIn] VarianceDecomposition[, 3] <- (Est@Est$vmod/Est@Est$vtotal)[-vBurnIn] VarianceDecomposition[, 4] <- (Est@Est$vtvp/Est@Est$vtotal)[-vBurnIn] ## d) vobs vobs <- VarianceDecomposition[, "vobs"] ## e) vcoef vcoef <- VarianceDecomposition[, "vcoef"] ## f) vmod and vtvp pdf(paste(sDir_Figures, "Figure5.pdf", sep = ""), width = 15, height = 10) layout(matrix(1:6, 3, 2), heights = c(rep(2, 2), 2.2, rep(2, 2), 2.2)) ## a) vDeltaHat par(mar = c(0, 4, 0.1, 2)) vLim <- c(0.9, 1) plot(vDates, vDeltaHat, type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ylim = vLim) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, vDeltaHat, col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) axis(4, at = mean(vLim), paste("(a)"), las = 1, tick = FALSE, hadj = 0.5, cex.axis = 1.5) axis(2, at = seq(0.9, 1, 0.02), labels = c("0.90", "0.92", "0.94", "0.96", "0.98", "1.00"), las = 1, cex.axis = 1.5) ## b) vsize DMS par(mar = c(0, 4, 0, 2)) vLim <- c(0, 10) plot(vDates, vsizeDMS, type = "n", xaxt = "n", xlab = "", yaxt = "n", ylab = "", ylim = vLim) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, vsizeDMS, col = "black", lwd = 2, lty = 1) xblocks(USRecessions == 1, col = rgb[4]) axis(4, at = mean(vLim), paste("(b)"), las = 1, tick = FALSE, hadj = 0.5, cex.axis = 1.5) axis(2, at = seq(0, 10, 2), labels = c("0.00", "2.00", "4.00", "6.00", "8.00", "10.00"), las = 1, cex.axis = 1.5) ## c) vhighmpTop01 par(mar = c(2, 4, 0, 2)) vLim <- c(0, 0.21) plot(vDates, vhighmpTop01, type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ylim = vLim) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, vhighmpTop01, col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) axis(4, at = mean(vLim), paste("(c)"), las = 1, tick = FALSE, hadj = 0.5, cex.axis = 1.5) axis(2, at = seq(0, 0.20, length.out = 6), labels = c("0.00", "0.04", "0.08", "0.12", "0.16", "0.20"), las = 1, cex.axis = 1.5) axis.Date(1, at = seq(min(vDates), max(vDates), "year"), cex.axis = 1.5) axis.Date(1, at = seq(min(vDates), max(vDates), "quarter"), labels = FALSE, tcl = -0.2) ## d) vobs par(mar = c(0, 4, 0.1, 2)) vLim <- c(0, 1) plot(vDates, vobs, type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ylim = vLim) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, vobs, col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) axis(4, at = mean(vLim), paste("(d)"), las = 1, tick = FALSE, hadj = 0.5, cex.axis = 1.5) axis(2, at = seq(0, 1, 0.2), labels = c("0.00", "0.20", "0.40", "0.60", "0.80", "1.00"), las = 1, cex.axis = 1.5) ## e) vcoef par(mar = c(0, 4, 0, 2)) vLim <- c(0, 1) plot(vDates, vcoef, type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ylim = vLim) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, vcoef, col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) axis(4, at = mean(vLim), paste("(e)"), las = 1, tick = FALSE, hadj = 0.5, cex.axis = 1.5) axis(2, at = seq(0, 1, 0.2), labels = c("0.00", "0.20", "0.40", "0.60", "0.80", "1.00"), las = 1, cex.axis = 1.5) ## f) vmod and vtvp par(mar = c(2, 4, 0, 2)) vmod_vtvp = VarianceDecomposition[, c("vmod", "vtvp")] vLim <- c(0, 1) plot(vDates, vmod_vtvp[, 1], type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ylim = vLim) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, vmod_vtvp[, "vmod"], col = 1, lwd = 2, lty = 1) lines(vDates, vmod_vtvp[, "vtvp"], col = 2, lwd = 2, lty = 2) xblocks(USRecessions == 1, col = rgb[4]) axis(4, at = mean(vLim), paste("(f)"), las = 1, tick = FALSE, hadj = 0.5, cex.axis = 1.5) axis(2, at = seq(0, 1, 0.2), labels = c("0.00", "0.20", "0.40", "0.60", "0.80", "1.00"), las = 1, cex.axis = 1.5) legend("topleft", legend = c("Mod", "TVP"), col = 1:2, lty = 1:2, lwd = 2, horiz = FALSE, bg = "white", cex = 2) axis.Date(1, at = seq(min(vDates), max(vDates), "year"), cex.axis = 1.5) axis.Date(1, at = seq(min(vDates), max(vDates), "quarter"), labels = FALSE, tcl = -0.2) dev.off() ## Figure 7 load(file = paste(sDir_h1, "Model_3_Z01.RData", sep = "")) load(file = paste(sDir_h1, "Model_3_Z20.RData", sep = "")) load(file = paste(sDir_h1, "Model_3_ZTo2.RData", sep = "")) load(file = paste(sDir_h1, "Model_3_Z206.RData", sep = "")) load(file = paste(sDir_h1, "Model0.RData", sep = "")) load(file = paste(sDir_h1, "Model1.RData", sep = "")) load(file = paste(sDir_h1, "Model3_4.RData", sep = "")) load(file = paste(sDir_h1, "Model5_6.RData", sep = "")) load(file = paste(sDir_h1, "Model7.RData", sep = "")) M0 <- pred.like(Model0, iBurnPeriod = iBurnPeriod) mLPDF <- matrix(0, length(M0), 4, dimnames = list(index(M0), c("M1", "M3", "M5", "M7"))) mLPDF[, "M1"] <- pred.like(Model1, iBurnPeriod = iBurnPeriod) mLPDF[, "M3"] <- pred.like(Model3_4, iBurnPeriod = iBurnPeriod) mLPDF[, "M5"] <- pred.like(Model5_6, iBurnPeriod = iBurnPeriod) mLPDF[, "M7"] <- pred.like(Model7, iBurnPeriod = iBurnPeriod) mLPDF <- apply(mLPDF, 2, function(x, M0) cumsum(x - M0), M0 = M0) mESize <- matrix(0, length(M0), 4, dimnames = list(index(M0), c("Z01", "Z20", "ZTo2", "Z206"))) mESize[, "Z01"] <- as.data.frame(Model_3_Z01, which = "vsize_DMS", iBurnPeriod = iBurnPeriod) mESize[, "Z20"] <- as.data.frame(Model_3_Z20, which = "vsize_DMS", iBurnPeriod = iBurnPeriod) mESize[, "ZTo2"] <- as.data.frame(Model_3_ZTo2, which = "vsize_DMS", iBurnPeriod = iBurnPeriod) mESize[, "Z206"] <- as.data.frame(Model_3_Z206, which = "vsize_DMS", iBurnPeriod = iBurnPeriod) vDates <- as.Date(index(M0)) pdf(paste(sDir_Figures, "Figure6.pdf", sep = ""), width = 15, height = 10) layout(matrix(1:8, 4, 2), heights = c(rep(2, 3), 2.2, rep(2, 3), 2.2)) ## LPDF vLim <- c(-18, 26) par(mar = c(0, 4.5, 0.1, 2)) plot(vDates, mLPDF[, 1], type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ylim = vLim) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, mLPDF[, 1], col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) abline(h = 0, lwd = 2, lty = 2, col = "orange") axis(4, at = mean(vLim), paste("(", letters[1], ")", sep = ""), las = 1, tick = FALSE, hadj = 0.5, cex.axis = 1.5) axis(2, at = seq(-18, 24, length.out = 6), labels = c("-18.00", "-9.60", "-1.20", "7.20", "15.60", "24.00"), las = 1, cex.axis = 1.5) par(mar = c(0, 4.5, 0, 2)) plot(vDates, mLPDF[, 2], type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ylim = vLim) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, mLPDF[, 2], col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) abline(h = 0, lwd = 2, lty = 2, col = "orange") axis(4, at = mean(vLim), paste("(", letters[2], ")", sep = ""), las = 1, tick = FALSE, hadj = 0.5, cex.axis = 1.5) axis(2, at = seq(-18, 24, length.out = 6), labels = c("-18.00", "-9.60", "-1.20", "7.20", "15.60", "24.00"), las = 1, cex.axis = 1.5) par(mar = c(0, 4.5, 0, 2)) plot(vDates, mLPDF[, 3], type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ylim = vLim) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, mLPDF[, 3], col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) abline(h = 0, lwd = 2, lty = 2, col = "orange") axis(4, at = mean(vLim), paste("(", letters[3], ")", sep = ""), las = 1, tick = FALSE, hadj = 0.5, cex.axis = 1.5) axis(2, at = seq(-18, 24, length.out = 6), labels = c("-18.00", "-9.60", "-1.20", "7.20", "15.60", "24.00"), las = 1, cex.axis = 1.5) par(mar = c(2, 4.5, 0, 2)) plot(vDates, mLPDF[, 4], type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ylim = vLim) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, mLPDF[, 4], col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) abline(h = 0, lwd = 2, lty = 2, col = "orange") axis(4, at = mean(vLim), paste("(", letters[4], ")", sep = ""), las = 1, tick = FALSE, hadj = 0.5, cex.axis = 1.5) axis(2, at = seq(-18, 24, length.out = 6), labels = c("-18.00", "-9.60", "-1.20", "7.20", "15.60", "24.00"), las = 1, cex.axis = 1.5) axis.Date(1, at = seq(min(vDates), max(vDates), "year"), cex.axis = 1.5) axis.Date(1, at = seq(min(vDates), max(vDates), "quarter"), labels = FALSE, tcl = -0.2) ## ESIZEDMS par(mar = c(0, 4.5, 0.1, 2)) vLim <- c(0, 11) plot(vDates, mESize[, 1], type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ylim = vLim) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, mESize[, 1], col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) axis(4, at = mean(vLim), paste("(", letters[5], ")", sep = ""), las = 1, tick = FALSE, hadj = 0.5, cex.axis = 1.5) axis(2, at = seq(0, 10, length.out = 6), labels = c("0.00", "2.00", "4.00", "6.00", "8.00", "10.00"), las = 1, cex.axis = 1.5) par(mar = c(0, 4.5, 0, 2)) plot(vDates, mESize[, 2], type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ylim = vLim) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, mESize[, 2], col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) axis(4, at = mean(vLim), paste("(", letters[6], ")", sep = ""), las = 1, tick = FALSE, hadj = 0.5, cex.axis = 1.5) axis(2, at = seq(0, 10, length.out = 6), labels = c("0.00", "2.00", "4.00", "6.00", "8.00", "10.00"), las = 1, cex.axis = 1.5) par(mar = c(0, 4.5, 0, 2)) plot(vDates, mESize[, 3], type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ylim = vLim) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, mESize[, 3], col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) axis(4, at = mean(vLim), paste("(", letters[7], ")", sep = ""), las = 1, tick = FALSE, hadj = 0.5, cex.axis = 1.5) axis(2, at = seq(0, 10, length.out = 6), labels = c("0.00", "2.00", "4.00", "6.00", "8.00", "10.00"), las = 1, cex.axis = 1.5) par(mar = c(2, 4.5, 0, 2)) plot(vDates, mESize[, 4], type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ylim = vLim) grid(nx = 10, ny = 10, col = "gray", lty = "dotted") lines(vDates, mESize[, 4], col = "black", lwd = 2) xblocks(USRecessions == 1, col = rgb[4]) axis(4, at = mean(vLim), paste("(", letters[8], ")", sep = ""), las = 1, tick = FALSE, hadj = 0.5, cex.axis = 1.5) axis(2, at = seq(0, 10, length.out = 6), labels = c("0.00", "2.00", "4.00", "6.00", "8.00", "10.00"), las = 1, cex.axis = 1.5) axis.Date(1, at = seq(min(vDates), max(vDates), "year"), cex.axis = 1.5) axis.Date(1, at = seq(min(vDates), max(vDates), "quarter"), labels = FALSE, tcl = -0.2) dev.off() ############################################## ## Appendix B ############################################## library("eDMA") data("SimData", package = "eDMA") newData <- c(NA, -0.07, 1.61, -2.07, 0.17, -0.80) SimData <- rbind(SimData, newData) Fit <- DMA(y ~ x2 + x3 + x4 + x5 + x6, data = SimData, vDelta = seq(0.9, 1.0, 0.01)) getLastForecast(Fit)