\section{introduction}

Factor models are a special kind of latent-variable models that are widely used in economics and other disciplines of research.
Due to the ubiquitous dependence across high-dimensional economic variables, it is appealing for explaining the variation of the high-dimensional economic measurements from only a small number of latent common factors.
Well-known examples include the arbitrage pricing theory \citep{Ross1976}, rank of demand systems \citep{Lewbel1991}, multiple-factor models \citep{FF1993}, components analysis in economic activities  \citep{GH1999}, and diffusion index \citep{SW2002}.
Successful applications rely heavily on the ability to correctly specify the number of factors, which is usually unknown in practice. Therefore, determining the number of factors is an essential step in applying factor models.

 One approach to this end is to utilize the eigen-structure of data matrix; see for example  \cite{F1977}, \cite{Schott1994},  \cite{YW2003}, \cite{Onatski2010} and \cite{LL2016}, amongst  others.
However, these methods require strong assumptions on the separation of eigenvalues,  and thus most of them are only applicable to fixed dimensional data.
More popular approaches are cross-validation (CV) and information criteria (IC).
{For factor model, several cross-validation methods have been proposed}; see, for example, \cite{Wold1978} and \cite{EK1982}, amongst  others.
\cite{BKSK2008} comprehensively reviewed these CV-based approaches, and concluded that most were not statistically consistent.
In addition, they recommended an alternative method based on the expectation maximization (EM) algorithm, which turned out to be extremely  computationally expensive.
As an alternative, the ICs are  more computationally effective and consistency can be guaranteed. As far as we know, \cite{CD1997} was the first paper that used ICs to determine the  number of factors.
Subsequently, other methods based on the Akaike information criterion (AIC) and Bayesian information criterion (BIC) \citep{SW1998,FLHR2000, BaiNg2002, LLS2017}  were proposed for different settings of sample size $n$ and dimension of variable, $p$.
When the variance of the idiosyncratic component is relatively small, these methods are able to provide very efficient estimates of the number of factors \citep{BaiNg2002}.
However, the IC-based approaches are usually not fully data-driven since their penalties depend on predetermined tuning parameters.
Even though several guides are proposed to mitigate the effect of  tuning parameters, the performances of these methods are not stable when the signal-to-noise ratio is relatively small \citep{Onatski2010}, i.e., the variation of the idiosyncratic component is relatively large, or when the idiosyncratic component has heavy tails in distribution.
Unfortunately, both of which are common in financial data.

In this paper, we propose to estimate the number of factors by a computationally-efficient CV method.
Note that the conventional CV methods, e.g., those used in the linear regression, only leave one or several \emph{observations} out.
In contrast, factor models involve a matrix, where both observations (or rows) and variables (or columns) play similar roles mathematically in the analysis.

Thus, we propose a double cross-validation (DCV) method that applies CV first to the \emph{observations} and then to the \emph{variables} in the {observations}.
Theoretically, we show that the method is consistent for high-dimensional data and allows dependency between the idiosyncratic errors.
The method  is thus applicable to approximate factor models.
In addition, computationally, the second CV can be easily calculated by the linear regression error using the full data \citep{Shao1993}, while the first CV can be done in a $K$-fold manner.
Simulation studies show that the proposed approach performs satisfactorily even when the idiosyncratic  error has relatively big variance or homoscedasticity, or heavy tail,  which are especially relevant for economic data.

The rest of this paper is organized as follows. Section~\ref{sec:method} reviews the basic notation of factor analysis.
Section~\ref{DCVsection} presents the DCV approach.
Section~\ref{sec:theory} establishes statistical consistency of the approach, while the proofs of those properties are given in the Appendix.
Sections~\ref{sec:simulation}--\ref{sec:application} present a set of simulation studies and  an empirical application to assess the finite sample performance of DCV. Concluding remarks are given in Section~\ref{sec:discussion}.

\section{Factor model and its  factors and  loadings}\label{sec:method}

