The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
82,123 characters
\fontsize{12}{14pt plus.8pt minus .6pt}\selectfont \vspace{0.8pc}
\begin{center}
\large\bf
Quasi-maximum likelihood estimation of break point in high-dimensional factor models
\end{center}
\vspace{.4cm} \centerline{ Jiangtao Duan\textsuperscript{1}, Jushan Bai\textsuperscript{2}, Xu Han\textsuperscript{3} } \vspace{.4cm} \centerline{\it
\textsuperscript{1}Northeast Normal University, \textsuperscript{2}Columbia University and \textsuperscript{3}City University of Hong Kong}
\vspace{.55cm} \fontsize{9}{11.5pt plus.8pt minus
.6pt}\selectfont
\begin{quotation}
\noindent {\it Abstract:}
This paper estimates the break point for large-dimensional factor models with a single structural break in
factor loadings at a common unknown date. We propose a quasi-maximum likelihood (QML) estimator of
the change point based on the second moments of factors, which are estimated by a single principal component analysis.
We show that the QML estimator is consistent for the true break point when the covariance matrix of the
pre- or post-break factor loading (or both) is singular. Consistency here means that the deviation of the estimated break date
from the actual break date $k_0$ converges to zero as the sample size grows. This is a much stronger result than the break fraction
$\hat k/T$ being $T$-consistent (super-consistent) for $k_0/T$. Also, singularity occurs for most types of
structural changes, except for a rotational change. Even for a notational change, the QML estimator is still $T$-consistent
in terms of the break fraction. Simulation results confirm the theoretical properties of this estimator, and in fact QML significantly outperforms existing estimators for change points in factor models. Finally we apply the method to estimate the break points in a U.S. macroeconomic dataset and a stock return dataset.
\vspace{9pt}
\noindent {\it Key words and phrases:}
Structural break,
High-dimensional factor models,
Factor loadings\par
\end{quotation}\par
\vspace{-1em}
\section{Introduction}
Large factor models assume that a few factors can capture the common driving forces of a large number of economic variables.
Although factor models are useful, practitioners have to be cautious about the potential structural changes. For example, either the number of factors or the factor loadings may change over time. This
concern is empirically relevant because parameter instability is pervasive in large-scale panel data.
So far, many methods have been developed to test structural breaks in factor models (e.g., \cite{Stock2008}, \cite{Breitung2011}, and \cite{Chen2014}). The rejection of the null hypothesis of no
structural change leads to the subsequent issues of how to estimate the change point, determine the numbers of pre- and post-break factors, and estimate the factor space. \cite{Chen2015} considers a
least-squares estimator of the break point and proves the consistency of the estimated break fraction (i.e., the break date $k$ divided by the full time series $T$, $\frac{k}{T}$).
\cite{Cheng2016} propose a shrinkage method to obtain a consistent estimator of the break fraction.
\cite{Baltagi2017} develop a least-squares estimator of the change point based on the second moments of the estimated pseudo-factors and show that the estimation error of the proposed estimator is
$O_p(1)$, which indicates the consistency of the estimated break fraction. A few recent studies also explore a consistent estimation of break points, which is technically more challenging.
\cite{Ma_Su2018} develop an adaptive fused group Lasso method to consistently estimate all break points under a multibreak setup.
\cite{Barigozzi2018} propose a method based on wavelet transformations to consistently estimate the number and locations of break points in the common and idiosyncratic components.
\cite{Bai2017} establish the consistency of the least-squares estimator of the break point in large factor models when factor loadings are subjected to a structural break and the size of the break is
shrinking as the sample size increases. Although the estimators proposed in these studies are consistent under certain assumptions, the simulation results show that they perform poorly when (1) the
number of factors changes after the break or (2) the loading matrix undergoes a rotational type of change.
According to the factor model literature, a factor model with a break in factor loadings is observationally equivalent to that with constant loadings and possibly more pseudo-factors (e.g.,
\cite{Han2015} and \cite{Bai2016}). Thus, the estimation of the change point of factor loadings can be converted into that of the change point of the second moment of the pseudo-factors. We propose a
quasi-maximum likelihood (QML) method to estimate the break point based on the second moment of the estimated pseudo-factors; therefore, the number of original factors is not required to be known for
computing our estimator. First, we estimate the number of pseudo-factors (defined as the factors in the equivalent representation that ignores the break), and then estimate the pre- and
post-break second moment matrices of the estimated pseudo-factors for all possible sample splits. The structural break date is estimated by minimizing the QML function among all possible split points.
This paper makes the following contributions to the literature. First, we establish the consistency of the QML break point estimator if the break leads to more pseudo-factors than the original pre- or
post-break factors. This occurs when the break augments the factor space or in the presence of disappearing or emerging factors. Under these circumstances, the covariance matrix of loadings on the pre-
or post-break pseudo-factors is singular, which is the key condition to establish the consistency of our QML estimator. To the best of our knowledge, this is the first study that links the consistency
of the break point estimator to the singularity of covariance matrices of loadings on pre- and post-break pseudo-factors. In addition, we prove that the difference between the estimated and true change
points is stochastically bounded when both pre- and post-break loadings on the pseudo-factors have nonsingular covariance matrices. In this case, the loading matrix only undergoes a rotational change,
and both the numbers of pre- and post-break original factors are equal to the number of pseudo-factors.
The aforementioned singularity leads to a technical challenge of analyzing the asymptotic property. The singular population covariance matrix of the pre(post)-break loadings has a zero determinant,
whose logarithm is undefined. To resolve this issue, we show that the estimated covariance matrices have nonzero determinants and a well-defined inverse for any given sample size, by obtaining the
convergence rate of the lower bound of their smallest eigenvalues. This ensures that the objective function based on the estimated covariance is appropriately defined in any finite sample.
Our second major contribution is that the QML method allows a change in the number of factors. Namely, it allows for disappearing or emerging factors after the break. This is an advantage over the
methods developed by \cite{Ma_Su2018} and \cite{Bai2017}, who assume that the number of factors remains constant after the break. Our simulation result indicates that the estimator proposed by
\cite{Bai2017} is inconsistent when some factors disappear and the remaining factors have time-invariant loadings. \cite{Baltagi2017} allow a change in the number of factors; however, their estimation
error was only stochastically bounded. In contrast, our QML estimator remains consistent under a varying number of factors.
Finally, the QML method has a substantial computational advantage over the estimators that iteratively implement high-dimensional principal component analysis (PCA). For example, the estimator proposed
by \cite{Bai2017} runs PCA for pre- and post-split sample covariance matrices for all possible split points. In comparison, our QML runs PCA for the entire sample only once, and thus, is computationally
more efficient, especially in large samples.
The rest of this paper is organized as follows. Section 2 introduces the factor model with a single break on the factor loading matrix and describes the QML estimator for the break date. Section 3
presents the assumptions made for this model. Section 4 presents the consistency and asymptotic distribution of the QLM estimator for the break date.
Section 5 investigates the finite-sample properties of the QML estimator through simulations. Section 6 implements the proposed method to estimate the break points in a monthly macroeconomic dataset of
the United States and a dataset of weekly stock returns of Nasdaq 100 components. Section 7 concludes the study.
The following notations will be 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}$, its MP inverse as $\mathbb{A}^{-}$, its $i$-th singular value as $\sigma_{i}(\mathbb{A})$, and its adjoint matrix as $\mathbb{A}^{\#}$ when $m=n$. Let $\mathrm{Proj}(\mathbb{A}|\mathbb{Z})$ denote the projection of matrix $\mathbb{A}$ onto the columns of matrix $\mathbb{Z}$. For a real number $x$, $[x]$ represents the integer part of $x$.
\vspace{-1em}
\section{Model and estimator}
Let us consider the following factor model with a common break at $k_0$ in the factor loadings for $i=1,\cdots,N$:
\begin{eqnarray}
x_{it}=
\begin{cases}
\lambda_{i1}f_{t}+e_{it} & for~~ t=1,2,\cdots,k_0(T) \cr
\lambda_{i2}f_{t}+e_{it} & for~~ t=k_0(T)+1,\cdots,T,
\end{cases}
\label{model_1}
\end{eqnarray}
where $f_t$ is an $r-$dimensional vector of unobserved common factors; $r$ is the number of pseudo-factors; $k_0(T)$ is the unknown break date; $\lambda_{i1}$ and $\lambda_{i2}$ are the pre- and
post-break factor loadings, respectively; and $e_{it}$ is the error term allowed to have serial and cross-sectional dependence as well as heteroskedasticity.
$\tau_0\in (0,1)$ is a fixed constant and $[x]$ represents the integer part of $x$. For notational simplicity, hereinafter, we suppress the dependence of $k_0$ on $T$. Note that the
dimension of $f_t$ is the same as that of the pseudo-factors (to be defined soon) instead of the original underlying factors. This formulation simplifies the representation of various types of breaks in a unified framework, which
will be clarified in the examples below.
In vector form, model (\ref{model_1}) can be expressed as
\begin{eqnarray}
x_{t}=
\begin{cases}
\Lambda_{1}f_{t}+e_{t} & for~~ t=1,2,\cdots,k_0 \cr
\Lambda_{2}f_{t}+e_{t} & for~~ t=k_0+1,\cdots,T,
\end{cases}
\label{model_vector}
\end{eqnarray}
where $x_t=[ x_{1t},\cdots,x_{Nt} ]^{'}$, $e_t=[ e_{1t},\cdots,e_{Nt} ]^{'}$, $\Lambda_1=[ \lambda_{11},\cdots,\lambda_{N1} ]^{'}$, and $\Lambda_2=[ \lambda_{12},\cdots,\lambda_{N2} ]^{'}$.
For any $k=1,\cdots,T-1$, we define
$$X_{k}^{(1)}=[x_1,\cdots,x_k]^{'},X_{k}^{(2)}=[x_{k+1},\cdots,x_T]^{'},$$
$$F_{k}^{(1)}=[f_1,\cdots,f_k]^{'},F_k^{(2)}=[f_{k+1},\cdots,f_T],$$
$$\mbox{\boldmath $e$}_{k}^{(1)}=[e_1,\cdots,e_k]^{'},\mbox{\boldmath $e$}_k^{(2)}=[e_{k+1},\cdots,e_T],$$
where the subscript $k$ denotes the date at which the sample is to be split, and the superscripts $(1)$ and $(2)$ denote the pre- and post-$k$ data, respectively. We rewrite (\ref{model_vector}) using
the following matrix representation:
\begin{eqnarray}
\left[
\begin{array}{cccccccccc}
X_{k_0}^{(1)}\\
X_{k_0}^{(2)}\\
\end{array}
\right]
&=&\left[
\begin{matrix}
F_{k_0}^{(1)}\Lambda_1^{'}\\
F_{k_0}^{(2)}\Lambda_2^{'}
\end{matrix}
\right]+
\left[
\begin{array}{cccccccccc}
e_{k_0}^{(1)}\\
e_{k_0}^{(2)}\\
\end{array}
\right]
=\left[
\begin{matrix}
F_{k_0}^{(1)}(\Lambda B)^{'}\\
F_{k_0}^{(2)}(\Lambda C)^{'}
\end{matrix}
\right]+
\left[
\begin{array}{cccccccccc}
e_{k_0}^{(1)}\\
e_{k_0}^{(2)}\\
\end{array}
\right],\nonumber \\
&=&\left[
\begin{matrix}
F_{k_0}^{(1)}B^{'}\\
F_{k_0}^{(2)}C^{'}\\
\end{matrix}
\right]
\Lambda^{'}+
\left[
\begin{array}{cccccccccc}
e_{k_0}^{(1)}\\
e_{k_0}^{(2)}\\
\end{array}
\right], \nonumber\\
&=&G\Lambda^{'}+E,
\label{Baltagi}
\end{eqnarray}
where $F_{k_0}^{(1)}$ and $F_{k_0}^{(2)}$ have dimensions $k_0\times r$ and $(T-k_0)\times r$, respectively, and $\Lambda$ is an $N\times r$ matrix with full column rank. The pre- and post-break
loadings are modeled as $\Lambda_1=\Lambda B$ and $\Lambda_2=\Lambda C$, respectively, where $B$ and $C$ are some $r\times r$ matrices. Both $\Lambda_1$ and $\Lambda_2$ have dimension $N\times r$.
In this model, $r_1=rank(B)\leq r$ and $r_2=rank(C)\leq r$ denote the numbers of \emph{original factors} before and after the break, respectively. We refer to $G$ in (\ref{Baltagi}) as the \emph{pseudo-factors} because the last line of (\ref{Baltagi}) provides an observationally equivalent
representation without a change in the loadings matrix $\Lambda$. In other words, if the break is ignored in the estimation process, then the factors being estimated by a full-sample PCA are actually
the pseudo-factors $G$ in (\ref{Baltagi}).
It is well known that the break can augment the factor space; thus, $r_1 \leq r$ and $r_2 \leq r$, with $rank(G)=r$.
Our representation in (\ref{Baltagi}) allows for changes in the factor loadings
and the number of factors. Below, several examples are provided
to illustrate that the pseudo-factor representation in (\ref{Baltagi}) is general
enough to cover three types of breaks.
\textbf{Type 1}. Both $B$ and $C$ are singular. In this case, the
number of original factors is strictly less than that of the pseudo-factors
both before and after the break (i.e., $r_{1}<r$ and $r_{2}<r$).
This means that the structural break in the factor loadings augments
the dimension of the factor space. Let us consider the following example.
Example (1): Let $\mathbb{F}_{k_{0}}^{(1)}$$(k_{0}\times r_{1})$
and $\mathbb{F}_{k_{0}}^{(2)}$$((T-k_{0})\times r_{2})$ denote the
original factors before and after the break, respectively, and $\Theta_{1}$
and $\Theta_{2}$ denote the pre- and post-break loadings on these
factors. Thus, this model can be represented and transformed as \begin{eqnarray}
\left[\begin{array}{c}
X_{k_{0}}^{(1)}\\
X_{k_{0}}^{(2)}\end{array}\right] & = & \left[\begin{array}{c}
\mathbb{F}_{k_{0}}^{(1)}\Theta_{1}^{\prime}\\
\mathbb{F}_{k_{0}}^{(2)}\Theta_{2}^{\prime}\end{array}\right]+e=\left[\begin{array}{cc}
\mathbb{F}_{k_{0}}^{(1)} & 0\\
0 & \mathbb{F}_{k_{0}}^{(2)}\end{array}\right]\left[\begin{array}{c}
\Theta_{1}^{\prime}\\
\Theta_{2}^{\prime}\end{array}\right]+e\nonumber \\
& = & \left[\begin{array}{c}
[\mathbb{F}_{k_{0}}^{(1)}\;\vdots\;*]B^{\prime}\\
{}[*\;\vdots\;\mathbb{F}_{k_{0}}^{(2)}]C^{\prime}\end{array}\right]\Lambda^{\prime}+e=\underbrace{\left[\begin{array}{c}
F_{k_{0}}^{(1)}B^{\prime}\\
F_{k_{0}}^{(2)}C^{\prime}\end{array}\right]}_{G}\Lambda^{\prime}+e,\label{eq:ex1}\end{eqnarray}
where $\Lambda=[\Theta_{1},\Theta_{2}]$, $B=diag(I_{r_{1}},0_{r_{2}\times r_{2}})$,
$C=diag(0_{r_{1}\times r_{1}},I_{r_{2}})$, $F_{k_{0}}^{(1)}=[\mathbb{F}_{k_{0}}^{(1)}\;\vdots\;*]$,
$F_{k_{0}}^{(2)}=[*\;\vdots\;\mathbb{F}_{k_{0}}^{(2)}]$, and the asterisk
denotes some unidentified numbers such that all rows in $F_{k_{0}}^{(1)}$
and $F_{k_{0}}^{(2)}$ have the same variance (to satisfy Assumption
\ref{factors} in Section 3). (Note that the asterisk entries are cancelled due to multiplication by zero in $B$ and $C$.) In the special case of $r_{1}=r_{2}$, $\Lambda$
is of full rank $2r_{1}$ (i.e., the dimension of the pseudo-factor
space is twice that of the original factor space) if the shift in
the loading matrix $\Theta_{2}-\Theta_{1}$ is linearly independent
of $\Theta_{1}$. We refer to this special case as the shift type
of change, because the augmentation of the factor space is induced
by a linearly independent shift in the loading matrix. Hence, Type
1 covers the shift type of change.
\textbf{Type 2}. Only $B$ or $C$ is singular. In this case, emerging or disappearing factors are present in the model. Let us consider the following
example of disappearing factors.
Example (2): Without loss of generality, let us assume that $r_{2}<r_{1}$ and
$\Theta_{2}$ is equal to the first $r_{2}$ columns of $\Theta_{1}$;
thus, the last $r_{1}-r_{2}$ factors disappear after the break. Therefore,
we can obtain the pseudo-factors by using the following transformation
from the original factors $\mathbb{F}$: \begin{eqnarray}
\left[\begin{array}{c}
X_{k_{0}}^{(1)}\\
X_{k_{0}}^{(2)}\end{array}\right] & = & \left[\begin{array}{c}
\mathbb{F}_{k_{0}}^{(1)}\Theta_{1}^{\prime}\\
\mathbb{F}_{k_{0}}^{(2)}\Theta_{2}^{\prime}\end{array}\right]+e=\left[\begin{array}{c}
\mathbb{F}_{k_{0}}^{(1)}\Theta_{1}^{\prime}\\
{}[\mathbb{F}_{k_{0}}^{(2)}\;\vdots\;*]C^{\prime}\Theta_{1}^{\prime}\end{array}\right]+e\nonumber \\
& = & \left[\begin{array}{c}
F_{k_{0}}^{(1)}\\
F_{k_{0}}^{(2)}C^{\prime}\end{array}\right]\Theta_{1}^{\prime}+e=\underbrace{\left[\begin{array}{c}
F_{k_{0}}^{(1)}\\
F_{k_{0}}^{(2)}C^{\prime}\end{array}\right]}_{G}\Lambda^{\prime}+e,\label{eq:ex2}\end{eqnarray}
where $F_{k_{0}}^{(1)}=\mathbb{F}_{k_{0}}^{(1)}$, $F_{k_{0}}^{(2)}=[\mathbb{F}_{k_{0}}^{(2)}\;\vdots\;*]$,
$C=diag(I_{r_{2}},0_{(r_{1}-r_{2})\times(r_{1}-r_{2})})$, $\Lambda=\Theta_{1}$,
and the asterisk is defined in a similar manner to that in \eqref{eq:ex1}.
In this example, $B=I_{r_{1}}$, $r=r_{1}$, and $r_{2}=\mathrm{rank}(C)<r$.
Symmetrically, if $B$ is singular and $C=I_{r_{2}}$, then $r_{2}=r$
and $r_{1}=\mathrm{rank}(B)<r$, which means that certain factors emerge after the break point. Type 2 changes are important in empirical
analysis. Please refer to \cite{Mcalinn2018} for empirical evidence
regarding the varying number of factors in the U.S. macroeconomic dataset. For
Types 1 and 2, we obtain a significant result that $P(\hat{k}\ifmmode\begingroup\defbold{bold}
\text{\ifx\math@versionbold\bfseries\fi\textminus}\endgroup\else\textminus\fik_{0}=0)\to1$
as $N,T\to\infty$.
\footnote{Technically, Types 1 and 2 can be combined into one type that involves
singularity, which renders our QML estimator consistent. We consider
Type 2 separately to emphasize the case of emerging and disappearing
factors.
}
\textbf{Type 3}. Both $B$ and $C$ are nonsingular. In this case,
the loadings on the original factors undergo a rotational change,
and the dimension of the original factors is the same as that of the pseudo-factors.
Example (3): Let us assume that $r_{2}=r_{1}$ and $\Theta_{2}=\Theta_{1}C$
for a nonsingular matrix $C$. The model with the original factors
$\mathbb{F}$ can be transformed into the following pseudo-factor representation: \begin{eqnarray}
\left[\begin{array}{c}
X_{k_{0}}^{(1)}\\
X_{k_{0}}^{(2)}\end{array}\right] & = & \left[\begin{array}{c}
\mathbb{F}_{k_{0}}^{(1)}\Theta_{1}^{\prime}\\
\mathbb{F}_{k_{0}}^{(2)}\Theta_{2}^{\prime}\end{array}\right]+e=\left[\begin{array}{c}
\mathbb{F}_{k_{0}}^{(1)}\Theta_{1}^{\prime}\\
\mathbb{F}_{k_{0}}^{(2)}C^{\prime}\Theta_{1}^{\prime}\end{array}\right]+e\nonumber \\
& = & \left[\begin{array}{c}
F_{k_{0}}^{(1)}\\
F_{k_{0}}^{(2)}C^{\prime}\end{array}\right]\Theta_{1}^{\prime}+e=G\Lambda^{\prime}+e,\label{eq:ex3}\end{eqnarray}
where $F_{k_{0}}^{(1)}=\mathbb{F}_{k_{0}}^{(1)}$, $F_{k_{0}}^{(2)}=\mathbb{F}_{k_{0}}^{(2)}$,
and $\Lambda=\Theta_{1}$. In this example, $B=I_{r_{1}}$ and $r=r_{1}=r_{2}$,
and the factor dimension remains constant. In the observationally
equivalent pseudo-factor representation, the loading is time-invariant
and the original post-break factors $\mathbb{F}_{k_{0}}^{(2)}$ are
rotated by $C$. We refer to this as the rotation type of change.
The above examples show that a factor model with any of these three
types of change can be unified and reformulated by the representation
in (\ref{Baltagi}) with pseudo-factors. This representation controls the break
type by varying the settings for $B$ and $C$, and thus, is convenient
for our theoretical analysis.
\cite{Bai2017} rule out the rotation type of change because the break date is not identifiable by minimizing the sum of squared residuals.
\cite{Baltagi2017} allow changes in the number of factors and rotation type of change; however, the difference between their estimator and the true break point is only stochastically bounded (i.e.,
their estimator is not consistent).
Ma and Su's (2018) setup requires $r_1 = r_2$; thus, Type 2 is ruled out under their assumptions. Our simulation result shows that Ma and Su's estimator does not perform well under rotational changes
(Type 3), whereas our QML method can handle changes in all three types discussed above. We obtain a significant result that $\hat{k}-k_0=O_p(1)$ if both $B$ and $C$ are of full rank (i.e.,
Type 3) and $\hat{k}-k_0=o_p(1)$ if $B$ or $C$, or both, is singular (i.e., Type 1 and Type 2).
In this paper, we consider the QML estimator of the break date for model (\ref{Baltagi}):
\begin{eqnarray}
&&\hat{k}=\arg\min_{[\tau_1 T] \leq k \leq [\tau_2 T]} U_{NT}(k),
\label{obj_fun}
\end{eqnarray}
where $[\tau_1 T]$ and $[\tau_2 T]$ denote the prior lower and upper bounds for the real break point $k_0$ with $\tau_1,\tau_2\in (0,1)$ and $\tau_1 \leq \tau_0 \leq \tau_2$. The QML objective function
$U_{NT}(k)$ is equal to
\begin{eqnarray}
&&U_{NT}(k)=k\log( \det( \hat{\Sigma}_1 ) )+(T-k)\log( \det( \hat{\Sigma}_{2} ) ),
\label{obj}
\end{eqnarray}
where $\hat{\Sigma}_1$ and $\hat{\Sigma}_{2}$ are defined as
\begin{eqnarray}
\hat{\Sigma}_1&=&\frac{1}{k}\sum\limits_{t=1}^k \hat{g}_t \hat{g}_t^{'},\nonumber\\
\hat{\Sigma}_{2}&=&\frac{1}{T-k}\sum\limits_{t=k+1}^T \hat{g}_t \hat{g}_t^{'},
\label{Sig}
\end{eqnarray}
and $\hat{g}_t$ is the PCA estimator of $g_t$ (i.e., the transpose of the $t$-th row of $G$). We define $\Sigma_{G,1}=E(g_t g_t^{'})$ for $t\leq k_0$, $\Sigma_{G,2}=E(g_t g_t^{'})$ for $t>k_0$, and $\Sigma_G=\tau_0 \Sigma_{G,1}+(1-\tau_0)\Sigma_{G,2}$. We define $\Sigma_\Lambda$ as the covariance matrix of $\Lambda$. The PCA estimator $\hat{g}_t$ is asymptotically close to $H^{'}g_t$ for a rotation matrix $H$, and $H \xrightarrow{p} H_0=\Sigma_{\Lambda}^{1/2}\Phi V^{-1/2}$ as $(N,T)\rightarrow \infty$, where $V$ and
$\Phi$ are the eigenvalue and eigenvector matrices of $\Sigma_{\Lambda}^{1/2}\Sigma_G \Sigma_{\Lambda}^{1/2}$, respectively. Evidently, the second moment of $H_0 g_t$ shares the same change point as
that of $g_t$. Therefore, we proceed to estimate the pre- and post-break second moments of $g_t$ by using the estimated factors $\hat{g}_t$, and then use (\ref{obj_fun}) to obtain the QML break point
estimator $\hat{k}_{QML}$. Similar QML objective functions have been used for multivariate time series with observed data (e.g., \cite{Bai2000}).
\vspace{-1em}
\section{Assumptions}
\vspace{-1em}
In this section, we state the assumptions made for establishing the consistency and asymptotic distribution of the QML estimator.
\noindent \begin{assum}\label{factors}
(i) $E\left\|f_t \right\|^4<M<\infty$, $E(f_tf_t^{'})=\Sigma_F$, where $\Sigma_F$ is positive definite, and
$\frac{1}{k_0}\sum_{t=1}^{k_0}f_tf_t^{'}\xrightarrow{p}\Sigma_F,\frac{1}{T-k_0}\sum_{t=k_0+1}^{T}f_tf_t^{'}\xrightarrow{p}\Sigma_F$;
(ii) There exists $d>0$ such that $\left\|\Delta\right\|\geq d>0$, where $\Delta =B\Sigma_FB^{'}- C\Sigma_FC^{'}$ and $B,C$ are $r\times r$ matrices.
\end{assum}
\noindent \begin{assum}\label{Factor_Loadings}
$\left\| \lambda_{\ell i} \right\|\leq \bar{\lambda}<\infty$ for $\ell=1,2$, $i=1,\cdots,N$, $\left\| \frac{1}{N}\Lambda^{'}\Lambda-\Sigma_{\Lambda} \right\|\rightarrow 0$ for some $r\times r$ positive
definite matrix $\Sigma_{\Lambda}$.
\end{assum}
\noindent \begin{assum}\label{Depen_and_Hetero}
There exists a positive constant $M<\infty$ such that
\begin{itemize}
\item[(i)] $E(e_{it})=0$ and $E|e_{it}|^8\leq M$ for all $i=1,\cdots,N$ and $t=1,\cdots,T$;
\item[(ii)] $E(\frac{e_s^{'}e_t}{N})=E(N^{-1}\sum_{i=1}^Ne_{is}e_{it})=\gamma_N(s,t)$ and $\sum_{s=1}^{T}|\gamma_N(s,t)|\leq M$ for every $t\leq T$;
\item[(iii)] $E(e_{it}e_{jt})=\tau_{ij,t}$ with $|\tau_{ij,t} |<\tau_{ij}$ for some $\tau_{ij}$ and for all $t=1,\cdots,T$ and $\sum_{j=1}^{N}|\tau_{ij}|\leq M$ for every $i\leq N$;
\item[(iv)] $E(e_{it}e_{js})=\tau_{ij,ts}$,
\begin{equation*} \frac{1}{NT}\sum\limits_{i,j,t,s=1}|\tau_{ij,ts}|\leq M;
\end{equation*}
\item[(v)] For every $(s,t)$, $E\left| N^{-1/2}\sum_{i=1}^{N}(e_{is}e_{it}-E[e_{is}e_{it}]) \right|^4\leq M$.
\end{itemize}
\end{assum}
\noindent \begin{assum}\label{Weak_Dependence}
There exists a positive constant $M<\infty$ such that
\begin{eqnarray*}
E(\frac{1}{N}\sum\limits_{i=1}^{N}\left\| \frac{1}{\sqrt{k_0}}\sum\limits_{t=1}^{k_0}f_te_{it} \right\|^2)&\leq& M,\\
E(\frac{1}{N}\sum\limits_{i=1}^{N}\left\| \frac{1}{\sqrt{T-k_0}}\sum\limits_{t=k_0+1}^{T}f_te_{it}\right\|^2)&\leq& M.
\end{eqnarray*}
\end{assum}
\noindent \begin{assum}\label{eigenvalues}
The eigenvalues of $\Sigma_G\Sigma_\Lambda$ are distinct.
\end{assum}
\noindent \begin{assum}\label{Hajek-Renyi}
Let us define $\epsilon_t=f_tf_t^{'}-\Sigma_F$. According to the data-generating process (DGP) of factors, the H\'{a}jek-R\'{e}nyi inequality applies to the processes $\{\epsilon_t,t=1,\cdots,k_0\}$,
$\{\epsilon_t,t=k_0,\cdots,1\}$, $\{\epsilon_t,t=k_0+1,\cdots,T\}$, and $\{\epsilon_t,t=T,\cdots,k_0+1\}$.
\end{assum}
\begin{remark}
Using the H\'{a}jek-R\'{e}nyi equality on $\epsilon_t$, we can ensure that $\max\limits_{k_0<k\leq [\tau_2 T]}\| \frac{1}{T-k} \sum\limits_{t=k+1}^{T} f_tf_t^{'}-\Sigma_F\|=O_p(\frac{1}{\sqrt{T}})$ in
Lemma \ref{differ} and $\max\limits_{[\tau_1 T]\leq k<k_0}\| \frac{1}{k_0-k} \sum\limits_{t=k+1}^{k_0} g_tg_t^{'} \|=O_p(1), \max\limits_{k_0<k\leq [\tau_2 T]}\| \frac{1}{k-k_0}
\sum\limits_{t=k_0+1}^{k} g_tg_t^{'} \|=O_p(1)$ in Lemmas \ref{differ} and \ref{differ2}.
\label{Hajek_Renyi_Generalized inequality}
\end{remark}
\noindent \begin{assum}\label{error_sup}
There exists an $M<\infty$ such that
(i) For each $s=1,\cdots,T$,
\begin{eqnarray*}
E(\max_{k<k_0}\frac{1}{k_0-k}\sum_{t=k+1}^{k_0}|\frac{1}{\sqrt{N}}\sum_{i=1}^{N}[e_{is}e_{it}-E(e_{is}e_{it})]|^2)&\leq& M,\\
E(\max_{k> k_0}\frac{1}{k-k_0}\sum_{t=k_0+1}^{k}|\frac{1}{\sqrt{N}}\sum_{i=1}^{N}[e_{is}e_{it}-E(e_{is}e_{it})]|^2)&\leq& M;
\end{eqnarray*}
(ii)
\begin{eqnarray*}
E(\max_{k<k_0}\frac{1}{k_0-k}\sum_{t=k+1}^{k_0}\left\|\frac{1}{\sqrt{N}}\sum_{i=1}^{N} \lambda_i e_{it} \right\|^2)&\leq &M,\\
E(\max_{k>k_0}\frac{1}{k_0-k}\sum_{t=k_0+1}^{k}\left\|\frac{1}{\sqrt{N}}\sum_{i=1}^{N} \lambda_i e_{it} \right\|^2)&\leq &M.
\end{eqnarray*}
\end{assum}
\noindent \begin{assum}\label{Central_Limit}
There exists an $M<\infty$ such that for all values of $N$ and $T$,
(i) for each $t$,
\begin{eqnarray*}
E\left(\max_{1 \leq k<k_0}\frac{1}{k_0-k}\sum_{t=k+1}^{k_0} \left\| \frac{1}{\sqrt{NT}}\sum_{s=1}^{T}\sum_{i=1}^Nf_s[e_{is}e_{it}-E(e_{is}e_{it})] \right\|^2\right)&\leq& M,\\
E\left(\max_{k_0<k\leq T}\frac{1}{k-k_0}\sum_{t=k_0+1}^{k} \left\| \frac{1}{\sqrt{NT}}\sum_{s=1}^{T}\sum_{i=1}^Nf_s[e_{is}e_{it}-E(e_{is}e_{it})] \right\|^2\right)&\leq& M;
\end{eqnarray*}
(ii) the $r\times r$ matrix satisfies
\begin{eqnarray*}
E\left\| \frac{1}{\sqrt{NT}}\sum_{t=1}^{T}\sum_{i=1}^Nf_t\lambda_i^{'}e_{it} \right\|^2\leq M.
\end{eqnarray*}
\end{assum}
\vspace{-1em}
\section{Asymptotic properties of the QML estimator}
\vspace{-1em}
In this section, we derive the asymptotic properties of the QML estimator for various breaks. In the literature of structural breaks for a fixed-dimensional time series, conventional break point estimators, such as the least-squares (LS) estimator of \cite{Bai1997} or the QML estimator of
\cite{Qu_Perron2007}, are usually inconsistent. The estimation error of these conventional estimators is $O_p(1)$ when the break size is fixed. To reach consistency, the cross-sectional dimension of the
time series must be large (e.g., \cite{Bai2010} and \cite{Kim2011}).
Recall that the observationally equivalent representation in (\ref{Baltagi}) has time-invariant loadings and varying pseudo-factors. Hence, our problem converges to estimating the break point in the
$r$-dimensional time series $g_t$, where $r$ is fixed. Theorems \ref{bound_theorem} and \ref{distribution_theorem} below show that, for rotational breaks (Type 3), the convergence rate and limiting
distribution are similar to those available in the literature. However, for Type 1 and 2 breaks, Theorem \ref{consistency} derives a much more significant result than that available in the literature,
according to which our QML estimator is consistent even if our $g_t$ has only a fixed cross-sectional dimension $r$.
\noindent
\begin{theorem}
Under Assumptions \ref{factors}--\ref{Central_Limit}, when both $B$ and $C$ are of full rank, $\hat{k}-k_0=O_p(1)$.
\label{bound_theorem}
\end{theorem}
This theorem implies that the difference between the QML estimator and the true change point is stochastically bounded in model (\ref{eq:ex3}). Although the estimation errors of both \cite{Baltagi2017} and our QML
methods are bounded, the QML estimator has much better finite sample properties. To confirm this theoretical result, we conduct a simulation where the factor loadings have a rotational change (see
DGP 1.B in Section 5). Table \ref{rotation_full_rank} presents the MAEs and RMSEs of different estimators. The simulation result shows that the QML estimators have much smaller MAEs and RMSEs than other
methods. In addition, $\hat{k}$ does not collapse to $k_0$, leading to a nondegenerate distribution. We will state the limiting distribution in Theorem \ref{distribution_theorem}. Nevertheless, this
theorem shows that the break point can be appropriately estimated because $\hat{\tau}=\hat{k}/T$ is still consistent for $\tau_0$.
To make an inference regarding the change point when both $B$ and $C$ are of full rank, we derive the limiting distribution of $\hat{k}$. Let us define
\begin{eqnarray*}
\xi_t & = & H_{0}^{'}g_{t}g_{t}^{'}H_{0}-\Sigma_{1}\text{ for }t\leq k_{0},\\
\xi_t & = & H_{0}^{'}g_{t}g_{t}^{'}H_{0}-\Sigma_{2}\text{ for }t>k_{0},
\end{eqnarray*}
where $\Sigma_{1}=H_{0}^{'}\Sigma_{G,1}H_{0}$ and $\Sigma_{2}=H_{0}^{'}\Sigma_{G,2}H_{0}$
are the pre- and post-breaks of $H_{0}^{'}E(g_{t}g_{t}^{'})H_{0}$. The limiting distribution of $\hat{k}$ is given by the following theorem:
\noindent
\begin{theorem}
Under Assumptions \ref{factors}--\ref{Central_Limit}, when both $B$ and $C$ are of full rank,
\begin{eqnarray*}
\hat{k}-k_0\xrightarrow{d} \arg\min\limits_\ell W(\ell),
\end{eqnarray*}
where
\begin{eqnarray*}
&&W(\ell)=\sum\limits_{t=k_0+\ell}^{k_0-1}tr((\Sigma_2^{-1}-\Sigma_1^{-1})\xi_t)-\left( tr(\Sigma_1\Sigma_2^{-1})-r-\log|\Sigma_1\Sigma_2^{-1}| \right)\ell\\
&&\text{for }\ell=-1,-2,\cdots,\\
&&W(\ell)=0\text{ for }\ell=0,\\
&&W(\ell)=\sum\limits_{t=k_0+1}^{k_0+\ell}tr((\Sigma_1^{-1}-\Sigma_2^{-1})\xi_t)+\left( tr(\Sigma_1^{-1}\Sigma_2)-r-\log|\Sigma_1^{-1}\Sigma_2| \right)\ell\\
&&\text{for }\ell=1,2,\cdots.
\end{eqnarray*}
\label{distribution_theorem}
\end{theorem}
This result shows that the limiting distribution depends on $\xi_t$.
If $\xi_t$ is independent over time, then $W(\ell)$ is a two-sided random walk. If $f_t$ is stationary, then $\xi_t$ is stationary in each regime.
Here, the limiting distribution of the estimated break date is dependent on the generation processes of the unobserved factors, and thus, cannot be directly used to construct a confidence interval for a
true break point.
\cite{Bai2017} propose a bootstrap method to construct a confidence interval for $k_0$ when the change in the factor loading matrix shrinks as $N\to \infty$.
However, their bootstrap procedure lacks robustness in the cross-sectional correlation in the error terms.
In the current setup, the break magnitude $\left\| \Sigma_2-\Sigma_1 \right\|$ is fixed and we leave the case of shrinking break magnitude as a future topic.
Next, we establish a much stronger result than that available in the literature, which states that the QML estimator remains consistent when $B$ or $C$, or both, is singular. We make the following
additional assumptions.
\noindent \begin{assum}\label{as-.LeeL}
With probability approaching one (w.p.a.1), the following inequalities hold:
\begin{align*}
&0<\underline{c}\le \min_{[\tau_{1}T]\le k\le k_{0}}\rho_{j}\left(\frac{1}{Nk}\sum_{t=1}^{k}\Lambda^{\prime}e_{t}e_{t}^{\prime}\Lambda\right),\\
&0<\underline{c}\le \min_{k_{0}\le k\le[\tau_{2}T]}\rho_{j}\left(\frac{1}{N(T-k)}\sum_{t=k+1}^{T}\Lambda^{\prime}e_{t}e_{t}^{\prime}\Lambda\right),\,\text{for }j=1,\cdots,r;\\
&\rho_{1}\left(\frac{1}{NT}\sum_{t=1}^{T}\Lambda^{\prime}e_{t}e_{t}^{\prime}\Lambda\right)\le \overline{c}<+\infty,
\end{align*}
as $N,T\to\infty$, where $\underline{c}$ and $\overline{c}$ are some constants.
\noindent \end{assum}
\noindent \begin{assum}\label{invariance}
\begin{align*}
\max_{[\tau_{1}T]\le k\le k_{0}}\left\Vert \frac{1}{\sqrt{Nk}}\sum_{t=1}^{k}\sum_{i=1}^{N}f_{t}e_{it}\lambda_{i}^{\prime}\right\Vert & =O_{p}(1),\\
\max_{k_{0}\le k\le[\tau_{2}T]}\left\Vert \frac{1}{\sqrt{N(T-k)}}\sum_{t=k+1}^{T}\sum_{i=1}^{N}f_{t}e_{it}\lambda_{i}^{\prime}\right\Vert & =O_{p}(1).\end{align*}
\end{assum}\vspace{2em}
Assumption \ref{as-.LeeL} is useful to derive the lower bound of the smallest eigenvalue of $\hat{\Sigma}_1$ (or $\hat{\Sigma}_2$) if $B$ (or $C$) is a singular matrix. Assumption \ref{invariance} strengthens Assumption \ref{Central_Limit}(ii), which is similar to Assumption F2 of \cite{Bai2003}. Note that the summation $\sum\limits_{t=1}^k$ in Assumptions \ref{as-.LeeL}-\ref{invariance} involves a positive fraction of observations over time since the lower bound of $k$ is $\tau_1 T$, with $\tau_1\in (0,1)$.
Also, as the log of matrix determinant is involved in the QML function, a natural problem is that the
log determinant of a singular population covariance matrix is undefined when $B$ or $C$, or both, is singular. Fortunately, the determinants of $\hat{\Sigma}_1=\frac{1}{k}\sum\limits_{t=1}^k \hat{g}_t \hat{g}_t^{'}$ and
$\hat{\Sigma}_2=\frac{1}{T-k}\sum\limits_{t=k+1}^T \hat{g}_t \hat{g}_t^{'}$ are small but not equal to zero in finite samples, when $\Sigma_1$ and $\Sigma_2$ are
singular matrices. The following proposition develops a lower bound for the smallest
eigenvalues of $\hat{\Sigma}_1$ and $\hat{\Sigma}_2$.
\begin{proposition}\label{low_bound} Under Assumptions \ref{factors}--\ref{invariance}, for $k\ge k_{0}$
and $k\le[\tau_{2}T]$, if $C$ is singular and $\sqrt{N}/T\to0$
as $N,T\to\infty$, then there exist constants
$c_{U}\ge c_{L}>0$ such that \begin{align*}
P\left(\min_{k\in[k_{0},[\tau_{2}T]]}\rho_{j}(\hat{\Sigma}_{2})\ge\frac{c_{L}}{N}\right) & \to1,\\
P\left(\max_{k\in[k_{0},[\tau_{2}T]]}\rho_{j}(\hat{\Sigma}_{2})\le\frac{c_{U}}{N}\right) & \to1,\end{align*}
for $j=r_{2}+1,...,r$.
\end{proposition}
In proposition \ref{low_bound}, the lower bound of the smallest eigenvalue of the estimated sample covariance matrix $\hat{\Sigma}_2$ is $c_L/N$ for a constant $c_L>0$ w.p.a.1. A similar lower bound for the smallest eigenvalue of $\hat{\Sigma}_1$ can be obtained when $B$ is singular under the same assumptions.
This ensures a lower bound for the determinants of the estimated sample covariance matrices. Proposition \ref{low_bound} provides a useful tool to establish the consistency of our QML estimator.
Although this technical result is a byproduct in our analysis, we believe that it is of independent interest and useful in other contexts.
\noindent \begin{assum}\label{B_C_full_rank_project}
(i) $[B,C]$ is of full row rank.
(ii) $C^{\#}Bf_{k_0}\neq 0$ when $r-1=r_2>0$; and $B^{\#}Cf_{k_0+1}\neq 0$ when $r-1=r_1>0$, where $\mathbb{A}^{\#}$ denotes the adjoint matrix for a singular matrix $\mathbb{A}$.
(iii) $\|Bf_{k_{0}}-\mathrm{Proj}(Bf_{k_{0}}|C)\|\ge d>0$ when $r-r_2\geq 2$ or $r_{2}=0$; and $\|Cf_{k_{0}+1}-\mathrm{Proj}(Cf_{k_{0}+1}|B)\|\ge d>0$ when $r-r_1\geq 2$ or $r_{1}=0$, where
$\mathrm{Proj}(\mathbb{A}|\mathbb{Z})$ denotes the projection
of $\mathbb{A}$ onto the columns of $\mathbb{Z}$, and $d$ is a constant.
\end{assum}
Assumption \ref{B_C_full_rank_project}(i) implies that $\Sigma_G$ is positive definite.$\footnote{Since \begin{eqnarray*}
rank(\Sigma_G)&=&rank\left(\left[
\begin{array}{cccccccccc}
\sqrt{\tau_0}B,\sqrt{1-\tau_0}C
\end{array}
\right]
diag\left(\Sigma_F,\Sigma_F\right)
\left[
\begin{array}{cccccccccc}
\sqrt{\tau_0}B,\sqrt{1-\tau_0}C
\end{array}
\right]^{'}\right)
=rank\left(\left[
\begin{array}{cccccccccc}
\sqrt{\tau_0}B,\sqrt{1-\tau_0}C
\end{array}
\right]\right)\\
&=&rank\left(
[B,C]diag\left( \sqrt{\tau_0}I_r, \sqrt{1-\tau_0}I_r \right) \right)=rank\left([B,C]\right)
\end{eqnarray*}
and $1<\tau_0<1$, Assumption \ref{B_C_full_rank_project} implies that $\Sigma_G$ is a positive definite matrix.
}
$ Assumptions \ref{B_C_full_rank_project}(ii) implies that $B^{\#}C\neq 0$ when $r-1=r_1>0$, and $C^{\#}B\neq 0$ when $r-1=r_2>0$. It also excludes the possibility that $f_{k_0}$ and $f_{k_0+1}$ are in the null space of $C^{\#}B$ and $B^{\#}C$, respectively. Similarly, Assumption \ref{B_C_full_rank_project}(iii) rules out the cases that $Bf_{k_0}$ lies in the column space of $C$ when $r-r_2\geq 2$ or $r_2=0$ and that $Cf_{k_0+1}$ lies in the column space of $B$ when $r-r_1\geq 2$ or $r_1=0$.\footnote{Note that $r_2=0$ means $C=0$, so $B$ has to be nonsingular by Assumption \ref{B_C_full_rank_project}(i). Thus, Assumption \ref{B_C_full_rank_project}(iii) implies that $f_{k_0}\neq 0$ when $C=0$.}
Assumption \ref{B_C_full_rank_project} is used to establish Lemma \ref{differ2}, which is useful for validating the consistency result that $Prob(\hat{k}-k = 0) \to 1$ in the proof of Theorem
\ref{consistency}. It ensures that the value of the objective function becomes larger even if $\hat{k}$ slightly deviates from the true break point in large samples. Assumption \ref{B_C_full_rank_project} is flexible enough to allow various data generating processes for $f_t$. For example, if $f_{k_{0}}$ and $f_{k_{0}+1}$ have continuous probability distribution functions, then Assumptions \ref{B_C_full_rank_project}(ii)-(iii) just exclude a zero probability event since $C^{\#}B$ and $B^{\#}C$ are not equal to zero.
It is remarkable that existing estimators such as \cite{Baltagi2017} and \cite{Bai2017} are not consistent even if Assumptions \ref{as-.LeeL} -- \ref{B_C_full_rank_project} hold. In contrast, our QML estimator is shown to be consistent under these additional assumptions. The following theorem summarizes the result.
\begin{theorem}
Under Assumptions \ref{factors}--\ref{B_C_full_rank_project} and $\frac{N}{T}\to\kappa$, as
$N,T\to\infty$ for $0<\kappa<\infty$, when $B$ or $C$, or both, is singular, $Prob(\hat{k}-k = 0) \to 1.$
\label{consistency}
\end{theorem}
Theorem \ref{consistency} shows that the estimated change point converges to the true change point w.p.a.1 when $B$ or $C$, or both, is singular (Types 1 and 2 in Section 2). This result is much more
significant than that obtained by \cite{Baltagi2017}, who show that the distance between the estimated and true break dates is bounded for Types 1--3. Note that the case in which only $B$ (or $C$) is
singular corresponds to Type 2 with emerging (or disappearing) factors. Our QML estimator is consistent under this type of change, whereas \cite{Bai2017} and \cite{Ma_Su2018} rule out this type by
assumption. In empirical applications, the conditions of theorem \ref{consistency} are rather flexible and likely to hold and the consistency of the break date estimator is expected
in most economic data for the factor analysis.
\begin{remark}
An important contribution of Theorem \ref{consistency} is to link the consistency of the QML estimator with the singularity of the covariance matrices of the pre- or post-break factor loadings.
The singularity is generated by the special structure of the pseudo-factors $g_t$ shown in (\ref{eq:ex1}) and (\ref{eq:ex2}) in the presence of a structural change.
The PCA estimator $\hat{g}_t$ is consistent for $g_t$ (up to some rotation) for large $N$ and $T$, so the singularity structure is maintained in $\hat{\Sigma}_1$ and $\hat{\Sigma}_2$ and hence contributes to the consistency of our QML estimator. The result in Theorem \ref{consistency} is in contrast to
conventional break point estimators, which only have $O_p(1)$ estimation errors in multivariate time series with a small cross-sectional
dimension (e.g., \cite{Bai1997}; \cite{Qu_Perron2007}). Although our $\hat{g}_t$ has a fixed dimension, the divergence rate of the objective function depends on $N$. $\footnote{This is because the convergence rate of the smallest eigenvalue of $\hat{\Sigma}_2$ is $N^{-1}$ for $k = k_0$ when $C$ is singular. See Proposition 1.}$
In other words, our QML estimator still implicitly utilizes the information in the large cross-sectional dimension, which is the source of our consistency.
\label{consistency_remark}
\end{remark}
\begin{remark}
The conditions that $B$ or $C$, or both, is singular and $\frac{N}{T}\rightarrow \kappa\in (0,\infty)$ are likely to hold in many economic datasets for factor
analysis.
If both $B$ and $C$ are singular, the break occurs such that the number of pseudo-factors in the entire factor model is larger than that of the factors in the pre- and post-break subsamples. This can happen when the factor loadings undergo a shift type of change, as discussed in Example (1) for Type 1 changes.
If $B$ is of full rank and $C$ is singular, some factors become irrelevant, and thus, the loading coefficients attached to these disappearing factors become zero. For example, in the momentum portfolio,
some risks are not part of the firm's long-run structure as only sorting based on recent returns works; the reward is high but disappears within less than a year.
If $B$ is singular and $C$ is of full rank, some factors emerge after the break date, increasing the dimension of the post-break factor space. For example, changes in the technology or policy may
produce certain new factors.
\label{consistency_remark2}
\end{remark}
\begin{remark}
Theorem \ref{consistency} indicates that $U_{NT}(k)$ can be minimized to consistently estimate $k_0$. The intuition for this is that $U_{NT}(k)-U_{NT}(k_0)$ is always larger than zero, even if $k$
deviates only slightly from the true break point $k_0$, so that $\hat{k}$ must be equal to $k_0$ to minimize $U_{NT}(k)-U_{NT}(k_0)$.
For example, in Type 1, when both $B$ and $C$ are singular for $k<k_0$, we can decompose $\hat{\Sigma}_2$ as $\hat{\Sigma}_2=\frac{1}{T-k}\sum\limits_{t=k+1}^{k_0} \hat{g}_t
\hat{g}_t^{'}+\frac{1}{T-k}\sum\limits_{t=k_0+1}^T \hat{g}_t \hat{g}_t^{'}$, and the term $\frac{1}{T-k}\sum\limits_{t=k+1}^{k_0} \hat{g}_t \hat{g}_t^{'}$ results in a larger determinant of $\hat{\Sigma}_2$ than that of $\hat{\Sigma}_2^{0}=\frac{1}{T-k_0}\sum\limits_{t=k_0+1}^T \hat{g}_t \hat{g}_t^{'}$.
By symmetry, we obtain a similar result for $k>k_0$. (See Lemmas \ref{differ} and \ref{differ2} for more technical details.) Thus, $U_{NT}(k)-U_{NT}(k_0)>0$ w.p.a.1 as $N,T\to\infty$ if $k\neq k_0$.
\label{singular_remark}
\end{remark}
\begin{remark}
With the QML estimator, we do not need to know the numbers of original factors $r_1$ and $r_2$ before and after the break point, but only the number of pseudo-factors in the entire sample.
\cite{Bai2017} and \cite{Ma_Su2018} require knowledge of the number of original factors, which is much more difficult to estimate due to the augmented factor space resulting from the break. In practice,
the number of pseudo-factors is much easier to estimate by using one of a number of estimators, such as the information criteria developed by \cite{Bai2002}.
\label{pseudo_true_fators}
\end{remark}
\vspace{-1em}
\section{Simulation}
\vspace{-1em}
In this section, we consider DGPs corresponding to Types 1--3 to evaluate the finite sample performance of the QML estimator. We compare the QML estimator with three other estimators. As shown below,
$\hat{k}_{BKW}$ is the estimator proposed by Baltagi, Kao, and Wang (2017, BKW hereafter); $\hat{k}_{BHS}$ is the estimator proposed by Bai, Han, and Shi (2020, BHS hereafter); $\hat{k}_{MS}$ is the
estimator proposed by Ma and Su (2018, MS hereafter); and $\hat{k}_{QML}$ is the QML estimator. \cite{Barigozzi2018} develops a change point estimator using wavelet transformation, which exhibits
similar performance to that of the estimator proposed by \cite{Ma_Su2018}. Hence, the comparison with the estimator proposed by \cite{Barigozzi2018} is not reported here, but the result is available
upon request.
The DGP roughly follows BKW, which can be used to examine various elements that may affect the finite sample performance of the estimators, and we use this DGP for model (\ref{Baltagi}).
We calculate the root mean square error (RMSE) and mean absolute error (MAE) of these change point estimators $\hat{k}_{BKW}$, $\hat{k}_{BHS}$, and $\hat{k}_{QML}$, and each experiment is repeated 1000
times, where RMSE$=\sqrt{\frac{1}{1000}\sum\limits_{s=1}^{1000}(\hat{k}_s-k_0)^2}$ and MAE$=\frac{1}{1000}\sum\limits_{s=1}^{1000}|\hat{k}_s-k_0|$.
When $T$ is small, there is a possibility that Ma and Su's (2018) method detects no break or multiple breaks; thus, the definition of the estimation error for a single break point in such cases is not
straightforward. For a comparison, we compute the RMSE and MAE of the MS estimator by only using the results obtained by the MS estimator when it successfully detects a single break.
As the computation of $\hat{k}_{BHS}$ and $\hat{k}_{MS}$ requires the number of original factors and that of $\hat{k}_{BKW}$ and $\hat{k}_{QML}$ requires the number of pseudo-factors, we set
$\hat{r}=r_0$ for $\hat{k}_{BHS}$ and $\hat{k}_{MS}$ and $\hat{r}=r$ for $\hat{k}_{QML}$ and $\hat{k}_{BKW}$, where $r_0$ is the number of original factors and $r$ is the number of pseudo-factors.
We generate factors and idiosyncratic errors using a DGP similar to that of BKW. Each factor is generated by the following AR(1) process:
\begin{eqnarray*}
f_{tp}=\rho f_{t-1,p}+u_{t,p},\quad for\quad t=2,\cdots,T;\quad p=1,\cdots,r_0,
\end{eqnarray*}
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 i.i.d. $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
\begin{eqnarray*}
e_{i,t}=\alpha e_{i,t-1}+v_{i,t},\quad for\quad i=1,\cdots,N\quad t=2,\cdots,T,
\end{eqnarray*}
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 set $r_0=3$ and $k_0=T/2$.
We consider the following DGPs for factor loadings and investigate the performance of the QML estimator for the three types of breaks discussed in Section 2.
\textbf{DGP 1.A} We first consider the case in which $C$ is singular, and set $C=[1,0,0;0,1,0;0,0,0]$. This setup aims to model (\ref{eq:ex2}). In the pre-break regime, all elements of $\lambda_{i,1}$
are i.i.d. $N(0,\frac{1}{r_0^2}I_{r_0})$ across $i$. In the post-break regime, $\Lambda_2=(\lambda_{1,2},\cdots,\lambda_{N,2})^{'}=\Lambda_1C$. This case corresponds to a Type 2 change with a
disappearing factor. The number of pseudo-factors is the same as $r_0$, so $r=3$, and the numbers of pre- and post-break factors are 3 and rank$(C)=2$, respectively.
Table \ref{rotation_singular} lists the RMSEs and MAEs of three estimators for different values of $(\rho, \alpha, \beta)$.
In all cases, $\hat{k}_{QML}$ has much smaller MAEs and RMSEs than $\hat{k}_{BKW}$ and $\hat{k}_{BHS}$.
Moreover, the MAEs and RMSEs of $\hat{k}_{QML}$ tend to decrease as $N$ and $T$ increase. This confirms the consistency of $\hat{k}_{QML}$ established in Theorem \ref{consistency}. In addition, the
RMSEs and MAEs of $\hat{k}_{BKW}$ do not converge to zero as $N$ and $T$ increase, which confirms that $\hat{k}_{BKW}$ has a stochastically bounded estimation error. $\hat{k}_{BHS}$ does not appear to
be consistent when a factor disappears after the break. Moreover, a larger AR(1) coefficient $\rho$ tends to deteriorate the performance of $\hat{k}_{BKW}$, but does not have much impact on our QML
estimator.
\textbf{DGP 1.B} We next consider the case in which $C$ is of full rank. We set $C$ as a lower triangular matrix. The diagonal elements are equal to $0.5$, $1.5$, and $2.5$, and the elements below these
diagonal elements are i.i.d. and drawn from a standard normal distribution. Under this DGP, we have $r = r_0$.
Table \ref{rotation_full_rank} reports the performance of three estimators for different values of $(\rho, \alpha, \beta)$.
In all cases, $\hat{k}_{BKW}$ and $\hat{k}_{QML}$ appear to have stochastically bounded estimation errors, which confirms Theorem 1 of BKW and Theorem \ref{bound_theorem} of this paper.
Both $\hat{k}_{QML}$ and $\hat{k}_{BKW}$ are inconsistent under this DGP; however, under all settings, our QML estimator tends to have much smaller RMSEs and MAEs than the estimator of BKW.
The MAEs and RMSEs of $\hat{k}_{BHS}$ appear to increase with the sample size; thus, the BHS method cannot handle this case.
\textbf{DGP 1.C} In this case, we set $C=[1,0,0;2,1,0;3,2,m]$ and $m\in \{1,0.8,0.5,0.1,0\}$. As $m$ decreases to zero, the matrix $C$ changes from full rank to singular. We still consider serial
correlation in factors and serial correlation and cross-sectional dependence in idiosyncratic errors simultaneously with $N=100, T=100$. Table \ref{full_rank_to_singular} shows that the MAEs and RMSEs
of $\hat{k}_{QML}$ monotonically decrease with $m$, which confirms our findings in Theorems \ref{bound_theorem} and \ref{consistency}.
In addition, the RMSEs and MAEs of $\hat{k}_{BKW}$ and $\hat{k}_{BHS}$ are much larger than those of $\hat{k}_{QML}$, and do not tend toward zero as $m$ decreases. For each value of $m$, the experiment
is repeated 10000 times to more accurately estimate and compare the RMSEs (MAEs) of our QML estimator across different values of $m$.
\textbf{DGP 1.D} This DGP considers a Type 1 break. In the first regime, the last elements of $\lambda_{i,1}$ are zeros for all $i$, and the first two elements of $\lambda_{i,1}$ are both i.i.d.
$N(0,\frac{1}{2}I_{r_0})$. In the second regime, $\lambda_{i,2}$ is i.i.d. $N(0,\frac{1}{3}I_{r_0})$ across $i$. As $\lambda_{i,1}$ and $\lambda_{i,2}$ are independent, the numbers of factors in the two
regimes are $r_1=2$ and $r_2=3$, respectively, and the number of pseudo-factors is $r=5$. Because the numbers of pre- or post-break factors are smaller than that of the pseudo-factors, both $\Sigma_1$ and
$\Sigma_2$ are singular matrices.
Table \ref{singular_pre_post_break} reports the MAEs and RMSEs of $\hat{k}_{QML}$, $\hat{k}_{BHS}$, and $\hat{k}_{BKW}$ under this DGP.
Table \ref{singular_pre_post_break} shows the performances of both $\hat{k}_{BHS}$ and our $\hat{k}_{QML}$. Their MAEs (RMSEs) are less than 0.05 (0.25) for all combinations of $N$, $T$,
$\rho$, $\alpha$, and $\beta$. Although $\hat{k}_{BHS}$ is consistent under this DGP, our QML estimator still has smaller RMSEs than $\hat{k}_{BHS}$ in most cases reported in Table
\ref{singular_pre_post_break}. In addition, $\hat{k}_{BKW}$ performs better under this DGP than DGPs 1.A--1.C. However, its estimation error is much larger than that of our QML estimator. This is not
surprising because $\hat{k}_{BKW}$ is not consistent. Finally, a larger AR(1) coefficient $\rho$ tends to yield a larger bias for $\hat{k}_{BKW}$, but does not have much effect on the performances of
$\hat{k}_{BHS}$ and $\hat{k}_{QML}$.
In summary, Tables \ref{rotation_singular} and \ref{rotation_full_rank} show that the QML estimator performs much better than $\hat{k}_{BHS}$ under Type 2 and 3 breaks, which are ruled out under the
assumptions of \cite{Bai2017}. Table \ref{singular_pre_post_break} shows that the QML estimator often slightly outperforms $\hat{k}_{BHS}$, even though the latter is known to be consistent and has excellent finite-sample performance under Type
1 breaks. Note that the strength of the BHS method is the consistent estimation of break point for Type 1 break, especially when the size of the break is shrinking as the sample size increases. Under settings with a shrinking break size, the QML method will lose its power because the dimension
of $G$ (determined by the IC criterion in \cite{Bai2002}) will not be augmented, which means that the singularity does not show up in the covariance if breaks are small enough.
\begin{table}[H]
\caption{Simulated mean absolute errors (MAEs) and root mean squared errors (RMSEs) of $\hat{k}_{BKW}$, $\hat{k}_{BHS}$, and $\hat{k}_{QML}$ under DGP 1.A.}
\centering
\label{rotation_singular}
\begin{center}
\begin{tabular}{l l r r r r r r} \hline
$N,T$ & & \multicolumn{2}{c}{$\hat{k}_{BKW}$} & \multicolumn{2}{c}{$\hat{k}_{BHS}$} & \multicolumn{2}{c}{$\hat{k}_{QML}$} \\
& & MAE & RMSE & MAE & RMSE & MAE & RMSE \\ \hline
& & & &$\rho=0$ &$\alpha=0$ &$\beta=0$ & \\\hline
100,100 & &6.3130 &8.9546 &5.4600 &7.7325 &1.6070 &2.9293 \\
100,200 & &7.0230 &11.9053 &7.9580 &12.4801 &1.2990 &2.3206 \\
200,200 & &5.6730 &9.9774 &6.7150 &10.8610 &0.7960 &1.5218 \\
200,500 & &4.6940 &8.5732 &10.0960 &17.9778 &0.7340 &1.3799 \\
500,500 & &4.4580 &8.5789 &8.6770 &15.6509 &0.3890 &0.8597 \\\hline
& & & &$\rho=0.7$ &$\alpha=0$ &$\beta=0$ & \\\hline
100,100 & &9.7200 &12.0612 &4.5670 &6.9270 &1.3570 &2.7592 \\
100,200 & &14.3410 &19.5941 &7.0110 &11.1559 &1.0470 &2.2070 \\
200,200 & &13.6260 &19.1151 &6.7760 &10.9099 &0.5840 &1.2394 \\
200,500 & &15.4880 &27.5716 &10.5450 &18.7350 &0.5190 &1.1406 \\
500,500 & &16.9890 &29.5463 &8.2030 &15.1581 &0.3210 &0.7944 \\\hline
& & & &$\rho=0$ &$\alpha=0.3$ &$\beta=0$ & \\\hline
100,100 & &6.5060 &9.1533 &6.1520 &8.6248 &2.3740 &4.0635 \\
100,200 & &7.5490 &12.4416 &8.7150 &13.4473 &1.6920 &3.1464 \\
200,200 & &6.2890 &10.8337 &8.4910 &13.2894 &1.0230 &1.9409 \\
200,500 & &5.1220 &10.1068 &11.3960 &19.4945 &0.8110 &1.5156 \\
500,500 & &4.7580 &9.5055 &10.3660 &18.7453 &0.4570 &0.9407 \\\hline
& & & &$\rho=0$ &$\alpha=0$ &$\beta=0.3$ & \\\hline
100,100 & &6.6620 &9.2573 &4.7300 &6.9593 &1.7580 &3.1183 \\
100,200 & &7.8200 &12.5561 &6.1740 &10.2069 &1.4930 &2.6943 \\
200,200 & &6.4500 &10.9881 &5.8020 &9.7340 &0.7480 &1.4276 \\
200,500 & &4.9340 &10.3110 &5.9390 &10.6041 &0.7020 &1.3900 \\
500,500 & &4.0550 &7.7718 &5.8820 &11.0830 &0.3660 &0.8567 \\\hline
& & & &$\rho=0.7$ &$\alpha=0.3$ &$\beta=0.3$ & \\\hline
100,100 & &9.9510 &12.3063 &5.4430 &7.6969 &1.8080 &3.4531 \\
100,200 & &14.2890 &19.5804 &7.1810 &11.7141 &1.3250 &2.5367 \\
200,200 & &14.8820 &20.3572 &7.3080 &11.8072 &0.7450 &1.6592 \\
200,500 & &17.0330 &29.3210 &9.2010 &17.0803 &0.6680 &1.3461 \\
500,500 & &14.7130 &26.3587 &10.3800 &19.2727 &0.3540 &0.8331 \\\hline
\end{tabular}
\end{center}
\end{table}
\begin{table}[H]
\caption{ Simulated mean absolute errors (MAEs) and root mean squared errors (RMSEs) of $\hat{k}_{BKW}$, $\hat{k}_{BHS}$, and $\hat{k}_{QML}$ under DGP 1.B.}
\centering
\label{rotation_full_rank}
\begin{center}
\begin{tabular}{l l r r r r r r} \hline
$N,T$ & & \multicolumn{2}{c}{$\hat{k}_{BKW}$} & \multicolumn{2}{c}{$\hat{k}_{BHS}$} & \multicolumn{2}{c}{$\hat{k}_{QML}$} \\
& & MAE & RMSE & MAE & RMSE & MAE & RMSE \\ \hline
& & & &$\rho=0$, &$\alpha=0$, &$\beta=0$ & \\\hline
100,100 & &4.1610 &6.6934 &8.7430 &11.0347 &1.2180 &2.3259 \\
100,200 & &4.4450 &8.4477 &18.5660 &22.9913 &0.9960 &1.8799 \\
200,200 & &4.9160 &8.9420 &19.4440 &23.6923 &0.9060 &1.7082 \\
200,500 & &4.4530 &8.8368 &49.3330 &59.4865 &0.9130 &1.7085 \\
500,500 & &3.9420 &7.2061 &51.9270 &61.5507 &0.8370 &1.5959 \\\hline
& & & &$\rho=0.7$ &$\alpha=0$ &$\beta=0$ & \\\hline
100,100 & &6.4570 &9.4427 &10.1710 &12.3371 &1.9460 &3.7691 \\
100,200 & &9.1750 &14.8115 &21.3380 &25.1834 &1.8480 &3.6362 \\
200,200 & &9.6310 &15.0080 &21.5560 &25.2723 &1.7850 &3.5901 \\
200,500 & &11.4150 &21.3302 &51.9850 &61.9028 &1.6750 &3.4218 \\
500,500 & &9.5430 &18.4598 &53.6060 &62.7128 &1.6490 &3.5501 \\\hline
& & & &$\rho=0$ &$\alpha=0.3$ &$\beta=0$ & \\\hline
100,100 & &3.9840 &6.4778 &7.9990 &10.5485 &1.0910 &2.1824 \\
100,200 & &4.6820 &8.6151 &17.6010 &22.5002 &1.0360 &1.9432 \\
200,200 & &4.6350 &8.4454 &21.9190 &26.0996 &0.8770 &1.7306 \\
200,500 & &4.2690 &8.2870 &50.1790 &61.5307 &0.8600 &1.6474 \\
500,500 & &4.2040 &8.3094 &54.8050 &64.8615 &0.8040 &1.5492 \\\hline
& & & &$\rho=0$ &$\alpha=0$ &$\beta=0.3$ & \\\hline
100,100 & &4.3220 &6.9244 &7.7560 &10.0601 &1.0510 &1.9802 \\
100,200 & &4.7150 &8.6248 &14.5640 &19.4154 &0.9830 &1.8571 \\
200,200 & &4.5300 &8.2421 &18.7950 &23.1307 &0.9090 &1.7587 \\
200,500 & &3.9080 &7.3553 &42.6850 &54.8098 &0.8900 &1.6199 \\
500,500 & &4.3570 &8.5140 &49.6030 &59.6292 &0.8250 &1.6843 \\\hline
& & & &$\rho=0.7$ &$\alpha=0.3$ &$\beta=0.3$ & \\\hline
100,100 & &6.6990 &9.6327 &9.1750 &11.4037 &2.0750 &3.9735 \\
100,200 & &9.4990 &15.0852 &18.7450 &23.3085 &2.0590 &4.5305 \\
200,200 & &9.2240 &14.6721 &20.4670 &24.5054 &1.8140 &3.8021 \\
200,500 & &12.8110 &23.1517 &51.0760 &61.1890 &1.7200 &3.5844 \\
500,500 & &10.0590 &19.2453 &52.5400 &62.2628 &1.7000 &3.6521 \\\hline
\end{tabular}
\end{center}
\end{table}
\begin{table}[H]
\caption{ Simulated mean absolute errors (MAEs) and root mean squared errors (RMSEs) of $\hat{k}_{BKW}$, $\hat{k}_{BHS}$, and $\hat{k}_{QML}$ under DGP 1.C with $N=100,T=100$ among 10000 replications.}
\centering
\label{full_rank_to_singular}
\begin{center}
\begin{tabular}{l l r r r r r r} \hline
$m$ & & \multicolumn{2}{c}{$\hat{k}_{BKW}$} & \multicolumn{2}{c}{$\hat{k}_{BHS}$} & \multicolumn{2}{c}{$\hat{k}_{QML}$} \\
& & MAE & RMSE & MAE & RMSE & MAE & RMSE \\ \hline
& & & &$\rho=0$ &$\alpha=0$ &$\beta=0$ & \\\hline
1 & &3.9228 &6.4437 &7.2579 &9.7141 &0.6562 &1.2903 \\
0.8 & &3.9425 &6.4624 &6.6330 &9.1145 &0.6348 &1.2559 \\
0.5 & &3.7847 &6.2319 &5.4950 &7.9789 &0.5420 &1.0814 \\
0.1 & &3.8469 &6.2895 &4.6050 &6.9212 &0.5093 &1.0568 \\
0 & &3.8310 &6.2414 &4.4915 &6.8352 &0.4969 &1.0315 \\\hline
& & & &$\rho=0.7$ &$\alpha=0$ &$\beta=0$ & \\\hline
1 & &6.0404 &9.0280 &9.3131 &11.5733 &0.9478 &2.0680 \\
0.8 & &6.0063 &9.0017 &8.5168 &10.9192 &0.8547 &1.8960 \\
0.5 & &5.9803 &8.9390 &6.5127 &9.0641 &0.6925 &1.5752 \\
0.1 & &5.9300 &8.8833 &4.6894 &7.0610 &0.5178 &1.2335 \\
0 & &6.0197 &8.9771 &4.5440 &6.9049 &0.5070 &1.2057 \\\hline
& & & &$\rho=0$ &$\alpha=0.3$ &$\beta=0$ & \\\hline
1 & &3.8349 &6.2423 &7.1824 &9.7359 &0.6727 &1.3234 \\
0.8 & &3.8234 &6.2331 &6.6551 &9.2338 &0.6535 &1.2963 \\
0.5 & &3.8345 &6.3110 &5.8040 &8.3371 &0.6152 &1.2362 \\
0.1 & &3.9127 &6.4083 &5.0645 &7.4846 &0.5895 &1.1644 \\
0 & &3.9188 &6.4124 &4.9815 &7.3974 &0.5813 &1.1551 \\\hline
& & & &$\rho=0$ &$\alpha=0$ &$\beta=0.3$ & \\\hline
1 & &3.8250 &6.3150 &6.2535 &8.7224 &0.6622 &1.3039 \\
0.8 & &3.8135 &6.2932 &5.6808 &8.1438 &0.6259 &1.2379 \\
0.5 & &3.8253 &6.3061 &4.6189 &6.9328 &0.5619 &1.1171 \\
0.1 & &3.9120 &6.4147 &3.9299 &6.0820 &0.5424 &1.0949 \\
0 & &3.8176 &6.2881 &3.8963 &6.0564 &0.5199 &1.0497 \\\hline
& & & &$\rho=0.7$ &$\alpha=0.3$ &$\beta=0.3$ & \\\hline
1 & &6.0745 &9.0347 &8.0648 &10.5304 &1.0515 &2.2669 \\
0.8 & &6.0041 &8.9542 &7.3126 &9.8433 &0.9338 &2.0173 \\
0.5 & &6.0519 &9.0124 &5.8471 &8.4490 &0.7798 &1.7537 \\
0.1 & &6.0120 &8.9694 &4.6376 &7.1447 &0.6100 &1.4401 \\
0 & &6.0379 &8.9861 &4.5336 &7.0337 &0.5850 &1.3509 \\\hline
\end{tabular}
\end{center}
\end{table}
\begin{table}[H]
\caption{Simulated mean absolute errors (MAEs) and root mean squared errors (RMSEs) of $\hat{k}_{BKW}$, $\hat{k}_{BHS}$, and $\hat{k}_{QML}$ under DGP 1.D.}
\centering
\label{singular_pre_post_break}
\begin{center}
\begin{tabular}{l l r r r r r r} \hline
$N,T$ & & \multicolumn{2}{c}{$\hat{k}_{BKW}$} & \multicolumn{2}{c}{$\hat{k}_{BHS}$} & \multicolumn{2}{c}{$\hat{k}_{QML}$} \\
& & MAE & RMSE & MAE & RMSE & MAE & RMSE \\ \hline
& & & &$\rho=0$, &$\alpha=0$, &$\beta=0$ & \\\hline
100,100 & &0.4330 &1.3494 &0.0370 &0.1975 &0.0260 &0.1673 \\
100,200 & &0.3380 &1.0900 &0.0300 &0.1732 &0.0240 &0.1549 \\
200,200 & &0.2780 &0.7668 &0.0180 &0.1342 &0.0130 &0.1140 \\
200,500 & &0.2850 &0.8155 &0.0070 &0.0837 &0.0100 &0.1000 \\\hline
& & & &$\rho=0.7$ &$\alpha=0$ &$\beta=0$ & \\\hline
100,100 & &1.8760 &4.8750 &0.0120 &0.1095 &0.0110 &0.1049 \\
100,200 & &1.1140 &4.0007 &0.0150 &0.1225 &0.0110 &0.1140 \\
200,200 & &0.8700 &3.5000 &0.0050 &0.0707 &0.0020 &0.0447 \\
200,500 & &0.4070 &1.3435 &0.0030 &0.0548 &0.0010 &0.0316 \\\hline
& & & &$\rho=0$ &$\alpha=0.3$ &$\beta=0$ & \\\hline
100,100 & &0.4400 &1.4519 &0.0450 &0.2302 &0.0410 &0.2258 \\
100,200 & &0.3590 &1.3802 &0.0440 &0.2145 &0.0340 &0.1897 \\
200,200 & &0.3080 &0.8438 &0.0150 &0.1225 &0.0140 &0.1265 \\
200,500 & &0.2150 &0.6656 &0.0160 &0.1265 &0.0120 &0.1095 \\\hline
& & & &$\rho=0$ &$\alpha=0$ &$\beta=0.3$ & \\\hline
100,100 & &0.3710 &1.0747 &0.0380 &0.1949 &0.0360 &0.1897 \\
100,200 & &0.2850 &0.7918 &0.0340 &0.1897 &0.0220 &0.1483 \\
200,200 & &0.3150 &0.8972 &0.0100 &0.1000 &0.0110 &0.1049 \\
200,500 & &0.2380 &0.6885 &0.0120 &0.1183 &0.0050 &0.0707 \\\hline
& & & &$\rho=0.7$ &$\alpha=0.3$ &$\beta=0.3$ & \\\hline
100,100 & &1.9420 &4.8557 &0.0260 &0.1612 &0.0180 &0.1414 \\
100,200 & &1.0170 &3.6438 &0.0220 &0.1549 &0.0090 &0.0949 \\
200,200 & &0.9750 &3.9242 &0.0060 &0.0775 &0.0080 &0.0894 \\
200,500 & &0.6390 &2.5879 &0.0080 &0.0894 &0.0050 &0.0707 \\\hline
\end{tabular}
\end{center}
\end{table}
Tables \ref{rotation_singular_correct_Pro}--\ref{singular_pre_post_correct_Pro} present the probabilities of the correct estimation of the break date.
The results are consistent with those displayed in Tables \ref{rotation_singular}--\ref{singular_pre_post_break}: the QML estimator $\hat{k}_{QML}$ can detect the true break date with higher
probabilities than others regardless of the values of $(\rho,\alpha,\beta)$.
The MS method sometimes detects more than one or no break; hence, we only compute its probability of correctly estimating $k_0$ under the condition that it detects a single break. The probabilities of a
correct estimation of the QML method increase with the sample sizes $N$ and $T$ in Tables \ref{rotation_singular_correct_Pro}, \ref{rotation_full_rank_correct_Pro}, and
\ref{singular_pre_post_correct_Pro}.
Table \ref{full_rank_to_singular_correct_Pro} shows that the probabilities of correct estimation of the QML estimators increase as $m$ decreases. A smaller $m$ means that $C$ is closer to a singular
matrix. Table \ref{full_rank_to_singular_correct_Pro} is consistent with Table \ref{full_rank_to_singular}, and confirms Theorems \ref{bound_theorem} and \ref{consistency}. To explore in more detail the
effect of changes in $m$ on the QML estimator, we vary the value of $m$ using finer grids and find a similar pattern to that shown in Table \ref{full_rank_to_singular_correct_Pro}. The results are
reported in the supplementary appendix.
Figures \ref{1A_NT100} and \ref{1A_NT500} show the frequency of the estimated change points under DGP 1.A for $N=100,T=100$ and $N=500,T=500$ for 1000 replications. According to these figures, the QML
estimators exhibit the highest frequency around the true break under different settings. When we increase the $(N,T)$ value from $100$ to $500$, the frequency at the true break point increases and the
simulated distribution becomes tighter. This indicates that the QML estimators are highly likely to identify the true break point. This is consistent with our theory. However, the other three methods
are found to have much larger variation and substantially lower probabilities to correctly estimate the break point. Thus, the QML estimators are advantageous in this case.
Moreover, the simulation result indicates that for a sample size exceeding $N = 5000, T = 1000$, the probabilities of correctly estimating the QML estimator exceed $90\%$.
Recall that BKW and QML only have $O_p(1)$ estimation errors under DGP 1.B. However, Table \ref{rotation_full_rank_correct_Pro} shows that in all cases, the probabilities of correct estimation by the
QML estimator are much higher than those of correct estimation by the BKW estimator
Apparently, the BHS and MS methods cannot accurately estimate the true break point in this case.
Figures \ref{1B_NT100} and \ref{1B_NT500} show the distributions of the estimated change points under (1.B) for $N=100,T=100$ and $N=500,T=500$,
indicating that BHS and MS cannot handle rotational changes.
Although the estimation errors of BKW and QML are bounded under all settings, the QML estimators have a much tighter distribution around the true break point.
\begin{table}[H]
\caption{Probability of correct estimation under DGP 1.A.}
\centering
\label{rotation_singular_correct_Pro}
\begin{center}
\begin{tabular}{l l r r r r r r} \hline
$N,T$ & & \multicolumn{1}{c}{$\hat{k}_{BKW}$} & \multicolumn{1}{c}{$\hat{k}_{BHS}$} & \multicolumn{1}{c}{$\hat{k}_{MS}$} & \multicolumn{1}{c}{$\hat{k}_{QML}$} \\
& & & & & & \\ \hline
& &$\rho=0$ &$\alpha=0$ &$\beta=0$ & \\\hline
100,100& &0.1530 &0.1440 &0.1626 &0.4220 \\
100,200& &0.1920 &0.1510 &0.1863 &0.4370 \\
200,200& &0.2340 &0.1780 &0.1307 &0.5680 \\
200,500& &0.2540 &0.2030 &0.2020 &0.5780 \\
500,500& &0.2990 &0.2100 &0.2123 &0.7290 \\\hline
& &$\rho=0.7$ &$\alpha=0$ &$\beta=0$ & \\\hline
100,100& &0.1050 &0.2050 &0.2329 &0.5290 \\
100,200& &0.1250 &0.1850 &0.1779 &0.5510 \\
200,200& &0.1390 &0.1920 &0.1898 &0.6660 \\
200,500& &0.1750 &0.1890 &0.2031 &0.6940 \\
500,500& &0.2100 &0.2420 &0.2306 &0.7810 \\\hline
& &$\rho=0$ &$\alpha=0.3$ &$\beta=0$ & \\\hline
100,100& &0.1790 &0.1300 &0.1072 &0.3280 \\
100,200& &0.1850 &0.1380 &0.1897 &0.4090 \\
200,200& &0.2260 &0.1650 &0.1931 &0.5320 \\
200,500& &0.2530 &0.1730 &0.1845 &0.5650 \\
500,500& &0.2750 &0.1920 &0.1964 &0.6880 \\\hline
& &$\rho=0$ &$\alpha=0$ &$\beta=0.3$ & \\\hline
100,100& &0.1480 &0.1700 &0.1956 &0.3840 \\
100,200& &0.1730 &0.1810 &0.1847 &0.4210 \\
200,200& &0.2240 &0.2110 &0.2069 &0.5700 \\
200,500& &0.2770 &0.2250 &0.2370 &0.5930 \\
500,500& &0.3220 &0.2790 &0.2790 &0.7500 \\\hline
& &$\rho=0.7$ &$\alpha=0.3$ &$\beta=0.3$ & \\\hline
100,100& &0.1070 &0.1510 &0.1739 &0.4670 \\
100,200& &0.1210 &0.1860 &0.2157 &0.5030 \\
200,200& &0.1370 &0.1820 &0.2072 &0.6360 \\
200,500& &0.1670 &0.2180 &0.2149 &0.6520 \\
500,500& &0.1900 &0.2510 &0.2427 &0.7640 \\\hline
\end{tabular}
\end{center}
\end{table}
\begin{table}[H]
\caption{Probability of correct estimation under DGP 1.B.}
\centering
\label{rotation_full_rank_correct_Pro}
\begin{center}
\begin{tabular}{l l r r r r r r} \hline
$N,T$ & & \multicolumn{1}{c}{$\hat{k}_{BKW}$} & \multicolumn{1}{c}{$\hat{k}_{BHS}$} & \multicolumn{1}{c}{$\hat{k}_{MS}$} & \multicolumn{1}{c}{$\hat{k}_{QML}$} \\
& & & & & & \\ \hline
& &$\rho=0$ &$\alpha=0$ &$\beta=0$ & \\\hline
100,100& &0.2760 &0.0690 &0.0769 &0.4790 \\
100,200& &0.2920 &0.0540 &0.0362 &0.5180 \\
200,200& &0.2720 &0.0320 &0.0655 &0.5270 \\
200,500& &0.3110 &0.0140 &0.0091 &0.5340 \\
500,500& &0.2960 &0.0100 &0.0123 &0.5580 \\\hline
& &$\rho=0.7$ &$\alpha=0$ &$\beta=0$ & \\\hline
100,100& &0.2710 &0.0640 &0.0909 &0.4540 \\
100,200& &0.2500 &0.0270 &0.0398 &0.4530 \\
200,200& &0.2180 &0.0160 &0.0200 &0.4790 \\
200,500& &0.2370 &0.0120 &0.0144 &0.4970 \\
500,500& &0.2450 &0.0080 &0.0080 &0.5090 \\\hline
& &$\rho=0$ &$\alpha=0.3$ &$\beta=0$ & \\\hline
100,100& &0.3050 &0.1060 &0.1163 &0.5180 \\
100,200& &0.2930 &0.0740 &0.0989 &0.5020 \\
200,200& &0.2890 &0.0390 &0.0496 &0.5540 \\
200,500& &0.3000 &0.0230 &0.0328 &0.5630 \\
500,500& &0.3090 &0.0090 &0.0125 &0.5780 \\\hline
& &$\rho=0$ &$\alpha=0$ &$\beta=0.3$ & \\\hline
100,100& &0.2740 &0.0880 &0.1458 &0.5000 \\
100,200& &0.2970 &0.0650 &0.0692 &0.5220 \\
200,200& &0.2870 &0.0390 &0.0338 &0.5390 \\
200,500& &0.3100 &0.0300 &0.0320 &0.5290 \\
500,500& &0.2940 &0.0120 &0.0123 &0.5810 \\\hline
& &$\rho=0.7$ &$\alpha=0.3$ &$\beta=0.3$ & \\\hline
100,100& &0.2210 &0.1000 &0.1524 &0.4330 \\
100,200& &0.2400 &0.0610 &0.0763 &0.4640 \\
200,200& &0.2370 &0.0490 &0.0538 &0.4810 \\
200,500& &0.2230 &0.0230 &0.0218 &0.4770 \\
500,500& &0.2420 &0.0160 &0.0207 &0.5100 \\\hline
\end{tabular}
\end{center}
\end{table}
\begin{table}[H]
\caption{Probability of correct estimation under DGP 1.C with $N=100,T=100$.}
\centering
\label{full_rank_to_singular_correct_Pro}
\begin{center}
\begin{tabular}{l l r r r r r r} \hline
$m$ & & \multicolumn{1}{c}{$\hat{k}_{BKW}$} & \multicolumn{1}{c}{$\hat{k}_{BHS}$} & \multicolumn{1}{c}{$\hat{k}_{MS}$} & \multicolumn{1}{c}{$\hat{k}_{QML}$} \\
& & & & & & \\ \hline
& &$\rho=0$ &$\alpha=0$ &$\beta=0$ & \\\hline
1 & &0.3044 &0.1075 &0.1188 &0.6079 \\
0.8& &0.3033 &0.1252 &0.1389 &0.6153 \\
0.5& &0.3009 &0.1736 &0.1904 &0.6467 \\
0.1& &0.2976 &0.1998 &0.2014 &0.6680 \\
0 & &0.2977 &0.2031 &0.2192 &0.6705 \\\hline
& &$\rho=0.7$ &$\alpha=0$ &$\beta=0$ & \\\hline
1 & &0.2841 &0.0781 &0.1040 &0.6051 \\
0.8& &0.2871 &0.0975 &0.1135 &0.6254 \\
0.5& &0.2896 &0.1532 &0.1620 &0.6641 \\
0.1& &0.2876 &0.2073 &0.2369 &0.7131 \\
0 & &0.2848 &0.2194 &0.2297 &0.7154 \\\hline
& &$\rho=0$ &$\alpha=0.3$ &$\beta=0$ & \\\hline
1 & &0.2973 &0.1226 &0.1442 &0.6063 \\
0.8& &0.2981 &0.1416 &0.1648 &0.6134 \\
0.5& &0.2988 &0.1641 &0.1730 &0.6219 \\
0.1& &0.2993 &0.1828 &0.1954 &0.6316 \\
0 & &0.2988 &0.1860 &0.1995 &0.6342 \\\hline
& &$\rho=0$ &$\alpha=0$ &$\beta=0.3$ & \\\hline
1 & &0.3018 &0.1344 &0.1399 &0.6075 \\
0.8& &0.3016 &0.1549 &0.1679 &0.6164 \\
0.5& &0.3044 &0.1927 &0.2078 &0.6383 \\
0.1& &0.3009 &0.2141 &0.2093 &0.6461 \\
0 & &0.3036 &0.2211 &0.2314 &0.6566 \\\hline
& &$\rho=0.7$ &$\alpha=0.3$ &$\beta=0.3$ & \\\hline
1 & &0.2821 &0.1519 &0.1739 &0.5921 \\
0.8& &0.2844 &0.1710 &0.1964 &0.6082 \\
0.5& &0.2843 &0.2327 &0.2402 &0.6496 \\
0.1& &0.2850 &0.2889 &0.2911 &0.6898 \\
0 & &0.2868 &0.2951 &0.2966 &0.6951 \\\hline
\end{tabular}
\end{center}
\end{table}
\begin{table}[H]
\caption{Probability of correct estimation under DGP 1.D.}
\centering
\label{singular_pre_post_correct_Pro}
\begin{center}
\begin{tabular}{l l r r r r r r} \hline
$N,T$ & & \multicolumn{1}{c}{$\hat{k}_{BKW}$} & \multicolumn{1}{c}{$\hat{k}_{BHS}$} & \multicolumn{1}{c}{$\hat{k}_{MS}$} & \multicolumn{1}{c}{$\hat{k}_{QML}$} \\
& & & & & & \\ \hline
& &$\rho=0$ &$\alpha=0$ &$\beta=0$ & \\\hline
100,100& &0.7960 &0.9640 &0.9553 &0.9750 \\
100,200& &0.8200 &0.9700 &0.9700 &0.9760 \\
200,200& &0.8160 &0.9820 &0.9841 &0.9870 \\
200,500& &0.8260 &0.9930 &0.9930 &0.9900 \\\hline
& &$\rho=0.7$ &$\alpha=0$ &$\beta=0$ & \\\hline
100,100& &0.7000 &0.9880 &0.9864 &0.9890 \\
100,200& &0.7540 &0.9850 &0.9859 &0.9900 \\
200,200& &0.7950 &0.9950 &0.9949 &0.9980 \\
200,500& &0.8220 &0.9970 &0.9970 &0.9990 \\\hline
& &$\rho=0$ &$\alpha=0.3$ &$\beta=0$ & \\\hline
100,100& &0.8020 &0.9580 &0.9563 &0.9630 \\
100,200& &0.8170 &0.9570 &0.9589 &0.9670 \\
200,200& &0.8140 &0.9850 &0.9842 &0.9870 \\
200,500& &0.8510 &0.9840 &0.9840 &0.9880 \\\hline
& &$\rho=0$ &$\alpha=0$ &$\beta=0.3$ & \\\hline
100,100& &0.7910 &0.9620 &0.9671 &0.9640 \\
100,200& &0.8090 &0.9670 &0.9674 &0.9780 \\
200,200& &0.8150 &0.9900 &0.9904 &0.9890 \\
200,500& &0.8330 &0.9890 &0.9890 &0.9950 \\\hline
& &$\rho=0.7$ &$\alpha=0.3$ &$\beta=0.3$ & \\\hline
100,100& &0.6670 &0.9740 &0.9766 &0.9830 \\
100,200& &0.7670 &0.9790 &0.9801 &0.9910 \\
200,200& &0.7910 &0.9940 &0.9940 &0.9920 \\
200,500& &0.7900 &0.9920 &0.9920 &0.9950 \\\hline
\end{tabular}
\end{center}
\end{table}
\begin{figure}[htbp]
\centering
\subfigure[$(\rho,\alpha,\beta)=(0.7,0,0)$]{
\includegraphics[width=7.0cm]{corr_100_N_T_100.eps}
}
\quad
\subfigure[$(\rho,\alpha,\beta)=(0,0.3,0)$]{
\includegraphics[width=7.0cm]{corr_010_N_T_100.eps}
}
\quad
\subfigure[$(\rho,\alpha,\beta)=(0,0,0.3)$]{
\includegraphics[width=7.0cm]{corr_001_N_T_100.eps}
}
\quad
\subfigure[$(\rho,\alpha,\beta)=(0.7,0.3,0.3)$]{
\includegraphics[width=7.0cm]{corr_111_N_T_100.eps}
}
\caption{Plots of the frequency of the estimated break points among 1000 replications for DGP 1.A and $N=100,T=100$.}
\label{1A_NT100}
\end{figure}
\begin{figure}[htbp]
\centering
\subfigure[$(\rho,\alpha,\beta)=(0.7,0,0)$]{
\includegraphics[width=7.0cm]{corr_100_N_T_500.eps}
}
\quad
\subfigure[$(\rho,\alpha,\beta)=(0,0.3,0)$]{
\includegraphics[width=7.0cm]{corr_010_N_T_500.eps}
}
\quad
\subfigure[$(\rho,\alpha,\beta)=(0,0,0.3)$]{
\includegraphics[width=7.0cm]{corr_001_N_T_500.eps}
}
\quad
\subfigure[$(\rho,\alpha,\beta)=(0.7,0.3,0.3)$]{
\includegraphics[width=7.0cm]{corr_111_N_T_500.eps}
}
\caption{Plots of the frequency of the estimated break points among 1000 replications for DGP 1.A and $N=500,T=500$.}
\label{1A_NT500}
\end{figure}
\begin{figure}[htbp]
\centering
\subfigure[$(\rho,\alpha,\beta)=(0.7,0,0)$]{
\includegraphics[width=7.0cm]{full_rank_corr_100_N_T_100.eps}
}
\quad
\subfigure[$(\rho,\alpha,\beta)=(0,0.3,0)$]{
\includegraphics[width=7.0cm]{full_rank_corr_010_N_T_100.eps}
}
\quad
\subfigure[$(\rho,\alpha,\beta)=(0,0,0.3)$]{
\includegraphics[width=7.0cm]{full_rank_corr_001_N_T_100.eps}
}
\quad
\subfigure[$(\rho,\alpha,\beta)=(0.7,0.3,0.3)$]{
\includegraphics[width=7.0cm]{full_rank_corr_111_N_T_100.eps}
}
\caption{Plots of the frequency of the estimated break points among 1000 replications for DGP 1.B and $N=100,T=100$.}
\label{1B_NT100}
\end{figure}
\begin{figure}[htbp]
\centering
\subfigure[$(\rho,\alpha,\beta)=(0.7,0,0)$]{
\includegraphics[width=7.0cm]{full_rank_corr_100_N_T_500.eps}
}
\quad
\subfigure[$(\rho,\alpha,\beta)=(0,0.3,0)$]{
\includegraphics[width=7.0cm]{full_rank_corr_010_N_T_500.eps}
}
\quad
\subfigure[$(\rho,\alpha,\beta)=(0,0,0.3)$]{
\includegraphics[width=7.0cm]{full_rank_corr_001_N_T_500.eps}
}
\quad
\subfigure[$(\rho,\alpha,\beta)=(0.7,0.3,0.3)$]{
\includegraphics[width=7.0cm]{full_rank_corr_111_N_T_500.eps}
}
\caption{Plots of the frequency of the estimated break points among 1000 replications for DGP 1.B and $N=500,T=500$.}
\label{1B_NT500}
\end{figure}
\vspace{-1em}
\section{Empirical Application}
\vspace{-1em}
\subsection{Macroeconomic data}
In the first empirical application, we apply our proposed method to a U.S. macroeconomic dataset (\cite{Stock2012}) to detect
the possible structural breaks in the underlying factor model. We use the dataset adopted by \cite{Cheng2016}, which comprises monthly observations of 102 U.S. macroeconomic variables. The sample begins
after the Great Moderation and ranges from 1985:01 to 2013:01 $(T = 337)$. Following \cite{Bai2017}, we focus on the subsample period between 2001:12 and 2013:01 $(T=134,N=102)$ because the complete
data may have multiple breaks.
\cite{Cheng2016} find that 2007:12 is a single-break date, and that the pre-break and post-break subsamples have one factor and two or three factors, respectively. Following \cite{Cheng2016},
\cite{Bai2017} also set the number of factors equal to one and two for the pre- and post-break subsamples, respectively. Then, they implement the LS estimation and obtain the estimated break point
$\hat{k}=2008:12$.
To implement our QML method, we first use Bai and Ng's information criterion IC1 and determine three pseudo-factors in the complete sample. Based on this result, we compute our QML estimator and obtain
2007:07 as the estimated break point, using which we split the sample into pre- and post-break subsamples. IC1 of \cite{Bai2002} detects two pre-break and three post-break factors. Based on the numbers of pre- and post-break factors and that of pseudo-factors, we can conclude that a new factor emerges after the break, so the QML estimator is consistent based on Theorem \ref{consistency}.
\subsection{Stock data}
The second empirical application uses the weekly rate of return for Nasdaq 100 Index from April 18, 2019, to October 1, 2020. As all companies have data starting from April 18, 2019, we choose that as
the start date. Traditionally, the index is limited to 100 common-stock issues, with only one issue allowed per issuer. Now, the index is limited to 100 issuers, some of which may have multiple issues as
index components. The current index has 103 components, representing 100 issuers, four of which are from China: Baidu, JD.com, Ctrip, and NetEase. Thus, the sample size is $T=76$ and $N=103$. As IC1 and
IC2 of \cite{Bai2002}, the methods proposed by \cite{Onatski2010}, \cite{Ahn_Horenstein2013}, and \cite{Fan2020} yield different numbers of pseudo-factors for the sample, we use different number of
factors $r=2,3,4,5,6,7$ to estimate the break date by using the QML method, and find that the estimated break date always falls in the week of February 20, 2020. This result agrees with that obtained
using the method developed by \cite{Baltagi2017}. In fact, the stock market began to fall sharply in the week of February 20, 2020, and two weeks later, the circuit breaker was triggered and U.S. stock market trading halted for a couple of times. Thus, the
factor loading matrix appears to have changed in the early days of the epidemic.
\vspace{-1em}
\section{Conclusions}
\vspace{-1em}
We study the QML method for estimating the break point in high-dimensional factor models with a single structural change. We consider three types of changes and develop an asymptotic theory for the
QML estimator.
We show that the QML estimator is consistent when the covariance matrices of the pre- or post-break factor loadings, or both, are singular.
In addition, the estimation error of the QML estimator is $O_p(1)$ when there is a rotation type of change in the factor loading matrix. We also derive the limiting distribution of the estimated break point in this case.
Moreover, our QML estimator is computationally easy and fast because the eigendecomposition is conducted only once.
The simulation results validate the suitable performance of the QML estimator.
We use the proposed method to estimate the break point for U.S. macroeconomic data and stocks data.
\newpage