Inside the black box: Neural network-based real-time prediction of US
  recessions
Long short-term memory (LSTM) and gated recurrent unit (GRU) are used to
model US recessions from 1967 to 2021. Their predictive performances are
compared to those of the traditional linear models. The out-of-sample
performance suggests the application of LSTM and GRU in recession forecasting,
especially for longer-term forecasts. The Shapley additive explanations (SHAP)
method is applied to both groups of models. The SHAP-based different weight
assignments imply the capability of these types of neural networks to capture
the business cycle asymmetries and nonlinearities. The SHAP method delivers key
recession indicators, such as the S&P 500 index for short-term forecasting up
to 3 months and the term spread for longer-term forecasting up to 12 months.
These findings are robust against other interpretation methods, such as the
local interpretable model-agnostic explanations (LIME) and the marginal
effects.
Seulki Chung
# Inside the black box: Neural network-based real-time prediction of US recessions

###### Abstract

Feedforward neural network (FFN) and two specific types of recurrent neural network, long short-term memory (LSTM) and gated recurrent unit (GRU), are used for modeling US recessions in the period from 1967 to 2021. The estimated models are then employed to conduct real-time predictions of the Great Recession and the Covid-19 recession in US. Their predictive performances are compared to those of the traditional linear models, the logistic regression model both with and without the ridge penalty. The out-of-sample performance suggests the application of LSTM and GRU in the area of recession forecasting, especially for the long-term forecasting tasks. They outperform other types of models across 5 forecasting horizons with respect to different types of statistical performance metrics. Shapley additive explanations (SHAP) method is applied to the fitted GRUs across different forecasting horizons to gain insight into the feature importance. The evaluation of predictor importance differs between the GRU and ridge logistic regression models, as reflected in the variable order determined by SHAP values. When considering the top 5 predictors, key indicators such as the S&P 500 index, real GDP, and private residential fixed investment consistently appear for short-term forecasts (up to 3 months). In contrast, for longer-term predictions (6 months or more), the term spread and producer price index become more prominent. These findings are supported by both local interpretable model-agnostic explanations (LIME) and marginal effects.

keywords: Forecasting, Recession, Business cycle, LSTM, GRU, SHAP, LIME

_JEL:_ C01, C45, C51, C52, C53, C71, E37
## 1 Introduction

Recession forecasting has been a longstanding challenge for policymakers and market practitioners, as it enables them to make timely decisions that could mitigate the impact of a recession. However, due to the intricate and interconnected nature of modern economics, this task has proven to be difficult and has had limited success. Traditional methods such as probit or logit regression models, along with their variations, have been widely employed to address this forecasting issue. Conversely, in recent decades, machine learning techniques and artificial neural networks have become increasingly popular among economists and have been applied to certain macroeconomic forecasting problems. Nevertheless, it has yet to be demonstrated whether these modern approaches outperform traditional statistical models when it comes to predicting recessions.

Since the pioneering work of Mitchell and Burns (1938), which identified 21 variables as potential economic indicators for business cycles out of a larger set, researchers have been engaged in selecting indicators and developing theoretical frameworks or predictive models to link these indicators with business cycles. However, most studies have primarily focused on linear frameworks. Many of them have utilized probit regression models or their extensions to generate recession forecasts. Estrella and Mishkin (1996) and Estrella and Mishkin (1998) compiled a combination of financial and macroeconomic variables and conducted recession forecasting using a probit framework. Their findings revealed that stock prices exhibit greater short-term predictability, while the slope of the yield curve performs better for longer-term predictions. Wright (2006) demonstrated that probit models incorporating the federal funds rate and term spread as predictors outperform models solely relying on the term spread. Dueker (1997) and Dueker (2002) expanded the standard probit model by incorporating Markov regime switching within the probit framework, allowing for variation in coefficients. Chauvet and Potter (2005) introduced several specifications to the probit model that accounted for different business cycle dependencies and autocorrelation errors, concluding that their more complicated extensions improved the accuracy of recession forecasts. Fornari and Lemke (2010) as well as Nyberg (2014) incorporated vector autoregression components into the probit model to capture the endogenous dynamics of the predictors.

While probit regression and its extensions continue to be widely used in business cycle forecasting, there has been a growing interest in exploring the forecasting capabilities of nonlinear models, including machine learning methods. This interest stems from the recognition that macroeconomic
data often exhibit nonlinear patterns. Empirical evidence supports this notion, such as the work of Tiao and Tsay (1994), who demonstrate that a threshold autoregressive model outperforms a linear autoregressive model for predicting GDP growth. Maasoumi et al. (1994) examine multiple macroeconomic time series and confirm their nonlinear nature. Puglia and Tucker (2021) highlight the attractiveness of machine learning methods as an alternative to probit regression and its extensions, noting that probit methods typically require additional parameters for flexibility, whereas flexibility is inherent in machine learning methods. Stock and Watson (1998) compare the forecasting performance of 49 univariate linear and nonlinear models across 215 macroeconomic time series. They find that some of the nonlinear models perform poorly compared to linear models. Jaditz et al. (1998) explore the use of nearest neighbor regression models for forecasting industrial production but observe only marginal improvements in predictive performance. Vishwanathan and Murty (2002) present an iterative algorithm for support vector machines in classification problems. Ng (2014) applies a tree ensemble classifier to a large panel of predictors. Fornaro (2016) combines a Bayesian methodology with a shrinkage prior within the probit framework to predict recessions using extensive sets of predictors. More recently, Holopainen and Sarlin (2017), Bluwstein et al. (2020), and Vrontos et al. (2021) employ various machine learning methods for economic event forecasting. Vrontos et al. (2021) provide empirical evidence supporting the application of machine learning over traditional econometric techniques in the context of recession prediction.

Neural networks, as a subset of machine learning, have garnered significant attention and application in fields like finance, primarily due to their ability to establish flexible mappings between variables and exhibit high pattern recognition capabilities (Zhang et al. (1998)). Neural networks, being nonlinear and nonparametric models, can approximate almost any functional form accurately, as stated by the universal approximation theorem (Hornik et al. (1989)). Consequently, they are valuable modeling tools, particularly when there is limited prior knowledge about the appropriate functional relationships. However, the utilization of neural networks in macroeconomic studies has been relatively limited due to the small sample sizes and low-frequency nature of macroeconomic data. Swanson and White (1997) compare artificial neural networks with linear models in terms of predictive performance for nine macroeconomic variables, revealing only marginal improvements in forecast accuracy. Moshiri and Cameron (2000) apply neural networks to inflation forecasting using a dataset of 300 observations spanning 25 years of monthly data. Tkacz (2001) compares mul
tivariate neural networks with linear and univariate models, finding minor forecast improvements in the short term but more pronounced benefits for longer horizons, such as a one-year forecast. Qi (2001) employs a simple feed-forward neural network to predict US recessions using a range of financial and economic indicators, identifying some indicators as useful for prediction. More recently, Puglia and Tucker (2021) compare neural networks with probit regression in forecasting US recessions using the term spread and other macro-financial variables, finding little difference between the models when evaluated. Wang et al. (2022) employ a specific type of recurrent neural network, namely a Bi-LSTM with autoencoder, along with other machine learning models to predict the beginning and end of economic recessions in the US. Their results suggest that the Bi-LSTM with autoencoder is the most accurate model.

The novelty of this paper is twofold: Firstly, it focuses on two special types of recurrent neural network models, LSTM and GRU, which address the limitations of standard recurrent neural networks related to the exploding and vanishing gradient problems. The performance of these models is compared to simple feed-forward neural networks that suffer from the key limitation of having to specify the temporal dependence upfront in the design of the model and traditional linear models in the context of recession forecasting. Secondly, the paper applies the SHAP (SHapley Additive exPlanations) method to GRU to explore the importance of features for different forecast horizons. In the context of recession forecasting, Puglia and Tucker (2021) and Delgado et al. (2022) also use the SHAP method to decompose recession forecasts, but they are applied to other models than LSTM and GRU. The three main findings can be summarized as follows: Firstly, the out-of-sample performance strongly supports the application of LSTM and GRU in the area of recession forecasting, especially for the long-term forecasting tasks. They outperform other types of models across 5 forecasting horizons with respect to different types of statistical performance metrics. Secondly, GRU and ridge logistic regression models differ in assessing predictor importance, evident in the variable order based on SHAP values. Lastly, while the leading predictors for GRU and ridge logistic regression models slightly differ, key indicators like S&P 500 index, real GDP, and private residential fixed investment consistently emerge for short-term predictions (up to 3 months). For longer-term forecasts (6 months or more), the term spread and producer price index take precedence. These results are corroborated by local interpretable model-agnostic explanations (LIME) and marginal effects, respectively.
The remainder of the paper is as follows. Section 2 explains the data used. Section 3 describes the models and forecast evaluation metrics and outlines the research methodology. Section 4 reports the main results, and Section 5 concludes.

## 2 Data

