Estimation
The main function of the BGVAR package is its bgvar function. The unique feature of this toolbox is that we use Bayesian shrinkage priors with optionally stochastic volatility to estimate the country models in the GVAR. In its current version, three priors for the country VARs are implemented:
- Non-conjugate Minnesota prior (
MN, Litterman 1986; Koop and Korobilis 2010) - Stochastic Search Variable Selection prior (
SSVS, George, Sun, and Ni 2008) - Normal-Gamma prior (
NG, Huber and Feldkircher 2019)
The first two priors are described in more detail in Crespo Cuaresma, Feldkircher, and Huber (2016). For a more technical description of the Normal-Gamma prior see Huber and Feldkircher (2019) and for an application in the GVAR context Feldkircher and Siklos (2019). For the variances we can assume homoskedasticity or time variation (stochastic volatility). For the latter, the library relies on the stochvol package of Kastner (2016).
We start with estimating our toy model using the NG prior, the reduced eerData data set and the adjusted W.trade0012 weight matrix:
model.1<-bgvar(Data=eerData,
W=W.trade0012,
saves=100,
burns=100,
plag=1,
prior="NG",
hyperpara=NULL,
SV=TRUE,
thin=1,
trend=TRUE,
h=0,
save.country.store=FALSE,
multithread=FALSE,
eigen=1.05
)
The default prior specification in bgvar is to use the NG prior with stochastic volatility and one lag for both the endogenous and weakly exogenous variables (plag=1). In general, due to its high cross-sectional dimension, the GVAR can allow for very complex univariate dynamics and it might thus not be necessary to increase the lag length considerably as in a standard VAR (Burriel and Galesi 2018).
The setting hyperpara=NULL implies that we use the standard hyperparameter specification for the NG prior; see the helpfiles for more details. The slots saves and burns denote the posterior draws and the burn-in draws (i.e., the draws that are discarded). To ensure that the MCMC estimation has converged, a high-number of burn-ins is recommended (say 15,000 to 30,000). Saving the full set of posterior draws can eat up a lot of storage. To reduce this, we can use a thinning interval which stores only a thin\(^{th}\) draw of the global posterior output. For example, with thin=10 and saves=5000 posterior draws, the amount of MCMC draws stored is 500. TREND=TRUE implies that the model is estimated using a trend. Note that regardless of the trend specification, each equation always automatically includes an intercept term. To speed up computation, it is possible to set multithread=TRUE, which invokes parallel computing using the doParallel package. The package then attempts to detect the number of available CPU cores which it can use for parallel computing; ideally this number would be equal to the number of entities \(N\). In addition to storing the global posterior, we could also be interested in inspecting the output of the \(N\) country models in more detail. To do so, we could set save.country.store=TRUE allows to save the whole posterior distribution of the single country models, in the case we would like to inspect them. Due to storage reasons, the default is set to FALSE and only the posterior medians of the single country models are reported. Note that even in case save.country.store=FALSE, the whole posterior distribution of the is reported.
With have estimated the above model with stochastic volatility (SV=TRUE). There are several reasons why one may want to let the residual variances change over time. First and foremost, most time periods used in macroeconometrics are nowadays rather volatile including severe recessions. Hence accounting for time variation might improve the fit of the model (Primiceri 2005; Sims and Zha 2006; Dovern, Feldkircher, and Huber 2016; Huber 2016). Second, the specification implemented in the toolbox nests the homoskedastic case. It is thus a good choice to start with the more general case when first confronting the model with the data. For structural analysis such as the calculation of impulse responses, we take the variance covariance matrix with the median volatilities (over the sample period) on its diagonal. If we want to look at the volatilities of the first equation (y) in the euro area country model, we can type:
To discard explosive draws, we can compute the eigenvalues of the reduced form of the global model, written in its companion form. Unfortunately, this can only be done once the single models have been estimated and stacked together (and hence not directly built into the MCMC algorithm for the country models). To discard draws that lead to higher eigenvalues than 1.05 set eigen=1.05. We can look at the 10 largest eigenvalues by typing:
## [1] 1.0082184 0.9637162 1.0415077 1.0226959 1.0122765 1.0189337 1.0249688
## [8] 1.0476392 1.0022260 1.0105412Last, we have used the default option h=0, which implies that we use the full sample period to estimate the GVAR. For the purpose of forecast evaluation, h could be specified to a positive number, which then would imply that the last h observations are reserved as a hold-out sample and not used to estimate the model.
Model output and diagnostic checks
Having estimated the model, we can have summarize the outcome in various ways.
First, we can use theprint method
## ---------------------------------------------------------------------------------------
## Model Info:
## Prior: NG
## Nr. of lags: 1
## Nr. of posterior draws: 100/1=100
## Size of GVAR object: 1.1 Mb
## Trimming leads to 90 (90%) stable draws out of 100 total draws
## ---------------------------------------------------------------------------------------
## Model specification:
## $EA
## [1] "y, Dp, rer, stir, ltir, tb, y*, Dp*, rer*, stir*, ltir*, tb*, poil**, trend"
##
## $US
## [1] "y, Dp, rer, stir, ltir, tb, poil, y*, Dp*, rer*, stir*, ltir*, tb*, trend"
##
## $RU
## [1] "y, Dp, rer, stir, tb, y*, Dp*, rer*, stir*, ltir*, tb*, poil**, trend"
This just prints the submitted arguments of the bgvar object along with the model specification for each unit. The asterisks indicate weakly exogenous variables, double asterisks exogenous variables and variables without asterisks the endogenous variables per unit.
The summary method is a more enhanced way to analyze the output. It computes descriptive statistics like convergence properties of the MCMC chain, serial autocorrelation in the errors and the average pairwise autocorrelation of cross-unit residuals.
## ---------------------------------------------------------------------------------------
## Model Info:
## Prior: NG
## Nr. of lags: 1
## Nr. of posterior draws: 100/1=100
## Number of stable posterior draws: 90
## Number of countries: 3
## ---------------------------------------------------------------------------------------
## Convergence diagnostics
## Geweke statistic: 137 out of 1620 variables' z-values exceed the 1.96 threshold (8.46%)
## ---------------------------------------------------------------------------------------
## Global Likelihood: 4420.61
## F-test, first order serial autocorrelation of cross-country residuals
## Summary statistics:
## # p-values in %
## >0.1 "7" "38.89%"
## 0.05-0.1 "1" "5.56%"
## 0.01-0.05 "3" "16.67%"
## <0.01 "7" "38.89%"
## --------------------------------------------------------------------------------------
## Average pairwise cross-country correlation of country model residuals
## Summary statistics:
## y Dp rer stir ltir tb
## <0.1 "0 (0%)" "3 (100%)" "3 (100%)" "2 (66.67%)" "2 (100%)" "3 (100%)"
## 0.1-0.2 "0 (0%)" "0 (0%)" "0 (0%)" "1 (33.33%)" "0 (0%)" "0 (0%)"
## 0.2-0.5 "3 (100%)" "0 (0%)" "0 (0%)" "0 (0%)" "0 (0%)" "0 (0%)"
## >0.5 "0 (0%)" "0 (0%)" "0 (0%)" "0 (0%)" "0 (0%)" "0 (0%)"
## --------------------------------------------------------------------------------------
We can now have a closer look at the information provided by summary. The header contains some basic information about the prior used to estimate the model, how many lags, posterior draws and countries. The next line shows Geweke’s CD statistic, which is calculated using the coda package. Geweke’s CD assesses practical convergence of the MCMC algorithm. In a nutshell, the diagnostic is based on a test for equality of the means of the first and last part of a Markov chain (by default we use the first 10% and the last 50%). If the samples are drawn from the stationary distribution of the chain, the two means are equal and Geweke’s statistic has an asymptotically standard normal distribution.
The test statistic is a standard Z-score: the difference between the two sample means divided by its estimated standard error. The standard error is estimated from the spectral density at zero and so takes into account any autocorrelation. The test statistic shows that only a small fraction of all coefficients did not convergence. Increasing the number of burn ins can help decreasing this number further. The statistic can also be calculated by typing conv.diag(model.1).
The next model statistic is the likelihood of the global model. This statistic can be used for model comparison. Next and to assess, whether there is first order serial autocorrelation present, we provide the results of a simple F-test. The table shows the share of p-values that fall into different significance categories. Since the null hypothesis is that of no serial correlation, we would like to have as many large (\(>0.1\)) p-values as possible. The statistics show that already with one lag, serial correlation is modest in most equations’ residuals. This could be the case since we have estimated the unit models with stochastic volatility. To further decrease serial correlation in the errors, one could increase the number of lags via plag.
The last part of the summary output contains a statistic of cross-unit correlation of (posterior median) residuals. One assumption of the GVAR framework is that of negligible, cross-unit correlation of the residuals. Significant correlations prohibit structural and spillover analysis (Dees et al. 2007). In this example, correlation is reasonably small.
Some other useful methods the bgvar toolbox offers contain the coef (or coefficients as its alias) methods to extract the \(k \times k \times plag\) matrix of reduced form coefficients of the global model. Via the vcov command, we can access the global variance covariance matrix and the logLik function allows us to gather the global log likelihood (as provided by the summary command).
Last, we can have a quick look at the in-sample fit using either the posterior median of the country models’ residuals (global=FALSE) or those of the global solution of the GVAR (global=TRUE). The in-sample fit can also be extracted by using fitted().
Here, we show the in-sample fit of the euro area model (global=FALSE).
In-sample fit for euro area variables
We can estimate the model with two further priors on the unit models, the SSVS prior and the Minnesota prior. To give a concrete example, the SSVS prior can be invoked by typing:
model.ssvs.1<-bgvar(Data=eerData,
W=W.trade0012,
saves=100,
burns=100,
plag=1,
prior="SSVS",
hyperpara=NULL,
SV=TRUE,
thin=1,
trend=TRUE,
h=0,
save.country.store=FALSE,
multithread=FALSE,
eigen=1.05
)One feature of the SSVS prior is that it allows to look at the posterior inclusion probabilities to gauge the importance of particular variables. For example, we can have a look at the PIPs of the euro area model by typing:
## y Dp rer stir ltir tb
## y_lag1 1.00 0.17 0.28 0.06 0.09 0.69
## Dp_lag1 0.06 0.24 0.08 0.27 0.51 0.10
## rer_lag1 0.21 0.49 1.00 1.00 0.16 0.50
## stir_lag1 0.15 0.22 0.34 1.00 0.47 0.94
## ltir_lag1 0.13 1.00 0.34 0.42 1.00 0.96
## tb_lag1 0.62 0.30 0.85 0.08 0.36 1.00
## y* 1.00 0.17 0.18 0.09 0.08 0.12
## Dp* 0.15 0.31 0.17 0.20 0.16 0.18
## rer* 0.34 0.38 1.00 0.20 0.06 0.46
## stir* 0.08 0.92 0.13 0.36 0.09 0.13
## ltir* 0.48 0.16 0.19 0.96 0.04 0.24
## tb* 0.14 1.00 0.31 0.42 0.15 0.10
## poil** 0.40 1.00 0.24 0.09 0.30 0.28
## y*_lag1 0.61 0.05 0.02 0.11 0.12 0.49
## Dp*_lag1 0.54 0.23 0.41 0.15 0.71 0.15
## rer*_lag1 0.21 0.16 0.92 0.05 0.28 0.42
## stir*_lag1 0.20 0.26 0.11 0.53 0.19 0.04
## ltir*_lag1 1.00 0.76 0.35 1.00 0.17 0.28
## tb*_lag1 0.35 1.00 1.00 0.51 0.10 0.86
## poil**_lag1 0.85 0.53 0.13 0.10 0.07 0.27
## cons 0.40 0.27 1.00 0.06 0.10 1.00
## trend 0.30 0.91 0.04 0.03 0.08 0.97This prints in the columns, the equations in the EA country model and in the rows the included variables. The example shows that besides other variables, the trade balance (tb) is an important determinant of the real exchange rate (rer).
We can also have a look at the average of the PIPs across all units:
## y Dp rer stir ltir tb poil
## y_lag1 1.0000000 0.19666667 0.28666667 0.35666667 0.075 0.6400000 0.04
## Dp_lag1 0.2833333 0.69000000 0.11000000 0.17000000 0.365 0.3100000 0.93
## rer_lag1 0.1733333 0.31333333 1.00000000 0.52333333 0.295 0.7566667 0.17
## stir_lag1 0.7100000 0.16000000 0.20333333 1.00000000 0.705 0.5600000 0.13
## ltir_lag1 0.2800000 0.54000000 0.41000000 0.35000000 1.000 0.6450000 0.15
## tb_lag1 0.2766667 0.29666667 0.39000000 0.25333333 0.255 1.0000000 0.09
## y* 0.4933333 0.13000000 0.42000000 0.17333333 0.100 0.4033333 0.13
## Dp* 0.2600000 0.50000000 0.40000000 0.09333333 0.325 0.3433333 1.00
## rer* 0.2333333 0.32666667 0.79000000 0.21333333 0.335 0.3466667 0.25
## stir* 0.3366667 0.39333333 0.14333333 0.23333333 0.235 0.3600000 0.43
## ltir* 0.2366667 0.24333333 0.10666667 0.43666667 0.170 0.1800000 0.12
## tb* 0.2466667 0.47666667 0.26333333 0.34333333 0.210 0.2400000 0.14
## poil** 0.7000000 0.55000000 0.20000000 0.28000000 0.300 0.1900000 NaN
## y*_lag1 0.2966667 0.05666667 0.41666667 0.34333333 0.115 0.3033333 0.03
## Dp*_lag1 0.2666667 0.23333333 0.34333333 0.10333333 0.495 0.4100000 0.14
## rer*_lag1 0.6266667 0.19000000 0.62333333 0.06666667 0.475 0.2900000 0.56
## stir*_lag1 0.5033333 0.20666667 0.07333333 0.49333333 0.330 0.2733333 0.31
## ltir*_lag1 0.4300000 0.39666667 0.29333333 0.40666667 0.145 0.2500000 0.20
## tb*_lag1 0.3200000 0.46333333 0.53000000 0.30666667 0.225 0.4566667 0.13
## poil**_lag1 0.4500000 0.35000000 0.25500000 0.19500000 0.070 0.1750000 NaN
## cons 0.3066667 0.18000000 0.64000000 0.18666667 0.130 0.7833333 0.05
## trend 0.3133333 0.42333333 0.36666667 0.51000000 0.105 0.4800000 0.95
## poil_lag1 0.3500000 0.26000000 0.09000000 1.00000000 0.140 0.3500000 1.00This shows that the same determinants for the real exchange rate appear also as important regressors in other country models.
Different specifications of the model
In this section we explore different specifications of the structure of the GVAR model. Other, specification choices that relate more to the time series properties of the data, such as specifying different lags and priors are left for the reader to explore. We will use the SSVS prior and judge the different specifications by examining the posterior inclusion probabilities.
As a first modification, we could use different weights for different variable classes as proposed in Eickmeier and Ng (2015). For example we could use financial weights to construct weakly exogenous variables of financial factors and trade weights for real variables.
The eerData set provides us with a list of different weight matrices that are described in the help files.
Now we specify the sets of variables to be weighted:
eerData2<-eerData
variable.list<-list();variable.list$real<-c("y","Dp","tb");variable.list$fin<-c("stir","ltir","rer")We can then re-estimate the model:
# weights for first variable set tradeW.0012, for second finW0711
model.ssvs.2<-bgvar(Data=eerData2,
W=W.list[c("tradeW.0012","finW0711")],
plag=1,
saves=100,
burns=100,
prior="SSVS",
SV=TRUE,
thin=1,
hyperpara=NULL,
eigen=TRUE,
variable.list=variable.list,
OE.weights=NULL,
Wex.restr=NULL,
trend=TRUE,
save.country.store=FALSE
)Another specification would be to include a foreign variable only when its domestic counterpart is missing. For example, when working with nominal bilateral exchange rates we probably do not want to include also its weighted average (which corresponds to something like an effective exchange rate). Using the previous model we could place a restriction on long-term interest rates, since those are not contained in all country models.
# does include ltir* only when ltir is missing domestically
model.ssvs.3<-bgvar(Data=eerData,
W=W.trade0012,
plag=1,
saves=100,
burns=100,
prior="SSVS",
SV=TRUE,
thin=1,
hyperpara=NULL,
eigen=TRUE,
variable.list=NULL,
OE.weights=NULL,
Wex.restr="ltir",
trend=TRUE,
save.country.store=FALSE
)## ---------------------------------------------------------------------------------------
## Model Info:
## Prior: SSVS
## Nr. of lags: 1
## Nr. of posterior draws: 100/1=100
## Size of GVAR object: 1 Mb
## Trimming leads to 80 (80%) stable draws out of 100 total draws
## ---------------------------------------------------------------------------------------
## Model specification:
## $EA
## [1] "y, Dp, rer, stir, ltir, tb, y*, Dp*, rer*, stir*, tb*, poil**, trend"
##
## $US
## [1] "y, Dp, rer, stir, ltir, tb, poil, y*, Dp*, rer*, stir*, tb*, trend"
##
## $RU
## [1] "y, Dp, rer, stir, tb, y*, Dp*, rer*, stir*, ltir*, tb*, poil**, trend"Last, we could also use a different specification of oil prices in the model. Currently, the oil price is determined endogenously within the US model. Alternatively, one could set up an own standing oil price model with additional variables that feeds the oil price back into the other economies as exogenous variable (Mohaddes and Raissi 2019).
The model structure would then look something like in the Figure below:
“GVAR with oil prices modeled separately.”
For that purpose we have to remove oil prices from the US model and attach them to a separate slot in the data list. This slot has to have its own country label. We use ‘OC’ for “oil country”.
eerData2$OC<-eerData$US[,c("poil"),drop=FALSE] # move oil prices into own slot
eerData2$US<-eerData$US[,c("y","Dp", "rer" , "stir", "ltir","tb")] # exclude it from US m odel
Now we have to specify a list object that we label OC.weights. The list has to consist of three slots with the following names weights, variables and exo:
OC.weights<-list()
OC.weights$weights<-rep(1/3, 3)
names(OC.weights$weights)<-names(eerData2)[1:3] # last one is OC model, hence only until 3
OC.weights$variables<-c(colnames(eerData2$OC),"y") # first entry, endog. variables, second entry weighted average of y from the other countries to proxy demand
OC.weights$exo<-"poil"
The first slot, weights, should be a vector of weights that sum up to unity. In the example above, we simply use \(1/N\), other weights could include purchasing power parities (PPP). The weights are used to aggregate specific variables that in turn enter the oil model as weakly exogenous. The second slot, variables, should specify the names of the endogenous and weakly exogenous variables that are used in the model. In the oil price example, we include the oil price (poil) as endogenous variable (not contained in any other country model); a weighted average using weights of output (y) to proxy world demand is included as weakly exogenous variable. Last, we specify via exo which one of the endogenous variables of the oil price model are fed back into the other country models. In this example we specify poil. Last, we include OC.weights in a further list called OE.weights (other entity weights):
where the list entry has to has the same name as the oil country model, in our example OC. Note that specifying OE.weights in principle allows to include further “extra” country models by simply adding additional list slots.
Now we can re-estimate the model by specifying OE.weights, which is a list with its entries labelled exactly like the additional country models. The structure is designed to also allow for more than one other entity (henceforth labelled OE).
model.ssvs.4<-bgvar(Data=eerData2,
W=W.trade0012,
plag=1,
saves=100,
burns=100,
prior="SSVS",
SV=TRUE,
thin=1,
hyperpara=NULL,
variable.list=NULL,
OE.weights=OE.weights,
trend=TRUE,
save.country.store=FALSE
)and can compare the results of the four models by e.g., looking at the average PIPs.
aux1<-model.ssvs.1$cc.results$PIP$PIP.avg;aux1<-aux1[-nrow(aux1),1:6]
aux2<-model.ssvs.2$cc.results$PIP$PIP.avg;aux2<-aux2[-nrow(aux2),1:6]
aux3<-model.ssvs.3$cc.results$PIP$PIP.avg;aux3<-aux3[-nrow(aux3),1:6]
aux4<-model.ssvs.4$cc.results$PIP$PIP.avg;aux4<-aux4[-nrow(aux4),1:6]heatmap(aux1,Rowv=NA,Colv=NA, main="Model 1")
heatmap(aux2,Rowv=NA,Colv=NA, main="Model 2")
heatmap(aux3,Rowv=NA,Colv=NA, main="Model 3")
heatmap(aux4,Rowv=NA,Colv=NA, main="Model 4")
Heatmaps of PIPs.
We could also compare the models based on their fit, the likelihood, information criteria such as the DIC, residual properties or their forecasting performance.
