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.
87,452 characters · 9 sections · 95 citation commands
Taxonomy and Estimation of Multiple Breakpoints in High-Dimensional Factor Models
\def\spacingset#1{ {#1}} \spacingset{1}
\defcitealias{Duan2023}{Duan, Bai and Han}
\defcitealias{Baltagi2021}{Baltagi, Kao and Wang}
\if11 \fi
\if01 {
} \fi
{\it Keywords:} Quasi-Maximum Likelihood, Singularity, Consistency, Information criterion.
\spacingset{1.6}
Factor models have become a cornerstone in the analysis of large-scale economic and financial datasets characterized by both high cross-sectional and time dimensions. The models' versatility has led to widespread applications across economics and statistics, including the study of policy shocks (Stock2016, Stock2016; Han2025, Han2025), macroeconomic dynamics (Boivin2006, Boivin2006), asset pricing (Giglio2021, Giglio2021; Giglio2025, Giglio2025), portfolio theory (ross1976, ross1976; chamberlain1983, chamberlain1983; connor1986, connor1986, connor1988; Ding2021, Ding2021), causal inference ( Gobillon2016, Gobillon2016; Xu2017, Xu2017), and high-dimensional statistical machine learning (Fan2021, Fan2021). However, despite their wide applicability, real-world data are rarely stable over time. Structural breaks are ubiquitous, arising from policy interventions, financial crises, technological progress, or global events such as the COVID-19 pandemic. Importantly, they often embody the causal impact of interventions and shocks---policies may reshape economic structures, innovations may alter production and financial linkages, and major crises may transform underlying relationships. Identifying and modeling these breaks is therefore a powerful tool for causal inference, policy evaluation, and understanding how economic systems evolve in response to major events.
This paper develops a quasi-maximum likelihood (QML) approach for estimating breakpoints in high-dimensional factor models subject to multiple structural changes. We develop a unified theory that classifies structural breaks in factor loadings into two categories: singular changes and rotational changes. Our asymptotic results show that the proposed estimator precisely recovers the true breakpoints under singular changes, while for rotational changes the estimation error remains stochastically bounded. In addition, we propose an information criterion that consistently identifies the true number of breaks with probability approaching one as the sample size grows.
Previous studies (e.g., Han2015, Han2015) have shown that a single structural break in the factor loading matrix admits an observationally equivalent representation with time-invariant loadings and an expanded set of pseudo-factors. This expansion induces singularity in the covariance matrices of the pseudo-factors from either the pre-break or post-break subsample. \citetalias{Duan2023} (Duan2023, hereafter DBH) exploited this property and demonstrated that such singularity is key to establishing the consistency of the QML estimator in the single-break setting. However, the DBH approach is restricted by its single-break assumption. In comparison, our framework allows for multiple breaks and establishes a necessary and sufficient condition for break taxonomy. This is a critical advancement, as singularity alone no longer suffices to distinguish break types in multi-break settings. Moreover, we develop a novel proof strategy to overcome the technical challenges unique to multi-break setups and propose an information criterion that consistently selects the number of breaks.
Partly due to the simultaneous existence of both singular and rotational breaks, extending the analysis from a single-break to a multi-break setting presents several significant and non-intuitive challenges. First, unlike the single-break case, merely checking the singularity of the pseudo-factor covariance matrix within a given regime is insufficient to characterize the type of break in a multi-break environment. A break can be rotational even if neither the preceding nor the succeeding regime exhibits full column rank, a phenomenon that can occur due to singular breaks at other time points. Thus, a primary challenge lies in rigorously defining rotational and singular changes in factor loadings within a multi-break framework. Second, because multiple breakpoints are jointly estimated, the estimation errors of these individual breakpoint estimators become inherently interdependent, posing substantial technical challenges for asymptotic analysis. Third, given that the QML estimator exhibits distinct convergence rates for singular and rotational breaks, consistently determining the true number of breaks without prior knowledge of their specific types is challenging.
We address these challenges and make the following contributions. First, we establish a necessary and sufficient condition to precisely distinguish between singular and rotational structural breaks in factor loadings within a multi-break framework. Specifically, for each breakpoint, break types are identified by comparing the numbers of effective factors in the pre-break and post-break regimes with the number of pseudo-factors in the regime formed by pooling these two regimes. Second, to tackle the interdependence of estimation errors mentioned previously, we analyze the asymptotic properties of the $j$-th breakpoint by “conceptually” overfitting the remaining true breakpoints. This strategy requires us to compare the eigenvalues of the sample covariance matrix of the pseudo-factors in a regime against those from a subsample of that same regime. Third, leveraging this crucial result, we prove that the QML estimator identifies the true breakpoint with probability tending to one in large samples for singular breaks, and achieves $O_{p}(1)$ estimation errors for rotational breaks. Finally, we analyze how underfitting and overfitting the number of breaks affects the value of our QML objective function under both types of breaks. This analysis enables us to design an effective penalty function for our information criterion, ensuring its consistent selection of the true number of breaks without requiring prior knowledge of break types.
Moreover, we perform a series of Monte Carlo simulations to assess the finite sample properties of our estimator. The results demonstrate that our method can accurately estimate the break dates, even in the presence of small sample sizes, and that our information criterion consistently identifies the correct number of breaks. We further apply our approach to the FRED-MD dataset (McCracken2016, McCracken2016), detecting five breakpoints over the period from January 1959 to July 2024.
{\bf Links to the literature.} A growing body of literature investigates tests for structural breaks in high-dimensional factor models (e.g., Chen2014, Chen2014; Han2015, Han2015; Yamamoto2015, Yamamoto2015; Bai2024, Bai2024; fan2024tests, fan2024tests; Fan2024, Fan2024; Peng2025, Peng2025). Our study differs by focusing on the estimation of breakpoints rather than testing for their presence. Moreover, traditional literature on structural break estimation primarily addresses time series with small cross-sectional dimensions (e.g., Bai1998, Bai1998; Bai2000, Bai2000; Qu2007, Qu2007; Chan2014, Chan2014, Chan2021). In contrast, our factor model framework operates in a setting where both the cross-sectional ($N$) and time ($T$) dimensions are large.
Within this large $N$, large $T$ context, prior work has explored structural breaks in various settings. Qian2016, for instance, applied group fused Lasso to panel models with structural breaks, and Li2016 extended this to panel data with interactive fixed effects. While these studies considered large panels, their focus differed from ours, as they concentrated on breakpoints in regression coefficients rather than the high-dimensional factor loading matrix. For high-dimensional time series, Bai2010 developed a least squares (LS) estimator for breakpoints in the mean and variance, and Li2024 introduced a two-way MOSUM procedure with asymptotic properties for change-point estimators in trend functions. Our framework distinguishes itself from these by specifically addressing breakpoints in the loading coefficients of unobserved factors.
Regarding breakpoint estimation in factor models, earlier literature established the break fraction consistency, specifically that $\hat{\tau}-\tau_{0}=o_{p}(1)$, where $\hat{\tau}$ is the estimated breakpoint $\hat{k}$ divided by sample size $T$, and $\tau_{0}$ is the true fraction (e.g., Chen2015, Chen2015; Cheng2016, Cheng2016). More recent studies have provided theoretical results on the detailed rate of the error bound for the estimated breakpoint. For example, Barigozzi2018 proposed a wavelet-based method to estimate both the number and locations of structural breaks in common and idiosyncratic components. Barigozzi2024 further developed a MOSUM procedure for detecting multiple change points in high-dimensional factor models. Both studies showed that the estimation error of the breakpoint diverges slowly with sample size, implying a consistently estimated break fraction. Similarly, Baltagi2017 (Baltagi2017, Baltagi2021) developed an estimator that compares subsample second moments of factors estimated via full-sample principal components. Their breakpoint estimator exhibits a stochastically bounded estimation error, which also implies a consistently estimated break fraction.
In contrast to these approaches, our QML estimator achieves not only fraction-consistency but also point-consistency (sometimes referred to as super-consistency in the literature), meaning the estimated breakpoints equal the true ones with probability approaching one in large samples. To our knowledge, only a few studies have established point-consistency for breakpoint estimators of the loading matrix in factor models. Bai2020 demonstrated that the LS estimator is point-consistent under both large and small break settings. However, their framework is limited to a single break, whereas our methodology accommodates multiple breaks. Ma2018 proposed a fused Lasso estimator that is consistent for multiple breakpoints in factor models, under the assumption of a constant number of factors. Our work considers a more general and empirically relevant setting by allowing the number of effective factors to vary over time. Additionally, Ma2023 adopted a group fused Lasso approach to estimate changes in regression coefficients of one factor on the remaining factors. Their fused Lasso estimator exhibits stochastically bounded errors, which are then refined by LS to achieve point-consistency. Nonetheless, their method does not accommodate certain types of rotational changes (e.g., when the post-break loading matrix is rescaled by a constant), a feature that our Type 2 break in Section 3 does allow.
{\bf Plan of the Paper.} The rest of this paper is organized as follows. Section 2 introduces the factor model with multiple breaks on the factor loading matrix and describes the QML breakpoint estimator. Section 3 categorizes the two types of breakpoints. Section 4 outlines the assumptions for this model and presents the asymptotic theory. Section 5 investigates the finite-sample properties of the QML estimators through simulations. Section 6 provides an empirical study. Section 7 concludes.
{\bf Notation.} The following notation is used throughout the paper. Let $\rho_{i}(\mathbb{B})$ denote the $i$-th eigenvalue of an $n\times n$ symmetric matrix $\mathbb{B}$, and $\rho_{1}(\mathbb{B})\geq\rho_{2}(\mathbb{B})\geq\cdots\geq\rho_{n}(\mathbb{B})$. For an $m\times n$ real matrix $\mathbb{A}$, we denote its Frobenius norm as $\|\mathbb{A}\|=[tr(\mathbb{A}\mathbb{A}^{'})]^{1/2}$, and its adjoint matrix as $\mathbb{A}^{\#}$ when $m=n$. Let $\mathrm{Proj}(\mathbb{A}|\mathbb{Z})$ denote the projection of columns of matrix $\mathbb{A}$ onto the columns of matrix $\mathbb{Z}$. For $a,\ b>0$, $a\gtrsim b$ indicates that there exists a constant $c>0$ such that $a\geq cb$. We write $a\asymp b$ to mean that there exist constants $c_{1},c_{2}>0$ such that $c_{1}b\leq a\leq c_{2}b$. For a real number $x$, $[x]$ represents the integer part of $x$.
We consider the following factor model with $m_{0}$ common breaks at $k_{1}^{0}(T)<k_{2}^{0}(T)<\cdots<k_{m_{0}}^{0}(T)$ in the factor loadings for $i=1,\cdots,N$:
where $f_t$ is an $r$-dimensional vector of unobserved common factors; $r$ is the number of pseudo-factors in the entire sample; $k_j^0(T),j=1,\cdots,m_0$ are unknown break dates and $\min_{\{j=0,1,\cdots,m_0\}}k_{j+1}^0(T)-k_{j}^0(T)\geq \eta T$ for a constant $0<\eta<1$, with $k_{0}^0(T)=0$ and $k_{m_0+1}^0(T)=T$; $\lambda_{i,j}$ are $r\times1$ factor loading vectors for unit $i$ in the $j$th regime $t\in[k_{j-1}^0(T)+1,\,k_j^0(T)]$ for $j=1,\cdots,m_0+1$; and $e_{it}$ is the error term allowed to have serial and cross-sectional dependence as well as heteroskedasticity, and both $N$ and $T$ tend to infinity. $\tau_j^0=k_j^0(T)/T\in (0,1),j=1,\cdots,m_0$ are break fractions and are fixed constants. For notational simplicity, hereinafter, we suppress the dependence of $k_j^0$ on $T$. Note that the dimension of $f_{t}$ is the same as that of the pseudo-factors (defined below).
In vector form, model ((ref)) can be expressed as
where $x_{t}=[x_{1,t},\cdots,x_{N,t}]^{'}$, $e_{t}=[e_{1,t},\cdots,e_{N,t}]^{'}$, $\Lambda_{j}=[\lambda_{1,j},\cdots,\lambda_{N,j}]^{'},j=1,\cdots,m_{0}+1$.
For any $j=1,\cdots,m_0+1$, we define $$X_{j}=[x_{k_{j-1}^0+1},\cdots,x_{k_j^0}]^{'},~~F_{j}=[f_{k_{j-1}^0+1},\cdots,f_{k_j^0}]^{'},~~\mbox{\boldmath $e$}_{j}=[e_{k_{j-1}^0+1},\cdots,e_{k_j^0}]^{'}.$$ We rewrite ((ref)) using the following matrix representation:
where $F_j,j=1,\cdots,m_{0}+1$ has dimensions $(k_{j}^{0}-k_{j-1}^{0})\times r$ with full column rank for $j=1,\cdots,m_0+1$ and $\Lambda$ is an $N\times r$ matrix with full column rank. The loadings are modeled as $\Lambda_{j}=\Lambda B_{j}$, $j=1,\cdots,m_{0}+1$, where $B_{j}$ is an $r\times r$ matrix. Each $\Lambda_{j}$ has dimension $N\times r$. In this model, $$r_{j}=rank(B_{j})\leq r$$ denote the numbers of effective factors in regime $j$, i.e., the rank of the common component in regime $[k_{j-1}^0+1,k_j^0]$. In ((ref)), we define $$G=[g_1,\cdots,g_T]'$$ as the pseudo-factors because the final expression in ((ref)) provides an observationally equivalent representation, maintaining the structure of the loadings matrix \( \Lambda \) unchanged. More precisely, when the break is omitted during the estimation, the factors derived from a full-sample PCA correspond to the pseudo-factors \( G \) presented in ((ref)). It is well established that the presence of breaks can expand the factor space; consequently, for \( j = 1, \dots, m_{0}+1 \), we have \( r_j \leq r \), where \( \text{rank}(G) =\text{rank}([B_1, B_2, \dots, B_{m_0 + 1}])= r \). This representation allows for changes in factor loadings and the number of factors. Concrete constructions of $\Lambda$, $B_j$, $r$, and $r_j$ under different breakpoint types are given in Examples 1--3 of Section 3. We differentiate between effective factors and pseudo-factors. Effective factors are the factors required to capture common variation within a single, stable regime. In contrast, pseudo-factors are used to describe the factors (whose dimension may exceed that of effective factors due to breaks) within either the full sample (i.e., $G$) or a subsample consisting of consecutive regimes.
In what follows, we utilize the QML method to estimate the break dates for Model ((ref)), assuming the number of breaks is known. An information criterion for determining the number of breaks is also proposed in Section 4.
Let $\hat{G}=(\hat{g}_{1},...,\hat{g}_{T})^{\prime}$ denote the full-sample PCA estimator for $G$ subject to the normalization condition that $G^{'}G/T=I_r$ and $\Lambda^{'}\Lambda$ being diagonal. The number of pseudo-factors can be estimated following Bai2002. \footnote{Other methods for determining the number of factors include the eigenvalue ratio and growth ratio (Ahn2013, Ahn2013), a similar method using adjacent eigenvalue ratios proposed by LamYao2012, the empirical distribution estimator (Onatski2010, Onatski2010), and the bridge estimator (Caner2014, Caner2014).} For a given set of splitting points $(k_1,\cdots,k_{m_0})$, the QML objective function is expressed as follows:
where \[ \hat{\Sigma}(k_{\ell-1},k_{\ell})=\frac{1}{k_{\ell}-k_{\ell-1}}\sum\limits _{t=k_{\ell-1}+1}^{k_{\ell}}\hat{g}_{t}\hat{g}_{t}^{'},\] with $k_{0}=0$ and $k_{m_0+1}=T$. The QML estimator of the breakpoints for model ((ref)) is defined as follows:
where $K_\eta=\{(k_1,\cdots,k_{m_0}): 0=k_0<k_1<\cdots<k_{m_0}<k_{m_0+1}=T, k_j-k_{j-1} \geq T\eta,j=1,\cdots,m_0+1\}$ and $\eta \in (0,1)$ is a small positive constant.
The estimator can be viewed as a two-step procedure: first estimate $\hat G$ by full-sample PCA, and then locate the breakpoints by minimizing the QML objective function. Unlike classical methods, our log-determinant objective exploits the near-singularity of $\hat{\Sigma}(k_{\ell-1}, k_{\ell})$ when the segment boundaries approach true Type 1 breakpoints. This sensitivity to the rank-deficient structure enables the faster convergence rates established in Theorem 2, outperforming existing methods such as Barigozzi2018, Barigozzi2024, Baltagi2017, Baltagi2021.
Once the breakpoints are estimated, one can apply split-sample PCA to consistently estimate the loadings and factors for each regime. Using data within the interval $[\hat{k}_{j-1}+1, \hat{k}_j]$, one can recover the regime-specific $r_j$-dimensional factors $\hat{f}_t$ and loadings $\hat{\Lambda}_j$, subject to a local $r_j \times r_j$ rotation $\hat{H}_j$.
In the next section, we classify breakpoints into two types based on the rank properties of the loading matrices across regimes.
The representation in ((ref)) is convenient for theoretical analysis because one can control break types by setting $B_j$ for $j=1,\cdots,m_0+1$. To further illustrate the generality of model ((ref)), we present two types of breakpoint structures within the combined time interval \([k_{j-1}^0+1, k_{j+1}^0]\) (combined regime) by varying the ranks of \( B_j \) and \( B_{j+1} \). We define the number of pseudo-factors in the combined regime $[k_{j-1}^{0}+1,k_{j+1}^{0}]$ as \footnote{This is because
}\[ r_{j,j+1}=\text{rank}([B_{j}F_{j}',B_{j+1}F_{j+1}']')=\text{rank}([B_{j},B_{j+1}]).\] Next, we show that there are two types of breakpoints.
Type 1. $r_{j,j+1}>\min\{r_{j,}r_{j+1}\}$. In this case, the number of pseudo-factors in the combined regime $[k_{j-1}^{0}+1,k_{j+1}^{0}]$ exceeds the number of effective factors in at least one of the two regimes: $[k_{j-1}^{0}+1,k_{j}^{0}]$ and $[k_{j}^{0}+1,k_{j+1}^{0}]$. We refer to this as {\bf a singular change}.
If $r_{j,j+1}>\max\{r_{j},r_{j+1}\}$, the dimension of the factor space over the combined regime $[k_{j-1}^{0}+1,k_{j+1}^{0}]$ increases due to the structural break. This occurs when the factor loading matrices from regime $[k_{j-1}^{0}+1,k_{j}^{0}]$ and regime $[k_{j}^{0}+1,k_{j+1}^{0}]$ jointly span a larger space than either individual matrix. This space augmentation is analogous to a regression where a break in coefficients increases the number of effective regressors; specifically, the interaction of a break dummy with original regressors leads to an expanded regressor set.
Type 2. $r_{j,j+1}=\min\{r_j,r_{j+1}\}$. In this case, the number of pseudo-factors in the combined regime \([k_{j-1}^0+1, k_{j+1}^0]\) is the same as the number of effective factors in each regime, and we refer to this as {\bf a rotational change}. Specifically, there exists a nonsingular $r\times r$ matrix $\bf R$ such that
then
which satisfies $r_{j,j+1}=r_j=r_{j+1}=rank(B_j)$. When there is one break ($m_0=1$) and $B_1$ is nonsingular, it is analogous to the rotational change for the single break setup.
These two types of breaks have distinct economic interpretations. Rotational (Type 2) breaks can be interpreted as observationally equivalent to a change in the second moment of the existing factors, while leaving the loading space unchanged. In contrast, singular (Type 1) breaks represent “heterogeneous” shifts where the column space of the loading matrix changes. This corresponds to more profound structural transformations, including the emergence of new factors or the disappearance of existing ones.
The preceding discussion naturally raises the question of how to distinguish between the two types of changes in multi-breakpoint settings, particularly when the singularity of $B_{j}$ and $B_{j+1}$ (the DBH rule in a single-break setting) is no longer sufficient for distinction. The answer is to use the definition $r_{j}=r_{j+1}=r_{j,j+1}$ to determine whether the $j$-th breakpoint is rotational. Evidently, if Equation ((ref)) holds, then $r_{j}=r_{j+1}=r_{j,j+1}$ must also hold, implying that the $j$-th breakpoint is rotational. Conversely, another natural question arises: does $r_{j}=r_{j+1}=r_{j,j+1}$ always imply the existence of a nonsingular $R$ matrix satisfying Equation ((ref))? This is affirmed by Proposition (ref), which demonstrates that Equation ((ref)) is a necessary and sufficient condition for defining a Type 2 break. This proposition also explains why Type 2 is named a “rotational change.”
Proposition (ref) clarifies that the condition $B_{j}{\bf R}=B_{j+1}$ is equivalent to the definition of a Type 2 breakpoint (i.e., $r_{j}=r_{j+1}=r_{j,j+1}$). In practice, researchers do not need to directly verify the existence of ${\bf R}$ such that $B_{j}{\bf R}=B_{j+1}$ is satisfied. Instead, based on the conclusion of Proposition (ref), the breakpoint type can be readily determined by simply comparing the magnitudes of $r_{j}$, $r_{j+1}$, and $r_{j,j+1}$. If these quantities are equal, it indicates a Type 2 change; otherwise, it is a Type 1 change.
For Type 1 breakpoints, we will show that the QML estimators converge to the true breakpoints with probability approaching one. In contrast, for Type 2 breakpoints, the difference between QML estimators and the true change points is stochastically bounded. The specific asymptotic theory is detailed in Section 4.
The assumptions are as follows:
Assumptions (ref)--(ref) are standard in the factor model literature. Assumption (ref) (i) is similar to Assumption A in Bai2003 and (ii) imposes restrictions on $B_j$ and $B_{j-1}$ to ensure the identification of the breakpoint. Assumption (ref) (ii) rules out the case that $B_j=-B_{j-1}$ since the objective function ((ref)) cannot identify the difference between $\hat \Sigma_{j-1}$ and $\hat \Sigma_{j}$. Our analysis focuses on the case of fixed-magnitude breaks to facilitate the development of the statistical properties of near-zero eigenvalues established in Proposition 2. Assumption (ref) (iii) implies that $\Sigma_G$ is positive definite. Assumption (ref) is similar to Assumption B of Bai2003. Assumption (ref) allows for weakly correlated error terms in both time and cross-sectional dimensions. Assumption (ref) means that the factors and idiosyncratic errors are allowed to be weakly dependent within each regime. Assumption (ref) corresponds to Assumption G in Bai2003. Assumption (ref) corresponds to Assumption 7 of Baltagi2021, which allows the H\'{a}jek--R\'{e}nyi inequality to be applicable to the second moment process of the factors.
Assumptions (ref) and (ref) are similar to Assumptions 7 and 8 of DBH. The terms within the squared norm in Assumption (ref) are standard in the factor model literature and are typically $O_p(1)$ (e.g., Assumptions C(5), F(1), and F(3) in Bai2003). Our assumption strengthens these traditional conditions by requiring that the averages of these squared terms remain uniformly $O_p(1)$ over sub-intervals with length exceeding $(\log T)^2$. We provide sufficient conditions for this uniformity in Section C of the Supplementary Appendix. Assumption (ref) concerns the partial sums of the mean-zero processes $Z_{N,t}^{(1)}$ and $Z_{N,t}^{(2)}$ over candidate sub-intervals $(s, k]$, where $Z_{N,t}^{(1)}=\frac{1}{\sqrt N}\sum_{i=1}^{N}\lambda_i f_t'e_{it}$ and $Z_{N,t}^{(2)}=\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{N}\lambda_i\lambda_j'\bigl(e_{it}e_{jt}-E[e_{it}e_{jt}]\bigr)$. For $k-s \ge \eta T$, the required $O_p(1)$ bound is a direct consequence of a Functional Central Limit Theorem. For shorter intervals where $k-s \ge (\log T)^2$, the relevant upper envelope involves a $\log\log$ factor, consistent with the Law of the Iterated Logarithm. Detailed derivations for a benchmark example illustrating the feasibility of Assumption (ref) are provided in Section C of the Supplementary Appendix.
Assumption (ref) extends Assumption F3 in Bai2003 and ensures that $\rho_{r}(\hat{\Sigma}(k,s))\asymp\frac{1}{N}$ w.p.a.1 if the interval $[k+1,s]$ belongs to $j$-th regime and $B_{j}$ is singular. Assumption (ref) is critical when $B_j$ or $B_{j+1}$ is singular. Suppose $B_{j+1}$ is singular with $r_{j+1} < r$. Under the true partition, the pseudo-factor covariance matrix for $(k_j^0, k_{j+1}^0]$ has rank $r_{j+1}$. Shifting the split point from $k_j^0$ to $k_j^0-1$ mistakenly incorporates $B_j f_{k_j^0}$ from regime $j$ into the misaligned segment $(k_j^0-1, k_{j+1}^0]$. Assumption (ref) requires $B_j f_{k_j^0} \notin \mathrm{Col}(B_{j+1})$, ensuring that the misaligned segment's rank exceeds $r_{j+1}$. In finite samples, where estimated factors are used to compute $\hat{\Sigma}(k_{\ell-1}, k_\ell)$, this rank inflation translates into an increase in the objective function, facilitating precise breakpoint identification. Conversely, if $B_j f_{k_j^0} \in \mathrm{Col}(B_{j+1})$ (e.g., $f_{k_j^0} = 0$), misspecifying the breakpoint as $k_j^0-1$ fails to increase the rank of segment $(k_j^0-1, k_{j+1}^0]$, in which case the breakpoint cannot be precisely estimated. The case of $B_{j}$ being singular is similar by symmetry.
In this subsection, we establish the asymptotic properties of the QML estimators. The analysis of these properties is significantly more challenging in a multi-break setting than in a single-break setting. As discussed in the Introduction, the estimation error of any estimated breakpoint is generally interdependent with that of others. This interdependence is further complicated when both Type 1 and Type 2 breaks are present. Furthermore, since the true endpoints of a regime are unknown, an estimated segment $[\hat{k}_{j}+1,\hat{k}_{j+1}]$ may have no overlap with the true interval $[k_{j}^{0}+1,k_{j+1}^{0}]$. This contrasts sharply with the single-break setting, where the estimated segments are guaranteed to overlap with the true regimes because $k_{0}^{0}=0$ and $k_{m_{0}+1}^{0}=k_{2}^{0}=T$ are known.
To address these challenges, we introduce a strategy of “conceptually” overfitting the number of breaks to isolate the $j$-th breakpoint of interest. Specifically, assuming $\hat{k}_{j}\ne k_{j}^{0}$, we conceptually insert all true breakpoints, except for the $j$-th one, into the existing set of estimated breakpoints. After inserting these breakpoints, all resulting segments (derived from the estimated and inserted breakpoints) contain no true breakpoint, except for one segment that contains the $j$-th true breakpoint. This enables us to derive an upper bound for $|\hat{k}_{j}-k_{j}^{0}|$.
While this strategy helps isolate the target breakpoint and analyze its asymptotic properties, the conceptually overfitting strategy introduces two technical challenges. First, an inserted true breakpoint may be very close to an estimated breakpoint, leading to a segment with fewer than $r$ observations. In such a segment, the pseudo-factors will have an exactly singular covariance matrix. We address this issue by merging such a short segment with an adjacent regime whose length is necessarily proportional to $T$. More details are provided in the text preceding the proof of Theorem 1 in the Appendix. Second, overfitting the breaks inevitably cuts a true regime into sub-regimes. When the pseudo-factors in the $j$-th regime are not of full column rank (i.e., $B_{j}$ is singular), it is crucial and challenging to compare the nearly zero eigenvalues of the sample covariance matrices of pseudo-factors in the true regimes versus those in the sub-regimes created by conceptually overfitting. Proposition (ref) provides a detailed rate for the difference between these nearly-zero eigenvalues. This proposition plays a crucial role in establishing the point-consistency of the QML estimator and the selection consistency of our proposed information criterion.
To compare these nearly-zero eigenvalues, we consider the case of conceptually overfitting a point $\tilde{k}_{j}$ within the regime $[k_{j-1}^{0}+R,\ k_{j}^{0}-R]$, where $R=(\log T)^{2}$. Recall that $\hat{\Sigma}(k,s)$ denotes the estimated sample covariance matrix of $\hat{g}_t$ for $k+1 \leq t \leq s$. Specifically, \[ \hat{\Sigma}(k_{j-1}^{0},k_{j}^{0})=\hat{G}_{j}^{\prime}\hat{G}_{j}/(k_{j}^{0}-k_{j-1}^{0}),\ \ \hat{\Sigma}(k_{j-1}^{0},\tilde{k}_{j})=\hat{G}_{j,1}^{\prime}\hat{G}_{j,1}/(\tilde{k}_{j}-k_{j-1}^{0}),\ \ \hat{\Sigma}(\tilde{k}_{j},k_{j}^{0})=\hat{G}_{j,2}^{\prime}\hat{G}_{j,2}/(k_{j}^{0}-\tilde{k}_{j}),\] where $\hat{G}_{j}=[\hat{g}_{k_{j-1}^{0}+1},...,\hat{g}_{k_{j}^{0}}]'$, $\hat{G}_{j,1}=[\hat{g}_{k_{j-1}^{0}+1},...,\hat{g}_{\tilde{k}_{j}}]'$, $\hat{G}_{j,2}=[\hat{g}_{\tilde{k}_{j}+1},...,\hat{g}_{k_{j}^{0}}]'$. The following proposition then establishes the convergence rate of the difference between the eigenvalues of $\hat{\Sigma}(k_{j-1}^{0},k_{j}^{0})$ and $\hat{\Sigma}(k_{j-1}^{0},\tilde{k}_{j})$, which is crucial for the subsequent theoretical analysis when $B_{j}$ is singular.
In Proposition (ref), the convergence rate is determined by $N$ and the length of the subsample, $\tilde{k}_{j}-k_{j-1}^{0}$ (or $k_{j}^{0}-\tilde{k}_{j}$). Notably, the distance (in terms of Frobenius norm) between $\hat{\Sigma}(k_{j-1}^{0},k_{j}^{0})$ and $\hat{\Sigma}(k_{j-1}^{0},\tilde{k}_{j})$ is no smaller than the order $T^{-1/2}$, while the small eigenvalues of $\hat{\Sigma}(k_{j-1}^{0},k_{j}^{0})$ converge to zero at a rate of $T^{-1}$. This disparity implies that we cannot obtain the rate of $\rho_{\ell}\left(\hat{\Sigma}(k_{j-1}^{0},k_{j}^{0})\right)-\rho_{\ell}\left(\hat{\Sigma}(k_{j-1}^{0},\tilde{k}_{j})\right)$ for $\ell=r_{j}+1,\cdots,r$ using the continuous mapping theorem. This is because the distance between these two matrices is actually much larger than the magnitude of their smallest eigenvalues. To address this, Lemmas 4--6 in the Appendix develop essential technical tools that facilitate the proof of the detailed rate for $\rho_{\ell}\left(\hat{\Sigma}(k_{j-1}^{0},k_{j}^{0})\right)-\rho_{\ell}\left(\hat{\Sigma}(k_{j-1}^{0},\tilde{k}_{j})\right)$ as established in Proposition (ref). These lemmas may also be of independent interest.
Theorem (ref) provides an initial upper bound on the estimation error of $\hat{k}_{j}$ for both types of breaks. Part (i) indicates that the estimation error for Type 1 breaks is bounded by $\sqrt{T}$ with high probability in large samples. Moreover, Part (ii) shows that for a Type 2 break, the upper bound is $T^{\alpha}$, where $\alpha$ is some constant slightly larger than $1/2$. These results collectively imply fraction-consistency, i.e., $\hat \tau_{j}-\tau_{j}^{0}=o_{p}(1)$, where $\hat \tau_{j}=\hat{k}_{j}/T$ is the estimated break fraction. Building upon Theorem (ref), the subsequent theorem further refines the convergence rate of $\hat{k}_{j}-k_{j}^{0}$.
Theorem (ref) (i) demonstrates that the estimated change points converge to the true breakpoints w.p.a.1 for a Type 1 break. This result is stronger than that of Baltagi2021, who established that the distance between estimated and true breakpoints is $O_{p}(1)$ for Type 1 changes. It is important to note that Type 1 breaks incorporate scenarios of emergence and disappearance of factors, which correspond to $r_{j,j+1}=\max\{r_{j},r_{j+1}\}$. Our QML estimator proves to be consistent when the number of factors changes across regimes, whereas Ma2018 excluded this case by assumption. Furthermore, Theorem (ref) (ii) suggests that the discrepancy between the QML estimators and the true change points is stochastically bounded for Type 2 changes.
{\bf Determining the number of breaks}
We now discuss the choice of the number of breaks $m$, which is an important issue when the objective function is used in practice. We propose selecting $m$ to minimize the following information criterion:
where $(\hat{k}_1, \dots, \hat{k}_m)$ are obtained by solving ((ref)) for each given $m$, and $\hat{\rho}$ is the radius of the estimated AR(1) coefficient matrix, obtained by fitting a VAR(1) model to $\hat{g}_{t}$, although the actual $g_t$ is not required to be an AR(1) process.
Note that, for any fixed $m$, the penalty term in Eq. ((ref)) is constant with respect to the breakpoint locations. Hence, the joint search over $m$ and the breakpoint locations is numerically equivalent to the following two-step procedure: first, for each $m$, determine the optimal break locations by minimizing the QML objective function in Eq. ((ref)); second, substitute these minimizers into Eq. ((ref)) and select the value of $m$ that yields the smallest criterion value. In practice, for each candidate $m$, the minimization of $U_{NT}$ is implemented by the dynamic programming algorithm in Section 2.2.
Let $m_{max}$ be a bounded integer such that $m_{0}<m_{max}$ and $\hat{m}=\arg\min_{0\leq m\leq m_{max}}IC(m)$.
In this section, we consider DGPs corresponding to different break types to evaluate the finite sample performance of the QML estimator. We compare the QML estimator with two other estimators. As shown below, $\hat{k}_{BKW}$ is the estimator proposed by \citetalias{Baltagi2021} (Baltagi2021, BKW hereafter); $\hat{k}_{MS}$ is the estimator proposed by Ma2018 (Ma2018, MS hereafter); and $\hat{k}_{QML}$ is the QML estimator. We explore the performance of these estimators \footnote{We also evaluate the estimators of Ma2023, which perform similarly to those of Ma2018. Due to space constraints, a comprehensive comparison of these results is provided in Appendix D of the supplementary material.} under various DGPs, incorporating different types of changes. We calculate the RMSE of these change point estimators $\hat{k}_{BKW}$, $\hat{k}_{MS}$, and $\hat{k}_{QML}$, and each experiment is repeated 1000 times, where RMSE$_{i}=\sqrt{\frac{1}{1000}\sum\limits _{s=1}^{1000}(\hat{k}_{i,s}-k_{i}^{0})^{2}}$ for $i=1,\cdots,m_{0}$. Note that the MS method sometimes detects more or fewer than $m_{0}$ breaks. To ensure comparability, the RMSE for the MS estimator is calculated exclusively from instances where the method correctly identifies the true number of breaks.
Each factor is generated by the following AR(1) process:
where $u_{t}=(u_{t,1},\cdots,u_{t,r_{0}})^{'}$ is i.i.d. $N(0,I_{r_{0}})$ for $t=2,\cdots,T$ and $f_{1}=(f_{1,1},\cdots,f_{1,r_{0}})^{'}$ is drawn from $N(0,\frac{1}{1-\rho^{2}}I_{r_{0}})$. The scalar $\rho$ captures the serial correlation of factors, and the idiosyncratic errors are generated by
where $v_{t}=(v_{1,t},\cdots,v_{N,t})^{'}$ is i.i.d. $N(0,\Omega)$ for $t=2,\cdots,T$ and $e_{1}=(e_{1,1},\cdots,e_{N,1})^{'}$ is $N(0,\frac{1}{(1-\alpha^{2})}\Omega)$. The scalar $\alpha$ captures the serial correlation of the idiosyncratic errors, and $\Omega$ is generated as $\Omega_{ij}=\beta^{|i-j|}$ so that $\beta$ captures the degree of cross-sectional dependence of the idiosyncratic errors. In addition, $u_{t}$ and $v_{t}$ are mutually independent for all values of $t$. We consider the following DGPs for factor loadings and investigate the performance of the QML estimator for the different types of breaks discussed in Section 3.
DGP 1. We first consider the case in which $m_0=2$ and $r_0=3$. We set $k_1^0=[0.3T]$ and $k_2^0=[0.7T]$, where $[\cdot]$ denotes the rounding operation. This means that there are two breakpoints and three regimes.
DGP 1.A In this setup, we adopt the DGP from MS. Each of the three regimes has $r_{0}$ effective factors, i.e., $r_{1}=r_{2}=r_{3}=3$. The factor loading matrices on the effective factors are generated as follows: \[ \breve{\Lambda}_{1}=(\breve{\lambda}_{1,1},\ldots,\breve{\lambda}_{N,1})',\quad\breve{\lambda}_{i,1}\;\overset{\text{i.i.d.}}{\sim}\; N\bigl((0.5b,\,0.5b,\,0.5b)^{\prime},\,\tfrac{1}{r_{0}}I_{r_{0}}\bigr),\] \[ \breve{\Lambda}_{2}=(\breve{\lambda}_{1,2},\ldots,\breve{\lambda}_{N,2})',\quad\breve{\lambda}_{i,2}\;\overset{\text{i.i.d.}}{\sim}\; N\bigl((b,\, b,\, b)^{\prime},\,\tfrac{1}{r_{0}}I_{r_{0}}\bigr),\] \[ \breve{\Lambda}_{3}=(\breve{\lambda}_{1,3},\ldots,\breve{\lambda}_{N,3})',\quad\breve{\lambda}_{i,3}\;\overset{\text{i.i.d.}}{\sim}\; N\bigl((1.5b,\,1.5b,\,1.5b)^{\prime},\,\tfrac{1}{r_{0}}I_{r_{0}}\bigr),\] The $\Lambda$ matrix in ((ref)) can be written as $\Lambda\;=\;\bigl[\breve{\Lambda}_{1},\breve{\Lambda}_{2},\breve{\Lambda}_{3}\bigr]_{N\times3r_{0}}$. We define the three $3r_{0}\times 3r_{0}$ block-selector matrices\[ B_{1}=\operatorname{diag}\!\bigl(I_{r_{0}},\,0_{r_{0}\times r_{0}},\,0_{r_{0}\times r_{0}}\bigr),\quad B_{2}=\operatorname{diag}\!\bigl(0_{r_{0}\times r_{0}},\,I_{r_{0}},\,0_{r_{0}\times r_{0}}\bigr),\quad B_{3}=\operatorname{diag}\!\bigl(0_{r_{0}\times r_{0}},\,0_{r_{0}\times r_{0}},\,I_{r_{0}}\bigr), \] so that in regime $j$ we set $\Lambda_{j}=\Lambda\, B_{j}$, which yields three effective factors per regime and a total of $r=3r_{0}=9$ pseudo-factors over the full sample. This corresponds to a Type 1 breakpoint scenario with $r_{j,j+1}=r_{j}+r_{j+1}$, where loadings in adjacent regimes are independent. By varying $b$ we can control the magnitude of the break.
DGP 1.B We set $\lambda_{i}$ to be i.i.d. $N(0,\frac{1}{r}I_{r})$ across $i$, and define $\Lambda=(\lambda_{1},\ldots,\lambda_{N})'$. In the first regime, $\Lambda_{1}=(\lambda_{1,1},\cdots,\lambda_{N,1})^{'}=\Lambda B_{1}$ with $B_{1}=[1,0,0;0,1,0;0,0,0]$. In the second regime, $\Lambda_{2}=(\lambda_{1,2},\cdots,\lambda_{N,2})^{'}=\Lambda B_{2}$ with $B_{2}=[1,0,0;0,0,0;0,0,1]$. In the third regime, $\Lambda_{3}=(\lambda_{1,3},\cdots,\lambda_{N,3})^{'}=\Lambda B_{3}$ with $B_{3}=[0,0,0;0,1,0;0,0,1]$. For the full sample, the number of pseudo-factors is $r=3$. In each of the three regimes, the number of effective factors is 2, as $\text{rank}(B_{1})=\text{rank}(B_{2})=\text{rank}(B_{3})=2$. This setting corresponds to a Type 1 break in Section 3 with $r_{j,j+1}>\max\{r_{j},r_{j+1}\}$.
DGP 1.C The generation of $\Lambda$ remains the same as in DGP 1.B. For $B_j$, we set $B_{1}=[1,0,0;0,1,0;0,0,1]$, $B_{2}=[1,0,0;0,1,0;0,0,0]$ and $B_{3}=[0,0,0;0,0,0;0,0,1]$. The number of pseudo-factors for the full sample is $r=3$, with the numbers of effective factors in the three regimes being 3, 2, and 1, respectively. In the first regime, all three factors are present. In the second regime, only the first two factors are present and in the third regime, only the third factor is present.
DGP 1.D The generation of $\Lambda$ remains the same as in DGP 1.B. Let $\Lambda_{1}=\Lambda$, $\Lambda_{2}=2\Lambda$ and $\Lambda_{3}=\Lambda$. In this case, the number of factors remains constant across regimes, and the changes in factor loading are purely proportional, corresponding to a setting with two Type 2 breaks.
DGP 1.E In this setup, we consider a case involving both types of changes. The generation of $\Lambda$ remains the same as in DGP 1.B. We set $B_{1}=[1,0,0;0,1,0;0,0,0]$, $B_{2}=[2,x,y;0,2,z;0,0,0]$ where $x$, $y$, $z$ are i.i.d. $N(0,1)$ and $B_{3}=[0,0,0;0,0,0;0,0,1]$. By construction, there exists a nonsingular matrix ${\bf R}$ such that $B_{1}{\bf R}=B_{2}$, so the first breakpoint is Type 2 with $r=3>r_{1,2}=r_{1}=r_{2}=2$. The second breakpoint, however, is a Type 1 change and is responsible for the rank deficiency of the pseudo-factors in the first two regimes.
Since the computation of $\hat{k}_{MS}$ requires the number of effective factors ($r_{0}$), we use $\hat{r}=r_{0}$ for it. Moreover, as $\hat{k}_{QML}$ and $\hat{k}_{BKW}$ require the number of pseudo-factors ($r$), we set $\hat{r}=r$ for them. This allows us to focus on the estimation errors of the breakpoint estimators themselves, rather than errors resulting from misspecifying the number of factors.
Table (ref) reports the RMSEs of the estimated break dates under DGP 1.A--1.E for $(\rho,\alpha,\beta)=(0.7,0.3,0.3)$. For DGP (ref).A (both $b=1$ and $b=0$), the QML estimator yields very small RMSE in small samples, and the RMSE quickly decreases as $(N,T)$ increase; BKW and MS also improve with sample size but are uniformly less accurate than QML. In DGP (ref).B, QML again achieves the smallest RMSE and continues to improve with $(N,T)$, whereas BKW exhibits non-vanishing (stochastically bounded) errors as the sample size grows; MS improves but remains less accurate than QML. In DGP (ref).C, where factors disappear and the factor dimension changes across regimes, both competing methods deteriorate---especially MS, which relies on a constant factor dimension---while QML remains stable and its RMSE decreases with $(N,T)$. For DGP (ref).D (rotational change), QML delivers small and bounded estimation errors, consistent with Theorem (ref)(ii), and it clearly outperforms BKW; MS again performs poorly and does not improve with sample size. Finally, for the mixed design in DGP (ref).E (Type 2 followed by Type 1), QML yields bounded RMSE for the Type 2 breakpoint and near-zero RMSE for the Type 1 breakpoint, whereas both BKW and MS have substantially larger errors. Overall, the simulation evidence is consistent with Theorem (ref): QML is consistent under Type 1 changes and attains $O_p(1)$ estimation errors under Type 2 changes. Additional results for other values of $(\rho,\alpha,\beta)$ are reported in Appendix D of the supplementary material.
Finally, we assess the precision of our information criterion ((ref)) in selecting the number of breaks under DGPs 1.A--E and 2 (with $m_{max}=5$). Table (ref) reports the percentage of correct detection of the number of breaks under DGP 1.A--E (all featuring two breaks) using the information criterion ((ref)). It can be seen that the proposed information criterion can approximately correctly identify the number of breakpoints as $N$ and $T$ increase.
DGP 2. In addition, DGP 2 considers scenarios where the true number of breaks, $m_{0}$, ranges from 0 to 4. The factors are generated by ((ref)) and we set $r_{0}=2$. For each regime, the $N\times2$ factor loading matrix has elements independently drawn from $N(0,1/r_{0})$. This implies that the number of pseudo-factors, $r$, is $(m_{0}+1)r_{0}$.
Before estimating the factors and factor loadings, we use Bai2002's (Bai2002) information criterion to determine the number of factors for the full samples with $r_{max}=12$. Notably, if no breakpoint exists, the criterion function in ((ref)) is simply zero, as the sample covariance of the factors becomes an identity matrix. As Table (ref) illustrates, the probability of our proposed information criterion accurately selecting the true number of breakpoints approaches one as the sample size increases. This finding further validates our Theorem (ref).
{6pt} {4pt} {6pt}
In this section we apply our method to detect breakpoints for the FRED-MD (Federal Reserve Economic Data - Monthly Data) data set. The FRED database is maintained by the research division of the Federal Reserve Bank of St. Louis, and is publicly accessible and updated in real time. Our analysis utilizes a sample of 126 unbalanced monthly time series, covering the period from January 1959 to July 2024. These series are categorized into eight groups: (1) output and income, (2) labor market, (3) housing, (4) consumption, orders, and inventories, (5) money and credit, (6) interest rates and exchange rates, (7) prices, and (8) stock market. Following the literature, we replace the missing values using the expectation-maximization algorithm and transform the data into approximately stationary series with the outliers replaced, so that we obtain a total of $T=785$ monthly observations for each macroeconomic variable. One can refer to McCracken2016 for the detailed data description.
In selecting the number of breaks, we set the minimum segment length to $h=20$ months and restrict candidate break dates to the monthly grid. We consider $m\in\{0,1,\ldots,8\}$ as the candidate numbers of breaks. For each $m$, we minimize the QML objective over all admissible partitions using the dynamic programming algorithm and then select $\hat m$ by minimizing Eq. ((ref)). In our empirical application, $N=126$, $r=7$ (selected by Bai--Ng IC2 on the full sample), and $\hat\rho=0.9574$, so the penalty term is numerically \[ \log(\min\{N,T\})\,r^2(1+|\hat\rho|)=\log(126)\times 7^2\times (1+0.9574)\approx 463.86. \] Hence the penalty used in the information criterion is $463.86\times m$. No additional tuning constant is used. The selected model yields five breakpoints: January 1969, January 1983, June 2008, March 2010, and February 2020. These breakpoints correspond to significant economic events. January 1969 marked a surge in U.S. inflation, primarily driven by increased spending on the Vietnam War and the “Great Society” programs. January 1983 roughly coincides with the onset of the Great Moderation, a period characterized by reduced macroeconomic volatility. June 2008 signifies the early stages of the Global Financial Crisis, as the collapse of the housing market began to impact the broader financial sector. March 2010 is approximately the end of the Great Recession for the U.S. economy. February 2020 represents the initial economic disruptions stemming from the COVID-19 pandemic, which rapidly led to a sharp global downturn.
Using the estimated breakpoints, we split the sample into six regimes and apply Bai--Ng IC2 within each regime, obtaining 2, 5, 7, 3, 4, and 7 factors, respectively. Under our classification rule, the variation in the estimated numbers of factors across adjacent regimes suggests that all identified breaks are Type 1 changes. For clarity, the regime with three factors is 2008:07--2010:03, whereas the following regime 2010:04--2020:02 contains four factors. Since principal components are identified only up to rotation, the economic meaning of individual factors is not always clear, and higher-order components are particularly difficult to interpret. We therefore focus on a more stable feature of the data: the concentration of explained variance in the leading principal components.
This concentration rises substantially in the crisis regimes. In 2008:07--2010:03, the first principal component explains about $34.7\%$ of total variance and the first three explain about $59.4\%$. In 2020:03--2024:07, the first principal component alone explains about $52.6\%$ of total variance. By contrast, in the non-crisis regimes 1983:02--2008:06 and 2010:04--2020:02, the first principal component explains only about $10.0\%$ and $13.2\%$, respectively. This indicates that crisis periods are dominated by a small number of strong common shocks, whereas common variation is more diffuse in normal times.
The two crisis episodes also differ in the composition of the leading common variation. During 2008:07--2010:03, the leading common variation is more closely associated with financial and price-related categories, including money and credit, interest rates and exchange rates, and prices. During 2020:03--2024:07, the leading common variation is much more strongly concentrated in output and income, as well as labor-market variables. This difference is consistent with the distinct nature of the Global Financial Crisis and the COVID-19 shock. Detailed regime-wise explained-variance results are reported in Appendix E of the supplementary material.
This paper introduces a quasi-maximum likelihood (QML) method for estimating multiple breakpoints in high-dimensional factor models. We analyze two distinct types of structural changes and develop corresponding asymptotic theories for the QML estimators. For Type 1 changes, we establish the consistency of the QML estimator, demonstrating that it converges to the true breakpoint with probability approaching one as the sample size increases. For Type 2 changes, we show that the distance between the QML estimators and the true breakpoints is stochastically bounded (fractionally consistent).
The proposed QML method is straightforward to implement and computationally efficient, requiring only a single full-sample eigendecomposition. Furthermore, we demonstrate that the number of breaks can be consistently estimated using our information criterion. Extensive simulation results confirm the strong finite-sample performance of our approach. Finally, an application to the FRED-MD dataset underscores the method's practical utility in empirical macroeconomic analysis.