Bayesian estimation of large dimensional time varying VARs using copulas

Bayesian estimation of large dimensional time varying VARs using copulas

Introduction

Following the seminal contributions by sims80 and litterman1986, Vector AutoRegression (VAR) models and their variants are now ubiquitously applied to multivariate time series, as a flexible alternative to structural models. There is now a huge body of literature on both theory and applications: useful surveys are provided inter alia by watson94 and lutkepohl91. \newline

Although in their standard form VARs already offer a relatively flexible modelling approach, extensions have been considered to accommodate time variation. This may occur in the coefficients of the conditional mean equations (see e.g. doan; canova; sims93; stock96; and cogley), so affording a flexible alternative to models with abrupt breaks, known as the Time-Varying Parameter VAR (TVP-VAR). Time variation has also been considered in the covariance matrix of the error term, thereby allowing for time-varying heteroskedasticity.

Following seminal papers by uhlig, cogley and primiceri, recent examples where the assumption of homoskedasticity has been relaxed include koop13 and koop19, who attempt to reduce the dimensionality issue essentially by imposing a factor structure onto the volatilities - see also clark2011, carriero2015 , clark2015 and carriero2012. In a recent landmark paper, carriero2016large propose a far less restrictive set-up, which allows for fully Bayesian inference without imposing restrictions on the form of the heteroskedasticity.

Across the extensive literature on multivariate models, virtually all studies have one issue in common: the dimensionality of the model and the computational burden it brings about. On the one hand, unless the number of variables involved in the model is relatively large, omitted variable bias may impair the forecasting ability of the model (see giannone2006). carriero2016large make a compelling case for the superior performance of large dimensional VARs. On the other hand, computational difficulties may arise when there are a large number of variables and, more substantially, over-parameterisation can occur. Thus, the literature has focused on finding techniques that allow for the estimation of large VARs: see banbura for an excellent review of the various approaches which have been proposed.

In the case of homoskedastic VARs, the dimensionality issue can be handled by the choice of appropriate (conjugate) prior distributions, as shown by banbura who successfully apply their technique to the estimation of a VAR with 130 variables. Conversely, in the case of a heteroskedastic VAR, this is no longer possible, and the computational burden cannot be resolved through the choice of an appropriate prior. As explained in carriero2016large, heteroskedasticity invalidates the so-called \textquotedblleft symmetry\textquotedblright across equations that characterises homoskedastic VARs. Such property entails that a homoskedastic VAR is a SUR model where the regressors are the same across all equations; in turn, this entails that the covariance matrix of the VAR coefficients has a Kronecker structure, which makes estimation much simpler than if one had to deal with large matrices that do not have a simplifying structure. Few contributions consider estimation of a large VAR with time variation in both the coefficients of the conditional mean equations and of the covariance matrix of the error term. koop13 and koop19 propose an estimation technique for large, possibly heteroskedastic TVP-VARs, which essentially relies on the Kalman filter. However, their approach is not fully Bayesian, and it is liable to mis-specification issues if the assumed model for coefficient variation is not correct, also being, in practice, limited in dealing with the dimensionality issues - we refer however to a recent contribution by venditti which, through a non-parametric approach, solves these issues. carriero2016large solve the problem of fully Bayesian estimation of large VARs with heterskedasticity, by proposing a new estimation algorithm which is shown to perform very well in out-of-sample forecasting and also in structural analysis. However, their paper does not consider the presence of time varying coefficients in the VAR specification.

Proposed methodology and main contribution of this paper

