Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
79,211 characters · 12 sections · 77 citation commands
A new algorithm for structural restrictions in Bayesian vector autoregressions\@thefnmark\@footnotetextI would like to thank without implicating Christiane Baumeister, Martin Bruns, Fabio Canova, Filippo Ferroni, Luca Gambetti, Toru Kitagawa, Gary Koop, Michele Lenza, Laura Liu, Emanuel M\"onch, Alberto Musso, Serena Ng, Michele Piffer, Davide Pettenuzzo, Francesco Ravazzolo, Frank Schorfheide, Maximilian Schr\"oder, Christian Schumacher, Leif Anders Thorsrud, John Tsoukalas and Harald Uhlig, as well as seminar and conference participants, for useful discussions and comments. Any remaining errors should solely be attributed to the author. MATLAB code that replicates the Monte Carlo and empirical results of this paper is available on https://sites.google.com/site/dimitriskorobilis/https://sites.google.com/site/dimitriskorobilis/. Correspondence: Professor of Econometrics, Adam Smith Business School, University of Glasgow, 40 University Avenue, Glasgow, G12 8QQ, UK; email: mailto:[email removed]@googlemail.com.
\thispagestyle{empty}
\doublespacing
\setcounter{page}{1}
This paper proposes a new Bayesian Markov chain Monte Carlo (MCMC) algorithm for joint estimation of parameters of reduced-form vector autoregressions (VARs) and associated sign restrictions for structural identification. The main idea is to allow the reduced-form VAR disturbances to have a static factor model structure. By doing so, sign and other restrictions can be incorporated via straightforward parametric prior distributions, and the factors can be interpreted as the structural VAR (SVAR) disturbances. A new, computationally efficient, algorithm is able to jointly sample VAR parameters and identification restrictions. The implication of this feature is that the parameter estimates and the fit of the VAR depend on, and interact with, the identification restrictions the researcher has in mind. Existing reduced-form VAR approaches typically follow a two-step procedure in which an estimate of the VAR covariance matrix is obtained in the first step, and some identification scheme that seems plausible to the researcher is imposed in a second step.\footnote{Consider a VAR covariance matrix estimate $\mathbf{\widehat{\Omega}}$, structural identification simply boils down to finding a matrix $\mathbf{A}$ such that $\mathbf{A}\mathbf{A}^{\prime} = \mathbf{\widehat{\Omega}}$. There are infinite such matrices that satisfy this relationship, therefore, it is required to impose some zero or other restrictions on $\mathbf{A}$. However, these restrictions can never be (statistically) tested since $\mathbf{\widehat{\Omega}}$ is fixed and the data likelihood remains unchanged no matter what the researcher thinks restrictions in $\mathbf{A}$ should be.} In the proposed modeling approach, different identification schemes result to different VAR parameter estimates and general model fit. The benefit of the new approach is that the researcher can treat parameter estimation and identification as a joint problem, which is a very advantageous approach towards inference, due to the fact that either the “true” VAR parameters or the “true” structural restrictions are never known. By extracting unobserved factors from VAR disturbances, sign and zero structural restrictions become, respectively, inequality and zero parametric restrictions in the associated factor loadings matrix.
The sign restrictions approach to identification has become very popular in applied work compared to traditional identification methods such as exclusion restrictions; see kilianlutkepohl2017 for a detailed review of this literature. The main feature of popular Bayesian algorithms for inference in sign restrictions, such as RubioRamirezetal2010, is that they rely on rejection sampling schemes (also known as accept/reject algorithms) in order to search for matrices that satisfy the desired restrictions. If restrictions are tight, as it would be the case in models with many variables and many shocks, rejection sampling results in constantly rejecting draws. In contrast, the Gibbs sampler proposed in this paper allows to sample contemporaneous structural matrices from their conditional posterior and these samples are always accepted. The benefits of the new approach are demonstrated by revisiting the empirical results in Furlanettoetal2017 by using a single 15-variable VAR with many shocks, instead of many 5-variable VARs for identifying a few shocks at a time (which is what these authors do due to the computational constraints of the RubioRamirezetal2010 algorithm they adopt). As a rough indication of the computational efficiency of the new algorithm, I find that obtaining 5,000 uncorrelated samples from the benchmark six-variable VAR of Furlanettoetal2017 using the new algorithm takes less than five minutes; using the original RubioRamirezetal2010 algorithm that Furlanettoetal2017 adopted in order to produce their results, it takes roughly four hours to obtain 2,000 draws that satisfy the same restrictions. Empirically, using the same set of sign restrictions, the two algorithms produce comparable results as measured by shapes and magnitudes of impulse response functions.
The new algorithm for structural inference in VARs shares some similarities, from a computational perspective, with the SVAR approach of BaumeisterHamilton2015. These authors estimate a joint model for structural restrictions and parameter estimation. However, the need to derive a reasonably simple algorithm algorithm for inference, means that these authors integrate out autoregressive and variance parameters from the joint posterior of parameters and structural restrictions using natural conjugate priors. The result is a Metropolis-Hastings algorithm that is also of the accept-reject form and cannot scale up easily to very high dimensions.\footnote{To be exact, if a candidate sample is not accepted as a sample from the true posterior, then the immediately previous accepted sample values are used. RePEc:diw:diwwpp:dp1796 propose a more efficient Dynamic Striated Metropolis Hastings algorithm that builds on importance sampling proposals. While this algorithm is appropriate for high-dimensional models, in the case of the VAR it will still hit a computational bottleneck at much lower dimensions than the Gibbs sampler proposed in this paper.} Most importantly, the adoption of a natural conjugate prior means that it is not possible to extend the BaumeisterHamilton2015 methods with empirically relevant time series features, such as heteroskedasticity or structural breaks. In contrast, the algorithm proposed here builds on a standard reduced-form Bayesian VAR that is easy to work with and extend to large dimensions, nonlinear or asymmetric shocks, nonlinear parameters (e.g. stochastic volatility), and numerous other interesting features; see KoopKorobilis2010 for a thorough review of Bayesian tools for inference in reduced-form VARs.
The idea to decompose the VAR disturbances into common factors is also related to numerous other modelling approaches. Gorodnichenko2005 specified an identical VAR model with reduced-rank decomposition of the disturbance term. The purpose of specifying the VAR that way was to replace standard block diagonal restrictions in VARs BernankeMihov1992 with a more parsimonious identification scheme that imposes less (possibly unreasonable) zero restrictions. More recently, MatthesSchwartzman2019 specify a closely related VAR model in order to identify the structural impact of sectoral dynamics on GDP. Their identification is via a factor structure on the residuals that has the additional assumption of allowing for correlation within industries but no correlation across industries.
In a different strand of the VAR literature, StockWatson2005 specify a more general factor-augmented VAR (FAVAR) and discuss in detail how various identification schemes fit in this setting. They note StockWatson2005 that the sign restrictions identification scheme proposed by Uhlig2005 also fits the FAVAR framework. An application of this idea can be found in AhmadiUhlig2016. From a modeling point of view, the model I propose in this paper can be viewed as a special case of the AhmadiUhlig2016 FAVAR model. However, the specification I use has completely different implications both algorithmically and in terms of inference. AhmadiUhlig2016 project a large vector of observable macroeconomic variables into a smaller vector of factors and they model VAR dynamics only on these factors. This means that there is some loss of information (not all macro variables are explained well by the factors) and the statistical fit of the factors determines the contribution of each structural shock on each macroeconomic variable. Additionally, the autoregressive dynamics of the large macro dataset is represented only by the autoregressive dynamics of the smaller vector of factors. This modeling choice means that, inevitably, the FAVAR is unable to capture richer patterns of propagation of structural shocks to observed macroeconomic variables. In contrast, in this paper all observable macroeconomic variables are endogenous in the VAR and the sole role of the factors is to represent structural shocks. Additionally, the algorithm derived here is computationally simpler as it relies on posterior formulas for linear regression models, instead of building on more demanding simulation smoothing techniques, as is the case with the FAVAR (see AhmadiUhlig2016 and Bernankeetal2005).
The next Section introduces the new methodology and associated Gibbs sampler algorithm for inference, and it outlines the key components that help speed up and stabilize (numerically) posterior sampling in high dimensions. In Sections 3 I undertake several important exercises using synthetic datasets, in order to test both the computational features of the new algorithm as well as shed light on how joint inference on parameters and structural restriction is implemented. In Section 4 the algorithm is applied to the issue of measuring the impact of a financial shock to the macroeconomy. Section 5 concludes the paper.
The starting point is the reduced-form vector autoregression
where $\mathbf{y}_{t}$ is a $\left( n \times 1 \right)$ vector of observed variables, $\mathbf{x}_{t} = \left( 1,\mathbf{y}_{t-1}^{\prime},...,\mathbf{y}_{t-p}^{\prime} \right)^{\prime}$ a $\left( k \times 1 \right)$ vector (with $k=np+1$) containing a constant and $p$ lags of $\mathbf{y}$, $\mathbf{\Phi}$ is an $(n \times k)$ matrix of coefficients, and $\bm{\varepsilon}_{t}$ a $\left( n \times 1 \right)$ vector of disturbances distributed as $N\left( \mathbf{0}_{n \times 1},\mathbf{\Omega} \right)$ with $\mathbf{\Omega}$ an $n \times n$ covariance matrix. The structural VAR (SVAR) form associated with the reduced-form model in (ref) is
where $\mathbf{B} = \mathbf{A}\mathbf{\Phi}$, $\bm{u}_{t} = \mathbf{A}\bm{\varepsilon}_{t}$ and $cov(\bm{u}_{t}) = \mathbf{D}$, with $\mathbf{D}$ an $n \times n$ diagonal matrix which, sometimes, is normalized to be the identity matrix. The SVAR form can be obtained by means of a decomposition of the reduced-form covariance matrix of the form $\mathbf{A}\mathbf{\Omega}\mathbf{A}^{\prime} = \mathbf{D} $ where both sides of equation (ref) are left-multiplied with the $n \times n$ matrix $\mathbf{A}$. This decomposition is unique when $\mathbf{A}$ is a lower triangular matrix (typically with a unit diagonal, unless we do the normalization $\mathbf{D}=\mathbf{I}$), but it has infinite solutions for a full matrix $\mathbf{A}$.
I begin by building on fundamental ideas introduced in the factor model literature, as applied to empirical problems in macroeconomics: a few common forces (which in a structural setting we desire to identify as “primitive shocks”; see Ramey2016) are driving the set of reduced-form shocks in the system of $n$ endogenous variables. In order to materialize this idea, the reduced-form VAR disturbances in equation (ref) are decomposed using the following static factor model specification
where $\mathbf{\Lambda}$ is an $n \times r$ matrix of factor loadings, $\mathbf{f}_{t}$ is an $r \times 1$ vector of factors, and $\mathbf{v}_{t}$ is an $n \times 1$ vector of idiosyncratic shocks. While the $n$ shocks in $\bm{\varepsilon}_{t}$ are decomposed into $r + n$ shocks, only the $r$ common shocks in $\mathbf{f}_{t}$ are considered structural and the $n$ shocks in $\mathbf{v}_{t}$ are simply nuisance shocks (e.g. due to measurement or expectations error, incomplete information etc). The assumption here is that $n$ is large and that $r \leq n$, and not necessarily $r \ll n$, as is typically assumed in the factor literature. In line with the exact factor model literature, let $\mathbf{v}_{t} \overset{\text{i.i.d}}\sim N \left(\mathbf{0}_{n \times 1}, \mathbf{\Sigma} \right)$, with $\mathbf{\Sigma}$ an $n \times n$ diagonal matrix. Additionally, let $\mathbf{f}_{t} \sim N \left( \mathbf{0}_{r \times 1}, \mathbf{I}_{r}\right)$, such that the conditional covariance matrix of $\bm{\varepsilon}_{t}$ is now of the form
This factor model decomposition of $\mathbf{\Omega}$ reveals that, as long as $\mathbf{\Sigma}$ is diagonal, identification via sign restrictions can be achieved by imposing the desired signs on $\mathbf{\Lambda}$. Similarly, zero restrictions simply correspond to setting the respective elements of $\mathbf{\Lambda}$ to zero.\footnote{If desired, several other restrictions can be incorporated in a straightforward way, such as ranking restrictions https://doi.org/10.3982/QE1277, restrictions on elasticities, or other restrictions on magnitudes of shocks. For example, if for shock $j$ variable $i$ should react with a larger magnitude than variable $k$, then we get the restriction $\mathbf{\Lambda}_{ij} > \mathbf{\Lambda}_{kj}$ which is fairly simple to incorporate within an MCMC sampling setting. The key early reference for inference in the Bayesian regression model with general inequality constraints is Geweke1996.} To see this, consider a reduced-rank SVAR representation of this model, which can be obtained by left-multiplying the reduced-form VAR model given by equations (ref) - (ref) with the generalized inverse of $\mathbf{\Lambda}$, as follows:
In the equation above, the SVAR matrix $\mathbf{A}_{1}$ is equivalent to the generalized inverse $\left( \mathbf{\Lambda}^{\prime} \mathbf{\Lambda} \right)^{-1}\mathbf{\Lambda}^{\prime}$. While $\mathbf{\Lambda}$ is not observed, assume that a consistent estimator of this parameter exists. Given that in the exact factor model formulation the $\mathbf{v}_{t}$ are uncorrelated, the Central Limit Theorem in Bai2003 suggests that for each $t$ and for $n\rightarrow \infty$ we have $\left( \mathbf{\Lambda}^{\prime} \mathbf{\Lambda} \right)^{-1}\mathbf{\Lambda}^{\prime}\mathbf{v}_{t} \rightarrow 0$ making this term asymptotically negligible. Therefore, it is justified to view $\mathbf{v}_{t}$ as a residual or a noise shock that carries no structural interpretation. At the same time, $f_{t}$ can be interpreted as a projection of the SVAR structural shocks $\mathbf{u}_{t}$ into $\mathbb{R}^{r}$.
Impulse response functions (IRFs) can be obtained via the vector moving average (VMA) representation of the VAR. In the particular case with $p=1$ lags (for notational simplicity) the VMA form becomes
where $\mathbf{\mu} = \left( I - \mathbf{\Phi}_{1} L \right)^{-1}\mathbf{\phi}_0$ with $L$ is the lag operator, and $\mathbf{\Psi}_{i} = \sum_{j=1}^{i} \mathbf{\Psi}_{j-i}\mathbf{\Phi}_{1}$ with $\mathbf{\Psi}_{0} = \mathbf{I}$ kilianlutkepohl2017. The impulse response on impact is
showing that parametric restrictions in $\mathbf{\Lambda}$ correspond to structural restrictions on impact IRFs. This result is true even when considering the effect of $\mathbf{v}_{t}$, since this term has a diagonal covariance matrix and does not affect contemparenous relationships in the variables $\mathbf{y}_{t}$.
Finally, similar to the algorithm in BaumeisterHamilton2015, the proposed specification is efficient only for static sign restrictions. This is because dynamic restrictions are nonlinear and cannot be represented (in a straightforward way) as equivalent parametric inequality restrictions in the linear VAR parameters. Therefore, dynamic restrictions would require to turn to less efficient sampling schemes similar to Ariasetal2018. In practice, there is little consensus in economic theory about the signs of structural impulse responses at longer horizons CanovaPaustian2011, and for that reason the vast majority of empirical papers impose restrictions only on impact kilianlutkepohl2017. Nevertheless, as Uhlig2017 notes, it can be quite useful to have the option to impose sign restrictions in the longer-run responses of macroeconomic variables to shocks. Within the context of the proposed methodology, this issue can be addressed if equation (ref) is specified as a dynamic factor model instead of a static factor model. However, this more general modeling assumption would require to rely on filtering and smoothing sampling steps that would, in turn, lead to a completely different estimation algorithm compared to the algorithm presented in this paper. As a result, I exclusively focus here on impact (contemporaneous) structural restrictions in small and large VARs. Extending to the case of long-horizon sign restrictions is conceptually feasible, but it is left for future research.
The previous discussion established that the factor decomposition of the VAR disturbances projects the $n$ shocks into $r$ structural plus $n$ nuisance shocks. Therefore, the first identification issue relates to being able to separate the common component ($\mathbf{\Lambda}\mathbf{f}_{t}$) from the idiosyncratic one. Notice that the original VAR covariance matrix $\mathbf{\Omega}$ has $n(n+1)/2$ free elements, while the right-hand side of the factor decomposition in (ref) has $nr + n$ free parameters. Therefore, the first condition is that $n(n+1)/2 \ge nr + n$ or that $r \le (n-1)/2$, which implies that the common component will always be identified even if the factors (structural shocks) are not identified. This condition implies that in a 19-variable VAR a reasonable number of nine factors/shocks can be estimated. Next, restrictions are required for uniquely identifying the factors, which are also structural shocks. Following anderson1956 an additional $r(r-1)/2$ restrictions are needed in order to deal with the rotation problem. This is due to the fact that rotating $\mathbf{\Lambda}$ and $\mathbf{f}_{t}$ using an orthogonal matrix $\mathbf{P}$ leads to an observationally equivalent solution, that is,
where $ \widetilde{\mathbf{\Lambda}} = \mathbf{\Lambda} \mathbf{P}$ and $\widetilde{\mathbf{f}}_{t} = \mathbf{P}^{\prime} \mathbf{f}_{t}$.
From an estimation perspective, the above equation shows that sampling of unique values of the factors in equation (ref) cannot be achieved without these additional $r(r-1)/2$ restrictions. From a structural perspective, this condition shows that the same number of restrictions is required for identification of the structural factors (shocks) in the VMA form of equation (ref). Consequently, under the assumption that $\mathbf{f}_{t} \sim N(\mathbf{0},\mathbf{I})$, placing restrictions on $\mathbf{\Lambda}$ ensures both unique estimation of factors and identification of the structural model. For instance, anderson1956 show that setting to zero all elements of $\mathbf{\Lambda}$ above the main diagonal achieves the $r(r-1)/2$ restrictions required for identification of the factors. Additionally, the diagonal elements of $\mathbf{\Lambda}$ can be normalized to be nonegative, such that the sign of the factors is also always identified. However, as shown in detail in StockWatson2005 numerous other identifying assumptions can be used in structural factor models.
In this paper interest lies in structural identification via sign (and possibly some zero) restrictions. However, many other restrictions can be incorporated in a straightforward way. For example, a researcher may want to impose that the impact response of variable $i$ to shock $j$ is lower in magnitude than the response of variable $k$ to the same shock. Due to the fact that impact restrictions on the IRFs are equivalent to imposing restrictions to $\mathbf{\Lambda}$ (see equation (ref)), this magnitude restriction can be represented as $\mathbf{\Lambda}_{ij}<\mathbf{\Lambda}_{kj}$. In the proposed model, a large class of desired structural restrictions can be represented using parametric inequalities that are imposed upon estimation of $\mathbf{\Lambda}$. Following Geweke1996, these parametric inequality restrictions can be sampled efficiently using the Gibbs sampler. The next subsection derives such a Gibbs sampler algorithm, with particular focus on ensuring computational efficiency in high-dimensions.
Posterior sampling in the reduced-form VAR with factor structure in the residuals is straightforward due to the fact that posterior conditional distributions have standard forms. To see this, write the model using a single equation for convenience
Assume that all sign and zero restrictions in $\mathbf{\Lambda}$ are collected into a matrix $\mathbf{S}$, with entries $+1$ for positive signs, $-1$ for negative signs, $0$ for zero restrictions, and a missing value for no restriction (this case is denoted in this paper using the symbol NA, and in the code using the MATLAB value $NaN$). The priors for the VAR parameters are of the form
for $i=1,...,n$, $j=1,...,r$, where $\mathbf{\Phi}_{i}$ is the $i^{th}$ row of $\mathbf{\Phi}$, $\sigma_{i}^{2}$ is the $i^{th}$ diagonal element of the matrix $\mathbf{\Sigma}$, and $\delta_{0} \left( \mathbf{\Lambda}_{ij}\right)$ is the Dirac delta function for $\mathbf{\Lambda}_{ij}$ at zero (i.e. a point mass function with all mass concentrated at zero).
The joint posterior of the parameters is the distribution $p\left(\bm{\Phi}, \bm{\Lambda},\bm{F},\bm{\Sigma} \vert \bm{y}\right) \equiv p\left(\lbrace\bm{\phi}_{i}\rbrace_{i=1}^{n},\lbrace \mathbf{\Lambda}_{i} \rbrace_{i=1}^{n},\lbrace \mathbf{f}_{t}\rbrace_{t=1}^{T}, \lbrace \sigma_{i} \rbrace_{i=1}^{n} \vert \mathbf{y} \right)$, which by Bayes theorem is the product of the normal likelihood function implied by equation (ref) and the prior distributions presented above. This product is a complicated function making sampling from the joint posterior infeasible. However, the conditional posteriors are tractable and trivial to derive in this linear model. Therefore, Bayesian inference in this VAR breaks down to sequentially sampling from the following conditional posterior distributions \\ Factor-based structural restrictions (FSR) algorithm
Step 1 is efficient because autoregressive coefficients can be sampled equation-by-equation. Further speed enhancements can be achieved by using the sampling methodology of Bhattacharya2016 in the case where the prior covariance matrices $\underline{\mathbf{V}}_{i}$ are diagonal. This is the case here, as $\underline{\mathbf{V}}_{i}$ is diagonal and its elements follow the horseshoe hierarchical structure of Carvalhoetal2010 which has the form
This prior belongs to the class of local-global shrinkage priors, that is a prior that penalizes the likelihood and shrinks coefficients towards zero. In this prior, $\psi_{i,j}$ is the local shrinkage parameter for each scalar coefficient $\bm{\phi}_{i,j}$ while $\tau_{i}$ is the global shrinkage parameter pertaining to all coefficients in equation $i$. Unlike the popular (in macroeconomic VARs) Minnesota prior that typically requires subjective tuning KoopKorobilis2010, the horseshoe prior is tuning-free as the local and global shrinkage parameters have their own distributions and are, thus, updated by information in the data. The horseshoe prior has established posterior consistency properties when used in a variety of high-dimensional settings\footnote{See online Appendix for more details and citations.} making it an appropriate choice for penalized estimation in both smaller and higher dimensional VARs. Sampling from the truncated Normal posteriors in step 2 of the algorithm above can be done using the recent contribution of Botev2017.\footnote{In an ideal world one would want to sample the vector $\mathbf{\Lambda}_{i}$ in one step and unconditionally from the factors $\mathbf{f}_{t}$, in order to reduce correlation in the MCMC chain. In practice, I follow Geweke1996 and sample each element $\mathbf{\Lambda}_{ij}$ conditional on $\mathbf{\Lambda}_{ij}$, as this conditioning allow for the algorithm to be extended more easily (e.g. if structural breaks or time-varying loadings are required). This comes at the cost of thinning the Gibbs chain by a factor of 100, that is, every 100-th sample is stored in order to ensure posterior samples are uncorrelated. In all results in this paper the Gibbs chain runs for 550,000 iterations where the first 50,000 iterations are discarded and from the final 500,000 iterations I store every 100th iteration, leaving a total of 5,000 samples from the parameter posterior for inference.} Finally, sampling of the factors $\mathbf{f}_{t}$ for each $t=1,...,T$ is fairly fast for monthly or quarterly macroeconomic data, as is sampling of the scalar parameters $\sigma_{i}^{2}$ for each $i=1,...,n$. Computational details and further discussion on the excellent properties of the horseshoe prior are provided in the online Appendix.
Given that the horseshoe prior requires no subjective tuning, there are only a handful of prior hyperparameters in the whole VAR that need to be elicited by the researcher. The parameter affecting primarily structural identification, is the choice of the prior variance $\underline{h}_{ij}$, regardless of whether the associated loading parameter $\mathbf{\Lambda}_{ij}$ should be truncated (sign restricted) or not. Due to the fact that the coefficients $\mathbf{\Lambda}$ enter the VAR in a linear way, it turns out that a fairly large value of $\underline{h}_{ij}$ implies a diffuse prior that does not impact the posterior asymptotically. Therefore, for the remainder of this paper, I set $\underline{h}_{ij} = 100$ which is a fairly diffuse choice for the typical scale of variables encountered in macroeconomic VARs. Finally, I set $\rho_{i}=1$ and $\kappa_{i}=0.01$ for all $i$, leading to a prior mean of 0.16 and substantially large prior variance, reflecting the belief that the scale of the idiosyncratic variances $\sigma_{i}^{2}$ should be small and most of variability in the VAR disturbances should come from the common component $\mathbf{\Lambda} \mathbf{f}_{t}$. Model fit in this VAR can be assessed with the Deviance Information Criterion (DIC) of Spiegelhalteretal2002, due to its simplicity of implementation. The formula and justification for the use of the DIC is provided also in the online Appendix. What suffices to remember is that the DIC has the same interpretation as any other information criterion, that is, lower values signify better fit. The DIC can be used to test any kind of parametric restrictions that are of interest, whether it pertains to lag length selection, number of structural shocks, linearities vs nonlinearities, and so on.
The most important computational aspect of the new algorithm is that samples from the restricted $\mathbf{\Lambda}$ matrices are always accepted, making it very efficient in high-dimensions. This feature is in contrast with a large class of accept/reject algorithms used especially for sign restrictions in VARs; see for example the Bayesian algorithms of RubioRamirezetal2010 and BaumeisterHamilton2015 as well as the algorithm of OuliarisPagan2017. In RubioRamirezetal2010, for example, one has to first obtain posterior samples from the VAR covariance matrix $\mathbf{\Omega}$ and then rotate its Cholesky factor $\mathbf{P}$ using randomly generated orthogonal matrices $\mathbf{Q}$. If the random rotation $\mathbf{H} = \mathbf{P}\mathbf{Q}$ satisfies the required sign restrictions then $\mathbf{H}$ is a draw from the desired matrix of contemporaneous structural relationships. Inevitably, such an accept/reject scheme is deemed to fail in high-dimensions, when the desired restrictions may be so tight that no sample of $\mathbf{H}$ can be accepted by random chance, see also the discussion in Section 13.6.4 of kilianlutkepohl2017. Many authors express the belief that when the accept/reject algorithm results to a high rejection rate, this is evidence that identification is sharp.\footnote{See for example “principle 7” in Uhlig2017 and the corresponding discussion.} However, this premise is not testable in a statistical sense, and can be misleading especially in high dimensions: a high rejection rate could either be because of the researcher imposing too many restrictions, or because of imposing restrictions that simply do not comply with the evidence in the data.
In this section, the properties of the new algorithm are explored using artificially generated data. The core exercise involves generating multivariate time series from a data generating process (DGP) that fully matches equation (ref), and estimating parameters and impulse response functions based on time series generated from this DGP. I first implement this experiment assuming that a correctly specified model is estimated using the artificial data. Subsequently, various cases of misspecification errors during the estimation process are considered -- that is, I estimate models that do not perfectly match the correct DGP.
The DGP is of the form
The DGP parameters $\widehat{\mathbf{\Phi}}$, $\widehat{\mathbf{\Lambda}}$, $\widehat{\mathbf{\Sigma}}$ are based on estimates of a VAR on real data. First, monthly data on 14 monthly macroeconomic variables are collected\footnote{The variables are: 1) real GDP, 2) GDP deflator, 3) federal funds rate, 4) commodity price index, 5) total reserves, 6) nonborrowed reserves, 7) S&P 500, 8) M1 , 9) unemployment rate , 10) industrial production , 11) employment , 12) CPI , 13) core CPI , 14) core PCE. More details on these variables is provided in the online Appendix.} for the US over the period 1965M1 - 2007M12, providing $\widehat{T}=516$ observations.\footnote{In practice, I generate $\widehat{T} + 1000$ observations and discard the first 1000 observations in order to remove dependence to the initial values of the generated time series process.} At a second step, an estimate $\widehat{\mathbf{\Phi}}$ is obtained by applying OLS to an unrestricted VAR($12$) estimated with these 14 observed US variables. The third step is to obtain the first $r=3$ principal components of these OLS residuals, and store the estimate $\widehat{\mathbf{\Lambda}}$ using OLS in a regression between the VAR residuals and their principal components. Finally, the residuals from this latter regression provide the elements of the diagonal matrix $\widehat{\mathbf{\Sigma}}$, by means of equation-by-equation application of the usual least squares formula for the variance.
While it is not possible, or even interesting, to display all estimates $\widehat{\mathbf{\Phi}}$ used as input in the DGP, it is instead interesting to look at the estimates $\widehat{\mathbf{\Lambda}}$ obtained using the procedure described above. This is because both the signs and the magnitudes of the implied IRFs in the true DGP will be affected by those estimates. Panel (A) of (ref) shows the OLS estimates, where the diagonal is normalized to be one, by dividing each element in the $m^{th}$ column of $\widehat{\mathbf{\Lambda}}$ with the original value of its $m^{th}$ element, $m=1,2,3$. While this matrix is the outcome of using real data and applying simple principal components followed by OLS estimates (which carry no structural restrictions), looking at the signs of the loadings of the first three variables (output, inflation, interest rate) allows for the classification of the three pseudo-shocks as aggregate supply, aggregate demand, and monetary policy, respectively. The estimated magnitudes, of course, are not necessarily economically meaningful. For example, in the first column a shock of $1\%$ in GDP increases inflation by $1.49\%$, which is probably not a representative magnitude for a true aggregate supply shock. Nevertheless, this is an exercise where the main aim is to check the numerical precision of the new algorithm, so the estimates in panel (A) of (ref) are perfectly valid inputs for a DGP. Finally, panel (B) of (ref) shows the sign restrictions imposed on $\mathbf{\Lambda}$, that is the matrix $\mathbf{S}$ introduced in equation (ref). These comply with the signs imposed in the DGP, and in 11 instances no sign restrictions are imposed (these entries are denoted as NA). The choice of which signs are known in $\mathbf{S}$ is random, and the next subsection looks at varying assumptions about how many sign restrictions are known to the researcher.
\vskip 1cm
For estimation purposes five different scenarios are assumed: one correctly specified case and four misspecified cases. These are denoted as C1-C5, and are defined as follows:
500 datasets of size $T=516$ are generated and posterior mean estimates of all parameters, IRFs and DICs from all five cases above are obtained. Results presented next are based on the distribution of the posterior means over these 500 generated datasets.
Before evaluating precision of estimates over the Monte Carlo iterations, it is important to first evaluate general model fit using the DIC. (ref) shows the value of the deviance information criterion attained by the estimates of the VAR model in each of the five specification cases. Because case C2 refers to a VAR with $n=8$, it is impossible to directly compare it with the other four cases that assume $n=14$.\footnote{Information criteria can only be used to compare models with the same dependent variable $\mathbf{y}$.} For that reason I present two DIC metrics, a full one based on all $n=14$ variables (with no value available for C2) and a reduced DIC which is the same formula evaluated only in the first eight VAR equations (which are common to all five cases). These are labeled in (ref) as $DIC_{14}$ and $DIC_{8}$, respectively. According to both subsets of criteria, the correctly specified estimated model, case C1, is the best one as it attains the lowest DIC value. Interestingly, the case where an additional fourth shock is incorrectly estimated (C5), doesn't seem to harm estimation accuracy; at least not as much as the case of estimating one less shock (C4). This result makes sense because shocks in the proposed VAR are equivalent to factors. By far the worst type of misspecification seems to be the one related to the lag-length. This is a characteristic of the VAR model rather than a “problem” with the specific algorithm or prior. As long as the true DGP has $p=12$ important lags, estimating the VAR with $p=2$ provides a huge loss of structurally meaningful information. In contrast, reducing the VAR from $n=14$ (which is the truth in the DGP) to $n=8$ as in case C2 is less harmful for the general fit of the model. This is probably because the missing six variables are additional measures of output and prices, that do not offer more information compared to the real GDP and GDP deflator variables that are included in both the eight- and fourteen-variable VARs.
\vskip 0.5cm
Next, estimation accuracy of the proposed algorithm has to be evaluated. Since the main focus of sign restrictions algorithms is on impulse response analysis, I compare precision of the estimated impulse response functions using all generated datasets. IRFs are combinations of all VAR parameters $\mathbf{\Phi},\mathbf{\Lambda},\mathbf{f},\mathbf{\Sigma}$, therefore comparing their precision provides a convenient summary of overall estimation precision in a VAR model. (ref) shows the responses of the first three variables in the VAR to the three identified pseudo-shocks, in the correctly specified case (C1). Green solid lines are medians over the posterior IRFs in the 500 estimated VARs using an equal number of artificial datasets. Shaded areas show the 90% probability bands of these IRFs. Finally, black dashed lines show the true IRFs implied by the parameters that are fed into the DGP. The 90% bands always include the true IRF, which suggests that estimation precision is satisfactory. The online Appendix shows identical graphs for the four misspecified cases C2-C5. These graphs become a visual confirmation of the numerical results in (ref), that is, case C5 quite precisely captures the path of the true IRFs, while case C3 results in the largest estimation errors. \vskip 0.5cm
The next Section makes clear that in the context of the empirical application in Furlanettoetal2017, the new algorithm is multiple times faster than the algorithm of RubioRamirezetal2010 in a six-variable VAR with five identified shocks. Nevertheless, it would be interesting to use artificial data in order to provide more thorough evidence on how fast the factor sign restrictions algorithm is, and how large a VAR it can scale to. For that reason artificial data are generated from the same DGP described in equations (ref)-(ref) for various values of the key parameters that affect the dimensionality of the VAR, namely $T$, $n$ and $r$. Due to the fact that this exercise pushes the VAR dimension $n$ to very large values, I fix $p=1$ in order to be able to ensure that the VAR process in the DGP is always stationary, and generation of explosive data is excluded. For the purposes of this exercise I set $\mathbf{\Phi} = 0.9\mathbf{I}_{n}$, $\mathbf{\Lambda}_{ij} \sim U(-1,1)$ and $\Sigma_{i} \sim U(0,1)$, for all $i=1,...,n$ and $j=1,...,r$. During estimation $nk$ sign restrictions are imposed, simply by obtaining the signs of the randomly generated matrix $\mathbf{\Lambda}$.\footnote{The purpose of this exercise is not to estimate meaningful restrictions, rather just to measure times. In this case, I impose the maximum number of restrictions possible on $\mathbf{\Lambda}$ in order to test the new algorithm in a worst-case scenario where all $nk$ of its elements are restricted and have to be generated from a truncated Normal posterior.}
(ref) shows the average, over 10 Monte Carlo iterations, machine time in minutes (defined as the total estimation time in seconds divided by 60 and then rounded to the nearest integer) needed to obtain $10,000$ draws from the posterior of all parameters after discarding $2,000$ draws (hence, $12,000$ draws in total). These results show that in a huge-dimensional VAR with $n=100$ series, $T=500$ observations, and $r=20$ shocks, it only takes 25 minutes to obtain $10,000$ draws from all parameter matrices, including the $1000$ sign-restricted elements in $\mathbf{\Lambda}$. For the smaller model with $n=15$ -- which is already much larger than the vast majority of models considered in the sign restrictions literature -- it only takes less than five minutes to obtain the same number of draws when $T=500$, and only one minute when $T=200$. These fantastic timings justify the choice to focus on carefully developing a Gibbs sampler that is computationally efficient.\footnote{The Gibbs sampler typically loses efficiency when there is high correlation in the samples from the posterior. In the online Appendix I show that, in order to draw $\mathbf{\Lambda}_{ij}$ from univariate (instead of the intractable multivariate) truncated Normal conditionals, we need to condition on $\mathbf{\Lambda}_{-ij}$, i.e. the set of all elements of $\mathbf{\Lambda}$ excluding the $ij^{th}$. This conditioning increases correlation relative to sampling directly the full matrix $\mathbf{\Lambda}$. However, inefficiency factors for the Gibbs sampler in the linear factor-VAR specification are still quite low (MCMC diagnostic results are available upon request). Additionally, given the ability of the algorithm to obtain quickly tens of thousands of draws from the posterior, concerns about possible correlation of draws can be alleviated by doing “thinning” -- i.e. the procedure of storing only every $\rho^{th}$ sample from the posterior, where is $\rho$ is the order of the highest significant autocorrelation in the chain.}
The results above are based on code written in MATLAB2019b and run in a personal computer with Intel Core i7 8700K, tuned at 4.9Ghz, and 32GB of RAM. Note that the Gibbs sampler algorithm iterates over each VAR equation independently and, thus, significant speed improvements can be achieved by taking advantage of parallel processing abilities of modern computers and high-performance clusters (HPCs). In MATLAB this is as simple as replacing for loops with parfor loops. Therefore, the algorithm indeed allows the estimation of arbitrarily large VAR models, as it is claimed in the Introduction.
In practical situations, the only issue that might inhibit the performance of the algorithm (and any Monte Carlo-based algorithm, to that effect) is the fact that in very large dimensions we may be sampling parameters $\mathbf{\Phi}$ in a region of the posterior that implies nonstationarity of the VAR. In order to make sense out of impulse response functions, forecast error variance decompositions, historical decompositions etc, we need to make sure we maintain only samples from the posterior which are stationary. For that reason it is important to stress that, throughout my experiments, the horseshoe prior does a great job (especially relative to a subjectively chosen Minnesota prior) in shrinking the coefficients $\mathbf{\Phi}$ towards a more numerically stable region of their posterior, where the VAR model is stationary.
In this section I revisit the empirical exercise in Furlanettoetal2017, who aim to measure various financial shocks to the US economy.\footnote{The online Appendix provides the results of an additional numerical exercise (measuring optimism shocks) that builds on Ariasetal2018.} Given computational restrictions, due to their use of the RubioRamirezetal2010 accept/reject algorithm, Furlanettoetal2017 end up estimating a series of smaller VARs in order to sequentially measure and label interesting financial shocks, such as uncertainty, credit and housing. Before illustrating how to use the new algorithm to collectively measure all these shocks in one high-dimensional data setting, I first replicate their benchmark results using a smaller VAR. That way, the new algorithm can be contrasted against the output of the RubioRamirezetal2010 algorithm, something that is not computationally feasible in the large VAR case.
Among all VAR specifications they use in their work, Furlanettoetal2017 specify a baseline VAR specification with $p=5$ lags, using data on real GDP, consumer prices, interest rate, investment-to-output ratio, stock prices, and the external finance premium.\footnote{The external finance premium is defined as the spread between yields on Baa rated bonds and the federal funds rate. Notice that the three variables that are not already expressed as rate, ratio, or spread (i.e. GDP, consumer prices, and stock prices), are transformed only using logarithms of the levels and not growth rates. Also note that these authors use a noninformative (uniform) prior, while I use the shrinkage horseshoe prior described in Section 2.} All data are for the 1985Q1 - 2013Q2 period. The online Appendix provides exact details of all series and transformations used, which in this case they are identical to those reported in Table 11 of Furlanettoetal2017. Five shocks in total are identified by the authors using the six-variable baseline VAR. The names of the shocks and the associated sign restrictions adopted are shown in (ref). The first four shocks are standard macro-related shocks, and of interest is the fifth shock which is a generic financial sector shock.
(ref) shows the effects of a financial shock identified as a shock that causes GDP, consumer prices, stock prices, interest rate and the investment/output ratio to react positively contemporaneously. The sign of the spread is left unrestricted. Panel (a) replicates the impulse responses also shown in Figure 1 of Furlanettoetal2017, produced using the algorithm of RubioRamirezetal2010. Panel (b) shows the same responses produced by application of the new algorithm for sign restrictions. The responses on impact in both panels are of almost identical magnitude, showing that the new algorithm produces sensible results. Any observed differences in the propagation of the impulse responses in subsequent periods, especially for GDP, prices and investment/output ratio, is due to the different estimates of the autoregressive coefficients (Furlanettoetal2017 use noninformative priors).\footnote{The online Appendix replicates (ref) by exchanging the horseshoe prior for a diffusing (flat) prior. This can be done by simply droping the local-global hyperparameters of the horseshoe and, instead, setting $\mathbf{V}_{i} = c \times \mathbf{I}$ for $c \rightarrow \infty$. In this case, the normal prior becomes locally uniform on the parameter support. (ref) in the online Appendix reveals that in the case of this noninformative prior the shapes of the shocks between the two algorithms become identical. However, the error bands in the new algorithm are still sharper. This is because Furlanettoetal2017 sample the VAR covariance matrix from the standard inverse Wishart posterior KoopKorobilis2010, while the new algorithm samples the covariance matrix from the more parsimonious factor model.} Following up on the discussion in the previous section, it takes roughly four hours to obtain 2000 draws from the Furlanettoetal2017 using their MATLAB code and exact numerical settings based on the RubioRamirezetal2010 algorithm. In contrast, using the same PC\footnote{Specifications of the PC are reported in (ref).} it takes 20 minutes to obtain $600,000$ draws from the proposed Gibbs sampler (where out of these $600,000$ draws I discard $100,000$ and then save every $100^{th}$ draw, leading to $5,000$ draws from the posterior of VAR parameters and impulse response functions). Similar conclusions can be made for all other shocks in the system (supply, demand, monetary, investment), where impulse responses are qualitatively similar. Plots for these shocks are provided in the online Appendix.
The new algorithm relies on joint estimation of parameters and identification restrictions. Therefore, it could be argued that the qualitatively similar results in (ref) are an artifact as they can be very sensitive to the structural identification restrictions imposed. However, this is not the case and the algorithm works well in various different scenarios. I consider the following identification restrictions in $\mathbf{\Lambda}$
The first case is the one also plotted in panel (b) of (ref). The second case is used as a means of showing that the algorithm is able to incorporate the case with as many factors as variables, and is not sensitive to the motivating assumption that only a few shocks drive the VAR. This motivating assumption seems reasonable in larger systems (see next subsection), but what if a researcher is interested in smaller systems with as many shocks as variables? The third case allows to find out to what extend identification of the financial shock is affected by the restrictions in the remaining four shocks. Since shocks are identical factors, the aim is to find how estimates of the fifth factor are affected by assumptions in the first four factors. Finally, case four simply removes any information in the first four shocks and simply estimates a model with one shock. Identifying a single shock of interest is a very popular practice in empirical papers that rely on the RubioRamirezetal2010 algorithm, as it allows for faster inference (less accept/reject algorithmic steps) and results remain quantitatively unchanged. In contrast, the current algorithm is affected by the assumptions on the number of restrictions. In the fourth case using one shock means that only one factor is used for estimation, which in turn means that estimates of the VAR parameters (covariance matrix, and coefficients of lagged variables) will be affected.
(ref) presents the impulse responses from all four cases. They are qualitatively and quantitatively identical. There are only a couple of differences when moving from the models with five shocks (factors) to the model with only one factor (panel (d) of the figure). In the latter case, the $68\%$ bands of the IRFs of prices are a bit wider, and the median impact response of the spread is twice as large (in absolute value) as the other three cases. Additionally, the curvature of the IRFs of GDP, interest rate, investment/output ratio, and stock prices, over the first 10 periods following the shock, is more pronounced in panel (d) relative to the other three panels. Therefore, (ref) suggests that the number of shocks is more important than identification assumptions made in shocks other than the shock of interest.
In the previous Section, using artificially generated data, it was suggested that it can be hurtful to estimate less shocks compared to the true number of shocks in the DGP. In this case, when estimating $r=1$ shock does not distort the IRFs substantially, as from a statistical point of view one factor can be sufficient for a small, six-variable VAR. In order to find out if this is truly the case, (ref) shows the DIC values attained by each of the four cases. The worst case is the one where six shocks are identified from the six series -- this corresponds to an overparametrized and unnecessary (from a statistical point of view) factor decomposition. The case with $r=1$ has the second highest DIC value, meaning that the true number of factors (again in a statistical sense) is larger than one and smaller than six. Surprisingly, Case 3 which is builds on the benchmark Case 1 but lifts all sign restrictions in the first four shocks, is the one that fits the best from a statistical point of view.\footnote{Notice that in this case the IRFs of supply, demand, monetary, and investment will be flat around zero -- as a matter of fact these should be named as Shock 1, Shock 2, Shock 3, Shock 4 exactly because there are no sign restrictions or identifying assumptions. As it was argued in the previous Sections, even though these four shocks are not structurally identified, the common component $\mathbf{\Lambda}\mathbf{f}_{t}$ is identified and covariance matrix estimation is feasible. Despite the fact that these shocks are not identified, the fifth shock (financial) is identified based on its restrictions and it is not affected by the lack of identification of the first four shocks/factors.} This DIC value suggests that identifying supply, demand, monetary and investment shocks, using sign restrictions based on economic theory, is statistically inferior to estimating four generic shocks with no sign restrictions. However, the model with the best statistical fit is not necessarily the best model for structural analysis Bernankeetal2005, therefore, in this case the benchmark Case 1 should be preferred as it has a good fit and at the same time allows for identification of important macroeconomic shocks.
We next proceed to demonstrate how the new algorithm can estimate one, large-dimensional system in order to measure in one setting all the financial shocks that Furlanettoetal2017 identify. The larger VAR that these authors specify has seven variables and six shocks: aggregate supply, aggregate demand, investment, housing, uncertainty, and credit. These authors do not identify a monetary shock using this larger VAR, possibly due to computational concerns. Here we attempt to use all available variables in Furlanettoetal2017 to identify seven shocks, that is, the six shocks just listed plus a monetary shock. We also use additional measures of output, consumer prices, stock prices, interest rate, and credit spread, in order to enhance identification. We end up with a 15-variable VAR with $p=5$ on the following variables: 1) real GDP; 2) prices (GDP deflator); 3) interest rate (3-month Tbill); 4) investment to output ratio; 5) stock prices (real S&P500 prices); 6) credit spread (Baa minus Fed funds rate); 7) credit to real estate value ratio; 8) excess bond premium (EBP); 9) EBP to VIX ratio; 10) mortgage rate (30-year rates); 11) employment; 12) Federal funds rate; 13) core CPI; 14) stock prices 2 (real DJIA prices); and 15) credit spread 2 (“GZ” spread). The online Appendix has detailed definitions of these variables, transformations used, and sources.
(ref) shows the signs imposed on each of the 15 variables in order to identify each of the seven structural shocks. This is a large matrix of restrictions, but the new algorithm can handle computationally the task of drawing $600,000$ samples from the posterior of all parameters (including the structural matrix of contemporaneous shocks) in a matter of minutes. As it was the case with the baseline VAR above, out of these $600,000$ draws $100,000$ are discarded and every $100^{th}$ sample is stored, resulting in $5,000$ samples used to produce numerical results from this large model. The horseshoe prior also has a crucial role in the estimation of this model, as we have $1140$ parameters in $\mathbf{\Phi}$ and only $114$ observations for each of the $15$ endogenous variables.
(ref) shows the impulse responses of the 15 endogenous variables to a credit shock. The green lines are posterior medians, and the shaded areas 68% bands. The magnitudes and shapes of the IRFs are consistent with the ones reported in Furlanettoetal2017, despite the fact that in the case of variables such as GDP the IRFs are strongly different from zero. The most interesting feature of this figure is the effect of a credit shock on the two credit spread variables we used in the same VAR. Furlanettoetal2017 use these spreads (plus an additional third spread we haven't included here) one at a time in their VAR in order to assess robustness of their results. These authors do not impose sign restrictions on the credit spread and they find that in their baseline specification this tends to be negative. In the large VAR case, the first credit spread variable has a strong negative contemporaneous response before subsequently moving to positive territory, while the second credit variable does not have a contemporaneous response different from zero and in subsequent period reacts positively. Such results show the important avenues for identifying various structural shocks that the new algorithm opens up: by using large information sets we can have the ability to identify several structural shocks in one setting, thus, making comparisons and testing of structural hypotheses more transparent. The online Appendix provides further results for this 15-variable VAR.
This paper outlines a new algorithm based on a VAR methodology that fully utilizes the interpretability and parsimony of factor models. In particular, the novel element of the proposed approach is the formulation of reduced-form VAR disturbances using a common factor structure, and the derivation of an algorithm that allows for efficient sampling of sign-restricted decompositions of the VAR covariance matrix. The new algorithm can handle VARs with possibly 100 or more variables and it provides sensible numerical results compared to the algorithm of RubioRamirezetal2010 -- despite the fact that the two algorithms rely on different modeling assumptions and are not directly comparable.\footnote{Additional numerical results are provided in the online Appendix.} Therefore, the new algorithm can be seen as a useful tool in the toolbox of modern macroeconomists, that complements existing algorithms and at the same time opens up new avenues for empirical research using large-scale VAR models.
\addcontentsline{toc}{section}{References}
\onehalfspacing
\vskip 1cm