Factor Analysis for High-Dimensional Time Series with Change Point

TL;DR

This paper develops methods for detecting and estimating change points in high-dimensional time series with latent factors, allowing for strong cross-sectional dependence and using self-normalization for robust testing.

Contribution

It introduces consistent estimators for factor loadings and change points in high-dimensional series with dependence, advancing change-point analysis in complex data.

Findings

Proposed estimators are consistent under strong dependence.

Self-normalization improves change-point test robustness.

Numerical experiments validate the methods.

Abstract

We consider change-point latent factor models for high-dimensional time series, where a structural break may exist in the underlying factor structure. In particular, we propose consistent estimators for factor loading spaces before and after the change point, and the problem of estimating the change-point location is also considered. Compared with existing results on change-point factor analysis of high-dimensional time series, a distinguished feature of the current paper is that our results allow strong cross-sectional dependence in the noise process. To accommodate the unknown degree of cross-sectional dependence strength, we propose to use self-normalization to pivotalize the change-point test statistic. Numerical experiments including a Monte Carlo simulation study and a real data application are presented to illustrate the proposed methods.

Tables

Table 1. Table 1: Empirical sizes of tests with different levels of factor strength in Section 4.1

		$\alpha =10\%$		$\alpha =5\%$		$\alpha =1\%$
$\delta$		0	0.25	0	0.25	0	0.25
n=400	p=20	0.112	0.138	0.064	0.072	0.014	0.012
n=400	p=40	0.127	0.149	0.069	0.083	0.013	0.020
n=400	p=100	0.123	0.112	0.067	0.059	0.018	0.016
n=1000	p=20	0.116	0.122	0.058	0.057	0.009	0.011
n=1000	p=40	0.103	0.130	0.058	0.072	0.012	0.012
n=1000	p=100	0.115	0.137	0.068	0.071	0.014	0.019

Table 2. Table 2: Empirical powers of tests under different settings in Section 4.2

		$\alpha =10\%$				$\alpha =5\%$				$\alpha =1\%$
Setting		SS	SW	WS	WW	SS	SW	WS	WW	SS	SW	WS	WW
n=400	p=20	0.975	0.973	0.966	0.969	0.941	0.929	0.932	0.943	0.820	0.787	0.815	0.810
n=400	p=40	0.979	0.970	0.955	0.969	0.944	0.936	0.915	0.934	0.820	0.804	0.783	0.799
n=400	p=100	0.960	0.968	0.975	0.977	0.931	0.924	0.949	0.940	0.805	0.806	0.811	0.821
n=1000	p=20	0.996	0.996	0.998	0.996	0.990	0.989	0.992	0.992	0.956	0.949	0.958	0.959
n=1000	p=40	0.998	0.996	0.994	0.997	0.989	0.987	0.985	0.987	0.935	0.956	0.946	0.950
n=1000	p=100	0.999	0.990	0.997	0.999	0.996	0.995	0.992	0.992	0.960	0.962	0.960	0.962

Table 3. Table 3: Average estimation error | γ ^ − γ 0 | ^ 𝛾 subscript 𝛾 0 |\widehat{\gamma}-\gamma_{0}| in Section 4.2

$n$	$n=400$			$n=1000$
$p$	$20$	$40$	$100$	$20$	$40$	$100$
SS	0.035	0.039	0.040	0.015	0.018	0.018
SW	0.051	0.066	0.083	0.023	0.029	0.039
WS	0.054	0.060	0.081	0.023	0.027	0.038
WW	0.053	0.060	0.071	0.024	0.028	0.034

Table 4. Table 4: Average estimation error 𝒟 ​ { ℳ ​ [ 𝐐 ^ i ​ ( r ^ ) ] , ℳ ​ ( 𝐐 i ) } 𝒟 ℳ delimited-[] subscript ^ 𝐐 𝑖 ^ 𝑟 ℳ subscript 𝐐 𝑖 \mbox{$\mathcal{D}$}\{{\cal M}[\widehat{{\mathbf{Q}}}_{i}(\widehat{r})],\,{\cal M}({\mathbf{Q}}_{i})\} in Section 4.2

	$n$	$n=400$			$n=1000$
	$p$	$20$	$40$	$100$	$20$	$40$	$100$
Setting SS	${\delta }_1=0$	0.059	0.055	0.053	0.033	0.032	0.031
	${\delta }_2=0$	0.059	0.057	0.055	0.033	0.032	0.031
Setting SW	${\delta }_1=0$	0.058	0.056	0.053	0.033	0.032	0.031
	${\delta }_2=0.25$	0.095	0.099	0.113	0.049	0.053	0.058
Setting WS	${\delta }_1=0.25$	0.093	0.094	0.110	0.048	0.052	0.058
	${\delta }_2=0$	0.060	0.056	0.053	0.034	0.032	0.032
Setting WW	${\delta }_1=0.25$	0.089	0.095	0.103	0.049	0.052	0.057
	${\delta }_2=0.25$	0.093	0.094	0.105	0.050	0.052	0.057

Table 5. Table 5: Size comparison among different tests with α = 5 % 𝛼 percent 5 \alpha=5\% for Example 1 in Section 4.3

$\delta$		Our method		Wald-cdg		LM-hi		Wald-hi
$n$	$p$	0	0.25	0	0.25	0	0.25	0	0.25
$100$	$20$	0.086	0.100	0.238	0.203	0.016	0.017	0.118	0.130
$100$	$40$	0.110	0.148	0.267	0.203	0.004	0.010	0.161	0.136
$100$	$100$	0.119	0.204	0.235	0.210	0.006	0.008	0.161	0.141
$400$	$20$	0.074	0.071	0.072	0.092	0.028	0.023	0.045	0.042
$400$	$40$	0.086	0.073	0.072	0.092	0.029	0.031	0.061	0.045
$400$	$100$	0.099	0.077	0.084	0.093	0.023	0.026	0.050	0.050

Table 6. Table 7: Comparison among different models for real data analysis

Method	$p$-value for white noise test	RSS
Ours	0.075	308238
Wald-cdg	0.000	160476
LM-hi	0.070	1482874
Wald-hi	0.030	1156676

Equations

$${\mathbf{y}}_{t}={\mathbf{A}}{\mathbf{x}}_{t}+\mbox{\boldmath$\varepsilon$}_{t},$$