This paper proposes a copula-based Bayesian estimation methodology for large TVP-VARs with heteroskedasticity. Similarly to carriero2016large, our estimators are fully Bayesian, thus allowing for the computation of the uncertainty around all estimators. \newline Full details of our approach are in Section (ref); here, we give a heuristic preview of our methodology. Given a multivariate model of (possibly) very large dimension $n$, we reduce it into $n$ univariate models, which are more easily handled. In order to recover the cross-dependence among equations, we use a copula-like term. In consequence, the likelihood function in our system factors into the likelihoods of the individual autoregressive models, plus the likelihood of the copula term. Thereby we are able to obtain a posterior where each set of parameters (the $n$ sets corresponding to the $n$ equations, and the set corresponding to the copula) is, conditional on the sample, independently distributed of the other sets. Therefore, from a computational point of view, each univariate problem is dealt with separately, which greatly reduces the computational burden. In this respect, our idea of breaking down the multivariate estimation problem into separate univariate problems is similar to the approach for fixed parameter, homoskedastic VARs developed in korobilis2019adaptive, although in our case we allow for time variation in the conditional mean and variance. The use of copulas to model dependence has also been considered by the literature in a Bayesian context (see e.g. gruber2015 and gruber2018), including non-parametric Bayesian analysis (we refer, inter alia, to the contributions by rodriguez2008bayesian, taddy2010autoregressive, nieto2012time, di2013simple, bassetti2014beta and nieto2016bayesian). \\ Our approach allows for great flexibility in the specification of the univariate models. For example, in the simplest version of our methodology, each series is modelled as an AR(1) model. However, given that more sophisticated model selection tools may be desirable to construct the univariate models, we develop an approach based on a model selection technique known as Bayesian compression (see dunson). Moreover, whilst the focus of this paper in on TVP-VARs with heteroskedasticity, our approach can be also used to estimate other multivariate models such as e.g. VARMAs and Multivariate Stochastic Volatility models (MSV). In particular, in appendix we carry out an empirical exercise using VARMAs, to illustrate the computational convenience of our approach. Also, in another contribution (itt2019), we apply our methodology to large MSV models for financial variables. Here, we illustrate our methodology by estimating a large TVP-VAR model with possible heteroskedasticity, using the same data as koop13.

The remainder of the paper is organised as follows. Our methodology is spelt out in Section (ref). The empirical application is in Section (ref); we also report a further application to VARMAs in Appendix (ref). Section (ref) concludes.

Methodology

We begin by introducing the main model and some notation. We consider the TVP-VAR($p$)

where $\mathbf{y}_{t}$ is an $n\times 1$ vector and $\mathbf{u}_{t}$ is a zero mean, Gaussian process with possibly time varying variance - we discuss the specification of the second moment later on. Model ( (ref)) could be extended to consider e.g. exogenous regressors, latent regressors such as common factors, or deterministics such as a constant, (linear or nonlinear) trends and seasonal dummies. Also, ((ref) ) could also have an MA($q$) structure, in the spirit of chan2013; or it could have no autoregressive structure at all, and only time varying heteroskedasticity as in the case of creal. We prefer to focus on a simpler specification, so that the discussion is not overshadowed by the algebra. Similarly, the assumption that $\mathbf{u}_{t}$ is Gaussian is made only for simplicity. Note that, even in this simple set-up, the number of parameters grows rapidly with $p$ and $n$, whence the dimensionality issue.

Theory: the univariate equations, the copula and the likelihood function

The univariate equations

In the context of ((ref)), we consider the following univariate AR($p$) models

for $1\leq i\leq n$ with $u_{i,t}=e^{h_{i,t}/2}u_{i,t}^{*}$, where $ u_{i,t}^{\ast }\sim i.i.d.N(0,1)$ and

with $e_{i,t}\sim i.i.d.N(0,\delta _{i})$. As noted above, ((ref)) can be extended and/or modified to incorporate e.g. a different number of lags $p_{i}$ for each unit, an MA($q_{i}$) component, exogenous regressors, deterministics, (conditional or unconditional) heteroskedasticity, etc.. Similarly, ((ref)) could be replaced e.g. by a GARCH specification to allow for conditional heteroskedasticity (see also the discussion in uhlig on the relative merits of possible specifications for time heteroskedasticity).

Given that in ((ref)) $y_{i,t}$ is predicted using only its own past, this may lead to a loss of predictive accuracy. A possible alternative would be to use the Bayesian compression algorithm developed in dunson. In particular, we consider the specification