Prior studies in business cycle forecasting often rely on macroeconomic indicators, which are subject to revisions after their initial estimates. Stark and Croushore (2002) demonstrate that the accuracy of forecasts is influenced by using the most up-to-date data instead of real-time data. Therefore, when comparing forecasts from new models to benchmark forecasts, it is crucial to ensure that the comparisons are based on real-time data. In my research, the focus is primarily on assessing the real-time predictability of the Great Recession and the Covid-19 recession in the United States. This necessitates working with real-time data since the information that is available in hindsight was not accessible prior to the recessions. To evaluate the predictability of the recessions, I utilize the same data that actually was available to real-time forecasters for out-of-sample forecasting.

The dataset employed for prediction consists of 194 real-time vintages of macroeconomic and financial market variables, covering the period from February 1967 to October 2021. The out-of-sample forecasting commences in November 2006, utilizing the real-time vintage of data. Taking into account previous studies like Vrontos et al. (2021) and the availability of real-time data, a set of 25 predictors is chosen. A detailed list with descriptions of these variables can be found in Table 1.

The selected predictors cover a wide range of categories, including output, income, prices, the labor market, the housing market, money and credit, and the financial market. The data frequency varies from daily to quarterly. To facilitate model estimation, higher frequency data is aggregated into monthly data using the mean. Quarterly frequency variables are transformed into monthly equivalents using natural cubic spline interpolation, following the methodology outlined in Vrontos et al. (2021). Specifically, at each month, all available data up to that point are used to calculate the interpolating cubic spline. This spline curve is then utilized to generate monthly frequency data points between the given data points.

The majority of monthly data snapshots for the variables are obtained from ALFRED (Archival
Federal Reserve Economic Data). However, there are a few selected predictors for which the real-time data prior to 2013 is not available in ALFRED. The first set of variables, including real personal income excluding transfer receipts and real manufacturing and trade sales, along with total non-farm payroll employment and the industrial production index, are used by Chauvet and Piger (2008) to identify business cycle dates in real-time. The real-time data for these series, provided by Jeremy Piger on his website, ends in August 2013. However, it can be easily extended beyond 2013 by using the data in ALFRED. The second set of variables consists of real M1 and M2 money stock, for which the earliest available real-time data in ALFRED is from January 2014.

\begin{table}
\begin{tabular}{l l l l l l} \hline Nr. & Predictive variable & Abbreviation & Category & Transformation & Frequency \\ \hline
1 & Average hourly earnings of production and nonsupervisory employees & AHETPI & Income & Percent change & Monthly \\
2 & Average weekly hours of production and nonsupervisory employees & AWHNNOAG & Labor market & percent change & Monthly \\
3 & Moody’s BAA yield & BAA & Money and credit & First-order difference & Monthly \\
4 & Moody’s BAA yield relative to 10-Year treasury yield & BAA10YM & Money and credit & First-order difference & Monthly \\
5 & Real manufacturing and trade industries sales & CMRMTSPL & Output & Log growth rate & Monthly \\
6 & Corporate profits after tax & CP & Income & Log growth rate & Quarterly \\
7 & Real disposable personal income & DSPIC96 & Income & Log growth rate & Monthly \\
8 & Effective federal funds rate & FEDFUNDS & Financial market & First-order difference & Monthly \\
9 & Real gross domestic product & GDPC1 & Output & Log growth rate & Quarterly \\
10 & Privately-owned housing units started & HOUST & Housing market & Log growth rate & Monthly \\
11 & Industrial production index & INDPRO & Output & Log growth rate & Monthly \\
12 & Real M1 money stock & M1REAL & Money and credit & First-order difference & Monthly \\
13 & Real M2 money stock & M2REAL & Money and credit & First-order difference & Monthly \\
14 & Non-farm payroll total & PAYEMS & Labor market & Log growth rate & Monthly \\
15 & Real personal consumption expenditures & PCEC96 & Prices & Log growth rate & Monthly \\
16 & Privately-owned housing units permitted & PERMIT & Housing market & Log growth rate & Monthly \\
17 & Producer price index by all commodities & PPIACO & Prices & Log growth rate & Monthly \\
18 & Private residential fixed investment & PRFI & Housing market & Log growth rate & Quarterly \\
19 & S\&P 500 index & SP500 & Financial market & Log growth rate & Daily \\
20 & 3-month treasury bill rate & TB3MS & Financial market & First-order difference & Monthly \\
21 & Term spread - 5-year treasury yield minus 3-month treasury bill rate & T5Y3MM & Financial market & First-order difference & Monthly \\
22 & Consumer Sentiment - University of Michigan & UMCSENT & Prices & Log growth rate & Monthly \\
23 & Unemployment rate & UNRATE & Labor market & First-order difference & Monthly \\
24 & Producer price index by commodity: final demand: finished goods & WPSFD49207 & Prices & Log growth rate & Monthly \\
25 & Real personal income excluding current transfer receipts & W875RX1 & Income & Log growth rate & Monthly \\ \hline \end{tabular}
\end{table}
Table 1: Overview of predictors
To address the absence of real-time vintages for real M1 and M2 stock before this date, nominal M1 and M2 stock are adjusted for inflation using the consumer price index, which has real-time vintages readily accessible in ALFRED.

To identify recession periods in the United States, I rely on the business cycle expansion and contraction dates determined by the National Bureau of Economic Research (NBER). NBER is widely regarded as the standard reference for US business cycles in the existing literature. In this study, recession months are defined as the period following the peak and continuing until the trough, while all other months are considered periods of economic expansion. The earliest available NBER recession indicator vintage in ALFRED is from September 2014. For the monthly vintages preceding that date, I manually collect and construct them based on the official announcements made by the NBER business cycle dating committee1. However, it's worth noting that one major practical challenge with NBER business cycle dates is that they are often announced with significant publication delays.