For a $p$-dimensional random vector ${\bm x}=(X_1, X_2,...,X_p)^\top $, the factor model assumes the following \emph{bilinear} structure
\begin{equation*}
{\bm x}={\bm L^0{\bm f}^0}+{\bm e},
\end{equation*}
where ${\bm L}^0\in\mathbb{R}^{p\times d_0}$ is the loading matrix, ${\bm f}^0=(f_{1}, f_{2}, ... , f_{d_0})^\top $ is  the vector of $ \ d_0 \  $ common factors. Together ${\bm L^0{\bm f}^0}$ is called the common component of the data,  and ${\bm e}=(e_{1},e_{2},...,e_{p})^\top $ the idiosyncratic component. Let $\boldsymbol\Sigma_{\bm e}=\mathbb{E}[\bm e\bm e^\top]$, $\boldsymbol\Sigma_{\bm f_0}=\mathbb{E}[{\bm f_0\bm f_0^\top}]$, and $\boldsymbol\Sigma_{\bm x}=\mathbb{E}[\bm x\bm x^\top]$. Under this model, the variation of the $p$-dimensional variable is mainly generated by a small number of common factors.
In the statistical literature, the idiosyncratic component are usually assumed to be independent of the common component, and have diagonal covariance matrix, in which case {$\boldsymbol\Sigma_{\bm x}={\bm L}^0\boldsymbol\Sigma_{\bm f_0}{{\bm L}^{0}}^\top +\boldsymbol\Sigma_{\bm e}$}. However, such independence assumption and diagonal structure are usually not realistic in financial data or macroeconomic data where factor analysis are often applied. In this paper, we consider the approximate factor model where the idiosyncratic components can be weakly dependent.


Let  ${\bm X}=({\bm x}^\top _1,{\bm x}^\top _2,...,{\bm x}^\top _n)^\top $ be random samples of  $n$ observations on $ p$ variables. We have the following matrix form for the factor model:
\begin{equation}  \label{FM}
{\bm X}={\bm F}^0{{\bm L}^{0}}^\top +{\bm E},
\end{equation}
with ${\bm F}^0=({\bm f}_1^0,{\bm f}_2^0,...,{\bm f}_n^0)^\top $ and  ${\bm E}=({\bm e}^\top _1,{\bm e}^\top _2,...,{\bm e}^\top _n)^\top $.
According to the eigen-decomposition, ${\bm X^\top \bm X}$ has the  expression
\begin{equation*}
{\bm X^\top \bm X}=\boldsymbol\Phi\boldsymbol\Lambda\boldsymbol\Phi^\top,
\end{equation*}
where $\boldsymbol\Lambda$ is a $p\times p$  diagonal matrix  with diagonal  entries being  the eigenvalues of $\bm X^\top\bm X$, $\boldsymbol\Phi=(\boldsymbol\phi_1,\boldsymbol\phi_2,...,\boldsymbol\phi_n)\in\mathbb{R}^{n\times n}$ is the eigenvector matrix of ${\bm X^\top \bm X}$ such that $\boldsymbol\Phi^\top\boldsymbol\Phi={\bm I}_p$. Then, with working number of factors $d$, the estimators of factor loadings and factors  are respectively
\begin{equation}\label{FL1}
\widetilde{\bm L}^d=\sqrt{p}(\boldsymbol\phi_1,\boldsymbol\phi_2,...,\boldsymbol\phi_d) \ \  \mbox{and} \ \  \widetilde{\bm F}^d=\frac{1}{{p}}{\bm X}\widetilde{\bm L}^d.
\end{equation}
Denote by $\text{Tr}(\bm A)$ the trace of matrix $\bm A$. It is easy to see that
\begin{equation*} \label{LF}
(\widetilde{\bm F}^{d},\widetilde{\bm L}^{d})=\underset{\substack{{\bm F}^d\in\mathbb{R}^{n\times d},{\bm L}^d\in\mathbb{R}^{p\times d}\\{{\bm L}^{d}}^\top {\bm L}^d/p={\bm I}_d}}{\text{argmin}}\{\text{Tr}[({\bm X}-{\bm F}^d{{\bm L}^{d}}^\top )({\bm X}-{\bm F}^d{{\bm L}^{d}}^\top )^\top ]\}.
\end{equation*}

\section{Double Cross Validation for Number of Factors}\label{DCVsection}