with $u_{i,t}$ still satisfying ((ref)). As in ((ref)), $ z_{i,t}^{\left( 2\right) }$ is a subset of the regressors in each equation of the unrestricted VAR (say $\widetilde{z}_{i,t}$). However, in the case of ((ref)), the vector $z_{i,t}^{\left( 2\right) }$ can include lags of $ y_{i,t}$ and also lags of $y_{j,t}$ for $j\neq i$. In order to select the components of $z_{i,t}^{\left( 2\right) }$, dunson suggest the following technique. Let $z_{i,t}^{\left( 2\right) }=\Phi \widetilde{z}_{i,t} $ with $\Phi $ a $p\times np$ matrix whose entries are defined as

and $\phi $ and $p$ are drawn uniformly from $\left( 0.1,1\right) $ and $ \left\{ 1,...,p^{\max }\right\} $, with $p^{\max }$ chosen so that the marginal likelihood has a global peak. The matrix $\Phi $ is then normalised via the Gram-Schmidt orthonormalisation - see dunson for details.

The copula

We now introduce the copula term to model dependence among the univariate equations. Letting $X$ denote a continuous $k$-dimensional random variable whose density is given by $f\left( x\right) $, it holds that

where $f_{j}\left( x_{j}\right) $ is the density of the $j$-th coordinate of $X$, $v_{j}=F_{j}\left( x_{j}\right) =\int_{-\infty }^{x_{j}}f_{j}\left( u\right) du$, and $c^{\ast }\left( v_{1},...,v_{k}\right) $ is the copula density (which is unique since $X$ is continuous). This result is known as Sklar's theorem (see sklar1 and sklar2; see also the book by nelsen for a comprehensive introduction to copulas). Equation ((ref)) can equivalently be written as

The likelihood function

We now turn to specifying the likelihood. Henceforth, we use $z_{i,t}$ as short-hand for both $z_{i,t}^{\left( 1\right) }$ and $z_{i,t}^{\left( 2\right) }$; $\beta _{i,t}$ for both $\beta _{i,t}^{\left( 1\right) }$ and $ \beta _{i,t}^{\left( 2\right) }$ in ((ref)) and ((ref)) respectively. We assume the following law of motion

with $\epsilon _{i,t}\sim i.i.d.N\left( 0,\Sigma _{i}\right) $, independent across $i$. We point out that, in ((ref)), we do not impose the typical random walk model for the time-varying parameters (see e.g. koop13), which makes our set-up more general. For simplicity, we do not allow for time variation in any other parameter (i.e., we do not allow for the copula parameters, or the coefficients in ((ref)), to vary over time).

Let $b_{i}=\left( \alpha_{i},\gamma _{i},\delta _{i}\right) $. Then, the marginal density of $y_{i,t}$ conditional on $z_{i,t}$ can be denoted as $f_{i}\left( y_{i,t}|z_{i,t};\beta _{i,t},b_{i}\right) $ (we omit dependence on $A_{\beta ,i}$, $\Sigma _{i}$ and $\beta _{i,0}$ for short). Then, by ((ref)), it holds that

having defined $v_{i,t}=\int_{-\infty }^{y_{i,t}}f_{i}\left( u|z_{i,t};\beta _{i,t},b_{i}\right) du$, with $f_{i}\left( u |z_{i,t};\beta _{i,t},b_{i}\right) $ denoting the density of $y_{i,t}$ conditional on $ z_{i,t}$.

Although Sklar's theorem ensures that the copula density $c^{\ast }\left( v_{1,t},...,v_{n,t}\right) $ is unique, it does not provide an expression for it. One possibility would be to consider a non-parametric copula, and we refer to scaillet2002, ibragimov and chen2007nonparametric, and the references therein, for the relevant theory in a time series context. In this paper, we choose a different set-up. In particular, we consider a (parametric) Gaussian mixture copula model (GMCM henceforth; see tewari) model, viz.

