{'fraction_non_alphanumeric': 0.02710210146824157, 'fraction_numerical': 0.01052000574125514, 'mean_word_length': 4.875390644377201, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Machine Learning, the Treasury Yield Curve and Recession Forecasting

We use machine learning methods to examine the power of Treasury term spreads and other financial market and macroeconomic variables to forecast US recessions, vis-à-vis probit regression. In particular we propose a novel strategy for conducting cross-validation on classifiers trained with macro/financial panel data of low frequency and compare the results to those obtained from standard k-folds cross-validation. Consistent with the existing literature we find that, in the time series setting, forecast accuracy estimates derived from k-folds are biased optimistically, and cross-validation strategies which eliminate data “peeking” produce lower, and perhaps more realistic, estimates of forecast accuracy. More strikingly, we also document rank reversal of probit, Random Forest, XGBoost, LightGBM, neural network and support-vector machine classifier forecast performance over the two cross-validation methodologies. That is, while a k-folds cross-validation indicates that the forecast accuracy of tree methods dominates that of neural networks, which in turn dominates that of probit regression, the more conservative cross-validation strategy we propose indicates the exact opposite, and that probit regression should be preferred over machine learning methods, at least in the context of the present problem. This latter result stands in contrast to a growing body of literature demonstrating that machine learning methods outperform many alternative classification algorithms and we discuss some possible reasons for our result. We also discuss techniques for conducting statistical inference on machine learning classifiers using Cochrane’s Q and McNemar’s tests; and use the SHapley Additive exPlanations (SHAP) framework to decompose US recession forecasts and analyze feature importance across business cycles.


Introduction
It is well documented that an inverted Treasury yield curve is a strong signal of recession in the United States. In May 2018, the term spread between the 3-month Treasury bill discount and 10-year on-the-run yield-to-maturity -one of several common measures of the yield curve's slope -narrowed to less than one percentage point for the first time since the Great Recession. The spread continued to decline and in May 2019 the yield curve inverted; that is to say, the 10-year yield-to-maturity fell below level of the 3-month bill discount. At time of writing (October 2019), the yield curve's slope is 11 basis points.
The long-standing literature on forecasting recessions from financial market and macroeconomic data has used probit classification as its workhorse method. The predictive power of term spreads to forecast recession during the Volcker Era was particularly strong and the results of a probit regression of Treasury term spreads against an indicator of recession over a one-to four-quarter horizon will bear that out. Studies dating to the early 1990's were first to make this case formally.
Virtually every study of this topic following those first published in the 1990's has used some version of the probit method to investigate the problem, each extending the original studies in a new direction. Some researchers have sought out new data to augment the feature space, and credit or more broadly financial conditions appear to be an important input now, particularly during the Great Moderation. Other researchers, rather than seek out new features, have extended the probit method itself to test and control for a variety of the time-series properties of terms spreads and other inputs over recession.
In addressing the shortcomings of the standard probit method, all previous efforts essentially contend with one simple fact, which is that the probit framework is rigid. That rigidity can be credited to the probit link function. In order to classify data, the probit method attempts to draw a hyperplane through the feature space by minimizing a cost function. In doing so, the optimization must trade off goodness of fit in some areas of the feature space (by fitting the hyperplane separator well in those regions) against poorness of fit in other areas.
Machine learning methods -rooted in statistics and computer science -are an attractive alternative to the probit classification methods traditionally used in time-series macroeconometrics, particularly those that extend probit regression in complex ways (so as to address dynamic coefficient variation, structural breaks and serial correlation, etc.). Whereas probit methods must typically introduce additional parameters in order to create model flexibility, with machine learning flexibility is more often an innate feature of the methods, saving the researcher the task of reworking model specifications or the econometric approach when results are not satisfactory. That said, machine learning methods have the tendency to over-fit data, unless controlled sufficiently via a set of hyperparameters, which is not a problem that arises when probit methods are applied. Thus, the flexibility of machine learning methods also presents the researcher with new difficulties, which are managed via a bias-variance tradeoff.
In this paper we attempt a survey of machine learning methods and investigate their performance relative to probit methods in forecasting US recessions from financial market and macroeconomic data. Specifically we study artificial neural networks, support-vector machines and the Random Forest, XGBoost and LightGBM algorithms (the latter three being presently some of the most popular tree ensemble methods in use). In service of this goal, we propose a novel strategy for conducting cross-validation on classifiers trained with macro/financial panel data of monthly frequency and compare the results to those obtained from standard k-folds cross-validation. We also discuss techniques for conducting statistical inference on machine learning classifiers using Cochrane's Q and McNemar's tests; and use the SHapley Additive exPlanations (SHAP) framework of Lundberg and Lee (2017) to decompose US recession forecasts and analyze feature importance across business cycles.
In a slight preview of our results we find that, consistent with established results, forecast accuracy estimates derived from k-folds cross-validation in time series settings are biased optimistically. Furthermore, we confirm that more conservative cross-validation strategies which eliminate data "peeking" produce lower, and perhaps more realistic, estimates of forecast accuracy. More strikingly, we also document rank reversal of probit, Random Forest, XGBoost, LightGBM, neural network and support-vector machine classifier performance over the two cross-validation methodologies we use. That is, while a k-folds cross-validation indicates that the forecast accuracy of tree methods dominates that of neural networks, which in turn dominates that of probit regression, the more conservative crossvalidation strategy we propose indicates the exact opposite, and that probit regression should be preferred over machine learning methods, at least in the context of the present problem. This latter result stands in contrast to a growing body of literature demonstrating that machine learning methods outperform many alternative classification algorithms and we discuss some possible reasons for our result.
This paper is structured as follows. In the next section we survey the literature on forecasting recession from financial market data using probit methods, as well as the related machine learning literature. In Section III we describe the data used throughout the study. In Section IV we describe how the features in the data set are combined into a nest of models, to which classification algorithms (probit regression, neural networks, support-vector machines, Random Forest, XGBoost and LightGBM) are applied, and in Section V we describe the classification algorithms themselves. In Section VI we describe a novel cross-validation strategy for hyperparameter selection and forecast performance estimation, as well as our approach to classifier comparison and statistical inference. In Section VII we discuss the methods used to calculate feature importances. In Section VIII we present the results of cross-validation and compare classification methods. Those results are followed in Section IX by a discussion of the some of the salient features and economic implications of our findings, interpretation of the machine learning classifier outputs and forecast attribution using the SHAP framework. In Section X we conclude.

II. Review of the Literature
Prior to the 1990's, econometric research examining the predictive power of financial variables was focused on forecasting macroeconomic outcomes such as real output and inflation. Much research since then -and that most relevant to the current paper -has examined the predictive power of financial variables to forecast future states of recession, rather than continuous measures of activity. Mishkin (1996, 1998) are early and oft-cited works that use financial and macroeconomic variables in a probit framework to forecast recession. They find that "stock prices are useful with one-to three-quarter horizons, as are some well-known macroeconomic indicators," for predicting recession. For longer horizons, however, they conclude that "the slope of the yield curve emerges as the clear individual choice and typically performs better by itself out-of-sample than in conjunction with other variables." Indeed, for the recessions occurring during and before the Volcker Era, negative terms spreads turned out to be a very strong signal of impending recession.
The slope of the yield curve largely failed to predict the 1990-91 recession, however, or at least not as strongly as it had those before in the 1970's and 1980's. Dueker (1997Dueker ( , 2002 uses Markov switching in the probit framework to allow for coefficient variation and also investigates issues surrounding the application of probit methods to time-series data. He finds that, while it is important to allow for dynamic serial correlation of term spreads in the probit framework, "allowance for general coefficient variation is not particularly significant at horizons less than one year." Just prior to the onset of the 2001 recession, Chauvet and Potter (2001) extended the probit method to investigate the instability of the term spread's predictive power and the possible existence of structural breaks. Using Bayesian techniques and several specifications of the probit method that variously allow for business-cycle dependence and autocorrelated errors, they find that "allowing for business-cycle-specific innovation variance and an autoregressive component has a much better insample fit than the original model of [Estrella and Mishkin]." They also find that recession forecasts from their more complicated extensions of the probit framework "are very different from the ones obtained from the standard probit specification." As the Great Recession loomed, and as the yield curve was beginning to invert again for first time since the 2001 recession, Wright (2006) re-examined the standard probit specification, looking more closely at the policy instruments and other financial variables of interest to the Federal Open Market Committee (FOMC). Including the level of the effective federal funds rate, the real federal funds rate (inflation) and a measure of term premiums in his analysis, he finds that adding the effective federal funds rate to a measure of the term spread improves predictive power and that "there is more information in the shape of the yield curve about the likely odds of recession than that provided by term spreads alone." Around that same time, King, Levin and Perli (2007) embed credit spreads in the probit framework along with term spreads and find that predictive ability is improved, and type I error is reduced. 1 They also incorporate Bayesian model averaging into the analysis, in another break with the standard probit framework, and find that "optimal (Bayesian) model combination strongly dominates simple averaging of model forecasts in predicting recessions." More recently, Fornani and Lemke (2010) extend the probit approach by endogenizing the dynamics of the regressors using a VAR and study the US, Germany and Japan. Liu and Moench (2016) use the receiver operating curve to assess predictive performance of a number of previously proposed variables. They find that, while "the Treasury term spread has the highest predictive power at horizons four to six quarters ahead, adding lagged observations of the term spread significantly improves the predictability of recessions at shorter horizons." Favara, Gilchrist, Lewis and Zakrajsek (2016) decompose credit spreads and show that their power to predict recession is contained in a measure of investor risk appetite called the Excess Bond Premium (EBP), which is a feature used in this paper. Finally, Johansson and Meldrum (2018) use the principal components of the yield curve, as well as a measure of term premiums to predict recession and Engstrom and Sharpe (2018) investigate the forecast power of nearterm forward term spreads.
The methods used in the research described to this point (probit regression, Markov switching, Bayesian techniques, etc.) are well established in the field of econometrics. In contrast, the methods that we use in this paper have roots in statistics and computer science and -though used in industrial applications such as for credit scoring, actuarial modeling, online marketing, demand forecasting, etc.have only in recent years found application in macroeconometric analysis. Fornaro (2016) uses large panels of predictors (several to hundreds) and adds to the probit framework a Bayesian methodology with a shrinkage prior for the parameters to predict recession. Ng (2014) dispenses with the probit framework altogether and applies a tree ensemble classifier 2 to a panel of 132 real and financial features and their lags to do so. In contrast to Ng, in this paper we study a small panel of just 9 features and investigate the recession forecasting power of a broad cross-section of machine learning classifiers vis-à-vis probit methods. As previously mentioned, the methods we explore are artificial neural networks, support-vector machines, and the Random Forest, XGBoost and LightGBM algorithms. Finally, Holopainen and Sarlin (2017) use many of the machine learning methods used in this paper (and more) for the purpose of creating an early-warning/crisis detection mechanism for the Euro area. After conducting a horse race between the methods, their focus turns to model aggregation and ensemble voting techniques for producing forecasts from multiple classifiers. They also investigate the statistical significance of their results by means of bootstrap procedures.
In contrast to Holopainen and Sarlin's work, in this paper we propose a very conservative nested time-series cross-validation procedure and explore strategies for contending with the time-series properties of macro panel data containing multiple structure breaks. Furthermore, we conduct statistical inference and classifier comparison by means of joint omnibus (Cochrane's Q) and pairwise post hoc (McNemar's) tests, rather than via bootstrap methods. We also apply the SHapley Additive exPlanations framework of Lundberg and Lee (2017) to investigate feature importance and decompose the forecasts of our classifiers.
The Random Forest method is a stable, well-studied tree-based ensemble algorithm used for both classification and regression. Its origins can be traced as far back as the 1960's, when Morgan and Sonquist (1963) published the first paper on decision tree analysis, describing their AID (automatic interaction and detection) algorithm. In the early 1980's, Breiman et al (1984) published the CART (Classification and Regression Trees) algorithm and reignited interest in decision trees. Decision tree classifiers suffer from several shortcomings however. Notably they are prone to over-fit and are sensitive to small perturbations in data, which render it difficult to balance bias and variance with them in supervised learning applications.
The Random Forest method proposed by Breiman (2001) corrects for these faults by introducing randomness into the decision tree algorithm and then aggregating the results of many decision trees. Each decision tree in a random forest is fit with a random subset of features from the pool of predictors and is trained on data sampled with replacement from the training set, a process known as bagging (bootstrap aggregating). Breiman shows that the Random Forest is robust to over-fitting and limits the bias traded against variance in supervised learning applications in ways that decision trees cannot.
The XGBoost (eXtreme Gradient Boosting) algorithm of Chen et al (2016) is a recent example of a general gradient boosting algorithm (to contrast with bagging algorithms like Random Forest). Gradient boosting methods were first proposed by Friedman (2001Friedman ( , 2002 and Mason et al (1999aMason et al ( , 1999b, and can also be used for both classification and regression. Tree-based boosting classifiers build decision trees sequentially, re-weighting misclassified instances on each iteration, and employ gradient descent to minimize a loss function. The AdaBoost algorithm of Freund and Schapire (1999) is an early example of a boosting algorithm that has found widespread use, but now it is but one of many and a growing number of such algorithms. XGBoost has been shown to possess a number of advantages over previous alternatives, notably speed and scalability.
Like XGBoost, LightGBM is a tree-based gradient boosting algorithm developed at Microsoft DMTK by Ke et al (2017) and open-sourced for public use. Though newer than XGBoost and not as widely used or tested, it has garnered interest in machine learning applications for its speed and, because its hyperparameter interface is very similar to XGBoost's, ease of use.
Artificial neural networks trace their history the 1940's. Early developments include Rosenblatt's perceptron (1958), a simplified mathematical model of how biological neurons operate, which was subsequently implemented in computer hardware. Though developments towards today's neural network architecture continued through the 1960's, interest cooled considerably after Minsky and Papert (1969) proved the limits of the direction research on the topic had taken. During a relatively dormant period of work in the area, Werbos applied backpropagation (1974) to neural networks, which along with gradient descent forms the basis of modern-day training methods. Research interest gradually picked up again in the 1980's with the advent of distributed, parallel computing, or connectionism as it was called at the time. Further development of the technology continued through the rest of the 20 th century, but it has been the growing availability of commodity computing hardware over the past two decades that has made rapid advancement and widespread commercial application possible. Many open-source software packages for developing and training neural networks now exist, including scikit-learn, Tensorflow (Google) and PyTorch (Facebook) among others. Support-vector machines were originally invented by Vapnik and Chervonenkis (1963) and grew in popularity after Boser, Guyon and Vapnik (1992) proposed the kernel trick for creating non-linear classifiers. Cortes and Vapnik (1995) introduced the soft margin classifier for handling noisy data. For a time these methods enjoyed a high degree of popularity among machine-learning practitioners, but more recently they have been surpassed in power by neural networks for many applications.
The empirical performance of machine learning methods on problems outside of the fields of finance and economics are well-documented. Fernandez-Delgado et al. (2014) studied 179 different binary classifiers on the 121 datasets in the UCI Machine Learning Data Repository. 3 They found that, on average, Random Forest methods achieved the highest accuracy among all other families of classifiers, followed by support-vector machines and neural networks, and that in many cases the differences at the top of the list where not statistically significant. A similar study by Wainer (2016), using more robust hyperparameter searches, found Random Forest and gradient boosting to be two of the three top performing methods for classification, with the differences between the two methods also not statistically significant. Neither study utilized the (more recent) XGBoost or LightGBM packages specifically, but the former is currently the dominant statistical package in online data science competitions hosted by Kaggle 4 .

III. Data
The primary focus of this paper is to compare probit regression -used in most of the existing literature on forecasting recessions from financial market and macroeconomic data -to machine learning classifier algorithms. Our purpose is not to identify new features for forecasting recessions. As such, most of the features we have chosen for the analysis (6 of 9) are fairly standard in the literature and our results for models estimated via probit methods should be fairly well known. That said, we have experimented with two or three features that have not received much use in the published literature on this topic.
For a measure of terms spreads, or yield curve slope (Slope hereafter), we use the monthly average of the 10-year Treasury spot yield of Gurkaynak, Sack, and Wright (2006) less the monthly average of the 3-month Treasury bill discount in the Federal Reserve's H.15 series. Consistent with Estrella and Mishkin (1998), we use the 3-month log difference of end-of-month S&P 500 index values (SP500 hereafter) to capture equity market developments. For a measure of credit market conditions, we use the Excess Bond Premium (EBP hereafter) of Gilchrist and Zakrajsek (2012), which is a monthly series. As in Wright (2006), the end-of-month values of the effective federal funds rate (FF) is used as the short-term rate of interest. In addition, we include the Federal Reserve Bank of Chicago's National Financial Conditions Index (NFCI hereafter) as a measure of financial conditions more broad than that captured by the Excess bond premium and to our knowledge this is its first use in the literature.
In addition to these financial market inputs, we augment our feature space with two macroeconomic/business conditions indexes and a survey-based measure of recession probability. The 3-month log difference of the Conference Board's Leading Economic Index (LEI hereafter), previously studied by Estrella and Mishkin (1998) and Dueker (1997) in the context of the present problem, is used in this study to capture macroeconomic conditions writ large. The highly correlated monthly average of the Aruoba-Diebold-Scotti business conditions index (ADS hereafter), which is comprised of several of the same macroeconomic data series 5 , is also used 6 . For our survey-based measure, we use the cumulative probability of recession over four quarters calculated from the Federal Reserve Bank of Philadelphia's Survey of Professional Forecasters 7 (SPF hereafter). 4 https://www.kaggle.com/ 5 Specifically, the index is comprised of (seasonally adjusted) weekly initial jobless claims, monthly payroll employment, industrial production, personal income less transfer payments, manufacturing trade and sales and quarterly real GDP. 6 Note that the ADS index -like the financial variables -is available in real-time, while the LEI and SPF are only available at a lag. This motivates our use of the ADS index, despite its high correlation (77%) with changes in the LEI. We must credit Viktors Stebunovs with introducing us to this feature and its use in the context of recession forecasting. 7 That is, the cumulated probability of recession assuming independence between quarters calculated from RECESS2 through RECESS5 Finally, we add a measure of term premiums to the analysis, with one small twist. While Wright (2006) and Johansson and Meldrum (2018) included measures of term premiums in their analysis (Cochrane and Piazzesi's return forecasting factor and the 10-year term premium of Kim and Wright, respectively), here we use changes in term premiums, rather than levels. Specifically we use the 6month percentage point change in the 5-year term premium of Adrian, Crump and Moench (ACM, 2013). We have chosen to use changes in term premiums (hereafter ACM5D) due to the well documented fact that the most widely used term structure "models produce different estimates for the levels of the term premia, but broadly agree on the trends and dynamics" (Cohen et al, 2018). We have chosen to use the ACM model specifically due to the widely held view that term premium estimates derived from the model are less persistent than alternative models (and so, in our view, perhaps more able to capture changes in the outlook for recession in a timely manner).
The recession indicator used in the analysis is the same as that used in most of the existing literature as well. For any given month, it is defined as true if any of the following twelve months falls within a recession, as defined by the National Bureau of Economic Research (NBER), and is false otherwise. 8 All data is monthly, and covers the period from January 1972 9 to June 2019 (558 months or observations). Figure 1 below summarizes the nine data series. Shaded regions indicate recessions. All series are monthly. 3m/10y Slope is the monthly average of the 10-year Treasury spot yield of Gurkaynak, Sack, and Wright (GSW) less the monthly average of the 3-month Treasury bill discount in the Federal Reserve's H.15 series. Excess Bond Premium (EBP) is that of Gilchrist and Zakrajsek (2012). National Financial Conditions Index (NFCI) comes from the Federal Reserve Bank of Chicago. ADS Index is the monthly average of the Aruoba-Diebold-Scotti Business Conditions Index. LEI is the 3-month log-difference of the Conference Board's Leading Economic Index. S&P500 is the 3-month log-difference of the equity index. The effective federal funds rate is month-end. SPF is the cumulative 4-quarter probability of recession calculated from the Federal Reserve Bank of Philadelphia's Survey of Professional Forecasters RECESS series. ACM 5-year term premiums is a month-end value from the model of Adrian, Crump and Moench. The 6-month change in the ACM 5-year term premium of Adrian, Crump and Moench is a first difference in levels of average daily values. Shaded regions indicate recessions. Source: S&P 500 data from Bloomberg Finance LP, Bloomberg Per Security Data License.
The unconditional probability of the indicator across the entire dataset is 25.4% while the unconditional probability of recession in any given month over the entire dataset is 14.0%. The static moments, monthly volatility 10 , and the end-of-sample (June 2019) values of each series are summarized in Table 1 In order to see more clearly the signal in the noise, Figure 2 shows the same data when each series is truncated at the level of 50% probability in a univariate probit regression against the recession indicator 11 , then rebased to zero. The last row of Table 1 lists the 50% probability level for each feature. Recessions prior the mid-1980's were preceded by periods of strong inversion in term spreads and high short-term interest rates, while recessions after 1991 have been preceded by high levels of EBP (which is closely related to adverse credit conditions and/or credit investor risk aversion). 12 The NFCI is partially correlated with the EBP, but fails to predict the recessions of 1991 and 2000. The ADS, LEI and SPF features appear to capture the 1991 recession more strongly than did the financial variables. Given the small, sparse nature of the data set, these facts likely explain most of its predictive power. 11 For the slope regression, this level is 34 basis points (i.e. yields curves flatter than 34 basis points signal greater than 50% probability of recession in the following 12 months.) For the EBP regression, this level is 61 basis points (i.e. EBP greater than 61 basis points signals greater than 50% probability of recession in the following 12 months). For the ADS regression, this level is -0.81 (i.e. an index level less than -0.81 signals greater than 50% probability of recession in the following 12 months). And so on… 12 Estrella and Mishkin (1996) noted that "restrictive monetary policy probably induced the 1973-75, 1980 and 1981-82 recessions but it played a much smaller role in the 1990-91 recession." The latter could perhaps also be said of the 2001 recession and Great Recession as well.

IV. Models
Let us begin the present discussion by clarifying terminology. In machine learning studies generally, the concepts of a model, an algorithm and a classifier can be difficult to distinguish and the terms are often used interchangeably. Furthermore, there is a great deal of precedent in the fields of economics and finance that does not conform to conventions used in the disperse fields from which machine learning methods are derived. Here we use the following nomenclature: • Model -A specification or collection of features (possibly transformed) representing the underlying economic system (e.g. a linear combination of variables, such as in the right hand side of a linear regression specification, or just a loose collection of features to be fed as one into a machine learning algorithm). Also called a hypothesis.
• Algorithm -A strategy or set of instructions for approximating a target function using a model specification and data observations. In the present case, the target function is a recession forecast probability (or equivalently one that maps a recession label to a data point). Examples of classifier algorithms would include probit regression, neural networks or any other machine learning or econometric method. Also called a classifier algorithm or a learning algorithm. For practical purposes, algorithms exist as software packages to which researchers are able to apply models and data for the purpose of creating classifiers.
• Classifier -The result of applying a model and data to a classifier algorithm. A classifier is an object that resides in computer memory after estimation is complete. In econometric terms, it is a sample from the space defined by a model and an algorithm over the data population. In order to compare the performance of two models or two algorithms, two classifiers must be estimated and the comparison must be made through the classifier samples. In contrast to the previous two concepts, which are abstract concepts for the present purposes, a classifier is a tangible machine that an econometrician can work on and study.
In this paper we study a large number of nested models, which is the subject of the present section. Denote the feature values in month t as , , , , 500 , , , , 5 . Furthermore, denote the response indicator, which takes a value of 1 if there is an NBER-dated recession in any of the twelve months following month t and is 0 otherwise, as +1, +12 . Finally, denote Φ() as the probit link function (the standard normal cumulative distribution), () as neural network classifier algorithm over a collection of features, () as a support-vector machine algorithm over a collection of features, () as a Random Forest binary classifier algorithm over a collection of features, () as an XGBoost binary classifier algorithm over a collection of features, and () as a LightGBM binary classifier algorithm over a collection of features.
The first batch of nine models are univariate in each of the nine features. In the case of the Slope model, the probability of the indicator conditional on the model features for each classifier (i.e. the Slope classifiers to be estimated from the Slope model using the learning algorithms) is given by: Note that a constant is included when the probit method is used.
In the second batch of eight bivariate models, the Slope feature is combined with each of the other eight features. In the case of the Slope/LEI pair, the probability of the indicator conditional on the features for each classifier (i.e. the Slope/LEI classifiers to be estimated from the Slope/LEI model using the learning algorithms) is given by: The third batch of tri-variate models combines Slope with two other variables that scored highly among the batch of 2-feature models. We continue in this fashion, forming batches of 4-, 5-, 6-, 7-and 8-feature models using the features appearing in the highest-performing models from the previous batch and combining with each of the remaining features. Due to computing constraints it is not feasible to estimate classifiers on models with all combinations of 1 to 9 features, so some subjective judgement is required between batches to identify those that are likely to perform well or otherwise be of interest to the analysis.
In total we compare 106 models, each estimated with all six algorithms using two crossvalidation strategies, producing a total of 1272 classifiers. The models are listed in the Appendix in Table  13. Many of the models have been studied previously in the literature, particularly those that combine some of the Slope, LEI, EBP, SP500, FF and SPF features.

V. Algorithms
In this section, we provide a more detailed look at the machine learning methods studied in this paper, and provide a more rigorous mathematical treatment of the training and prediction processes. Specifically, we cover the Random Forest, boosting (XGBoost and LightGBM), support-vector machine, and neural network algorithms. That said, we only seek to provide a greater intuition for these so-called "black box" methods, and encourage the interested reader to explore the cited papers more fully. In order to keep mathematical notation consistent, for the remainder of this section denote as the design matrix of dimension × where each of the columns is a feature vector with observations. Similarly, define = ( 1 , 2 , … , ) as the binary response vector. The classical binary response is typically coded as ∈ {0,1}; however, for the mathematical formulations of boosting and support vector machines, each ∈ {−1,1}. Finally, let ( ′ , ′ ) be an arbitrary loss function with some inputs ′ , ′ . Further notation will be introduced as needed.
As previously introduced, the Random Forest, formulated by Breiman (2001), is an ensemble of many decision trees, each trained over random bootstrap samples from the original data. Bootstrap resampling of training data is a process known as bagging (as in bootstrap aggregating) in the machine learning literature. In addition to bagging, the Random Forest method employs another process called feature bagging to further decorrelate the trees. Feature bagging is the selection of a random subset without replacement of features to be considered in the fitting process for each decision tree. That is, only bagged features are considered at each split in the tree. Training of the decision tree (i.e. determining the splits in the tree) is accomplished via the Classification and Regression Tree (CART) algorithm (Breiman et al 1984). Illustrating the algorithm is most easily done through the following pseudocode, where we define ( ′ , ′ ) to represent a decision tree trained on some inputs ′ , ′ that optimizes some loss function ( ′ , ′ ): and are taken to be free parameters in this algorithm and are tuned via cross-validation. There are additional regulation parameters that can be used to further tune the algorithm and control the bias-variance tradeoff (for instance, to control how deep a tree may be grown or how many nodes a tree may contain) but details on these other parameters are omitted for the sake of brevity 13 .
The loss function most commonly employed by the Random Forest and other methods using the CART algorithm is Gini impurity, which measures the chance of a randomly chosen element being misclassified, given the distributions of classes in the problem. Each split in the tree ( , ) maximizes the amount of Gini impurity removed from the classification task over the variables considered. Other common choices of loss function are information gain, from information theory, for classification and variance reduction for regression.
Next, we present the general boosting framework first introduced by Friedman (2001Friedman ( , 2002 and Mason et al (1999aMason et al ( , 1999b to describe the XGBoost and LightGBM methods. Boosting employs an iterative gradient descent algorithm to approximate a function, ( ), that minimizes expected loss. This approximation to ( ), denoted by � ( ), is a weighted sum of many decision trees such that � ( ) minimizes the empirical risk. In contrast to bagging in the Random Forest, the boosting method, at a general iteration step , fits a new decision tree on the vector of "pseudo-residuals" in a greedy 14 fashion rather than on the response . More specifically, define the "pseudo-residuals" as: so that at step , the boosting tree fits ( , � ) rather than the bagged tree ( , ) in the Random Forest algorithm. In other words, boosting methods train on the error of the previous iteration, so as to pass more weight on to misclassified examples and less to correctly classified examples. Again, using pseudocode to describe this procedure from Friedman (2001): Algorithm 2: General boosting framework Create a base classifier to initialize the fitting process � 0 ( ) = ( , ) For s =1 to S do

End for
Classify new observations using � ( ) is taken to be a free parameter in this algorithm to be tuned during cross-validation. In most implementations there will exist other free parameters can be used to further tune the algorithm (for instance, to control the depth of trees or the maximum number of nodes in any tree). Similarly to Random Forests, boosting is adaptive to a variety of loss functions. The choice of function depends on the problem at hand -classification or regression -and on the utility of the econometrician. In this study, we utilize binary logistic loss for boosting methods.
The Random Forest and boosting methods share the free parameters and , which control the number of trees grown in the algorithms, as well as the depth and number of nodes in any particular tree. Each individual tree grown can be described as a weak learner in that the generalization of any tree is often poor. However, combining many of these weak learners in an ensemble is what drives the predictive power of the methods. Our choice of is the solution to a bias-variance tradeoff problem in classical statistics, with the added wrinkle that we must also consider computational feasibility when working with most boosting algorithms. A small may increase computational efficiency and bias, while keeping prediction variance low. Conversely, a large may increase the variance and computational load, but reduces bias. The following section presents more information on the various parameters encountered in the algorithms and our cross-validation strategy for tuning them.
Moving on we present the support-vector machine method originally proposed by Vapnik and Chervonenkis (1963). A support-vector machine classifier attempts to construct a pair of parallel hyperplanes to separate data with the largest possible margin between classes. For illustrative purposes, we consider a linear support-vector machine in 2-dimensions with data that is perfectly separable (or hard-margin, to contrast with soft-margin). In this simple case, there are many lines that define a margin separating the classes in the data; however the support-vector machine finds those which maximize the distance to the nearest data point in each class. Consider two parallel lines defined by: As shown in the figure below by the two dashed black lines, if the parameters of the lines are chosen appropriately, a margin between the classes in the data is formed. The distance between these lines is given by:

‖��⃗‖
To maximize the distance between the lines, minimize ‖��⃗‖ subject to the constraint (��⃗ • ⃗ − ) ≥ 1 (i.e. so that no data points fall inside the margin). The solution to this problem defines the classifier ⃗ → sgn(��⃗ • ⃗ − ). In this simple case, the support-vector machine is not that different from a probit regression, but the support-vector machine can be modified to be more flexible. For example, for data that is not separable, the hinge loss function can be introduced (soft-margin). For non-linear boundaries, kernel functions can be introduced. In our application, we search over 2 nd , 3 rd , and 4 th order polynomial kernel functions during cross-validation.
To conclude this section, we turn to the artificial neural network classifier. Artificial neural networks approximate non-linear regression and classification functions through a collection of nodes organized into multiple layers. The figure below depicts a simple example of a neural network classifier with three input vectors, two hidden layers with four nodes each and a single output node.

Figure 4 -Illustration of a feed-forward neural network.
This is a simple feed-forward neural network shown for illustrative purposes. It contains two hidden layers, with four nodes each. Data from the leftmost layer, the input layer, is passed into the hidden layers, which transform the data. Values from the hidden layers are passed to the rightmost layer, the output layer, which is an indicator of class probability in classification problems. Bias (constant) nodes may also be employed, but are omitted in the figure.
The first layer, the input layer, passes the data vectors 1 , 2 , 3 into hidden layer 1. Each node in hidden layer 1 uses a weighting vector (displayed as edges or arrows into the node) to linearly transform the information from the input layer, 0 + 1 1 + 2 2 + 3 3 , and then employs a (possibly) non-linear activation function f to construct the layer outputs at each node: Common choices of the activation function include the hyperbolic tangent function and the logistic function amongst others. In this study the rectified linear unit (ReLU) is used. This procedure continues to the right for each node in each hidden layer. The last layer of the neural network, the output layer, transforms the values from the final hidden layer into the output values of the network, which in the binary classification setting is a single node denoting the indicator class probability.
Like the neural network depicted above, in this study we use a relatively simple architecture, consisting of just two hidden layers. We make no attempt to cross-validate over the number of nodes in each layer (or over more exotic architectures such as recurrent or convolution networks or those employing other features such as dropout). The number of nodes used in each layer is a stepwise function of the number of features in the model being trained, as follows (neglecting the bias nodes):

# of Model Features, k
# Nodes, Hidden Layer 1 # Nodes, Hidden Layer 2 0<k≤3 9 5 3<k≤6 12 6 k≥7 18 9 Neural networks are trained on data (i.e. the weights are chosen) using backpropagation and gradient descent. Backpropagation is an iterative method that exploits the chain rule to update node weights in a neural networks so as to minimize a loss function. In this study cross-entropy loss is targeted.

VI. Classifiers and Cross-Validation
In this section we outline our approach to 1) selecting hyperparameters for the algorithms, 2) estimating the out-of-sample forecast performance of each classifier and 3) comparing classifiers across models and algorithms. Ultimately our goal is to determine what, if any, value the machine learning methods add to our ability to forecast recession vis-à-vis the probit method. The data set presents several challenges that motivate our overall strategy for doing so. Those challenges can be summarized as follows: -The data set is serially correlated, which violates the i.i.d. assumption required to use many of the cross-validation methods most popular in the machine learning literature, such as k-folds, leave-p-out, etc. Furthermore, because the data is time-series and not cross-sectional in nature, attempting k-folds cross validation on the data set would result in "data peeking" and overly optimistic estimations of forecast performance.
-The data set likely contains one if not multiple structural breaks (delineated by the Volcker Fed, the Great Moderation and the Financial Crisis perhaps), further violating the i.i.d. assumption and also making standard time-series cross-validation methods such as sliding or expanding windows problematic.
-The indicator (recession) is unbalanced in the data set (unconditional mean of ~25%) and occurs temporally in clusters of varying length at the end of each business cycle, further complicating the use of standard time-series cross-validation methods such as sliding or expanding windows.
-The data set is relatively short (558 observations) and sparse over the joint feature space.
Our strategy for contending with these features of the data is to implement a 4-fold nested time-series (NTS) cross-validation to estimate out-of-sample forecast performance for each combination of the 106 models and 6 algorithms under study. Originally designed by Varma and Simon (2006) for use on small datasets to address the unique difficulties they pose, according to Raschka (2018) nested crossvalidation "shows a low bias in practice where reserving data for independent test sets is not feasible." We augment the standard nested cross-validation strategy to incorporate several features that make it more amenable to conducting time-series analysis on the small macro/financial panel dataset under study in this paper, and to address the serial correlation, indicator imbalance and structural breaks that are present (hence the name nested time-series cross-validation). In particular, we overlay standard nested cross validation with an expanding window, so as to respect the time ordering of the data and prevent future leakage. We add one wrinkle to this feature in that, rather than forward chaining the window over a single data point or a fixed-size block of data points, we forward chain the outer loop of the NTS cross validation over business cycles.
The NTS cross-validation procedure we use is most easily described by means of the following pseudo-code algorithm: To be more concrete, consider the top panel of Figure 5, where we plot real GDP growth over our sample period, indicating recession periods in gray above the x-axis. Colored regions below the xaxis demarcate the individual business cycles that comprise the data; each colored region begins in a period of economic growth and ends when the NBER recession period following the period of growth terminates. In the top panel, Real GDP growth is plotted over the sample period. Gray shaded regions above the x-axis indicate recession periods. Colored regions below the x-axis demarcate the business cycles that comprise the data set. In the lower panel, the correspondence between the business cycles and the data folds used in cross-validation is shown.
The bottom panel of Figure 5 shows each iteration of the NTS cross-validation outer loop in a row. Within each iteration of the outer loop, a stratified/shuffled 3-fold cross-validation is conducted over each training set (labeled Train) to determine optimal hyperparameters. The optimal hyperparameters are then used to estimate a classifier over the full training set, and the outer loop iteration is scored by calculating forecast performance (e.g. accuracy) over the test set (labeled Test) using the estimated classifier. When all 4 iterations of the outer loop are complete, the model/algorithm pair is scored by the average out-of-sample forecast performance over the 4 iterations.
The first three business cycles are used for training in the first iteration of the outer loop due to their relatively short durations, and to ensure a large enough sample for classifier estimation. It is our opinion that the division of the data and the NTS cross-validation strategy described above represents the best balance that can be achieved between 1) size of the training/test sets on one hand, 2) sample size of the outer loop forecast accuracy estimate and 3) indicator imbalance in the data.
Two main alternatives to this cross-validation strategies exist, but when applied to time series data they become problematic, in that they do not always respect the time ordering of data. The first alternative is a standard k-folds cross-validation strategy, originally developed for use on cross-sectional data that obeys an i.i.d. assumption. It is widely recognized that this strategy is not appropriate for time series data and results is overly optimistic estimates of forecast performance. The second alternative is a hybrid approach, somewhat between a standard k-folds cross-validation strategy and the NTS crossvalidation strategy described above. While it is generally accepted that the time-ordering of data must be respected when estimating classifier forecast performance, it is less clear whether it is necessary to respect the time-ordering of data when determining the optimal hyperparameters of the classifier algorithm. A strategy in which optimal hyperparameters are chosen through a strategy such as k-folds cross-validation on the entire data set, but then forecast performance is estimated using a rolling or expanding window strategy and the optimal hyperparameters determined by k-folds, would respect the time-ordering of data during the latter half of the exercise, but violate it during the former. Holopainen and Sarlin (2016) implement a hybrid strategy similar to this and find ultimately that "conventional statistical approaches are outperformed by more advanced machine learning methods." In this paper we've taken a more conservative approach in using the NTS cross validation strategy, precluding any violation of the time-ordering of data during the study. In a later section we will discuss the possible implications of our choice, and we present results from NTS cross-validation and k-folds cross-validation side-by-side for comparison.
With regards to the forecast performance metric, we use accuracy as the target value over which hyperparameters are selected and out-of-sample forecast performance is measured. Our choice of forecast accuracy as a scoring metric within the NTS cross validation procedure is motivated by several considerations. Though we would have preferred to use average precision (AP) of the precisionrecall curve to score forecast performance, the last business cycle in the sample does not yet contain a positive instance of the indicator (recession), 15 and so AP is ill-conditioned in that business cycle during NTS cross-validation. We also considered area under the receiver operating curve (ROCAUC). Like AP, however, ROCAUC is also ill-conditioned when a positive instance of the indicator does not exist and in any case Davis and Goadrich (2006) show that, in this setting, the area under the receiver operating characteristic curve does not give a complete picture of a classification method's performance 16 . Since accuracy and Brier score are well-behaved over the last business cycle, they remain as the best candidates for scoring forecast performance. For the present analysis we have chosen to use accuracy as the target metric.
The grid search that is used in the inner loop of the NTS cross-validation procedure is the principle means by which variance is traded against bias for the machine learning algorithms. For example, if the tree-based methods (Random Forest, XGBoost, LightGBM) are allowed to grow too deeply, or if neural networks with a large number of neurons and hidden layers are used, the methods are prone to overfit data, resulting in high forecast variance, which is perhaps the main difficulty that arises in working with the algorithms.
For each of the tree-based algorithms, we search over a subset of about five of the hyperparameters exposed. Since each of the algorithms has a different set of hyperparameters it is difficult to draw direct comparisons between them, but generally there are a few main parameters used in each to control overfitting and balance speed of estimation against accuracy. The ensemble methods' hyperparameter "number of iterations" or "number of trees" dictates the number of weak learners created. Setting this number too low can cause models to not capture the signal in the data, but setting this number too high will almost certainly result in an infeasible computation time. For the Random Forest and XGBoost algorithms, the "maximum depth" hyperparameter (exposed in all of the major implementations of these algorithms) is a strong control. LightGBM also has a maximum depth hyperparameter, but a separate parameter controlling the maximum number of leaves in an individual tree is more readily able to control overfitting. For the 1-feature models, maximum depth is limited to be at most 3 in the grid search; for the 2-feature models it is limited to 4; while for all other models it is limited to 5. While this restriction may seem arbitrary, we found its interpretation and symmetry to be pleasing. 17 As a practical matter, we also found the restriction to be very necessary in light of our very small and sparse dataset.
Where the neural network and the support-vector machines are concerned, the situation is greatly simplified. We were able to obtain satisfactory results using neural networks with a simple two hidden-layer architecture of roughly 12 and 6 neurons respectively (slightly less for the models with just a couple features, and slightly more for models with seven or eight, as shown in Table 2). We made no attempt to optimize over the neural network architecture, though it is possible our results could be improved by doing so. For support-vector machines, the kernel function search was limited to 2 nd , 3 rd and 4 th order polynomials. Finally, note that, since the probit method does not have hyperparameters, grid-search cross-validation is not necessary and the inner loop can be forgone, which renders our strategy equivalent to time-series expanding windows (over business cycles) cross-validation in this special case.
After applying the NTS cross-validation strategy to our models, algorithms and data, we will be left with 636 estimates of forecast performance across all combinations of 106 models and 6 algorithms, from which we hope to choose those that work best for the problem at hand. Accounting for estimation uncertainty is critical to this exercise. Unfortunately this is a topic of research that is still developing in the machine learning literature and there does not yet exist a generally accepted approach to statistical inference or classifier comparison. In this paper we rely heavily on the strategies suggested in Raschka (2018).
The first stage of our model comparison amounts to what is effectively an omnibus F-test over our classifiers. Rather than group all of our classifiers into a single test, however, we group them once by model and again by algorithm, to identify the dimensions along which large, statistically significant differences in forecast performance may exist, and those along which they do not. More specifically, we use Cochran's Q Test, which tests the null hypothesis that there is no difference between the grouped classifier forecast accuracy estimates. Consider the group of six classifiers for a single model x (the 1feature Slope model, for instance). The null hypothesis 0 is: where , is the forecast accuracy of the classifier estimated on model x, algorithm y. If the = 6 classifiers for the model do not differ in forecast performance, then Cochran's Q statistic for model x, defined below, will be distributed approximately χ 2 with − 1 degrees of freedom: where = 558 is the total number of observations in the data set; is the number of classifiers out of m that correctly classified the j th observation in the data set; , , is 1 if the classifier for model x, algorithm y correctly identified the j th observation and 0 otherwise; and: Similarly we could group the classifiers by algorithm (neural networks for instance) and conduct the omnibus test in that dimension. In this case the null hypothesis 0 is: where again , is the forecast accuracy of the classifier estimated on model x, algorithm y. If the = 106 classifiers for the algorithm do not differ in forecast performance, then Cochran's Q statistic for algorithm y, defined below, will be distributed approximately χ 2 with − 1 degrees of freedom: and all other quantities are the same as for .
Note that in the case of , the grouping of classifiers by sixes keeps the individual tests reasonably small. In effect, the hypothesis that tests is along the lines of "does it matter which algorithm is used to estimate this model?" If the null hypothesis is rejected, it suggests that at least one classifier is greatly outperforming or underperforming the others, and that further investigation is warranted. In the case of , however, grouping 106 classifiers at a time is fairly pointless. In a sample so large it is almost certain that one classifier will perform differently than the rest. In order to address this issue (that is, to pose a more reasonable hypothesis) some a priori means of first filtering among the 106 models is required before applying Cochran's Q test along that dimension. In a later section we will discuss the means by which we select down from 106 to just 8 classifiers before applying the test.
If the null hypothesis is rejected during one of the omnibus tests, we have reason to conclude that a statistically significant difference in forecast performance exists among the classifiers. In the second stage of our model comparison, we conduct pairwise post-hoc tests, with adjustments for multiple comparisons, to determine where the differences occur. In particular, we use pairwise McNemar (1947) or Chi-square within-subjects tests between all classifiers in the grouping. To do so, we first construct 2x2 contingency tables for each pairing of classifiers in the group, shown below: In the contingency table, A is the number of observations out of n=558 that both classifier 1 and 2 classified correctly. B is the number of observations that classifier 1 classified correctly and classifier 2 identified incorrectly, and vice-versa for C. D is the number of observations that both classified incorrectly. Note that + + + = and that the accuracy of classifier 1 and 2 is: Note also that the differences between the classifiers are captured by B and C, so these are the natural quantities of interest to the problem at hand. If the null hypothesis is that the probabilities of B and C are equal (i.e. both classifiers have the same rate of error), then McNemar's statistic, defined below, will be distributed χ 2 with one degree of freedom.
Edwards (1948) proposed a continuity correction to this statistic, which becomes: If the null hypothesis is rejected after application of McNemar's test to a pairing of classifiers, we may conclude that the classifier with the higher forecast accuracy estimate outperforms the other and that the difference is statistically significant.

VII. Feature Importances
After the forecast performances of all models and algorithms have been determined in NTS cross validation and weighed against one another by means of the omnibus Cochran's Q and pairwise McNemar tests presented in the previous sections, it is natural to ask which of the nine features under study contributed to the overall performance of the classifier and how. In this section we outline a few methods for doing so.
At this point of the analysis we dispense with NTS cross-validation and calculate our feature importance metrics in-sample. That is, we use the entire data set as a training sample, use k-folds crossvalidation to determine optimal hyperparameters, then estimate classifiers using the optimal hyperparameters, and then calculate feature importances by means of the estimated classifiers. It is not clear whether this is a valid approach. We know that the (out-of-sample) forecasts from NTS crossvalidation differ greatly from in-sample forecasts, and forecasts estimates are an important input to some of the metrics we use (e.g. SHAP values) as will be shown below. We also know that it is preferable to calculate some importance metrics (e.g. permutation importance or mean decrease in accuracy) out-of-sample, but it is generally common practice to calculate them in-sample. Our reason for calculating feature importance in-sample is mainly practical. First, it would require a great deal of new machinery to integrate the metrics into the NTS cross-validation strategy we propose. Second, even if that were theoretically possible (which is not yet clear) we would have serious issues of computational feasibility to contend with. As such, we move forward using in-sample metrics from this point and leave those considerations to future work.
The feature importance metrics that we consider fall into two broad categories: those that are specific to only one or a few algorithms and those that can be applied broadly to all algorithms. In the probit literature referenced in Section II, it is common practice to report the sensitivity or marginal effects of an estimated model as a measure of variable importance. This is an example of the former type of feature importance, since it relies on the parametric nature of the probit method and cannot be applied directly to the machine learning methods we study. The sensitivity metric has many variations, but in essence all of them are partial derivatives of the output probability with respect to the inputs at some point in the feature space. One commonly used variation of the metric calculates these partial derivatives at each point in the dataset used for estimation and averages over all observations. That is, the sensitivity of the probit classifier to the k th feature f k is: where k=1…m; m is the total number of features in the model; α k , k=1…m are the estimated coefficients; Obs i is the set of feature values in observation i=1…n; and n is the total number of observations. Other variations may use the partial derivative evaluated at some other point in the feature space, such as the unconditional mean of the features.
A more useful importance metric for the present study is the permutation importance, or mean decrease in accuracy (MDA) metric. Originally proposed by Breiman (2001), the calculation of this metric uses bootstrap methods, allowing it to be applied broadly to any classifier. The algorithm is most easily described via the pseudo-code below: For each iteration in the outer loop of the algorithm, a fraction of the observations (perhaps 70%) are sampled from the data to create a bootstrap sample. Baseline accuracy for the sample set is then calculated and recorded, before entering the inner loop of the algorithm. On each iteration of the inner loop, a single feature's column is permuted to break the association between the feature and indicator in the data. Accuracy of the classifier on the bootstrap/permuted data set is then calculated, and the decrease in accuracy from the baseline level that results is recorded. When all iterations (perhaps 10 or 100, depending on sample size) of the outer loop are complete, the sampled decreases in accuracy for each feature are averaged over the number of outer loop iterations to produce the mean decreases in accuracy for each feature. Features with large MDA's indicate that classifier accuracy will suffer if removed from the model and re-estimated, and those with low MDA's indicate that classifier accuracy will not change greatly if removed.
An argument can be made that a classifier's MDA should not be calculated on the same data used to train the classifier, but rather a test or validation set that is held out. Due to the many data issues described already, and because we will have documented forecast performance prior to calculating these feature importance metrics, we calculate MDA on our entire data sample. While this may be cause for criticism, most of the recession forecasting literature uses in-sample data when calculating probit classifier sensitivities as well. We also feel that, while absolute forecast performance estimates will likely be overstated if the time ordering of the data is not respected during cross validation (necessitating the proper use of test sets and cross-validation methods that preclude data peeking) it is not clear that relative feature importances (which are probably more germane to this discussion than absolute feature importances) will be so greatly impacted if calculated using in-sample data.
Data issues aside, the feature importance metrics described so far -probit sensitivities and permutation importances -belong to a class of so-called global feature attribution methods, in that they score the value of a feature to a classifier as a whole through a single summary metric. At the core of the SHAP framework, the Shapley value is a concept from coalitional game theory that prescribes how gains or "payouts" resulting from a cooperative game should be divided amongst players in proportion to their contribution to the outcome. In the present case, the players are the features in the dataset and the "payout" is the forecast probability output by the classifier given the observations at time t 18 . More technically, consider the set f of m features in the dataset at time t, and a classifier prediction function v that maps the feature values at t to a probability forecast ( ) ∶ 2 → Ρ . If s is a coalition or subset of features ⊂ , then v(s) describes the worth of coalition s, or the total "payout" (i.e. probability) that the members of s can expect to receive by working together.
The Shapley value suggests a way to distribute the total forecast probability to the grand coalition of all features f after considering all possible coalitions and their individual worths. The Shapley value, which in the present case is the amount of "probability" apportioned to feature i in f at t, is given by: The second term inside the summation can be interpreted as the amount of payout fairly credited to feature i if i were to join coalition s. The summation is conducted over all subsets s of f not containing feature i, or rather all permutations of coalitions s that i does not belong to but could join. In effect, feature i's total payout is the average of its payouts over all those permutations.
The Shapley value is considered to be a fair attribution of the total probability in the sense that it possesses a few desirable properties. Those pertinent to the present discussion are summarized as follows: • Efficiency -The Shapley values for all features in the set f add up to the total payout v.
That is ∑ • Null Player -If a feature cannot contribute to any coalition to which it is not a member, than its Shapley value is zero. That is, if ( ∪ { }) = ( ) for every coalition ⊆ \{ } then ( ) = 0.
Though these ideas form the basis of the SHAP framework, direct implementation of Shapley values in a machine learning setting is generally not feasible, since it may require the classifier to be retrained on all the permutations of ⊆ \{ } for each feature i in the data set (as well as all observations t). In order to apply Shapley values in practice, an approximation is required. The SHAP framework unifies several local Shapley value approximation methods under a more general conditional expectation function of an abstract machine learning model. In general, each of these local approximation methods maps to a single class of machine learning algorithms. The SHAP methods relevant to the present discussion and the algorithms they map to are: • TreeSHAP -SHAP value approximation for tree ensembles (Lundberg et al, 2018). We apply TreeSHAP to the Random Forest, XGBoost and LightGBM classifiers in this study.
• DeepSHAP -SHAP value approximations for deep neural networks. This method wraps DeepLIFT (Shrikumar et al, 2017). Though we study neural network classifiers in this paper, the DeepSHAP framework will not integrate with our implementation 19 , and so we use the more generally applicable KernelSHAP for neural network SHAP value decomposition.
• KernelSHAP -SHAP value approximations for general functions. This method uses LIME (Reibero et al, 2016) to locally approximate any classifier function. We apply KernelSHAP to the probit, support-vector machine and neural network classifiers in this paper (i.e. all except the tree ensembles).

VIII. Results
We begin by presenting the results of NTS cross-validation across all models and all algorithms. As mentioned in the previous section, we have chosen accuracy to score the models. Table 14 in the appendix shows accuracy estimates for each model and classifier with standard errors, as well as the mean accuracies conditional on the model. Rows are sorted in order of decreasing mean accuracy conditional the model. The top 5 individual classifiers are highlighted in dark orange. The "best-inclass" nest of models (to be defined shortly) is highlighted in light orange. Table 14 can be summarized and more easily ascertained in Figure 6 below, which scatterplots the estimates of classifier forecast accuracy against the number of features in its model. The points are color-coded by algorithm. The dashed lines connect the means of the forecast accuracy estimates conditional on algorithm and number of features.  Table 14.

Much of the information in
The first column of Table 3 below repeats the bottom row of Table 14, which is the mean forecast accuracy, conditional on algorithm, across all 106 models, with conditional standard errors in parenthesis. The second column weights the conditional mean and standard errors calculation (Cochran 1977) to control for the varying sample size within each tier of models with the same number of features (because some algorithms perform better with a few number of features, and others may perform better on many features). The results are generally consistent between the two methods of calculation. As a first order consideration, probit performs the best, followed generally by neural networks and/or LightGBM, then XGBoost, Random Forest and/or support-vector machines.   Table 16 in the Appendix reports the means of the estimated forecast accuracies, conditional on algorithm and number of features. They are plotted below as well (and in Figure 6 as dotted lines).  Table 16. A few features stand out in the figures above and their table equivalents in the appendix: 1) Before considering the statistical significance of the results, the probit methods seem to outperform the other algorithms generally across all models, at least up through 6 features. Neural networks appear to be a close second over many models, though LightGBM is sometimes competitive. XGBoost, Random Forest and support-vector machines generally round out the bottom of the ranking, in that order.
2) The best single model generally was the 7-feature Slope/LEI/ACM5D/NFCI/EBP/SPF/SP500 model. This model also nests many of best-performing 2-, 3-, 4-, 5-and 6-feature models as well, as highlighted in light orange in Table 14. The best single classifier was the 6-feature Slope/ADS/NFCI/EBP/ACM5D/SPF probit classifier, highlighted in dark orange in Table 14. Two of the top five classifiers were probit classifiers and two were neural networks.
3) Among the 1-feature models, the forecast performance of LEI, NFCI, SPF and ADS exceeded that of Slope, in that order. Slope outperformed EBP, SP500, ACM5D and FF, which is consistent with the literature.
4) The 2-Feature Slope/LEI model is ranked 12 th in Table 14, and the 3-feature Slope/LEI/ACM5D is ranked 7 th , both higher than many models with far more features. This highlights the high degree of forecast power that Slope and LEI contain.
5) The worst performing models generally include the effective federal funds rate. 20 Even when combined with other features that together perform well, the federal funds rate appears to diminish the power of a model, particularly if estimated by probit regression.
6) Accuracy increases generally in the number of features used, up to 7 features, then declines. It seems that 4 features is better than 1, but it is not clear that 7 is necessarily better than 4, after taking into account classifier variance. This topic will be addressed in more detail shortly.
Next we compare the results of NTS cross-validation to 5-fold cross validation. Stratification is used when forming folds, due to the class imbalance in the indicator series and to ensure that the training and testing sets contain similar proportions of positive and negative instances. Furthermore, shuffling is employed to contend with the existence of structural breaks in the data. 21 Table 15 in the appendix shows the results of 5-fold cross-validation for each model and classifier, as well as the averages conditional on the model. Rows are sorted in order of decreasing mean accuracy conditional the model. The top 5 individual classifiers are highlighted in dark orange. The "best-in-class" nest of models from NTS cross-validation (which again will be defined shortly) is highlighted in light orange, for comparison.
For the purposes of comparing the two cross-validation strategies, we have included one 12feature and one 20-feature model, in addition to those described in Section IV. The 12-feature model begins with the nine features under study and adds the natural rate of interest (r-star) of Laubach and Williams (2003), GARCH(1,1) conditional volatility of West Texas Intermediate (WTI) spot crude oil prices and GARCH(1,1) conditional volatility of the Federal Reserve's Trade Weighted U.S. Dollar Index. All three of these features have very low recession forecasting power and essentially amount to noise in the algorithms. The 20-feature model further adds 1) the near-term forward term spread of Engstrom and Sharpe (2018), the 6-month percentage point change in the 10-year term premium as well as the levels of the 5-year and 10-year term premiums of Adrian, Crump and Moench (ACM, 2013), the Treasury 3-month bill discount, the real Federal Funds rate (as defined in Wright, 2006) and the default and "GZ" spread of Gilchrist and Zakrajsek (2012), all of which also either amount to noise in the algorithms or are highly correlated with some of the nine features under study. Our purpose for doing this is mainly to investigate the robustness of both the algorithms and the cross-validation methodologies to noisy or redundant (co-linear) features.
The results are displayed in Figure 8 below in two panels. The panel on the left shows the results from NTS cross validation, and like Figure 6 it scatterplots the estimates of classifier forecast accuracy against the number of features in each model, color-coding points by algorithm and showing the means of the forecast accuracy estimates conditional on algorithm and number of features as dashed lines. The panel on the right shows the results from 5-fold cross-validation in the same format.  Table 14 and Table 15. A few conclusions may be drawn from the figure: -The accuracy estimates from k-folds cross-validation are higher than the estimates from NTS cross-validation. This seems to hold no matter the criteria used to make this assessment. The highs are higher, the lows are higher, and the means and medians are higher all along the curves shown. This is consistent with the existing literature and here we have it is very clear visual form.
-The accuracy estimates from NTS cross-validation, conditional on algorithm and number of features have much higher variance than k-folds cross-validation. This calls into question whether the k-folds cross-validated classifiers, which have been allowed to "peek" at future data before forecasting a sample that is held out, are truly able to forecast recession or have merely memorized the joint time-series of recession and the features.
-Both cross-validation strategies suggest that there is value in using up to 4 or 5 features and possibly up to 7 or 8 to forecast recession. The NTS cross-validation suggests that the addition 3 to 11 more features, beyond the 9 under study, causes forecast accuracy to deteriorate for all algorithms. The k-folds cross-validation suggests that, while there is perhaps little to be gained by adding 3 to 11 more features, there is little to be lost either. With the exception of the probit method, all of the algorithms are able to handle a large number of features fairly well, according to the k-folds cross-validation.
Another salient feature of this comparison between cross-validation methods is more difficult to see in the figure above. Similar to Table 3, the first column of Table 4 below repeats the bottom row of  Table 15, which is the mean forecast accuracy, conditional on algorithm, across all 106 models, with conditional standard errors in parenthesis. The second column weights the conditional mean and standard error calculation to control for the varying sample size. Though the results are generally consistent between the two methods of calculation, they stand in stark contrast to the results in Table 3, in that they indicate rank reversal of probit, neural network and the tree ensemble classifier forecast performance over the two cross-validation methodologies. That is, the NTS cross-validation shows that probit regression is generally the best of all methods, followed by neural networks and the tree ensembles, while the k-folds cross-validation indicates the exact opposite. Under 5-fold cross validation, LightGBM, Random Forest and XGBoost perform best, in that order, followed by support vector machines, neural networks and probit regression. We take that as some evidence, however incomplete, that caution should be used when assessing the performance of machine learning algorithms vis-à-vis probit regression when data peeking or future leakage occurs during cross-validation.  We may now ask why in particular the probit method has outperformed the other methods in NTS cross-validation. Upon inspection of the outer loop scores from the nested-cross validation algorithm, it appears that the outperformance occurs broadly over all outer loop iterations. The table below averages the outer scores for each iteration in the NTS cross validation across all of the classifiers all 106 models and for each of the algorithms. The figure following shows the same information in chart format.     Table 5.