$${\mathbf{y}}_{t}=\left\{\begin{array}[]{ll}{\mathbf{A}}_{1}{\mathbf{x}}_{t,1}+\mbox{\boldmath$\varepsilon$}_{t},&\mathrm{if}\ t\leq r_{0};\\ {\mathbf{A}}_{2}{\mathbf{x}}_{t,2}+\mbox{\boldmath$\varepsilon$}_{t},&\mathrm{if}\ t>r_{0},\end{array}\right.$$

$$\mathcal{D}\{\mathcal{M}({\mathbf{S}}_1),\mathcal{M}({\mathbf{S}}_2)\}={\left\{1-\frac{\mathrm{tr}({\mathbf{O}}_1{\mathbf{O}}_1^{\prime }{\mathbf{O}}_2{\mathbf{O}}_2^{\prime })}{\min (q_1,q_2)}\right\}}^{1/2}.$$

$$\mathcal{D}\{\mathcal{M}({\mathbf{S}}_1),\mathcal{M}({\mathbf{S}}_2)\}={\left\{1-\frac{\mathrm{tr}({\mathbf{O}}_1{\mathbf{O}}_1^{\prime }{\mathbf{O}}_2{\mathbf{O}}_2^{\prime })}{q}\right\}}^{1/2},$$

$$\Vert {\mathbf{A}}_i{\Vert }_2^2\asymp \Vert {\mathbf{A}}_i{\Vert }_{\min }^2\asymp p^{1-{\delta }_i}$$

$${\mathbf{\Sigma }}_{x,i}(h,\gamma )=\frac{1}{n}\sum _{t=1}^n\mathrm{E}\left[{\mathbf{x}}_{t,i}{\mathbf{x}}_{t+h,i}^{\prime }{I}_{t,i}(\gamma )I_{t+h,i}(\gamma )\right],{\mathbf{\Sigma }}_{y,i}(h,\gamma )=\frac{1}{n}\sum _{t=1}^n\mathrm{E}\left[{\mathbf{y}}_t{\mathbf{y}}_{t+h}^{\prime }{I}_{t,i}(\gamma )I_{t+h,i}(\gamma )\right],$$

$${\mathbf{\Sigma }}_{y,1}(h,\gamma )={\mathbf{A}}_1{\mathbf{\Sigma }}_{x,1}(h,\gamma ){\mathbf{A}}_1^{\prime },\text{mbox}if\gamma \le {\gamma }_0,$$

$${\mathbf{\Sigma }}_{y,2}(h,\gamma )={\mathbf{A}}_2{\mathbf{\Sigma }}_{x,2}(h,\gamma ){\mathbf{A}}_2^{\prime },\text{mbox}if\gamma >{\gamma }_0.$$

$${\mathbf{M}}_i(\gamma )=\sum _{h=1}^{h_0}{\mathbf{\Sigma }}_{y,i}(h,\gamma ){\mathbf{\Sigma }}_{y,i}(h,\gamma {)}^{\prime },$$

$${\mathbf{M}}_i(\gamma )=\sum _{h=1}^{h_0}{\mathbf{A}}_i\left[{\mathbf{\Sigma }}_{x,i}(h,\gamma ){\mathbf{A}}_i^{\prime }{\mathbf{A}}_i{\mathbf{\Sigma }}_{x,i}(h,\gamma {)}^{\prime }\right]{\mathbf{A}}_i^{\prime }$$

$${\mathbf{Q}}_i(\gamma )=({\mathbf{q}}_{i,1}(\gamma ),\dots ,{\mathbf{q}}_{i,k_i}(\gamma )),$$

$${\hat{\mathbf{\Sigma }}}_{y,i}(h,\gamma )=\frac{1}{n}\sum _{t=1}^n\left({\mathbf{y}}_t{\mathbf{y}}_{t+h}^{\prime }{I}_{t,i}(\gamma )I_{t+h,i}(\gamma )\right),{\hat{\mathbf{M}}}_i(\gamma )=\sum _{h=1}^{h_0}{\hat{\mathbf{\Sigma }}}_{y,i}(h,\gamma ){\hat{\mathbf{\Sigma }}}_{y,i}(h,\gamma {)}^{\prime },$$

$${\hat{\mathbf{Q}}}_i(\gamma )=({\hat{\mathbf{q}}}_{i,1}(\gamma ),\dots ,{\hat{\mathbf{q}}}_{i,k_i}(\gamma )),$$

$$\mathcal{D}\{\mathcal{M}[{\hat{\mathbf{Q}}}_1(\gamma )],\mathcal{M}({\mathbf{Q}}_1)\}=O_p(p^{{\delta }_1}{n}^{-1/2})$$

$$\mathcal{D}\{\mathcal{M}[{\hat{\mathbf{Q}}}_2(\gamma )],\mathcal{M}({\mathbf{Q}}_2)\}=O_p(p^{{\delta }_2}{n}^{-1/2})$$

$$G(\gamma)=\sum_{i=1}^{2}g_{i}(\gamma),\quad g_{i}(\gamma)=\Big{\|}{{\mathbf{B}}_{i}(\eta_{i})}^{\prime}\,{\mathbf{M}}_{i}(\gamma)\,{\mathbf{B}}_{i}(\eta_{i})\Big{\|}_{2},$$

$$G(\gamma_{0})=\sum_{i=1}^{2}\Big{\|}{{\mathbf{B}}_{i}}^{\prime}\,{\mathbf{M}}_{i}(\gamma_{0})\,{\mathbf{B}}_{i}\Big{\|}_{2}=\sum_{i=1}^{2}\Big{\|}\sum_{h=1}^{h_{0}}\left\{{{\mathbf{B}}_{i}}^{\prime}\,{\mathbf{A}}_{i}\left[\boldsymbol{\Sigma}_{x,i}(h,\gamma_{0}){\mathbf{A}}_{i}^{\prime}{\mathbf{A}}_{i}\boldsymbol{\Sigma}_{x,i}(h,\gamma_{0})\right]{\mathbf{A}}_{i}^{\prime}\,{\mathbf{B}}_{i}\right\}\Big{\|}_{2}=0.$$

$${\hat{\mathbf{B}}}_i({\eta }_i)=({\hat{\mathbf{q}}}_{i,k_i+1}({\eta }_i),\dots ,{\hat{\mathbf{q}}}_{i,p}({\eta }_i)),$$

$$\widehat{G}(\gamma)=\sum_{i=1}^{2}\Big{\|}\widehat{{\mathbf{B}}}_{i}(\eta_{i})^{\prime}\,\widehat{{\mathbf{M}}}_{i}(\gamma)\,\widehat{{\mathbf{B}}}_{i}(\eta_{i})\Big{\|}_{2}.$$

$$\hat{\gamma }={\mathrm{argmin}}_{\gamma \in \{0,\frac{1}{n},\dots ,1\}\cap ({\eta }_1,{\eta }_2)}\hat{G}(\gamma ),$$

$$P(\hat{\gamma }<{\gamma }_0-ϵ)\le \frac{C{p}^{{\delta }_1}}{ϵn^{1/2}},P(\hat{\gamma }>{\gamma }_0+ϵ)\le \frac{C{p}^{{\delta }_2}}{ϵn^{1/2}},$$

$${\hat{\mathbf{Q}}}_i(\hat{r})=({\hat{\mathbf{q}}}_{i,1}(\hat{r}),\dots ,{\hat{\mathbf{q}}}_{i,k_i}(\hat{r})),$$

$$\mathcal{D}\{\mathcal{M}[{\hat{\mathbf{Q}}}_i(\hat{\gamma })],\mathcal{M}({\mathbf{Q}}_i)\}=O_p(p^{{\delta }_i}{n}^{-1/2})$$

$${\hat{k}}_i={\mathrm{argmin}}_{1\le k\le p/2}\frac{{\hat{\lambda }}_{i,k+1}({\eta }_i)}{{\hat{\lambda }}_{i,k}({\eta }_i)},$$

$$\mathrm{var}({\mathbf{b}}^{\prime }{\mathbf{y}}_t)={\mathbf{b}}^{\prime }\mathbf{\Sigma }\mathbf{b},t\le r_0,$$

$$\mathrm{var}({\mathbf{b}}^{\prime}{\mathbf{y}}_{t})=\mathrm{var}({\mathbf{b}}^{\prime}{\mathbf{A}}_{2}{\mathbf{x}}_{t,2}+{\mathbf{b}}^{\prime}\mbox{\boldmath$\varepsilon$}_{t})={\mathbf{b}}^{\prime}{\mathbf{A}}_{2}\mathrm{cov}({\mathbf{x}}_{t,2}){\mathbf{A}}_{2}^{\prime}{\mathbf{b}}+{\mathbf{b}}^{\prime}\boldsymbol{\Sigma}{\mathbf{b}},\quad t>r_{0}.$$

$$H_0:\mathrm{var}({\mathbf{b}}^{\prime }{\mathbf{y}}_1)=\cdots =\mathrm{var}({\mathbf{b}}^{\prime }{\mathbf{y}}_n)$$

$${\hat{\nu }}_{i,j}=(j-i{)}^{-1}\left(\sum _{t=i}^j{\mathbf{b}}^{\prime }{\mathbf{y}}_t{\mathbf{y}}_t^{\prime }\mathbf{b}-\sum _{t=i}^j{\mathbf{b}}^{\prime }{\mathbf{y}}_t\sum _{t=i}^j{\mathbf{y}}_t^{\prime }\mathbf{b}\right),$$

$$T_n=\underset{n{\eta }_1<r<n{\eta }_2}{\max }\frac{\{r(n-r)({\hat{\nu }}_{1,r}-{\hat{\nu }}_{r+1,n}){\}}^2}{n^2{V}_r},$$

$$V_r=\frac{1}{n}\left[\sum _{i=1}^r{\left\{\frac{i(r-i)({\hat{\nu }}_{1,i}-{\hat{\nu }}_{i+1,r})}{r}\right\}}^2+\sum _{i=r+1}^n{\left\{\frac{(i-r-1)(n-i+1)({\hat{\nu }}_{r+1,i-1}-{\hat{\nu }}_{i,n})}{n-r}\right\}}^2\right]$$

Full text

KEYWORDS: Factor model; high-dimensional time series; change point detection; strong cross-sectional dependence.

1 Introduction

High-dimensional time series has been emerging as a common and important data type in applications from a number of disciplines, including climate science, economics, finance, medical science, and telecommunication engineering among others. Although numerous statistical methods and their associated theory have been developed for the modeling and inference of time series data, existing results mostly focused on the univariate or finite-dimensional multivariate case. The problem of extending existing results developed under low-dimensional settings to handle high-dimensional time series, however, is typically nontrivial and requires significant innovations. For example, when the dimension is larger than the length of the observed time series, the commonly used autoregressive moving-average (ARMA) model in its conventional form may face a serious identification problem as commented by Lam et al. (2011). To handle the phenomenon of high dimensionality, one typically resorts to certain sparsity-type conditions for the purpose of dimension reduction. For example, when considering vector autoregressive (VAR) models in the high-dimensional setting, one typically need to assume that the coefficient matrices are sparse in a suitable sense in order to obtain their meaningful estimators; see for example Basu and Michailidis (2015), Davis et al. (2016) and references therein for research in this direction.

Unlike the aforementioned sparse VAR approach that aims at extending existing parametric time series models to their sparse high-dimensional counterparts, a more commonly used approach in the literature for modeling high-dimensional time series is through the use of factor models; see for example Chamberlain and Rothschild (1983) and Stock and Watson (1998), Bai and Ng (2002), and Forni et al. (2004) among others. In the aforementioned works, it was assumed that most of the variations in the observed high-dimensional time series can be explained by a set of common factors, and as a result such factor models cannot be used to capture strong cross-sectional dependence. In addition, the common component suffers from the identifiability issue when the dimension is finite; see for example Bai and Ng (2002). To alleviate these problems, Lam et al. (2011) proposed an alternative type of factor models which has become more and more popular in the last decade. In the model of Lam et al. (2011), the common factors are viewed as the force that drives all the dynamics and are used to explain the serial dependence in the data. Under this setting, the noise process can accommodate the strong cross-sectional dependence and is white. Lam et al. (2011) proposed an approach based on autocovariance matrices of the observed process at nonzero lags for loading space estimation. This method is also applicable to non-stationary processes, processes with uncorrelated or endogenous regressors, and matrix-valued process; see for example Liu and Chen (2016), Chang et al. (2017), Wang et al. (2019), and Liu and Chen (2019). Compared with the approach of Bai and Ng (2002) and Forni et al. (2004), the factor model by Lam et al. (2011) not only ensures that the common component is identifiable but also captures all serial dependence of data which enables us to build forecasting models after dimension reduction. However, existing results along the line of Lam et al. (2011) were generally developed under the assumption that the underlying factor structure remains the same over the whole time period, while recent empirical applications reveal that the factor structure tends to exhibit structural breaks at a certain point, either due to a sudden market change or in response to some unpredictable reasons. This motivates us to consider factor models with possible change points.

The problem of change-point analysis for factor models has been an active area of research. For this, Breitung and Eickmeier (2011) considered testing the existence of change in loadings when the noise process is either independent or autoregressive with a finite order. Chen et al. (2014) came up with a Lagrange multiplier test and a Wald test for detecting the break by regressing one of the factors estimated by PCA on the remaining estimated factors. Han and Inoue (2015) testes the break point through the second moments of the estimated factors. Barigozzi et al. (2018) constructed a cumulative sum (CUSUM) test using wavelet coefficients for detecting multiple structural breaks when the noise sequence is Gaussian. Besides testing-based methods, Chen (2015) and Baltagi et al. (2017) considered using methods based on least squares to estimate the change-point location. Ma and Su (2018) proposed an adaptive fused group lasso approach to estimate the multiple change-point location. However, existing results in this direction were typically developed for models with the so-called idiosyncratic noise, where no strong cross-sectional dependence is allowed in the noise. In the current paper, we follow settings by Lam et al. (2011) where the common factor drives all the dynamics and explain serial dependence of the observed process, and consider the situation when the strength of cross-sectional dependence of the noise sequence is unknown and potentially strong. As discussed above, one advantage of this setting is that we can fully extract the dependence of data through common factors for future forecasting model building if needed. The aforementioned paper considered the problem of estimating the factor loading by assuming that there is no change point, and we shall here focus on the change-point case.

The remaining of the paper is organized as follows. Section 2 introduces the change-point factor model and considers its associated estimation problem, including estimating factor loading spaces before and after the change point and the location of the change point. The asymptotic properties of the proposed estimators are also investigated. Section 3.1 proposes a self-normalized approach to testing the existence of change point, which has a pivotalized asymptotic distribution regardless of whether the noise process has weak or strong cross-sectional dependence. Details on its practical implementation are discussed in Section 3.2. Numerical experiments, including a Monte Carlo simulation study and a real data application, are presented in Sections 4 and 5 respectively to illustrate the proposed methods. Section 6 concludes the paper. Technical proofs are deferred to the Appendix.

2 Change-Point Factor Model and Its Estimation

In this section, we will introduce the our change-point factor model and study its associated estimation problem, including the loading space before and after the change point and the change point location.

2.1 Change-Point Factor Model

Suppose we observe a $p$-dimensional time series ${\mathbf{y}}_{t}$, $t=1,\ldots,n$, according to a factor model, then one can write

[TABLE]

where $({\mathbf{x}}_{t})$ is a latent factor process whose dimension $k_{0}$ is typically much smaller than $p$, ${\mathbf{A}}\in\mathbb{R}^{p\times k_{0}}$ is the associated loading matrix, and (\mbox{\boldmath\varepsilon}_{t}) denotes the white noise process. The latent factor model (1) has been widely used in the literature for dimension reduction of high-dimensional time series; see for example Bai and Ng (2002), Lam et al. (2011), Lam and Yao (2012), Chang et al. (2015) and references therein. It also relates to the generalized dynamic factor model of Forni et al. (2005) in which the latent factor process $({\mathbf{x}}_{t})$ is assumed to follow a low-dimensional autoregressive model. We shall here consider the general setting where $({\mathbf{x}}_{t})$ is not necessarily an autoregressive process. In model (1), the factor loading structure remains the same over the whole sampling period, while in many applications a structural break may occur due to various reasons. For this, we consider the change-point factor model

[TABLE]

where ${\mathbf{x}}_{t,i}\in\mathbb{R}^{k_{i}}$, $i=1,2$, represents the underlying latent factor before and after the change point whose location is denoted by $r_{0}$, and ${\mathbf{A}}_{1}$ and ${\mathbf{A}}_{2}$ are the associated loading matrices. Recently, there have been efforts in studying the change-point factor model by incorporating certain beliefs on the change-point mechanism into the analysis. For example, Liu and Chen (2016) modeled the change-point mechanism by a finite-state hidden Markov chain, and thus structural breaks occur when there is a regime switching in the hidden state variable. On the other hand, Liu and Chen (2019) considered using a threshold variable to model the change-point mechanism, where the threshold variable is assumed to be $\alpha$-mixing and observable up to a small number of unknown parameters. Instead of introducing a Markov chain process or an additional threshold variable, we in the current paper focus on the change-point factor model (2) which uses time to naturally divide the observed process into homogenous pieces before and after the change point.

In practice, one does not observe $({\mathbf{x}}_{t,i})$ nor (\mbox{\boldmath\varepsilon}_{t}) but only $({\mathbf{y}}_{t})$, and thus the loading matrices in (2) are not fully identifiable. To be more specific, for $i=1,2$, one can replace $({\mathbf{A}}_{i},{\mathbf{x}}_{t,i})$ in (2) by $({\mathbf{A}}_{i}{\mathbf{U}}_{i},{\mathbf{U}}^{-1}_{i}{\mathbf{x}}_{t})$ for any $k_{i}\times k_{i}$ non-singular matrix ${\mathbf{U}}_{i}$. However, the space spanned by the columns of the loading matrix ${\mathbf{A}}_{i}$, denoted by ${\cal M}({\mathbf{A}}_{i})$, is always uniquely defined; see for example the discussions in Lam et al. (2011), Chang et al. (2015), Liu and Chen (2016) and Liu and Chen (2019). As a result, when estimating the factor model (2), we shall focus on the loading space ${\cal M}({\mathbf{A}}_{i})$ instead of the loading matrix ${\mathbf{A}}_{i}$ itself.

We shall here introduce some notations. For a matrix ${\mathbf{H}}$, we use $\|{\mathbf{H}}\|_{F}$ and $\|{\mathbf{H}}\|_{2}$ to denote its Frobenius and L-2 norms respectively. In addition, we use ${\rm tr}({\mathbf{H}})$ for the trace, $\sigma_{i}({\mathbf{H}})$ for the $i$-th largest singular value, and $\|{\mathbf{H}}\|_{\min}$ for the square root of minimum nonzero eigenvalue of ${\mathbf{H}}^{\prime}{\mathbf{H}}$. Also, we write $a\asymp b$ if $a=O(b)$ and $b=O(a)$, and we use $\lfloor{x}\rfloor$ and $\lceil{x}\rceil$ to denote the largest previous and smallest following integers of $x$.

2.2 Estimating the Loading Spaces

Before introducing our estimators for loading spaces in the change-point setting, we first need a measure to quantify the distance between two linear spaces which can then be used to assess the statistical performance our estimators. For this, let ${\mathbf{S}}_{1}$ and ${\mathbf{S}}_{2}$ be full rank matrices in $\mathbb{R}^{p\times q_{1}}$ and $\mathbb{R}^{p\times q_{2}}$ respectively with $\max(q_{1},q_{2})\leq p$. Denote ${\mathbf{O}}_{i}$ the matrix whose columns form an orthonormal basis of ${\cal M}({\mathbf{S}}_{i})$ for $i=1,2$, then the distance between column spaces of ${\mathbf{S}}_{1}$ and ${\mathbf{S}}_{2}$ can be measured by

[TABLE]

The distance measure (3) was first introduced in Liu and Chen (2019), and is a quantity between 0 and 1. In particular, it equals to 0 if ${\cal M}({\mathbf{S}}_{1})\in{\cal M}({\mathbf{S}}_{2})$ or ${\cal M}({\mathbf{S}}_{2})\in{\cal M}({\mathbf{S}}_{1})$, and equals to 1 if ${\cal M}({\mathbf{S}}_{1})$ and ${\cal M}({\mathbf{S}}_{2})$ are orthogonal. For the special case when $q_{1}=q_{2}=q$, the two spaces ${\mathbf{S}}_{1}$ and ${\mathbf{S}}_{2}$ have the same dimension, and the distance measure (3) reduces to

[TABLE]

which was used in Chang et al. (2015) and Liu and Chen (2016). Since the number of factors is usually unknown in practice and may be estimated in a nonperfect way, we shall in the current paper use the generalized version in (3) to measure the distance between two linear spaces.

For factor models in high-dimensional cases, it is common to assume that the number of factors is fixed but the squared L-2 norm of the $p\times k_{i}$ loading matrix ${\mathbf{A}}_{i}$ grows with the dimension $p$ (Bai and Ng, 2002; Doz et al., 2011). The growth rate is called the strength of the factors in Lam et al. (2011), Lam and Yao (2012), Chang et al. (2015), and Liu and Chen (2019+). For $i=1,2$, assume that

[TABLE]

for some $0\leq\delta_{i}\leq 1$, then factors in regime $i$ are said to be strong if $\delta_{i}=0$ and weak if $\delta_{i}\in(0,1)$. The strength of factors measures the relative growth rate of the amount of information carried by the observed process ${\mathbf{y}}_{t}$ about the common factors ${\mathbf{x}}_{t}$ as $p$ increases, with respect to the growth rate of the amount of noise process in regime $i$. It can be seen from our theoretical results below that the factor strength plays an important role in the estimation efficiency.

Let $\gamma_{0}=r_{0}/n$, and we shall now introduce the estimation procedure for loading spaces when a tentative break date $\gamma$ is given. For this, we separate the data into two subsets before and after $\gamma$, namely $t\leq\lfloor\gamma n\rfloor$ and $t>\lfloor\gamma n\rfloor$, and propose a change-point generalization of the estimator of Lam et al. (2011). To be more specific, define the generalized second cross moment matrices of ${\mathbf{x}}_{t,i}$ and ${\mathbf{y}}_{t}$ with lag $h$ in each regime $i$ as

[TABLE]

where $I_{t,i}(h,\gamma)$ is the indicator function for regime $i$ satisfying $I_{t,1}(\gamma)=1$ if $1\leq t\leq\lfloor\gamma n\rfloor$, $I_{t,2}(h,\gamma)=1$ if $\lfloor\gamma n\rfloor<t\leq n$, and zero otherwise. With the white noise assumption, $\boldsymbol{\Sigma}_{x,i}(h,\gamma)$ and $\boldsymbol{\Sigma}_{y,i}(h,\gamma)$ satisfy

[TABLE]

and

[TABLE]

For a pre-determined positive integer $h_{0}$, let

[TABLE]

then we can see that

[TABLE]

holds for $i=1$ if $\gamma\leq\gamma_{0}$ and for $i=2$ if $\gamma\geq\gamma_{0}$. Therefore, ${\mathbf{M}}_{i}(\gamma)$ is a symmetric non-negative definite matrix sandwiched by ${\mathbf{A}}_{i}$ and ${\mathbf{A}}_{i}^{\prime}$ when $\gamma$ is in regime $i$. If there exists at least one nonzero $h$ such that $1\leq h\leq h_{0}$ and $\boldsymbol{\Sigma}_{x,i}(h,\gamma)$ is full rank, then ${\mathbf{M}}_{i}(\gamma)$ has rank $k_{i}$. Its eigenspace corresponding to the nonzero eigenvalues is ${\cal M}({\mathbf{A}}_{i})$. Hence, ${\cal M}({\mathbf{A}}_{i})$ can be estimated through an eigen-decomposition of the sample version of ${\mathbf{M}}_{i}$. In particular, let ${\mathbf{q}}_{i,k}(\gamma)$ be the unit eigenvector of ${\mathbf{M}}_{i}(\gamma)$ corresponding to the $k$-th largest nonzero eigenvalue for $k=1,\ldots,k_{i}$. Define

[TABLE]

where $k_{i}$ is the number of factors. For simplicity, in the rest of this paper we use ${\mathbf{Q}}_{i}$ and ${\mathbf{M}}_{i}$ to denote ${\mathbf{Q}}_{i}(\gamma_{0})$ and ${\mathbf{M}}_{i}(\gamma_{0})$ for $i=1,2$.

To estimate the loading spaces in the change-point setting, we shall use the sample version of the above quantities and define

[TABLE]

for $i=1,2$. Let $\widehat{\lambda}_{i,1}\geq\widehat{\lambda}_{i,2}\geq\ldots\geq\widehat{\lambda}_{i,p}$ be the $p$ eigenvalues of $\widehat{{\mathbf{M}}}_{i}(\gamma)$ and $\widehat{{\mathbf{q}}}_{i,1}(\gamma),\widehat{{\mathbf{q}}}_{i,2}(\gamma),\ldots,\widehat{{\mathbf{q}}}_{i,p}(\gamma)$ be the set of corresponding orthonormal eigenvectors with $\mathbf{1}^{\prime}\widehat{{\mathbf{q}}}_{i,j}(\gamma)>0$. Define

[TABLE]

then ${\cal M}({\mathbf{A}}_{i})$ can be estimated by ${\cal M}(\widehat{{\mathbf{Q}}}_{i}(\gamma))$. When $\gamma$ is in regime $i$, similar to Theorem 1 in Lam et al. (2011), we can show that $\widehat{{\mathbf{Q}}}_{i}(\gamma)$ provides a consistent estimator for the loading space under mild conditions.

Proposition 1**.**

Assume Conditions 1–8 in Appendix A.1. If $p^{\delta_{\max}}n^{-1/2}=o(1)$, then as $n,p\to\infty$, with true $k_{1}$ and $k_{2}$, we have

[TABLE]

when $\gamma\leq\gamma_{0}$, and

[TABLE]

when $\gamma\geq\gamma_{0}$, where $\delta_{\max}=\max(\delta_{1},\delta_{2})$.

We remark that when $\delta_{i}=0$, the estimator ${\cal M}[\widehat{{\mathbf{Q}}}_{i}(\gamma)]$ converges to ${\cal M}({\mathbf{Q}}_{i})$ at the rate of $n^{-1/2}$, and thus the curse of dimensionality does not exist. When the factors in regime $i$ are weak, however, the convergence rate is slower and the noise process distorts the information on the latent factor; see for example Lam et al. (2011).

We shall here also provide a discussion about the situation when $\gamma$ does not fall in regimes $i$. Without loss of generality, we illustrate by using the case when $\gamma\geq\gamma_{0}$. In this case, ${\mathbf{M}}_{2}(\gamma)$ is sandwiched by ${\mathbf{A}}_{2}$, and $\widehat{{\mathbf{Q}}}_{2}(\gamma)$ is a reasonable estimate for ${\cal M}({\mathbf{A}}_{2})$. However, misclassification does occur for the estimation of regime 1. In particular, since data points with $\lfloor\gamma n\rfloor<t\leq r_{0}$ from regime 2 are included in the calculation of ${{\mathbf{M}}}_{1}(\gamma)$, it is no longer sandwiched by ${\mathbf{A}}_{1}$. However, when $\gamma$ is sufficiently close $\gamma_{0}$, one can show that the misclassification effect becomes negligible on the estimation, and the estimated space using the sample version of ${\mathbf{M}}_{1}(\gamma)$ will continue to be consistent; see the results in Section 2.3.

2.3 Estimating the Change Point Location

In our estimation approach, we assume that the change point does not happen in the boundary area, namely there exists $0<\eta_{1}<\eta_{2}<1$ such that $\gamma_{0}\in(\eta_{1},\eta_{2})$. Let $\gamma$ be a hypothesized change point location, then we can use it to split the data into two subsets, namely the one with $t\leq\lfloor\gamma n\rfloor$ and the other with $t>\lfloor\gamma n\rfloor$, and we define

[TABLE]

where ${\mathbf{B}}_{i}$ is a $p\times(p-k_{i})$ matrix for which $({\mathbf{Q}}_{i},{\mathbf{B}}_{i})$ forms a $p\times p$ orthonormal matrix with ${\mathbf{Q}}_{i}^{\prime}{\mathbf{B}}_{i}=\mathbf{0}$ and ${\mathbf{B}}_{i}^{\prime}{\mathbf{B}}_{i}={\mathbf{I}}_{p-k_{i}}$. In this case, ${{\cal M}}({\mathbf{B}}_{i})$ represents the orthogonal complement space of ${{\cal M}}({\mathbf{Q}}_{i})$ for $i=1,2$. Note that although ${\mathbf{B}}_{i}$ is not uniquely defined and subject to any orthogonal transformation, $g_{i}(\gamma)$ is invariant under such transformations. If we project the cross moment matrices $\{\boldsymbol{\Sigma}_{y,i}(h,\gamma),\ h=\pm 1,\ldots,\pm h_{0}\}$ onto ${\cal M}({\mathbf{B}}_{i})$, then by (6) we can see that $G(\gamma)$ measures the squared norm of the projections. If $\gamma=\gamma_{0}$ is correctly specified, then the data in the two regimes identified by $\gamma$ do belong to the correct regimes. Since ${\mathbf{M}}_{i}(\gamma_{0})={\mathbf{M}}_{i}$, $i=1,2$, by the definition of ${\mathbf{M}}_{i}$ in (5), we have

[TABLE]

If $\gamma\neq\gamma_{0}$, the data are not correctly separated and at least one of the subsets contains data from both regimes. The following proposition shows that under mild conditions $G(\gamma)>0$ for $\gamma\neq\gamma_{0}$.

Proposition 2**.**

Under Conditions 1–8, $G(\gamma)>0$ if $\gamma\neq\gamma_{0}$.

Since $\gamma_{0}\in(\eta_{1},\eta_{2})$, we can use data corresponding to $t\leq\lfloor\eta_{1}n\rfloor$ and $t>\lfloor\eta_{2}n\rfloor$ to get consistent estimates for ${\cal M}({\mathbf{B}}_{1})$ and ${\cal M}({\mathbf{B}}_{2})$. Specifically, for $i=1,2$, we estimate ${\mathbf{B}}_{i}$ by

[TABLE]

where $\widehat{{\mathbf{q}}}_{i,k}(\eta_{i})$ is the unit eigenvector of $\widehat{{\mathbf{M}}}_{i}(\eta_{i})$ corresponding to the $k$-th largest eigenvalue. Then the sample version of the objective function $G(\gamma)$ is given by

[TABLE]

We propose to estimate the change point location $\gamma_{0}$ by

[TABLE]

whose asymptotic property is provided in Theorem 1.

Theorem 1**.**

Assume Conditions 1–8 in Appendix A.1. If $p^{\delta_{\max}}n^{-1/2}=o(1)$, then for any $\epsilon>0$, with true $k_{1}$ and $k_{2}$ we have

[TABLE]

as $n,p\to\infty$.

By Theorem 1, the proposed estimator $\widehat{\gamma}$ in (11) for the change point location is consistent under mild conditions. It also reveals that the estimation performance can depend critically on the strength of factors in both regimes. In particular, if the factors are strong in both regimes ($\delta_{1}=\delta_{2}=0$), then the estimation is immune to the curse of dimensionality. On the other hand, if factors are weak in one regime, then the resulting estimator can become less efficient as $p$ increases. When factors have different levels of strengths before and after the break, the probability that the $\widehat{\gamma}_{0}$ falls in the weaker regime is larger but the estimation precision in the stronger regime is better. As a result, the overall rate of convergence of $\widehat{\gamma}$ depends on the strength of the weaker regime.

By plugging $\widehat{\gamma}$ into the estimation procedure described in Section 2.2, we obtain the final estimators for ${\cal M}({\mathbf{A}}_{i})$. To be more specific, let $\widehat{r}=\lfloor n\widehat{\gamma}\rfloor$ and

[TABLE]

where $\widehat{{\mathbf{q}}}_{i,k}(\widehat{r})$ is the unit eigenvector of $\widehat{{\mathbf{M}}}_{i}(\widehat{r})$ corresponding to its $k$-th largest eigenvalue. Theorem 2 provides the asymptotic property of the estimated loading spaces when the estimated break date is used.

Theorem 2**.**

Assume Conditions 1–8 in Appendix A.1. If $p^{\delta_{\max}}n^{-1/2}=o(1)$, then as $n,p\to\infty$, with true $k_{1}$ and $k_{2}$, we have

[TABLE]

for $i=1,2$.

By Theorem 2, the convergence rate of the associated loading space estimators is the same as that in Proposition 1 when the true change point location is known. Compared with the results in Liu and Chen (2019) which used a threshold variable to split the data, the ‘helping effect’ disappears and asymptotically there is no interaction between regimes. In practice, the number of factors is typically unknown but can be estimated through a similar eigenvalue ratio estimator as that used in Lam et al. (2011), namely

[TABLE]

where $\widehat{\lambda}_{i,k}(\eta_{i})$ is the $k$-th largest eigenvalue of $\widehat{{\mathbf{M}}}_{i}(\eta_{i})$, $i=1,2$.

3 Determining the Existence of Change Point

3.1 A Self-Normalized Change-Point Test

Although the change-point factor model (2) is able to capture potential structural breaks in the loading space, it can be unnecessarily complicated when there is actually no change point. We shall here consider the problem of determining the existence of change points in the high-dimensional factor model (2). The problem has been studied by Chen et al. (2014) using linear regression of estimated factors, and by Ma and Su (2018) using fused Lasso penalization on block estimators; see also reference therein. However, the aforementioned results both require the noise vector to have a sparse cross-sectional dependence structure; see for example Assumption 4 in Chen et al. (2014) and Assumption A3 in Ma and Su (2018). In addition, the fused Lasso approach of Ma and Su (2018) does not provide a $p$-value which can be an informative measure for quantifying the amount of statistical evidence in favor of the change-point model. To handle noise processes with possibly non-sparse cross-sectional dependence structures, the main challenge is to deal with the unknown scale parameter that associates with the strength of cross-sectional dependence. For this, we propose to adopt the idea of self-normalization (Lobato, 2001; Shao, 2010), which is capable of adaptively handling unknown scale parameters caused by different dependence strengths; see for example Shao (2011), Bai et al. (2016) and Taqqu and Zhang (2019).

Assume that the noise process (\mbox{\boldmath\varepsilon}_{t}) has a constant covariance structure $\boldsymbol{\Sigma}_{t}\equiv\boldsymbol{\Sigma}$, then for any vector ${\mathbf{b}}\in{\cal M}({\mathbf{B}}_{1})$, by (2) and Condition 3 in Appendix A.1 we have

[TABLE]

and

[TABLE]

Therefore, one can determine the existence of change points by testing for structural breaks in the variance of the transformed sequence $({\mathbf{b}}^{\prime}{\mathbf{y}}_{t})$, namely

[TABLE]

versus the alternative that there exists a time point $r_{0}$ such that $\mathrm{var}({\mathbf{b}}^{\prime}{\mathbf{y}}_{r_{0}})\neq\mathrm{var}({\mathbf{b}}^{\prime}{\mathbf{y}}_{r_{0}+1})$. Compared with traditional change-point testing problems, the key difference here is that the scale of the random variable sequence $({\mathbf{b}}^{\prime}{\mathbf{y}}_{t})$ can grow to infinity at an unknown rate depending on the strength of cross-sectional dependence. This is particularly due to the fact that the dimension of the covariance matrix $\boldsymbol{\Sigma}$ is $p\times p$, where $p$ can grow to infinity with the sample size $n$. Under the idiosyncratic error assumption of Chen et al. (2014) and Ma and Su (2018), the covariance matrix $\boldsymbol{\Sigma}$ has a sparse structure by which one can show that $\mathrm{var}({\mathbf{b}}^{\prime}{\mathbf{y}}_{t})=O(\|{\mathbf{b}}\|_{2}^{2})$. However, for situations with non-sparse cross-sectional dependence structures, $\mathrm{var}({\mathbf{b}}^{\prime}{\mathbf{y}}_{t})$ can grow to infinity even when $\|{\mathbf{b}}\|_{2}$ is bounded, and the rate can depend on the strength of the underlying cross-sectional dependence, which is typically unknown in practice. As a result, this can pose a challenge when one seeks an appropriate normalization for the test statistic. For this, we consider the following self-normalization approach.

Let

[TABLE]

we consider the test statistic

[TABLE]

where

[TABLE]

is the self-normalizer that pivotalizes the asymptotic distribution of the test statistic (15) to make it free of any nuisance scale parameter.

Theorem 3**.**

Assume Condition 9 in Appendix A.1. Then under the null hypothesis of no change point, we have as $n,p\to\infty$,

[TABLE]

where

[TABLE]

3.2 A Data-Driven Choice

Although the hypothesis test in (13) is able to detect the change point with any ${\mathbf{b}}\in{\cal M}({\mathbf{B}}_{1})$, we shall here consider a data-driven choice that aims to optimize the power performance. Intuitively, ${\mathbf{b}}$ should be chosen to be far away from the linear space ${\cal M}({\mathbf{B}}_{2})$, which can be achieved by

[TABLE]

As a result, we propose to choose ${\mathbf{b}}$ by

[TABLE]

When ${\mathbf{u}}$ is the unit right singular vector of ${\mathbf{B}}_{2}^{\prime}{\mathbf{B}}_{1}$ that corresponds to the largest singular value, ${\mathbf{b}}$ will be in the column space of ${\mathbf{B}}_{1}$ and the distance of ${\mathbf{b}}$ and ${{\cal M}}({\mathbf{B}}_{2})$ will then be maximized. Let $\widehat{{\mathbf{u}}}(\eta_{1},\eta_{2})$ be the unit right singular vector of $\widehat{{\mathbf{B}}}_{2}^{\prime}(\eta_{2})\widehat{{\mathbf{B}}}_{1}(\eta_{1})$ that corresponds to the largest singular value, the associated test statistic with the data-driven choice can then be obtained by replacing ${\mathbf{b}}{\mathbf{y}}_{t}$ with ${\mathbf{B}}_{1}\widehat{{\mathbf{u}}}(\eta_{1},\eta_{2})$ in (14) and (15). When the degrees of factor strength are different in the two regimes, the data from the one where factors are stronger will provide more information. We can further incorporate this into the data-driven choice of ${\mathbf{b}}$. To be more specific, if $\|{\widehat{{\mathbf{M}}}}_{2}(\eta_{2})\|_{2}>\|{\widehat{{\mathbf{M}}}}_{1}(\eta_{1})\|_{2}$, then factors after the potential change point are likely to be stronger, in which case we can choose ${\mathbf{b}}\in{{\cal M}}({{\mathbf{B}}_{2}})$ that maximizes the distance from ${{\cal M}}({\mathbf{B}}_{1})$ by letting ${\mathbf{b}}={\mathbf{B}}_{2}{\mathbf{v}}$, where ${\mathbf{v}}$ is the right singular vector of ${\mathbf{B}}_{1}^{\prime}{\mathbf{B}}_{2}$ corresponding to the largest singular value.

4 Simulation

In this section, we present a Monte Carlo simulation study to examine the finite sample performance of the proposed inference procedure. For this, we generate entries of the loading matrix ${\mathbf{A}}_{i}$ as independent sample from the uniform distribution on $[-p^{-\delta_{i}/2},p^{-\delta_{i}/2}]$. As a result, the factor strength of ${\mathbf{A}}_{i}$ is characterized by $\delta_{i}$. The associated noise process is generated as a Gaussian process whose covariance matrix has 1 on the diagonal and $\rho_{e}$ for all off-diagonal entries. For simplicity, we use $h_{0}=1$ throughout our numerical examples. For each setting, we generate 1000 realizations, and the results are summarized in the following subsections.

4.1 The Case with No Change Point

We first consider the case with no change point and examine the empirical size performance of the self-normalized change-point test proposed in Section 3. For this, we consider both the setting with strong factors ($\delta=0$) and the setting with weak factors ($\delta=0.25$). Three common factors drive the time series, and the factor process is set to be three independent AR(1) processes with $N(0,4)$ noises process and AR coefficients 0.9, -0.7, and 0.8. $\rho_{e}=0.5$ for the noise process. $\eta_{1}=0.1$ and $\eta_{2}=0.9$. The results are summarized in Table 1, from which we can see that the empirical sizes of the proposed test are reasonably close to their nominal levels $\alpha=10\%,5\%,1\%$ when factors are strong. When $\delta=0.25$, our tests are slightly oversized, and this is because we only utilize information obtained before time $\lfloor\eta_{1}n\rfloor$ and after time $\lfloor\eta_{2}n\rfloor$ to estimate ${\mathbf{b}}$ in (14).

4.2 The Change-Point Case

We now consider the change-point case and examine if the self-normalized change-point test proposed in Section 3 can successfully detect the existence of the change point. For this, we follow the data generating mechanism described in Section 4.1 for generating the factor and noise processes, and set $k_{1}=k_{2}=3$. Let the change point location $\gamma_{0}=0.5$, and we consider four different scenarios on the factor strength, namely SS ($\delta_{1}=\delta_{2}=0$) in which strong factors are used before and after the change point, SW ($\delta_{1}=0$ and $\delta_{2}=0.25$) in which strong factors are used before the change point and weak factors after, WS ($\delta_{1}=0.25$ and $\delta_{2}=0$) in which weak factors are used before the change point and strong factors after, and WW ($\delta_{1}=\delta_{2}=0.25$) in which weak factors are used before and after the change point. The results are summarized in Table LABEL:power_ex2, from which we can see that the proposed test performs reasonably well as it successfully identifies the change point with high probabilities. In addition, the performance seems to improve in general when the sample increases.

After the change point existence has been identified by the test, we shall apply the methods proposed in Section 2 to estimate the change point location and the loading spaces in two regimes. Figure 1 provides the histograms of our change point location estimator $\widehat{\gamma}$ under different settings, from which we can observe the followings. First, when the factor strength is weak in at least one regime, before or after the change point, the estimation efficiency in that weak regime can suffer from the increase in dimension. In contrast, the estimation efficiency in the strong regime does not seem to be affected by the curse of dimensionality. This is in line with the results in Theorem 1; see also the discussions thereafter. Second, it can be seen from the top panels in Figure 1 that, when the factor strengths before and after the change point are different, namely settings SW and WS, the estimation bias, though asymptotically negligible, is more likely to be toward the regime with weaker factors. In particular, when the factor strength after the change point is weaker as in the SW setting, then it is more likely to overestimate $\gamma_{0}$. On the other hand, if the factor strength before the change point is weaker as in the WS setting, then it is more likely to underestimate $\gamma_{0}$. We also provide in Table 3 a summary of the estimation error $|\widehat{\gamma}-\gamma_{0}|$. The estimation error for the loading spaces are summarized in Table 4, from which we can see that the estimation procedure proposed in Section 2.2 performs reasonably well under all the considered settings.

4.3 A Comparison with Existing Results

In this subsection, we will compare the performance of our change point detection procedure with Chen et al. (2014) and Han and Inoue (2015) who both focus on single change point detection in factor models but do not allow strong cross-sectional dependence in noise process. Set $\eta_{1}=0.2$ and $\eta_{2}=0.8$. Data are simulated with processes described in Section 4.1 and in Section 4.2 with $\rho_{e}=0.95$ to compare the sizes and powers of different tests, respectively. From Table 6 and Table 6, we can see that our test controls the sizes better than Wald test proposed by Chen et al. (2014) and is more powerful than LM and Wald tests by Han and Inoue (2015) when the strong cross-sectional exists in noise process.

5 Real data analysis

We applied our method to the Stock-Watson data (Stock and Watson, 1998, 2005), containing 132 U.S. monthly economic indicators from March 1960 to December 2003, with $n=526$ and $p=132$. The data include real output and income, employment, real retail, manufacturing and trade sales, consumption, interest rates, price index and other economic indicators. Stock and Watson (2005) provided more detailed information about this data set and transformations needed for stationarity before analysis.

Set $h_{0}=1$, $\eta_{1}=0.1$, $\eta_{2}=0.9$, and level of significance $\alpha=5\%$. We applied the wild binary segmentation method (Fryzlewicz, 2014) to check if there are multiple change points. The result shows that during this period there is only one change point.

Left and right panels in Figure 5 demonstrate the ratio of eigenvalues of $\widehat{{\mathbf{M}}}(\eta_{1})$ and $\widehat{{\mathbf{M}}}(\eta_{2})$, respectively, and it reaches its minimum values at 1 and 2, which implies that $\widehat{k}_{1}=1$ and $\widehat{k}_{2}=2$. Using methods described in Section 2.3, we have $\widehat{\gamma}=0.635$.

It means that the dynamics of economic indicators experienced once a permanent structural change around November 1988 possibly due to saving and loan crisis (Curry and Shibut, 2000). Before the change, there was one common factor, while after the change, there was two common factors driving 132 economic indicators.

Extending the method described in Lam et al. (2011), we can estimate the residuals of factor models with

[TABLE]

We compare our method with Wald test by Chen et al. (2014), LM and Wald tests by Han and Inoue (2015) through white noise test (Chang et al., 2017) and residual sum of squares (RSS). All the three tests reject the null and agree that there is a structural break in data. With the estimated break dates, we fix the number of factors before and after break at 1 and 2 for comparison and calculate the residual series. Table 7 reports the $p$-values associated with white noise test for the residual series and RSS. Our model is validated by the white noise test with $p$-value 0.07. It confirms that after extracting the common factors, there is no serial dependence left. The residuals series from Wald test by Chen et al. (2014) and Wald test by Han and Inoue (2015) fail the white noise test. Although $p$-value for LM test by Han and Inoue (2015) is greater than 0.05, Han and Inoue (2015) yields a larger RSS. In sum, our method can successfully capture all the serial dynamics of the data, and lower the residual sum of squares comparing to others.

6 Conclusion

Factor models have been frequently used in the study of high-dimensional time series. Their associated change-point analyses were mostly studied under the framework of Bai and Ng (2002) and Forni et al. (2004). However, such type of factor models suffer from the identifiability issue in the finite dimensional case and cannot be used to explain the phenomenon of strong cross-sectional dependence in noise process. To alleviate the aforementioned problems, Lam et al. (2011) proposed an alternative framework for factor analysis of high-dimensional time series. In the current paper, we develop statistical methods and their associated theory for change-point analysis of the recently proposed factor model of Lam et al. (2011). In the presence of a change point, we propose consistent statistical estimators for the change point location and factor loading spaces before and after the change point, and provide their explicit convergence rates. In particular, our results reveal that the resulting estimators can have very different asymptotic behaviors in response to the phenomenon of high dimensionality depending on the factor strength. To be more specific, for processes with strong factor strength, the convergence rate will not be affected by the dimension, while the curse of dimensionality can be observed for processes with weak factor strength; see the discussion in Section 2 and our numerical results in Section 4. We also propose a self-normalized test for determining the change-point existence, which, due to the self-normalization, can adaptively handle the case with strong and weak cross-sectional dependence. It can be seen from our numerical results in Section 4 that the proposed change-point test performs reasonably well in terms of both the size and power.

Appendix A.1. Regularity Conditions

Define

[TABLE]

[TABLE]

and three intervals

[TABLE]

For any $0\leq c_{1}<c_{2}\leq 1$ and $c_{1}$ and $c_{2}$ are from the same one of the three intervals, $I_{1}$, $I_{2}$ or $I_{3}$, let

[TABLE]

The following regularity conditions are needed for theoretical properties.

Condition 1. Let ${\cal F}_{\ell}^{j}$ be the $\sigma$-field generated by $\{({\mathbf{x}}_{t,1},{\mathbf{x}}_{t,2}):\ell\leq t\leq j\}$ The latent process $\{{\mathbf{x}}_{t,1},{\mathbf{x}}_{t,2}\}$ is $\alpha$-mixing with mixing coefficients satisfying

[TABLE]

for $i=1,2$ and some $\zeta>2$, where $\alpha(t)=\sup_{j}\sup_{A\in{\mathcal{F}}_{-\infty}^{j},B\in{\mathcal{F}}_{j+t}^{\infty}}|P(A\cap B)-P(A)P(B)|$.

Condition 2. For any $i=1,2$, $j=1,\ldots,k_{i}$, and $t=1,\ldots,n$, $E(|x_{t,i,j}|^{4\zeta})<\sigma_{x}^{4\zeta}$, where $x_{t,i,j}$ is the $j$-th element of ${\mathbf{x}}_{t,i}$, $\sigma_{x}>0$ is a constant, and $\zeta$ is given in Condition 1.

Condition 3. \{\mbox{\boldmath\varepsilon}_{t}\} is a white noise process. \mbox{\boldmath\varepsilon}_{t} and $\{{\mathbf{x}}_{t,1},{\mathbf{x}}_{t,2}\}$ are uncorrelated given ${\cal F}_{-\infty}^{t-1}$. Each element of $\boldsymbol{\Sigma}$ remains bounded by a positive constant $\sigma_{\varepsilon}^{2}$ as $p$ increases to infinity.

Condition 4. For $i=1,2$, there exists a constant $\delta_{i}\in[0,1]$ such that $\|{\mathbf{A}}_{i}\|_{2}^{2}\asymp\|{\mathbf{A}}_{i}\|^{2}_{\min}\asymp p^{1-\delta_{i}}$, as $p$ goes to infinity.

Condition 5. $\gamma_{0}\in(\eta_{1},\eta_{2})$. For any $\gamma\in[\eta_{1},\eta_{2}]$, there exists an integer $h_{i}\in[1,h_{0}]$ such that $\boldsymbol{\Sigma}_{x,i}(h_{i},\gamma)$ and $\boldsymbol{\Sigma}_{x,i}(h_{i},\gamma)$ are full rank, and $\|\boldsymbol{\Sigma}_{x,i}(h_{i},\gamma)\|_{\min}$ is uniformly bounded above 0, for $i=1,2$.

Condition 6. ${\mathbf{M}}_{i}(\gamma)$ admits $k_{i}$ distinct positive eigenvalues, for $\gamma\in[\eta_{1},\eta_{2}]$, $i=1,2$.

Condition 7. There exists a positive constant $d$ such that ${\cal D}({\cal M}({\mathbf{Q}}_{1}),{\cal M}({\mathbf{Q}}_{2}))>d$ as $n$ and $p$ go to infinity.

Condition 8. For any $\gamma\in(\eta_{1},\gamma_{0})$, there exists an integer $h_{1}^{*}\in[1,h_{0}]$ such that $\boldsymbol{\Gamma}_{x}(h_{1}^{*},\gamma,\gamma_{0}-h_{1}^{*}/n)$ is full rank. For any $\gamma\in(\gamma_{0},\eta_{2})$, there exists an integer $h_{2}^{*}\in[1,h_{0}]$ such that $\boldsymbol{\Gamma}_{x}(h_{2}^{*},\gamma_{0},\gamma)$ is full rank. The minimum singular values of these two matrices mentioned are uniformly bounded above $u_{0}>0$.

Let

[TABLE]

and we use $\rightsquigarrow$ to denote the weak convergence of processes (Van der Vaart and Wellner, 1996).

Condition 9. There exists a deterministic sequence $(\sigma_{n,p})$ such that

[TABLE]

where $\mathcal{D}=\{(s,t),\ 0\leq s\leq t\leq 1\}$ and $\{\mathbb{B}(t)\}_{t\in[0,1]}$ is a standard Brownian motion.

Appendix A.2 Proofs and Lemmas

In this section, we mainly focus on the mathematical proofs for before the break and when $\epsilon>0$. The results for after the break or $\epsilon<0$ are included, but most of proofs are omitted since they are quite similar. For any fixed $\epsilon\neq 0$, there exists a positive integer $N$ such that when $n\geq N$, $|\epsilon|>(h+1)/n$, therefore, for Lemmas 2-7, we only consider when $|\epsilon|>(|h|+1)/n$. In addition, the model is not distinguishable for all values between $[k/n,(k+1)/n]$ as the break point, so for simplicity we treat $\epsilon n$ as an integer in the proofs. We use $Cs$ to denote generic uniformly positive constants which only depend on the parameters.

Lemma 1**.**

For $0\leq c_{1}<c_{2}\leq 1$, and $c_{1}$ and $c_{2}$ are from the same one of the three intervals, $I_{1}$, $I_{2}$ or $I_{3}$, let

[TABLE]

Under Conditions 1 and 2, for any $h\in[1,h_{0}]$, it holds that

[TABLE]

[TABLE]

where $\alpha=\sum_{t=1}^{\infty}\alpha(t)^{1-2/\zeta}$, and $k_{\max}=\max\{k_{1},k_{2}\}$.

Proof: Let $a_{q,\ell}$ and $\widehat{a}_{q,\ell}$ be the $(q,\ell)$-th entry in $\boldsymbol{\Gamma}_{x}(h,c_{1},c_{2})$ and $\widehat{\boldsymbol{\Gamma}}_{x}(h,c_{1},c_{2})$ respectively. By Condition 2 and Jensen’s inequality we know that ${\rm E}(x_{t,i,j}^{2})<\sigma_{x}^{2}$ and ${\rm E}(x_{t,i,j}^{4})<\sigma_{x}^{4}$, for $i=1,2$, $j=1,\ldots,k_{i}$, and $t=1,\ldots,n$. Let $x_{t,q}$ be the $q$-th entry in ${\mathbf{x}}_{t}$. We have ${\rm E}(x_{t,j}^{2})<\sigma_{x}^{2}$ and ${\rm E}(x_{t,j}^{4})<\sigma_{x}^{4}$, for $j=1,\ldots,k_{1}$ when $t\leq\gamma_{0}n$, for $j=1,\ldots,k_{2}$ when $t>\gamma_{0}n$. By Cauchy-Schwarts inequality,

[TABLE]

It follows $\|\boldsymbol{\Gamma}_{x}(h,c_{1},c_{2})\|_{2}^{2}\leq\|\boldsymbol{\Gamma}_{x}(h,c_{1},c_{2})\|^{2}_{F}\leq k_{\max}^{2}\sigma_{x}^{4}$.

By Proposition 2.5 in Fan and Yao (2003),

[TABLE]

Thus, ${\rm E}\|\widehat{\boldsymbol{\Gamma}}_{x}(h,c_{1},c_{2})-{\boldsymbol{\Gamma}}_{x}(h,c_{1},c_{2})\|_{2}^{2}\leq{\rm E}\|\widehat{\boldsymbol{\Gamma}}_{x}(h,c_{1},c_{2})-{\boldsymbol{\Gamma}}_{x}(h,c_{1},c_{2})\|_{F}^{2}\leq(3h+8\alpha)k_{\max}^{2}\sigma_{x}^{4}/N(c_{1},c_{2})$.

Lemma 2**.**

Under Conditions 1-4 and 6, for $\epsilon\in(-\gamma_{0},1-\gamma_{0})$ and $|\epsilon|>(h+1)/n$, it holds that

[TABLE]

where $\sigma_{0}=\max\{\sigma_{x},\sigma_{\varepsilon},1\}$, and $a_{1}>1$ satisfies $\|{\mathbf{A}}_{i}\|_{2}\leq a_{1}p^{1/2-\delta_{i}/2}$, for i=1,2.

Proof: When $\epsilon>0$,

[TABLE]

Condition 4 implies that there exists a positive constant $a_{0}$ such that $\|{\mathbf{A}}_{i}\|_{2}\leq a_{0}p^{1/2-\delta_{i}/2}$ for $i=1,2$. Let $a_{1}=\max\{a_{0},1\}$. By Lemma 1 and Condition 3, we have

[TABLE]

Since ${\mathbf{x}}_{t}$ and \mbox{\boldmath\varepsilon}_{t} are independent,

[TABLE]

and

[TABLE]

where $\varepsilon_{t,j}$ is the $j$-th entry in \mbox{\boldmath\varepsilon}_{t}. On the other hand,

[TABLE]

Together with (16) and (17) we have

[TABLE]

where $\sigma_{0}=\max\{\sigma_{x},\sigma_{\varepsilon},1\}$. Similarly, we can show that

[TABLE]

Hence,

[TABLE]

When $\epsilon<0$, it can be proven in a similar fashion.

Lemma 3**.**

Under Conditions 1-4 and 6, for $\epsilon\in(-\gamma_{0},1-\gamma_{0})$ and $|\epsilon|>(h+1)/n$, it holds that

[TABLE]

and

[TABLE]

Proof: By the definition of $\boldsymbol{\Sigma}_{y,i}(h,\gamma)$ and Lemma 1, when $\epsilon>(h+1)/n$, we have

[TABLE]

and

[TABLE]

Lemma 4**.**

Under Conditions 1-4 and 6, for $\epsilon\in(-\gamma_{0},1-\gamma_{0})$and $|\epsilon|>(h+1)/n$, it holds that

[TABLE]

and

[TABLE]

Proof: When $\epsilon>(h+1)/n$, by Lemma 1

[TABLE]

and

[TABLE]

Lemma 5**.**

Under Conditions 4 and 7, it holds that

[TABLE]

where $a_{2}$ is a positive constant such that $a_{2}p^{1/2-\delta_{i}/2}\leq\|{\mathbf{A}}_{i}\|_{\min}$ for $i=1,2$, and $\tau=\min\{k_{2}/k_{1},k_{1}/k_{2}\}$.

Proof: Note that

[TABLE]

On the other hand,

[TABLE]

Hence, \mbox{tr}({\mathbf{Q}}_{2}^{\prime}{\mathbf{B}}_{1}{\mathbf{B}}_{1}^{\prime}{\mathbf{Q}}_{2})=k_{\min}[\mbox{\mathcal{D}}({\cal M}({\mathbf{Q}}_{1}),{\cal M}({\mathbf{Q}}_{2}))]^{2}+k_{2}-k_{\min}. Then we have $\|{\mathbf{B}}_{1}^{\prime}{\mathbf{Q}}_{2}\|_{2}^{2}\geq\mbox{tr}({\mathbf{Q}}_{2}^{\prime}{\mathbf{B}}_{1}{\mathbf{B}}_{1}^{\prime}{\mathbf{Q}}_{2})/k_{2}\geq\tau d^{2}$. Condition 4 implies that there exists a positive constant $a_{2}$ such that $\|{\mathbf{A}}_{2}\|_{\min}\geq a_{2}p^{1/2-\delta_{i}/2}$. It follows

[TABLE]

Lemma 6**.**

Under Conditions 1-8, we have $G(\gamma_{0})=0$ and for $\epsilon\in(\eta_{1}-\gamma_{0},\eta_{2}-\gamma_{0})$ and $|\epsilon|>(h+1)/n$,

[TABLE]

Proof: Under Condition 8, by Lemmas 1 and 5, and Theorem 6 in Merikoski and Kumar (2004), when $\epsilon<-(h+1)/n$,

[TABLE]

when $\epsilon>(h+1)/n$,

[TABLE]

From Lemma 4, we have

[TABLE]

Lemma 7**.**

Under Conditions 1-6, if $p^{\delta_{\max}}n^{-1/2}=o(1)$, with true $k_{1}$ and $k_{2}$, as $n,p,\to\infty$, we have

[TABLE]

where $C$ is a generic uniformly positive constant which only depends on the parameters.

Proof: Let $Y_{t}=x_{t,i,q}x_{t+h,i,\ell}-{\rm E}(x_{t,i,q}x_{t+h,i,\ell})$. Condition 2 indicates that there exists a positive constant $\sigma_{y}$ such that ${\rm E}(|Y_{t}^{2\zeta}|)<\sigma_{y}^{2\zeta}$. For any $0\leq c_{1}<c_{2}\leq 1$, by Cauchy-Schwartz inequality, Proposition 2.5 in Fan and Yao (2003) and Lemma 1,

[TABLE]

It follows

[TABLE]

Thus we have

[TABLE]

where $C_{1}$,$C_{2}$ and $C_{3}$ are positive constants and depend only on the parameters.

Hence, with Lemmas 2 and 3

[TABLE]

Following the proof of Theorem 1 in Lam et al. (2011), we can reach the conclusion.

Lemma 8**.**

*Under Conditions 1-8, for $\epsilon\in[-\gamma_{0},1-\gamma_{0}]$, it holds that *

[TABLE]

Proof: By the definition of ${\mathbf{M}}_{i}(\eta_{i})$, we can see that ${\cal M}({\mathbf{B}}_{i})={\cal M}({\mathbf{B}}_{i}(\eta_{i}))$. It implies that there exists an orthogonal $(p-k_{i})\times(p-k_{i})$ matrix ${\mathbf{R}}_{i}$ such that ${\mathbf{B}}_{i}={\mathbf{B}}_{i}(\eta_{i}){\mathbf{R}}_{i}$.

[TABLE]

By the definition of $\widehat{G}(\gamma)$ we have,

[TABLE]

By Lemmas 2-4 and 7,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

From (20), it follows,

[TABLE]

Proof of Proposition 1. We can apply Theorem 1 in Lam et al. (2011) to obtain the conclusions.

Proof of Theorem 1. Since $G(r)\geq 0$ and $G(r_{0})=0,$ for any fixed $\epsilon>(h+1)/n$, it follows that

[TABLE]

By Lemma 6, Lemma 8, and Chebyshev’s inequality, if $p^{\max}n^{-1/2}=o(1)$ and $n$ is large enough, when $\widehat{r}>r_{0}+\epsilon$, we have

[TABLE]

Hence, there exists a constant $C$ such that

[TABLE]

Proof of Theorem 2. When $\widehat{\gamma}>\gamma_{0}$, from Proposition 1, it follows

[TABLE]

Now we start to investigate the asymptotic properties of ${\cal M}[\widehat{{\mathbf{Q}}}_{1}(\widehat{\gamma})]$ when $\widehat{\gamma}>\gamma_{0}$.

Lemmas 2-4 imply that

[TABLE]

Under Conditions 2 and 4, it follows from Lemma 1

[TABLE]

Hence,

[TABLE]

If $p^{\delta_{\max}}n^{-1/2}=o(1)$, together with (25), we have

[TABLE]

Theorem 1 tells us if $\widehat{r}>r_{0}$, $|\widehat{r}-r_{0}|=O_{p}(p^{\delta_{2}}n^{-1/2})$. Therefore,

[TABLE]

Under Condition 5, by Theorem 9 in Merikoski and Kumar (2004), we can see that

[TABLE]

Following the proof of Theorem 2 in Liu and Chen (2016), we have

[TABLE]

as $n,p\to\infty$, when $\widehat{r}>r_{0}$.

The conclusions for $\widehat{r}<r_{0}$ can be proven in a similar way.

Proof of Theorem 3. The key idea is to show that the test statistic (15) can be approximated by a functional of $\{(\lfloor nt\rfloor-\lfloor ns\rfloor\vee 1+1)(\widehat{\nu}_{\lfloor ns\rfloor\vee 1,\lfloor nt\rfloor}^{\star}-{\mathbf{b}}^{\prime}\boldsymbol{\Sigma}{\mathbf{b}})\}_{(s,t)}$ with a negligible difference under the null hypothesis of no change point. For this, by Condition 10 in Appendix A.1, we have

[TABLE]

and

[TABLE]

Therefore, the left hand sides of both (26) and (27) are of order $O_{p}(1)$. Let

[TABLE]

then we have

[TABLE]

and thus by (26) and (27),

[TABLE]

By condition 10 in Appendix A.1,

[TABLE]

we have

[TABLE]

Note that one can write

[TABLE]

the result follows by the continuous mapping theorem.