The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
94,315 characters
Divide-and-Conquer: A Distributed Hierarchical Factor Approach to Modeling Large-Scale Time Series Data
\title{Divide-and-Conquer:
A Distributed Hierarchical Factor Approach to Modeling Large-Scale Time Series Data}
\author{
Zhaoxing Gao$^1$ and Ruey S. Tsay$^2$ \\
$^1$Department of Mathematics, Lehigh University\\
$^2$Booth School of Business, University of Chicago
}
\date{}
\maketitle
\begin{abstract}
This paper proposes a hierarchical approximate-factor approach to analyzing high-dimensional,
large-scale heterogeneous time series data using distributed computing. The new method employs a multiple-fold dimension reduction procedure using Principal Component Analysis (PCA) and shows great promises for modeling large-scale data that cannot be stored nor analyzed by a single machine.
Each computer at the basic level performs a PCA to extract common
factors among the time series assigned to it and transfers those factors to one and only
one node of the second level. Each 2nd-level computer collects the common factors from its subordinates and performs another PCA
to select the 2nd-level common factors. This process is repeated until the central server is reached, which collects
common factors from its direct subordinates and performs a final PCA to
select the global common factors. The noise terms of the 2nd-level approximate factor model
are the unique common factors of the 1st-level clusters.
We focus on the case of 2 levels in our theoretical derivations, but the idea can easily be
generalized to any finite number of hierarchies. We discuss
some clustering methods when the group memberships are unknown
and introduce a new diffusion index approach to forecasting. We further extend the analysis
to unit-root nonstationary time series. Asymptotic properties of the proposed method are derived for the diverging dimension of the data in each computing unit
and the sample size $T$. We use both simulated data and real examples to assess the performance of the proposed method in finite samples, and compare our method with the commonly used ones in the literature concerning the forecastability of extracted factors.
\end{abstract}
\noindent {\sl Keywords}: Distributed Factor Model, Large-scale Time Series, Principal Component Analysis, Commonality, Heterogeneity
\newpage
\section{Introduction}
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 \cite{ross1976}, \cite{StockWatson_1989}, \cite{famafrench1993,famafrench2015}, and \cite{StockWatson_1998}. Macroeconomic variables of multiple countries are often found to have comovements, see \cite{GregoryHead1999} and \cite{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 \cite{chamberlain1983}, \cite{StockWatson_1989,StockWatson_2002a,StockWatson_2002b}, \cite{BaiNg_Econometrica_2002}, and \cite{Bai_Econometrica_2003}; the dynamic factor models in \cite{forni2000,forni2005}; the factor models with dynamically dependent factors in \cite{lamyao2012} and \cite{gaotsay2018a,gaotsay2020a, gaotsay2020b}, among others. Other dimension reduction techniques for time series with a low dimensional structure are also available,
see \cite{BoxTiao_1977}, \cite{TiaoTsay_1989}, \cite{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 \cite{zhang2013} and \cite{chenxie2014}, the distributed PCA in \cite{qu2002} and \cite{fan2019} with a horizontal partition regime, where each server contains a small portion of the data of all subjects,
and another distributed PCA in \cite{kargupta2001}, \cite{lietal2011}, and
\cite{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.
\begin{figure}
\begin{center}
{\includegraphics[width=0.95\textwidth]{pca3er.png}}
\caption{A 3-layer sketch of distributed computing framework for the proposed hierarchical factor modeling of large-scale time series data.}\label{fig0}
\end{center}
\end{figure}
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{fig0}, 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{fig0}, 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 \cite{fan2019},
and it is also different from the one in \cite{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 \cite{lamyao2012} and \cite{ahn2013} and the BIC-type criterion in \cite{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 \cite{Moench2013}, the dynamic factor model in \cite{alonso2020} with cluster structure, and the panel data model in \cite{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 \cite{Moench2013} and \cite{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 \cite{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{sec2} and
study the theoretical properties of the proposed model and its associated
estimates in Section \ref{sec3}. The proposed method is extended to unit-root nonstationary time series data in Section \ref{sec4}. Section \ref{sec5} 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{sec6}, and
Section \ref{sec7} 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)$.
\section{Models and Methodology}\label{sec2}
\subsection{Setting}
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:
\begin{equation}\label{m-factor}
{\mathbf y}_{i,t}={\mathbf A}_i{\mathbf g}_{i,t}+{\mathbf e}_{i,t},\quad i=1,...,m,
\end{equation}
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
\begin{equation}\label{global-f}
{\mathbf g}_t={\mathbf B}{\mathbf f}_t+{\mathbf u}_t,
\end{equation}
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{m-factor}) and (\ref{global-f}) together,
we obtain the following matrix form of the proposed 2-level model:
\begin{equation}\label{full:m}
{\mathbf y}_t={\mathbf A}{\mathbf B}{\mathbf f}_t+{\mathbf A}{\mathbf u}_t+{\mathbf e}_t,
\end{equation}
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{global-f}) 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 \cite{andobai2017} and \cite{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{m-factor}) and (\ref{global-f}) are not uniquely identified. We follow the restrictions used in \cite{BaiNg_Econometrica_2002} and \cite{fan2013} and assume that
\begin{equation}\label{constr}
\textnormal{Cov}({\mathbf g}_{i,t})={\mathbf I}_{r_i},\,\,\textnormal{Cov}({\mathbf f}_t)={\mathbf I}_r,\,\, {\mathbf A}_i'{\mathbf A}_i\,\,\text{and}\,\, {\mathbf B}'{\mathbf B}\,\,\text{are diagonal},\,\, i=1,...,m.
\end{equation}
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.
\subsection{Estimation of Distributed Principal Components}
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:
\begin{equation}\label{est-g}
(\widehat{\mathbf A}_i,\widehat{\mathbf g}_{i,t})=\arg\min_{{\mathbf A}_i\in \mathbb{R}^{p_i\times r_i}, {\mathbf g}_{i,t}\in \mathbb{R}^{r_i}}\sum_{t=1}^T\|{\mathbf y}_{i,t}-{\mathbf A}_i{\mathbf g}_{i,t}\|_2^2,
\end{equation}
subject to the normalization conditions:
\begin{equation}\label{norm-g}
\frac{1}{T}\sum_{t=1}^T{\mathbf g}_{i,t}{\mathbf g}_{i,t}'={\mathbf I}_{r_i},\,\,\text{and}\,\, {\mathbf A}_i'{\mathbf A}_i\,\,\text{is diagonal}.
\end{equation}
Similarly to \cite{BaiNg_Econometrica_2002} and \cite{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{est-g}) 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
\begin{equation}\label{est-c}
(\widehat{\mathbf B},\widehat{\mathbf F})=\arg\min_{{\mathbf B}\in \mathbb{R}^{k\times r},{\mathbf F}\in \mathbb{R}^{T\times r}}\|\widehat{\mathbf G}-{\mathbf F}{\mathbf B}'\|_F^2,
\end{equation}
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 \cite{BaiNg_Econometrica_2002}. It estimates $r_i$ by
\begin{equation}\label{ri}
\widehat r_i=\arg\min_{0\leq K_i\leq K}\log\left(\frac{1}{p_iT}\|{\mathbf Y}_i-\widehat{\mathbf G}_{i,K_i}\widehat{\mathbf G}_{i,K_i}'{\mathbf Y}_i/T\|_F^2\right)+K_ig(T,p_i),
\end{equation}
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 \cite{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{ri}), we also adopt the ratio-based method in \cite{lamyao2012} and \cite{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
\begin{equation}\label{est-r}
\widehat r=\arg\min_{1\leq l\leq R}{\frac{\widehat\lambda_{l+1}}{\widehat\lambda_{l}}},
\end{equation}
where we take $R=[k_m/2]$ as suggested by \cite{lamyao2012}. Alternatively, we may estimate the number of global factors by selecting a threshold $0<c_0<1$ such that
\begin{equation}\label{var-exp}
\widehat r= \arg \min_j\left\{\frac{\sum_{l=1}^j\widehat\lambda_{l}}{\sum_{l=1}^{k_m}\widehat\lambda_{l}}\geq c_0\right\},
\end{equation}
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\%$).
\subsection{Clustering Time Series If Needed}
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, \cite{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, \cite{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 \texttt{1234} in \texttt{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
\begin{equation*}\label{pca:v}
{\mathbf P}_5=\left(\begin{array}{rrrrr}
0.451 &-0.332& 0.613 &0.093 & 0.550\\
-0.456 & 0.310 & 0.476& -0.670 & 0.144\\
-0.533 & 0.001 & 0.483 & 0.661 & -0.213\\
-0.147 & -0.833 & 0.084& -0.316& -0.423\\
-0.532& -0.318& -0.396 & 0.078 & 0.673
\end{array}\right),
\end{equation*}
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:
\begin{equation*}\label{pca:vc}
{\mathbf P}_5^*=\left(\begin{array}{rrrrr}
0.999& -0.043 & 0.011 & 0.010& -0.002\\
-0.024 & -0.468 & 0.521& -0.311& -0.642\\
-0.025 &-0.604 &0.243 &-0.195 & 0.733\\
0.001& -0.299 &-0.767& -0.550 &-0.139\\
-0.029& -0.570& -0.284 & 0.750 &-0.177\\
\end{array}\right).
\end{equation*}
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.
\subsection{Prediction}
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
\begin{equation}
\label{df-id}
\widehat x_{t+h}=\widehat\alpha_h+\widehat\boldsymbol{\beta}_h(L)\widehat{\mathbf f}_t+\widehat\boldsymbol{\gamma}_h(L)\widehat{\mathbf u}_{i,t}+\widehat\phi_h(L)x_{t},
\end{equation}
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{df-id}), we may include all the group-specific factors and using group penalized method to estimate the coefficients such as the group Lasso in \cite{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.
\section{Theoretical Properties}\label{sec3}
We begin with some assumptions for Models (\ref{m-factor}) to (\ref{full:m}), where only $\{{\mathbf y}_{it}|1\leq i\leq m,1\leq t\leq T\}$ are observable. Similar ones are also given in \cite{BaiNg_Econometrica_2002} and \cite{fan2013}. All proofs of the theorems are given in
the Appendix.
\begin{assumption}\label{asm1}
All the eigenvalues of $p_i^{-1}{\mathbf A}_i'{\mathbf A}_i$ and $k_m^{-1}{\mathbf B}'{\mathbf B}$ are bounded away from both $0$ and $\infty$, as $p_i, k_m\rightarrow\infty$.
\end{assumption}
\begin{assumption}\label{asm2}
The process $\{{\mathbf y}_{it},{\mathbf g}_{it}, {\mathbf e}_{i,t},{\mathbf f}_{t}\}$ is stationary with zero mean, and $\alpha$-mixing with the mixing coefficient $\alpha_p(k)\leq \exp(-Ck^{\gamma})$ for some positive constants
$C$ and $\gamma>0$, where
\[\alpha_p(k)=\sup_{i}\sup_{A\in\mathcal{F}_{-\infty}^i,B\in \mathcal{F}_{i+k}^\infty}|P(A\cap B)-P(A)P(B)|,
\]
and $\mathcal{F}_i^j$ is the $\sigma$-field generated by $\{({\mathbf y}_{it},{\mathbf g}_{it}, {\mathbf e}_{i,t},{\mathbf f}_{t}):i\leq t\leq j\}$.
\end{assumption}
\begin{assumption}\label{asm3}
(i) There exist constants $c_1,c_2>0$ such that $\lambda_{\min}(\boldsymbol{\Sigma}_{e_i})>c_1$, $\lambda_{\min}(\boldsymbol{\Sigma}_{u})>c_1$, $\|\boldsymbol{\Sigma}_{e_i}\|_1<c_2$, and $\|\boldsymbol{\Sigma}_u\|_1<c_2$, where $\boldsymbol{\Sigma}_{e_i}=\textnormal{Cov}({\mathbf e}_{i,t})$ and $\boldsymbol{\Sigma}_u=\textnormal{Cov}({\mathbf u}_t)$;\\
(ii) There exists a constant $C>0$ such that ${\mathbf g}_{i,t}$, ${\mathbf f}_t$, ${\mathbf e}_{i,t}$, and ${\mathbf u}_t$ are sub-exponentially distributed in the sense that
\[P(|{\mathbf v}'({\mathbf x}_t-E({\mathbf x}_t))|>x)\leq C\exp(-Cx),\,\,\text{for}\,\, x>0,\]
where $\|{\mathbf v}\|_2=1$ and ${\mathbf x}_t$ is any vector of ${\mathbf g}_{i,t}$, ${\mathbf f}_t$, ${\mathbf e}_{i,t}$, and ${\mathbf u}_t$ .
\end{assumption}
\begin{assumption}\label{asm4}
There exists an $M>0$ such that, for all $l, k\leq p_i$, $j,\iota\leq k_m$,
$s, t\leq T$, $1\leq i\leq m$,\\
(i) $\|{\mathbf a}_{i,l}\|_\infty<M$ and $\|{\mathbf b}_{j}\|_\infty<M$, where ${\mathbf a}_{i,l}$ and ${\mathbf b}_j$ are, respectively, the $l$-th and the $j$-th rows of ${\mathbf A}_i$ and ${\mathbf B}$; \\
(ii) $E\{p_i^{-1/2}[{\mathbf e}_{i,s}'{\mathbf e}_{i,t}-E({\mathbf e}_{i,s}'{\mathbf e}_{i,t})]\}^4<M$, $E\{k_m^{-1/2}[{\mathbf u}_{s}'{\mathbf u}_{t}-E({\mathbf u}_{s}'{\mathbf u}_{t})]\}^4<M$;\\
(iii) $E\|p_i^{-1/2}\sum_{l=1}^{p_i}{\mathbf a}_{i,l}e_{i,lt}\|_2^4<M$ and $E\|k_m^{-1/2}\sum_{j=1}^{k_m}{\mathbf b}_{j}u_{jt}\|_2^4<M$;\\
(iv) $\sum_{l=1}^{p_i}|E(e_{i,lt}e_{i,kt})|<M$ and $\sum_{j=1}^{k_m}|E(u_{jt}u_{\iota t})|<M$;\\
(v) $E\|\frac{1}{\sqrt{p_iT}}\sum_{t=1}^T\sum_{l=1}^{p_i}{\mathbf g}_{i,t}[e_{i,ls}e_{i,lt}-E(e_{i,ls}e_{i,lt})]\|_2^2<M$ and $E\|\frac{1}{\sqrt{k_mT}}\sum_{t=1}^T\sum_{j=1}^{k_m}{\mathbf f}_{i,t}[u_{js}u_{jt}-E(u_{js}u_{jt})]\|_2^2<M$.
\end{assumption}
Assumption \ref{asm1} implies that the factors used are pervasive as those in \cite{BaiNg_Econometrica_2002} and \cite{fan2013}. Assumption \ref{asm2} is standard for dependent random processes, and it implies the temporal dependence conditions in \cite{BaiNg_Econometrica_2002} and \cite{Bai_Econometrica_2003}. See \cite{gaoetal2017} for a theoretical justification of the mixing conditions for VAR models. Assumption \ref{asm3}(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{est-r}), and (ii) implies that the factors and the idiosyncratic terms all have exponential tails, which is stronger than the moment conditions in \cite{BaiNg_Econometrica_2002}. Assumption \ref{asm4} controls the cross-sectional dependence of the idiosyncratic terms.
Under Assumptions \ref{asm1}--\ref{asm3} and \ref{asm4}(i)-(iii), all the assumptions used in Bai and Ng (2002) are satisfied. See also \cite{fan2013}. Assumptions \ref{asm4}(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 \cite{Bai_Econometrica_2003},
and the proof of Theorem 1 in \cite{BaiNg_joe_2006} is valid under these assumptions.
\begin{theorem}\label{tm1}
Suppose Assumptions \ref{asm1} to \ref{asm4} hold and $r_i$ and $r$ are known. Let $\widehat{\mathbf g}_{i,t}$ and $\widehat{\mathbf f}_{t}$ be the principal components estimators of ${\mathbf g}_{i,t}$ and ${\mathbf f}_t$. Then, as $p_i, T\rightarrow\infty$, there exist full-rank matrices ${\mathbf H}_i\in \mathbb{R}^{r_i\times r_i}$ and ${\mathbf K}\in \mathbb{R}^{r\times r}$ such that
\begin{equation}\label{rate:g}
\frac{1}{T}\sum_{t=1}^T\|\widehat{\mathbf g}_{i,t}-{\mathbf H}_i'{\mathbf g}_{i,t}\|_2^2=O_p\left(\frac{1}{p_i}+\frac{1}{T}\right),\quad i=1,...,m
\end{equation}
and
\begin{equation}\label{rate:f}
\frac{1}{T}\sum_{t=1}^T\|\widehat{\mathbf f}_{t}-{\mathbf K}'{\mathbf f}_{t}\|_2^2=O_p\left(\frac{1}{k_m}+\frac{1}{T}+\frac{m}{k_m\underline{p}}+\frac{m^2}{{T^2}}+\frac{m^2}{{\underline{p}^2}}\right),
\end{equation}
where $\underline{p}=\min\{p_1,...,p_m\}$.
Furthermore, under the same assumptions, we have
\begin{equation}\label{max:rate1}
\max_{1\leq t\leq T }\|\widehat{\mathbf g}_{i,t}-{\mathbf H}_i'{\mathbf g}_{i,t}\|_2=O_p\left(\frac{1}{\sqrt{T}}+\frac{T^{1/4}}{\sqrt{p_i}}\right),
\end{equation}
and
\begin{equation}\label{max:rate2}
\max_{1\leq t\leq T }\|\widehat{\mathbf f}_{t}-{\mathbf K}'{\mathbf f}_{t}\|_2=O_p\left\{\sqrt{\frac{m}{k_m}}\log(T)\left[\frac{1}{\sqrt{T}}+\frac{1}{\sqrt{\underline{p}}}+\frac{m}{T}+\frac{m}{\underline{p}}\right]+\sqrt{\frac{m}{k_m}}\frac{T^{1/4}}{\sqrt{\underline{p}}}+\frac{T^{1/4}}{\sqrt{k_m}}\right\}.
\end{equation}
\end{theorem}
Some remarks of Theorem~\ref{tm1} are in order.
\begin{remark}\label{rm1}
(i) Theorem 1 implies that the convergence rates of the second step factor modeling depend on the number of clusters $m$, which is comparable with the rate of $k_m$ if the number of common factors in each group is small and fixed;\\
(ii) For the consistency of the estimated factors in the sense of (\ref{rate:g}) and (\ref{rate:f}) as $p_i,T\rightarrow\infty$, we require that $m\rightarrow\infty$ but $m=o(T,\underline{p})$, which is a reasonable setting as the number of clusters is not extremely large in practice;\\
(iii) From (\ref{max:rate1}) and (\ref{max:rate2}), the uniform convergence can be achieved so long as $T=o(\underline{p}^2, m^2)$ and $m=o(T,\underline{p})$ as $p_i,T\rightarrow\infty$ under the assumption that $m=O(k_m)$. The conditions for the uniform convergence in (\ref{max:rate1}) and (\ref{max:rate2}) are stronger than those of (\ref{rate:g}) and (\ref{rate:f}) as they involve the maximums over the length of the time series $T$. However, the uniform rates do not provide a consistency pattern as $T\rightarrow\infty$, but $p_i$ and $m$ are large and fixed. We consider the discrepancies between the estimated loading matrices and factors with the true ones in Theorem \ref{tm3} below, which shows a consistency pattern as $T\rightarrow\infty$; \\
(iv) By \cite{BaiNg_joe_2006}, under the condition that all the eigenvalues of $p_i^{-1}{\mathbf A}_i'{\mathbf A}_i$ and $k_m^{-1}{\mathbf B}'{\mathbf B}$ are distinct, we can show that ${\mathbf H}_i\rightarrow_p {\mathbf I}_{r_i}$ and ${\mathbf K}\rightarrow_p {\mathbf I}_r$ up to a change of signs in the diagonals. We do not pursue this issue here.
\end{remark}
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.
\begin{theorem}\label{tm2}
Suppose Assumptions \ref{asm1} to \ref{asm4} hold, and $r_i$ and $r$ are known. Then, for $1\leq i\leq m$,
\begin{equation}\label{max-loading1}
\max_{1\leq l\leq p_i}\|\widehat{\mathbf a}_{i,l}-{\mathbf H}_i'{\mathbf a}_{i,l}\|_2=O_p\left(\sqrt{\frac{\log(p_i)}{T}}+\frac{1}{\sqrt{p_i}}\right),
\end{equation}
and
\begin{equation}\label{max-loading2}
\max_{1\leq l\leq p_i}\|\widehat{\mathbf c}_{i,l}-{\mathbf K}'{\mathbf c}_{i,l}\|_2=O_p\left(\sqrt{\frac{1}{k_m}+\frac{m}{k_m\underline{p}}+\frac{m^2}{{T^2}}+\frac{m^2}{{\underline{p}^2}}}+\sqrt{\frac{\log(p_i)}{T}}\right),
\end{equation}
where $\underline{p} = \min\{p_1,\ldots,p_m\}$. Consequently,
\begin{equation}\label{max:gt}
\max_{l\leq p_i,t\leq T}\|\widehat{\mathbf a}_{i,l}'\widehat{\mathbf g}_{i,t}-{\mathbf a}_{i,l}'{\mathbf g}_{i,t}\|_2=O_p\left(\log(T)\sqrt{\frac{\log(p_i)}{T}}+\frac{T^{1/4}}{\sqrt{p_i}}\right),
\end{equation}
and
\begin{equation}\label{max:ft}
\max_{l\leq p_i, t\leq T}\|\widehat{\mathbf c}_{i,l}'\widehat{\mathbf f}_{i,t}-{\mathbf c}_{i,l}'{\mathbf f}_{i,t}\|_2=O_p\left(\log(T)w_T+\sqrt{\frac{m}{k_m}}\frac{T^{1/4}}{\sqrt{\underline{p}}}+\frac{T^{1/4}}{\sqrt{k_m}}\right),
\end{equation}
where $w_{T}=\sqrt{\frac{m}{k_m\underline{p}}+\frac{m^2}{{T^2}}+\frac{m^2}{{\underline{p}^2}}}+\sqrt{\frac{\log(p_i)}{T}}$.
\end{theorem}
\begin{remark}
(i) Theorem 2 implies that, the uniform convergences of the estimated factor loading matrices and the extracted factors depend on the number of groups together with the dimension $p_i$ and the sample size $T$. Under the conditions specified in Remark \ref{rm1}, we can also achieve the consistency in Theorem \ref{tm2};\\
(ii) The convergence rate of $\|\widehat{\mathbf a}_{i,l}'\widehat{\mathbf u}_{i,t}-{\mathbf a}_{i,l}'{\mathbf u}_{i,t}\|_2$ follows from (\ref{max:gt}) and (\ref{max:ft}), and we do not repeat it here.
\end{remark}
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 \cite{gaotsay2020a,gaotsay2020b} for details.
To this end, we adopt the discrepancy measure used by
\cite{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
\begin{equation}
D({\bf Q}_1,{\bf
Q}_2)=\sqrt{1-\frac{1}{d}\textrm{tr}({\bf Q}_1{\bf Q}_1'{\bf
Q}_2{\bf Q}_2')}.\label{eq:D}
\end{equation}
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 \cite{stewart-sun1990}. The squared Frobenius norm of the so-called $\sin\boldsymbol{\Theta}$ distance is
\begin{equation}\label{sin:t}
\|\sin\boldsymbol{\Theta}({\mathbf Q}_1,{\mathbf Q}_2)\|_F^2:=\sum_{i=1}^d\sin^2(\theta_i)=\sum_{i=1}^d(1-\sigma_i^2)=d\{D({\bf Q}_1,{\bf
Q}_2)\}^2,
\end{equation}
which implies that the above two measures are equivalent for finite $d$.
Therefore, we shall use only the distance in (\ref{eq:D}) 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$.
\begin{theorem}\label{tm3}
Under the assumptions of Theorem \ref{tm1}, for $1\leq i\leq m$, we have
\begin{equation}\label{DA:DAB}
D(\widehat{\mathbf A}_i,{\mathbf A}_i)=O_p(T^{-1/2})\quad \text{and}\quad D(\widehat{\mathbf A}\widehat{\mathbf B},{\mathbf A}{\mathbf B})=O_p\left(mT^{-1/2}+\underline{p}^{-1}\bar{p}T^{-1/2}+\underline{p}^{-1}\right).
\end{equation}
Furthermore,
\[p_i^{-1/2}\|\widehat{\mathbf A}_i\widehat{\mathbf g}_{i,t}-{\mathbf A}_i{\mathbf g}_{i,t}\|_2=O_p\left(p_i^{-1/2}+T^{-1/2}\right),\]
and
\[p^{-1/2}\|\widehat{\mathbf A}\widehat{\mathbf B}\widehat{\mathbf f}_t-{\mathbf A}{\mathbf B}{\mathbf f}_t\|_2=O_p\left(p^{-1/2}k_m^{1/2}\delta_T+p^{-1/2}\right),\]
where $\bar{p}=\max\{p_1,...,p_m\}$, $\underline{p}=\min\{p_1,...,p_m\}$, and $\delta_T=mT^{-1/2}+\underline{p}^{-1}\bar{p}T^{-1/2}+\underline{p}^{-1}$.
\end{theorem}
Theorem \ref{tm3} explores a different aspect of the consistency from that in Theorem \ref{tm2}. The uniform rates in Theorem \ref{tm2} do not imply the convergence pattern as $T\rightarrow\infty$ for fixed but large $p_i$, while the average consistencies in Theorem \ref{tm3} 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{rate:g}) and (\ref{rate:f}) so long as $m=o(\sqrt{T},\underline{p})$. The uniform consistency in (\ref{max:rate2}) 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{tm3} can be relaxed, because we can apply the Bernstein-type inequality in, for example, \cite{merlevede2011}, to reduce the term $m$ in (\ref{DA:DAB}) 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.
\begin{theorem}\label{tm4}
Under the assumptions of Theorem \ref{tm1}, $P(\widehat r_i=r_i)\rightarrow 1$ and $P(\widehat r=r)\rightarrow 1$, as $p_i, T\rightarrow\infty$, where $\widehat r_i$ and $\widehat r$ are defined in (\ref{ri}) and (\ref{est-r}), respectively.
\end{theorem}
By Theorem \ref{tm4}, we can replace $r_i$ and $r$ in Theorems \ref{tm1} and \ref{tm2} by $\widehat r_i$ and $\widehat r$, respectively, and the results therein continue to hold. See \cite{gaotsay2020b} for a similar argument. For the discrepancy measure defined in (\ref{eq:D}) and those considered in Theorem \ref{tm3}, we take $\widehat{\mathbf A}_i$ for example and modify (\ref{eq:D}) 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 \cite{gaotsay2020b} for a similar argument.
\section{Extension to Unit-Root Nonstationary Time Series}\label{sec4}
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{global-f}) 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 \cite{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
\begin{equation}\label{non:cond}
{\mathbf A}_i'{\mathbf A}_i={\mathbf I}_{r_i}\quad\text{and}\quad{\mathbf B}'{\mathbf B}={\mathbf I}_{r}.
\end{equation}
The above identification conditions and Assumption \ref{asm1} 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{asm3}, by Theorems 1 and 2 of \cite{penaponcela2006}, we have
\begin{equation}\label{eigen}
0\leq {\mathbf a}_{i,j}'\widehat\boldsymbol{\Sigma}_{y_i}{\mathbf a}_{i,j}\leq\left\{ \begin{array}{cc}
CT&\text{if},\, E({\mathbf a}_{i,j}'{\mathbf y}_{i,t})=0,\\
CT^2&\text{if}\,\, E({\mathbf a}_{i,j}'{\mathbf y}_{it})\neq 0,
\end{array}\right.
\end{equation}
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{non:cond}) 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 \cite{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
\begin{equation}\label{cov:gs}
\widetilde\boldsymbol{\Sigma}_{g}=\frac{1}{T}\sum_{t=1}^T(\widehat{\mathbf g}_t^*-\bar{\widehat{\mathbf g}}^*)(\widehat{\mathbf g}_t^*-\bar{\widehat{\mathbf g}}^*)',
\end{equation}
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 \cite{gaotsay2020a}. Specifically, for some prescribed integer $\bar{k}>0$, define
\begin{equation}\label{auto:cor}
S_{i,l}(\bar{k})=\sum_{k=1}^{\bar{k}}|\widehat\rho_{i,l}(k)|,\,\, \text{and}\,\,S_j(\bar{k})=\sum_{k=1}^{\bar{k}}|\widehat\rho_j(k)|,
\end{equation}
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})'$.
\begin{assumption}\label{asm5}
(i) The process $\{{\mathbf w}_t,{\mathbf e}_{i,t},{\mathbf u}_t\}$ is $\alpha$-mixing with the mixing coefficient satisfying the condition $\alpha_p(k)\leq \exp(-ck^{\gamma})$ for some constants $c>0$ and $\gamma$ given in Assumption \ref{asm2}; \\
(ii) ${\mathbf w}_t,{\mathbf e}_{i,t}$ and ${\mathbf u}_t$ are sub-exponentially distributed as given in Assumption 3(ii);\\
(iii) There exists a Gaussian process ${\mathbf W}({\mathbf t})=(W_1(t_1),...,W_{r}(t_{r}))'$ such that as $T\rightarrow \infty$,
${\mathbf S}_T^{r}({\mathbf t})\overset{J_1}{\Longrightarrow}{\mathbf W}({\mathbf t})\,\, \text{on}\,\, D_{r_1}[0,1],$
where $\overset{J_1}{\Longrightarrow}$ denotes weak convergence under the Skorokhod $J_1$ topology, and ${\mathbf W}({\bf 1})$ has a positive definite covariance matrix $\boldsymbol{\Omega}=[\sigma_{ij}]$.
\end{assumption}
\begin{assumption}\label{asm6}
There exists a constat $\delta$ with $0<\delta<1/2$ such that
\[\frac{1}{T^{1+\delta}}\sum_{t=1}^T f_{j,t} u_{k,t}=o_p(1)\quad \text{and}\quad \frac{1}{T^{1+\delta}}\sum_{t=1}^Tf_{j,t}e_{i,\iota t}=o_p(1),\]
uniformly, for $1\leq j\leq r$, $1\leq k\leq k_m$, $1\leq i\leq m$, and\, $1\leq \iota\leq p_i$.
\end{assumption}
Assumptions \ref{asm5}(i)-(ii) are used to replace the ones in Assumptions \ref{asm2} and \ref{asm3}(ii) and Assumption \ref{asm5}(iii) is standard for unit-root processes. Assumption \ref{asm6} is not strong and it can be established under the setting of \cite{Stock1987}, where we can assume ${\mathbf y}_{i,t}$ and ${\mathbf g}_t$ have similar structures as the one in (2.4) of \cite{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{tm1} and \ref{tm2} continue to hold with a modified normalization in the proofs; see also \cite{bai2004} for a similar argument regarding the theoretical results with stochastic trend factors and the stationary ones. We do not repeat them here.
\begin{theorem}\label{tm5}
Suppose Assumptions \ref{asm1}, \ref{asm3}(i), \ref{asm4}(i) to (iv), \ref{asm5}, and \ref{asm6} hold, and $r_i$ and $r$ are known, Then, for $1\leq i\leq m$,
we have
\begin{equation}\label{DA:DAB-1}
D(\widehat{\mathbf A}_i,{\mathbf A}_i)=O_p\left(\frac{1}{T^{1-\delta}}\right)\quad \text{and}\quad D(\widehat{\mathbf A}\widehat{\mathbf B},{\mathbf A}{\mathbf B})=O_p\left(\frac{m}{T^{1-\delta}}+\frac{\bar{p}}{T^{1-\delta}\underline{p}}\right).
\end{equation}
Furthermore,
\[p_i^{-1/2}\|\widehat{\mathbf A}_i\widehat{\mathbf g}_{i,t}-{\mathbf A}_i{\mathbf g}_{i,t}\|_2=O_p\left(p_i^{-1/2}+T^{-1/2+\delta}\right),\]
and
\[p^{-1/2}\|\widehat{\mathbf A}\widehat{\mathbf B}\widehat{\mathbf f}_t-{\mathbf A}{\mathbf B}{\mathbf f}_t\|_2=O_p\left(p^{-1/2}k_m^{1/2}\eta_T+p^{-1/2}\right),\]
where $\bar{p}=\max\{p_1,...,p_m\}$, $\underline{p}=\min\{p_1,...,p_m\}$, and $\eta_T=mT^{-1/2+\delta}+\underline{p}^{-1}\bar{p}T^{-1/2+\delta}$.
\end{theorem}
\begin{remark}
(i) The convergence rate of $\widehat{\mathbf A}_i$ is informative. In fact, we can drop the small constant $\delta$ if we make an assumption similar to the distribution results in \cite{Stock1987} in which the sample statistics in Assumption \ref{asm6} are stochastically bounded with $\delta=0$. The resulting rate is $T$, which is in agreement with that of the estimated cointegrating vector in the literature;\\
(ii) The second rate in (\ref{DA:DAB-1}) depends on the number of nodes $m$ and the ratio between the largest and the smallest group dimensions. The constant $m$ can be reduced to a smaller one as we applied the Bonferroni's inequality to the $m$ groups in the proof. If $m$ is reduced and the dimensions of the groups are comparable, then the convergence rate will be similar to that of $\widehat{\mathbf A}_i$.
\end{remark}
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 \cite{gaotsay2020a} and the details are omitted.
\section{Incorporating the Horizontal Partition}\label{sec5}
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{fig0}, 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{sec2} 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{est-c}), 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{tm3},
we have the following corollary regarding the consistency of the loading matrices.
\begin{corollary}
If the conditions of Theorem \ref{tm1} hold for the $j$-th segment of the data. Then, for $1\leq i\leq m$, and $1\leq j\leq n$, we have
\begin{equation}\label{DA:DAB-2}
D(\widehat{\mathbf A}_{i,j},{\mathbf A}_i)=O_p(T_j^{-1/2})\quad \text{and}\quad D(\widehat{\mathbf A}_j^n\widehat{\mathbf B}_j^n,{\mathbf A}{\mathbf B})=O_p(mT_j^{-1/2}+\underline{p}^{-1}\bar{p}T_j^{-1/2}+\underline{p}^{-1}).
\end{equation}
\end{corollary}
In fact, all the theoretical properties in Theorems \ref{tm1} to \ref{tm3} 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 \cite{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 \cite{fan2019}, the theoretical properties
are much more involved as the data are time series with dynamically dependent structures and the theoretical results in \cite{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 \cite{fan2019} is no longer valid under the
current framework. We leave details of this approach for future research.
\section{Numerical Results}\label{sec6}
We use simulation and a real example to assess the performance of the proposed analysis in finite samples.
\subsection{Simulation}
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{eq:D}) 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
\begin{equation}\label{dmeasure}
\bar{D}({\mathbf Q}_1,{\mathbf Q}_2)=\sqrt{1-
\frac{1}{\max{(d_1,d_2)}}\textrm{tr}({\mathbf P}_1{\mathbf P}_2)},
\end{equation}
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{eq:D}).
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.
\noindent {\bf Example 1.} Consider the models in (\ref{m-factor}) and (\ref{global-f}) 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 \texttt{R} to \texttt{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{est-r}) 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{Table1} 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 \cite{BaiNg_Econometrica_2002} and \cite{lamyao2012}.
\begin{table}
\caption{Empirical probabilities $\min_{1\leq i\leq m}EP(\widehat{r}_i=r_i)$,
denoted by $\min_{1\leq i\leq m}EP_i$, and $EP(\widehat r=r)$ of various $(p_i,m,T)$ configurations
for the model of Example 1 with $r_i=5$ and $r=3$, where $p_i$ and $T$ are, respectively, the dimension and the sample size of each group. $500$ iterations are used.}
\label{Table1}
\begin{center}
\setlength{\abovecaptionskip}{0pt}
\setlength{\belowcaptionskip}{3pt}
\begin{tabular}{c|ccccccccc}
\hline
\hline
&&\multicolumn{2}{c}{$T=400$}&&\multicolumn{2}{c}{$T=800$}&&\multicolumn{2}{c}{$T=1200$}\\
$m=10$&&${\min\atop{1\leq i\leq m}}EP_i$&$EP(\widehat r=r)$&&${\min\atop{1\leq i\leq m}}EP_i$&$EP(\widehat r=r)$&&${\min\atop{1\leq i\leq m}}EP_i$&$EP(\widehat r=r)$\\
\cline{1-2}\cline{3-4}\cline{6-7}\cline{9-10}
$p_i=20$&&1&1&&1&1&&1&1\\
$p_i=30$&&1&1&&1&1&&1&1\\
$p_i=40$&&1&1&&1&1&&1&1\\
$p_i=50$&&1&1&&1&1&&1&1\\
\hline
\hline
&&\multicolumn{2}{c}{$T=400$}&&\multicolumn{2}{c}{$T=800$}&&\multicolumn{2}{c}{$T=1200$}\\
$m=20$&&${\min\atop{1\leq i\leq m}}EP_i$&$EP(\widehat r=r)$&&${\min\atop{1\leq i\leq m}}EP_i$&$EP(\widehat r=r)$&&${\min\atop{1\leq i\leq m}}EP_i$&$EP(\widehat r=r)$\\
\cline{1-2}\cline{3-4}\cline{6-7}\cline{9-10}
$p_i=20$&&1&1&&1&1&&1&1\\
$p_i=30$&&1&1&&1&1&&1&1\\
$p_i=40$&&1&1&&1&1&&1&1\\
$p_i=50$&&1&1&&1&1&&1&1\\
\hline
\hline
\end{tabular}
\end{center}
\end{table}
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
\begin{equation}\label{da:am}
\bar{D}_{1}(\widehat{\mathbf A},{\mathbf A})=\frac{1}{m}\sum_{i=1}^m \bar{D}(\widehat{\mathbf A}_i{\mathbf A}_i)\,\,\text{and}\,\,\bar{D}_2(\widehat{\mathbf A},{\mathbf A})=\max_{1\leq i\leq m}\bar{D}(\widehat{\mathbf A}_i,{\mathbf A}_i),
\end{equation}
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{fig1}. From Figure \ref{fig1}(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
\begin{equation}\label{dab:m}
\bar{D}_3(\widehat{\mathbf A}\widehat{\mathbf B},{\mathbf A}{\mathbf B})=\frac{1}{m}\sum_{i=1}^m\bar{D}(\widehat{\mathbf A}_i\widehat{\mathbf B}_i,{\mathbf A}_i{\mathbf B}_i),
\end{equation}
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{fig2}. 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.
\begin{figure}
\begin{center}
\subfigure[]{\includegraphics[width=0.49\textwidth]{ADA-m10.pdf}}
\subfigure[]{\includegraphics[width=0.49\textwidth]{MDA-m10.pdf}}
\caption{(a) Boxplots of the average discrepancy $\bar{D}_1(\widehat{\mathbf A},{\mathbf A})$ when $r_i=5$ and $m=10$ of Example 1; (b) Boxplots of the maximum discrepancy $\bar{D}_2(\widehat{\mathbf A},{\mathbf A})$ when $r_i=5$ and $m=10$ of Example 1, where $\bar{D}_1(\cdot,\cdot)$ and $\bar{D}_2(\cdot,\cdot)$ are defined in (\ref{da:am}). The sample sizes used are $400,\ 800,\ 1200$, and the results are based on $500$ iterations.}\label{fig1}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
{\includegraphics[width=0.6\textwidth]{DAB-m10.pdf}}
\caption{Boxplots of $\bar{D}_3(\widehat{\mathbf A}\widehat{\mathbf B},{\mathbf A}{\mathbf B})$ with $r_i=5$, $r=3$, and $m=10$ of Example 1, where $\bar{D}_3(\widehat{\mathbf A}\widehat{\mathbf B},{\mathbf A}{\mathbf B})$ is defined in (\ref{dab:m}). The sample sizes used are $T=400,\ 800,\ 1200$,
and the number of iterations is $500$.}\label{fig2}
\end{center}
\end{figure}
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
\begin{equation}\label{rmse:f}
\text{RMSE}_1=(\frac{1}{p_imT}\sum_{t=1}^T\|\widehat{\mathbf A}\widehat{\mathbf g}_t-{\mathbf A}{\mathbf g}_t\|_2^2)^{1/2},\,\,\text{RMSE}_2=(\frac{1}{p_imT}\sum_{t=1}^T\|\widehat{\mathbf A}\widehat{\mathbf B}\widehat{\mathbf f}_t-{\mathbf A}{\mathbf B}{\mathbf f}_t\|_2^2)^{1/2},
\end{equation}
which quantify the accuracy in estimating the group factors and the global factors, respectively. Boxplots of the RMSEs are shown in Figure \ref{fig3}. 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{tm3}.
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.
\begin{figure}
\begin{center}
\subfigure[]{\includegraphics[width=0.49\textwidth]{AG-m10.pdf}}
\subfigure[]{\includegraphics[width=0.49\textwidth]{ABF-m10.pdf}}
\caption{(a) Boxplots of RMSE$_1$ with $r_i=5$ and $m=10$ in Example 1; (b) Boxplots of RMSE$_2$ with $r_i=5$ and $m=10$ in Example 1. The sample sizes used are $400,\ 800,\ 1200$, and the results are based on $500$ iterations.}\label{fig3}
\end{center}
\end{figure}
Next, we apply the BIC-type criterion IC1 of Section 2.2 as suggested by \cite{BaiNg_Econometrica_2002} and the ratio-based method in \cite{ahn2013} and \cite{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{fig-r}. From the plots, we see that,
overlooking the hierarchical factor structure, the BIC criterion of \cite{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{Table1}, our hierarchical framework works well to identify both numbers of common factors.
\begin{figure}
\begin{center}
\subfigure[]{\includegraphics[width=0.49\textwidth]{SWR.pdf}}
\subfigure[]{\includegraphics[width=0.49\textwidth]{RTR.pdf}}
\caption{(a) Boxplots of the estimated number of factors using the full data by the IC1 criterion in \cite{BaiNg_Econometrica_2002} when $r_i=5$, $r=3$, and $m=10$ in Example 1; (b) Boxplots of the estimated number of factors using the full data by the ratio-based method in \cite{ahn2013} and \cite{lamyao2012} when $r_i=5$, $r=3$, and $m=10$ in Example 1. The sample sizes used are $400,\ 800,\ 1200$, and the results are based on $500$ iterations.}\label{fig-r}
\end{center}
\end{figure}
\vskip 0.5cm
\noindent{\bf Example 2.} In this example, we consider the models in (\ref{m-factor}) and (\ref{global-f}) 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{auto:cor}) 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{Table2-1} 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{Table2-1}, 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{Table2-1}. Therefore, we do not repeat them here.
\begin{table}
\caption{Empirical probabilities $\min_{1\leq i\leq m}EP(\widehat{r}_i=r_i)$,
denoted by $\min_{1\leq i\leq m}EP_i$, and $EP(\widehat r=r)$ of various $(p_i,T)$ configurations
for the model of Example 2 with $m=10$, $r_i=3$ and $r=5$, where $p_i$ and $T$ are, respectively, the dimension and the sample size of each group. $500$ iterations are used.}
\label{Table2-1}
\begin{center}
\setlength{\abovecaptionskip}{0pt}
\setlength{\belowcaptionskip}{3pt}
\begin{tabular}{c|ccccccccc}
\hline
\hline
&&\multicolumn{2}{c}{$T=400$}&&\multicolumn{2}{c}{$T=800$}&&\multicolumn{2}{c}{$T=1200$}\\
$m=10$&&${\min\atop{1\leq i\leq m}}EP_i$&$EP(\widehat r=r)$&&${\min\atop{1\leq i\leq m}}EP_i$&$EP(\widehat r=r)$&&${\min\atop{1\leq i\leq m}}EP_i$&$EP(\widehat r=r)$\\
\cline{1-2}\cline{3-4}\cline{6-7}\cline{9-10}
$p_i=20$&&0.982&0.980&&1&1&&1&1\\
$p_i=30$&&0.984&0.984&&1&1&&1&1\\
$p_i=40$&&0.992&0.982&&1&1&&1&1\\
$p_i=50$&&0.996&0.994&&1&1&&1&1\\
\hline
\hline
\end{tabular}
\end{center}
\end{table}
\subsection{Real Data Analysis}
In this section, we demonstrate the application of the proposed method using U.S. monthly macroeconomic data with 157 time series and 480 observations.
\noindent{\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 \cite{StockWatson_2002a}. The description of the original data set consisting of 215 series was given in \cite{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 \cite{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 \cite{StockWatson_2002a} and the series of 12 groups are shown in Figure \ref{fig4} with $m=12$.
\begin{figure}
\begin{center}
{\includegraphics[width=0.8\textwidth]{group-data.pdf}}
\caption{Time plots of the differenced series in 12 categories (or groups) from February 1959 to December 1998.
They are (1) Real output and income (Out), (2) Employment and hours (EMP), (3) Retail and manufacturing trade (RTS), (4) Consumption (PCE), (5) Housing starts and sales (HSS), (6) Real inventories and inventory-sales ratios (Inv), (7) Orders and unfilled orders (Ord), (8) Stock prices (Spr), (9) exchange rates (EXR), (10) Interest rates (Int), (11) Money and credit quantity aggregates (Mon), (12) Price indexes (Pri). }\label{fig4}
\end{center}
\end{figure}
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{fig5}. 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{fig6}. 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{fig7}. These factors capture the common cross-sectional dependence between the group factors of series in Figure~\ref{fig5}.
\begin{figure}
\begin{center}
{\includegraphics[width=0.85\textwidth]{group-factors.pdf}}
\caption{Time series plots of group common factors of 12 categories from February 1959 to December 1998. They are (1) Real output and income (Out, $r_1=7$), (2) Employment and hours (EMP, $r_2=9$), (3) Retail and manufacturing trade (RTS, $r_3=3$), (4) Consumption (PCE, $r_4=2$), (5) Housing starts and sales (HSS, $r_5=3$), (6) Real inventories and inventory-sales ratios (Inv, $r_6=4$), (7) Orders and unfilled orders (Ord, $r_7=6$ ), (8) Stock prices (Spr, $r_8=3$), (9) exchange rates (EXR, $r_9=2$,), (10) Interest rates (Int, $r_{10}=7$), (11) Money and credit quantity aggregates (Mon, $r_{11}=3$), (12) Price indexes (Pri, $r_{12}=7$). }\label{fig5}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
{\includegraphics[width=0.9\textwidth]{global-ratio.pdf}}
\caption{(a) The first 26 eigenvalues of $\widehat{\mathbf G}'\widehat{\mathbf G}$; (b) Ratios of consecutive eigenvalues
of $\widehat{\mathbf G}'\widehat{\mathbf G}$ in Example 3.}\label{fig6}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
{\includegraphics[width=0.9\textwidth]{global-factors.pdf}}
\caption{(a) Time series plot of the first global factor process $\widehat f_{1,t}$ of Example 3; (b) Time series plot of the second global factor process $\widehat f_{2,t}$ of Example 3.}\label{fig7}
\end{center}
\end{figure}
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 \cite{BaiNg_Econometrica_2002} to select the number of factors, we find that $\widehat r$ increases as we increase the upper bound $K$ in (\ref{ri}), and it eventually reaches a value greater than $50$, so we adopt the ratio-based method in \cite{ahn2013} and \cite{lamyao2012}, as specified in a similar fashion as (\ref{est-r}). 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{fig8}. 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, \cite{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 \cite{StockWatson_2002a} found that the simple diffusion index forecasts excluding the lags are even better than that with the lagged variables. Table \ref{Table2} 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{df-id}). They are (1) 1glb: including the first global factor in Figure~\ref{fig7}(a) as a predictor; (2) 2glb: including the two global factors in Figure~\ref{fig7}(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 \cite{StockWatson_2002a}. From Table \ref{Table2}, 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 \cite{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.
\begin{figure}
\begin{center}
{\includegraphics[width=0.9\textwidth]{ratio-full.pdf}}
\caption{(a) The first 25 eigenvalues of $\widehat\boldsymbol{\Sigma}_y$; (b) The plot of ratios of consecutive eigenvalues
of $\widehat\boldsymbol{\Sigma}_y$ in Example 3.}\label{fig8}
\end{center}
\end{figure}
\begin{table}[h]
\footnotesize{
\caption{The 1-step to 4-step ahead forecast errors. Five methods employed are according to
factors used in (\ref{df-id}). They are (1) 1glb: including the first global factor in Figure~\ref{fig7}(a) as a predictor; (2) 2glb: including two global factors in Figure~\ref{fig7}(a,b) as predictors; (3)
1glb+spf: including
the frist global factor and the associated group-specific factors ${\mathbf u}_{i,t}$ of $x_{i,t}$; (4)
2glb+spf: including two global factors and the associated group-specific factors ${\mathbf u}_{i,t}$ of $x_{i,t}$; (5) SW: including the 3-dimensional factor process obtained from the full data as that in \cite{StockWatson_2002a}.
Boldface numbers denote the smallest one for a given variable and step. In the table, msmtq$^*$ and
lpnag$*$ denote, respectively, msmtq$/10^7$ and lpnag$/10^4$.}
\label{Table2}
\begin{center}
\setlength{\abovecaptionskip}{0pt}
\setlength{\belowcaptionskip}{3pt}
\begin{tabular}{c|c|cccccccc}
\hline
\hline
Step&Method&ip&gmyxpq&msmtq$^*$ &lpnag$^*$ &punew&gmdx&puxf&pwfsa\\
\hline
1&1glb&0.315&{\bf 268.7}&2.8&2.1&0.050&0.024&0.054&0.212\\
&2glb&0.322&270.9&2.9&2.3&0.050&0.024&0.054&{\bf 0.194}\\
&1glb+spf&{\bf 0.278}&281.2&2.7&{\bf 1.7}&{\bf 0.040}&{\bf 0.020}&{\bf 0.045}&0.198\\
&2glb+spf&0.284&283.8&{\bf 2.6}&{\bf 1.7}&{\bf 0.040}&{\bf 0.020}&{\bf 0.045}&0.200\\
&SW&0.346&279.0&2.7&2.8&0.052&0.024&0.055&0.201\\
\hline
2&1glb&0.305&272.8&2.8&2.2&0.051&0.024&0.054&0.212\\
&2glb&0.285&273.6&2.6&2.2&0.051&0.024&0.056&{\bf 0.203}\\
&1glb+spf&0.301&291.9&2.7&{\bf 1.6}&0.050&{\bf 0.023}&{\bf 0.050}&0.228\\
&2glb+spf&{\bf 0.290}&292.8&{\bf 2.5}&{\bf 1.6}&{\bf 0.049}&{\bf 0.023}&{\bf 0.050}&0.225\\
&SW&0.338&{\bf 272.1}&2.8&2.4&0.058&0.027&0.062&0.218\\
\hline
3&1glb&0.286&272.9&{\bf 2.6}&1.8&0.051&0.024&0.055&0.212\\
&2glb&{\bf 0.269}&262.9&{\bf 2.6}&1.7&0.053&0.024&0.057&{\bf 0.209}\\
&1glb+spf&0.305&268.6&2.7&{\bf 1.4}&0.044&{\bf 0.020}&0.048&0.215\\
&2glb+spf&0.290&{\bf 256.8}&{\bf 2.6}&{\bf 1.4}&{\bf 0.043}&{\bf 0.020}&{\bf 0.046}&0.213\\
&SW&0.333&284.0&2.7&2.5&0.059&0.027&0.065&0.215\\
\hline
4&1glb&0.322&260.2&2.8&2.3&0.051&0.025&0.055&{\bf 0.212}\\
&2glb&{\bf 0.315}&{\bf 259.5}&{\bf 2.7}&2.3&0.052&0.025&0.056&0.215\\
&1glb+spf&0.341&267.2&2.8&{\bf 1.8}&0.046&{\bf 0.021}&{\bf 0.049}&0.240\\
&2glb+spf&0.333&262.6&2.8&{\bf 1.8}&{\bf 0.045}&{\bf 0.021}&0.050&0.236\\
&SW&0.332&263.5&2.8&2.7&0.060&0.027&0.064&0.220\\
\hline
\hline
\end{tabular}
\end{center}}
\end{table}
\section{Concluding Remarks}\label{sec7}
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.