where: $\left\{ p_{g}\right\} _{g=1}^{G}$ (such that $\sum_{g=1}^{G}p_{g}=1$ and $p_{1}<...<p_{G}$) is a set of weights and $f_{N}\left( \cdot |\mu _{g},\Omega _{g}\right) $ is the density of an $n$-variate Gaussian random variable with mean $\mu _{g}$ and covariance matrix $\Omega _{g}$. In ((ref)), we use the short-hand notation $\alpha =\left( \left( p_{1},...,p_{G}\right) ^{\prime },\mu _{1}^{\prime },...,\mu _{G}^{\prime },\left( vech\left( \Omega _{1}\right) \right) ^{\prime },...,\left( vech\left( \Omega _{G}\right) \right) ^{\prime }\right) ^{\prime }$.

Finally, let now $\beta _{t}=\left( \beta _{1,t}^{\prime },...,\beta _{n,t}^{\prime }\right) ^{\prime }$, $\omega =\left( b_{1}^{\prime },...,b_{n}^{\prime }\right) ^{\prime }$, $A_{\beta }=\left\{ A_{\beta ,1},...,A_{\beta ,n}\right\} $, $\Sigma =\left\{ \Sigma _{1},...,\Sigma _{n}\right\} $, and, for short,

Putting everything together, the resulting likelihood function (conditional on the initial observations $\left\{ \mathbf{y}_{t}\right\} _{t=1}^{p}$) is given by

where we have now emphasized the dependence of the marginal densities on $ A_{\beta ,i}$, $\Sigma _{i}$ and $\beta _{i,0}$.

It follows that

which indicates that maximisation with respect to each unit $i$ can be carried out separately, like maximisation with respect to $\alpha $.

Dimension reduction and estimation

Equation ((ref)) indicates that the likelihood can be factored into $n+1$ independent problems. We choose the prior

where: $p\left( \alpha \right) $ and $p\left( \omega _{i,0}\right) $ are flat priors (in the latter, coefficients are restricted to be non-negative); $p\left( \Sigma _{i}\right) \propto \left\vert \Sigma _{i}\right\vert ^{-\left( n+1\right) /2}$ as a standard non-informative prior; finally, $p\left( A_{\beta ,i}\right) $ and $p\left( \beta _{i,0}\right) $ are Gaussian priors and we discuss them in details in Section (ref). Thus, by construction, $p\left( \theta \right) $ can also be factorised into $n+1$ independent problems.

Hence, the posterior is given by

Note that, based on ((ref)) and ((ref)), the posterior again factorises into separate posteriors for each unit-specific set of parameter. This entails, as discussed in the introduction, that the estimation of the TVP-VAR with possible heteroskedasticity\ can be decomposed into $n+1$ estimation problems that can be carried out in parallel, independently of each other. We point out that this result holds as long as the prior on $\alpha $, $p\left( \alpha \right) $, is independent of the other parameters; conversely, the prior on the other parameters can have a hierarchical structure, so that ((ref)) might alternatively be written as

Then, by standard arguments, ((ref)) would become

Prior to discussing estimation, some considerations on the potential for dimensionality reduction are in order. Despite the presence of the copula term, the number of parameters in $\theta $ is still proportional to $n^{2}$ , which does not fully resolve the challenge represented by dimensionality in a large VAR. More specifically, from ((ref)), it is apparent that, when estimating $\mu _{g}$, the number of parameters to be estimated is $Gn$; conversely, the matrices $\Omega _{g}$ contain each $ \frac{n\left( n+1\right) }{2}$ elements and this is where the dimensionality issue arises from. In order to attenuate this problem, in Section (ref) we consider two ways of restricting the $\Omega _{g}$s, which both reduce the number of free parameters in the copula to being proportional to $n$ as opposed to $n^{2}$.

Univariate regressions estimation - $\beta_{i,t}$ and $b_{i}$

Each equation ((ref)) and ((ref)) is a regression (or, if specification ((ref)) is indeed chosen, an autoregression) with time varying parameters and stochastic volatility. Thus, we use the approach by kim1998 to estimate $\beta_{i,t}$ and $b_{i}$ (and the other hyperparameters).\\ More precisely, note that ((ref)) entails that

Thus, conditional on $\beta_{i,t}$s, we have

