Time Series Forecasting with Many Predictors

1. Introduction

In the current data-rich environment, many predictors are often available in time series forecasting. However, conventional methods have encountered serious difficulties in exploiting the information contained in such high-dimensional data. In particular, the curse of dimensionality often leads to unreliable forecasts when the conventional methods are applied naively. To exploit the information in the high-dimensional predictors, the use of latent factors has emerged among the first successful approaches. See, for instance, Pena and Box [1], Stock and Watson [2,3], Lam et al. [4], and Chapter 6 of Tsay [5]. These factor-based methods are also widely used in the econometric analysis of high-dimensional time series [e.g. 6,7,8,9].

Nevertheless, existing factor-based methods may still be inadequate, and several issues have been constantly observed in practice. For example, weak factors are prevalent in real-world data and extracting factors using all predictors may not be optimal [10]. In addition, some factor models may still contain many parameters when the factor dimension is not small. On the other hand, Kim and Swanson [11] reported, from an extensive experiment, that combining shrinkage methods with factor-based approaches can yield more accurate forecasts among many targeted macroeconomic time series. This implies that there might exist some irrelevant predictors in an empirical application. The recently proposed scaled principal component analysis of Huang et al. [12] and the supervised dynamic PCA (sdPCA) of Gao and Tsay [13] partially remedy this issue by constructing factors in a supervised fashion, so that the effects of irrelevant predictors are reduced. Limited experience indicates that when the number of predictors far exceed the sample size or when the predictors are highly correlated, the performance of such supervised approaches can still be compromised.

To efficiently extract predictive factors from high-dimensional data, in this paper, we propose a new method that blends variable selection in factor estimation. The proposed method uses the group orthogonal greedy algorithm (GOGA, Chan et al. 14) to screen variables that are relevant to the prediction problem at hand, and apply the sdPCA of Gao and Tsay [13] to extract factors from the selected variables. The orthogonal greedy algorithm, as well as its variant GOGA for grouped predictors, have been employed for variable selection for high-dimensional linear models, especially with dependent data and highly correlated predictors [14,15,16]. Since both GOGA and sdPCA serve to minimize the effects of noisy irrelevant variables, the combined procedure, which we call GO-sdPCA, can effectively construct highly predictive factors.

