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.
85,984 characters · 24 sections · 59 citation commands
Kernel Three Pass Regression Filter
In recent years, the surge in high-dimensional datasets across fields like economics has ushered in new opportunities and challenges. A paramount issue is the `curse of dimensionality,' which undermines the effectiveness of traditional finite-dimensional estimation methods. Most modeling techniques applied to high-dimensional data assume the existence of a low-dimensional structure that effectively summarizes the data. One stylized feature of high-dimensional economic datasets is the presence of high and pervasive collinearity among variables, leading researchers to posit a data-generating process that assumes all variables are a function of a few latent factors. This formulation is commonly referred to as the factor model. A vast amount of literature focuses on using this latent factor structure for forecasting applications. A typical example is found in diffusion index models (stock2002diff), where latent factors are derived from a high-dimensional set of variables using Principal Components (hereafter, PC) method. These factors are subsequently utilized to forecast a target variable. A limitation of this PC-based factor estimation is its unsupervised nature, i.e., no information from the target variable is incorporated.
Since the primary goal is to forecast a target rather than estimate the underlying factor structure, introducing a degree of supervision can be beneficial. This can help filter out irrelevant information from the predictor set, thus enhancing the predictive accuracy. This can be done in different ways: using soft and hard thresholding methods to remove predictors with no predictive content, as in bai2008, or assigning varying weights to predictors based on their predictive capabilities for the target (see, for example, scaledpca), or estimate the subset of factors that exhibit predictive power for the target rather than the complete set of factors that drive the predictors, as in kelly.
The aforementioned models, though offering a supervised alternative to the unsupervised nature of Principal Components, nevertheless rest on the convenient assumption of linearity. However, as underscored in ml_jae, non-linearity often characterizes many predictive relationships in economics, particularly over extended time horizons and within data-rich environments.
Various approaches have been proposed to integrate non-linearity into factor models. For instance, squared principal components (PCs) or principal component squared ($PC^2$) as seen in bai2008, sufficient forecasting by fanSDR, the kernel trick to estimate factors (Varlam) among others. However, these approaches have limited supervision in the prediction process, if any. For example, fanSDR estimates factors through an unsupervised method (PC) and then derives sufficient indices using these PCs. Similarly, Varlam essentially applies kernel PCA (an unsupervised method) to estimate the set of factors driving a higher-dimensional space obtained by lifting the set of predictors through the kernel method. In bai2008, a very particular form of non-linearity (quadratic) is examined, which is somewhat ad hoc. Although they employ thresholding methods to reduce predictors to a smaller set, their screening method, however, may still encounter challenges in filtering relevant factors within this subset, leading to inefficient forecasts.
Our paper incorporates both non-linearity and supervision by introducing a novel kernel three-pass regression filter. The approach applies the three-pass regression filter of kelly to a non-linear transformation of the predictor set. While we adopt the lifting idea used in Varlam, we depart from unsupervised techniques such as kernel PCA and instead employ a supervised procedure to estimate factors obtained from non-linear transformations of the predictors that are relevant for the target variable. Our results align with the findings in ml_jae; over long horizons, where non-linear predictive relations become more pronounced, methods that capture non-linearity, such as ours, consistently outperform their linear counterparts. However, other non-linear approaches do not perform as well. Although they display the expected pattern of improving relative to linear models at longer horizons, the gains are substantially smaller. We attribute the superior performance of our approach to the combination of non-linearity and supervision. Importantly, the non-linear gains at longer horizons do not come at the cost of short-horizon underperformance; while it may not always be the best at very short horizons, it remains close to the best and in several cases even surpasses competing methods.
The paper proceeds as follows. Section 2 provides a brief introduction to Kernel methods. Section 3 introduces our estimator and discusses its similarity with the estimator of kelly. We also list a set of assumptions that ensure the theoretical properties of our estimators, which are given in the subsequent section 4. We present our empirical results in sections 5 and 6 and conclude in section 7. Mathematical proofs and implementation details are given in the appendix.
We use $\boldsymbol{y}$ to denote the $T \times 1$ vector of the target variable, i.e. \[ \boldsymbol{y} = (y_h, y_{h+1}, \dots, y_{T+h})'. \] We have $N$ predictors with $T$ observations for each predictor. The cross section of predictors at time $t$ is given by the $N \times 1$ vector $\boldsymbol{x}_t$, while the temporal observations of predictor $i$ are given by the vector $\mathbf{x}_i \in \mathbb{R}^T$. Stacking the predictors yields the $T \times N$ matrix \[ \boldsymbol{X} = (\boldsymbol{x}_1, \boldsymbol{x}_2, \dots, \boldsymbol{x}_T)' = (\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N), \] so that the $t$-th row of $\boldsymbol{X}$ is the cross section $\boldsymbol{x}_t' \in \mathbb{R}^{1 \times N}$, and the $i$-th column contains the temporal observations $\mathbf{x}_i \in \mathbb{R}^{T \times 1}$ of predictor $i$.
We also have $L < \infty$ proxies, collected in the $T \times L$ matrix \[ \boldsymbol{Z} = (\boldsymbol{z}_1, \boldsymbol{z}_2, \ldots, \boldsymbol{z}_T)', \] where each $\boldsymbol{z}_t \in \mathbb{R}^{L}$ is the vector of all $L$ proxies observed at time $t$.
We use $\varphi$ to denote a transformation $\varphi:\mathbb{R}^N\to\mathbb{R}^M$, which transforms each $\boldsymbol{x}_t$, at time $t$, into an $M\times 1$ vector $\varphi(\boldsymbol{x}_t)$. Applying $\varphi$ to the vector of predictors at each time point and stacking the resulting vectors yields the transformed predictor matrix, \[ \varphi(\boldsymbol{X}) = \big(\varphi(\boldsymbol{x}_1), \varphi(\boldsymbol{x}_2), \dots, \varphi(\boldsymbol{x}_T)\big)' = \big(\varphi_1(\mathbf{x}), \varphi_2(\mathbf{x}), \dots, \varphi_M(\mathbf{x})\big)\in\mathbb{R}^{T\times M}, \] where the $t$-th row of $\varphi(\boldsymbol{X})$ is the transformed cross section $\varphi(\boldsymbol{x}_t)'\in\mathbb{R}^{1\times M}$ at time $t$, and the $j$-th column $\varphi_j(\mathbf{x})\in\mathbb{R}^{T\times 1}$ represents the vector of temporal values of the $j$-th transformed feature across times $\{1, \dots, T\}$. $\boldsymbol{I}_L$ denotes an identity matrix of dimesnion $L$. We use $\boldsymbol{J}_T$ to denote the demeaning matrix, i.e., $\boldsymbol{J}_T \equiv \boldsymbol{I}_T - \frac{1}{T}\mathbf{\iota}_T\mathbf{\iota}_T'$, and for conformable matrices $\boldsymbol{U}$ and $\boldsymbol{V}$ we define $ \boldsymbol{S}_{UV} \equiv \boldsymbol{U}' \boldsymbol{J}_T \boldsymbol{V}. $. Stochastic orders are denoted by the usual ${O}_{p}$ and ${o}_{p}$. For a matrix, $\boldsymbol{O}_{p}$ and $\boldsymbol{o}_{p}$ denotes the element wise stochastic order, i.e., a matrix is said to be $\boldsymbol{O}_{p}(1)$ or $\boldsymbol{o}_{p}(1)$ if all it's elements are ${O}_{p}(1)$ or ${o}_{p}(1)$ respectively.
We consider a {transformation} \(\varphi: \mathbb{R}^N \to \mathbb{R}^M\) that maps the original set of predictors \(\boldsymbol{x}_t\) into a vector \(\varphi(\boldsymbol{x}_t)\) in a (typically higher- or potentially infinite-) dimensional space. These transformed vectors are often called features. Stacking the transformed predictors across all times yields the matrix \[ \varphi(\boldsymbol{X}) = \big(\varphi(\boldsymbol{x}_1), \varphi(\boldsymbol{x}_2), \dots, \varphi(\boldsymbol{x}_T)\big)' \in \mathbb{R}^{T \times M}. \]
Methods such as principal components and 3PRF, depend on the predictors only through the \(T \times T\) matrix of dot products \(\boldsymbol{X}\boldsymbol{X}'\). Extending these methods to the transformed set of predictors would require computing the Gram matrix \(\varphi(\boldsymbol{X}) \varphi(\boldsymbol{X})'\), which can be computationally expensive when \(M\) is large or infeasible if the transformed features lie in an infinite-dimensional space.
To address this, we introduce a kernel function $\mathcal{K}(\boldsymbol{x}_t, \boldsymbol{x}_s)$, which computes the inner product of the transformed predictors indirectly as
This approach, known as the `kernel trick', allows us to compute the inner products in the feature space without explicitly evaluating the transformation $\varphi$. Stacking these evaluations across all times yields the kernel (Gram) matrix \[ \mathcal{K}(\boldsymbol{X}, \boldsymbol{X}') = \varphi(\boldsymbol{X}) \varphi(\boldsymbol{X})' \in \mathbb{R}^{T \times T}, \] whose $(t,s)$-th entry corresponds to $\mathcal{K}(\boldsymbol{x}_t, \boldsymbol{x}_s)$ as defined in (ref). A function $\mathcal{K}$ is called a valid kernel if, for any collection of inputs $\{\boldsymbol{x}_1, \dots, \boldsymbol{x}_T\}$, the resulting matrix $\mathcal{K}(\boldsymbol{X}, \boldsymbol{X}')$ is symmetric and positive semi-definite. By Mercer’s theorem (see Appendix (ref)), every such valid kernel is associated with a Hilbert space $\mathcal{H}$ and a mapping $\varphi: \mathbb{R}^N \to \mathcal{H}$ such that the equality in (ref) holds. In other words, a valid kernel can always be interpreted as the inner product between transformed observations in some feature space.\footnote{In the Online Appendix-B.1, we illustrate how different kernel functions correspond to the inner products of transformed inputs.}
Working in a transformed space provides a substantial advantage: many nonlinear relationships among the original predictors can become linear in the appropriately transformed set of predictors. As a stylized example, consider the illustration in Figure (ref). We generate two variables $X$ and $Y$ from the uniform distribution $U[-2,2]$ and define a binary variable $z$ as \[ z =
\] In the original predictor space (Figure (ref), left), the two classes of $z$ cannot be separated by a linear boundary. However, after transforming the inputs to $\varphi_1(X) = \sin^2(X)$ and $\varphi_2(Y) = \cos^2(Y)$, the classes become linearly separable in the transformed feature space (Figure (ref), right), where blue and red points represent the two classes. Although this example involves only a two-dimensional input $\boldsymbol{W} = (X, Y)$ transformed into a two-dimensional feature space, in practice $\varphi(\boldsymbol{X})$ can be high- or even infinite-dimensional. Such transformations typically introduce a large set of nonlinear basis functions of the original predictors, making it more likely that an underlying nonlinear relationship in the data becomes linearly representable in the transformed space.
We delineate the three regression passes that we use to construct our forecast. The first two passes, as explained below, are not feasible in practice, whilst the eventual closed-form solutions are. Nonetheless, these steps offer valuable insights into the underlying process of our estimator and elucidate its similarity to the well-known linear three-pass filter proposed by kelly.
Below, we list the data generation process for the transformed predictor set ($\varphi(\boldsymbol{X})$), the target ($\boldsymbol{y}$), and the proxies employed for supervision ($\boldsymbol{Z}$). Given the data structure, it is easy to explain why this supervised methodology is effective in estimating the target relevant factors.
The mapping $\varphi$ transforms the original $N$-dimensional predictors $\boldsymbol{x}_t$ into an $M$-dimensional set of features, where $M$ depends on the input dimension $N$. In general, for any (fixed) valid kernel, $M$ increases with $N$, meaning that adding more predictors weakly increases the dimension of the (implied) transformed feature space for a given Kernel. For example, under a second-degree polynomial kernel, the transformed feature vector includes all original linear terms, all squared terms, and all pairwise products of the predictors, so that the total number of features is $M = N + N(N+1)/2$, which grows with $N$.\footnote{See the Online Appendix for more details.}
Assumption (ref) endows the transformed set of predictors with a factor structure. The coefficient vector $\boldsymbol{\beta}_f$ corresponds to the relevant factors $\boldsymbol{f}_t$ in the main predictive equation for $y_{t+h}$; being nonzero, it indicates that these factors directly influence the outcome, whereas the coefficients on the irrelevant factors $\boldsymbol{g}_t$ are zero. The dimensions $K_f$, $K_g$, and $L$ are all fixed and remain small relative to $N$.
The factor loading matrix for the transformed predictors $\boldsymbol{\Phi}$ has dimension $M \times K$, mapping the $K$ latent factors to the $M$ transformed features. Individual factor loadings $\boldsymbol{\phi}_k \in \mathbb{R}^M$ correspond to the weights of factor $k$ across the $M$ transformed features, with $\boldsymbol{\Phi} = [\boldsymbol{\phi}_1, \dots, \boldsymbol{\phi}_K]$. The submatrix $\boldsymbol{\Phi}_f$ is $M \times K_f$ and contains the loadings on the relevant factors, while $\boldsymbol{\Phi}_g$ is $M \times K_g$ and contains the loadings on the irrelevant factors. Similarly, the loading matrix for the proxies, $\boldsymbol{\Lambda}$, has dimension $L \times K$, mapping the $K$ latent factors to the $L$ proxy variables. Its submatrices $\boldsymbol{\Lambda}_f$ and $\boldsymbol{\Lambda}_g$ have dimensions $L \times K_f$ and $L \times K_g$, respectively, and are defined analogously, capturing the contribution of the relevant and irrelevant factors to each proxy.
The proxy variables $\boldsymbol{z}_t$ are generated by the underlying factors together with idiosyncratic proxy noise. Overall, the framework parallels kelly, with the key distinction that the factor structure is applied to transformed (and potentially nonlinear) features rather than the raw predictor space.
Our estimator is characterized by three infeasible three-passes, which are summarized in Table (ref).
These three passes rely on the fact that the correlation between the transformed $\varphi(\boldsymbol{X})$ and the proxies is only due to target relevant factors. Therefore, pass 1 of the regression asymptotically yields a rotation of the relevant-factor loadings of the $j^{th}$ predictor. Cross-sectional covariance between these loadings and the predictors, across $t$, is solely affected by the target relevant factor(s). Hence, pass 2 of this process traces the factor(s) out as a slope parameter. The last pass involves regressing the target variable on the estimated factor(s). Although these three passes offer valuable insights into the mechanics of our process, they are infeasible in practice due to the unavailability of the transformed inputs $\varphi(\boldsymbol{X})$. This is where the kernel trick proves to be useful. To see this, we note that factor(s), their predictive coefficients, and the forecast can be expressed in closed form as below,
The estimated factor(s) :
The estimated coefficient(s) of the factor(s) :
Finally, the estimated target :
These expressions are obtained by simply replacing $\boldsymbol{X}$ by $\varphi(\boldsymbol{X})$ in the three-pass regression filter of kelly.\footnote{The only difference is that the matrix $\boldsymbol{J}_N$, used in kelly to demean along the cross-sectional dimension, is absent here. By omitting the intercept in the regression of $\boldsymbol{X}$, we avoid the need for this demeaning step. This is justified because we work with the kernel function rather than explicit transformed features, so the transformation is implicit and passes 1 and 2 are infeasible in practice. Such an approach is standard in kernel-based methods like kernel PCA, where demeaning is performed along the temporal (sample) dimension but not along the transformed predictor dimension.} As evident from the expression of $\hat{\boldsymbol{F}}^{\prime}$, the filtration process applied on the transformed predictor space results in a favorable scenario where the eventual estimate of the factor(s) depends upon the transformed predictors only through their dot products in the transformed space. This holds true for $\hat{\boldsymbol{\beta}}$ and $\hat{\boldsymbol{y}}$ as well.
This inner product can be computed using a suitable kernel function. Alternatively, it can be inferred that employing a positive semidefinite (psd) kernel function to calculate dot products in these derived expressions is akin to executing the three-pass filter process on the transformed set of predictor(s), which, according to Mercer's theorem, are guaranteed to exist.
The key insight is that the variation in the predictors that is relevant for forecasting $y_{t+h}$ may lie on a low-dimensional, nonlinear manifold in the original predictor space. After applying a kernel-induced transformation $\varphi(\cdot)$, this nonlinear structure becomes much closer to linear in the high-dimensional feature space spanned by the transformed features.
By assuming that $\varphi(x_t)$ admits an approximate factor structure, we are effectively stating that the predictive manifold becomes well-approximated by a low-dimensional linear subspace in the transformed feature space. Therefore, while the conditional expectation $E[y_{t+h}\mid x_t]$ may be highly nonlinear in the original $N$-dimensional predictors, the same object becomes linear in a small number of relevant latent factors $f_t$ after the transformation $\varphi(\cdot)$.
This delivers the same idea as kelly: forecasting is driven by a few latent factors that summarize the relevant predictive variation. The only difference is that, rather than extracting these factors from the raw predictors, we extract them after a kernel-induced nonlinear transformation, which makes nonlinear predictive structure amenable to linear factor methods.
The Kernel three-pass regression, like the linear 3PRF, relies on the availability of suitable proxies. kelly show that such proxies can always be constructed using the target variable $\boldsymbol{y}$. That process is explained in Table-(ref).
Assumption (ref) lays out the factor structure of our model. Below, we delineate a set of additional assumptions under which our model delivers consistent forecasts.
Assumption (ref).(ref) imposes regularity on the latent factors by requiring their covariance matrix to converge to a well-behaved, non-stochastic limit, denoted $\boldsymbol{\Delta}_F$, a $K \times K$ matrix. Assumption (ref).(ref), adapted from kelly, imposes an analogous condition on the factor loadings, treating $\mathcal{P}$ as a non-stochastic $K \times K$ matrix representing their limiting cross-sectional covariance. Likewise, $\boldsymbol{\Delta}_\omega$ is a non-stochastic $L \times L$ matrix capturing the population covariance of the proxy noise. Since we assume a factor structure on the transformed space instead of the original predictor space, the normalization in various terms features $M$ and not $N$, where $M$ is the dimension of our transformed space. Assumptions (ref).(ref)-(ref).(ref), borrowed from kelly, impose regularity on various error processes.
Assumption (ref).(ref)-(ref).(ref) allow various forms of weak cross-sectional and temporal dependence between the idiosyncratic components of the transformed predictors. These assumptions characterize our `Approximate' factor model. The terminology of approximate, as opposed to a strict factor model, alludes to the allowance of these weak correlations, as outlined by chamberlain. These assumptions are standard in the literature; see bai2003. Assumption (ref).(ref)-(ref).(ref) are either borrowed from or are weaker versions of Assumptions in kelly. They are reasonable because each of them involves a product of orthogonal series.
Assumption (ref) characterizes the covariance matrices defined in Assumption (ref). This is a normalization assumption and is common in factor model literature, see kelly Assumption 5. It pertains to the non-identifiability of the true factor(s). It is well known that only the vector space spanned by the factor(s) can be consistently estimated but not the factor themselves. Imposing some normalization condition for the uniqueness of solution(s) is common in literature.
Assumption (ref) outlines the role and structure of the proxies in our framework. The proxy loading matrix is of the form $\boldsymbol{\Lambda} = [\boldsymbol{\Lambda}_f \;\; \mathbf{0}]$, so that the proxies load only on the relevant factors and have zero loading on the irrelevant factors ($\boldsymbol{\Lambda}_g = \mathbf{0}$). Non-singularity of $\boldsymbol{\Lambda}_f$ ensures that the proxies are linearly independent and none are redundant. Implicitly, this also requires that the number of proxies $L$ matches the number of relevant factors $K_f$, so that the relevant factor space is fully spanned by the proxies. In other words, the common component of the proxies captures all variation in the relevant factors, while none of the proxy variation is due to irrelevant factors.
We show that our estimated forecast converges to the infeasible best in probability. To show the same, we prove some intermediate results. All the proofs are in the appendix. \\[2mm] Define $\delta_{M T} \equiv min\{ \sqrt{M}, \sqrt{T}\}$. Define $\boldsymbol{H}_{f} \equiv \hat{\boldsymbol{F}}_{A} \hat{\boldsymbol{F}}_{B}^{-1} \mathbf{\Lambda}_f\boldsymbol{\Delta}_{f}$ where, $\hat{\boldsymbol{F}}_{A} = T^{-1} \boldsymbol{Z}^{\prime} \boldsymbol{J}_{T} \boldsymbol{Z}$ and $\hat{\boldsymbol{F}}_{B} =M^{-1} T^{-2} \boldsymbol{Z}^{\prime} \boldsymbol{J}_{T} \varphi(\boldsymbol{X}) \varphi(\boldsymbol{X}^{\prime}) \boldsymbol{J}_{T} \boldsymbol{Z}$, and $\boldsymbol{\Delta}_{f}$ denotes the submatrix of $\boldsymbol{\Delta}_F$ formed by its first $K_f$ rows and first $K_f$ columns.
This theorem establishes the estimated factor(s) convergence to the true relevant factors up to a rotation. It is well known in the literature on factor models\footnote{This feature of inherent unidentifiability has been emphasized in bai2003, kelly among other papers. The normalization imposed in assumption (ref) is done to handle this issue.}, that true underlying factor(s) are not identifiable; we instead estimate a rotated version of the true factors, which preserves their span.
Define $\boldsymbol{G}_{\beta} \equiv \hat{\boldsymbol{\beta}}_{1} ^{-1} {\hat{\boldsymbol{\beta}}_{2}} \left[\boldsymbol{\Lambda}_f \boldsymbol{\Delta}_{f}^3 \boldsymbol{\Lambda}_f^{\prime}\right]^{-1} \boldsymbol{\Lambda}_f \boldsymbol{\Delta}_{f}^2 $, where\\ $\hat{\boldsymbol{\beta}}_{1}=\hat{\boldsymbol{F}}_{A}$ and $\hat{\boldsymbol{\beta}}_{2}=\hat{\boldsymbol{F}}_{B}$, and $\boldsymbol{\Delta}_{f}$ is defined in the paragraph preceding Theorem (ref).
This theorem establishes the convergence of the predictive coefficients to a rotation of the true coefficients. Just like in the case of factor(s), true coefficients are not identifiable and we instead estimate their rotation. The rates established in Theorems (ref) and (ref) differ from the rates established in kelly and the reason is that our definition of rotation matrices $\boldsymbol{H}_f$ and $\boldsymbol{G}_{\beta}$ are different from kelly. (See Remark (ref)).
where $\mathbb{E}_{t} y_{t+h} = \beta_0+\boldsymbol{\beta}^{\prime} \boldsymbol{F}_t = \beta_0+\boldsymbol{\beta}_f^{\prime} \boldsymbol{f}_t $.
Combining Theorem (ref) and (ref), the convergence of $\hat{y}_{t+h}$ follows directly. Our proof, unlike kelly uses the convergence results for the estimated factor(s) and coefficients to obtain this result.
We apply our proposed method to real-world applications, focusing on forecasting time series variables across various economic domains such as national income, finance, labor, housing, prices, etc. To assess the performance of our approach, we conduct comparative analyses against competitive methods, employing the out-of-sample $R^2$ performance metric as a benchmark. The out-of-sample $R^2$ metric is computed as: \[ R^2 = 1 - \frac{\sum_{i \in \text{test-data}} (y_i - \hat{y}_i)^2}{\sum_{i \in \text{test-data}} (y_i - \bar{y}_{\text{train}})^2} \] It computes the out-of-sample proximity of our forecast $\hat{y}$ with the target ($y$) relative to a historical mean $(\bar{y})$; a positive value indicates that the forecast is better than the historical mean. Detailed explanations of performance metrics computation are provided in the Online Appendix-B.2.
This section provides an introduction to our competing methods and then a high-level summary of the empirical applications.
We compare our method against eight different forecasting approaches. The first is the principal components (PC) regression proposed by stock2002diff, denoted as PC in our performance tables. We also consider PC-Squared (PC-Sq) and Squared-PC (Sq-PC) from bai2008, kernel PCA (kPCA) [kpca], our linear counterpart \textit{3PRF} (kelly), and the autoregressive (AR) model with optimal lag $p$.\footnote{We evaluate several choices of $p$ and report the best-performing one. This represents the ideal performance for AR models in a setting where the true lag is known; in practice, $p$ cannot be chosen based on maximum out-of-sample $R^2$.}
The Diffusion Index (DI) serves as a benchmark forecasting method that combines AR and PCA. We select the number of principal components using the eigenvalue ratio test of ahn2013eigenvalue and run several DIs for different AR($p$), retaining the one with the $p$ that achieves the highest out-of-sample $R^2$. In this sense, we assess the upper bound of DI performance, since in practice the optimal $p$ cannot be chosen based on out-of-sample results.
Finally, we consider the Autoencoder (AE), a modern neural network method that learns a latent-space or manifold representation (AE_hinton). The crucial hyperparameter, the number of latent dimensions, is chosen via cross-validation, consistent with the hyperparameter tuning applied to the other methods. AE is conceptually the closest neural-network-based competitor to our approach, as it also learns a latent representation of the data. Some of these methods require hyperparameter tuning to achieve optimal results; we perform such tuning as described in subsection (ref).
As discussed in Section (ref), different kernel functions correspond to different $\varphi(\cdot)$. We use the Gaussian\footnote{The Gaussian kernel is mathematically equivalent to the commonly used RBF kernel in some papers.} kernel throughout our applications.
\[ K_\sigma(\mathbf{x}, \mathbf{x'}) = \exp \left( -\frac{\|\mathbf{x} - \mathbf{x'}\|^2}{2\sigma^2} \right) \]
where \(\sigma\) is a hyperparameter determined via cross-validation. The use of this kernel is justified by its strong performance in macroeconomic forecasting, as documented by rbf_ref1, k_ridge, and Varlam.
We begin by discussing the data sources, data transformations for stationary series, and the hyperparameter tuning procedure. The first empirical application illustrates the working of our method using theory-guided proxies. We then provide visual evidence showing how our method captures the nonlinearities present in the target variable compared with its nearest competitor. Next, we present the out-of-sample forecasting performance for several key variables of interest across a broad range of economic categories. Finally, we compare our method against others in a comprehensive forecasting exercise using all available series in the dataset as a target.
We utilize the quarterly macroeconomic dataset, FRED-QD. It spans the time period 1959-2023. This dataset encompasses a comprehensive array of more than 250 variables\footnote{There are 176 series for which we can find no missing values from 1959 to 2023, therefore our number of time series used are 176 in this analysis.}, including macroeconomic (such as GDP, Consumption, and Investment), financial, labor market, housing, and industrial and manufacturing variables. We present a tabulation comprising the mnemonic codes and details of the variables in the FRED-QD dataset alongside their counterparts in the Stock-Watson dataset in Online Appendix-B.4 for the series we forecast in this section.
Estimation exercises involving nonstationary time series pose significant challenges. Nonstationarity is ubiquitous in economic and financial data. Nonstationary variables lack a defined population mean, and the sample standard deviation tends to diverge as the number of observations increases, see Onatski_Wang for a more detailed discussion. Generally, researchers address this by manually examining each series to identify necessary transformations before computing principal components. hamilton_tranformation offers an improved method for transforming the predictors to achieve stationarity. We use their method to make our data stationarity.
Scholars in the literature often employ sample periods devoid of structural breaks. fan2023latent notes that “There exist significant structural breaks for many variables around the year of financial crisis in 2008 which makes our data non-stationary even after performing the suggested transformations”. Therefore, our study focuses on the stationary period spanning from 1965 to 2007\footnote{Another indirect advantage of the choice of this sample period is that it gives us the number of samples less than the number of predictors ($T<N$), hence a truly high-dimensional scenario to test our method in.}. For robustness checks, we examine three other sample periods. First, 1965–2019, which captures the Global Financial Crisis (GFC), allowing us to assess performance during high-volatility periods, as the GFC falls within the test sample under our 70% training split. Second, 1965–2023, the full sample\footnote{Although our data begins in 1959, we exclude early observations to accommodate AR(p) models with sufficiently large lags.}, encompassing both high- and low-volatility periods, including the GFC, COVID-19, and the Great Moderation. Third, 1984–2007, which corresponds to the Great Moderation and enables the evaluation of performance in a low-volatility environment. We found that the results are qualitatively similar across all sub-sample periods and convey the same story as in the main text (The results are tabulated in Online Appendix B.9). These findings indicate that our method is robust to both high-volatility (turbulent) and low-volatility (Great Moderation) periods.
In our main analysis, the sample period from 1964:Q4 to 2007:Q1 comprises $T=170$ observations (periods) and $N=176$ variables (predictors). While the data is initially available for around 250 series, those with missing values are removed; this leaves us with a total of 176 series. The model training and hyperparameter tuning are conducted within a rolling window framework, utilizing 70% of the total observations as the width of the rolling window. We observe qualitatively similar performance across varying window widths (50%, 60% of total data).
Our methodology incorporates the kernel as a fundamental element of the estimation process. The kernel function includes a hyperparameter that necessitates optimization. Concurrently, a similar hyperparameter requires tuning in the context of a competitor method, namely kernel PCA. Thus, we employ an identical tuning procedure for both methodologies. We use Gaussian kernel, which relies on a single hyperparameter, denoted as $\sigma$, for our specific applications. We partition the data into two folds and conduct cross-validation to determine the optimal tuning parameter. Further elaboration on the algorithm employed for this purpose is provided in the Online Appendix-B.5. We need to choose an optimal lag $p$ for the autoregressive (AR) and Diffusion Index (DI) methods, which we have already discussed. For the Autoencoder (AE) method, the key tuning parameters include the number of layers, layer size, and the dimension of the latent space, of which the latent space dimension is the most critical. We conduct experiments by varying the number of layers and layer sizes, and determine the optimal latent space dimension through cross-validation. Furthermore, a similar cross-validation procedure is applied to determine other hyperparameters in competing methods to ensure a fair comparison.\\[2mm] Among our competitive methodologies, where factors are computed as PCs, we are required to specify the number of factors. To address this, we employ the eigenvalue ratio test method proposed by ahn2013eigenvalue. This method computes the ratio of each eigenvalue to its predecessor and selects the number of principal components corresponding to the index where this ratio attains its maximum value. We employ a single factor throughout our analyses in both the 3PRF and kernel 3PRF models. This choice is often prudent within the 3PRF setting, as elucidated by kelly, who highlight instances where a single factor can effectively represent a multi-factor system. When factors exhibit the same variances\footnote{See appendix section A.7.2 in kelly}, a single proxy achieves optimal performance, and even when variances are not identical but closely aligned, one factor estimated through a single auto proxy typically explains a significant portion of the variation\footnote{See simulations in appendix A.7.3 in kelly}, rendering residual variation minimal. While we assessed the performance of our estimator with varying numbers of factors, we consistently observed that a single factor predominates; thus, we report results based on this configuration.
The primary objective of this subsection is to establish the viability of theory-guided proxies in forecasting using our method. We also compare with the linear benchmark, i.e., 3PRF of kelly. A more extensive performance evaluation will be presented in subsequent subsections, where we used the auto-proxy method discussed in table-(ref) to construct forecasts using Kernel 3PRF.
Consumption and investment are well-known drivers of GDP, and they likely capture the underlying latent factors that drive future economic activity. Using these variables as theory-guided proxies allows us to extract the relevant predictive components from the broader cross section of macroeconomic predictors. We construct GDP forecasts using these variables as proxies and report the results in Table-(ref).
This exercise proves the efficacy of theory-guided proxies. Our method outperforms the nearest competitive method kelly.
We also reproduce the theory-guided proxy example from kelly, focusing on one-period-ahead inflation forecasts. The target variable is the change in the price level ($\Delta$Price Level), and the dynamic quantity theory of money relates future inflation to current growth in real output and money supply:
\[ \frac{\Delta(\text{Money supply}) \times \Delta(\text{Velocity of money})}{\Delta (\text{Real Product})} = \Delta (\text{Price level}) \]
Table (ref) reports the results. Using both $\Delta$GDP and $\Delta$Money Supply as proxies captures the inflation-relevant variation more effectively than using either alone.
The results indicate that the theory-guided proxies successfully capture the key inflation dynamics, achieving performance comparable to the linear benchmark 3PRF. Importantly, this also demonstrates that in settings where the relationship is approximately linear, the kernel 3PRF can adapt and produce forecasts of similar quality, as appropriate tuning of the hyperparameters allows it to effectively capture linear effects. It is important to emphasize again that this analysis focuses on one-step-ahead forecasts, which are not the primary strength of our methodology. The purpose of presenting these results is solely to demonstrate the workings of the procedure through the theory-guided proxies.
\FloatBarrier
To visually demonstrate the enhanced performance of kernel 3PRF compared to its linear counterpart, we provide comparative performance plots across four distinct types of economic series spanning various domains: macroeconomic series (Exports), price series (GDP Deflator), manufacturing series (Industrial Production), and financial series (S&P 500 Index) in figure-(ref) and (ref). Plots of all other series on different forecast horizons are given in the Online Appendix-B.6.
\FloatBarrier
An astute economic decision, such as monetary policy formulation, hinges upon well-informed anticipations of future trends in macroeconomic and financial data. Consequently, forecasting macroeconomic variables emerges as a pivotal pursuit for economists. Quoting Federal Reserve of New York's website, kim_macro_forecasts notes, “In formulating the nation’s monetary policy, the Federal Reserve considers a number of factors, including the economic and financial indicators which follow, as well as the anecdotal reports compiled in the Beige Book. Real Gross Domestic Product (GDP); Consumer Price Index (CPI); Nonfarm Payroll Employment Housing Starts; Industrial Production/Capacity Utilization; Retail Sales; Business Sales and Inventories; Advance Durable Goods Shipments, New Orders and Unfilled Orders; Lightweight Vehicle Sales; Yield on 10-year Treasury Bond; S&P 500 Stock Index; M2”. We, therefore, aim to forecast some of these crucial indicators in this paper. We compare the performance of our model against the competitors. This section forecasts seven macro series: GDP, Consumption, Investment, Exports, Imports, Fixed Investment, and Industrial Production (Final).
To present the results in an organized manner, we create two tables. In Table-(ref), we display the forecasting performance for three series: GDP, Consumption, and Investment, which we informally refer to as `Group-I'. Table-(ref) presents a comparative analysis of forecasting performance for `Group-II' macro variables\footnote{ Variables' FRED-QD code and description can be found in Online Appendix-B.4}: Exports, Imports, Fixed Investments, and Industrial Production (Final Index). As defined earlier in the text, the reported numbers in the tables represent out-of-sample $R^2$ values across various forecast horizons ranging from one period ahead to twelve periods ahead.\\[2mm] Results highlight a secular observation that among various unsupervised forecasting methodologies—PC, Squared-PC, PC-Squared, and non-linear unsupervised approaches such as kernel PCA—none exhibit superior performance compared to our proposed method across any forecast horizon for the seven series under consideration. While the supervised linear forecasting model 3PRF demonstrates improved performance relative to the unsupervised techniques, it still falls short of our non-linear supervised approach. Notably, the autoregressive (AR) model and its variant, the Diffusion Index (DI), emerge as the sole contenders capable of surpassing our method at shorter horizons, albeit only marginally and for a few series. The neural-network-based competitor, the Autoencoder (AE), does not outperform our method at any forecast horizon—a result that is not surprising, as neural-network-based methods are known to require large sample sizes for effective training, which is not the case here.
In summary, our method significantly outperforms all competitors across longer horizons and competes well with AR(p) and DI in the short run. Therefore, it serves as a dependable and preferred forecasting framework across all forecast horizons in macroeconomic prediction tasks.
\FloatBarrier
This analysis aims to forecast key labor market and price variables. Within the labor market category, we focus on unemployment rates and total non-farm employment (Nonfarm Emp). We examine the GDP Deflator and the Consumer Price Index (CPI) for price variables. The GDP Deflator offers insights into overall inflation at the macroeconomic level, while the CPI captures inflation experienced by consumers at a more disaggregated level. The results of this analysis are summarized in Table (ref).
While the major story remains of the forecast performance results are qualitatively similar to those of aggregate macroeconomic series forecasting. We want to highlight the following: in the Unemployment Rate and GDP Deflator, our method beats AR(p) and DI method even in the short-run, barring $h=1$, and in the remaining two variables, our method is the best except initial two horizons.
\FloatBarrier
We evaluate the relative performance of our method across several key indicators: Privately Owned Housing Starts (HStart), Privately Owned Housing Starts in the Western Census region (HStart-W), GS-1 (Treasury Bills), GS-10 (Treasury Notes), and the S&P 500 Index. The first two indicators pertain to the housing market, while the latter three belong to the financial market. These financial variables are listed in ascending order of volatility.
We observe that the linear AR(p) and DI models perform exceptionally well for housing variables, with our method closely trailing them while outperforming all other competitors. In contrast, for financial variables, our method decisively surpasses all competing approaches as the forecast horizon extends to $h \geq 3$.
\FloatBarrier
We observe that the AR(p) model and its closely related method, DI, perform particularly well in the short run. It is well-established in the literature that AR(p) models are difficult to outperform over shorter forecast horizons. Possibly for this reason, some forecasting approaches evaluate performance after accounting for the variation captured by AR(p). For instance, our closest competitor kelly effectively forecasts the innovation of an AR(4) process for each variable, thereby “partialing out” the predictive information contained in the target variable’s own lags.
Following a similar approach, we also conduct an exercise in which we partial out AR(4) variation and subsequently compare the forecasting performance of all methods. We find that the qualitative patterns remain consistent, and in some cases are even more pronounced. The detailed results are presented in Online Appendix B.8.
To enhance the robustness of our empirical analysis, we conducted comparative assessments of our method against competing methods across all 176 series within our dataset. This entailed selecting each series as the target and repeating the comparative analysis for every series in our dataset.
Our investigation encompasses the comparative performance of models across a total of $176\times 8 =1408$ target-horizon combinations. The results of these comparisons indicate the percentage of instances where a particular method demonstrated superior performance to all other competitors. For example, if a method emerged as the best performer in 704 out of 1408 combinations, it would be represented by a value of 50 in the table. Essentially, we list the relative frequency of the occurrence of the best performance of a given method.
While the preceding frequency comparisons provide insight into the number of times each method proved superior to others, they do not measure the extent to which the best-performing method surpassed its nearest competitor. In other words, while method A may marginally outperform method B on one forecast horizon, method B might exhibit a considerable advantage over method A on another horizon. Then, the aforementioned frequency comparison may not depict the full picture. To account for this, we introduce a notion of `Tolerance' level. We call a method `best' under tolerance level $\epsilon$ if the out-of-sample $R^2$ of a method is within $\epsilon$ percentage lower than the best method's performance \footnote{For example, if the AR model is the best for a of series $y_\ell$ and horizon $h_0$ with a $R^2=0.60$. For tolerance=5, another method will also be considered `best' if its $R^2 \geq 0.60(1-5/100)=0.57$}. Therefore, for a non-zero tolerance, it is possible to have multiple `best' methods. \\[2mm] The “All Horizons” set of rows summarizes all 1408 comparisons, i.e., encompassing all horizons and all series. Recognizing that forecast objectives may vary in time horizon, we scrutinize comparative performances in short- and long-run contexts. The “Short-run" rows incorporate horizons $h={1,2,3,4}$, comprising 708 (calculated as $176 \times 4$) combinations, while the “Long-run" row includes horizons $h={6,8,10,12}$, similarly amounting to 708 combinations. Additionally, the portion labeled as “Excluding AR" excludes the auto-regressive AR(p) and closely related Diffusion Index (DI) methods and compares the remaining methods across all 1408 combinations. For more granular analysis, we report comparative performance numbers for each forecast horizon $h$. These numbers are reported in the Online Appendix-B.7.
It is important to note that multiple `best' methods may exist for a non-zero tolerance level, resulting in the sums of rows (in Table 5 in Online Appendix-B.7) exceeding 100 percent. However, for a tolerance level of zero, the rows sum to 100 percent.
\FloatBarrier
We present the results in table-(ref). The findings presented above yield several noteworthy observations. First, it is evident that unsupervised forecasting techniques, including PC, Squared-PC, PC-squared, Autoencoder, and kernel PCA, exhibit inferior performance across the majority of scenarios when compared to our method. Second, our method, kernel 3PRF, demonstrates unequivocal superiority in longer-horizon forecasting endeavors. Third, our method is unequivocally superior across all horizons when autoregressive (AR) and closely related DI methods are excluded. Fourth, even in the short run, our method is the single best-performing approach under strict (zero-tolerance) evaluations of absolute performance. Moreover, its performance remains comparable to the leading short-run benchmarks, AR(p) and DI, when we allow for some tolerance level.
Taken together, by demonstrating consistent performance across all series in our dataset, we conclude that our method can be regarded as a reliable forecasting tool across all forecast horizons.
\FloatBarrier
Building upon the three-pass regression filter by kelly, we introduce a new forecasting method, kernel three-pass regression filter. Through extensive empirical exercises, we show that this approach holds promise as a dependable forecasting tool. Improved performance can be attributed to two noteworthy features of our method. First, it integrates non-linear relationships by transforming input data into a higher-dimensional space, encapsulating its non-linear functions. Second, it operates as a supervised method, effectively filtering out and discarding irrelevant factors while predicting the target variable.