The model is linear in $h_{i,t}$. It is well known (see kim1998) that using a Quasi-Maximum Likelihood estimator under the assumption that $\ln u_{i,t}^{*2}$ is normal results in poor small-sample properties; thus, we follow the approach suggested by kim1998. In particular, we approximate the distribution of $\ln u_{i,t}^{*2}$ using a mixture of normals with seven components. Thence, for each $i$, $h_{i,t}$s is sampled at once using the Kalman filter. In turn, conditionally on $h_{i,t}$, the model for $y_{i,t}$ has a linear state space representation in terms of $\beta_{i,t}$s. Therefore, for each $i$, we draw the entire vector $\beta_{i,t}$ at once, using again the Kalman filter.\footnote{We point out that an alternative approach is to use the Gibbs sampler to draw from the conditional posterior distribution $p\left(\beta_{i,t}|\{h_{i,\tau},\tau\neq t\},\{h_{i,t}\},\{\mathbf{y}\}_{t=1}^{T}\right) $ but this approach, although simpler, results in slower convergence and higher autocorrelation in MCMC draws.}

Copula density estimation - $\protect\alpha $

As is typical with copula models, we first obtain an estimate of the univariate densities $f_{i}\left( y_{i,t}|z_{i,t};A_{\beta ,i},\Sigma _{i},\beta _{i,0},b_{i},\beta _{i,t}\right) $. We then obtain the probability integral transforms, $v_{i,t}$, and use these as data to estimate $\alpha $.\footnote{ This procedure can be viewed as \textquotedblleft two-stage\textquotedblright Bayesian, as opposed to a \textquotedblleft full-information\textquotedblright Bayesian estimator (see also creal). Whilst this could, in principle, be carried out by modifying the MCMC algorithm, it adds to the computational complexity of the estimation; further, we have tried to use it in some of our empirical applications, but results were - if anything - marginally worse than with the proposed two-step procedure which we study here.}

The dimensionality issue can be further addressed by imposing some restrictions on $\alpha $. We discuss two possible approaches (denoted as S1 and S2), where the priors employed are flat.

Dimension reduction: strategy S1

The first dimension reduction strategy is based on a recursive model for the $\Omega _{g}$s:

having initialised ((ref)) by leaving $\Omega _{1}$ unrestricted (and thus setting $V_{1}=0$). In ((ref)), $h_{g}$ is a scalar to be estimated, and the idiosyncratic shock $V_{g}$ is restricted to be $V_{g}=diag\left\{ v_{g,1},...,v_{g,n}\right\} $.

Consequently, the number of parameters associated to the copula is $\left( G-1\right) \left( n+1\right) $, which therefore grows linearly, as opposed to quadratically, with $n$.

Dimension reduction: strategy S2

Another possible dimension reduction approach is intimately related to Principal Components (we refer to humphreys for a full treatment, which we briefly summarize here), and to the Bayesian compression literature (dunson). We again leave $\Omega _{1}$ unrestricted, and model the $\Omega _{g}$s, for $2\leq g\leq G$, as

In ((ref)), $D_{g}=diag\left\{ d_{g,1},...,d_{g,n}\right\} $ and $Q_{0,g}$ is an $n\times k$ matrix. \newline We make no attempt to estimate $Q_{0,g}$. Instead, we randomly generate the elements of $Q_{0,g}$, say $\left\{ Q_{0,g}\right\} _{i,j}$, $1\leq i\leq n$ , $1\leq j\leq k$, as independent of each other with $\left\{ Q_{0,g}\right\} _{i,j}\sim N\left( 0,q_{g}^{2}\right) $ for a total of $ 1,000,000$ iterations, choosing the specification which maximises the log marginal likelihood. Thus, the only parameters that need to be estimated are $q_{g}$ and $\left\{ d_{g,1},...,d_{g,n}\right\} $. Under these restrictions, the number of parameters is $\left( G-1\right) \left( n+1\right) $ - i.e., the same as in S1.

Sampling from the posterior: the MALA\ algorithm