Outer Loop
In the figure above, it is clear that the outperformance of the probit method occurs broadly across iterations (and hence business cycles). In a later section we discuss some possible reasons for this result, and how it may be interpreted. For now, this concludes our presentation of results from kfolds cross-validation, and from this point forward we consider only the results of NTS cross-validation unless otherwise noted, assuming these to be the more conservative results where recession forecasting accuracy is concerned. Next we turn to matters of statistical inference, before concluding this section.
After considering all of these results, and in light of the discussion in section VI on the omnibus Cochrane's Q and pairwise McNemar's test and the need to pre-filter on our 106 models in order to make those exercises feasible, at this juncture we will appeal to intuition and limit all further analysis to what we'll call the eight "best-in-class" nest of models, mentioned above and listed in Table 6. On the grounds that these are essentially best-in-class for any given number of features, all subsequent analysis will be conducted on these models -listed below for clarity -only. The other 98 models we will consider to be either inferior or at least not statistically different from the eight chosen. Model  1  Slope  2  Slope/LEI  3  Slope/LEI/ACM5D  4 Slope/LEI/NFCI/ACM5D 5 Slope/LEI/NFCI/EBP/ACM5D 6 Slope/LEI/NFCI/EBP/ACM5D/SP500 7 Slope/LEI/NFCI/EBP/ACM5D/SP500/SPF 8 Slope/LEI/NFCI/EBP/ACM5D/SP500/SPF/FF Table 6 -"Best-in-Class" Nested Family Models. All subsequent analysis will be conducted on classifiers estimated from the models listed.