Footnote 1: [https://www.nber.org/research/business-cycle-dating/business-cycle-dating-committee-announcements](https://www.nber.org/research/business-cycle-dating/business-cycle-dating-committee-announcements)

Table 2 presents the peak and trough dates of the US business cycle from 1980 to 2021, along with their corresponding announcement dates. The publication lags for the most recent six recessions range from 5 to 21 months, with troughs being identified later than peaks on average. While the NBER business cycle dates remain unchanged once finalized, the presence of publication lags complicates the creation of real-time versions of the NBER recession indicator. To address this, the NBER recession indicator in ALFRED is constructed under the assumption that the previous state remains unchanged until a new turning point is officially announced.

## 3 Econometric methodology

In this section, the focus is on the technical aspects of the models used to predict the two most recent recessions in the United States. The paper specifically explores neural networks, which are complex nonlinear models composed of interconnected nodes arranged in multiple layers. These networks have the ability to approximate any linear or nonlinear continuous functions, as stated by the universal approximation theorem Hornik et al. (1989), given that the network is wide or deep enough. In this paper, three types of neural networks are employed for recession forecasting: feedforward neural network (FNN), long short-term memory (LSTM), and gated recurrent unit
(GRU). Section 3.1 provides a detailed description of these models. Model specifications, including estimation and prediction techniques, are discussed in Section 3.2. Furthermore, Section 3.3 introduces statistical measures that effectively evaluate the performance of the predictions.

### Neural Networks

#### 3.1.1 Feedforward neural network

The feed-forward neural network (FFN) is a widely used and straightforward type of artificial neural network. It operates by transmitting information in a unidirectional manner, where data flows from the input layer to the output layer without feedback loops. Figure 1 illustrates, based on the NN-SVG2, a three-layer FFN designed for binary classification. It comprises an input layer with eight nodes, a hidden layer with four nodes, and an output layer with one node. The top

\begin{table}
\begin{tabular}{l l l l} \hline \hline Date & Type & Duration & Announcement \\ \hline
1980:01 & Peak & 6 & 1980:06(+5) \\
1980:07 & Trough & 12 & 1981:07(+12) \\
1981:07 & Peak & 16 & 1982:01(+6) \\
1982:11 & Trough & 92 & 1983:07(+8) \\
1990:07 & Peak & 8 & 1991:04(+9) \\
1991:03 & Trough & 120 & 1992:12(+21) \\
2001:03 & Peak & 8 & 2001:11(+8) \\
2001:11 & Trough & 73 & 2003:07(+20) \\
2007:12 & Peak & 18 & 2008:12(+12) \\
2009:06 & Trough & 128 & 2010:09(+15) \\
2020:02 & Peak & 2 & 2020:06(+4) \\
2020:04 & Trough & ongoing & 2021:07(+15) \\ \hline \hline \end{tabular} The table reports NBER business cycle dates, including the type of cycle (contraction or expansion), the duration in months, and the time of announcement. The data covers the period from 1980 to 2021, and publication lags in months are indicated in parentheses.

\end{table}
Table 2: US Business Cycle dates
nodes in both the input and hidden layers serve as bias nodes. In this configuration, data from the input layer is passed through the hidden layer, which transforms the data. The values obtained from the hidden layer are then forwarded to the output layer, which translates them into desired outputs based on the problem at hand. In the case of binary classification, the last node in the output layer utilizes the sigmoid activation function, producing a value between 0 and 1. This value represents the probability of an event occurring.

Figure 1: The architecture of a feed-forward neural network
Depending on the number of explanatory variables in the data and of nodes in the hidden layer the unknown underlying function \(f\) for an output node can be written as

\[f(X)=g_{2}\Bigg{[}\alpha_{0}+\sum_{j=1}^{k}\alpha_{j}g_{1}\Bigg{(}\beta_{0j}+\sum _{i=1}^{n}\beta_{ij}x_{i}\Bigg{)}\Bigg{]}+\epsilon, \tag{1}\]

where \(n\) is the number of predictors, \(k\) is the number of units in the hidden layer, \(g_{1}\) is the activation function in the hidden layer, \(g_{2}\) is the activation function in the output layer, \(\beta_{ij}\) and \(\alpha_{j}\) represent the weight parameters from the input to the hidden layer and from the hidden to the output layer, respectively, and \(\epsilon\) is the error term. The weight parameters are estimated by the backpropagation process that seeks to repeatedly update the weights until convergence based on the derivatives of the cost function with respect to input and the hidden layers.

#### 3.1.2 Long short-term memory

FFNs are restricted to one-way signal flow, meaning that there is no feedback mechanism where the output of a layer can influence the same layer. Consequently, FFNs lack the ability to capture temporal dependencies in time series data. Although time series data can be fed into an FFN by incorporating additional input units representing previous time points, the main limitation lies in the fixed dimensionality of inputs and outputs. In other words, the precise length of temporal dependence must be predetermined, which is often unknown in real-world scenarios. This is where recurrent neural networks (RNNs) come into play. RNNs establish connections that form cycles, allowing for feedback loops where data can be fed back into the input before being forwarded again. This feedback loop enables RNNs to maintain an internal state or memory to process sequences of inputs. Theoretically, RNNs can retain all information over time and handle long-term dependencies. However, they face two computational challenges. Firstly, as input sequences grow longer, the backpropagation process relies heavily on the chain rule, which may lead to vanishing gradients. If any gradient approaches zero, all other gradients will diminish exponentially fast due to the multiplicative nature of the chain rule. This phenomenon, known as the vanishing gradient problem, prevents effective learning in the model. Secondly, depending on the length of input sequences, the gradient of the loss function can become excessively large and result in numerical instability, referred to as the exploding gradient problem.

Hochreiter and Schmidhuber (1997) introduced the long short-term memory (LSTM) architecture and a corresponding learning algorithm to address the challenges of error back-flow and long
term dependency in recurrent networks. Figure 2 illustrates a basic LSTM neural network comprising an LSTM layer and an output layer. The input data is represented by a three-dimensional tensor, with dimensions for batch size, number of features, and sequence length. The LSTM layer consists of LSTM cells, with each cell responsible for information retrieval at a specific time point in a time series. These cells contain hidden units comprising special nodes or gates that are pre-determined to process the information. Figure 3 provides a closer look at the hidden units of two commonly used variants of recurrent neural networks (RNNs): LSTM and gated recurrent unit (GRU). For both figure 2 and 3, the graphics of the hidden units are adapted from Chris Colah's

Figure 2: The architecture of a long short-term memory network
blog3. The LSTM unit on the left features three gates (forget, input, and output gates) denoted by red dotted lines, along with a cell state, which is a crucial distinction between RNNs and LSTM networks. The cell state functions like a conveyor belt that extends across the entire sequence, facilitating the retention of information over long periods. This mechanism closely resembles the long-term memory function of the human brain. The three gates regulate the flow of information to and from the cell state, enhancing the overall functionality of the network.

Footnote 3: [https://colah.github.io/posts/2015-08-Understanding-LSTMs/](https://colah.github.io/posts/2015-08-Understanding-LSTMs/)

Upon the encounter of information from the previous time step \(h_{t-1}\) with new information \(x_{t}\), the forget gate evaluates which information from the cell state \(C_{t-1}\) should be discarded. Formally, this can be described as follows:

\[f_{t}=\sigma(W_{f}\cdot[h_{t-1}\;\;x_{t}]+b_{f}) \tag{2}\]

Figure 3: Hidden units of a LSTM and a GRU network with specialized gates
Equation (1) presents a concise representation of the mathematical operations within the forget gate. The terms enclosed in brackets represent the linear operations within the activation function, introducing nonlinearity. The input vector \([h_{t-1}\;x_{t}]\) is multiplied by the weight matrix \(W_{f}\) and combined with the bias vector \(b_{f}\), which are then passed through the activation function, typically a sigmoid function. This ensures that the output values fall within the range of 0 to 1. Later, these values will be multiplied by each corresponding element in the previous cell state \(C_{t-1}\), determining the proportion of old information to be retained in the new cell state.

The subsequent step involves identifying the portion of new information deemed valuable for storage in the cell state. This process occurs in two stages. Firstly, a node utilizing the hyperbolic tangent function proposes a vector of new candidate values for the cell state, denoted as \(\tilde{C}_{t}\):

\[\tilde{C}_{t}=\tanh(W_{C}\cdot[h_{t-1}\;x_{t}]+b_{C}). \tag{3}\]

Secondly, the input gate regulates the magnitude of the update, determining which values of \(\tilde{C}_{t}\) will be stored in the cell state \(C_{t}\):

\[i_{t}=\sigma(W_{i}\cdot[h_{t-1}\;x_{t}]+b_{i}). \tag{4}\]

Similarly, the input gate generates values ranging from 0 to 1, which, when multiplied by \(\tilde{C}_{t}\), determine the proportion of the new candidate that should be incorporated into the cell state:

\[C_{t}=f_{t}\times C_{t-1}+i_{t}\times\tilde{C}_{t}/ \tag{5}\]

Equation (5) provides a summary of the update process. The new cell state \(C_{t}\) is obtained by combining the previous state \(C_{t-1}\) multiplied by the forget gate \(f_{t}\), which discards irrelevant information from the old state, and the new candidate values \(\tilde{C}_{t}\) scaled by the input gate \(i_{t}\), retaining only the most important information from the new candidate state.

The resulting new hidden state \(h_{t}\) is determined by the current cell state \(C_{t}\) and the output gate \(o_{t}\):

\[o_{t}=\sigma(W_{o}\cdot[h_{t-1}\;x_{t}]+b_{o}). \tag{6}\]
The output gate in equation (6) utilizes a sigmoid activation function to determine the portion of the new cell state that will be outputted:

\[h_{t}=o_{t}\times\tanh(C_{t}). \tag{7}\]

The cell state undergoes the hyperbolic tangent function and is subsequently multiplied by the output of the output gate, yielding the new hidden state \(h_{t}\). Therefore, \(h_{t}\) selectively captures the most pertinent information from the cell state, representing the short-term memory, while the cell state retains the long-term memory.

#### 3.1.3 Gated recurrent unit

In Figure 3, the right panel shows a gated recurrent unit (GRU) which differs from the LSTM in terms of its gate structure and quantity. GRU, introduced by Cho et al. (2014), aims to simplify the LSTM by combining two gates into a single gate and merging the cell state with the hidden state. Unlike LSTM, GRU eliminates the output gate and integrates the functions of the forget and input gates into a single gate known as the update gate. This update gate processes new information based on the previous hidden state, consequently updating the hidden state.

Initially, the previous hidden state \(h_{t-1}\) and the current input \(x_{t}\) are introduced and passed through the reset gate, which determines the extent to which past information should be disregarded:

\[r_{t}=\sigma(W_{r}\cdot[h_{t-1}\;\;x_{t}]+b_{r}). \tag{8}\]

As the sigmoid activation function outputs values between 0 and 1, the resulting values can be interpreted as the proportion of past information to retain. This value, denoted as \(r_{t}\), is multiplied by \(h_{t-1}\) and passed through the hyperbolic tangent (tanh) activation function. This process yields a vector of new candidate values, represented as \(\tilde{h}_{t}\):

\[\tilde{h}_{t}=\tanh(W_{h}\cdot[r_{t}\cdot h_{t-1}\;\;x_{t}]+b_{h}). \tag{9}\]

The update gate governs the decision of which information to update in the next step:
\[z_{t}=\sigma(W_{z}\cdot[h_{t-1}\ x_{t}]+b_{z}). \tag{10}\]

Finally, the new hidden state \(h_{t}\) is then computed as a weighted average of the previous hidden state \(h_{t-1}\) and the new candidate values \(\tilde{h}_{t}\):

\[h_{t}=(1-z_{t})\cdot h_{t-1}+z_{t}\cdot\tilde{h}_{t}. \tag{11}\]

### Model specifications

Before estimating the model, it is necessary to prepare the chosen type of neural network. Figure 4 illustrates a simplified workflow chart outlining the process of forecasting using neural networks. Recession forecasting using neural networks involves five steps: Data must first undergo preprocessing before being fed into the model. Depending on the specific type of neural networks the network architecture is chosen. Hyperparameters are then fine-tuned using the blocked form of cross-validation. The second and third steps can be combined, treating the number of layers and units within those layers as hyperparameters. Once the data is prepared and the optimal set of hyperparameters is determined, the model is ready to be estimated and generate predictions.

### Data preprocessing

The first step involves preprocessing the data before it can be utilized by the model. This preprocessing stage comprises multiple smaller steps. Initially, each series is evaluated and transformed individually to ensure stationarity. Zhang et al. (2020) argue that while advanced optimization algorithms like RMSProp and Adam are effective for non-stationary data, they incorporate historical gradients in the update calculations, which can result in a lack of relevance to past information when dealing with non-stationary data where distributions change over time. Therefore, it is

Figure 4: Workflow chart of the forecasting process
prudent to use stationary time series data. The second step involves seasonally adjusting and standardizing all explanatory variables. Lastly, a generator function is employed to generate data batches throughout the estimation process. Optimization algorithms based on stochastic gradient descent often utilize smaller batch sizes or minibatches. Bottou (2010) and Ge et al. (2015) highlight that minibatches offer the advantage of avoiding memory loss and introducing sufficient noise into each gradient update, aiding in escaping saddle points or local minima while achieving faster convergence. Additionally, the LSTM and GRU networks require a three-dimensional tensor as input, consisting of the batch size, number of features, and sequence length as axes. To accommodate LSTM and GRU networks, the generator function transforms the data into a suitable three-dimensional tensor format.

#### Network architecture

In Section 3.1, three neural network models are presented, each based on a different architecture. The number of nodes or units in the input layer is determined by the number of features. Depending on the desired number of time points to consider, the input dimension increases accordingly. For example, if the last 12 months of data from 25 predictors are to be used for future predictions, the input dimension would increase from 25 to \(25\times 12=300\). Given the relatively small dataset with a maximum of 657 observations across 25 explanatory variables, I allow a maximum of two hidden layers. The number of hidden layers and the number of nodes in those layers are treated as hyperparameters that should be optimized during the cross-validation process, which will be discussed in more detail later. Both the input and hidden layers use the rectified linear unit as the activation function for their units. The output layer consists of a single unit that uses the sigmoid activation function, producing a number between 0 and 1. This number can be interpreted as the probability of a recession.

The LSTM neural network architecture consists of as many layers as the number of time steps to look back, along with an output layer. In the context of recession forecasting, I choose to examine the temporal variations of the predictors over the last 12 months, or 1 year, which corresponds to 12 lagged variables of the same predictor. For each time point from \(t-1\) to \(t-12\), there are 12 LSTM layers, each with a same number of hidden units ranging from 16 to 64. The first unit in the first LSTM layer is linked to the first unit in the second LSTM layer, and the same holds
for the rest of the LSTM layers. This approach ensures that each chain of units explores different feature dimensions at different time points. Additionally, the network may require a second chain of layers, leading to a stacked LSTM. The optimal number of stacked layers is determined by the hyperparameter optimization process. The output layer in LSTM is the same as in the FFN. The GRU follows the same architectural structure as LSTM but differs in the structure and functionality of the hidden units within the GRU layer.

#### Hyperparameter optimization

Deep neural networks with multiple hidden layers and numerous parameters have the capacity to learn complex relationships. However, when training data is limited, some of these relationships may be the outcome of sampling noise and might not be present in real test data. This gives rise to overfitting, which transpires when a model fails to generalize from observed data to new data. Overfitting is a commonly acknowledged challenge in the supervised machine learning framework, and various methods employing different strategies have been devised to mitigate its effects and prevent the model from overfitting. Ying (2019) categorizes methods to address overfitting into four groups: network reduction, data augmentation, early stopping, and regularization.

Network reduction involves reducing the depth and width of a neural network by decreasing the number of layers and units, thereby reducing the total number of parameters to estimate. This approach helps the model focus on capturing essential patterns in the training data, leading to a better generalization of test data. To implement network reduction, I limit the network to a maximum of two hidden layers, and choose a relatively small number of units in these layers. Additionally, through cross-validation, the model selects the optimal depth and width that minimizes validation loss, ensuring good generalization on unseen validation data.

Complex models with numerous parameters require a substantial amount of data to distinguish meaningful patterns from noise. Therefore, expanding the training data is an effective approach to enhance the model's generalizability. Data augmentation involves not only acquiring more training data but also employing techniques directly applied to the existing data, such as random sampling or reshuffling. I am unable to address this issue due to two reasons. Firstly, I use the same data as in Chung (2022) to compare the performance of neural networks with traditional statistical models like logistic regression. Therefore, data augmentation is not suitable as it would alter the original
dataset. Secondly, the temporal dependence of the time series data across multiple variables is a crucial aspect that needs to be explored. In this context, data augmentation strategies such as random sampling or reshuffling are less desirable.

Instead, early stopping is implemented as a strategy to halt the training process when the gap between training loss and validation loss begins to widen. The weights are then set to the values corresponding to the smallest gap, preventing the model from excessively memorizing the training data. This helps in improving the model's performance on unseen data. To incorporate early stopping, I include a call during the estimation process that stops the iteration and restores the best weights when the validation loss does not decrease for 5 consecutive epochs during cross-validation and for 10 consecutive epochs during the final estimation, where an epoch means one complete pass of the training dataset through the algorithm.

Regularization is a technique used to mitigate the impact of less important features in a model. There are two commonly applied categories of regularization methods in neural networks. The first category involves adding a penalty term or regularizer to the loss function. Two well-known types of regularizers, namely \(L1\) and \(L2\), are commonly used, similar to their application in penalized regression. \(L1\) regularization, known as Lasso regression (Tibshirani (1996)), assigns zero weights to unimportant features, effectively removing them from the model. This ensures that only influential features with significant effects on the variable of interest are retained. On the other hand, \(L2\) regularization, known as Ridge regression (Hoerl and Kennard (1970)), assigns lower weights to unimportant features rather than discarding them completely. This approach aims to extract as much relevant information as possible while controlling model complexity and reducing overfitting. In my model, I incorporate \(L2\) regularization to address these concerns, as Chung (2022) demonstrated the superiority of Ridge regression over Lasso regression in the context of recession forecasting. The second category of regularization methods involves using a technique called dropout. Dropout was introduced by Srivastava et al. (2014) and involves randomly dropping units and their connections from the neural network during training. A certain percentage of hidden units are randomly dropped to create a thinned network, which is then trained using stochastic gradient descent. After training, the dropped units are restored, and the process is repeated. Dropout can be seen as approximating the effects of averaging the estimates over multiple smaller networks, while also preventing overfitting. I leverage the dropout technique in my model,
treating the dropout percentage as a hyperparameter to be optimized through cross-validation.

Neural networks have hyperparameters that need to be determined prior to training, in addition to the parameters that are estimated during training. These hyperparameters include the number of layers, units in each layer, activation function, batch size, learning rate, dropout percentage, and more. The values of these hyperparameters must be specified before training can begin. There are various methods to search for the best combination of hyperparameters, known as hyperparameter optimization. One simple approach is manual search, where the researcher selects their own set of hyperparameters based on existing literature and evaluates the model's performance on a validation dataset. This process is repeated with different hyperparameter settings until the best combination is found. However, manual search can be time-consuming and does not guarantee finding the optimal hyperparameters. Grid search is another approach, where a predetermined set of values is specified for each hyperparameter. Every possible combination of values is then evaluated, resulting in a large number of trials. However, as indicated by Bellman (1961), the number of possible combinations grows exponentially with the number of hyperparameters, leading to the curse of dimensionality. This can make grid search impractical for models with many hyperparameters. A more effective approach, demonstrated by Bergstra and Bengio (2012), is random search. Random search involves randomly selecting values from predetermined ranges for each hyperparameter. The advantage of random search is that it can identify good or even better models within a smaller fraction of computation time compared to manual or grid search methods. Following this approach, I perform random search using predefined ranges of values for the hyperparameters. The specific ranges used for random search in the optimization process are reported in Table 3.

Figure 5 illustrates the form of cross-validation, which I employ to determine the optimal set of hyperparameters. Each block of the training set is divided into two sections at every iteration, with the validation set always following the training split. This approach ensures that the natural order of observations is maintained within each block. For each block, the size of the validation set remains the same, while the training set gets larger as the validation set of the prior block is added to it. The number of iterations in the cross-validation setup depends on the size of each block. If the block size is small, the number of blocks and iterations increases. However, selecting an excessively small block size may result in certain blocks having insufficient or no recession data, leading to sampling bias. Considering the information presented in Table 1, which indicates the
longest period between two recessions in US history prior to the Great Recession as 128 months, I set the length of each validation block to be 128 months to ensure that each validation set contains data from at least one recession. I let the length \(l\) of the first training block be equal to \(l=N_{total}-128\times(\lfloor\frac{n_{total}}{128}\rfloor-1)\), where \(N_{total}\) denotes the total number of observations in a data vintage. This guarantees that the first training block is always larger than the validation block. Subsequently, I employ an expanding window approach, increasing the length of the training set by 128 months at each step while maintaining a validation block of the same length (128 months). This methodology ensures that scarce recession data are fully utilized.

For each combination of hyperparameters, I calculate the average validation loss across all cross-validations and store it. The set of hyperparameters with the lowest average validation loss is selected as the optimal choice. According to mathematical calculations by Bergstra and Bengio (2012), using random search with 60 trials is likely to yield a set of hyperparameters that falls within the top 5% interval around the optimal solution with a 95% probability. Therefore, I conduct 60 trials each time I train the model for new predictions. Once the optimal combination of hyperparameters is determined, I split the entire dataset into training and validation parts in the same proportion, maintaining the chronological order of the data. I then perform a final estimation and validation of the model to improve the accuracy of the weights and check for any

\begin{table}
\begin{tabular}{l|l} \hline Type & Range \\ \hline \# hidden layer & \(\in(1,2)\) \\ \hline \# unit & \(\in(16,32,64)\) \\ \hline batch size & \(\in(16,32,64)\) \\ \hline dropout & \(\in(0,0.1,0.2,0.3,0.4,0.5)\) \\ \hline recurrent dropout & \(\in(0,0.1,0.2,0.3,0.4,0.5)\) \\ \hline weight decay & \(\in(0,0.1,0.2)\) \\ \hline learning rate & \(\in(0.01,0.001)\) \\ \hline \end{tabular} The table provides the ranges of possible values for different hyperparameters of the neural networks used in the analysis. The optimal combination of these hyperparameters is determined using the time series cross-validation technique during the hyperparameter optimization process.

\end{table}
Table 3: Tuning hyperparameters
signs of overfitting.

#### Model estimation

The initial dataset used to train the model consists of real-time data from October 2005. It contains the latest available data in October 2005 which stretches from February 1967 to August 2005. For each month in the forecasting period, a dataset is constructed containing the NBER recession indicator and other relevant time series data. This dataset is used to train the model and make predictions for the following month. The forecasting process assumes that predictions are made on the first day of each month. The forecaster incorporates all available data up to the end of the previous month to predict for the current month and beyond, based on the forecast horizon. Most of the data for a specific month is typically published by the last day of the following month. As a result, there is a lag of at least two months (depending on the forecast horizon) between the

Figure 5: Cross-validation for time series
latest available data and the month being predicted. Occasionally, data may have time delays or missing values. In such cases, a \(k\)-nearest-neighbor imputation method is used during the data pre-processing stage. Experimenting with alternative methods, such as the bagged tree imputation algorithm, does not affect the results significantly, although the latter is known for its superior performance on account of higher computational costs according to existing literature (Stekhoven and Buhlmann (2011)).

Once the hyperparameters are optimized, the neural network is prepared to fit the training data. The data is fed into the model, processed, in the manner of equation (1) in case of a FFN, and transformed into an output. This output is then compared to the actual output, and the error signals are propagated backward through the network to adjust the weights accordingly. This iterative process is known as backpropagation and continues until a specific condition is satisfied. Although the concept appears straightforward, the estimation process involves extensive mathematical computations. The parameter values in equation (1) are estimated to minimize the following binary cross-entropy loss function:

\[\text{Loss}=\frac{1}{n}\sum_{i=1}^{N}-[y_{i}\cdot\ln(p_{i})+(1-y_{i})\cdot\ln( 1-p_{i})]. \tag{12}\]

In most cases, there is no analytical solution available for minimizing such problems, and therefore the parameters need to be estimated numerically. To perform this estimation, I utilize the Adam optimizer, which was introduced by Kingma and Ba (2014). The Adam optimizer is specifically designed for optimizing stochastic loss functions using first-order gradients. It combines the strengths of two popular optimization methods, AdaGrad (Duchi et al. (2011)) and RMSProp (Tieleman and Hinton (2012)). The Adam optimizer offers computational efficiency and requires minimal memory resources.

#### Prediction

To fully utilize the capacity of neural networks, I incorporate all 12 monthly lags of each predictor over the past year. De Veaux and Ungar (1994) argue that due to the overparameterization of neural networks, individual weights associated with multicollinearity become less influential. The final output of a neural network is a result of various combinations of activation functions that
involve interactions among predictors, making the impact of multicollinearity typically insignificant. Moreover, the backpropagation algorithm used in neural networks does not require inverting matrices, which can be problematic when there is perfect or severe multicollinearity.

Since the effectiveness of predictive variables may vary across different forecast horizons, I explore five distinct windows: nowcasting, immediate-term, short-term, medium-term, and long-term. These windows correspond to predicting the current month, one month ahead, three months ahead, six months ahead, and twelve months ahead, respectively. Accounting for the two-month publication lag, the forecast horizons represent 2-, 3-, 5-, 8-, and 14-steps-ahead forecasts using the latest set of predictive variables. For example, suppose it is November 2006, and I have data available up to September 2006, which was published in October 2006. If I aim to conduct nowcasting and forecast the recession for the current month, I use all available data up to September 2006 to predict for November 2006, resulting in a two-steps-ahead forecast. Similarly, for other forecast horizons, two extra months or steps need to be considered.

After the models generate probability predictions, they are transformed into monthly zero-one indicators that represent the state of the economy, such as recession or boom, based on predefined cutpoints. The conventional approach is to use a fixed threshold, typically set at 0.5. Under this approach, if a probability prediction is equal to or greater than 0.5, it is classified as a recession, and vice versa. However, it is often mentioned that business cycles exhibit asymmetry, and therefore a 50% threshold may not be optimal for linear modeling methods. Berge and Jorda (2011) propose an optimal cutpoint within the range of 0.3 to 0.6 based on smoothed state probabilities estimated by Chauvet and Piger (2008). Ng (2014) suggests using different thresholds for different forecast horizons, ranging from 0.3 to 0.44. Vrontos et al. (2021) adopt a fixed threshold of 0.33 for classification and evaluation purposes. In the case of neural networks, they are capable of capturing the asymmetry and other nonlinear characteristics of business cycles. Therefore, I opt to follow the traditional approach and use a fixed threshold of 0.5.

The out-of-sample forecasts are produced in a quasi-recursive fashion to capture dynamic structural changes in the data. To manage computational resources, a forecast window of 12 months is selected, during which the estimated parameters remain constant. This means that once a neural network model is built, it is used to predict the next 12 months using the latest available data. Afterward, the observed data for those 12 months is added to the existing sample, and the neu
ral network model is reestimated using this augmented dataset. The updated parameters are then used to forecast the following 12 months, and so on. Although this recursive reestimation approach is time-consuming and costly, it mimics the real-time process of generating predictions using the most up-to-date data available. To clarify this procedure, I focus on the nowcasting forecasts. The out-of-sample period begins in November 2006, and we want to predict for that month. The available data in November 2006 goes up until September 2006. The training and validation set consists of 476 monthly observations from February 1967 to September 2006. Using the optimized hyperparameters obtained through blocked cross-validation, I compute the nowcasting forecast for November 2006. The predicted value for this month is then compared to the actual value from the January 2007 vintage. The same model uses the most recent data available each month to generate predictions for the next 12 months. For the next iteration, data of the last 12 months are added to the training and validation set. The model is reestimated and used to predict for November 2007 and onwards. This expanding window procedure is repeated, gradually increasing the size of the training and validation set, to generate out-of-sample forecasts for the testing period spanning from November 2006 to October 2021. The same approach is adopted for the other forecasting horizons.

### Performance evaluation

The estimated models produce two types of forecasts. Firstly, the models generate probability predictions. In the second step, these probabilities are transformed into binary point predictions using a specified threshold. While the second type may not always be necessary, I include it to ensure the forecast is more easily interpretable for the end-users. The performance of both the probability predictions and point predictions is evaluated using various statistical measures. With the exception of metrics specifically related to probability predictions, the evaluation is based on the output of a contingency table, commonly known as a confusion matrix as illustrated in Table 4. This matrix provides a framework for assessing the accuracy of the predictions.

The confusion matrix is composed of elements that represent the count of observations belonging to different categories. True positives (TP) and true negatives (TN) indicate the correct classification of positive and negative outcomes, respectively. False positives (FP) occur when the prediction is positive while the actual value is negative. Conversely, false negatives (FN) arise when the prediction is negative while the actual value is positive. By utilizing these counts, various
performance evaluation measures are derived to assess the accuracy of the predictions.

The performance evaluation of the forecasting model is assessed using the metrics listed in Table 5. The ROC (Receiver Operating Characteristics) curve displays the entire set of possible combinations of true positive rates \(TPR(c)=\frac{TP(c)}{\#\:Actual\:Positive}\) and false positive rates \(FPR(c)=\frac{FP(c)}{\#\:Actual\:Negative}\) for some cutpoint \(c\in(0,1)\) that maps the predicted probability to a binary cate

\begin{table}
\begin{tabular}{|c|c|c|c|} \hline \multicolumn{2}{|c|}{} & \multicolumn{2}{c|}{Predicted} \\ \cline{3-4} \multicolumn{2}{|c|}{} & \multicolumn{2}{c|}{Positive} & \multicolumn{1}{c|}{Negative} \\ \hline \multirow{3}{*}{Actual} & Positive & TP & FN \\  & Negative & FP & TN \\ \hline \end{tabular} The table exhibits an example of the confusion matrix containing elements representing the count of observations belonging to each category. True positives (TP) and true negatives (TN) indicate the correct classification of positive and negative outcomes, respectively. False positives (FP) occur when the prediction is positive, but the actual value is negative. Conversely, false negatives (FN) arise when the prediction is negative, but the actual value is positive.

\end{table}
Table 4: Confusion Matrix

\begin{table}
\begin{tabular}{|c|l|c|} \hline \hline Type & Metric & Formula \\ \hline \multirow{2}{*}{Probability} & Area under the ROC curve & \(\int_{0}^{1}ROC(c)\,dc\) \\ \cline{2-3}  & Area under the PR curve & \(\int_{0}^{1}PR(c)\,dc\) \\ \hline \hline \multirow{5}{*}{Pediction} & Sensitivity & \(\frac{TP}{TP+FN}\) \\ \cline{2-3}  & Specificity & \(\frac{TN}{FP+TN}\) \\ \cline{2-3}  & Precision & \(\frac{TP}{TP+FP}\) \\ \cline{2-3}  & Balanced accuracy & \(\frac{Sensitivity+Specificity}{2}\) \\ \cline{2-3}  & Matthews correlation coefficient & \(\frac{TP\times TN-FP\times FN}{\sqrt{(TP+FP)(TP+FN)(TN+FP)(TN+FN)}}\) \\ \cline{2-3}  & \(F_{1}\)-Score & \(\frac{2}{Sensitivity+\frac{1}{Precision}}\) \\ \hline \end{tabular} The table reports the metrics used for the performance evaluation of the forecasting model. They are divided into two groups, depending on which type of predictions they refer to.

\end{table}
Table 5: Performance evaluation metrics
gory. Similarly, the PR (Precision-Recall) curve plots the complete set of possible combinations of precision and recall for \(c\in(0,1)\), where recall is a synonym for sensitivity. Both the area under the ROC curve and the PR curve increase with the underlying metrics for a given cutpoint. By aggregating over the entire set of cutpoints, these curves assess the overall predictive ability of the forecasting model, regardless of the specific cutpoint used. Tharwat (2021) provides a comprehensive overview of these metrics, while Chung (2022) briefly discusses their strengths and weaknesses, particularly in the context of imbalanced binary classification problems.

## 4 Empirical results

Table 6 presents the out-of-sample forecast performance of the models for various forecast periods. A carefully selected set of statistical metrics is used to assess the performance from different perspectives. The metrics AUROC (Area Under the Receiver Operating Characteristic Curve) and AUPRC (Area Under the Precision-Recall Curve) are employed to measure the models' raw performance in the sense that these metrics do not require any threshold to convert probabilities into binary outcomes, making them suitable for evaluating the models' pure performance. Another set of metrics is derived from a contingency table, which necessitates the use of thresholds to convert predicted probabilities into binary outcomes. The performance values of the logistic and ridge regression models, obtained from Chung (2022), are included in the analysis and can be directly compared to the neural network models in this study since they are based on the same data and forecast horizons.

In the nowcasting setup, there is little difference in predictive performance between the types of logistic regression models and of neural network models. Logistic regression models are not outperforming other models in any considered metric, but the difference is smaller compared to other forecast horizons. The ridge regression model performs similarly to the feed-forward neural network (FFN). Although the FFN performs slightly better, the differences are less than 0.1 for most of the metrics. The LSTM and GRU models show significantly better performance in terms of MCC and \(F_{1}\)-Score, mainly due to their high precision values. For instance, the LSTM correctly identifies 90% of recession months, while 75% of its positive predictions are accurate. In contrast, the ridge regression model has a precision of 51.5%, indicating that the LSTM accurately predicts recessions without generating too many false alarms. The AUPRC values also highlight
\begin{table}
\begin{tabular}{l c c c c c c c c} \hline Method & AUROC & AUPRC & BAcc & MCC & \(F_{1}\)-Score & Sensitivity & Specificity & Precision \\ \hline \multicolumn{10}{l}{Panel A: nowcasting setup} \\ \hline Logit & 0.853 & 0.365 & 0.815 & 0.506 & 0.556 & 0.750 & 0.881 & 0.441 \\ Ridge & 0.920 & 0.529 & 0.875 & 0.609 & 0.642 & 0.850 & 0.900 & 0.515 \\ FFN & 0.917 & 0.642 & 0.906 & 0.668 & 0.692 & 0.900 & 0.912 & 0.563 \\ LSTM & 0.899 & 0.754 & 0.931 & 0.797 & 0.818 & 0.900 & 0.962 & 0.750 \\ GRU & 0.890 & 0.837 & 0.928 & 0.778 & 0.800 & 0.900 & 0.956 & 0.720 \\ \hline \multicolumn{10}{l}{Panel B: immediate-term setup} \\ \hline Logit & 0.671 & 0.243 & 0.668 & 0.265 & 0.357 & 0.500 & 0.837 & 0.278 \\ Ridge & 0.931 & 0.555 & 0.894 & 0.693 & 0.723 & 0.850 & 0.937 & 0.630 \\ FFN & 0.860 & 0.652 & 0.825 & 0.623 & 0.667 & 0.700 & 0.950 & 0.636 \\ LSTM & 0.847 & 0.514 & 0.890 & 0.607 & 0.632 & 0.900 & 0.881 & 0.487 \\ GRU & 0.896 & 0.831 & 0.944 & 0.887 & 0.900 & 0.900 & 0.987 & 0.900 \\ \hline \multicolumn{10}{l}{Panel C: short-term setup} \\ \hline Logit & 0.524 & 0.170 & 0.591 & 0.165 & 0.267 & 0.300 & 0.881 & 0.240 \\ Ridge & 0.928 & 0.522 & 0.856 & 0.601 & 0.640 & 0.800 & 0.912 & 0.533 \\ FFN & 0.872 & 0.800 & 0.841 & 0.732 & 0.757 & 0.700 & 0.981 & 0.824 \\ LSTM & 0.858 & 0.395 & 0.784 & 0.508 & 0.565 & 0.650 & 0.918 & 0.500 \\ GRU & 0.885 & 0.544 & 0.931 & 0.797 & 0.818 & 0.900 & 0.962 & 0.750 \\ \hline \multicolumn{10}{l}{Panel D: medium-term setup} \\ \hline Logit & 0.632 & 0.174 & 0.565 & 0.097 & 0.226 & 0.350 & 0.780 & 0.167 \\ Ridge & 0.910 & 0.423 & 0.818 & 0.517 & 0.566 & 0.750 & 0.887 & 0.455 \\ FFN & 0.795 & 0.496 & 0.731 & 0.510 & 0.556 & 0.500 & 0.962 & 0.625 \\ LSTM & 0.853 & 0.519 & 0.903 & 0.655 & 0.679 & 0.900 & 0.906 & 0.546 \\ GRU & 0.731 & 0.250 & 0.647 & 0.266 & 0.356 & 0.400 & 0.893 & 0.320 \\ \hline \multicolumn{10}{l}{Panel E: long-term setup} \\ \hline Logit & 0.737 & 0.250 & 0.668 & 0.297 & 0.383 & 0.450 & 0.887 & 0.333 \\ Ridge & 0.818 & 0.388 & 0.712 & 0.353 & 0.431 & 0.550 & 0.874 & 0.355 \\ FFN & 0.873 & 0.391 & 0.640 & 0.294 & 0.368 & 0.350 & 0.931 & 0.389 \\ LSTM & 0.862 & 0.371 & 0.725 & 0.396 & 0.468 & 0.550 & 0.899 & 0.407 \\ GRU & 0.909 & 0.674 & 0.837 & 0.707 & 0.737 & 0.700 & 0.975 & 0.778 \\ \hline \end{tabular} The table reports the performance evaluation measures of forecasts obtained by logistic regression models, ridge regression models, feed-forward (FFN), long short-term memory (LSTM), and gated recurrent unit (GRU) neural networks over different time horizons for the out-of-sample period, November 2006 to October 2021: Panel (A) presents the nowcasts. Panel (B), (C), (D), and (E) display the 1-months-ahead forecasts, the 3-months-ahead forecasts, the 6-months-ahead forecasts, and the 12-month-ahead forecasts, respectively.

\end{table}
Table 6: Performance evaluation measures: A real-time assessment
the superiority of LSTM and GRU models over the others, with GRU achieving an AUPRC of 0.837 compared to 0.642 for FFN and 0.529 for ridge regression. Davis and Goadrich (2006) argue that precision-recall curves give a more informative picture of an algorithm's performance, when dealing with highly imbalanced dataset. In fact, the ratios of the number of months in boom to the number of months in recession in 194 real-time datasets range from 4.5 to 6.7 with a mean value of 6, which means that the period of booms is in average 6 times longer than that of recessions. The percentage of recessions in the datasets fluctuates around 0.2 across the vintages. Although not extremely skewed, the datasets are certainly unbalanced, and the difference in AUPRC values suggests that LSTM and GRU models handle recessions better in the presence of class imbalance.

Moving to the immediate-term setup, logistic regression models perform significantly worse, with MCC and \(F_{1}\)-Score approximately 40% lower than before. Conversely, the other models show similar or slightly improved performance. Among them, the GRU exhibits the highest performance across various summarizing metrics such as balanced accuracy, MCC, and \(F_{1}\)-Score.

In the short-term setup, although the ridge regression model demonstrates better predictive performance than LSTM, on average, neural network models, particularly GRU, outperform the other models. GRU maintains a 90% accuracy in classifying recession months with higher specificity and precision, resulting in higher values of balanced accuracy, MCC, and \(F_{1}\)-Score.

In the medium-term setup, LSTM performs slightly better than the other models, while GRU performs poorly, and thus the difference between the two groups is minimal. However, LSTM achieves 90% sensitivity and over 90% specificity, while maintaining precision above 0.5. The overall average performance of the models is the lowest among different forecast horizons.

In the long-term setup, GRU surpasses the other models by a significant margin, with an \(F_{1}\)-Score of 0.737 compared to 0.423 for the ridge regression model. LSTM accurately predicts 70% of recession months, with 77.8% of its positive predictions being correct.

Neural network models, especially LSTM and GRU, show a significant improvement in predictive performance compared to logistic regression models. Across all metrics considered, neural network models consistently outperform logistic regression models by a significant margin, except for the nowcasting setup. This superiority is evident in both pure performance metrics such as AUROC and AUPRC, as well as threshold-based point prediction performances. The disparity is particularly prominent in AUPRC, which is a more reliable indicator when dealing with imbal
anced data, highlighting the neural network models' ability to handle class imbalance effectively. Despite logistic regression models employing more advanced techniques to adjust thresholds in the presence of class imbalance while neural network models use a fixed threshold of 0.5, the neural network models consistently outperform logistic regression models. The performance difference diminishes when comparing them to ridge regression models. Ridge regression models exhibit higher AUROC values compared to other model specifications, except for the long-term setup. However, in terms of AUPRC, neural networks again prove to be a better choice for capturing class imbalance. On average, neural network models perform better than ridge regression models in terms of these metrics. The gap widens further when considering only the two recurrent neural network models. Depending on the forecast horizons, either LSTM or GRU may lead the competition, with occasional substantial differences between them and the other models.

The ROC and PR curves of the models for five different forecast horizons are presented in Figures 6 and 7, respectively. Panel B and C of the ROC curves highlight the noticeable difference in forecast performance between the standard logistic regression model (represented by the black line) and the other models. In Panel D and E, although the lines are closer together, there is still some discernible gap between each of the models. Overall, the ridge regression, LSTM, and GRU models exhibit superior performance across various forecasting horizons. When dealing with class imbalance, it is beneficial to complement ROC curves with PR curves since the latter focuses on precision. Precision refers to the fraction of true positives among all predicted positives. PR curves are specifically designed for detecting rare events and determine if the classifier can accurately classify most instances of the minority class. Figure 7 demonstrates that, on average, the green, blue, and purple lines lie above the red line and, even more so, above the black line. This finding supports the argument that neural network models, compared to linear models, are more capable of handling class imbalance in the data.

Figure 8 displays the out-of-sample forecasts in probabilities from November 2006 to September 2021, accompanied by shaded areas representing the Great Recession and the Covid-19 recession. Across all panels of Figure 8, similar patterns emerge: the graphs start at high levels before the Great Recession, either remain high or increase during that period, decrease afterward to a range between 0 and 0.25 for an extended period, and then exhibit significant fluctuations around the Covid-19 recession. The models successfully indicate economic and financial downturns during the
The figure illustrates receiver operating characteristic (ROC) curves for the out-of-sample period from November 2006 to October 2021: Panel (A) shows the ROC curves for the nowcasting forecast horizon, comparing the logistic regression model, ridge regression model, FFN, LSTM, and GRU models. Panels (B), (C), (D), and (E) present the ROC curves for the immediate-term (1-months-ahead), short-term (3-months-ahead), medium-term (6-months-ahead), and long-term (12-months-ahead) forecast horizons, respectively.

Figure 6: Receiver Operating Characteristic (ROC) curves
The figure depicts precision recall (PR) curves for the out-of-sample period from November 2006 to October 2021: Panel (A) shows the PR curves for the nowcasting forecast horizon, comparing the logistic regression model, ridge regression model, FFN, LSTM, and GRU models. Panels (B), (C), (D), and (E) present the ROC curves for the immediate-term (1-months-ahead), short-term (3-months-ahead), medium-term (6-months-ahead), and long-term (12-months-ahead) forecast horizons, respectively.

Figure 7: Precision Recall (PR) curves
The figure presents out-of-sample recession probabilities for the period from November 2006 to October 2021:: Panel (A) displays the predicted recession probabilities for the nowcasting forecast horizon, comparing the logistic regression model, ridge regression model, FFN, LSTM, and GRU models. Panels (B), (C), (D), and (E) show the predicted recession probabilities for the immediate-term (1-months-ahead), short-term (3-months-ahead), medium-term (6-months-ahead), and long-term (12-months-ahead) forecast horizons, respectively. The grey shaded areas depict NBER recession months.

Figure 8: Out-of-sample recession probabilities
Great Recession, but an external event like the Covid-19 recession remains largely unpredictable. Depending on the forecast horizon, once the model learns from the impact of the Covid-19 pandemic on the economy, it struggles to interpret this abrupt change and generates alternating forecasts.

## 5 Feature importance

To gain a better understanding of the underlying factors driving the results, it is important to identify the main variables or features that have the most significant impact on the output. While linear models can utilize estimated coefficients or marginal effects to measure feature importance, more complex models like neural networks require different approaches. Molnar (2020) provides a comprehensive list of methods that enable the interpretation of complex machine learning models. One such method is Local Interpretable Model-Agnostic Explanations (LIME), which suggests the use of local models to estimate effects. For an instance of interest, LIME constructs a new dataset by altering samples feature-wise, drawing from a normal distribution with sample mean and sample standard deviation, and obtains corresponding predictions from the black box model. This modified dataset is then used to train an interpretable model, such as Lasso or a decision tree, which is weighted by the proximity of the sampled instances to the instance of interest. The interpretable model serves as an approximation of the black box model's predictions at a local level, even though it may not accurately represent the global approximation.

Another solution stems from cooperative game theory, particularly the Shapley value introduced by Shapley (1953). The Shapley value method involves assigning payouts to players based on their contribution to the overall payout. This cooperative framework resembles a game where players form coalitions and receive profits based on their cooperation. In the context of model predictions and their interpretability, the game pertains to the prediction task, and the payout represents the difference between the actual prediction for that data point and the average prediction across all data points. The players correspond to the feature values of the data point, working together to achieve the payout, which signifies predicting a specific value. To calculate the Shapley value for a specific feature value, all possible coalitions of feature values, excluding the feature of interest, are formed for each data point. The values of features outside a coalition are substituted with random values of those features from data to generate a prediction. The predictions are computed for each coalition, both with and without the feature value of interest. The difference between
these predictions represents the marginal contribution of the feature value for that coalition. This computation is repeated for all possible coalitions, and the average of the marginal contributions across all coalitions yields the Shapley value for that feature value. Finally, the Shapley values for a feature can be averaged across data points to assess the relative importance of features in comparison to each other.

The Shapley value stands out as the sole explanation method supported by a robust theory that satisfies the axioms of efficiency, symmetry, dummy, and additivity to allow for fair distribution of the contributions among the features. Conversely, methods like LIME rely on the assumption of local linear behavior of the machine learning model, but lack a theoretical justification for why this approach is effective. Lundberg and Lee (2017) propose an alternative estimation method, the Shapley additive explanations (SHAP), which uses a kernel-based estimation approach for Shapley values (KernelSHAP). SHAP introduces an innovative aspect by presenting the Shapley value explanation in the form of an additive feature attribution method, which can be viewed as a linear model. This perspective establishes a connection between LIME and Shapley values. In this context, SHAP defines the explanation model as follows:

\[f(c^{\prime})=\alpha_{0}+\sum_{j=1}^{K}\alpha_{j}c^{\prime}_{j}\]

In the provided equation, the explanation model is denoted as \(f\), the coalition vector as \(c^{\prime}\in\{0,1\}^{K}\), the maximum coalition size as \(K\), and the Shapley value for a feature \(j\) represented by \(\alpha_{j}\). The coalition vector \(c\) indicates the presence or absence of feature values, with a value of 1 representing presence and 0 indicating absence.

KernelSHAP is a method that estimates the contributions of individual feature values to the predictions by first generating a random coalition \(c^{\prime}\in\{0,1\}^{K}\) consisting of K members by flipping a coin multiple times until we obtain a sequence of 0's and 1's. The sampled coalition of size \(K\) is then used as a data point for the regression model. In this model, the target is the prediction for a coalition, whereby in case of \(c_{k}=0\), the absent feature value is substituted with random feature values from the data. This process is repeated for each data point. Then Shapley compliant weights are computed according the the SHAP kernel that is proposed by Lundberg and Lee (2017). Finally, we fit a weighted linear regression model on the modified data, and the estimated coefficients of the model are the Shapley values.
The main distinction from LIME lies in how instances are weighted in the regression model. In LIME, the weighting is determined based on the proximity of instances to the original instance. The closer an instance is, the higher its weight in LIME. On the other hand, SHAP assigns weights to sampled instances based on the weight they would receive in the Shapley value estimation. Small coalitions (with fewer 1's) and large coalitions (with many 1's) receive the highest weights in SHAP. The underlying intuition is that we gain the most knowledge about individual features when we can study their effects in isolation. For a coalition consisting of a single feature, we can observe the isolated main effect of that feature on the prediction. When a coalition includes all features except one, we can learn about the total effect of that particular feature, including both the main effect and feature interactions. However, if a coalition comprises half the features, it provides limited insight into the contribution of an individual feature due to the numerous possible coalitions with half of the features.

The SHAP values for logistic regression models are computed by the KernelSHAP method using KernelExplainer from SHAP package in python, whereas the SHAP values for deep learning models are approximated by an enhanced version of the DeepLIFT (Deep Learning Important FeaTures) algorithm introduced by Shrikumar et al. (2017) using DeepExplainer from the same package. DeepLIFT is a method used to analyze the output prediction of a neural network for a specific input. It achieves this by backpropagating the contributions of all neurons in the network to each feature of the input. DeepLIFT compares the activation of each neuron to its reference activation and assigns contribution scores based on the difference between them. Unlike other approaches, DeepLIFT has the ability to uncover dependencies that may be overlooked. It can also consider positive and negative contributions separately. Additionally, the scores can be efficiently computed in a single backward pass.

SHAP method generates a matrix of Shapley values that contains one row per data instance and twelve columns per feature for each time step. These Shapley values can be used to create global explanations by averaging the absolute Shapley values per feature across the twelve time steps and data points. This process is repeated for the entire out-of-sample period. The results for GRUs are presented in Figure 9, which depicts the distribution of average absolute SHAP values across features. The features are ranked based on their medians, providing insights into their importance for different forecasting horizons.
Figure 9: SHAP values of the predictors: GRU
The figure shows the boxplots of the absolute average SHAP values of the predictors in descending order based on ridge models: Panel (A) to (E) display absolute average SHAP values of the predictors related to the nowcasting, immediate-term (1-month-ahead), short-term (3-months-ahead), medium-term (6-months-ahead), and long-term (12-months-ahead) forecast horizon, respectively

Figure 10: SHAP values of the predictors: ridge logistic regression
The ranking of features may vary over time, but it offers an initial understanding of the key drivers influencing the prediction results for each forecast horizon. Depending on the forecasting horizon, variables associated with financial conditions such as S&P 500 index and term spread, and macroeconomic variables related to GDP, inflation, and house market, which are commonly regarded as important recession indicators, consistently rank high (Estrella and Mishkin (1998)). These findings provide valuable insights into the black box and suggest that neural networks effectively capture the effects of significant recession indicators.

Figure 10 presents the SHAP results for ridge logistic regression models. There are three major differences, among others, between the results in Figure 9 and 10. Firstly, the order of variables differs for each forecasting horizon between the two types of models. While some variables appear in the leading groups for both models, their specific order varies significantly. For instance, the money supply M2 is influential in the ridge model but less prominent in the GRU model. Spearman's rank correlations between the variable series for the two model types range from -0.24 to 0.10, indicating weak correlation. This suggests that GRU and Ridge models assign different weights to the variables. This difference could be attributed to the neural network's ability to capture nonlinear relationships, allowing for varied weighting of the variables based on different patterns. Secondly, the SHAP values of the GRU models exhibit a more even distribution among the variables compared to the ridge model. This is reflected in the smaller variation in medians between the variables. The distinction becomes particularly pronounced in the short-term setting (up to 3 months), where the ridge model heavily depends on the S&P 500 index. Additionally, the ridge model exhibits numerous outliers in the SHAP values, indicating its high sensitivity to changes in macroeconomic and financial conditions in SHAP value estimation. In contrast, this highlights the GRU model as a more robust modelling framework. Finally, the within variation of the GRU-based SHAP values for some variables are significantly larger than those of the others. Certain variables demonstrate little change in their SHAP values, while others exhibit significant variation. This phenomenon is less prominent in the outcomes of the ridge model. The variances are more uniformly spread, evident in the roughly equal sizes of boxes across variables. This could be viewed as another indication of the GRU's more resilient modeling framework. The GRU is capable of adapting to shifts in economic conditions, adjusting the weights of specific variables accordingly. In contrast, the ridge model has relatively fixed variable orders, with limited flexibility for variation.
Figure 11: Marginal effects of the predictors: ridge logistic regression
The figure presents the boxplots of the GRU-based absolute average LIME values of the predictors in descending order: Panel (A) to (E) display absolute average LIME values of the predictors related to the nowcasting, immediate-term (1-month-ahead), short-term (3-months-ahead), medium-term (6-months-ahead), and long-term (12-months-ahead) forecast horizon, respectively.

Figure 12: LIME values of the predictors: GRU
The marginal effect analysis is another option to evaluate the feature importance of linear models. Since deep learning models can not be interpreted in terms of marginal effects, I compute the average marginal effects of the ridge model for the out-of-sample period and plot them in Figure 11. The Spearman's rank correlations between the variable series for the ridge model, evaluated using either SHAP values or marginal effects, range from 0.22 to 0.51. The magnitudes of the correlations are not especially large to conclude that both metrics are the same, but the previous findings, especially the first one about the leading recession indicators, also apply to the Ridge model evaluated by the marginal effects. The averaged LIME values in Figure 12 which may be regarded as the marginal effect version of the GRU model reaffirm the main findings discussed above. The only minor distinction is that it places greater emphasis on the term spread as the comprehensive recession indicator, both in short-term and long-term contexts.

The SHAP method allows for a detailed analysis of the impact of a single variable on predictions. In order to anticipate recessions as early as possible, I focus on the long-term forecasting horizon

Figure 13: Dependence plot: Term Spread
(12 months) and consider the term spread as the most relevant and influential variable to examine its effects on recession predictions for the GRU model. Figure 13 presents the dependence plot of the term spread using the latest available data from October 2021. The x-axis represents the standardized value of the feature (term spread) while the y-axis represents the corresponding SHAP value. Notably, there is a cutoff point around 0.2, above which a higher value of the term spread has a negative impact on the predicted recession probability, while a value below the cutoff positively influences the predicted probability. Additionally, the farther the value deviates from the cutoff, the stronger the effect of the term spread on the predicted probability becomes. This corroborates the traditional interpretation of the term spread in relation to the probability of a recession according to Estrella and Mishkin (1996) and indicates that advanced neural network models like GRU can also uncover and utilize this relationship.

## 6 Conclusion

This research focuses on examining the real-time predictability of neural network models, compared to linear models, with respect to the recent two US recessions, the Great Recession and the Covid-19 recession. Three different neural network models are trained and updated throughout the study: a standard feed-forward neural network model, as well as two recurrent neural network models (LSTM and GRU) designed to capture temporal dependencies in time series data. The performance of these models is evaluated using out-of-sample forecasts and compared to standard and ridge logistic regression models. Additionally, the SHAP method is utilized to rank the predictors based on their importance for each forecast horizon, providing initial insights into the most influential predictors. The results are then compared to SHAP results obtained from the ridge logistic regression model. The main finding are validated using both the LIME approach and the marginal effect method. For in-depth analysis, the most influential variable for the long-term forecast horizon, which is the term spread, is chosen to investigate the impact of the variable on the recession probability.

This paper introduces two main contributions. Firstly, it focuses on LSTM and GRU, specialized recurrent neural network models that address issues like exploding and vanishing gradients in standard recurrent neural networks. The performance of these models is compared with simple feed-forward neural networks, which face the challenge of specifying temporal dependence in
advance, and with traditional linear models in the context of predicting recessions. Secondly, the paper employs the SHAP method to assess the importance of features in GRU for different forecast horizons. The three main findings are as follows: Firstly, LSTM and GRU demonstrate strong out-of-sample performance in recession forecasting, particularly excelling in long-term predictions across five forecasting horizons based on various statistical performance metrics. Secondly, there are differences in how GRU and ridge logistic regression models evaluate the importance of predictors. The variable order differs for each forecasting horizon between GRU and Ridge models. While some variables are significant in both models, their specific ranking varies notably. This suggests that GRU and Ridge models assign different weights to variables, potentially due to the neural network's capacity to capture nonlinear relationships and assign varied weights based on distinct patterns. Furthermore, the SHAP values of GRU models show a more balanced distribution among variables compared to the ridge model. This is evident in the smaller variation in medians between variables. The difference is particularly noticeable in the short-term scenario (up to 3 months), where the ridge model heavily relies on the S&P 500 index. Moreover, the ridge model displays numerous outliers in SHAP values, indicating high sensitivity to changes in macroeconomic and financial conditions. In contrast, this underscores the GRU model's robustness. The within variation of GRU-based SHAP values for some variables is significantly larger than for others. Certain variables show minimal change in SHAP values, while others exhibit considerable variation. This contrast is less pronounced in the ridge model outcomes, where variances are more evenly distributed, reflected in similar box sizes across variables. This suggests that the GRU has a more adaptable modeling framework, capable of adjusting variable weights in response to shifts in economic conditions. Conversely, the ridge model has relatively fixed variable orders, with limited flexibility for variation. Lastly, although the primary predictors for GRU and ridge logistic regression models show slight differences, key indicators such as the S&P 500 index, real GDP, and private residential fixed investment consistently play a significant role in short-term predictions (up to 3 months). For longer-term forecasts (6 months or more), the term spread and producer price index become more prominent. These results are supported by LIME and marginal effects, respectively.