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.
94,315 characters · 13 sections · 100 citation commands
Divide-and-Conquer: A Distributed Hierarchical Factor Approach to Modeling Large-Scale Time Series Data
{\sl Keywords}: Distributed Factor Model, Large-scale Time Series, Principal Component Analysis, Commonality, Heterogeneity
Multivariate time series with low-dimensional factor structures have been widely used in many scientific fields, including economics, finance, and statistics, among others. Asset returns in finance are often modeled as functions of a small number of factors, see ross1976, StockWatson_1989, famafrench1993,famafrench2015, and StockWatson_1998. Macroeconomic variables of multiple countries are often found to have comovements, see GregoryHead1999 and forni2000. For large-scale data, factor models not only provide an effective way to extract useful common information, but also reveal their low-dimensional structure when the dimension is high. Indeed, various factor models have been developed in the literature for multivariate and, recently, high-dimensional time series. See the approximate factor models in chamberlain1983, StockWatson_1989,StockWatson_2002a,StockWatson_2002b, BaiNg_Econometrica_2002, and Bai_Econometrica_2003; the dynamic factor models in forni2000,forni2005; the factor models with dynamically dependent factors in lamyao2012 and gaotsay2018a,gaotsay2020a, gaotsay2020b, among others. Other dimension reduction techniques for time series with a low dimensional structure are also available, see BoxTiao_1977, TiaoTsay_1989, Tsay_2014, and the references therein. However, all available methods require the accessibility to the entire data set of interest by a single computer, and such a requirement has become increasingly infeasible with the rapid developments in information science and technology.
A fundamental challenge in modeling large-scale, high-dimensional time series is the storage of big data with huge sample size along the time horizon and a large number of series from the spatial domain. In many applications, we need to analyze internet data containing billions, or even trillions, of data points across spatial dimension and time horizon, making even a linear pass of the whole data set infeasible. In addition, such data are often highly dynamically dependent and may not be stored in a centralized machine. For example, stock returns of different exchanges and different countries may be stored in different databases, macroeconomic data are often collected and stored in different government agencies, and the high-frequency spatio-temporal data, such as the real-time PM$_{2.5}$ indexes, are usually monitored and collected by local agencies across a country. Consequently, it has become increasingly difficult to analyze them due to the communication cost, privacy concern, and data security, among others. Therefore, a distributed statistical method and computing architecture is needed. The fundamental approach to store and process such data is to divide and conquer with the basic idea being to partition a large problem into many tractable sub-problems. The sub-problems are tackled in parallel by different processing units. Intermediate results from individual machines are then combined to be further processed on the central server to yield the final output. Many distributed methods following this approach have been developed during the past decades. See the distributed regression methods in zhang2013 and chenxie2014, the distributed PCA in qu2002 and fan2019 with a horizontal partition regime, where each server contains a small portion of the data of all subjects, and another distributed PCA in kargupta2001, lietal2011, and BertrandMoonen2014 with a vertical partition regime, where each server contains all the data of a small subset of the subjects. However, most of the methods mentioned above were developed for independent and identically distributed data, and they cannot tackle high-dimensional time series directly. In addition, the statistical guarantees are absent in some of the aforementioned literature.
The goal of this paper is to propose a hierarchical factor framework for modeling large-scale time series data. When data are stored across multiple machines, or data centers, the traditional factor modeling requires analyzing the whole data set in a central server, resulting in excess computational cost caused by data transfer and communication. Our proposed framework can overcome this difficulty efficiently. To illustrate, we sketch a 3-layer hierarchical factor modeling framework in Figure (ref), where the first layer has $n_1+\cdots+n_m$ processing units with $n_i$ of them connected to the $i$-th node in the second layer, and the $m$ nodes of the second layer are connected to the central server, which forms the 3rd layer. Suppose the marketplaces submit orders to the central server and expect to receive its responses later. Here the relevant data are either stored in the first-layer units, or assigned to them by the central server based on the capacities of the units and/or the problem under study. Each 1st-layer unit performs a PCA to extract common factors from the data assigned to it and sends those common factors to its connected node in the 2nd layer. After collecting the common factors sent to it, each node of the 2nd layer performs, again, a PCA to obtain the 2nd-level common factors. It then transfers those extracted common factors to the central server, and the central server performs yet another PCA to extract the global common factors using the 2nd-level common factors it received. Clearly, the idea of hierarchical factor analysis can easily be extended to more layers. In this way, the proposed scheme is similar to that of deep neural networks (DNN) with multiple layers, where each computing unit or node can be treated as a neuron. However, our objective is different from that of DNN, because here each neuron only computes linear combinations of the factors it receives from its direct subordinates via PCA.
We focus on the scenario of 2 levels in our methodological formulations and theoretical derivations because it serves as a building-block for the general case, and the idea can easily be generalized to any finite number of hierarchies. Specifically, we consider the $m$ computing nodes and the central server of Figure (ref), in which each node contains all observations of a small subset of the $p$-dimensional time series. For a $p$-dimensional time series and $m$ servers, we first apply PCA to the $p_i$-dimensional time series stored on the $i$-th node to obtain some common factors for the $i$-th group (or cluster), which are called {\em within group factors}. The central server then collects all within group factors to form a new vector time series, and performs another PCA to select higher-level common factors. The resulting common factors serve as the global factors (or {\em between group factors}) that capture the cross-sectional information across all the groups, and the idiosyncratic terms of the 2nd-level PCA would be the {\em group-specific factors} of the corresponding 1st-level group. This approach involves a two-fold dimension reduction procedure using PCA, and is in line with the distributed computing with a vertical partition regime. Thus, it is different from the one in fan2019, and it is also different from the one in kargupta2001, which only focuses on a single-stage PCA to reduce the dimension of each group without providing rigorous statistical guarantees. Asymptotic properties of the proposed method are established for diverging group dimension $p_i$ and the sample size $T$. We use both simulated data and real examples to assess the performance of the proposed method in finite samples. In the presence of group structure and assuming that the data set can be analyzed by a single machine, we apply the ratio-based method in lamyao2012 and ahn2013 and the BIC-type criterion in BaiNg_Econometrica_2002 to the whole data set to select the number of common factors. We find that the former tends to find the number of global factors, whereas the latter selects the total number of group specific factors, as the number of common factors of the whole data. We also compare our method with the commonly used ones in the literature concerning the forecastability of the extracted factors, and find that our proposed new diffusion-index model outperforms the traditional one in forecasting eight commonly employed U.S. economic time series.
From a statistical perspective, it is worth mentioning that our resulting model is similar to the dynamic hierarchical factor models in Moench2013, the dynamic factor model in alonso2020 with cluster structure, and the panel data model in andobai2017 with predictors and factor structures. However, their methods can only handle stationary time series, and require storing the entire data on a single machine, resulting in high computational cost for large data sets. In addition, the statistical guarantees remain absent in Moench2013 and alonso2020. Our proposed method is different from those of the aforementioned literature, and its statistical properties are established for both stationary and unit-root nonstatinary time series data under some reasonable assumptions. Furthermore, when the data are in small-scale and stored in a single machine, our method provides an alternative way to modeling panel time series data with global and group-specific factor structures if a clustering method is available. We discuss this issue in Section 2.3 and propose a variation-based clustering method to distribute the data to different nodes if the group memberships are not available. This method for distributed computing is of independent interest by itself.
The contributions of this paper are multi-fold. First, the proposed distributed framework accommodates the large-scale features of modern big data and can also be used to analyze independent and identically distributed data with parallel computing. The procedure is easy to implement and addresses two important features of the data: Heterogeneity and Commonality, where the heterogeneity is preserved in the group specific factors of the first PCA, and the commonality of different groups is accounted for by the global factors obtained in the subsequent analysis. Second, unlike most of the distributed and hierarchical modeling procedures available in the literature, we established statistical guarantees of the proposed method provided that the number of machines used is not extremely large. Third, the framework is extended to unit-root nonstationary time series to extract the common stochastic trends of the data. We also briefly discuss ways to incorporate horizontal partition regime into the proposed vertical regime for analyzing extremely large-scale data. Finally, we propose a new diffusion-index model for time series forecasting. This new model makes use of the global common factors and the group specific factors as predictors, and it is able to embrace the empirical findings in BoivinNg_joe_2006 that factors extracted from a smaller number of series might yield satisfactory or even better forecasting results.
The rest of the paper is organized as follows. We introduce the proposed model and estimation methodology in Section (ref) and study the theoretical properties of the proposed model and its associated estimates in Section (ref). The proposed method is extended to unit-root nonstationary time series data in Section (ref). Section (ref) briefly discusses ways to combine the horizontal and vertical partitions together when the sample size $T$ is extremely large. Numerical studies with both simulated and real data sets are given in Section (ref), and Section (ref) provides some concluding remarks. All technical proofs are given in an Appendix. Throughout the article, we use the following notation. For a $p\times 1$ vector ${\mathbf u}=(u_1,..., u_p)',$ $||{\mathbf u}||_2 =\|{\mathbf u}'\|_2= (\sum_{i=1}^{p} u_i^2)^{1/2} $ is the Euclidean norm, $\|{\mathbf u}\|_\infty=\max_{1\leq i\leq p}|u_i|$ is the $\ell_\infty$-norm, and ${\mathbf I}_p$ denotes a $p\times p$ identity matrix. For a matrix ${\mathbf H}=(h_{ij})$, $\|{\mathbf H}\|_1=\max_j\sum_i|h_{ij}|$, $|{\mathbf H}|_\infty=\max_{i,j}|h_{ij}|$, $\|{\mathbf H} \|_F=\sqrt{\sum_{i,j}h_{ij}^2}$ is the Frobenius norm, $\|{\mathbf H} \|_2=\sqrt{\lambda_{\max} ({\mathbf H}' {\mathbf H} ) }$ is the operator norm, where $\lambda_{\max} (\cdot) $ denotes for the largest eigenvalue of a matrix, and $\|{\mathbf H}\|_{\min}$ is the square root of the minimum non-zero eigenvalue of ${\mathbf H}'{\mathbf H}$. The superscript ${'}$ denotes the transpose of a vector or matrix. We also use the notation $a\asymp b$ to denote $a=O(b)$ and $b=O(a)$.
Consider a 2-level framework. Let ${\mathbf y}_t= ({\mathbf y}_{1,t}',{\mathbf y}_{2,t}'...,{\mathbf y}_{m, t}')'$ be an observable $p$-dimensional time series consisting of $m$ groups with ${\mathbf y}_{i,t}=(y_{i1,t},...,y_{i p_i, t})'\in \mathbb{R}^{p_i}$, for $1\leq i\leq m$, stored on $m$ machines or servers, and $p_1+p_2+...+p_m=p$. We consider the case of large $p$, but slowly diverging $m$. If there are more than two layers, $\{{\mathbf y}_{i,t},t=1,...,T\}$ may consist of stacked common factors from its associated groups in the previous layer. For simplicity, we assume that $E({\mathbf y}_{i,t})={\bf 0}$ and $\{{\mathbf y}_{i,t}| t=1,...,T\}$ admits a latent factor structure:
where ${\mathbf g}_{i,t}=(g_{i1,t},...,g_{ir_i,t})'$ is a $r_i$-dimensional latent factor process for the $i$-th group ${\mathbf y}_{i,t}$, ${\mathbf A}_i$ is the associated factor loading matrix, and ${\mathbf e}_{it}$ is the idiosyncratic term of the $i$-th group, which is independent of the group common factors ${\mathbf g}_{it}$. Letting ${\mathbf g}_t:=({\mathbf g}_{1,t}',...,{\mathbf g}_{m,t}')'$ be the collection of all group common factors, which is a $k_m$-dimensional latent vector with $k_m=r_1+...+r_m$, we assume that
where ${\mathbf f}_t=(f_{1,t},...,f_{r,t})'$ is an $r$-dimensional global factor process capturing the cross-sectional dependence across the $m$ groups via their group factors, ${\mathbf B}\in \mathbb{R}^{k_m\times r}$ is the associated loading matrix, and ${\mathbf u}_t$ = $({\mathbf u}_{1t}',\ldots,{\mathbf u}_{mt}')'$ is the idiosyncratic term of the stacked process ${\mathbf g}_t$ with ${\mathbf u}_{it}$ constituting the {\em group specific factors} of the $i$-th group, which is orthogonal to the global factors ${\mathbf f}_t$. Putting Equations ((ref)) and ((ref)) together, we obtain the following matrix form of the proposed 2-level model:
where ${\mathbf A}=\mbox{diag}({\mathbf A}_1,...,{\mathbf A}_m)$ is a $p\times k_m$ block-diagonal matrix, ${\mathbf f}_t$ is the global common factor process, and ${\mathbf u}_t$ consists of $m$ vectors of group specific factor processes. For the identification issue, we assume ${\mathbf f}_t$, ${\mathbf u}_t$, and ${\mathbf e}_t$ are mutually independent of each other. The decomposition in ((ref)) indicates that the global common factors come from the cross-sectional dependence between the group factors, and the ${\mathbf u}_t$ process, which captures less cross-sectional dependence, would be the specific factors uniquely representing the individual groups.
We mention that andobai2017 and alonso2020 studied panel data with global and group-specific factors. However, our factors and the associated loading matrices are different from those of the aforementioned two papers. Most importantly, the methods used in the prior papers require the entire data set to estimate the global and group-specific factors, even when the group membership is known.
It is well known that the factors and their associated loading matrices in Models ((ref)) and ((ref)) are not uniquely identified. We follow the restrictions used in BaiNg_Econometrica_2002 and fan2013 and assume that
Our goal here is to estimate ${\mathbf A}_i$ and ${\mathbf B}$, and to extract the factor processes ${\mathbf g}_{it}$ and ${\mathbf f}_t$. We also study asymptotic properties of the proposed estimates.
In this section, we discuss an estimation method based on a given realization $\{{\mathbf y}_t,|t=1,...,T\}$. We start with the case of knowing $r_i$ and $r$, where $r_i$ is the number of common factors of the $i$-th group and $r$ is the number of global common factors. The specifications of $r_i$ and $r$ are discussed shortly after. For each group (or machine), we use the least squares (LS) method to estimate ${\mathbf A}_i$ and ${\mathbf g}_{i,t}$ as follows:
subject to the normalization conditions:
Similarly to BaiNg_Econometrica_2002 and Bai_Econometrica_2003, let ${\mathbf Y}_i=({\mathbf y}_{i,1},...,{\mathbf y}_{i,T})'$ and ${\mathbf G}_i=({\mathbf g}_{i,t},...,{\mathbf g}_{i,T})'$. The Principal Components (PC) estimation is equivalent to the LS method in ((ref)) that minimizes $\mbox{tr}[({\mathbf Y}_i-{\mathbf G}_i{\mathbf A}_i')'({\mathbf Y}_i-{\mathbf G}_i{\mathbf A}_i')]$ subject to the normalizations ${\mathbf G}_i'{\mathbf G}_i/T={\mathbf I}_{r_i}$ and ${\mathbf A}_i'{\mathbf A}_i$ being diagonal. Consider the eigenvalue-eigenvector analysis of the $T\times T$ matrix ${\mathbf Y}_i{\mathbf Y}_i'$. The estimator for ${\mathbf G}_i$, denoted by $\widehat{\mathbf G}_i=(\widehat{\mathbf g}_{i,1},...,\widehat{\mathbf g}_{i,T})'$, is $\sqrt{T}$ times the eigenvectors corresponding to the $r_i$ largest eigenvalues, in decreasing order. Then, $\widehat{\mathbf A}_i={\mathbf Y}_i'\widehat{\mathbf G}_i/T$ is the corresponding factor loading matrix. The common component matrix ${\mathbf G}_i{\mathbf A}_i'$ is estimated by $\widehat{\mathbf G}_i\widehat{\mathbf A}_i'$.
Turn to the estimations of ${\mathbf f}_t$ and ${\mathbf B}$. Let $\widehat{\mathbf G}=(\widehat{\mathbf G}_1,...,\widehat{\mathbf G}_m)$ and ${\mathbf F}=({\mathbf f}_1,...,{\mathbf f}_T)'$. We apply, again, the PC estimation method so that
subject to the constraints ${\mathbf F}'{\mathbf F}/T={\mathbf I}_r$ and ${\mathbf B}'{\mathbf B}$ being diagonal. Therefore, $\widehat{\mathbf F}$ consists of the eigenvectors (multiplied by $\sqrt{T}$) associated with the $r$ largest eigenvalues, in decreasing order, of the $T\times T$ matrix $\widehat{\mathbf G}\widehat{\mathbf G}'$. Then, $\widehat{\mathbf B}=\widehat{\mathbf G}'\widehat{\mathbf F}/T$ and, hence, the group specific factors $\widehat{\mathbf U}=(\widehat{\mathbf u}_1,...,\widehat{\mathbf u}_T)'=\widehat{\mathbf G}-\widehat{\mathbf F}\widehat{\mathbf B}'$, where ${\mathbf u}_t=({\mathbf u}_{1,t}',...,{\mathbf u}_{m,t}')'$ with ${\mathbf u}_{i,t}$ being the group-specific factor process associated with the $i$-th group.
Next, we discuss the determination of the number of factors in each group and the global one. For the estimation of $r_i$, we apply the well-known method in BaiNg_Econometrica_2002. It estimates $r_i$ by
where $K$ is a prescribed upper bound, $\widehat{\mathbf G}_{i,K_i}$ is a $T\times K_{i}$ matrix whose rows are $\sqrt{T}$ times the eigenvectors corresponding to the $K_i$ largest eigenvalues of the $T\times T$ matrix ${\mathbf Y}_i{\mathbf Y}_i'$, and $g(T,p_i)$ is a penalty function of $(p_i,T)$ such that $g(T,p_i)=o(1)$ and $\min\{p_i,T\}g(T,p)\rightarrow\infty$. Two examples of $g(T,p_i)$ suggested by BaiNg_Econometrica_2002 are IC1 and IC2, respectively, \[g(T,p_i)=\frac{p_i+T}{p_iT}\log\left(\frac{p_iT}{p_i+T}\right)\quad \text{and}\quad g(T,p_i)=\frac{p_i+T}{p_iT}\log(\min\{p_i,T\}).\] For the estimation of $r$, besides the information criterion as that in ((ref)), we also adopt the ratio-based method in lamyao2012 and ahn2013 with some modifications. Let $\widehat\lambda_1\geq ...\geq \widehat\lambda_{k_m}$ be the $k_m$ eigenvalues of $\widehat{\mathbf G}'\widehat{\mathbf G}$, we estimate $r$ by
where we take $R=[k_m/2]$ as suggested by lamyao2012. Alternatively, we may estimate the number of global factors by selecting a threshold $0<c_0<1$ such that
where the ratio denotes the fraction of the variances explained by the first $j$ global factors. This estimate is easy to implement and could be practically useful for a prescribed $c_0$ (e.g., $c_0=85\%$).
In this section, we consider the case in which the data need to be distributed to different machines by the central server according to certain characteristics of the observed time series, or the data are in small scale and stored on a single machine, but require clustering. In other words, we consider the case that the group memberships are not available. To apply the proposed method, we need to infer the group memberships of the series under study, which can be done via a clustering analysis.
Let ${\mathbf y}_t=(y_{1t},y_{2t},...,y_{pt})'$ be a $p$-dimensional time series of interest. There are many time-series clustering methods available in the literature, especially when the data can be stored in a single machine, and we briefly mention some of the commonly used ones. The first method is based on some selected features. Specifically, let ${\mathbf S}=\{{\mathbf s}_1,...,{\mathbf s}_p\}$, where ${\mathbf s}_i$ is a vector of selected features of the time series $\{y_{it},t=1,...,T\}$. For example, ${\mathbf s}_i$ may contain the first few lags of the sample autocorrealtions or the sample partial autocorrelations of $y_{it}$. Clustering can then be based on a chosen distance measure between ${\mathbf s}_i$ and ${\mathbf s}_j$. See, for instance, Caiado2006 and the references therein. The second method is based on a distance measure between the two series $y_{it}$ and $y_{jt}$ directly. For time series $\{y_{it}\}_{t=1}^T$ and $\{y_{jt}\}_{t=1}^T$, one can choose a distance measure $d_{ij}$, which can be the Euclidean or the Manhattan distance, to define the similarities, and the smaller the distance the closer the two series. The $k$-means clustering method can then be used. $d_{ij}$ can also be the cross-sectional correlations between the two series, and a large correlation implies that the two series can be grouped together. See, for instance, alonso2019 and the references therein.
When the dimension $p$ is sufficiently large so that the data cannot be stored in a single node, the above two classes of clustering methods have advantages and disadvantages. For clustering methods based on selected features: (a) There is no need to store all data in a single machine, but one needs to assume that the selected feature vectors ${\mathbf S}$ can be stored in a single machine; (b) One may also need to transfer data between many computing units based on the clustering results obtained by the central server, which incurs communication costs. For clustering methods based on similarity measures: (a) One needs to compute the distance between all pairs of time series, which also requires data transfer in and out of a single computing unit; (b) One also needs to transfer data to other machines based on the clustering results obtained by the central server. The communication costs can be extremely high due to the computation of the distances between all possible pairs of the time series.
In view of the above discussion, we leverage the advantages of feature-based procedures to propose a variation clustering procedure, which not only can distribute the series to different machines efficiently, but also avoids requiring to transfer the data to and out of the center server. The procedure is also motivated by the fact that PCA is not scale invariant and often allocates heavier weights to series with larger sample variances. For example, we set the seed to 1234 in R and simulated a 5-dimensional time series of a VAR(1) model with $500$ observations, where the coefficient matrix $\boldsymbol{\Phi}$ is generated by drawing elements of a $5\times 5$ matrix ${\mathbf M}$ randomly from $U(-2,2)$, then properly normalized to ensure stationarity, and the noises are independently generated from the standard multivariate normal distribution. Specifically, the VAR(1) coefficient matrix used is $\boldsymbol{\Phi} = 0.9{\mathbf M}/\|{\mathbf M}\|_2$. The variances of the five series are $(0.1484,0.1482,0.1574,0.1366,0.1582)$ so that they have similar variances. We apply PCA to these five time series and obtain the principal vectors below
where the $i$-th column corresponds to the $i$-th principal direction in decreasing order of the captured variances. We now multiply the first series by 10 so that its variance becomes $14.84$. We apply PCA, again, and obtain the following principal vectors:
Here the dynamical dependence between the five time series remains the same, but PCA assigns rather different weights to the first series in estimating the first principal component, as shown by the first column in ${\mathbf P}_5$ and ${\mathbf P}_5^*$. This phenomenon may lead to complications in selecting and interpreting the resulting PCs and, hence, the common factors. To mitigate the scaling effect in applying PCA, let $\widetilde{\mathbf y}_t$ be the original $p$-dimensional series with unknown group memberships and $s_i^2$ be the sample variance of the $i$-th component, we require that $s_i^2$ be transferred to the central server, if needed. The central server then sorts the sample variances in the decreasing order as $s_{(1)}^2\geq s_{(2)}^2\geq...\geq s_{(p)}^2$. Denote the re-arranged series by ${\mathbf y}_t$ and we divide the $p$ series into $m$ groups of roughly equal sizes. That is, we set $p_i= [ p/m]$ for $i=1,...,m-1$, and $p_m=p-(m-1)[ p/m]$, where $[p/m]$ is the integer part of $p/m$. The first $p_1$ series of ${\mathbf y}_t$ are assigned to the first node, the next $p_2$ series to the 2nd node and so on until the last $p_m$ series are assigned to the $m$-th node. In this way, the observed time series are assigned to groups based on their sample variances. There are several advantages of this variation-based clustering method. First, it can mitigate the scaling effect affecting the performance of PCA. Second, it does not require transferring the entire time series data to the center server and it is easy to implement, even when both the dimension $p$ and sample size $T$ are reasonably large.
There are other efficient clustering methods available in the literature to assign individual time series to the $m$ nodes. For instance, some prior (or domain) knowledge may be available for clustering. The stock returns can naturally be clustered using their Standard Industrial Classification Codes, macroeconomic time series using their industrial characteristics, and spatio-temporal time series using their geographical locations. An important issue to consider is the cost of computing and data transfer in the assignment.
We introduce next a new variant of the Diffusion Index (DI) model to predict a selected series $x_t$ based on all the components of ${\mathbf y}_t$. In some applications, $x_t$ may be associated with a specific group of ${\mathbf y}_t$. For example, in macroeconomic forecasting, $x_t$ can be classified into a group of ${\mathbf y}_t$ according to its industrial sector; in spatio-temporal data, $x_t$ often belongs to a cluster of ${\mathbf y}_t$ according to its geographic features. Assume that $x_t$ belongs to the $i$-th group. Then, with the estimated global factors $\widehat{\mathbf f}_t$ and the group-specific ones $\widehat{\mathbf u}_{i,t}$, the $h$-step ahead prediction of $x_t$ at the forecast origin $t$ is
where $\widehat\boldsymbol{\beta}_h(L)=\widehat\boldsymbol{\beta}_0'+\widehat\boldsymbol{\beta}_1'L+...+\widehat\boldsymbol{\beta}_{q_1}'L^{q_1}$, $\widehat\boldsymbol{\gamma}_h(L)=\widehat\boldsymbol{\gamma}_0'+\widehat\boldsymbol{\gamma}_1'L+...+\widehat\boldsymbol{\gamma}_{q_2}'L^{q_2}$, and $\widehat\phi_h(L)=\widehat\phi_0+\widehat\phi_1L+...+\widehat\phi_{q_3}L^{q_3}$ are matrix or scalar polynomials of order $q_1, q_2$ and $q_3$, respectively, with $L$ being the lag operator, $\widehat\alpha_h, \widehat\boldsymbol{\beta}_i,\widehat\boldsymbol{\gamma}_j$, and $\widehat \phi_k$ are parameter estimates based on the available data $\{x_1,...,x_t\}$, $\{\widehat{\mathbf f}_1,...,\widehat{\mathbf f}_t\}$ and $\{\widehat{\mathbf u}_{i,1},...,\widehat{\mathbf u}_{i,t}\}$ using, for example, the ordinary LS method. On the other hand, if we are not certain about which group-specific factor process to use in ((ref)), we may include all the group-specific factors and using group penalized method to estimate the coefficients such as the group Lasso in yuanlin2006. Here the cross-validation method for choosing the tuning parameters should be replaced by a rolling-window method for time series data due to its dynamic structure. We do not pursue further any detail here.
We begin with some assumptions for Models ((ref)) to ((ref)), where only $\{{\mathbf y}_{it}|1\leq i\leq m,1\leq t\leq T\}$ are observable. Similar ones are also given in BaiNg_Econometrica_2002 and fan2013. All proofs of the theorems are given in the Appendix.
Assumption (ref) implies that the factors used are pervasive as those in BaiNg_Econometrica_2002 and fan2013. Assumption (ref) is standard for dependent random processes, and it implies the temporal dependence conditions in BaiNg_Econometrica_2002 and Bai_Econometrica_2003. See gaoetal2017 for a theoretical justification of the mixing conditions for VAR models. Assumption (ref)(i) guarantees that the covariances of the idiosyncratic terms are well defined, which will be used to show the consistency of the ratio-based method in ((ref)), and (ii) implies that the factors and the idiosyncratic terms all have exponential tails, which is stronger than the moment conditions in BaiNg_Econometrica_2002. Assumption (ref) controls the cross-sectional dependence of the idiosyncratic terms. Under Assumptions (ref)--(ref) and (ref)(i)-(iii), all the assumptions used in Bai and Ng (2002) are satisfied. See also fan2013. Assumptions (ref)(iv)-(v) are imposed to control the cross-sectional dependence of the idiosyncratic terms as that of Assumption E(2) and Assumption F(1) in Bai_Econometrica_2003, and the proof of Theorem 1 in BaiNg_joe_2006 is valid under these assumptions.
Some remarks of Theorem (ref) are in order.
Partition ${\mathbf B}=({\mathbf B}_1',...,{\mathbf B}_m')'$ and let ${\mathbf C}_i={\mathbf A}_i{\mathbf B}_i$. Denote the $l$-th loading of ${\mathbf A}_i$ and ${\mathbf C}_i$ by ${\mathbf a}_{i,l}$ and ${\mathbf c}_{i,l}$, respectively. We have the following theorem on the consistency of the estimated factor loading matrices.
Note that ${\mathbf A}_i$ and ${\mathbf B}$ cannot be uniquely determined due to the identification issue. But the linear spaces spanned by their columns, denoted by $\mathcal{M}({\mathbf A}_i)$ and $\mathcal{M}({\mathbf B})$, respectively, can be uniquely estimated. See gaotsay2020a,gaotsay2020b for details. To this end, we adopt the discrepancy measure used by panyao2008, which is a well-defined distance measure on a quotient matrix space as shown therein: for two $p\times d$ semi-orthogonal matrices ${\bf Q}_1$ and ${\bf Q}_2$ satisfying the condition ${\bf Q}_1'{\bf Q}_1={\bf Q}_2'{\bf Q}_2={\mathbf I}_{d}$, the difference between the two linear spaces $\mathcal{M}({\bf Q}_1)$ and $\mathcal{M}({\bf Q}_2)$ is measured by
Note that $D({\bf Q}_1,{\bf Q}_2) \in [0,1].$ It is equal to $0$ if and only if $\mathcal{M}({\bf Q}_1)=\mathcal{M}({\bf Q}_2)$, and to $1$ if and only if $\mathcal{M}({\bf Q}_1)\perp \mathcal{M}({\bf Q}_2)$. In addition, if we denote the singular values of ${\mathbf Q}_1'{\mathbf Q}_2$ by $\{\sigma_i\}_{i=1}^d$ in descending order, then the principal angles between $\mathcal{M}({\mathbf Q}_1)$ and $\mathcal{M}({\mathbf Q}_2)$, denoted by $\boldsymbol{\Theta}({\mathbf Q}_1,{\mathbf Q}_2)=\mbox{diag}(\theta_1,...,\theta_d)$, are defined as $\mbox{diag}\{\cos^{-1}(\sigma_1),...,\cos^{-1}(\sigma_d)\}$; see, for example, Theorem I.5.5 of stewart-sun1990. The squared Frobenius norm of the so-called $\sin\boldsymbol{\Theta}$ distance is
which implies that the above two measures are equivalent for finite $d$. Therefore, we shall use only the distance in ((ref)) to present our theoretical results. The following theorem establishes the consistency of the estimated loading matrix $\widehat{\mathbf A}_i$, the matrix $\widehat{\mathbf B}$, and the extracted common factors $\widehat{\mathbf A}_i\widehat{\mathbf g}_{i,t}$ and $\widehat{\mathbf A}_i\widehat{\mathbf B}_i\widehat{\mathbf f}_t$.
Theorem (ref) explores a different aspect of the consistency from that in Theorem (ref). The uniform rates in Theorem (ref) do not imply the convergence pattern as $T\rightarrow\infty$ for fixed but large $p_i$, while the average consistencies in Theorem (ref) provide the asymptotic behavior of the estimated loading matrices and the extracted factors, and they are indeed consistent as the sample size $T\rightarrow\infty$, $\underline{p} \rightarrow \infty$, $m=o(\sqrt{T})$, and the clusters are of similar sizes. The condition for $m$ is different from that in Theorem 1, but we can achieve the consistencies in ((ref)) and ((ref)) so long as $m=o(\sqrt{T},\underline{p})$. The uniform consistency in ((ref)) can still be achieved under this setting so long as $T^{1/4}/\sqrt{k_m}\rightarrow\infty$ if some $r_i$ diverges slowly. On the other hand, the conditions of $m$ in Theorem (ref) can be relaxed, because we can apply the Bernstein-type inequality in, for example, merlevede2011, to reduce the term $m$ in ((ref)) to an exponent of $\log(m)$ under sightly stronger conditions such that $\max_{1\leq i\leq m}D(\widehat{\mathbf A}_i,{\mathbf A}_i)=O_p(\sqrt{\log(m)/T})$, instead of $O_p(mT^{-1/2})$. We omit the details here to save space.
Next, we discuss the selections of the number of factors as they are unknown in practice.
By Theorem (ref), we can replace $r_i$ and $r$ in Theorems (ref) and (ref) by $\widehat r_i$ and $\widehat r$, respectively, and the results therein continue to hold. See gaotsay2020b for a similar argument. For the discrepancy measure defined in ((ref)) and those considered in Theorem (ref), we take $\widehat{\mathbf A}_i$ for example and modify ((ref)) as \[\widetilde D(\widehat{\mathbf A}_i,{\mathbf A}_i)=\sqrt{1-\frac{1}{\max(\widehat r_i,r_i)}\textrm{tr}({\widehat{\mathbf A}}_i{\widehat{\mathbf A}}_i'{{\mathbf A}}_i{{\mathbf A}}_i')},\] which takes into account the fact that the dimensions of $\widehat{\mathbf A}_i$ and ${\mathbf A}_i$ may be different. See gaotsay2020b for a similar argument.
In this section, we consider the case when the series ${\mathbf y}_{i,t}$ is an $I(1)$ process, for $1\leq i\leq m$, and, hence, extend the proposed method to unit-root nonstationary time series. To do so, we assume that ${\mathbf f}_t$ in ((ref)) is an $I(1)$ process that captures the global stochastic trends, and ${\mathbf e}_{i,t}$ and ${\mathbf u}_t$ are stationary components. Consequently, each ${\mathbf y}_{i,t}$ is an $I(1)$ process and ${\mathbf g}_{i,t}$ is the group common stochastic trend of the $i$-th group. This setup is natural and is similar to the nonstationary factor models of penaponcela2006, in which the dimension is fixed and ${\mathbf u}_{i,t}$ is the stationary common factor of the $i$-th group. However, we do not need the idiosyncratic terms ${\mathbf e}_{i,t}$ to be white noises. Under the assumption, the first PCA conducted in each group is to seek the group common stochastic trends, and the second one tries to find the global common stochastic trends and, hence, allows for the existence of co-integration relationships between the group common factors. On the other hand, as the dimension of ${\mathbf g}_t$ is relatively large, we consider the case that its cointegration rank is low. Therefore, we assume that $r_i\leq r$ in the sequel, implying that the number of global common trends can be relatively large.
Although the estimation procedures used to obtain the common stochastic trends and their associated loading matrices remain valid for unit-root series in finite samples, except for the methods to determine the number of stochastic trends, we introduce another set of identification conditions to maintain the dynamic patterns of all unit-root common factors. Specifically, in the estimation, we assume
The above identification conditions and Assumption (ref) imply that an additional strength of order $p_i^{1/2}$ (and $k_m^{1/2}$) is imposed on ${\mathbf g}_{i,t}$ (and ${\mathbf f}_{t}$). The estimation of ${\mathbf A}_i$ is based on the sample covariance $$\widehat\boldsymbol{\Sigma}_{y_i}=\frac{1}{T}\sum_{t=1}^T({\mathbf y}_{i,t}-\bar{{\mathbf y}}_i)({\mathbf y}_{i,t}-\bar{{\mathbf y}}_i)',$$ where $\bar{{\mathbf y}}_i=T^{-1}\sum_{t=1}^T{\mathbf y}_{i,t}$. It is a scaled version of ${\mathbf Y}_i'{\mathbf Y}_i$ in Section 2 if the data are centered first. Under Assumption (ref), by Theorems 1 and 2 of penaponcela2006, we have
for any column vector ${\mathbf a}_{i,j}$ of ${\mathbf A}_i$ and some constant $0<C<\infty$. The rate depends on whether there is a drift or not in the unit-root series. On the other hand, the constraints in ((ref)) imply that the eigenvalues of $\widehat\boldsymbol{\Sigma}_{g_i}=T^{-1}\sum_{t=1}^T({\mathbf g}_{it}-\bar{{\mathbf g}}_i)({\mathbf g}_{it}-\bar{{\mathbf g}}_i)'$ is bounded below by $Cp_iT$ or $Cp_iT^2$ depending on whether $E({\mathbf f}_t)= \mathbf{0}$ or not. We assume in the section that $E({\mathbf f}_t)=\mathbf{0}$ for simplicity.
We first assume that $r_i$ and $r$ are known. Similarly to the argument in Harris1997 that the $I(1)$ series components of ${\mathbf g}_{i,t}$ are the linear combinations of ${\mathbf y}_{i,t}$ that yield explosive variances, and the ones with minimum sample variances would be the “most stationary" components. Therefore, the estimated loading matrix $\widehat{\mathbf A}_i$ consists of the eigenvectors of $\widehat\boldsymbol{\Sigma}_{y_i}$ corresponding to the $r_i$ largest eigenvalues. The estimated factors are $\widehat{\mathbf g}_{i,t}=\widehat{\mathbf A}_i'{\mathbf y}_{i,t}$. Although the way to obtain $\widehat{\mathbf A}$ and $\widehat{\mathbf g}_{i,t}$ is different from the method of Section 2, but the extracted factor term $\widehat{\mathbf A}_i\widehat{\mathbf g}_{i,t}$ remains the same due to the identification issue.
Next, define the normalized group common factors $\widehat{\mathbf g}_{i,t}^*=\widehat{\mathbf g}_{i,t}/\sqrt{p_i}$ and the stacked group factors $\widehat{\mathbf g}_t^*=(\widehat{\mathbf g}_{1,t}^{*'},...,\widehat{\mathbf g}_{m,t}^{*'})'$. The estimation of ${\mathbf B}$ is based on the following sample covariance matrix
where $\bar{\widehat{\mathbf g}}^*=T^{-1}\sum_{t=1}^T\widehat{\mathbf g}_t^*$. Consequently, the estimator $\widehat{\mathbf B}$ consists of the $r$ eigenvectors of $\widetilde\boldsymbol{\Sigma}_{g}$ associated with the $r$ largest eigenvalues. The resulting estimated global stochastic trend is $\widehat{\mathbf f}_t=\widehat{\mathbf B}'\widehat{\mathbf g}_t^*$.
To estimate the number of common trends, we adopt the auto-correlation based method in gaotsay2020a. Specifically, for some prescribed integer $\bar{k}>0$, define
where $\widehat\rho_{i,l}(k)$ is the lag-$k$ sample autocorrelation function (ACF) of the principal components $\widehat g_{i,lt}$, for $1\leq i\leq m$ and $1\leq l\leq p_i$, and $\widehat\rho_j(k)$ is the corresponding one of $\widehat f_{j,t}$, for $1\leq j\leq k_m$. If $\widehat g_{i,lt}$ (or $\widehat f_{j,t}$) is stationary, then under some mild conditions, $\widehat\rho_{i,l}(k)$ (or $\widehat\rho_j(k)$) decays to zero exponentially as $k$ increases, and $\lim_{\bar{k}\rightarrow\infty} S_{i,l}(\bar{k})<\infty$ (or $\lim_{\bar{k}\rightarrow\infty} S_{j}(\bar{k})<\infty$) as $T\rightarrow\infty$. If $\widehat g_{i,lt}$ (or $\widehat f_{j,t}$) is unit-root nonstationary, then $\widehat\rho_{i,l}(k)\rightarrow 1$ (or $\widehat\rho_j(k)\rightarrow 1$), and $\lim_{\bar{k}\rightarrow\infty} S_{i,l}(\bar{k})=\infty$ (or $\lim_{\bar{k}\rightarrow\infty} S_{j}(\bar{k})=\infty$) as $T\rightarrow\infty$. Therefore, for the $i$-th group, we start with $l=1$. If the average of the absolute sample ACFs $S_{i,l}(\bar{k})/\bar{k}\geq \delta_0$ for some constant $0< \delta_0<1$, then $\widehat g_{i,lt}$ has a unit root and we increase $l$ by $1$ to repeat the detecting process. This detecting process is continued until $S_{i,l}(\bar{k})/\bar{k}< \delta_0$ or $l=p_i$. If $S_{i,l}(\bar{k})/\bar{k}\geq \delta_0$ for all $l$, then $\widehat r_i=p_i$; otherwise, we denote $\widehat r_i=l-1$. The estimate $\widehat r$ of the number of global common trends can similarly be obtained.
To study the theoretical properties of the estimators of the unit-root processes, we make some additional assumptions. Let ${\mathbf w}_t=(w_{1t},...,w_{rt})'={\mathbf f}_t-{\mathbf f}_{t-1}$ with $E({\mathbf w}_t)=\bf 0$, and define ${\mathbf S}_T^r({\mathbf t})=(S_T^1(t_1),...,S_T^r(t_r))'=\left(\frac{1}{\sqrt{n}}\sum_{s=1}^{[Tt_1]}w_{1s},...,\frac{1}{\sqrt{T}}\sum_{s=1}^{[Tt_{r}]}w_{rs}\right)'$, where $0< t_1<...<t_r\leq 1$ are constants and ${\mathbf t}=(t_1,...,t_{r})'$.
Assumptions (ref)(i)-(ii) are used to replace the ones in Assumptions (ref) and (ref)(ii) and Assumption (ref)(iii) is standard for unit-root processes. Assumption (ref) is not strong and it can be established under the setting of Stock1987, where we can assume ${\mathbf y}_{i,t}$ and ${\mathbf g}_t$ have similar structures as the one in (2.4) of Stock1987 and there are cointegration relationships among them.
We now state the convergence of the linear spaces spanned by the columns of the estimated loading matrices. Similar results as those in Theorems (ref) and (ref) continue to hold with a modified normalization in the proofs; see also bai2004 for a similar argument regarding the theoretical results with stochastic trend factors and the stationary ones. We do not repeat them here.
The consistency of the estimated numbers of the group stochastic trends and the global trends using the average of the absolute autocorrelations can be shown in a similar way as that in gaotsay2020a and the details are omitted.
In this section, we briefly discuss the extreme case in which the sample size $T$ is too large for a single machine to store the entire data of the $p_i$ time series, leading to the need to partition the data of those $p_i$ series. This situation is referred to as the horizontal-partition regime in the {\em Introduction}. Again, for simplicity, we focus on the 2-level case in Figure (ref), but assume that additional $n$ new machines are affiliated with each of the $m$ nodes, where $n$ depends implicitly on the sample size $T$, the number of nodes $m$ and the dimension $p_i$. We assume that the sample size $T$ is partitioned, in time order, as $T=T_1+\cdots+T_n$ such that the $T_j$ observations of the $p_i$-dimensional series can be analyzed by a single machine, for $i=1,\ldots,m$ and $j=1,\ldots,n$. Thus, we partition the $p_i$ dimensional series into $n$ blocks as $\{{\mathbf y}_{i,jt}|t=1,...,T_j; j=1,\ldots, n\}$, for $i=1,\ldots, m$, and store the data $\{{\mathbf y}_{i,jt}|t=1,\ldots,T_j\}$ in the $j$-th new machine affiliated with the $i$-th node. Assuming further that all the numbers of common factors $r_i$ and $r$ are known, we introduce next two approaches to estimating the common factors and their associated loading matrices.
The first approach applies the proposed method of Section (ref) to each sub-sample $\{{\mathbf y}_{i,jt},t=1,...,T_j\}$ and obtains the estimated loading matrix and factors, say, $\widehat{\mathbf A}_{i,j}$ and $\{\widehat{\mathbf g}_{i,jt},t=1,...,T_j\}$, for $j=1,\ldots, n$ and $i=1,\ldots,m$. Then we collect all the factors resulted from the $j$-th machine associated with the $i$-th node and form $\widehat{\mathbf G}_{i,j}=(\widehat{\mathbf g}_{i,j1},...,\widehat{\mathbf g}_{i,jT_j})'$. Define $\widehat{\mathbf G}_j^n=(\widehat{\mathbf G}_{1,j},...,\widehat{\mathbf G}_{m,j})$, which consists of all the group factors for the $j$-th segment, where the subscript $n$ signifies the partition of $n$ blocks. Next, we apply PCA again to $\widehat{\mathbf G}_j^n$ as in ((ref)), and denote the estimated loading matrix and factors by $\widehat{\mathbf B}_j^n$ and $\{\widehat{\mathbf f}_{j,t}: t=1,...,T_j\}$, respectively, which correspond to the $j$-th segment. Consequently, the estimated group common factors and global factors are $\{\widehat{\mathbf g}_{i,11,},...,\widehat{\mathbf g}_{i,1T_1},...,\widehat{\mathbf g}_{i,n1},...,\widehat{\mathbf g}_{i,nT_n}|i=1,...,m\}$ and $\{{\mathbf f}_{1,1},...,{\mathbf f}_{1,T_1},...,{\mathbf f}_{n,1},...,{\mathbf f}_{nT_n}\}$, respectively.
Let $\widehat{\mathbf A}_j^n=\mbox{diag}(\widehat{\mathbf A}_{1,j},...,\widehat{\mathbf A}_{m,j})$ be the block diagonal matrix consisting of all the estimated loading matrices from the $j$-th segment. Then, $\widehat{\mathbf A}_j^n$ is an estimate of ${\mathbf A}$ and $\widehat{\mathbf B}_j^n$ is an estimate of ${\mathbf B}$ based on the $j$-th segment. By Theorem (ref), we have the following corollary regarding the consistency of the loading matrices.
In fact, all the theoretical properties in Theorems (ref) to (ref) remain valid for the $j$-th segment so long as we replace $T$ therein by $T_j$. A key consideration for using each data segment to estimate ${\mathbf A}$ and ${\mathbf B}$ is that, in real applications with extremely long time series, the stationarity of ${\mathbf y}_{j,t}$ may fail. If this occurs, our first approach can extract more accurately the common factors using $\widehat{\mathbf A}_{i,j}$ for the $j$-th segment rather than using a globally estimated loading matrix $\widehat{\mathbf A}_i$.
The second approach is similar to that in fan2019 of the horizontal-partition regime. We adopt the identification conditions that ${\mathbf A}_i'{\mathbf A}_i={\mathbf I}_{r_i}$ and ${\mathbf B}'{\mathbf B}={\mathbf I}_r$. Letting $\boldsymbol{\Sigma}_{y_i,j}=\textnormal{Cov}({\mathbf y}_{i,jt})$ and $\widehat\boldsymbol{\Sigma}_{y_i,j}$ be its sample estimate, we denote $\widehat{\mathbf A}_{i,j}$ the matrix consisting of the eigenvectors of $\widehat\boldsymbol{\Sigma}_{y_i,j}$ associated with the $r_i$ largest eigenvalues. The estimated loading $\widehat{\mathbf A}_i$ for the entire data set of the $p_i$-dimensional time series is the matrix consisting of the eigenvectors associated with the $r_i$ largest eigenvalues of \[\widetilde\boldsymbol{\Sigma}_{y_i}^n=\frac{1}{n}\sum_{j=1}^n\widehat{\mathbf A}_{i,j}\widehat{\mathbf A}_{i,j}',\] which is similar to the model averaging technique in machine learning, and we compute the factors $\{\widehat{\mathbf g}_{i,t},t=1,...,T\}$ by making use of the final estimator $\widehat{\mathbf A}_i$ for the $i$-th group.
Next, we partition $\{\widehat{\mathbf g}_{i,t},t=1,...,T\}$ into $n$ segments and stack all the $j$-th segment associated with the $i$-th group to obtain $\widehat{\mathbf G}_j^n$ as before. Similarly to the way we obtain $\widehat{\mathbf A}_i$, letting $\widehat{\mathbf B}_j^n$ be the matrix consisting of the eigenvectors of the $r$ largest eigenvalues of $\widehat{\mathbf G}_j^n{'}\widehat{\mathbf G}_j^n/T_j$, we perform eigenvalue-eigenvector analysis on \[\widetilde\boldsymbol{\Sigma}_{g}^{n}=\frac{1}{n}\sum_{j=1}^n\widehat{\mathbf B}_j^n\widehat{\mathbf B}_j^n{'},\] and the final estimated loading matrix $\widehat{\mathbf B}$ consists of the eigenvectors associated with the $r$ largest eigenvalues of $\widetilde\boldsymbol{\Sigma}_{g}^{n}$.
Although the second procedure is similar to that in fan2019, the theoretical properties are much more involved as the data are time series with dynamically dependent structures and the theoretical results in fan2019 are only valid for i.i.d. data. Furthermore, the data involved in the second-step PCA may not be weakly dependent as $\widehat{\mathbf A}_i$ contains information across the entire time span. Therefore, the bias-variance argument in fan2019 is no longer valid under the current framework. We leave details of this approach for future research.
We use simulation and a real example to assess the performance of the proposed analysis in finite samples.
We study the finite-sample properties of the proposed methodology under the scenarios when $p$ is both small and large. As the dimensions of $\widehat{\mathbf A}_i$ and ${\mathbf A}_i$ are not necessarily the same, and ${\mathbf A}_i$ is not necessarily an orthogonal matrix, we extend the discrepancy measure in Equation ((ref)) to a more general one given below. Let ${\mathbf Q}_i$ be a $p\times d_i$ matrix with rank$({\mathbf Q}_i) = d_i$, and ${\mathbf P}_i = {\mathbf Q}_i({\mathbf Q}_i'{\mathbf Q}_i)^{-1} {\mathbf Q}_i'$, $i=1,2$, where we use the Moore-Penrose generalized inverse. Define
then $\bar{D} \in [0,1]$. Furthermore, $\bar{D}({\mathbf Q}_1,{\mathbf Q}_2)=0$ if and only if either $\mathcal{M}({\mathbf Q}_1)\subset \mathcal{M}({\mathbf Q}_2)$ or $\mathcal{M}({\mathbf Q}_2)\subset \mathcal{M}({\mathbf Q}_1)$, and it is 1 if and only if $\mathcal{M}({\mathbf Q}_1) \perp \mathcal{M}({\mathbf Q}_2)$. When $d_1 = d_2=d$ and ${\mathbf Q}_i'{\mathbf Q}_i= {\mathbf I}_d$, $\bar{D}({\mathbf Q}_1,{\mathbf Q}_2)$ reduces to $D({\mathbf Q}_1,{\mathbf Q}_2)$ of Equation ((ref)).
In the simulation study, we mainly focus on the numerical performance of the proposed method in the vertical partition regime for stationary and unit-root non-stationary time series data, because, in a horizontal partition regime, each block will produce similar results as those shown below, and we omit the details here to save space.
{\bf Example 1.} Consider the models in ((ref)) and ((ref)) with the global common factors following the VAR(1) model \[{\mathbf f}_t=\boldsymbol{\Phi}{\mathbf f}_{t-1}+\boldsymbol{\eta}_t,\] where $\boldsymbol{\eta}_t$ is a white noise process and $\boldsymbol{\Phi}$ is a diagonal matrix. We set the number of groups $m=10,\ 20$. For each $m$, the true number of global factors is $r=3$, the numbers of group factors are $r_1=...=r_m=5$, the group dimensions are $p_1=...=p_m=20,30,40$, or 50, and the sample size is $T=400,\ 800$, or 1200. For the coefficient matrices, we first set the seed in R to 1234 and the elements of ${\mathbf B}$ and ${\mathbf A}_{i}$ are drawn independently from $U(-2,2)$, and the diagonal elements of $\boldsymbol{\Phi}$ are drawn independently from $U(0.5,0.9)$. For each realization of ${\mathbf y}_t$, $\boldsymbol{\eta}_t \sim N(0,{\mathbf I}_{r})$, ${\mathbf u}_t\sim N(0,{\mathbf I}_{k_m})$, and ${\mathbf e}_{i,t}\sim N(0,{\mathbf I}_{p_i})$. We use $500$ replications for each $(p_i,m,T)$ configuration, where all the $p_i$'s are set to be equal.
We first consider the performance of estimating the number of factors, and use the IC1 criterion with $K=[p_i/3]$ to determine $r_i$ and the ratio-based method in ((ref)) with $R=[k_m/2]$ to estimate $r$ . The empirical probabilities of $\min_{1\leq i\leq m}P(\widehat r_i=r_i)$ and $P(\widehat r=r)$ are reported in Table (ref) for $m=10$ and $20$. From the table, we see that the proposed methods can successfully estimate the numbers of group factors using the BIC and that of the global factors using the ratio-based procedure. This is understandable because the factors used are all strong ones. Similar results can also be found in BaiNg_Econometrica_2002 and lamyao2012.
Next, we study the estimation accuracy of the loading matrices. We only show the results for $m=10$ since the results are similar for $m$ = 20. Define
which measure the average and the maximum estimation errors between the linear spaces spanned by $\widehat{\mathbf A}_i$ and ${\mathbf A}_i$, respectively. The boxplots of $\bar{D}_1(\widehat{\mathbf A},{\mathbf A})$ and $\bar{D}_2(\widehat{\mathbf A},{\mathbf A})$ are presented in Figure (ref). From Figure (ref)(a)--(b), we see that, for each $p_i$, the discrepancy decreases as the sample size increases and this is in agreement with our asymptotic theory. For the estimation errors of the loading matrix associated with the global factors, we further define
which measures the average discrepancy between $\widehat{\mathbf A}_i\widehat{\mathbf B}_i$ and ${\mathbf A}_i{\mathbf B}_i$ in each group. The boxplots of $\bar{D}_3(\widehat{\mathbf A}\widehat{\mathbf B},{\mathbf A}{\mathbf B})$ are shown in Figure (ref). From the figure, we also see that, for each $p_i$, the discrepancy between the estimated global loading matrix and the true ones decreases as the sample size increases and the result is also in agreement with our theorems.
Finally, we study the estimation errors of the extracted factors. We only consider the case when $m=10$ since the results are similar for $m$ = 20. For each $(p_i,T)$, define the root-mean-squared error (RMSE) as
which quantify the accuracy in estimating the group factors and the global factors, respectively. Boxplots of the RMSEs are shown in Figure (ref). From the plot, we see clearly that, as the sample size increases, the RMSE decreases for a given $p_i$, which is consistent with the results of Theorem (ref). Also, as expected from the same theorem, the RMSE decreases with the dimension $p_i$. Overall, our proposed method works well for various settings of moderately large dimensions and sample sizes.
Next, we apply the BIC-type criterion IC1 of Section 2.2 as suggested by BaiNg_Econometrica_2002 and the ratio-based method in ahn2013 and lamyao2012 to the entire data set. Our goal here is to illustrate the impact of overlooking the underlying factor structure of the model in data analysis. The boxplots of the estimated numbers of common factors by these two methods are shown in Figure (ref). From the plots, we see that, overlooking the hierarchical factor structure, the BIC criterion of BaiNg_Econometrica_2002 is able to identify the total number of the group factors while ignoring the number of global common factors, and the ratio-based method tends to identify the number of global factors while ignoring the total number of group factors. On the other hand, as shown in Table (ref), our hierarchical framework works well to identify both numbers of common factors.
\vskip 0.5cm
{\bf Example 2.} In this example, we consider the models in ((ref)) and ((ref)) with global common factors following the random walk model \[{\mathbf f}_t={\mathbf f}_{t-1}+\boldsymbol{\eta}_t,\] where $\boldsymbol{\eta}_t$ is a white noise process. We set $r_i=3$ and $r=5$, and everything else is generated in the same way as that of Example 1.
To study the performance of estimating the number of stochastic trends, we use the the autocorrelation-based method in ((ref)) with $\delta_0=0.3$ and $\bar{k}=20$ to estimate $r_i$ and $r$. The empirical probabilities of $\min_{1\leq i\leq m}P(\widehat r_i=r_i)$ and $P(\widehat r=r)$ are reported in Table (ref) for $m=10$. Similar results are also found for the other choice of $m$, and we omit the detail to save space. From Table (ref), we see that the proposed method works well in selecting the numbers of group common stochastic trends and the global ones. Overall, as the sample size increases, the proposed method shows improved performance in estimating the numbers of common factors, which is in agreement with the asymptotic theory. In addition, the performance of the proposed method seems to be stable or even improved as the dimension increases, showing evidence of the `blessing of dimensionality' in high-dimensional PCA estimation. The boxplots of the empirical discrepancies between the estimated loading matrices and the true ones are similar to those in Example 1 because the estimated numbers of factors are quite accurate according to Table (ref). Therefore, we do not repeat them here.
In this section, we demonstrate the application of the proposed method using U.S. monthly macroeconomic data with 157 time series and 480 observations.
{\bf Example 3.} Consider the data set of U.S. monthly macroeconomic time series available from January 1959 to December 1998, which was used in StockWatson_2002a. The description of the original data set consisting of 215 series was given in StockWatson_2002a, but we eliminate the series with missing data and the series in the categories of average hourly earnings and miscellaneous, because all, but one or two, series in these two categories contain missing values. The resulting data set consists of 157 series with $480$ observations. To achieve stationarity, we take the first difference of all series and, hence, we employ the transformed data with $p=157$ and $T=479$.
The series can roughly be classified into $12$ groups (or clusters) as indicated in StockWatson_2002a: (1) Real output and income (Out, series 1--21), (2) Employment and hours (EMP, series 22--48), (3) Retails and manufacturing trade (RTS, series 49--57), (4) Consumption (PCE, series 58--62), (5) Housing starts and sales (HSS, 53--69), (6) Real inventories and inventory-sales ratios (Inv, series 70--80), (7) Orders and unfilled orders (Ord, series 81--96), (8) Stock prices (Spr, series 97--103), (9) exchange rates (EXR, series 104--109), (10) Interest rates (Int, series 110--128), (11) Money and credit quantity aggregates (Mon, series 129--137), (12) Price indexes (Pri, series 138--157). Details are given in StockWatson_2002a and the series of 12 groups are shown in Figure (ref) with $m=12$.
We first applied our proposed method to each group and found that the numbers of common factors in each group are (1) $r_1=7$, (2) $r_2=9$, (3) $r_3=3$, (4) $r_4=2$, (5) $r_5=3$, (6) $r_6=4$, (7) $r_7=6$, (8) $r_8=3$, (9) $r_9=2$, (10) $r_{10}=7$, (11) $r_{11}=3$, and (12) $r_{12}=7$. The total number of group factors is $k_m=56$. The time plots of all group factors are shown in Figure (ref). We then performed the second stage factor analysis. The first $26$ eigenvalues ($\widehat\lambda_i$) of $\widehat{\mathbf G}'\widehat{\mathbf G}$ and the associated eigenvalue ratios are presented in Figure (ref). From the plots, we see that the sharpest drop occurs at the first ratio, which implies that there is one global common factor and, hence, $\widehat r=1$. For ease in illustration, we show the first two estimated global factors in Figure (ref). These factors capture the common cross-sectional dependence between the group factors of series in Figure (ref).
Next, we compare the forecasting ability of our global and group-specific factors with the ones obtained using the whole $157$ series jointly. Observations of the last ten years are used for out-of-sample forecasting. When using the information criterion method of BaiNg_Econometrica_2002 to select the number of factors, we find that $\widehat r$ increases as we increase the upper bound $K$ in ((ref)), and it eventually reaches a value greater than $50$, so we adopt the ratio-based method in ahn2013 and lamyao2012, as specified in a similar fashion as ((ref)). The first 25 eigenvalues, $\widehat\mu_i$, of the covariance of the full data $\widehat\boldsymbol{\Sigma}_y$ and the ratios of the consecutive eigenvalues are given in Figure (ref). The sharpest drop of the ratios occurs at $\widehat r_{full}=3$. For comparisons of the forecasting ability of different factors, we investigate the 1-step to 4-step ahead forecasts of eight commonly studied economic indicators ($x_{i,t}, 1\leq i\leq 8$). They are the industrial production (ip), real personal income less transfers (gmyxpq), real manufacturing trade and sales (msmtq), number of employees on non-agricultural payrolls (lpnag), the consumer price index (CPI) (punew), the personal consumption expenditure deflator (gmdc), the consumer price index less food and energy (puxf), and the producer price index for finished goods (pwfsa); see, also, StockWatson_2002a.
To see the forecastability of the extracted factors, we only regress $x_{i,t+h}$ on $\widehat{\mathbf f}_t$ and $\widehat{\mathbf u}_{i,t}$ without including the lagged variable of $x_{i,t}$, because StockWatson_2002a found that the simple diffusion index forecasts excluding the lags are even better than that with the lagged variables. Table (ref) reports the $1$-step to $4$-step ahead mean squared forecast errors for the eight selected series described above. We adopt five models according to the included factors in ((ref)). They are (1) 1glb: including the first global factor in Figure (ref)(a) as a predictor; (2) 2glb: including the two global factors in Figure (ref)(a,b) as predictors; (3) (1glb+spf): including the first global factor and the associated group-specific factors ${\mathbf u}_{i,t}$ of $x_{i,t}$; (4) 2glb+spf: including the two global factors and the associated group-specific factors ${\mathbf u}_{i,t}$ of $x_{i,t}$ (2glb+spf); (5) SW: including the 3-dimensional factor process obtained from the full data set as that in StockWatson_2002a. From Table (ref), we see that the smallest error was never achieved by the SW method except the 2-step ahead prediction for gmyxpq. In addition, the forecastability of the first global factor identified by our ratio-based method performs better than that of the factors recovered using the full data set by the method of StockWatson_2002a, except for the 1-step ahead prediction of msmtq; the method of 2glb outperforms SW except for the 1-step ahead forecast of msmtq, and the 2-step ahead forecast for gmyxpq; the method of 1glb+spf outperforms SW except for the 1-step ahead forecast for gmyxpq, the 2-step ahead forecasts of gmyxpq and pwfsa, and the 4-step ahead forecasts for ip, gmyxpq, and pwfsa; the method of 2glb+spf outperforms SW except for the 1-step ahead forecast of gmyxpq, the 2-step ahead forecast for gmyxpq and pwfsa, and the 4-step ahead forecasts for ip and pwfsa. Overall, the factors extracted by the proposed method tend to have higher forecastability than those obtained by the full data set in predicting the selected eight economic variables.
This article introduced a distributed factor approach to modeling large-scale stationary and unit-root non-stationary time series data. In the modern era of data science, large data sets are usually scattered across different locations such that it is extremely difficult to fuse or aggregate them due to the issues of communication cost, privacy concern, data security, and ownership, among others. The proposed divide-and-conquer method can overcome this difficulty by performing the factor modeling in parallel, and the resulting factors are of lower dimensions and easier to transfer, and the privacy can also be protected. The numerical results suggest that the factors produced by our proposed method have higher forecasting ability than those obtained by applying the traditional PCA to the full data set. The proposed method thus provides another option in the toolkit for the practitioners who are interested in out-of-sample forecasting with large-scale data sets. We also briefly discussed the analysis of extremely large-scale time series data for which both the horizontal and vertical partitions are needed. This topic deserves a careful further investigation.