# Features
The forecast accuracy estimates from NTS cross-validation for these models are shown below in Figure 10. This figure is exactly the same as Figure 6, except that it isolates the eight models of present interest.  Table 6, estimated using nested time-series cross-validation, is scattered against the number of features in its model above. Points are color-coded by algorithm.
Finally we consider formal tests of significance of our results. To this point we have appealed to intuition or left the matter unspoken, but the omnibus Cochrane's Q and pairwise McNemar tests suggested by Raschka (2018) and described in Section VI provide a framework for incorporating classifier uncertainty into the analysis. In Table 7 below we first present Cochrane's Q for joint tests on each of the models in Table 6 across algorithms. That is, we conduct 8 omnibus tests on the test set forecasts (outer loop/out-of-sample) from NTS cross-validation, one on each of the best-in-class models model and each test containing 6 classifiers corresponding to the 6 learning algorithms under study.  Table 6 and over the six classifiers. **p<0.01, *p<0.05.

Model
The first obvious point to draw out of Table 7 is that the six classifiers for each of the 2-, 3-, 5and 8-feature models are not statistically different from one another. At almost any reasonable threshold we fail to reject the null hypothesis. This is an interesting result in light of the considerable attention that variations on the probit method have received in the literature over the past two decades, particularly those that extend the method on 2-and 3-variable models to adjust for the timeseries properties of the data. It suggests that, though we may throw a great deal technological firepower at the problem, there are limits to the amount of information or forecasting power that can be extracted from a monthly panel of 2 or 3 macro/financial series and that more data in the time-series dimension may be what's really necessary to improve forecasts. The null hypothesis for the 1-, 4-, 6, and 7-feature models is rejected at the 1% level. For these four models, pairwise McNemar tests are required to determine where the differences can be found, which we will present shortly.
In addition to running the omnibus tests for each of the 8 models over the algorithms, we also run them for each of the 6 algorithms over the 8 models. In Table 8 below we present Cochrane's Q tests on each of the algorithms, across the models in Table 6. In this case we reject the null hypothesis in all cases. This should come as little surprise however, as it is tantamount to saying that the 7-feature model always performs better than the 1-feature model, regardless of the algorithm used, which is more or less evident in Figure 6 even before testing for significance. Turning to the pairwise McNemar tests by model, Table 9 below displays four matrices of pairwise p-values, one matrix for each of the four models for which the null hypothesis was rejected in Table 7 after running joint Cochrane's Q tests. Cells highlighted in light orange indicate that the null hypothesis of no difference in error rates between the classifiers is rejected at the 5% level, dark orange indicates rejection at the 1% threshold. A few points can be drawn from the table:

