update to 0.11.0

kthohr · Jun 24, 2018 · 703d5e8 · 703d5e8
1 parent e5fb9d2
commit 703d5e8
Show file tree

Hide file tree

Showing 18 changed files with 599 additions and 41 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: BMR
 Type: Package
 Title: Bayesian Macroeconometrics in R
-Version: 0.10.0
-Date: 2018-03-08
+Version: 0.11.0
+Date: 2018-06-24
 Author: Keith O'Hara <[email protected]>
 Maintainer: Keith O'Hara <[email protected]>
 Description: A package for estimating Bayesian macroeconometric models.
@@ -17,7 +17,8 @@ Imports:
     Rcpp,
     grid
 NeedsCompilation: yes
-RcppModules: bvarm_module, bvars_module, bvarcnw_module, bvarinw_module,
+RcppModules: bvarm_module, bvars_module,
+             bvarcnw_module, bvarinw_module,
              bvartvp_module, cvar_module,
              gensys_module, uhlig_module, 
              dsge_gensys_module, dsge_uhlig_module, 

diff --git a/NAMESPACE b/NAMESPACE
@@ -17,7 +17,7 @@ exportClasses("Rcpp_uhlig")
 exportClasses("Rcpp_dsge_gensys")
 exportClasses("Rcpp_dsgevar_gensys")
 
-export(FEVD,IRF,forecast,mode_check,prior,gtsplot)
+export(FEVD,IRF,IRFcomp,forecast,mode_check,prior,gtsplot)
 
 exportPattern("^[[:alpha:]]+")
 

diff --git a/R/IRF.R b/R/IRF.R
@@ -138,14 +138,14 @@ IRF.Rcpp_dsgevar_uhlig <- function(obj,periods,var_names=NULL,percentiles=c(.05,
         }
     }
 
-    rm("irf_temp")
-
-    irf_tess <- apply(irf_tess,c(3,1,2),sort)
-    irf_tess <- aperm(irf_tess,c(2,3,1,4))
+    irf_tess_1 <- apply(irf_tess,c(3,1,2),sort) # fix annoying bug
+    irf_tess <- aperm(irf_tess_1,c(2,3,1,4))
+
+    rm("irf_temp","irf_tess_1")
 
-    irf_upper <- round(percentiles[3]*n_draws)
+    irf_upper <- min(round(percentiles[3]*n_draws),n_draws)
     irf_mid <- round(percentiles[2]*n_draws)
-    irf_lower <- round(percentiles[1]*n_draws)
+    irf_lower <- max(round(percentiles[1]*n_draws),1)
 
     #
 
