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.
113,810 characters · 23 sections · 88 citation commands
Macroeconomic Forecasting and Machine Learning
\thispagestyle{empty}
JEL Codes: C22, C52, C53, C55. \\ Keywords: Deep Neural Network, Quantile Regression, Unemployment at Risk, Regularization, Big Data, Out-of-Sample Validation. \vskip3cm \onehalfspacing \setcounter{page}{1}
Forecasting has undergone a profound transformation in the 21st century, driven by advancements in methodology, computational power, and data availability. The origins of this transformation can be traced back to the late 1990s, when the field of economics began grappling with the challenges and opportunities presented by what is now termed “Big Data”—the availability of large datasets with numerous predictors. This period marked the emergence of systematic efforts to develop tools capable of addressing the high-dimensional nature of these datasets with seminal contributions by Frank Diebold, Mario Forni, Marc Hallin, Marco Lippi, Lucrezia Reichlin, Jim Stock, and Mark Watson Reichlin1998, Watson1998, Diebold1998. They laid the foundation for a new era of forecasting, as presented at the World Meeting of the Econometric Society in the summer of 2000 DemolEtal2017,Diebold2021.
Since these early contributions, the field of macroeconomic forecasting has experienced significant progress. Many of these advances were developed within macroeconometrics, while others originated independently in different disciplines, often leading to valuable cross-fertilization. Taken together, these developments fall under the broader umbrella of machine learning (ML), which has become one of the most influential ideas of the 21st century, reshaping empirical modeling across science, industry, and society. Motivated by the success of these methods in other fields, economists have in some cases rediscovered and re-emphasized key ideas, giving them renewed momentum and sparking broader interest across the profession. At times, however, the connections between modern machine learning and the macroeconomic forecasting tradition are overlooked. This can create confusion about what is genuinely new: some economists, familiar with the history of macroeconomic forecasting, may see ML as “reinventing the wheel,” while others may regard it as an entirely new paradigm. This disconnect can hinder communication, mutual understanding, and trust, even within the same field, by obscuring both novelty and continuity in today's forecast innovations.
At a high level, three key ideas now underpin modern forecasting approaches: using high-dimensional data with appropriate regularization, adopting rigorous out-of-sample procedures for both testing and validation, and incorporating nonlinearities. Most of these elements have deep roots in the theory and practice of macroeconomic forecasting. In addition, over the past decade, there has been growing attention to predicting the full distribution of risks, reflecting increased interest in tail events and macroeconomic uncertainty.
In this paper, we take stock of these developments and synthesize them into a unified framework for predicting macroeconomic risk. In summary, our results show that the model can effectively capture emerging macroeconomic vulnerabilities—especially those related to downside risks—by extracting information from high-dimensional data through a flexible and general predictive framework. We find little variation in upside risk but substantial time variation in downside risk. Thanks to validation-based selection and appropriate regularization, the high dimensionality of the data and the complexity of deep neural networks do not lead to overfitting issues. While shrinkage proves essential in this context, allowing for nonlinearities provides only limited gains in predictive accuracy.
The rest of the paper is organized as follows. Section (ref) provides a summary of the evolution of the literature on macroeconomic forecasting. Section (ref) introduces the problem formally. In Section (ref), we apply the framework to forecasting the unemployment rate one month ahead. Section (ref) discusses the results, and Section (ref) concludes.
In this paper, we use the term “Big Data” as it is understood in macroeconomics: a setting in which model complexity is high relative to the number of observations DemolEtal2017, GiannoneEtal2021. In such environments, even simple linear models can overfit, especially when the number of predictors approaches or exceeds the number of observations. In-sample performance often fails to translate into reliable out-of-sample forecasting accuracy, as the high parameter-to-observation ratio leads to poor generalization. Two major strategies have emerged to address these challenges: dimensionality reduction via principal components and regularization through shrinkage. SW2002, SW2005JBES, ForniEtal2000, ForniEtal2005, and West2003 introduced dynamic factor models for large datasets and showed that these models provide an effective way to summarize the information in many predictors. For a comparative overview of dynamic factor models in macroeconomic forecasting, see SW2006, DagostinoGiannone2012, SW2016, and DozFuleky2020.
In parallel, shrinkage methods such as ridge regression, Lasso, and spike-and-slab regressions were developed to penalize complexity and mitigate overfitting. DemolEtal2008 introduced shrinkage methods in the context of forecasting with many predictors. They showed that ridge regression (which uses an $\ell_2$ penalty) is asymptotically equivalent to dynamic factor models when predictors are strongly correlated, a common feature of macroeconomic datasets. This literature and the connection between penalty functions and informative priors have led to the widespread adoption of Bayesian inference in large-dimensional systems, as Bayesian inference provides both computational tractability and principled uncertainty quantification, particularly through the development of Large Bayesian VARs BanburaEtal2010,Koop2013JAE,GiannoneEtal2015.
DemolEtal2008 also found that Lasso (which uses an $\ell_1$ penalty) achieves similar predictive performance, but the selection of predictors is unstable. BianchiEtal2022 studied elastic net (combining of $\ell_1$ and $\ell_2$ penalties) while GiannoneEtal2021 considered spike-and-slab priors (combining $\ell_0$ and $\ell_2$ penalties). Sparsity-inducing methods have gained significant traction in statistics and have been successfully applied across a variety of disciplines, from seismography and medical imaging to astronomy and signal processing SantosaSymes1986, DaubechiesDefriseDeMol2004, Tibshirani1996, BickelEtal2009, HastieTibshiraniWainwright2015. In economics, these methods have proven effective in financial applications—for instance, in constructing sparse and stable portfolios BrodieEtAl2009—as well as in micro-econometric settings, particularly for high-dimensional treatment effect estimation and instrument selection BelloniEtal2013. However, in macroeconomic forecasting, the use of sparsity-based techniques has remained relatively limited. Researchers have instead relied more heavily on ridge regression and factor model approaches. This reflects the fact that macroeconomic datasets typically exhibit strong co-movement among predictors, which makes sparsity less relevant and favors dense representations that capture common underlying dynamics DemolEtal2017, GiannoneEtal2021.
During the 2000s, most forecasting applications remained grounded in linear models, given their empirical success and interpretability. Nevertheless, researchers explored nonlinear methods. Notable early contributions include White1988,white1989learning,DieboldNason1990,SwansonWhite1997,HaefkeEtal1997,White2006, who employed neural networks for macroeconomic forecasting. Yet empirical evidence during that period suggested that nonlinearities yielded limited improvements over linear benchmarks HevaviEtal2004. White1988 and DieboldNason1990 found similar results for predictions of stock prices and exchange rates, respectively. Combined with the computational burden of estimating nonlinear models in high dimensions, this led to a continued emphasis on linear approaches. Recent increases in computational power and algorithmic efficiency have revitalized interest in nonlinear methods, especially deep learning. Models such as regression trees and deep neural networks (DNNs) can now be estimated and applied at scale, enabling richer approximations of complex functional relationships LenzaEtal2023,ClarkEtal2023,CarrieroEtal2024,HauzenbergerEtal2025.
Another major development was the growing role of out-of-sample (OOS) validation. While OOS evaluation has long been a cornerstone of econometric forecast comparison and hypothesis testing, it became especially prominent in high-dimensional settings. Seminal contributions by DieboldMariano1995, West1996, White2000, ClarkMcCraken2001,GiacominiKomunjer2005, GiacominiWhite2006, AmisanoGiacomini2007, ClarkWest2007, GiacominiKomunjer2005, and RossiSekhposyan2019 established OOS validation as a key tool for evaluating forecast accuracy ClarkMcCracken2013,GiacominiRossi2013,Komunjer2013. The shift toward OOS reasoning occurred in tandem with the Big Data revolution, where model complexity made in-sample fit unreliable. A key early contribution was SW1996, SW1999, who pioneered the systematic use of OOS model comparison in high-dimensional macroeconomic forecasting. In this context, in-sample criteria often fail to capture genuine predictive performance. Machine learning (ML) has pushed this logic further by embedding out-of-sample evaluation into every stage of model development. In modern ML workflows, performance on held-out data is not just used for testing, but also for selecting hyperparameters and designing model architectures. The logic of out-of-sample evaluation for hyperparameter tuning and model selection is deeply rooted in classical tools such as cross-validation and hierarchical Bayesian inference. These approaches offer a principled framework for regularization that naturally incorporates out-of-sample validation into the estimation process. For a discussion of the connection between hierarchical Bayes and out-of-sample validation for model selection in a Big Data setting, see GiannoneEtal2015, GiannoneEtal2021.
While earlier work in macroeconomic forecasting focused primarily on point predictions, the past decade has seen growing attention to forecasting the full distribution of risks. AdrianEtal2019 show that macroeconomic risk—especially on the downside—varies substantially over time and can be predicted using financial and macroeconomic indicators. Their original analysis focused on U.S. GDP growth and introduced the concept of growth-at-risk to summarize the evolving distribution of future outcomes. This perspective has since been extended to other countries and macroeconomic variables, including inflation and unemployment Figueres_Jarocinski_2020, AABG2021, kiley, Amburgey2023, boyarchenko2023outlook, Chernis2023JAE,Loria2024Inflation. Although this literature often incorporates Big Data, typically by using indexes or factors extracted from large panels of macroeconomic and financial variables, it has not fully exploited recent advances in machine learning, particularly the flexible incorporation of high-dimensional predictors, nonlinearities, and rigorous out-of-sample validation. Several recent studies represent important steps toward a more systematic integration of these components. For example, LenzaEtal2023 use quantile regression forests to forecast the distribution of euro area inflation, capturing nonlinear interactions between inflation and financial conditions. Hengge2019 integrates multiple measures of financial conditions and economic and political uncertainty to forecast economic growth. The IMF2024 applies deep neural networks to assess the predictive power of financial conditions and uncertainty for macroeconomic risk. furceri2024global use many macroeconomic and political factors to predict both the central trajectory of public debt and the uncertainty surrounding it across more than one hundred countries. There is also a long tradition of forecasting risk in finance: ForesiPeracchiJasa2005 forecast the cumulative distribution function of excess returns conditional on a broad set of predictors and CrumpEtal2024 combine density forecasts of stock returns to improve overall forecast accuracy. Out-of-sample evaluation beyond point forecasts also has a long tradition in macroeconomics. AmisanoGiacomini2007 study the evaluation of the accuracy of density forecasts. DieboldEtal1999,RossiSekhposyan2019 develop methods to assess the calibration of density forecasts. GiacominiKomunjer2005 consider the accuracy of quantile predictions Komunjer2013. JoreEtal2010,ConflittiEtal2015 study how to combine density forecasts to improve out-of-sample accuracy.
Our contribution builds on these advancements by systematically combining all key elements—regularization, nonlinearity, out-of-sample validation and testing, and full predictive distribution—into a unified high-dimensional forecasting framework. We adopt a nested modeling approach that not only integrates these components, but also allows us to assess the relative importance of each in shaping predictive performance.
The forecasting exercise is deliberately stylized and does not capture the full complexity of real-time macroeconomic data, nor significant nonlinearities. While we carefully avoid look-ahead bias by using a recursive design that mimics real-time information availability during model training, validation, and testing, we do not incorporate several key real-time features. In particular, our analysis omits the non-synchronicity of data releases, mixed-frequency observations, and data revisions—dimensions emphasized in the nowcasting literature GRS2008,BanburaEtal2013,Luciani,CascaldiEtal2024. Nonetheless, the core lessons from our analysis remain relevant, as nowcasting models often rely on similar modeling principles and estimation strategies, particularly the use of regularization to handle big data.
We are interested in estimating the following object: \[ P(Y_{T+h} < y \mid X_T = x_T) = F_{Y_{T+h} \mid X_T = x_T}(y) \] The function \( F_{Y_{T+h} \mid X_T = x_T}(y) \) is the conditional cumulative distribution function (CDF) of \( Y_{T+h}\) given the predictors \( X_T = x_T \). Here, \( X_T \) is a \( k \)-dimensional vector of predictors, with \( x_T \) being the observed realization of \( X_T \), and \( k \) is large, meaning the model must handle a high-dimensional feature space.
Quantile prediction has a long tradition in econometrics and statistics, where it has been used extensively to model conditional distributions beyond the mean Komunjer2013.
Instead of estimating the cumulative distribution function (CDF) directly,\footnote{The direct estimation of the CDF can be performed using distribution regression, as developed by ForesiPeracchiJasa2005.} we estimate its inverse—the quantile function—defined as \[ Q_{Y_{T+h} \mid X_T = x_T}(\tau) = y \quad \text{where} \quad \tau = F_{Y_{T+h} \mid X_T = x_T}(y). \] That is, the \(\tau\)-th conditional quantile of \(Y_{T+h}\) given \(X_T = x_T\) is the value \(y\) such that the probability of observing an outcome below \(y\) is \(\tau\).
Quantiles can be equivalently characterized as the solution to an optimization problem. Specifically, the \(\tau\)-th quantile choose the optimal \(q\) that minimizes an asymmetric linear loss function known as the quantile loss, or pinball loss:
This loss function penalizes under-predictions and over-predictions differently depending on the chosen quantile $\tau$. For \(\tau = 0.5\), the pinball loss function is symmetric: it assigns equal weight to positive and negative errors, and the resulting V-shape has equal slopes on both sides. However, when \(\tau \neq 0.5\), the function becomes tilted, introducing asymmetry in the penalization of forecast errors. When \(\tau < 0.5\), the loss function assigns more weight to negative errors—that is, when the predicted quantile \(q\) exceeds the actual outcome \(y_{T+h}\). The slope on the left side of the V is steeper, encouraging the model to shift the quantile downward so that a smaller proportion of observations fall below the estimated quantile line. When \(\tau > 0.5\), positive errors—when the predicted quantile is below the actual value—are penalized more. The right-hand slope becomes steeper, shifting the estimated quantile upward so that a larger proportion of observations fall below it. These asymmetries are clearly visualized in Figure (ref), which plots the pinball loss for \(\tau = 0.05\), \(\tau = 0.5\), and \(\tau = 0.95\), respectively. The loss is also known as lin-lin loss since it is linear on each side of the origin, with potentially different slopes ChristoffersenDiebold1997. It also goes by the names check loss or tick loss, due to its resemblance to a check mark when asymmetric, with a steeper slope on one side.
This tilting of the loss function reflects the fundamental goal of quantile estimation: to find the value \(q\) such that approximately \(\tau\) of the conditional distribution lies below it. As Figure (ref) shows, the shape of the loss function adjusts with \(\tau\), becoming increasingly asymmetric as the quantile moves away from the median.
The case \(\tau = 0.5\) can play a special role since it corresponds to the mean absolute deviation, a robust alternative to the squared loss. In contrast to the quadratic loss, which increases the weight of large errors quadratically, the pinball loss grows linearly with the size of the error. This makes it more robust to outliers and skewed distributions. Historically, the quadratic loss (i.e., squared error) has been widely used due to its analytical convenience in settings like OLS or ridge regression. However, with modern computational resources, estimating quantiles via the pinball loss is no longer burdensome. In our analysis, we use the median prediction for point forecasts and report results under both the \(\ell_1\) and \(\ell_2\) losses for robustness.
The quantile function is estimated by minimizing this loss. We treat the quantile function \( Q_{Y_{T+h} \mid X_T = x_T}(\tau) \) for each quantile level \( \tau \) as a general function of the predictors \( x_T \). We approximate this quantile function in a parametric form, where it depends on the predictors \(x_T\), on some parameters \(\theta, \psi \), denoted as $$ Q_{Y_{T+h} \mid X_T = x_T}(\tau) \approx f_{\psi,\tau}(x_T; \theta).$$
This function will be estimated using a deep neural network. The parameters $\psi$ specify the architecture and $\theta$ the trainable parameters. As described later, to regularize the estimation, we will augment the quantile loss with a penalty on the model's parameters.
The computation of the quantiles was once challenging computationally, up to a few decades ago, and there is a large literature using linear programming to compute conditional quantiles. Nowadays, we are able to solve these problems in high dimensions and also for non-linear relationships.
This section introduces Deep Neural Networks (DNNs) and their application to nonlinear prediction problems. We begin by outlining the general structure of a DNN, which consists of an input layer, one or more hidden layers, and an output layer. The core mechanism involves stacking multiple linear transformations interleaved with nonlinear `activation' functions\footnote{Mimicking the thresholding behavior in biological neurons}. For simplicity and clarity, we then illustrate specific cases with zero, one, and two hidden layers. These cases help relate DNNs to familiar econometric models such as linear regression and reduced-rank regression. Finally, we discuss the choice of activation function, focusing on the Leaky ReLU, a widely used and flexible nonlinear transformation.
Deep Neural Networks (DNNs), which extend classical Multi-Layer Perceptrons (MLPs), have become the foundation of modern machine learning. Originally introduced to overcome the limitations of shallow architectures, DNNs rose to prominence as advances in computation and optimization made training deep architectures feasible. Although the concept of MLPs dates back to the 1950s and 1960s---originating from early models like Rosenblatt's perceptron---the breakthrough came in the 1980s with the introduction of the backpropagation algorithm. The term "deep" entered common use in the mid-2000s, when researchers showed that networks with many hidden layers could be trained effectively, notably in work by rumelhart1986learning and hinton2006fast.
Graphically, the structure of a DNN is illustrated in Figure (ref). The figure shows a typical architecture consisting of an input layer, two hidden layers, and an output layer. Here, the input layer takes 130 features (which might represent macroeconomic indicators or other predictors). The hidden layers, assemble the inputs into the layer and apply nonlinear transformations, enabling the network to learn complex representations and shown to be `universal' (continuous) function approximators hornik1989multilayer. The output layer produces a scalar prediction. We define each of these components mathematically below.
Let \( x_t \in \mathbb{R}^k \) denote the input vector at time \( t \), consisting of \( k \) features. A deep neural network maps this input into an output \( y_t \in \mathbb{R}^n \) by passing it through a sequence of transformations. The architecture consists of an input layer, a series of hidden layers, and an output layer. Each hidden layer applies a non-linear activation function to a linear transformation of its input.
The input layer simply takes the feature vector \( x_t \) and feeds it into the network: \[ x^{(0)} = x_t \in \mathbb{R}^k. \] This vector represents the raw input data, for instance, a set of macroeconomic indicators or other predictors.
Each hidden layer performs two operations: a linear transformation followed by a non-linear function, the `activation' function. Let \( x^{(i)} \in \mathbb{R}^{r_i} \) denote the output of the \( i \)-th hidden layer ($i>0$). The transformation is given by: \[ x^{(i)} = a_i(W^{(i)} x^{(i-1)} + b^{(i)}), \] where:
The final layer takes the output of the last hidden layer and maps it to a prediction. For quantile regression, we assume a linear output layer: \[ Q_{Y_{t+h} \mid X_t = x_t}(\tau) \approx \gamma^\top x^{(d)} + c, \] where \( x^{(d)} \) is the scalar output of the last hidden layer, \( \gamma \) is the weight vector, and \( c \) is a constant.
While the general structure of a DNN is straightforward, the recursive notation used to describe multiple layers can appear cumbersome. To build intuition and relate the architecture to familiar forecasting models, we now present a sequence of examples. We begin with the simplest case involving only the output layer, and then consider models with one and two hidden layers.
Example 0: If the network has no hidden layers (\( d = 0 \)), the model reduces to: \[ Q_{Y_{t+h} \mid X_t = x_t}(\tau) \approx \gamma^\top x_t + c, \] which is the familiar linear quantile regression model. If a penalty is added on \( \gamma \), the model corresponds to a ridge-regularized quantile regression.
Example 1: With one hidden layer of dimension \( r_1 \), the network becomes: \[ x^{(1)} = a_1(W^{(1)} x_t + b^{(1)}), \] \[ Q_{Y_{t+h} \mid X_t = x_t}(\tau) \approx \gamma^\top x^{(1)} + c. \] If the activation function \( a_1 \) is linear, this corresponds to a reduced-rank quantile regression.
This model might resemble unsupervised dimensionality reduction like PCA. However, in the latter case, the projection coefficients \( W^{(1)} \) are chosen to maximize the explained variance of the predictors. Instead, here the projection is optimized directly to minimize prediction error. This makes the approach closer to reduced-rank regression than to unsupervised dimensionality reduction.
Example 2: With two hidden layers, the transformations become: \[ x^{(1)} = a_1(W^{(1)} x_t + b^{(1)}), \] \[ x^{(2)} = a_2(W^{(2)} x^{(1)} + b^{(2)}), \] \[ Q_{Y_{t+h} \mid X_t = x_t}(\tau) \approx \gamma^\top x^{(2)} + c. \] Again, if both activation functions are linear, the model reduces to a multi-layer generalization of reduced-rank regression.
The activation function introduces nonlinearity into the model, allowing the network to approximate complex relationships. The list of activation functions is long, from piecewise linear like the ReLU (Rectifier Linear Unit) to leaky ReLU, to sigmoidal, hyperbolic tangent, and so on. We use the piecewise linear Leaky ReLU:
\[ a(x; \alpha) =
\] where the hyperparameter \( \alpha \in [0,1] \) controls the slope on the negative side.
This function spans a range of behaviors:
When the activation function is linear ($\alpha=1$), the model itself becomes entirely linear. If all coefficients are left unrestricted, the number of layers becomes irrelevant: models with zero, one, two, or more layers will be observationally equivalent, and all reduce to standard linear quantile regression. However, in what follows, we will introduce regularization by placing a penalty on the parameters. In this case, while the model remains linear when ($\alpha=1$), different numbers of layers will generate different regularization structures, because the penalty depends on how the parameters are distributed across layers. As a result, models with more layers will shrink differently, even if they remain linear in functional form.
Different from the Leaky ReLU, this activation function will saturate the output for both large negative inputs and large positive ones. It also has the advantage of non-vanishing gradients.
Since the models we consider are very densely parameterized, we add a penalty to regularize the problem. Even without any hidden layers—i.e., with only the output layer—we already have more parameters than predictors, which in our setting is a large number. This high dimensionality is typical in economic applications, where the number of predictors often rivals or exceeds the sample size. In such situations, complex models are prone to overfitting: they may fit the training data well but perform poorly out of sample, particularly in real-time forecasting. A widely adopted solution to this problem is to introduce a penalty that shrinks the individual parameters towards zero, thereby controlling model complexity and improving generalization. Regularization through penalization has a long history in mathematics and statistics Tikhonov1963,hoerl1970ridge. It was introduced and studied in the context of forecasting with Big Data by DemolEtal2008. For a general discussion, see DemolEtal2017.
Specifically, we use a quadratic penalty, which corresponds to applying an \( \ell_2 \)-norm to all weights in the network: \[ \text{Pen}(\theta, \psi, \lambda ) = \lambda \left( \sum_{i=1}^{d}\left(\|W^{(i)}\|_F^2 \right) + \|\gamma\|_2^2 \right) \]
This regularization strategy is closely related to ridge regression. In fact, when the model includes only a single linear layer and no hidden layers, this formulation reduces to a ridge-regularized quantile regression model.
Ridge regression is known to be intimately related to data compression and regularization through factor structures, as shown in DemolEtal2008,DemolEtal2024. Similarly, one could consider alternative penalties. Popular choices include sparsity-inducing penalties such as the \( \ell_1 \)-norm (Lasso), the \( \ell_0 \)-norm (subset selection), or combinations such as elastic nets or spike-and-slab priors. We choose to retain the \( \ell_2 \)-norm in our setup because there is neither strong theoretical motivation nor empirical evidence in favor of sparse forecasting models. In particular, GiannoneEtal2021 show that shrinkage---rather than sparsity---is the key to improving forecasting performance in high-dimensional settings.
More broadly, as in any forecasting model, there is a fundamental trade-off between bias and variance. The degree of penalization, governed by the hyperparameter \( \lambda \), determines where the model lies along this trade-off. At one extreme, when \( \lambda \rightarrow \infty \), all parameters are shrunk to zero, and the forecast collapses to the empirical unconditional quantile---resulting in zero variance but high bias, since all covariate information is discarded. At the other extreme, when \( \lambda = 0 \), there is no shrinkage, and the model fits both signal and noise---yielding low bias but potentially high variance. Intermediate values of \( \lambda \) balance this trade-off. In our implementation, the optimal degree of shrinkage is selected through out-of-sample validation, as discussed in the following section.
In the discussion above, we implicitly distinguish between two sets of parameters:
The distinction between parameters and hyperparameters is related to hierarchical Bayesian model GiannoneEtal2015,GiannoneEtal2021. In our empirical analysis, hyperparameters are estimated or selected once based on out-of-sample performances over a validation sample and held fixed thereafter, while parameters, continuously re-estimated as new data becomes available (`online training'). This distinction is particularly important when we perform recursive estimation or online learning, as we do in our empirical implementation.
We begin by describing the data and how we split the sample; then we will discuss the estimation of parameters over a recursive training sample, hyperparameters over a validation sample, and finally evaluate model's performance over a test sample.
We use the well-established dataset originally developed for macroeconomic forecasting of the US economy using a large number of predictors. The dataset was constructed by SW2002 when they introduced principal component regression in the context of Big Data. The dataset is a reference for all subsequent work on Big Data and machine learning, including ForniEtal2005 and DagostinoGiannone2012 with factor models, DemolEtal2008 and GiannoneEtal2021 with shrinkage, BanburaEtal2010 with large vector auto-regressions, Coulombe2022 and CarrieroEtal2024 with DNN.
The dataset consists of 130 predictors, including various monthly macroeconomic indicators, labor market variables, house prices, consumer and producer prices, money, credit, and asset prices. The sample ranges from February 1960 to January 2024, and all the data has been normalized.\footnote{The dataset has been updated by McCrackenNg2016 and is available through the Federal Reserve Economic Data (FRED).}
Most of the empirical analysis focuses on the one-month-ahead forecast of the change in the unemployment rate. To assess the robustness of the results, we also consider a longer horizon (one year ahead) and the one-month-ahead growth rate of industrial production.
We conduct a simulated out-of-sample analysis by estimating the model recursively starting in January 1980. As detailed below, we split the out-of-sample period into two parts: we use the first part for selecting the architecture and the hyperparameters (validation), and the second part for testing the model's performance (testing).
To exploit the time series structure of the sample and use all available information, instead of splitting the data once into fixed training, validation, and test sets, we use an expanding training sample.
After each iteration, the window expands forward by one step, and the model is retrained on the new training window. For each step, we train the model on a data window from $1$ to $T$ and then test or validate the model on a future observation $T+h$. This expanding approach ensures the model continuously learns from past information. Alternatively, a rolling window with fixed size could be used.
Based on the trained predictions, we use a validation sample $T+h=T_1+1,\ldots,T_2$ to estimate the hyperparameters, and the rest of the sample $T+h=T_2+1,\ldots,T_3$ to evaluate the real-time model performance.
The validation and test samples are defined in terms of the forecast date $T+h$. As a result, the corresponding end of the estimation sample $T$ changes with the horizon. For the baseline case $h=1$, we set $T_1$ as January 1980, $T_2$ as January 2000, and $T_3$ as December 2024. The first forecast we evaluate is made in January 1980 ($T=T_1$) to predict February 1980 ($T_1+1$). The last forecast is made in November 1999 ($T=T_2-1$) to predict December 1999 ($T_2$).
In the robustness analysis, we also consider $h=12$. In that case, the first forecast is made in January 1980 to predict January 1981 ($T_1+12$). The last forecast is made in December 2023 ($T=T_3-12$) to predict December 2024 ($T_3$).
We approximate the quantiles of the predictive density as follows: \[ Q_{Y_{T+h} \mid X_T = x_T}(\tau) \approx f_{\psi,\tau}(x_T; \theta). \]
We estimate the parameters \(\theta\) for a given set of hyperparameters \(\lambda\) over recursive samples: \[ \hat \theta^{(T)}_{\lambda,\psi} = \arg \min_{\theta} \sum_{t+h=h+1}^{T} L_{\tau}(y_{t+h}, f_{\psi,\tau}(x_t; \theta)) + \text{Pen}(\theta, \psi, \lambda) \]
The out-of-sample predictions are: \[ \hat Q^{( \psi,\lambda)}_{Y_{T+h} \mid X_T = x_T}(\tau) = f_{\psi,\tau}(x_T; \hat \theta^{(T)}_{\lambda,\psi}) \]
Out-of-sample means that the estimation of \(\theta\) does not use observations \(t > T\). Since the number of parameters \(\theta\) is large, estimation is prone to overfitting. We mitigate this with limited architectural complexity and regularization. After each training step, the model is evaluated on the next point in time. This recursive estimation preserves the temporal structure of the data and ensures that no future information is used when training the model. In the test period, hyperparameters are fixed, and parameters are updated recursively.
To save computation time, we initialize each new optimization using the estimates from the previous step: e.g., \(\hat \theta^{(T+1)}_{\lambda}\) is initialized with \(\hat \theta^{(T)}_{\lambda}\). In machine learning terminology, this is called warm-starting from a checkpoint. In the initial iteration, we use 500 epochs; in subsequent steps, we reduce to 100 epochs. The expanding window ensures that the model evolves gradually, helping maintaining prediction quality despite the warm-starting strategy.
We train the models using four Nvidia Quadro RTX 8000 GPUs and the PyTorch library for architecture design and gradient descent. Stochastic gradient descent produces locally optimal solutions, but due to its randomness, results may vary slightly depending on the seed. We tested sensitivity to different random seeds and found the results to be robust.
Out-of-sample is used throughout. We split into validation, where we optimize hyperparameters and test where we freeze those hyperparameters and evaluate completely ut-of-sample.
As discussed above, the predictions \[ \hat Q^{( \psi,\lambda)}_{Y_{T+h} \mid X_T = x_T}(\tau)= f_{\psi,\tau}\left(x_T; \hat \theta^{(T)}_{\lambda,\psi} \right) \] still depend on the hyperparameters. We estimate these in the validation sample by minimizing the out-of-sample prediction loss: \[ \hat \lambda,\hat \psi = \arg \min_{\lambda,\psi } \sum_{T+h=T_1+h}^{T_2} L_{\tau}\left(y_{T+h}, \hat Q^{( \psi,\lambda)}_{Y_{T+h} \mid X_T = x_T}(\tau) \right) \]
We use the same loss function as in the training sample—the pinball loss—ensuring that the evaluation criterion aligns with the estimation target. This is consistent with the guidance in GiacominiKomunjer2005, who emphasize that out-of-sample evaluation of quantile forecasts should be conducted using the same tick loss used for estimation. While their focus is on comparing forecast accuracy across models or relative to a benchmark, the same principle applies to selecting hyperparameters and architectures.
Below is the set of hyperparameters we have trained in the validation period (Table (ref)). When the number of nonlinear layers is zero, the model reduces to a linear regression, and there is no need to choose hidden dimensions or activation functions.\footnote{Therefore, the number of validation jobs is approximately $ 3\times3\times 4\times 40 \times 8 = 11520 $.}
The performances in the validation sample are not genuine out-of-sample forecasts: they are subject to a look-ahead bias since the hyperparameters were selected to optimize performance over this very same evaluation period. To obtain an honest assessment of real-time performance, one must fix the hyperparameters based on the validation sample and then evaluate the model in a new test sample—never seen by the econometrician during either parameter or hyperparameter estimation. For this reason we evaluate model performance over the test sample \( T = T_2+1, \ldots, T_3 \), using the hyperparameters \(\hat{\lambda}, \hat{\psi}\) selected in the validation period. At each forecast time \(T\), parameters \(\theta_{\hat{\lambda}, \hat{\psi}}\) are re-estimated using data up to \(T\).
For \(T > T_2\), these predictions are fully real-time: \[ \hat Q_{Y_{T+h} \mid X_T = x_T}(\tau)= g_{\hat \psi,\tau}\left(x_T; \hat \theta^{(T)}_{\hat \lambda, \hat \psi}\right) \]
and the accuracy will be measured as \[ \hat \lambda,\hat \psi = \arg \min_{\lambda,\psi } \sum_{T+h=T_2+h}^{T_3} L_{\tau}\left(y_{T+h}, \hat Q^{( \psi,\lambda)}_{Y_{T+h} \mid X_T = x_T}(\tau) \right) \]
The difference in accuracy between validation and testing is informative and usually small if the hyperparameter space is not too large and if there are no major structural changes in the data.
Summing up, we split the out-of-sample in two parts, a validation and a test. We use the validation sample to optimize the hyperparameters, namely the architecture $\psi$ and the shrinkage parameter $\lambda$. Summarizing the outcome of this optimization directly is challenging because the hyperparameters involve multiple dimensions that interact in complex ways. For example, deeper networks with more neurons typically require stronger penalization to avoid overfitting. In the following section we summarize all models in a single dimension that represent complexity.
A simple way to understand model complexity is through the in-sample variance of the forecasts, following the insights of GiannoneEtal2015. Intuitively, more complex models fit more closely to the idiosyncrasies of the training data and therefore produce forecasts that fluctuate more in sample. Conversely, heavily penalized or highly constrained models generate smoother forecasts with lower variance.
A natural benchmark is the naïve model with full shrinkage, which produces constant forecasts equal to the unconditional quantiles. This corresponds to a simple linear regression on a constant and contains no information from the predictors. We use these naïve forecasts as the reference point against which more complex models are evaluated.
To formalize this idea, for any architecture $\psi$ and penalty parameter $\lambda$, we generate a sequence of fitted forecasts over the initial estimation window $t+h=h+1,\ldots,T_1$: \[ \widetilde Q^{\psi,\lambda}_{t+h}(\tau)=f_{\psi,\tau}\big(x_t;\,\hat\theta^{(T_1)}_{\psi,\lambda}\big), \qquad t+h=h+1,\ldots,T_1, \] where $\hat\theta^{(T_1)}_{\psi,\lambda}$ is estimated using all data available up to $T_1$. We use the tilde notation to emphasize that these are in-sample fitted forecasts, in contrast to the recursive forecasts $\widehat Q$ used in the validation and testing periods.
The mean fitted forecast is \[ \bar Q^{(0)}_{\psi,\lambda} = \frac{1}{T_1-h}\sum_{t+h=h+1}^{T_1}\widetilde Q^{\psi,\lambda}_{t+h}(\tau), \] and the forecast variance is \[ \operatorname{Var}_0(\psi,\lambda) = \frac{1}{T_1-h}\sum_{t+h=h+1}^{T_1}\Big(\widetilde Q^{\psi,\lambda}_{t+h}(\tau)-\bar Q^{(0)}_{\psi,\lambda}\Big)^2. \]
For each architecture, the unpenalized model $(\lambda=0)$ delivers the maximum forecast variance. Under full shrinkage ($\lambda \to \infty$), all weights are set to zero and predictions collapse to the unconditional quantiles. In this case, fitted forecasts are the same across all architectures, so we drop the index $\psi$ from the notation. The variance of these flat forecasts is \[ \operatorname{Var}_0^{\text{flat}} = \frac{1}{T_1-h}\sum_{t+h=h+1}^{T_1}\Big(\widetilde Q^{\text{naive}}_{t+h}(\tau)-\bar Q^{\text{naive}}\Big)^2, \qquad \bar Q^{\text{naive}}=\frac{1}{T_1-h}\sum_{t+h=h+1}^{T_1}\widetilde Q^{\text{naive}}_{t+h}(\tau). \]
We then define the complexity index \[ r(\psi,\lambda) = \frac{\operatorname{Var}_0(\psi,\lambda)}{\operatorname{Var}_0^{\max}} \;\;\in [0,1], \] where $\operatorname{Var}_0^{\max}=\operatorname{Var}_0(\psi,0)$ is the maximum variance achieved by the unpenalized model. Since the denominator corresponds to the variance of the unconditional quantile forecasts, it is the same for all $\psi$, and we therefore omit the architecture index in the notation.
In this definition, $r=0$ corresponds to flat, constant forecasts and the largest $r$ corresponds to the most complex, typically the most volatile, forecasts achievable with the model considered. Given the large set of predictors, the maximum will be close to 1, meaning that in sample it is possible to achieve an almost perfect fit.
The complexity parameter $r$ provides a natural one-dimensional summary of the bias--variance trade-off. Models with low complexity ($r$ near zero) are highly shrunk and produce very smooth forecasts. In this region, variance is low but bias is high, since useful predictive information from the data is largely discarded. By contrast, models with high complexity ($r$ near one) produce forecasts that closely track the training sample. These forecasts have low bias, because they use the predictors aggressively, but high variance, because they are sensitive to noise and prone to overfitting. Optimal forecast performance is typically achieved at intermediate values of $r$, where the model captures the systematic signal in the data without being dominated by idiosyncratic noise.
To operationalize this construction, we define a discrete grid of target complexity values, \[ \mathcal R = \{0.0, 0.1, 0.2, \ldots, 1.0\}. \] For each $r \in \mathcal R$, we search over a set $S$ of candidate hyperparameters consisting of architectures $\psi$ and penalties $\lambda \in \{0\} \cup \{10^{-3},10^{-2.9},\ldots,10^2\}$ on a logarithmic scale. For each pair $(\psi,\lambda)$, we compute the implied complexity $r(\psi,\lambda)$. We then select the combination that produces a complexity index closest to the grid point: \[ (\psi_r,\lambda_r) \in \arg\min_{(\psi,\lambda)\in S} \big|\, r - r(\psi,\lambda)\,\big|. \]
This mapping compresses the high-dimensional hyperparameter space into a single, interpretable index. Each target value $r$ on the grid corresponds to a representative model $(\psi_r,\lambda_r)$, and performance can be evaluated as a function of $r$ rather than of individual hyperparameters. In this sense, our approach generalizes the linear case studied by DemolEtal2008 and BanburaEtal2010, where shrinkage is the only dimension of regularization and the mapping between $\lambda$ and forecast variance is one-to-one. In the nonlinear setting considered here, by contrast, multiple combinations of architectures and shrinkage values can achieve a similar complexity index. To avoid reporting all possible combinations, we select one representative pair $(\psi_r,\lambda_r)$ for each grid point and summarize results only in terms of the effective complexity index $r$. The selected $(\psi_r,\lambda_r)$ should be seen as one representative among the set of models that yield the same complexity. The simplifying assumption is that models with similar complexity also display similar predictive performance, so that selecting hyperparameters in terms of $(\psi_r,\lambda_r)$ is approximately equivalent to selecting directly across $r$. Extensive checks over the full hyperparameter space confirm that this is a reasonable approximation, and that little is lost by summarizing performance in terms of the complexity index alone.
Formally, the model used for forecasting is based on the triplet $(\hat r,\hat\lambda_{\hat r},\hat\psi_{\hat r})$. The complexity level $\hat r$ is selected in the validation sample by minimizing the quantile loss, \[ \hat r = \arg\min_{r } \sum_{T+h=T_1+h}^{T_2} L_{\tau}\left(y_{T+h}, \widehat Q^{( \psi_r,\lambda_r)}_{Y_{T+h} \mid X_T = x_T}(\tau) \right), \] where $L_\tau$ denotes the quantile loss function. For each candidate value of $r$, the pair $(\hat\lambda_r,\hat\psi_r)$ is chosen among the set of hyperparameters that achieve the in-sample complexity closest to $r$. Thus the final choice $(\hat r,\hat\lambda_{\hat r},\hat\psi_{\hat r})$ should be interpreted as representative of this set. Empirical results based on the simplified procedure are very similar to those obtained from the full minimization over $(\lambda,\psi)$ discussed in the previous section, \[ \hat \lambda,\hat \psi = \arg \min_{\lambda,\psi } \sum_{T+h=T_1+h}^{T_2} L_{\tau}\left(y_{T+h}, \widehat Q^{( \psi,\lambda)}_{Y_{T+h} \mid X_T = x_T}(\tau) \right). \] This empirical similarity allows us to present the results as a function of $r$ alone, rather than in terms of individual hyperparameters.\footnote{This equivalence between the simplified procedure based on $r$ and the full minimization over $(\lambda,\psi)$ is an empirical feature of our data and model. In other contexts, or with different architectures and datasets, the correspondence might be weaker, and the full hyperparameter search could be necessary.}
We start our analysis by focusing on one-month-ahead forecasts of the unemployment rate. Later we will discuss robustness for other variables and horizons.
In Table (ref), we report the average loss over the evaluation sample for each level of complexity, which we expressed in terms of the implied variance ratio. A variance ratio of zero implies full shrinkage (constant prediction), while a ratio of one corresponds to no shrinkage and maximum complexity. Losses are reported relative to the unconditional quantile benchmark, so the top row equals one by construction.
This highlights the fundamental bias-variance trade-off in statistical learning. More complexity allows the model to capture more structure from the data but increases variance. More shrinkage stabilizes predictions but may lead to bias. The optimal point minimizes expected loss by balancing these effects.
For lower quantiles, the optimal complexity ratio around \(0.0{-}0.2\); for upper quantiles, it is closer to \(0.2{-}0.4\). The loss is lower for upper quantiles, suggesting greater predictability in the upper tail (i.e., rising unemployment). The need for more complexity in this region suggests that the predictors contain useful information. Conversely, for lower quantiles, optimal performance is achieved with limited complexity, suggesting limited information in the predictors for this region. In this case, variance reduction dominates.
For the upper quantiles, the model delivers clear gains relative to the unconditional quantile, thanks to its ability to extract predictive signals from the large information set through a flexible architecture. The more predictive content there is in the data, the broader the range of complexity ratios for which the model outperforms the naïve benchmark. However, even in cases of strong predictability, the danger of overfitting is always present. For instance, when the variance ratio exceeds 0.8, predictive losses increase sharply, and all forecasts deteriorate below the performance of the unconditional model. This underscores the critical role of out-of-sample validation in selecting the hyperparameters that govern model complexity.
Optimal complexity is fund to correspond to a range of variance ratios around 0.1 to 0.3 across quantiles.
We compare the performances in the test sample in order to study the stability of the performances between the validation and the testing sample. Comparing table for validation and testing we see very similar patters, in particular similar improvements over relevant ranges. This highlights the stability of the model performances across the two sample. Specifically, we see that in the test period the optimal complexity ratio is in a similar range (0.1 to 0.3)
This is further clarified by looking at predictions over time. We report predicted quantiles in Figure (ref), alongside the actual unemployment rate. The chart shows forecasts over the entire out-of-sample period, including both validation and testing. The vertical dotted line separates the two. Recall that forecasts in the testing period are fully real-time, since parameter estimates are recursive and hyperparameters are held fixed at the values selected during the validation period.
The central panel shows forecasts based on the architecture and hyperparameters selected in the validation sample. Forecasts in the validation period are not truly out-of-sample because hyperparameters are optimized using that sample. Thus, some overfitting is possible. However, because the number of hyperparameters is small relative to the sample size, the risk is limited.
To illustrate the role of complexity, we also show two extremes. The left panel imposes infinite shrinkage (variance ratio = 0), collapsing predictions to the unconditional quantiles. The right panel applies no shrinkage (variance ratio = 1), resulting in excessively volatile forecasts, indicating severe overfitting.
These three panels clearly display the bias-variance trade-off. At one extreme, excessive shrinkage eliminates all information, yielding stable but biased forecasts. At the other extreme, the absence of restrictions allow to incorporates all information exploiting the richeness of the unconstraned model, yielding low bias but high variance. The middle panel strikes an optimal balance. Since validation is done out-of-sample, this trade-off is clearly observable.
The middle panel shows that the model captures rising downside risk when it emerges. It also reflects business cycle asymmetry: downside risks increase, while upside risks (lower quantiles of unemployment) remain relatively flat. Notably, the similarity between forecasts in validation and testing periods indicates that overfitting due to hyperparameter selection is limited. This highlights the value of using out-of-sample performance as the criterion for model selection---a core innovation of machine learning.
In summary, the model detects time-varying labor market risk, especially on the downside. It identifies emerging vulnerabilities using a large set of macro indicators.
Importantly, all forecasts are generated in a fully automated, hands-off-the-wheel fashion. No manual selection of predictors or model tuning is used. All decisions---architecture, shrinkage, and learning rate---are guided by out-of-sample validation. Remarkably, this approach reproduces key insights from earlier studies AdrianEtal2019,kiley,boyarchenko2023outlook, but with no handcrafted inputs. Previous studies used two-step procedures with expert-selected factors; here, all features are learned from data.
Through careful regularization and model selection, the DNN handles a very flexible architecture without overfitting. We attribute this robustness to the combination of shrinkage and the DNN's ability to operate on low-dimensional manifolds, as discussed in GhigliazzaEtal2020.
We will now study the role of nonlinearities, and then assess robustness by considering different horizons and target variables. Specifically, we examine one-month-ahead unemployment, twelve-month-ahead unemployment, and one-month-ahead industrial production growth. These exercises demonstrate the flexibility and robustness of quantile estimation via DNNs.
In general, the model allows for complex non-linear interactions. However, one should check that a linear model is insufficient to match the problem's complexity, which, in this case is a penalized linear quantile regression, albeit over-parametrized.
In the left panel of Figure (ref), we report the predictions from the full model, where the slope of the Leaky ReLU is treated as a hyperparameter and chosen (along with others) to optimize performance in the validation sample. This is the same model as the one shown in the middle panel of Figure (ref). The linear activation case (\( \alpha = 1 \)) is among the candidate values included in this optimization.
In the right panel, we report the corresponding predictions when we impose a linear activation function, holding all other settings fixed. It is evident from the chart that the predictions are nearly identical. This suggests that, in this case, nonlinearities do not play a substantial role in improving predictive accuracy.
It is important to emphasize that our architecture allows for flexible nonlinearities, and we select among them using out-of-sample validation. Therefore, the finding that linearity performs just as well is a genuine result, not a consequence of model mis-specification or limitation.
In Table (ref) we report the losses of the linear model in the test and validation sample. For the convenience of the reader we report again the results for the full model, to ease comparability.
We conclude that the model is able to detect macroeconomic vulnerabilities as they emerge. However, in this setting, the empirical evidence indicates that nonlinear transformations are not essential for forecasting accuracy.
Forecasts based on the pinball loss function are generally more robust than those based on the quadratic loss. This difference stems from how these loss functions treat extreme observations. The quadratic loss assigns disproportionate weight to large errors, making it highly sensitive to outliers. In contrast, the pinball loss increases linearly with the error, thus downweighting extreme values and enhancing robustness.
This distinction is well-known in the context of estimating location parameters: the median (which minimizes the pinball loss for \( \tau = 0.5 \)) is more robust than the mean (which minimizes the quadratic loss). The left panel of Figure (ref) illustrates this point. For a given observed outcome, large deviations from the predicted value have a much larger impact on the quadratic loss than on the pinball loss.
The right panel of Figure (ref) provides supporting evidence from our empirical application. After the COVID-19 outbreak, predictions based on the pinball loss remain stable and plausible, while predictions based on the quadratic loss become erratic and unstable. This underscores the empirical value of quantile-based forecasting, particularly in turbulent periods.
The baseline exercise has focused on forecasting unemployment one month ahead. Interestingly, the main conclusions are valid for longer horizons and other target variables. We explored many combinations and found that these insights consistently hold. Comprehensive results are omitted for brevity but are available upon request. As illustrative examples, we report forecasts for the year-over-year change in unemployment one year ahead (Figure (ref)).
The model captures upside macroeconomic risk at longer horizons, detecting predictable movements in the upper quantiles of unemployment. This confirms that the general approach is effective even when forecasting risk further into the future.
We also report the forecast of the monthly growth rate of industrial production one month ahead (Figure (ref)).
The model detects downside macroeconomic risk by identifying predictable declines in the lower quantiles of industrial production growth.
In both cases, the forecasts generated by the full DNN and those with linear activation functions are nearly identical. The generality of this empirical result echoes CarrieroEtal2024, who use Time-Series Large Language Models (TS-LLMs) for macroeconomic forecasting with Big Data. These models adapt transformer-based large language models, originally developed for text, to time-series forecasting. Using a similar dataset, they find that TS-LLMs do not outperform linear factor models, which are closely related to principal components DozEtal2012, or large Bayesian VARs, which are closely related to linear ridge regression BanburaEtal2010, BanburaEtal2015.
In summary, the results above show that the model is able to capture business cycle vulnerabilities as they emerge by exploiting the information contained in Big Data through a highly flexible and general predictive model. There are no major variations in upside macroeconomic risk, while substantial variation is observed on the downside. Methodologically, our conclusion is that machine learning models, when properly regularized and validated, can robustly uncover predictive structures in high-dimensional macroeconomic data, even when their functional complexity is not ultimately required. We have shown that using out-of-sample performance to guide architecture and hyperparameter selection is essential. Otherwise, the signal in the data would be lost in the noise due to overfitting. Thanks to this validation-based selection process, the complexity of the DNN does not lead to overfitting. In this context, shrinkage is crucial. Allowing for nonlinearities does not provide substantial gains in predictive accuracy, in this case.
Throughout the paper, we focused on partial models, in which the objective was to predict one variable at a time. The natural extension is to full models, where all variables are predicted jointly. In the linear setting, the literature has built up developing large Bayesian vector autoregressions (VARs) that leverage regularization techniques with a direct connection to informative prior specifications, and shown that they deliver accurate point and density forecasts while offering both computational tractability and theoretical grounding for inference in high dimension BanburaEtal2010,Koop2013JAE,GiannoneEtal2015. Related contributions incorporated stochastic volatility to further enhance density forecasting performance Clark2011JBES,Chan2020JBES. In the more general setting of this paper, however, the extension was far from trivial: the multivariate counterpart of quantile predictions posed conceptual challenges, since there is no natural ordering in multiple dimensions, and the simulation of nonlinear models to recover longer horizon forecasts required particular care. Recent work by AdrianEtal2019 explored nonparametric multivariate distribution regression models, but their approach was limited to a handful of variables and did not easily extend to high-dimensional settings. We therefore leave the development of such high-dimensional nonlinear multivariate models for future research. At the same time, the results of this paper provided a useful guide in that direction: in particular, the finding that linear models performed well at the quantile level offered a natural starting point for developing tractable high-dimensional multivariate frameworks. Such extensions would be crucial for advancing risk assessment and policy analysis, as understanding the joint distribution of many variables is essential for monitoring vulnerabilities and the buildup of risks in the economy.
While this paper focuses on forecasting macroeconomic risks in the United States, the framework developed here—evaluating the role of regularization through shrinkage, model architecture choices such as depth and activation functions, and the use of out-of-sample validation—can be applied more broadly. Future work could extend this analysis to other economic domains, including finance, microeconomic applications, and international macroeconomics, as well as to different data structures such as cross-sectional and panel data, in the spirit of GiannoneEtal2021.