Algorithm
-The results for the 1-feature Slope model suggest that that probit and XGBoost classifiers outperform the other four algorithms. Furthermore, the probit and XGBoost are not distinguishable from each other and, likewise, the other four algorithms are not distinguishable from each other. This is consistent with the results shown in Figure 10.
-The results for the 4-feature models suggests that LightGBM underperformed generally, but the other algorithms are indistinguishable, also consistent with casual inspection of Figure 10.
-The results for the 6-feature model suggests that the neural network outperformed all other algorithms, and that the probit method underperformed most, also consistent with casual inspection of Figure 10.
-The results for the 7-feature model suggests that the neural network and support-vector machine outperformed XGBoost and LightGBM, but the Random Forest and probit were indistinguishable from all, consistent with Figure 10.   Table 7 are shown. **p<0.01, *p<0.05. Table 17 in the appendix displays six matrices of pairwise p-values, one matrix for each of the algorithms for which the null hypothesis was rejected in Table 8 after running joint Cochrane's Q tests. Cells highlighted in dark orange indicate that the null hypothesis of no difference in error rates between the classifiers is rejected at the 1% level, light orange for 5%. The results for LightGBM are representative of the general pattern in the data, which is that, within a given classifier, the 1-feature model is statistically different than all other models. The differences stop there however and we fail to reject the null hypothesis of no difference between classifiers with 2 or more features. This is consistent the results shown in Figure 10, which shows great forecast accuracy improvement if LEI is added to the Slope model, but no consistent marginal benefit to adding more features. 22 This again highlights the explanatory power contained in just these two features, and the limitations posed by the sparse dataset.
This concludes the discussion on tests of significance. In closing this section a few feature importances are presented. Recall that, having estimated the forecast accuracies of the classifiers, we will dispense with the NTS cross-validation from this point forward use classifiers that have been trained on all of the data, so all results are in-sample. 23 Probit sensitivities are shown in Table 10 below. As would be expected, the Slope coefficient is strongly significant, and the marginal effect is consistent with previously published estimates, as is LEI. NFCI is also very significant, and absorbs much of the explanatory power of EBP, which is a new result. ACM5D is modestly significant while SP500, SPF and FF add little to the regression after considering the other features. The R 2 of the regressions grow rapidly as features are added to the univariate model, but beyond five features there less to be gained, consistent with many of our previous results. The next table displays the permutation importances of the features for each classifier for only the 7-feature "best-in-class" model, with standard errors in parenthesis. Notice that not all algorithms rank Slope highest along this measure, including the probit method. In general, however, the rank ordering is similar to that of the probit sensitivities for this model. NFCI is ranked generally higher than EBP, and in some cases NFCI is ranked 2 nd . At the bottom of the list the permutation importances are generally low and there is noticeably less agreement and statistical significance among the results. For purposes of comparison to the global feature importance metrics just presented, the mean absolute SHAP values for the 7-feature best-in-class model are shown below for each classifier. After decomposing the classifier outputs at the level of the observation, the metrics are aggregated up to a global level to get a sense of feature importance across time. Note that the values for LightGBM and XGBoost are in units of log-odds, while the others are in percent. As with the MDA metric, LEI is ranked highly and often above the Slope measure. NFCI largely outranks EBP, which outranks ACM5D. Overall the results are consistent with the probit sensitivities and MDAs. In the next section we will discuss the decomposition of the SHAP values over observations in more depth. To conclude this section, our results show that -across all 106 models studied -probit classifiers outperform machine learning classifiers for the problem at hand when NTS cross-validation (which precluded data peeking) is used. When considering only the eight "best-in-class" subset of models, the probit method stands out, at least when four or fewer features are used, as shown in Figure 10, statistical significance notwithstanding. This stands in contrast to a body of research documenting the superior performance of the machine learning methods investigated here to other algorithms, including probit regression (Fernandez-Delgado et al. [2014], Wainer [2016]). That said, the differences among classifier performances across algorithms for the 2-, 3-, 5-and 8-feature best-in-class models are not statistically significant. Furthermore, for all of the classifier algorithms tested, there is little to distinguish the 2-through 8-feature models from one another, and it appears that there is not much to be gained by including more features in the models, after using Slope and LEI.