Sampling from ((ref)) can be done along similar lines as in the case of a fixed coefficient VARs, but with the complications arising from $\beta _{t}$ being time-varying. We use the Metropolis Adjusted Langevin (MALA) algorithm by roberts1998 (see also girolami), which is likely to be more efficient than an ordinary Random Walk Metropolis algorithm in light of the large dimensionality of $\theta $.

In order to illustrate the algorithm, we begin by defining the matrix

computed at a generic value $\widetilde{\theta }$. The likelihood $L\left( \left\{ \mathbf{y}_{t}\right\} _{t=p+1}^{T}|\theta \right) $ is differentiable up to any order, within the whole parameter space, due to the normality assumption; thus, by the Schwarz Lemma, $G\left( \widetilde{\theta }\right) $ is symmetric for any $\widetilde{\theta }$ within the parameter space.

Based on the definitions above, the resampling scheme is as follows:

We now discuss the proposal density. In ((ref)), the scale parameter $\lambda $ is discussed later on, and the mean $m\left( \theta _{k}\right) $ is given by

where \textquotedblleft $\nabla $\textquotedblright\ refers to the gradient, which is computed with respect to $\theta $ and then specialised in the value $\theta _{k}$ (we use the same notation as nemeth). In ((ref)), the main difficulty is the computation of

Assuming - as is typical, see nemeth - that $\nabla \ln p\left( \theta _{k}\right) $ is known, this boils down to estimating $\nabla \ln L\left( \left\{ \mathbf{y}_{t}\right\} _{t=1}^{T}|\theta _{k}\right) $. Note that by Fisher's identity (cappe), it holds that

where $E_{\left\{ \beta _{t}\right\} _{t=1}^{T}}$ denotes expectation taken with respect to $p\left( \left\{ \beta _{t}\right\} _{t=1}^{T}|\left\{ \mathbf{y}_{t}\right\} _{t=1}^{T}\right) $, with

We carry out the estimation of $\nabla \ln L\left( \left\{ \mathbf{y} _{t}\right\} _{t=1}^{T}|\theta _{k}\right) $ using the Rao-Blackwellised estimator proposed in nemeth, as described below.

The output of the algorithm is the estimate $\nabla \ln \widehat{L}\left( \left\{ \mathbf{y}_{t}\right\} _{t=1}^{T}|\theta _{k}\right) $, which can then be plugged in ((ref)). As indicated by nemeth, the algorithm also readily affords the computation of other important quantities such as the predictive likelihood, etc.

Empirical application

In this section, we illustrate our methodology by applying it to the estimation of a (large) TVP-VAR with heteroskedasticity. We use the same data as koop13, namely $n=25$ US macroeconomic variables (see Table (ref)) running from 1959:Q1 to 2010:Q2. The focus of our exercise is the prediction of three series: inflation, GDP and interest rate. Given that all series are transformed into first differences in order to ensure stationarity, our model predicts the percentage change in inflation (the second log difference of CPI), GDP growth (the log difference of real GDP) and the change in the interest rate (the difference of the Fed funds rate). To ensure a fair comparison with koop13, we have also demeaned all variables and then standardised them (we use the standard deviation calculated from 1959Q1 through 1969Q4). The forecasting horizon is 1970:Q1 till 2010:Q2.

Results are in Tables (ref)-(ref), where we report the relative Mean Squared Forecast Errors (MSFE) when using the various VAR specifications to predict GDP, inflation and interest rate (respectively). The numbers in the tables are the MSFE relative to the TVP-VAR-DMA model in koop13, which is therefore our baseline model.

