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 Martin
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,
draws=100,
burnin=100,
plag=1,
prior="NG",
hyperpara=NULL,
SV=TRUE,
thin=1,
trend=TRUE,
hold.out=0,
eigen=1
)
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.
Other standard specifications that should be submitted by the user
comprise the number of posterior draws (draws) and burn-ins
(burnin, 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
draws=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.
Expert users might want to take further adjustments. These have to be
provided via a list (expert). For example, to speed up
computation, it is possible to invoke parallel computing in
R. The number of available cpu cores can be specified via
cores. Ideally this number is equal to the number of units
\(N\)
(expert=list(cores=N)). Based on the user’s operating
system, the package then either uses parLapply (Windows
platform) or mclapply (non-Windows platform) to invoke
parallel computing. If cores=NULL, the unit models are
estimated subsequently in a loop (via R’s
lapply function). To use other / own apply functions, pass
them on via the argument applyfun. As another example, we
might be interested in inspecting the output of the \(N\) country models in more detail. To do
so, we could provide expert=list(save.country.store=TRUE),
which allows to save the whole posterior distribution of each unit /
country model. 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 this case, the whole
posterior distribution of the global model is stored.
We 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] 0.9708296 0.9929929 0.9636729 0.9788546 0.9914594 0.9959264 0.9927437
## [8] 0.9905900 0.9920747 0.9683688
Last, 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 summarize the outcome in various ways.
First, we can use theprint method
## ---------------------------------------------------------------------------
## Model Info:
## Prior: Normal-Gamma prior (NG)
## Number of lags for endogenous variables: 1
## Number of lags for weakly exogenous variables: 1
## Number of posterior draws: 100/1=100
## Size of GVAR object: 0.5 Mb
## Trimming leads to 36 (36%) stable draws out of 100 total draws.
## ---------------------------------------------------------------------------
## Model specification:
##
## EA: y, Dp, rer, stir, ltir, tb, y*, Dp*, rer*, stir*, ltir*, tb*, poil**, trend
## US: y, Dp, rer, stir, ltir, tb, poil, y*, Dp*, rer*, stir*, ltir*, tb*, trend
## RU: 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
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: Normal-Gamma prior (NG)
## Number of lags for endogenous variables: 1
## Number of lags for weakly exogenous variables: 1
## Number of posterior draws: 100/1=100
## Number of stable posterior draws: 36
## Number of cross-sectional units: 3
## ---------------------------------------------------------------------------
## Convergence diagnostics
## Geweke statistic:
## 90 out of 360 variables' z-values exceed the 1.96 threshold (25%).
## ---------------------------------------------------------------------------
## F-test, first order serial autocorrelation of cross-unit residuals
## Summary statistics:
## ========= ========== ======
## \ # p-values in %
## ========= ========== ======
## >0.1 6 33.33%
## 0.05-0.1 0 0%
## 0.01-0.05 4 22.22%
## <0.01 8 44.44%
## ========= ========== ======
## ---------------------------------------------------------------------------
## Average pairwise cross-unit correlation of unit-model residuals
## Summary statistics:
## ======= ======== ======== ========== ========== ======== ========
## \ y Dp rer stir ltir tb
## ======= ======== ======== ========== ========== ======== ========
## <0.1 0 (0%) 3 (100%) 2 (66.67%) 1 (33.33%) 2 (100%) 3 (100%)
## 0.1-0.2 0 (0%) 0 (0%) 0 (0%) 2 (66.67%) 0 (0%) 0 (0%)
## 0.2-0.5 3 (100%) 0 (0%) 1 (33.33%) 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 output 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 fraction 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).
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,
draws=100,
burnin=100,
plag=1,
prior="SSVS",
hyperpara=NULL,
SV=TRUE,
thin=1,
Ex=NULL,
trend=TRUE,
expert=list(save.shrink.store=TRUE),
hold.out=0,
eigen=1,
verbose=TRUE
)
One feature of the SSVS prior is that it allows us to look at the
posterior inclusion probabilities to gauge the importance of particular
variables. Per default, bgvar does not save the
volatilities of the coefficients to save memory. If we set
expert=list(save.shrink.store=TRUE) to TRUE
(default is FALSE) then those probabilities are saved and
posterior inclusion probabilities (PIPs) are computed. 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.53 0.32 0.10 0.10 0.46
## Dp_lag1 0.30 0.19 0.17 0.14 0.31 0.26
## rer_lag1 0.56 0.26 1.00 1.00 0.06 0.53
## stir_lag1 1.00 0.42 0.29 1.00 0.12 1.00
## ltir_lag1 1.00 0.67 0.08 0.21 1.00 1.00
## tb_lag1 0.05 0.67 0.34 0.20 0.15 1.00
## y* 0.24 0.48 0.14 0.09 0.07 0.69
## Dp* 0.06 0.72 0.33 0.26 0.23 0.37
## rer* 0.10 0.56 1.00 0.15 0.19 0.32
## stir* 0.20 0.80 0.09 0.30 0.12 0.09
## ltir* 0.22 0.34 0.35 1.00 0.09 0.28
## tb* 0.25 0.68 0.31 0.37 0.10 0.27
## poil** 0.46 1.00 0.09 0.14 0.17 0.71
## y*_lag1 0.84 0.22 0.06 0.03 0.04 0.26
## Dp*_lag1 0.68 0.49 0.24 0.09 0.36 0.20
## rer*_lag1 0.46 0.53 1.00 0.04 0.25 0.28
## stir*_lag1 0.23 0.26 0.04 0.41 0.08 0.27
## ltir*_lag1 1.00 0.26 0.25 1.00 0.80 0.39
## tb*_lag1 0.17 0.50 0.97 0.69 0.04 1.00
## poil**_lag1 0.44 0.99 0.14 0.17 0.14 0.61
## cons 0.25 0.32 1.00 0.27 0.07 1.00
## trend 1.00 0.59 0.23 0.04 0.26 1.00
The equations in the EA country model can be read column-wise with the
rows representing the associated explanatory 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.3900000 0.26000000 0.4066667 0.120 0.4000000 0.02
## Dp_lag1 0.3966667 0.4700000 0.19333333 0.1266667 0.425 0.2333333 0.22
## rer_lag1 0.3033333 0.4766667 1.00000000 0.4966667 0.245 0.5400000 0.97
## stir_lag1 0.8666667 0.2000000 0.21333333 1.0000000 0.290 0.4300000 0.19
## ltir_lag1 0.6150000 0.3500000 0.14000000 0.2100000 1.000 0.6250000 0.29
## tb_lag1 0.2733333 0.6133333 0.44000000 0.2166667 0.110 1.0000000 0.53
## y* 0.3566667 0.2266667 0.38666667 0.2500000 0.095 0.3866667 0.92
## Dp* 0.2100000 0.6466667 0.22000000 0.1600000 0.310 0.7866667 0.90
## rer* 0.3333333 0.3566667 0.85666667 0.1633333 0.570 0.2266667 0.17
## stir* 0.2333333 0.3633333 0.06666667 0.1700000 0.510 0.1533333 0.47
## ltir* 0.1300000 0.5000000 0.19333333 0.4900000 0.150 0.2000000 0.07
## tb* 0.4600000 0.5900000 0.25333333 0.2200000 0.160 0.3400000 0.83
## poil** 0.7300000 0.5300000 0.16500000 0.1500000 0.170 0.8550000 NaN
## y*_lag1 0.3900000 0.1800000 0.37000000 0.3733333 0.050 0.2000000 0.91
## Dp*_lag1 0.3700000 0.2366667 0.32333333 0.1166667 0.260 0.3900000 0.08
## rer*_lag1 0.4366667 0.3300000 0.80666667 0.0900000 0.605 0.2033333 0.19
## stir*_lag1 0.4766667 0.2066667 0.08000000 0.4866667 0.390 0.3100000 0.11
## ltir*_lag1 0.4900000 0.1966667 0.19666667 0.4200000 0.540 0.2266667 0.30
## tb*_lag1 0.2366667 0.3633333 0.48666667 0.3300000 0.195 0.4833333 0.28
## poil**_lag1 0.2600000 0.5400000 0.29000000 0.1600000 0.140 0.8050000 NaN
## cons 0.2733333 0.1800000 0.73000000 0.3333333 0.190 0.4500000 0.36
## trend 0.3800000 0.2433333 0.15333333 0.4033333 0.270 0.4000000 0.81
## poil_lag1 0.1300000 0.1100000 0.09000000 1.0000000 0.570 0.6000000 1.00This shows that the same determinants for the real exchange rate appear 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:
variable.list<-list();variable.list$real<-c("y","Dp","tb");variable.list$fin<-c("stir","ltir","rer")We can then re-estimate the model and hand over the
variable.list via the argument expert:
# weights for first variable set tradeW.0012, for second finW0711
model.ssvs.2<-bgvar(Data=eerData,
W=W.list[c("tradeW.0012","finW0711")],
plag=1,
draws=100,
burnin=100,
prior="SSVS",
SV=TRUE,
eigen=1,
expert=list(variable.list=variable.list,save.shrink.store=TRUE),
trend=TRUE
)
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 an exclusion
restriction on foreign long-term interest rates using
Wex.restr which is again handed over via
expert. The following includes foreign long-term rates only
in those country models where no domestic long-term rates are available:
# does include ltir* only when ltir is missing domestically
model.ssvs.3<-bgvar(Data=eerData,
W=W.trade0012,
plag=1,
draws=100,
burnin=100,
prior="SSVS",
SV=TRUE,
eigen=1,
expert=list(Wex.restr="ltir",save.shrink.store=TRUE),
trend=TRUE,
)## ---------------------------------------------------------------------------
## Model Info:
## Prior: Stochastic Search Variable Selection prior (SSVS)
## Number of lags for endogenous variables: 1
## Number of lags for weakly exogenous variables: 1
## Number of posterior draws: 100/1=100
## Size of GVAR object: 0.4 Mb
## Trimming leads to 18 (18%) stable draws out of 100 total draws.
## ---------------------------------------------------------------------------
## Model specification:
##
## EA: y, Dp, rer, stir, ltir, tb, y*, Dp*, rer*, stir*, tb*, poil**, trend
## US: y, Dp, rer, stir, ltir, tb, poil, y*, Dp*, rer*, stir*, tb*, trend
## RU: y, Dp, rer, stir, tb, y*, Dp*, rer*, stir*, tb*, poil**, trendLast, 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 (Kamiar Mohaddes and Raissi 2019).
The model structure would then look something like in the Figure below:
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<-eerData
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 OC model.
In the oil price example, we include the oil price (poil)
as an endogenous variable (not contained in any other country model) and
a weighted average using weights of output (y)
to proxy world demand as weakly exogenous variable. Next, we specify via
exo which one of the endogenous variables of the OC model
are fed back into the other country models. In this example we specify
poil. Last, we put all this information in a further list
called OE.weights (other entity weights). This is done to
allow for multiple other entity models (i.e., an oil price model, a
joint monetary union model, etc.). It is important that the list entry
has the same name as the other entity model, in our example
OC.
Now we can re-estimate the model where we pass on
OE.weights via the expert argument.
model.ssvs.4<-bgvar(Data=eerData2,
W=W.trade0012,
plag=1,
draws=100,
burnin=100,
prior="SSVS",
SV=TRUE,
expert=list(OE.weights=OE.weights,save.shrink.store=TRUE),
trend=TRUE
)##
## BGVAR incurred an error when estimating the models.
## See original error message:## Error occured in countrymodel: OC. Please check.## Error in detail:## Error in gamma_store[, , count] <- gamma :
## number of items to replace is not a multiple of replacement length## Error in bvar$Y: $ operator is invalid for atomic vectorsand 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]## Error in eval(expr, envir, enclos): object 'model.ssvs.4' not found## Error in eval(expr, envir, enclos): object 'aux4' not foundheatmap(aux1,Rowv=NA,Colv=NA, main="Model 1", cex.main=2, cex.axis=1.7)
heatmap(aux2,Rowv=NA,Colv=NA, main="Model 2", cex.main=2, cex.axis=1.7)
heatmap(aux3,Rowv=NA,Colv=NA, main="Model 3", cex.main=2, cex.axis=1.7)
heatmap(aux4,Rowv=NA,Colv=NA, main="Model 4", cex.main=2, cex.axis=1.7)## Error in heatmap(aux4, Rowv = NA, Colv = NA, main = "Model 4", cex.main = 2, : object 'aux4' not found
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.
