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.
63,072 characters · 11 sections · 77 citation commands
The Kernel Trick for Nonlinear Factor Modeling
\spacingset{1.45}
Over the past century, factor models have become an integral part of multivariate analysis and high-dimensional statistics, and have had a substantial effect on a number of different fields, including psychology (Thompson1938), biology (Hirzel2002) and economics (Chamberlain). In economics and finance, applications range from portfolio optimization (Fama1992) and covariance estimation (Fan2013) to forecasting (StockWatson2002a). The application that we consider is macroeconomic forecasting, where various forms of factor analysis have become state-of-the-art techniques for prediction.
The general idea behind factor analysis consists of determining a few latent variables, or factors, that drive the dependence of the entire outcomes. Factors are designed to capture the common dynamics in a large panel of data. This feature is crucial in the context of increasing availability of macroeconomic time series coupled with the inability of standard econometric methods to handle many variables. While classical econometric tools break down in such data-rich or “big data" environments, factor models help to compress a large amount of the available information into a few factors, turning the curse of dimensionality into a blessing.
Factor analysis possesses several attractive properties that justify a large amount of literature in its support. First, it effectively handles large dimensions thereby enhancing forecast accuracy in such regimes. This was demonstrated in StockWatson2002a and StockWatson2002b who used so-called diffusion indexes, or factors, in forecasting models when dealing with a large number of predictors. More recently, KimSwanson2018 find that factor augmented models nearly always outperform a wide range of big data and machine learning models in terms of predictive power. Second, due to its conceptual simplicity, this methodology found its use beyond academic research. For example, the Federal Reserve Bank of Chicago constructs the Chicago Fed National Activity Index (CFNAI) simply as the first principal component of a large number of time series. Third, factor analysis aligns naturally with the dynamic equilibrium theories as well as the stylized fact of Sargent77 of a small number of variables explaining most of the fluctuations in macroeconomic time series. And finally, factor estimates can also be used to provide efficient instruments for augmenting vector autoregressions (VARs) (Bernanke05) to assist in tracing structural shocks.
Factors, however, are not observable and need to be estimated. Two classical estimation strategies rely on either intertemporal or contemporaneous smoothing. The former casts the model into a state-space representation and estimates it by the maximum likelihood via the Kalman filter. The disadvantages of this approach are that it requires parametric assumptions and that it quickly becomes computationally infeasible as the number of predictor series grows \footnote{Interestingly, there has been some evidence to the contrary, see Doz2012}. Contemporaneous smoothing is a more predominant and computationally simpler way based on principal component analysis (PCA) (Pearson), nonparametric least-squares approach for estimating factors.
Forecasts are obtained via a two-step procedure. First, factors estimates are derived from the set of available time series by one of the two methods described above. Once the factors are estimated, run a linear autoregression of the variable of interest onto factor estimates and observed covariates (e.g. lagged values of the dependent variable).
We analyze a factor model that is high-dimensional, static and approximate. High-dimensional framework (BaiNg2002), as opposed to classical framework (Anderson1984), allows both time and cross-section dimensions to grow. Static models do not explicitly model time-dependence of factors contrary to more general dynamic counterparts (Forni2000a). Approximate factor structure (Chamberlain) is more flexible compared with a strict version (Ross1976) as it imposes milder assumptions on the idiosyncratic component.
Factor analysis is closely related to PCA, although the two are not the same (Jolliffe1986). It is, however, well documented that the two are asymptotically equivalent under suitable conditions (see the pervasiveness assumption in Fan2013 for a recent treatment). There are several results on consistency of PCA estimators of factors (Connor1986, StockWatson2002a, Bai2006 among others) for various forms of factor models. One of the most relevant of the results is established in BaiNg2002 who derive convergence rates of such estimators for an approximate static factor model of large dimensions.
Despite its widespread use, factor modeling, and diffusion index forecasting methodology in particular, is still fundamentally limited to linear framework. Over the past two decades, leading researchers noted multiple times (e.g. see StockWatson2002b, Bai2008, StockWatson2012, Hansen2015) that further forecast improvements “will need to come from models with nonlinearities and/or time variation" and that “nonlinear factor-augmented regression should be considered" for forecasting. There have been attempts to incorporate time dependence (see, for example, DelNegro2008, Mikkelsen2015 and Stevanovic2016 among others) as well as some work documenting the superiority of nonlinear models in time series context (Terasvirta1994, Giovannetti, KimSwanson2014). However, the literature on addressing nonlinearity within the factor modeling or diffusion index methodology framework, which arguably remains to be the state-of-the-art technique for macroeconomic prediction (Stevanovich2019), is scarce.
One of the first such attempts is Amemiya2001 who assume nonlinear errors-in-variables parametrization of a factor model, which is not designed for forecasting applications. The most prominent work with focus on a prediction exercise is Bai2008. They either augment the set of predictor time series with their squares and apply standard PC to the augmented set, or use squares of principal components obtained from the original (non-augmented) set. Another closely related work is Exterkate who substitute the linear second step with kernel ridge regression and discover that this leads to more accurate forecasts of the key economic indicators.
This study adds to the scarce literature on nonlinear factor models. Specifically, the factors are allowed to capture nontrivial functions of predictors. To circumvent the computational difficulties associated with such novelties, we use the kernel trick, or kernel method (Hofmann2008), which is discussed in the next section within the diffusion index methodology context.
The rest of the paper is organized as follows. Subsections 2.1 and 2.2 of Section 2 review the diffusion index methodology and the kernel trick; subsections 2.3 and 2.4 derive the estimators and provide the theoretical guarantees. Section 3 outlines the forecasting models, describes the data and provides the empirical results. Section 4 concludes and discusses possible extensions. All proofs are given in the Appendix.
Notation. For a vector $\mathrm{v}\in\mathbb{R}^{d}$, we write its $i$-th element as $v_i$. The corresponding $\ell_p$ norm is $\left\lVert\mathrm{v}\right\rVert_p = \big(\sum_{i=1}^d |v_i|^p\big)^{1/p}$, which is a norm for $1\le p \le \infty$. An inner product between two vectors of the same dimension is $\left\langle\mathrm{v_i},\mathrm{v_j}\right\rangle = \mathrm{v_i}'\mathrm{v_j}$. For a matrix $A\in\mathbb{R}^{m \times d}$, we write its $(i,j)$-th entry as $\{A\}_{ij} = a_{ij}$ and denote its $i$-th row (transposed) and $j$-th column as column vectors $A_{i\cdot}$ and $A_{\cdot j}$ respectively. Its singular values are $\sigma_1(A) \ge \sigma_2(A) \ge \ldots \ge \sigma_q(A)$, where $q=\min(m,d)$. The spectral norm is $\left\lVertA\right\rVert_2 = \underset{\mathrm{v}\ne0}{\max} \frac{\left\lVertA\mathrm{v}\right\rVert_2}{\left\lVert\mathrm{v}\right\rVert_2} = \sigma_1(A)$. The $\ell_1$ norm is $\left\lVertA\right\rVert_1 = \underset{1\le j\le d}{\max} \sum_{i=1}^m |u_{ij}|$ and $\ell_\infty$ norm is $\left\lVertA\right\rVert_\infty = \underset{1\le i\le m}{\max} \sum_{j=1}^d |u_{ij}|$. The Frobenius norm is $\left\lVertA\right\rVert_F = \sqrt{\langle A,A \rangle} = \sqrt{tr(A'A)} = \sqrt{\sum_{i=1}^q \sigma_{i}^2(A)}$. For a symmetric matrix $W\in\mathbb{R}^{d\times d}$ with eigenvalues $\lambda_1(W) \ge \lambda_2(W) \ge \ldots \ge \lambda_d(W)$, define $eig_r(W)\in\mathbb{R}^{d\times r}$ to be a matrix stacking $r\le d$ normalized eigenvectors in the order corresponding to $\lambda_1(W),\ldots, \lambda_r(W)$. Finally, for a sequence of random variables $\{X_n\}_{n=1}^\infty$ and a sequence of real nonnegative numbers $\{a_n\}_{n=1}^\infty$, denote $X_n = O_\mathbb{P}(a_n)$ if $\forall \epsilon >0, \exists M,N>0$ such that $\forall n>N, \; \mathbb{P}(\left\lvertX_n/a_n\right\rvert\ge M) < \epsilon$; and denote $X_n = o_\mathbb{P}(a_n)$ if $\forall \epsilon >0, \; \underset{n\to\infty}{\lim} \mathbb{P}(\left\lvertX_n/a_n\right\rvert\ge \epsilon) = 0$. Finally, let $\mathbf{1}_{1/T}$ be a $T\times T$ matrix of ones divided by $T$.
Our goal is to accurately forecast a scalar variable $Y_t$, given a $T\times N$ data matrix $X$ with $t$th row $X_t'$, or $X_{t\cdot}'$. Both the number of observations $T$ and the number of series $N$ are typically large.
Consider the following baseline model, known as a Diffusion Index (DI) model:
Equation (ref) is a linear forecasting model, where $Y_{t+h}$ is the value of the target variable $h$ periods in the future, $F_t$ is the vector of $r$ factors at time $t$, $W_t$ is a vector of $p$ observed covariates (e.g. an intercept and lags of $Y_{t+h}$), $\epsilon_{t+h}$ is a disturbance term. Equation (ref) specifies the factor model, where $X_t$ is vector of $N$ candidate predictor series, $\Lambda$ is a loading matrix for $r$ common driving forces in $F_t$, $e_t$ is an idiosyncratic disturbance; and $t=1,\ldots,T$. The latter equation can be rewritten in matrix form
where $X = [X_1,\ldots,X_T]'$ and $F = [F_1,\ldots,F_T]'$. Throughout the paper it is assumed that all series are weakly stationary, while variables in $X$ have been standardized, meaning that each variable is separately demeaned and set to have unit $\ell_2$-norm.
If the above set of equations is augmented with transition equations for $F_t$, we obtain a dynamic factor model which is estimated by the Kalman filter as discussed above. Let us instead focus on a nonparametric estimation approach as suggested in StockWatson2002a. The goal at first stage is to solve
where the restrictions are in place for identifying the unique solution (up to a column sign change). It is well known that the estimator of factor loadings $\widehat{\Lambda}$ is given by the $r$ eigenvectors associated with largest eigenvalues of $X'X$, while $\widehat{F} = X\widehat{\Lambda}$. This estimator $\widehat{F}$ is equivalent to principal component (PC) scores derived from the matrix $X$. Once we have an estimate of $F$, the second stage involves least squares estimation of equation (ref) with $F$ substituted with its estimate.
It is clear that the standard PC estimator reduces the dimensionality of $X$ linearly: $F_t$ represents the projection of $X_t$ onto $r$ eigenvector directions exhibiting the most variation. However, if there is a nonlinearity in $X$, that is if the true lower dimensional representation is a nonlinear submanifold in the original space, such linear projections will be inaccurate. There are several ways to take into account a possible nonlinearity. For example, Bai2008 propose a squared principal components (SPC) procedure, which applies the standard PCA algorithm to the matrix $X$ augmented by its square, that is $[X,X^2]$. Although this procedure supposedly leads to additional forecasting gains, it is limited by the second-order features of the data.
Other nonlinear dimension reduction techniques include Laplacian eigenmaps (Belkin2001), Locally-Linear Embeddings (Roweis2000), Isomaps (Tenenbaum2000) and a number of others (Hastie1989, Mark1991, Hinton2008). In this paper we use the approach that applies the kernel trick to the standard PCA, so\hypcalled kernel PCA (kPCA) (Scholkopf1999). This algorithm can be shown to contain a number of widely used dimensionality reduction methods, including the ones listed above (Hofmann2008). While it permits modeling a set of nonlinearities rich enough for successful applications in nontrivial pattern recognition tasks such as face recognition (Kim2002), the algorithm does not involve any iterative optimization.
In this subsection we review the kernel trick methodology and illustrate its usefulness. The kernel method implicitly maps the original data nonlinearly into a high\hypdimensional space, known as a feature space, $\varphi(\cdot): \mathcal{X} \rightarrow \mathcal{F}$. This space can in fact be infinite\hypdimensional which would seemingly prohibit any calculations. However, the trick is precisely in avoiding such calculations. The focus is instead on similarities between any two transformed data points $\varphi(\mathrm{x_i})$ and $\varphi(\mathrm{x_j})$ in the feature space\footnote{Formally, the feature space is thus a Hilbert space, that is a vector space with a dot product defined on it.} as measured by $\varphi(\mathrm{x_i})'\varphi(\mathrm{x_j})$, calculating which would at first sight require the knowledge of the functional form of $\varphi(\cdot)$. The solution is to use a kernel function $k(\cdot,\cdot): \mathcal{X} \times \mathcal{X} \rightarrow \mathbb{R}$, which would output the inner product in the feature space without ever requiring the explicit functional form of $\varphi(\cdot)$. Moreover, a valid kernel function guarantees the existence of a feature mapping $\varphi(\cdot)$ although its analytic form may be unknown. The only requirement for this is positive-definiteness of the kernel function (Mercer's condition, see (ref)), specifically, $\int \int f(\mathrm{x_i}) k(\mathrm{x_i}, \mathrm{x_j}) f(\mathrm{x_j}) d\mathrm{x_i} d\mathrm{x_j} \ge 0$, for any square-integrable function $f(\cdot)$.
This kernel function forms a Gram matrix $K$, known as a Kernel matrix, elements of which are inner products between transformed training examples, that is $\{K\}_{ij} = k(\mathrm{x_i},\mathrm{x_j}) = \varphi(\mathrm{x_i})'\varphi(\mathrm{x_j})$. What makes the kernel trick useful is the fact that many models can be written exclusively in terms of dot products between data points. For example, the ridge regression coefficient estimator can be formulated as $X'(XX' + \lambda I_T)^{-1}\mathrm{y}$\footnote{This can be derived by solving the dual of the ridge least squares optimization problem, however a simpler approach would be to apply the matrix identity $(B'C^{-1}B + A^{-1})^{-1}B'C^{-1} = AB'(BAB'+C)^{-1}$ to the usual ridge estimator $(X'X + \lambda I_N)^{-1}X'\mathrm{y}$.}, and hence the prediction for a test example $\mathrm{x_*}$ is $\hat{y} = \mathrm{x_*}' X'(XX' + \lambda I_T)^{-1}\mathrm{y}$, where the dependence is exclusively on inner products between the covariates. This property allows us to apply the kernel method by substituting dot products between original variables with their nonlinear kernel evaluations, that is dot products between transformed variables. Hence, an alternative form is $ k_*'(K + \lambda I_T)^{-1}\mathrm{y}$, where $k_*$ with $\{k_*\}_i = \{\varphi(\mathrm{x}_*)'\varphi(X)'\}_{i}$ is a vector of similarities between the test example and training examples in the feature space. In terms of the time complexity, the algorithm needs to invert a $T\times T$ matrix instead of inverting an $N\times N$ matrix.
The key advantage of the kernel method is that it effectively permits using a linear model in a high-dimensional nonlinear space, which amounts to applying a nonlinear technique in the original space. As a toy example, consider a classification problem shown in Figure (ref), where the true function separating the two classes is a circle of radius .5 around the origin. The left panel depicts observations from two classes which are not linearly separable in the original two-dimensional space. Applying a simple polynomial kernel of degree $2$, $k(\mathrm{x_i}, \mathrm{x_j)} = (\mathrm{x_i}'\mathrm{x_j})^2$, implicitly corresponds to working in the feature space depicted on the right panel, since for $\varphi(\mathrm{x_i}) = (x_{i1}^2, \; \sqrt{2}x_{i1}x_{i2}, \; x_{i2}^2)'$ we have $\varphi(\mathrm{x_i})'\varphi(\mathrm{x_j}) = (\mathrm{x_i}'\mathrm{x_j})^2$. In this toy example, a linear classifier could perfectly separate the observations in the right panel of Figure (ref).
While all valid kernel functions are guaranteed to have the corresponding feature space, in many cases it is implicit and infinite-dimensional, as, for instance, for the radial basis function (RBF) kernel $k(\mathrm{x_i}, \mathrm{x_j}) = e^{-\gamma\left\lVert\mathrm{x_i} - \mathrm{x_j}\right\rVert_2^2}$ (see (ref)). Again, luckily, the knowledge of the feature mapping is not required.
Suppose there is a nonlinear function $\varphi(\cdot): \mathbb{R}^N \rightarrow \mathbb{R}^M$, where $M \gg N$ is very large (often infinitely large), mapping each observation to a high-dimensional feature space, $X_t \rightarrow \varphi(X_t)$. For now we consider $M$ to be finite for simplicity of exposition, so the original $T\times N$ data matrix $X$ can be represented as a $T\times M$ matrix $\Phi = [\varphi(X_1),\ldots,\varphi(X_T)]'$ in the transformed space, which may not be observable. Infinite-dimensional case induces several complications and is considered later.
Both the original $X$ and its transformation $\Phi$ are assumed to be demeaned. The latter requirement is simple to incorporate in the kernel matrix despite the mapping being unobserved. Specifically, supposing the original (non-demeaned) transformation is $\widetilde{\Phi}$, the kernel associated with demeaned features is
where $\widetilde{K} = \widetilde{\Phi} \widetilde{\Phi}'$ is based on the original $\widetilde{\Phi}$.
Our modeling of nonlinearity is through the feature mapping $\varphi(\cdot)$. This function replaces the original variables of interest in equation (ref) with their transformations,
where the subscript $\varphi$ indicates the association with the transformation. By stacking these into a $T\times M$ matrix $\Phi$ we can rewrite the minimization problem (ref) as
Note that solving this directly through the eigendecomposition of $\Phi'\Phi$ is generally infeasible, since $\Phi'\Phi$ is $M\times M$ dimensional. Even if the dimension $M$ was not prohibitive, the map $\varphi(\cdot)$ is unknown for interesting problems rendering any computation dependent on $\Phi$ or $\Phi'\Phi$ alone impossible. Fortunately, it is possible to reformulate this problem in terms of the $T\times T$ Gram matrix $K = \Phi\Phi'$.
While we are assuming $M$ to be prohibitively large but finite, the following decomposition generalizes to infinite dimensions. Starting from the “infeasible" eigendecomposition of the unknown covariance matrix of $\Phi$
where the eigenvalues $\lambda_i^c = \lambda_i(\frac{\Phi'\Phi}{T})$ satisfy $\lambda_1^c \ge \lambda_2^c \ge \ldots \ge \lambda_T^c$ and $\lambda_j^c=0$ for $j > T$ (assuming $M\ge T$) and $V^{[i]}$ is an $M$-dimensional eigenvector associated with the $i$th eigenvalue $\lambda_i^c$.
The key is to observe that each $V^{[i]}_\varphi$ can be expressed as a linear combination of features
where $A^{[i]} = \frac{\Phi V^{[i]}_\varphi}{\lambda_i^c T} = \left[ \alpha_1^{[i]}, \; \ldots \;, \alpha_T^{[i]} \right]'$ is a vector of weights which is determined next. Plugging this back into (ref) yields
Finally, premultiplying equation (ref) to the left by $\Phi$ and removing $K=\Phi\Phi'$ from both sides we obtain
hence the $i$th vector of weights $A^{[i]}$ corresponds to an eigenvector of a finite-dimensional Gram matrix $K$ associated with the $i$th largest eigenvalue $\lambda_i(\frac{K}{T}) = \lambda_i^c$ with $\lambda_j(\frac{K}{T})=0$ for $j > T$.
Notice that while solving the eigenvalue problem of $\frac{K}{T}$ allows to compute $A^{[i]}$, we are still unable to obtain the vector $V^{[i]} = \Phi'A^{[i]}$ since $\Phi$ may not be known. However, the main object of interest is recoverable: to calculate principal component projections, we project the data onto (unknown) eigenspace,
Stacking estimated factors corresponding to the first $r$ eigenvalues, define a $T\times r$ matrix $\widehat{F}_\varphi = \left[\widehat{F}^{[1]}_\varphi, \; \ldots\;, \widehat{F}^{[r]}_\varphi \right]$, where the subindex $r$ is dropped to simplify the notation. We refer to the factors constructed this way as kernel factors.
A similar alternative solution that only involves the Gram matrix has been known in econometrics since at least Connor1993. In particular, for a given $r$ the optimization problem in (ref) with the identification constraints $T^{-1}F_\varphi'F_\varphi = I_r$ and diagonal $\Lambda_\varphi'\Lambda_\varphi$, has the solution $\widetilde{F}_\varphi = \sqrt{T} eig_r(\Phi \Phi') = \sqrt{T} A_r$, where $A_r = \left[\widehat{A}^{[1]}, \; \ldots\;, \widehat{A}^{[r]} \right]$. Hence, the kPCA estimator is equivalent to the latter premultiplied by $\frac{\Phi \Phi'}{\sqrt{T}} = \frac{K}{\sqrt{T}}$. Note, however, that both estimators yield the same predictions when passed to the main forecasting equation (ref) as they have identical column spaces. This idea is summarized in the following proposition.
Importantly, we also establish that certain commonly used kernels allow the kernel factor estimator to incorporate the usual PC estimator. The following proposition demonstrates that RBF and sigmoid kernels allow to nest (a constant multiple of) the PC estimator for limiting values of the hyperparameter.
Adjustment by $\gamma^{-1}$ is necessary as the entries of $\widehat{F}_\varphi = K eig_r(K)$ approach $0$ as $\gamma \to 0$. This is due to the nature of $K = \widetilde{K} - \mathbf{1}_{1/T}\widetilde{K} - \widetilde{K}\mathbf{1}_{1/T} + \mathbf{1}_{1/T}\widetilde{K}\mathbf{1}_{1/T}$, where $k_{ij}$ converge to $0$ for both kernels for all $i,j=1,\ldots,T$. Luckily $\gamma^{-1}K$ has a nice non-degenerate limiting value (shown in the proof) that ensures the limit in proposition (ref) holds.
Proposition (ref) states that the (properly scaled) kernel factor matrix converges pointwise to its PCA analog (up to a sign flip) as the value of the hyperparameter $\gamma$ nears zero. That is, in the limit the two factor estimators are constant multiples of one another and hence produce the same forecasts. This ability to mimic the linear estimator is important for certain applications. For example, in macroeconomic forecasting it is notoriously difficult to beat linear models in a short horizon prediction exercise.
Figure (ref) illustrates the implicit procedure for obtaining $r$ kernel factors. The selected kernel function induces nonlinearity $\varphi(\cdot)$ on each element of the input layer, $N$-dimensional observations $X_{1\cdot}, \ldots, X_{T\cdot}$. Next, pairwise similarities between high\hypdimensional vectors $\varphi(X_{1\cdot}), \ldots, \varphi(X_{T\cdot})$ are computed, with $k(X_{i\cdot},\cdot) = \Big[ \varphi(X_{i\cdot})'\varphi(X_{1\cdot}),\allowbreak \; \ldots \;, \allowbreak \varphi(X_{i\cdot})'\varphi(X_{T\cdot}) \Big]'$. Finally, each factor is obtained as a linear combination of these inner products with the weights given as a solution to the eigenvalue problem discussed above. Of course, the kernel PCA algorithm does not explicitly nonlinearize the inputs as the kernel trick allows us to immediately calculate similarities and avoid the expensive high-dimensional computation. However, as mentioned earlier, the existence of such implicit nonlinearities is guaranteed by Mercer's theorem ((ref)). As opposed to standard feedforward neural networks, there is no iterative training involved and no peril of being trapped in local optima. On the other hand, it has been noted that certain kernels, including RBF and sigmoid, allow extracting features of the same type as the ones extracted by neural networks (Scholkopf1999). The two necessary steps involve evaluation of similarities in the kernel matrix and solving its eigenvalue problem; the complexity is thus dependent only on the sample size.
Depending on the choice of the kernel, the induced feature space could be either finite or infinite. A polynomial kernel considered earlier generates a finite-dimensional feature space, consisting of a set of polynomial functions over the inputs. In this simple case, the eigenspace associated with the kernel factor estimator in equation (ref) can generally be consistently estimated within the framework of Bai2003. Particularly, proposition (ref) and the following theorem in Bai2003 immediately imply $\sqrt{M}$-consistency of $\widetilde{F}_\varphi$.
For the model associated with equation (ref):
Assumption A: There exists a constant $c_1 < \infty$ independent of $M$ and $T$, such that
Part (a) is standard in factor model literature. Part (b) ensures pervasiveness of factors in the sense that each factor has non-negligible effect on the variability of covariates. This is crucial for asymptotic identification of the common and idiosyncratic components. Parts (c$\cdot$) partially permit time-series and cross-section dependence in the idiosyncratic component, as well as heteroskedasticity in both dimensions. Possible correlation of $e_{\varphi,it}$ across $i$ sets up the model to have an approximate factor structure. Part (d) allows weak dependence between factors and idiosyncratic errors and (e) is for identification.
The following theorem establishes consistency of kernel factors, for kernels inducing finite nonlinearities, in a large-dimensional framework.
Some comments are in order. First, the theorem states that the squared differences between the proposed factor estimator and (a rotation of) the true factor vanish as $M,T \to \infty$. While true factors themselves are not identifiable unless additional assumptions are imposed (see BaiNg13), identification of the latent space spanned by factors is just as good as exact identification for forecasting purposes. Second, since for any given number of original variables $N$ the dimension of the transformed space $M$ is fixed and finite, the growth in $M$ is only possible through $N$ and thus the limit on $M$ implies one on $N$. Third, this result does not imply uniform convergence in $t$. Lastly, the results suggest the possibility of $\sqrt{T}$-consistent estimation of the forecasting equation with respect to its conditional mean.
Unfortunately, this result does not generalize to the most interesting kernels (e.g. RBF) inducing infinite-dimensional Hilbert spaces, rendering the traditional approach unsuitable for establishing theoretical properties. Hence we turn to a functional analytic framework which allows rigorous treatment of infinite-dimensional spaces and which has been the classical framework for analyzing statistical properties of functional PCA, and kernel PCA in particular, in machine learning literature.
As will be shown later, it turns out that we can still show that the estimator concentrates around its population counterpart. We briefly present the necessary terms for understanding this result, without aiming to be exhaustive. A sufficiently detailed introduction to the analysis in Hilbert spaces can be found in Blanchard06, while a classical reference for operator perturbation theory is Kato52.
One of the first major investigations of the statistical properties of kernel PCA can be found in Shawe-Taylor02. The study provides concentration bounds on the sum of eigenvalues of the kernel matrix towards that of corresponding (infinite-dimensional) kernel operators. This permits to characterize the accuracy of kPCA in terms of the reconstruction error, that is the ability to preserve the information about a high-dimensional input in low dimensions. This is of significant interest for certain types of applications, such as pattern recognition. Blanchard06 further extend the results of the aforementioned study and improve the bounds on eigenvalues using tools of perturbation theory.
Although the theoretical discussion and the mathematical approaches developed in this literature are extremely valuable, our interest is not in kPCA's ability to reconstruct a given observation. The kernel factor estimator only reduces the dimensionality and passes it to the next stage, without ever going through the reconstruction phase. As the dimensionality is reduced by projecting observations onto the eigenspace -- the space spanned by eigenfunctions of the true covariance operator with largest eigenvalues -- the interest is in convergence of empirical eigenfunctions towards the true counterparts. Importantly, the proximity of eigenvalues does not guarantee that underlying eigenspaces will also be close.
We now briefly introduce the technical background for understanding our result. Let $\mathcal{H}$ be an inner product space, that is a linear vector space endowed with an inner product $\left\langle\cdot,\cdot\right\rangle_\mathcal{H}$, commonly denoted together as $(\mathcal{H}, \left\langle\cdot,\cdot\right\rangle_\mathcal{H})$. An inner product space $(\mathcal{H}, \left\langle\cdot,\cdot\right\rangle_\mathcal{H})$ is a Hilbert space if and only if the norm induced by the inner product $\left\lVert\cdot\right\rVert_\mathcal{H} = \left\langle\cdot. \cdot\right\rangle_\mathcal{H}^{1/2}$ is complete\footnote{Specifically, the limits of all Cauchy sequences of functions must be in the Hilbert space.}. The results below rely on this norm, so we drop the subscript $\mathcal{H}$ to simplify the notation.
Let $X$ be a random variable taking values on a general probability space $\mathcal{X}$. A function $k: \mathcal{X} \times \mathcal{X} \to \mathbb{R}$ is a kernel if there exists a real-valued Hilbert space and a measurable feature mapping $\phi: \mathcal{X} \to \mathcal{H}$ such that $\forall x,x' \in \mathcal{X}, \; k(x,x') = \left\langle\phi(x), \phi(x')\right\rangle.$ A function $k: \mathcal{X} \times \mathcal{X} \to \mathbb{R}$ is a reproducing kernel of $\mathcal{H}$ and $\mathcal{H}$ is a reproducing kernel Hilbert space (RKHS), if $k$ satisfies (i) $\forall x \in \mathcal{X}, k(\cdot,x) \in \mathcal{X}$, and the reproducing property (ii) $\forall x \in \mathcal{X}, \forall f \in \mathcal{H}, \left\langlef,k(\cdot,x)\right\rangle = f(x)$. An important result from Aronszajn50 guarantees the existence of a unique RKHS for every positive definite $k$.
Assume that $\operatorname*{\mathbb{E}} \phi(X) = 0$ and $\operatorname*{\mathbb{E}}\left\lVert\phi(X)\right\rVert^2<\infty$. A unique covariance operator on $\phi(X)$, $\Sigma = \operatorname*{\mathbb{E}} \phi(X) \otimes \phi(X)$, satisfying $\left\langleg,\Sigma h\right\rangle_\mathcal{H} = \operatorname*{\mathbb{E}} \left\langleh,\phi(X)\right\rangle_\mathcal{H} \left\langleg,\phi(X)\right\rangle_\mathcal{H},\; \forall g,h \in \mathcal{H}$, always exists (Theorem 2.1 in Blanchard06) and is a positive, self-adjoint trace-class operator. Denote $\widehat{\Sigma} = \frac{1}{T} \sum_{i=1}^{T}\phi(X_i)\otimes \phi(X_i)$ to be its empirical counterpart. Finally, an orthogonal projector in $\mathcal{H}$ onto a closed subspace $V$ is an operator $\Pi_V$ such that $\Pi_V^2 = \Pi_V$ and $\Pi_V = \Pi_V^*$.
We now lay out two lemmas which lead to the result. Lemma (ref) bounds the difference between the true and empirical covariance operators.
Note that both sigmoid and RBF kernels are bounded and hence satisfy the requirement of the above lemma. Lemma (ref) is an operator perturbation theory result and is adapted from Koltchinskii2000 and Zwald:
Finally, the following theorem bounds the difference between empirical and true eigenvectors.
Some comments are in order. First, note that we require the eigenvalues to be distinct, a well-known restriction (similar to $sin(\theta)$ theorem of Davis70), since it is impossible to identify eigenspaces with the same eigenvalues. Second, Theorem (ref) suggests that the eigenspace estimated by kernel PCA will concentrate close to the true eigenspace. Our kernel factor estimator simply projects onto that eigenspace and hence the precision is expected to increase as $T\to\infty$. Third, this rate does not address the case when variables exhibit dependence, although the exercise in the next section is indicative of some form of concentration, which suggests that the theoretical assumptions might be too conservative. Lastly, it may be possible to obtain a sharper bound since the proof relies on crude inequalities (e.g. triangle inequality).
This subsection discusses specific forms of equations (ref) and (ref) that are used for forecasting. Autoregressive Diffusion Index (ARDI) model is specified as
where superscript $h$ indicates dependence on the time horizon. Note, $P_t^h, M_t^h, K_t^h$ are the number of lags of the target variable, the number of lags of factors, the number of factors respectively. These three parameters are estimated simultaneously for each time period and time horizon using BIC. Since the true factors are unknown, we instead plug in the estimates from the factor equation, which is discussed next.
Three factor equation specifications are considered,
Factors in equation (ref) are estimated by PCA. Equation (ref) is similar, replacing the left-hand side with an augmented vector $X_{*,t}=[X_t,X_t^2]$. This procedure was dubbed as squared principal components (SPC) in Bai2008. Finally, the last equation applies nonlinearity induced by the selected kernel and is estimated by kPCA.
Forecasts using PCA and SPC are produced in three steps. First, we extract three factors from the set transformed and standardized predictors using one of the two methods. Second, three parameters are determined according to BIC for each out-of-sample forecasting period and each prediction horizon: the number of lags of the target variable $P^h_t$, the number of lags of factors $M^h_t$, the number of factors $K^h_t$. Third, the forecasting equation is estimated by least squares and forecasts are produced. The procedure for predicting with kPCA is similar, except there is an additional step where the value of the hyperparameter is specified, and the estimation is instead made in accordance with Algorithm (ref).
The algorithm starts by computing similarities via a given kernel. This induces transformed observation that need to be demeaned according to equation (ref). As in ordinary PCA, one then needs to compute the eigenvectors. The difference here is that our object is a Gram matrix, not a covariance matrix. The kernel factors are then simply derived as projections of the kernel matrix onto $r$ eigenvectors associated with the largest eigenvalues.
As an empirical investigation, we examine whether using kernel factors leads to improved performance in forecasting several key macroeconomic indicators. We use a large dataset from FRED-MD (McCracken), which has become one of the classical datasets for empirical analysis of big data. Its latest release consists of 128 monthly US variables running from $1959:01$ through $2020:04$, 736 observations in total. Following previous studies, we set $1960:01$ as the first sample, leaving 724 observations. Since the models presented in this study require stationary series, each of the variables undergo a transformation to achieve stationarity. The decision on a particular form of transformation is generally dependent on the outcome of a unit root test, which is known to lack power in finite samples. So instead, following McCracken, all interest and unemployment are assumed to be $I(1)$, while price indexes are assumed to be $I(2)$. The transformations applied to each series are described in supplemental materials.
We aim to predict a single time series from this dataset by utilizing the remaining variables. The series to be predicted include 8 variables characterizing different aspects of the economy. Specifically, we take one series from each of the eight variable “groups" in the dataset. The summary is provided in Table \hyperref[app:series]{A1} in Appendix.
Forecasts are constructed for $h=1,3,6,9,12,18,24$ months ahead with a rolling time window, the size of which is taken to be $120-h$. Thus, the pseudo-out-of-sample forecast evaluation period is $1970:01$ to $2020:04$, which is 604 months. We estimate 6 variants of autoregressive diffusion index models. The first model, taken as a benchmark, is a classical ARDI with PC estimates. Several studies have documented a strong performance of this model (see for example, Stevanovich2019). The second and third take SPC and so-called PC-squared ($PC^2$) estimates (Bai2008) respectively, where the latter is identical to the first model with squares of factor estimates added in the forecasting equation. The remaining models are based on kPCA estimates with three different kernels: a sigmoid $k(\mathrm{x_i}, \mathrm{x_j}) = \tanh(\gamma(\mathrm{x_i}'\mathrm{x_j})+1)$, a radial basis function (RBF) $k(\mathrm{x_i}, \mathrm{x_j}) = e^{-\gamma\left\lVert\mathrm{x_i} - \mathrm{x_j}\right\rVert^2}$ and a quadratic\footnote{Polynomial kernels of lower and higher order demonstrated poor forecasting ability and are not included.} polynomial (poly(2)) kernel $k(\mathrm{x_i}, \mathrm{x_j}) = (\mathrm{x_i}'\mathrm{x_j}+1)^d, \; d=2$.
Optimal parameters for each model at each step, $P^h_t,M^h_t,K^h_t$ and the kernel hyperparameter, are determined within the rolling window period, that is our setup only permits the information set that would be available at the moment of making a prediction. Specifically, $P^h_t,M^h_t,K^h_t$ are selected by BIC (maximum value allowed for each is set equal to $6$) for each out-of-sample period, while $\gamma$ is determined over a grid of values by so-called time series cross-validation. Specifically, we consecutively predict the latest 5 available observations and select the hyperparameter that minimizes the average error. The standard cross-validation may not theoretically be fully adequate due to the presence of serial correlation in the data and several approaches were suggested to correct it (Racine00).
The main empirical findings are presented in Table (ref). The subtitles of series indicate the names of dependent variables, while each value in the table represents the ratio of out-of-sample MSPE of a given estimation method to out-of-sample MSPE of the autoregression augmented by diffusion indexes estimated by PCA. The results range for $8$ variables across $7$ different forecast horizons. Our results are reproducible: in supplemental materials we provide the script written in Python $3.6$ that generates all results within a few hours.
Some comments are in order. First, sigmoid and RBF kernel approaches do lead to improved forecasting accuracy, especially at medium- and long-term time horizons. One of these two approaches dominates others in nearly 95% of the cases considered. Additionally, Diebold-Mariano test (Diebold2002) of equal predictive ability suggests that only these two methods are capable of often significantly outperforming the baseline PCA autoregression across a range of variables and horizons, while never being dominated. The kernel method is least advantageous for one-step-ahead forecasting. The phenomenon that linearity is hard to beat in a very short horizon is rather well known in the literature. Luckily, as was shown in Proposition (ref), kPCA is capable of mimicking a linear PCA by adjusting the kernel hyperparameter closer to 0, which often leads to the parity of the two methods in near-term forecasting. While the gains are not pronounced at $h=1$, they become apparent at longer horizons. Results vary across variables, but the improvement is prevailing at medium-term horizons and is uniform in one-year and longer predictions. The superiority exhibited by kPCA in many cases is remarkable for macroeconomic forecasting literature.
Second, both SPC and PC$^2$ perform substantially worse than a simple PCA. This result contradicts to Bai2008, but is consistent with a recent empirical comparison of Exterkate. Besides, PC$^2$ is extremely unreliable for long-term predictions. Similar to SPC, poly(2) kernel seeks to model the second-order features of the data and, as a result, often performs on par with SPC. As pointed out by the referee, the additional squared terms might often be leading to inefficiency.
Ultimately, note that kPCA's computational complexity is dependent on the number of time periods for estimation $120-h$, making kPCA slightly faster in this particular exercise. Most importantly, kPCA's advantage would grow in a macroeconomic setting, where “bigger" data (i.e. larger $N$) is becoming the norm.
In this study we have introduced a nonlinear extension of factor modeling based on the kernel method. Although our exposition mainly focused on a feature mapping $\varphi(\cdot)$ enforcing nonlinearity, it is also convenient to think of this approach as kernel smoothing in an inner product space. That is, the kernel factor estimator implicitly relies on the weighted distances between original observations. This alternative viewpoint presumes that analyzing the variation in the inner product space, rather than the original space, may be more beneficial. This idea had a profound impact on machine learning and pattern recognition fields, especially as regards to support vector machines (SVMs). By using a positive definite kernel, one can be very flexible in the original space while effectively retaining the simplicity of the linear case in the high-dimensional feature space.
We have demonstrated that constructing factor estimates nonlinearly can be beneficial for macroeconomic forecasting. Specifically, the nonlinearity induced by the sigmoid and RBF kernels leads to considerable gains at medium- and long-term time horizons. This gain in performance comes at no substantial sacrifice, the algorithm remains scalable and computationally fast.
There are several possible extensions. First, it is interesting to see how the performance would change have we pre-selected the variables (targeting) before reducing the dimensionality. As shown in Bai2008 and Bulligan2015 this generally leads to better precision. Second, the forecasting accuracy can be compared with other nonlinear dimension reduction techniques mentioned earlier, such as autoencoders. For the latter, however, one must be aware of the possibility of implicit overfitting by tuning the network architecture. This is not an issue in the current framework as there are a lot fewer parameters to specify. Third, the static factors considered here could possibly be extended to dynamic factors (Forni2000a), by explicitly incorporating the time domain, or “efficient" factors, by weighing observations by the inverse of the estimated variance.