@@ -350,6 +350,336 @@ IRF.Rcpp_dsgevar_uhlig <- function(obj,periods,var_names=NULL,percentiles=c(.05,
     return=list(plot_vals=plot_vals)
 }
 
+IRFcomp <- function(obj_1, obj_2, periods=10, cumulative=FALSE, cumul_inds=NULL, var_names=NULL, percentiles=c(.05,.50,.95), 
+                    which_shock=NULL, which_response=NULL, shocks_row_order=TRUE, save=FALSE, save_format=c("pdf","eps"), 
+                    save_title=NULL, height=13, width=13)
+{
+    #
+
+    if (periods <= 0) {
+        stop("error: need periods to be > 0")
+    }
+
+    #
+
+    if (obj_1$M != obj_2$M) {
+        stop("error: models need to be of the same data dimensions")
+    }
+
+    #
+
+    M <- obj_1$M
+    n_draws <- dim(obj_1$beta_draws)[3]
+
+    irf_temp <- obj_1$IRF(periods)$irf_vals
+
+    # put the IRFs in a tesseract-type format
+
+    irf_tess <- array(NA,dim=c(M,M,periods,n_draws))
+
+    for (i in 1:n_draws) {
+        irf_tess[,,,i] <- irf_temp[,,((i-1)*periods+1):(i*periods)]
+    }
+
+    if (cumulative)
+    {
+        if (is.null(cumul_inds)) {
+            cumul_inds <- 1:M
+        }
+
+        for (i in 1:n_draws) {
+            irfs_draw_i <- irf_tess[,,,i]
+
+            for (ll in cumul_inds) {
+                for (jj in 2:periods) {
+                    irfs_draw_i[ll,,jj] <- irfs_draw_i[ll,,jj] + irfs_draw_i[ll,,jj-1]
+                }
+            }
+
+            irf_tess[,,,i] <- irfs_draw_i
+        }
+    }
+
+    irf_tess_1 <- apply(irf_tess,c(3,1,2),sort) # fix annoying bug
+    irf_tess <- aperm(irf_tess_1,c(2,3,1,4))
+
+    rm("irf_temp","irf_tess_1")
+
+    irf_upper <- min(round(percentiles[3]*n_draws),n_draws)
+    irf_mid <- round(percentiles[2]*n_draws)
+    irf_lower <- max(round(percentiles[1]*n_draws),1)
+
+    #
+
+    plot_vals_1 <- array(NA,dim=c(periods,4,M,M))
+    IRFPData <- 0
+
+    for (i in 1:M) {
+        for (k in 1:M) {
+            IRFPData <- data.frame(irf_tess[,k,irf_lower,i],irf_tess[,k,irf_mid,i],irf_tess[,k,irf_upper,i],1:(periods))
+            IRFPData <- as.matrix(IRFPData)
+            plot_vals_1[,,k,i] <- IRFPData
+        }
+    }
+
+    #
+    # Model 2
+
+    M <- obj_2$M
+    n_draws <- dim(obj_2$beta_draws)[3]
+
+    irf_temp <- obj_2$IRF(periods)$irf_vals
+
+    # put the IRFs in a tesseract-type format
+
+    irf_tess <- array(NA,dim=c(M,M,periods,n_draws))
+
+    for (i in 1:n_draws) {
+        irf_tess[,,,i] <- irf_temp[,,((i-1)*periods+1):(i*periods)]
+    }
+
+    if (cumulative)
+    {
+        if (is.null(cumul_inds)) {
+            cumul_inds <- 1:M
+        }
+
+        for (i in 1:n_draws) {
+            irfs_draw_i <- irf_tess[,,,i]
+
+            for (ll in cumul_inds) {
+                for (jj in 2:periods) {
+                    irfs_draw_i[ll,,jj] <- irfs_draw_i[ll,,jj] + irfs_draw_i[ll,,jj-1]
+                }
+            }
+
+            irf_tess[,,,i] <- irfs_draw_i
+        }
+    }
+
+    irf_tess_1 <- apply(irf_tess,c(3,1,2),sort) # fix annoying bug
+    irf_tess <- aperm(irf_tess_1,c(2,3,1,4))
+
+    rm("irf_temp","irf_tess_1")
+
+    irf_upper <- min(round(percentiles[3]*n_draws),n_draws)
+    irf_mid <- round(percentiles[2]*n_draws)
+    irf_lower <- max(round(percentiles[1]*n_draws),1)
+
+    #
+
+    plot_vals_2 <- array(NA,dim=c(periods,4,M,M))
+    IRFPData <- 0
+
+    for (i in 1:M) {
+        for (k in 1:M) {
+            IRFPData <- data.frame(irf_tess[,k,irf_lower,i],irf_tess[,k,irf_mid,i],irf_tess[,k,irf_upper,i],1:(periods))
+            IRFPData <- as.matrix(IRFPData)
+            plot_vals_2[,,k,i] <- IRFPData
+        }
+    }
+
+    #
+
+    vplayout <- function(x,y){viewport(layout.pos.row=x, layout.pos.col=y)}
+
+    #
+
+    if (is.null(which_shock)) {
+        which_shock <- 1:M
+    }
+
+    if (is.null(which_response)) {
+        which_response <- 1:M
+    }
+
+    n_response <- length(which_response)
+    n_shocks   <- length(which_shock)
+
+    if (class(var_names) != "character") {
+        var_names <- character(length=M)
+        for (i in 1:M) {  
+            var_names[i] <- paste("VAR",i,sep="")
+        }
+    }
+
+    #
+    # plot IRFs
+
+    save_format <- match.arg(save_format)
+
+    IRFL <- IRFM <- IRFU <- Time <- NULL # CRAN check workaround
+
+    if (n_response == M && n_shocks == M) {
+
+        if (save==TRUE)
+        {
+            if(class(dev.list()) != "NULL"){dev.off()}
+
+            save_name <- ""
+            if (!is.null(save_title)) {
+                save_name <- paste(save_title,".",save_format,sep="")
+            } else {
+                save_name <- paste("IRFs.",save_format,sep="")
+            }
+
+            if (save_format=="eps") {
+                cairo_ps(filename=save_name,height=height,width=width)
+            } else {
+                cairo_pdf(filename=save_name,height=height,width=width)
+            }
+        }
+
+        grid.newpage()
+        pushViewport(viewport(layout=grid.layout(M,M)))
+
+        for (i in 1:M) {
+            for (k in 1:M) {
+                NameResponse <- var_names[k]
+                NameImpulse  <- var_names[i]
+
+                IRFDF_1 <- plot_vals_1[,,k,i]
+                #IRFDF_1 <- data.frame(IRFDF_1)
+
+                IRFDF_2 <- plot_vals_2[,,k,i]
+                IRFDF <- data.frame(IRFDF_1,IRFDF_2)
+
+                #colnames(IRFDF_1) <- c("IRFL","IRFM","IRFU","Time")
+                colnames(IRFDF) <- c("IRFL1","IRFM1","IRFU1","Time1","IRFL2","IRFM2","IRFU2","Time2")
+
+                #
+
+                gg1 <- ggplot(data=(IRFDF),aes(x=Time1)) + xlab("") + ylab(paste("Shock from ",NameImpulse," to", NameResponse)) + geom_hline(yintercept=0)
+                gg2 <- gg1 + geom_line(aes(y=IRFM1),color="blue",size=2) + geom_line(aes(y=IRFL1),color="blue",size=1,linetype=4) + geom_line(aes(y=IRFU1),color="blue",size=1,linetype=4)
+                gg3 <- gg2 + geom_line(aes(y=IRFM2),color="darkred",size=2) + geom_line(aes(y=IRFL2),color="darkred",size=1,linetype=4) + geom_line(aes(y=IRFU2),color="darkred",size=1,linetype=4)
+                gg4 <- gg3 + theme(panel.background = element_rect(fill='white', colour='grey5')) + theme(panel.grid.major = element_line(colour = 'grey89'))
+
+                if (shocks_row_order) {
+                    print(gg4,vp = vplayout(i,k))
+                } else {
+                    print(gg4,vp = vplayout(k,i))
+                }
+
+                #
+
+                Sys.sleep(0.3)
+            }
+        }
+
+        if (save==TRUE) {dev.off()}
+
+    } else {
+
+        if (n_response < 4) {
+            MR <- n_response; MC <- 1
+        } else if (n_response == 4) {
+            MR <- 2; MC <-2
+        } else if (n_response > 4 && n_response < 7) {
+            MR <- 3; MC <- 2
+        } else if (n_response > 6 && n_response < 10) {
+            MR <- 3; MC <- 3
+        } else if (n_response > 9 && n_response < 13) {
+            MR <- 4; MC <- 3
+        }else if (n_response > 12 && n_response < 17) {
+            MR <- 4; MC <- 4
+        }else if (n_response > 17 && n_response < 21) {
+            MR <- 5; MC <- 4
+        }else if (n_response > 20 && n_response < 26) {
+            MR <- 5; MC <- 5
+        } else if (n_response > 25 && n_response < 31) {
+            MR <- 5; MC <- 6
+        } else if (n_response > 30 && n_response < 37) {
+            MR <- 6; MC <- 6
+        } else {
+            stop("You have too many IRFs to plot!")
+        }
+
+        #
+
+        for (i in which_shock) {
+            plot_ind_r <- 1
+            plot_ind_c <- 1
+
+            if (save==TRUE) {
+                if(class(dev.list()) != "NULL"){dev.off()}
+
+                if (n_shocks==1) {
+                    cairo_pdf(filename="IRFs.pdf",height=height,width=width)
+                } else {
+                    SaveIRF <- paste(var_names[i],"_Shock",".pdf",sep="")
+                    cairo_pdf(filename=SaveIRF,height=height,width=width)
+                }
+            }
+
+            if (save==TRUE)
+            {
+                if(class(dev.list()) != "NULL"){dev.off()}
+
+                save_name <- ""
+                if (n_shocks==1) {
+                    if (!is.null(save_title)) {
+                        save_name <- paste(save_title,".",save_format,sep="")
+                    } else {
+                        save_name <- paste("IRFs.",save_format,sep="")
+                    }
+                } else {
+                    save_name <- paste(var_names[i],"_Shock",".",save_format,sep="")
+                }
+
+                if (save_format=="eps") {
+                    cairo_ps(filename=save_name,height=height,width=width)
+                } else {
+                    cairo_pdf(filename=save_name,height=height,width=width)
+                }
+            }
+
+            grid.newpage()
+            pushViewport(viewport(layout=grid.layout(MR,MC)))
+
+            for (k in which_response) {
+                NameResponse <- var_names[k]
+                NameImpulse  <- var_names[i]
+
+                IRFDF_1 <- plot_vals_1[,,k,i]
+                #IRFDF_1 <- data.frame(IRFDF_1)
+
+                IRFDF_2 <- plot_vals_2[,,k,i]
+                IRFDF <- data.frame(IRFDF_1,IRFDF_2)
+
+                #colnames(IRFDF_1) <- c("IRFL","IRFM","IRFU","Time")
+                colnames(IRFDF) <- c("IRFL1","IRFM1","IRFU1","Time1","IRFL2","IRFM2","IRFU2","Time2")
+
+                #
+
+                gg1 <- ggplot(data=(IRFDF),aes(x=Time1)) + xlab("") + ylab(paste("Shock from ", NameImpulse," to", NameResponse)) + geom_hline(yintercept=0)
+                gg2 <- gg1 + geom_line(aes(y=IRFM1),color="blue",size=2) + geom_line(aes(y=IRFL1),color="blue",size=1,linetype=4) + geom_line(aes(y=IRFU1),color="blue",size=1,linetype=4)
+                gg3 <- gg2 + geom_line(aes(y=IRFM2),color="darkred",size=2) + geom_line(aes(y=IRFL2),color="darkred",size=1,linetype=4) + geom_line(aes(y=IRFU2),color="darkred",size=1,linetype=4)
+                gg4 <- gg3 + theme(panel.background = element_rect(fill='white', colour='grey5')) + theme(panel.grid.major = element_line(colour = 'grey89'))
+
+                print(gg4,vp = vplayout(plot_ind_r,plot_ind_c))
+
+                Sys.sleep(0.3)
+
+                #
+
+                plot_ind_c <- plot_ind_c + 1
+
+                if (plot_ind_c > MC) {
+                    plot_ind_c <- 1
+                    plot_ind_r <- plot_ind_r + 1
+                }
+            }
+        }
+
+        if(save==TRUE){dev.off()}
+
+    }
+
+    #
+
+    return=list(plot_vals_1=plot_vals_1,plot_vals_2=plot_vals_2)
+}
+
 .irf_bvartvp <- function(obj, periods=10, which_irfs=NULL, var_names=NULL, percentiles=c(.05,.50,.95), 
                          which_shock=NULL, which_response=NULL, save=FALSE, height=13, width=13)
 {