The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
76,960 characters
Large Bayesian VARs for Binary and Censored Variables
\maketitle
\vspace*{-4.5em}
\normalsize
\begin{center}
\begin{minipage}{.49\textwidth}
\large\centering Joshua C.C. \MakeUppercase{Chan}\\[0.25em]
\textit{Purdue University}
\end{minipage}
\begin{minipage}{.49\textwidth}
\large\centering Michael \MakeUppercase{Pfarrhofer}\\[0.25em]
\textit{WU Vienna}
\end{minipage}
\vspace*{2em}
{\large June 2025}
\end{center}
\vspace*{1em}
\doublespacing
\begin{center}
\begin{minipage}{0.85\textwidth}
\noindent\small
We extend the standard VAR to jointly model the dynamics of binary, censored and continuous variables, and develop an efficient estimation approach that scales well to high-dimensional settings. In an out-of-sample forecasting exercise, we show that the proposed VARs forecast recessions and short-term interest rates well. We demonstrate the utility of the proposed framework using a wide rage of empirical applications, including conditional forecasting and a structural analysis that examines the dynamic effects of a financial shock on recession probabilities.\\
{\sffamily \textbf{JEL}}: C34, C35, C53, E32, E47 \\
{\sffamily \textbf{KEYWORDS}}: macroeconomic forecasting, effective lower bound, financial shocks, shadow rate, recession
\end{minipage}
\end{center}
\singlespacing\vfill\noindent{\footnotesize\textit{Contact}: Michael Pfarrhofer ([email removed]), Department of Economics, WU Vienna University of Economics and Business.}
\thispagestyle{empty}\newpage
\doublespacing\normalsize
\section{Introduction}
Since the seminal work of \citet{Sims80} and \citet{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 \citep[e.g.,][]{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 \citet{chib1992bayes} and \citet{AC93}; see also \citet{anceschi2023bayesian} for a review of related approaches. More recently, \citet{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 \citet{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 \citet{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 \citet{dueker2006business}, \citet{dueker2010forecasting} and \citet{mccracken2022binary}. \citet{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, \citet{johannsen2021time} propose a multivariate unobserved components model for nominal interest rates that are bounded below by their ELB. Alternatively, \citet{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{s:framework} introduces the econometric framework and develops the estimation approach. Section~\ref{s:forecasting} 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{s:application} presents several applications to illustrate how the proposed framework can be used in different settings, including conditional forecasting and structural analysis. Finally, Section~\ref{s:conclusion} concludes and outlines some possible future research directions.
\section{Econometric Framework} \label{s:framework}
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.
\subsection{VAR for mixed variables} \label{sec:VAR}
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:
\begin{itemize}
\item $n_b$ binary variables $\tilde{y}_{b,it} \in \{0,1\}$ in $\tilde{\bm{y}}_{b,t} = (\tilde{y}_{b,1t},\hdots,\tilde{y}_{b,n_bt})'$;
\item $n_c$ censored variables $\tilde{y}_{c,it} \in [l_i,\infty)$ in $\tilde{\bm{y}}_{c,t} = (\tilde{y}_{c,1t},\hdots,\tilde{y}_{c,n_ct})'$ with known thresholds $l_i$ (for ease of exposition we take $l_i=0$; the proposed framework can be easily adapted to cases with non-zero thresholds, or when observations are censored from above);
\item $n_u$ unrestricted (continuous) variables $y_{u,it} \in \mathbb{R}$ in $\bm{y}_{u,t} = (y_{u,1t},\hdots,y_{u,n_ut})'$.
\end{itemize}
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
\begin{align}
\tilde{y}_{b,it} &= \mathbb{I}(y_{b,it} > 0), \quad \text{for } i = 1,\hdots,n_b,\label{eq:linkbinary}\\
\tilde{y}_{c,it} &= y_{c,it}\mathbb{I}(y_{c,it} > 0) = \max(y_{c,it},0), \quad \text{for } i = 1,\hdots,n_c,\label{eq:linkcensored}
\end{align}
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:
\begin{equation}
\bm{y}_t = \bm{a} + \bm{A}_1 \bm{y}_{t-1} + \hdots + \bm{A}_P \bm{y}_{t-P} + \bm{\epsilon}_t, \quad \bm{\epsilon}_t \sim \mathcal{N}(\bm{0}_n, \bm{\Sigma}_t),\label{eq:VAR}
\end{equation}
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 \citet{CCM16}, heavy-tailed error distributions considered in \citet{chan20}, an outlier component proposed in \citet{SW16}, and a volatility break model for fitting post COVID-19 pandemic data developed in \citet{LP22}; see also \citet{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 \citet{dueker2005dynamic} that contains a single binary variable; when $n_b = 0$, it then becomes what has been termed a \textit{simple} shadow rate VAR by \citet{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. \citet{zhang2015bayesian} consider some of these extensions in a static setting; \citet{PPZ25} propose a related framework for count data.}
\subsection{The conditional distribution of the latent variables }\label{sec:lik}
If $\bm{y}_{b,t}$ and $\bm{y}_{c,t}$ were observed, it would be relatively straightforward to estimate the VAR in~\eqref{eq:VAR} even when $n$ is large---one can simply employ any recently developed efficient posterior sampler designed for large VARs with complete data \citep[e.g.,][]{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 \citet{chan2023high} for drawing missing data designed for conditionally Gaussian state space models to our setting. Specifically, \citet{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 \texttt{l} and \texttt{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$:
\begin{itemize}
\item \textbf{Missing} values (labeled \texttt{m}, irrespective of variable type). Let $n_{b,t}^{\text{\texttt{m}}}$ denote the number of binary variables that are missing at time $t$; similarly define $ n_{c,t}^{\text{\texttt{m}}}$ and $ n_{u,t}^{\text{\texttt{m}}}$. Let $\bm{y}_{\text{\texttt{m}}, t}$ represent the $(n_{b,t}^{\text{\texttt{m}}} + n_{c,t}^{\text{\texttt{m}}} + n_{u,t}^{\text{\texttt{m}}})\times1$ subvector of $\bm{y}_t$ consisting of all missing values at time $t$. Then, define the corresponding selection matrix $\bm{s}_{\text{\texttt{m}},t}'$ of size $(n_{b,t}^{\text{\texttt{m}}} + n_{c,t}^{\text{\texttt{m}}} + n_{u,t}^{\text{\texttt{m}}})\times n$ such that $\bm{s}_{\text{\texttt{m}},t}'\bm{y}_t = \bm{y}_{\text{\texttt{m}},t}$.
\item \textbf{Restricted} values (labeled \texttt{r}, for observed binary variables and observations where censoring applies). Here, we use the selection matrix $\bm{s}_{\text{\texttt{r}},t}'$ of size $(n_{b,t}^{\text{\texttt{r}}} + n_{c,t}^{\text{\texttt{r}}})\times n$ to select the \textit{observed} binary and truncated variables, such that $\bm{s}_{\text{\texttt{r}},t}'\bm{y}_t = \bm{y}_{\text{\texttt{r}},t}$.
\item \textbf{Observed} unrestricted values (labeled \texttt{o}). Again, we define the selection matrix $\bm{s}_{\text{\texttt{o}},t}'$ of size $(n_{c,t}^{\text{\texttt{o}}} + n_{u,t}^{\text{\texttt{o}}})\times n$ such that $\bm{s}_{\text{\texttt{o}},t}'\bm{y}_t = \bm{y}_{\text{\texttt{o}},t}$.
\end{itemize}
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}}}$:
\begin{equation}
\bm{y}_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},\quad \bm{y}_t = \bm{s}_{\text{\texttt{l}},t} \bm{y}_{\text{\texttt{l}},t} + \bm{s}_{\text{\texttt{o}},t} \bm{y}_{\text{\texttt{o}},t},\label{eq:decompY}
\end{equation}
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. \eqref{eq:decompY}, 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
\begin{equation} \label{eq:y_partition}
\bm{y} = \bm{S}_{\text{\texttt{l}}} \bm{y}_{\text{\texttt{l}}} + \bm{S}_{\text{\texttt{o}}} \bm{y}_{\text{\texttt{o}}}.
\end{equation}
Following \citet{chan2023high}, we encode the VAR dynamics in the $Tn \times (T+P)n$ matrix
\begin{equation*}
\bm{M} = \begin{bmatrix}
-\bm{A}_p & \cdots & -\bm{A}_1 & \bm{I}_n & \bm{0}_{n\times n} & \cdots & \cdots & \bm{0}_{n\times n}\\
\bm{0}_{n\times n} & -\bm{A}_p & \cdots & -\bm{A}_1 & \bm{I}_n & \bm{0}_{n\times n} & \cdots & \bm{0}_{n\times n}\\
\vdots & \ddots & \ddots & & \ddots & \ddots & \ddots & \vdots\\
\bm{0}_{n\times n} & \cdots & \bm{0}_{n\times n} & -\bm{A}_p & \cdots & -\bm{A}_1 & \bm{I}_n & \bm{0}_{n\times n}\\
\bm{0}_{n\times n} & \cdots & \cdots & \bm{0}_{n\times n} & -\bm{A}_p & \cdots & -\bm{A}_1 & \bm{I}_n
\end{bmatrix},
\end{equation*}
so that Eq. \eqref{eq:VAR}, for $t = 1,\ldots,T,$ can then be written compactly as:
\begin{equation}
\bm{M}\bm{y} = \bm{m} + \bm{\epsilon}, \quad \bm{\epsilon}\sim\mathcal{N}(\bm{0}_{Tn},\bm{\Omega}),\label{eq:static}
\end{equation}
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. \eqref{eq:y_partition} into Eq. \eqref{eq:static}, 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 \citet{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:
\begin{equation}
\bm{V}_{\text{\texttt{l}}} = (\bm{H}_{\text{\texttt{l}}}'\bm{\Omega}^{-1}\bm{H}_{\text{\texttt{l}}})^{-1}, \quad \bm{\mu}_{\text{\texttt{l}}} = \bm{V}_{\text{\texttt{l}}}\bm{H}_{\text{\texttt{l}}}'\bm{\Omega}^{-1}(\bm{m} - \bm{H}_{\text{\texttt{o}}}\bm{y}_{\text{\texttt{o}}}).\label{eq:poststates}
\end{equation}
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:
\begin{equation}
\bm{y}_{\text{\texttt{l}}}\:\vert\:(\tilde{\bm{y}}_b,\tilde{\bm{y}}_c,\bm{y}_{\text{\texttt{o}}},\bullet) \sim \mathcal{N}(\bm{\mu}_{\text{\texttt{l}}},\bm{V}_{\text{\texttt{l}}}), \quad \text{subject to} \quad \bm{R}\bm{y}_{\text{\texttt{l}}} \geq \bm{0}, \label{eq:truncGauss}
\end{equation}
where $\bm{R}$ is a diagonal matrix incorporating the inequality restrictions in Eqs. (\ref{eq:linkbinary}) and (\ref{eq:linkcensored}). 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 \citep[e.g.,][]{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{eq:truncGauss}) 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 \citep[see, e.g.,][]{pakman2014exact} to sample from the high-dimensional truncated normal distribution, as discussed in the next section.
\subsection{Priors, posteriors and sampling algorithms}
In this section, we outline the priors on the model parameters and provide details of key sampling steps in the posterior simulation.
\noindent\textsc{\textbf{Sampling the latent states}}. As discussed in Section \ref{sec:lik}, 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{eq:truncGauss}). Here we implement an HMC step, and specifically versions referred to as harmonic HMC or Zigzag HMC \citep[see][]{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 \citet{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 \citet{nishimura2020discontinuous,nishimura2021hamiltonian,zhang2022hdtg} for more discussion.
Our default implementation uses Zigzag HMC because we may compute the mean in Eq. \eqref{eq:poststates} 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{s:simulations} provides some simulation results that show that the proposed algorithm works well.
\noindent\textsc{\textbf{Restrictions on the covariance matrix}}. As noted in Section~\ref{sec:VAR}, 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 \citep{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:
\begin{equation}
p(\bm{\Sigma},\bm{d} \:\vert\:\bm{y}, \bullet) \propto p(\bm{\Sigma},\bm{d}) p(\bm{y}\:\vert\:\bullet)
= \text{det}(\bm{D})^{(n-1)/2} p_{\mathcal{W}^{-1}}(f(\bm{\Sigma},\bm{d}))p(\bm{y}\:\vert\:\bullet).\label{eq:condlikiW}
\end{equation}
This distribution is not of a known form, and we follow the parameter-expanded Metropolis-Hastings (PX-MH) algorithm proposed in \citet{zhang2020parameter} to sample from it.\footnote{\citet{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.}
\noindent\textsc{\textbf{Priors 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 \citep[HS,][]{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{sec:VAR}, 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 \citet{SW16}, which is used in a VAR setting for modeling COVID-19 outliers in \citet{carriero2021addressing}. In particular, we assume that:
\begin{equation*}
o_t = \begin{cases}
1 & \text{with probability} \quad 1 - \mathfrak{p}\\
\mathcal{U}(2,10) & \text{with probability} \quad \mathfrak{p},
\end{cases}
\end{equation*}
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 \citet{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.
\subsection{Forecasts, scenario analysis and structural inference}\label{sec:forecasting}
The VAR$(P)$ process for the latent and observed variables, as defined in Eq. (\ref{eq:VAR}), 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:
\begin{equation}
p(\bm{y}_{\tau+1:\tau+h}|\bm{\Xi},\mathcal{I}_{1:\tau}) \sim\mathcal{N}(\bm{\mu}_{y,\tau h},\bm{V}_{y,\tau h}),\label{eq:preddist}
\end{equation}
with known moments that can be derived from considerations analogous to Eq. (\ref{eq:static}) in Section \ref{sec:lik}; see also \citet{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{eq:linkbinary}) and (\ref{eq:linkcensored}). 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 \citet{mccracken2022binary}. We may compute this probability for each sweep of our sampler which provides us with draws from the corresponding posterior.
\noindent\textsc{\textbf{Conditional 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{eq:preddist}) 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 \citep[e.g.,][]{antolin2021structural} can be directly applied using the latent representation given in Eq. (\ref{eq:VAR}).
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 \citet{chan2024conditional} for computational details.
\noindent\textsc{\textbf{Generalized 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{eq:linkbinary}) and (\ref{eq:linkcensored}) require the use of generalized IRFs \citep[GIRFs, see][]{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 \citet{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{eq:preddist}) 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{eq:linkcensored}) 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 \textit{average partial effect}. By contrast, one may also compute a \textit{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
\section{Out-of-Sample Forecasting} \label{s:forecasting}
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 \citet{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 \citet{giusto2017identifying} and \citet{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 ({\sffamily C}), small ({\sffamily S}), medium ({\sffamily M}), and large ({\sffamily L}). Table \ref{tab:data} 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 \citet{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 \citet{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.
\begin{table}[t]
\caption{List of variables used in the forecasting exercise. \textit{Notes}: Variable codes, descriptions and transformation: (0) no transformation $h(x_t) = x_t$; (1) logarithms $h(x_t) = \log(x_t)$; (2) log-differences $h(x_t) = 100\cdot\log(x_t/x_{t-1})$. The check marks in the remaining columns indicate inclusion in the respective information set.}\label{tab:data}
\centering\footnotesize
\begin{tabular}{llcllll}
\toprule
\textbf{Code} & \textbf{Description} & $h(x_t)$ & {\sffamily C} & {\sffamily S} & {\sffamily M} & {\sffamily L} \\
\midrule
NBERREC & NBER Recession Indicator & 0 & \checkmark & \checkmark & \checkmark & \checkmark \\
FEDFUNDS & Effective Federal Funds Rate & 0 & \checkmark & \checkmark & \checkmark & \checkmark \\
TB3MS & 3-Month Treasury Bill & 0 & & \checkmark & \checkmark & \checkmark \\
TB6MS & 6-Month Treasury Bill & 0 & & \checkmark & \checkmark & \checkmark \\
GS1 & 1-Year Treasury Rate & 0 & & \checkmark & \checkmark & \checkmark \\
GS5 & 5-Year Treasury Rate & 0 & & \checkmark & \checkmark & \checkmark \\
GS10 & 10-Year Treasury Rate & 0 & & \checkmark & \checkmark & \checkmark \\
TS10G1 & Term-spread 10-Year minus 1-Year Treasury Rate & 0 & \checkmark & \checkmark & \checkmark & \checkmark \\
PAYEMS & All Employees: Total non-farm & 2 & \checkmark & \checkmark & \checkmark & \checkmark \\
INDPRO & Industrial Production Index & 2 & \checkmark & \checkmark & \checkmark & \checkmark \\
CMRMTSPLx & Real Manuf. and Trade Industries Sales & 2 & \checkmark & \checkmark & \checkmark & \checkmark \\
W875RX1 & Real personal income ex transfer receipts & 2 & \checkmark & \checkmark & \checkmark & \checkmark \\
DPCERA3M086SBEA & Real personal consumption expenditures & 2 & & & & \checkmark \\
CUMFNS & Capacity Utilization: Manufacturing & 1 & & & & \checkmark \\
UNRATE & Civilian Unemployment Rate & 0 & & & & \checkmark \\
CES0600000007 & Avg Weekly Hours: Goods-Producing & 1 & & & & \checkmark \\
CES0600000008 & Avg Hourly Earnings: Goods-Producing & 2 & & & & \checkmark \\
PCEPI & Personal Consumption Expenditure (PCE) Prices & 2 & & & & \checkmark \\
CPIAUCSL & Consumer Price Index (CPI): All Items & 2 & & & & \checkmark \\
HOUST & Housing Starts: Total New Privately Owned & 1 & & & & \checkmark \\
SP500 & S\&Ps Common Stock Price Index & 2 & & & \checkmark & \checkmark \\
EXUSUKx & US/UK Foreign Exchange Rate & 2 & & & \checkmark & \checkmark \\
OILPRICEx & Crude Oil, spliced WTI and Cushing & 2 & & & \checkmark & \checkmark \\
BAA & Moody's Seasoned Baa Corporate Bond Yield & 0 & & & \checkmark & \checkmark \\
EBP & Excess Bond Premium & 0 & & & \checkmark & \checkmark \\
\bottomrule
\end{tabular}
\end{table}
Our most general model is called the {\sffamily ProToVAR}, because it features both the binary target variable (addressed by the {\sffamily Pro}bit part) and censored variables (the {\sffamily To}bit part). Another specification includes the binary target variable but treats all other variables as continuous. It thus only features the {\sffamily Pro}bit part, and we label it the {\sffamily 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 {\sffamily hom}) and the other with the outlier component that captures heteroskedasticity (labeled {\sffamily 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.}
\subsection{Recession forecasts}
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{fig:recprob}. 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 {\sffamily ProVAR-hom} specification estimated using the core dataset coincides with the preferred model of \citet{mccracken2022binary}.
\begin{figure}[ht]
\centering
\includegraphics[width=\textwidth]{figures/auc.pdf}
\caption{Area under ROC curve (AUC). Higher values indicate superior performance of the respective specification. Rows refer to subsamples, colored squares indicate the best performing specification by row.}
\label{fig:recprob}
\end{figure}
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, {\sffamily ProToVAR-het} performs best for dating recessions; {\sffamily ProVAR-het} forecasts well for horizons less than a year; and {\sffamily ProVAR-hom} tends to be the best model for horizons more than a year. But for post 2000 subsample, {\sffamily ProToVAR-het} and {\sffamily ProVAR-het} tend to outperform their homoskedastic counterparts.
Figure \ref{fig:prob-dist} plots the out-of-sample recession probabilities across multiple horizons from {\sffamily 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.
\begin{figure}[ht]
\centering
\includegraphics[width=\linewidth]{figures/fc_probs.pdf}
\caption{Recession probabilities from {\sffamily ProToVAR-het} for horizons $h= -6, 0, 12$ estimated using the core (upper panel) and large datasets (lower panel). Predictive median alongside 68/90 percent credible sets.}
\label{fig:prob-dist}
\end{figure}
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.
\subsection{Interest rate forecasts}
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 {\sffamily VAR}. The second model treats the interest rates as censored variables, and is referred to as {\sffamily ToVAR}. We use the continuous ranked probability score \citep[CRPS][]{gneiting2007strictly} as a density forecast metric.
Table \ref{tab:crps_ir} presents the CRPSs for different models ({\sffamily VAR}, {\sffamily ProVAR}, {\sffamily ToVAR} and {\sffamily ProToVAR}), with and without heteroskedasticity ({\sffamily het} and {\sffamily hom}), estimated using different datesets (small, medium and large). We use the homoskedastic VAR estimated with the large dataset ({\sffamily VAR-hom-L}) as the benchmark model. The values in the row associated with {\sffamily 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.
\begin{table}[t]
\caption{Predictive accuracy for interest rates as measured by average continuous ranked probability scores (CRPSs) for the post 2000 subsample. Values are ratios benchmarked relative to {\sffamily VAR-hom-L} (the gray-shaded row shows raw losses), best performing specification for each horizon in bold.}\label{tab:crps_ir}
\includegraphics[width = \textwidth]{tables/crps_ir.pdf}
\end{table}
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 {\sffamily 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 {\sffamily VAR} models sometimes produce longer-horizon forecasts that are below the ELB. Moreover, it is interesting to note that {\sffamily ProToVAR-het} models consistently outperform their corresponding {\sffamily 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, {\sffamily 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{app:results}. 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, {\sffamily ProToVAR-het} is always competitive, and in several cases even yields improvements in predictive accuracy for these variables as well.
\FloatBarrier
\section{Applications} \label{s:application}
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.
\subsection{Recession probabilities and shadow rates}
A key output from the proposed VAR is the latent states associated with the binary and censored variables. For example, from the {\sffamily ProToVAR-het-L} model in Section~\ref{s:forecasting}, 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{fig:recession_full-sample} 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{eq:linkbinary}) to $\tilde{y}_{b,it} = \mathbb{I}(y_{b,it} < 0)$. As discussed in more detail in \citet{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.
\begin{figure}[ht]
\centering
\includegraphics[width=\linewidth]{figures/ProToVARrecession.pdf}
\caption{Estimates of the latent states associated with the binary NBER recession indicator and the corresponding recession probabilities for the sampling period from 1959M2 to 2022M12.}
\label{fig:recession_full-sample}
\end{figure}
\begin{figure}[ht]
\centering
\includegraphics[width=\linewidth]{figures/ProToVARshadowrates.pdf}
\caption{Estimates of the partially latent states associated with the interest rates subject to truncation at the ELB from {\sffamily ProToVAR-het-L}.}
\label{fig:shadowrates}
\end{figure}
Figure \ref{fig:shadowrates} 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 \citet{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
\subsection{Conditional forecasts and scenario analysis}
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 \citep{waggoner1999conditional, banbura2015conditional,antolin2021structural}. For example, \citet{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 {\sffamily 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 \citet{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{fig:cf_H2019restricted}, 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.
\begin{figure}[ht]
\centering
\begin{subfigure}{\textwidth}
\centering
\caption{Hard restrictions}
\includegraphics[width=\linewidth]{figures/cf/cf_restr_2019-12-01.pdf}
\end{subfigure}
\begin{subfigure}{\textwidth}
\centering
\caption{Conditional recession forecast}
\includegraphics[width=\linewidth]{figures/cf/cf_restr_2019-12-01_NBERREC.pdf}
\end{subfigure}
\caption{Conditional forecasts and restrictions for the \texttt{2019-12-01} forecast origin. Solid dots indicate realized values.}
\label{fig:cf_H2019restricted}
\end{figure}
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 {\sffamily 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{fig:cf_S2023restricted} 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.
\begin{figure}[ht]
\centering
\begin{subfigure}{\textwidth}
\centering
\caption{Soft restrictions}
\includegraphics[width=\linewidth]{figures/cf/cf_restr_2023-12-01.pdf}
\end{subfigure}
\begin{subfigure}{\textwidth}
\centering
\caption{Conditional recession forecast}
\includegraphics[width=\linewidth]{figures/cf/cf_restr_2023-12-01_NBERREC.pdf}
\end{subfigure}
\caption{Conditional forecasts and restrictions for the \texttt{2023-12-01} forecast origin. Solid dots indicate realized values.}
\label{fig:cf_S2023restricted}
\end{figure}
\FloatBarrier
\subsection{Financial shocks in the US}
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{s:forecasting} to estimate the model {\sffamily 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 \citet{gilchrist2012credit}. This approach yields the impact vector $\bm{\delta}$ discussed in Section \ref{sec:forecasting}, 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.}
\begin{figure}[ht]
\centering
\includegraphics[width=\linewidth]{figures/irf/default_irf.pdf}
\caption{Impulse response functions for selected latent and observed variables. Cumulated responses for differenced variables. Posterior median (solid black) alongside 68/90 percent posterior credible sets (gray shades).}
\label{fig:default_irf}
\end{figure}
Figure \ref{fig:default_irf} presents the IRFs of selected variables to a one SD shock from the linear VAR specification of Eq. (\ref{eq:VAR}). These variables are selected based on their economic importance, and roughly correspond to those shown in \citet{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{fig:recession_full-sample} 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{fig:peak_irf}, 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{fig:recession_full-sample}). 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).
\begin{figure}
\centering
\includegraphics[width=\linewidth]{figures/irf/peak_irf.pdf}
\caption{Peak generalized impulse response functions for recession probabilities over time, peak horizon $h$ in shades of blue and red. Dots mark the posterior median alongside 68 percent credible set error bars. Gray shades indicate NBER recessions.}
\label{fig:peak_irf}
\end{figure}
Figure \ref{fig:prob_girf} 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.
\begin{figure}
\centering
\includegraphics[width=\linewidth]{figures/irf/prob_girf.pdf}
\caption{Generalized impulse response functions for recession probabilities. The left panels plot the posterior median response for each period (``Period medians''). The panel ``Periods'' shows the responses for selected periods alongside 68 percent credible sets, the rightmost panel ``Partial FX'' shows partial effects averaged across subsets of periods.}
\label{fig:prob_girf}
\end{figure}
\FloatBarrier
\section{Concluding remarks and future research} \label{s:conclusion}
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, \citet{zhang2015bayesian} develop a model for ordinal and categorical variables in a static setting, whereas \citet{PPZ25} construct a VAR for count data. It would be useful to generalize the proposed framework to incorporate these approaches.
{\setstretch{0.87}
\small\addcontentsline{toc}{section}{References}
\bibliographystyle{custom.bst}
\bibliography{lit}}\normalsize\clearpage\doublespacing