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.
70,800 characters · 14 sections · 69 citation commands
{8pt} {8pt}
In recent years, factor-augmented forecasting regression (FAR) models have received much attention in economic and finance following the seminal work of Stock01122002. They show that principal components extracted from a large set of predictors can substantially enhance forecasting accuracy of economic variables. Building on this, baing2006 establish the theory for the estimated parameters of the factor-augmented regressions to enable statistical inference. Following their work, researchers have conducted many theoretial and applied work for FAR. For instance, on the theoretical side, BaiAndNg2008 develop variable selection methods using thresholding rules to refine factor construction; XuCheng2015 propose a forecast combination approach within the FAR framework, assigning weights to individual models based on the Mallows and leave-h-out cross-validation criteria. Recently, Hannadige2024 propose the mixture-FAR method to incorporate both stationary and nonstationary factors. fanjianqing2024 propose FAR-NN model, which uses a diversified projection matrix to linearly compress observed variables to get a surrogate for latent factors and fits the regression function via deep ReLU networks. chenfanzhu2024 extend FAR from vector covariates to matrix covariates and introduce a matrix factor model structure. Other theoretical advances were made, among others, by wangcuili2015, PerronMcCracken, Karabiyik2018, Djogbenou2020 and TuWang2025modelselection. On the empirical side, Dijk2016 use factors extracted from a large set of macroeconomic variables to predict monthly U.S. excess stock returns. KellySeth2015 show that the factor-augmented approach can improve the conventional predictive regression forecasts in out-of-sample tests. We can also see the usefulness of factor-augmented models for forecasting finance and economics data from the discussions in Aastveit2015, Swanson2020 and Qiu2020.
It should be noticed that the existing statistical theory has primarily focused on factor-augmented forecasting regression with continuous response variables baing2006,XuCheng2015,Hannadige2024. However, an increasing number of studies have shown that binary response variable models have been widely applied in the fields of economics and finance, such as the prediction of bank failure audrino2018predicting, kim2022predicting, corporate bankruptcy, default, and financial distress li2020predicting,Charalambous2023, credit default liuzhangxiong2024 and business circlebybee2024business.
Despite the large number of potential applications, to the best of our knowledge the theoretical framework of the factor-augmented model with binary response variables is undeveloped. In this paper, we aim to fill this gap in the literature. Specifically, we consider the model:\footnote{Recent studies have also explored the incorporation of non-linear features and time-varying coefficients into factor-augmented regression models. For example, yancheng2022 introduce threshold effects into the factor-augmented framework; massacci2024 examine the forecasting performance of such models in the presence of structural breaks; Lietal2020 propose a nonlinear functional-coefficient model where regression coefficients vary smoothly with an index variable; and \citet*{chenhong2024} develop a time-varying forecast combination approach for factor-augmented regressions with smooth structural changes. In this paper, we focus on addressing the challenge of binary response variables within the factor-augmented framework, and we leave the exploration of more flexible model structures to future research.}
where $t=1,2,\cdots,T-h$, $h\geq0$ is the forecast horizon, $y_{t+h}$ is the binary dependent variable, $w_{t}$ is a $p\times1$ vector of observed regressors, $f_{t}$ is a $d_0\times1$ vector of unobserved factors, $(\beta_0,\beta_w^\top,\beta_f^\top)^\top$ are unknown parameters to be estimated, and $\epsilon_{t+h}$ is the idiosyncratic error term. In addition, we assume that both $p$ and $d_0$ are known and define $k=p+d_0+1$ as the total number of observable and unobservable factors. Suppose that we observe a large number of predictors $x_{it}$ that contain information about $f_{t}$, which is captured by a factor model as follows:
where $\lambda_{i}$ is the factor loadings and $e_{it}$ is another idiosyncratic error.
From a modelling standpoint, our framework is closely related to those of baing2006 and Hannadige2024, but with a key distinction that we consider a binary response variable. Unlike the continuous‐response framework, the binary nature of our dependent variable induces two fundamental deficiencies in their least‐squares approaches: (i) fitted probabilities can fall outside the unit interval, and (ii) the regression errors are inherently heteroskedastic. To remedy these shortcomings, we adopt a maximum‐likelihood framework and establish the resulting estimators’ consistency and asymptotic normality. Moreover, we relax some of the original assumptions by introducing $\alpha$-mixing conditions to capture the temporal dependence in the variables. Our model also shares structural similarities with that of gao2023. Both studies consider the binary response variable and MLE method for estimation. However, gao2023 focus on a panel data model in which the factors and model parameters are jointly estimated by MLE. Their specification relies on a fixed, low‐dimensional set of regressors and thus cannot accommodate a large number of covariates. In contrast, our approach applies principal component analysis (PCA) to a growing number of observed predictors, thereby enabling consistent forecasting in high‐dimensional settings.
From a practical perspective, our model is well-suited for classification tasks commonly encountered in economics and finance. Typical examples include applying high-dimensional features or large-scale factor structures for binary outcome prediction tasks, such as corporate bankruptcy forecasting in financial markets Jones2017highdimensional,JonesStewartJohnstone2017. Additionally, the model can effectively incorporate a wide range of high-dimensional variables---such as financial ratios, macroeconomic indicators, and transactional data---for credit risk assessment and credit scoring Laitinen1999,Masten2021,Silva2022,Liuyanghuangmazheng2024. In our empirical analysis, we focus on forecasting U.S. recessions. Using a large monthly macroeconomic dataset spanning from January 1960 to December 2024, we show that our proposed model consistently outperforms conventional Probit regression across both in-sample and out-of-sample exercises. In particular, our model delivers superior predictive accuracy at short- and medium-term horizons, effectively identifying major recessions such as the 2008 subprime crisis and the 2020 COVID-19 pandemic. These results demonstrate the practical relevance and applicability of the proposed model and estimation method, highlighting the value of incorporating high-dimensional information through latent factors in binary prediction tasks involving macroeconomic turning points.
In summary, the main contributions of the paper are as follows. (i) We develop a binary response factor-augmented forecasting regression that extends the original FAR framework baing2006,Hannadige2024 to accommodate binary response variables, thereby broadening its applicability to classification-oriented forecasting problems. (ii) We propose a maximum likelihood estimation method to estimate the unknown parameters and establish the corresponding asymptotic properties. In addition, different from baing2006, our theoretical results are derived under $\alpha$-mixing condition of $\{w_t, f_t, x_{it},\epsilon_{t},e_{it}\}$, which represents a more general dependence structure FanYao,Gao2007. (iii) We examine the finite-sample performance of the estimators through simulation studies and further use the proposed model to forecast U.S. recessions, highlighting its real-world relevance.
The structure of this paper is as follows. Section \hyperref[section-methods]{2} introduces the proposed model, develops the estimation methodology, and establishes the asymptotic theory. Section \hyperref[simulation]{3} conducts simulation experiments to examine the finite-sample performance of the theoretical findings. Section \hyperref[a_case_study]{4} provides an empirical study to demonstrate the model's practical utility by applying it to forecast U.S. recessions. Section \hyperref[conclusion]{5} concludes. The proof of the main theorems is provided in the appendix.
Before proceeding further, we introduce some mathematical symbols that will be used repeatedly throughout the paper. $\| \cdot \|$ denotes the Euclidean norm of a vector or the Frobenius norm of a matrix; for a vector or a matrix $A$, $A^\top$ denotes the transpose of $A$; $\mathrm{diag}(A_1, \ldots, A_n)$ represents a block-diagonal matrix with its diagonal elements are $A_1,\cdots, A_n$, and $I_m$ denotes an identity matrix of dimension $m$. In addition, $\xrightarrow{P}$ and $\xrightarrow{D}$ denote convergence in probability and in distribution, respectively. For a square matrix $M$, $\rho_{\min}(M)$ stands for its smallest eigenvalue; $a \wedge b = \min\{a, b\}$ and $a \vee b = \max\{a, b\}$.
In this section, we introduce our proposed estimation method and then establish the corresponding asymptotic theory.
Before proceeding further, we first introduce some additional notations. Let $\beta=(\beta_0,\beta_w^\top,\beta_f^\top)^\top$, $z_t=(1, w_t^\top, f_t^\top)^\top$, $X=(X_{1},\ldots, X_{N})^\top$ with $X_i=(x_{i1},\ldots, x_{iT})^\top$, $F = (f_1,\ldots, f_T)^\top$ is a $T\times d_0$ matrix consisting of factor regressors, $\Lambda=(\lambda_1,\lambda_2,\cdots, \lambda_N)^\top$ is a $N\times d_0$ matrix and $E$ denotes a $T\times N$ error matrix containing the elements $\{e_{it}\}$. With these notations, Model ((ref)) and ((ref)) can be rewritten as
Noteworthily, we assume the true number of factors $d_0$ is known in this step of estimation and will formally discuss the factor number selection issue in Section (ref).
In this paper, our object of interest is to estimate the unknown parameter vector $\beta$ in model ((ref)) and the conditional probability $\mathbb{P}(y_{t+h}=y \vert z_{t})$. To this end, we introduce a three-step estimation procedure, which is given as follows.
{\bf Step 1}: PCA estimation of unknown factors. We follow the methodology of bai2003 and apply principal component analysis to estimate the unknown factors. Let $\tilde{F}=(\tilde f_1,\ldots, \tilde f_T)^\top$ denote the estimated factor matrix. Specifically,
where $V_{NT}$ is a $d_0\times d_0$ diagonal matrix that contains the first $d_0$ largest eigenvalues of $\frac{1}{NT}XX^{\top}$ in descending order, $\Tilde{F}$ is a $T\times d_0$ matrix that contains the corresponding eigenvectors and satisfies $ \frac{1}{T}\tilde{F}^\top\tilde{F}=I_{d_0}$.
{\bf Step 2}: MLE estimation of $\beta$. Simple algebra yields
for $y=0$ or $1$, where $\Phi_{\epsilon}(\cdot)$ is the CDF of $\epsilon_{t}$. Drawing upon this fact, by replacing the unknown factors with their estimators, we obtain the factor augmented log-likelihood function:
where $\tilde{z_{t}}=(1,w_{t}^{\top},\tilde{f}_{t}^{\top})^{\top}$.
Therefore, we can construct the maximum likelihood estimator of $\beta$:
where we define $\widehat{\beta}=(\widehat{\beta}_0,\widehat{\beta}_{w}^{\top},\widehat{\beta}_f^{\top})^{\top}$.
{\bf Step 3}: Construction of predicted probabilities. Given the PCA estimates \(\tilde{F}\) and the MLE \(\widehat{\beta}\), we compute the predicted probability of the binary outcome: \[ \widehat{\mathbb{P}}(y_{t+h} = 1 \mid \tilde{z}_{t}) = \Phi_{\epsilon}(\widehat{\beta}^\top \tilde{z}_{t}). \]
In the first step of our procedure, we extract factor estimates from a large set of continuous predictor variables via the PCA method. Alternatively, if some of the predictors in $x_{it}$ are discrete, we can apply the MLE approach of wang2020 for nonlinear factor models to construct factor proxies. Specifically, we may consider the following nonlinear factor structure: \[ \mathbb{P}(x_{it}=x|\lambda_i^\top f_t)=g_{it}(x,\lambda_i^\top f_t), \] where \(g_{it}(x,\cdot)\) is some known joint probability mass function of $x_{it}$ given the factors and loadings. The factor estimators can then be obtained by maximizing the following log-likelihood function: \[ (\tilde{\Lambda}, \tilde{F})=\arg\max_{\Lambda,F}\sum_{i=1}^N \sum_{t=1}^T\log g_{it}(x_{it}, \lambda_i^\top f_t), \] Under suitable regularity conditions, wang2020 show that these maximum likelihood estimators achieve asymptotic properties comparable to those of the PCA estimators in linear factor models. In this regard, our proposed method can be extended to the case with discrete predictor variables by replacing the PCA step in the first stage with this nonlinear MLE-based factor estimation approach.
In the next subsection, we establish the asymptotic theory for $\widehat{\beta}$.
To facilitate the theoretical developments, we further define some notations. Let $E_t = (e_{1t},\ldots, e_{Nt})^\top$, $H=(\Lambda^\top\Lambda/N)(F^\top\tilde{F}/T)V_{NT}^{-1}$, $f_t^\circ=H^\top f_t$, $z^{\circ}_{t}=(1, w_{t}^\top, (H^\top f_{t})^\top)^\top$ and $\beta^{\circ}=(\beta_0,\beta_w^{\top},\beta_f^\top(H^\top)^{-1})^{\top}$. It is straightforward to see that $z_{t}^{\circ}=\text{diag}(I_{p+1},H^\top)z_t$, $\beta^\circ=\text{diag}(I_{p+1},H^{-1})\beta$ and $ \beta^\top z_t=\beta^{\circ\top} z_t^{\circ}$. Additionally, $M$ represents a sufficiently large constant that could be different at each appearance.
\hyperref[assumption_1]{Assumption 1} imposes $\alpha$-mixing conditions to control the serial and cross-sectional dependence among the variables. Such conditions are standard in the literature on factor models (see, for example, Assumption 3.1 in gao2023 and Assumption 1.(i) of GAO2020329). In addition, we assume strict exogeneity by requiring that the variables $w_t, \lambda_i$, and $f_t$ are independent of the random errors. It is worth noting that dependence between the observable and unobservable predictors is permitted.
With \hyperref[assumption_1]{Assumption 1}, we establish the consistency of the factor estimators obtained in the first step of our estimation procedure.
Lemma (ref) follows directly from $\frac{1}{T} \sum_{t=1}^{T} \left\| \tilde{f}_t - H^\top f_t \right\|^2=O_P\left(\frac{1}{N\wedge T}\right)$, which characterizes the rate of convergence for the factor estimators obtained via the PCA method. This rate is consistent with the seminal result of bai2003 and extends their framework to settings with $\alpha$-mixing dependence. Notably, Lemma (ref) provides the theoretical foundation for our augmented estimation of the coefficient $\beta$ in the second step of our procedure.
For notational convenience, let \( u_t = \beta^\top z_t \) and \( \widehat{u}_t = \widehat{\beta}^\top \tilde{z}_t \). Additionally, define the log-likelihood contribution at time \( t \) as $l_t(x) = y_{t+h} \log \Phi_\epsilon(x) + (1 - y_{t+h}) \log (1 - \Phi_\epsilon(x))$, where \( \Phi_\epsilon(\cdot) \) is the distribution function of the error term \( \epsilon \). Let \( l_t'(x) \), \( l_t''(x) \), and \( l_t'''(x) \) denote the first, second, and third derivatives of \( l_t(x) \) with respect to \( x \), respectively.
\hyperref[assumption_2]{Assumption 2} imposes standard conditions for maximum likelihood estimation. Specifically, \hyperref[assumption_2]{Assumption 2.1} requires that the innovation terms are identically distributed with a known distribution. This assumption can be relaxed by allowing the error distribution to be unknown, in which case its PDF and CDF can be estimated nonparametrically. \hyperref[assumption_2]{Assumption 2.2} imposes uniform bounds on the support of \(\beta^\top z_t\) with probability approaching one. Assuming bounded support is conventional in the MLE literature for nonlinear models wang2020. Additionally, \hyperref[assumption_2]{Assumption 2.3} and \hyperref[assumption_2]{2.4} regulate the derivatives of the error PDF and the individual log‐likelihood function, ensuring the existence and nonsingularity of the Hessian matrix.
With \hyperref[assumption_1]{Assumptions 1} and \hyperref[assumption_2]{2}, we establish the consistency of the augmented maximum likelihood estimators in the following theorem.
As established in \hyperref[Theorem_2.1_consistency]{Theorem 1}, \(\widehat{\beta}_0\) and \(\widehat{\beta}_w\) are consistent estimators of the true coefficients \(\beta_0\) and \(\beta_w\), respectively. Moreover, \(\widehat{\beta}_f\) consistently estimates \(\beta_f\) up to a rotation matrix \(H^{-1}\), which is a consequence of using the PCA method to recover the factors in the first stage. This indeterminacy is well known in the factor model literature bai2003, bai2009panel. By Theorem 1 of baing2002, one obtains \(\|H\| = O_P(1)\). Under additional identification conditions on the true factors (see Conditions 2.1–2.3 in BAI201318), the rotation matrix \(H\) converges in probability to the identity at rate $H - I \;=\; O_P\left(\frac{1}{\sqrt{N}\wedge \sqrt{T}}\right)$.
Before deriving the asymptotic distribution of \(\widehat{\beta}\), we impose additional conditions ensuring the positive definiteness of its asymptotic covariance matrix.
The positive‐definiteness requirements on \(\Sigma_\beta\) and \(\Omega_\beta\) in \hyperref[assumption_3]{Assumption 3} follow directly under the i.i.d.\ case. We now state the asymptotic distribution of \(\widehat{\beta}\) in the following theorem.
In \hyperref[Theorem_2.2_asymptotic]{Theorem 2}, we derive a central limit theorem for the augmented maximum likelihood estimator \(\widehat{\beta}\). The requirement \(\sqrt{T}/N \to 0\) controls the asymptotic order of the bias terms, ensuring that \(\widehat{\beta}\) achieves the optimal \(\sqrt{T}\)-rate of convergence. A similar condition appears frequently in the factor model literature for panel data (see, for example, Theorem 1 of baing2006).
In the forecasting regression literature, the asymptotic behaviour of fitted values is frequently examined. For instance, Theorem 3 of baing2006 establishes the limiting distribution for predicted outcome values in linear factor‐augmented regressions. In contrast, in binary response models, the discrete nature of \(y_{t+h}\) precludes directly forecasting the outcome itself. Accordingly, we focus on the asymptotic behaviour of the predicted probability \(\widehat{\mathbb{P}}(y_{t+h}=1\mid \tilde{z}_{t})=\Phi_{\epsilon}(\widehat{\beta}^\top \tilde z_t)\).
\hyperref[Theorem_2.3_asymptotic]{Theorem 3} demonstrates that the proposed factor‐augmented method consistently predicts the probability that the outcome falls into one (and consequently the complementary zero category) at the rate \( O_P\left(\frac{1}{\sqrt{N}\wedge\sqrt{T}}\right)\). Under some stronger regularity conditions, one can also derive the asymptotic distribution of this predicted probability. We leave this extension for future research.
For inference, deriving analytic estimators of the asymptotic covariance matrix \(\Omega_\beta\) becomes challenging in the presence of serial dependence. Accordingly, we follow higgins2025inference and employ the moving‐block bootstrap (MBB) to approximate the sampling distribution of our maximum‐likelihood estimators. Originally proposed by kunsch1989jackknife for the dependent time series data, the MBB has recently been adapted by higgins2025inference to accommodate panel‐data settings. This procedure randomly draws \(L\) overlapping time blocks with each of length \(q = \lfloor (T-h)/L\rfloor\) from the full sample and applies each block to all cross‐sectional units, where $\lfloor m\rfloor$ denotes the largest integer below $m$. Without loss of generality, we assume $T-h=qL$. We consider the following bootstrap procedure.
{\bf Step 1:} We draw $s_1, \ldots,s_L $ independently from the uniform discrete distribution on $\{0,\ldots,T-q-h\}$ and define \[ y_{(l-1)q+m+h}^\ast= y_{s_l+m+h}^\ast,\, w_{(l-1)q+m}^\ast= w_{s_l+m}^\ast,\, x_{i,(l-1)q+m}^\ast= x_{i,s_l+m}^\ast, \] for $l=1,\ldots,L$ and $m = 1,\ldots,q$.
{\bf Step 2:} Using the bootstrap sample \(\{y^{*}_{t+h},\,w^{*}_t,\,x^{*}_{i,t}\}\), re‐estimate the model by the augmented maximum‐likelihood procedure described in Section (ref), yielding the bootstrap estimator \(\widehat\beta^{*}\).
{\bf Step 3:} Repeat Steps 1 and 2 a large number of times to generate the empirical bootstrap distribution of \(\widehat\beta^{*}\).
Under standard regularity conditions, the bootstrap distribution of \(\widehat\beta^{*}\) consistently approximates the sampling distribution of the original estimator \(\widehat\beta\). We omit the formal proof, as it follows directly from existing results in the MBB literature.
In this section, we describe our approach to determining the number of factors. Selecting the factor dimension is a well‐studied problem in the factor‐model literature. For instance, baing2002 propose estimating the factor count by minimizing several information criteria, while ahn2013eigenvalue advocate choosing the number of factors as the maximizer of the ratio of adjacent eigenvalues. In this paper, we follow baing2002 and employ the information criterion method to select the factor dimension.
Specifically, we consider an information criterion of the following form:
where $\tilde{\lambda}_i^{(d)}$ and $ \tilde{F}_t^{(d)}$ represent the loading and factor estimators, respectively, when the number of factors is specified as $d$, and $C_{NT}=\frac{NT}{N+T}$.
Let \( d_{\max} \) be a prespecified upper bound on the number of factors. The true number of factors is then estimated as
Using arguments similar to those in the proof of Theorem 2 and Corollary 1 of baing2002, we can establish the consistency of this information criterion method: $\lim_{N,T\to\infty} \mathbb{P}(\widehat{d}=d_0)=1$, under \hyperref[assumption_1]{Assumption 1}.
Throughout the simulation and empirical analyses, we utilize this information criterion to select the number of factors whenever it is unknown.
In this section, we conduct simulations to examine the finite-sample performance of our proposed estimator. Specifically, we consider the following data generating process (DGP):
where $\beta_0=-2$, $\beta_f=(1,1)^\top$ and $\beta_w=(1,1)^\top$. The observable factors $w_t=(w_{1t}, w_{2t})^\prime$ are generated from uniform distributions $U(0,2)$ and $U(-3,3)$, respectively, and both the unobservable factors $f_t=(f_{1t}, f_{2t})^\top$ follow AR(1) processes:
for $t\geq1$, with $\rho_i=0.8^i$ for $i=1,2$ and $(\kappa_{1t},\kappa_{2t})$ are independently generated from the standard normal distribution. The initial values $f_{10}$ and $f_{20}$ are independently drawn from the uniform distribution $U(0,2)$.
Additionally, the DGP for observable predictors \{$x_{it}$\} is given by \[ x_{it}=\lambda_{i1} f_{1t}+\lambda_{i2} f_{2t}+\gamma_{it}, \] where the loadings are independently generated from the uniform distribution $U(0,6)$, and $\gamma_{it}$ is following the standard normal distribution. We denote the factor estimators as $(\tilde{f}_{1t},\tilde{f}_{2t})^\top$, where $t=1,2,...,T$.
We consider two distributions for the innovations $\epsilon_{t+1}$, corresponding to light-tailed and heavy-tailed distributions, respectively. In each example, we further vary the time-series correlation of the error term, covering the i.i.d. and autocorrelation scenarios to examine the robustness of the proposed method.
Example 1. Three normal distributions are considered in this example:
Example 2. We adopt DGPs analogous to those in {\bf Example 1}, replacing every normal distribution with the logistic distribution with zero mean and unit variance.
For each simulation dataset, we employ the method outlined in Section 2 to construct the estimators and predicted probabilities. We vary the values of $N$ and $T$ from $\{100,200,300\}$ and $\{100, 200,400\}$, respectively. After repeating the procedure for $R=500$ times, we compute the root mean squared errors (RMSEs) to evaluate the finite-sample estimation accuracy of the proposed method: \[ {\rm RMSE}_{all}=\sqrt{\frac{1}{R}\sum_{i=1}^{R}\| \widehat{\beta}_i-\beta^\circ\|^2},\, \] where $\widehat{\beta}_i$ denotes the estimated coefficient vector in the $i$th replication, and $\beta^\circ$ is the rotated coefficeints as defined in Section (ref). Additionally, We also report RMSE values for each individual coefficient, including RMSE$_{\text{cons}}$, RMSE$_{f_1}$, RMSE$_{f_2}$, RMSE$_{w_1}$, and RMSE$_{w_2}$.
We compute the area under the receiver operating characteristic curve (AUC) to assess predictive performance. The AUC measures the model's ability to distinguish between the two outcome classes (0 vs 1) across different threshold values. The receiver operating characteristic curve (ROC) curve maps the true positive rate
against the false positive rate
across all thresholds $0 \leq \xi \leq 1$, forming a function over the unit square $[0,1] \times [0,1]$. Here, $p_t$ denotes the predicted conditional probability of the target variable. And a higher AUC indicates stronger predictive discrimination. Specifically, for each replication, we calculate the AUC based on the estimated probabilities and the true binary outcomes. The simulation results are then summarized by reporting the mean, median, and standard deviation of the AUC across all replications.
Table (ref) and Table (ref) present the simulation results for the two Examples. The findings indicate that our method delivers robust performance under both light-tailed and heavy-tailed error distributions. The RMSE values decrease steadily as the sample size \(T\) increases. This pattern supports the \(\sqrt{T}\)-consistency of the augmented maximum likelihood estimators established in Theorem (ref).
Tables (ref) and (ref) summarize the mean, median, and standard deviation of the resulting AUCs. Overall, the proposed model exhibits high predictive accuracy, with both the mean and median AUC values exceeding 0.89 across simulations under various settings. Comparisons between Tables (ref) and (ref) reveal that variations in the distribution of idiosyncratic error terms have limited influence on the AUC results.
Recession forecasting is a critical task undertaken by many economic institutions. Predicting whether the economy will enter an expansion or recession in the coming months or years provides valuable insights for policymakers, investors, and households. For instance, government agencies can adjust fiscal policies based on projected recovery timelines, while central banks can refine monetary strategies in anticipation of future business cycle shifts. Seminal studies, such as estrella1998predicting,chauvet2005forecasting and kauppi2008predicting, demonstrate the predictive power of yield curve inversions and interest rate spread for U.S. recessions, using dynamic Probit models. An expanding body of research continues to identify additional predictive variables, including sentiment indicators christiansen2014forecasting, economic uncertainty measures ercolani2020forecasting, and credit market variables ng2012forecasting,ponka2017role. Meanwhile, researchers have increasingly focused on employing high-dimensional information and more advanced methods to predict U.S. recessions fornaro2016forecasting,berge2015predicting,liu2016predicts,vrontos2021modeling,proano2014predicting. Building on this literature, we conduct an empirical study to evaluate the proposed binary FAR model in predicting U.S. recession periods.
Following the same definition of recessions as estrella2006yield, we use the NBER defined business cycle expansion and contraction dates to determine US recessions. In particular, the first month following a peak month defines the first recession month and the last month of a through defines the last recession month. Let $y_t$ denote the state of the business cycle as determined by the NBER, where $y_t=1$ denotes that month $t$ is an NBER-defined recession and $y_t=0$ is an expansion. That is,
A crucial limitation arises from the NBER's ex post dating procedure, which introduces substantial reporting delays. We follow chauvet2005forecasting,kauppi2008predicting,ng2012forecasting and christiansen2014forecasting to assume that the recession indicator is announced with a three month delay, as they argue that the recent advances in real-time business cycle dating let researchers make reasonable assumptions on the state of the economy even if the NBER is yet to announce its classification. Specifically, we estimate model parameters using data measured at time $t - d - q$ (rather than $t - q$) for prediction construction, where $d=3$ denotes the assumed publication lag of the recession indicator to account for the announcement delay.
The NBER business cycle dates are publicly available and are standard in the business cycle literature.\footnote{The business cycle dates are available at \url{www.nber.org.}} We use a large, monthly frequency, macroeconomic dataset from which we estimate the monthly macroeconomic factors: 121 variables included in the FRED-MD database described in mccracken2016fred. The series represent broad categories of macroeconomic and financial time series including: output and income; labor market; housing; consumption, orders, and inventories; money and credit; interest and exchange rates; prices; and stock market. For detailed descriptions of these macroeconomic datasets, please refer to Michael W. McCracken's website\footnote{\url{https://research.stlouisfed.org/econ/mccracken/fred-databases/}}. We follow mccracken2016fred, using the transform precodure they offer to ensure the stationarity of the dataset. Our sample period spans from January 1960 to December 2024, providing a total of $T=780$ monthly observations. The sample period includes nine recession periods, with the most recent occurring during the COVID-19 pandemic from March to April 2020.
The $IC_2$ criterion baing2002 indicates that the optimal number of factors is 8. To better understand how these unobservable common factors relate to observable variables, we identify which series they best explain. Following mccracken2016fred, we regress each $i$-th series in the transformed dataset on a set of $k$ stationary common factors. This yields the adjusted $R_i^2(r)$ for each series $i$, where $r$ ranges from $1$ to $k$. Since the factors are orthogonal and ordered by decreasing eigenvalues, we can determine the incremental explanatory power of factor $r$ for series $i$ using the “marginal $R$-squared" ($mR^2$). This is calculated as $mR^2_i(r) = R^2_i(r) - R^2_i(r-1)$ for $r \ge 2$, with $mR^2_i(1) = R_i^2(1)$. The average importance of factor $r$ across all series is then $mR^2(r) = \frac{1}{N}\sum_{i=1}^N mR_i^2(r)$.
Table (ref) presents the $mR^2(r)$ values and lists the five series with the highest $mR^2_i(r)$ for each factor $r$. Typically, the observable series with the highest $mR^2_i(r)$ for a given factor are similar and belong to the same category. This allows us to assign economic meaning to each common factor based on its corresponding observable series. Furthermore, we select the observable series with the highest $mR^2(r)$ within its category to serve as an observable representative for each macroeconomic factor, which we refer to as the “observable factor" (bolded in Table (ref)). These eight factors collectively explain 52.3% of the total variation across all series. They are categorized as: two output factors, a prices factor, three interest factors, a labor factor, and a stock market factor. The corresponding observable factors are: IP: Manufacturing (SIC) index (IPMANSICS), CPI index: Commodities (CPIAUCSL), Baa-FF spread (BAAFFM), 1-year treasury rate (GS1), 5-year treasury spread (T5YFFM), average weekly working hours: manufacturing (AWHMAN), Real Personal Income (\textbf{RPI}), and S&P 500 index (\textbf{S&P 500}).
To have a clear view of the performance of the proposed binary FAR model, we consider two models in our empirical study.
In our empirical study, $w_t$ contains the eight observable factors in section (ref), that is, $w_t=(IPMANSICS_t,CPIAUCSL_t,BAAFFM_t,GS1_t,T5YFFM_t,AWHMAN_t,RPI_t,$ $S\&P 500_t)^\top$. All series are transformed to be stationary, for example, $S\&P 500_t$ denotes the log return of S&P 500 index. Given the limitations of traditional Probit models in handling high-dimensional settings, we assess the merits of the binary FAR model-based on latent factors—by contrasting its performance with that of a Probit regression relying on observable proxies for the same factors. Our analysis considers $h=1, 3, 6, 9, 12$, corresponding to 1-month, 1-quarter, 2-quarter, 3-quarter and 1-year-ahead predictions, respectively.
In this section, we use the full sample from 1960:M1 to 2024:M12 to evaluate the in-sample fitting performance of binary FAR and traditional Probit regression. Two goodness-of-fit measures are considered, the first is the area under the receiver operating characteristic curve (AUC), as introduced in section (ref). The second measure is the pseudo-$R^2$ of estrella1998new, which is a counterpart of the coefficient of determination ($R^2$) designed for binary response models, given by
where $\log L_u$ and $\log L_c$ are the maximum values of the constrained and unconstrained log-likelihood functions respectively, and $T$ is the sample size. This measure takes on values between 0 and 1, and can be interpreted in the same way as the coefficient of determination in the usual linear predictive regression model.
Table (ref) presents the in-sample performance of the binary FAR and Probit regression models across multiple forecast horizons ($h = 1, 3, 6, 9, 12$). Both models shows strong discriminatory power, as measured by the AUC, and moderate explanatory power, as indicated by the Pseudo-$R^2$. Across all horizons, binary FAR consistently outperforms Probit regression, with the performance gap being especially pronounced at shorter horizons. Figure (ref) plots the corresponding in-sample ROC curves, which confirm the findings in Table (ref). A curve lying above the 45$^\circ$ diagonal indicates classification performance better than random guessing. Notably, the binary FAR model demonstrates superior classification performance, while the Probit regression model, which relies solely on observable factors, shows weaker discriminatory capacity.
Figure (ref) contains the in-sample implied recession probability for binary FAR and Probit regression with NBER defined recession dates plotted in grey shading. The figure demonstrates that binary FAR consistently generates higher fitted probabilities during recession periods compared to Probit regression, particularly at shorter horizons. This superior performance is especially evident during the early 1990s and early 2000s recessions.
To sum up, the in-sample results demonstrate that the binary FAR model, by leveraging high-dimensional information in latent factors, provides significantly richer information about future business cycle periods than the conventional low-dimensional Probit regression.
It is well-known that in-sample results may not reliably reflect out-of-sample predictive accuracy. In order to evaluate the predictive accuracy of our models, we conduct an out-of-sample analysis. In particular, the full sample is split into an initial estimation period from 1960:M1 to 1999:M12 and an out-of-sample forecasting period from 2000:M1 to 2024:M12, yielding a total of $T = 300$ out-of-sample forecasts. The out-of-sample period includes the three major U.S. recessions: the dot-com crisis in 2001, the subprime crisis in 2008 and the recent COVID-19 pandemic in 2020. We employ an expanding window approach, whereby the estimation sample grows by one observation with each step forward in time. Since the Pseudo-$R^2$ is an in-sample measure, we use the out-of-sample AUC as the sole criterion to evaluate predictive performance.
Table (ref) reports the out-of-sample area under the ROC curve (AUC) for the binary FAR model and the traditional Probit regression across multiple forecast horizons ($ h = 1, 3, 6, 9, 12 $). The binary FAR model consistently outperforms the Probit regression, achieving AUC values above 0.88 at all horizons and near-perfect classshowsification performance at short-term horizons (0.982 at $ h=1 $ and 0.981 at $ h=3 $). In contrast, the Probit regression shows reasonable short-term predictive accuracy but suffers a substantial decline in performance as the forecast horizon increases, with AUC values falling to 0.778 at $ h=9$ and 0.787 at $h=12$.
Figures (ref) display the out-of-sample ROC curves for the binary FAR and Probit models across multiple horizons, confirming the findings from Table (ref). The binary FAR model consistently shows strong predictive accuracy, with curves for $h=1$ and $h=3$ closely approaching the top-left corner, and solid performance maintained even at longer horizons. In contrast, the Probit model exhibits weaker classification ability, particularly at $h=9$ and $h=12$, where the curves flatten toward the diagonal. These findings underscore the forecasting advantage of incorporating high-dimensional latent information through the binary FAR framework. By exploiting a richer information set embedded in the latent factors, the binary FAR model captures more persistent and informative signals related to future business cycle phases, particularly at medium- to long-term horizons.
Figure (ref) shows the out-of-sample recession probability forecasts for binary FAR and Probit regression with NBER defined recession dates plotted in grey shading. At shorter horizons ($h = 1, 3$), the binary FAR model effectively forecasts both the Dot-com crisis in 2001 and the COVID-19 pandemic in 2020. In contrast, the Probit regression fails to accurately capture the onset of the Dot-com crisis and exhibits two distinct probability peaks around 2020, resulting in a fragmented signal for the COVID-19 recession. At the longer horizon ($h = 12$), the binary FAR model continues to identify the COVID-19 episode, though with a lag of one to two months, whereas the Probit model produces a more diffuse and ambiguous forecast for the same period.
In summary, the out-of-sample results reinforce the superiority of the binary FAR model over the traditional Probit regression. By leveraging high-dimensional latent factors, binary FAR incorporates a richer and more comprehensive information set than the Probit model, which relies solely on a limited set of observable predictors. This informational advantage translates into consistently higher predictive accuracy across all horizons. Notably, the performance gap between the two models is even more pronounced in the out-of-sample setting than in-sample, highlighting the robustness and generalizability of the binary FAR model in capturing business cycle dynamics beyond the estimation sample.
In this paper, we have introduced a factor-augmented model with a binary response variable by allowing for both the cross-sectional dimension and the time dimension to diverge. We have developed an MLE-based procedure to estimate regression parameters and established the corresponding asymptotic theories. From a practical perspective, the framework can be applied to conduct bank failure prediction, credit default prediction, etc. Following estrella1998predicting,kauppi2008predicting,christiansen2014forecasting, in the empirical study, we focus on forecasting U.S. recessions. The proposed binary FAR model consistently outperforms conventional Probit regression across both in-sample and out-of-sample exercises, demonstrating the practical relevance and applicability of the proposed model and estimation method. Additionally, we conduct intensive simulation studies to examine the finite sample performance of the newly proposed model and methodology.
For future work, one may consider an extension to the binary response factor-augmented regression model with functional coefficients, time-varying parameters, structure breaks, or threshold effects gu2021autoencoder,yancheng2022,chenhong2024 to capture instable relationships between the predictors and outcome variables.