Broadly speaking, results show that our methodology affords good forecasting ability especially for shorter horizons; a notable exception is the poor performance of the TVP-VAR for GDP, although using strategies S1 and S2 yields a marked improvement. Indeed, there is no clearly superior model, although the results seem to make a case for heteroskedastic VARs (nonetheless, homoskedastic VARs with GMCM show very good results). In general, using GMCM (and determining the $G$) works better than restricting $ G$ to 1 (as could be expected). Similarly, reducing the dimensionality of the copula model with either strategy S1 or S2 generally improves forecasting ability. Although Bayesian compression works well, it does not seem to yield a uniformly superior predictive performance than the univariate models proposed in equation ((ref)). As a final point, a distinctive feature of koop13 is that the authors propose to use \textquotedblleft forgetting factors\textquotedblright (a procedure not dissimilar to an exponentially weighted moving average); thus, they avoid estimating the covariance matrix of the VAR and the covariance matrix of the time-varying coefficients. In our case, we are dealing with univariate models, and therefore we do not have to estimate covariance matrices.

Prior sensitivity

We have carried out a further exercise to explore the sensitivity of our methodology to the choice of the (main) priors on $A_{\beta,i}$ and $\beta_{i,0}$. We point out that - in this contribution - the main focus is not so much the choice of the prior but the copula-based dimensionality reduction. Indeed, we propose flat priors in general, although of course some parameters undergo nonlinear transformations which invalidates this argument (see the classical reference by jeffreys for an early treatment of the issue). Hence the importance of at least validating the choice of our priors through sensitivity analysis.

We begin by describing the benchmark prior. For each element in the vector $\left(vec\left( A_{\beta,i}\right)', \beta_{i,0}' \right)'$, we have chosen the prior $N(\overline{b},\overline{s}_{b}^{2})$, independent across elements. As far as the copula functions are concerned, recall ((ref)). We have used both dimension reduction strategies S1 and S2, with:

where

and

again independent for $1\leq g\leq G$. Finally, in ((ref)), we have used $ \Omega _{g}=C_{g}C_{g}^{\prime }$, with :

for $2\leq g\leq G$, where we have defined $c_{g}=vech\left( C_{g}\right) $. We have set the priors parameters as follows:

In our analysis, we have used $1,000$ different priors by sampling randomly from ((ref))-((ref)), given the parameters defined in ((ref)). For each prior, we have used MCMC sampling, employing $ 10,000$ iterations starting from the posterior moments delivered by the benchmark prior. Note that we have not examined sensitivity with respect to other priors, which are anyway rather diffuse.

In order to compare our results against the TVP-VAR-DMA model in koop13, we have computed the relative MSFEs as above. The sampling distributions of these relative MSFEs are reported in Figures (ref)-(ref). As can be seen, strategy S2 seems to deliver the best results, both in terms of the mode of the sampling distribution, and the dispersion around it.

Conclusions

Our paper has developed an alternative methodology for the direct estimation of large TVP-VARs with possible heteroskedasticity. The original multivariate model is decomposed into $n$ simpler models, whose interactions are modelled separately through a copula. We use the GMCM copula, whose good performance in our context is in line with the conclusions of other papers (see e.g. geweke2007 and villani). In principle, however, it would be possible to use also other copula specifications; given that considering this approaches goes beyond the scope of our paper (and the GMCM copula did not pose any particular runtime issues), this is an area for future research.

Our empirical applications (see also the estimation of VARMAs in Appendix (ref)) show that our approach is computationally more convenient than directly estimating multivariate models. In addition, our results also show excellent goodness of fit and predictive ability. We note that, when reducing the original multivariate model into $n$ separate models, it is not necessary to impose a pure AR(1) structure in which each series is predicted using solely its own lags. Indeed, we also consider a different model reduction strategy based on Bayesian compression. However, we found that even univariate, simple AR(1) models afford good forecasting ability. These considerations support the conclusion that the use of copulas, particularly in high dimension, is advantageous in that the copula manages to capture features of the data that the original, standard multivariate models are likely to miss. Thus, our contribution may also be viewed as a complement to the recent advances in the Bayesian analysis of large VARs, such as the ones developed in banbura and glp, where - instead of using copulas - new, more sophisticated priors are proposed as a way to deal with large VARs.

We point out that our applications mainly focus on \textquotedblleft reduced form\textquotedblright\ examples, as can be seen by the emphasis on forecasting ability. We conjecture however that, in light of its excellent performance, our technique could also be employed in the context of more structural applications. This issue is currently under examination by the authors.