It is worth pointing out that the novelty of our method lies in a key step that successfully combines variable selection with factor-based methods. Indeed, variable selection techniques commonly employed such as the Lasso of Tibshirani [17] or the OGA of Ing and Lai [16] tend to select a sparse representation of the data. These methods will select only a few variables when many variables are driven by common factors and thereby highly correlated. For factor estimation, this can be undesirable since factor recovery typically benefits from the many variables that are loaded on the shared factors [see, e.g. [4,13]. This issue can be even more pronounced when the factors are strong, in which case the factors can be accurately recovered by employing many predictors. To circumvent this predicament, we propose to use a “peeling” technique, which repeatedly apply GOGA to the data after previously selected variables are dropped from the set of candidate predictors. In this way, our method can select more variables in the model and discriminate between relevant and irrelevant predictors. In Section 4, we intoduce the proposed method in more details. To better motivate the proposed method, we briefly review the orthogonal greedy algorithm (OGA) and its variants in Section 2. We also review the high-dimensional information criteria that are often used to balance the model complexity and the model fit along OGA iterations. The sdPCA is reviewed in Section 3, where we particularly emphasize the method’s design principles.

We use simulation studies and some empirical analyses to examine the performance of the proposed GO-sdPCA approach. In addition to comparing with some factor-based methods, such as the diffusion index approach of Stock and Watson [2,3], the sdPCA, and the time series factor model of Lam et al. [4], we also compare GO-sdPCA against the Lasso and the random forests of Breiman [18], both of which are widely-used and versatile tools in machine learning. We found that GO-sdPCA improves upon most of factor-based approaches and fares well with selected competing methods, even when the number of relevant variables is comparable with the sample size. The simulation results are discussed in Section 5.

In Section 6, we apply the proposed method to two real datasets. The first dataset is the well-known FRED-MD macroeconomic data [9]. Time series forecasting has been a vital topic both in econometric research and in policy-making pertaining to this data set. The other dataset contains the hourly particulate matters (PM_2.5) measurements in Taiwan. The PM_2.5 data play an important role in the environmental studies as well as informing public health policies. However, since its measurements are taken over time and space, the dataset is naturally of high dimension and presents serious challenges to forecast. Our findings demonstrate that the proposed method consistently offers more accurate forecasts than the competing methods, validating its practical utility.

2. Orthogonal Greedy Algorithm

In this section, we briefly review the (group) orthogonal greedy algorithm. The OGA is an iterative algorithm that sequentially chooses predictors to form a regression model. Theoretically grounded in approximation theory [19], the OGA is also easy to implement computationally. In the statistical learning and high-dimensional linear regression literature, Barron et al. [20] and Ing and Lai [16] analyzed its convergence rate as well as variable screening capability. Its generalization to grouped predictors was studied by Chan et al. [14] with application to threshold autoregressive time series models.

Throughout this paper, we denote by $\mathbf{y}={(y_1,\dots ,y_n)}^{\top }$ the data of the response variable we are interested in, where n is the sample size. We also have data for the p predictors, $\{{\mathbf{x}}_{\left(j\right)}:j=1,2,\dots ,p\}$, where ${\mathbf{x}}_{\left(j\right)}={(x_{1,j},\dots ,x_{n,j})}^{\top }$. The OGA is defined as follows. Starting with ${\widehat{J}}_0=\varnothing$ and ${\mathbf{u}}^{\left(0\right)}=\mathbf{y}={(y_1,y_2,\dots ,y_n)}^{\top }$, the OGA computes, at iteration $k=1,2,\dots$,

$$\begin{array}{cc}\hfill {\widehat{j}}_k& =arg\underset{1\le j\le p}{min}{\parallel {\mathbf{u}}^{(k-1)}-{\mathbf{x}}_{\left(j\right)}{\left({\mathbf{x}}_{\left(j\right)}^{\top }{\mathbf{x}}_{\left(j\right)}\right)}^{-1}{\mathbf{x}}_{\left(j\right)}^{\top }{\mathbf{u}}^{(k-1)}\parallel }^2,\hfill \end{array}$$

(1)

and updates

$$\begin{array}{cc}\hfill {\widehat{J}}_k& ={\widehat{J}}_{k-1}\cup \left\{{\widehat{j}}_k\right\},\hfill \\ \hfill {\mathbf{u}}^{\left(k\right)}& =(\mathbf{I}-{\mathbf{H}}_{\left(k\right)})\mathbf{y}.\hfill \end{array}$$

where ${\mathbf{H}}_{\left(k\right)}$ is the orthogonal projection matrix associated with the linear space spanned by $\{{\mathbf{x}}_{\left(j\right)}:j\in {\widehat{J}}_k\}$. Clearly, the OGA sequentially selects variable to include in the model; the set ${\widehat{J}}_k$ denotes the index set corresponding to the predictors already selected at iteration k. Intuitively, in each iteration OGA selects the variable that best explains the current residuals. There are numerical schemes to speed up the computation of the residuals such as using sequential orthogonalization. We refer to Ing and Lai [16] and Chan et al. [14] for details.

For the purpose of this paper, we will make use of a slight generalization of the OGA that deals with grouped predictors. Chan et al. [14] studied the group OGA (GOGA) and applied the method to estimate threshold time series models. Consider the j-th predictor, ${\left\{x_{t,j}\right\}}_{t=1}^n$. Instead of being a scalar predictor, it contains $d_j$ component predictors, $\{(x_{t,j,1},\dots ,x_{t,j,d_j}):t=1,2,\dots ,n\}$, thereby forming a “group” predictor. Then, by substituting

$$\begin{array}{c}\hfill {\mathbf{x}}_{\left(j\right)}=\left(\begin{array}{cccc}{x}_{1,j,1}& x_{1,j,2}& \dots & x_{1,j,d_j}\\ x_{2,j,1}& x_{2,j,2}& \dots & x_{2,j,d_j}\\ ⋮& ⋮& \ddots & ⋮\\ x_{n,j,1}& x_{n,j,2}& \dots & x_{n,j,d_j}\end{array}\right)\end{array}$$

in (1), the procedure becomes the GOGA.

After $K_n$ iterations, the GOGA selects the (group) predictors corresponding to ${\widehat{J}}_{K_n}=\{{\widehat{j}}_1,{\widehat{j}}_2,\dots ,{\widehat{j}}_{K_n}\}$. Bias-variance trade-off manifests when selecting the number of predictors to be included, as a large $K_n$ may lead to over-fitting. To select a desirable model complexity, Ing [15] suggests using the high-dimensional Akaike information criterion (HDAIC). The HDAIC of the model at iteration k is defined as

$$\begin{array}{c}\hfill \mathrm{HDAIC}\left({\widehat{J}}_k\right)=\left(1+C\frac{klogp}{n}\right){\widehat{\sigma }}_{\left(k\right)}^2,\end{array}$$

(2)

where ${\widehat{\sigma }}_{\left(k\right)}^2=n^{-1}{\parallel {\mathbf{u}}^{\left(k\right)}\parallel }^2$, and C is a constant to be tuned. Then the model selected by HDAIC is

$$\begin{array}{c}\hfill {\widehat{J}}_{\widehat{k}}=\{{\widehat{j}}_1,\dots ,{\widehat{j}}_{\widehat{k}}\},\mathrm{where} \widehat{k}=arg\underset{1\le k\le K_n}{min}\mathrm{HDAIC}\left({\widehat{J}}_k\right).\end{array}$$

(3)

Theoretically, Ing [15] proved that the resulting model selected by OGA+HDAIC adapts to the underlying sparsity structure. In particular, consider the high-dimensional regression model,

$$\begin{array}{c}\hfill y_t=\sum _{j=1}^p{\beta }_j{x}_{t,j}+{ϵ}_t.\end{array}$$

Then the conditional mean squared prediction error of OGA+HDAIC is of rate

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\left(\frac{logp}{n}\right)}^{1-1/2\gamma },\hfill & \mathrm{if}\sum _{j\in J}\left|{\beta }_j\right|\le D_1{\left(\sum _{j\in J}{\beta }_j^2\right)}^{(\gamma -1)/(2\gamma -1)}\mathrm{for} \mathrm{all}J\hfill \\ \frac{lognlogp}{n},\hfill & \mathrm{if}\sum _{j\in J}|{\beta }_j|\le D_2\underset{j\in J}{max}\left|{\beta }_j\right|\mathrm{for} \mathrm{all}J\hfill \\ \frac{k_{0}logp}{n},\hfill & \mathrm{if}\underset{j:{\beta }_j\ne 0}{min}\left|{\beta }_j\right|\ge D_3,\hfill \end{array}\right.\end{array}$$

