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.
72,164 characters · 13 sections · 41 citation commands
Structural Change Detection in High-Dimensional Transformed Factor Models via Canonical Correlation Analysis
{\sl Keywords}: Structural change; High-dimensional time series; Transformed factor model; Canonical correlation analysis; Factor number estimation
High-dimensional time series are now routinely encountered in economics, finance, environmental studies, and many other empirical fields, where a large number of variables are observed over time and often exhibit both temporal dynamics and strong cross-sectional dependence. Although classical statistical models provide a foundation for time-series analysis, many of them were developed primarily for low-dimensional settings. For instance, the vector autoregressive moving-average (VARMA) model is a canonical parametric framework for multivariate time-series analysis. In high-dimensional settings, however, unrestricted VARMA models quickly become overparameterized and may suffer from identification difficulties. One way to mitigate these difficulties is to impose simplified dynamic structures, such as the scalar component model (SCM) proposed by TiaoTsay_1989, or sparsity and regularization restrictions; see, for example, ShojaieMichailidis_2010, HanLuLiu_2015, and Davis2016.
An alternative and widely used strategy is dimension reduction, which reduces both the statistical and computational complexity of high-dimensional time-series modeling. Prominent dimension-reduction approaches include the canonical correlation analysis (CCA) framework of BoxTiao_1977 and principal component analysis (PCA)-based methods developed in, for example, StockWatson_2002. Factor models have become especially influential because they represent temporal and cross-sectional dependence through a small number of unobserved common factors; see, among others, ChamberainRothschild_1983, BaiNg_Econometrica_2002, bai2003, StockWatson_2005, panyao2008, LamYaoBathia_Biometrika_2011, lamyao2012, and changguoyao2015. Within this literature, LamYaoBathia_Biometrika_2011 and lamyao2012 estimate factor loadings and latent factors through the eigen-decomposition of sample cross-autocovariance matrices. By contrast, gaotsay_2019, building on TiaoTsay_1989, proposed a transformed factor approach that uses CCA of time-lagged moment matrices to separate dynamically dependent common factors from serially uncorrelated idiosyncratic components. A key feature of this framework is that, after a nonsingular linear transformation of the original high-dimensional series, the dynamically dependent factor component can be distinguished from the noise component through nonzero canonical correlations between the current observations and their lagged values.
Most classical high-dimensional factor models assume that the loading matrix is time invariant. In empirical applications, however, the underlying data-generating mechanism may undergo structural changes that alter model parameters. For example, major policy adjustments, financial crises, technological changes, or shifts in market conditions may induce abrupt changes in factor loadings, thereby changing the effects of latent factors on observed variables. Ignoring such unobserved structural changes may lead to substantial estimation bias, spurious inflation in the estimated number of factors, and degraded forecasting or inferential performance. Consequently, detecting structural changes is important for reliable factor modeling and subsequent empirical analysis.
A growing literature has investigated structural changes in high-dimensional factor models; see, among others, BreEick_2011, Chenetal_2014, LiuChen_D_2016, and MaSu_2018. LiuChen_D_2016 proposed a regime-switching dynamic factor model in which regime transitions are governed by an unobserved Markov chain, so that the transition times effectively serve as change points. Subsequently, BaiHanShi_2020 developed a least-squares-based estimation framework for change-point detection in high-dimensional time series. Their method assumes that a piecewise factor model with a structural change point can be recast as an equivalent pseudo-factor model without an explicit change point. Their procedure first estimates the pseudo-factor sequence and then selects the time point that minimizes the residual sum of squares. LiuChen_D_2020 further studied a threshold-variable approach to identifying change points within the factor estimation procedure. Building on this line of research, LiuZhang_2022 proposed a data-driven projection method for sequentially detecting change points in latent factor models. Their method projects a second-order sample cross-moment matrix, constructed under a candidate partition, onto the estimated orthogonal noise space associated with the true change-point partition, and then minimizes the squared norm of the projected matrix.
Estimating the number of latent factors is another fundamental model-selection problem in high-dimensional time-series analysis. Onatski_2010 and Onatski_2012 proposed thresholding methods for determining the number of factors based on the empirical spectral distribution of the sample covariance or autocovariance matrix. lamyao2012 estimated the number of factors by maximizing the ratio of adjacent eigenvalues derived from the cross-autocovariance matrix. Although eigenvalue-ratio methods are computationally convenient, they may be sensitive to dominant factors and can suffer from a masking effect, leading to substantial underestimation of the true number of factors. To address this issue, Wu_2016 proposed applying a monotone nonlinear transformation, such as the standard normal cumulative distribution function, to the eigenvalues before forming adjacent ratios. Such a transformation reduces the scale disparity between strong and weak factors and mitigates the masking effect induced by dominant factors. Following this idea, Xia_TCR_2017 introduced the transformed contribution ratio (TCR) method, which incorporates cumulative contribution information to improve finite-sample accuracy. In the transformed factor setting, gaotsay_2019 estimated the number of factors through a sequential hypothesis test for zero canonical correlations and established the asymptotic chi-square distribution of the proposed test statistic under the null hypothesis.
In the presence of structural change, however, factor-number estimation and change-point estimation are intrinsically connected. A misspecified break location may mix observations from different regimes and thereby distort the estimated factor dimension. Conversely, an inaccurate factor-number estimate may affect the construction of the noise subspace and distort the change-point criterion. This circular dependence is particularly relevant in transformed factor models, where the dynamic signal is identified through canonical correlations and the noise space is determined by the number of zero canonical correlations. Therefore, a change-point procedure for transformed factor models should jointly account for the unknown break location and the regime-specific factor numbers.
Motivated by the transformed factor framework of gaotsay_2019 and the structural change literature for high-dimensional factor models, this paper develops a CCA-based procedure for estimating a single structural change point in high-dimensional transformed factor models. The proposed method exploits the low-rank canonical-correlation structure induced by dynamically dependent common factors. At the true change point, the regime-specific canonical correlation matrix has a signal subspace associated with the common factors and a complementary noise subspace associated with zero canonical correlations. We therefore construct an eigenvalue-based criterion that measures the residual canonical dependence in the estimated noise subspace. The population counterpart of this criterion is minimized at the true change point, which provides the basis for identification and estimation.
The contribution of this paper is threefold. First, we introduce a canonical-correlation noise-space criterion for structural change detection in transformed factor models. Unlike methods based directly on covariance or cross-moment matrices, the proposed criterion targets the dynamic dependence structure captured by canonical correlations and is therefore tailored to the transformed factor framework. Second, we develop an alternating iterative estimation (AIE) procedure to handle the mutual dependence between the unknown change point and the regime-specific factor numbers. The procedure alternates between factor-number estimation and change-point estimation until convergence, thereby reducing the sensitivity of the change-point estimator to initial factor-number choices. Third, we establish the asymptotic properties of the proposed estimators under suitable mixing and moment conditions. The resulting convergence rates explicitly reflect the effects of factor strength, cross-sectional dimension, and sample size. Monte Carlo experiments and empirical applications further illustrate the finite-sample performance and practical usefulness of the proposed method.
The remainder of this paper is organized as follows. Section (ref) introduces the single-change-point transformed factor model and develops the proposed CCA-based change-point estimation procedure. It also presents the alternating iterative estimation algorithm for jointly estimating the change point and the regime-specific factor numbers. Section (ref) establishes the theoretical properties of the proposed estimators. Section (ref) reports Monte Carlo simulation results. Section (ref) applies the proposed method to intraday stock returns and daily temperature series. Section (ref) concludes. Technical proofs are collected in the Appendix.
We conclude this section by introducing notation used throughout the paper. For a matrix ${\mathbf H}=(h_{ij})$, $||{\mathbf H}||_2=\sqrt{\lambda_{\max} ({\mathbf H}^{T} {\mathbf H} ) }$ and $||{\mathbf H}||_F=\sqrt{\mbox{tr}({\mathbf H}^{T} {\mathbf H})}$ denote, respectively, the $L_2$ norm and the Frobenius norm, where $\lambda_{\max} ({\mathbf H}) $ is the largest eigenvalue of ${\mathbf H}$, and $\mbox{tr}(\cdot)$ denotes the trace operator. $\text{Diag}[{\mathbf H}_{1},{\mathbf H}_{2},\ldots]$ denotes a block diagonal matrix, and $\mathcal{M}({\mathbf H})$ denotes the column space spanned by the matrix ${\mathbf H}$. $\sigma_{\min}^{+}$ denotes the minimum positive singular value. The superscript $T$ denotes the transpose of a vector or matrix. In addition, we use the notation $a \asymp b$ to denote $a=O(b)$ and $b=O(a)$. Finally, $\left \lfloor x \right \rfloor $ and $\left \lceil x \right \rceil $ denote the greatest integer less than or equal to $x$ and the smallest integer greater than or equal to $x$, respectively.
Let $\boldsymbol{Y}=({\mathbf y}_{1},{\mathbf y}_{2},\ldots,{\mathbf y}_{T})$ be a $p$-dimensional time series, where each observation at time $t$ is a column vector ${\mathbf y}_{t}=(y_{1t},y_{2t},\ldots,y_{pt})^{T}$. We allow the latent factor structure to undergo a single structural change point, which partitions the sample into two temporal regimes. For any $1\le t \le T$, we consider a structural transformed factor model with a single change point:
Here, ${\mathbf f}_{t}^{(i)}\in\mathbb R^{r_i}$, $i=1,2$, denotes the latent common factor vector in the pre- and post-change regimes, respectively, and the positive integer $k_0$ is the true but unknown location of the change point. Furthermore, $\mbox{\boldmath$\varepsilon$}_{t}^{(i)} \in \mathbb R^{v_i}$ denotes a serially uncorrelated idiosyncratic component with mean zero and covariance matrix $\textnormal{Cov}(\mbox{\boldmath$\varepsilon$}_{t}^{(i)})$. For $i=1,2$, the matrices $\widetilde{\mathbf L}_{1}^{(i)}\in \mathbb R^{p\times r_i}$ denote the unknown regime-specific loading matrices associated with the common factors, whereas $\widetilde{\mathbf L}_{2}^{(i)}\in \mathbb R^{p\times v_i}$ are the corresponding loading matrices associated with the idiosyncratic components.
When $r_i+v_i=p$ and the matrix $[\widetilde{\mathbf L}_{1}^{(i)},\widetilde{\mathbf L}_{2}^{(i)}]$ is nonsingular, Model (ref) can be regarded as a transformed factor model. In particular, by defining ${\mathbf T}_i=[\widetilde{\mathbf L}_{1}^{(i)},\widetilde{\mathbf L}_{2}^{(i)}]^{-1}$, we have ${\mathbf T}_i{\mathbf y}_t=( {\mathbf f}_t^{(i)}{^T}, \mbox{\boldmath$\varepsilon$}_t^{(i)}{^T})^T$ within regime $i$. Hence, after the regime-specific linear transformation, the observed vector is decomposed into a signal component and an innovation component. This formulation is consistent with the scalar component model of TiaoTsay_1989 and gaotsay_2019. For each $i=1,2$, we assume ${\mathbf f}_{t}^{(i)}$ and $\mbox{\boldmath$\varepsilon$}_{t}^{(i)}$ are mutually independent. To model the temporal dependence, we further assume that ${\mathbf f}_{t}^{(i)}$ follows a stationary VAR($d$) process,
where $\{{\mathbf u}_{t}^{(i)}\in\mathbb{R}^{r_i},t=1,2,\ldots\}$ is the innovation series with diagonal covariance matrix, and $\boldsymbol{\Phi}_{j}^{(i)}$ denotes the autoregressive coefficient matrix at lag $j$.
Model ((ref)) can be equivalently reformulated into the following compact form:
where $\widetilde{\mathbf L}^{(i)}=(\widetilde{\mathbf L}_1^{(i)},\widetilde{\mathbf L}_2^{(i)})$ is a nonsingular $p \times p$ real-valued matrix, and $\boldsymbol{\eta}_{t}^{(i)}=[{\mathbf f}_{t}^{(i)^T},\mbox{\boldmath$\varepsilon$}_{t}^{(i)^T}]^{T}$ is a $p$-dimensional vector whose covariance matrix is block diagonal, namely \[\mathrm{Diag}[\textnormal{Cov}({\mathbf f}_t^{(i)}),\textnormal{Cov}(\mbox{\boldmath$\varepsilon$}_t^{(i)})].\] Note that any nonsingular linear transformation of the original series ${\mathbf y}_{t}$ does not alter the canonical correlation between ${\mathbf y}_{t}$ and its lagged values; consequently, the representation in $(\ref{cpf-transform-simple})$ is not unique. In light of this rotational non-uniqueness, we consider the following standardized transformed factor model with a change point:
where $\boldsymbol{\Sigma}_{\mathbf y}=\textnormal{Cov}({\mathbf y}_t)$ is the covariance matrix used for standardization, and ${\mathbf L}^{(i)}=({\mathbf L}_{1}^{(i)},{\mathbf L}_{2}^{(i)})=\boldsymbol{\Sigma}_{{\mathbf y}}^{-1/2} \widetilde{\mathbf L}^{(i)} $ is a $p\times p$ orthogonal matrix satisfying ${\mathbf L}^{(i)}{\mathbf L}^{(i)^{T}}={\mathbf L}^{(i)^{T}}{\mathbf L}^{(i)}={\mathbf I}_{p}$. Let $\gamma_1$ and $\gamma_2$ with $0<\gamma_1<\gamma_2<1$ denote the lower and upper truncation constants for the change-point search interval, respectively, so that the effective search-window length is $n=\left \lfloor \gamma_{2}T \right \rfloor-\left \lfloor \gamma_{1}T \right \rfloor$. For any candidate change point $k$ within the interval $ [\left \lfloor \gamma_{1}T \right \rfloor,\left \lfloor \gamma_{2}T \right \rfloor]$, let $\mathbb{I}_{t,1}(k)$ and $\mathbb{I}_{t,2}(k)$ denote the corresponding indicator functions such that $\mathbb{I}_{t,1}(k)=1$ if $1\le t \le k$ and $\mathbb{I}_{t,2}(k)=1$ if $k+1\le t \le T$. For notational convenience, we use \(\mathbb{I}_{t,i}(k)\) to indicate the segment induced by a candidate split \(k\). All covariance and cross-covariance matrices indexed by \(i\) and \(k\) are understood as segment-restricted quantities computed within the corresponding segment. Thus, for example, \[ \textnormal{Cov}({\mathbf y}_t\mathbb{I}_{t,i}(k)) \equiv \textnormal{Cov}({\mathbf y}_t\mid \mathbb{I}_{t,i}(k)=1), \] and the same convention applies to cross-covariance matrices such as \(\textnormal{Cov}({\mathbf y}_t\mathbb{I}_{t,i}(k),{\mathbf y}_{t,m}\mathbb{I}_{t,i}(k))\) below. Under the partition induced by the candidate change point $k$, the latent factors can be explicitly recovered as follows:
where \(\boldsymbol{\Sigma}_{{\mathbf y},i}(k)=\textnormal{Cov}({\mathbf y}_t\mathbb{I}_{t,i}(k))\) denotes the segment-restricted covariance matrix of \({\mathbf y}_t\) in the \(i\)-th segment induced by the candidate split \(k\), following the notational convention above.
Let ${\mathbf y}_{t,m}=({\mathbf y}_{t-1}^T ,{\mathbf y}_{t-2}^T,\ldots,{\mathbf y}_{t-m}^T)^T$ be the $mp$-dimensional vector of past $m$ lagged observations of ${\mathbf y}_t$ where the integer m satisfies $1\le m \le t-1$. We then define the cross-covariance matrix $\boldsymbol{\Sigma}_{{\mathbf y}{\mathbf y}_m,i}(k)=\textnormal{Cov}({\mathbf y}_{t}\mathbb{I}_{t,i}(k),{\mathbf y}_{t,m}\mathbb{I}_{t,i}(k))$ and the autocovariance matrix of the lagged vector $\boldsymbol{\Sigma}_{{\mathbf y}_m,i}(k)=\textnormal{Cov}({\mathbf y}_{t,m}\mathbb{I}_{t,i}(k))$. To estimate the parameters $\{{\mathbf L}^{(1)},{\mathbf L}^{(2)},{\mathbf f}_{t}^{(1)},{\mathbf f}_{t}^{(2)},k_0\}$ under model ((ref)), we compare the canonical correlation structures between ${\mathbf y}_t$ and its lagged vector ${\mathbf y}_{t,m}$ across the two regimes induced by the candidate change point $k$. From a statistical perspective, applying CCA to maximize the linear association between ${\mathbf y}_t$ and ${\mathbf y}_{t,m}$ is equivalent to eigen-decomposition on the following target matrix:
Let $\boldsymbol{\eta}_{t,m}^{(i)}=(\boldsymbol{\eta}_{t-1}^{(i)\ ^T} ,\boldsymbol{\eta}_{t-2}^{(i)\ ^T},...,\boldsymbol{\eta}_{t-m}^{(i)\ ^T})^T$. We define $\boldsymbol{\Sigma}_{\boldsymbol{\eta}\boldsymbol{\eta}_m,i}(k)$, $\boldsymbol{\Sigma}_{\boldsymbol{\eta},i}(k)$ and $\boldsymbol{\Sigma}_{\boldsymbol{\eta}_m,i}(k)$ analogously as the covariance matrices of the corresponding stochastic vectors. When $k=k_0$, the matrix admits the following representation:
Under the assumption that the idiosyncratic components are serially uncorrelated, it follows that $\textnormal{Cov}(\mbox{\boldmath$\varepsilon$}_{t}^{(i)}\mathbb{I}_{t,i}(k),\boldsymbol{\eta}_{t,m}^{(i)}\mathbb{I}_{t,i}(k))=0$ for $i=1,2$. Therefore, we have the simplified form of the matrix ${\mathbf M}_{i}(k_0)$:
where $\boldsymbol{\Sigma}_{{\mathbf f}\boldsymbol{\eta}_m,i}(k_0)=\textnormal{Cov}({\mathbf f}_{t}^{(i)}\mathbb{I}_{t,i}(k_0),\boldsymbol{\eta}_{t,m}^{(i)}\mathbb{I}_{t,i}(k_0))$.
From ((ref)), at the true change point $k_0$, the CCA matrix ${\mathbf M}_i(k_0)$ is driven solely by the common factor component, while the idiosyncratic component contributes no serial dependence under the serial-uncorrelated noise assumption. As a consequence, the rank of ${\mathbf M}_i(k_0)$ is determined by the factor dimension $r_i$, and the remaining eigenvalues associated with the orthogonal noise space are theoretically equal to zero. This structure provides the foundation for change-point detection: when the candidate point coincides with the true change-point location, the residual canonical dependence lying in the estimated noise subspace is minimized.
We first consider the case in which the regime-specific factor numbers $r_i$ in both regimes are known. For each regime $i \in \{1,2\}$, let $\lambda_{j}^{(i)}(k)$ for $j=1,2,\ldots,p,$ denote the $j$-th largest eigenvalue of ${\mathbf M}_i(k)$, and let ${\mathbf b}_{j}^{(i)}(k)$ be the corresponding unit eigenvector. Define the signal-subspace eigenvector matrix ${\mathbf P}^{(i)}(k)=({\mathbf b}_{1}^{(i)}(k),{\mathbf b}_{2}^{(i)}(k),..., {\mathbf b}_{r_i}^{(i)}(k))$ and the complementary noise-subspace eigenvector matrix ${\mathbf Q}^{(i)}(k)=({\mathbf b}_{r_{i}+1}^{(i)}(k),{\mathbf b}_{r_{i}+2}^{(i)}(k),\ldots,{\mathbf b}_{p}^{(i)}(k))$, such that $({\mathbf P}^{(i)}(k),{\mathbf Q}^{(i)}(k))$ forms a complete orthogonal basis of $\mathbb{R}^p$. From ((ref)), it is straightforward to verify that the factor loading space $\mathcal{M}({\mathbf L}_1^{(i)})$ is equivalent to the eigenspace $\mathcal{M}({\mathbf P}^{(i)}(k_0))$. Therefore, strict orthogonality holds between the loading space and the noise subspace, i.e., $\mathcal{M}({\mathbf L}_1^{(i)})\ \bot \ \mathcal{M}({\mathbf Q}^{(i)}(k_0))$. For any $i=1,2$, we obtain the following relation:
Furthermore, we have $\mbox{rank}({\mathbf M}_i(k_0))=r_i$, which implies that all subsequent eigenvalues are theoretically equal to zero, that is, $\lambda_{l}^{(i)}(k_0)=0$ for any $l=r_i+1,r_i+2,\ldots,p$. Thus, equation ((ref)) admits an alternative spectral representation:
Motivated by the projection of a second-order cross-moment matrix onto its orthogonal noise subspace in LiuZhang_2022, we define the following criterion for estimating the change point:
where the criterion function $G(k) \ge 0$ for all admissible $k$, and equality holds if and only if the candidate point matches the true change point $k=k_0$, as implied by ((ref)). Since ${\mathbf Q}^{(i)}(k)$ is constructed from the eigendecomposition of ${\mathbf M}_i(k)$ under the candidate partition, the criterion in ((ref)) is equivalently a candidate-specific residual eigenvalue criterion. It measures the deviation of ${\mathbf M}_i(k)$ from the ideal rank-$r_i$ canonical correlation structure through the magnitude of its noise-space eigenvalues. The normalization by $\|{\mathbf M}_{i}(k)\|_2$ is introduced to reduce the impact of regime-dependent scale variation across candidate partitions. In this way, the criterion function focuses on the relative magnitude of the residual component in the estimated noise subspace, rather than on the overall size of the canonical correlation matrix itself.
By spectral decomposition, the proposed criterion function admits a simpler representation in terms of eigenvalue ratios:
By ((ref)), it follows immediately that $G(k_0)=0$. The condition $G(k_0)=0$ characterizes the ideal low-rank canonical correlation structure under the true partition. To see the identification mechanism more clearly, consider a candidate split $k \neq k_0$. If $k < k_0$, the second segment contains observations generated from both regimes; if $k > k_0$, the first segment is similarly contaminated by observations from the two regimes. In either case, the regime-specific CCA matrix computed under the candidate partition is no longer generated by a single homogeneous transformed factor structure. Provided that the factor loading spaces, or the corresponding dynamic canonical correlation structures, differ sufficiently across the two regimes, such a mixed segment induces additional nonzero canonical-correlation components outside the regime-specific signal space. Consequently, the residual eigenvalues in the candidate-specific noise subspace become positive, leading to $G(k)>G(k_0)$. This explains why the proposed criterion identifies the true change point at the population level.
Furthermore, replacing the matrix $L_2$ norm in ((ref)) by the Frobenius norm yields,
By definition of the Frobenius norm, ((ref)) can be rewritten as the ratio of cumulative squared noise eigenvalues to cumulative squared total eigenvalues, as shown in ((ref)),
The corresponding change-point estimator under the Frobenius norm is defined analogously. The two criterion functions emphasize different aspects of the residual structure in the noise subspace. The $L_2$-norm-based criterion is more sensitive to the largest deviation from the ideal low-rank factor structure, whereas the Frobenius-norm-based criterion captures the cumulative contribution of all residual noise-space eigen-components. For this reason, both versions are retained and compared in the subsequent simulation study.
In empirical implementation, we compute the following sample covariance matrices before and after the candidate change point $k$, \[ \widehat\boldsymbol{\Sigma}_{{\mathbf y},1}(k)=\frac{1}{k}\sum_{t=1}^{k}({\mathbf y}_{t}-\bar{\mathbf y}^{(1)})({\mathbf y}_{t}-\bar{\mathbf y}^{(1)})^T,\quad\text{and}\quad \widehat\boldsymbol{\Sigma}_{{\mathbf y},2}(k)=\frac{1}{T-k}\sum_{t=k+1}^{T}({\mathbf y}_{t}-\bar{\mathbf y}^{(2)})({\mathbf y}_{t}-\bar{\mathbf y}^{(2)})^T, \] where $\bar{\mathbf y}^{(1)}=\dfrac{1}{k}\sum_{t=1}^{k}{\mathbf y}_{t}$ and $\bar{\mathbf y}^{(2)}=\dfrac{1}{T-k}\sum_{t=k+1}^{T}{\mathbf y}_{t}$. Analogous definitions are used for other sample covariance matrices. It follows that
For each $i=1,2$ and $j=1,2,\ldots,p$, let $\widehat\lambda_j^{(i)}(k)$ denote the estimate of the $j$-th largest eigenvalue of $\widehat{\mathbf M}_{i}(k)$, and let $\widehat{\mathbf b}_{j}^{(i)}(k)$ be the associated unit eigenvector. Then, if the number of factors $r_i$ is known, the estimator of the factor loading matrix in regime $i$ is given by $\widehat{\mathbf L}_1^{(i)}(k)=(\widehat{\mathbf b}_{1}^{(i)}(k),\widehat{\mathbf b}_{2}^{(i)}(k),...,\widehat{\mathbf b}_{r_i}^{(i)}(k))$. Finally, the factor vector can be estimated by ((ref)) as follows:
Next, we need to estimate the change point. Based on the preceding discussion, the sample criterion functions under the $L_2$ norm and the Frobenius norm are given by ((ref)) and ((ref)), respectively:
Since $G(k_0)=0$ characterizes the true change point, we estimate its location by minimizing the empirical criterion function over the admissible search domain:
Similarly, the change-point estimator under the Frobenius norm is defined by:
A key identifying condition in the proposed procedure is that the idiosyncratic components are serially uncorrelated. This assumption ensures that, at the true change point $k_0$, the residual canonical correlation matrix in the noise subspace has exact zero eigenvalues, which forms the basis of the criterion function $G(k)$. This assumption is stronger than the weak-dependence conditions often imposed in approximate factor models, where idiosyncratic components may exhibit limited temporal or cross-sectional dependence. Its role here is to deliver a sharp zero-eigenvalue structure in the population criterion, which in turn yields a tractable identification argument for the change point.
Estimating the regime-specific factor numbers $r_1$ and $r_2$ is an important component of the proposed procedure. We next present the one-step estimation method. For $i=1,2$, we can estimate the number of zero canonical correlations $v_i$ (and hence $r_i$) by testing the null hypothesis $H_0:\lambda_{p-v_{i}+1}=...=\lambda_p=0$ and $\lambda_{p-v_i}\neq0$ versus the alternative hypothesis $H_a:\lambda_{p-v_i}=0$. Let $n_1=\left \lfloor \gamma_{1}T \right \rfloor$ and $n_2=T-\left \lfloor \gamma_{2}T \right \rfloor$, which represent the sample sizes of the two boundary regions that exclude the possible change-point neighborhood, respectively. For simplicity, let $\widehat\lambda_{i,j}(\gamma_i)$ be the $j$-th largest eigenvalue of $\widehat{\mathbf M}_i(\gamma_i)$, where $\widehat{\mathbf M}_i(\gamma_i)$ is the canonical correlation matrix computed from the boundary interval $[1,\left \lfloor\gamma_1T\right \rfloor]$ or $[\left \lfloor\gamma_2T\right \rfloor,T]$, respectively. A test statistic for testing the hypothesis of zero canonical correlations is given by
where the statistic $S_{n_i}(v_i)$ asymptotically follows the chi-square distribution with degrees of freedom $df_i=v_i[(m-1)p+v_i]$ under the null hypothesis. The above procedure conducts the test sequentially from $v_i=1$ until the null hypothesis is accepted, leading to the estimator
where $\chi^2_{\alpha}(df_i)$ is the $\alpha$-quantile of the chi-square distribution.
Furthermore, we consider the transformed contribution ratio (TCR) estimator of Xia_TCR_2017 for comparison.
where $V_{i,j}=\sum_{s=j+1}^{p}\widehat\lambda_{i,s}(\gamma_i)$.
The use of boundary data serves as a stable starting point for factor-number estimation. Since the two boundary segments are far away from the central search region for the change point, they are less likely to be contaminated by observations near the structural change point and therefore more likely to reflect relatively homogeneous factor regimes. This feature makes them suitable for constructing initial estimators of the factor dimension. Nevertheless, relying only on boundary data does not fully exploit the information contained in the entire sample. Therefore, the one-step method for estimating the number of factors described above has certain limitations. Conversely, using a larger portion of the sample requires reasonably accurate knowledge of the change-point location, creating an intrinsic circular dependency between change-point detection and factor-number estimation. To address this issue, we propose the Alternating Iterative Estimation (AIE) algorithm, as shown in Algorithm (ref), where \(\mathcal F(a,b)\) denotes a generic factor-number estimation rule based on the two subsamples \([1,a]\) and \([b,T]\). The AIE algorithm alternates between change-point estimation and factor-number selection until the two estimates stabilize. We first obtain initial estimates of the factor numbers from the boundary data and then estimate the change point based on these preliminary values. Given the updated change-point estimate, the sample is repartitioned and the factor numbers in the two regimes are re-estimated. Repeating these two steps allows the procedure to gradually reduce the circular dependence between the unknown change-point location and the unknown factor numbers.
In this AIE algorithm, the initial estimators $\widehat{r}^{(0)}_1$ and $\widehat{r}^{(0)}_2$ are obtained by the same method (e.g., TCR or HT) using the boundary data, and the initial change-point estimate $\widehat{k}^{(0)}$ is then computed by treating $\widehat{r}^{(0)}_1$ and $\widehat{r}^{(0)}_2$ as the number of factors. The iteration terminates when the successive changes in the factor-number estimates and the change-point estimate fall below the thresholds ${{\varepsilon}}_2$ and ${{\varepsilon}}_1$, respectively, or when the number of iterations reaches $K_{\max}$. Upon termination, the resulting estimators are reported.
We next establish the theoretical properties of the proposed estimators. The analysis is conducted under an $\alpha$-mixing condition on the joint process $\{({\mathbf y}_t,{\mathbf f}_t^{(1)},{\mathbf f}_t^{(2)}):1\leq t\leq T\}$. The corresponding mixing coefficients are defined as follows:
where $\mathcal{F}_i^j$ is the $\sigma$-field generated by $\{({\mathbf y}_t,{\mathbf f}_t^{(1)},{\mathbf f}_t^{(2)}):i\leq t\leq j\}$.
Assumption (ref) is a standard regularity condition for weakly dependent time series, and its validity for VAR-type models is discussed in gaoetal2017. For each $i=1,2$, the moment conditions $\mathbb{E}(y_{jt})=0$ and $\mathbb{E}|y_{jt}|^{2\xi}\leq c$ in Assumption (ref) can be justified under suitable conditions on the series $\boldsymbol{\eta}_t^{(i)}$ and the corresponding loading matrix $\widetilde{\mathbf L}^{(i)}$ in model ((ref)). For instance, these conditions hold if $\mathbb{E}(\eta_{jt}^{(i)})=0$, $\mathbb{E}|\eta_{jt}^{(i)}|^{2\xi}< \infty$, and $\|\widetilde{\mathbf L}^{(i)}\|_{\infty}<\infty$, where $\|\cdot\|_{\infty}$ denotes the row norm of a matrix.
Under the normalization condition $\textnormal{Cov}(\boldsymbol{\eta}_t^{(i)})={\mathbf I}_p$ in Assumption (ref) and the stationarity condition in Assumption (ref), it follows as a mathematical corollary that the eigenvalues of $\boldsymbol{\Sigma}_{{\mathbf f}{\boldsymbol{\eta}_m,i}}(k_0)$ and $\boldsymbol{\Sigma}_{\boldsymbol{\eta}_m,i}(k_0)$ are bounded in Assumption (ref). The constants $\kappa_1$ and $\kappa_2$ of Assumption (ref) control the strength of the matrix $\widetilde{\mathbf L}^{(1)}$ and $\widetilde{\mathbf L}^{(2)}$, respectively. Assumption (ref) guarantees the identifiability of the eigenspaces associated with the positive eigenvalues of ${\mathbf M}_i(k_0)$ for $i=1,2$. This allows Theorem 2 below to establish the convergence rates of the estimator for ${\mathbf L}_1^{(i)}$ directly. Assumption (ref) requires the factor loading spaces before and after the change point to be sufficiently separated, thereby ensuring the identifiability of the structural change point.
Under a suitable orthogonal transformation of ${{\mathbf L}}_{j}^{(i)}$, Theorem (ref) establishes the convergence rate of the estimated loading matrices at the true change point when the numbers of factors are known. For each $i=1,2$, the choice of the semi-orthonormal loading matrix ${\mathbf L}_1^{(i)}$ in Model (ref) is generally not unique. Consequently, we measure the estimation error in the column space $\mathcal{M}({\mathbf L}_1^{(i)})$ rather than in the matrix ${\mathbf L}_1^{(i)}$ itself. This is because the column space $\mathcal{M}({{\bf L}}^{(i)}_1)$ is uniquely determined and remains invariant under different choices of ${\mathbf L}_1^{(i)}$. We therefore define the estimated and true factor loading spaces as $\mathcal{M}(\widehat{{\bf L}}^{(i)}_1)$ and $\mathcal{M}({{\bf L}}^{(i)}_1)$, respectively. The discrepancy measure proposed by panyao2008 is defined as
Note that $\widetilde{D}\left(\mathcal{M}(\widehat{{\bf L}}^{(i)}_1),\mathcal{M}({\bf L}^{(i)}_1)\right)$ takes values in $[0,1]$. In particular, it is equal to 0 if $\mathcal{M}(\widehat{{\bf L}}^{(i)}_1)=\mathcal{M}({\bf L}^{(i)}_1)$ and is equal to 1 if $\mathcal{M}(\widehat{{\bf L}}^{(i)}_1)$ and $\mathcal{M}({\bf L}^{(i)}_1)$ are orthogonal. Theorem (ref) below establishes the convergence of $\widetilde{D}\left(\mathcal{M}(\widehat{{\bf L}}^{(i)}_1),\mathcal{M}({\bf L}^{(i)}_1)\right)$ as $p$ and $T$ go to infinity.
Theorem (ref) shows that, when $p$ is fixed, the estimator $\mathcal{M}(\widehat{\bf L}_{1}^{(i)})$ converges to the true loading space $\mathcal{M}({\bf L}_{1}^{(i)})$ at the standard rate of $T^{-1/2}$. Furthermore, under the additional condition $\kappa_1\asymp\kappa_2\asymp p^{l}$ for some $0<l<1$ in Assumption (ref), the convergence rate can be refined to $O_p(p^{1-l}T^{-1/2})$. This makes explicit how the interplay between factor strength and dimensionality influences loading-space recovery.
This section investigates the finite-sample behavior of the proposed estimators through Monte Carlo experiments. The simulation design is chosen to reflect the theoretical features discussed in Section 3, with particular attention to the effects of sample size and cross-sectional dimension. We consider both the benchmark case in which the regime-specific factor numbers are known and the more realistic case in which they must be estimated from the data. The data-generating process follows Model ((ref)). We set the truncation parameters to $\gamma_1=0.1$ and $\gamma_2=0.9$. The true number of factors is set to $r_1=r_2=3$. We consider cross-sectional dimensions $p\in\{10,15,30,50\}$ and sample sizes $T\in\{200,500,1000\}$. For regimes $i=1,2$, the idiosyncratic components are generated from $\mbox{\boldmath$\varepsilon$}_{t}^{(i)}\sim \mathcal{N}({\mathbf 0},{\mathbf I}_{p-r_i})$. Within each regime, the latent factors ${\mathbf f}_{t}^{(i)}$ are generated from a VAR(1) process of the form ${\mathbf f}_{t}^{(i)} = \boldsymbol{\Phi}^{(i)} {\mathbf f}_{t-1}^{(i)} + {\mathbf u}_{t}^{(i)}$. The diagonal entries of the coefficient matrix are independently drawn from $U(0.2,0.9)$, and the innovation term ${\mathbf u}_{t}^{(i)} \sim \mathcal{N}({\mathbf 0},4{\mathbf I}_{r_i})$. We focus on the setting in which the factor strengths remain unchanged across the change point. The entries of the loading matrices $\widetilde{\mathbf L}^{(i)}$ are independently drawn from $U(-1,1)$. For each parameter configuration, we conduct 1,000 Monte Carlo replications to evaluate the performance of change-point estimation and loading-space recovery.
To investigate the estimation performance of the proposed procedure, we first consider the benchmark case in which the factor numbers are known, namely $r_1=r_2=3$. Table (ref) reports the average normalized absolute change-point estimation errors, $|\widehat{k}-k_0|/T$. As illustrated in Table (ref), for a fixed dimension $p$, the average normalized error $|\widehat{k}-k_0|/T$ decreases as the sample size $T$ increases. Conversely, for a fixed $T$, the estimation error tends to increase with the dimension $p$. The two norm-based criteria deliver broadly similar results, although their relative performance differs slightly across specific combinations of $p$ and $T$. For relatively small sample sizes ($T=200,500$), the estimation error exhibits an overall increasing trend as $p$ grows under both norm choices. This suggests that dimensionality has a more pronounced effect in smaller samples, whereas its impact becomes less substantial as the sample size increases. The two norm-based estimators yield broadly comparable performance, with only small differences manifesting under some combinations of $p$ and $T$. Overall, the simulation evidence supports the theoretical findings in Theorem (ref) and indicates that the proposed CCA-based procedure becomes more accurate in larger samples.
We then examine the discrepancy between the true loading space and its estimate under the same settings for $T$ and $p$ as in the previous analysis. In this simulation study, we set ${\bf L}^{(i)}_1=\widehat\boldsymbol{\Sigma}_{\mathbf y}^{-\frac{1}{2}}\widetilde{\bf L}^{(i)}_1$ for $i=1,2$, where $\widehat\boldsymbol{\Sigma}_{\mathbf y}$ denotes the sample covariance matrix of ${\mathbf y}_t$. Since $\widetilde{\bf L}^{(i)}$ is not constructed to be strictly orthogonal in our simulation design, we modify the discrepancy measure ((ref)) accordingly and use the criterion defined in ((ref)). Let ${\mathbf H}_j$ be a $p\times r_j$ matrix with $\mbox{rank}({\mathbf H}_j)=r_j$ for $j=1,2$. Define
This measure guarantees that $\bar{D}(\mathcal{M}({\mathbf H}_1),\mathcal{M}({\mathbf H}_2))=0$ if and only if one subspace is contained in the other (e.g. $\mathcal{M}({\mathbf H}_1)\subset \mathcal{M}({\mathbf H}_2)$), and $\bar{D}(\mathcal{M}({\mathbf H}_1),\mathcal{M}({\mathbf H}_2))=1$ if and only if $\mathcal{M}({\mathbf H}_1) \perp \mathcal{M}({\mathbf H}_2)$. Under the above definition, Table (ref) reports the average discrepancy $\frac{1}{2}\sum_{i=1}^{2}\bar{D}(\mathcal{M}(\widehat{{\mathbf L}}^{(i)}_{1}),\mathcal{M}({\mathbf L}^{(i)}_{1}))$ for the case where the factor numbers are known and equal to $r_1=r_2=3$. Table (ref) shows that the subspace discrepancy decreases as $T$ increases for fixed $p$; in contrast, for fixed $T$, loading-space estimation becomes more difficult as the dimension $p$ increases. This behavior also agrees with the theoretical analysis, since loading-space recovery relies on accurate estimation of the leading eigenspaces of the canonical correlation matrix, a task that becomes more challenging in relatively high-dimensional settings. First, the subspace is a high-dimensional object, and finite-sample noise can produce substantial discrepancy even when the change point is accurately estimated. Second, the difficulty is exacerbated by the simulation design, particularly when the factor strengths or eigen-gaps are small, which naturally leads to larger deviations in the estimated subspace. Importantly, these discrepancies do not undermine the practical utility of the proposed method for change-point detection. The criterion function relies on minimizing residual canonical correlation in the noise subspace, which remains effective even if the estimated loading spaces deviate from the true ones.
We subsequently evaluate change-point estimation performance when $r_{1}$ and $r_{2}$ are treated as unknown and must be estimated from the data. We compare the one-step procedure with the proposed AIE algorithm for $T\in \{200, 500, 1000\}$, using both the Transformed Contribution Ratio (TCR) and Hypothesis Testing (HT) methods to estimate the unknown factor numbers $r_1$ and $r_2$. As detailed in Table (ref), for a fixed sample size and under the $L_2$-norm criterion, the HT approach generally yields smaller change-point estimation errors than the TCR approach. By contrast, under the Frobenius-norm criterion, the TCR method performs better than the HT method. Across most specifications, the AIE procedure improves upon the corresponding one-step estimator, suggesting that iterative updating helps mitigate the dependence between change-point estimation and factor-number estimation. The improvement is consistent with the purpose of the iterative scheme. By re-estimating the factor numbers after updating the break location, the AIE procedure uses information from more homogeneous regime-specific subsamples and reduces the dependence on the initial boundary-based estimates. The results suggest that the iterative update can be beneficial when the factor numbers are unknown. Taken together, the results in Table (ref) suggest using HT with AIE under the $L_2$-norm criterion, and TCR with the AIE algorithm under the Frobenius-norm criterion.
Figure (ref) presents box plots of the estimators $\widehat{r}_{1}$ and $\widehat{r}_{2}$ under different methods, constraining the settings to $T \in \{500, 1000\}$, $p=10$, and $r_1=r_2=3$. The plots indicate that the one-step estimators exhibit substantial variability and noticeable bias, although such dispersion gradually diminishes as $T$ increases. By contrast, the AIE algorithm produces factor-number estimates that are much more tightly concentrated around the true values, regardless of whether TCR or HT is used. In addition to reducing bias, the boxplots indicate that the AIE procedure substantially improves the stability of factor-number estimation, as reflected by the noticeably smaller dispersion of the resulting estimators.
Finally, Table (ref) presents the discrepancy measure between the estimated and true factor loading spaces when the factor numbers are unknown. As $T$ increases, the discrepancy measure associated with the AIE algorithm remains comparable to that of the one-step estimator, indicating that the iterative estimation of the change point and factor numbers does not materially deteriorate loading-space estimation. Note that the loading-space discrepancies are relatively smaller than those in Table (ref). {This phenomenon can be explained by the fact that the estimated numbers of factors tend to exceed the true ones in this experiment. Consequently, the column spaces of the estimated loading matrices are enlarged and may provide a closer approximation to the true loading spaces, resulting in smaller loading-space errors.}
This section applies the proposed method to two empirical datasets. The first application uses intraday log excess returns of S&P 500 firms, and the second uses daily temperature records from U.S. monitoring stations. These applications are designed to illustrate how the proposed CCA-based procedure identifies interpretable changes in latent dynamic dependence structures in both financial and spatiotemporal panels.
{\bf Example 1 (Stock Returns).} The first application uses intraday log excess returns of S&P 500 stocks from Pelger_2020. The dataset covers the period from January 2004 to December 2016 and contains a balanced panel of 332 firms that are continuously observed over the full sample period. For illustration, we select 12 firms with relatively strong pairwise correlations and apply the proposed procedure to this subset. The correlation heatmap of these 12 firms is presented in Figure (ref). This subset is used mainly to provide a transparent visualization of the detected change point.
For this dataset ($p=12,T=3213$), we set the maximum lag order to $m=2$ and restrict the change-point search interval to $[0.2, 0.8]$. Figure (ref) plots the intraday return series of the 12 selected S&P 500 firms, with the vertical red line marking the estimated change point on November 21, 2008. Around the estimated change point, the return series exhibit a pronounced increase in volatility and more frequent extreme movements, particularly among real estate and housing-related firms such as BXP, SPG, UDR, HCP, PHM, KBH, and LEN, while the remaining stocks also display notable shifts in the magnitude of fluctuations. These synchronous changes across a broad set of firms indicate that the detected change point is unlikely to be driven by idiosyncratic shocks alone.
The estimated number of factors increases from $\widehat r_1=2$ before the change point to $\widehat r_2=4$ afterward. This increase suggests that the dimension of the latent dependence structure expanded during the crisis period. Economically, the estimated date lies within a period of severe financial stress in the United States, when concerns about major financial institutions, including Citigroup, intensified and shocks were transmitted from the financial sector to the broader economy. The result is therefore consistent with a reconfiguration of systematic risk following the Lehman Brothers collapse.
{\bf Example 2 (U.S. Temperatures).} We apply the proposed method to a 12-dimensional spatiotemporal dataset of daily surface temperatures. The data are obtained from the NOAA GSOD archive (2023); this dataset contains daily temperature records (in Fahrenheit) from January 1 to December 31, 2023, for 12 monitoring stations in the contiguous United States. Their geographic locations are shown in Figure (ref).
These monitoring stations are distributed across 11 states in the contiguous United States. More specifically, the stations cover a broad range of climatic regions across the Western, Midwestern, Eastern, and Southern United States. This broad geographic coverage provides reasonable spatial representation for the temperature dynamics under study. The raw temperature series for the individual stations are shown in Figure (ref). Notably, the series exhibit substantially greater variability from January to May than from June to August. This visual pattern suggests the possible presence of a change point during 2023.
After applying Seasonal-Trend decomposition using Loess (STL) to the original series, we set the maximum lag order to $m=2$ and restrict the change-point search interval to $[0.2, 0.8]$ in view of typical seasonal meteorological patterns. The estimated change point occurs on April 2, 2023. This estimate is consistent with the preliminary visual impression from the series and is marked in Figure (ref). It may be interpreted as reflecting a seasonal transition from the cold season to the warm season over mid-latitude North America. The estimation results indicate a reduction in factor numbers across the two regimes, with $\widehat{r}_{1}=4$ factors in the more heterogeneous cold-season regime and $\widehat{r}_{2}=2$ factors in the comparatively more homogeneous warm-season regime. The sign pattern of the first factor loadings changes little across the estimated change point, suggesting a broadly synchronous temperature anomaly pattern across locations, with the main difference lying in the loading magnitudes across the two regimes. The second post-change factor appears to capture a latitudinal contrast in temperature anomalies, roughly separated around the 42nd parallel north.
This paper proposes a CCA-based approach to the estimation of a single change point in high-dimensional transformed factor models. The method uses the low-rank structure of regime-specific canonical correlation matrices to construct an eigenvalue-based criterion for identifying changes in the latent factor structure. When the regime-specific factor numbers are unknown, an alternating iterative estimation procedure is introduced to update the factor numbers and the change point sequentially.
The theoretical analysis establishes convergence properties of the proposed estimators under mixing and moment conditions. The results show explicitly how the estimation accuracy depends on factor strength, sample size, and cross-sectional dimension. Monte Carlo experiments provide numerical evidence consistent with these theoretical findings. The empirical applications to intraday stock returns and daily temperature series further illustrate that the method can identify interpretable changes in latent dependence structures.
Several extensions remain possible. One direction is to allow for multiple change points in transformed factor models. Another is to relax the serial uncorrelatedness assumption on the idiosyncratic component and study settings with weakly dependent idiosyncratic errors. These extensions would broaden the applicability of the proposed framework to more general high-dimensional time-series environments.
The datasets used in the real-data analysis of this paper are all publicly available from open-access repositories. Their sources are listed below for the two empirical examples. Example 1 uses the intraday stock returns dataset, which is publicly available at \url{https://mpelger.people.stanford.edu/data-and-code}. Example 2 uses the global summary of the day dataset, which is accessible from the National Centers for Environmental Information (NCEI) archive at \url{https://www.ncei.noaa.gov/data/global-summary-of-the-day/archive/}. All datasets can be downloaded and used for research purposes subject to the respective access and usage policies of the hosting repositories.
{CRediT:} Lei Jia: Conceptualization, Methodology, Formal analysis, Investigation, Software, Validation, Data curation, Writing -- original draft, Visualization. Shouri Hu: Conceptualization, Methodology, Supervision, Writing -- review & editing. Zhaoxing Gao: Conceptualization, Methodology, Supervision, Project administration, Writing -- review & editing, Funding acquisition.
{No potential conflict of interest was reported by the authors.}
{This work was partially supported by the National Natural Science Foundation of China (NSFC) under Grant Nos. 72573029, U23A2064, 12201558, and 12401336.}