Features/Algorithm
We believe this mixed, but generally consistent evidence that machine learning classifiers underperform the probit method to be consequence of the relatively small, sparse data that has been used. As described in Section VI, the data used, though standard in the literature and the best available for the task at hand, is of monthly frequency and extends back to only 1972. Furthermore, the data is serially correlated and contains multiple structural breaks, and the binary indicator we are attempting to forecast is unbalanced and occurs temporally in clusters. The machine methods that we are using, however, have been developed and optimized for application on much larger data sets and much of the existing literature on their use presumes independence and identical distribution of data.
The challenges posed by the dataset are likely to be common to many applications of machine learning methods to macroeconometric problems. Though machine learning methods and their power against large datasets have advanced greatly in recent years, time-series macroeconomic data is still accrued at the same rate it always has been: slowly over the course of business cycles and generations. And even if the collection of macroeconomic data too is accelerating in many ways, our ability to backfill newer sets with data from even ten years ago, let alone the 20 th century is questionable at best. In the next section we discuss the implications of the small size and sparseness of the data set vis-à-vis machine learning methods and what can be gleaned from these algorithms about recession forecasting (inferior performance notwithstanding). Even if the machine learning classifiers do not outperform the probit method on the basis of forecast accuracy, there is still much to be learned by attempting to interpret their output.