Our basic idea to specify the number of factor $d_0$, is based on the prediction error in a cross-validated manner. Write the model element-wise as
\begin{equation}\label{reg1}
{x}_{is}={{\bm f}_i^0}^\top {\bm l}^0_s+{{e}}_{is}, \ \ \ \ i=1, ..., n, \ \  s=1,2,...,p.
\end{equation}
CV needs to make prediction of $ {x}_{is} $ based on other elements except for itself.

Because neither $ {{\bm f}_i^0} $ nor $  {\bm l}^0_s $ in \eqref{reg1} is observable, we implement the prediction by two-stage fitting: leave-\emph{observation}-out  and leave-\emph{variable}-out.

The first stage is to estimate $  {\bm l}^0_s $ from the data by the $K$-fold CV. Divide the rows of $ {\bm X} $ into $K$ folds, $ 1< K \le n  $. Denote them by $M_1, ..., M_K$. Let $n_k=\#M_k$ be the number of elements in $M_k$, and ${\bm X}_{-M_k}$ be the sub-matrix of $ {\bm X} $ with  rows in $ M_k $ being removed. Apply  the eigen-decomposition to ${\bm X}_{-M_k}$ to obtain corresponding matrices $\widetilde{\bm F}^{k,d}$ and $\widetilde{\bm L}^{k,d}$ as in \eqref{FL1}.
Here, we use the superscript to highlight the fact that the estimator is obtained from the data with the rows in $M_k$ removed, with working number of factors $d$. Similar to \cite{BaiNg2002}, we rescale the estimated factor loading matrix and let
\begin{equation*}
\widehat{\bm L}^{k,d}= (\frac{1}{n-n_k}{\bm X}^\top_{-M_k}{{\bm X}_{-M_k}}) \widetilde{\bm L}^{k,d} = (\widehat{\bm l}^{k,d}_1, ..., \widehat{\bm l}^{k,d}_p).
\end{equation*}