for some $\gamma \ge 1$ and positive constants $D_1,D_2,D_3$. See Theorem 3.1 of Ing [15] for details. Note that these rates are minimax optimal, and, are automatically achieved by OGA+HDAIC without knowledge about the sparsity pattern. In this paper, we will leverage this property to select good predictors in constructing forecasts. In this way, the variables used in our method is carefully selected in a supervised fashion, which is more effective than employing all predictors in the high-dimensional data.

3. Supervised Dynamic PCA

Next, we review the supervised dynamic PCA (sdPCA) proposed by Gao and Tsay [13] for forecasting. The sdPCA is a factor-based forecasting approach, which took a major role in the literature of time series forecasting. Tailored to incorporate dynamic information, the sdPCA has been shown to outperform some existing factor-based approaches such as the diffusion index model of Stock and Watson [2,3], and the scaled PCA of Huang et al. [12].

We first outline the sdPCA procedure. Given a forecast horizon h, the sdPCA first constructs intermediate predictions using each individual predictor and its lagged values. For instance, one can regress $y_{t+h}$ on $x_{t,j},x_{t-1,j},\dots ,x_{t-q_2+1,j}$, where $q_2\in \mathbb{N}$ is a user-specified lag value, and obtain the fitted values

$$\begin{array}{c}\hfill {\widehat{\mu }}_j+\sum _{k=0}^{q_2-1}{\widehat{\gamma }}_{j,k}{x}_{t-k,j}:={\widehat{\mu }}_j+{\widehat{x}}_{t,j},\end{array}$$

(4)

where ${\widehat{\mu }}_j$ is the intercept estimate and ${\widehat{\gamma }}_{0,1},\dots ,{\widehat{\gamma }}_{q_2-1,j}$ are the slope estimates. For different j’s, the number of lags $q_2$ used in the regression can differ. For instance, $q_2$ can be selected by the BIC. Then, with the constructed intermediate predictions, ${\widehat{x}}_{t,j}$’s, the sdPCA apply PCA to estimate a lower dimensional factor ${\widehat{\mathbf{g}}}_t\in {\mathbb{R}}^r$, where $r<p$. Therefore, the data to which the PCA applies, ${\widehat{x}}_{t,1},\dots ,{\widehat{x}}_{t,p}$, are in the same unit (as $y_t$), which scales the variables according to their predictive power.

Finally, one may employ a linear model for the predictive equation,

$$\begin{array}{c}\hfill y_{t+h}\sim \widehat{\alpha }+{\widehat{\mathit{\beta }}}^{\top }{\widehat{\mathbf{g}}}_t,\end{array}$$

(5)

where $\widehat{\alpha }$ and $\widehat{\mathit{\beta }}$ are intercept and slope estimates respectively, and ∼ signifies that we run the linear regression while the underlying relationship of $y_{t+h}$ and ${\widehat{\mathbf{g}}}_t$ may not be exactly linear. Gao and Tsay [13] also suggests to use Lasso to estimate (5) instead of the usual linear regression if the number of common factors is large. Additionally, one can also include some lags of ${\widehat{\mathbf{g}}}_t$ is (5) and let Lasso perform variable selection.

The sdPCA has several advantages over the conventional PCA methods. First, the PCA is not scale-invariant. On the contrary, the sdPCA constructs predictors that are in the same unit, which naturally scales the predictors according to their predictive capabilities. Second, instead of performing PCA directly on contemporaneous data ${\mathbf{x}}_t$, the sdPCA sources from the lagged information in ${\mathbf{x}}_{t-l}$, $l=0,1,\dots ,q_2-1$, where ${\mathbf{x}}_t={(x_{t,1},\dots ,x_{t,p})}^{\top }$. For the conventional PCA to use the lagged information, one often needs to augment the data by appending the lagged variables so that PCA is performed on $({\mathbf{x}}_t,\dots ,{\mathbf{x}}_{t-q_2+1})$, which leads to even higher dimensionality. Lastly, the usual PCA is performed in an unsupervised fashion, whereas the sdPCA constructs the factors in a supervised fashion. While the usual principal component directions are not necessarily predictive of the response, factors extracted by sdPCA can potentially yield better forecasts. In fact, Gao and Tsay [13] have showed that sdPCA has a lower mean square forecasting error than the approaches of Stock and Watson [2,3] and Huang et al. [12].

In this paper, we employ the sdPCA to capitalize on the aforementioned properties. However, for noisy high-dimensional data, the performance of sdPCA may be severely compromised. Hence, careful dimension reduction before applying the sdPCA is desirable. In the next section, we describe the proposed procedure, which combines GOGA and HDAIC with sdPCA to improve the accuracy of prediction.