IX. Interpretation
In the analysis that follows, we present some salient details of our results. Again we limit the discussion to just the eight models listed in Table 6. As in the presentation on feature importance, we use classifiers that have been trained on all of the data, so all results are in-sample.

1-Feature Slope Model
Beginning with the 1-feature Slope model, Figure 11 shows the model-implied probability of the indicator as a function of the Slope variable for each of the algorithms. Going forward we will refer to this as the decision function, for the sake of simplicity. Recall that, after application of the McNemar test, it was shown that the XGBoost and probit algorithms outperformed the other algorithms for this model. Note that all algorithms place the 50% decision threshold at nearly the same place in the narrow range of 2bps to 36 bps and only 31 of 558 observations (5.6% of the sample) fall in this range.
There are some noticeable differences between the classifiers in Figure 11. First, the machine learning classifiers closely track each other across the Slope domain, but differ from the probit classifier in important ways. As Slope falls from high values towards 1 percentage point (pp), the machine learning probabilities remain stable at less than 20% generally. As Slope crosses the 1 pp threshold and continues to fall further, however, the machine learning probabilities rise quickly over a range of 1.5 pps, reaching about 85% at a Slope of -0.5 pp. (The probability of the indicator in the data, conditional on Slope < -0.5 pps is 100%, as can be seen in the scatter of red and black dots in Figure 11 to the left of that level.) In contrast, the probit probabilities rise gradually over the full range of Slope displayed, in a manner determined strictly by the probit link function. Additionally, the machine learning probabilities generally do not fall below 10%, while the probit probabilities approach zero asymptotically for high values of Slope. Given that the probability of the indicator in the data conditional on Slope > 1.5 pps is 11.5%, it would appear that the machine learning estimates and in particular the ensemble tree probabilities are about right.
These results suggest that the flexibility of the machine learning methods allows them to capture the empirical distribution of the recession indicator in ways that the probit classifier cannot. Finally note that the neural network classifier combines the best qualities of all methods: while capturing the non-linear nature of the decision function in ways that only the machine learning methods can, the neural network decision boundary is a smooth curve like the probit decision function with asymptotic properties more appealing than the support-vector machine, and places the 50% threshold at nearly the same point in the Slope domain. The policy implications of all these facts, when attempting to judge recession probabilities from probit regression should be clear. In the introduction, we pointed out that, at time of writing, the spread between the 3-month Treasury bill discount and the Treasury 10-year note yield-to-maturity was 11 bps. Though our data set ends in June 2019 and we do not show forecast values past that date in Figure 1, at time of writing the value of the term spread used here and in Figure 11 is about 14bps. According to the machine learning classifiers (which together forecast recession in the following 12months with 58% probability at this level of Slope), a further decrease of just 50bps (the equivalent of two hikes of the short term interest rate) would signal a coming recession with over 85% probability.
The probit method on the other hand (which places the probability of recession at 55% at the current level of Slope) suggests there remains much more room for term spreads to run (about 161 bps or roughly six interest rate hikes) before that high level of probability is reached. To the extent that policymakers use a 1-feature Slope model for the purpose of forecasting recessions, attention to this result is warranted in our view.

