############################################################################### ### Replication Material for ### BGVAR: Bayesian Global Vector Autoregressions with Shrinkage Priors in R ### Maximilian Boeck, Martin Feldkircher, and Florian Huber ### Journal of Statistical Software ### ### created by Maximilian Boeck 30/04/2022 ### created with Sweave ### ############################################################################### ### R code from vignette source 'jss4147.Rnw' ################################################### ### code chunk number 1: preliminaries ################################################### options(prompt = "R> ", continue = "+ ", width = 76, useFancyQuotes = FALSE) library("MASS") ################################################### ### code chunk number 2: main (eval = FALSE) ################################################### ## model.1 <- bgvar(Data, W, plag=1, draws=5000, burnin=5000, prior="NG", ## SV=TRUE, hold.out=0, thin=1, hyperpara=NULL, eigen=TRUE, ## Ex=NULL, trend=FALSE, expert=NULL, verbose=TRUE) ################################################### ### code chunk number 3: start ################################################### library("BGVAR") data("pesaranData", package = "BGVAR") set.seed(571) ################################################### ### code chunk number 4: model.1 ################################################### model.1 <- bgvar(Data = pesaranData, W = W.8016) ################################################### ### code chunk number 5: model.1.summary ################################################### summary(model.1) # geweke: 19.6% ################################################### ### code chunk number 6: insample ################################################### yfit <- fitted(model.1) plot(model.1, global=FALSE, resp="US", cex.axis=1.2) ################################################### ### code chunk number 7: forecasting ################################################### model1.ssvs <- bgvar(Data=pesaranData,W=W.8016,plag=1,hold.out=8,thin=2, draws=1000,burnin=1000, prior="SSVS") model1.ng <- bgvar(Data=pesaranData,W=W.8016,plag=1,hold.out=8,thin=2, draws=1000,burnin=1000, prior="NG") model1.mn <- bgvar(Data=pesaranData,W=W.8016,plag=1,hold.out=8,thin=2, draws=1000,burnin=1000, prior="MN") model1.hs <- bgvar(Data=pesaranData,W=W.8016,plag=1,hold.out=8,thin=2, draws=1000,burnin=1000,prior="HS") model2.ssvs <- bgvar(Data=pesaranData,W=W.8016,plag=2,hold.out=8,thin=2, draws=1000,burnin=1000,prior="SSVS") model2.ng <- bgvar(Data=pesaranData,W=W.8016,plag=2,hold.out=8,thin=2, draws=1000,burnin=1000, prior="NG") model2.mn <- bgvar(Data=pesaranData,W=W.8016,plag=2,hold.out=8,thin=2, draws=1000,burnin=1000, prior="MN") model2.hs <- bgvar(Data=pesaranData,W=W.8016,plag=2,hold.out=8,thin=2, draws=1000,burnin=1000,prior="HS") ################################################### ### code chunk number 8: model.eval ################################################### fcast <- predict(model1.ssvs, n.ahead=8) lps.ssvs<-sum(lps(fcast)) ################################################### ### code chunk number 9: model.eval2 ################################################### fcast1.ng <- predict(model1.ng, n.ahead=8);lps.ng<-sum(lps(fcast1.ng)) fcast1.mn <- predict(model1.mn, n.ahead=8);lps.mn<-sum(lps(fcast1.mn)) fcast1.hs <- predict(model1.hs, n.ahead=8);lps.hs<-sum(lps(fcast1.hs)) fcast2.ssvs <- predict(model2.ssvs, n.ahead=8);lps.ssvs2<-sum(lps(fcast2.ssvs)) fcast2.ng <- predict(model2.ng, n.ahead=8);lps.ng2<-sum(lps(fcast2.ng)) fcast2.mn <- predict(model2.mn, n.ahead=8);lps.mn2<-sum(lps(fcast2.mn)) fcast2.hs <- predict(model2.hs, n.ahead=8);lps.hs2<-sum(lps(fcast2.hs)) dic1.ng<-dic(model1.ng);dic2.ng<-dic(model2.ng) dic1.mn<-dic(model1.mn);dic2.mn<-dic(model2.mn) dic1.ssvs<-dic(model1.ng);dic2.ssvs<-dic(model2.ssvs) dic1.hs<-dic(model1.hs);dic2.hs<-dic(model2.hs) ################################################### ### code chunk number 10: summary.table1 ################################################### library(xtable) summarize.lps<-rbind(c(lps.ssvs,lps.ng,lps.mn,lps.hs), c(lps.ssvs2,lps.ng2,lps.mn2,lps.hs2) ) colnames(summarize.lps)<-c("SSVS","NG","MN","HS") rownames(summarize.lps)<-c("p=1","p=2") xtable(summarize.lps, digits = 2, caption = "LPS scores",label="tbl:lps") ################################################### ### code chunk number 11: summary.table2 ################################################### summarize.dic<-rbind(c(dic1.ssvs,dic1.ng,dic1.mn,dic1.hs), c(dic2.ssvs,dic2.ng,dic2.mn,dic2.hs) ) colnames(summarize.dic)<-c("SSVS","NG","MN","HS") rownames(summarize.dic)<-c("p=1","p=2") rel.dic<-summarize.dic/summarize.dic["p=1","NG"] xtable(rel.dic, digits = 2, caption = "DIC relative to NG (p=1)",label="tbl:dic") ################################################### ### code chunk number 12: free.up.memory1 ################################################### rm(model1.hs,model1.mn,model1.ng,model1.ssvs,model2.hs,model2.mn,model2.ng,model2.ssvs) rm(fcast,fcast1.ng,fcast1.mn,fcast1.hs,fcast2.ssvs,fcast2.ng,fcast2.hs,fcast2.mn) rm(dic1.hs,dic1.mn,dic1.ng,dic1.ssvs,dic2.hs,dic2.mn,dic2.ng,dic2.ssvs) rm(lps.hs,lps.mn,lps.ng,lps.ssvs,lps.hs2,lps.ng2,lps.mn2,lps.ssvs2) rm(summarize.lps,summarize.dic) gc() ################################################### ### code chunk number 13: gfevd_allgemein ################################################### gfevd.1=gfevd(model.1,n.ahead=24,running=TRUE)$FEVD ################################################### ### code chunk number 14: gfevd.CA ################################################### idx<-which(grepl("CA.",dimnames(gfevd.1)[[2]])) own<-colSums(gfevd.1["CA.y",idx,]) foreign<-colSums(gfevd.1["CA.y",-idx,]) ################################################### ### code chunk number 15: Canada_barplot ################################################### par(mfrow=c(1,1), mar=c(3.1, 4.1, 3.1, 8.1), xpd=TRUE) barplot(t(cbind(own,foreign)),legend.text =c("own","foreign"), args.legend = list(bty="n",x=40,x.intersp=0.1)) box() ################################################### ### code chunk number 16: gfevd_CA ################################################### idx<-which(grepl(".y",dimnames(gfevd.1)[[2]])) # get index of output variables xx<-dimnames(gfevd.1)[[2]][idx] xx<-substr(xx,12,13) # get country names shares<-matrix(0,ncol=2,nrow=length(xx)) colnames(shares)<-c("Impact","Long-run");rownames(shares)<-xx for(i in 1:length(xx)){ idx<-which(grepl(xx[i],dimnames(gfevd.1)[[2]])) own<-colSums(gfevd.1[paste0(xx[i],".y"),idx,]) shares[i,]<-own[c(1,24)] } ################################################### ### code chunk number 17: Canada_longrun ################################################### par(mfrow=c(2,1),mar=c(3,3,2,2)) barplot(sort(shares[,"Impact"]),las=2,ylim=c(0,1.1),main="h=0");box() barplot(sort(shares[,"Long-run"]),las=2,ylim=c(0,1.1),main="h=24");box() ################################################### ### code chunk number 18: do_cholesky ################################################### shockinfo_chol <- get_shockinfo("chol") shockinfo_chol$shock <- "US.eq" shockinfo_chol$scale <- -1 irf_chol <- irf(model.1, n.ahead = 24, shockinfo = shockinfo_chol) ################################################### ### code chunk number 19: chol_plot1 ################################################### plot(irf_chol, resp=c("US.eq","US.y")) ################################################### ### code chunk number 20: chol_plot2 ################################################### plot(irf_chol, resp=c("CH.y","TR.y","CA.y","DE.y")) ################################################### ### code chunk number 21: free.up.memory2 ################################################### rm(irf_chol) ################################################### ### code chunk number 22: do_sign_restrictions ################################################### shockinfo_sign <- get_shockinfo("sign") shockinfo_sign <- add_shockinfo(shockinfo_sign, shock="US.eq", restriction=c("US.eq","US.y"), sign=c(">",">"), horizon=c(4,4), prob=1, scale=-1) irf_sign <- irf(model.1, n.ahead = 24, shockinfo=shockinfo_sign) ################################################### ### code chunk number 23: sign_plot ################################################### plot(irf_sign, resp=c("US.eq","US.y"), shock="US.eq") ################################################### ### code chunk number 24: free.up.memory3 ################################################### rm(irf_sign)