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.
92,416 characters
Modeling High-Dimensional Unit-Root Time Series
\title{Modeling High-Dimensional Unit-Root Time Series}
\author{
Zhaoxing Gao, Department of Mathematics, Lehigh University \\
and Ruey S. Tsay,
Booth School of Business, University of Chicago
}
\maketitle
\begin{abstract}
This paper proposes a new procedure to build factor models for high-dimensional unit-root time
series by postulating that a $p$-dimensional unit-root process is a nonsingular linear transformation of a set of unit-root processes, a set of stationary common factors, which are dynamically dependent, and some idiosyncratic white noise components. For the stationary components,
we assume that the factor process captures the temporal-dependence and the idiosyncratic white noise series explains, jointly with the factors, the cross-sectional dependence.
The estimation of nonsingular linear loading spaces is carried out in two steps.
First, we use an eigenanalysis of a nonnegative definite matrix of the data to separate the unit-root processes from the stationary ones and a modified method to specify the number of unit roots. We then employ another eigenanalysis and a projected principal component analysis to identify the stationary common factors and the white noise series.
We propose a new procedure to specify the number of white noise series and, hence, the
number of stationary common factors, establish asymptotic properties of the proposed method for both fixed and diverging $p$ as the sample size $n$ increases, and
use simulation and a real example to demonstrate
the performance of the proposed method in finite samples. We also compare our method with some commonly used ones in the literature regarding the forecast ability of the extracted factors and find that the proposed method performs well in out-of-sample forecasting of a 508-dimensional PM$_{2.5}$ series in Taiwan.
\end{abstract}
\noindent {\sl Keywords}: Common factor, Cointegration, Eigenanalysis,
Factor model, High-dimensional time series, Unit root.
\newpage
\section{Introduction}
High-dimensional data are common in many scientific fields including biology, business, economics,
and environmental studies. In many applications, the data consist naturally of high-dimensional time series and exhibit characteristics of unit-root nonstationarity. For instance,
the monthly consumer price indexes of European countries tend to exhibit upward trends
associated with inflation.
In theory, the vector autoregressive integrated moving-average (VARIMA) models can be used to
analyze such a high-dimensional time series, but they often encounter the difficulties of
high-dimensional co-integration testing, over-parametrization, and lack of identifiability in real applications. See, for instance, \cite{johansen2002} for the first difficulty and
\cite{TiaoTsay_1989}, \cite{lutkepohl2006}, \cite{Tsay_2014}, and the references therein, for the latter. Therefore, detecting the number of unit roots and
dimension reduction become a necessity in analyzing high-dimensional unit-root time series.
Various methods have been developed in the literature for
dimension reduction or structural specification of multivariate
time series analysis, including the scalar component models of \cite{TiaoTsay_1989}, the LASSO regularization in VAR models by \cite{ShojaieMichailidis_2010} and \cite{SongBickel_2011}, the sparse VAR model via partial spectral coherence in \cite{Davis2012}, and the factor modeling by \cite{BaiNg_Econometrica_2002}, \cite{StockWatson_2002a, StockWatson_2005}, \cite{forni2005} and \cite{lamyao2012}, among others. However, most of the
aforementioned studies focus
on stationary processes. For unit-root time series, cointegration is often used
to account for the common trends and to avoid non-invertibility induced by over-differencing.
See \cite{engle-granger1987}, \cite{johansen1988, johansen1991}, \cite{Tsay_2014}, and the
references therein. But the cointegration rank of a multiple time series is unknown in applications,
and many approaches have been proposed to estimate the unknown rank from data, starting from \cite{engle-granger1987} and the popular likelihood ratio (LR) test in \cite{johansen1988,johansen1991} with a parametric
integrated VAR setting, to \cite{saikkonen2000} and \cite{aznar2002}.
As discussed in \cite{johansen2002}, the conventional co-integration tests may fare poorly when the
dimension of the time seres is high. Yet there exist many applications involving high-dimensional
time series. For example, \cite{engel-etal2015} contemplated the possibility of determining the cointegration rank of a system of seventeen OECD exchange rates. \cite{banerjee-etal2004} emphasized the importance of testing for no cross-sectional cointegration in panel cointegration analysis, and the cross-sectional dimension of modern macroeconomic panel can easily be as large as several hundreds. Therefore, the complexity of the dynamical dependence in high-dimensional unit-root time series requires further investigation, especially extracting dynamic information from such data plays an important role in modeling and forecasting large serially dependent data.
This article provides a new approach to analyze high-dimensional unit-root time series from a factor modeling perspective. Like \cite{zhang-etal2019}, we assume that a $p$-dimensional time series
is a nonsingular linear transformation of some common unit-root processes and a stationary vector process. However, in contrast to \cite{zhang-etal2019}, but in agreement with \cite{bai2004},
we postulate that the number of unit roots is relatively small.
Due to the nature of a unit-root process, we assume the unit-root factors contribute to both the cross-sectional and temporal dependencies of the data, which are different from those in \cite{bai2004}.
To further reduce the dimensionality, we assume the stationary vector process is a nonsingular linear transformation of certain common stationary factors, which are dynamically dependent, and a vector idiosyncratic white noise series. In other words, for the
stationary part, we assume the common factors capture all the non-trivial dynamics of the data, but the cross-sectional dependence may be explained by both the common factors and the idiosyncratic components. This is different from the traditional factor analysis where factors capture most of the cross-sectional dependence, while the idiosyncratic terms may contain some non-trivial temporal dependence; see, for example, \cite{BaiNg_Econometrica_2002}. Under the entertained model,
the common factors explain all the dynamic dependencies of the stationary component of the data and the idiosyncratic white noises may be contemporaneously correlated with each other
so that they also contribute to the cross-sectional dependence between the series.
Therefore, the idiosyncratic white noises of the
proposed model is also different from the orthogonal factor models in, for example, \cite{maxwell1977} that the idiosyncratic terms (or specific factors therein) tend to be combinations of measurement errors and disturbances that are uniquely associated with the individual variables. Our approach generalizes that of \cite{gaotsay2019a, gaotsay2019b} by allowing the number of stationary factors to diverge with the sample size.
To summarize, under the proposed model, a $p$-dimensional time series
is a nonsingular linear transformation of certain unit-root common trends, some stationary common factors which are dynamically dependent, and a white noise idiosyncratic
process. The proposed model is an extension of the work of \cite{zhang-etal2019} and \cite{gaotsay2019b}, and is in line with the framework of \cite{TiaoTsay_1989}
because any finite-order $p$-dimensional
VARMA time series can always
be written as a nonsingular linear transformation of $p$ scalar component models (SCM) via
canonical correlation analyses of constructed vector time series.
See \cite{TiaoTsay_1989} and Section 2.1 below. The only difference is that we
focus on white noise series, which are SCM of order (0,0), and do not consider specifically individual
SCM of order beyond (0,0).
\cite{penaponcela2006} also considered a multiple time series model driven by some common unit-root and some stationary factors when the dimension is finite. This paper also marks an extension of their approach to high-dimensional nonstationary factor modeling with certain model structure.
Although the maximum likelihood method is more efficient than
the principal components method in applying the traditional factor model, it is not feasible for the
proposed model because our factors not only explain the variance of the data, but also capture the dynamic dependencies, and the covariance matrix of the idiosyncratic term may not be
diagonal as in \cite{baili2012}.
Instead, similar to \cite{po1988}, \cite{robinson2002}, \cite{penaponcela2006}, and \cite{zhang-etal2019}, we employ methods based on eigenanalysis. We first estimate the number of unit-root factors (or equivalently the cointegration rank) and extract them from the data by an eigenanalysis of a nonnegative definite matrix, which is a function of the sample covariance and lagged autocovariance matrices of the data. The nonnegative definite matrix used in this paper is different from the sample covariance used in \cite{bai2004} and the fixed lag sample covariance or autocovariances used in \cite{penaponcela2006}, because we adopt a combination of the covariance and some lagged autocovariances together to capture simultaneously the cross-sectional and temporal dependence of the data. In addition, we propose to use an average of the absolute autocorrelations of the transformed components to identify the number of unit-root series. The absolute value can avoid the impact of sign changes in the autocorrelations introduced by the stationary part embedded in the individual unit-root component. Limited simulation studies suggest that, although the performance of the proposed method is comparable with the one in \cite{zhang-etal2019} when the sample size is sufficiently large, the use of absolute autocorrelations can improve the accuracy in estimating the number of unit-root processes, especially when the sample size is small. An accurate specification of the number of unit-root factors can provide more accurate information for the second eigenanalysis, from which the number of stationary
common factors is identified. Specifically, to estimate the number of stationary common factors in
the second eigenanalysis, we apply the method of \cite{gaotsay2019b} to the transformed data which are orthogonal to the unit-root components. Since eigenanalysis
is mainly based on spectral decomposition of a nonnegative definite matrix, the resulting
components are arranged according to the amount of variabilities explained. The ordering
of the components thus provides no information concerning the temporal dependence in each
principal component. Consequently, to specify the number of white noises and, hence, the number of
stationary common factors, we propose to reorder the transformed components of
eigenanalysis according to their $p$-values (in ascending order) of the Ljung-Box statistic in
testing their serial correlations. Limited experience shows that this re-ordering
procedure is helpful in detecting the number of white noise series when the dimension is large
and the sample size is small.
Under the proposed framework, the dimension of the idiosyncratic white noise may go to infinity
and some largest eigenvalues of the covariance matrix of the white noise may also diverge.
We refer to the latter case as prominent noise effect and apply the projected
principal component analysis of \cite{gaotsay2019b} to mitigate the effect of such prominent noises
in estimating the stationary common factors. Consequently, our proposed method
can successfully separate the nonstationary unit-root processes, the stationary common factors,
and the idiosyncratic white noise components.
In estimating the number of unit-root series (or equivalently the cointegration rank), we could allow the dimension $p$ to grow as fast as the sample size $n$. This relaxes the constraint of
\cite{zhang-etal2019} and the error-correction factor models of \cite{tu-etal2019} that $p$ can
grow at most with $\sqrt{n}$. Asymptotic properties of the proposed method are established for both fixed $p$ and diverging $p$ as the sample size $n$ tends to infinity.
The rest of the paper is organized as follows. We introduce the proposed model and
estimation methodology in Section 2.
In Section 3, we study the theoretical properties of the proposed model and its associated
estimates.
Section 4 illustrates the performance of the proposed model using both
simulated and real data sets, including analysis of a 508-dimensional time series of
PM$_{2.5}$ measurements in Taiwan with 744 observations.
Section 5 provides some discussions and concluding remarks.
All technical proofs are relegated to an Appendix.
Throughout the article,
we use the following notation: $||{\mathbf u}||_2 = (\sum_{i=1}^{p} u_i^2)^{1/2} $
is the Euclidean norm of a $p$-dimensional vector
${\mathbf u}=(u_1,..., u_p)'$, $\|{\mathbf u}\|_\infty=\max_i |u_i|$, and ${\mathbf I}_k$ denotes the $k\times k$ identity matrix. For a matrix ${\mathbf H}=[h_{ij}]$, $\|{\mathbf H}
\|_2=\sqrt{\lambda_{\max} ({\mathbf H}' {\mathbf H} ) }$ is the operator norm, where
$\lambda_{\max} (\cdot) $ denotes 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.
Finally, we use the notation $a\asymp b$ to denote $a=O(b)$ and $b=O(a)$.
\section{The Proposed Methodology}
\subsection{The Setting and Framework}
Let ${\mathbf y}_t=(y_{1t},...,y_{pt})'$ be a $p$-dimensional $I(1)$ time series process.
We assume ${\mathbf y}_t$ is observable and admits the following latent structure
\begin{equation}\label{int:eq}
{\mathbf y}_t={\mathbf L}\left[\begin{array}{c}
{\mathbf f}_{1t}\\
{\mathbf f}_{2t}\\
\mbox{\boldmath$\varepsilon$}_t
\end{array}\right]=[{\mathbf L}_1,{\mathbf L}_2,{\mathbf L}_3]\left[\begin{array}{c}
{\mathbf f}_{1t}\\
{\mathbf f}_{2t}\\
\mbox{\boldmath$\varepsilon$}_t
\end{array}\right]={\mathbf L}_1{\mathbf f}_{1t}+{\mathbf L}_2{\mathbf f}_{2t}+{\mathbf L}_3\mbox{\boldmath$\varepsilon$}_t,
\end{equation}
where ${\mathbf L}\in R^{p\times p}$ is a full rank loading matrix, ${\mathbf f}_{1t}=(f_{1,1t},\ldots,f_{1,r_1t})'$
is an $r_1$-dimensional $I(1)$ process, ${\mathbf f}_{2t}=(f_{2,1t},\ldots,f_{2,r_2t})'$ is an $r_2$-dimensional stationary process, and $\mbox{\boldmath$\varepsilon$}_t=({{\varepsilon}}_{1t},\ldots,{{\varepsilon}}_{vt})$ is a $v$-dimensional white noise series with $v=p-r$, $r=r_1+r_2$, and $r_i \geq 0$. For meaningful dimension reduction, we assume $r_1$ is a relatively small and fixed integer as that in \cite{bai2004}, and $r_2$ can be either fixed or slowly growing with the dimension $p$, which extends the results in \cite{gaotsay2019b}.
In addition, we also assume that ${\mathbf f}_{2t}$ and $\mbox{\boldmath$\varepsilon$}_t$ are
independent of each other with $\textnormal{Cov}({\mathbf f}_{2t})={\mathbf I}_{r_2}$ and $\textnormal{Cov}(\mbox{\boldmath$\varepsilon$}_t)={\mathbf I}_v$,
and no linear combination of ${\mathbf f}_{1t}$ is a stationary process and no linear combination
of ${\mathbf f}_{2t}$ is a white noise. In theory, \textnormal{Cov}(${\mathbf f}_{1t}$) is time-varying
because ${\mathbf f}_{1t}$ consists of unit-root processes, but its sample version may assume
an identity matrix when the sample size is given and the processes are assumed to
start at $t$ = 0 with fixed starting values.
The decomposition of (\ref{int:eq}) is general and in line with the framework
of \cite{TiaoTsay_1989}.
Under the scalar component models of \cite{TiaoTsay_1989}, any finite-order
$p$-dimensional VARIMA($p,d,q$) process consists of $p$ scalar component models (SCM)
of orders $(p_i,q_i)$ with $p_i+d_i \leq p+d$ and
$q_i \leq q$. The SCM and the associated transformation matrix can be obtained via
canonical correlation analyses between constructed random vectors of ${\mathbf y}_t$
and its lagged variables. Model (\ref{int:eq}) can be considered as consisting of
three classes of SCM, namely, the unit-root processes, the processes of SCM(0,0), and
the process of stationary SCM($p_i,q_i$) with $p_i+q_i > 0$. Specifically,
let ${\mathbf V}'={\mathbf L}^{-1}$ and ${\mathbf V}=({\mathbf V}_1,{\mathbf V}_2,{\mathbf V}_3)$ with ${\mathbf V}_1\in R^{p\times r_1}$, ${\mathbf V}_2\in R^{p\times r_2}$ and ${\mathbf V}_3\in R^{p\times v}$. Model (\ref{int:eq}) is a transformation model employed in \cite{TiaoTsay_1989}.
Another way to see the generality of Model (\ref{int:eq}) is to employ the canonical correlation
analysis of \cite{BoxTiao_1977} involving a VARIMA process ${\mathbf y}_t$ and its
1-step ahead prediction at time index $t-1$. Under such an analysis, $({\mathbf f}_{1t}', {\mathbf f}_{2t}', \mbox{\boldmath$\varepsilon$}_t')'$ are the canonical variates whose correlations with the past values are in descending order, and the correlation that is close to unity corresponds to a unit-root nonstationary series, and those close to zero are close to a white noise; see also \cite{penaponcela2006}. As a result, we consider the top canonical variates ${\mathbf f}_{1t}$ a unit-root process, the bottom $\mbox{\boldmath$\varepsilon$}_t$ a white noise, and the ones in the middle are dynamically dependent stationary processes. Consequently, the columns of $({\mathbf V}_2,{\mathbf V}_3)$ are the cointegrating vectors of ${\mathbf y}_t$.
For more details, we refer interested readers to \cite{BoxTiao_1977} and \cite{TiaoTsay_1989} for general discussions.
To illustrate the identification issue of Model (\ref{int:eq}) and to provide a concrete analysis,
we let ${\mathbf g}_{2t}=({\mathbf f}_{2t}',\mbox{\boldmath$\varepsilon$}_t')'$ and ${\mathbf G}_2=[{\mathbf L}_2,{\mathbf L}_3]$, and rewrite
Model (\ref{int:eq}) as
\begin{equation}\label{2:eq}
{\mathbf y}_t=[{\mathbf L}_1,{\mathbf G}_2]\left[\begin{array}{c}
{\mathbf f}_{1t}\\
{\mathbf g}_{2t}
\end{array}\right],
\end{equation}
where ${\mathbf g}_{2t}$ is a $(p-r_1)$-dimensional stationary process.
Note that Model (\ref{2:eq}) is not uniquely defined, as $[{\mathbf L}_1,{\mathbf G}_2]$ and $({\mathbf f}_{1t}',{\mathbf g}_{2t}')'$ can be replaced by $[{\mathbf L}_1,{\mathbf G}_2]{\mathbf H}^{-1}$ and ${\mathbf H}({\mathbf f}_{1t}',{\mathbf g}_{2t}')'$, respectively, for any invertible ${\mathbf H}$ with the form
\begin{equation}\label{H}
{\mathbf H}=\left[\begin{array}{cc}
{\mathbf H}_{11}&{\mathbf H}_{12}\\
\bf 0&{\mathbf H}_{22}
\end{array}\right],
\end{equation}
where ${\mathbf H}_{11}$ and ${\mathbf H}_{22}$ are square matrices of sizes $(p-r_1)$ and $r_1$, respectively. In other words, the nonstationary components can include any linear combinations of the stationary ones. However, for any nonorthogonal invertible matrix $[{\mathbf L}_1,{\mathbf G}_2]$, we always have the decomposition $[{\mathbf L}_1,{\mathbf G}_2]={\mathbf A}{\mathbf U}$, where ${\mathbf A}$ is orthonormal and ${\mathbf U}$ is upper-triangular, and we may replace $[{\mathbf L}_1,{\mathbf G}_2]$ and $({\mathbf f}_{1t}',{\mathbf g}_{2t}')'$ by ${\mathbf A}$ and ${\mathbf U}({\mathbf f}_{1t}',{\mathbf g}_{2t}')'$ without altering the structure of the model. Let ${\mathbf A}=[{\mathbf A}_1,{\mathbf A}_2]$ and $({\mathbf x}_{1t}',{\mathbf x}_{2t}')'={\mathbf U}({\mathbf f}_{1t}',{\mathbf g}_{2t}')'$. There is no loss of generality in assuming that
\begin{equation}\label{w:eq}
{\mathbf y}_t={\mathbf A}\left[\begin{array}{c}
{\mathbf x}_{1t}\\
{\mathbf x}_{2t}
\end{array}\right]=[{\mathbf A}_1,{\mathbf A}_2]\left[\begin{array}{c}
{\mathbf x}_{1t}\\
{\mathbf x}_{2t}
\end{array}\right],
\end{equation}
where ${\mathbf A}$ is an orthonormal matrix, ${\mathbf x}_{1t}$ is an $r_1$-dimensional $I(1)$ process, and ${\mathbf x}_{2t}$ is a $(p-r_1)$-dimensional stationary process. Therefore, ${\mathbf x}_{1t}={\mathbf A}_1'{\mathbf y}_t$ and ${\mathbf x}_{2t}={\mathbf A}_2'{\mathbf y}_t$. For any ${\mathbf H}$ in the form of (\ref{H}) to be orthonormal, we can show that ${\mathbf H}$ is a block-orthonormal matrix. Thus, Model (\ref{w:eq}) is still not identifiable and ${\mathbf A}_1$ and ${\mathbf A}_2$ cannot be uniquely defined. However, the linear spaces spanned by the columns of ${\mathbf A}_1$ and ${\mathbf A}_2$, denoted by $\mathcal{M}({\mathbf A}_1)$ and $\mathcal{M}({{\mathbf A}_2})$, can be uniquely defined.
To proceed with the proposed dimension reduction procedure,
noting that $({\mathbf x}_{1t}',{\mathbf x}_{2t}')'={\mathbf U}({\mathbf f}_{1t}',{\mathbf g}_{2t}')'$ symbolically for an upper triangular matrix ${\mathbf U}$, we further assume that
\begin{equation}\label{st:eq}
{\mathbf x}_{2t}={\mathbf U}_{22}\left[\begin{array}{c}
{\mathbf f}_{2t}\\
\mbox{\boldmath$\varepsilon$}_t
\end{array}\right]={\mathbf U}_{22,1}{\mathbf f}_{2t}+{\mathbf U}_{22,2}\mbox{\boldmath$\varepsilon$}_t,
\end{equation}
where ${\mathbf U}_{22}=[{\mathbf U}_{22,1},{\mathbf U}_{22,2}]$ is the lower diagonal block of ${\mathbf U}$ in the form of (\ref{H}). Given Models (\ref{w:eq}) and (\ref{st:eq}), we estimate $r_1$, $r_2$, the linear spaces $\mathcal{M}({\mathbf A}_1)$, $\mathcal{M}({\mathbf A}_2)$, and $\mathcal{M}({\mathbf A}_2{\mathbf U}_{22,1})$ as well as
recover the processes ${\mathbf x}_{1t}$ and ${\mathbf f}_{2t}$.
\subsection{Estimation Methods}
To introduce estimation, we assume first that $r_1$ and $r_2$ are known and estimate ${\mathbf A}_1$, ${\mathbf A}_2$, and ${\mathbf U}_{22,1}$, or equivalently the linear spaces spanned by their columns. The estimation
of $r_1$ and $r_2$ is given in Section 2.3 below. We start with the case that $p$ is finite.
For $k\geq 0$, we define the lag-$k$ sample covariance matrix of ${\mathbf y}_t$ as
\begin{equation}\label{sigy}
\widehat\boldsymbol{\Sigma}_y(k)=\frac{1}{n}\sum_{t=k+1}^n({\mathbf y}_{t}-\bar{{\mathbf y}})({\mathbf y}_{t-k}-\bar{{\mathbf y}})',\quad\bar{{\mathbf y}}=\frac{1}{n}\sum_{t=1}^n{\mathbf y}_t.
\end{equation}
For any ${\mathbf a}_1\in\mathcal{M}({\mathbf A}_1)$ and ${\mathbf a}_2\in\mathcal{M}({\mathbf A}_2)$, ${\mathbf a}_1'\widehat\boldsymbol{\Sigma}_y(k){\mathbf a}_1$ is the lag-$k$ sample autocovariance
of the $I(1)$ process ${\mathbf a}_1'{\mathbf y}_t$, and ${\mathbf a}_2'\widehat\boldsymbol{\Sigma}_y(k){\mathbf a}_2$ is that of the weakly stationary univariate time series ${\mathbf a}_2'{\mathbf y}_t$. When $p$ is finite, it is not hard to see that ${\mathbf a}_2'\widehat\boldsymbol{\Sigma}_y(k){\mathbf a}_2$ converges to a finite constant almost surely under some mild conditions, and with probability tending to one that
\begin{equation}\label{auto:order}
0<{\mathbf a}_1'\widehat\boldsymbol{\Sigma}_y(k){\mathbf a}_1\leq \left\{\begin{array}{ll} Cn & \mbox{if $E({\mathbf a}_1'{\mathbf y}_t) = 0$}, \\ Cn^2 & \mbox{if $E({\mathbf a}_1'{\mathbf y}_t) \neq 0$,}\end{array}\right.
\end{equation}
for some constant $0< C<\infty$. See Theorems 1 and 2 of \cite{penaponcela2006}.
Hence, the $r_1$ directions in the space $\mathcal{M}({\mathbf A}_1)$ make ${\mathbf a}_1'\widehat\boldsymbol{\Sigma}_y(k){\mathbf a}_1$ as large as possible for all $k\geq 0$.
\subsubsection{Estimation of Unit-Root Processes}
Similar to several works in the literature, e.g., \cite{lamyao2012}, we combine the information over different lags of ${\mathbf y}_t$ and define
\begin{equation}\label{wy}
\widehat{\mathbf M}_1=\sum_{k=0}^{k_0}\widehat\boldsymbol{\Sigma}_y(k)\widehat\boldsymbol{\Sigma}_y(k)',
\end{equation}
where $k_0\geq 1$ is a prespecified fixed integer. We use the product $\widehat\boldsymbol{\Sigma}_y(k)\widehat\boldsymbol{\Sigma}_y(k)'$ instead of $\widehat\boldsymbol{\Sigma}_y(k)$ to ensure that each term in the sum of (\ref{wy}) is nonnegative definite, and there is no information cancellation over different lags. Unlike \cite{bai2004} and \cite{penaponcela2006}, $\widehat{\mathbf M}_1$ of (\ref{wy}) combines the sample covariance and lag autocovariances together to focus on both the cross-sectional and the dynamic dependence of the data simultaneously. Limited experience suggests that a relatively small $k_0$ is sufficient in providing useful information concerning the model structure of ${\mathbf y}_t$, because, for the stationary component ${\mathbf x}_{2t}$, cross-correlation matrices decay to zero exponentially as $k$ increases. Also, the choice of $k_0$ seems to be not sensitive. See, for instance, the simulation results in Section 4. It can be shown that the $r_1$ largest eigenvalues of $\widehat{\mathbf M}_1$ are at least of order $n^2$, while the other $(p-r_1)$ eigenvalues are $O_p(1)$. Hence, $\mathcal{M}({\mathbf A}_1)$ can be estimated by the linear space spanned by the $r_1$ eigenvectors of $\widehat{\mathbf M}_1$ corresponding to the $r_1$ largest eigenvalues, and $\mathcal{M}({\mathbf A}_2)$ can be estimated by that spanned by the $(p-r_1)$ eigenvectors of $\widehat{\mathbf M}_1$ corresponding to the $(p-r_1)$ smallest eigenvalues.
Let $(\widehat{\mathbf a}_{1,1},...,\widehat{\mathbf a}_{1,r_1},\widehat{\mathbf a}_{2,1},...,\widehat{\mathbf a}_{2,p-r_1})$ be the orthonormal eigenvectors of $\widehat{\mathbf M}_1$ corresponding to the eigenvalues arranged in descending order. Define $\widehat{\mathbf A}_1=(\widehat{\mathbf a}_{1,1},...,\widehat{\mathbf a}_{1,r_1})$ and $\widehat{\mathbf A}_2=(\widehat{\mathbf a}_{2,1},...,\widehat{\mathbf a}_{2,p-r_1})$, the estimated $\widehat{\mathbf x}_{1t}$ and $\widehat{\mathbf x}_{2t}$ are given by
\begin{equation}\label{xhat}
\widehat{\mathbf x}_{1t}=\widehat{\mathbf A}_1'{\mathbf y}_t\quad\text{and}\quad\widehat{\mathbf x}_{2t}=\widehat{\mathbf A}_2'{\mathbf y}_t.
\end{equation}
Then, $\mathcal{M}(\widehat{{\mathbf A}}_1)$ and $\mathcal{M}(\widehat{{\mathbf A}}_2)$, the linear spaces spanned by the eigenvectors of $\widehat{\mathbf M}_1$, are consistent estimators for $\mathcal{M}({\mathbf A}_1)$ and $\mathcal{M}({\mathbf A}_2)$, respectively. See Theorem 1 below for details.
When $p$ is diverging, it is reasonable to consider strengths of the factors ${\mathbf x}_{1t}$ and ${\mathbf x}_{2t}$, and the strengths from the columns of ${\mathbf L}_1$ and ${\mathbf G}_2$. See the discussion in \cite{gaotsay2019b} for details. For simplicity, we introduce a parameter $\delta\in[0,1)$ such that the nonzero singular values of ${\mathbf L}_1$ and ${\mathbf L}_2$, and a few largest singular values
of ${\mathbf L}_3$ are of order $p^{(1-\delta)/2}$. It is not hard to see that the order of ${\mathbf a}_2'\widehat\boldsymbol{\Sigma}_y(k){\mathbf a}_2$ is $O_p(p^{1-\delta})$ under some mild conditions, and with probability tending to one that $0<{\mathbf a}_1'\widehat\boldsymbol{\Sigma}_y(k){\mathbf a}_1\leq Cp^{1-\delta}n$ or $Cp^{1-\delta}n^2$ for some $0<C<\infty$, depending on whether $E({\mathbf a}_1'{\mathbf y}_t)=0$ or not. Therefore, $\mathcal{M}({\mathbf A}_1)$ can still be estimated by the linear space spanned by the $r_1$ eigenvectors of $\widehat{\mathbf M}_1$ corresponding to the $r_1$ largest eigenvalues, and $\mathcal{M}({\mathbf A}_2)$ can be estimated by that spanned by the remaining $(p-r_1)$ eigenvectors of $\widehat{\mathbf M}_1$ corresponding to the $(p-r_1)$ smallest eigenvalues. The estimators for ${\mathbf x}_{1t}$ and ${\mathbf x}_{2t}$ are the same as those in (\ref{xhat}).
\subsubsection{Estimation of Stationary Common Factors}
Turn to the estimation of ${\mathbf U}_{22,1}$ and ${\mathbf f}_{2t}$. Because ${\mathbf U}_{22,1}$ is not uniquely defined and we can replace $({\mathbf A}_{2},{\mathbf U}_{22,1})$ by $({\mathbf A}_{2}{\mathbf H}_1',{\mathbf H}_1{\mathbf U}_{22,1})$ or replace $({\mathbf U}_{22,1},{\mathbf f}_{2t})$ by $({\mathbf U}_{22,1}{\mathbf H}_2',{\mathbf H}_2{\mathbf f}_{2t})$ for any orthonormal matrices ${\mathbf H}_1\in R^{(p-r_1)\times (p-r_1)}$ and ${\mathbf H}_2\in R^{r_2\times r_2}$ without altering the relations. Therefore, only $\mathcal{M}({\mathbf A}_{2}{\mathbf U}_{22,1})$ can be estimated.
We decompose ${\mathbf U}_{22,1}$ and ${\mathbf U}_{22,2}$ as ${\mathbf U}_{22,1}={\mathbf U}_1{\mathbf Q}_1$ and ${\mathbf U}_{22,2}={\mathbf U}_2{\mathbf Q}_2$, where ${\mathbf U}_{1}$ and ${\mathbf U}_2$ are two half orthonormal matrices in the sense that ${\mathbf U}_1'{\mathbf U}_1={\mathbf I}_{r_2}$ and ${\mathbf U}_2'{\mathbf U}_2={\mathbf I}_{p-v}$. Then, Model (\ref{st:eq}) can be written as
\begin{equation}\label{re:st}
{\mathbf x}_{2t}={\mathbf U}_1{\mathbf z}_{2t}+{\mathbf U}_2{\mathbf e}_t,
\end{equation}
where ${\mathbf z}_{2t}={\mathbf Q}_1{\mathbf f}_{2t}$ and ${\mathbf e}_t={\mathbf Q}_2\mbox{\boldmath$\varepsilon$}_t$. Equivalently, we estimate
$\mathcal{M}({\mathbf A}_2{\mathbf U}_1)$ and recover ${\mathbf z}_{2t}$.
First, we can apply the method of \cite{gaotsay2019b} to estimate the linear space $\mathcal{M}({\mathbf U}_1)$ and recover ${\mathbf z}_{2t}$. Specifically, letting $\widehat{\mathbf x}_{2t}=\widehat{\mathbf A}_2'{\mathbf y}_t$, we
estimate $\mathcal{M}({\mathbf U}_1)$ using the transformed data $\widehat{\mathbf x}_{2t}$ because
$\mathcal{M}(\widehat{\mathbf A}_2)$ is a consistent estimator for $\mathcal{M}({\mathbf A}_2)$. Define
\begin{equation}\label{M2}
\widehat{\mathbf M}_2=\sum_{j=1}^{j_0}\widetilde\boldsymbol{\Sigma}_2(j)\widetilde\boldsymbol{\Sigma}_2(j)',
\end{equation}
where $\widetilde\boldsymbol{\Sigma}_2(j)$ is the lag-$j$ sample autocovariance matrix of $\widehat{\mathbf x}_{2t}$
and $j_0$ is a prespecified and fixed positive integer. Assume that ${\mathbf V}_1$ and ${\mathbf V}_2$ are the orthogonal complements of ${\mathbf U}_1$ and ${\mathbf U}_2$, respectively. From (\ref{re:st}), we see that ${\mathbf V}_1'{\mathbf x}_{2t}={\mathbf V}_1'{\mathbf U}_2{\mathbf e}_t$ and ${\mathbf V}_2'{\mathbf x}_{2t}={\mathbf V}_2'{\mathbf U}_1{\mathbf z}_{2t}$ are uncorrelated with each other, and therefore, ${\mathbf V}_2$ consists of the eigenvectors corresponding to the $r_2$ zero eigenvalues of ${\mathbf S}:=\boldsymbol{\Sigma}_2{\mathbf V}_1{\mathbf V}_1'\boldsymbol{\Sigma}_2$, where $\boldsymbol{\Sigma}_2=\textnormal{Var}({\mathbf x}_{2t})$. This leads to the following projected PCA method. Let the columns of $[\widehat{\mathbf U}_1,\widehat{\mathbf V}_1]$ be the orthonormal eigenvectors of $\widehat{\mathbf M}_2$ corresponding to the eigenvalues arranged in descending order, where $\widehat{\mathbf U}_{1}$ contains the first $r_2$ columns and $\widehat{\mathbf V}_{1}$ contains the
remaining $(p-v)$ columns, and $\widetilde\boldsymbol{\Sigma}_2$ be the sample covariance of $\widehat{\mathbf x}_{2t}$. Then the estimation of ${\mathbf V}_2$ is based on the eigenanalysis of
\begin{equation}\label{Spca}
\widehat{\mathbf S}=\widetilde\boldsymbol{\Sigma}_2\widehat{\mathbf V}_1\widehat{\mathbf V}_1'\widetilde\boldsymbol{\Sigma}_2,
\end{equation}
which is a projected PCA and $\widehat{\mathbf S}$ is an estimator of ${\mathbf S}$. That is, we project the data $\widehat{\mathbf x}_{2t}$ onto the direction of $\widehat{\mathbf V}_1$, then perform the PCA between $\widehat{\mathbf x}_{2t}$ and the projected coordinates.
When $p$ is finite, assume $\widehat{\mathbf V}_2$ contains the last $r_2$ eigenvectors corresponding the the smallest $r_2$ eigenvalues of $\widehat{\mathbf S}$, then from the relation ${\mathbf V}_2'{\mathbf x}_{2t}={\mathbf V}_2'{\mathbf U}_1{\mathbf z}_{2t}$, the recovered ${\mathbf z}_{2t}$ process is
\begin{equation}\label{z2hat}
\widehat{\mathbf z}_{2t}=(\widehat{\mathbf V}_2'\widehat{\mathbf U}_1)^{-1}\widehat{\mathbf V}_2'\widehat{\mathbf x}_{2t}.
\end{equation}
When $p$ is large, similar as that in \cite{gaotsay2019b}, we assume the largest $K$ eigenvalues of the covariance of ${\mathbf U}_2{\mathbf e}_t$ are diverging, and therefore, the largest $K$
eigenvalues of $\widehat{\mathbf S}$ are also diverging. Suppose $\widehat{\mathbf V}_2^*$ contains $(p-r_1-K)$ eigenvectors of $\widehat{\mathbf S}$ corresponding to its $(p-r_1-K)$ smallest eigenvalues, we choose $\widehat{\mathbf V}_2$ as $\widehat{\mathbf V}_2=\widehat{\mathbf V}_2^*\widehat{\mathbf R}$, where $\widehat{\mathbf R}=[\widehat{\mathbf r}_1,...,\widehat{\mathbf r}_{r_2}]\in R^{(p-r_1-K)\times r_2}$ with $\widehat{\mathbf r}_i$ being the eigenvector associated with the $i$-th largest eigenvalues of $\widehat{\mathbf V}_2^*{'}\widehat{\mathbf U}_1\widehat{\mathbf U}_1'\widehat{\mathbf V}_{2}^*$. This choice of estimator guarantees that the matrix $(\widehat{\mathbf V}_2'\widehat{\mathbf U}_1)^{-1}$ behaves well in recovering the factor $\widehat{\mathbf z}_{2t}$. The resulting estimator $\widehat{\mathbf z}_{2t}$ is the same as that in (\ref{z2hat}).
\subsubsection{Prediction}
With $\widehat{\mathbf A}_1$, $\widehat{\mathbf A}_2$, $\widehat{\mathbf U}_1$, and the estimated nonstationary factor process $\widehat{\mathbf x}_{1t}$ and the stationary one $\widehat{\mathbf z}_{2t}$, we compute the $h$-step ahead prediction of the ${\mathbf y}_t$ series using the formula $\widehat{\mathbf y}_{n+h}=\widehat{\mathbf A}_1\widehat{\mathbf x}_{1,n+h}+\widehat{\mathbf A}_2\widehat{\mathbf U}_1\widehat{\mathbf z}_{2,n+h}$, where $\widehat{\mathbf x}_{1,n+h}$ and $\widehat{\mathbf z}_{2,n+h}$ are $h$-step ahead forecasts for ${\mathbf x}_{1t}$ and ${\mathbf z}_{2t}$ based on the estimated past values $\{\widehat{\mathbf x}_{11},...,\widehat{\mathbf x}_{1n}\}$ and $\{\widehat{\mathbf z}_{21},...,\widehat{\mathbf z}_{2n}\}$, respectively. This can be done, for example, by fitting a VAR model to $\{\widehat{\mathbf x}_{11},...,\widehat{\mathbf x}_{1n}\}$ and $\{\widehat{\mathbf z}_{21},...,\widehat{\mathbf z}_{2n}\}$, respectively. In addition, we may also adopt the factor-augmented error correction model in \cite{banerjee-etal2014} to do forecasting using the nonstationary factors $\widehat{\mathbf x}_{1t}$ if we are interested in some particular components of ${\mathbf y}_t$.
\subsection{Determination of $r_1$ and $r_2$}
Turn to the estimation of $r_1$ and $r_2$, which are unknown in practice.
We begin with the estimation of the number of unit roots $r_1$.
Note that the components of $\widehat{\mathbf x}_t=\widehat{\mathbf A}'{\mathbf y}_t=(\widehat x_{1t},\ldots,\widehat x_{pt})'$, defined in (\ref{xhat}), are arranged according to the descending order of eigenvalues of $\widehat{\mathbf M}_1$. Therefore, the order of the components reflects inversely the closeness to stationarity of the component series, with $\{\widehat x_{pt}\}$ being most likely stationary, and $\{\widehat x_{1t}\}$ being most likely an $I(1)$ process. Based on this observation, we consider a modified
method to estimate the number of nonstationary components $r_1$.
\cite{zhang-etal2019} use the average of sample autocorrelations of each transformed
component to determine its stationarity. Our proposed method makes two important modifications,
because (a) limited simulation studies show that some transformed stationary components also have large autocorrelations when the dimension is high, especially at the lower-order lags,
and (b) the stationary part included in a unit-root process may change the sign
of its autocorrelations. Our modified method is given below.
Let $k_1=1$ and $k_j=k_1+(j-1)l$ for some constant $l\geq 1$. Define
\begin{equation}\label{sm}
S_i(l,m)=\sum_{j=1}^m|\widehat\rho_i(k_j)|,
\end{equation}
where $\widehat\rho_i(j)$ is the lag-$j$ sample autocorrelation function (ACF) of $\widehat x_{it}$.
If $\widehat x_{it}$
is stationary, then under some mild conditions, $\widehat{\rho}_i(k)$ decays exponentially as $k$ increases. Therefore, $\lim_{m\rightarrow\infty} S_i(l,m)<\infty$ in probability. On the other hand,
if $\widehat x_{it}$ is unit-root nonstationary, $\widehat\rho_i(k)\rightarrow 1$ in probability for any fixed $k$ as $n\rightarrow\infty$. Therefore, $\lim_{m\rightarrow\infty} S_i(l,m)=\infty$ as $n\rightarrow
\infty$.
We use a gap of size $l\geq 1$ in (\ref{sm}) to reduce the effect of high-dimensionality on the
lower-order sample ACFs of transformed stationary components, and employ the absolute ACFs
to avoid the effect of sign changes in ACFs due to the existence of stationary part embedded in the unit-root components. Using the statistic $S_i(l,m)$, we propose the following thresholding
procedure to estimate $r_1$.
We start with $i=1$, if the average of absolute sample autocorrelations
$S_i(l,m)/m\geq c_0$ for some constant $0<c_0<1$, the $\widehat{x}_{it}$ series has a unit root
and increase $i$ by 1 to repeat the detecting process. This detecting process is continued
until $S_i(l,m)/m< c_0$ or $i=p$. If $S_i(l,m) \geq c_0$ for all $i$, then $\widehat{r}_1=p$; otherwise,
an estimator of $r_1$ is $\widehat r_1=i-1$. In our numerical experiments of Section 4,
we set $c_0=0.3$, $l=3$ and $m=10$, and the estimator $\widehat r_1$ performs very well. Consequently, $\widehat{\mathbf A}_1$ and $\widehat{\mathbf A}_2$ in (\ref{xhat}) can be replaced by $\widehat{\mathbf A}_1=(\widehat{\mathbf a}_{1,1},...,\widehat{\mathbf a}_{1,\widehat r_1})$ and $\widehat{\mathbf A}_2=(\widehat{\mathbf a}_{2,1},...,\widehat{\mathbf a}_{2,p-\widehat r_1})$.
Next, turn to the estimation of $r_2$, which is the number of stationary common factors.
Because Model (\ref{re:st}) is the same as (2.2) in \cite{gaotsay2019b} and $\mathcal{M}(\widehat{\mathbf A}_2)$ is a consistent estimator for $\mathcal{M}({\mathbf A}_2)$, we apply the white noise test procedure there to the transformed data $\widehat{\mathbf x}_{2t}$ to estimate $v$ and use
$\widehat r_2 = p-\widehat r_1 - \widehat v$. If the dimension $p$ is small (say less than 10),
we use a bottom-up procedure to determine the number of
white noise components. Specifically, we test the null hypothesis that $\widehat\xi_{it}$ is a white noise starting with $i=p-\widehat r_1$ using, for example, the well-known Ljung-Box statistic $Q(m)$. Clearly, this testing process can only last until $i=1$. If all transformed series $\widehat \xi_{it}$ are white noise, then $\widehat r_2=0$ and $\widehat v=p-\widehat r_1$. In general, if $\widehat \xi_{it}$ is not a white noise series but $\widehat \xi_{jt}$ are for $j=i+1,...,p$, then $\widehat r_2=i$ and $\widehat v=p-\widehat r_1-i$, and we have $\widehat{\mathbf W}=[\widehat{\mathbf U}_1,\widehat{\mathbf V}_1]$, where $\widehat{\mathbf U}_1\in R^{(p-\widehat r_1)\times \widehat r_2}$ and $\widehat{\mathbf V}_1\in R^{(p-\widehat r_1)\times \widehat v}$.
For large $p$, we propose a modification before
conducting the high-dimensional white noise test.
Specifically, let $\widehat {\mathbf W}$ be the matrix of eigenvectors (in the decreasing order of eigenvalues) of the sample matrix $\widehat {\mathbf M}_2$ in Equation (\ref{M2}) and $\widehat\boldsymbol{\xi}_t=\widehat{\mathbf W}'\widehat{\mathbf A}_2'{\mathbf y}_t=(\widehat\xi_{1t},...,\widehat\xi_{p-\widehat r_1,t})'$ be the transformed series.
Note that the ordering of components in $\widehat\boldsymbol{\xi}_t$ is based on the eigenvalues which do not have
any specific association with the temporal dependence of its components. Yet our goal is
to detect the number of white noise components. It is then reasonable to re-order the
components of $\widehat\boldsymbol{\xi}_t$ based on the $p$-values of the Ljung-Box statistics of testing
zero serial correlations in the components of $\widehat\xi_t$.
This re-ordering step enables us to conduct the
high-dimensional white noise test more efficiently starting with all $p-\widehat r_1$ components.
For simplicity, let $\widehat\boldsymbol{\xi}^*_t$ be the re-ordered series of $\widehat\boldsymbol{\xi}_t$, which is of
dimension $p-\widehat r_1$, and $\widehat {\mathbf W}^*$ be the corresponding matrix of eigenvectors.
Our high-dimensional white-noise test starts with
the null hypothesis that $\widehat\boldsymbol{\xi}^*_t$ is
white noise. If the null hypothesis is rejected, we remove the first component from $\widehat\boldsymbol{\xi}^*_t$
and repeat the testing procedure. If all null hypotheses cannot be rejected, then the number of
white noises is $\widehat v = p-\widehat r_1$. In general, $\widehat v$ is the dimension of the
subset of $\widehat\boldsymbol{\xi}^*_t$ for which the white-noise hypothesis cannot be rejected.
The number of stationary common factors is then $\widehat r_2 = p-\widehat r_1-\widehat v$.
Once $\widehat v$ and $\widehat r_2$ are estimated, $\widehat {\mathbf U}_1$ and $\widehat{\mathbf V}_1$ can be estimated
accordingly using columns of $\widehat{\mathbf W}^*$. Simulation studies in Section 4.2 suggest that the performance of detecting the number of stationary factors in both small and large dimensions benefit from this reordering procedure, especially when the sample size is small.
\section{Theoretical Properties}
In this section, we investigate some asymptotic theory for the estimation methods used in the paper. Starting with the assumption that the number of common factors $r_1$ and $r_2$ are known, we divide the derivations into two cases depending on the dimension $p$. The case of estimated $r_1$ and $r_2$ is discussed later. In the derivation, we adopt the discrepancy measure used by
\cite{panyao2008}: for two $p\times r$ half orthogonal
matrices ${\bf H}_1$ and ${\bf H}_2$ satisfying the condition ${\bf
H}_1'{\bf H}_1={\bf H}_2'{\bf H}_2={\mathbf I}_{r}$, the difference
between the two linear spaces $\mathcal{M}({\bf H}_1)$ and
$\mathcal{M}({\bf H}_2)$ is measured by
\begin{equation}
D({\bf H}_1,{\bf
H}_2)=\sqrt{1-\frac{1}{r}\textrm{tr}({\bf H}_1{\bf H}_1'{\bf
H}_2{\bf H}_2')}.\label{eq:D}
\end{equation}
Note that $D({\bf H}_1,{\bf H}_2) \in [0,1].$
It is equal to $0$ if and only if
$\mathcal{M}({\bf H}_1)=\mathcal{M}({\bf H}_2)$, and to $1$ if and
only if $\mathcal{M}({\bf H}_1)\perp \mathcal{M}({\bf H}_2)$. Denote the $I(1)$ factors by ${\mathbf x}_{1t}=(x_{1,1t},\ldots,x_{1,r_1t})'$ and define ${\mathbf w}_t=(w_{1t}, \ldots, w_{r_1t})'$, where $w_{it}\equiv \nabla x_{1,it}$, which is $I(0)$ for $1\leq i\leq r_1$.
For simplicity, we assume $\mathrm{E}{\mathbf w}_t=\bf 0$ and denote the vector of partial sums
of the components of ${\mathbf w}_t$ by
${\mathbf S}_n^{r_1}({\mathbf t})\equiv (S_n^1(t_1),...,S_n^{r_1}(t_{r_1}))'=\left(\frac{1}{\sqrt{n}}\sum_{s=1}^{[nt_1]}w_{1s},...,\frac{1}{\sqrt{n}}\sum_{s=1}^{[nt_{r_1}]}w_{r_1s}\right)',$
where $0< t_1<...<t_p\leq 1$ are constants and ${\mathbf t}=(t_1,...,t_{r_1})'$.
We always assume that the processes ${\mathbf w}_t$ and $({\mathbf x}_{2t},{\mathbf f}_{2t})$ are weakly stationary
and $\alpha$-mixing, that is, their
mixing coefficients $\alpha_p(k)\rightarrow 0$ as $k\rightarrow\infty$, 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 w}_t:i\leq t\leq j\}$ or $\{({\mathbf x}_{2t},{\mathbf f}_{2t}):i\leq t\leq j\}$. We use $c$ or $C$ as a finite positive generic constant below.
\subsection{Asymptotic Properties When $p$ is Fixed, but $n\rightarrow\infty$}
We begin with some assumptions.
\begin{assumption}
(i) The process $\{{\mathbf w}_t\}$ is $\alpha$-mixing with the mixing coefficient satisfying the condition $\alpha_p(k)\leq \exp(-cn^{\gamma_1})$ for some constants $c>0$ and $\gamma_1\in(0,1]$. \\
(ii) For any $x>0$, $\sup_{t>0,1\leq i\leq r_1}P(|w_{it}|>x)\leq c\exp(-cx^{\gamma_2})$ for constants $c>0$, $\gamma_2\in(0,2]$.\\
(iii) There exists a Gaussian process ${\mathbf W}({\mathbf t})=(W_1(t_1),...,W_{r_1}(t_{r_1}))'$ such that as $n\rightarrow \infty$,
${\mathbf S}_n^{r_1}({\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 Skorokhod $J_1$ topology, and ${\mathbf W}({\bf 1})$ has a positive definite covariance matrix $\boldsymbol{\Omega}=[\sigma_{ij}]$.
\end{assumption}
\begin{assumption}
The process $\{({\mathbf x}_{2t},{\mathbf f}_{2t})\}$ is $\alpha$-mixing with the mixing coefficient satisfying the condition $\sum_{k=1}^\infty\alpha_p(k)^{1-2/\gamma}<\infty$ for some $\gamma>2$.
\end{assumption}
\begin{assumption}
(i) $E|f_{2,it}|^{2\gamma}<C$ and $E|{\varepsilon}_{jt}|^{2\gamma}<C$ for $1\leq i\leq r_2$ and $1\leq j\leq v$, where $\gamma$ is given in Assumption 2.\\
(ii) For any ${\mathbf h}_1\in R^{r_1}$ and ${\mathbf h}_2\in R^{v}$ with $\|{\mathbf h}_1\|_2=\|{\mathbf h}_2\|_2=1$, there exists a constant $C>0$ such that
$P(|{\mathbf h}_1'{\mathbf f}_{2t}|>x)\leq 2\exp(-Cx)$ and $P(|{\mathbf h}_2'\mbox{\boldmath$\varepsilon$}_t|>x)\leq 2\exp(-Cx)$ for any $x>0$.
\end{assumption}
The mixing conditions in Assumptions 1(i) and 2 are standard for dependent random processes. See \cite{gaoetal2017} for a theoretical justification for VAR models. Assumption 1(ii) implies that the tail probabilities of linear combinations of ${\mathbf w}_t$ decay exponentially fast, which is used to establish the Bernstein-type inequality of (\ref{bern}) together with Assumption 1(i) by Theorem 1 of \cite{merlevede2011}. The restrictions of $\gamma_1$ and $\gamma_2$ are introduced only for the presentation convenience and large $\gamma_1$ and $\gamma_2$ would only make the conditions therein stronger.
Assumption 1(iii) is also used in \cite{zhang-etal2019}. The conditions in Assumption 3(i) implies that $E|y_{it}|^{2\gamma}<C$ under the setting that $p$ is fixed, and Assumption 3(ii) suggests that ${\mathbf f}_{2t}$ and $\mbox{\boldmath$\varepsilon$}_t$ are sub-exponential, which is a larger class of distributions than sub-Gaussian. The following theorem establishes the consistency of the estimated loading matrices $\widehat{\mathbf A}_1$, $\widehat{\mathbf A}_2$, $\widehat{\mathbf A}_2\widehat{\mathbf U}_1$, its orthonormal complement $\widehat{\mathbf A}_2\widehat{\mathbf V}_1$, the matrix $\widehat{\mathbf A}_2\widehat{\mathbf V}_2$, and the extracted common factors $\widehat{\mathbf A}_1\widehat{\mathbf x}_{1t}$ and $\widehat{\mathbf A}_2\widehat{\mathbf U}_1\widehat{\mathbf z}_{2t}$.
\begin{theorem}\label{tm1}
Suppose Assumptions 1--3 hold and $r_1$ and $r_2$ are known and fixed. Then, for fixed $p$,
\begin{equation}\label{Ai}
D(\mathcal{M}(\widehat{\mathbf A}_i),\mathcal{M}({\mathbf A}_i))=O_p(n^{-1}),\,\, \text{for}\,\, i=1,2,
\end{equation}
\begin{equation}\label{A2U1}
D(\mathcal{M}(\widehat{\mathbf A}_2\widehat{\mathbf U}_{1}),\mathcal{M}({\mathbf A}_2{\mathbf U}_{1}))=O_p(n^{-1/2}), D(\mathcal{M}(\widehat{\mathbf A}_2\widehat{\mathbf V}_{1}),\mathcal{M}({\mathbf A}_2{\mathbf V}_{1}))=O_p(n^{-1/2}),
\end{equation}
\begin{equation}\label{A2V2}
D(\mathcal{M}(\widehat{\mathbf A}_2\widehat{\mathbf V}_{2}),\mathcal{M}({\mathbf A}_2{\mathbf V}_{2}))=O_p(n^{-1/2}).
\end{equation}
As a result,
\begin{equation}\label{ext:ft}
\|\widehat{\mathbf A}_1\widehat{\mathbf A}_1'{\mathbf y}_t-{\mathbf A}_1{\mathbf x}_{1t}\|_2=O_p(n^{-1/2})\quad\text{and}\quad \|\widehat{\mathbf A}_2\widehat{\mathbf U}_1\widehat{\mathbf z}_{2t}-{\mathbf A}_2{\mathbf U}_1{\mathbf z}_{2t}\|_2=O_p(n^{-1/2}).
\end{equation}
\end{theorem}
Similar results are also obtained in \cite{zhang-etal2019} when each component of ${\mathbf y}_t$ is $I(1)$.
The standard $\sqrt{n}$-consistency is still achieved because the super-consistency of ${\mathbf A}_i$ does not affect the second step asymptotically.
If we further assume that the largest $r_2$ eigenvalues of ${\mathbf M}_2$ are distinct,
${\mathbf U}_1$ can be uniquely defined and estimated if we ignore the signs of the columns, as described in the illustrations of Assumption 3 of \cite{gaotsay2019b} and Theorem 1 therein. Nevertheless, we just estimate the corresponding linear space spanned by its columns, and this will not alter the uniquenesses of $\widehat{\mathbf A}_1\widehat{\mathbf A}_1'{\mathbf y}_t$ and $\widehat{\mathbf A}_2\widehat{\mathbf U}_1\widehat{\mathbf z}_{2t}$ since any orthonormal rotations of the columns of $\widehat{\mathbf A}_1$ or $\widehat{\mathbf A}_2\widehat{\mathbf U}_1$ will be canceled out by its inverse associated with the rest matrix product, as illustrated by \cite{lamyao2012} for equation (2.3) therein.
The next theorem states that the proposed method in Section 2.3 can estimate $r_1$ and $r_2$ consistently.
\begin{theorem}
Under Assumptions 1--3, $ P(\widehat r_1=r_1)\rightarrow 1$ and $P(\widehat r_2=r_2)\rightarrow 1$ as $n\rightarrow\infty$.
\end{theorem}
\subsection{Asymptotic Properties When $n\rightarrow\infty$ and $p\rightarrow\infty$}
{We consider the case when $p\rightarrow\infty$ with $p=O(n^c)$ for some constant $c$ which will be specified in the theorems below.}
\begin{assumption}
There exists a constant $\delta\in [0,1)$ such that $\sigma_1({\mathbf L}_1)\asymp$ ... $\asymp\sigma_{r_1}({\mathbf L}_1)\asymp p^{(1-\delta)/2}$, $\sigma_1({\mathbf U}_{22,1})\asymp$ ... $\asymp\sigma_{r_2}({\mathbf U}_{22,1})\asymp p^{(1-\delta)/2}$, $\sigma_1({\mathbf U}_{22,2})\asymp$ ... $\asymp\sigma_{K}({\mathbf U}_{22,2})\asymp p^{(1-\delta)/2}$, and $\sigma_{K+1}({\mathbf U}_{22,2})\asymp$ ... $\asymp\sigma_{v}({\mathbf U}_{22,2})\asymp 1$, where $0\leq K< v$ is an integer.
\end{assumption}
\begin{assumption}
(i) The process $\{p^{-(1-\delta)/2}{\mathbf w}_t\}$ is $\alpha$-mixing with the mixing coefficient satisfying the condition $\alpha_p(k)\leq \exp(-cn^{\gamma_1})$ for some constants $c>0$ and $\gamma_1\in(0,1]$. \\
(ii) For any $x>0$, $\sup_{t>0,1\leq i\leq r_1}P(|p^{-(1-\delta)/2}w_{it}|>x)\leq c\exp(-cx^{\gamma_2})$ for $c>0$, $\gamma_2\in(0,2]$.\\
(iii) There exists a Gaussian process ${\mathbf W}({\mathbf t})=(W_1(t_1),...,W_{r_1}(t_{r_1}))'$ such that as $n\rightarrow \infty$,
$p^{-\frac{1-\delta}{2}}{\mathbf S}_n^{r_1}({\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 Skorohod $J_1$ topology, and ${\mathbf W}({\bf 1})$ has a positive definite covariance matrix $\boldsymbol{\Omega}=[\sigma_{ij}]$.
\end{assumption}
\begin{assumption}
(i) For $\gamma$ given in Assumption 2, any ${\mathbf h}\in R^{v}$ and $0<c_h<\infty$ with $\|{\mathbf h}\|_2=c_h$, $E|{\mathbf h}'\mbox{\boldmath$\varepsilon$}_t|^{2\gamma}<\infty$; (ii) $\sigma_{\min}({\mathbf R}'{\mathbf V}_2^{*}{'}{\mathbf U}_1)\geq C_3$ for some constant $C_3>0$ and some half orthogonal matrix ${\mathbf R}\in R^{(p-r_1-K)\times r_2}$ satisfying ${\mathbf R}'{\mathbf R}={\mathbf I}_{r_2}$, where $\sigma_{\min}$ denotes the minimum non-zero singular value of a matrix.
\end{assumption}
Assumption 4 is similar to Assumption 5 of \cite{gaotsay2019b} to quantify the strength of the factors and the noises, and Assumption 6 is the same as the one therein. Assumption 5 is similar to Assumption 1 in Section 3.1 and we take the strength of the nonstationary factors into account according to Assumption 4 and the decomposition of ${\mathbf L}$ in Section 2.1.
\begin{theorem}\label{tm3}
Suppose Assumptions 2--6 hold and $r_1$ is known and fixed. If $r_2=o(\min(p^\delta,p^{1-\delta}))$ and $p=o(n^{1/(1+\delta)})$, then
\begin{equation}\label{Ai:p}
D(\mathcal{M}(\widehat{\mathbf A}_i),\mathcal{M}({\mathbf A}_i))=O_p(p^{1/2}n^{-1}),\,\, \text{for}\,\, i=1,2.
\end{equation}
Furthermore,
\begin{equation}\label{A2U1:p}
D(\mathcal{M}(\widehat{\mathbf A}_2\widehat{\mathbf U}_{1}),\mathcal{M}({\mathbf A}_2{\mathbf U}_{1}))=O_p(\frac{p^{(1+\delta)/2}}{n^{1/2}}), D(\mathcal{M}(\widehat{\mathbf A}_2\widehat{\mathbf V}_{1}),\mathcal{M}({\mathbf A}_2{\mathbf V}_{1}))
=O_p(\frac{p^{(1+\delta)/2}}{n^{1/2}}),
\end{equation}
\begin{equation}\label{A2V2:p}
D(\mathcal{M}(\widehat{\mathbf A}_2\widehat{\mathbf V}_{2}^{*}),\mathcal{M}({\mathbf A}_2{\mathbf V}_{2}^{*}))=O_p(p^{(1+\delta)/2}n^{-1/2}).
\end{equation}
Consequently,
\begin{equation}\label{ext:ft:1}
p^{-1/2}\|\widehat{\mathbf A}_1\widehat{\mathbf A}_1'{\mathbf y}_t-{\mathbf A}_1{\mathbf x}_{1t}\|_2=O_p(p^{(1-\delta)/2}n^{-1/2}),
\end{equation}
\begin{equation}\label{ext:ft:2}
p^{-1/2}\|\widehat{\mathbf A}_2\widehat{\mathbf U}_1\widehat{\mathbf z}_{2t}-{\mathbf A}_2{\mathbf U}_1{\mathbf z}_{2t}\|_2=O_p(p^{1/2}r_2^{1/2}n^{-1/2}+p^{-1/2}).
\end{equation}
\end{theorem}
From Theorem~\ref{tm3}, the convergence rate of $\mathcal{M}(\widehat{\mathbf A}_1)$ and $\mathcal{M}(\widehat{\mathbf A}_2)$ does not depend on the strength $\delta$ and coincides with that in Theorem 3.3 of \cite{zhang-etal2019} if assuming $r_1$ is small. Under this assumption and even with a slowly growing $r_2$, we can handle the cases when the dimension $p=o(n^{1/(1+\delta)})$ which is higher than the maximal rate of $p=o(n^{1/2-\tau})$ in \cite{zhang-etal2019} for some $0<\tau<1/2$. In addition, if $p=o(n^{1/(1+\delta)})$, all the remaining estimators are all consistent as shown in Theorem~\ref{tm3}. The next theorem states the consistency of the estimated $\widehat r_1$ and $\widehat r_2$ with additional constraints on $p$.
\begin{theorem}\label{tm4}
Let Assumptions 2--6 hold, $r_1$ be fixed and $r_2=o(\min(p^\delta,p^{1-\delta}))$.
(i) if $p=o(n^{1/(1+\delta)})$, $ P(\widehat r_1=r_1)\rightarrow 1$ as $n\rightarrow\infty$; and
(ii) if $p^{1+\delta/2}n^{-1/2}\log(np)=o(1)$, $P(\widehat r_2=r_2)\rightarrow 1$ as $n\rightarrow\infty$.
\end{theorem}
\section{Numerical Properties}
In this section, we study the finite-sample performance of the proposed methodology under the scenarios when $p$ is both small and large.
As the dimensions of $\widehat{\mathbf A}_1$ and ${\mathbf A}_1$ are not necessarily the same,
and ${\mathbf L}_1$ is not an orthogonal matrix in general, we first extend the discrepancy measure
in Equation (\ref{eq:D}) to a more general form below. Let ${\mathbf H}_i$ be a
$p\times d_i$ matrix with rank$({\mathbf H}_i) = d_i$, and ${\mathbf P}_i =
{\mathbf H}_i({\mathbf H}_i'{\mathbf H}_i)^{-1} {\mathbf H}_i'$, $i=1,2$. Define
\begin{equation}\label{dmeasure}
\bar{D}({\mathbf H}_1,{\mathbf H}_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 H}_1,{\mathbf H}_2)=0$ if and only if
either $\mathcal{M}({\mathbf H}_1)\subset \mathcal{M}({\mathbf H}_2)$ or
$\mathcal{M}({\mathbf H}_2)\subset \mathcal{M}({\mathbf H}_1)$, and it is 1 if and only if
$\mathcal{M}({\mathbf H}_1) \perp \mathcal{M}({\mathbf H}_2)$.
If $d_1 = d_2=r$ and ${\mathbf H}_i'{\mathbf H}_i= {\mathbf I}_r$, then
$\bar{D}({\mathbf H}_1,{\mathbf H}_2)$ reduces to that in Equation (\ref{eq:D}).
We only present the simulation results
for $k_0=2$ and $j_0=2$ in Equations (\ref{wy}) and (\ref{M2}), respectively,
to save space because other choices of $k_0$ and $j_0$ produce similar results.\\
\subsection{Simulation}
{\bf Example 1.} Consider the models in (\ref{w:eq}) and (\ref{st:eq}) with common factors following
\[{\mathbf x}_{1,t}={\mathbf x}_{1,t-1}+\boldsymbol{\eta}_{1,t}\quad\text{and}\quad {\mathbf f}_{2,t}=\boldsymbol{\Phi} {\mathbf f}_{2,t-1}+\boldsymbol{\eta}_{2,t},\]
where $\boldsymbol{\eta}_{1,t}$ and $\boldsymbol{\eta}_{2,t}$ are independent white noise processes. We set the true numbers of factors $r_1=2$ and $r_2=2$, the dimension $p=6, 10, 15, 20$, and the sample size $n=200, 500, 1000, 1500, 2000$. For each dimension $p$, we set the seed value to \texttt{1234} in \texttt{R} and generate first a matrix ${\mathbf M}\in R^{p\times p}$ with each element independently drawn from $U(-2,2)$, then use the method and package of \cite{hoff2009} to simulate a random orthonormal matrix ${\mathbf A}$ of Model (\ref{w:eq}) that follows a matrix-variate von Mises-Fisher distribution for a given matrix ${\mathbf M}$. The elements of ${\mathbf U}_{22,1}$ and ${\mathbf U}_{22,2}$ are drawn independently from $U(-1,1)$, and the elements of ${\mathbf U}_{22,2}$ are then divided by $\sqrt{p}$ to balance the accumulated variances of $f_{2,it}$ and ${\varepsilon}_{it}$ for each component of ${\mathbf x}_{2,t}$. $\boldsymbol{\Phi}$ is a diagonal matrix with its diagonal elements being drawn independently from $U(0.5,0.9)$. For each realization, $\mbox{\boldmath$\varepsilon$}_t\sim N(0,{\mathbf I}_v)$, $\boldsymbol{\eta}_{1,t}\sim N(0,{\mathbf I}_{r_1})$ and $\boldsymbol{\eta}_{2,t}\sim N(0,{\mathbf I}_{r_2})$, we use $500$ replications for each $(p,n)$ configuration. For the estimation of $\widehat r_1$, we use `a' to denote using the average of sample autocorrelations and
`a$^*$' the average of the absolute autocorrelations as described in Section 2.3. We denote `w$^*$' and `w' the method of white noise test with and without the reordering procedure
discussed in Section 2.3, respectively. `a$^*$w$^*$' is the estimation of $r_2$ using `w$^*$' based on the estimated $r_1$ by `a$^*$', and `aw' is similarly defined.
We first study the performance of estimating the numbers of unit-root and stationary factors. We use $l=3$ and $m=10$ in (\ref{sm}) and $c_0=0.3$ to determine $r_1$. Since $p$ is relatively small compared to the sample size $n$, for each iteration, we use Ljung-Box test statistics $Q(m)$ with
$m$ = 10 to determine the number of common factors $r_2$.
The empirical probabilities $P(\widehat r_1=r_1)$, $P(\widehat r_2=r_2)$, and $P(\widehat r_1+\widehat r_2=r_1+r_2)$ are reported in Table~\ref{Table1} for different methods. From the table, we see that, for a given $p$, the performance of all methods improves as the sample size increases.
To detect the number of unit roots, the proposed method that uses average of absolute autocorrelations performs better than the one without using the absolute values when the sample size $n$ is not large. Also, the proposed method with reordering procedure also outperforms the one without reordering, especially when $n$ is relatively small. On the other hand, for a given $n$ and method, the empirical probability of correct specification
decreases as $p$ increases, especially for the determination of $r_2$. This is understandable
because it is harder to determine the correct number of stationary factors when the dimension increases and the errors in the white noise testing procedure accumulates. In addition, when the sample size is small and the dimension is relatively high (e.g., $n=200$ and
$p > 10$), the empirical probabilities $P(\widehat r_1=r_1)$ and $P(\widehat r_2=r_2)$ are not satisfactory, but the total number of factors ($r_1+ r_2$) can still be estimated reasonably well if the dimension is low, say $p \leq 10$. Finally, for large $n$ (say, $n \geq 500$) and the low dimension $p$ considered in the simulation, the impact
of modifications discussed in Section 2.3 is small. This is not surprising as the serial dependence
of the stationary models used is positive and decays quickly.
Overall, the estimation of $r_1$ performs well when the sample size is sufficiently large,
and the Ljung-Box test to determine $r_2$ works well for the
case of small dimension (e.g., $p\leq 10$). However, when $p$ is larger (e.g., $p=15, 20$), the Ljung-Box test statistic tends to overestimate the number of stationary factors, implying that we can still keep sufficient information of the original process ${\mathbf y}_t$.
To illustrate the consistency in estimating the loading matrices of the proposed methods that use absolute autocorrelations and the reordering procedure in white noise testing, we present the boxplots of $\bar{D}(\widehat{\mathbf A}_1,{\mathbf A}_1)$, $\bar{D}(\widehat{\mathbf A}_2,{\mathbf A}_2)$ and $\bar{D}(\widehat{\mathbf A}_2\widehat{\mathbf U}_1,{\mathbf A}_2{\mathbf U}_{22,1})$ in Figure~\ref{fig1}(a), (b) and (c), respectively, where $\bar{D}(\cdot,\cdot)$ is defined in (\ref{dmeasure}). From Figure~\ref{fig1}, for each dimension $p$, the discrepancy decreases as the sample size increases and this is in agreement with the theory.
The plots also show that, as expected,
the mean discrepancy increases as the dimension $p$ increases.
\begin{table}
\caption{Empirical probabilities $P(\widehat{r}_1=r)$, $P(\widehat{r}_2=r_2)$ and $P(\widehat r_1+\widehat r_2 =r)$ of various $(p,n)$ configurations
for the model of Example 1 with $r_1=2$, $r_2=2$ and $r=r_1+r_2=4$, where $p$ and $n$ are the dimension and the sample size, respectively. `a' denotes the method of using the average of sample autocorrelations and
`a$^*$' the average of absolute sample autocorrelations, `w$^*$' and `w' the method of white noise test with and without the reordering procedure of Section 2.3, respectively. `a$^*$w$^*$' is the estimation of $r_2$ using `w$^*$' based on the estimated $r_1$ by `a$^*$', and `aw' is similarly defined. $500$ iterations are used.}
\label{Table1}
\begin{center}
\setlength{\abovecaptionskip}{0pt}
\setlength{\belowcaptionskip}{3pt}
\scriptsize
\begin{tabular}{c|c|r|rrrrr}
\hline
&&&\multicolumn{5}{c}{$n$}\\
$p$&EP&Methods&$200$&$500$&$1000$&$1500$&$3000$\\
\hline
&$P(\widehat r_1=r_1)$&a$^*$(a)&0.874(0.682)&1(0.996)&1(1)&1(1)&1(1)\\
$6$&$P(\widehat r_2=r_2)$& a$^*$w$^*$(aw)&0.788(0.604)&0.902(0.898)&0.906(0.906)&0.908(0.908)&0.914(0.914)\\
&$P(\widehat r_1+\widehat r_2=r)$&a$^*$w$^*$(aw)&0.908(0.904)&0.902(0.902)&0.906(0.906)&0.908(0.906)&0.914(0.914)\\
\hline
&$P(\widehat r_1=r_1)$&a$^*$(a)&0.844(0.690)&1(1)&1(1)&1(1)&1(1)\\
$10$&$P(\widehat r_2=r_2)$&a$^*$w$^*$(aw)&0.606(0.512)&0.740(0.740)&0.732(0.732)&0.726(0.726)&0.762(0.762)\\
&$P(\widehat r_1+\widehat r_2=r)$&a$^*$w$^*$(aw)&0.716(0.724)&0.740(0.740)&0.732(0.732)&0.726(0.726)&0.762(0.762)\\
\hline
&$P(\widehat r_1=r_1)$&a$^*$(a)&0.780(0.656)&0.996(0.994)&1(1)&1(1)&1(1)\\
$15$&$P(\widehat r_2=r_2)$&a$^*$w$^*$(aw)&0.420(0.338)&0.544(0.540)&0.586(0.586)&0.592(0.592)&0.562(0.562)\\
&$P(\widehat r_1+\widehat r_2=r)$&a$^*$w$^*$(aw)&0.524(0.520)&0.544(0.544)&0.586(0.586)&0.592(0.592)&0.562(0.562)\\
\hline
&$P(\widehat r_1=r_1)$&a$^*$(a)&0.678(0.526)&0.988(0.982)&1(1)&1(1)&1(1)\\
$20$&$P(\widehat r_2=r_2)$&a$^*$w$^*$(aw)&0.286(0.226)&0.390(0.390)&0.420(0.420)&0.434(0.434)&0.482(0.482)\\
&$P(\widehat r_1+\widehat r_2=r)$&a$^*$w$^*$(aw)&0.406(0.396)&0.398(0.400)&0.420(0.420)&0.434(0.434)&0.482(0.482)\\
\hline
\end{tabular}
\end{center}
\end{table}
\begin{figure}
\begin{center}
\subfigure[]{\includegraphics[width=0.325\textwidth]{da1-psm-n.pdf}}
\subfigure[]{\includegraphics[width=0.325\textwidth]{da2-psm-n.pdf}}
\subfigure[]{\includegraphics[width=0.325\textwidth]{dau-psm-n.pdf}}
\caption{Estimation results of Example 1 when $p$ is relatively small and
$r_1=r_2=2$: (a) Boxplots of $\bar{D}(\widehat{\mathbf A}_1,{\mathbf A}_1)$; (b) Boxplots of $\bar{D}(\widehat{\mathbf A}_2,{\mathbf A}_2)$; (c) Boxplots of $\bar{D}(\widehat{\mathbf A}_2\widehat{\mathbf U}_1,{\mathbf A}_2{\mathbf U}_{22,1})$. The sample sizes used are $200, 500, 1000, 1500, 3000$, and the results are based on $500$ iterations.}\label{fig1}
\end{center}
\end{figure}
Next, for each $(p,n)$ configuration, we study the root-mean-square errors (RMSEs):
\begin{equation}\label{rmse:psm}
\text{RMSE}_1=[\frac{1}{n}\sum_{t=1}^n\|\widehat{\mathbf A}_1\widehat{\mathbf x}_t-{\mathbf A}_1{\mathbf x}_{1t}\|_2^2]^{1/2},\text{RMSE}_2=[\frac{1}{n}\sum_{t=1}^n\|\widehat{\mathbf A}_2\widehat{\mathbf U}_1\widehat{\mathbf z}_{2t}-{\mathbf A}_2{\mathbf U}_{22,1}{\mathbf f}_{2t}\|_2^2]^{1/2},
\end{equation}
which quantify the accuracy in estimating the common factor processes.
Boxplots of the RMSE$_1$ and RMSE$_2$ are shown in Figure~\ref{fig2}(a)-(b), respectively. From the plots, we observe a clear pattern that, as the sample size increases, the RMSEs decrease for a given $p$, which is consistent with the results of Theorem \ref{tm1}.
Overall, the one-by-one testing procedure works well when the dimension is small, and the RMSEs decrease when the sample size increases, even though the performance of the
Ljung-Box test may deteriorate due to the overestimation of the number of the
stationary common factors for higher dimension $p$.
\begin{figure}
\begin{center}
\subfigure[]{\includegraphics[width=0.45\textwidth]{dx1-psm-n.pdf}}
\subfigure[]{\includegraphics[width=0.45\textwidth]{dx2-psm-n.pdf}}
\caption{Estimation accuracies of Example 1
when $p$ is relatively small and $r_1=r_2=2$:
(a) Boxplots of the RMSE$_1$ defined in (\ref{rmse:psm}); (b) Boxplots of the RMSE$_2$ defined in (\ref{rmse:psm}). The sample sizes used are $200, 500, 1000, 1500, 3000$, and the results are based on $500$ iterations.}\label{fig2}
\end{center}
\end{figure}
\noindent{\bf Example 2.} In this example, we consider Models (\ref{w:eq}) and (\ref{st:eq}) with ${\mathbf x}_{1t}$ and ${\mathbf f}_{2t}$
being the same as those of Example 1. We set the true numbers of factors $r_1=4$, $r_2=6$ and the number of prominent components of the noise covariance $K=2$ as defined in Assumption 4. The dimensions
used are $p=50, 100, 300, 500$, and the sample sizes are $n=300, 500, 1000, 1500, 2000$. We consider two scenarios for $\delta$ defined in Assumption 4: $\delta=0$ and $\delta=0.5$. For each setting, since it is time consuming to simulate a random orthonormal matrix when the dimension is high, instead we first simulate a matrix ${\mathbf M}\in R^{p\times p}$ with each element drawn from $U(-2,2)$ using the same seed value as Example 1, then perform the SVD on ${\mathbf M}$, and choose the columns of ${\mathbf A}$ as the left singular vectors of ${\mathbf M}$ multiplied by $p^{(1-\delta)/2}$. The elements of ${\mathbf U}_{22,1}$ and ${\mathbf U}_{22,2}$ are drawn independently from $U(-1,1)$, then we divide ${\mathbf U}_{22,1}$ by $p^{\delta/2}$, the first $K$ columns of ${\mathbf U}_{22,2}$ by $p^{\delta/2}$ and the remaining $v-K$ columns by $p$ to satisfy Assumption 4.
$\boldsymbol{\Phi}$, $\mbox{\boldmath$\varepsilon$}_t$, $\boldsymbol{\eta}_{1,t}$ and $\boldsymbol{\eta}_{2,t}$ are drawn similarly as those of Example 1. We also use $500$ replications in each experiment and denote the methods by `a', `a$^*$', `w', and w$^*$ as in Example 1.
We first study the performance of $S_i(l,m)$ in (\ref{sm}) to estimate the number of unit-root series and that of high-dimensional white-noise tests of \cite{gaotsay2019b} to estimate the number of stationary common factors with and without the re-ordering modification of Section 2.3.
The choices of $m$, $c_0$ and $l$ are the same as those in Example 1. For simplicity, we only present the results of the $T(m)$ statistics there to estimate $r_2$, and
the results for the other white-noise test are similar.
When $(p-\widehat r_1)\geq n$, we only keep the upper ${{\varepsilon}} n$ components of $\widehat\boldsymbol{\xi}_t=\widehat{\mathbf W}'\widehat{\mathbf A}_2'{\mathbf y}_t$ with ${{\varepsilon}}=0.75$ in white-noise testing. The results are reported in Table~\ref{Table2}.
From the table, we see that for each setting of $\delta$ and fixed $p$, the performance of all methods improves as the sample size increases except for some cases in which
$P(\widehat r_2=r_2)$ of $n=500$ is smaller than that of $n=300$ for $p=500$. This is understandable since the actual dimension used in the white noise testing for $n=300$ is $0.75n=225$, while $(500-\widehat r_1)$ for $n=500$ is much larger. For the estimation of $r_1$, the proposed method using absolute autocorrelations fares better when the sample size $n$ is small, and it is comparable with the one without using the absolute values when $n$ is large for each $p$. When estimating $r_2$, the proposed method consisting the reordering procedure (a$^*$w$^*$) outperforms
the other in most cases when $n$ is relatively small ($n=300, 500$)
and the two methods are comparable with each other when the sample size is large.
As expected, when the factor strength is
strong ($\delta = 0$), most procedures work better. The table also shows that
the white noise test needs further improvement in selecting the number of stationary
common factor $r_2$ when $p \geq n$, but the reordering procedure provides substantial
improvements over the one without re-ordering.
Overall, the performance of the proposed modeling procedure is satisfactory
when the sample size is larger than the dimension.
\begin{table}
\caption{Empirical probabilities $P(\widehat{r}_1=r)$, $P(\widehat{r}_2=r_2)$ and $P(\widehat r_1+\widehat r_2 =r)$ of various $(p,n)$ configurations
for the model of Example 2 with $r_1=4$, $r_2=6$ and $r=r_1+r_2=10$, where $p$ and $n$ are the dimension and the sample size, respectively. `a' denotes the method of using the average of autocorrelations and
`a$^*$' the average of absolute autocorrelations, `w$^*$' and `w' denote the method of white noise test with and without the reordering procedure of Section 2.3. `a$^*$w$^*$' is the estimation of $r_2$ using `w$^*$' based on the estimated $r_1$ by `a$^*$', and `aw' is similarly defined. $500$ iterations are used.}
\label{Table2}
\begin{center}
\setlength{\abovecaptionskip}{0pt}
\setlength{\belowcaptionskip}{2pt}
\tiny
\begin{tabular}{c|c|c|r|rrrrr}
\hline
&&&&\multicolumn{5}{c}{$n$}\\
$\delta$&$p$&EP&Methods&$300$&$500$&$1000$&$1500$&$2000$\\
\hline
$0$&&$P(\widehat r_1=r_1)$&a$^*$(a)&0.798(0.412)&0.968(0.926)&0.992(0.992)&0.996(0.996)&0.998(0.998)\\
&$50$&$P(\widehat r_2=r_2)$&a$^*$w$^*$(aw)&0.626(0.380)&0.892(0.846)&0.928(0.928)&0.898(0.898)&0.924(0.924)\\
&&$P(\widehat r_1+\widehat r_2=r)$&a$^*$w$^*$(aw)&0.770(0.778)&0.920(0.920)&0.936(0.936)&0.902(0.902)&0.926(0.926)\\
\cline{2-9}
&&$P(\widehat r_1=r_1)$&a$^*$(a)&0.740(0.442)&0.930(0.928)&0.982(0.984)&0.998(0.998)&0.988(0.988)\\
&$100$&$P(\widehat r_2=r_2)$&a$^*$w$^*$(aw)&0.634(0.408)&0.864(0.856)&0.898(0.900)&0.930(0.930)&0.904(0.904)\\
&&$P(\widehat r_1+\widehat r_2=r)$&a$^*$w$^*$(aw)&0.832(0.832)&0.924(0.924)&0.916(0.916)&0.932(0.932)&0.916(0.916)\\
\cline{2-9}
&&$P(\widehat r_1=r_1)$&a$^*$(a)&0.748(0.450)&0.892(0.922)&0.962(0.980)&0.970(0.974)&0.982(0.982)\\
&$300$&$P(\widehat r_2=r_2)$&a$^*$w$^*$(aw)&0.336(0.154)&0.804(0.830)&0.874(0.894)&0.902(0.908)&0.896(0.896)\\
&&$P(\widehat r_1+\widehat r_2=r)$&a$^*$w$^*$(aw)&0.370(0.234)&0.896(0.896)&0.910(0.910)&0.926(0.926)&0.910(0.910)\\
\cline{2-9}
&&$P(\widehat r_1=r_1)$&a$^*$(a)&0.690(0.446)&0.916(0.926)&0.980(0.986)&0.962(0.962)&0.984(0.984)\\
&$500$&$P(\widehat r_2=r_2)$&a$^*$w$^*$(aw)&0.526(0.408)&0.764(0.172)&0.906(0.912)&0.892(0.894)&0.924(0.922)\\
&&$P(\widehat r_1+\widehat r_2=r)$&a$^*$w$^*$(aw)&0.702(0.706)&0.824(0.188)&0.926(0.926)&0.928(0.928)&0.936(0.936)\\
\hline
$0.5$&&$P(\widehat r_1=r_1)$&a$^*$(a)&0.786(0.440)&0.962(0.950)&0.996(0.996)&0.996(0.996)&0.996(0.996)\\
&$50$&$P(\widehat r_2=r_2)$&a$^*$w$^*$(aw)&0.596(0.368)&0.880(0.868)&0.930(0.930)&0.928(0.928)&0.886(0.886)\\
&&$P(\widehat r_1+\widehat r_2=r)$&a$^*$w$^*$(aw)&0.728(0.830)&0.914(0.916)&0.934(0.934)&0.932(0.922)&0.890(0.890)\\
\cline{2-9}
&&$P(\widehat r_1=r_1)$&a$^*$(a)&0.732(0.438)&0.936(0.930)&0.990(0.994)&0.994(0.994)&0.990(0.990)\\
&$100$&$P(\widehat r_2=r_2)$&a$^*$w$^*$(aw)&0.436(0.300)&0.818(0.818)&0.912(0.916)&0.932(0.932)&0.902(0.902)\\
&&$P(\widehat r_1+\widehat r_2=r)$&a$^*$w$^*$(aw)&0.594(0.612)&0.866(0.882)&0.922(0.922)&0.938(0.938)&0.912(0.912)\\
\cline{2-9}
&&$P(\widehat r_1=r_1)$&a$^*$(a)&0.764(0.450)&0.904(0.920)&0.968(0.980)&0.968(0.974)&0.982(0.982)\\
&$300$&$P(\widehat r_2=r_2)$&a$^*$w$^*$(aw)&0.128(0.012)&0.474(0.146)&0.796(0.708)&0.888(0.892)&0.896(0.892)\\
&&$P(\widehat r_1+\widehat r_2=r)$&a$^*$w$^*$(aw)&0.144(0.018)&0.536(0.144)&0.820(0.720)&0.914(0.910)&0.914(0.910)\\
\cline{2-9}
&&$P(\widehat r_1=r_1)$&a$^*$(a)&0.708(0.450)&0.914(0.916)&0.982(0.986)&0.958(0.958)&0.982(0.982)\\
&$500$&$P(\widehat r_2=r_2)$&a$^*$w$^*$(aw)&0.264(0.008)&0.262(0.014)&0.628(0.318)&0.788(0.686)&0.880(0.852)\\
&&$P(\widehat r_1+\widehat r_2=r)$&a$^*$w$^*$(aw)&0.380(0.012)&0.274(0.010)&0.634(0.312)&0.822(0.710)&0.888(0.858)\\
\hline
\end{tabular}
\end{center}
\end{table}
To demonstrate the advantages of using white-noise tests to select the number of stationary
common factors with dynamic dependencies, we compare it with ratio-based method in \cite{lamyao2012}, which defines
$\widehat r_2=\arg\min_{1\leq j\leq R}\left\{\widehat\lambda_{j+1}/\widehat\lambda_j\right\}$, where $\widehat\lambda_1,...,\widehat\lambda_{p-\widehat r_1}$ denote the eigenvalues of $\widehat{\mathbf M}_2$ in (\ref{M2}) and we choose $R=(p-\widehat r_1)/2$ suggested in their paper. Figure~\ref{fig3a} presents the boxplots of $\widehat r_2$ for $p=100$ and $300$ with different sample sizes.
The selection of $\widehat r_2$ centers around 8, instead of the true value 6,
indicating that the ratio-based method may fail to identify the correct number of stationary factors $r_2$ if the noise effect is diverging, while the high-dimensional white noise test continues to work well as shown in Table~\ref{Table2}.
\begin{figure}
\begin{center}
{\includegraphics[width=0.6\textwidth]{ly-ratio-2d0.pdf}}
\caption{Boxplots of Lam and Yao (2012)'s ratio method in selecting the number of
stationary common factors when $r_1=4$, $r_2=6$, $K=2$. Left: $p=100$;
Right: $p=300$ and the results are based on $500$ iterations.}\label{fig3a}
\end{center}
\end{figure}
Next, we study the accuracy of the estimated loading matrices as that of Example 1
using the methods a$^*$ and a$^*$w$^*$.
The boxplots of $\bar{D}(\widehat{\mathbf A}_1,{\mathbf A}_1)$, $\bar{D}(\widehat{\mathbf A}_2,{\mathbf A}_2)$ and $\bar{D}(\widehat{\mathbf A}_2\widehat{\mathbf U}_1,{\mathbf A}_2{\mathbf U}_{22,1})$ are shown in Figure~\ref{fig3}(a)--(c), respectively. Similar patterns are also obtained for the estimation of other matrices so we omit them here.
From Figure~\ref{fig3},
there is a clear indication that the estimation
accuracy of the loading matrix improves as the sample size increases even
for moderately large $p$,
which is in line with the asymptotic theory. The results also confirm that the proposed $S_i(l,m)$ and the white noise
test that select $\widehat r_1$ and $\widehat r_2$, respectively, perform reasonably well even for large $p$.
\begin{figure}
\begin{center}
\subfigure[]{\includegraphics[width=0.325\textwidth]{da1-plg-nd0.pdf}}
\subfigure[]{\includegraphics[width=0.325\textwidth]{da2-plg-nd0.pdf}}
\subfigure[]{\includegraphics[width=0.325\textwidth]{dau-plg-nd0.pdf}}
\caption{Estimation results of Example 2 when $p$ is relatively large, $r_1=4$, $r_2=6$, and
$K=2$; (a) Boxplots of $\bar{D}(\widehat{\mathbf A}_1,{\mathbf A}_1)$; (b) Boxplots of $\bar{D}(\widehat{\mathbf A}_2,{\mathbf A}_2)$;
(c) Boxplots of $\bar{D}(\widehat{\mathbf A}_2\widehat{\mathbf U}_1,{\mathbf A}_2{\mathbf U}_{22,1})$. The sample sizes used are $300, 500, 1000, 1500, 3000$, and the results are based on $500$ iterations.}\label{fig3}
\end{center}
\end{figure}
For each $(p,n)$, we further define the RMSEs for high dimension as:{\small
\begin{equation}\label{rmse:plg}
\text{RMSE}_3=[\frac{1}{np}\sum_{t=1}^n\|\widehat{\mathbf A}_1\widehat{\mathbf x}_t-{\mathbf A}_1{\mathbf x}_{1t}\|_2^2]^{1/2},\text{RMSE}_4=[\frac{1}{np}\sum_{t=1}^n\|\widehat{\mathbf A}_2\widehat{\mathbf U}_1\widehat{\mathbf z}_{2t}-{\mathbf A}_2{\mathbf U}_{22,1}{\mathbf f}_{2t}\|_2^2]^{1/2},
\end{equation} }
which quantify the accuracy in estimating the common factor processes.
When the dimension is moderately large, the
boxplots of the RMSE$_3$ and RMSE$_4$ are shown in Figure~\ref{fig4}(a)-(b), respectively. From the plots, similar to Example 1, we see a clear pattern that, as the sample size increases, the RMSEs decrease for a given $p$, which is consistent with the results of Theorem \ref{tm3}.
Overall, the proposed method works well even for moderately high dimensions.
This is especially so when the sample size is greater than the dimension.
\begin{figure}
\begin{center}
\subfigure[]{\includegraphics[width=0.48\textwidth]{dx1-plg-nd0.pdf}}
\subfigure[]{\includegraphics[width=0.48\textwidth]{dx2-plg-nd0.pdf}}
\caption{Estimation accuracies of Example 2 when $p$ is relatively large,
$r_1=4$, $r_2=6$ and $K=2$: (a) Boxplots of the RMSE$_3$ defined in (\ref{rmse:plg});
(b) Boxplots of the RMSE$_4$ defined in (\ref{rmse:plg}). The sample sizes used are $300, 500, 1000, 1500, 2000$, and the results are based on $500$ iterations.}\label{fig4}
\end{center}
\end{figure}
\subsection{Real Data Analysis}
{\bf Example 3.} Consider the hourly measurements of PM$_{2.5}$ collected by
Airboxes at 508 locations in Taiwan. The locations of the 508 locations are shown in Figure~\ref{TWNlocations}, which mainly consist of three clusters, signifying the major cities
of Taiwan, namely Taipei, Taichung, Tainan, and Kaohsiung. The latter two cities are adjacent and
form a large cluster.
The isolated location outside of Taiwan
denotes part of the Orchid Island of Taiwan.
We apply our proposed method to the hourly measurements of PM$_{2.5}$ for March 2017
with a total of 744 observations. The 508 time series are shown in Figure~\ref{PM25plots}.
We first applied the method of Section 2.3 with $k_0=2$ and found that $\widehat r_1=3$,
i.e., there are 3 unit-root
processes. When applying the method, we choose $c_0=0.3$, $m=30$, and $l=3$. Similar results are also obtained by varying the values of $c_0$, $m$, and $l$ for several choices.
For example, when $(c_0,m,l)=(0.3,15,2)$, $(0.2,15,3)$, and $(0.2,8,4)$, we obtain that $\widehat r_1=3$, 3, and 1, respectively.
It is true that these choices are subjective, but they are not unique to our empirical study.
See, for instance, \cite{bai2004} in selecting the number of factors for different choices of $kmax$.
However, the key message of the detection is that there exist some unit-root common trends in
the hourly PM$_{2.5}$ measurements.
The three recovered unit-root factors and their sample ACFs are shown in Figure~\ref{PM25unitroot}.
From the plots, we see that these three unit-root factors capture most of the trends in the
original data of Figure~\ref{PM25plots} and, as expected, their sample ACFs
decay slowly. The extracted 505-dimensional stationary process is shown in
Figure~\ref{fig8}(a).
\begin{figure}
\begin{center}
{\includegraphics[width=0.8\textwidth]{loc-508.pdf}}
\caption{Locations (latitude vs. longitude) of the 508 AirBoxes across Taiwan of Example 3.}\label{TWNlocations}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
{\includegraphics[width=0.8\textwidth]{pmdata508.pdf}}
\caption{Time plots of hourly measurements of PM$_{2.5}$ of 508 locations across Taiwan in March 2017 of Example 3.}\label{PM25plots}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
{\includegraphics[width=0.62\textwidth]{act-unit.pdf}}
\caption{Time plots of the 3 estimated unit-root common trends by the proposed method and their sample ACFs of Example 3.}\label{PM25unitroot}
\end{center}
\end{figure}
Next, we applied the white noise testing method of Section 2.3
with $k_0 = j_0 = 2$ for Equations (\ref{wy}) and (\ref{M2}) and found
that $\widehat r_2=256$, i.e., there are 256 stationary common factors and the remaining
249 components are white noises.
If we set the parameters $k_0=j_0=1$ and $k_0=j_0=3$, we obtain that $\widehat r_2=292$ and $281$, respectively. But using these factors produces similar results below and, for simplicity,
we only report the analysis with $k_0=j_0=2$. To obtain the extracted stationary factors, by the projected PCA in \cite{gaotsay2019b}, we first examine the eigenvalues of the sample covariance matrix $\widehat{\mathbf S}$ defined in (\ref{Spca}).
From Figure~\ref{fig7}(a), we see that there is an eigenvalue of the covariance matrix of the white noise components that is much larger than the others. Therefore, we choose $K=1$ and the recovered stationary factors and the white noise components are shown in Figure~\ref{fig8}(b) and (c), respectively. From Figure~\ref{fig8}, we see that the 256 stationary common factors in part (b) capture most of the nontrivial dynamic dependencies of the stationary components $\widehat{\mathbf x}_{2t}$ shown in part (a), and the remaining 249 white noise series capture little dynamic information of $\widehat{\mathbf x}_{2t}$.
Consequently, for the 508-dimensional time series of PM$_{2.5}$ measurements, our proposed method recovers 3 unit-root series of common trends, 256 stationary common factors with
non-trivial dynamic dependence, and
249 white noise series which capture part of the contemporaneous variability of the data.
\begin{figure}
\begin{center}
{\includegraphics[width=0.62\textwidth]{ppca-eigen.pdf}}
\caption{(a) The first 20 eigenvalues of $\widehat{\mathbf S}$; (b) The plot of ratios of consecutive eigenvalues of $\widehat{\mathbf S}$ .}\label{fig7}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
{\includegraphics[width=0.6\textwidth]{decom.pdf}}
\caption{(a) Time plots of the recovered 505-dimensional stationary process $\widehat{\mathbf x}_{2t}$; (b) Time plots of the recovered 256 stationary common factors $\widehat{\mathbf z}_{2t}$; (c) Time plots of the 249 extracted white noise series.}\label{fig8}
\end{center}
\end{figure}
Next we examine and compare the forecasting performance of the extracted factors via
the proposed method with different methods available in the literature.
We estimate the models using the data in the time span $[1,\tau]$ with $\tau=600,...,744-h$ for the $h$-step ahead forecasts, i.e., we use the last 6 days of March 2017 for out-of-sample forecasting.
First, we applied the stationary factor approach of \cite{BaiNg_Econometrica_2002} and the nonstationary one in \cite{bai2004}, and found that the numbers of common factors are 20 and 1, respectively. We also set the number of nonstationary factors in \cite{zhang-etal2019} (denoted by ZRY) to $\widehat{r}_1=3$.
Second, following the discussion of Section 2.2, we fit a VARIMA(1,1,0) model to the nonstationary factors $\widehat{\mathbf x}_{1t}$ and scalar AR(1) models to the stationary common factors $\widehat{\mathbf z}_{2t}$, and denote this approach by GT. Since the dimension of the extracted stationary common factor
is high, we employ scalar AR models to provide a quick and simple approximation. We also fit a univariate ARIMA(1,1,0) model to the factor process extracted by \cite{bai2004} and a VAR(1) model to the 20-dimensional factors extracted by \cite{BaiNg_Econometrica_2002}, where the insignificant parameters in the coefficient matrix are removed. We denote these two approaches
by B-2004 and BN-2002, respectively. As a benchmark model, we also employ scalar AR(1)
models to the differenced PM$_{2.5}$ series and denote the result by DFAR.
We compute the $h$-step ahead predictions of ${\mathbf y}_t$ using the predictions of factors and the associated factor loadings. The forecast error is defined as
\begin{equation}\label{fe}
\text{FE}_h=\frac{1}{144-h+1}\sum_{\tau=600}^{744-h}E(\tau,h)\quad\text{with}\quad E(\tau,h)=\frac{1}{\sqrt{p}}\|\widehat{\mathbf y}_{\tau+h}-{\mathbf y}_{\tau+h}\|_2,
\end{equation}
where $p=508$.
Table \ref{Table4.21} reports the 1-step to 4-step ahead forecast errors of Equation (\ref{fe})
for the various models considered. The smallest forecast error of each step is shown in boldface. From the table, we see that our proposed method is capable of producing accurate forecasts and the associated forecast errors based on the extracted stationary and nonstationary factors by our method are smaller than the factors extracted by other methods. In addition, we note that the ZRY approach
and the stationary approach of BN-2002 can also produce relatively small forecast errors,
but the single nonstationary factor process recovered by \cite{bai2004} is not able to make accurate predictions. The benchmark approach DFAR produces results similar to that of B-2004,
which are not satisfactory. For further illustration, the pointwise 1-step ahead forecast errors
of various methods are shown in Figure~\ref{fig9}. From the
plot, we see that the proposed approach outperforms the other four methods in most time periods
of the forecasting subsample.
\begin{table}[h]
\caption{The 1-step to 4-step ahead out-of-sample forecast errors. GT denotes the proposed method, ‘BN-2002’ denotes the stationary approach of \cite{BaiNg_Econometrica_2002}, B-2004 is the nonstationary method of \cite{bai2004}, ZRY denotes the approach of \cite{zhang-etal2019}, and DFAR is the
scalar AR approach to the differenced PM$_{2.5}$ data.
Boldface numbers denote the smallest one for a given model.}
\label{Table4.21}
\begin{center}
\setlength{\abovecaptionskip}{0pt}
\setlength{\belowcaptionskip}{3pt}
\begin{tabular}{c|cccccc}
\hline
&&\multicolumn{5}{c}{Methods}\\
\cline{3-7}
Step&>&ZRY&B-2004&BN-2002&DFAR\\
\hline
1 & & {\bf 6.25} & 11.29 & 45.77&7.79&45.80 \\
2& &{\bf 8.58}&12.03&45.93&9.68&45.95\\
3& &{\bf 10.18}&12.82&46.08&11.12&50.33\\
4& &{\bf 11.50}&13.67&46.25&12.26&46.27\\
\hline
\end{tabular}
\end{center}
\end{table}
\begin{figure}
\begin{center}
{\includegraphics[width=0.65\textwidth]{forecast.pdf}}
\caption{Time plots of the 1-step ahead pointwise forecast errors using various models of Example 3. GT denotes the proposed method, ‘BN-2002’ denotes the stationary approach of \cite{BaiNg_Econometrica_2002}, B-2004 is the nonstationary method of \cite{bai2004}, ZRY denotes the approach in \cite{zhang-etal2019}, and DFAR is the scalar AR approach to the differenced PM$_{2.5}$ data. }\label{fig9}
\end{center}
\end{figure}
We then compare the forecast ability of extracted factors by different approaches and the benchmark method by adopting the asymptotic test in \cite{diebold-1995} to test for equal predictive ability without a rigorous argument. The null hypothesis of interest is that the approaches considered have equal predictive ability, and the alternative is that our proposed method performs better than the others in predicting the PM$_{2.5}$ series.
Table~\ref{Table43} reports the testing results, where we use the method in \cite{andrews1991} to calculate the long-run covariance matrix (L-COV). In the table, (GT)--(ZRY), (GT)--(BN-2002), (GT)--(B-2004), and (GT)--(DFAR) denote the comparisons between our approach and that in ZRY, BN-2002, B-2004, and DFAR, respectively, and the test statistic is the difference between the forecast errors of the two methods involved, as illustrated in \cite{diebold-1995}. From Table~\ref{Table43}, we see that all $p$-values of the test statistics are
small and less than 0.05, indicating that the factors
extracted by our method have better forecast ability than those extracted by
other four methods in predicting all components of the PM$_{2.5}$ measurements. The results in Tables~\ref{Table4.21} and \ref{Table43} are understandable since the factors extracted by our method capture most of the dynamic dependence of the data by eliminating the white noise effect. These extracted common factors are helpful in predictions.
\begin{table}
\caption{Testing for equal predictive ability of different methods, using the asymptotic test of \cite{diebold-1995}. GT denotes the proposed method, ‘BN-2002’ denotes the stationary approach in \cite{BaiNg_Econometrica_2002}, B-2004 is the nonstationary method in \cite{bai2004}, ZRY denotes the approach in \cite{zhang-etal2019}, and DFAR is the componentwise AR approach to the differenced data.}
\label{Table43}
\begin{center}
\setlength{\abovecaptionskip}{0pt}
\setlength{\belowcaptionskip}{3pt}
\begin{tabular}{c|cccccc}
\hline
Step&&(GT)--(ZRY)&(GT)--(BN-2002)&&(GT)--(B-2004)&(GT)--(DFAR)\\
\hline
1 & L-COV&0.161 & 0.015 & &279.5& 266.1 \\
& $p$-value& $\approx 0$& $\approx 0$& &$0.01$ &$0.01$ \\
\hline
2&L-COV&0.154&0.055&&61.48&61.64\\
&$p$-value&$\approx 0$&$\approx 0$&&$\approx 0$&$\approx 0$\\
\hline
3&L-COV&0.263&0.086&&71.14&71.17\\
&$p$-value&$\approx 0$&$\approx 0$&&$\approx 0$&$\approx 0$\\
\hline
4&L-COV&0.273&0.131&&70.20&70.23\\
&$p$-value&$\approx 0$&$0.017$&&$\approx 0$&$\approx 0$\\
\hline
\end{tabular}
\end{center}
\end{table}
To further compare the forecast ability of the nonstationary factors extracted by the proposed method and the one in \cite{bai2004}, we employ the factor-augmented error correction models (FECM) in \cite{banerjee-etal2014}. As the factors in \cite{banerjee-etal2014} are extracted by the method of \cite{bai2004}, we denote the
approach of FECM of \cite{banerjee-etal2014} by Bai-BMM, and our method by GT-BMM. Since the method in \cite{bai2004} identifies 1 factor and the proposed method estimates 3, we include 1 to 3 factors in the FECM defined in Equation (8) of \cite{banerjee-etal2014} with $q=1$. That is, we include the lag-1 of the differenced augmented vectors in the error correction model and denote it
by FECM(1). We study the forecast of the PM$_{2.5}$ in Daliao district of Kaohsiung, which is one of the most polluted areas in Southern Taiwan. Table~\ref{Table4.22} presents the 1-step to 4-step ahead root mean squared forecast errors (RMSFE) defined as
\begin{equation}\label{rmfe}
\text{RMSFE}_h=\left(\frac{1}{144-h+1}\sum_{\tau=600}^{744-h}(\widehat y_{i,\tau+h}-y_{i,\tau+h})^2\right)^{1/2},
\end{equation}
where $y_{i,t}$ is the observed PM$_{2.5}$ data at Daliao district. From Table~\ref{Table4.22}, we see that our nonstationary factors perform as good as the ones in \cite{bai2004} using
the FECM of \cite{banerjee-etal2014}, and our factors perform slightly better when we adopt
all 3 detected unit-root common factors.
In conclusion, for the PM$_{2.5}$ data in Taiwan, the forecasting performance of the extracted factors seems to favor the proposed approach when we adopt the factor-based modeling and make predictions using the associated loadings. The nonstationary factors extracted by our method also fare well in comparison with the one by \cite{bai2004} when using the FECM approach in \cite{banerjee-etal2014}. Finally, we emphasize that the proposed model is different from the ones in \cite{bai2004} and \cite{StockWatson_2002a} because it explores a different aspect of the data. The proposed method is intended as another option for modeling high-dimensional unit-root time series.
\begin{table}[h]
\caption{The 1-step to 4-step ahead root mean squared forecast errors. GT-BMM denotes the FECM using the factors of proposed method, Bai-BMM is the FECM approach using the factors in \cite{bai2004}.}
\label{Table4.22}
\begin{center}
\setlength{\abovecaptionskip}{0pt}
\setlength{\belowcaptionskip}{3pt}\footnotesize
\begin{tabular}{ccccccccccc}
\hline
&&&\multicolumn{2}{c}{1 factor}&&\multicolumn{2}{c}{2 factors}&&\multicolumn{2}{c}{3 factors}\\
\cline{1-2}\cline{4-5}\cline{7-8}\cline{10-11}
Step&Model&>-BMM&Bai-BMM&>-BMM&Bai-BMM&>-BMM&Bai-BMM\\
\cline{1-2}\cline{4-5}\cline{7-8}\cline{10-11}
1&FECM(1)&&9.08&9.07&&14.45&14.40&&14.28&15.48\\
2&FECM(1)&&13.32&13.26&&18.32&18.33&&18.97&20.58\\
3&FECM(1)&&15.79&15.70&&19.37&19.36&&20.30&21.20\\
4&FECM(1)&&17.22&17.15&&19.98&19.95&&20.75&21.33\\
\hline
\end{tabular}
\end{center}
\end{table}
\section{Concluding Remarks}
This paper proposed a new approach to analyze high-dimensional unit-root time series from a factor modeling perspective. The proposed approach marks an extension of the work of \cite{zhang-etal2019}, \cite{gaotsay2019a,gaotsay2019b}, and
\cite{penaponcela2006}, and is in line with the frameworks of \cite{TiaoTsay_1989}
and \cite{BoxTiao_1977}.
Our approach uses modified methods of \cite{zhang-etal2019} to separate
the unit-root series from the stationary components via an eigenanalysis of certain
nonnegative definite matrix, which combines variance and lagged autocovariance matrices together.
To further reduce the dimension of the stationary components, we extend the
method of \cite{gaotsay2019b} with improvements in high-dimensional white noise
testing. The proposed approach is easy to implement for high-dimensional time series
and empirical results show that
it can effectively extract the unit root series and stationary common factors from complex data.
In addition, the extracted common trends and stationary common factors are useful in
out of sample predictions. Thus, the proposed approach expands the toolkits for
analyzing high-dimensional unit-root time series.