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.
115,499 characters · 17 sections · 78 citation commands
Real-time Inflation Forecasting Using Non-linear Dimension Reduction Techniques
\thispagestyle{empty}
\onehalfspacing
Inflation expectations are used as crucial inputs for economic decision making in central banks such as the European Central Bank (ECB) and the US Federal Reserve (Fed). Given current and expected inflation, economic agents decide on how much to consume, save and invest. In addition, measures of inflation expectations are often employed to estimate the slope of the Phillips curve, infer the output gap or the natural rate of interest. Hence, being able to accurately predict inflation is key for designing and implementing appropriate monetary policies in a forward looking manner.
Although the literature on modeling inflation is voluminous and the efforts invested considerable, predicting inflation remains a difficult task and simple univariate models are still difficult to beat sw07inflation. The recent literature, however, has shown that using large datasets stock2002macroeconomic and/or sophisticated models koop2007estimation, kk2012, giannone2013, koop2013var, clark2015macroeconomic, CCK2018, JaroLenza2018 has the potential to improve upon simpler benchmarks.
These studies often exploit information from huge datasets. This is commonly achieved by extracting a relatively small number of principal components (PCs) and including them in a second stage regression model stock2002macroeconomic. While this approach performs well empirically and yields consistent estimators for the latent factors, it fails to capture non-linear relations in the dataset. In the presence of non-linearities, using simple PCs potentially reduces predictive accuracy by ignoring important features of the data. Some studies deal with this issue by using flexible factor models which allow for non-linearities in the data. Bai-Ng2008 use targeted predictors coupled with quadratic principal components and show that allowing for non-linearities yields non-trivial improvements in predictive accuracy for inflation. This suggests that non-linearities (of a known form) are present in US macroeconomic datasets which are commonly employed for inflation forecasting. More recently, pelger2021state propose a flexible state-dependent factor model and apply this method to US bond yields and stock returns. Using this non-linear and non-parametric technique yields results which differ from linear, PC-based models by extracting significantly more information from the data.
One additional assumption commonly made is that the relationship between inflation and the latent factors is constant. For longer time series which feature multiple structural breaks this assumption is a strong one and may be deleterious for predictive accuracy. Several recent papers deal with this issue by using time-varying parameter (TVP) regressions which, in addition, allow for heteroscedasticity through stochastic volatility (SV) models koop2007estimation, kk2012, giannone2013, belmonte2014hierarchical, clark2015macroeconomic, JaroLenza2018, korobilis2019high.
Investigating whether allowing for non-linearities in the compression stage pays off for inflation forecasting is the key objective of the present paper. Building on recent advances in machine learning gallant-white1992, McNelis2004Forecasting, exterkate2016nonlinear, Chakraborty2017ML, heaton2017, Mullainathan2017ML, Polson2018, Coulombe2019ML, Kelly2018AE, medeiros2019forecasting, we adopt several non-linear dimension reduction techniques. The resulting latent factors are then linked to inflation in a second stage regression. To investigate whether there exists a relationship between non-linear factor estimation and flexible modeling of the predictive inflation equation, we introduce dynamic regression models that allow for TVPs and SV. Since the inclusion of a relatively large number of latent factors can still imply a considerable number of parameters (and this problem is even more severe in the TVP regression case), we rely on state-of-the-art shrinkage techniques.
From an empirical standpoint it is necessary to investigate how these dimension reduction techniques perform over time and during different business cycle phases. We show this by carrying out a thorough real-time forecasting experiment for the US. Our forecasting application uses monthly real-time datasets mccracken2016fred and includes a battery of well established models commonly used in central banks and other policy institutions to forecast inflation. These include simple benchmarks as well as more elaborate models such as the specification proposed in stock2002macroeconomic.
Our results show that non-linear dimension reduction techniques yield forecasts that are highly competitive to (and in fact often better than) the ones obtained from using linear methods based on PCs. In terms of one-month-ahead forecasts we find that models based on the Autoencoder yield point and density forecasts which are more precise than the ones obtained from other sophisticated non-linear dimension reduction techniques as well as traditional methods based on PCs. When the focus is on one-quarter-ahead forecasts we find that non-linear variants of PCs perform best. This performance, however, is not homogeneous over time and some of the models do better than others during different stages of the business cycle. In a brief discussion, we also analyze how our set of models performs during the COVID-19 pandemic.
These findings give rise to the second contribution of our paper. Since we observe that more sophisticated non-linear dimension reduction methods outperform simpler techniques during recessions, we combine the different models using dynamic model averaging RafteryDMA, koop2013var. We show that combining our proposed set of models with a variety of standard forecasting models yields predictive densities which are very close to the single best performing model in overall terms. Since the set of models we consider is huge, this indicates that using model and forecast averaging successfully controls for model uncertainty.
The remainder of this paper is structured as follows. Section (ref) discusses our proposed set of dimension reduction techniques. Section (ref) introduces the econometric modeling environment that we use to forecast inflation. Section (ref) first provides some in-sample features, then discusses the results of the forecasting horse race and finally presents our findings based on forecast averaging. The last section summarizes and concludes the paper. The Online Appendix provides further details on the econometric techniques as well as the data and additional empirical results.
Suppose that we are interested in predicting inflation using a large number of $K$ regressors that we store in a $T \times K$ matrix $\bm X = (\bm x_1, \dots, \bm x_T)'$, where $\bm x_t$ denotes a $K$-dimensional vector of observations at time $t$. If $K$ is large relative to $T$, estimation of an unrestricted model that uses all columns in $\bm X$ quickly becomes cumbersome and overfitting issues arise. As a solution, dimension reduction techniques are commonly employed stock2002macroeconomic, bernanke2005measuring. These methods strike a balance between model fit and parsimony. At a very general level, the key idea is to introduce a function $f$ that takes the matrix $\bm X$ as input and yields a lower dimensional representation $\bm Z = f(\bm X) = (\bm z_1, \dots, \bm z_T)'$, which is of dimension $T \times q$, as output. The critical assumption to achieve parsimony is that $q \ll K$. The latent factors in $\bm Z$ are then linked to inflation through a dynamic regression model (see Section (ref)).
The function $f: \mathbb{R}^{T \times K} \to \mathbb{R}^{T \times q}$ is typically assumed to be linear with the most prominent example being PCs. In this paper, we will consider several choices of $f$ that range from linear to highly non-linear (such as manifold learning as well as deep learning) specifications. We subsequently analyze how these different specifications impact inflation forecasting accuracy. In the following sub-sections, we briefly discuss the different techniques and refer to the original papers for additional information.
We start our discussion by considering principal component analysis (PCA). Minor alterations of the standard PCA approach allow for introducing non-linearities in two ways. First, we can introduce a non-linear function $g$ that maps the covariates onto a matrix $\bm W = g(\bm X)$. Second, we could alter the sample covariance matrix (the kernel) with a function $h$: $\bm \kappa = h(\bm W' \bm W)$. Both $\bm W$ and $\bm \kappa$ form the two main ingredients of a general PCA reducing the dimension to $q$, as outlined below Schoelkopf1998.
Independent of the functional form of $g$ and $h$, we obtain PCs by performing a truncated singular value decomposition (SVD) of the transformed sample covariance matrix $\bm \kappa$. Conditional on the first $q$ eigenvalues, the resulting factor matrix $\bm Z$ is of dimension $T \times q$. These PCs, for appropriate $q$, explain the vast majority of variation in $\bm X$. In the following, the relationship between the PCs and $\bm X$ is:
with $\bm \Lambda(\bm \kappa)$ being the truncated $K \times q$ eigenvector matrix of $\bm \kappa$ stock2002macroeconomic. Notice that this is always conditional on deciding on a suitable number $q$ of PCs. The number of factors is a crucial parameter that strongly influences predictive accuracy and inference Bai-Ng-2002. In our empirical work, we consider a small ($q=5$), moderate ($q=15$), and large ($q=30$) number of PCs.
By varying the functional form of $g$ and $h$ we are now able to discuss the first set of linear and non-linear dimension reduction techniques belonging to the class of PCA:
Diffusion maps, originally proposed in coifman2005diff and coifman2006diffusion, are another set of non-linear dimension reduction techniques that retain local interactions between data points in the presence of substantial non-linearities in the data.\footnote{For an application to astronomical spectra, see richards2009diff.} The local interactions are preserved by introducing a random walk process.
The random walk captures the notion that moving between similar data points is more likely than moving to points which are less similar. We assume that the weight function which determines the strength of the relationship between $\bm x_{\bullet i}$ to $\bm x_{\bullet j}$ is given by
where $||\bm x_{\bullet i} - \bm x_{\bullet j}||$ denotes the Euclidean distance between $\bm x_{\bullet i}$ and $\bm x_{\bullet j}$ and $c_2$ is a tuning parameter set such that ${w}(\bm x_{\bullet i}, \bm x_{\bullet j})$ is close to zero except for $\bm x_{\bullet i} \approx \bm x_{\bullet j}$. Here, $c_2$ is determined by the median distance of the $k$-nearest neighbors of $\bm x_{\bullet i}$ as suggested by Zelnik2005diff. The number of $k$ is approximated using the algorithm suggested by angerer2016destiny.
The probability of moving from $\bm x_{\bullet i}$ to $\bm x_{\bullet j}$ is then simply obtained by normalizing:
This probability tends to be small except for the situation where $\bm x_{\bullet i}$ and $\bm x_{\bullet j}$ are similar to each other. As a result, the probability that the random walk moves from $\bm x_{\bullet i}$ to $\bm x_{\bullet j}$ will be large if they are equal but rather small if both covariates differ strongly.
Let $\bm P$ denote a transition matrix of dimension $K \times K$ with $(i,j)$th element given by $p_{i \to j}$. The probability of moving from $\bm x_{\bullet i}$ to $\bm x_{\bullet j}$ in $n=1,2, \dots$ steps is then simply the matrix power of $\bm P^n$, with typical element denoted by $p_{i \to j}^n$. Using a biorthogonal spectral decomposition of $\bm P^n$ yields:
with $\psi_s$ and $\phi_s$ denoting left and right eigenvectors of $\bm P$, respectively. The corresponding eigenvalues are given by $\lambda_s$.
We then proceed by computing the so-called diffusion distance as follows:
with $p_0$ being a normalizing factor that measures the proportion the random walk spends at $\bm x_{\bullet j}$. This measure turns out to be robust with respect to noise and outliers. coifman2006diffusion show that
This allows us to introduce the family of diffusion maps from $\mathbb{R}^{K} \to \mathbb{R}^{q}$ given by:
The distance matrix can then be approximated as:
Intuitively, this equation states that we now approximate diffusion distances in $\mathbb{R}^K$ through the Euclidian distance between $\bm \Xi_n(\bm x_{\bullet i})$ and $\bm \Xi_n (\bm x_{\bullet j})$. This discussion implies that we have to choose $n$ and $q$ and we do this by setting $q=\{5,15,30\}$ according to our approach with either a small, moderate or large number of factors and $n=T$, the number of time periods. The algorithm in our application is implemented using the R packages diffusionMap and destiny diffusionMap, angerer2016destiny.
Locally linear embeddings (LLE) have been introduced by roweis2000lle. Intuitively, the LLE algorithm maps a high dimensional input dataset $\bm X$ into a lower dimensional space while preserving the neighborhood structure. This implies that points which are close to each other in the original space are also close to each other in the transformed space.
The LLE algorithm is based on the assumption that each $\bm x_{\bullet i}$ is sampled from some underlying manifold. If this manifold is well defined, each $\bm x_{\bullet i}$ and its neighbors $\bm x_{\bullet j}$ are located close to a locally linear patch of this manifold. One consequence is that each $\bm x_{\bullet i}$ can be reconstructed from its neighbors $\bm x_{\bullet j}$ with $j \neq i$, conditional on suitably chosen linear coefficients. This reconstruction, however, will be corrupted by measurement errors. roweis2000lle introduce a cost function to quantify these errors:
with $\omega_{ij}$ denoting the $(i,j)$th element of a weight matrix $\bm \Omega$. This cost function is then minimized subject to the constraint that each $\bm x_{\bullet i}$ is reconstructed only from its neighbors. This implies that $\omega_{ij} = 0$ if $\bm x_{\bullet j}$ is not a neighbor of $\bm x_{\bullet i}$. The second constraint is that the matrix $\bm \Omega$ is row-stochastic, i.e., the rows sum to one. Conditional on these two restrictions, the cost function can be minimized by solving a least squares problem.
To make this algorithm operational we need to define our notion of neighbors. In the following, we will use the $k$-nearest neighbors in terms of the Euclidean distance. We choose the number of neighbors by applying the algorithm proposed by kayo2006, which automatically determines the optimal number for $k$. The $q$ latent factors in $\bm Z$, with typical $i$th column $\bm z_{\bullet i}$, are then obtained by minimizing:
which implies a quadratic form in $\bm z_t$. Subject to suitable constraints, this problem can be easily solved by computing:
and finding the $q+1$ eigenvectors of $\bm M$ associated with the $q+1$ smallest eigenvalues. The bottom eigenvector is then discarded to arrive at $q$ factors. For our application, we use the R package lle LLE.
Isometric Feature Mapping (ISOMAP) is one of the earliest methods developed in the category of manifold learning algorithms. Introduced by tenenbaum2000, the ISOMAP algorithm determines the geodesic distance on the manifold and uses multidimensional scaling to come up with a low number of factors describing the underlying dataset. Originally, ISOMAP was constructed for applications in visual perception and image recognition. In economics and finance, some recent papers highlight its usefulness Ribeiro2008, lin2011isomap, orsenigo2013isomap, Zime2014.
The algorithm consists of three steps. In the first step, a dissimilarity index that measures the distance between data points is computed. These distances are then used to identify neighboring points on the manifold. In the second step, the algorithm estimates the geodesic distance between the data points as shortest path distances. In the third step, metric scaling is performed by applying classical multidimensional scaling (MDS) to the matrix of distances. For the dissimilarity transformation, we determine the distance between point $i$ and $j$ by the Manhattan index $d_{ij} = \sum_k |x_{ki} - x_{kj}|$ and collect those points where $i$ is one of the $k$-nearest neighbors of $j$ in a dissimilarity matrix. For our empirical application, we again choose the number of neighbors by applying the algorithm proposed by kayo2006 and use the implementation in the R package vegan vegan.
The described non-linear transformation of the dataset enables the identification of a non-linear structure hidden in a high-dimensional dataset and maps it to a lower dimension. Instead of pairwise Euclidean distances, ISOMAP uses the geodesic distances on the manifold and compresses information under consideration of the global structure.
Deep learning algorithms are characterized by not only non-linearly converting input to output but also representing the input itself in a transformed way. This is called representation learning in the sense that representations of the data are expressed in terms of other, simpler representations before mapping the data input to output values.
One tool which performs representation of itself as well as representation to output is the Autoencoder (AE). The first step is accomplished by the encoder function, which maps an input to an internal representation. The second part, which maps the encoded (transformed) data to the output, is called the decoder function. Their ability to extract factors, which explain a large fraction of the variability in the observed data, in a non-linear manner makes deep learners a powerful tool complementing the range of commonly used dimension reduction techniques Goodfellow-et-al-2016. andreini2020deep, for example, embed a dynamic Autoencoder structure in a dynamic factor model and show that it yields a good now- and forecasting performance for US GDP. In their paper, they allow for additional flexibility by simultaneously estimating the non-linear latent factors and the parameters. In empirical finance, heaton2017, Polson2018 and Kelly2018AE find that the application of these methods is beneficial to predict asset returns.
Based on deep learning techniques, we propose obtaining hierarchical predictors $\bm Z$ by applying a number of $l \in \{1, \dots ,L\}$ non-linear transformations to $\bm X$. These transformations are called hidden layers with $L$ giving the depth of our architecture and $f$ denoting an univariate activation function.\footnote{In principle, $f$ can vary over the different layers.} More specifically, in each layer, activation functions (non-linearly) transform the inputs (which are the outputs of the previous layer). A common choice, which we adopt, is the hyperbolic tangent (tanh) given by
We apply this function element-wise to the entries of $\bm X$. Using tanh activation functions is justified by its strong empirical properties identified in recent studies such as Saxe2019deep and andreini2020deep.
The structure of our deep learning algorithm can be represented in form of a composition of univariate semi-affine functions given by
and $\hat{\bm X}^{(0)} = \bm X$ for $l=0$. Here, $\bm W^{(l)}$ denotes a weighting matrix of dimension $N_{l-1} \times N_{l}$ (with $N_l$ being the number of neurons in layer $l$), $\bm b_{l}$ is a $N_l \times 1$ bias vector and $\bm \iota_T$ is a $T \times 1$ vector of ones.
The output of the network is then obtained by setting:
Notice that if we set $N_L=q (\ll K)$, we achieve dimension reduction and the output of the network is a (non-linearily) compressed version of the input dataset. In principle, what we have just described constitutes the encoding part of the Autoencoder. If we are interested in recovering the original dataset $\bm X$ we simply have to add additional layers characterized by increasing numbers of neurons until we reach $N_{L+j} = K$ for $j=1,2,\dots$.
The optimal sets of $\hat{\bm W} = (\hat{\bm W}^{(1)},\dots,\hat{\bm W}^{(L)})$ and $\hat{\bm b} = (\hat{\bm b}_1,\dots,\hat{\bm b}_L)$ are obtained by computing a loss function, most commonly the mean squared error of the in-sample fit. The complexity of the neural network is determined by choosing the number of hidden layers $L$ and the number of neurons in each layer $N_l$. We perform our forecasting exercise with different sets of tuning parameters and choose one, three, five, and eight hidden layers with the number of neurons evenly being downsized to the desired number of factors.
For the loss function and the optimization algorithm we stick to common choices in the literature and use the mean squared error loss function and the Adaptive Moment Estimation (ADAM). We repeat the optimization procedure in $100$ epochs on at least $84$ batches which corresponds to the average duration of a business cycle in the US.\footnote{The average duration of a business cycle was determined using data provided by The National Bureau of Economic Research on business cycle expansions and recessions.} This implies that we train the algorithm in each epoch with a partition of the original data set of at least the length of one business cycle. To capture the dynamics of the different cycles present in the data the optimization procedure needs to be repeated in a reasonably high number of epochs. We find that the algorithm converges quickly and setting the number of epochs to $100$ is sufficient.
We employ the R interface to keras keras, a high-level neural networks API and widely used package for implementing deep learning models.
In the following, we introduce the predictive regression that links our target variable, inflation in consumer prices, to $\bm Z$ and other observed factors. Following stock1999forecasting, inflation is specified such that:
with $\text{CPI}_{t+h}$ denoting the consumer price index in period $t+h$.
In the empirical application we set $h \in \{1, 3\}$. $y_{t+h}$ is then modeled using a dynamic regression model:
where $\bm \beta_{t+h}$ is a vector of TVPs associated with $M (=q + p)$ covariates denoted by $\bm d_t$ and $\sigma^2_{t+h}$ is a time-varying error variance. $\bm d_t$ might include the latent factors extracted from the various methods discussed in the previous sub-section, lags of inflation, an intercept term or other covariates which are not compressed.
Following much of the literature Taylor1982stochvol, bkk, kg2014, kastner2014ancillarity, stock2016core, chan2017stochastic, huber2020 we assume that the TVPs and the error variances evolve according to independent stochastic processes:
with $\mu_h$ denoting the conditional mean of the log-volatility, $\rho_h$ its persistence parameter and $\vartheta_h^2$ the error variance of $\log \sigma_{t+h}^2$. The matrix $\bm V$ is an $M \times M$-dimensional variance-covariance matrix with $\bm V = \text{diag}(v_1^2, \dots, v_M^2)$ and $v_j^2$ being the process innovation variance that determines the amount of time-variation in $\bm \beta_{t+h}$. This setup implies that the TVPs are assumed to follow a random walk process while the log-volatilities evolve according to an AR(1) process.
The model described by Eq. ((ref)) and Eq. ((ref)) is a flexible state space model that encompasses a wide range of models commonly used for forecasting inflation. For instance, if we set $\bm V = \bm 0_M$ and $\vartheta^2=0$, we obtain a constant parameter model with homoscedastic errors. If $\bm V$ is instead a full $M \times M$ matrix but of reduced-rank, we obtain the model proposed in chan2020reducing. If $\bm d_t$ includes the lags of inflation and (lagged) PCs, we obtain a model closely related to the one used in stock2002macroeconomic. If we set $d_t =1$ and allow for TVPs, we obtain a specification similar to the unobserved components stochastic volatility model successfully adopted in stock1999forecasting. A plethora of other models can be identified by appropriately choosing $\bm d_t$, $\bm V$ and $\vartheta^2$. This flexibility, however, calls for model selection. We select appropriate submodels by using Bayesian methods for estimation and forecasting. These techniques are further discussed in Section (ref) of the Online Appendix and allow for data-based shrinkage towards simpler nested alternatives.
In our empirical application we consider the popular FRED-MD database. This dataset is publicly accessible and available in real-time. The monthly data vintages ensure that we only use information that would have been available at the time a given forecast is being produced. A detailed description of the databases can be found in mccracken2016fred. To achieve approximate stationarity we transform the dataset as outlined in Section (ref) of the Online Appendix. Furthermore, each time series is standardized to have sample mean zero and unit sample variance prior to using the non-linear dimension reduction techniques.
Our US dataset includes 105 monthly variables that span the period from 1963:01 to 2021:01. The forecasting design relies on a rolling window, as justified in clark2011, that initially ranges from 1980:01 to 1999:12. For each month of the hold-out sample, which starts in 2000:01 and ends in 2019:12, we compute the $h$-month-ahead predictive distribution for each model (for $h \in \{1, 3\}$), keeping the length of the estimation sample fixed at $240$ observations (i.e., a rolling window of $20$ years).\footnote{In addition to our baseline sample ending in 2019:12, we present the results of our forecasting exercise including observations covering the COVID-19 pandemic (2020:01 to 2020:08) in Sub-section (ref). Since the pandemic caused severe outliers in our dataset, including those periods helps to test the forecasting performance of our models during turbulent times.} For these periods we contrast each forecast with the realization of inflation in the vintage one-quarter-ahead, following the evaluation approach of chan2017stochastic. As most data revisions take place in the first quarter while afterwards the vintages remain relatively unchanged croushore2011frontiers, pfarrhofer2020forecasts, we make sure that realized inflation is not subject to revisions anymore.
One key limitation is that all methods are specified conditionally on $\bm d_t$ and thus implicitly on the specific function $f$ used to move from $\bm X$ to $\bm Z$. Another key objective of this paper is to control for uncertainty with respect to $f$ by using dynamic model averaging techniques. For obtaining predictive combinations, we use the first $24$ observations of our hold-out sample. The remaining periods (i.e., ranging from 2002:01 to 2019:12) then constitute our evaluation sample and the respective predictions are again contrasted to the one-quarter-ahead vintage of inflation.
In terms of competing models we can classify the specifications along two dimensions:
In this sub-section we analyze bivariate correlations between the factors, obtained from using different dimension reduction techniques, and the variables in our dataset as well as inflation. These correlations provide some information on the specific factor dynamics and (with caution) on how to interpret the factors in $\bm Z$ from a structural perspective.\footnote{The estimates of the factor are considerably more difficult to interpret. Nevertheless, to provide some intuition on how the factors for the best performing specifications evolve over time, see (ref) in the Online Appendix.} The recent literature crawford2018approx, crawford2019approx, joseph2019parametric advocates using linear approximations or Shapley values to improve interpretability of these highly non-linear models. In this paper, we opt for a simple correlation-based approach given the large amount of competing dimension reduction techniques and the fact that for some of these the different techniques work better than for other methods.
(ref) is a heatmap of the correlations with rows denoting the different covariates in $\bm X$ and columns representing the different dimension reduction techniques. These correlations are averages across the factors (in case that $q>1$) and, since we include several lags of the input dataset, are also averaged across the lags.
The figure suggests for most dimension reduction techniques that the factors are correlated with housing quantities (PERMIT and HOUST alongside their sub-components) as well as interest rate spreads. Some variables which measure real activity (such as industrial production and several of its components) also display comparatively large correlations with the factors. In some cases, these correlations are positive whereas in other cases, correlations are negative. In both instances, however, the absolute magnitudes are similar. The three exceptions from this rather general pattern are diffusion maps as well as PCA quadratic and squared. In this case, the corresponding columns indicate lower correlations.
Averaging over the factors, as done in (ref), potentially masks important features of individual factors. Next we ask whether there are relevant differences by analyzing the correlations between each $\bm z_j~(j=1, \dots, q)$ and each column of $\bm X$. For brevity, we focus on a specific model that performs extraordinarily well in terms of density forecasts: the Autoencoder with a single hidden layer and $30$ factors. (ref) shows, for each factor, the five variables which display the largest absolute correlation. The variables in the rows are a union over the sets of top-five variables for each factor. This figure shows that several factors display quite similar correlation patterns. For all of them, housing quantities are either positively or negatively correlated (with similar magnitudes). Apart from that, and in consistence with the findings discussed above, we observe that financial market variables (such as interest rate spreads) show up frequently for several factors. Only very few factors depart from this overall pattern. In the case of factors 9, 22, 23 and 24 we find low correlations with housing and much stronger correlations with financial markets. In fact, factor 9 is closely tracking the credit (BAAFFM) and term spreads (e.g., T10YFFM).
These two heatmaps provide a rough overview on what variables drive the factors. Next, we ask whether we can construct models based on including variables which display the strongest correlations with the factors. This approach can be interpreted as a simple selection device which takes non-linearities in the input dataset implicitly into account. Since the heatmap is based on full-sample results and we are interested in using these small-scale models for out-of-sample forecasting we compute the correlation for each point in our hold-out period. In summary, the variables which frequently show up across hold-out periods and dimension reduction techniques are:
These variables are also the ones which display high correlations to the factors in (ref) and are included in the large-scale ARX model. Here, it is worth stressing that there seems to be appreciable heterogeneity with respect to dimension reduction methods. Most of them generate factors that are highly correlated with real activity and housing measures as well as interest rates and other stock market variables. Interestingly, when we focus on the second group we observe that the factors arising from using PCA squared (and to a somewhat lesser extent PCA quadratic) are heavily related to labor market measures. Average correlations with prices (i.e., CUSR0000SA0L5) are small for most techniques (with PCA quadratic yielding the largest correlations of around $0.3-0.4$). Some methods also yield factors that are strongly correlated to money stocks and reserves (e.g., diffusion maps). (ref) of the Online Appendix provides a much more detailed picture on the precise variables used to build the small-scale models.
Next, we ask whether the factors are correlated to inflation. (ref) shows the correlation with inflation averaged across the number of factors for each dimension reduction techniques as well as the minimum and maximum value (across these factors) in parentheses. To assess whether these correlations differ over time, we divide our sample into expansionary and recessionary periods.\footnote{Recessions are defined by using the the business cycle classification of the National Bureau of Economic Research (NBER).} Since the COVID-19 pandemic marks an extraordinary period in our sample, we also compute the correlations for 2020 only and include it at the bottom of (ref).
For the full sample as well as during expansions, we find that the factors obtained from using the linear variants of PCA display comparatively higher correlations relative to the other dimension reduction techniques (with some of the factors featuring a correlation of close to $0.2$). In recessions and the pandemic, these correlations increase substantially to reach average correlations close to $0.3$ (with the factor displaying the maximum correlation being strongly related to inflation, with values of around $0.6$). The non-linear dimension reduction techniques yield strong correlations during turbulent times (i.e., recessions and the pandemic). This is not surprising since these methods tend to work well if there are strong deviations from linearity (which mostly occurs in recessions). Such a feature can be easily demonstrated by considering PCA squared. In normal times, the factors will be centered around zero and typically display little variation. But in recessions the link function implies that larger changes will dominate the shape of the factors and imply pronounced movements which could be helpful for predicting turning points in inflation.
We now consider point and density forecasting performance of the different models and dimension reduction techniques. The forecast performance is evaluated through log predictive likelihoods (LPLs) for density forecasts and root mean squared errors (RMSEs) for point forecasts. Superior models are those with high scores in terms of LPL and low values in terms of RMSE. We benchmark all models relative to the autoregressive (AR) model with constant parameters and the Minnesota prior. The first entry in the tables gives the actual value of the LPL (cumulated over the hold-out sample) with actual RMSEs in parentheses (averaged over the hold-out sample) for our benchmark model. The remaining entries are differences in LPLs with relative RMSEs in parentheses. We mark statistically significant results according to the diebold1995dmtest test at the one, five and ten percent significance levels with one, two and three asterisks, respectively.
Starting with the one-month-ahead horizon, (ref) depicts the inflation forecasting results. This table suggests that, in terms of density forecasts, using dimension reduction techniques (both linear and non-linear) improves predictions substantially. These improvements arise not only relative to the AR benchmark but also related to the large AR models with additional exogenous regressors. For some models, these improvements are sizable, irrespective of the regression specification (i.e., whether we use a constant parameter or a TVP model). Especially the Autoencoder with one and five layers sharply improves upon the benchmark (and all the remaining competitors) by large margins. Moreover, it yields statistically significant improvements at the one percent level. A similar story emerges when we focus on point forecasts. Non-linear dimension reductions help slightly. Relative RMSEs are smaller but close to one for most models. Again, the Autoencoder works well and yields RMSEs which are, across regression specifications, almost 20 percent lower than the ones from the benchmark. It is worth emphasizing that PCA squared also yields highly competitive point forecasts which are statistically significant according to the diebold1995dmtest test.
When we compare model performance across regression specifications and focus on the Minnesota-type priors, we find that constant parameter models work quite well if non-linear dimension reduction techniques such as the Autoencoder are adopted. With a single exception (diffusion maps), introducing TVPs does not pay off and yields density forecasts which are slightly more imprecise than the ones obtained from their time-invariant counterparts. If we use a horseshoe prior, this result somewhat reverses (with the caveat that the models coupled with the horseshoe sometimes yield weaker inflation forecasts than the benchmark). Here, we observe that introducing TVPs often improves log predictive likelihoods relative to the constant parameter model with the same prior.
Summing up this discussion, we observe that the Autoencoder yields favorable point and density forecasts, irrespective of the prior and regression specification chosen. This strong performance of the Autoencoder, however, depends on the number of layers as well as the number of factors. The literature huang2003learning, heaton2008introduction suggests that the number of hidden layers should increase with the complexity of the dataset. Our results, however, suggest the opposite. For a typical US macroeconomic dataset the forecast performance of the Autoencoder seems to be strongest when a single hidden layer coupled with a large number of factors is used.
Next, we inspect the longer forecast horizon in greater detail. (ref) depicts the forecast performance of all competitors for the one-quarter-ahead horizon. The table indicates that several non-linear dimension reduction techniques (most notably the Autoencoder, PCA quadratic and PCA squared) clearly outperform the autoregressive benchmark as well as models based on linear PCs. The improvements relative to linear PCs are sizable and statistically significant (especially for squared and quadratic PCs). For this horizon, the large ARX models also exhibit excellent forecasting properties.
Zooming into the different approaches to dimension reduction reveals that PCA quadratic with TVPs yields highest LPLs. Among the different dimension reduction techniques, both PCA quadratic and squared stand out and improve appreciably against all other dimension reduction techniques. The Autoencoder also provides density forecasts which are highly competitive.
When we consider point forecasts a similar picture arises. Here we find that models which do well in terms of LPLs also yield precise point forecasts. The single best performing specification, however, is PCA squared with five factors, constant parameters and a horseshoe prior. This model improves upon the AR benchmark by around 24 percent. Notice, however, that the same model but with TVPs also yields forecasts which are 24 percent more precise than the ones of the benchmark. This suggests that the horseshoe shrinks the TVPs close to zero and the corresponding point estimators are almost identical. Since the LPLs differ, the remaining time-variation in the coefficients mainly impacts the LPLs through the predictive variance.
(ref) provides a summary of the best performing models for the one-month and the one-quarter-ahead forecasts, respectively. Moreover, we also assess how a model which includes the five variables displaying the highest average correlations to the factors performs (see discussion in Sub-section (ref)). These models are labeled small ARX in the table. The results suggest that for one-month-ahead forecasts, replacing the latent factors arising from the Autoencoder with observed variables that display a high correlation to the factors does not pay off in terms of point and density forecasts. Across all regression specifications, the RMSEs are higher and the LPLs lower. When we consider one-quarter-ahead forecasts, however, this strategy seems to work much better. In this case, smaller models with covariates selected based on their correlations to the factors obtained by using PCA squared and PCA quadratic yield density predictions which are close to the ones obtained from exploiting all available data. Notice, however, that for point forecasts, the performance of the small models is inferior.
Before proceeding to the next sub-section we briefly discuss two important issues. First, it is worth stressing that the factors used in this forecasting exercise are extracted from the full set of variables in $\bm X$. In (ref) and (ref) of the Online Appendix we divide the dataset into slow- and fast-moving variables bernanke2005measuring and extract the latent factors from these partitioned datasets exclusively. The main results based on extracting the factors from the full dataset remain in place: for one-month-ahead forecasts we find the Autoencoder to perform particularly well whereas for one-quarter-ahead predictions PCA squared and quadratic yield accurate forecast densities.
Second, for one of our best performing models (the Autoencoder with one hidden layer) forecasting performance changes sharply when the number of factors is changed. This raises the question on how the relationship between the number of factors and forecast performance is. In (ref) in the Online Appendix we show two graphs that discuss how point and density forecasting performance change with the number of factors. In this exercise we find that the largest jumps in predictive accuracy is found when increasing the number of factors from 17 to 24 and again from 29 to 30 in terms of LPLs and from 18 to 26 and 29 to 30 in terms of RMSE.
The results based on RMSEs and LPLs provide information on relative forecasting performance. In the next step, we ask whether the different methods and models we propose yield predictive distribution which are better calibrated. To this end, we consider the normalized forecast errors obtained through the probability integral transform (PIT). If a model is correctly specified the PITs are iid uniformly distributed and the respective standardized forecast errors should be iid normally distributed. Departures from the standard Gaussian distribution allow us to inspect along what dimensions the model is poorly calibrated. For instance, if the variance of the normalized forecast error is too small (i.e., below one) this is evidence that the predictive distribution is too wide (i.e., too many predictions are in the tails) while values greater than one indicate that the variance is too tight (i.e., the tails are not adequately represented).
(ref) shows the results for the one-month-ahead normalized forecast errors.\footnote{The results for one-quarter-ahead forecasts are provided in Section (ref) of the Online Appendix.} In principle, we observe that the mean across methods is close to zero (with some few exceptions such as PCA quadratic for $q=15$). Nevertheless, these differences are never statistically significantly different from zero. Considering the variances shows that most models yield forecast distributions which seem to be slightly too narrow (with variances exceeding one). The asterisks indicate whether the variances are significantly different from one. For some few models, this is the case (especially if we assume constancy of the parameters) but if we allow for TVPs there are only a handful of cases left. This, however, strongly depends on the shrinkage prior adopted. Turning to the autocorrelation of the normalized shocks reveals that these are mostly close to zero and never statistically significantly different from zero.
Comparing sophisticated to simple dimension reduction methods suggests no discernible differences in model calibration. In principle, approaches based on linear PCs yield normalized forecast errors with similar statistical properties than the ones obtained from using more sophisticated dimension reduction techniques.
The discussion above might mask important differences in calibration of different parts of the predictive distribution. We now turn to a deeper analysis of the one-month-ahead predictive distribution of the two best performing models vis-\'{a}-vis the benchmark: the Autoencoder 1l $(q=30)$ and PCA squared $(q=5)$. This analysis is based on visual inspection of the normalized forecast errors (left panel of (ref)), a histogram of the PITs (middle panel of (ref)) and the visual diagnostic of the empirical cumulative density function proposed in rossi2019alternative (right panel of (ref)). Recall that, under correct specification, the PITs should be iid uniformly and the normalized forecast errors should be iid standard normally distributed, respectively.
The left panel of the figure indicates that for both models under consideration, normalized forecast errors are centered on zero, display little serial correlation and a variance close to one (with the Autoencoder generating slightly more spread out normalized forecast errors). In some periods, normalized forecast errors depart significantly from the standard normal distribution (i.e., the corresponding observations lie outside the 95% confidence intervals). But in general, and for both models (and the benchmark), model calibration seems to be adequate. Next, we focus on the histogram in the middle panel of (ref) (which includes 95% confidence intervals). From this figure, we learn that both models are well calibrated with some tendency to overestimate the upper tail risk. Finally, considering the right panel shows that all models appear to be well calibrated, with most observations being clustered around the 45 degree lines and not a single observation being outside the 95% confidence intervals.
To the detriment of linear modeling techniques, the COVID-19 pandemic caused severe outliers for several of the time series we include in our dataset. Following the recent literature huber2020nowcasting, chkmp2021, coulombe2021covid19 which advocates using non-linear and non-parametric modeling techniques in turbulent times, we briefly investigate whether the non-linear dimension reduction techniques proposed in this paper yield more precise inflation forecasts during the pandemic.
(ref) depicts the differences in LPLs for the period 2020:01 to 2020:08. For illustrative purposes, we only consider the models with 30 factors.\footnote{The findings for the other factors are very similar and available from the corresponding author upon request.}
The figure provides a few interesting insights. First, we observe that in March 2020, models based on the Autoencoder improve upon the benchmark, irrespective of the prior and regression specification adopted. This finding is less pronounced for the other techniques in the constant parameter case. Comparing the performance of the constant parameter and the TVP regression models reveals that, irrespective of the prior, allowing for time variation in the parameters improves density forecasts during the pandemic. This finding is consistent with findings in, e.g., huber2020nowcasting, who show that flexible models improve upon linear models during the pandemic due to increases in the predictive variance.
In the previous sub-section and Section (ref) of the Online Appendix we provide some evidence that model performance varies considerably over time (see (ref)). The key implication is that non-linear compression techniques (and time-varying parameters) might be useful during turbulent times whereas forecast evidence is less pronounced in normal times. In this sub-section, we ask whether combining models in a dynamic manner further improves predictive accuracy.
After having obtained the predictive densities of $y_{t+h}$ for the different dimension reduction techniques and model specifications, the goal is to exploit the advantages of both linear and non-linear approaches. This is achieved by combining models in a model pool such that better performing models over certain periods receive larger weights while inferior models are subsequently down-weighted. The literature on forecast combinations suggests several different weighting schemes, ranging from simply averaging over all models hendry2004pooling, Hall2007combi, clark2010averaging, berg2015combination to estimating weights based on the models' performances according to the minimization of an objective or loss function timmermann2006combi, Hall2007combi, Geweke2011pools, conflitti2015combination, pettenuzzo2016combi or according to the posterior probabilities of the predictive densities RafteryDMA, kk2012, beckmann2020exchange. More recent approaches set up separate state space models which assume sophisticated law of motions for the weights associated with each predictive distribution billio2013combi, pettenuzzo2016combi, mcalinn2019bma. These approaches, while being elegant and having the advantage of incorporating all available sources of uncertainty (i.e., also control for estimation uncertainty in the weights), are computationally cumbersome if the number of models to be combined is large.
Since our model space is huge, we use computationally efficient approximations to dynamically combine models. Our approach builds on combining predictive densities according to their posterior probabilities. This is referred to as Bayesian model averaging (BMA). The resulting weights are capable of reflecting the predictive power of each model for the respective periods. Dynamic model averaging (DMA), as specified by RafteryDMA, extends the approach by adding a discount (or forgetting) factor to control for a model's forecasting performance in the recent past. The `recent past' is determined by the discount factor, with higher values attaching greater importance to past forecasting performances of the model and lower values gradually ignoring results of past predictive densities. Similar to RafteryDMA, kk2012 and beckmann2020exchange, we apply DMA to combine the predictive densities of our various models. These methods do not require computationally intensive MCMC or sequential Monte Carlo techniques and are thus fast and easy to implement.
DMA works as follows. Let $\bm \varrho_{t+h|t} = (\varrho_{t+h|t,1}, \dots, \varrho_{t+h|t,J})'$ denote a set of weights for $J$ competing models at time $t+h$ given all available information up to time $t$. These (horizon-specific) weights vary over time and depend on the recent predictive performance of the model according to:
where $p_j(y_{t+h}|y_{1:t})$ denotes the $h$-month-ahead predictive distribution of model $j$ evaluated at $y_{t+h}$ and $\delta \in (0, 1]$ denotes a forgetting factor close to one. Intuitively speaking, the first equation is a prediction of the weights based on all available information up to time $t$ while the second equation shows how the weights get updated if new data flows in.
In our empirical work we set $\delta = 0.97$.\footnote{koop2013var find robust results over the interval [0.95,1]. After optimizing over this set of parameter values we choose $\delta = 0.97$.} Notice that if $\delta=1$, we obtain standard BMA weights while $\delta=0$ would imply that the weights depend exclusively on the forecasting performance in the last period.
Weights obtained by combining models according to their predictive power convey useful information about the adequacy of each model over time. In order to get a comprehensive picture of the effects of different model modifications, we combine our models and model specifications in various ways.
(ref) presents the forecasting results when we use DMA to combine models. Again, all models are benchmarked to the AR model with constant parameters and the Minnesota prior. The first row depicts the relative performance of the single best performing model for the chosen time horizon.
The table can be understood as follows. Each entry includes all dimension reduction techniques. The rows define whether the model space includes all factors $q \in \{5, 15, 30\}$ or whether we combine models with a fixed number of factors exclusively. The columns refer to model spaces which include only constant parameter, time-varying parameter or both specifications in the respective model pool. Since we also discriminate between two competing priors we consider model weights conditioning on either the horseshoe or the Minnesota prior or average across both prior specifications (the first upper part of the table with $\{\text{HS, MIN}\}$).
Across the two forecast horizons considered, we find pronounced accuracy improvements for point and density forecasts relative to the AR model. When we benchmark the different combination strategies to the single best performing model we find no accuracy gains for both horizons. Differences in terms of LPLs are, however, rather small. This suggests that while the best performing model (i.e., a constant parameter regression with factors obtained through the Autoencoder) is hard to beat, one can effectively reduce model and specification uncertainty and thus obtain competitive forecasts without the need to rely on a single model.
Comparing whether restricting the model a priori improves predictions yields mixed insights. For the one-month-ahead horizon we observe that pooling over models which use our variant of the Minnesota prior yields more favorable forecasts as compared to a pooling strategy which uses both priors or the horseshoe only. When we pool over constant and TVP regressions we find small decreases in predictive accuracy relative to a model pool which only includes constant parameter regressions.
Turning to one-quarter-ahead forecasts yields a similar picture. Using a large pool of models generally leads to slightly less precise forecasts. For higher-order forecasts our results suggest that pooling models which use the horseshoe yields higher LPLs. When we compare the different regression specifications we find that integrating out uncertainty with respect to whether parameters should be time-varying yields forecasts which are very similar to the strategy that only pools over constant parameter models.
In general, the differences in predictive performance across the DMA-based averaging schemes are small. Hence, as a general suggestion we can recommend applying DMA and using the most exhaustive model space available (i.e., including both priors, the different number of factors and TVP and constant parameter regressions).
To investigate which model receives substantial posterior weight over time, (ref) depicts the weights associated with the one-month-ahead LPLs over the hold-out period. Panel (a) displays the weight placed on models that allow for TVP, panel (b) shows the weight attached to the different number of factors and panel (c) shows the weight attached to each model. These weights are obtained by using the full model space (i.e., that includes both priors, TVP and constant parameter regressions and all number of factors). The weight placed on TVP specifications, for instance, is then simply obtained by summing up the weights associated with the different models that feature TVPs.
Starting with the top panel of the figure, we observe that during the beginning of the sample, appreciable model weight is placed on constant parameter models. In the mid of 2006, this changes and DMA places increasing posterior mass on models that allow for time-variation in the parameters. In the period from the beginning of 2007 to the onset of the financial crisis, we see that the weight on TVP models somewhat decreases. During the financial crisis, we again experience a pronounced increase in posterior weight towards TVP regression. In that period, constant parameter models only play a limited role in forming inflation forecasts. With some few exceptions, the remainder of the hold-out period is characterized by evenly distributed posterior mass across constant and TVP regressions.
The middle panel of (ref) shows that DMA places increasing posterior mass on models with a large number of factors during the period prior to the global financial crisis. During the recession and the immediate period afterwards we observe that models with a different number of factors obtain substantial model weight. This suggests that if the number of factors is small, our dimension reduction techniques soak up information which might be useful during a recession (such as sharp changes in $\bm X$). If the number of factors becomes large this information is extracted as well but (potentially) reflected by more factors that display changes which are less pronounced. In the period after the global financial crisis we again find a large number of factors retrieving substantial posterior weight.
The bottom panel (panel (c)) of (ref) provides information on how much weight is allocated to models that exploit non-linear dimension reduction techniques. This figure corroborates our full sample findings: the Autoencoder performs extremely well and dominates our pool of models. Notice, however, that this statement is not true during the global financial crisis. During that period we observe that models based on PCA squared and PCA quadratic feature large weights. We also find that linear techniques (PCA linear) and other non-linear techniques (PCA with a Gaussian kernel, LLE, ISOMAP, diffusion maps) retreive almost no posterior weight over time.
Summing up this discussion we find that the single best performing model (the Autoencoder) is hard to beat when we dynamically combine models. However, this comparison is, to some extent, unfair since the researcher does not have this information at her disposal. Hence, combining models helps to integrate out this uncertainty by producing forecasts which are close to the single best performing model but, at the cost of higher computational costs, without the necessity of knowing the strongest single model specification.
In macroeconomics, the vast majority of researchers compress information using linear methods such as principal components to efficiently summarize information embodied in huge datasets in forecasting applications. Machine learning techniques describing large datasets with relatively few latent factors have gained relevance in the last years in various areas. In this paper, we have shown that using such approaches potentially improves real-time inflation forecasts for a wide range of competing model specifications. Our findings indicate that point forecasts of simpler models are hard to beat (especially at the one-month-ahead horizon). For density forecasts, however, we find that more sophisticated modeling techniques that rely on non-linear dimension reduction do particularly well. Among all the techniques considered, our results suggest that the Autoencoder, a particular form of a deep neural network, produces the most precise inflation forecasts (both in terms of point and density predictions). The large battery of competing models gives rise to substantial model uncertainty. We address this issue by using dynamic model averaging to dynamically weight different models, dimension reduction methods and priors. Doing so yields forecasts which are almost as accurate as the ones obtained from the single best performing models.
{\setstretch{0.85} \addcontentsline{toc}{section}{References} }
\setcounter{page}{1}