2-Feature Slope/LEI Model
Next we take up the 2-feature Slope/LEI model. Recall that, though the probit classifier outperformed during NTS cross validation, the classifiers were not found to be statistically different across algorithms for this model. Also recall the Slope/LEI model performed very highly among the 126 models in the overall study, and better than many 3-, 4-and 5-feature models, suggesting this pair of features contains a high degree of explanatory power.
In Figure 12 the Slope and LEI data is shown in a scatter plot, along with the probit boundaries for 25%, 50% and 75% probability. The points in the scatter are red if the indicator variable is true and black otherwise and all points are labeled with the year in which they were observed. Data from the Volcker and pre-Volcker recessions mainly concentrate to the left of the LEI axis, and data from the Great Recession mainly concentrates in the lower right of the fourth quadrant, with the 1991 and 2001 recessions in between. Probability of a recession in the following 12-months as a function of 3m/10y slope and 3-month log-difference of the Conference Board's Leading Economic Index is displayed above. Indicators code the occurrence of an NBER recession in the next 12 months in the sample data and are shaded black (0) and red (1).
The results from the probit classifier are relatively well documented in the literature. Before moving on to the machine learning classifiers, we will point out a few features in the plot. First we might question whether a linear separator, regardless of whether it is estimated by probit regression, logistic regression, support-vector machines or any other method is the appropriate separator for this data. The data is sparse, and while admittedly a linear separator performs well near the center of the mass around the 1991 and 2001 recessions, it is less clear that a linear separator is appropriate in the regions where Slope is negative or LEI below 1% or so. Certainly there is no theoretical justification for extending the probit line to infinity in both directions. The justification for doing so at this point is merely a practical one, since our sample includes very few realizations in the tail regions of the joint distribution.
As an alternative to a linear separator, we posit there exists a critical level of Slope below which the level of LEI ceases to matter. If that were so, we would consider the recession probability everywhere in the left half-plane defined by that critical Slope level to be very high. Conversely, if macro conditions are severe enough -as in the Great Recession -it could be argued that there ceases to be any level of Slope steep enough (or perhaps any monetary policy accommodative enough) to avert recession. As such, we would consider the probability of recession everywhere in the lower half-plane defined by that critical LEI level to be very high as well.
Next we move to the Random Forest decision function in Figure 13 to consider what the ensemble tree methods have to say about these matters. Here we have plotted the decision function in detail, color coding the model probabilities. Red shading denotes degree of probability above 50% and gray shading denotes the degree of probability below 50%. The white region separating the red and gray is the region of 50% probability, and constitutes the decision boundary. The decision function has been overlaid with the contents of Figure 12 for comparison. Probability of a recession in the following 12months as a function of 3m/10y slope and 3-month log-difference of the Conference Board's Leading Economic Index is displayed above. Indicators code the occurrence of an NBER recession in the next 12 months in the sample data and are shaded black (0) and red (1). Results from a probit classifier are overlaid for comparison.
In the 1-feature Slope model, there was little to distinguish the decision boundary of the probit classifier from that of the ensemble methods; for all methods the boundary existed as a point in the feature space. Here, in two dimensions, the distinction becomes more obvious. While the probit method linearly separates the data by drawing a line (hyperplane) through the feature space, the decision boundary for the Random Forest classifier (and all tree ensembles) is a piece-wise continuous collection of axis-aligned line segments always intersecting at right angles to one another. In the limit that we allow the Random Forest algorithm to grow an infinity of trees with unlimited depth on a sample closely approximating the population, the decision boundary will asymptotically approach a smooth curve (or even a line if need be). Our sparse sample and the bias-variance tradeoff however have forced us to limit the depth of the trees, and computing constraints force us to grow a finite number of them. This in turn limits the number of line segments available to create a decision boundary and makes the discrete nature of its estimated form more observable. Figure 14 shows the decision function for the XGBoost algorithm. The power and propensity of these methods to overfit data becomes clear in this plot. That said, forecast variance seems to have been controlled enough during cross validation to obtain results that are not statistically different from any of the other Slope/LEI classifiers. Probability of a recession in the following 12-months as a function of 3m/10y slope and 3-month log-difference of the Conference Board's Leading Economic Index is displayed above. Indicators code the occurrence of an NBER recession in the next 12 months in the sample data and are shaded black (0) and red (1). Results from a probit classifier are overlaid for comparison. Figure 15 shows the decision function for the LightGBM algorithm. The decision boundaries, particularly that for LEI is very consistent across the three tree ensemble classifiers presented so far. Probability of a recession in the following 12-months as a function of 3m/10y slope and 3-month log-difference of the Conference Board's Leading Economic Index is displayed above. Indicators code the occurrence of an NBER recession in the next 12 months in the sample data and are shaded black (0) and red (1). Results from a probit classifier are overlaid for comparison. Figure 16 shows the decision function for the neural network algorithm. As in the 1-feature Slope model, this algorithm appears to combine the best qualities of all models. On the one hand it is a non-linear separator and captures the underlying phenomenon in ways that the probit method cannot, and on the other it has managed the bias-variance tradeoff well and produced a smooth decision function that is easy to interpret vis-à-vis the tree ensemble methods. Probability of a recession in the following 12months as a function of 3m/10y slope and 3-month log-difference of the Conference Board's Leading Economic Index is displayed above. Indicators code the occurrence of an NBER recession in the next 12 months in the sample data and are shaded black (0) and red (1). Results from a probit classifier are overlaid for comparison. Finally, Figure 17 shows the decision function for the support-vector machine algorithm. Note that a 2 nd order polynomial was chosen as the kernel function in cross-validation, and underlies these results. As in the 1-feature Slope model, its main flaw appears to be its asymptotic properties. In short, while the decision boundary for this algorithm displays many of the qualities of the neural network classifier, it retains a degree of rigidity similar to the probit method that (in conjunction with the requirement that the researcher choose a kernel function for this algorithm) is difficult to manage, particularly in a high number of dimensions or on sparse data. Probability of a recession in the following 12-months as a function of 3m/10y slope and 3-month log-difference of the Conference Board's Leading Economic Index is displayed above. Indicators code the occurrence of an NBER recession in the next 12 months in the sample data and are shaded black (0) and red (1). Results from a probit classifier are overlaid for comparison.
Where policy considerations are concerned, at time of writing and using the latest data available, the LEI (3-month log-difference) is -0.1%, placing the most recent observation to the right of the origin at (-0.1%, 11bps). Together the algorithms forecast an average probability of recession in the next 12 months of about 36% (standard error 3.3%), which is about 22 percentage points lower than the 1-feature Slope model forecast average across algorithms (58%, standard error 3.6%). Recall that the pairwise McNemar tests across models for each algorithm suggested that there exists a statistically meaningful difference between these 1-and 2-feature models. This highlights the value to including a macroeconomic indicator alongside term spreads in a recession forecasting model, regardless of which classifier algorithm is used.