In the second stage, we replace $ {\bm l}^0_s $ in \eqref{reg1} by $ \widehat{\bm l}^{k,d}_s $ and rewrite \eqref{reg1} as a regression model
\begin{equation}\label{reg2}
{x}_{is}={{\bm f}_i^0}^\top \widehat{\bm l}^{k,d}_s + {{e}}'_{is}, \ \ \ \ i \in M_k, \ \  s=1,2,...,p,
\end{equation}
where $ e'_{is} $ is the regression error in place of $ e_{is} $ due to the replacement. This time,  $\widehat{\bm l}^{k,d}_s, s=1,2,...,p$ are known and treated as the `regressors', but ${\bm f}^0_i$ is treated as  `regression coefficients' that need to be estimated.
By leaving $x_{is}$ out,  we estimate ${\bm f}^0_i$ by
\begin{equation*}
 \widehat {{\bm f}}_{i,s}^{k,d} = \underset{\bm f}{\text{argmin}} \sum_{t\neq s, t=1}^p ( {x}_{it} -\bm f^\top \widehat{\bm l}^{k,d}_t)^2 ,\  i \in M_k.
\end{equation*}
Thus we can predict $x_{is}$ by $(\widehat {{\bm f}}_{i,s}^{k,d})^\top \widehat{\bm l}^{k,d}_s $.
The average squared prediction error for $\bm x_i = (x_{i1}, ..., x_{ip})^\top$ is
\begin{equation*}
 V_i^{k,d} = \frac{1}{p}\sum\limits_{s=1}^p [x_{is} - (\widehat {{\bm f}}_{i,s}^{k,d})^\top \widehat{\bm l}^{k,d}_s]^2, \ \  i \in M_k.
\end{equation*}
The calculation of $  V_i^{k,d} $ can be much simplified  as shown in \cite{Shao1993}, i.e.,
\begin{equation*}
V_i^{k,d} =\frac{1}{p}\sum\limits_{s=1}^p(1-w^{k,d}_s)^{-2}(x_{is}-(\widehat{\bm f}_{i}^{k,d})^\top\widehat{\bm l}^{k,d}_s)^2,
\end{equation*}
where  $\widehat{\bm f}_i^{k,d}=(\widehat{\bm L}^{k,d\top}\widehat{\bm L}^{k,d})^{-1}\widehat{\bm L}^{k,d\top}\bm x_i$ is the conventional least squares estimator for \eqref{reg2} and $w^{k,d}_s$  is the $s$-th diagonal element of projection matrix $\mathcal{P}_{\widehat{\bm L}^{k,d}}=\widehat{\bm L}^{k,d}(\widehat{\bm L}^{k,d\top}\widehat{\bm L}^{k,d})^{-1} \times \widehat{\bm L}^{k,d\top}$.
Finally, consider the averaged  prediction error over all the elements
\begin{equation*}
DCV(d)=\frac{1}{n} \sum\limits_{k= 1}^K \sum\limits_{i \in M_k}^nV_i^{k,d}.
\end{equation*}
 Let $  d_{\max}$ be a fixed positive integer that is large enough such that  $  p> d_{\max} > d_0$.  The DCV estimator for the number of factors is given by
\begin{equation}\label{definek}
\widehat{d}_{DCV}=\underset{0\leq d\leq d_{\max}}{\text{argmin}}DCV(d).
\end{equation}


The above approach has similarity with the vanilla row-wise CV \citep{BKSK2008}, which also predicts ${\bm x}_{i}$  by two steps. Its first step is exactly the same as ours. However, in the second step, row-wise CV uses the full data instead of cross validation method,  leading to overfitting and inconsistency of the method. To fix this problem, \cite{EK1982}  suggested the element-wise cross validation. Although their methods solve the problem of overfitting, they are costly and involve immense computation.

Our DCV  offers a simpler solution. By utilizing CV in both the first stage of leaving ``observations'' out, and the second regression step of leaving   ``variables'' out, our method  guarantees that the prediction of each element does not use any information from itself, and thus the consistency is ensured as shown in the next section. In addition, the implement of $K$-fold CV can significantly reduce the computational complexity, especially when the number of rows is large. It is interesting to see that the $K$-fold CV in the first stage will also help in simplifying the calculation in the second stage. Because of this appealing nature, to further facilitate the calculation, we can transpose $ {\bm X} $ when $n<p$.






\section{Asymptotic consistency of the estimation}\label{sec:theory}
In this section, we investigate the consistency of our derived estimator.  We start with the following assumptions.
\begin{ass}\label{ass1}
The eigenvalues of ${\bm F}^{0\top}{\bm F}^0/n$ are bounded away from zero and infinity, and  that $\mathbb{E}\|{\bm f}^0_i\|^4<\infty$ for $ i=1, ..., n$.
\end{ass}
\begin{ass}\label{ass2} The eigenvalues of ${\bm L}^{0\top}{\bm L}^0/p$ are bounded away from zero and infinity, and that $\mathbb{E}\|{\bm l}^0_i\|^4<\infty$ for  $i=1, ..., p$.
 \end{ass}



\begin{ass}\label{ass3} There exists a constant $M<\infty$ such that
\begin{enumerate}
\item[1.]  $\mathbb{E} e_{it}=0$  and  $\mathbb{E} e_{it}^4\leq M$ for $1\leq i\leq n$ and $1\leq t\leq p$;

  \item[2.] $\gamma_p(i,j) = \mathbb{E}({\bm e}_i^\top{\bm e}_j/p)\leq M$  and $\frac{1}n\sum\limits_{j=1}^n\gamma^2_p(i,j)\leq M$ for all $i$;

 \item[3.] $\mathbb{E}({\bm e}_{is}{\bm e}_{it})=\tau_{st,i}$ with $\tau_{st,i}\leq |\tau_{st}|$ for some $\tau_{st}$ and for all $i$; in addition,  $\frac{1}p\sum\limits_{s=1}^p\sum\limits_{t=1}^p\tau_{st}$ $\leq M$;

 \item[4.] $\mathbb{E}({\bm e}_{is}{\bm e}_{jt})=\tau_{ij,st}$ and $\frac{1}{np}\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_{s=1}^p\sum\limits_{t=1}^p\tau_{ij,st}\leq M$;

\item[5.]  $\frac{1}{  n }\mathbb{E}|\sum\limits_{j=1}^n[e_{jt}e_{js}-\mathbb{E}(e_{jt}e_{js})]|^2\leq M$ for all $(s,t)$;

\item[6.]  $\|\frac{1}{\sqrt{n}}\sum\limits_{i=1}^n{\bm f}^{0}_{i}e_{is}\|^2\leq M$ for $1\leq s\leq p$;

\item[7.]  $\|\frac{1}{\sqrt{p}}\sum\limits_{s=1}^p{\bm l}^{0}_{s}e_{is}\|^2\leq M$ for $1\leq i\leq n$.

\end{enumerate}
\end{ass}


Assumptions \ref{ass1}--\ref{ass3} are commonly used in the  analysis of  factor models, for example, in \citet{BaiNg2002} and \citet{LLS2017}.  Assumptions \ref{ass1} and  \ref{ass2} together ensure that each factor plays a nontrivial role in contributing to the variation of ${\bm X}$. Unlike the strict factor model that assumes all entries of $\bm E$ are I.I.D., Assumption~\ref{ass3} allows weak dependency for elements of  $\bm E$, making our methods applicable to the approximate factor models.  These assumptions also indicate that  all the  eigenvalues of the covariance matrix of common components would dominate the eigenvalues corresponding to idiosyncratic components, which is crucial to the identifiability for the approximate factor models.

\begin{ass}\label{ass4}
$\sup\limits_k n_k=o(n^{1/3})$ and
 $K\sup\limits_kn_k/n< \infty$, where $n_k=\#M_k$ is the number of elements in $M_k$.

\end{ass}
Assumption \ref{ass4} is a weak condition on $K$-fold CV in the first stage. The number of elements in each fold can either be fixed or tend to infinity. In particular, this assumption includes the leave-one-out CV.
\edef\oldass{\the\numexpr\value{ass}+1}

\setcounter{ass}{0}
\renewcommand{\theass}{\oldass.\alph{ass}}

\begin{ass}\label{ass5a}  The idiosyncratic components satisfies
\begin{enumerate}
\item[1.] $\mathbb{E}[e_{it}^2\mid {\bm F}^0_{-M_k}{\bm L}^{0\top}]=\sigma^2$ for $i\in M_k$, where ${\bm F}^0_{-M_k}$ is a submatrix of ${\bm F}^0 $ with rows in $M_k$ being removed;

\item[2.] $\lambda_1({\bm E^\top}{\bm E}/n)<2 \sigma^2$.
\end{enumerate}

\end{ass}



\begin{ass}\label{ass5b} Let $p/n\to \rho\in(0,\infty)$.  There exists a $p\times p$ positive definite matrix $\boldsymbol\Sigma_p$ such that ${\bm e}_i=\boldsymbol\Sigma_p^{1/2}{\bm y}_i$, where ${\bm y}_i=( y_{i1},y_{i2},...,y_{ip})^\top$,
\begin{enumerate}
\item[1.] $\mathbb{E} y_{it}=0$, $\mathbb{E} y_{it}^4<\infty$,  and $y_{it}$ are I.I.D. for $1\leq i\leq n, 1\leq t\leq p$;

\item[2.]  ${\bm y}_i$ is independent of ${\bm F}^0_{-M_k}$;

 \item[3.] The spectral distribution of $\boldsymbol\Sigma_p$ is convergent, i.e., there exists a distribution function $H$ such that $\frac{1}{p}\sum\limits_{s=1}^p{\mathds{1}}(\lambda_s(\boldsymbol\Sigma_p)\leq t)\to H(t)$ as $p\to \infty$;

\item[4.] $2\text{\normalfont{diag}}(\boldsymbol\Sigma_p) - \boldsymbol\Sigma_p $ is  a positive definite matrix, where $\text{\normalfont{diag}}(\boldsymbol\Sigma_p)$ is a diagonal matrix that has the same diagonal elements as $\boldsymbol\Sigma_p$.
\end{enumerate}
\end{ass}
\let\theass\origtheass


Assumptions \ref{ass5a} and  \ref{ass5b} are two technical assumptions, either of which ensures the viability of the proposed method.  They also allow weak dependence among idiosyncratic components as well as the observations, which generalize the I.I.D. assumption required in the strict factor model.  In comparison, \cite{BaiNg2002} needs weaker assumptions on  the dependence than ours, possibly due to the relatively large penalty in their methods. Of course, this large penalty is at the cost of selecting inappropriately fewer factors as shown in our simulation study.






\begin{thm} \label{thm1}
Suppose  Assumptions \ref{ass1}--\ref{ass4} hold. If  either Assumption \ref{ass5a}  or \ref{ass5b} holds, then the DCV estimator \eqref{definek} is consistent as both $ n $ and $ p $ tend to infinity, i.e.,
\begin{equation*}
\lim\limits_{n,p\to\infty}Prob(\widehat{d}_{DCV}=d_0)=1.
\end{equation*}
\end{thm}

{\color{black}
\begin{rem}  Recent empirical findings suggest that the factor structure of economic data may change with  sample size and dimensionality of variable \citep{Onatski2010,  JLN2015,LLS2017}. In fact, Theorem \ref{thm1} can be easily extended to the case that $ d_0 $ varies slowly with both $ n $ and $ p $. {\color{red}DCV} is thus applicable to practical data with changing structure of factorization.
\end{rem}}

\section{Simulation Study}\label{sec:simulation}
In this section, we conduct simulation studies to compare the finite sample performances of the DCV with the other methods, including
the {panel} criterion $IC_1$ \citep{BaiNg2002} and $Ladle$ \citep{LL2016}.  
As shown in \cite{LL2016},  $Ladle$ estimator outperforms all the other eigen-structrue based methods in determining the number of factors, and is thus used  as a representative of those methods.
Similarly, for the IC-based methods, \cite{BaiNg2002} showed that their criterion outperforms other information criteria such as $AIC$ and $BIC$.
We do not consider other cross-validation methods as they are either inconsistent or suffer from excessive computational burden.
For our $DCV$, we consider both leave-one-out CV ($DCV_1$) and 10-fold CV ($DCV_{10}$) in the first stage.
Matlab codes for the calculations are available at \href{https://github.com/XianliZeng/Double-Cross-Validation}{https://github.com/XianliZeng/Double-Cross-Validation}.


\begin{figure*}[ht!]
\vspace{-0.3cm}
\centerline{\includegraphics[width=1.0\linewidth]{model1.eps} } \vspace{-0.5cm}
 \caption{\footnotesize The relative frequencies of correctly selecting the number of factors for model~\eqref{SM1} with independent Gaussian idiosyncratic errors (E1). In each panel, we use red line for $DCV$, black line for $DCV_{10}$, {blue} line for $IC_1$ and  green line for \textit{Ladle}.\vspace{-0.5cm}
}
\label{fig1}
\end{figure*}
\begin{figure*}[ht!]
\vspace{-0.3cm}
\centerline{\includegraphics[width=1.0\linewidth]{model2.eps} } \vspace{-0.5cm}
 \caption{\footnotesize The relative frequencies of correctly selecting the number of factors for  model \eqref{SM1} with independent $t$-distributed idiosyncratic errors (E2). In each panel, we use red line for $DCV$, black line for $DCV_{10}$, blue line for $IC_1$ and  green line for \textit{Ladle}.\vspace{-0.5cm}
}
\label{fig2}
\end{figure*}
\begin{figure*}[ht!]
\vspace{-0.3cm}
\centerline{\includegraphics[width=1.0\linewidth]{model3.eps} } \vspace{-0.5cm}
 \caption{\footnotesize The relative frequencies of correctly selecting the number of factors for   model~\eqref{SM1} with heteroskedastic idiosyncratic errors (E3). In each panel, we use red line for $DCV$, black line for $DCV_{10}$, blue line for $IC_1$ and  green line for \textit{Ladle}.\vspace{-0.5cm}
}
\label{fig3}
\end{figure*}

\begin{figure*}[ht!]
\vspace{-0.3cm}
\centerline{\includegraphics[width=1.0\linewidth]{model4.eps} } \vspace{-0.5cm}
 \caption{\footnotesize The relative frequencies of correctly selecting the number of factors for  model~\eqref{SM1} with serial correlated  idiosyncratic errors (E4). In each panel, we use red line for $DCV$, black line for $DCV_{10}$, blue line for $IC_1$ and  green line for \textit{Ladle}.\vspace{-0.5cm}
}
\label{fig4}
\end{figure*}

\begin{figure*}[ht!]
\vspace{-0.3cm}
\centerline{\includegraphics[width=1.0\linewidth]{model5.eps} } \vspace{-0.5cm}
 \caption{\footnotesize The relative frequencies of correctly selecting the number of factors for  model~\eqref{SM1} with cross-sectional correlated idiosyncratic errors (E5). In each panel, we use red line for $DCV$, black line for $DCV_{10}$, blue line for $IC_1$ and  green line for \textit{Ladle}.\vspace{-0.5cm}
}
\label{fig5}
\end{figure*}

We simulate data from a factor model with $d_0 = 5$:
\begin{equation}
 x_{is}=\sum_{j=1}^{5}f_{ij}l_{sj}+\sqrt{\theta}e_{is}; \quad f_{ij},l_{sj}\sim N(0,1)  \ \ i=1, ..., n, \ \  s=1,2,...,p.  \label{SM1}
\end{equation}
For the idiosyncratic errors $ e_{is} $, the following settings are considered:
\begin{description}
\item[] (E1) independent Gaussian:  $ e_{is}\sim N(0,1)$;
\item[] (E2) independent t-distributed:  $ e_{is}\sim t_3 $;
\item[] (E3) heteroskedastic:  $ e_{is} \sim N(0,\delta_s), $ $ \delta_s = 1 $ if $s$ is odd,  and 2 if $s$ is even;
\item[] (E4) serial correlated: $ e_{is}=0.3 e_{i,s-1}+\nu_{is} $,  $\nu_{is}\sim N(0,1)$;
\item[] (E5) cross-sectional correlated: $  e_{is}=\sum\limits_{j=-10}^{10} 0.15^j \nu_{i-j,s}$,  $\nu_{is}\sim N(0,1)$.
\end{description}

In all simulations, we select $d$ from $\{1,2,...,d_{\max}\}$ where $d_{\max}= 8 $ as in \cite{BaiNg2002}.
We consider the combinations of { $n$ and $p$} with $n=\{40,160,640\}$ and $p=\{30,90, 270\}$.
To see the effect of  the signal-to-noise ratio on the methods, $\theta$ varies in a range of intervals depending on different settings of idiosyncratic errors and different $ n $ and $p$.
These ranges are made clear in each panel of the figures below.
All comparisons  were evaluated based on 1000 replications.
Figures~\ref{fig1} to~\ref{fig5} summarize the simulation results, where the curves represent the frequencies of selecting the right number of factors at different noise level $\theta$.

Generally, the simulation results suggest that the performance of all the methods becomes better when $n$ and $p$ get larger, which lends support to the asymptotic consistency of the methods.
$Ladle$ usually shows better performance than $IC_1$, possibly due to its utilization of  both eigenvalues and eigenvectors.
Notably, $DCV_1$ and $DCV_{10}$ stand out and perform the best in most scenarios.
They possess a much higher frequency of correct selection when $\theta$ is large, implying their ability to select the right number of factors even when the {factors' contribution to the variation of the variables,  i.e., signal-to-noise ratio,} is low, which is particularly the situation in economic data.
Moreover, $DCV_1$ and $DCV_{10}$ present almost
identical performance for all the models, implying that the performance of the DCV is not sensitive to the choice of~$K$.

For different settings of the idiosyncratic errors, we have some additional observations.
Figure \ref{fig2} shows that $IC_1 $ is not so robust to the data with heavy-tail especially when the size of the data is small.
Figure \ref{fig3} implies that $IC_1 $ is also affected adversely by the heteroskedastic errors.
Note that both heavy tail and  heteroskedasticity are stylized facts in financial data.
In contrast, for these two types of data,  $DCV_1$ and $DCV_{10}$ still demonstrate very stable and efficient ability in selecting the number of factors.
In Figures \ref{fig4} and \ref{fig5}, where the idiosyncratic errors have serial or cross-sectional correlation,  DCV has inferior performance than $IC_1$ and $ Ladle $ when the variance of idiosyncratic components is small and  both $ n $ and $ p $ are small.
However, DCVs quickly gains advantage when either $ n $ or $ p $ increases.



\section{An empirical application: Fama-French  three-factor model}\label{sec:application}
\begin{figure*}[p]
\vspace{-0.5cm}
\centerline{Validation of three-factor model  based on value weighted return of 25 portfolios}
\centerline{\includegraphics[width=1\linewidth, height=0.56\linewidth]{figreal1.pdf}}

\centerline{Validation of three-factor model  based on value weighted return of 100 portfolios}
\centerline{\includegraphics[width=1\linewidth, height=0.56\linewidth]{figreal2.pdf}} \vspace{-0.5cm}
 \caption{\footnotesize The  frequencies of selected numbers of factors for the value weighted returns of 25 and 100 portfolios. Panels in each column represents one method; panels in the same rows are based on the same data in fixed period of years.
}
\label{fig6}
\end{figure*}

In asset pricing and portfolio management, the so-called three-factor model of \cite{FF1993} is widely applied to describe the stock returns.
It shows that the excess return of a stock or a portfolio ($R_{it}-{r_{ft}}$) can be satisfactorily explained by three common factors:  (1) excess return of the market portfolio ($R_{Mt}-r_{ft}$), (2) the outperformance of small-cap companies versus big-cap companies ($SMB$), and (3) the outperformance of companies with  high book to market ratio versus those with small book to market ratio ($HML$), i.e.,
\begin{equation*}
R_{it}-r_{ft}=\beta_1(R_{Mt}-r_{ft})+\beta_2SMB_t+\beta_{3}HML_t+\epsilon_{it},
\end{equation*}
where $ r_{ft} $ is the risk-free return.
Next,  we use the data provided by Professor  Kenneth R. French\footnote{http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data\_library.html} to mutually validate our methods and the three-factor model.
Stocks are divided into 6, 25 or 100 groups according to company's capitalization (market equity) and  Book-to-Market ratio, and portfolios are constructed with  value weighted daily returns or equally weighted daily returns in each group.
As those portfolios were compiled based on company's data at the end of June of every year, we set the first of July of every year as a starting date of a period.
It is understandable that the  portfolios change from year to year, thus we only consider periods of 1 to 3 years in our calculation.

Figure \ref{fig6} presents the frequencies of the numbers of factors detected by the $DCV_1$, $DCV_{10}$, $IC_1$ and $Ladle$ in two of the data sets with $d_{\max}$ set at $15$.
In most situations, all the  approaches identify three common factors for the returns of portfolios, which are in line with the \emph{three}-factor model of \cite{FF1993}.
In general,  DCV selects three factors most often amongst the approaches, showing its superior performance in practice.


\section{Concluding Remarks}\label{sec:discussion}
In the literature, both CV and IC based approaches were popularly used for determining the complexity of a model, such as  the number of factors.
As noticed in \cite{BKSK2008}, most CV approaches for factor models are not consistent due to their insufficient validation, while existing remedies for the consistency is very computationally expensive.
In contrast, ICs are consistent and easy to implement, but they depend on  predetermined penalty  functions \citep{BaiNg2002}.
Simulation results not reported here show that the penalty function in \cite{BaiNg2002} could  be modified to be more adaptive to the data.
In addition, ICs are unstable when the signal-to-noise ratio is relatively low \citep{Onatski2010}.

By validating the model twice, the proposed DCV method does not only ensure the consistency under mild conditions, but is also easy to implement. %
Because DCV is based on prediction error, it automatically selects the number of factors that balances the model complexity and stability.
Our simulation studies also demonstrate its superior efficiency over the existing methods most of the time.
The only exception is when there is strong serial dependence in the idiosyncratic errors.
However, this deficiency disappears when the sample size and dimension of the variables increase.
The advantages of our method are more  pronounced for data with relative large variation, heteroscedasticity  or heavy tails in the idiosyncratic errors. The method is thus particularly relevant for financial data.



\bibliographystyle{apalike}
