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.
101,432 characters · 14 sections · 107 citation commands
Large Bayesian VARs with Factor Stochastic Volatility: Identification, Order Invariance and Structural Analysis
\onehalfspacing
\thispagestyle{empty}
Bayesian vector autoregressions (VARs) with multivariate stochastic volatility, first developed in CS05 and primiceri05, are now the workhorse models in empirical macroeconomics. These multivariate stochastic volatility models, however, have the undesirable property that the implied likelihoods are not invariant to the order of the dependent variables.\footnote{This non-invariance problem is explicitly acknowledged and discussed in both CS05 and primiceri05. See also the discussion in CCM19.} This ordering issue has become an increasingly pertinent problem due to two prominent developments in the VAR literature. First, in the last two decades there has been a gradual departure from conventional recursive or zero identification restrictions to other more credible identification schemes---such as identification by sign restrictions Faust98,CD02, Uhlig05---that do not restrict the order of the variables. Despite this development, models of CS05 and primiceri05 continue to be used to first obtain reduced-form estimates, which are then taken as inputs in the subsequent structural analysis. Since the reduced-form estimates are not order invariant, the results from the structural analysis depend on the order of the variables in a subtle way, often without explicit recognition by the user.\footnote{The implications of this non-invariance problem for structural analysis have been illustrated in Bognanni18 and Hartwig2019.}
Second, following the seminal contributions by BGR10 and koop13, there is an increasing desire to use large VARs involving more than dozens of dependent variables for structural analysis. This development is partly motivated by the concern of informational deficiency of using a limited information set---by expanding the set of relevant variables, one can alleviate this concern HS91, LR93, LR94. However, unless there is a natural variable ordering (e.g., using recursive identification restrictions), the ordering issue becomes more severe as the number of ways to order the variables increases exponentially with the number of variables.
In view of these developments, we consider an alternative Bayesian VAR based on the factor stochastic volatility that is constructed to be invariant to the order of the dependent variables. Factor stochastic volatility models are commonly used for modeling high-dimensional financial data, but are less widely employed in empirical macroeconomics.\footnote{A notable exception is KH20, who use Bayesian VARs with factor stochastic volatility for macroeconomic forecasting. CCM18 consider a related multiplicative 2-factor stochastic volatility model to study the impact of macroeconomic and financial uncertainty.} In specifying a suitable factor stochastic volatility model, there is often a tension between identification and order invariance. On the one hand, one can identify the factors and the associated factor loadings by fixing the orientation of the factors GZ96,CNS06. But this identification strategy essentially fixes the order of the variables, and therefore the identified model is not order invariant CLS18. On the other hand, one could avoid fixing the orientation of the factors and obtain an order-invariant model, but it is unclear that the factors and the loadings are identified.\footnote{For example, Kastner19 does not impose any orientation restrictions on the factors, arguing that identification of the factor loadings is not necessary for his purpose of estimating the reduced-form covariance matrix.}
We solve this dilemma between achieving identification and order invariance by carefully teasing out a set of conditions strong enough for identification, yet they are weak enough that the model remains order invariant. More specifically, we construct a VAR in which the innovations have a factor structure, and both the factors and the idiosyncratic errors follow stochastic volatility processes. We first show that the likelihood implied by this model is invariant to the order of the dependent variables. We then discuss sufficient conditions for identification of the factors and the factor loadings, building upon the approach in SF01 and extending it to a more general setting in which both the factors and the idiosyncratic errors are heteroscedastic. Under mild regularity conditions, we show that the factor loadings under our setup are identified up to permutation and sign changes. Furthermore, with additional sign restrictions that satisfy a set of conditions, we show that the factor loadings and the associated factors are point-identified.
To determine the number of factors, we develop an estimator of the marginal likelihood based on an importance sampling approach to evaluate the observed-data or integrated likelihood. Through a series of Monte Carlo experiments, we show that our marginal likelihood estimator works well and is able to select the correct number of factors under a variety of settings.
We then discuss how our VAR with factor stochastic volatility (VAR-FSV) can be used for structural analysis. More specifically, we develop various structural analysis tools for VAR-FSV similar to those designed for standard structural VARs. In particular, we describe methods to construct structural impulse response functions, forecast error variance decompositions and historical decompositions. We demonstrate the methodology by revisiting the 6-variable VAR identified by a set of sign restrictions on the contemporaneous impact matrix considered in FRS19. We augment their system to a 20-variable VAR by including additional, seemingly relevant macroeconomic and financial variables, which helps alleviate the concern of informational deficiency. In addition, the impulse responses obtained using the VAR-FSV with the sign restrictions imposed are point-identified. Empirically, we show that by including the additional variables and sign restrictions, one can substantially sharpen inference.
Our paper is related to the recent work by Korobilis20, who uses a VAR with a factor error structure for structural analysis. His work is motivated by the computational challenge of imposing a large number of sign restrictions to obtain admissible draws using conventional accept-reject methods RWZ10. This computational hurdle has so far limited the use of sign restrictions to relatively small systems with at most half a dozen dependent variables.\footnote{Large VARs, on the other hand, are mostly identified using recursive or zero restrictions. See, for example, LSZ96, BGR10 and ER17.} Instead of using standard structural VARs, Korobilis20 assumes that the factors in his model play the role of structural shocks, and shows that in this case structural analysis can be done efficiently even when one imposes a large number of sign restrictions. His model, however, is homoscedastic, and consequently it is only set-identified. By contrast, in our VAR-FSV both the factors and the idiosyncratic errors follow stochastic volatility processes. This feature does not only accommodate the empirical finding that macroeconomic and financial variables typically exhibit time-varying volatility clark11, CR15, it also allows us to achieve point-identification of the factors and the factor loadings.
Our work also contributes to the recent literature on using heteroscedasticity to identify conventional structural VARs, including WD15, LLM10, HL14, BB20, Lewis21 and BPSS21. Our paper considers the alternative setting of a VAR with a factor stochastic volatility specification and establishes sufficient conditions for identification. One key advantage of using VAR-FSV for structural analysis, compared to structural VARs, is that under VAR-FSV it is computationally feasible to estimate large systems and impose a large number of sign restrictions.
Our work is also related to the growing literature on constructing multivariate stochastic volatility models that are order invariant. One approach is based on Wishart or inverse-Wishart processes; examples include PG2006, AM09, CDLS18 and SZ20. These models, however, are typically computationally intensive to estimate as the estimation involves drawing from non-standard high-dimensional distributions. As such, these models are generally not applicable to large datasets. An alternative approach is based on the common stochastic volatility models in CCM16 and chan20. Although these models are designed for large systems and can be estimated quickly, they are more restrictive since the time-varying error covariance matrix depends on a single stochastic volatility process---in particular, the error variances are always proportional to each other.
There are also order-invariant models that are based on the discounted Wishart process, such as those in Uhlig97, WH06 and Bognanni18. These models are convenient to estimate as they admit Kalman-filter type filtering and smoothing algorithms. The cost for this tractability, however, is that they are generally too tightly parameterized, and consequently, they tend to underperform in forecasting macroeconomic variables relative to standard stochastic volatility models such as CS05 and primiceri05 ARRS21. Lastly, the recent paper CKY21 extends the stochastic volatility model of CS05 by avoiding the use of Cholesky decomposition so that the extension is order-invariant. So far this reduced-form VAR is used for forecasting, and further research is needed to incorporate identification restrictions for structural analysis.
The rest of this paper is organized as follows. Section (ref) first introduces the VAR with factor stochastic volatility. Its theoretical properties, including order invariance and sufficient conditions for identification, are discussed in Section (ref). We then outline a posterior sampler and a marginal likelihood estimator for the model in Section (ref) and Section (ref), respectively. Next, Section (ref) develops various structural analysis tools for the VAR-FSV model, including algorithms to construct structural impulse response functions and to perform various decompositions. Then, Section (ref) presents Monte Carlo results to illustrate how well the marginal likelihood estimator works under a variety of settings. We next demonstrate the proposed methodology via a structural analysis with sign restrictions in Section (ref). Finally, Section (ref) concludes and discusses some future research directions.
In this section we outline a Bayesian VAR with factor stochastic volatility (FSV) and the associated prior distributions. To that end, let $\mathbf{y}_t$ be an $n\times 1 $ vector of dependent variables at time $t$. Then, for $t=1,\ldots, T$, consider the following VAR-FSV model:
where $\mathbf{f}_t = (f_{1,t},\ldots, f_{r,t})'$ denotes a $r\times 1$ vector of latent factors and $\mathbf{L}$ is an $n\times r$ matrix of factor loadings. Note also that $\mathbf{L}$ is unrestricted. The disturbances $\mathbf{u}_t^y$ and the latent factors $\mathbf{f}_t$ are assumed to be independent at all leads and lags. Moreover, they are specified as jointly Gaussian:
where $\boldsymbol \Sigma_t = \text{diag}(\text{e}^{h_{1,t}},\ldots, \text{e}^{h_{n,t}})$ and $ \boldsymbol \Omega_t = \text{diag}(\text{e}^{h_{n+1,t}},\ldots, \text{e}^{h_{n+r,t}})$ are diagonal matrices. For $t=2,\ldots, T$, the log-volatilities evolve as:
where we impose $|\phi_1|<1, \ldots |\phi_{n+r}|<1$ to ensure stationarity. Finally, the initial conditions follow the stationary distributions $h_{i,1} \sim \mathcal{N}(\mu_i,\sigma_i^2/(1-\phi_i^2)), i=1,\ldots, n,$ and $h_{n+j,1} \sim \mathcal{N}(0,\sigma_{n+j}^2/(1-\phi_{n+j}^2)), j=1,\ldots, r$. Note that the stationary distributions of the log-volatilities associated with the idiosyncratic errors have nonzero means, whereas the means of those associated with the factors are set to be zero for normalization.
To facilitate estimation, we rewrite the VAR in (ref) as
where $\mathbf{I}_n$ is the identity matrix of dimension $n$, $\otimes$ is the Kronecker product, $\boldsymbol \beta = \text{vec}([\mathbf{a}_0, \mathbf{A}_1, \ldots, \mathbf{A}_p]')$ and $\mathbf{x}_t = (1,\mathbf{y}_{t-1}',\ldots,\mathbf{y}_{t-p}')'$ is a $k \times 1$ vector of intercept and lagged values with $k=np+1$.
Next, we specify the prior distributions on the model parameters. Let $\boldsymbol \beta_i$ and $\mathbf{l}_i$ denote the VAR coefficients and the elements of $\mathbf{L}$ in the $i$-th equation, respectively, for $i=1,\ldots, n$. We assume the following independent priors on $\boldsymbol \beta_i$ and $\mathbf{l}_i$ for $i=1,\ldots, n$: \[ \boldsymbol \beta_i\sim\mathcal{N}(\boldsymbol \beta_{0,i},\mathbf{V}_{\boldsymbol \beta_i}),\quad \mathbf{l}_i\sim\mathcal{N}(\mathbf{l}_{0,i},\mathbf{V}_{\mathbf{l}_i}). \] We elicit the prior mean vector $\boldsymbol \beta_{0,i}$ and the prior covariance matrix $\mathbf{V}_{\boldsymbol \beta_i}$ similar to the Minnesota prior DLS84, litterman86, KK93. Specifically, for growth rates data, we set $\boldsymbol \beta_{0,i} = \mathbf{0}$ to shrink the VAR coefficients to zero. For level data, $\boldsymbol \beta_{0,i}$ is set to be zero as well except for the coefficient associated with the first own lag, which is set to be one. The prior covariance matrix $\mathbf{V}_{\boldsymbol \beta_i}$ is constructed so that it depends on two key hyperparameters, $\kappa_1$ and $\kappa_2$, that control respectively the overall shrinkage strength of `own' lags and `other' lags. For a more detailed discussion of the Minnesota prior, see, e.g., KK10, DNS11 or karlsson13.
Finally, for the parameters in the stochastic volatility equations, we assume the priors: \[ \mu_i \sim \mathcal{N}(\mu_{0,i},V_{\mu_i}), \; \phi_j\sim \mathcal{N}(\phi_{0,j},V_{\phi_j})1(|\phi_j|<1),\; \sigma_{j}^2 \sim \mathcal{IG}(\nu_{j},S_{j}), \] $i=1,\ldots, n$ and $j=1,\ldots, n+r$.
In this section we describe a few important properties of the VAR-FSV model specified in (ref)-(ref). First, the likelihood implied by the model is invariant to the order of the variables (after permuting the relevant parameters appropriately). To see that, let $\mathbf{P}$ be an $n\times n$ permutation matrix such that $\mathbf{P}\mathbf{P}' = \mathbf{P}'\mathbf{P} = \mathbf{I}_n$. For the $n$-variate Gaussian density $f_{\mathcal{N}}(\cdot ; \boldsymbol \mu,\boldsymbol \Sigma)$ with mean vector $\boldsymbol \mu$ and covariance matrix $\boldsymbol \Sigma$, it is easy to see that $f_{\mathcal{N}}(\mathbf{x} ; \boldsymbol \mu,\boldsymbol \Sigma) = f_{\mathcal{N}}(\mathbf{P}\mathbf{x} ; \mathbf{P}\boldsymbol \mu, \mathbf{P}\boldsymbol \Sigma\mathbf{P}')$.
Next, we derive an expression of the likelihood function. To that end, stack $\mathbf{h}_t^y = (h_{1,t},\ldots,h_{n,t})'$ and $\mathbf{h}_t^f = (h_{n+1,t},\ldots,h_{n+r,t})'$. We similarly define $\boldsymbol \phi_y, \boldsymbol \phi_f, \boldsymbol \sigma^{2}_y $ and $\boldsymbol \sigma^{2}_f$. In addition, we let $\mathbf{h}_t = (\mathbf{h}_t^{y'},\mathbf{h}_t^{f'})'$, $\boldsymbol \phi = (\boldsymbol \phi_y', \boldsymbol \phi_f')'$, $\boldsymbol \sigma^2 = (\boldsymbol \sigma^{2'}_y, \boldsymbol \sigma^{2'}_f)$ and $\boldsymbol \mu = (\mu_1,\ldots, \mu_n)'$. Then, the state equations (ref)-(ref) imply that the densities of $\mathbf{h}_t^y$ and $\mathbf{h}_t^f,$ for $t=2,\ldots, T$, are, respectively, \[ f_{\mathcal{N}}(\mathbf{h}_t^y; \boldsymbol \mu + \boldsymbol \phi_y\odot(\mathbf{h}_{t-1}^y-\boldsymbol \mu),\text{diag}(\boldsymbol \sigma^{2}_y)) \text{ and } f_{\mathcal{N}}(\mathbf{h}_t^f; \boldsymbol \phi_f\odot \mathbf{h}_{t-1}^f,\text{diag}( \boldsymbol \sigma^{2}_f)), \] where $\odot$ is the element-wise multiplication. Moreover, the initial conditions $\mathbf{h}_1^y$ and $\mathbf{h}_1^f$ have, respectively, the densities \[ f_{\mathcal{N}}(\mathbf{h}_1^y; \boldsymbol \mu,\text{diag}(\boldsymbol \sigma^{2}_y\oslash(\mathbf{1} - \boldsymbol \phi_y))), \text{ and } f_{\mathcal{N}}(\mathbf{h}_1^f; \mathbf{0},\text{diag}(\boldsymbol \sigma^{2}_f\oslash(\mathbf{1} - \boldsymbol \phi_f))), \] where $\oslash$ denotes the element-wise division.
Next, using the representation in (ref) and integrating out the factors, the density of $\mathbf{y}_t$ given the parameters and log-volatilities is $f_{\mathcal{N}}(\mathbf{y}_t; (\mathbf{I}_n \otimes \mathbf{x}_t')\boldsymbol \beta,\mathbf{L}\boldsymbol \Omega_t\mathbf{L}'+ \boldsymbol \Sigma_t)$. Stacking $\mathbf{y}=(\mathbf{y}_1',\ldots, \mathbf{y}_T')'$, the likelihood function, or more precisely the integrated or observed-data likelihood, can therefore be written as
Now, for an arbitrary permutation matrix $\mathbf{P}$, suppose we permute the order of the dependent variables $\widetilde{\mathbf{y}}_t = \mathbf{P}\mathbf{y}_t$ and the associated lagged values $\widetilde{\mathbf{x}}_t' = (1,(\mathbf{P}\mathbf{y}_{t-1})',\ldots, (\mathbf{P}\mathbf{y}_{t-p})') = \mathbf{x}_t'\mathbf{Q}'$, where $\mathbf{Q} = \text{diag}(1,\mathbf{I}_p\otimes\mathbf{P})$. We claim that the likelihood implied by the VAR-FSV model is invariant to the permutation $\mathbf{P}$ in the sense that \[ p(\mathbf{y} \,|\,\boldsymbol \beta,\mathbf{L},\boldsymbol \mu,\boldsymbol \phi,\boldsymbol \sigma^2) = p(\widetilde{\mathbf{y}} \,|\, \widetilde{\boldsymbol \beta},\widetilde{\mathbf{L}},\widetilde{\boldsymbol \mu},\widetilde{\boldsymbol \phi},\widetilde{\boldsymbol \sigma}^2), \] where $\widetilde{\mathbf{L}} = \mathbf{P}\mathbf{L}, \widetilde{\boldsymbol \beta} = (\mathbf{P}\otimes\mathbf{Q}) \boldsymbol \beta, \widetilde{\boldsymbol \mu} = \mathbf{P}\boldsymbol \mu, \widetilde{\boldsymbol \phi} = ((\mathbf{P}\boldsymbol \phi_y)', \boldsymbol \phi_f')'$ and $\widetilde{\boldsymbol \sigma}^{2} = ((\mathbf{P}\boldsymbol \sigma_y^2)', \boldsymbol \sigma_f^{2'})'$.\footnote{Note that the permuted vector $\widetilde{\boldsymbol \beta}$ consists of the VAR coefficients of the following system stacked by rows: \[ \widetilde{\mathbf{y}}_t = \widetilde{\mathbf{a}}_0 + \widetilde{\mathbf{A}}_1 \widetilde{\mathbf{y}}_{t-1} + \cdots + \widetilde{\mathbf{A}}_p \widetilde{\mathbf{y}}_{t-p} + \widetilde{\boldsymbol \varepsilon}_t, \] where $ \widetilde{\mathbf{a}}_0 = \mathbf{P}\mathbf{a}_0$ and $\widetilde{\mathbf{A}}_j = \mathbf{P}\mathbf{A}_j\mathbf{P}', j=1,\ldots,p.$ That is, \[ \widetilde{\boldsymbol \beta} = vec
= vec\left( \mathbf{Q}
\mathbf{P}'\right) = (\mathbf{P}\otimes\mathbf{Q})vec
= (\mathbf{P}\otimes\mathbf{Q})\boldsymbol \beta. \]} That is, we obtain the same likelihood value for any permutation of $\mathbf{y}_t$ if the lagged values and the parameters are permuted accordingly.
This claim of order invariance can be readily verified as follows. First, noting that
we therefore obtain \[ f_{\mathcal{N}}(\mathbf{y}_t; (\mathbf{I}_n \otimes \mathbf{x}_t')\boldsymbol \beta,\mathbf{L}\boldsymbol \Omega_t\mathbf{L}'+ \boldsymbol \Sigma_t) = f_{\mathcal{N}}(\widetilde{\mathbf{y}}_t; (\mathbf{I}_n \otimes \widetilde{\mathbf{x}}_t')\widetilde{\boldsymbol \beta}, \widetilde{\mathbf{L}}\boldsymbol \Omega_t\widetilde{\mathbf{L}}' + \widetilde{\boldsymbol \Sigma}_t), \] where $\widetilde{\boldsymbol \Sigma}_t = \mathbf{P}\boldsymbol \Sigma_t\mathbf{P}'$. Similarly, we also have
Since the Gaussian densities in (ref) are equal to their permuted counterparts, the integrand in $p(\widetilde{\mathbf{y}} \,|\, \widetilde{\boldsymbol \beta},\widetilde{\mathbf{L}},\widetilde{\boldsymbol \mu},\widetilde{\boldsymbol \phi},\widetilde{\boldsymbol \sigma}^2)$ is exactly the same as that in (ref). The only difference between the two integrals is the order of integration: $(\mathbf{h}_t^y,\mathbf{h}_t^f)$ versus $(\mathbf{P}\mathbf{h}_t^y,\mathbf{h}_t^f)$. But since the integral is finite, one can change the order of integration without changing the integral. Hence, the desired result follows. The following proposition summarizes this result.
Next, we discuss sufficient conditions for identification of the factor loadings and latent factors. We mainly follow the approach in SF01, but consider a more general setting in which the idiosyncratic errors $\mathbf{u}_t^y$ in (ref) are also heteroscedastic. First, note that it follows from (ref) and (ref) that the covariance matrix of $\boldsymbol \varepsilon_t$ is given by $\text{Var}[\boldsymbol \varepsilon_t \,|\, \boldsymbol \Omega_t,\boldsymbol \Sigma_t] = \mathbf{L}\boldsymbol \Omega_t\mathbf{L}'+ \boldsymbol \Sigma_t :=\boldsymbol \Gamma_t$. The covariance structure of any observationally equivalent model to (ref)--(ref) with the same number of factors must satisfy $\boldsymbol \Gamma_t = \mathbf{L}^* \boldsymbol \Omega_t^* \mathbf{L}^{*'} + \boldsymbol \Sigma_t^*$ for all $t$, where $ \mathbf{L}^*$ is $n\times r$ and $\boldsymbol \Omega_t^*$ is $r\times r$. Furthermore, for a square matrix $\mathbf{A}$ of dimension $m$, we define $\text{vecd}(\mathbf{A})$ to be the $m\times 1$ vector that stores its diagonal elements. Now, we consider the following assumptions that are used throughout the paper.
Assumption (ref) requires that no stochastic volatilities in the common factors can be expressed as a linear combination of other factor stochastic volatilities. Under our factor stochastic volatility model, this assumption is automatically satisfied. Assumption (ref) limits the extent of sparseness in the factor loadings matrix to ensure one can separately identify the common and the idiosyncratic components. This assumption can be traced back to AR56, and is widely adopted in the literature. Implicitly, it also requires that $r \leqslant (n-1)/2$. Since factor models are mostly applied to situations where the number of variables $n$ is much larger than the number of factors $r$, this is not a stringent condition.
With Assumptions 1 and 2, one can show that the factor stochastic volatility model specified in (ref)-(ref) is identified up to permutations and sign changes of the factors. The identification results are summarized in the following proposition.
We prove the proposition by adapting the results in AR56 and SF01 to our setting. The details are provided in Appendix A. Proposition (ref) contains two sets of identification results. First, it shows that the common and idiosyncratic variance components can be separately identified. Second, the common variance component is identified up to permutations and sign switches of the latent factors.
So far we have only considered the case where all $r$ factors are heteroskedastic. It turns out this is not necessary for identification of the common variance component. More generally, one can show that part of the factor loadings matrix is identified even when some of the factors are homoskedastic (their variances are normalized to be one). The following proposition summarizes such a partial identification result.
The proof of this proposition is given in Appendix A. The condition that $\text{diag}(\boldsymbol \Omega_{1t},1)$ satisfies Assumption 1---i.e., $(\text{vecd}(\boldsymbol \Omega_{1t})', 1)'$ are linearly independent for all $t$---requires all stochastic processes in $\text{vecd}(\boldsymbol \Omega_{1t})$ to be non-degenerate. (Otherwise those homoskedastic factors should be relocated to the homoskedastic part.) It is also worth noting that Proposition (ref) does not imply that for $r_1<r$, the common variance component is not identifiable. In fact, it turns out that the minimum number of heteroskedastic factors for identifying $\mathbf{L}$ (up to permutations and sign switches) is $ r_1 = r-1$. We summarize this result in the following corollary.
The reason why only $r_1 = r-1$ heteroskedastic factors are needed for identification is intuitive. Under the assumptions in Proposition (ref), when $r_1 = r-1$, only one element in $\boldsymbol \Omega_t$ is normalized to one; the remaining $r-1$ stochastic processes in $\boldsymbol \Omega_{1t}$ are linearly independent. Consequently, $\boldsymbol \Omega_t$ also satisfies Assumption 1. And Corollary (ref) follows from Proposition (ref).
For $r_1 < r-1$, part of the factor loadings matrix $\mathbf{L}$ is invariant under general orthogonal transformation. To see that, suppose $r_1 < r-1$, and hence $\mathbf{L}_2$ has at least $r_2 \geqslant 2$ columns. Let $\mathbf{R}_{f_2}$ be a $r_2\times r_2$ orthogonal matrix other than permutation $\mathbf{P}_{r_2}$ or reflection $\mathbf{P}_{\pm}$ such that $\mathbf{R}_{f_2}\mathbf{R}_{f_2}' = \mathbf{I}_{r_2}$.\footnote{The only one-dimensional orthogonal matrices are reflections, namely, $+1$ and $-1$. Hence, $r_2$ must be at least 2.} Then, we have $\mathbf{L}\boldsymbol \Omega_t\mathbf{L}' = \mathbf{L}_1\boldsymbol \Omega_{1t}\mathbf{L}_1' + \mathbf{L}_2\mathbf{L}_2' = \mathbf{L}_1\boldsymbol \Omega_{1t}\mathbf{L}_1' + \mathbf{L}_2\mathbf{R}_{f_2} \mathbf{R}_{f_2}'\mathbf{L}_2'.$ Hence, $\mathbf{L}^* = (\mathbf{L}_1,\mathbf{L}_2\mathbf{R}_{f_2})$ and $\boldsymbol \Omega^*_t = \boldsymbol \Omega_t$ form an observationally equivalent model.
For point-identification, one needs additional restrictions on $\mathbf{L}$ (or the latent factors) to pin down the unique permutation and sign configuration. As is common in macroeconomic analysis using VARs, sign restrictions implied by economic theory are often available to assist structural identification. For a recent contribution linking sign restrictions and factor models, see Korobilis20. Below we describe how we can incorporate sign restrictions to achieve point-identification. To that end, let $\mathbf{S}$ denote the $n\times r$ matrix that collects the corresponding restrictions on the factor loadings matrix $\mathbf{L}$. The entries of $\mathbf{S}$ can take four values: 1, $-1$, 0 and N/A, which denotes positive restriction, negative restriction, zero restriction and no restrictions, respectively. For example, if economic theory implies that $\mathbf{L}_{ij}>0 $, then $\mathbf{S}_{ij} = +1$; if there are no restrictions on $\mathbf{L}_{ij} $, then $\mathbf{S}_{ij} = \text{N/A}$.
Recall that under Assumptions 1-2, Proposition (ref) dictates that the factor loadings matrix $\mathbf{L}^*$ of any observationally equivalent model must be of the form $\mathbf{L}^* = \mathbf{L}\mathbf{P}$, where $\mathbf{P}$ is a product of a reflection and a permutation. To be observationally equivalent with the sign restrictions imposed---i.e., satisfying exactly the same sign restrictions---we must have $\mathbf{S}\mathbf{P} = \mathbf{S}$. Intuitively then, for point-identification of $\mathbf{L}$ there must be enough structure in $\mathbf{S}$ such that the only possible $\mathbf{P}$ is the identify matrix. Now, suppose that each column of $\mathbf{S}$ has at least one sign restriction and no columns are the same or negative of any other columns. These conditions are sufficient as they rule out any permutations or sign changes except the identity.
To see that the conditions are necessary, suppose there is a column that has no sign restrictions. Then, changing the sign of the associated column in $\mathbf{L}$ (and the associated rows in $\mathbf{f}_t$) would leave the model observationally equivalent. Next, suppose one column is the same or the negative of any other column, then we can permute (and change signs if necessary) the relevant columns to leave the model observationally equivalent. We summarize these results in the following corollary.
In this section we describe an efficient posterior sampler to estimate the model in (ref)--(ref) with signs or zero restrictions specified in $\mathbf{S}$. Below we note a few details in our implementation with the goal of improving speed and sampling efficiency. First, even though the factors $\mathbf{f}_1,\ldots, \mathbf{f}_T$ are conditionally independent given the data and other parameters, we sample them jointly in one step using the precision sampler of CJ09---instead of drawing them sequentially in a for-loop---to speed up the computations.
Second, since VARs tend to have a lot of parameters even for small and medium systems, we implement an equation-by-equation estimation approach similar in spirit to that in CCM19 to sample the VAR coefficients. Specifically, given the latent factors $\mathbf{f}$, the VAR becomes $n$ unrelated regressions, and one can sample the VAR coefficients equation by equation without any loss of efficiency. Third, with the sign restrictions imposed in $\mathbf{S}$, the full conditional distribution of the factor loadings in each equation becomes a truncated multivariate normal distribution. To sample from such a distribution, we use the algorithm in botev17 that is based on quadratic programming.
For notational convenience, stack $\mathbf{y} = (\mathbf{y}_1',\ldots, \mathbf{y}_T')',$ $\mathbf{f} = (\mathbf{f}_1',\ldots, \mathbf{f}_T')',$ $\mathbf{h} = (\mathbf{h}_1',\ldots, \mathbf{h}_{T}')'$ and $\boldsymbol \beta = (\boldsymbol \beta_1',\ldots, \boldsymbol \beta_n')'$. In addition, let $\mathbf{y}_{i,\boldsymbol{\cdot}} = (y_{i,1},\ldots, y_{i,T})'$ denote the vector of observed values for the $i$-th variable, $i=1,\ldots, n$. We similarly define $\mathbf{h}_{i,\boldsymbol{\cdot}} = (h_{i,1},\ldots, h_{i,T})', i=1,\ldots, n+r$. Then, posterior draws can be obtained by sampling sequentially from:
Step 1. As mentioned above, since the factors $\mathbf{f}_1,\ldots, \mathbf{f}_T$ are conditionally independent given other parameters, in principle one can sample each factor sequentially in a for-loop. Here, however, we vectorize all the operations and sample them jointly in one step to improve computational speed. More specifically, we first stack the VAR in (ref)-(ref) over $t=1,\ldots, T$ and write it as: \[ \mathbf{y} = \mathbf{X}\boldsymbol \beta + (\mathbf{I}_T\otimes \mathbf{L})\mathbf{f} + \mathbf{u}^y, \quad \mathbf{u}^y \sim \mathcal{N}(\mathbf{0},\boldsymbol \Sigma), \] where $\mathbf{X}$ is the matrix of intercepts and lagged values and $\boldsymbol \Sigma = \text{diag}(\boldsymbol \Sigma_1,\ldots, \boldsymbol \Sigma_T)$ with $\boldsymbol \Sigma_t = \text{diag}(\text{e}^{\mathbf{h}_t^y})$. In addition, it follows from (ref) that $(\mathbf{f} \,|\, \mathbf{h})\sim \mathcal{N}(\mathbf{0},\boldsymbol \Omega),$ where $\boldsymbol \Omega = \text{diag}(\boldsymbol \Omega_1,\ldots, \boldsymbol \Omega_T)$ with $\boldsymbol \Omega_t = \text{diag}(\text{e}^{\mathbf{h}_t^f})$.
Then, by standard linear regression results CKPT19, we have
where
Note that the precision matrix $\mathbf{K}_{\mathbf{f}}$ is a band matrix, i.e., it is sparse and all the nonzero entries are arranged along the diagonal bands above and below the main diagonal. As such, once can use the precision sampler of CJ09 to sample $\mathbf{f}$ efficiently.
Step 2. Next, we sample $\boldsymbol \beta$ and $\mathbf{L}$ jointly to improve sampling efficiency. Given the latent factors $\mathbf{f}$, the VAR becomes $n$ unrelated regressions and we can sample $\boldsymbol \beta$ and $\mathbf{L}$ equation by equation. Recall that $\mathbf{y}_{i,\boldsymbol{\cdot}} = (y_{i,1},\ldots, y_{i,T})'$ is defined to be the $T\times 1$ vector of observations for the $i$-th variable; and that $\boldsymbol \beta_i$ and $\mathbf{l}_i$ represent, respectively, the VAR coefficients and the factor loadings in the $i$-th equation. Then, the $i$-th equation of the VAR can be written as \[ \mathbf{y}_{i,\boldsymbol{\cdot}} = \mathbf{X}_i\boldsymbol \beta_i + \mathbf{F} \mathbf{l}_i + \mathbf{u}_{i,\boldsymbol{\cdot}}^y, \] where $\mathbf{F} = (\mathbf{f}_{1,\boldsymbol{\cdot}},\ldots, \mathbf{f}_{r,\boldsymbol{\cdot}})$ is the $T\times r$ matrix of factors with $\mathbf{f}_{i,\boldsymbol{\cdot}} = (f_{i,1},\ldots, f_{i,T})'$. The vector of disturbances $\mathbf{u}_{i,\boldsymbol{\cdot}}^y= (u_{i,1},\ldots, u_{i,T})'$ is distributed as $\mathcal{N}(\mathbf{0},\boldsymbol \Omega_{\mathbf{h}_{i,\boldsymbol{\cdot}}})$, where $\boldsymbol \Omega_{\mathbf{h}_{i,\boldsymbol{\cdot}}} =\text{diag}(\text{e}^{h_{i,1}},\ldots, \text{e}^{h_{i,T}})$.\footnote{Note that zero restrictions on $ \mathbf{l}_i$ can be easily handled by redefining $\mathbf{l}_i $ and $ \mathbf{F}$ appropriately. For example, if the first element of $\mathbf{l}_i $ is restricted to be zero, we can define $\widetilde{\mathbf{l}}_i $ to be the vector consisting of the second to $r$-th elements of $\mathbf{l}_i$ and $\widetilde{\mathbf{F}} = (\mathbf{f}_{2,\boldsymbol{\cdot}},\ldots, \mathbf{f}_{r,\boldsymbol{\cdot}})$. Then, we replace $\mathbf{F} \mathbf{l}_i$ by $\widetilde{\mathbf{F}}\widetilde{\mathbf{l}}_i$.} Letting $\boldsymbol \theta_i = (\boldsymbol \beta_i',\mathbf{l}_i')'$ and $\mathbf{Z}_i = (\mathbf{X}_i,\mathbf{F})$, we can further write the VAR systems as \[ \mathbf{y}_{i,\boldsymbol{\cdot}} =\mathbf{Z}_i\boldsymbol \theta_i + \mathbf{u}_{i,\boldsymbol{\cdot}}^y. \] Let $R_i\subset \mathbb{R}^r$ be the support of $\mathbf{l}_i$ defined by the sign restrictions specified in the $i$-th row of $\mathbf{S}$. Then, using standard linear regression results, we obtain: \[ (\boldsymbol \theta_i \,|\, \mathbf{y}_{i,\boldsymbol{\cdot}}, \mathbf{f}, \mathbf{h}_{i,\boldsymbol{\cdot}}) \sim \mathcal{N}(\widehat{\boldsymbol \theta}_i,\mathbf{K}_{\boldsymbol \theta_i}^{-1})1(\mathbf{l}_i\in R_i), \] where \[ \mathbf{K}_{\boldsymbol \theta_i} = \mathbf{V}_{\boldsymbol \theta_i}^{-1} + \mathbf{Z}_i'\boldsymbol \Omega_{\mathbf{h}_{i,\boldsymbol{\cdot}}}^{-1}\mathbf{Z}_i, \quad \widehat{\boldsymbol \theta}_i = \mathbf{K}_{\boldsymbol \theta_i}^{-1}(\mathbf{V}_{\boldsymbol \theta_i}^{-1}\boldsymbol \theta_{0,i} + \mathbf{Z}_i\boldsymbol \Omega_{\mathbf{h}_{i,\boldsymbol{\cdot}}}^{-1} \mathbf{y}_{i,\boldsymbol{\cdot}}) \] with $\mathbf{V}_{\boldsymbol \theta_i} = \text{diag}(\mathbf{V}_{\boldsymbol \beta_i},\mathbf{V}_{\mathbf{l}_i})$ and $\boldsymbol \theta_{0,i} = (\boldsymbol \beta_{0,i}',\mathbf{l}_{0,i}')'$. A draw from this truncated multivariate normal distribution can be obtained using the algorithm in botev17. The remaining steps are standard and we leave the details to Appendix B. Some simulation results are reported in Appendix D to show that the posterior sampler works well and the posterior estimates track the true values closely.
It is worth noting that Proposition (ref) only guarantees that the factors and factor loadings are identified up to permutations and sign changes. Hence, in practice one might encounter the so-called label-switching problem. One way to handle this issue is to post-process the posterior draws to sort them into the correct categories; see, e.g., KS19 for such an approach. In our empirical application we impose sign restrictions that satisfy Corollary (ref)---and consequently the factors and factor loadings are point-identified.
Next, we document the runtimes of estimating the VAR-FSV of different dimensions to assess how well the posterior sampler scales to higher dimensions. More specifically, we report in Table (ref) the computation times (in minutes) to obtain 10,000 posterior draws from the VAR-FSV of dimensions $n= 15, 30, 50$ and sample sizes $T=300, 800$. The posterior sampler is implemented in $\mathrm{M}\mathrm{{\scriptstyle ATLAB}}$ on a typical desktop with an Intel Core i5-9600 @3.10 GHz processor and 16 GB memory. It is evident from the table that even for high-dimensional applications with 50 variables, the VAR-FSV model with sign restrictions imposed on the factor loadings can be estimated fairly quickly.
This section first gives a brief overview on the theory of Bayesian model comparison via the marginal likelihood. Then, we introduce an algorithm to evaluate the likelihood, or more precisely the integrated likelihood marginal of the latent states, implied by the VAR-FSV model. Finally, we present an adaptive importance sampling algorithm to estimate the marginal likelihood under the VAR-FSV model.
Suppose we wish to compare a collection of models $\{M_{1},\ldots, M_{K} \}$, where each model $M_k$ is defined by a likelihood function $p(\mathbf{y}\,|\, \boldsymbol \theta_k, M_{k})$ and a prior on the model-specific parameter vector $\boldsymbol \theta_k$ denoted by $p(\boldsymbol \theta_k \,|\, M_k)$. The gold standard for Bayesian model comparison is the Bayes factor in favor of $M_i$ against $M_j$, defined as \[ \text{BF}_{ij} = \frac{p(\mathbf{y}\,|\, M_i)}{p(\mathbf{y}\,|\, M_j)}, \] where
is the {\em marginal likelihood\/} under model $M_k$, $k=i,j.$ This Bayes factor is related to the posterior odds ratio between the two models: \[ \frac{\mathbb P(M_i\,|\,\mathbf{y})}{\mathbb P(M_j\,|\,\mathbf{y})} = \frac{\mathbb P(M_i)}{\mathbb P(M_j)}\times \text{BF}_{ij}, \] where $\mathbb P(M_i)/\mathbb P(M_j)$ is the prior odds ratio. It if clear that if both models are equally probable a priori, i.e., $p(M_i) = p(M_j)$, then the posterior odds ratio between the two models is equal to the Bayes factor. Hence, the Bayes factor has a natural interpretation and is easy to understand. For example, under equal prior odds, if $\text{BF}_{ij} = 50$, then model $M_i$ is 50 times more likely than model $M_j$ given the data. For a more detailed discussion of the Bayes factor and its role in Bayesian model comparison, see koop03 or CKPT19. From here onwards we suppress the model indicator.
To estimate the marginal likelihood, we first present an efficient way to evaluate the likelihood, or more precisely the integrated likelihood marginal of the latent states, given in (ref) . For notational convenience, we rewrite (ref) as
where the conditional density of $\mathbf{y}$ given $\mathbf{h}$ but marginal of $\mathbf{f}$ has the explicit expression \[ p(\mathbf{y}\,|\,\boldsymbol \beta,\mathbf{L},\mathbf{h}) = (2\pi)^{-\frac{Tn}{2}}\prod_{t=1}^T|\mathbf{L}\boldsymbol \Omega_t\mathbf{L}'+\boldsymbol \Sigma_t|^{- \frac{1}{2}} \text{e}^{-\frac{1}{2}(\mathbf{y}_t - (\mathbf{I}_n\otimes\mathbf{x}_t')\boldsymbol \beta)'(\mathbf{L}\boldsymbol \Omega_t\mathbf{L}'+\boldsymbol \Sigma_t)^{-1} (\mathbf{y}_t - (\mathbf{I}_n\otimes\mathbf{x}_t')\boldsymbol \beta)}. \] The second term of the integrand, $p(\mathbf{h} \,|\, \boldsymbol \mu,\boldsymbol \phi,\boldsymbol \sigma^2)$, is a $T(n+r)$-variate Gaussian density implied by the state equations specified in (ref)-(ref). Its analytical expression is given in Appendix C, and in particular, its precision matrix is banded. Hence, both densities can be evaluated quickly.
The main difficulty in evaluating the integrated likelihood in (ref), however, is that it requires integrating out all the latent log-volatilities, which involves solving a $T(n+r)$-dimensional integral. In what follows, we adopt the importance sampling approach developed for time-varying parameter VARs in CE18 to our setting.\footnote{There is a long tradition of using importance sampling to evaluate the integrated likelihood of stochastic volatility models. Earlier papers, such as DK97, SP97, KH02, FW08, McCausland12, have focused mostly on univariate stochastic volatility models.} More specifically, given an importance sampling density $g$---that might depend on model parameters and the data---we evaluate the integrated likelihood via importance sampling:
where $\mathbf{h}^{(1)},\ldots, \mathbf{h}^{(R_1)}$ are independent draws from $g$.
The choice of the importance sampling density $g$ is vital as it determines the variance of the estimator. In general, we would like to use an importance sampling density so that it well approximates the integrand in (ref). Our particular choice is motivated by the observation that there is, in fact, a theoretical zero-variance importance sampling density---it is $p(\mathbf{h}\,|\,\mathbf{y},\boldsymbol \beta,\mathbf{L},\boldsymbol \mu,\boldsymbol \phi,\boldsymbol \sigma^2)$, the conditional posterior distribution of $\mathbf{h}$ given the other parameters but marginal of $\mathbf{f}$. In practice, however, this density cannot be used as an importance sampling density as it is non-standard (e.g., its normalizing constant is unknown and it is unclear how one can efficiently generate samples from this density). But this observation provides us guidance for selecting a good importance sampling density.
In particular, we aim to approximate this ideal importance sampling density using a Gaussian density. This is accomplished as follows. We first develop an expectation-maximization (EM) algorithm to locate the mode of $\log p(\mathbf{h}\,|\,\mathbf{y},\boldsymbol \beta,\mathbf{L},\boldsymbol \mu,\boldsymbol \phi,\boldsymbol \sigma^2)$, denoted as $\widehat{\mathbf{h}}$. Then, we obtain the negative Hessian of this log-density evaluated at the mode, denoted as $\mathbf{K}_{\mathbf{h}}$. The mode and the negative Hessian are then used, respectively, as the mean vector and precision matrix of the Gaussian approximation. That is, the importance sampling density is $\mathcal{N}(\widehat{\mathbf{h}},\mathbf{K}_{\mathbf{h}}^{-1}).$ We leave the technical details to Appendix C. Below we comment on a few computational details.
First, in the M-step of the EM algorithm, one needs to solve a $T(n+r)$-dimensional maximization problem, which is in general extremely computationally intensive. In our case, however, we are able to obtain analytical expressions of the gradient and the Hessian of the objective function (i.e., the Q-function), which allows us to implement the Newton-Raphson method. Furthermore, one can show that the Hessian is a) negative definite anywhere in $\mathbb{R}^{T(n+r)}$, and b) a band matrix. The former property guarantees rapid convergence of the Newton-Raphson method, while the latter property substantially speeds up the computations.
Second, to construct the importance sampling estimator in (ref), one needs to both evaluate and sample from the $T(n+r)$-dimensional Gaussian importance sampling density $M$ times. For very high-dimensional Gaussian densities, both operations are generally computational costly. For our Gaussian importance sampling density, however, we can show that its precision matrix is banded. As such, samples from this Gaussian density can be obtained quickly using the precision sampler in CJ09. Evaluation of the density can be done just as quickly. We summarize the evaluation of the integrated likelihood in Algorithm (ref).
Next, we discuss the marginal likelihood estimation of the VAR-FSV model using an adaptive importance sampling approach called the improved cross-entropy method. This method requires little explicit analysis from the user and is applicable to a wide variety of problems (in contrast to the importance sampling estimator of the integrated likelihood estimation presented in Algorithm (ref) that requires a lot of analysis). More specifically, suppose we wish to estimate the marginal likelihood $p(\mathbf{y})\equiv p(\mathbf{y}\,|\, M_k)$ given in (ref) using the following importance sampling estimator:
where $\boldsymbol \theta^{(1)},\ldots, \boldsymbol \theta^{(R_2)}$ are independent draws from the importance sampling density $g(\cdot)$. In particular, for our FSV model, $\boldsymbol \theta = \{\boldsymbol \beta,\mathbf{L},\boldsymbol \sigma^2,\boldsymbol \mu,\boldsymbol \phi\}$. While this importance sampling estimator in theory is unbiased and simulation consistent for any $g$---as long as it dominates $p(\mathbf{y}\,|\,\cdot)p(\cdot)$, i.e., $g(\boldsymbol \theta)=0\Rightarrow p(\mathbf{y}\,|\,\boldsymbol \theta)p(\boldsymbol \theta)=0$---in practice its performance heavily depends on the choice of $g$. Here we follow CE15 to use the improved cross-entropy method to construct $g$ optimally.\footnote{The original cross-entropy method was developed for rare-event simulation by rubinstein97, rubinstein99 using a multi-level procedure to construct the optimal importance sampling density. CE15 later show that this optimal importance sampling density can be obtained more accurately in one step using Markov chain Monte Carlo methods.}
To motivate the improved cross-entropy method, first note that the ideal zero-variance importance sampling density is the posterior density $p(\boldsymbol \theta \,|\,\mathbf{y})$. That is, if we use $g^*(\boldsymbol \theta) = p(\boldsymbol \theta\,|\,\mathbf{y}) = p(\mathbf{y}\,|\,\boldsymbol \theta)p(\boldsymbol \theta)/p(\mathbf{y})$ as the importance sampling density, then the associated estimator in (ref) has zero variance. Unfortunately, $g^*$ cannot be used in practice as its normalization constant is precisely the marginal likelihood, the unknown quantity we aim to estimate. This nevertheless provides a benchmark to construct an optimal importance sampling density. More specifically, we aim to find a density that is `close' to this benchmark $g^*$ that can be used as an importance sampling density.
To that end, consider a parametric family $\mathcal{G} = \{ g(\boldsymbol \theta;\mathbf{v}) \}$ indexed by the parameter vector $\mathbf{v}$. We then find the density $g(\boldsymbol \theta;\mathbf{v}^*)\in\mathcal{G}$ such that it is, in a well-defined sense, the `closest' to $g^*$. One convenient measure of closeness between densities is the Kullback-Leibler divergence or the cross-entropy distance. More precisely, for two density functions $g_1$ and $g_2$, the cross-entropy distance from $g_1$ to $g_2$ is defined as: \[ \mathcal{D}(g_1,g_2) = \int g_1(\mathbf{x})\log \frac{g_1(\mathbf{x})}{g_2(\mathbf{x})}\text{d}\mathbf{x}. \] Given this measure of closeness, we obtain the density $g(\cdot;\mathbf{v})\in\mathcal{G}$ such that $\mathcal{D}(g^*,g(\cdot;\mathbf{v}))$ is minimized, i.e., $\mathbf{v}_{\text{ce}}^* = \mathop{\rm argmin}_{\mathbf{v}}\mathcal{D}(g^*,g(\cdot;\mathbf{v}))$. It can be shown that solving the CE minimization problem is equivalent to finding \[ \mathbf{v}^*_{\text{ce}} = \mathop{\rm argmax}_{\mathbf{v}}\int p(\mathbf{y}\,|\,\boldsymbol \theta)p(\boldsymbol \theta)\log g(\boldsymbol \theta;\mathbf{v})\text{d}\boldsymbol \theta. \]
In general this optimization problem is difficult to solve analytically as it involves a high-dimensional integral. Instead, we consider its stochastic counterpart:
where $\boldsymbol \theta_1,\ldots, \boldsymbol \theta_M$ are posterior draws from $p(\boldsymbol \theta\,|\,\mathbf{y}) \propto p(\mathbf{y}\,|\,\boldsymbol \theta)p(\boldsymbol \theta)$. It is useful to note that $\widehat{\mathbf{v}}^*_{\text{ce}}$ is exactly the maximum likelihood estimate for $\mathbf{v}$ if we view $g(\boldsymbol \theta;\mathbf{v})$ as the likelihood function with parameter vector $\mathbf{v}$ and $\boldsymbol \theta_1,\ldots, \boldsymbol \theta_M$ as an observed sample. Since finding the maximum likelihood estimator is a standard problem, solving (ref) is typically easy. For example, analytical solutions to (ref) are available for the exponential family rk:ce.
Next, we discuss the choice of the parametric family $\mathcal{G}$. One convenient class of densities is one in which each member $g(\boldsymbol \theta ; \mathbf{v})$ is a product of probability densities, e.g., $g(\boldsymbol \theta; \mathbf{v}) = g(\boldsymbol \theta_1; \mathbf{v}_1)\times\cdots\times g(\boldsymbol \theta_B; \mathbf{v}_B)$, where $\boldsymbol \theta = \{\boldsymbol \theta_1,\ldots, \boldsymbol \theta_B\}$ and $\mathbf{v} = \{\mathbf{v}_1,\ldots,\mathbf{v}_B\}$. One main advantage of this choice is that we can then reduce the generally high-dimensional maximization problem (ref) into $B$ separate low-dimensional maximization problems. For example, for the our FSV model, we divide $\boldsymbol \theta = \{\boldsymbol \beta,\mathbf{L},\boldsymbol \sigma^2,\boldsymbol \mu,\boldsymbol \phi\}$ into 5 natural blocks, and consider the parametric family \[
\] where $g_{\mathcal{N}}$ and $g_{\mathcal{IG}}$ are, respectively, the Gaussian and the inverse-gamma densities. Given this choice of the parametric family, the maximization problem in (ref) can be readily solved (either analytically or using numerical optimization).
Given the optimal importance density, denoted as $g(\boldsymbol \beta,\mathbf{L},\boldsymbol \sigma^2,\boldsymbol \mu,\boldsymbol \phi; \mathbf{v}^*)$, we construct the following importance sampling estimator:
where $(\boldsymbol \beta^{(1)},\mathbf{L}^{(1)},\boldsymbol \sigma^{2(1)},\boldsymbol \mu^{(1)},\boldsymbol \phi^{(1)}), \ldots, (\boldsymbol \beta^{(R_2)},\mathbf{L}^{(R_2)},\boldsymbol \sigma^{2(R_2)},\boldsymbol \mu^{(R_2)},\boldsymbol \phi^{(R_2)})$ are independent draws from $g(\boldsymbol \beta,\mathbf{L},\boldsymbol \sigma^2,\boldsymbol \mu,\boldsymbol \phi; \mathbf{v}^*)$ and $p(\mathbf{y}\,|\,\boldsymbol \beta,\mathbf{L},\boldsymbol \sigma^2,\boldsymbol \mu,\boldsymbol \phi)$ is the integrated likelihood, which can be estimated using the estimator in (ref). We refer the readers to CE15 for a more thorough discussion of this adaptive importance sampling approach. We summarize the algorithm in Algorithm (ref).
Note that Algorithm (ref) has two nested importance sampling steps, and it falls within the importance sampling squared (IS$^2$) framework in \citet*{TSPK14}. We follow their recommendation to set $R_1$, the simulation size of the inner importance sampling loop (the importance sampling step for estimating the integrated likelihood), adaptively so that the variance of the log integrated likelihood is around 1. See also the discussion in \citet*{PdGK12}.
The VAR-FSV in (ref)-(ref) can be used to draw structural inference by employing standard tools such as impulse response functions, forecast error variance decompositions and historical decompositions. In particular, letting $\mathbf{A}(L) = \mathbf{I}_n - \mathbf{A}_1 L - \dotsm - \mathbf{A}_p L^p$, where $L$ is the lag operator, the representation
where $\boldsymbol \Phi(L) = \mathbf{A}(L)^{-1}$ is well-defined assuming $\det \mathbf{A}(z) \ne 0$ for all $|z| < 1, z\in \mathbb{C}$ (i.e., the process $\{\mathbf{y}_t : t \in \mathbb{Z}\}$ is covariance-stationary).
Although $\boldsymbol \varepsilon_t$ does not contain structural shocks (since its elements are correlated), the reduced-form representation (ref) can be matched to a hypothetical structural representation of the form
where $\mathbf{u}_t$ is a vector of structural shock, and hence, its elements are uncorrelated. Note that $\widetilde{\boldsymbol \Phi}_t(L)$ is time-varying because $\text{Var}(\boldsymbol \varepsilon_t) = \mathbf{L}\boldsymbol \Omega_t\mathbf{L}' + \boldsymbol \Sigma_t$ is time-varying; consequently, hypothetical structural representations that can be matched to (ref) will generally have time-varying parameters.
The standard structural VAR approach is to assume that (i) $\mathbf{u}_t$ is $n \times 1$ and (ii) $\boldsymbol \varepsilon_t = \widetilde{\boldsymbol \Phi}_{0,t} \mathbf{u}_t$, where $\widetilde{\boldsymbol \Phi}_{0,t}$ is a $n\times n$ constant matrix with $\operatorname{rank} \widetilde{\boldsymbol \Phi}_{0,t} = n$ and $\text{Var}(\mathbf{u}_t) = \mathbf{I}_n$. Then, $\widetilde{\boldsymbol \Phi}_t(L) = \boldsymbol \Phi(L)\widetilde{\boldsymbol \Phi}_{0,t}$ and $\widetilde{\boldsymbol \Phi}_{0,t}$ satisfies \[ \widetilde{\boldsymbol \Phi}_{0,t}\widetilde{\boldsymbol \Phi}_{0,t}' = \mathbf{L}\boldsymbol \Omega_t\mathbf{L}' + \boldsymbol \Sigma_t. \] In this case, identification of $\widetilde{\boldsymbol \Phi}_{0,t}$ requires additional restrictions since \[ \check{\boldsymbol \Phi}_{0,t}\check{\boldsymbol \Phi}_{0,t}' = \mathbf{L}\boldsymbol \Omega_t\mathbf{L}' + \boldsymbol \Sigma_t \] for all $\check{\boldsymbol \Phi}_{0,t} = \widetilde{\boldsymbol \Phi}_{0,t}\mathbf{R}_t$, given an arbitrary orthonormal matrix $\mathbf{R}_t$ (i.e. satisfying $\mathbf{R}_t'\mathbf{R}_t = \mathbf{R}_t\mathbf{R}_t' = \mathbf{I}_n$).
An alternative way to obtain structural inference in our settings---similar to Korobilis20---is to assume that $\widetilde{\boldsymbol \Phi}_t(L)$ is $n\times(r+n)$ and $\mathbf{u}_t$ is $(r+n)\times 1$, such that
where $\widetilde{\mathbf{f}}_t = \boldsymbol \Omega_t^{-\frac{1}{2}}\mathbf{f}_t$ and $\widetilde{\mathbf{u}}_t^y = \boldsymbol \Sigma_t^{-\frac{1}{2}} \mathbf{u}_t^y$. In this case, $\widetilde{\boldsymbol \Phi}_t(L) = \boldsymbol \Phi(L)\widetilde{\boldsymbol \Phi}_{0,t}$ is also $n \times (r+n)$, which in departure from standard structural VARs leads to a `short' system FGS19.
Identification of impulse response functions and forecast error variance decompositions in short systems is generally problematic PR22,CF22. However, in the formulation above, $\widetilde{\boldsymbol \Phi}_{0,t}$ is identified due to $\boldsymbol \Sigma_t$ being restricted to a diagonal matrix and $\mathbf{L}$ being identified by sign restrictions as described in Section (ref). Hence, impulse response functions and forecast error variance decompositions to all shocks in $\mathbf{u}_t$ are identified, even though $\mathbf{u}_t$ is generally {\em not recoverable\/} from past and future observations of $\mathbf{y}_t$, as defined in CJ21.
In our setting, the main interest lies in quantifying the effects of shocks in $\widetilde{\mathbf{f}}_t$, and therefore, sign restrictions on $\mathbf{L}$ play the role of endowing these shocks with economic meaning. The remaining elements in $\widetilde{\mathbf{u}}_t^y$ are not of direct interest and we do not treat them as economically meaningful shocks. Nevertheless, the restriction that $\boldsymbol \Sigma_t$ is diagonal plays a crucial role in the overall identification strategy along with an economically meaningful interpretation of $\widetilde{\mathbf{f}}_t$. We provide explicit expressions for computing impulse response functions and forecast error variance decompositions in Appendix F.
Finally, computing historical decompositions requires $\mathbf{f}_t$ and $\mathbf{u}_t^y$ (see Appendix F for details). The fact that $\mathbf{u}_t$ is not recoverable implies that historical decompositions may not be point identified. In a Bayesian setting, however, the posterior distribution of a historical decomposition at each horizon may still be constructed using draws from the posterior distribution of the VAR-FSV parameters together with draws of $\mathbf{f}_t$.
In the algorithm developed in Section (ref), draws of $\mathbf{f}_t$ are a by-product of simulation, while $\mathbf{u}_t^y$ is easily obtained as \[ \mathbf{u}_t^y = \mathbf{A}(L)\mathbf{y}_t - \boldsymbol \mu - \mathbf{L}\mathbf{f}_t, \] for each draw of $\boldsymbol \mu, \mathbf{A}_1, \dotsc, \mathbf{A}_p, \mathbf{L}, \mathbf{f}_t$. Therefore, draws from the posterior distribution of a HD are straightforward to compute.
In addition, $\mathbf{u}_t$ can be regarded as being recoverable in the limit as $n\longrightarrow \infty$ from the VAR residual $\boldsymbol \varepsilon_t$ (and therefore past and present $\mathbf{y}_t$) under a suitable assumption on the factor loadings $\mathbf{L}$. To see this, let $\mathbf{L}^+$ denote the Moore–Penrose inverse of $\mathbf{L}$. By Assumption (ref), $\operatorname{rank} \mathbf{L} = r$ and $\mathbf{L}^+ = (\mathbf{L}'\mathbf{L})^{-1}\mathbf{L}'$. It also follows that a right inverse (although not a Moore-Penrose inverse) of $\widetilde{\boldsymbol \Phi}_{0,t}$ is
In the factor model literature, a standard assumption bai2003,FGLR09 is that $n^{-1}\mathbf{L}'\mathbf{L} \longrightarrow \boldsymbol \Lambda$ as $n \longrightarrow \infty$, where $\boldsymbol \Lambda$ is a constant (strictly) positive-definite matrix. It implies that the factors $\mathbf{f}_t$ are {\em pervasive\/} in the sense that they significantly affect most of the variables {\em on impact\/}.\footnote{It is worth emphasising, however, that this {\em does not\/} imply $\mathbf{u}_t^y$ is a vector of {\em idiosyncratic\/} errors, as defined by FHLR00,FL01 in the context of generalised dynamic factor models. In particular, the overall effect of $\mathbf{u}_t^y$ on $\mathbf{y}_t$ is $\boldsymbol \Phi(L)\mathbf{u}_t^y$, which is generally pervasive, albeit with a delay.} An immediate consequence of the pervasiveness assumption, together the regularity condition that $\text{Var}(\varepsilon_{i,t}) < \infty$ for all $i = 1, \dotsc, n$, is $\mathbf{L}^+\boldsymbol \Sigma_t^\frac{1}{2}\widetilde{\mathbf{u}}_t^y \overset{m.s.}{\longrightarrow} \mathbf{0}$. Combining this result with (ref) yields
Consequently, $\mathbf{u}_t$ is recoverable from $\boldsymbol \varepsilon_t$ in the limit.\footnote{A more general result on recoverability with a fixed $n$ is given in CJ21.}
In this section we conduct a series of simulation experiments to assess the adequacy of using the proposed marginal likelihood estimator to determine the number of factors. More specifically, we generate datasets from the VAR-FSV in (ref)--(ref), but we change the error structure to $\boldsymbol \varepsilon_t = \mathbf{L} \mathbf{f}_t + \sqrt{\theta} \mathbf{u}_t^y$, where $\theta$ measures the signal-to-noise ratio, following bai2002determining. We set parameter values so that if $\theta=r$, the idiosyncratic component will then have the same variance as the common component. In particular, we generate $L_{ij} \sim \mathcal{N}(0, 1)$ for $i=1,\ldots,n$ and $j=1,\ldots,r$ and set $\mu_i=0$ for $i=1,\ldots,n,$ so that the log-volatility processes associated with the idiosyncratic errors have 0 unconditional mean.
The remaining parameters are generated as follows. The intercepts are drawn independently from the uniform distribution on the interval $(-10,10)$, i.e., $\mathcal{U}(-10, 10)$. For the VAR coefficients, the diagonal elements of the first VAR coefficient matrix are iid $\mathcal{U}(0,0.5)$ and the off-diagonal elements are from $\mathcal{U}(-0.2,0.2)$; all other elements of the $j$-th ($j > 1$) VAR coefficient matrices are iid $\mathcal{N}(0,0.1^2/j^2).$ Finally, for the log-volatility processes, we set $\phi_i=0.98$ and $\sigma_i^2 = 0.1^2$ for $i=1,\ldots,n+r$.
In this Monte Carlo study, we select the true number of factors $r \in \{1,3,5\}$ and $\theta \in \{1,3,5,10\}$; and we consider the number of variables $n \in \{15, 30\}$ and sample size $T \in \{300,500,800\}$. For each set of $(r,\theta, n, T)$, we generate 100 datasets. For each dataset, we estimate the VAR-FSV models with $r=1,\ldots, 6$ factors and compute the associated marginal likelihood values. For this Monte Carlo experiment, a total of 14,400 separate MCMCs and marginal likelihood estimation are run (24 settings $\times$ 6 factor models $\times$ 100 datasets). Among the 6 factor models for each dataset and parameter setting, we select the one with the largest marginal likelihood value. Table (ref) reports the selection frequency.
The Monte Carlo results show that the marginal likelihood estimator generally performs well in selecting the correct number of factors under a variety of settings. For example, for $n=15,$ $T=300,$ $r = 3$ and $\theta = 3$ (moderate signal-to-noise ratio), the marginal likelihood estimator is able to pick the correct number of factors 83% of the times; for the rest of the cases, the model with one fewer factor (13%) or one more factor (4%) is selected. In addition, as the sample size $T$ increases to 500, the selection frequency of the $r=3$ factor model increases to 97%. More generally, the selection frequency of the correct number of factors increases as the sample size $T$ increases for all cases considered, confirming that the marginal likelihood is a consistent model selection criterion.
We illustrate the proposed methodology by revisiting the structural analysis in FRS19 that is based on a standard structural VAR. More specifically, they use a 6-variable structural VAR to study the impacts of 5 structural shocks---demand, supply, monetary, investment and financial shocks---on a number of key economic variables, where these structural shocks are identified using sign restrictions on the contemporaneous impact matrix. The size of the VAR in their application is typical among empirical works that use sign restrictions for identification because of the computational burden.\footnote{For their 6-variable structural VAR, FRS19 report estimation time of about a week using a 12-core workstation.}
However, there are a number of reasons in favor of using a larger set of macroeconomic and financial variables. First, in practice the mapping from variables in an economic model to the data is often not unique. For example, as argued in LMW21, the economic variable inflation could be matched to data based on the CPI, PCE, or the GDP deflator, and it is not obvious which time series should be used. Instead of arbitrarily choosing one inflation measure, it is more appropriate to include multiple time series corresponding to the same economic variable in the analysis.
Second, one might be concerned about the problem of informational deficiency that arises from using a limited information set. More specifically, influential papers such as HS91 and LR93, LR94 have pointed out that when the econometrician considers a narrower set of variables than the economic agent, the underlying model used by the econometrician is non-fundamental. That is, current and past observations of the variables do not span the same space spanned by the structural shocks. As a consequence, structural shocks cannot be recovered from the model. A natural way to alleviate this concern of informational deficiency is to use a larger set of relevant variables gambetti21.
In view of these considerations, we augment the 6-variable VAR with additional macroeconomic and financial variables, and consider a 20-variable VAR with factor stochastic volatility identified using sign restrictions. There are two related papers that use large VARs to study the role of financial shocks in economic fluctuations. First, chan21 considers a 15-variable structural VAR with a new asymmetric conjugate prior to identify the financial shocks. Given the larger system and the large number of sign restrictions, estimation time is about a week to obtain 1,000 admissible draws using the algorithm of RWZ10. In contrast, the proposed approach takes less than a minute to obtain the same number of admissible draws, and is applicable to even larger systems. Second, Korobilis20 uses a 15-variable VAR with a factor error structure to identify the financial shocks, which can also be done quickly. The main advantage of our approach, however, is that the structural shocks obtained using our factor stochastic volatility model are point-identified, whereas they are only set-identified under a homoskedastic VAR. In practice, our approach can often provide sharper inference.
We use a dataset that consists of 20 US quarterly variables, which are constructed from raw time-series taken from from different sources, including the Federal Reserve Bank of Philadelphia and the FRED database at the Federal Reserve Bank of St. Louis. For easy comparison with the results in FRS19, we use the same sample period that spans from 1985:Q1 to 2013:Q2. The complete list of these time-series and their sources are given in Appendix E.
We include the same 6 variables used in the baseline model in FRS19, namely, real GDP, GDP deflator, 3-month treasury rate, ratio of private investment over output, S&P 500 index and a credit spread defined as the difference between Moody's baa corporate bond yield and the federal funds rate. In addition, we also include 14 additional macroeconomic and financial variables, such as the ratio of total credit over real estate value, labor market variables, mortgage rates, as well as other measures of inflation, interest rates and stock prices. These 20 variables are listed in Table (ref) and the details of the raw data are given in Appendix E.
In this section we re-examine the empirical application in FRS19 that identifies 5 structural shocks using a structural VAR with sign restrictions on the contemporaneous impact matrix. We first use the proposed VAR-FSV model to replicate their baseline results from a 6-variable VAR, but here we impose the sign restrictions on the factor loadings instead of the impact matrix. We then consider a larger VAR-FSV model with 20 variables to identify the same structural shocks.
Now, we first employ the same 6 variables and the associated sign restrictions used in FRS19, which are presented in the first six rows of Table (ref). The sign restrictions to identify the supply, demand, monetary, investment and financial shocks are exactly the same as in FRS19, and we refer the readers to their paper for more details. Here we only note that in order to distinguish investment and financial shocks from demand shocks, they are assumed to have different effects on the ratio of investment over output. In particular, positive investment and financial shocks have a positive effect on the ratio, motivating by the idea that investment and financial shocks create investment booms. By contrast, positive demand shocks reduce the ratio of investment over output, i.e., even though investment level could increase in response to demand shocks, it does not increase as much as other components of output.
We compute the impulse responses from the VAR-FSV with 5 factors, where the sign restrictions are imposed on the factor loadings. Since FRS19 use an improper/non-informative prior in their analysis, to make our results comparable, we consider a proper but relatively vague prior by setting $\kappa_1=\kappa_2=1$.\footnote{The variables are expressed in level. As such, the prior means of the first own lags are all set to be 1, whereas those of other VAR coefficients are set to be 0. In addition, the prior mean of $\mu_i$, the mean of the idiosyncratic log-volatility for the $i$-th variable, is set to be $\log(0.1\times \text{Var}(\mathbf{y}_{i,\cdot}))$. That is, a priori about 10% of the sample variance is attributed to idiosyncratic component. Finally, the prior variance $V_{\mu_{i}}$ is set to be 10 for $i=1,\ldots, n$.} We use the Gibbs sampler described in Section (ref) to obtain 50,000 posterior, storing every 10-th draw, after a burn-in period of 5,000.
Figure (ref) plots the impulse responses of the 6 variables to an one-standard-deviation financial shock. Despite the differences in methodology, the impulse responses are very similar to those given in FRS19. Consistent with the findings in FRS19, the results show that financial shocks have a substantial impact on output, stock prices and investment, but have a limited impact on inflation (measured by GDP deflator). Furthermore, even though the impact on the spread is unrestricted, we find that its reaction to financial shocks is significantly counter-cyclical. These results highlight one advantage of the proposed methodology: the median impulse responses from the VAR-FSV are very similar to those obtained using a standard structural VAR, but instead of using an accept-reject algorithm to obtain admissible draws, the sign restrictions can be easily incorporated in the estimation, and consequently, it can be done much faster.
Next, we augment the 6-variable VAR with 14 additional macroeconomic and financial variables. Many of these new variables are alternative data series corresponding to the same economic variable. For example, in addition to GDP deflator as prices, we also include CPI and PCE index as alternative measures. Similarly, Dow Jones Industrial Average and NASDAQ indexes are added as alternative measures of stock prices. Furthermore, other seemingly relevant variables, such as labor market and national accounts variables, are also included to alleviate the concern of informational deficiency. The additional variables and the corresponding sign restrictions are listed in rows 7-20 of Table (ref).
With $n=20$ variables and $r=5$ factors, this large VAR-FSV satisfies the condition that $r \leqslant (n-1)/2$. In addition, it is easy to verify that the sign restrictions given in Table (ref) satisfy the conditions in Corollary 2, and therefore the latent factors, which we interpret as structural shocks, are point-identified. Given the large number of variables, it is crucial to apply proper shrinkage on the VAR coefficients. Following the bulk of the literature CCM15, we set $\kappa_1 = 0.04$ and $\kappa_2 = 0.04^2$---i.e., the VAR coefficients associated with lags of other variables are shrunk more strongly to 0 than those on own lags. Again we obtain 50,000 posterior after a burn-in period of 5,000 to compute the impulse responses. The results are reported in Figure (ref).
The impulse responses from this 20-variable VAR-FSV are qualitatively similar to those from the smaller system, but the inference is much sharper. Specifically, the credible intervals of the 6 impulse response functions are much narrower, highlighting the benefits of incorporating more relevant information---more variables and sign restrictions as well as a more informative prior---to sharpen inference. For example, the credible intervals associated with the responses of investment and stock prices exclude zero for the first 32 quarters after the initial impact of a financial shock. This is in contrast to the much wider credible intervals from the 6-variable VAR (the median impulse responses of stock prices even become negative at longer horizons). The results from this large system therefore better highlight the impact of a positive financial shock, which FRS19 define as “a shock that generates an investment and a stock market boom."
Next, we plot in Figure (ref) the median impulse responses of the 6 variables from the remaining 4 structural shocks. These impulse responses are similar to those presented in Figure 3 in FRS19. In particular, we confirm that supply shocks generate large effects not only on output, but also on investment and stock prices. On the other hand, demand shocks have smaller effects on output, investment and stock prices, at least for short to medium horizons, but they are the main driver of prices. Finally, while we also find that monetary shocks have a protracted positive effect on output, their effects on stock prices are more subdued.
To quantify how much of the historical fluctuations in GDP and spread can be attributed to each of the structural shocks, we compute the historical decompositions of these two variables using the formulas derived in Appendix F. The results are reported in Figure (ref) and Figure (ref).
These historical fluctuations from the VAR-FSV are in line with those obtained using a standard structural VAR presented in FRS19. In particular, financial shocks play a large role in explaining the historical fluctuations in both GDP and spread, especially in the lead-up and aftermath of the Great Recession of 2007-2009.
Next, we quantify the amount of the prediction mean squared errors of 6 selected variables accounted for by each of the 5 structural shocks at different forecast horizons. More specifically, using the expressions developed in Appendix F, we compute the forecast error variance decompositions of the variables and the results are presented in Table (ref).
Overall, financial shocks play a large role in explaining the forecast error variances of the majority of the variables. The two exceptions are prices (measured by GDP deflator), which are mainly impacted by demand shocks, and stock prices (measured by S&P 500 index), which are mostly driven by supply shocks.
We have considered an order-invariant VAR with factor stochastic volatility and shown how the presence of multivariate stochastic volatility allows for statistical identification of the model. Furthermore, we have worked out sufficient conditions in terms of sign restrictions on the impact of the structural shocks for point-identification of the corresponding structural model. To estimate the proposed order-invariant VAR, we developed an efficient MCMC algorithm that can incorporate a large number of variables and sign restrictions. In an empirical application involving 20 macroeconomic and financial variables, we demonstrated the ability of our methods to produce more precise impulse responses compared to a medium-sized structural VAR.