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.
75,182 characters · 20 sections · 53 citation commands
Variable Selection in Macroeconomic Forecasting with Many Predictors
Recent development in information technology makes it possible to collect hundreds of economic variables in real time, with a reasonable cost. In such data-rich environment, using many economic predictors to forecast a few target variables has become a new trend in econometric research. In the last two decades, both theoretical and empirical works have been substantially built up on this direction, especially in the fields of macroeconomic forecasting SW2002a, BaiNg2008 and real-time now-casting GIANNONE2008,BANBURA2013. This new trend has also made practical impact – economic forecasts using many predictors are currently being produced by fiscal and monetary authorities in both the U.S. FREDMD,FREDQD and Europe ECB,REPEC2013. A key aspect of many-predictor forecasts is to impose suitable parsimonious structure on data so that the curse of dimensionality is circumvented and useful information can be extracted. There are two directions to accomplish this, which are based on two different assumptions about the economic data structure.
The first direction is factor augment (FA) approach. It has been found to produce superior forecasts over traditional methods such as AR and VAR, thus attained favor from both econometricians and practitioners. This approach assumes that many predictors are relevant to the target variable and they have a factor structure. Dynamic factor model is applied first to compress the information of predictors into a handful of estimated factors. Then, the factors are augmented to a linear forecasting equation for the target variable. The rationale behind this approach is that the common variation among many observed economic variables can be represented by a handful of unobserved factors, and disturbances to these factors correspond to the major aggregate shocks to the economy such as demand or supply shocks SW2006. This idea has a long tradition in macroeconomics . One example is SW2002a, which indicates the estimated factors can be interpreted as the diffusion indexes developed by NBER business cycle analysts to measure common movement in a set of macroeconomic variables.
Despite of the advantages of FA approach, there are a few drawbacks as well. First of all, it lacks explanation on the interrelationship among different economic variables, thus it cannot identify which predictors influence the target variable. Secondly, the estimated factors only capture the variation of major economic aggregates, but may lose information that is contained in a few predictors but beyond major economic aggregates. More importantly, the commonly used FA approach SW2002a,SW2002b does not take into account the target variable when estimating the factors, which means the factors used in the forecasting equation are the same no matter which target variable is being forecasted.
The second direction of many-predictor forecasts is to directly select the best predictors and their lagged values to carry out forecast, and we call it variable selection (VS) approach. This direction implies another rationale of the economic data -- given the selected predictors in the forecasting equation, all others have insignificant prediction power to the target variable anymore. To be noticed, it does not means the unpicked predictors are irrelevant or independent to the target variable. The forecasting equation derived from VS approach indicates that, conditional on the selected predictors, the remaining predictors have little prediction power on the target variable. VS is not a new topic, but it has not drawn as much attention as FA in economic forecasting. In contrast, VS has substantial development in other fields such as statistics, computer science, and bioinformatics, and impressive new methods and applications keep coming forward.
The first goal of this paper is to review several cutting-edge VS methods, and compare their performance with FA approach in the context of economic forecasting. One advantage of VS is its capability of interpreting the individual impact of each predictor on the target variable, including both direction and magnitude. This is helpful for understanding the interrelationship among different economic variables. More importantly, VS can select predictors that may contain useful information beyond the aggregate economic activity explained by the factors in FA approach. Only including the important predictors will avoid overfitting issue, thus enhance prediction power. In our empirical studies in Section (ref), we apply both FA and VS approaches to forecast three important macroeconomic variables -- Employment (EMP), Industrial Production (IP) and Consumer Price Index (CPI), and find that some VS methods achieve considerable improvement over FA approach for EMP and CPI. Also the relationship between the target variable and the selected predictor can be well explained by economic theories.
The second goal of this paper is to evaluate several groups of VS methods, including both classical procedures and cutting-edge algorithms, in terms of both variable selection accuracy and out-of-sample forecasting. Due to the huge body of VS literature, it is impossible to do an exhaustive review for VS methodologies. For this paper, we only focus on the high dimensional regime (dimension of predictors is larger than the number of observations), which is the case of economic forecasting. We pick the following four groups of methods: (1) classical procedure (forward selection); (2) $l_1$ regularization (adaptive LASSO); (3) gradient descent algorithms with sparsification (iterative hard thresholding and thresholding pursuit); and (4) a meta-heuristic algorithm called sequential Monte Carlo (SMC) proposed by Duan2019. We do not consider any machine learning algorithm such as random forest and neural network, due to their lack of interpretability. All these VS methods are applied in the framework of linear regression. Their performance in both variable selection and out-of-sample forecasting are examined through several simulation studies. The results show that, SMC, the most time-consuming algorithm, works the best across all simulation settings. Surprisingly, the performance of the classical forward selection matches up to the SMC and better than other advanced modern methods ($l_1$ regularization and gradient descent algorithms with sparsification).
The remainder of the paper is structured as follows. Section (ref) introduces the setting of economic forecasting with many predictors, and the implementation of FA and VS approaches. Section (ref) briefly reviews the four groups of VS methods, including their methodologies, advantages and disadvantages. Several simulation studies are carried out in Section (ref) to evaluate the performance of the four groups of VS methods. In Section (ref) we apply these VS methods on economic forecasting using the FRED-MD database (FREDMD), and compares their forecasting performance with that of FA approach. Conclusions and discussions are presented in Section (ref).
\paragraph{Notations} Throughout this paper, bold letters denote vectors, unbold letters denote scalar quantities and calligraphy letters denote matrices. $\mathbf{0}$ and $\mathbf{1}$ stand for a vector of zeros and ones respectively. For a $p$-dimensional vector $\boldsymbol{v}=(v_1,\dots,v_p)'$, we use $\left\Vert\boldsymbol{v}\right\Vert_0=\sum_{i=1}^p I(v_i\neq 0)$ with $I(.)$ being the indicator function, $\left\Vert\boldsymbol{v}\right\Vert_1=\sum_{i=1}^p|v_i|$ and $\left\Vert\boldsymbol{v}\right\Vert_2=\sqrt{\sum_{i=1}^p|v_i|^2}$ to denote the $l_0$ norm, $l_1$ norm and $l_2$ norm of $\boldsymbol{v}$ respectively. $\operatorname*{\text{supp}}(\boldsymbol{v})$ denotes the support of $\boldsymbol{v}$, i.e. $\operatorname*{\text{supp}}(\boldsymbol{v})=\{$indices of nonzero elements in $\boldsymbol{v}\}$. For a set $U$, $|U|$ denotes its cardinality, i.e. the number of elements in $U$. $\boldsymbol{v}_U=[v_i]_{i\in U}$ stands for a sub-vector of $\boldsymbol{v}$ whose indices of elements belong to $U$. For every iterative algorithm, we use superscript $(r)$ to stand for the $r$-th iteration.
In this section we first describe the setting of economic forecasting with many predictors, including notations and assumptions, then we outline FA and VS approaches under this framework.
We adopt the notations and assumptions per usual in economic forecasting literature SW2002a,SW2006,Ng2013. Let $y^h_{t+h}$ be the $h$-step ahead value of the variable to be forecasted. For example, in Section (ref) we consider forecasts of 1, 3, 6 and 12-month growth of the Employment (EMP). Let $EMP_t$ denote the value of EMP on month $t$. Then the $h$-month growth of EMP, at an annual rate, is
Let $\boldsymbol{Z}_t$ be the $n$-dimensional vector of predictor variables, which also includes the current value of the target variable $y_t$. In economic forecasting, both $y^h_{t+h}$ and $\boldsymbol{Z}_t$ are required to be stationary. This is accomplished by suitable preliminary transformations which are determined by a combination of statistical tests and expert judgment. For instance, unit root tests indicate that the logarithm of industry production (IP) series (denoted as $\{IP_t\}$) has a unit root. Therefore, its appropriate transformation is taking the log first difference, i.e. the corresponding predictor variable is $z_t=\log(IP_{t}) - \log(IP_{t-1})$. After transformation, each predictor $z_t$ is standardized to have mean zero and sample variance one. This standardization is required for FA approach and some VS methods.
FA approach assumes the predictor variables admit the following factor model representation:
where $\boldsymbol{F}_t$ is the $s \times 1$ latent factors, $\boldsymbol{\Lambda}$ is the $n\times s$ matrix of factor loadings, and $\boldsymbol{e}_t$ is the vector of idiosyncratic components satisfying $E(\boldsymbol{e}_t|\boldsymbol{F}_t) =\mathbf{0}$ and finite second moments. Here the latent factors $\boldsymbol{F}_t$ are estimated by the principle component analysis. SW2002a has proved that the principal component estimator is point-wise (for any date $t$) consistent and has limiting mean squared error (MSE) over all $t$ that converges to 0, under a suitable set of identifiability conditions. If some series contain missing values, the expectation-maximization (EM) algorithm described in SW2002a is utilized to estimate factors $\boldsymbol{F}_t$. After the factors have been estimated, the $h$-step ahead forecast is the linear projection of $y^h_{t+h}$ onto the $t$-dated factors, $y_t$ and their lagged values:
Here $q$ is the auto-regressive order, $m$ is the order of lagged factors and $\boldsymbol{f}_t$ is a vector of first $d$ factors in $\boldsymbol{F}_t$. In practice, $q$, $m$ and $d$ can be selected by some criteria, such as Schwarz’s BIC schwarz1978 and forward cross validation (FCV).
Note that factor model ((ref)) only includes the current value of predictor variables ($\boldsymbol{Z}_t$) without considering their lagged values. The historical information of $\boldsymbol{Z}_t$ are incorporated to forecasting through the lagged value of factors ($\boldsymbol{f}_{t-l}$, $l=1,\dots,m-1$). In our empirical study, we have 128 monthly economic time series, thus the vector $\boldsymbol{Z}_t$ for FA approach is 128-dimensional. However, the predictors used in VS approach are $\boldsymbol{X}_t=\left( \boldsymbol{Z}'_t,\boldsymbol{Z}'_{t-1},\dots,\boldsymbol{Z}'_{t-5}\right)'$, assuming that the lag order is five. Then the dimension of $\boldsymbol{X}_t$ is larger than the number of observations. In order to distinguish the predictors in these two different approaches, we use $\boldsymbol{Z}_t$ to denote the predictor variables in FA approach, and use $\boldsymbol{X}_t$ to denote the predictors in VS approach, respectively.
We apply various VS methods on economic forecasting through the following linear regression:
where $\boldsymbol{X}_t=(X_{1t},\dots,X_{pt})'$ is a $p-$dimensional vector of the predictors. In our real data application, $p=6n=768$. and $\boldsymbol{\beta}=(\beta_1,\dots,\beta_p)'$ is the vector of regression coefficients. The matrix form of ((ref)) is as follows:
Here we use $\boldsymbol{x}_i$ ($i=1,\dots,p$) to denote each column of the model matrix $\mathcal{X}$, i.e. $\mathcal{X}=\left [ \boldsymbol{x}_1,\dots,\boldsymbol{x}_p \right ] $. As mentioned in Subsection (ref), all $\boldsymbol{x}_i$'s are standardized to have mean zero and sample variance one.
The center part of the VS approach is the best subset problem with subset size $k$, which is given by the following optimization:
Here the $l_0$ norm of $\boldsymbol{\beta}$ (i.e. $\left\Vert\boldsymbol{\beta}\right\Vert_0$) counts the number of nonzeros in $\boldsymbol{\beta}$, which is bounded by $k$. Let $\hat{\boldsymbol{\beta}}_k$ be the optimal solution of ((ref)), then the support of $\hat{\boldsymbol{\beta}}_k$, denoted as $\hat{U}_k:=\left \{ i: \hat{\beta}_i\neq 0 \right \} $, is the best subset of predictors with size $k$. In practice, the subset size $k$ can be determined by AIC Akaike1974, BIC, FCV or other criteria.
The discrete nature of cardinality constraint ($\left\Vert\boldsymbol{\beta}\right\Vert_0\leq k$) poses a great difficulty in finding the global optimum. It requires comparison of all $\binom{p}{k}$ subsets of predictors, which is infeasible for large $p$. To the best of our knowledge, computing the optimal solution to problem ((ref)) is in general deemed as intractable. However, the last few decades have seen a flurry of activity in developing algorithms trying to solve ((ref)) at reasonable time cost, with associated optimality under certain conditions. In Section {(ref)}, we will review four groups of VS methods that try to obtain the good sub-optimal solution more efficiently.
As there is a vast literature on this topic, we present a selective overview. We select the following four types of VS methods: (1) classical greedy procedures, (2) $l_1$ regularization methods, (3) gradient descent algorithms with sparsification, and (4) meta-heuristic algorithms. The first two groups have already been investigated in many econometric literature, therefore their introduction are relatively concise. The last two groups are proposed in computer science and mathematical optimization but have not been seen wide adoption in economic forecasting. Thus these two groups will be introduced more elaborately. For each group, we mainly focus on the methods applied in our empirical study, presenting their ideas, advantages and disadvantages. The algorithm details can be found in Appendix (ref).
Classical VS procedures such as forward selection (FS), backward elimination (BE), and stepwise regression (SR) are available in many statistical software packages. These algorithms are greedy algorithms, which follow the heuristic of making the locally optimal choice at each iteration with the intent of finding a global optimum. For example, when adding a new predictor to the model, FS selects the one which maximize the decrement of sum of square errors (SSE) given the predictors already included in the model, until the model has $k$ predictors in total. As a result, such greedy search only examines a small portion of possible subsets of predictors and may be trapped in a local solution thus cannot guarantee to obtain the global optimum.
Due to the nature of greedy algorithms, the classical procedures seem to be inferior to modern VS algorithms such as LASSO. However, the latter also rely on certain assumptions to achieve global optimality, and these assumptions are hard to be verified in practice. In addition, FS and BE do not involve any tuning parameters, while more advanced algorithms are usually sensitive to the choice of tuning parameters and initial values. Thus, it is meaningful to compare classical greedy procedures with the modern algorithms in terms of their empirical performance. Moreover, another advantage of FS and BS is their fast update of model estimation when adding or deleting a predictor, which can avoid calculating the inverse of $\mathcal{X}'\mathcal{X}$. This fast updating algorithm and the detailed procedures of FS and BE are presented in Appendix (ref).
Since it is hard to directly solve the $l_0$ constraint optimization ((ref)), a group of methods pursue another direction of VS -- $l_1$ regularization, which is a convex relaxation of the $l_0$ constraint. One representative is LASSO Tibshirani1996:
As a surrogate for problem ((ref)), $l_1$ regularized methods produce a sparse estimation by shrinking many coefficients toward zero. There has been a large amount of work on this topic in terms of algorithms, theoretical properties and real world applications. Readers can refer to the books buhlmann2011statistics, hastie2009elements, wainwright_2019 and references therein. $l_1$ regularization methods enjoy several attractive properties. The first advantage of $l_1$ regularization, which is also an important reason to its popularity, is great computational efficiency. The problem ((ref)) is a convex quadratic optimization and there are several efficient algorithms to solve it. For example, pathwise coordinate optimization friedman2007 can compute the solution path at the same cost as a least squares calculation. Second, under some conditions it can be shown that LASSO can recover the true sparseness of $\boldsymbol{\beta}$ and deliver good prediction. However, the sufficient conditions to achieve this buhlmann2011statistics are related to the model matrix $\mathcal{X}$, thus are difficult to be verified in practice.
As argued in SCAD, the original LASSO ((ref)) leads to biased coefficient estimates. To address this shortcoming, several approaches such as adaptive LASSO adaptive and SCAD SCAD are proposed. It has been shown that their estimates possess the so-called oracle properties: (1) identifies the true set of predictors asymptotically; (2) has the optimal convergence rate. In our simulation study and real data application, we choose adaptive LASSO as a representative of $l_1$ regularization approaches, with objective function shown as follows:
Here, $w_i$ ($i=1,\dots,p$) is the penalty weight of $\beta_i$ that can be derived from $\sqrt{T}$-consistent estimator such as LASSO estimator ($w_i=1/|\hat{\beta}^{LASSO}_i|$) or ridge estimator ($w_i=1/|\hat{\beta}^{ridge}_i|$).
In spite of the good properties mentioned above, $l_1$ regularization is a relaxation of VS problem, thus do not provide provably optimal solution to ((ref)). In addition, as shown in our simulation studies, for very sparse cases, $l_1$ regularization tends to introduce a lot of spurious variables. As a consequence, it may be inferior to other types of methods in terms of identifying the important predictors.
In the fields of computer science and signal processing, a group of algorithms, referred to gradient decent with sparsification (GDS) algorithms, are developed to directly provide a good solution to the $l_0$-constraint optimization ((ref)). As the name suggests, these methods are extensions of gradient descent algorithms which impose sparsity on the parameter estimate. These algorithms include but not limited to iterative hard thresholding BLUMENSATH2009, compressive sampling matching pursuit CoSaMP2009, subspace pursuit subspace2009, hard thresholding pursuit HTP2011 and orthogonal matching pursuit with replacement orthogonal2011. The analysis of their theoretical properties in high dimensional regression setting can be found in IHTtheory2014, Bertsimas2016 and references therein. In general, this group of methods aim to minimize a loss function $f(\boldsymbol{\beta})$ subject to the $l_0$ constraint:
In our case, $f(\boldsymbol{\beta})$ equals to $\left\Vert\boldsymbol{y} - \beta_0\mathbf{1} - \mathcal{X}\boldsymbol{\beta}\right\Vert^2_2$, the SSE of linear regression. Without considering the $l_0$ constraint, gradient descend algorithm minimizes $f(\boldsymbol{\beta})$ iteratively by updating $\boldsymbol{\beta}$ as:
where superscript $(r)$ stands for the $r$th iteration, $\nabla f(\boldsymbol{\beta}^{(r)})=2\mathcal{X}'\mathcal{X}\boldsymbol{\beta}^{(r)}-2\mathcal{X}'\boldsymbol{y}$ is the gradient of $f(\boldsymbol{\beta})$ at $\boldsymbol{\beta}^{(r)}$, and $\eta$ is the step size. Since a $k$-sparse vector is desired, GDS-type algorithms modify the updating equation ((ref)) in some ways so that $\boldsymbol{\beta}^{(r+1)}$ becomes a $k$-sparse vector. In the following, we will introduce four GDS-type algorithms.
Iterative hard thresholding (IHT), also known as projected gradient descent, modifies the updating equation by adding a hard thresholding operator $H_k(.)$ to it:
For any vector $\boldsymbol{v}$, the hard thresholding operator $H_k(\boldsymbol{v})$ keeps the $k$ largest (in magnitude) elements of $\boldsymbol{v}$ and set the rests to zero. This operator is the simplest way to obtain a $k-$sparse vector. It has been also proven that for any arbitrary vector $\boldsymbol{v}$, $H_k(\boldsymbol{v})$ is the closest $k$-sparse vector to it in $l_2$ distance.
IHT directly converts $\boldsymbol{\beta}^{(r)} - \eta\nabla f(\boldsymbol{\beta}^{(r)})$ into a sparse vector, while some other methods, such as compressive sampling pursuit (CoSaMP) and subspace pursuit (SP), chase a good support of $\boldsymbol{\beta}$ based on the gradient descent and then finds the best fit within this support. Specifically, the CoSaMP algorithm iterates the following three-step scheme:
SP algorithm is similar to CoSaMP, except that it replaces $2k$ with $k$ in the first step of CoSaMP and modifies the third step:
Unlike IHT, CoSaMP and SP have one and two OLS calculations respectively. Since OLS is the most time consuming part in each iteration, these two methods have more computational costs than IHT within one iteration. However, on the other hand, OLS offers the best fit within the proposed support, which may make these two methods converge with less iterations.
HTP2011 combined the hard thresholding operator in IHT and the idea of pursuing a good support of $\boldsymbol{\beta}$ (CoSaMP and SP) into the following hard thresholding pursuit (HTP) algorithm:
The difference between HTP and CoSaMP/SP lies on their ways of proposing the support of $\boldsymbol{\beta}$. HTP uses the support of IHT result while CoSaMP and SP derives the support by the result of previous iteration $\boldsymbol{\beta}^{(r)}$ and the corresponding gradient $f(\boldsymbol{\beta}^{(r)})$.
In general, GDS-style algorithms are not limited to the aforementioned four algorithms, and there is no guarantee that one algorithm outperforms the others. In this paper, we select IHT and HTP as representatives of this group of methods. The initial input $\boldsymbol{\beta}^{(0)}$ should be a $k$-sparse vector, typically $\boldsymbol{\beta}^{(0)}=\mathbf{0}$. The iteration is stopped when the difference between $\boldsymbol{\beta}^{(r+1)}$ and $\boldsymbol{\beta}^{(r)}$ is small enough or the maximum number of iterations is reached.
Optimization ((ref)) is a combinatory optimization problem whose search space is discrete. Some meta-heuristic algorithms, such as simulated annealing SA1983,Cerny1985 and genetic algorithms goldberg2006genetic, are dedicated to solve such combinatory optimization, and become appealing to variable selection chatterjee1996genetic,brusco2014,brooks2003. Unlike greedy algorithms that always reject worse solutions, meta-heuristic algorithms accept worse solution with a probability to yield a more extensive research. In this subsection, we introduce one meta-heuristic algorithm proposed by Duan2019, called sequential Monte Carlo (SMC). This algorithm incorporates the idea of simulated annealing and Monte Carlo methods to considerably extend the searching space. Even the key idea involves Monte Carlo, SMC algorithm is different from the Bayesian variable selection methods that assume hierarchical distribution for the data and variables are selected by imposing some spike-and-slab priors on the model parameters. Indeed, SMC approach is not a Bayesian method. It does not require any distributional assumption and solves the $l_0$ constraint optimization ((ref)) in its original form.
For any nonempty set of indices $U$, let $\mathcal{X}_U=[\boldsymbol{x}_i]_{i\in U}$ be the sub-matrix of $\mathcal{X}$ whose columns belong to set $U$, and let $\boldsymbol{\beta}_U=[\beta_i]_{i\in U}$ be the corresponding sub-vector of coefficients. Then the $l_0$ constraint optimization ((ref)) is rewritten as a minimization problem:
where $\hat{\beta}_0$ and $\hat{\boldsymbol{\beta}}_U$ are the OLS estimates corresponding to $\mathcal{X}_U$. Let $\mathbb{U}_k$ be the collection of all the permutations of $k$ indices, i.e. $\mathbb{U}_k=\{U\subseteq\{1,\cdots,p\}: |U|=k\}$. Duan2019 assigned a discrete distribution function $f(U)$ on $\mathbb{U}_k$:
Then finding the global optimum of ((ref)) is done through generating a representative sample from this distribution. Obviously, among all permutations in $\mathbb{U}_k$, the optimal one corresponds to the peak of this distribution thus is more likely to be generated. Distribution $f$ is defined over permutations instead of combinations because permutations are easier to sample. Since this distribution has no tractable analytical solution and is multimodal, it is hard to directly construct a suitable proposal distribution for it. SMC applies the idea of distribution tempering that starts from an easy-to-sample distribution $f_0$, moves “smoothly" to the complex target distribution $f$ by composing a sequence of artificial intermediate distributions:
The distribution sequence $\{f_j\}_{j=1}^J$ is called “distribution-tempering bridge" and sequence $\{\gamma_j\}$ satisfies $0=\gamma_0<\gamma_1<\cdots<\gamma_J=1$. Clearly, $\gamma_0=0$ corresponds to the initial distribution $f_0$, and $\gamma_n=1$ corresponds to target distribution $f$. For each round $j$, the choice of $\gamma_j$ is self adapted by the algorithm, which makes the difference between $f_{j-1}$ and $f_j$ small enough so that the former distribution, $f_{j-1}$, can be a good proposal for the latter distribution, $f_j$.
The initial distribution, $f_0$, takes into consideration the prediction power of each individual predictor. Let $R^2_i$ be the $R^2$ of a single linear regression with the $i$-th predictor ($i=1,\dots,p$). The initial distribution, $f_0$, is the random sampling without replacement based on inclusion probability $q_i=R^2_i/\left(\sum_{i=1}^pR^2_i\right)$:
After obtaining the initial sample, denoted by $\Omega_0=\{U_{0,m}\}_{m=1}^M$ with $M$ being the sample size, we can generate a representative sample for the next distribution $f_1$ by following the three-step scheme:
By repeating this three-step scheme until $\gamma$ reaches 1, we finally arrive at a representative sample for the target distribution $f$. At last, a $k$-fold duplication technique Duan2015 is carried out to enlarge the sample size $k$ times, and the best subset of predictors is the $U$ with the smallest SSE in the final sample. The detailed algorithm can be found in Appendix (ref), and readers can refer to Duan2019 for its theoretical properties.
We present a variety of computational experiments to: (1) evaluate different types of VS methods in terms of variable selection accuracy and out-of-sample prediction; (2) compare FA approach to VS approach under the framework of time series forecasting. These two goals are addressed in Section (ref) and (ref) respectively.
We conduct three simulation studies to compare among different VS methods in terms of variable selection accuracy and out-of-sample prediction. For each group of VS methods, we select one or two representatives: (1) FS; (2) adaptive LASSO (adaLASSO); (3) IHT and HTP algorithms; (4) SMC algorithm. In our simulation studies, we investigate the $p>T$ case, thus BE is not applicable and excluded. In order to examine their performance, we consider five criteria:
Here, “True Positive" stands for the number of predictors that are correctly identified, “False Positive" is the number of predictors which are wrongly selected, and “False Negative" counts the number of predictors belonging to the true model but missed by the VS method. In this way, recall is the proportion of predictors in the true model that correctly identified, while precision is the proportion of the selected predictors that are truly significant. These two criteria quantify different aspects of the variable selection accuracy, and higher value indicates higher accuracy. However, these is a trade-off between precision and recall: as more predictors are selected, precision tends to increase while recall would decrease. The criteria, DC, can balance this trade-off and serves as an overall measurement of variable selection accuracy. Since this paper aims to apply VS methods on economic forecasting, we also include the MSPE as the criterion to evaluate out-of-sample prediction.
Simulation is carried out under the linear regression framework. The true model is:
The regressors are selected from $p$ potential predictors with $p\gg K$. All the potential predictors are generated from normal distribution with mean zero and variance one. Note that in this subsection we do not consider temporal dependence thus the observations are independent across time. In total, we construct three settings to evaluate VS methods under different scenarios. The first one adopts the simulation setting of Duan2019, which presents moderate high dimensional case ($p=900$ and $T=200$). The other two settings are cases of ultra-high dimension ($p=2000$, $T=50,\;100,\;200$), with setting 2 imposing independent structure and setting 3 imposing correlated structure among predictors.
All simulations are conducted with 100 repetitions. For FS, IHT, HTP and SMC, the tuning parameter is the the subset size $k$, while for adaLASSO, the tuning parameter is the penalty factor $\lambda$. Both of them are selected by 5-fold cross validation.
Table (ref) lists the five quantiles of the number of selected predictors by each VS method in setting 1, and the mean of $R^2$ in 100 repetitions. Recall that the true model contains 12 predictors by design. Clearly adaLASSO tends to pick too many predictors, with median of 38 and 34 under the two scenarios of theoretical $R^2$ respectively. Moreover, the range is much wider than those from other methods, reflecting its low stability in selecting predictors. This is not surprising since adaLASSO fails to distinguish between a zero and a nonzero coefficient when the corresponding predictors are correlated.
In contrast to adaLASSO, the other four methods are more conservative and tends to select less predictors than the true number 12. For example, FS yields a median of 8 predictors and covers the range from 5 to 13. Such under-selection is actually expected, because cross validation is a conservative way to find the correct number of predictors by avoiding over-fitting. In addition, the under-selection does not ruin the model fitting performance that is measured by in-sample $R^2$. The $R^2$ obtained by FS, SMC, IHT and HTP are close to the theoretical $R^2$. To better understand the reason of under-selection, we calculate the hit ratio (the proportion of successfully identifying the predictor in 100 repetitions) of each predictor in the true model, as shown in Table (ref). It appears that the predictors with small magnitudes of coefficient (0.1 and 0.4) would be excluded from the model due to their relatively low prediction power, especially when they are highly correlated with other predictors and the sample size ($T=200$) is not very large. It is also natural to see that the hit ratio is higher when the signal strength is higher (theoretical $R^2$=0.8) and the predictors are less correlated.
Table (ref) summarizes the results of the five criteria in setting 1. When the signal is strong (theoretical $R^2=0.8$), SMC, adaLASSO and FS are the best group in terms of prediction performance, which are significantly superior than IHT and HTP. As for the variable selection accuracy, SMC and FS outperforms the rests. Since adaLASSO suffers from excessive over-selection, its precision and DC are the lowest while its recall is the highest. HTP is a little better than IHT but still worse than SMC and FS. As for time cost, SMC is the most expensive algorithm, which takes about three hours to finish one simulation while others cost less than one minute. FS and adaLASSO are the most time-efficient, which only take two seconds or less to finish one simulation. In summary, the most sophisticated algorithm, SMC, is the best in terms of both variable selection accuracy and out-of-sample prediction but extremely time consuming. FS, the classical greedy algorithm, obtains very similar performance as SMC while it is much more computationally efficient. If we only focus on prediction, adaLASSO is one of the best choices given its great efficiency and promising prediction results.
In the scenario of theoretical $R^2=0.5$, the differences in prediction among all five methods are quite small, while adaLASSO outperforms the rests a little bit. As for variable selection accuracy, adaLASSO is still inferior in terms of DC and precision. SMC and FS are very close to each other and significantly better than IHT and HTP. This scenario asserts that, when the signal is not strong enough, adaLASSO is preferred for prediction while FS and SMC can provide more accurate variable selection results.
In setting 2 and 3, we only report the results of MSPE and DC in Table (ref) and (ref), and leave other criteria in Appendix (ref) (Table (ref) and (ref)). it can be found that all the VS methods work better in the independent case (setting 2) than in dependent case (setting 3). In addition, as sample size $T$ increases, their results become closer to that of the true model. In setting 2, $T=50$ is too small to successfully recover the true sparsity or obtain a satisfying prediction. When sample size increases to 100, SMC and FS have DC as 0.96 and 0.98 respectively which means they almost completely recover the true model, and their prediction are also quite close to the true model. While adaLASSO, IHT and HTP are still falling behind. When sample size increases to 200, except for adaLASSO which still over-selects predictors, all other VS methods are almost the same as the true model. In setting 3, since the potential predictors are correlated, VS methods can not provide satisfying performance until sample size increases to 200. Under this sample size, FS and SMC have the smallest prediction error and provide the most accurate models, followed by adaLASSO whose prediction is close to SMC and FS but performs worse in variable selection. In both settings, when the sample size is not large enough ($T=50$ in setting 2 and $T=50$ or 100 in setting 3), adaLASSO has the best prediction, followed by SMC, and SMC is superior than others in terms of identifying the true predictors.
We summarize our findings in these three simulation studies as follows. (1) When the signal is strong and the sample size is large enough, SMC and FS are the best among the five VS methods in terms of variable selection accuracy and prediction. However, SMC is very time consuming. The classical procedure FS is preferred if the time constraint is a concern. (2) adaLASSO is very good at prediction in all settings, especially when the signal is not strong or sample size is not large. However, it suffers from over-selection, and it is worse than others in terms of variable selection. (3) In general, IHT is slightly better than HTP, but both are significantly worse than SMC and FS. (4) When the sample size is too small or signal is not strong, none of these method can provide good performance, and their difference are insignificant.
In this simulation study, we investigate FA versus VS approaches under the framework of time series forecasting. The dimension of predictors and number of observations mimics the scales of our real data in Section (ref). The simulation setting is described as follows.
We generate four groups of predictors with group id as $i=1$, 2, 3 and 4. In each group, a series of latent factor $\{f_{it}\}$ is generated based on a AR(1) process, then 30 predictors are generated from this latent factor:
Clearly the series $\{x_{ijt}\}$ also follows a VAR(1) process with parameter $\phi_i$. Within one group, the temporal dependence is captured by the AR parameter $\phi_i$, while the cross-sectional correlation, cor$(x_{ijt},x_{ij't})$, is controlled by the relative magnitude between var($f_{it}$) and $\delta^2_{i2}$. By setting the values of $(\phi_i,\delta^2_{i1},\delta^2_{i2})$ as shown in Table {(ref)}, the four groups are constructed to have different temporal and cross-sectional dependence. Specifically, both temporal and cross-sectional dependence are strong in group 1. While for group 2, only the temporal dependence is strong but the cross-sectional dependence is weak. Group 3 has the opposite pattern to group 2, and in group 4, both temporal and cross-sectional dependence are weak.
For the response variable, we consider two generating mechanisms:
Setting1 indicates that $y_t$ is a sparse linear function of the predictors and their lagged values, which favors VS approach over FA. In contrast, setting 2 implies that $y_t$ is directly generated from the latent factors, thus FA is the correct approach while VS approach miss-specifies the model. In both settings, the theoretical $R^2$ is set to be 0.8.
In each setting, we generate 300 observations ($T=300$) and preserve the last 50 observations as test set to make prediction. We select FS, SMC, adaLASSO and IHT as representatives of VS approach and compare them with FA. For both FA and VS approaches, we do not include the lagged value of $y_t$ as predictor. Thus the models for FA and VS approaches are:
To predict $y_t$, the vector $\boldsymbol{X}_t$ in VS approach consists of $x_{ijt}$'s and their historical values up to lag 5, i.e. $\boldsymbol{X}_t=\{x_{11t},x_{11,t-1},\dots,x_{11,t-5},\dots,x_{4,30,t},x_{4,30,t-1},\dots,x_{4,30,t-5}\}$, with total of 720 predictors. While the vector $\boldsymbol{Z}_t$ in FA approach only contains $x_{ijt}$'s, i.e. $\boldsymbol{Z}_t=\{x_{11t},\dots,x_{4,30,t}\}$. The tuning parameters in VS methods are selected by FCV with validation set of size 50. For FA approach, we consider both BIC and FCV to determine the number of factors $d$ and the order of lagged factors $m$. The corresponding results are labeled as FA_BIC and FA_FCV respectively.
Tale (ref) shows the results of comparison between FA and VS approaches. Unsurprisingly, VS methods are uniformly better than FA in setting 1 and are inferior to FA in setting 2. For FA approach, there is little difference between the two criteria, BIC and FCV. While among the four VS methods, their comparisons are similar to those in the first simulation study: SMC and adaLASSO yield the most precise prediction followed by FS with no significant difference. As for variable selection accuracy, SMC and FS are the best while adaLASSO over-selects predictors. Based on this simulation study, we can conclude that, if only a handful of predictors are relevant to the response variable, VS methods are more advantageous than FA approach. However if there are many relevant predictors which possess a factor structure, FA approach is more suitable.
In this section, we apply both FA and VS approaches on forecasting of several important macroeconomic indices. By simulating real-time forecasting, we compare their empirical performance in real application. Our target indices are employment (EMP), industrial production index (IP) and consumer price index-all urban consumers (CPI). The predictors contains 128 economic variables in the FRED-MD dataset FREDMD, which is available online: \url{research.stlouisfed.org/econ/mccracken/fred-databases}.
We follow the framework of economic forecasting stated in Section (ref). Let $y^h_{t+h}$ be the $h$-month ahead value of the variable to be forecasted and $y_t$ be the corresponding $t$-dated value, which are required to be stationary. Following SW2002b and BaiNg2008, we define $y^h_{t+h}$ and $y_t$ as follows:
We carry out the transformation proposed by FREDMD to make the 128 economic variables in FRED-MD dataset stationary. Then each variable is standardized to have mean zero and variance one.
After the transformation and standardization, the 128 economic variables compose the input vector $\boldsymbol{Z}_t$ of factor model ((ref)) in the FA approach. As for the VS methods, the potential predictors contain both $\boldsymbol{Z}_t$ and its historical values up to 5 months, i.e. $\boldsymbol{X}_t=(\boldsymbol{Z}_t,\dots,\boldsymbol{Z}_{t-5})$, with dimension 768. Since both FA and VS approaches are under the framework of linear regression, their parameter estimation is sensitive to outliers. According to SW2002a,SW2002b, in a series $\{z_{it}\}$, an outlier is defined as an observation that deviates from the sample median by more than ten interquantile ranges. After identifying the outliers, FA approach replaces them with missing values (NA), and estimates the factors using EM algorithms to account for missing values. As for VS approach, we can either delete the series contaminated by outliers or impute the values by kalman filter. Therefore, we consider two strategies to account for outliers in the predictors:
For the target variable $y$, it does not contain any outlier, i.e. no value is beyond the ten interquantile ranges from the median. In our empirical results, generally FA and VS approaches perform slightly better under strategy 1. Thus we only report the results of strategy 1 in the main context, and put the results of strategy 2 in Table (ref) in Appendix (ref).
There are three tuning parameters in FA approach: number of factors $d$, lags of $y_t$, $q$, and lags of the factors $m$. We choose $(d, p, m)$ from $\{ 0 \le d \le 5, 0 \le q \le 6, 1 \le m \le 6 \}$ by FCV, where $d=0$ means no factor is included in the forecasting equation ((ref)) and $p=0$ means no auto-regressive term is included. The tuning parameter for VS is the number of selected predictors, which is also determined by FCV. In addition, we include univariate AR as benchmark model in the comparison:
where the AR order $\tilde{p}$ is selected from $0\leq \tilde{p} \leq 6$ by FCV.
The time span from 2015:1 to 2018:12 is reserved as test set to evaluate out-of-sample forecasting, with forecast horizon $h$ being 1, 3, 6 or 12 months. For both FA and VS approaches, the estimations and forecasts are conducted to simulate real-time forecasting using a rolling window scheme with window size 240. For example, if we want to forecast the $h$-month growth rate of EMP on 2015:1, The target variable is $y^h_{t+h}$ with $t\,$= 2015:1$\,-\,h$. Here 2015:1$\,-\,h$ stands for the month that is $h$ months before 2015:1. For example, if $h=3$, 2015:1$\,-\,h$ is 2014:10. To obtain its forecast, both FA and VS use data from 1995:1$\,-\,h$ to 2015:1$\,-\,h$ (240 months) to estimate model and select tuning parameters. The validation set for selecting tuning parameters is the sub-sample from 2011:1$\,-\,h$ to 2015:1$\,-\,h$ (48 months). When forecasting $y^h_{t+h}$ on 2015:2 ($t\,$= 2015:2$\,-\,h$), the tuning parameters are re-selected and the models are re-estimated using data from 1995:2$\,-\,h$ to 2015:2$\,-\,h$ (validation set is 2011:2$\,-\,h$ to 2015:2$\,-\,h$). In this way, we allow the parameter estimation and the optimal number of predictors and factors to change across $t$. For each pair of target variable and forecasting horizon $h$, a method produces 48 forecasts on the test set. We use MSPE of these 48 forecasts to evaluate the forecasting performance of this method. The ratio between MSPE of each method and that of the benchmark model ((ref)) is reported. Ratio less than one indicates the method has smaller forecasting error than the univariate AR.
The MSPE ratios are presented in Table (ref). For each forecast horizon $h$, we define the best group of results as those with MSPE ratios which do not exceeding 105% of the smallest one. These best ratios are marked in bold. For each target variable, the predictors\footnote{We only report the predictors selected with at least 12 times in the 48 rolling-window forecasts. The predictors with frequency less than 12 are viewed as lacking of systematic prediction power.} selected by the best VS method and their number of occurrences in the 48 forecasts are reported in Table (ref) in Appendix (ref).
For EMP, all the ratios are smaller than one, which manifests the usefulness of incorporating many predictors in forecasting. Among different approaches, SMC ranks the best for all the forecasting horizons. For $h=1$ and 3, the improvement of SMC over FA is not obvious, while for $h=6$ and 12, SMC reduces the MSPE of FA approach by 33% and 52% respectively.
The success of SMC implies EMP can be forecasted by only a few predictors. This point is also demonstrated in Table (ref), which shows that only two to five predictors are selected for each horizon. Specifically, $t$-dated EMP is selected in almost all 48 forecasts for all horizons, which implies the historical value has significant prediction power for EMP. For $h=6$ and 12, the linear forecasting model also includes two interest rates: 3-month treasury C minus FEDFUNDS and 3-month commercial paper minus FEDFUNDS. This leading effect of interest rate on EMP can be explained by economic theory. With low interest rate, consumers are more likely to consume now rather than wait for later. Low interest rate also drop the cost of borrowing to invest. Thus the increase in consumption and investment leads to higher demand for labor. As for one-month ahead forecast ($h=1$), except for the $t$-dated value of EMP, the model also includes M2 money stock, real personal consumption expenditures and two variables related to IP (IP:Fuel and IP:final products and nonindustrial supplies).
In terms of IP, neither FA nor VS approaches improve one-year-ahead forecasting ($h=12$) over AR, which implies other economic variables do not help long-term forecast of IP. While for short-term forecasts ($h=1$, 3 and 6), both FA and VS methods are better than AR, with VS method--FS--performs either the best or close to the best among all the methods. In addition, FS surpass FA by reducing MSPE by 18%, 6% and 6% for 1, 3 and 6-months ahead forecast respectively. Among the predictors selected by FS, IP:durable subcategory, S&P’s composite common stock:dividend yield and 3-month treasury C minus FEDFUNDS occurs the most often. Some other subcategories of IP are also frequently selected (more than 24 times). This leads us to conclude that, besides its own subcategories, the stock market and interest rate can forecast the movement of IP within six months.
The results of CPI show that, across all horizons, FA does not have significant advantage over AR, and even worse than it for one-year ahead forecast. This finding is consistent with the results in FREDMD. However, some of the VS methods do achieve substantial improvement over AR, especially the IHT algorithm. IHT performs the best among all the methods across all horizons, which reduces the MSPE of AR by 25%, 23% and 19% for 3, 6 and 12-months ahead forecasts respectively. Comparing FS with FA, the improvement is even more apparent: FS reduce the MSPE of FA by 23%, 27% and 36% for 3, 6, and 12-months ahead forecasts respectively.
The predictors selected by IHT are listed in Table (ref) in Appendix (ref). For short-term forecast ($h=1$, 3 and 6), only a handful of predictors are included in the model, and the predictors occurred most are real M2 money stock, $t$-dated CPI and its transportation subcategory, and oil price. This finding is consistent with economic theory and the composition of CPI. First, the money supply M2 has prediction power for CPI. This argument is supported by the classical quantity theory of money mill1965principles and evidenced by some empirical studies Bachmeier2005, BaiNg2008. Second, since CPI series has temporal dependence, it is natural to use current value ($t$-dated CPI) to forecast its future value. Third, transportation is the second largest category in CPI and is very sensitive to the oil price. Oil price also has direct impact on prices of many industrial materials which are the upstream prices of consumer price. In addition, oil price has great influence on other aspects of US economy such as stock market and investment. All these points make oil price an important index for CPI forecasting.
To summarize, some VS methods can provide better forecasts than FA approach, especially in the forecasting of EMP and CPI where SMC and IHT improve upon FA to a large extent. Moreover, for each target variable, the predictors selected by the best VS method are consistent with the underlying economic theory, which highlights the good interpretability of VS approach.
FA and VS approach indicate two different directions in economic forecasting. FA approach implies the target variable has many relevant predictors which can be explained by a few latent factors, while VS approach assumes only a handful of predictors have prediction powers on the target variable. Which approach is the best depends on the true data structure and the target variable to be forecasted. This paper aims to draw readers' attention on VS approach, which is less emphasized in the economic literature. In this paper, we introduce several cutting-edge VS algorithms to economic forecasting and compare to FA approach. It turns out for some target variables, VS approach is superior than FA approach. In particular, SMC significantly outperforms FA in forecasting of EMP, FS is superior than FA and AR for short-term forecasting of IP, and IHT is the best when predicting CPI. These methods also provide interpretable models which well explain the relationship between the target variable and the selected predictors. The second contribution of this paper lies on the overview and comparison among four different groups of VS methods. The last two groups, GDS-type algorithms and meta-heuristic algorithms, are popular in computer science but have not been widely used in economic forecasting. Several simulation studies are conducted to compare their prediction performance and variable selection accuracy. The most interesting finding is that, in all the simulation studies, the classical procedure FS works pretty well and sometimes even better than some advanced algorithms. Among all the reviewed VS methods, only the very time-consuming algorithm, SMC, is slightly but not significantly better than FS.
In economic forecasting, it is often the case that the underlying data structure is very complex. The relationships among different economic variables may vary in different time period. Therefore, it is unrealistic to expect any approach to be uniformly better than others. In the past two decades, several works have been developed to combine the ideas of FA and VS into one path, called supervised factor models (SFM), which includes targeted predictor approach BaiNg2008, supervised principle component analysis bair2006supervisedPCA and combining forecasts using principal components tu2019forecasting. These methods take into consideration the target variable when estimating the latent factors. For example, BaiNg2008 and tu2019forecasting first apply least angle regression to select relevant predictors and then construct factors only within the selected predictors. Most of existing methods only consider the commonly used $l_1$ regularization method, and little attention has been paid on $l_0$ constraint optimization algorithms such as the IHT and SMC algorithms. SFM under $l_0$ constraint optimization is an approach worthy for investigation and will be covered in a future paper.
No potential conflict of interest is declared by the authors.