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.
59,868 characters
Composite Quantile Factor Model
\author{Xiao Huang\thanks{The author is grateful for valuable comments from Zheng Fang, Daniel Henderson, Junsoo Lee, Xiaochun Liu, Whitney Newey, Alejandro Sanchez-Becerra, Ruoxuan Xiong, Zhaoguo Zhan, and seminar participants at the University of Alabama and the 3rd Georgia Econometrics Workshop at Emory University. The author also thanks the Office of Research at the Kennesaw State University for computation support. Correspondence address: Department of Economics, Finance, and Quantitative Analysis, Coles College of Business, Kennesaw State University, GA 30144, USA. Email: [email removed].}}
\title{\Large Composite Quantile Factor Model}
\date{ \today}
\maketitle
\doublespace
\begin{abstract}
This paper introduces the method of composite quantile factor model for factor analysis in high-dimensional panel data. We propose to estimate the factors and factor loadings across multiple quantiles of the data, allowing the estimates to better adapt to features of the data at different quantiles while still modeling the mean of the data. We develop the limiting distribution of the estimated factors and factor loadings, and an information criterion for consistent factor number selection is also discussed. Simulations show that the proposed estimator and the information criterion have good finite sample properties for several non-normal distributions under consideration. We also consider an empirical study on the factor analysis for $246$ quarterly macroeconomic variables. A companion R package \texttt{cqrfactor} is developed.
\end{abstract}
\bigskip
\textbf{JEL Classification}: C18, C21.
\bigskip
\textbf{Keywords}: Composite quantiles, factor analysis, panel data.
\newpage
\normalsize
\doublespace
\section{Introduction} \label{intro}
Factor model is a useful statistical tool to describe data with unobserved and systematic components. Detailed textbook discussions on classical factor analysis for data with a fixed or small number of variables can be found in \cite{lawleymaxwell1971factoranalysis,anderson2003introduction}. Following the important work in \cite{stockandwatson1998nber,stockandwatson2002JASA,stockandwatson2002jbes,baiandng2002ecma,bai2003ecma} on high-dimensional panel data, the research on factor analysis in panel data has been extended to many directions, and there now exists a large body of literature on factor analysis for high-dimensional panel data. See \cite{baiandwang2016are} for a review on recent developments in this burgeoning field.
At the heart of factor analysis for high-dimensional panel data is the use of principal component analysis (PCA) method. The estimates for factors are chosen to be the normalized eigenvectors of the sample covariance matrix of the data, based on which we can estimate the factor loadings using the least-squares (LS) method. These eigenvectors coincide with the solutions in PCA, and the procedure's simplicity greatly contributes to its wide popularity as a research tool in empirical macroeconomics and finance.
Despite its simplicity, the PCA-based procedure typically imposes some higher-order moment conditions on the factor and error terms (see, e.g., the assumptions in \cite{stockandwatson2002JASA,bai2003ecma}) in order to obtain desirable asymptotic results. However, the sample covariance matrix on which PCA operates may be irrelevant for data with infinite (or very large) variances and PCA becomes invalid. Weakening or even removing these conditions will be appealing since many data in applications such as finance are either heavy tailed or of unknown nature. Two approaches for robust factor analysis have emerged in recent literature. \cite{chenetal2021qfm} introduce the quantile factor models (QFM), where both the factors and factor loadings are quantile-dependent. At the quantile position $ 0.5 $, QFM estimator can be interpreted as the least absolute deviation (LAD) estimator that may be robust to certain error distributions. \cite{andoandbai2020jasa} provide a more general framework that adds a regression component with heterogeneous coefficients. By assuming the factors and errors follow a joint elliptical distribution, \cite{heetal2022JBESnomoment} propose a second approach to replacing the sample covariance matrix in PCA with the spatial Kendall's tau matrix that can handle non-normal distributions with large or infinite variances, a situation in which the standard PCA fails.
This paper studies another approach to robust factor analysis and we term it composite quantile factor models (CQFM). Our approach is inspired by the interesting work in \cite{zouandyuan2008cqr}. \cite{zouandyuan2008cqr} notice that, in a linear regression with infinite error variance, the parameter estimator will no longer have root-$ n $ consistency or asymptotic normality; a robust procedure such as LAD can be used but its relative efficiency to the LS estimator can be very small. The authors propose to estimate the regression coefficients by simultaneously minimizing the standard quantile regression objective function at multiple quantile positions and call this procedure composite quantile regression (CQR). \cite{zouandyuan2008cqr} demonstrate the good finite sample properties of the CQR estimator for several non-normal error distributions.
Although factor analysis is different from linear regression, they are intrinsically connected (see, e.g., \cite{stockandwatson1998nber} for the use of LS method in deriving the solution to the factor model). If CQR works in linear regression, we conjecture a variant of it will also work in factor model. This paper studies the extension of CQR to the estimation of factor model by simultaneously minimizing the objective function at multiple quantiles. The resulting estimates are shown to have good finite sample properties under various non-normal error distributions. Because the CQFM visits different quantiles of data during estimation, the estimated factors can usually pick up more skewness (and kurtosis) information about the data, often resulting a better fit of the model. It is important to point out that, although CQFM uses the method of quantile regression, its estimates are for the mean factor model. This sets our paper apart from the work in \cite{andoandbai2020jasa,chenetal2021qfm}, where the goal is to estimate parameters at a specific quantile position.
We make the following contributions to the growing literature on panel factor analysis. First, we introduce CQFM as a new method to perform factor analysis on the mean factor model that can capture features of data at different quantiles. Second, we develop the asymptotic distribution results for the estimated factors and factor loadings; an information criterion is also developed to consistently select the factor number. Third, we provide extensive simulation evidence to show that CQFM works well for several non-normal error distributions. A special case of CQFM is when one chooses to optimize the objective function at a single quantile position, say $ 0.5 $. This reduces CQFM to the QFM in \cite{chenetal2021qfm}. In several simulation examples, we demonstrate the advantage of CQFM as a result of using information at multiple quantiles. We also develop an R package \texttt{cqrfactor} that implements the CQFM method in this paper. The \texttt{cqrfactor} package can be downloaded from \url{https://github.com/xhuang20/cqrfactor}.
The rest of the paper is organized as follows. Section 2 sets up the objective function for CQFM and discusses the estimation procedure and the asymptotic results. Section 3 discusses the information criterion for the selection of factor numbers. Section 4 presents all simulation results. Section 5 applies CQFM method to the modeling of the quarterly macroeconomic data in \cite{mccrackenandng2020FREDQD}. Section 6 concludes. The online supplement contains all proofs, additional figures and tables.
\section{Model estimation and the asymptotic results}
\subsection{The model and the algorithm}
Let $ Y_{it} $ be the observation at time $ t $ for the $ i $th cross-section unit. Consider the following factor model:
\begin{equation} \label{eq: factor model}
Y_{it} = \lambda_{0i}' F_{0t} + \varepsilon_{it}, \text{ for } i = 1, \cdots, N, t = 1,\cdots,T,
\end{equation}
where $F_{0t}$ is an $r \times 1$ vector of factors, $\lambda_{0i}$ is an $r \times 1$ vector of factor loading, and $\varepsilon_{it}$ is the error term. Both $F_{0t}$ and $\lambda_{0i}$ are unobserved, and the goal is to estimate them jointly. We assume the number of factors ($r$) is known. An information criterion will be developed to estimate $r$ consistently in \Cref{sec:factor number}. Rewrite \cref{eq: factor model} in matrix form to have
\begin{equation} \label{eq: factor model matrix form}
\underset{T \times N}{Y \vphantom{\Lambda_0'}} = \underset{T \times r}{F_0 \vphantom{\Lambda_0'}} \underset{r \times N}{\Lambda_0'} + \underset{T \times N}{\varepsilon \vphantom{\Lambda_0'}}.
\end{equation}
where $Y_{it}$ and $\varepsilon_{it}$ are elements of $Y$ and $\varepsilon$, respectively, and
\begin{equation}
F_0 = \big[F_{01}, \cdots,F_{0t},\cdots, F_{0T}\big]', \qquad
\Lambda_0 = \big[\lambda_{01},\cdots, \lambda_{0i}, \cdots, \lambda_{0N}\big]'.
\end{equation}
Let $ \tau$ be a quantile position with $0 < \tau < 1$ and $b_{0\tau}$ be the $ 100\tau\% $ quantile of $\varepsilon_{it}$. For the quantile factor model, we seek the estimates for $b_{0\tau}$, $F_{0t}$, and $\lambda_{0i}$ that minimize the following QFM objective function:
\begin{equation} \label{eq: QFM obj}
\frac{1}{NT} \sum_{i=1}^{N} \sum_{t=1}^{T} \rho_{\tau}\left(Y_{it} - b_{\tau} - \lambda_i'F_t\right),
\end{equation}
where $\rho_{\tau}(u) = u(\tau - \mathbf{I}(u \leq 0))$ is the check function in quantile regression. The estimates from \cref{eq: QFM obj} depend on $\tau$ and can be written as $\hat{F}_t(\tau)$ and $\hat{\lambda}_i(\tau)$, an approach adopted in \cite{andoandbai2020jasa,chenetal2021qfm}.
In CQFM, instead of estimating the model at a single quantile position $\tau$, we estimate the model simultaneously at multiple quantiles by choosing a sequence of $K$ quantiles, $0 < \tau_1<\tau_2<\cdots<\tau_K<1$, and minimizing the following objection function:
\begin{equation} \label{eq: CQFM obj}
\frac{1}{NT} \sum_{k=1}^{K}\sum_{i=1}^{N}\sum_{t=1}^{T} \rho_{\tau_k}(Y_{it} -b_{\tau_{k}}- \lambda_i'F_t),
\end{equation}
where $b_{\tau_{k}}$ estimates $b_{0\tau_{k}}$, the $100\tau_{k}\%$ quantile of $\varepsilon_{it}$. Let $\hat{\lambda}_i$ and $\hat{F}_t$ be the estimators for $\lambda_{0i}$ and $F_{0t}$ in \cref{eq: CQFM obj}. Unlike the solutions to \cref{eq: QFM obj}, $\hat{\lambda}_i$ and $\hat{F}_t$ are not dependent on any specific quantile position, and they estimate parameters in the mean factor model. By minimizing the objective function across multiple quantiles, the estimators can adapt to data features at different quantiles while still giving estimates for the mean of the process. We usually select equally spaced quantiles with $\tau_{k} = \frac{k}{K+1}$ for $k = 1,2,\cdots,K$ and $K$ is an odd number such as $5$ or $7$. This will always include the $50\%$ quantile in estimation, but an even number of quantiles also works for \cref{eq: CQFM obj}. The number $K$ can be viewed as a tuning parameter of CQFM. In the special case of $K=1$, i.e., when a single quantile position is used in \cref{eq: CQFM obj}, CQFM reduces to the quantile factor model. Our R package can estimate both CQFM and QFM.
There is no closed-form solution to the minimization exercise in \cref{eq: CQFM obj}, $\hat{b}_{\tau_{k}}$, $\hat{\lambda}_i$, and $\hat{F}_t$ need to be obtained through an iterative algorithm. Since $\lambda_i$ and $F_t$ appear as a product in \cref{eq: CQFM obj}, they are not separately identifiable. We use the following normalization for identification purposes:
\begin{equation} \label{eq: normalization}
\begin{aligned}
\frac{1}{T} \sum_{t=1}^{T} \hat{F}_t \hat{F}_t' &= I_r, \quad \text{an identity matrix of dimension $r$},\\
\frac{1}{N} \sum_{i=1}^{N} \hat{\lambda}_i \hat{\lambda}_i' &= \Sigma_{\hat{\lambda}}, \quad \text{a diagonal matrix with decreasing diagonal elements}.
\end{aligned}
\end{equation}
We describe the steps of the algorithm below. Let $s$ denote the iteration step and the ranges of the subscripts $i,t,k$ are the same as those appear in \cref{eq: CQFM obj}.
\begin{enumerate}[nosep]
\item[] \textit{Step 1}. Choose a random starting value for $F_t^{(0)}$ for all $t$. Use $F_t^{(0)}$ to get the initial estimates for $\lambda_{i}^{(0)}$ and $b_{\tau_{k}}^{(0)}$.
\item[] \textit{Step 2}. Given $b_{\tau_{k}}^{(s-1)}$ and $\lambda_i^{(s-1)}$ , obtain $F_t^{(s)} $ that minimizes \cref{eq: CQFM obj}.
\item[] \textit{Step 3}. Given $b_{\tau_{k}}^{(s-1)}$ and $F_t^{(s)} $, obtain $\lambda_i^{(s)}$ that minimizes \cref{eq: CQFM obj}.
\item[] \textit{Step 4}. Given $\lambda_i^{(s)}$ and $F_t^{(s)} $, obtain $b_{\tau_{k}}^{(s)}$ that minimizes \cref{eq: CQFM obj}.
\item[] \textit{Step 5}. Repeat Steps 2 to 4 for $s = 1, 2, \cdots$ until estimates converge. Normalize the final solution according to \cref{eq: normalization}.
\end{enumerate}
A few remarks follow.
\begin{remark}
The minimization exercise in \cref{eq: CQFM obj} is non-convex in the parameters. However, it is convex, for example, when we solve $\lambda_i^{(s)}$ given $b_{\tau_{k}}^{(s-1)}$ and $F_t^{(s-1)}$. Our simulation experience indicates the final solution is not sensitive to the starting values of $F_t^{(0)}$. Our R package \texttt{cqrfactor} allows the user to supply different random seeds to initialize $F_t^{(0)}$, making it easy to check the solution's sensitivity to starting values. This iterative strategy is also used in several other papers such as \cite{bai2009ecma,chenetal2021qfm}.
\end{remark}
\begin{remark}
When solving \cref{eq: CQFM obj} in steps 2 to 4, we use the majorization-minimization (MM) algorithm for quantile regression described in \cite{hunterandlange2000JCGSMM}. This is one of the several popular methods for solving a quantile regression problem.
\end{remark}
\begin{remark}
In step 5, there is no unique way to define the convergence of the algorithm. Between steps $s$ and $s+1$, one can check the difference in the loss function \cref{eq: CQFM obj} to see if it is small enough; alternatively, one can check the (average) absolute change in parameter estimates between steps $s$ and $s+1$.
\end{remark}
\subsection{Asymptotic results of the estimators}
We make the following assumptions to derive the asymptotic results.
\begin{assumption} \label{asump: factor and facor loading}
The factors $F_{0t}$ are random with $\frac{1}{T}\sum_{t=1}^{T}F_{0t} \rightarrow E(F_{0t}) = 0$ and $\frac{1}{T}\sum_{t=1}^{T}F_{0t} F_{0t}' \rightarrow \Sigma_{F_0} = I_r$ as $T \rightarrow \infty$. The factor loadings have the limit $\frac{1}{N} \sum_{i=1}^{N} \lambda_{0i}\lambda_{0i}' \rightarrow \Sigma_{\lambda_0} $, a diagonal matrix with $\sigma_{ii} > \sigma_{jj} > 0$ if $i < j$.
\end{assumption}
\begin{assumption} \label{asump: error density}
The distribution of the error term $\varepsilon_{it}$ has an absolutely continuous cumulative function $F_{\varepsilon}$ with a continuous density function $f_{\varepsilon}$ that is uniformly bounded away from $0$ and $\infty$.
\end{assumption}
\begin{assumption} \label{asump: iid}
The error terms $\varepsilon_{it}$ are i.i.d. and are independent of the factors $F_{0t}$ across all $i$ and $t$.
\end{assumption}
\Cref{asump: factor and facor loading} is almost identical to \citet[Assumption~F1]{stockandwatson2002JASA}
and \citet[Assumption~1(i)]{chenetal2021qfm} and can help identify both $F_0$ and $\Lambda_0$. See \cite{baing2013pcrfactor} for a more detailed discussion on the identification in factor models. \Cref{asump: error density} is a standard one in quantile regression. This assumption is made for the unconditional distribution of $\varepsilon_{it}$. If we consider the conditional distribution of $\varepsilon_{it}$ given $F_{0t}$ in \Cref{asump: error density}, all expectations in the proof will be conditional. The i.i.d. requirement in \Cref{asump: iid} is strong. However, this assumption simplifies the presentation of the asymptotic results and allows us to make direct comparison of our asymptotic results to the cross-section regression result in \cite{zouandyuan2008cqr}; in addition, the simple form of our asymptotic results facilitates the efficiency comparison between the CQFM-based factors and the PCA-based factors in \cite{bai2003ecma} (see a remark following \Cref{thm: asymptotic distribution} for a discussion). We can modify \Cref{asump: iid} so that it is conditional on $F_{0t}$, and the asymptotic covariance in \Cref{thm: asymptotic distribution} will have the standard sandwich form. We discuss this in a remark following \Cref{thm: asymptotic distribution}. Our simulation section includes results for errors with heteroskedasticity and AR(1) structure, and CQFM continues to give good results especially when the sample size is large.
The following theorem gives the asymptotic distribution of the estimated factors and factor loadings.
\begin{theorem} \label{thm: asymptotic distribution}
Under \cref{asump: factor and facor loading,asump: error density,asump: iid}, the asymptotic distribution of $\sqrt{N}(\hat{F}_t - F_{0t})$ is $N(0, \Sigma_{\text{CQFM,}F})$ with
\begin{equation*}
\Sigma_{\text{CQFM,}F} = \frac{\sum_{k_1 = 1}^{K}\sum_{k_2=1}^{K} \min(\tau_{k_1}, \tau_{k_2}) (1 -\max(\tau_{k_1}, \tau_{k_2})}{\Big(\sum_{k=1}^{K} f_{\varepsilon}(b_{0\tau_{k}}) \Big)^2} \Sigma_{\lambda_0}^{-1};
\end{equation*}
the asymptotic distribution of $\sqrt{T}(\hat{\lambda}_{i} - \lambda_{0i})$ is $N(0, \Sigma_{\text{CQFM,}\lambda})$ with
\begin{equation*}
\Sigma_{\text{CQFM,}\lambda} = \frac{\sum_{k_1 = 1}^{K}\sum_{k_2=1}^{K} \min(\tau_{k_1}, \tau_{k_2}) (1 -\max(\tau_{k_1}, \tau_{k_2})}{\Big(\sum_{k=1}^{K} f_{\varepsilon}(b_{0\tau_{k}}) \Big)^2} \Sigma_{F_0}^{-1}.
\end{equation*}
\end{theorem}
The format of the limiting distributions in \Cref{thm: asymptotic distribution} resembles the result for linear regression coefficients in \citet[Theorem~2.1]{zouandyuan2008cqr}.
\begin{remark}
When $K=1$, CQFM reduces to the quantile factor model. Results in \Cref{thm: asymptotic distribution} are comparable to those in \cite{andoandbai2020jasa}. Use the asymptotic distribution for factor loadings as an example. At quantile position $\tau$, its asymptotic variance is
\begin{equation} \label{eq:AndoBai lambda var}
\text{ \citet[Theorem~2]{andoandbai2020jasa}: }\tau(1-\tau)\Gamma_{i,0,\tau}^{-1}V_{i,0,\tau}\Gamma_{i,0,\tau}^{-1},
\end{equation}
where both $ V_{i,0,\tau}$ and $ \Gamma_{i,0,\tau}$ are defined in their theorem and ``$V_{i,0,\tau}$" is similar to $\Sigma_{F_0}$ in \Cref{thm: asymptotic distribution}. This sandwich estimator for covariance is commonly found in other papers on quantile regression with panel data such as \cite{katoandgalvao2012panelquantile,galvaoandkato2016joesmoothedquantile,chenetal2021qfm}. In \Cref{thm: asymptotic distribution} with $K=1$, based on \cref{eq:lambda diff equation,eq:lambda distribution}, we have
\begin{equation} \label{eq:sigma lambda k=1}
\Sigma_{\text{CQFM,}\lambda} = \tau(1-\tau) \Big(f_{\varepsilon}(b_{0\tau}) \Sigma_{F_0}\Big)^{-1} \Sigma_{F_0} \Big(f_{\varepsilon}(b_{0\tau}), \Sigma_{F_0}\Big)^{-1} = \frac{\tau(1-\tau)}{f_{\varepsilon}(b_{0\tau})^2 }\Sigma_{F_0}^{-1},
\end{equation}
which matches the result in \cite{andoandbai2020jasa}. Our result is made simpler by the i.i.d. errors in \Cref{asump: iid} that allow us to separate $f_{\varepsilon}(b_{0\tau_{k}})$ from $\Sigma_{F_0}$ in the term $f_{\varepsilon}(b_{0\tau_{k}}) \Sigma_{F_0}$; other papers typically consider the distribution of $\varepsilon_{it}$ conditional on either some regressors or the factors, see, for example, the term ``$\Gamma_{i,0,\tau} =T^{-1} \sum_{t=1}^{T}g_{it}(0|\cdot)z_{it,0,\tau}z_{it,0,\tau}$'' in \citet[Theorem~2]{andoandbai2020jasa}, where the conditional density function $g_{it}(0|\cdot)$ cannot be taken out of the summation sign as $T \rightarrow \infty$. This simplification can also be found in \citet[Theorem~4.1]{koenker_2005} for the linear quantile regression with i.i.d. errors. If the distribution of $\varepsilon_{it}$ is conditional on $F_{0t}$, $\Sigma_{\text{CQFA,}\lambda}$ will have a format similar to \cref{eq:AndoBai lambda var}.
\end{remark}
\begin{remark}
If the true factor and factor loading, $F_{0t}$ and $\lambda_{0i}$, do not meet the normalization conditions in \cref{eq: normalization}, $\hat{F}_t$ and $\hat{\lambda}_{i}$ estimate a rotation of the corresponding true values. Our proof can be adapted to incorporate a rotation matrix. To simplify the presentation of the asymptotic results, we assume factors and loadings are identifiable under the normalization assumptions and omit the rotation matrix in \Cref{thm: asymptotic distribution}, similar to \cite{andoandbai2020jasa}.
\end{remark}
\begin{remark}
Although the asymptotic results in \Cref{thm: asymptotic distribution} are developed for the panel mean factor model while those in \cite{andoandbai2020jasa,chenetal2021qfm} are for panel quantile factor model, all proofs are related to techniques in quantile regression. \cite{andoandbai2020jasa} give a proof based on the uniform consistency of parameter estimates and higher-order moment conditions on the error term; \cite{chenetal2021qfm} derive the asymptotic results based on a smoothed quantile objective function by replacing the indicator function with a differentiable kernel function. In our proof, we replace the objective function with an asymptotic quadratic form of the parameters and solve $\hat{\lambda}_{i} - \lambda_{i}$ and $\hat{F}_t - F_{0t}$ directly from the first-order conditions, similar to the proof strategy in \cite{zouandyuan2008cqr} for CQR and \cite{koenker_2005} for quantile regression.
\end{remark}
\begin{remark}
To compare the relative efficiency between CQFM and PCA-based solutions, we compute the asymptotic relative efficiency (ARE) of CQFM relative to PCA --- the ratio of their asymptotic variances. Consider the estimator for $F_0$. In CQFM, its variance is given in \Cref{thm: asymptotic distribution}; for PCA-based factor analysis, the variance is given in \citet[Theorem~1(i)]{bai2003ecma}. We will simplify the variance expression for $\hat{F}_t$ in \citet[Theorem~1(i)]{bai2003ecma} to facilitate the comparison. The notation for factor estimator is ``$\tilde{F}_t$" in \cite{bai2003ecma}, while we use $\hat{F}_t^{\text{PCA}}$ to denote the same estimator. An $r \times r$ rotation matrix, $H = (\Lambda_0'\Lambda_0/N)(F_0'\hat{F}^{\text{PCA}}/T)V_{NT}^{-1}$, is introduced in \citet[p.~158]{bai2003ecma} to describe the indeterminacy of the solutions, where $V_{NT}$ is a diagonal matrix that contains the eigenvalues of $(NT)^{-1}YY'$. For our purpose, it will be desirable to set $H = I_r$ so that $\sqrt{N}(\hat{F}_t^{\text{PCA}} - H'F_{0t})$ in \citet[Theorem~1(i)]{bai2003ecma} becomes $\sqrt{N}(\hat{F}_t^{\text{PCA}} - F_{0t})$, matching the format in \Cref{thm: asymptotic distribution}. Replacing $F_0$ in $H$ with the estimator $\hat{F}^{\text{PCA}}$ and using the normalization $ F^{\text{PCA}\prime} F^{\text{PCA}} / T= I_r$, we obtain $V_{NT} = \Lambda_0'\Lambda_0/N \rightarrow \Sigma_{\lambda}$, where the convergence result follows Assumption B in \cite{bai2003ecma}. This result, combined with equation (7) in \citet[p.~150]{bai2003ecma}, suggests the variance of $\sqrt{N}(\hat{F}_t^{\text{PCA}} - H'F_{0t})$ in \citet[Theorem~1(i)]{bai2003ecma} can be written as $\sigma_{\varepsilon}^2 \Sigma_{\lambda}^{-1}$ if $\varepsilon_{it}$ is i.i.d. and independent of $F_{0t}$, where $\sigma_{\varepsilon}^2$ is the variance of $\varepsilon_{it}$ and is assumed to be a finite number. This result greatly simplifies the efficiency comparison between CQFM and PCA-based factor analysis. Define the ARE of CQFM relative to the PCA-based factor analysis as
\begin{equation} \label{eq:are factor}
\text{ARE}(K)_F = \frac{\sigma_{\varepsilon}^2\Big(\sum_{k=1}^{K} f_{\varepsilon}(b_{0\tau_{k}}) \Big)^2}{\sum_{k_1 = 1}^{K}\sum_{k_2=1}^{K} \min(\tau_{k_1}, \tau_{k_2}) (1 -\max(\tau_{k_1}, \tau_{k_2})},
\end{equation}
which is identical to equation (3.1) in \cite{zouandyuan2008cqr}. As a result, we can apply \cite[Theorem~3.1]{zouandyuan2008cqr} to show that the ARE for the factor estimator from CQFM in \cref{eq:are factor} has a relative efficiency of at least $0.7026$ with respect to that of the PCA-based factor analysis when $K \rightarrow \infty$. This result suggests that, compared to PCA, CQFM factors will have about $30\%$ efficiency loss in the worst scenario. This is a conservative theoretical result. In our simulations (see \Cref{tab:mse asym error} and \Cref{tab:mse sym error} in the supplement), the mean squared error (MSE) of the estimated component $\hat{\lambda}_{i}' \hat{F}_t$ are mostly smaller or much smaller than that of PCA-based estimate. Efficiency loss, if any, is small based on our simulation study.
\end{remark}
\begin{remark}
To compute $\Sigma_{\text{CQFM,}F} $ and $\Sigma_{\text{CQFM,}\lambda} $, we first note that quantities such as $\tau_{k_1}$ and $\tau_{k_2}$, along with $K$, are chosen beforehand by the researcher. $\Sigma_{\lambda_0}$ can be replaced with the normalized diagonal matrix $ \hat{\Lambda}'\hat{\Lambda}/N$ while $\Sigma_{F_0}$ is $I_r$. There is no unique way to estimate the density $f_{\varepsilon}(b_{0\tau_{k}})$ (and its inverse). Since $f_{\varepsilon}(b_{0\tau_{k}})$ is the density of $\varepsilon$ at $100\tau_{k}\%$ quantile and the estimate for $\hat{b}_{\tau_{k}}$ is available, we can first obtain the residuals $\hat{\varepsilon}_{it}$, and use a consistent nonparametric density estimator for the residuals to obtain $\hat{f}_{\varepsilon}(\hat{b}_{\tau_{k}})$. Because of the i.i.d. error assumption, this simple estimate for $\Sigma_{\text{CQFM,}F} $ and $\Sigma_{\text{CQFM,}\lambda} $ is always positive definite.
\end{remark}
\begin{remark}
In a quantile factor model, the conditional quantile at $\tau_{k}$ can be written as
\begin{equation} \label{eq:qfm}
Q_{Y_{it}}(\tau_{k}) = \lambda_{0i}(\tau_{k})' F_{0t}(\tau_{k}) + b_{0\tau_{k}},
\end{equation}
where the factors and factor loadings vary with $\tau_{k}$. But our factor model in \cref{eq: factor model}, when used inside \cref{eq: CQFM obj}, have constant factors and factor loadings across selected quantiles. This is not a misspecification since our goal is to estimate the $\lambda_{0i}$ and $F_{0t}$ in the mean of \cref{eq: factor model} but not $ \lambda_{0i}(\tau_{k})$ and $F_{0t}(\tau_{k})$ in \cref{eq:qfm}. Much like in a standard linear regression, in addition to the least-squares method, one can use the lasso, principal components regression, LAD, Huber loss regression, CQR, \textit{etc.}, for estimation, there are several ways to estimate the mean factor model, and CQFM is one of the alternatives. While still permitting the quantile model in \cref{eq:qfm} for the data, CQFM combines the mean factor model in \cref{eq: factor model} with the composite quantile loss in \cref{eq: CQFM obj}. A single quantile loss in \cref{eq: QFM obj} gives estimates that adapt to data at a particular quantile. By using multiple quantiles, CQFM is designed to give the mean estimates that can adapt to data at multiple quantiles. Our simulation results demonstrate that this approach works well for several examples of data with asymmetry, heteroskedasticity and time series correlation.
\end{remark}
\section{Factor number selection} \label{sec:factor number}
The number of factors is assumed to be known in \Cref{thm: asymptotic distribution}. We discuss the selection of factor number in this section. Since the important work in \cite{baiandng2002ecma} on consistent factor number selection, there has been continued development of new methods in the literature. See \cite{baing2007jbesprimitive,amengualandwatson2007JBES,hallinandliska2007jasa,onatski2009ecma,lamandyao2012AOS} for panel mean regression models and \cite{andoandbai2020jasa,chenetal2021qfm} for panel quantile regression models.
Denote $r$ the estimated number of factors. To work with \cref{eq: CQFM obj}, we propose the following information criterion (IC):
\begin{align} \label{eq:IC}
IC(r) &= \log \Bigg[\frac{1}{NT} \sum_{k}^{K}\sum_{i=1}^{N}\sum_{t=1}^{T} \rho_{\tau_k}(Y_{it} -\hat{b}_{\tau_{k}}(r)- \hat{\lambda}_i(r)'\hat{F}_t(r))\Bigg] \nonumber \\
&\quad + r \times q(N,T),
\end{align}
where
\begin{equation} \label{eq:qnt}
q(N,T) = \left(\frac{N+T}{NT}\right) \log\left(\frac{NT}{N+T}\right),
\end{equation}
and we use $\hat{b}_{\tau_{k}}(r), \hat{\lambda}_i(r)$ and $\hat{F}_t(r)$ to denote estimates based on $r$ number of factors. This information criterion is similar to the one used in \cite{andoandbai2020jasa} and $\textit{IC}_{p1}$ in \citet[p.~201]{baiandng2002ecma}. \Cref{thm:factor number} shows the consistency of $IC(r)$. Let $C_{NT} = \min(N,T)$.
\begin{theorem} \label{thm:factor number}
Under \Cref{asump: factor and facor loading,asump: error density,asump: iid}, as $N,T \rightarrow \infty$, if $q(N,T) \rightarrow 0$,
the information criterion in \cref{eq:IC} selects the number of factors consistently.
\end{theorem}
See the online supplement for the proof.
\begin{remark}
The condition for $q(N,T)$ in \Cref{thm:factor number} defines a class of penalty functions, and \cref{eq:qnt} is an example of possibly many other penalty functions. The IC with \cref{eq:qnt} works quite well for most of the simulation examples in our study. However, it fails when the error term follows a \textit{t} distribution with $1$ degree of freedom ($t_1$). In this case, we propose another $q(N,T)$ function that works well with the $t_1$ distribution
\begin{equation} \label{eq:qnt 2}
q(N,T) = \log \left(\log\left(\frac{NT}{N+T}\right)\right) \left(\frac{N+T}{NT}\right).
\end{equation}
This penalty function also meets the requirement for $q(N,T)$ in \Cref{thm:factor number}, but it converges to $0$ faster and, consequently, imposes less penalty than \cref{eq:qnt} . Its performance for the $t_1$ error distribution is reported in \Cref{tab:factor number sym error} in the online supplement.
\end{remark}
\section{Monte Carlo simulation}
In this section, we use Monte Carlo simulation to study the finite sample properties of the CQFM method. To compare CQFM to other methods, we use the R code in \cite{heetal2022JBESnomoment} to compute the robust two-step (RTS) factors and the matlab code in \cite{chenetal2021qfm} to compute the QFM factors at quantile position $0.5$ (QFM(0.5)) and also the estimated factor numbers. When space permits, we also add the PCA results.
The number of quantiles in CQFM is an additional tuning parameter, and we choose $K = 5$ for demonstration purposes, which corresponds to the quantiles of $0.17, 0.33, 0.5, 0.67$ and $0.83$. A convergence criterion of $ 10^{-3}$ is used in the MM algorithm.
\subsection{Data simulation}
Consider the following 3-factor data generating process (DGP):
\begin{equation*}
Y_{it} = \sum_{j=1}^{3} \lambda_{0i,j} F_{0t,j} + \varepsilon_{it},
\end{equation*}
where $F_{0t,1} = 0.8 F_{0t-1,1} + e_{1t}$, $F_{0t,2} = 0.5 F_{0t-1,2} + e_{2t}$, $F_{0t,3} = 0.2 F_{0t-1,3} + e_{3t}$, and both $e$ and $\lambda_{0i,j}$ are i.i.d. $N(0,1)$. This is identical to the DGP in \citet[Section~5.1]{chenetal2021qfm} except that we consider several asymmetrical i.i.d. errors. They are summarized in \Cref{tab:error_asym_description}. Let $\gamma_1$ and $\gamma_2$ be the skewness and excess kurtosis coefficient, respectively.
\begin{table}[htp]
\centering
\caption{Description of the $5$ asymmetric error distributions}
\begin{tabular}{ll}
\toprule
\text{Error distribution}& parameter setting \\
\midrule
1. skewed normal (\textit{sn}) & $\mu_{\varepsilon} = 0, \sigma_{\varepsilon}=1, \gamma_1=0.99$ \\
2. skewed \textit{t} & $\mu_{\varepsilon} = 0, \sigma_{\varepsilon}=1, \gamma_1=0.99, \gamma_2=3$ \\
3. asymmetric Laplace& location$=0$, scale$=0.5$, asymmetry$=4$\\
4. log-normal & $\mu = 0, \sigma = 1.5$ \\
5. mixture of skewed normal& $0.9\cdot sn(0,1,0.99) + 0.1\cdot sn(0,9,0.99)$ \\
\bottomrule
\end{tabular}
\label{tab:error_asym_description}
\end{table}
The R package \texttt{sn} is used to simulate the skewed normal and skewed \textit{t} distributions in \Cref{tab:error_asym_description}. If one specifies the skewness parameter directly, the \texttt{sn} package restricts $|\gamma_1| < 0.99527$; we set $\gamma_1 = 0.99$ for the first two error distributions. The asymmetric Laplace error term is generated using the \texttt{rlaplace} function in the R package \texttt{LaplacesDemon}. The three numbers, $0,0.5,4$ correspond to the location, scale, and kappa parameter in the \texttt{rlaplace} function in R. For the asymmetric Laplace distribution, a kappa value of $4$ implies a skewness of about $-1.99$. Next, we consider a more skewed log-normal distribution with mean and s.d. equal to $0$ and $1.5$, respectively. These are the parameter values for the log-normal density, which implies the error term $\varepsilon_{it}$ has its mean, s.d., and skewness equal to $ 3.08, 8.97$ and $33.47$. For both the asymmetric Laplace distribution and the log-normal distribution, we subtract the theoretical mean from the simulated errors so that all error terms have zero mean. Finally, we consider a mixture of skewed normal distribution, where $sn(0,1,0.99)$ and $sn(0,9,0.99)$ denote the skewed normal distribution with $\mu_{\varepsilon} = 0, \sigma_{\varepsilon}=1, \gamma_1=0.99$ and $\mu_{\varepsilon} = 0, \sigma_{\varepsilon}=3, \gamma_1=0.99$, respectively. The weights for the mixture normal are $0.9$ and $0.1$.
We consider five different sample sizes: $\left(N,T\right) = (50,100), (100,50), (100,200), (200,100)$ and $(300,300)$. For each error distribution and sample size, we report the value of an evaluation metric based on $100$ replications for every estimation method.
\Cref{tab:adj R2 sym error with pca,tab:mse sym error,tab:factor number sym error} in the online supplement report the results for $5$ symmetric error distributions, including $N(0,1)$, \textit{t} distribution with $1$ degree of freedom, \textit{etc}. \Cref{tab:adj R2 asym error with pca heterodasticity,tab:mse asym error heteroskedasticity,tab:factor number asym error heterskedasticity,tab:adj R2 asym ar1 error with pca,tab:mse asym ar1 error,tab:factor number asym ar1 error} report the results for the following heteroskedasticity and AR(1) asymmetric errors:
\begin{align}
\text{heteroskedasticity } & Y_{it} = \sum_{j=1}^{3} \lambda_{0i,j} F_{0t,j} + \left[2+\cos(2\pi\times\lambda_{0i,4}F_{0t,4})\right]\times\varepsilon_{it}, \label{eq:heter}\\
\text{AR(1) error } & Y_{it} = \sum_{j=1}^{3} \lambda_{0i,j} F_{0t,j} + \varepsilon_{it} \text{ with } \varepsilon_{it} = 0.5\varepsilon_{i,t-1} + u_{it}, \label{eq:ar1}
\end{align}
where $\lambda_{0i,4}$ and $F_{0t,4}$ are i.i.d. $N(0,1)$ and $\varepsilon_{it}$ and $u_{it}$ are asymmetric errors defined in \Cref{tab:error_asym_description}. CQFM is found to have good finite properties in these cases too.
\subsection{Estimation of the factor and factor loading}
Similar to \citet[Table~1]{chenetal2021qfm}, \Cref{tab:adj R2 asym error} reports the average adjusted $R^2$ from regressing $F_{0t,1}, F_{0t,2}$, and $F_{0t,3}$ on the $3$ estimated factors from the RTS, QFM(0.5), and CQFM methods. Results in \Cref{tab:adj R2 asym error} assess how well the estimated factors span the space spanned by the true factors. \Cref{tab:adj R2 asym error with pca} in the supplement expands \Cref{tab:adj R2 asym error} to include the PCA results.
\begin{table}[htp] \centering
\begin{center}
\caption{Adj. $R^2$ of regressing $3$ true factors on the estimated factors}
\label{tab:adj R2 asym error}
\begin{threeparttable}
\begin{tabular}{lrrrrrrrrr}
\cmidrule(lr){1-10}
\multirow{2}{*}{(T,N)} & $R^2_{1,\text{RTS}}$ & $R^2_{2,\text{RTS}}$ & $R^2_{3,\text{RTS}}$ & $R^2_{1,\text{QFM}}$ & $R^2_{2,\text{QFM}}$ & $R^2_{3,\text{QFM}}$ & $R^2_{1,\text{CQFM}}$ & $R^2_{2,\text{CQFM}}$ & $R^2_{3,\text{CQFM}}$ \\
\cmidrule(lr){2-4} \cmidrule(lr){5-7} \cmidrule(lr){8-10}
& \multicolumn{9}{c}{$\varepsilon_{it} \sim$ skewed normal} \\
\cmidrule(lr){2-10}
(50,100) & 0.9950 & 0.9914 & 0.9893 & 0.9915 & 0.9857 & 0.9821 & 0.9951 & 0.9917 & 0.9898 \\
(100,50) & 0.9909 & 0.9828 & 0.9786 & 0.9847 & 0.9712 & 0.9645 & 0.9911 & 0.9832 & 0.9790 \\
(100,200) & 0.9978 & 0.9961 & 0.9949 & 0.9960 & 0.9928 & 0.9908 & 0.9979 & 0.9962 & 0.9952 \\
(200,100) & 0.9959 & 0.9921 & 0.9898 & 0.9926 & 0.9856 & 0.9816 & 0.9961 & 0.9924 & 0.9903 \\
(300,300) & 0.9986 & 0.9975 & 0.9967 & 0.9974 & 0.9953 & 0.9938 & 0.9987 & 0.9976 & 0.9969 \\
& \multicolumn{9}{c}{$\varepsilon_{it} \sim$ skewed t} \\
\cmidrule(lr){2-10}
(50,100) & 0.9950 & 0.9916 & 0.9895 & 0.9936 & 0.9895 & 0.9866 & 0.9954 & 0.9923 & 0.9903 \\
(100,50) & 0.9908 & 0.9828 & 0.9790 & 0.9885 & 0.9783 & 0.9737 & 0.9914 & 0.9840 & 0.9804 \\
(100,200) & 0.9978 & 0.9960 & 0.9949 & 0.9972 & 0.9950 & 0.9936 & 0.9981 & 0.9964 & 0.9954 \\
(200,100) & 0.9959 & 0.9921 & 0.9899 & 0.9947 & 0.9900 & 0.9871 & 0.9962 & 0.9928 & 0.9908 \\
(300,300) & 0.9986 & 0.9975 & 0.9967 & 0.9983 & 0.9968 & 0.9958 & 0.9988 & 0.9978 & 0.9971 \\
& \multicolumn{9}{c}{$\varepsilon_{it} \sim$ asymmetric Laplace} \\
\cmidrule(lr){2-10}
(50,100) & 0.9569 & 0.9254 & 0.9085 & 0.9482 & 0.9108 & 0.8682 & 0.9759 & 0.9590 & 0.9518 \\
(100,50) & 0.9277 & 0.8698 & 0.8438 & 0.9151 & 0.8443 & 0.8062 & 0.9573 & 0.9221 & 0.9057 \\
(100,200) & 0.9826 & 0.9673 & 0.9577 & 0.9777 & 0.9566 & 0.9341 & 0.9911 & 0.9834 & 0.9784 \\
(200,100) & 0.9664 & 0.9375 & 0.9204 & 0.9571 & 0.9143 & 0.8900 & 0.9823 & 0.9663 & 0.9572 \\
(300,300) & 0.9891 & 0.9797 & 0.9739 & 0.9843 & 0.9704 & 0.9613 & 0.9946 & 0.9900 & 0.9874 \\
& \multicolumn{9}{c}{$\varepsilon_{it} \sim$ log-normal} \\
\cmidrule(lr){2-10}
(50,100) & 0.5958 & 0.3578 & 0.2543 & 0.9541 & 0.8272 & 0.6834 & 0.9889 & 0.9823 & 0.9769 \\
(100,50) & 0.5637 & 0.3286 & 0.2245 & 0.9216 & 0.7739 & 0.5421 & 0.9781 & 0.9598 & 0.9512 \\
(100,200) & 0.8045 & 0.5902 & 0.4469 & 0.9754 & 0.8402 & 0.5698 & 0.9964 & 0.9932 & 0.9914 \\
(200,100) & 0.7470 & 0.5382 & 0.4320 & 0.9687 & 0.8033 & 0.4895 & 0.9931 & 0.9863 & 0.9826 \\
(300,300) & 0.8950 & 0.7967 & 0.7229 & 0.9881 & 0.8448 & 0.4473 & 0.9980 & 0.9963 & 0.9952 \\
& \multicolumn{9}{c}{$\varepsilon_{it} \sim$ mixture of skewed normal} \\
\cmidrule(lr){2-10}
(50,100) & 0.9908 & 0.9851 & 0.9807 & 0.9899 & 0.9835 & 0.9788 & 0.9937 & 0.9900 & 0.9870 \\
(100,50) & 0.9829 & 0.9694 & 0.9609 & 0.9816 & 0.9670 & 0.9586 & 0.9880 & 0.9785 & 0.9731 \\
(100,200) & 0.9960 & 0.9929 & 0.9908 & 0.9954 & 0.9920 & 0.9893 & 0.9974 & 0.9954 & 0.9941 \\
(200,100) & 0.9926 & 0.9858 & 0.9818 & 0.9914 & 0.9837 & 0.9789 & 0.9952 & 0.9908 & 0.9880 \\
(300,300) & 0.9976 & 0.9955 & 0.9941 & 0.9971 & 0.9946 & 0.9929 & 0.9985 & 0.9972 & 0.9964 \\
\hline
\end{tabular}
\begin{tablenotes}[flushleft]
\item[] \textit{Notes}: Each number is the average of adjusted $R^2$ over $100$ replications of regressing one of the three true factors on the estimated factors based on the RTS, QFM(0.5), and CQFM method, respectively. We choose $\tau=0.5$ for the QFM method and $K=5$ for the CQFM method.
\end{tablenotes}
\end{threeparttable}
\end{center}
\end{table}
For the first two error distributions with small skewness in \Cref{tab:adj R2 asym error}, all three methods perform well. Their differences in the adjusted $R^2$ mostly appear in the third digit. Still, we see CQFM performs slightly better than RTS and QFM. In the case of asymmetric Laplace error, the difference between CQFM and the other two methods start to grow larger. For example, for the sample size $(50,100)$, the adj. $R^2$ associated with $F_{0t,3}$ is $0.8682$ for QFM, while it is $0.9518$ for CQFM. The log-normal error distribution poses the greatest challenge to the other methods, as the regression yields much lower adj. $R^2$s compared to those of CQFM, and increasing sample size from $(50,100)$ to $(300,300)$ does not seem to help.
To further investigate the accuracy of the estimates, we compute the mean squared error (MSE) of the estimated components and report them in \Cref{tab:mse asym error}. The MSE is defined as
\begin{equation*}
\text{MSE} = \frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T} \big(\lambda_{0i}'F_{0t} - \hat{\lambda}_i'\hat{F}_{t} \big)^2.
\end{equation*}
\Cref{tab:mse asym error} clearly indicates that CQFM always yields the smallest MSE for the $5$ error distributions. Depending on the error distribution, CQFM's reduction in MSE can be huge -- its MSE can be a fraction of that of PCA.
We can draw several conclusions based on \Cref{tab:adj R2 asym error,tab:mse asym error}. First, compared to QFM at a single quantile position $\tau = 0.5$, the higher adj. $R^2$ and smaller MSE of CQFM suggests that there is some benefit in performing the estimation at multiple quantile positions simultaneously. Second, CQFM continues to work well in cases such as the asymmetric Laplace and log-normal errors, implying that CQFM can be a useful alternative to PCA in certain cases.
The online supplement also includes the adj. $R^2$ and MSE results (\Cref{tab:adj R2 sym error with pca,tab:mse sym error}) for $5$ symmetric error distributions. Overall, CQFM continues to provide robust estimates.
\begin{table}[ht] \centering
\begin{center}
\caption{MSE under asymmetric errors}
\label{tab:mse asym error}
\begin{threeparttable}
\begin{tabular}{lrrrrrrrr}
\cmidrule(lr){1-9}
\multirow{2}{*}{(T,N)} & \multicolumn{4}{c}{$\varepsilon_{it} \sim$ skewed normal} & \multicolumn{4}{c}{$\varepsilon_{it} \sim$ skewed t} \\
& RTS & QFM & CQFM & PCA & RTS & QFM & CQFM & PCA \\
\cmidrule(lr){2-5} \cmidrule(lr){6-9}
(50,100) & 0.101 & 0.157 & 0.097 & 0.090 & 0.099 & 0.114 & 0.090 & 0.089\\
(100,50) & 0.092 & 0.155 & 0.088 & 0.089 & 0.092 & 0.114 & 0.084 & 0.089\\
(100,200) & 0.048 & 0.085 & 0.048 & 0.045 & 0.048 & 0.058 & 0.044 & 0.045\\
(200,100) & 0.046 & 0.084 & 0.043 & 0.045 & 0.046 & 0.058 & 0.041 & 0.045\\
(300,300) & 0.021 & 0.040 & 0.020 & 0.020 & 0.021 & 0.026 & 0.019 & 0.020\\
& \multicolumn{4}{c}{$\varepsilon_{it} \sim$ asymmetric Laplace} & \multicolumn{4}{c}{$\varepsilon_{it} \sim$ log-normal} \\
\cmidrule(lr){2-5} \cmidrule(lr){6-9}
(50,100) & 0.836 & 1.139 & 0.458 & 0.801 & 19.888 & 4.124 & 0.214 & 36.456 \\
(100,50) & 0.794 & 1.134 & 0.439 & 0.799 & 14.344 & 4.157 & 0.211 & 36.752\\
(100,200) & 0.390 & 0.747 & 0.198 & 0.379 & 6.727 & 4.294 & 0.084 & 26.414\\
(200,100) & 0.382 & 0.744 & 0.196 & 0.381 & 4.865 & 4.283 & 0.076 & 26.531\\
(300,300) & 0.167 & 0.440 & 0.080 & 0.165 & 1.726 & 4.377 & 0.038 & 17.464\\
& \multicolumn{4}{c}{$\varepsilon_{it} \sim$ mixture of skewed normal} & \\
\cmidrule(lr){2-5}
(50,100) & 0.176 & 0.182 & 0.120 & 0.161 &&&& \\
(100,50) & 0.168 & 0.182 & 0.113 & 0.164 &&&&\\
(100,200) & 0.086 & 0.097 & 0.058 & 0.081 &&&& \\
(200,100) & 0.083 & 0.097 & 0.053 & 0.081 &&&&\\
(300,300) & 0.037 & 0.046 & 0.025 & 0.036 &&&&\\
\hline
\end{tabular}
\begin{tablenotes}[flushleft]
\item[] \textit{Notes}: Each number is the average MSE over $100$ replications for the RTS, QFM(0.5), CQFM, and PCA method, respectively. We choose $\tau=0.5$ for the QFM method and $K=5$ for the CQFM method.
\end{tablenotes}
\end{threeparttable}
\end{center}
\end{table}
\subsection{Estimation of the factor number}
Next, we study the performance of the information criterion in \cref{eq:IC,eq:qnt}. \Cref{tab:factor number} reports the average estimated factor number and the frequency of correct factor number estimation. Both CQFM and PCA perform well in majority of the cases. For log-normal error, all methods fail in small sample. However, CQFM yields better results when sample size is large with $\text{Prob}(\hat{r} = 3) = 82\%$ in the case of $(300,300)$.
\begin{table}[htp] \centering
\begin{center}
\caption{Average estimated factor number and frequency of correct estimation}
\label{tab:factor number}
\begin{threeparttable}
\makebox[\linewidth] {
\begin{tabular}{lrrrrrr}
\toprule
(T,N) & QFM & CQFM & PCA & QFM & CQFM & PCA \\
\hline
& \multicolumn{3}{c}{avg. $\hat{r}$} & \multicolumn{3}{c}{$\text{Prob}(\hat{r} = 3)$} \\
\cmidrule(lr){2-4} \cmidrule(lr){5-7}
& \multicolumn{6}{c}{$\varepsilon_{it} \sim$ skewed normal} \\
\cmidrule(lr){2-7}
(50,100) & 2.52 & 3 & 3 & 0.61 & 1 & 1 \\
(100,50) & 2.55 & 3 & 3 & 0.6 & 1 & 1 \\
(100,200) & 2.94 & 3 & 3 & 0.95 & 1 & 1 \\
(200,100) & 2.92 & 3 & 3 & 0.92 & 1 & 1 \\
(300,300) & 3 & 3 & 3 & 1 & 1 & 1 \\
& \multicolumn{6}{c}{$\varepsilon_{it} \sim$ skewed t} \\
\cmidrule(lr){2-7}
(50,100) & 2.52 & 3 & 3 & 0.59 & 1 & 1 \\
(100,50) & 2.55 & 3 & 3 & 0.61 & 1 & 1 \\
(100,200) & 2.94 & 3 & 3 & 0.95 & 1 & 1 \\
(200,100) & 2.92 & 3 & 3 & 0.92 & 1 & 1 \\
(300,300) & 3 & 3 & 3 & 1 & 1 & 1 \\
& \multicolumn{6}{c}{$\varepsilon_{it} \sim$ asymmetric Laplace} \\
\cmidrule(lr){2-7}
(50,100) & 2.66 & 1.8 & 2.9 & 0.66 & 0.14 & 0.9 \\
(100,50) & 2.75 & 1.74 & 2.91 & 0.71 & 0.12 & 0.91 \\
(100,200) & 3.31 & 3 & 3 & 0.64 & 1 & 1 \\
(200,100) & 3.37 & 3 & 3 & 0.57 & 1 & 1 \\
(300,300) & 3.95 & 3 & 3 & 0.05 & 1 & 1 \\
& \multicolumn{6}{c}{$\varepsilon_{it} \sim$ log-normal} \\
\cmidrule(lr){2-7}
(50,100) & 2.71 & 1.21 & 3.51 & 0.57 & 0.02 & 0.16 \\
(100,50) & 2.95 & 1.21 & 3.49 & 0.56 & 0.01 & 0.2 \\
(100,200) & 3.46 & 2.42 & 3.38 & 0.42 & 0.27 & 0.16 \\
(200,100) & 3.57 & 2.3 & 3.21 & 0.37 & 0.28 & 0.23 \\
(300,300) & 4 & 3.19 & 3.83 & 0 & 0.82 & 0.21 \\
& \multicolumn{6}{c}{$\varepsilon_{it} \sim$ mixture of skewed normal} \\
\cmidrule(lr){2-7}
(50,100) & 2.54 & 3 & 3 & 0.6 & 1 & 1 \\
(100,50) & 2.61 & 3 & 3 & 0.64 & 1 & 1 \\
(100,200) & 2.95 & 3 & 3 & 0.96 & 1 & 1 \\
(200,100) & 2.92 & 3 & 3 & 0.92 & 1 & 1 \\
(300,300) & 3 & 3 & 3 & 1 & 1 & 1 \\
\bottomrule
\end{tabular}
}
\begin{tablenotes}[flushleft]
\item[] \textit{Notes}: To estimate $r$, we use the rank minimization method in \cite{chenetal2021qfm} for QFM at $\tau = 0.5$, the IC in \cref{eq:IC,eq:qnt} for CQFM, and the $IC_{p1}$ in \cite[p.~201]{baiandng2002ecma} for the PCA method. Avg. $\hat{r}$ is based on $100$ replications.
\end{tablenotes}
\end{threeparttable}
\end{center}
\end{table}
For the $5$ symmetric error distributions in \Cref{tab:factor number sym error}, our proposed IC with \cref{eq:qnt} continues to work well except for the $t_1$ error. In this case, the rank-based approach in \cite{chenetal2021qfm} gives good results when sample size is large. After standardizing the data and using \cref{eq:qnt 2} in \cref{eq:IC}, CQFM also gives satisfactory results when the sample size is large.
\section{Empirical application}
In this section, we use the quarterly macroeconomic data set, FRED-QD, in \cite{mccrackenandng2020FREDQD} to study the properties of CQFM factors. We use the version ``2023-06.csv", which contains $258$ quarterly observations from 1959/3/1 to 2023/3/1 for $246$ macroeconomic variables. The data link is: \url{https://research.stlouisfed.org/econ/mccracken/fred-databases/}. We use the matlab code in \cite{mccrackenandng2020FREDQD} to prepare the data, including transforming all variables to stationary time series based on the \texttt{tcode} in \cite{mccrackenandng2020FREDQD}, removing outliers, and using the EM algorithm to fill in missing values. The final data set has $255$ quarterly observations and $246$ variables ($T = 255, N = 246$).
The number of estimated factors varies across different methods. For example, the CQFM estimate is $1$ ($3$ if \cref{eq:qnt 2} is used); the rank minimization method in \cite{chenetal2021qfm} reports $4$ factors at $\tau = 0.5$; the $IC_{p1}$ and $IC_{p2}$ in \citet[p.~201]{baiandng2002ecma} give $12$ and $8$ factors, respectively. Since our focus is on the property of the estimated factor, we follow \cite{stockandwaston2012nberrecession} and choose the number $6$ across different methods. The scree plot in \Cref{fig:screeplot} reveals why the proposed IC with \cref{eq:qnt} selects only one factor: the first eigenvalue is $54$ and explains about $22\%$ of the variation in the (standardized) data while the second eigenvalue is $19$ and explains about $7.8\%$ of the variance in the data.
\begin{figure}[th!]
\centering
\includegraphics[width=1.0\textwidth,keepaspectratio=TRUE]{fig_factors1-3.pdf}
\caption{The first three CQFM and PCA factors from 1959/3/1 to 2023/3/1 }
\label{fig:factors1-3}
\end{figure}
\Cref{fig:factors1-3} plots the first three CQFM and PCA factors (factors 4 to 6 are plotted in \Cref{fig:factors4-6}). The first three factors from CQFM and PCA are very similar to each other. Despite this visual similarity, the estimated factors exhibit different moment properties. \Cref{tab:factor moments} summarizes the skewness and kurtosis of the six estimated factors from the four methods. We make a few observations. First, the CQFM-based factors tend to have larger skewness and kurtosis in the first few factors. This means, if the data have large skewness and/or kurtosis, the CQFM-based factors will likely give a better fit for the component ($\lambda_{0i}' F_{0t}$). Second, even if other methods such as PCA-based factors exhibit larger skewness and/or kurtosis in later factors -- for example, the 5th PCA factor exhibits larger skewness than CQFM, these larger value will unlikely be helpful in capturing the skewness and kurtosis in the data since it is typically the first few factors that determines the overall variability of the data. Third, compared to CQFM-based factors, the QFM-based factors exhibit less skewness and kurtosis, suggesting that estimation done at a single quantile position such as $\tau = 0.5$ may not be effective in capturing certain features of the data; the composite quantile approach is more effective in this regard. The small MSEs for CQFM in \Cref{tab:mse asym error} attest to the above arguments.
\begin{table}[t] \centering
\begin{center}
\caption{Skewness and kurtosis of the $6$ estimated factors}
\label{tab:factor moments}
\begin{threeparttable}
\makebox[\linewidth] {
\begin{tabular}{lrrrrrr}
\toprule
method & $\hat{F}_1$ & $\hat{F}_2$ & $\hat{F}_3$ & $\hat{F}_4$ & $\hat{F}_5$& $\hat{F}_6$ \\
\hline
& \multicolumn{6}{c}{skewness ($\gamma_1$)}\\
\cmidrule(lr){2-7}
\text{RTS} & 1.85 & -0.33 & -0.38 & -0.27 & -0.52 & 0.17 \\
\text{QFM} & -1.95 & 0.70 & -0.33 & 0.34 & 0.35 & 0.43 \\
\text{CQFM} & -2.51 & -0.93 & 0.59 & -0.52 & -0.02 & -0.16 \\
\text{PCA} & -2.23 & -0.76 & 0.51 & -0.18 & -0.12 & 0.00 \\
& \multicolumn{6}{c}{kurtosis ($\gamma_2$)}\\
\cmidrule(lr){2-7}
\text{RTS} & 14.55 & 4.43 & 13.57 & 2.76 & 3.68 & 3.47 \\
\text{QFM} & 17.62 & 6.54 & 5.25 & 2.80 & 3.08 & 7.50 \\
\text{CQFM} & 26.03 & 7.98 & 4.93 & 15.58 & 2.71 & 3.32 \\
\text{PCA} & 24.36 & 6.62 & 4.51 & 16.70 & 2.88 & 4.87 \\
\bottomrule
\end{tabular}
}
\begin{tablenotes}[flushleft]
\item[] \textit{Notes}: This table reports the skewness and kurtosis of the estimated six factors for different methods based on the FRED-QD data between 1959Q1 and 2023Q1. We focus on the magnitude of $\gamma_1$ since factors have sign indeterminacy.
\end{tablenotes}
\end{threeparttable}
\end{center}
\end{table}
Next, following \cite{stockandwatson2002jbes}, we use the diffusion indexes to forecast one-quarter-ahead macroeconomic variables. The forecasting function is given by
\begin{equation} \label{eq:diffusion_indexes_forecast}
y_{i,t+1} = \beta_i + \sum_{j=0}^{3} \beta_j y_{i,t - j} + \beta_F' \hat{F}_{t} + \epsilon_{i,t+1}, \text{ for } i = 1, \cdots, 246,
\end{equation}
where $y_{it}$ is the original data transformed according to the \texttt{tcode} in FRED-QD ``2023-06.csv" and $\hat{F}_t$ is the estimated $6$ factors at time $t$. The forecast period starts from $2000 Q1$ to $2023Q1$, a total of $93$ forecasts for each of the $246$ macroeconomic variables, and this forecast period covers three NBER-determined recessions, including the one induced by the recent pandemic. For each rolling forecast, we use a rolling window of $120$ quarters to estimate factors and the coefficients $\beta_j$ and $\beta_F$. We forecast the data that are transformed using the \texttt{tcode} in \cite{mccrackenandng2020FREDQD} and convert the forecast back to data in their original levels.
\begin{table}[t] \centering
\begin{center}
\caption{Forecast RMSE of GDP, Unemployment rate, and Inflation}
\label{tab:forecast_rmse}
\begin{threeparttable}
\makebox[\linewidth] {
\begin{tabular}{lrrrrr}
\toprule
variable & \multicolumn{1}{l}{RTS} & \multicolumn{1}{l}{QFM} & \multicolumn{1}{l}{CQFM} & \multicolumn{1}{l}{PCA} & \multicolumn{1}{l}{AR(4)} \\
\midrule
\texttt{GDPC1} & 322.131 & 303.498 & \textbf{258.293} & 355.526 & 299.965 \\
\texttt{UNRATE} & 1.227 & 1.110 & \textbf{1.093} & 1.246 & 1.203 \\
\texttt{CPIAUCSL} & \textbf{3.616} & 4.173 & 3.832 & 3.798 & 3.739 \\
\texttt{avg RMSE} & 8269.7 & 8302.2 & \textbf{8053.6} & 8578.6 & 8219.2 \\
\bottomrule
\end{tabular}
}
\begin{tablenotes}[flushleft]
\item[] \textit{Notes}: This table reports the average forecast RMSE over $93$ forecasts from $2000Q1$ to $2023Q1$. \texttt{GDPC1} is the real GDP in chained 2012 dollars; \texttt{UNRATE} is the civilian unemployment rate (percent); \texttt{CPIAUCSL} is the CPI for all urban consumers. Results for columns 1 to 4 are based on an AR(4) model with six factors as additional regressors. The last column reports the forecast RMSE of the AR(4) model with no augmented factors. The variable \texttt{avg RMSE} reports the average of RMSE for all the $246$ macroeconomic time series for each of the $5$ methods.
\end{tablenotes}
\end{threeparttable}
\end{center}
\end{table}
\Cref{tab:forecast_rmse} reports the forecast root-MSE (RMSE) for the three most common macroeconomic variables, real gross domestic product (\texttt{GDPC1}), civilian unemployment rate (\texttt{UNRATE}), and consumer price index for all urban consumers (\texttt{CPIAUCSL}) in the FRED-QD data set. CQFM gives good results, but its performance is not the best for the CPI data. It's also somewhat surprising that the AR(4) model can sometimes do better than factor-augmented methods. In the last row, we compute the average of RMSE over the $246$ macroeconomic variables for each of the $5$ models, and CQFM gives the smallest average RMSE. Notice that the results for the $4$ factor-based models are obtained by simply choosing $6$ factors without any additional tuning of the model. \cite{mccrackenandng2020FREDQD} consider $7$ factors, and, for the regression in \cref{eq:diffusion_indexes_forecast}, they try $2^7-1 = 127$ different combinations of the $7$ factors. Similar approach can also be used here to possibly improve the performance of the factor-based models. In addition, many other aspects of diffusion index modeling can be tuned to yield a favorable model, which includes, but not limited to, the number lags of the factors (we consider only $1$ in our regression), the forecast horizon (3-month, 6-month, one-year, \textit{etc}.), the inclusion of lag variables in the $Y$ matrix in \cref{eq: factor model matrix form} in factor analysis, the use of balanced panel data vs. unbalance panel data with EM-algorithm-generated data, the types of data transformation used, whether to split the data before and after a recession, among others. In the case of CQFM, we can also tune the parameter $K$ to possibly improve its performance. \Cref{tab:forecast_rmse} is a simple demonstration of the use of CQFM-based factors. A more comprehensive study is needed to further study the properties of different factor-based models.
\section{Conclusions}
In this paper, we develop the method of composite quantile factor model. We demonstrate in both simulations and an empirical study that, compared to PCA and several other methods, CQFM can be more effective in modeling asymmetric data due to its capability of adapting to data at multiple quantile positions. Asymptotic distributional theory and an information criterion for consistent factor number selection are also discussed. PCA-based method is popular for factor analysis, and CQFM will be a useful addition to a researcher's toolkit when handling non-normal data.
Many extensions of the current research are possible, and we give two examples that are highly relevant to data modeling. One is the creation of sparsity in CQFM. Adding penalty functions to \cref{eq: CQFM obj} gives
\begin{equation} \label{eq: CQFM obj with penalty}
\frac{1}{NT} \sum_{k}^{K}\sum_{i=1}^{N}\sum_{t=1}^{T} \rho_{\tau_k}(Y_{it} -b_{\tau_{k}}- \lambda_i'F_t) + \text{penalty}(F) + \text{penalty}(\Lambda).
\end{equation}
\cite{zouandyuan2008cqr} use the adaptive lasso in \cite{zou2006jasaadplasso} to induce sparsity in linear regression, and many other penalty functions are available for $F$ and $\Lambda$. The other example, following the work in \cite{bai2009ecma}, is to add a regression component to \cref{eq: CQFM obj} so that it becomes the panel data model with interactive fixed effects
\begin{equation} \label{eq: CQFM obj with regression}
\frac{1}{NT} \sum_{k}^{K}\sum_{i=1}^{N}\sum_{t=1}^{T} \rho_{\tau_k}(Y_{it} -b_{\tau_{k}}- X_{it}'\beta - \lambda_i'F_t),
\end{equation}
where $X_{it}$ is a vector of regressors. The model in \cref{eq: CQFM obj with regression} is a hybrid of CQR in \cite{zouandyuan2008cqr} and CQFM in the current paper. Given the good finite sample properties of CQR and CQFM under certain non-normal data, we expect estimators from \cref{eq: CQFM obj with regression} will also show some robustness to non-normal data. We leave these topics for future research.
\spacing{1.45}
\bibliographystyle{ecca}
\bibliography{reference}
\newpage
\setcounter{page}{1}
\spacing{1.42}