4. The Proposed GO-sdPCA

In this section, we introduce the proposed method, GO-sdPCA, which screens variables by the GOGA and HDAIC and then estimates factors by the sdPCA approach. In order to facilitate factor recovery, we apply a “peeling” technique to select more factor-relevant variables in the new procedure.

To tackle the difficulties encountered by the factor-based approaches when applied to high-dimensional data, our method begins with dimension reduction by employing the GOGA introduced in Section 2. Because of the serial dependence in the data, it is beneficial to select variables based not only on cross-sectional correlation but also on the lagged information. To this end, for each predictor ${\mathbf{x}}_j={(x_{1,j},\dots ,x_{n,j})}^{\top }$, we consider the group predictor,

$$\begin{array}{c}\hfill {\mathbf{x}}_{\left(j\right)}=\left(\begin{array}{cccc}{x}_{q_1,j}& x_{q_1-1,j}& \dots & x_{1,j}\\ x_{q_1+1,j}& x_{q_1,j}& \dots & x_{2,j}\\ ⋮& ⋮& \ddots & ⋮\\ x_{n-h,j}& x_{n-h-1,j}& \dots & x_{n-h-q_1+1,j}\end{array}\right),j=1,2,\dots ,p,\end{array}$$

(6)

where $q_1$ is a pre-specified integer for the number of lagged values to consider, and $h\in \mathbb{N}$ is the forecast horizon. Let

$$\begin{array}{c}\hfill \mathbf{y}={(y_{q_1+h},y_{q_1+h+1},\dots ,y_n)}^{\top }\end{array}$$

(7)

be the response vector. Then, we employ the GOGA algorithm in Section 2, with ${\mathbf{x}}_{\left(j\right)}$ and $\mathbf{y}$ substituted by (6) and (7), respectively. Since each group predictor already incorporates past information, the GOGA algorithm can select variables while accounting for the dynamic information.

We also employ the HDAIC (defined in (2)) after GOGA to prevent the inclusion of irrelevant variables. Since the above procedure depends on the number of GOGA iterations $K_n$, the HDAIC constant C in (2), the number of lags $q_1$ used in forming the group predictors, and the pool of candidate predictors $\mathcal{I}=\{1,2,\dots ,p\}$, we conveniently denote it as a set-valued function,

$$\begin{array}{c}\hfill \widehat{J}=\mathcal{A}(K_n,C,q_1,\mathcal{I}),\end{array}$$

(8)

where $\widehat{J}$ is the index set outputted by GOGA+HDAIC, as in (3). For completeness, we summarize the GOGA+HDAIC method introduced in Section (Section 2) in Algorithm 1.

Algorithm 1:GOGA+HDAIC ($\mathcal{A}$)

