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.
76,995 characters · 14 sections · 63 citation commands
Large Bayesian VARs for Binary and Censored Variables
\doublespacing
\singlespacing{\footnotesizeContact: Michael Pfarrhofer ([email removed]), Department of Economics, WU Vienna University of Economics and Business.}
\thispagestyle{empty} \doublespacing
Since the seminal work of Sims80 and DLS84, vector autoregressions (VARs) have become the workhorse models for structural analysis and macroeconomic forecasting in both academia and policy institutions. Since most macroeconomic variables are continuous and can take any real value, the literature has focused on developing flexible models and associated tools for real-valued time series. However, many important macroeconomic and financial indicators have restricted ranges or can only take discrete values---e.g., nominal interest rates are bounded below by their effective lower bound (ELB); events like financial crises, recessions, or sovereign debt defaults are binary---which make standard VARs ill-suited for modeling these time series. The growing popularity of large datasets that typically include a variety of these censored and binary variables makes this issue more acute.
We develop an econometric framework that extends standard VARs to jointly model a mix of binary, censored, and continuous variables. This provides a unified framework to study, for example, the impact of some structural shock or the effect of a particular scenario on the recession probability, while explicitly taking into account the ELB on short-term interest rates. In addition, motivated by the trend of using large datasets in applied work banbura2010large, koop2013forecasting, the proposed approach is designed to be scalable to high dimensions. At the same time, it can also accommodate empirically important features such as time-varying volatility and robustness to COVID outliers.
More specifically, our approach is based on data augmentation, which treats binary and censored observations as a missing data problem (e.g., the value of a binary variable is determined by the sign of an associated real-valued latent variable).\footnote{Bayesian approaches for modeling binary and censored variables can be traced back to the seminal papers by chib1992bayes and AC93; see also anceschi2023bayesian for a review of related approaches. More recently, zhang2015bayesian propose a data augmentation approach to analyze a combination of continuous, ordinal and categorical variables. These papers mostly focus on cross-sectional data or static settings.} Given the real-valued latent variables, all the machinery developed for standard VARs can be directly applied. The proposed approach thus makes it possible to take advantage of all the recent advances in VAR modeling, such as the development of adaptive hierarchical shrinkage priors and efficient estimation algorithms for heteroskedastic VARs.
Although the additional step of simulating the latent variables is conceptually simple, it is practically challenging in our high-dimensional, dynamic setting for two reasons. First, the latent variables are serially dependent due to the VAR setup. As a result, the dilemma is that while sampling them sequentially would induce high autocorrelations and poor mixing, sampling all of them jointly using conventional methods is computationally infeasible. Second, the latent variables are subject to a large number of inequality restrictions by construction, making the sampling problem even more difficult. To tackle these challenges, we adapt the approach of chan2023high, which is designed for drawing missing data in conditionally Gaussian state space models, to our setting. In particular, we show that the joint conditional posterior density of the latent variables is truncated Gaussian, where its inverse covariance or precision matrix is banded (i.e., it is sparse and the non-zero elements are arranged along a main diagonal band), and this feature can be exploited to vastly speed up the computations. We then implement a Hamiltonian Monte Carlo (HMC) step to sample from this high-dimensional truncated Gaussian distribution.
Our paper contributes to the relatively sparse literature on modeling censored and discrete-valued time series. The influential paper by dueker2005dynamic develops the so-called qualitative VAR, or QualVAR, that extends a standard VAR to include a single binary variable. This model is used for macroeconomic forecasting and scenario analysis in dueker2006business, dueker2010forecasting and mccracken2022binary. poon2024recessions introduce a multinomial logit model for jointly modeling several binary indicators, with an application of forecasting concurrent recessions in the US, UK and Euro Area. For censored variables, johannsen2021time propose a multivariate unobserved components model for nominal interest rates that are bounded below by their ELB. Alternatively, carriero2023shadow extend a standard VAR by modeling interest rates as censored observations of latent variables. This paper contributes to this literature by providing a unified framework that can simultaneously handle multiple binary and censored time series, and developing efficient and scalable sampling algorithms.
We demonstrate the usefulness of the proposed framework in a variety of empirical applications. In an out-of-sample forecasting exercise, we use datasets of different sizes to forecast recessions and short-term interest rates. The forecasting results show that using larger datasets---including a wide range of financial and labor market variables beyond those typically considered in earlier research---improves recession forecasts. Moreover, models that accommodate heteroskedasticity tend to outperform their homoskedastic counterparts, especially for the post 2000 subsample. Lastly, the forecasting results highlight the value of treating interest rates as censored variables: VARs that do so forecast especially well for short-term interest rates compared to those that do not.
We then investigate the dynamic evolution of the model-implied latent state variables corresponding to the NBER recession indicator and the short-term interest rates. The first set of state variables may be interpreted as a latent indicator of the business cycle and its turning points. The state variables corresponding to the interest rates provide a measure of the overall stance of monetary policy when the ELB is binding. In the next exercise, we demonstrate how one can use conditional forecasting to assess recession probabilities of various scenarios. We consider scenarios based on the 2020 Federal Reserve (Fed) Board stress test and the December 2023 economic projections from the Fed Board. Finally, we perform a structural analysis to investigate the effects of a financial shock on the recession probability, identified by recursive restrictions, while explicitly taking into account the ELB.
The rest of this article is organized as follows. Section (ref) introduces the econometric framework and develops the estimation approach. Section (ref) conducts an out-of-sample exercise to forecast recessions and short-term interest rates using a variety of models and datasets of different sizes. Section (ref) presents several applications to illustrate how the proposed framework can be used in different settings, including conditional forecasting and structural analysis. Finally, Section (ref) concludes and outlines some possible future research directions.
Many macroeconomic time-series have restricted supports that prevent the direct application of standard modeling tools such as VARs and the associated apparatus for forecasting and structural analysis. We aim to develop an econometric framework that is suitable for modeling a mix of binary, censored, and unrestricted variables, while allowing us to leverage all the machinery available for standard VARs. In particular, the proposed framework enables taking advantage of recent advances in adaptive hierarchical shrinkage priors and efficient estimation algorithms designed for high-dimensional VARs.
To set the stage, let $\tilde{\bm{y}}_t = (\tilde{\bm{y}}_{b,t}',\tilde{\bm{y}}_{c,t}',\bm{y}_{u,t}')'$, for $t=1,\hdots,T,$ be an $n \times 1$ vector, with $n = n_b + n_c + n_u$, comprised of:
Due to the restricted supports of the subvectors $\tilde{\bm{y}}_{b,t}$ and $\tilde{\bm{y}}_{c,t}$, it is inappropriate to model $\tilde{\bm{y}}_t$ using a VAR with errors drawn from a continuous distribution. A standard approach to handling discrete and censored variables is via data augmentation: introduce latent variables $\bm{y}_{b,t}\in\mathbb{R}^{n_b}$ and $\bm{y}_{c,t}\in\mathbb{R}^{n_c}$ that are linked to, respectively, the binary observations $\tilde{\bm{y}}_{b,t}$ and censored variables $\tilde{\bm{y}}_{c,t}$ through
where $\mathbb{I}(E)$ is an indicator function that takes the value 1 if the event $E$ is true and 0 otherwise. Next, we define $\bm{y}_t = (\bm{y}_{b,t}',\bm{y}_{c,t}',\bm{y}_{u,t}')'$ which collects the latent variables associated with the binary and censored variables, and the observed continuous variables. Then, this vector can be modeled as a VAR($P$) process:
where $\bm{A}_1, \ldots, \bm{A}_P $ are the $n\times n$ coefficient matrices and $\bm{\epsilon}_t$ is a zero mean Gaussian error with an $n \times n$ covariance matrix $\bm{\Sigma}_t$. Since only the signs of the binary variables $\tilde{\bm{y}}_{b,t}$ are identified, the diagonal elements of $\bm{\Sigma}_t$ corresponding to the variances of $\tilde{\bm{y}}_{b,t}$ need to be fixed (e.g., set equal to 1) for identification. This identification restriction makes it more difficult to specify a suitable time-varying model for $\bm{\Sigma}_t$, especially when $n$ is large. To strike the right balance between flexibility and computational feasibility, we assume that the time-varying covariance matrix takes the form $\bm{\Sigma}_t = o_t\bm{\Sigma}$, where $\bm{\Sigma}$ is an $n\times n$ positive-definite matrix and $o_t$ is a time-varying scalar. By specifying different processes for $o_t$, this setup includes a wide variety of commonly used heteroskedastic or non-Gaussian error structures designed for large VARs, such as the common stochastic volatility model of CCM16, heavy-tailed error distributions considered in chan20, an outlier component proposed in SW16, and a volatility break model for fitting post COVID-19 pandemic data developed in LP22; see also CQ25 for other examples. Naturally, it also nests the standard homoskedastic VAR with $o_t =1, t=1,\ldots, T$.
Next, we partition $\bm{\Sigma}$ into blocks associated with the variable types (binary, censored and unrestricted), e.g., $\bm{\Sigma}_{bb}, \bm{\Sigma}_{cc}$ and $ \bm{\Sigma}_{uu}$ are $n_b\times n_b$, $n_c\times n_c$, and $n_u\times n_u$ matrices on the (block) diagonal, respectively. The other submatrices that collect the covariances are conformally defined. We noted previously that the scale of $\tilde{\bm{y}}_{b,t}$ is not identifiable. Hence, for identification purposes, the diagonal elements of $\bm{\Sigma}_{bb}$ are set to 1.
To summarize, the proposed framework incorporates multivariate probit and tobit models to construct a VAR. It generalizes a variety of VARs previously considered in the literature. When $n_b = n_c = 0$, it reduces to the conventional VAR with continuous variables; when $n_c = 0$, it becomes a multivariate generalization of the so-called QualVAR of dueker2005dynamic that contains a single binary variable; when $n_b = 0$, it then becomes what has been termed a simple shadow rate VAR by carriero2023shadow in an application to infer interest rates at the ELB.\footnote{Further extensions of the proposed framework can be considered. For example, one could include ordinal and categorical variables by introducing additional latent variables. zhang2015bayesian consider some of these extensions in a static setting; PPZ25 propose a related framework for count data.}
If $\bm{y}_{b,t}$ and $\bm{y}_{c,t}$ were observed, it would be relatively straightforward to estimate the VAR in (ref) even when $n$ is large---one can simply employ any recently developed efficient posterior sampler designed for large VARs with complete data CCM19,CCCM22,chan23JE. Of course $\bm{y}_{b,t}$ and $\bm{y}_{c,t}$ are not observed, but this missing data problem can be tackled by augmenting the posterior sampler for complete data with an extra step of simulating the missing values from their joint conditional posterior distribution. Due to the modular nature of MCMC methods, all one needs is an extra block to draw the missing values.
However, while this additional step is conceptually simple, in practice it is challenging in our high-dimensional, dynamic setting for two reasons. First, unlike in static models where the missing data are serially independent given other model parameters and the observed data---one can therefore sample $(\bm{y}_{b,t}', \bm{y}_{c,t}')'$ for $t=1,\ldots, T$ one at a time sequentially without loss of efficiency---they are serially dependent in our setting because of the VAR structure. Consequently, while drawing $(\bm{y}_{b,t}', \bm{y}_{c,t}')'$ one at a time would induce high autocorrelations and poor mixing in the Markov chain, sampling all missing data jointly using conventional Kalman filter-based methods is computationally challenging, if not practically infeasible. Second, the conditional posterior distribution of the missing data is subject to a large number of inequality restrictions. For example, if $\tilde{y}_{b,it} = 0$, then the corresponding missing value should be negative, that is, $y_{b,it}< 0 $. This thus compounds the already difficult sampling problem.
To address these challenges, we adapt the algorithm of chan2023high for drawing missing data designed for conditionally Gaussian state space models to our setting. Specifically, chan2023high show that the conditional posterior density of the missing data given the model parameters and observed data is Gaussian. Moreover, its inverse covariance or precision matrix is banded---i.e., it is sparse and the non-zero elements are arranged along a main diagonal band---and this special structure can be exploited to vastly speed up the computations. We adapt the algorithm to our setting and extend it to incorporate a large number of inequality restrictions.
To that end, it is useful to divide $\bm{y}_t$ into two subvectors, $\bm{y}_{\text{\texttt{l}}}$ and $\bm{y}_{\text{\texttt{o}}}$, consisting of variables that are latent and observed, respectively (with subscripts l and o). Furthermore, latent variables arise from two cases: the first is from data augmentation when we are dealing with binary or censored variables, and the second occurs whenever observations are missing (e.g., at the beginning or end of the sample because of an unbalanced panel or release delay) irrespective of the variable type. For notational convenience, we begin by partitioning $\bm{y}_t$:
With these definitions, we thus have $n_{b,t} = n_{b,t}^{\text{\texttt{m}}} + n_{b,t}^{\text{\texttt{r}}}$, $n_{c,t} = n_{c,t}^{\text{\texttt{m}}} + n_{c,t}^{\text{\texttt{r}}} + n_{c,t}^{\text{\texttt{o}}}$, $n_{u,t} = n_{u,t}^{\text{\texttt{m}}} + n_{u,t}^{\text{\texttt{o}}}$. In addition, the total number of latent variables at time $t$ is $n_{t}^{\text{\texttt{l}}} = n_{b,t}^{\text{\texttt{m}}} + n_{b,t}^{\text{\texttt{r}}} + n_{c,t}^{\text{\texttt{m}}} + n_{c,t}^{\text{\texttt{r}}} + n_{u,t}^{\text{\texttt{m}}}$ and the number of observed variables is $n_{t}^{\text{\texttt{o}}} = n - n_{t}^{\text{\texttt{l}}} = n_{c,t}^{\text{\texttt{o}}} + n_{u,t}^{\text{\texttt{o}}}$. Finally, we can divide $\bm{y}_t$ into the subvectors $\bm{y}_{\text{\texttt{l}}}$ and $\bm{y}_{\text{\texttt{o}}}$:
with $\bm{s}_{\text{\texttt{l}},t} \bm{y}_{\text{\texttt{l}},t} = \bm{s}_{\text{\texttt{m}},t} \bm{y}_{\text{\texttt{m}},t} + \bm{s}_{\text{\texttt{r}},t} \bm{y}_{\text{\texttt{r}},t},$ where the selection matrix $\bm{s}_{\text{\texttt{l}},t}'$ of dimension $n_{t}^{\text{\texttt{l}}} \times n$ is defined by stacking the selection matrices $\bm{s}_{\text{\texttt{m}},t}'$ and $\bm{s}_{\text{\texttt{r}},t}'$.\footnote{To derive Eq. (ref), first note that $\bm{s}_t = [\bm{s}_{\text{\texttt{m}},t}, \bm{s}_{\text{\texttt{r}},t}, \bm{s}_{\text{\texttt{o}},t}] $ is an $n\times n$ permutation matrix with $\bm{s}_t\bm{s}_t' = \bm{I}_n$, where $\bm{I}_n$ is the $n$-dimensional identity matrix. By construction, we have $ \bm{s}_t'\bm{y}_t = [\bm{y}_{\text{\texttt{m}},t}', \bm{y}_{\text{\texttt{r}},t}', \bm{y}_{\text{\texttt{o}},t}']'.$ After pre-multiplying both sides by $\bm{s}_t$, we thus obtain $\bm{y}_t = [ \bm{s}_{\text{\texttt{m}},t}, \bm{s}_{\text{\texttt{r}},t}, \bm{s}_{\text{\texttt{o}},t}] [ \bm{y}_{\text{\texttt{m}},t}', \bm{y}_{\text{\texttt{r}},t}', \bm{y}_{\text{\texttt{o}},t}']' = \bm{s}_{\text{\texttt{m}},t} \bm{y}_{\text{\texttt{m}},t} + \bm{s}_{\text{\texttt{r}},t} \bm{y}_{\text{\texttt{r}},t} + \bm{s}_{\text{\texttt{o}},t} \bm{y}_{\text{\texttt{o}},t} = \underbrace{[ \bm{s}_{\text{\texttt{m}},t}, \bm{s}_{\text{\texttt{r}},t}]}_{\bm{s}_{\text{\texttt{l}},t}} [\bm{y}_{\text{\texttt{m}},t}', \bm{y}_{\text{\texttt{r}},t}']' + \bm{s}_{\text{\texttt{o}},t} \bm{y}_{\text{\texttt{o}},t}$. } For later reference, we also define the average numbers of latent and observed variables, $n^{\text{\texttt{l}}} = \sum_{t=1}^T n_{t}^{\text{\texttt{l}}}/T$ and $n^{\text{\texttt{o}}} = \sum_{t=1}^T n_{t}^{\text{\texttt{o}}}/T$.
The next step is to derive the joint density of $\bm{y} = (\bm{y}_{-P+1}',\ldots,\bm{y}_{0}',\bm{y}_{1}',\ldots,\bm{y}_{T}')'$, which is of dimension $(T+P)n\times 1$. Let $\bm{y}_{\text{\texttt{l}}} $ and $\bm{y}_{\text{\texttt{o}}}$ denote, respectively, the vectors of all latent variables and observed data. Furthermore, let $\text{bdiag}(\bullet)$ represent the operator that constructs a block diagonal matrix from its inputs. Next, we define the selection matrices $\bm{S}_{\text{\texttt{l}}} = \text{bdiag}(\bm{s}_{\text{\texttt{l}},(-P+1)},\hdots,\bm{s}_{\text{\texttt{l}},T})$ and $\bm{S}_{\text{\texttt{o}}} = \text{bdiag}(\bm{s}_{\text{\texttt{o}},(-P+1)},\hdots,\bm{s}_{\text{\texttt{o}},T})$, so that $\bm{S}_{\text{\texttt{l}}}'\bm{y} = \bm{y}_{\text{\texttt{l}}}$ and $\bm{S}_{\text{\texttt{o}}}'\bm{y} = \bm{y}_{\text{\texttt{o}}}$. It follows that
Following chan2023high, we encode the VAR dynamics in the $Tn \times (T+P)n$ matrix
so that Eq. (ref), for $t = 1,\ldots,T,$ can then be written compactly as:
where $\bm{\Omega} = (\text{diag}(o_1, \cdots, o_T) \otimes \bm{\Sigma})$ is the $Tn \times Tn$ error covariance matrix, $\otimes$ is the Kronecker product, $\bm{m} = (\bm{\iota}_T \otimes \bm{a})$ is the $Tn$-vector of intercepts, and $\bm{\iota}_T$ is a $T$-vector of ones.
Substituting the decomposition of $\bm{y}$ given in Eq. (ref) into Eq. (ref), we have $\bm{M}\bm{S}_{\text{\texttt{l}}} \bm{y}_{\text{\texttt{l}}} + \bm{M}\bm{S}_{\text{\texttt{o}}} \bm{y}_{\text{\texttt{o}}} = \bm{m} + \bm{\epsilon}$. Using the properties of the Gaussian distribution, one can derive the distribution of the latent variables $\bm{y}_{\text{\texttt{l}}}$ (missing values and restricted processes) conditional on the (unrestricted) observables $\bm{y}_{\text{\texttt{o}}}$. Specifically, defining $\bm{H}_{\text{\texttt{l}}} = \bm{M}\bm{S}_{\text{\texttt{l}}}$ and $\bm{H}_{\text{\texttt{o}}} = \bm{M}\bm{S}_{\text{\texttt{o}}}$, we obtain again a Gaussian $p(\bm{y}_{\text{\texttt{l}}}\:\vert\:\bm{y}_{\text{\texttt{o}}},\bullet) \propto \exp\left(-\frac{1}{2} (\bm{H}_{\text{\texttt{l}}}\bm{y}_{\text{\texttt{l}}} + \bm{H}_{\text{\texttt{o}}}\bm{y}_{\text{\texttt{o}}} - \bm{m})'\bm{\Omega}^{-1}(\bm{H}_{\text{\texttt{l}}}\bm{y}_{\text{\texttt{l}}} + \bm{H}_{\text{\texttt{o}}}\bm{y}_{\text{\texttt{o}}} - \bm{m})\right)$, where $\bullet$ stands for conditioning on all parameters. Following the derivations in chan2023high, one can show that $\bm{y}_{\text{\texttt{l}}}\:\vert\:(\bm{y}_{\text{\texttt{o}}},\bullet)\sim\mathcal{N}(\bm{\mu}_{\text{\texttt{l}}},\bm{V}_{\text{\texttt{l}}})$, with moments:
Note that the matrices $\bm{H}_{\text{\texttt{l}}}, \bm{H}_{\text{\texttt{o}}}$ and $\bm{\Omega}$ are all banded, and so is the inverse covariance or precision matrix $\bm{V}_{\text{\texttt{l}}}^{-1}$. This banded structure can be exploited to vastly speed up computations---efficient implementations of sparse and band matrix operations are available in standard statistical software.
Next, we further condition on the observed binary and censored variables $\tilde{\bm{y}}_b$ and $\tilde{\bm{y}}_c$ to obtain $p(\bm{y}_{\text{\texttt{l}}}\:\vert\:\tilde{\bm{y}}_b,\tilde{\bm{y}}_c,\bm{y}_{\text{\texttt{o}}},\bullet)$, which is a truncated multivariate normal distribution:
where $\bm{R}$ is a diagonal matrix incorporating the inequality restrictions in Eqs. ((ref)) and ((ref)). Specifically, if the $j$th element of $\bm{y}_{\text{\texttt{l}}}$ is a binary variable, $\bm{R}_{jj} = 2\tilde{y}_{b,it} - 1$; if it is a censored variable with $\tilde{y}_{c,it} = 0$, $\bm{R}_{jj} = -1$; and if it is a missing value, there are no associated restrictions, and we set $\bm{R}_{jj} = 0$. Obtaining samples from high-dimensional truncated normal distributions is, in general, challenging. While there is much recent progress botev2017normal, direct sampling approaches remain computationally prohibitive for our purpose that involves very high dimensions.\footnote{For missing values, the corresponding components of the multivariate normal distributions are unrestricted. In principle, one could sample from the distribution in Eq. ((ref)) in two steps: one samples from the lower dimensional truncated normal distribution $p(\bm{y}_{\text{\texttt{r}}} \:\vert\:\tilde{\bm{y}}_b,\tilde{\bm{y}}_c,\bm{y}_{\text{\texttt{o}}},\bullet)$ and then draws the unrestricted missing values from the normal distribution $p(\bm{y}_{\text{\texttt{m}}} \:\vert\:\bm{y}_{\text{\texttt{r}}},\bm{y}_{\text{\texttt{o}}},\bullet)$ conditional on $\bm{y}_{\text{\texttt{r}}}$. In typical applications, however, this two-step scheme does not result in significant computational gains, because the number of missing values is usually small compared to the number of restricted observations.} Instead, we implement variants of the HMC pakman2014exact to sample from the high-dimensional truncated normal distribution, as discussed in the next section.
In this section, we outline the priors on the model parameters and provide details of key sampling steps in the posterior simulation.
\noindentSampling the latent states. As discussed in Section (ref), a key step in our estimation approach is sampling the latent variables. Given the values of these latent variables, the estimation problem reduces to fitting a large VAR with complete data, and any efficient sampling method can then be used. We therefore focus on the step of simulating from the high-dimensional truncated normal distribution given in Eq. ((ref)). Here we implement an HMC step, and specifically versions referred to as harmonic HMC or Zigzag HMC nishimura2020discontinuous. They are efficient in terms of the low autocorrelations of the MCMC draws and the time required for simulations. The latter HMC is particularly attractive due to the band structure of the precision matrix of the normal distribution.
We provide a brief sketch of the relevant HMC algorithms below. For simplicity, we omit the conditioning arguments. To sample from the target distribution $p(\bm{y}_\text{\texttt{l}})$, we introduce the joint distribution $p(\bm{y}_\text{\texttt{l}},\bm{\xi}) = p(\bm{y}_\text{\texttt{l}})p(\bm{\xi})$, where $\bm{\xi}$ is called the momentum variable. Define $U(\bm{y}_\text{\texttt{l}}) = -\log(p(\bm{y}_\text{\texttt{l}}))$ and $K(\bm{\xi}) = -\log(p(\bm{\xi}))$, which are referred to as potential and kinetic energy, respectively; the total energy $H(\bm{y}_\text{\texttt{l}},\bm{\xi}) = U(\bm{y}_\text{\texttt{l}}) + K(\bm{\xi})$ is called a Hamiltonian.
We may use these quantities to simulate the solution of the Hamiltonian dynamics, which are governed by the differential equations $\text{d}\bm{y}_\text{\texttt{l}} / \text{d}s = \nabla K(\bm{\xi})$ and $\text{d}\bm{\xi} / \text{d}s = -\nabla U(\bm{y}_\text{\texttt{l}})$ with respect to the time $s$ for a set time duration $S$. The end state yields a proposal that can be used in a Metropolis-Hastings accept/reject step. Harmonic HMC chooses a Gaussian momentum distribution, which yields a closed-form solution when the target density is multivariate truncated normal, see pakman2014exact. Zigzag HMC can be implemented by choosing a Laplace momentum distribution, i.e., $p(\bm{\xi}) \propto \prod_{i=1}^{Tn^{\text{\texttt{l}}}} \exp\left(-|\xi_i|\right)$, where we obtain $\text{d}\bm{y}_\text{\texttt{l}} / \text{d}s = \text{sign}(\bm{\xi})$, which depends only on the sign of $\bm{\xi}$ but not its magnitude; this results in several favorable properties, see nishimura2020discontinuous,nishimura2021hamiltonian,zhang2022hdtg for more discussion.
Our default implementation uses Zigzag HMC because we may compute the mean in Eq. (ref) without explicitly obtaining the posterior covariance matrix via the inversion of $\bm{V}_\text{\texttt{l}}^{-1}$. The precision matrix is also sufficient to set up this algorithm. By contrast, harmonic HMC requires a Cholesky decomposition of either the covariance or precision matrix and Zigzag HMC is thus preferable in high-dimensional settings, at least from a computational viewpoint. Appendix (ref) provides some simulation results that show that the proposed algorithm works well.
\noindentRestrictions on the covariance matrix. As noted in Section (ref), the scale of $\tilde{\bm{y}}_{b,t}$ is not identifiable, and the diagonal elements of $\bm{\Sigma}_{bb}$---and hence the first $n_b$ diagonal elements of $\bm{\Sigma}$---are set to 1. These equality restrictions make the sampling of $\bm{\Sigma}$ more challenging. To address this difficulty, we follow the parameter expansion approach liu1999parameter, meng1999seeking by introducing $n_b$ auxiliary (non-identifiable) parameters $\bm{d} = (d_1,\ldots, d_{n_b})'$ and define $\bm{\Sigma}^{\ast} = \bm{D}^{1/2} \bm{\Sigma} \bm{D}^{1/2}$, where $\bm{D}=\text{diag}(\bm{d}', \bm{\iota}_{n_c}',\bm{\iota}_{n_u}')$. We denote this mapping from $(\bm{d},\bm{\Sigma})$ to $\bm{\Sigma}^{\ast}$ as $f$.\footnote{It is straightforward to see that $f$ is one-to-one and its inverse can be computed easily. Specifically, given $\bm{\Sigma}^{\ast}$, one can extract its first $n_b$ diagonal elements and set $\bm{d} = (\sigma^{\ast}_{11},\ldots, \sigma^{\ast}_{n_b n_b})'$. Then, one obtains $ \bm{\Sigma}$ via $ \bm{\Sigma} = \bm{D}^{-1/2}\bm{\Sigma}^{\ast}\bm{D}^{-1/2}$.} Given a proper prior on $\bm{d}$, the expanded joint posterior distribution is proper even though $\bm{d}$ is non-identifiable. The main advantage of this parameter expansion approach is that it is easier to sample from the expanded joint posterior distribution. An additional benefit is that the Markov chain produced by the expanded, non-identifiable model typically has better convergence properties compared to that from the original, identifiable model.
We assume an inverse Wishart prior on $\bm{\Sigma}^{\ast}$ with degree of freedom $n+2$ and scale matrix $\bm{I}_n$, i.e., $\bm{\Sigma}^{\ast} \sim \mathcal{W}^{-1}(n + 2, \bm{I}_n)$, and we denote its density function as $p_{\mathcal{W}^{-1}} (\bm{\Sigma}^{\ast})$. This prior induces a joint prior on $(\bm{\Sigma}, \bm{d})$. To derive the density of the induced prior, note that the Jacobian of transformation $f$ is $\text{det}(\bm{D})^{(n-1)/2}$. Hence, the induced density is given by $p(\bm{\Sigma},\bm{d}) = \text{det}(\bm{D})^{(n-1)/2} p_{\mathcal{W}^{-1}}(f(\bm{\Sigma},\bm{d}))$. Putting all this together, the joint conditional posterior has the form:
This distribution is not of a known form, and we follow the parameter-expanded Metropolis-Hastings (PX-MH) algorithm proposed in zhang2020parameter to sample from it.\footnote{dueker2005dynamic assumes an inverse Wishart prior directly on $\bm{\Sigma}$ and samples from the unrestricted posterior. The resulting draw is normalized ex post such that it aligns with the required restrictions.}
\noindentPriors on model parameters. We assume a global-local shrinkage prior on the VAR coefficients. Specifically, we collect $\bm{A} = (\bm{A}_1, \hdots, \bm{A}_P, \bm{a})'$ and vectorize $\bm{\beta} = \text{vec}(\bm{A})$. Each of the elements of $\bm{\beta}$ is assumed to follow an independent Gaussian prior: $\beta_j \sim \mathcal{N}(0, \tau^2\lambda_j^2)$, for $j = 1, \hdots, n^2 P + n.$ The prior variance is determined by two components, the global variance $\tau^2$ that is common to all VAR coefficients and the local component $\lambda_j^2$. In particular, we use the horseshoe prior carvalho2010horseshoe due to its favorable shrinkage properties and because there are no additional tuning parameters required. We assume $\tau\sim\mathcal{C}^{+}(0,1)$ that pushes all parameters in $\bm{\beta}$ towards the prior mean of $0$; $\lambda_j\sim\mathcal{C}^{+}(0,1)$ are local scalings that may offset heavy global shrinkage if the $j$th parameter is important. Note that we do not shrink the intercept terms hierarchically but instead impose a large prior variance.
As discussed in Section (ref), the error covariance matrix is assumed to take the form $\bm{\Sigma}_t = o_t\bm{\Sigma}$, which can accommodate many commonly used heteroskedastic or non-Gaussian covariance structures designed for large VARs. Here we adopt the heteroskedastic specification of SW16, which is used in a VAR setting for modeling COVID-19 outliers in carriero2021addressing. In particular, we assume that:
where $\mathcal{U}(2,10)$ is a discrete uniform distribution with support between $2$ and $10$ and $\mathfrak{p}\sim\mathcal{B}(a_\mathfrak{p},b_\mathfrak{p})$ is the probability associated with observing an outlier. Conditional on the latent states, updating most of the other parameters of our model is standard. In particular, the VAR coefficients are drawn equation-by-equation using a version of the efficient sampling approach of CCCM22.
Finally, we standardize the dataset (apart from the binary variables) before estimation to limit the influence of the scalings of the included variables. All relevant estimates and model outputs are subsequently scaled back to the original units of the variables ex post.
The VAR$(P)$ process for the latent and observed variables, as defined in Eq. ((ref)), can be used to compute multi-step forecasts using standard approaches. The joint distribution of the forecasts up to horizon $h$ made at time $\tau$, $\bm{y}_{\tau+1:\tau+h} = (\bm{y}_{\tau+1}',\hdots,\bm{y}_{\tau+h}')'$, conditional on all model parameters $\bm{\Xi}$ and the respective information set $\mathcal{I}_{1:\tau}$, is the Gaussian distribution:
with known moments that can be derived from considerations analogous to Eq. ((ref)) in Section (ref); see also chan2024conditional for explicit expressions. We may obtain a sample from this distribution in each iteration of our sampling algorithm and collect these draws from the predictive distribution $p(\bm{y}_{\tau+1:\tau+h}\:\vert\:\mathcal{I}_{1:\tau}) = \int p(\bm{y}_{\tau+1:\tau+h} \:\vert\:\mathcal{I}_{1:\tau}, \bm{\Xi}) p(\bm{\Xi} \:\vert\:\mathcal{I}_{1:\tau}) d \bm{\Xi}$. Forecasts for the continuous variables are directly available from the predictive distribution, while those for binary and censored variables can be obtained by applying the transformations in Eqs. ((ref)) and ((ref)). For example, one can compute the probability forecast associated with the binary variables via $\Pr(\tilde{y}_{b,i\tau+h} = 1) = \Phi(\mu_{y_b,i\tau+h} v_{y_b,i\tau+h}^{-1/2})$, where $ \Phi(\cdot)$ is the cumulative distribution of the standard normal distribution. This quantity is calculated by singling out the moments associated with the binary variable at time $\tau+h$ from $\bm{\mu}_{y,\tau h}$ and $\bm{V}_{y,\tau h}$, which we denote as $\mu_{y_b,i\tau+h}$ and $v_{y_b,i\tau+h}$; see also mccracken2022binary. We may compute this probability for each sweep of our sampler which provides us with draws from the corresponding posterior.
\noindentConditional forecasts. Next, we consider conditional forecasts implied by hard and/or soft constraints on the path of one or more of the endogenous variables. They can be obtained by conditioning the predictive distribution in Eq. ((ref)) on an additional information set $\mathcal{R}$, i.e., $p(\bm{y}_{\tau+1:\tau+h}\:\vert\:\mathcal{R}, \mathcal{I}_{1:\tau})$. We focus on conditioning-on-observables forecasting in our empirical work, but any standard VAR-based tool antolin2021structural can be directly applied using the latent representation given in Eq. ((ref)).
We consider two types of restrictions: hard and soft constraints, expressed as equality and inequality restrictions, respectively. Let $\bm{R}_{c}$ be an $r \times nh$ selection matrix that selects the restricted variables and horizons, and let $\bm{r}$ denote an $r$-vector of constants. Then, the set the hard restrictions can be written as $\mathcal{R} = \{\bm{R}_{c}\bm{y}_{\tau+1:\tau+h} = \bm{r}\}$. Similarly, if we let the $r$-vectors $\underline{\bm{r}}$ and $\overline{\bm{r}}$ denote, respectively, the lower and upper bounds, then the set of soft constraints can be written as $ \mathcal{R} = \{\underline{\bm{r}} \leq \bm{R}_{c} \bm{y}_{\tau+1:\tau+h} \leq \overline{\bm{r}}\}.$ Naturally, one can combine both hard and soft constraints. The conditional predictive distribution $p(\bm{y}_{\tau+1:\tau+h}\:\vert\:\mathcal{R}, \mathcal{I}_{1:\tau})$ is a multivariate truncated Gaussian distribution, but one can sample from this distribution efficiently. We refer the readers to chan2024conditional for computational details.
\noindentGeneralized impulse response functions. Impulse response functions (IRFs) for the latent representation are easy to compute using standard methods, but the nonlinearities arising from Eqs. ((ref)) and ((ref)) require the use of generalized IRFs koop1996impulse for the binary and truncated variables. The former, i.e., the default IRFs of the latent and observed variables in response to a shock, are time-invariant and symmetric for different shock signs and sizes as in any linear VAR model. By contrast, the GIRFs are needed because the configuration of the vector of endogenous variables at the time of the occurrence of the shock matters. For related discussions, see pfarrhofer2025scenario.
Suppose we have identified the contemporaneous effect of a structural shock on the vector of endogenous variables using tools from the literature on structural VARs (e.g., via timing/sign restrictions, or with internal/external instruments). We store this impact in an $n$-vector $\bm{\delta}$. We rely on Eq. ((ref)) to simulate the no-shock scenario first. However, different from before, we compute the moments of the predictive distribution $p(\bm{y}_{\tau+1:\tau+h}|\bm{\Xi},\mathcal{I}_{1:T})$ at each point in time $\tau$, conditional on the full information set $t = 1,\hdots,\tau,\hdots,T$. The counterfactual shock scenario can be obtained by adjusting the initial configuration that is projected forward when forecasting by the impact of the identified shock of interest at time $\tau$, such that $\bm{y}_{\tau|\delta} = \bm{y}_{\tau} + \bm{\delta}$. This also allows for simulating shocks of different signs and sizes, which potentially result in asymmetric GIRFs. The probability predictions for the shock and non-shock scenarios for the binary indicators can then be obtained via the cdf discussed above. The difference between the two yields the corresponding GIRF. Similarly, we can apply Eq. ((ref)) to obtain the truncated predictions, and compute the GIRF as the difference between the conditional and unconditional forecasts.
Computing the GIRF for each point in time relates to common summary interpretation schemes in the context of binary dependent variable models. Specifically, what we compute in this case, by varying the initial configurations $\bm{y}_{\tau|\delta}$ and obtaining estimates of the dynamic causal effects at each $\tau$, is usually processed further by taking a time average ex post. This yields an average partial effect. By contrast, one may also compute a partial effect at the average, which would imply synthetically constructing a single initial configuration $\overline{\bm{y}}_{\delta} = \sum_{\tau=1}^T \bm{y}_{\tau} + \bm{\delta}$. While the latter is computationally attractive, this may lead to unusual initial configurations $\overline{\bm{y}}_{\delta}$ that were never observed in reality. So although it is computationally more expensive, we opt for the former approach. This also implies that we can investigate potential time-variation (due to variation over initial conditions) in the GIRFs.
\FloatBarrier
We conduct an out-of-sample forecasting exercise to evaluate predictions of recession probabilities and forecasts of other macroeconomic and financial variables. Our full sample runs on a monthly frequency from 1959M2 until 2023M12 and uses a variety of variables taken from the well-known FRED-MD dataset described in mccracken2016fred. The initial training sample ends in 1979M12, and we simulate monthly updates of our predictions with an expanding window of observations. We do not use real-time vintages due to the relatively short sample this would entail, but our forecasting design respects the release calendar (i.e., it simulates the information set at each point in time that would have been available to forecasters with each release, and there are ragged edges).
To estimate recession probabilities, our binary target variable is the indicator of expansions and recessions (contractions) released by the \href{https://www.nber.org/research/data/us-business-cycle-expansions-and-contractions}{NBER Business Cycle Dating Committee}. This indicator does not have a structured publication schedule, because only business cycle turning points are announced. To simulate realistic information sets, we thus follow the precedent set by giusto2017identifying and mccracken2022binary. Specifically, after turning points, we treat the classification of the subsequent period as known; after troughs, we assume that the announcement delay for the following peak is shorter than 12 months. Following the dating of a peak, we treat only the first month of this recession as known. This irregular dating schedule implies that the window of horizons $h$ that needs to be backcasted ($h < 0$) or nowcasted $h = 0$ changes over time. That is, the number of out-of-sample observations we use to evaluate our models varies across horizons. For forecasts ($h > 0$), we consider up to two-year-ahead (24 months) predictions.
We consider four differently sized information sets: core ({ C}), small ({ S}), medium ({ M}), and large ({ L}). Table (ref) provides a list of the variables included in each of these specifications alongside transformation codes. We estimate each of the VARs with $P = 6$ lags. The core information set uses the key variables assessed by the NBER Business Cycle Dating Committee---non-farm payroll employment (PAYEMS), industrial production (INDPRO), real manufacturing and trade sales (CMRMTSPLx), real personal income excluding transfers (W875RX1)---and adds the federal funds rate (FEDFUNDS) and a term spread (TS10G1) as important financial predictors. This information set mirrors the preferred specification of mccracken2022binary.
The small information set adds interest rates at several maturities. Note that some of these interest rates have hit their ELB in our sample. Here, we follow the definition of carriero2023shadow, and consider all observations below a threshold of $0.25$ percent as constrained; this results in FEDFUNDS, TB3MS, TB6MS, and GS1 having at least one observation that is treated as censored (for the full sample). The medium information set then adds additional financial variables and oil prices, while the large information set adds mostly labor market variables and prices. For the excess bond premium (EBP), no observations are available before 1973. Our approach treats this as a partially missing data problem and imputes these observations for this part of the sample.
Our most general model is called the { ProToVAR}, because it features both the binary target variable (addressed by the { Pro}bit part) and censored variables (the { To}bit part). Another specification includes the binary target variable but treats all other variables as continuous. It thus only features the { Pro}bit part, and we label it the { ProVAR}. For a significant part of our holdout sample, these two specifications coincide because none of the interest rates hit the ELB. For each of these models, we consider two versions, one with homoskedastic errors (labeled { hom}) and the other with the outlier component that captures heteroskedasticity (labeled { het}). All specifications are equipped with the same prior setup.\footnote{For all of our empirical work, we cycle through our algorithm 12,000 times, discard the initial 3,000 draws as burn-in, and use every $3$rd of the remaining draws for inference; i.e., we retain 3,000 draws.}
To assess the predictive accuracy of the forecasts for the binary recession variable, we use a measure based on the receiver operating characteristic (ROC) curve. The ROC curve represents the false positive rate (FPR) and true positive rate (TPR) as a function of probability thresholds (i.e., on the unit square). A common evaluation metric for binary predictions is the area under the ROC curve (AUC). The larger the AUC, the better is the predictive performance of the respective model over the holdout. We compute the AUC for each model size as a function of the predicted horizon $h = -10, -9, \hdots, -1, 0, 1,\hdots, 24,$ for the full evaluation period (1980M1 to 2023M12) and two sub-periods (the “pre-COVID” sample ends in 2019M12, while the “Post 2000” sample begins in 2000M1).
The AUC values across the model specifications, information sets, subsamples, and forecast horizons are reported in Figure (ref). The rows correspond to evaluation periods (full, pre-COVID and post 2000) and the columns refer to different datasets (core, small, medium and large). For each sample period and forecast horizon, we indicate the best performing specification (combination of a model and dataset) using a colored square at the bottom of each panel. Finally, note that the homoskedastic { ProVAR-hom} specification estimated using the core dataset coincides with the preferred model of mccracken2022binary.
Our results show that while the core dataset appears to be sufficient for retrospectively dating recessions (i.e., backcasting with $h<0$) for the full sample, additional variables improve recession forecasts for all $h>0$ horizons. For example, models estimated using the large dataset with all 25 variables consistently outperform alternatives over short horizons (e.g., $1\leq h \leq 7$); adding interest rates and other financial variables also tend to improve recession forecasts over medium and long horizons. The results are even more stark if we consider the post 2000 subsample, in which models estimated using the medium and large datasets dominate. Across all model specifications, results are more mixed and no models dominate for all horizons. For the full sample, { ProToVAR-het} performs best for dating recessions; { ProVAR-het} forecasts well for horizons less than a year; and { ProVAR-hom} tends to be the best model for horizons more than a year. But for post 2000 subsample, { ProToVAR-het} and { ProVAR-het} tend to outperform their homoskedastic counterparts.
Figure (ref) plots the out-of-sample recession probabilities across multiple horizons from { ProToVAR-het}. Naturally, forecasting recession one-year-ahead ($h=12$) is inherently difficult, but the model is able to pick up some early warning signs, especially when the large dataset is used. The accuracy of the model improves markedly for nowcasting recessions ($h = 0$), and it is able to closely track the NBER recession dates.
Overall, two broad conclusions can be drawn from this forecasting exercise. First, our findings highlight the benefits of exploiting the information content in larger datasets that include a wide variety of financial and labor market variables in improving recession forecasts. Second, while no forecasting model strictly dominates, allowing for heteroskedasticity tends to improve forecast performance, especially for the post 2000 subsample.
Next, we report forecasting results for some key interest rates. We focus on the post 2000 subsample since the ELB became binding in this later part of the holdout sample. For comparison, we include two models that omit the binary recession indicator. The first is a standard VAR which treats all variables as continuous, which we denote as { VAR}. The second model treats the interest rates as censored variables, and is referred to as { ToVAR}. We use the continuous ranked probability score gneiting2007strictly as a density forecast metric.
Table (ref) presents the CRPSs for different models ({ VAR}, { ProVAR}, { ToVAR} and { ProToVAR}), with and without heteroskedasticity ({ het} and { hom}), estimated using different datesets (small, medium and large). We use the homoskedastic VAR estimated with the large dataset ({ VAR-hom-L}) as the benchmark model. The values in the row associated with { VAR-hom-L} are raw CRPSs, and all other values are relative to this benchmark as ratios. Values smaller than $1$ (blue shades) signal better relative performance; values larger than $1$ (red shades) indicate worse relative performance. The best performing model for each horizon and variable is marked in bold.
Our results show that models that treat interest rates as censored variables forecast well for short-term interest rates up to a year relative to { VAR} models that do not. They perform especially well for horizons longer than 1-month-ahead, as iterative forecasts from these models are guaranteed to be above the ELB, whereas { VAR} models sometimes produce longer-horizon forecasts that are below the ELB. Moreover, it is interesting to note that { ProToVAR-het} models consistently outperform their corresponding { ToVAR-het} models for forecasting short-term interest rates. This highlights the value of incorporating the binary recession indicator, which contains additional information beyond that is encoded in the many key macroeconomic variables. Furthermore, consistent with the recession forecasting results, models estimated using larger datasets tend to forecast better than those fitted using smaller datasets, and heteroskedastic models outperform their homoskedastic counterparts. Overall, { ProToVAR-het-L} tend to be the best model for forecasting interest rates up to a year. These results again underscore the importance of designing time-series models that respect the range of the target variables. In addition, they also highlight the empirical relevance of using a large set of variables and accommodating heteroskedasticity.
Results for the full sample as well as those for other (unrestricted) macroeconomic variables can be found in Appendix (ref). We note that the additional complexity due to including binary and censored variables in our VARs does not materially affect predictive losses of the forecasts of the unrestricted variables. Indeed, { ProToVAR-het} is always competitive, and in several cases even yields improvements in predictive accuracy for these variables as well.
\FloatBarrier
In this section we consider several applications to demonstrate the usefulness of the proposed framework. We first investigate the dynamic evolution of the model-based latent recession indicator and the shadow interest rates. We then illustrate how one can employ conditional forecasting to track the impact of two scenarios on key macroeconomic variables. The first relates to the 2020 Federal Reserve Board stress test and the second is based on the economic projections from the Federal Reserve Board that were made public in December 2023. Finally, we conduct a structural analysis to investigate the effects of a financial shock on the recession probability, identified by recursive restrictions, while explicitly taking the ELB into account.
A key output from the proposed VAR is the latent states associated with the binary and censored variables. For example, from the { ProToVAR-het-L} model in Section (ref), we obtain two sets of latent variables corresponding to the NBER recession indicator and the short-term interest rates. The first set of estimates can be interpreted as a latent indicator of the business cycle and its turning points; the latent variables associated with the interest rates, or the shadow rates, provide a measure of the overall stance of monetary policy when the ELB is binding.
Figure (ref) presents the dynamic evolution of the latent business cycle indicator and the corresponding recession probabilities. Note that our implementation yields a counter-cyclical indicator. A pro-cyclical indicator can be readily obtained by re-defining the binary observed variable or adjusting Eq. ((ref)) to $\tilde{y}_{b,it} = \mathbb{I}(y_{b,it} < 0)$. As discussed in more detail in dueker2010forecasting, the distance of the latent indicator from $0$ provides a measure of how far the business cycle is from its next turning point. This aspect will also be relevant in interpreting the time-variation in GIRFs for recession probabilities that we discuss below. By this metric, the most severe recessions in our sample were the 1973--1975 one associated with the 1973 oil crisis, and the Great Recession between 2007--2009. Interestingly, the indicator suggests that the US economy was farthest from a recession immediately after the pandemic-induced turmoil in early 2020. This is due to the rebound of the economy in late 2020 and early 2021, which followed the severe contraction associated with the outbreak of the pandemic.
Figure (ref) report the shadow rates that arise when the interest rates hit the truncation threshold, in this case the ELB. In our sample there are two episodes in which many short-term interest rates hit the ELB. The first episode begins during the Great Recession and lasts till around 2016; the second corresponds to the COVID recession. As discussed in carriero2023shadow, the shadow rates may be viewed as the monetary policy reactions to the variables included in the VAR. In particular, they capture the overall stance of monetary policy when the ELB is binding and the central bank implements alternative policies---e.g., balance sheet policies and forward guidance on future rates---to provide stimulus.
The two episodes are very different in terms of their depth and duration. At the onset of the Great Recession, the monetary policy gradually loosens and the shadow rate associated with Fed funds rate reaches about $-2\%$ in 2021. It then gradually increases over the next few years, but does not reach zero until around 2016. In contrast, during the COVID recession, many shadow rates quickly reach unprecedented magnitudes, reflecting the historically loose monetary policy to support the drastic decline in employment and output. But the monetary policy tightens shortly after, when faced with subsequent inflationary pressures, and interest rates are above zero again by 2022.
\FloatBarrier
Conditional forecasting is a popular tool for assessing the impact on a set of variables of interest under different assumptions of the future paths of some other variables waggoner1999conditional, banbura2015conditional,antolin2021structural. For example, mccracken2022binary use it to examine the effects of monetary policy interventions and increases in oil prices on the likelihood of a recession. Below we employ the { ProToVAR} model to investigate the recession probabilities of two sets of scenarios considered by the Federal Reserve Bank. An advantage of our approach is that it takes into account the effect of the ELB on short-term interest rates in assessing the recession probabilities.
In the first exercise, we consider the projections made for the baseline and adverse scenarios in the 2020 Federal Reserve Board stress test for the period from the beginning of 2020 to the end of 2022, similar to chan2024conditional. More specifically, we assume the CPI inflation (CPIAUCSL), the unemployment rate (UNRATE), and 10-year government bond rates (GS10) follow the projections specified under the baseline and adverse scenarios (hard constraints), and assess the likelihood of a recession. These projections are shown in the panel (a) in Figure (ref), alongside the actual observations (solid dots). We compute the conditional forecasts under the baseline and adverse scenarios, as well as the unconditional forecasts (for which no knowledge of the future path of any variable is assumed), using data up to 2019M12. Panel (b) displays the conditional forecasts of the latent recession indicator and the implied recession probabilities, alongside the unconditional forecasts.
The conditional forecasts under the baseline scenario are rather similar to the unconditional forecasts, but with narrower credible bands as expected, and both indicate that a recession is unlikely. In sharp contrast, under the adverse scenario, the conditional forecasts project an almost certain recession starting at the beginning of 2020 that will last until the second half of 2021. Perhaps not surprisingly, the responses under the adverse scenario are somewhat different from the actual path of the US economy took after the extreme and unexpected COVID-19 shocks that occurred in the first half of 2020. In particular, under the adverse scenario, the rise and fall of unemployment rate is more gradual and the recovery takes much longer in comparison.
In our second set of conditional forecasts, we consider the scenario based on the Federal Open Market Committee Projections, particularly the “Economic Projections of Federal Reserve Board Members and Federal Reserve Bank Presidents," which were made public on December 13, 2023. We estimate the model { ProToVAR} using data up to 2023M12, and we restrict the future unemployment rate (UNRATE) and the Fed funds rate (FEDFUNDS) from 2024 to 2028 to fall within the central tendency ranges provided in the economic projections (soft constraints).\footnote{The ranges for central tendency exclude the three highest and three lowest projections for each variable in each year. In addition, since these projections are made available in the annual frequency, we set the projected ranges as end-of-year values and interpolate their monthly paths linearly.}
Panel (a) in Figure (ref) depicts the ranges for the unemployment rate and the Fed funds rate, as well as the unconditional forecasts from the model. Panel (b) plots the conditional forecasts of the latent recession indicator and the implied recession probabilities. The economic projections for the two variables coincide with the unconditional forecasts until around 2025, but after that the economic projections remain substantially lower than the unconditional forecasts. For example, the economic projections for the unemployment rate are concentrated around 4% starting from 2025, but the unconditional forecasts rise gradually and reach over 5% at the end of the forecast horizon. This more robust outlook for the labor market under the economic projections is consistent with a very low probability of a recession from 2025 to 2028, compared to the unconditional recession probabilities that can reach as high as 0.25.
\FloatBarrier
For this application, we investigate the dynamic effects of a financial shock on the recession probability. In contrast to previous studies, we include the NBER recession indicator in our VAR and explicitly take into account the ELB on nominal interest rates. More specifically, we use the large information set described in Section (ref) to estimate the model { ProToVAR-het}. We order the macroeconomic variables above the excess bond premium (EBP), followed by the recession indicator and all financial variables. Then, we identify a financial shock with zero impact restrictions via a Cholesky decomposition of the reduced-form covariance matrix and consider the orthogonalized innovation in the EBP equation as the structural financial shock. This identification scheme is similar to that used in gilchrist2012credit. This approach yields the impact vector $\bm{\delta}$ discussed in Section (ref), and the generalized impulse response functions can be computed as discussed there. We consider a positive and negative one standard deviation (SD) shock.\footnote{We have experimented with up to $3$ SD shocks, but asymmetries were muted. Simulating even larger shocks (e.g., such as the one during the Great Recession) may potentially yield stronger effects and/or more asymmetries.}
Figure (ref) presents the IRFs of selected variables to a one SD shock from the linear VAR specification of Eq. ((ref)). These variables are selected based on their economic importance, and roughly correspond to those shown in carriero2023shadow. The responses of the variables are in line with those found in the preceding literature. For example, industrial production and employment contract following the financial shock, reaching their troughs in 1-2 years after impact. The magnitude of the response of the latent recession indicator (NBERREC) is noteworthy; the impact of a $1$ SD shock translates to a contemporaneous increase in this indicator by about $0.2$. Depending on the state of the business cycle indicator (see Figure (ref) in the previous subsection), this may lead to substantial effects on recession probabilities during specific economic episodes.
We investigate this claim empirically in Figure (ref), which shows the peak GIRF at each point in time for positive/negative shocks. Intuitively, in light of the discussion above, financial shocks matter most when the business cycle is already close to a turning point (i.e., just before/during/after recessions, see Figure (ref)). For example, a one SD shock in early 2008 increases the recession probability by about 8 basis points, whereas it has virtually no impact if it occurs in, e.g., 2010. It is also interesting to note that our estimates point towards mostly symmetric effects for shocks of different signs (and different sizes which we investigated in unreported results).
Figure (ref) zooms into specific time periods and assesses the dynamic responses of recession probabilities across horizons. This is done in three ways. First, we plot the median response month-by-month, for each point in time (“Period medians,” for subsets of about $20$ years at a time to improve readability). Second, we show the posterior medians and 68 percent credible intervals for selected periods (“Periods”). Third, we compute variants of time averages (i.e., average partial effects, “Partial FX”). Specifically, we consider the average across all periods and an average for recession episodes.\footnote{They are computed for each MCMC draw so that we obtain draws from the posterior of these average partial effects for any desired (sub)set of periods.} It is clear that the responses to the financial shock depend on the phase of the business cycle. In general, the responses during expansion tend to be muted compared to those during recession. In particular, the strongest responses occurred during the 2000 dot-com bubble burst, followed by the period just before the start of the Great Recession.
\FloatBarrier
We extended the standard VAR to jointly model binary, censored and continuous variables and developed an efficient estimation approach that can be used to fit large datasets. We showed that the proposed VARs forecast recession and short-term interest rates well. And we demonstrated the usefulness of the proposed framework in a range of empirical applications.
The proposed framework can be further extended to include ordinal data (e.g., sovereign credit ratings) and count data (e.g., number of bank failures) using data augmentation. For example, zhang2015bayesian develop a model for ordinal and categorical variables in a static setting, whereas PPZ25 construct a VAR for count data. It would be useful to generalize the proposed framework to incorporate these approaches.
{\setstretch{0.87} \addcontentsline{toc}{section}{References} }\doublespacing