5-Feature Slope/LEI/ACM5D/NFCI/EBP Model
For models with more than two features it becomes more difficult to interpret classifier output. Even in three dimensions it can be difficult to visualize decision boundaries as done above for the 1-and 2-feature models. In this section we show how the SHAP value framework can be used to assist in classifier interpretation, as well as to decompose recession forecasts among features in a model. The 5feature Slope/LEI/NFCI/EBP/ACM5D model is used for this part of the discussion, since it was not statistically different from the best-performing 7-feature model across classifiers but is more parsimonious.
In Figure 18 below, the SHAP values for the Slope feature are scatter-plotted against the feature values for each algorithm. High-order polynomials are fit to the scatter (except for XGBoost) to make the direction of the marginal effect more clear. Note that XGBoost and LightGBM SHAP values are in units of log odds (right-hand side) while the others are in units of probability. The direction of interaction is consistent for all algorithms (downward slope) and the probit and Random Forest classifiers are less sensitive to Slope (shallower slope) than the neural network and SVM classifiers. It is difficult to compare LightGBM and XGBoost directly to the other classifiers due to the different units in which they are expressed, but generally the scatters for these two classifiers are similar. All of these results are consistent with the permutation importances and mean absolute SHAP values report in Table  11 and Table 12. A similar plot for the ACM5D feature in the 5-feature Slope/LEI/NFCI/EBP/ACM5D model is shown below in Figure 19. While this feature has a lower mean absolute SHAP value than the Slope feature (as evidenced in Table 12 and confirmed below by the lower range displayed for the y-axis), all algorithms identify the same direction of interaction (upward slope). Though the ACM5D feature lacks the importance of Slope, LEI and NFCI in this model, this result indicates that rising term premiums are a feature of impending recession.    Using the SHAP value framework, we decompose these forecasts to attribute the predictions (and discrepancies) to specific features. Figure 22 below displays the SHAP values for the probit classifier on the 5-feature model as a stacked bar chart, and Figure 23 shows the same data for the neural network classifier. When SHAP values for an observation are summed with its expected (unconditional) model forecast (~25%) the results equal the classifier-implied recession probability for the observation, as displayed in Figure 20 and Figure 21.
Recall that both of these classifiers performed similarly in NTS cross-validation. Note the correspondence below between their SHAP decompositions across the sample. At the end of the sample, the recession signal emanating from the term spread (Slope) signal is canceled by an even mix of low term premiums (ACM5D), easy credit (EBP) and financial (NFCI) conditions. The neural network appears to view macroeconomic conditions (LEI) more favorably than the probit classifier, however.

Discussion Wrap-up
As mentioned earlier in the introduction, the flexibility of machine learning classifiers is both an attractive feature of the methodologies and a source of new difficulty. Managing the bias-variance tradeoff, particularly on a dataset as small and sparse as that used here, requires experience with the hyperparameters exposed by each individual algorithm. There is little science to the process of gridsearching over a hyperparameter space and, when confronted with constraints of time and finite computing resources, subjective judgment is required to identify the subset of hyperparameters and values to be searched over.
In our discussion of the 1-and 2-feature and 5-feature models above, we were able to show that results for the machine learning methods track closely and that the decision functions and SHAP values are well-behaved and appeal to economic intuition. We take this to suggest that our search for "optimal" hyperparameters was more or less successful and we have arrived at some kernel of truth in the data, despite their inferior performance in NTS cross validation.
The dataset used in this analysis is sparse, and it is doubtful any method is able to estimate recession probabilities in the tails of the joint feature space well. This fact needs to be considered when assessing any of the figures discussed here. It also highlights the importance of using visualization to supplement the more technical results of the previous section when assessing performance of any binary classifier.
Setting aside these caveats, we might still ask what can be learned from the present discussion that could not be learned from previous studies using probit methods. We summarize the answer to that question as follows: -In the discussion of the 1-feature Slope models it was shown that the machine learning decision functions are generally steeper and more "non-linear" than the probit decision function. Although perhaps not shown as clearly, this idea can be extended to the 2-and higher-feature models as well. This result should be considered when setting policy.
-Adding a leading indicator of macroeconomic conditions (such as the Conference Board's Leading Economic Indicators, or the Philadelphia Fed's ADS Index) to term spreads greatly improves forecast power at little cost, regardless of the classifier algorithm used for the purpose.
-The Chicago Fed's NFCI appears to capture some aspects of financial conditions that the Excess Bond Premium does not.
-Rising term premiums may be a subtle forecaster of recession.
We close this discussion with a few comments on the use of machine learning methods in macroeconometrics more generally. Many of the lessons learned here are likely applicable to other studies using the same or similar data.
As we mentioned earlier, the machine learning methods we've used are perhaps more appropriately applied to datasets larger than our panel of monthly time-series data. An online retailer deploying these algorithms to build recommender systems, or a credit card company using them to score applicants will collect every minute if not every second new data equivalent in size to the entire 50-year macro/financial panel we've used in this paper. It is our opinion that studies employing machine learning methods on new sources of "big" macro and micro data for which most machine learning algorithms are optimized represent an exciting new frontier in econometrics. We also believe that the vast majority of those studies using supervised classification algorithms on "big" data will find that machine learning methods do outperform probit methods along many dimensions, consistent with the growing literature on the subject.
For macroeconometric studies like ours, which employ more traditional monthly and quarterly macro and financial time-series data to (re-)explore parts of the established literature, we believe the challenges posed by the small size and sparse nature of those sets are considerable. In contrast to the large data sets employed by the online retailers and credit card companies of the world, which closely approximate a population and for which standard i.i.d. assumptions often hold, a 50-year time-series of monthly financial market and macroeconomic observations will exhibit not only serial correlation but also multiple structural breaks and, if categorical features like recession are involved, class imbalance as well. It is not clear (to the authors at least) that current technology offers all the tools necessary to contend with these problems.
Complicating our situation further is the fact that the major cross-validation methodologies (such as k-folds, leave-p-out, etc.) commonly used to estimate models using machine learning algorithms and manage the bias-variance tradeoff are also built for application on i.i.d. data, rather than serially correlated time-series data. Though expanding and sliding windows cross-validation are available for the latter type of data, and are likely work well on time-series of daily or higher frequency and many years or decades in length, in cases such as ours these methods presents difficulties.
All this notwithstanding, in managing the bias-variance tradeoff, most of our lessons learned so far are encapsulated in the cross-validation methodology we've employed. After many attempts with other strategies, the NTS cross-validation strategy that we have arrived at achieves, we believe, a good balance between training/test set size and validation set size given the constraints posed by serial correlation, class imbalance and structural breaks in the data. It is also sure to provide a very conservative estimate of forecast performance. While allowing hyperparameters to be chosen and outof-sample forecast error to be estimated simultaneously, the NTS cross-validation strategy: -Avoids the creation of a test set, to be held out entirely from training for the purpose of outof-sample forecast error estimation. Holding out the roughly 20-30% of sample data that is commonly needed for this purpose is very costly on small datasets and can result in much higher forecast variance if the time-series available for training is made 20-30% shorter. In cases such as ours, when a positive instance of the indicator does not exist in the most recent 20% of the data, it is perhaps fatal to the study.
-Mitigates class imbalance and problems that arise when sampling unbalanced categorical data, through the use of stratification. Some business cycles are long and some business cycles are short, and likewise some recessions are long and some recessions are short. Stratifying by business cycles ensures that all recession and non-recession periods in the data are sampled from in a balance manner.
-Mitigates some of the problems posed by sampling without replacement from data containing structural breaks, by chaining forward over entire business cycles, rather than data blocks of fixed size. If data blocks of fixed size are used, individual folds of data may end up containing zero instances of a positive indicator, or may bridge business cycles, structural breaks and recession periods awkwardly.
-Simplifies model diagnosis since the samples of forecast error (the outer loop scores) are mapped to specific business cycles. If the forecast error estimate of an outer loop iteration is poor or unusual, the cause is more readily linked back to the idiosyncrasies of the business cycle it is attempting to forecast.
-Avoids data peeking or future leakage, since it is always the case that forecast error is estimated from a test set consisting of observations occurring strictly after all those used to train the classifier.

X. Conclusion
Machine learning methods have gained widespread use and acceptance in binary classification problems outside of finance and economics. This paper has investigated their use in forecasting US recessions from term spreads and other financial market and macroeconomic data. To date, most studies examining the problem at hand have used probit methods exclusively.
It has been documented by others in previous work that Treasury term spreads and many of the other financial market and macroeconomic data used in this study (the Excess Bond Premium of Gilchrist and Zakrajsek, the Conference Board's Leading Economic Index and the Survey of Professional Forecasters, etc.) are forecasters of recession in the US. In this study we show that the Aruoba-Diebold-Scotti (ADS) business conditions index, 6-month changes in ACM 5-year term premiums and the Chicago Fed's NFCI index also have power to predict US recession.
Although the machine learning methods generally underperform probit regression in NTS cross validation, there is still something to be learned from them. We have shown that machine learning methods are able to capture important features of the joint empirical distribution of Treasury yields and other financial market and macroeconomic data over a recession indicator that probit methods cannot. In particular, machine learning methods, due to their flexibility, are able to capture the "non-linear" nature of those empirical distributions.
We have noted that the data set we have used is very sparse. This presents many challenges in navigating the bias-variance tradeoff when estimating models via ensemble methods. The algorithms we used are designed for application on much larger data sets and even then are prone to over-fit data. Still, in the discussion of our results, we were able to show visually that the results of the machine learning methods (the decision functions) exhibit many of the features of probability distributions and, for practical if not theoretical purposes, can be considered as such. Furthermore, the hyperparameters determined in grid search cross-validation have enabled us to produce results that generally conform to financial and economic intuition. As term spreads flatten, as the credit and financial conditions tighten or macro conditions deteriorate… recession probabilities rise across all of our models and methods. We also showed that the machine learning methods are able to identify characteristics of the empirical distribution of our features over recession that the probit method cannot.
In closing, we offer this analysis (technical, visual and qualitative) to suggest that machine learning methods be included among the standard tools now used to forecast US recession from financial and macroeconomic data. Definitions for all models used in the paper, named according to the features they contain. Each column collects models with an equal number of features (i.e. the first column contains all 1-feature models). Each of these models is estimated by all six classification algorithms studied in the paper.