Input: Number of maximum iterations $K_n$, HDAIC parameter C, number of lags $q_1$, candidate set $\mathcal{I}$ Initialization: ${\mathbf{u}}^{\left(0\right)}=\mathbf{y}$; selected index $\widehat{J}=\varnothing$ 1  for $k=1,2,\dots ,K_n$ do $\phantom{(}$ 2  Select $$\begin{array}{cc}\hfill {\widehat{j}}_k& =arg\underset{j\in \mathcal{I}}{min}{\parallel {\mathbf{u}}^{(k-1)}-{\mathbf{x}}_{\left(j\right)}{\left({\mathbf{x}}_{\left(j\right)}^{\top }{\mathbf{x}}_{\left(j\right)}\right)}^{-1}{\mathbf{x}}_{\left(j\right)}^{\top }{\mathbf{u}}^{(k-1)}\parallel }^2,\hfill \end{array}$$ where ${\mathbf{x}}_{\left(j\right)}$ is defined in (6). 3  Update $\widehat{J}\leftarrow \widehat{J}\cup \left\{{\widehat{j}}_k\right\}$ 4  Update ${\mathbf{u}}^{\left(k\right)}=(\mathbf{I}-{\mathbf{H}}_{\left(k\right)})\mathbf{y}$, where ${\mathbf{H}}_{\left(k\right)}$ is the projection matrix associated with $\{{\mathbf{x}}_{\left(j\right)}:j\in \widehat{J}\}$ 5  end 6  Choose, as in (3), $$\begin{array}{c}\hfill \widehat{m}=arg\underset{1\le m\le K}{min}\mathrm{HDAIC}\left(\{{\widehat{j}}_1,\dots ,{\widehat{j}}_m\}\right),\end{array}$$ where HDAIC is defined in (2) with C set to the inputted value. Output: selected indices $\{{\widehat{j}}_1,\dots ,{\widehat{j}}_{\widehat{m}}\}$

   Because the orthogonal residuals are used in the greedy search, GOGA tends to differentiate highly correlated predictors and only selects those predictors that explain distinct (close to orthogonal) directions of the response. However, in factor estimation, the relevant variables loaded on the common factors tend to be highly correlated and failing to employ these correlated predictors may lose some statistical efficiency for the inference of the underlying factors (the blessing of dimensionality, [4,13]). Therefore, to encourage GOGA to screen the factor-relevant, correlated predictors, we propose a technique, referred to as “peeling,” that allows careful inclusion of more variables. The idea is to repeatedly apply GOGA+HDAIC, with variables selected from previous implementations discarded from the candidate set. Hence, in each iteration, GOGA+HDAIC is forced to select a new set of variables that best predict $\mathbf{y}$. Formally, the peeling algorithm is described in Algorithm 2. Let M be the number of peeling iterations. In the following, we denote the peeling procedure, which depends on M, $K_n$, and C, by the set-valued function $\mathcal{P}$. Namely,

$$\begin{array}{c}\hfill \widehat{Q}=\mathcal{P}(M,K_n,C,q_1).\end{array}$$

Algorithm 2:Peeling ($\mathcal{P}$)

Input: Number of peeling iterations M, number of GOGA iterations $K_n$, HDAIC parameter C, number of lags $q_1$ in group predictors Initialization: $\widehat{Q}=\varnothing$, $\mathcal{I}=\{1,2,\dots ,p\}$ 1  for $m=1,2,\dots ,M$ do 2  Run GOGA+HDAIC $$\begin{array}{c}\hfill {\widehat{J}}^{\left(m\right)}=\mathcal{A}(K_n,C,q_1,\mathcal{I})\end{array}$$ 3  Update $\widehat{Q}\leftarrow \widehat{Q}\cup {\widehat{J}}^{\left(m\right)}$ 4  Discard selected variables from the candidate set $\mathcal{I}\leftarrow \mathcal{I}-\widehat{Q}$ 5  end Output: selected indices $\widehat{Q}$

It is worth noting that peeling is different from simply running GOGA for many iterations. Although GOGA also selects distinct variables along its iteration because of orthogonalization (that is, GOGA does not revisit any previously selected variable), peeling, which discards previously selected variables from the candidate set, would produce very different results. For high-dimensional data, running many iterations of GOGA will lead to extremely small residuals, which may no longer carry sufficient variations to detect the remaining relevant variables (with finite sample). On the other hand, peeling re-starts GOGA in every iteration, using $\mathbf{y}$ instead of the previous residuals in the GOGA algorithm. Since GOGA is re-started with a smaller pool of candidate variables, the residuals in peeling will not shrink to zero after many (potentially more than n) variables are already contained in $\widehat{Q}$. We also remark that the idea of peeling is akin to the random forest which uses randomly selected variables in building each tree. Hence, peeling can be used to detect the weak predictors in the case where a few of the predictors are highly predictive and many others are only weakly predictive [see, e.g. Chpt. 8 of [21].

After variable selection by peeling, sdPCA is employed to estimate the factors ${\widehat{\mathbf{f}}}_t$ from the selected variables. Let $q_2$ be the number of lags used in constructing the intermediate predictions (4). Finally, the predictive model

$$\begin{array}{c}\hfill y_{t+h}=\sum _{k=1}^{q_3}{\alpha }_k{y}_{t-q+1}+{\mathit{\beta }}^{\top }{\widehat{\mathbf{f}}}_t+{ϵ}_{t+h},\end{array}$$

(9)

is estimated by OLS or Lasso, where $q_3$ is the number of autoregressive variables. The above procedure, which combines GOGA, peeling, and sdPCA, is called GO-sdPCA.

In closing this section, we briefly discuss the selection of the number of lags. In practice, the number of lags $q_1$ used in GOGA, $q_2$ used in the sdPCA step, and $q_3$ in the predictive equation can all differ. If $q_1$ is large, we can detect the predictors whose distant lags are predictive of the response. However, each step of GOGA iteration will then consume more useful variations, so GOGA will consequently select a smaller number of grouped predictors, which can be counterproductive for the subsequent factor estimation. Thus it is advisable to use a smaller $q_1$ in the screening step, which leads to more predictors selected, and use a larger $q_2$ in the sdPCA step to extract information from the past.

5. Simulation Studies

In this section, we assess the finite-sample performance of the proposed GO-sdPCA method via simulation studies. Some existing factor-based methods are employed as benchmarks, such as the diffusion index approach of Stock and Watson [2,3], the time series factor model of Lam et al. [4], and the supervised dynamic PCA of Gao and Tsay [13]. These benchmarks are referred to as SW, LYB, and sdPCA, respectively.

After factor estimation, the predictive model

$$\begin{array}{c}\hfill y_{t+1}=\sum _{k=1}^q{\alpha }_k{y}_{t-q+1}+{\mathit{\beta }}^{\top }{\widehat{\mathbf{f}}}_t+{ϵ}_{t+1},\end{array}$$

(10)

is estimated by OLS, where q is an integer specified later, and ${\widehat{\mathbf{f}}}_t$ is constructed using different approaches. For the proposed GO-sdPCA (shorthanded as GsP* hereafter), we set $C=2$ in the HDAIC definition (2) and use 10 peeling iterations. In each peeling iteration, GOGA is applied with $q_1$ set to 2 for the group predictors. Then, in the sdPCA step, r factors are extracted with $q_2=q$ lags used in constructing the intermediate predictions in (4). To demonstrate the usefulness of the peeling technique, we also consider implementing GO-sdPCA by naively combining GOGA+HDAIC with sdPCA. That is, we only run one peeling iteration and the variables selected are exactly the ones selected by GOGA+HDAIC. This method is denoted as GsP in the sequel. Similarly, the forecasts of sdPCA are constructed by estimating (10) with the factors estimated as in Section 3. The time series factors of Lam et al. [4] are estimated using the eigenanalysis of a non-negative definite matrix computed from the autocovariance matrices at nonzero lags. In our implementation, q lags of past predictors are used by LYB in the eigenanalysis. Finally, SW follows that of Stock and Watson [2,3] and uses PCA to extract the factors.

In addition, we employ some commonly used alternatives, including the Lasso of Tibshirani [17] and the random forests of Breiman [18] as competing methods. The Lasso is a versatile tool for building sparse linear regression models, while the random forest (RF) excels in capturing non-linear relationships. Recently, Chi et al. [22] investigated the asymptotic consistency of RF for high-dimensional data. See also Saha et al. [23] for application of RF to dependent data. The tuning parameters for the Lasso are selected by the BIC as suggested by Medeiros and Mendes [24], whereas we adopt the hyper-parameters for RF as recommended by the randomForest [25] package in R. Therefore, about one-third of the candidate variables is randomly selected at each split. For both methods, q lags of the dependent variable and the predictors are used for fair comparison.

5.1. Simulation Designs and Results

In the simulations, we consider three data generating processes (DGP) to generate the synthetic data. Throughout the experiments, we use one-step-ahead forecasts ($h=1$) where each method makes a forecast for $y_{n+1}$, which is not in the training sample. The root mean squared forecast errors, averaged over 500 Monte-Carlo simulations, are used for comparing different approaches.

DGP 1. 

Let ${\mathbf{f}}_t{\sim }_{\mathrm{i}.\mathrm{i}.\mathrm{d}.}N(0,{\mathbf{I}}_{r_{\mathrm{DGP}}})$, where $r_{\mathrm{DGP}}\in \mathbb{N}$ is the number of underlying factors. The predictors ${\mathbf{x}}_t\in {\mathbb{R}}^p$ are generated by

$$\begin{array}{c}\hfill {\mathbf{x}}_t=\mathbf{B}{\mathbf{f}}_t+2{\mathbf{\delta }}_t,\end{array}$$

where $\left\{{\mathbf{\delta }}_t\right\}$ are independent p-dimensional t-distributed random vectors with independent components and five degrees of freedom, and $\mathbf{B}\in {\mathbb{R}}^{p\times r_{\mathrm{DGP}}}$ has independent $\mathrm{Unif}(-2,2)$ entries, with $p-s$ rows randomly set to zero. That is, $\mathbf{B}$ only has s nonzero rows. Randomly generate ${\mathit{\beta }}_1={({\beta }_{1,1},\dots ,{\beta }_{r_{\mathrm{DGP}},1})}^{\top }$ and ${\mathit{\beta }}_2={({\beta }_{1,2},\dots ,{\beta }_{r_{\mathrm{DGP}},2})}^{\top }$ via ${\beta }_{j,1}\sim \mathrm{Unif}(1.0,2.5)$ and ${\beta }_{j,2}\sim \mathrm{Unif}(-2.0,-0.8)$. Finally,

$$\begin{array}{c}\hfill y_t=0.6y_{t-1}+0.2y_{t-2}+{\mathit{\beta }}_1^{\top }{\mathbf{f}}_{t-1}+{\mathit{\beta }}_2^{\top }{\mathbf{f}}_{t-2}+{ϵ}_t,\end{array}$$

where $\left\{{ϵ}_t\right\}$ are independent standard Gaussian.

In this DGP, the parameter s dictates both the number of relevant variables and the strength in recovering the factors. The larger s is, the stronger the factors are. Factor strength plays a critical role in factor recovery [2,3,4]. In practice, s is seldom known. Therefore, we will consider the cases $s\in \{0.25n,0.5n,0.75n\}$, where n is the sample size, to see whether the methods adapt well to various levels of factor strength.

Table 1 reports the root mean squared forecast error (RMSFE) averaged over 500 Monte Carlo simulations. Across all sparsity levels $s\in \{0.25n,0.5n,0.75n\}$, the proposed GsP* delivered the most accurate forecasts and the sdPCA ranks the second. This suggests that the proposed peeling procedure and GOGA improves accuracy in forecasting. The GsP, which naively combines GOGA+HDAIC with sdPCA, shows limited forecasting capabilities with RMSFE being higher than those of the Lasso. This indicates, again, the peeling strategy markedly improved the forecasting performance by selecting more factor-relevant variables. The other forecasting methods, including SW, LYB and RF, seem to suffer from the effect of employing many irrelevant variables in the high-dimensional data. We remark that DGP 1 is essentially a sparse model because $r_{DGP}$ used is relatively small. Therefore, it is not surprising to see that Lasso fares reasonably well.

DGP 2. 

In this DGP, ${\mathbf{x}}_t$ is generated via a vector MA(1) model:

$$\begin{array}{c}\hfill {\mathbf{x}}_t={\mathbf{\delta }}_t+0.8\mathbf{B}{\mathbf{\delta }}_{t-1},\end{array}$$

where $\mathbf{B}$ is a randomly drawn $p\times p$ matrix of rank $r_{\mathrm{DGP}}$. The coefficients ${\mathit{\beta }}_1={({\beta }_{1,1},\dots ,{\beta }_{p,1})}^{\top }$ and ${\mathit{\beta }}_2={({\beta }_{1,2},\dots ,{\beta }_{p,2})}^{\top }$ are randomly generated via ${\beta }_{j,1}\sim \mathrm{U}(1.0,3.0)$ and ${\beta }_{j,2}\sim \mathrm{U}(-2.5,-0.5)$. But, for a set of random indices J with cardinality equal to s (i.e., $♯\left(J\right)=s$), we set ${\beta }_{k,1}={\beta }_{k,2}=0$ for $k\notin J$. In this way, ${\mathit{\beta }}_1$ and ${\mathit{\beta }}_2$ share the same nonzero coordinates. Finally,

$$\begin{array}{c}\hfill y_t=0.6y_{t-1}+0.2y_{t-2}+{\mathit{\beta }}_1^{\top }{\mathbf{x}}_{t-1}+{\mathit{\beta }}_2^{\top }{\mathbf{x}}_{t-2}+{ϵ}_t.\end{array}$$

In this example, the relevant predictors have a direct impact on the response, instead of through any common factors. Additionally, when s is large, it is very difficult to recover the regression coefficients because of the lack of sparsity. Therefore, DGP 2 fits neither the factor model nor the sparse linear model frameworks. Nevertheless, the covariance matrix of ${\mathbf{x}}_t$ has a special structure. Observe that

$$\begin{array}{c}\hfill \mathbb{E}\left({\mathbf{x}}_t{\mathbf{x}}_t^{\top }\right)={\mathbf{I}}_p+0.64\mathbf{B}{\mathbf{B}}^{\top }.\end{array}$$

That is, the covariance matrix of ${\mathbf{x}}_t$ has the form of a spiked matrix.

Table 2 reports the RMSFEs of the competing methods considered. The sdPCA fares the best as it yields the smallest RMSFE, especially when s is small. This reflects that the sdPCA can better utilize the spiked covariance structure by selecting contributions from relevant predictors. However, the performance of GsP* and sdPCA converge when $s=0.75n$, and both of them compare favorably against the other alternatives. As expected, for non-sparse models, Lasso and random forest do not work well.

DGP 3. 

In this example, the predictors follow a VAR model:

$$\begin{array}{c}\hfill {\mathbf{x}}_t=\mathbf{B}{\mathbf{x}}_{t-1}+{\mathbf{\delta }}_t,\end{array}$$

where $\left\{{\mathbf{\delta }}_t\right\}$ are independent p-dimensional standard Gaussian processes. Let $\tilde{\mathbf{B}}$ be a randomly generated $p\times p$ matrix of rank $r_{\mathrm{DGP}}$. Then the AR coefficient matrix is constructed as

$$\begin{array}{c}\hfill \mathbf{B}=\frac{\tilde{\mathbf{B}}}{1.05\parallel \tilde{\mathbf{B}}\parallel },\end{array}$$

where $\parallel ·\parallel$ denotes the operator norm. The target variable is generated via

$$\begin{array}{c}\hfill y_t=0.5y_{t-1}+{\mathit{\beta }}^{\top }{\mathbf{x}}_{t-1}+{ϵ}_t,\end{array}$$

where $\left\{{ϵ}_t\right\}$ are independent standard Gaussian variates and $\mathit{\beta }={({\beta }_1,\dots ,{\beta }_p)}^{\top }$ with

$$\begin{array}{c}\hfill {\beta }_j=\left\{\begin{array}{cc}{(-1)}^j{u}_j,& 1\le j\le s\\ 0,& \mathrm{otherwise}\end{array}\right.,\end{array}$$

in which $u_j{\sim }_{\mathrm{i}.\mathrm{i}.\mathrm{d}.}\mathrm{U}(0.1,3.0)$.

In this DGP, the predictors, following a high-dimensional VAR(1) model, exhibit complex dynamics and correlations. Without simplifying structures such as latent factors, spiked covariance matrices, and sparsity, forecasting becomes difficult. Table 3 reports the RMSFEs of various competing methods. With small s, which corresponds to sparse models, Lasso can suitably choose the predictors and yield relatively accurate forecasts. However, as s increases, its RMSFE quickly increases and becomes similar to those of the factor-based approaches. The factor-based approaches as well as the RF performed similarly under this DGP. This example also demonstrates that the sdPCA method may encounter loss in prediction accuracy if the number of selected common factors is under-specified; see the case of $r_{\mathrm{GDP}}=15$ and $r=10$.

Overall, our simulation studies show that no forecasting method always dominates the competing methods used in the study, but the proposed GsP* procedure can be effective in some cases.

6. Empirical Examples

In this section, we apply the proposed method to two real data sets. The first data set is the U.S. macroeconomic data and the other consists of Taiwan particulate matter (PM_2.5) measurements. Forecasting plays an essential role in the applications pertaining to these two data sets, despite that both of them contain high-dimensional predictors. Moreover, these two datasets have different characteristics, which enable us to examine the forecasting performance of various forecasting methods available in the literature. In addition to the proposed GO-sdPCA approach, the competing methods used in the simulation studies, such as sdPCA, SW, LYB, Lasso, and RF, are also employed in this section as benchmarks.

For both datasets, we consider the rolling-window h-step-ahead forecasting. Let $\left\{y_t\right\}$ be the target variable of interest. For predicting $y_{t+h}$, the factor-based methods use the predictive equation

$$\begin{array}{c}\hfill {\widehat{y}}_{t+h}=\sum _{k=1}^q{\alpha }_k{y}_{t-q+1}+{\mathit{\beta }}^{\top }{\widehat{\mathbf{f}}}_t+{ϵ}_{t+h},\end{array}$$

where ${\widehat{\mathbf{f}}}_t$ is the vector of estimated common factors by different methods. For Lasso and RF, $y_t,y_{t-1}\dots ,y_{t-q+1},{\mathbf{x}}_t,\dots ,{\mathbf{x}}_{t-q+1}$ are used as potential predictors.

6.1. U.S. Macroeconomic Data

First, we consider the U.S. monthly macroeconomic data from January 1973 to June 2019. The data are from the FRED-MD database maintained by St. Louis Federal Reserve at https://research.stlouisfed.org/econ/mccracken/fred-databases/. We transform the time series to stationarity according to McCracken and Ng [9], and, after discarding some variables containing missing values, there are 125 macroeconomic time series. Among them, we focus on predicting (1) industrial production index, (2) unemployment rate, (3) CPI-All, and (4) Real manufacturing and trade industries sales (M&T sales). Time plots of these four target series after transformations are depicted in Figure 1.

For this data set, we consider $h=1$ and different combinations of q and r (the number of factors extracted). Since q lags of the predictors are used, there are $125q$ total predictors, which exceeds the sample size in each window. The last twenty years of data (240 time periods) are reserved for testing. In addition to the root mean squared forecasting errors, we also report the p-values of the Diebold-Mariano test [26,27] against the alternative hypothesis that the proposed GsP* procedure is more accurate.

Table 4, Table 5, Table 6 and Table 7 present the results. The proposed GO-sdPCA achieved the lowest RMSFE for three of the targeted series: industrial production index, unemployment rate, and CPI. It notably outperformed the factor-based alternatives sdPCA, SW, and LYB with high confidence in forecasting these three series. In addition, no other methods demonstrated such consistent effectiveness across the targeted series.

6.2. Particulate Matters in Taiwan

We next consider the data of hourly PM_2.5 measurements in Taiwan during March of 2017. The data are sourced from the SLBDD package [28] in R. Each series in the data set represents measurements (in micrograms per cubic meters, $\mu \mathrm{g}/{\mathrm{m}}^3$) taken by a novel device known as the AirBox. After an initial examination of these series, we remove series 29 and 70 because these series are near identically zero except at a few time points, offering no useful variations. Among the 516 series in the data set, we choose series 101, 201, 301, 401, as target series. Figure 2 depicts their time plots.

We consider both one-step-ahead ($h=1$) and two-step-ahead ($h=2$) forecasts, and employ $q\in \{2,3\}$ lags and $r\in \{2,4,6\}$ factors in forecasting. The last ten days of data (240 time periods) are reserved for out-of-sample testing. Table 8, Table 9, Table 10 and Table 11 present the results for $h=1$. The proposed GO-sdPCA ranks as the most predictive method in terms of RMSFE for two of the four targeted series: Series 101 and 201. It outperforms the Lasso in forecasting all four targeted series. The performance of LYB is also noteworthy. It ranks among the best two methods for all four targeted series. The DM test, on the other hand, indicates the difference in forecasting accuracy is not statistically significant. Contrary to the case in the macroeconomic data, the proposed GsP* only outperformed the factor-based methods, including the sdPCA, with some weak confidence. For the two-step-ahead forecasts, for which the results are reported in Table 12, Table 13, Table 14 and Table 15, the performance of various methods is more entangled, with the RF consistently ranked among the top two methods for three of the four targeted series. With some weaker confidence, GsP* is the best performing factor-based method for two of the four targeted series. Again, the DM test fails to separate significantly various forecasting methods. We believe this result might be caused by the substantial uncertainty in the PM_2.5 measurements.

In sum, for $h=1$, GsP* is effective in forecasting PM_2.5 data as well as the U.S. macroeconomic data. This implies the proposed procedure is able to capture highly predictive factors across various applications. On the other hand, for $h=2$, the dynamic nature of the data may be much more involved. The performance of the proposed method becomes similar to those of the other forecasting methods.

7. Discussion and Concluding Remarks

In this work, we proposed a novel method for time series forecasting when many predictors are available. The rationale behind our method is to mine the possibly many (compared to the sample size) factor-relevant predictors while reducing the effect of the many irrelevant variables in the high-dimensional data. The results in the simulation studies and the empirical applications suggest that the proposed method is useful for improving upon both the factor-based methods and methods for sparse linear models such as the Lasso. Finally, we remark that the theoretical investigation of the peeling technique, a key ingredient in our method, would be an interesting topic for future research.