Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
92,416 characters · 15 sections · 120 citation commands
Modeling High-Dimensional Unit-Root Time Series
{\sl Keywords}: Common factor, Cointegration, Eigenanalysis, Factor model, High-dimensional time series, Unit root.
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, johansen2002 for the first difficulty and TiaoTsay_1989, lutkepohl2006, 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 TiaoTsay_1989, the LASSO regularization in VAR models by ShojaieMichailidis_2010 and SongBickel_2011, the sparse VAR model via partial spectral coherence in Davis2012, and the factor modeling by BaiNg_Econometrica_2002, StockWatson_2002a, StockWatson_2005, forni2005 and 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 engle-granger1987, johansen1988, johansen1991, 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 engle-granger1987 and the popular likelihood ratio (LR) test in johansen1988,johansen1991 with a parametric integrated VAR setting, to saikkonen2000 and aznar2002. As discussed in 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, engel-etal2015 contemplated the possibility of determining the cointegration rank of a system of seventeen OECD exchange rates. 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 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 zhang-etal2019, but in agreement with 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 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, 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, 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 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 zhang-etal2019 and gaotsay2019b, and is in line with the framework of 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 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). 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 baili2012. Instead, similar to po1988, robinson2002, penaponcela2006, and 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 bai2004 and the fixed lag sample covariance or autocovariances used in 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 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 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 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 zhang-etal2019 and the error-correction factor models of 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)$.
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
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 bai2004, and $r_2$ can be either fixed or slowly growing with the dimension $p$, which extends the results in 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, 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)) is general and in line with the framework of TiaoTsay_1989. Under the scalar component models of 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)) 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)) is a transformation model employed in TiaoTsay_1989.
Another way to see the generality of Model ((ref)) is to employ the canonical correlation analysis of 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 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 BoxTiao_1977 and TiaoTsay_1989 for general discussions.
To illustrate the identification issue of Model ((ref)) 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)) as
where ${\mathbf g}_{2t}$ is a $(p-r_1)$-dimensional stationary process. Note that Model ((ref)) 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
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
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)) to be orthonormal, we can show that ${\mathbf H}$ is a block-orthonormal matrix. Thus, Model ((ref)) 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
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)). Given Models ((ref)) and ((ref)), 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}$.
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
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
for some constant $0< C<\infty$. See Theorems 1 and 2 of 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$.
Similar to several works in the literature, e.g., lamyao2012, we combine the information over different lags of ${\mathbf y}_t$ and define
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)) is nonnegative definite, and there is no information cancellation over different lags. Unlike bai2004 and penaponcela2006, $\widehat{\mathbf M}_1$ of ((ref)) 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
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 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)).
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)) can be written as
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 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
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)), 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
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
When $p$ is large, similar as that in 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)).
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 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$.
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)), 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$. 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
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)) 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)) 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)) is the same as (2.2) in 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)) 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.
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 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
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.
We begin with some assumptions.
The mixing conditions in Assumptions 1(i) and 2 are standard for dependent random processes. See 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)) together with Assumption 1(i) by Theorem 1 of 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 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}$.
Similar results are also obtained in 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 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 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.
{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.}
Assumption 4 is similar to Assumption 5 of 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.
From Theorem (ref), 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 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 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). The next theorem states the consistency of the estimated $\widehat r_1$ and $\widehat r_2$ with additional constraints on $p$.
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)) 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
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)). We only present the simulation results for $k_0=2$ and $j_0=2$ in Equations ((ref)) and ((ref)), respectively, to save space because other choices of $k_0$ and $j_0$ produce similar results.\\
{\bf Example 1.} Consider the models in ((ref)) and ((ref)) 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 1234 in 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 hoff2009 to simulate a random orthonormal matrix ${\mathbf A}$ of Model ((ref)) 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)) 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) 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)(a), (b) and (c), respectively, where $\bar{D}(\cdot,\cdot)$ is defined in ((ref)). From Figure (ref), 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.
Next, for each $(p,n)$ configuration, we study the root-mean-square errors (RMSEs):
which quantify the accuracy in estimating the common factor processes. Boxplots of the RMSE$_1$ and RMSE$_2$ are shown in Figure (ref)(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). 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$.
{\bf Example 2.} In this example, we consider Models ((ref)) and ((ref)) 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)) to estimate the number of unit-root series and that of high-dimensional white-noise tests of 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). 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.
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 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)) and we choose $R=(p-\widehat r_1)/2$ suggested in their paper. Figure (ref) 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).
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)(a)--(c), respectively. Similar patterns are also obtained for the estimation of other matrices so we omit them here. From Figure (ref), 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$.
For each $(p,n)$, we further define the RMSEs for high dimension as:{
} 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)(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). Overall, the proposed method works well even for moderately high dimensions. This is especially so when the sample size is greater than the dimension.
{\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), 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).
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, 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). From the plots, we see that these three unit-root factors capture most of the trends in the original data of Figure (ref) and, as expected, their sample ACFs decay slowly. The extracted 505-dimensional stationary process is shown in Figure (ref)(a).
Next, we applied the white noise testing method of Section 2.3 with $k_0 = j_0 = 2$ for Equations ((ref)) and ((ref)) 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 gaotsay2019b, we first examine the eigenvalues of the sample covariance matrix $\widehat{\mathbf S}$ defined in ((ref)). From Figure (ref)(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)(b) and (c), respectively. From Figure (ref), 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.
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 BaiNg_Econometrica_2002 and the nonstationary one in bai2004, and found that the numbers of common factors are 20 and 1, respectively. We also set the number of nonstationary factors in 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 bai2004 and a VAR(1) model to the 20-dimensional factors extracted by 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
where $p=508$. Table (ref) reports the 1-step to 4-step ahead forecast errors of Equation ((ref)) 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 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). From the plot, we see that the proposed approach outperforms the other four methods in most time periods of the forecasting subsample.
We then compare the forecast ability of extracted factors by different approaches and the benchmark method by adopting the asymptotic test in 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) reports the testing results, where we use the method in 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 diebold-1995. From Table (ref), 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) and (ref) 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.
To further compare the forecast ability of the nonstationary factors extracted by the proposed method and the one in bai2004, we employ the factor-augmented error correction models (FECM) in banerjee-etal2014. As the factors in banerjee-etal2014 are extracted by the method of bai2004, we denote the approach of FECM of banerjee-etal2014 by Bai-BMM, and our method by GT-BMM. Since the method in bai2004 identifies 1 factor and the proposed method estimates 3, we include 1 to 3 factors in the FECM defined in Equation (8) of 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) presents the 1-step to 4-step ahead root mean squared forecast errors (RMSFE) defined as
where $y_{i,t}$ is the observed PM$_{2.5}$ data at Daliao district. From Table (ref), we see that our nonstationary factors perform as good as the ones in bai2004 using the FECM of 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 bai2004 when using the FECM approach in banerjee-etal2014. Finally, we emphasize that the proposed model is different from the ones in bai2004 and 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.
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 zhang-etal2019, gaotsay2019a,gaotsay2019b, and penaponcela2006, and is in line with the frameworks of TiaoTsay_1989 and BoxTiao_1977. Our approach uses modified methods of 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 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.