Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
66,947 characters · 9 sections · 124 citation commands
Predicting Default Probabilities for Stress Tests: A Comparison of Models
Stress tests have a twenty-year history as tools for micro- and macroprudential supervision and are now used regularly by financial institutions and those who supervise them. The aim of such tests is to assess institutions' resilience to adverse financial and economic developments, as well as to contribute to the overall assessment of systemic risk in the financial system. In order to assess the behavior of financial institutions to stress, a multi-year macrofinancial scenario needs to be designed. It captures relevant systemic risks and the materialization thereof to generate stress in the system.
The present study focuses on the use of macrofinancial linkage models to translate the country-level scenario into bank-level risk parameters, which are an essential input to every stress test. These so-called satellite models thus provide scenario-conditional forecasts for the probabilities of default (PD). What is key in these models is the proper detection of signals from an array of global variables which are not directly linked to the financial health of companies.
Outside the context of stress testing, such credit risk models have a long-standing history keeton:1987,wilson:1998 and still represent a very active research area in which various models in different setups have been proposed. Specifically, the literature covers linear models (see, e.g., aver:2008; bofondi:2011), vector autoregressive (VAR) models (see, e.g., gambera:2000; pesaran:2006; castren2010), panel models (see, e.g., pesola:2001; castro:2013), latent factor models (see, e.g., koopman:2005; kerbl:2011), quantile regressions schechtman:2012 and machine learning models jacobs:2018. Interestingly, besides jacobs:2018, in the machine learning literature, there seems to be little to no coverage of credit risk models that estimate or predict PDs. Even in very recent literature surveys that specifically cover specifically machine learning in banking risk management leo:2019 and machine learning in credit risk breeden:2020, there seems to be just a few papers in the approximate vicinity of this topic.
The current state-of-the art satellite model for PD translation is Bayesian model averaging (BMA) raftery:1995. It has a long track record as being a reliable tool for generating scenario-conditional projections for credit risk and is being adopted by more and more central banks and institutions. The inherent advantage of BMA is the explicit tackling of model uncertainty by operating on a large pool of competing models which are weighted by their predictive performance and combined to one final model. However, with easier access to sophisticated regression approaches provided by open source programming languages such as R R:2020, Julia bezanson:2017 or Python vanrossum:2009 and the advent of new predictive models in the field of machine learning, the question arises if there are other models that could deliver better results.
The aim of this paper is to conduct a systematic forecast comparison with a large number of different regression models to find the best performing credit risk satellite model. The winning model is evaluated for the ability to precisely forecast default probabilities conditional on a standardized set of macroeconomic variables as provided to financial institutions by the esrb:2020 for the EU-wide banking sector stress test. We implemented a total of 43 models, which can be assigned to 9 categories, ranging from conventional statistical models to more recent machine learning methods. We tried to encompass as many models as possible with a proven track record in forecasting linear and non-linear relationships and which were readily available within an R library R:2020. Additionally, we also combine the models with different forecast combination approaches to further gauge their potential accuracy. For the purpose of this paper, we implement a framework that allows us to conduct this comparison with a standardized data set for all models, to tune the respective hyperparameters for each model and to cross-validate the results based on recursive pseudo out-of-sample forecasts.
There are only two papers in the literature that are related to our work. papadopoulos:2016 created a composite satellite model for stress testing by weighting candidate models from the full space of all possible variable combinations. grundke:2019 used a combination of BMA to select relevant variables and OLS to regress the selected variables onto a credit default index. They further analyzed different modifications to various steps in their estimation framework and assessed the different outcomes in terms of out-of-sample default rate forecasts. Both papers focus on the proper identification of the model given one estimation technique but lack comparisons across different procedures.
Our paper contributes to the literature in the following ways: To the best of our knowledge, it represents the first systematic model comparison in the field of credit risk and stress testing. First and foremost, we deliver insights on a large number of potential credit risk satellite models across a wide range of modelling techniques. Since it is not evident which modelling technique will achieve the best results, we take a naive approach by testing many different models without prejudice on how well they will fare. Our results indicate that there are better performing models than the current state-of-the art model (i.e. BMA), which have not been mentioned in the literature until now. Second, due to the large variation in models, we are able to shed light on the potential positive effects of machine learning models and thus extend the scarce literature in this area. Moreover, the use of different forecast combination techniques across different sets of models shows that not only can these techniques help to hedge against model uncertainty, but they can also enhance the predictive accuracy.
The remainder of the paper is structured as follows. In Section (ref) our estimation and evaluation framework are explained in detail, including the used data, models, hyperparameter tuning and performance measures. Thereafter, the results are shown in Section (ref), which also presents a deep dive into the winning model. Section (ref) concludes the paper.
This section outlines the setup of our forecasting exercise. First, we present the underlying data set and the applied transformations. Second, we give a short summary of the 43 models within the 9 overarching categories. Third, as many models need a prior setting of parameters, we discuss the process of tuning the respective hyperparameters. Fourth, we briefly discuss the measure used to evaluate the forecasting performance.
The data involved in the modelling exercise refer to Austria and include as dependent variable a measure for the probability of default for the non-financial corporate sector and macroeconomic and financial data as independent variables. The deployed default probabilities are based on Moody’s KMV Expected Default Frequency (EDF) measure, available at quarterly frequency from 2002Q2 to 2019Q2 ($T = 1,\dots,69$). The EDFs have been successfully deployed in multiple credit risk models (see, e.g. alves:2005; castren2010; gross:2019). In adherence to the IMF's approach of credit risk modelling in their Financial Sector Assessment Programs (FSAP), we also take the provided mean PD measures and not the median (see imf:2020). Hence, the measure relates to the average default probability across non-financial corporations.
The independent covariates are based on the variables within the macroeconomic scenario as designed for the EU-wide stress tests by the European Systemic Risk Board esrb:2020. These variables are real GDP growth (GDP), the unemployment rate (UNE), the inflation rate (INF), real estate price growth (RRE), stock price growth (EQP), exchange rates (EXR) and short-term (\textit{STR}) and long-term interest rates (\textit{LTR}). However, the actual scenario figures are not needed for the estimation as we conduct pseudo out-of-sample forecasts for which the actual default probability time series is needed. Hence, the scenario serves as guidance for the choice of covariates, but the models are evaluated on the basis of historical figures.
Although the range of the time series is limited by the availability of the default probability, it still includes non-linear events such as the financial crisis of 2008 and the European sovereign debt crisis of 2011, which is an important feature for credit risk satellite models, as they must be able to estimate and predict such structural breaks. We get data on real GDP, Harmonized Index of Consumer Prices (HICP) and the unemployment rate from Eurostat. The GDP figure is seasonally and calendar adjusted and transformed to year-on-year (YoY) growth rates to fit the figures used by the ESRB. Similarly, the HICP is also transformed to YoY growth rates to match the ECB definition of inflation ecb:HICP. The real estate prices and EUR/USD exchange rate are sourced from the ECB's Statistical Data Warehouse (SDW). The house prices cover all new and existing residential properties across all dwelling types and are also transformed into YoY growth rates. Finally, we take 10-year government bond yields as long-term interest rates, 3-month Euribor as short-term interest rates and the equity price index for Austria, the ATX, from the OECD database. The equity price index is also transformed into YoY growth rates.
A first descriptive analysis of the variables and potential correlations among them can be drawn from Figure (ref). The chart depicts the variables without further transformations and the grey shaded areas mark the two previous crisis periods -- the GFC (2008 Q2 - 2009 Q2) and the European sovereign debt crisis (2011 Q1 - 2013 Q1). The sample starts exactly at the peak of the crisis that unfolded in the aftermath of the burst of the dotcom bubble in 2000, the uncertainty triggered by the 9/11 attacks and the very volatile years after the introduction of the euro as a new currency from 2000 - 2001. The economic shockwave led to a significant increase in expected and actual corporate defaults around the globe altman:2003.
In the case of Austria, the default probability starts with a value of 5.3% and decmidrules until the start of the financial crisis, when it reaches the sample maximum of 6%. The pattern during the European sovereign debt crisis is not as clear, even though there is a small increase in the EDFs towards the end of 2012. The movements of the macro and financial variables behave as expected during the downturns. In both periods, we see significant drops in GDP growth and inflation, an uptick in the unemployment rate, large negative distortions on the stock market, devaluations of the euro vis-a-vis the US dollar and an increase in real estate prices reflecting a flight to safe investments, and we also see the ECB reacting to these events by pushing the short-term interest rate and indirectly the long-term interest rates towards zero and beyond. Even though the co-movements of the variables within the structural breaks seem to be going in sensible directions (e.g. GDP down, PD up), the credit risk satellite models in our setup will need to be able to pick up correct signals in calmer periods. In the majority of European countries, the same co-movement patterns have been observed in the last two decades. Therefore, our conclusions based on Austria can be generalized to other regions.
As a last transformative step, we need to make sure that the variables have optimal properties and behave well in predictions. First, to ensure that the point forecasts of the default probability is bounded between a 0% and 100% interval, we apply the following logit transformation for the regression and calculation of the performance measures,
Second, we analyze the trend and cyclical components of the variables and perform a series of unit root tests. The seasonal decomposition by loess robert:1990 shows that all variables, except GDP (since it has already been adjusted), exhibit a form of cyclicality, which we remove in due course. In order to get a clear picture of the stationarity of the variables we deploy the augmented Dickey-Fuller (ADF) test dickey:1979, the Elliott, Rothenberg & Stock (ERS) test elliott:1996, the Phillips-Perron (PP) test phillips:1988 and the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test kwiatkowski:1992. Without further transformations, all variables would suffer from unit roots. Hence, by taking the first difference for each variable, the test statistics in each unit root test indicate robust level and trend stationarity.
In this section we give a short overview of, and introduction to, the adopted credit risk satellite models. In order not to artificially expand the paper, the formal definitions of the models are not mentioned below. The interested reader is referred to the publications cited in each paragraph. The selection of the models was based, on the one hand, on the desire to cover as many models as possible with a wide array of different features. In doing so, we want to extend the existing model space in the literature from mainly linear models to non-linear, data-driven models with a special focus on regularization. Especially the latter point will turn out to be very important when analyzing the results. On the other hand, as the implementation of so many models is time-intensive, we focus our attention on proven models that are readily available, using open source computing environments, such as R R:2020.\footnote{For a compact overview of the models and R libraries see Table (ref) in the annex.}
In total, we have implemented 43 models which are placed in 9 overarching groups to give the reader a better overview of the models at hand. All satellites use the same basic model structure in which the dependent variable $\mathbf{y} = (y_1,\dots,y_T)'$ is described as a function of contemporaneous and lagged predictors $\mathbf{X} = (x_1,\dots,x_T)'$, such that
The actual approximation of $f()$ is dependent on the respective model in our comparison. Notice that the equation is ordered by the lags, i.e. $p = 0, 1,\dots P$. This fact is important for certain stepwise regressions and it achieves overall better prediction results than with the equation being sorted by the variables. We set $P=4$ as a standard choice for quarterly data. We follow gross:2019 by forcing a "closed" lag structure without gaps between the first and fourth lag for the initial equation. gross:2019 found in their analysis that this type of structure appeared to be more meaningful and robust. However, models with regularization or variable selection are not constrained in their choice of the proper equation.
We do deviate from some papers in the literature on credit risk satellite models in that we do not include an autoregressive lag of the dependent variable. When we tested both options, nearly all models showed a worse performance with the autoregressive lag than without, while the overall ranking across the models remained stable. Moreover, we did not want to increase the already large set of contemporaneous and lagged regressors ($n = 40$) in combination with the limited sample length ($T = 64$).
(Generalized) Linear Models
The first class of competitors are (generalized) linear models with a proven track record in a number of disciplines. We start with naive benchmarks based on a standard ordinary least squares (OLS) regression (short name: lm) and its robust alternative based on the M-estimator (rlm) as defined in huber:1992.
Nevertheless, due to the high number of variables relative to the sample, overfitting in the context of predictions is problematic. Thus, the next group of models uses regularization techniques to reduce the number of covariates. More precisely, we implement a forward selection algorithm (lmfs) hocking:1976 and a least-angle regression (lars) by efron:2004.\footnote{The corresponding backward selection procedure yields exactly the same best model subset. Using the Akaike Information Criterion (AIC) instead of the BIC, diminishes the out-sample performance as too many predictors are chosen as being relevant.}
Shrinkage estimators deviate from the stepwise approach as they penalize the coefficients to reduce multicollinearity in the equation. Here we chose the ridge estimator (ridge) (tikhonov:1943; phillips:1962) and, as the ridge cannot yield sparse equations because no covariates are dropped due to the regularization procedure, the least absolute shrinkage and selection operator (lasso) by tibshirani:1996.
We further implement three refinements to the lasso estimator that try to deal with certain shortcomings of the original estimator. First, if one takes the temporal structure of the data into account, we get the fused lasso (flasso) by tibshirani:2005, which adds a second penalty to the differences of the coefficients. Second, meinshausen:2006 showed that using only one penalty factor implies an inherent conflict between model selection and shrinkage estimation, leading to many noise variables in the final variable set. Thus, meinshausen:2007 introduced the relaxed lasso (relaxo) with a new parameter controlling the applied shrinkage. Third, lasso introduces a potential bias for large coefficients fan:2001, which can be tackled by adding weights to the regularization term -- i.e. the adaptive lasso (adalasso) by zou:2006.
Building on the strengths and weaknesses of the ridge and lasso estimator and combining both penalties, zou:2005 introduced the elastic net (glmnet). This algorithm includes a complexity parameter controlling the strength of the regularization and a mixing parameter between ridge and lasso regression.
The next three models also belong in the class of shrinkage estimators, yet they do not alter the coefficients but the covariates themselves. Specifically, they assume that the variables can be described as a linear combination of a reduced set of factors and loadings. These are principal component regressions (pcr), independent component regression (icr) by comon:1994 and partial least squares (pls) by wold:2001.
Until now the models assumed that the estimated parameters are fixed and driven by an unknown underlying data-generating process. With the upcoming models we want to venture into Bayesian territory and thereby assume that the parameters are random variables following certain distributions for which we can apply prior knowledge. As before, starting with a simple alteration of the linear model in which the coefficients follow a Student-t distribution (bayesglm) gelman:2008. Additionally, we introduce shrinkage via the Bayesian ridge regression (bridge) and the spike-and-slab (spikeslab) prior ishwaran:2005.
Unsurprisingly, there is also a Bayesian alternative for the lasso estimator (blasso) introduced by park:2008. The results are very similar to the original lasso algorithm, but the Bayesian treatment has the advantage that the penalty factor has not to be determined via cross-validation but can be implicitly derived in a fully Bayesian fashion. To further reduce the time-intensive and computationally demanding calculations, cai:2011 combined an empirical Bayes (EB) method with lasso to create the empirical Bayesian lasso (eblasso).
As a last subgroup of the Bayesian models, we introduce "global-local" shrinkage estimators that, as the name suggests, introduce a global shrinkage parameter pushing the coefficients uniformly towards the origin, while the local parameter allows for coefficient-specific deviations. Although there are many different priors to choose from, we settled for the Dirichlet–Laplace prior (dlbayes) by bhattacharya:2015, the horseshoe prior (horseshoe) by carvalho:2010 and the extended horseshoe prior by bhadra:2017 named horseshoe+ (horseshoePlus), which exhibits significant improvements in the case of "ultra-sparse" signals (i.e. nearly all coefficients are zero).
Now coming to the last subgroup of the (generalized) linear models: ensembles of linear models. We deploy a straightforward gradient boosting algorithm (bstlm), which is based on linear models as weak learners.
Model Averaging
Variable selection, either in the way of regularization or shrinkage, can lead to an over-correction due to too many variables being penalized and thus biased estimates of the remaining covariates and too narrow confidence intervals as the inherent model uncertainty is not taken into account by focusing on only one equation lukacs:2010. Therefore, we introduce three model averaging models that can overcome these issues by combining multiple equations of the same base model. First, Mallow’s model averaging (mma) by hansen:2007, which forms the final model by weighting multiple nested models based on minimized mean squared forecasting errors. Second, jackknife model averaging (jma), introduced by hansen:2012, calculates the weights by minimizing the leave-one-out (LOO) cross-validated residuals. Third, we introduce the current state-of-the-art credit risk satellite model: Bayesian model averaging (bma) raftery:1995. With this technique, the weights are based on Bayes’ theorem and the posterior model probabilities thereof.
Exponential Smoothing
Forecasts based on exponential smoothing have a long and successful track record going back to the late 1950s brown:1957. We implement a simple exponential smoothing (es) algorithm with additive errors and no trend or seasonal component based on hyndman:2008.\footnote{The model type for error, trend and seasonality -- in our case ANN -- is implicitly chosen by the sample-size-corrected AIC. This is inline with our data transformations, which include de-seasonalizing and differencing. The exogenous variables are also chosen based on this criterion.} As the choice of model types is crucial for the performance of the model, svetunkov:2016 introduced complex exponential smoothing (ces), which avoids the artificial distinction of a time series in level, trend, and seasonality.
(Generalized) Additive Models
The class of generalized additive models (GAMs) marks the point from which the presented models become more non-parametric and data-driven in their estimation routines. A feature which is often attributed to machine learning models.
GAMs allow more flexibility compared to a standard linear model due to the built-in smoothing function. Additive models were originally proposed by friedman:1981, and the first model in this class will be the one which they proposed in their seminal paper, namely Projection Pursuit Regression (ppr). Furthermore, we implement boosted smoothing spline (bstpline), which utilizes the same gradient boosting algorithm as the linear version (i.e. bstlm). Lastly, a boosted generalized additive model (boostgam) as outlined by schmid:2008 is implemented, combining the features of the first and second models.
Multivariate Adaptive Regression Splines
Another regression technique, which was proposed by Jerome Friedman, is the multivariate adaptive regression spline (mars) friedman:1991. Somewhat similarly to (generalized) additive models, the setup is now fully non-parametric, consisting of linear combinations of hinge functions. This extension of linear models allows the implicit modelling of non-linearities and interactions between variables.
Support Vector Machines
Support vector machines (SVM) have a long-standing history in providing robust classifications for linear and non-linear problems. The work of drucker:1997 introduced the concept of support vectors to regression problems. For our model comparison, we utilize a specific SVM model called a L2-regularized L1-loss support vector regression (svr), which uses a linear kernel and, as the name states, performs ridge regularization on the covariates hsieh:2008.
The second model in this class is a relevance vector machine (RVM) which has the same functional form as SVMs but uses Bayesian inference to estimate the equation tipping:2001. The advantages are that, compared to SVMs, RVMs are highly sparse and no prior cross-validation is needed to tune the cost function. We deploy the model with the standard Gaussian radial basis kernel (rvm), which delivered the best predictive performance.
Gaussian Process
A Gaussian process (GP) regression williams:1996 is a Bayesian kernel machine which bridges the gap between Bayesian linear models or spline models and SVMs. A GP builds on a covariance matrix which is used to assess how much information contiguous observations convey about each other. By applying a kernel matrix, in our case a polynomial kernel (gp), the estimation is extended into a non-linear realm karatzoglou:2006.
Tree-Based Models
Similar to the class of (generalized) linear models, tree-based models also contain a multitude of different approaches, ranging from simple to more complex. The first subclass are rule-based tree structures such as classification and regression trees (cart) breiman:1984 that recursively split the data set to make predictions about the outcome variable. As a further advancement in this field, hothorn:2006 introduced conditional inference trees (ctree) which tackle the issues of overfitting and selection bias prevalent in CART by introducing permutation tests during the partitioning.
In the realm of tree-based models, ensemble learning methods are more common than in the case of linear models. One of the most well-known ensemble models is random forests (rf) by breiman:2001. The method combines multiple decision trees which are trained on different subsampled parts of the data and with random subsets of variables. The results of all trees are averaged across to form the final prediction (i.e. bootstrap aggregation or bagging).\footnote{We want to note that the standard implementation of bagging in the used R library ranger is not sensitive to time series data as it assumes iid data. Currently, there exists no available library that includes bootstrap methods for dependent data.} The low bias in the results comes at the price of large variance. In order to tackle these issues, geurts:2006 suggested extremely randomized trees (ert) which use the whole data set for each tree instead of subsampling and randomize the splitting rule for each node.
Another popular modelling technique is gradient boosting friedman:2001, which is again an ensemble method, but in contrast to random forests, gradient boosted trees (or gradient boosting machines, GBM) are built sequentially. Each new tree improves the shortcomings of former trees, combining the results along the way. We opt for a standard Gaussian loss function (bsttree). It may come as no surprise that there also exists a Bayesian version of GBM called Bayesian additive regression trees (bart) established by chipman:2010. Instead of combing the trees via a learning rate, an iterative backfitting Markov chain Monte Carlo (MCMC) algorithm is used. As we will see in Section (ref), this flexible setup turns out to be the winning model.
The last competitor in this class of models is called node harvest meinshausen:2010 and settles itself between easy-to-understand regression trees and more accurate ensembles like random forests (enstree). The model delivers sparse results by initially creating random nodes and then finding suitable weights for each node based on an empirical loss function.
Neural Networks
The last class of competitor models is one of the earliest and most commonly used techniques in machine learning: neural networks mcculloch:1943. As a first model, we deploy a deep neural network (nn) with three hidden layers using resilient backpropagation with weight backtracking riedmiller:1994. The last model will be again Bayesian, namely a Bayesian regularized neural network (brnn). This model fits a two-layer neural network as described in foresee:1997. In contrast to classical neural networks, Bayesian networks are graphical models in which each node represents a variable with probabilistic relationships among them.
The increasing complexity of the models outlined above is also reflected in the number of parameters that need to be set before a model can be estimated. These hyperparameters control complexity and are thus crucial ingredients to the overall outcome of the model. Although many machine learning libraries provide default values for most parameters, olson:2018 showed that tuned hyperparameters can significantly reduce the variance compared to the out-of-the-box values.
There are different methods to find the most suitable set of hyperparameters (see feurer:2019 for an overview), again with different layers of complexity. We choose a model-free, non-black-box method for this task: grid search. After specifying a set of values for each parameter, grid search evaluates each set combination. The drawback of grid search is the possible large number of combinations that must be evaluated. However, to tackle this problem, we combine grid search with expert judgment. Specifically, on the one hand, we conduct research on the proper parameter space for each model and, on the other hand, we manually fine tune the grid to reduce the computational burden. In combination with a highly efficient implementation by the caret library kuhn:2020, the whole process stays feasible and very transparent.\footnote{We implement our own grid search for models which are not implemented in caret.}
For the performance measure we follow hyndman:2006 and stay within the realm of the well-known scale-dependent measures. The used data set stays the same across the models and we thus do not need to take the possible different data characteristics in consideration. Therefore, we choose the mean absolute error (MAE) as our indicator on how well a model predicts the default probability. Since some models exhibit an unstable forecasting behavior, we disregarded the more widely used mean square error (MSE) or root mean square error (RMSE) as they are more sensitive to outliers armstrong:2001.\footnote{For completeness' sake, we also implemented the RMSE, mean absolute percentage error (MAPE) and mean absolute scaled error (MASE). However, the tuning parameters and model rankings are nearly identical.}
Finally, in order to generalize the parameters for different sample lengths and time periods, a cross-validation strategy is introduced. More precisely, we apply a rolling-origin evaluation, starting with an initial training set of $T_{t,1} = 1,\dots,41$ and a fixed test (or holdout set) of $T_h = 12$, representing the 12 quarters we want to forecast and base our model comparison upon. In each iteration of the cross-validation, the training set is extended by one observation, leading to a total of 12 estimation rounds and a final training set of $T_{t,12} = 1,\dots,52$.\footnote{Due to the limited data availability and the large set of predictors, we do not deploy a fixed, moving window.} The best performing sets of hyperparameters are saved and used in the following estimations. The fairly large initial training set is justified by our sizeable set of $n = 40$ predictors, which allow us to estimate each model without a constraint on their ability to restrict the set of predictors with regularization. As a result, to obtain the hyperparameters for 37 models, around 15,000 estimations have to be carried out.\footnote{An overview of the tuned hyperparameters per model can be found in Table (ref).}
For the same reasons as outlined above, the estimation of the 43 models is carried out by deploying the same cross-validation strategy. The only difference is that the initial training set is $T_{t,1} = 1,\dots,4$. Hence, we provide for the initial estimation only one full year of data to the models. In combination with the large set of variables, this enables us to gauge the performance of the models under such extreme overfitting conditions. Using such as small data set also reflects reality, in which data availability is often limited to incomplete data or only a few recent observations. Given these circumstances, we would assume that especially regularized models and machine learning models, which are often advertised for their ability to handle such cases well bzdok:2018, will fare well in the comparison.
As with the hyperparameter tuning, the MAE is again used as our main evaluation factor. After each cross-validation window, the pseudo out-of-sample forecast is computed and compared to the actual PD time series. We thereby follow the idea outlined by hastie:2017 that for the purpose of model comparisons based on conditional forecasts, we do not need to focus on the causal relationship between the variables but identify statistical dependencies that are stable over time. Thus, this section will not present any results on the estimated coefficients, but only focus on the forecasting performance.
Before we come to the evaluation of the models in terms of their predictive accuracy, we need to make sure the comparison is as fair as possible. Specifically, if models exhibit certain instabilities during the estimation or deliver no proper output at all, combining these predictions in one final score can lead to distorting effects on the overall ranking of the models.\footnote{In the majority of cases the instabilities occurred due to over-regularization effects, which led to flat forecasts or algorithms not being estimable in the case of $T \ll n$.} Thus, if 25% of the predictions of a model fall in such a category, the whole model is dropped for further processing. As a consequence, 21 models have been dropped in due process -- 16 of which due to estimation instabilities and 5 could not estimate Eq. (ref) in the case of $T \ll n$. Due to this process, the multivariate adaptive regression splines model has been dropped and thus the whole category, leaving 22 models in eight categories.
Figure (ref) plots the evolution of the out-of-sample performance criterion across all cross-validation windows. The dashed vertical lines indicate the overfitting threshold after which sample $T$ is larger than the number of covariates $n$. At first glance, the figure indicates a significant variation of prediction accuracy across model categories, but also within the groups. Another key takeaway, without diving into the details, is that nearly all models are able to improve and stabilize forecast accuracy once the overfitting threshold is overcome. This is especially true for models which have no built-in feature to treat the overfitting problem. However, even models with such features, like BMA, struggle at first until a certain length of the training set is reached. This confirms our expectation that there is indeed a large variation in results across models and that thus a proper model comparison is needed to gain more insight into the driving factors.
Indeed, there is more to be learned from Figure (ref). The group of (generalized) linear models still contains most competitors (a total of eleven). For the first 35 iterations, the models depict a somewhat similar, erratic forecasting behavior. However, from then onwards, there is a clear separation of the majority of the models towards a MAE region of 0.2, whereas the Bayesian generalized linear model shows an increasingly worse performance. In contrast, (generalized) additive models attain a more stable forecasting pattern earlier in the process. The worst predictions are obtained from exponential smoothing and the Gaussian process, in both we see a large increase in forecasting error at around the 20th iteration. Interestingly, the overall best performance stems from the tree-based models. Specifically, classical random forests and Bayesian additive regression trees are able to accurately pick up the correct signals even in a remarkable overfitting setting, thereby keeping a very stable profile across the iterations. Lastly, the categories of model averaging, support vector machines and neural networks show a somewhat similar pattern of irregularities for the first 30 to 35 iterations and stable performance afterwards. At first glance, it seems that the state-of-the-art credit risk satellite model BMA (red line) already has strong competition.
After descriptively reviewing the evolution of the out-of-sample MAE it would be difficult to state which model performed best and what the actual difference between the models would be. In order to gain more insight into the ranking of the results, we could resort to parametric statistical tests (F-test or t-test) or to their non-parametric alternatives (Friedman test friedman:1937 or Wilcoxon signed-ranks test wilcoxon:1945). However, all of these allow only a pairwise comparison, which would mean 462 comparisons per cross-validation iteration and hence 22,638 in total. Moreover, the parametric tests rely on strong assumptions about the distribution of prediction errors, which seem to be at least partly violated for most models, as can be seen in Figure (ref). Another popular method to compare the accuracy of forecast methods is the pairwise Diebold-Mariano test diebold1995 and its multivariate alternative mariano:2012. However, even with the multivariate test, we would still be left with one test per cross-validation iteration, i.e. 49 tests, if we were to combine all models at once. Even without conducting the test, we can already assume that the equal predictive accuracy (EPA) hypothesis would not hold on such a diverse set of prediction models and thus refrain from implementing the tests.
However, the question of how to compare multiple forecasting models across multiple cross-validation runs has been tackled before. Most notably, koning:2005 ranked the results of the M3 competition using two non-parametric tests: multiple comparisons with the mean (ANOM) and multiple comparisons with the best (MCB). On the one hand, the new testing strategy allowed them to provide new insights on the statistical significance of the comparative model performance. On the other hand, both tests have received criticism due to the binary nature of the conclusion why certain models performed better or worse. Thus, demsar:2006 introduced a third test as an alternative based on the Nemenyi test nemenyi:1963. The test ranks the performance of each model across the various data sets or cross-validation iterations, averages over the ranks and produces confidence bounds. From the (non-)overlapping confidence bounds one can deduce if the models are statistically different from each other.\footnote{We are aware of the possible drawbacks of using null hypothesis significance testing (see, e.g., benavoli:2017 for an overview) and the Bayesian alternatives that could help with such issues (see, e.g. calvo:2018). However, to the best of our knowledge, there is currently no (publicly available) Bayesian hypothesis testing strategy that can handle multiple models and multiple cross-validation results at once.}
Figure (ref) depicts the results of the Regression for Multiple Comparison with the Best (RMCB), which is the regression-based version of the Nemenyi test introduced by svetunkov:2020. The test constructs a simple linear regression with the model ranks as dummy variables and uses the estimated coefficients and confidence intervals to determine the performance differences. The main difference between the Nemenyi test and RMCB is the underlying critical distance. For the former it is a atudentized range distribution, while the latter uses a Student's t-distribution. This leads to narrower confidence bounds for the RMCB, which can be helpful with smaller sample sizes.
The models on the x-axis in Figure (ref) are sorted by the mean rank which they achieved across the cross-validation, depicted by the points in the figure, while the vertical lines reflect the confidence bounds. Thus, this test allows us to reveal the winning model: Bayesian additive regression trees (bart). However, as indicated by the dashed line, there are seven more models that are not statistically different from the winning model on a 5% significance level. Particularly, BART is closely followed by the spike-and-slab prior (spikeslab) and Random Forests (rf). The following section will provide a deep dive into BART and give more details on the estimation and results.
Within these eight best performing models, two belong to the group of tree-based models (bart, rf), four are (generalized) linear models (spikeslab, icr, lasso, pcr), one is a neural network (\textbf{nn}) and one is a support/relevance vector machine (\textbf{rvm}). Given the fairly long cross-validation period in which overfitting prevailed (36 out of 49 iterations), it is remarkable that the tree-based models, neural network and relevance vector machine are able to provide such accurate forecasts without the need for regularization. In contrast, it comes as no surprise that all of the linear models use some form of regularization to tackle the overfitting issue.
Lastly, the red indicator depicts again the current state-of-the-art satellite model among central bankers. Within our framework and given these results, we can deduce that Bayesian model averaging (bma) is significantly worse than the first eight models, i.e. the winning group. Hence, for the use case of Austrian corporate default probabilities, these eight models would, on average, deliver more precise forecasts compared to BMA.
From the 43 implemented models we started out with and the 22 models that remained, Bayesian additive regression trees (BART) by chipman:2010 turned out victorious. The following paragraphs will introduce the model in more detail, give insights into why the model worked that well and provide detailed results.
The framework consists of two parts, a sum-of-trees model and regularization priors on the parameters that constraint each tree. BART approximates the function $f$ from Eq. (ref) by
whereas $N$ binary trees $\mathcal{T}_{j}$ are used, each $j^{th}$ tree with a vector of $\mathbf{m}_{j} = (\mu_{j1},\dots,\mu_{jb_j})'$ terminal nodes and $b_j$ leaves.
The second part imposes a set of regularization priors over the grown trees $p(\mathcal{T}_{j})$, the model parameters $p(\mu_{ij}|\mathcal{T}_{j})$ and error variance $p(\sigma)$. These priors ensure that no individual tree is too influential in the sum of trees. Specifically, we want to highlight the prior on the terminal node parameter $\mu_{ij}$, representing the effect of a tree,
The prior ensures that the tree parameters are shrunken towards zero, thereby constraining each weak learner. The variable $k$ is the prior probability that $E(\mathbf{y}|\mathbf{X})$ is within the range of $\mathbf{y}$. The prior standard deviation $\sigma_{\mu}$ is related to the gradient boosting shrinkage parameter of friedman:2001, which also balances the effect of each tree. For more details on the model and inference, we refer the interested reader to the original paper.
The first step in our forecasting framework is the tuning of the hyperparameters. BART, as many other machine learning models, offers the possibility to tune all prior hyperparameters. However, chipman:2010 also specify out-of-the-box parameters that work well on a range of different data sets. More precisely, one can tune the number of trees $N$, the base ($\alpha$) and power ($\beta$) parameter for the prior $p(\mathcal{T}_{j})$, the above-mentioned variable $k$ for $p(\mu_{ij}|\mathcal{T}_{j})$ and the parameters for the inverse chi-square distribution on $p(\sigma)$, $\nu$ and $q$.
As outlined in Section (ref), we use a grid of potential hyperparameters and the combinations thereof, centered around the default values by chipman:2010. Nevertheless, the default values returned the best prediction accuracy, and we will thus not go into the details of the tuning process. There is one exception for which we diverted from the default value, namely the number of trees $N$. chipman:2010 state a default value of $N=200$ as a larger number of trees increases BART's representation flexibility and thus predictive capabilities. However, they also state that BART can be used for variable selection when the number of trees is reduced. The more trees are grown, the more irrelevant covariates are mixed with relevant ones, diminishing its selection effectiveness. When the number of trees is reduced, BART endogenously picks the more relevant variables. Given our large overfitting period in the training sample, a lower number of tress ($N=50$) achieved the best results.
Until now, for the purpose of a unified model comparison, we concentrated on point forecasts and the deviation from the true default probability. This was a deliberate choice in order to focus on the prediction accuracy of each model and to keep the set of results manageable. Nevertheless, the uncertainty surrounding the point forecasts is at least as important as the forecast itself. Thus, Figure (ref) shows the 12-step ahead forecast from the BART model including prediction intervals. The forecast covers the whole length of the training sample, starting after the initial four quarters as the first cross-validation iteration. The solid blue line is the estimated posterior mean, while the dark blue shaded area represents the 80% prediction interval and the light blue area the 95% interval. The solid black line is the actual PD time series. Plotting all cross-validation results (49 in total) would lead to many overlapping points forecasts and indistinguishable prediction intervals. We thus decided to only show five non-overlapping predictions, which nonetheless span all forecasted quarters. The vertical dotted lines indicate the predicted region. Especially in this aspect we can gauge the Bayesian inference as the predictive distribution can simply be calculated from the posterior draws, thereby incorporating the inherent parameter uncertainty.
Given our long forecasting horizon, there are some deviations between the solid blue and solid black line. Especially in the first two segments -- from around 2004 to 2007 and 2007 to 2010 -- BART was not able to pick up proper signals from the data to detect the surge in default probabilities in the first quarter of 2006 and during the financial crisis of 2007 and 2008. In the subsequent segments, the accuracy of the model increases as the blue line starts to trace the peaks and troughs of the actual PD. The prediction interval keeps a reasonable width over the whole time span. However, the uncertainty around the estimates is undeniable and would need special attention when being used as a credit risk satellite model.
Overall, the non-parametric Bayesian approach coupled with the adaptive weak learners seems to be a very sensitive, capable model that worked well in our setting. Given that many machine learning models depend critically on the chosen set of hyperparameters, the performance of BART with default setting is quite remarkable. The success of BART can be seen in the growing literature in various fields and the technical extensions that have been proposed. Particularly, BART with heteroscedastic errors pratola:2020, BART in (non-linear) VAR settings huber:2020b, multinomial logistic regression via BART murray:2017 and survival models sparapani:2016. A more detailed overview of the technical extensions can be found in hill:2020.
As we have seen in last section, even the best models are fraught with uncertainty regarding the point prediction. This inherent uncertainty regarding econometric models is undisputed in the literature and relates to unknown data generating processes, misspecified models and a generally complex reality the models try to replicate. In order to hedge against such uncertainties, bates:1969 introduced in their seminal paper the concept of forecast combinations. The idea is straight-forward: there is no one true model, there are only different approximations of the data generating process. The underlying models have their own strength and weaknesses, which, when combined, should yield an overall better forecast. Even though such forecast combinations can be hard to interpret in terms of marginal effectiveness of the coefficients, we can again take up the point made by hastie:2017 that especially with comparisons among prediction models, a stable statistical dependency outweighs the underlying causal relationships.
In the last five decades since bates:1969, a wide range of combination methods have been suggested. For the purpose of this paper, we focus on three groups: simple combination methods, regression-based combination methods and eigenvector-based combination methods. For the simple methods, we choose the naive average, which weights all models equally (AVG) and the Newbold/Granger method newbold:1974, which calculates the weights from the estimated mean squared prediction error (MSPE) matrix (NG). The second group is still based on linear functions of the individual forecasts, but the weights are determined using a constrained least squares (CLS) regression. Finally, the standard eigenvector-based approach by hsiao:2014 uses, unlike newbold:1974, a normalization condition that leads to an unconstrained minimum of the MSPE (SEA). These methods (and more) have been implemented by weiss:2018b.\footnote{Besides a wide array of combination method, the implementation by weiss:2018b also allows a dynamic version of combinations, which is related to the idea of time series cross-validation. However, in our setting the normal static version achieved better results than the dynamic version.}
Moreover, we conduct the analysis with three different scenarios: first, in a naive approach we combine all 22 models; second, we only use the winning eight models which have been determined by the Nemenyi test in Figure (ref); third, we combine the best model of each category, as outlined in section (ref).\footnote{Using all 22 models with the Newbold/Granger and the eigenvector-based method leads to mathematical problems and is thus excluded.}
Figure (ref) shows, similar to Figure (ref), the evolution of the mean absolute errors (MAE) across the cross-validation iterations. The four panels reflect the forecast combinations group outlined above with the three different model scenarios -- all models (green line), top eight winning models (purple line) and the best models of each category (orange line). In order to properly frame the results of the forecast combinations, we additionally plot the performance of BART (blue line) and BMA (red line). First of all, the sole comparison between BMA and BART again emphasizes the significant increase in terms of predictive accuracy that BART delivers. However, while the simple average method (AVG) is not able to fully outperform BART, the constrained least squares version (CLS) is able to combine the underlying models in a way that beats BART across nearly all cross-validation iterations. The difference becomes even clearer in the case of the Newbold/Granger method (NG) and the standard eigenvector-based approach (SEA). While the eigenvector approach depicts a somewhat erratic behavior, driven by the combined models, the Newbold/Granger approach delivers a significant improvement across all cross-validation iterations.
In order to get a clearer picture of the results, Table (ref) shows the average out-of-sample MAE and ranks for BART, BMA and the forecast combinations, calculated over the cross-validation iterations. The numbers in parentheses indicate the average forecasting accuracy and rank relative to BART. On average, BART is able to beat two of the three simple average scenarios and is nearly on par with the top eight scenarios. However, all other combination methods are able to improve the average prediction accuracy significantly. The best accuracy is achieved by the Newbold/Granger method using the group-wise best performing models. Although this setup includes rather badly performing models, such as exponential smoothing and the Gaussian process, the Newbold/Granger approach is able to calculate the weights in a way that extract only the positive features from the underlying models.
Since the Great Financial Crisis, the use of stress tests as a tool for assessing the resilience of financial institutions to adverse financial and economic developments has increased significantly. One key part in such exercises is the translation of macroeconomic variables into default probabilities for credit risk by using macrofinancial linkage models. A key requirement for such models is that they should be able to properly detect signals from a wide array of macroeconomic variables in combination with a mostly short data sample.
The current state-of-the art satellite model for PD translation is Bayesian model averaging (BMA) raftery:1995. It has a long track record as being a reliable tool for generating scenario-conditional projections for credit risk and is being adopted by more and more central banks and institutions. However, with the easier access to regression models and the advent of new predictive models in the field of machine learning, the question arises if there are other models that could deliver better results.
The aim of this paper is to conduct a systematic forecast comparison with a large number of different regression models to find the best performing credit risk satellite model. The best model is evaluated for the ability to precisely forecast default probabilities conditional on a standardized set of macroeconomic variables as provided to financial institutions by the esrb:2020 for the EU-wide banking sector stress test. We implement a total of 43 models which we assigned to 9 categories, ranging from conventional statistical models to more recent machine learning methods. Additionally, we combine the models with different forecast combination approaches to further gauge their potential accuracy. For the purpose of this paper, we implement a framework that allows us to conduct this comparison with a standardized data set for all models, to tune the respective hyperparameters for each model and to cross-validate the results based on recursive pseudo out-of-sample forecasts. The data used in the modelling exercise refer to Austria and include as dependent variable a measure for the probability of default for the non-financial corporate sector and macroeconomic and financial data as independent variables.
Our results indicate that there are indeed better performing models than the current state-of-the art model. Specifically, a group of eight models significantly outperforms BMA in terms of their capability to forecast default probabilities. Five of these eight models belong to the category of machine learning models. Among these models, there is also the overall winner of the model comparison: Bayesian additive regression trees (BART) by chipman:2010. The combination of a flexible sum-of-trees model and well calibrated regularization priors gives BART the advantage needed to outperform all other models, especially in situations where overfitting is prevalent. Additionally, given that the majority of winning models has not been explicitly covered in the literature yet, our comparison sheds light on potential other credit risk models to be further investigated. We specifically highlight the advantages of machine learning models in the context of default probability prediction and more generally their applicability in high dimensions where overfitting prevails. Lastly, as most forecast combinations even outperform BART, we show that simple combination techniques can help to further hedge against model uncertainty and boost predictive accuracy.