Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
95,837 characters · 19 sections · 82 citation commands
Quantile Factor Models
\affil[1]{School of Economics, Shanghai University of Finance and Economics, [email removed] } \affil[2]{Department of Economics, Universidad Carlos III de Madrid, [email removed]} \affil[3]{Department of Economics, Universidad Carlos III de Madrid, [email removed]}
\thispagestyle{empty}
\affil[1]{School of Economics, Shanghai University of Finance and Economics, [email removed] } \affil[2]{Department of Economics, Universidad Carlos III de Madrid, [email removed]} \affil[3]{Department of Economics, Universidad Carlos III de Madrid, [email removed]}
{\ }\setcounter{page}{1} Following the key contributions by ross1976arbitrage, chamberlain1983arbitrage and connor1986performance to the theory of approximate factor models (AFM henceforth) in the context of asset pricing, the analysis and applications of this class of models have proliferated thereafter. As it is well known, AFM imply that a panel $\{X_{it}\}$ of $N$ variables (units), each with $T$ observations, has the representation $X_{it}=\lambda _{i}^{\prime }f_{t}+\epsilon _{it}$, where $\lambda _{i}=[\lambda _{i1},..,\lambda _{ir}]^{\prime }$ and $f_{t}$ $=[f_{t1},..,f_{tr}]^{ \prime }$ are $r\times 1$ vectors of factor loadings and common factors, respectively, with $r\ll N$, and $\{\epsilon _{it}\}$ are zero-mean weakly dependent idiosyncratic disturbances which are uncorrelated with the factors.
The fact that it is easy to construct theories involving common factors, at least in a narrative version, together with the availability of fairly straightforward estimation procedures for AFM--- e.g. via Principal Components Analysis (PCA), has led to their extensive use in many fields of economics.\footnote{ Early applications of AFM abound in Aggregation Theory, Consumer Theory, Business Cycle Analysis, Finance, Monetary Economics, and Monitoring and Forecasting; see,\ inter alia, bai2003inferential, bai2008large, stock2011dynamic. More recently, the characterization of cross-sectional dependence among error terms in Panel Data has relied on the use of a finite number of unobserved common factors which originate from economy-wide shocks affecting units with different intensities (loadings). Interactive fixed-effects models can be easily estimated by PCA (see bai2009panel) or by common correlated effects (see pesaran2006estimation), and there are even generalizations of these techniques for nonlinear panel single-index models (see chen2018panel). Lastly, the surge of Big Data and Machine Learning technologies has made factor models a key tool for dimension reduction and predictive analytics when using very large datasets (see athey2019machine for a survey).}
Inspired by the generalization of linear regression to quantile regression (QR) models, our starting point in this paper is to claim that the standard regression interpretation of static AFM as linear conditional mean models of $X_{it}$ given $ f_{t}$ (i.e. $\mathbb{E}(X_{it}|f_{t})=\lambda _{i}^{\prime }f_{t}$), entails two possibly restrictive features. On the one hand, PCA does not capture hidden factors that may shift characteristics (moments or quantiles) of the distribution of $X_{it}$ other than its mean. On the other hand, neither the loadings $\lambda _{i}$ nor the factors $f_{t}$ are allowed to vary across the distributional characteristics of each unit in the panel.
Highlighting these limitations, a growing literature in empirical finance has been documenting a much more pronounced co-movement of financial asset returns at the lower part than at the rest of their distributions. In particular, doi:10.1080/07350015.2018.1505631 reject the null hypothesis of a Gaussian copula when analyzing the cross-sectional dependence among monthly returns on individual US stocks, for which they find nonlinear tail dependence, co-skeweness and co-kurtosis. Likewise, ando2018quantile (\textcolor{blue}{AB 2020}, hereafter) report that the common factor structures explaining the asset return distributions in global financial markets since the subprime crisis are different in the lower and the upper tails. In international finance, maravalle2018changes find that, while global and regional factors were key in explaining fluctuations of European government bond yields between the Great Moderation and the onset of the Great Recession (lower tails of the distributions of bond yields), their role decreased afterwards, giving way to country-specific factors as the main driving forces.\footnote{The diminishing role of the former factors could result from quantitative easing and the increase in yields (upper tail of the distribution), while the emergence of the latter factors could be due to the financial fragmentation in a number of vulnerable euro-area countries during sovereign debt crisis.} On the macro side, adrian2019vulnerable find that, while the estimated lower conditional quantiles of the distribution of future GDP growth in the US exhibit strong dependence on current financial conditions, the upper quantiles are stable over time. Lastly, in micro theory, de2019dynamic have recently extended the expected utility model of rational behavior to quantile utility preferences, where e.g. factor structures determining hedonic pricing of consumption goods or financial stocks may exhibit large differences across quantiles.
A simple way of illustrating the above limitations of the conventional formulation of AFM\ is to consider the factor structure in a location-scale shift model with the following Data Generating Process (DGP): $ X_{it}=\alpha _{i}f_{1t}+\eta_{i}f_{2t}\epsilon _{it}$, with $f_{1t}\neq f_{2t}$ (both are scalars), $\eta_i,f_{2t}>0$ and $\mathbb{E}(\epsilon _{it})=0$. The first factor ($f_{1t}$) shifts location, whereas the second factor ($f_{2t}$) shifts the scale and therefore governs the volatility of shocks to $X_{it}$. This model has been proposed by herskovic2016common to empirically document a strong co-movement of the volatilities in the idiosyncratic component of individual stock returns and firm-level cash flows.\footnote{ This DGP is further discussed in subsection 2.2 below, where we present a larger set of illustrative models as examples of potential DGPs for $X_{it}$. Notice that the simplifying assumption of a known number of factors in this specific example is later relaxed.} Such a DGP can be rewritten in QR format as $X_{it}=\lambda _{i}^{\prime }(\tau )f_{t}+u_{it}(\tau )$, with $0<\tau <1 $, $\lambda _{i}(\tau )=[\alpha _{i},\eta_{i}\mathsf{Q} _{\epsilon}(\tau )]^{\prime }$, where $\mathsf{Q} _{\epsilon}(\tau )$ represents the quantile function of $\epsilon _{it}$, $f_{t}=[f_{1t},f_{2t}]^{\prime }$, $u_{it}(\tau )=\eta_{i}f_{2t}[\epsilon _{it}-\mathsf{Q} _{\epsilon}(\tau )]$, and the conditional quantile $\mathsf{Q}_{u_{it}(\tau)}[\tau|f_{t}]=0$. \footnote{ Throughout the paper we use $\mathsf{Q}_{W}[\tau |Z]$ to denote the conditional quantile of $W$ given $Z$.} PCA will only extract the location-shifting factor $ f_{1t}$ in this model, but it will fail to capture the scale-shifting factor $f_{2t}$ and the quantile-dependent loadings $\lambda _{i}(\tau )$ in its QR representation. As will be explained below, our estimation procedure allows to estimate the space spanned by $f_{1t}$ and $f_{2t}$.\footnote{Given that $f_{1t}$ can be consistently estimated by PCA, it is also feasible to separate $f_{2t}$ from their joint space.} Also notice that, when the distribution of $\epsilon _{it}$ is symmetric, then $f_{t}$ can be considered as being quantile dependent, i.e. $f_{t}(\tau)$, since $f_{t}(\tau )=f_{1t}$ for $\tau =0.5$, and $f_{t}(\tau )=[f_{1t},f_{2t}]^{\prime }$ for $\tau \neq 0.5$. Together with the remaining examples discussed below, this means that the general class of models to be considered in the sequel would be one where the loadings, factors and the number of factors are all allowed to be quantile-dependent objects, namely, $\lambda _{i}(\tau )$, $f_{t}(\tau )$ and $r(\tau)$ for $\tau \in (0,1)$. In what follows, we coin this class of models Quantile Factor Models (QFM, hereafter), whose detailed definition is provided in Section 2 below.
That said, our goal in this paper is to develop a common factor methodology for QFM which is flexible enough to capture the quantile-dependent objects that standard AFM tools fail to recover. To do so, we analyze their estimation and inference, including selection criteria for the number of factors at each quantile $\tau $. Put succinctly, QFM could be thought of as capturing the same type of flexible generalization that QR techniques represent for linear regression models.
To help understand how this new methodology works, we start by proposing an estimation approach for the quantile-dependent objects in QFM, labeled Quantile Factor Analysis (QFA, henceforth). Our QFA estimation procedure relies on the minimization of the standard check function in QR (instead of the conventional quadratic loss function used in AFM) to estimate jointly the common factors $f_{t}(\tau )$ and the loadings $\lambda _{i}(\tau )$ at a given quantile $\tau $, once the number of factors has been selected. However, since the objective function for QFM is not convex in the relevant parameters, we introduce an iterative QR algorithm which yields estimators of the quantile-dependent objects. We then derive their average rates of convergence, and propose two consistent selection criteria (one based on information criteria and another on rank minimization) for the number of factors at each $\tau$. In addition, we establish asymptotic normality for QFA estimators based on smoothed QR (see e.g., horowitz1998bootstrap and galvao2016smoothed). Moreover, given that QFA estimation captures all quantile-shifting factors (including those affecting the means of observed variables), our asymptotic results and the proposed selection criteria provide a natural way to differentiate AFM from QFM.
In sum, the key contributions of this paper to the literature on factor models can be summarized as follows:
\noindentRelated literature
There is a recent literature that attempts to make the AFM setup more flexible. For example, su2017time allow for the factor loadings to be time-varying and pelger2018interpretable admit these loadings to be state dependent. chen2009nonlinear provide a theory for nonlinear PCA, where they favor sieve estimation to retrieve nonlinear factors. Finally, gorodnichenko2017level propose an algorithm to estimate level and volatility factors simultaneously. Different from these studies, our approach to modelling nonlinearities in factor models is through the conditional quantiles of the observed data.
On top of this, there is an emerging literature on heterogeneous panel quantile models with factor structures, especially in financial economics. The main idea is that a few unobservable factors explain co-movements of asset return distributions in a large range of asset returns observed at high frequencies, as in stock markets. In parallel and independent research, there have been two papers related to ours. First, ma2017estimation propose estimation and inference procedures in semiparametric quantile factor models. In these models, factor loadings/betas are smooth functions of a small number of observables under the assumption that the included factors all have non zero mean. Then, sieve techniques are used to obtain preliminary estimation of these functions for each time period. Finally, the factor structure is imposed in a sequential fashion to estimate the factor returns by GLS under weak conditions on cross-sectional and temporal dependence. We depart from these authors in that we do not need to assume the loadings to depend on observables and, foremost, in that not only loadings but also factors are quantile-dependent objects in our setup. Second, in a closely related paper, \textcolor{blue}{AB (2020)} use a similar setup to ours, where the unobservable factor structure is also allowed to be quantile dependent. These authors use Bayesian MCMC and frequentist estimation approaches, the latter building upon our proposed iterative procedure, as duly acknowledged in their paper. However, we differ from \textcolor{blue}{AB (2020)} in several respects which make our QFA approach valuable: (i) our assumptions are less restrictive, since we rely on properties of the density, as in QR, while \textcolor{blue}{AB (2020)} needs all the moments of the idiosyncratic errors to exist, (ii) our proofs of the main results are different from theirs, and (iii) our rank-minimization selection criterion to estimate the number of factors is novel, behaves well in finite samples and is computationally more efficient than the information criteria-based method.
Finally, it is noteworthy that the illustrative location-scale shift model above, where $f_{1t}$ $\neq f_{2t}$, is behind a current line of research in asset pricing which has been coined the \textquotedblleft idiosyncratic volatility puzzle\textquotedblright\ by ang2006cross. This approach focuses on the co-movements in the idiosyncratic volatilities of a panel of asset returns, and consists of applying PCA to (or taking cross-sectional averages of) the squared residuals, once the mean (PCA) factors have been removed from the original variables (a procedure labeled PCA-SQ hereafter).\footnote{ See, e.g., barigozzi2016generalized, herskovic2016common and renault2016apt. Notice that the volatility co-movement does not arise from omitted factors in the AFM but from assuming a genuine\ factor structure in the idiosyncratic volatility processes.} For example, this technique would be valid for our illustrative example above. Yet, while the QFA approach is able to recover the whole QFM structure for more general DGPs than the previous model (see subsection 2.2), PCA-SQ fails to do so. It will also fail when the idiosyncratic errors do not have bounded eighth moments. Hence, to the best of our knowledge, our QFA approach becomes the first estimation procedure capable of dealing with these issues.
\noindentStructure of the Paper
The rest of the paper is organized as follows. Section 2 defines QFM and provides a list of simple illustrative examples where the new QFM methodology applies. In Section 3, we present the QFA estimator and its computational algorithm, establish the average rates of convergence of the quantile-dependent factors and factor loadings, and propose two consistent selection criteria to choose the number of factors at each quantile. Section 4 introduces a kernel-smoothed version of the QFA estimators to derive their asymptotic distributions. Section 5 contains some Monte Carlo simulation results to evaluate the performance in finite samples of our estimation procedures relative to other alternative approaches with different assumptions about the idiosyncratic error terms. Section 6 considers three empirical applications using three large panel datasets, where we document the relevance of extra factors in causal analysis, forecasting and economic interpretation of common factors. Finally, Section 7 concludes and suggests several avenues for further research. Proofs of the main results are collected in the Online Appendix.
\noindentNotations
The Frobenius norm is denoted as $\Vert \cdot \Vert $. For a matrix A with real eigenvalues, $\rho _{j}(A)$ denotes the $j$th largest eigenvalue. Following vanweak, the symbol $\lesssim $ means \textquotedblleft left side bounded by a positive constant times the right side\textquotedblright\ (the symbol $\gtrsim $ is defined similarly), and $ D(\cdot ,g,\mathcal{G})$ denotes the packing number of space $\mathcal{G}$ endowed with semimetric $g$.
This section starts by introducing the main definitions to be used throughout the paper. Next, we show how to derive the QFM representation of several illustrative DGPs exhibiting different factor structures.
Suppose that the observed variable $X_{it},$ with $i=1,2,..,N$ and $ t=1,2,...,T$, has the following QFM structure at some $\tau \in (0,1)$:
where the common factors $f_{t}(\tau )$ is a $r(\tau )\times 1$ vector of unobservable random variables, $\lambda _{i}(\tau )$ is a $r(\tau )\times 1$ vector of non-random factor loadings with $r(\tau)\ll N$. Note that in the QFM defined above, the factors, the loadings, and the number of factors are all allowed to be quantile-dependent.
Alternatively, the above equation implies that
where the quantile-dependent idiosyncratic error $u_{it}(\tau )$ is assumed to satisfy the following quantile restrictions:
In this section we provide a few illustrative examples of how QFMs can be derived from different specifications of location-scale shift models and related ones. The goal of these simple illustrations is to show instances where the standard AFM methodology may fail to capture the full factor structure, therefore requiring the use of the alternative QFM approach.
Not surprisingly, the standard AFM methodology based on PCA only works in Example 1, insofar as the idiosyncratic errors satisfy certain moment conditions. In the remaining examples, PCA will only yield consistent estimates of those factors shifting the locations; however (except in Example 2), it will fail to capture those extra factors which shift quantiles other than the means, or their corresponding quantile-varying loadings. In the sequel, QFA is therefore proposed as a new estimation procedure to estimate both sets of quantile-dependent objects in QFM.
To simplify the notations, we suppress hereafter the dependence of $ f_{t}(\tau ),\lambda _{i}(\tau ),r(\tau )$ and $u_{it}(\tau )$ on $\tau $, so that the QFM in (ref) is rewritten as:
where $\lambda _{i},f_{t}\in \mathbb{R}^{r}$. Suppose that we have a sample of observations $\{X_{it}\}$ generated by (ref) for $i=1,\ldots ,N,$ and $t=1,\ldots ,T$, where the realized values of $\{f_{t}\}$ are $ \{f_{0t}\} $ and the true values of $\{\lambda _{i}\}$ are $\{\lambda _{0i}\} $. We take a fixed-effects approach by treating $\{\lambda _{0i}\}$ and $ \{f_{0t}\}$ as parameters to be estimated, and our asymptotic analysis is conditional on $\{f_{0t}\}$. In Section 3.1, we consider the estimation of $\{\lambda _{0i}\}$ and $\{f_{0t}\}$ while $r$ is assumed to be known. Finally, Section 3.2 deals with the estimation of $r$ for each quantile.
It is well known in the literature on factor models that $\{\lambda _{0i}\}$ and $\{f_{0t}\}$ cannot be separately identified without imposing normalizations (see \textcolor{blue}{BN 2002}). Without loss of generality, we choose the following normalizations:
Let $M=(N+T)r$, $\theta =(\lambda _{1}^{\prime },\ldots ,\lambda _{N}^{\prime },f_{1}^{\prime },\ldots ,f_{T}^{\prime })^{\prime }$, and $ \theta _{0}=(\lambda _{01}^{\prime },\ldots ,\lambda _{0N}^{\prime },f_{01}^{\prime },\ldots ,f_{0T}^{\prime })^{\prime }$ denotes the vector of true parameters, where we also suppress the dependence of $\theta $ and $ \theta _{0}$ on $M$ to save notation. Let $\mathcal{A},\mathcal{F} \subset \mathbb{R}^{r}$ and define:
Further, define:
where $\rho _{\tau }(u)=(\tau -\mathbf{1}\{u\leq 0\})u$ is the check function. The QFA estimator of $\theta _{0}$ is defined as:
It is obvious that the way in which our estimator is related to the PCA estimator studied by \textcolor{blue}{BN (2002)} and bai2003inferential is analogous to how QR is related to standard least-squares regressions. However, unlike bai2003inferential's PCA estimator, our estimator $\hat{\theta}$ does not yield an analytical closed form. This makes it difficult not only to find a computational algorithm that would yield the estimator, but also the analysis of its asymptotic properties. In the sequel, we introduce a computational algorithm called iterative quantile regression (IQR, hereafter) that can effectively find the stationary points of the object function. In parallel, Theorem 1 shows that $\hat{\theta}$ achieves the same convergence rate as the PCA estimators for AFM.
To describe the algorithm, let $\Lambda =(\lambda _{1},\ldots ,\lambda _{N})^{\prime }$, $F=(f_{1},\ldots ,f_{T})^{\prime }$, and define the following averages:
Note that we have $\mathbb{M}_{NT}(\theta )=N^{-1}\sum_{i=1}^{N}\mathbb{M} _{i,T}(\lambda _{i},F)=T^{-1}\sum_{t=1}^{T}\mathbb{M}_{t,N}(\Lambda ,f_{t})$ . The main difficulty in finding the global minimum of $\mathbb{M}_{NT}$ is that this object function is not convex in $\theta $. However, for given $F$ , $\mathbb{M}_{i,T}(\lambda ,F)$ happens to be convex in $\lambda $ for each $i$ and likewise, for given $\Lambda $, $\mathbb{M}_{t,N}(\Lambda ,f)$ is convex in $f$ for each $t$. Thus, both optimization problems can be efficiently solved by various linear programming methods (see Chapter 6 of koenker2005quantile). Based on this observation, we propose the following iterative procedure:
{Iterative quantile regression (IQR): }\newline {Step 1: }Choose random starting parameters: $F^{(0)}$. \newline {Step 2: }Given $F^{(l-1)}$, solve $\lambda _{i}^{(l-1)}=\operatorname*{arg\min} _{\lambda }\mathbb{M}_{i,T}(\lambda ,F^{(l-1)})$ for $i=1,\ldots ,N$; given $ \Lambda ^{(l-1)}$, solve $f_{t}^{(l)}=\operatorname*{arg\min}_{f}\mathbb{M}_{t,N}(\Lambda ^{(l-1)},f)$ for $t=1,\ldots ,T$. \newline {Step 3: }For $l=1,\ldots ,L$, iterate the second step until $ \mathbb{M}_{NT}(\theta ^{(L)})$ is close to $\mathbb{M}_{NT}(\theta ^{(L-1)}) $, where $\theta ^{(l)}=(\text{vech}(\Lambda ^{(l)})^{\prime }, \text{vech}(F^{(l)})^{\prime })^{\prime }$. \newline {Step 4: }Normalize $\Lambda ^{(L)}$ and $F^{(L)}$ so that they satisfy the normalizations in (ref).
To see the connection between the IQR algorithm and the PCA estimator of bai2003inferential, suppose that $r=1$, and replace the check function in the IQR algorithm by the least-squares loss function. Then, it is easy to show that the second step of the algorithm above yields $\Lambda ^{(l-1)}=(X'F^{ (l-1)})/\Vert F^{(l-1)}\Vert ^{2}$ and $ F^{(l)}=(X\Lambda ^{(l-1)})/\Vert \Lambda ^{(l-1)}\Vert ^{2}=XX'F^{(l-1)}/C_{l-1}$, where $X$ is the $T\times N$ matrix with elements $ \{X_{it}\}$, and $C_{l}=\Vert F^{(l)}\Vert ^{2}\cdot \Vert \Lambda ^{(l)}\Vert ^{2}$. Thus, with proper normalizations at each step, the iterative procedure is equivalent to the well-known power method of hotelling1933analysis, and the sequence $F^{(0)},F^{(1)},\ldots $ will converge to the eigenvector associated with the largest eigenvalue of $XX^{\prime }$. In the more general case $r>1$, if we replace the check function in the IQR algorithm by the least-squares loss function and normalize $F^{(l-1)},\Lambda^{(l-1)}$ to satisfy (ref) at step 2, it can be shown that the above iterative procedure is similar to the method of orthogonal iteration (see Section 7.3.2 of golub2013matrix) for calculating the eigenvectors associated with the $r$ largest eigenvalues of $XX'$, which is the PCA estimator of bai2003inferential. Therefore, the IQR algorithm and its corresponding QFA estimator can be viewed as an extension of PCA to QFM.
Similar algorithms have been proposed in the machine learning literature to reduce the dimensions for binary data, where the check function is replaced by some smooth nonlinear link functions, e.g. collins2002generalization. However, unlike PCA, whether such methods guarantee finding the global minimum remains an important open question. Nonetheless, in all of our Monte Carlo simulations we found that the QFA estimators of the factors using the IQR algorithm always converge to the space of the true factors, which is somewhat reassuring in this respect.
To prove the consistency of the QFA estimator $\hat{\theta}$, we make the following assumptions:
Assumptions 1 (i) is essentially the strong factors assumption that is standard in the literature (see Assumption B of bai2003inferential). The requirement that $\sigma_1,\ldots,\sigma_r$ are distinct is similar to Assumption G of bai2003inferential, which is a convenient assumption to order the factors. Assumptions 1 (ii) and (iii) are similar to (C1) and (C2) in \textcolor{blue}{AB (2020)}, except that we do not require moments of $u_{it}$ to exist. Also notice that Assumption (iii), which allows for both cross-sectional and time series heteroskedasticity, requires the idiosyncratic errors to be mutually independent. This stems from the use of Hoeffding's inequality in the proofs of some results, which provides a sub-Gaussian tail bound for the sum of bounded independent random variables. There have been attempts to relax this assumption (see Remark 1.4 below) but it is difficult to characterize the minimal set of conditions that the error terms should satisfy to achieve the sub-Gaussian inequality required in our proofs. Notice, however, that in exchange for the independence assumption, we can dispense with the bounded moment conditions in the idiosyncratic terms, whose violation would render PCA invalid. At any rate, in sub-section 5.2 we run some Monte Carlo simulation on the performance of our QFA estimation when error terms are allowed to exhibit mild cross-sectional and serial dependence in order to check the robustness of our results to these features.
Write $\hat{\Lambda}=(\hat{\lambda}_{1},\ldots ,\hat{\lambda}_{N})^{\prime }$ , $\Lambda _{0}=(\lambda _{01},\ldots ,\lambda _{0N})^{\prime }$, $\hat{F}=( \hat{f}_{1},\ldots ,\hat{f}_{T})^{\prime }$, $F_{0}=(f_{01},\ldots ,f_{0T})^{\prime }$, and let $L_{NT}=\min \{\sqrt{N},\sqrt{T}\}$. The following theorem provides the average rate of convergence of $\hat{\Lambda}$ and $\hat{F}$.
The sign matrix $\mathsf{S}$ appears in the above result due the intrinsic sign indeterminacy of factors and loadings -- that is, the factor structure remains unchanged if a factor and its loading are both multiplied by $-1$ (e.g., see Theorem 1.b of stock2002forecasting for a similar result).
Remark 1.1: Since our proof strategy is substantially different from that of \textcolor{blue}{BN (2002)}, we briefly sketch the main ideas underlying our proof here. To facilitate the discussion, for any $\theta _{a},\theta _{b}\in \Theta ^{r}$ define the semimetric $d$ by:
and let
The semimetric $d$ plays an important role in our asymptotic analysis. We first show that $d(\hat{\theta},\theta _{0})=o_{P}(1)$. Next, it can be shown that:
and that for sufficiently small $\delta >0$,
where $\Theta ^{r}(\delta )=\{\theta \in \Theta ^{r}:d(\theta ,\theta _{0})\leq \delta \}$. Intuitively, the above two inequalities and $d(\hat{ \theta},\theta _{0})=o_{P}(1)$ imply that $d^{2}(\hat{\theta},\theta _{0})\lesssim d(\hat{\theta},\theta _{0})/L_{NT}$, or $d(\hat{\theta},\theta _{0})\lesssim L_{NT}^{-1}$. Then, the desired results follow from the fact that $\Vert \hat{\Lambda}-\Lambda _{0}\mathsf{S}\Vert /\sqrt{N}+\Vert \hat{F} -F_{0}\mathsf{S}\Vert /\sqrt{T}\lesssim d(\hat{\theta},\theta _{0})$.
Inequality (ref) follows easily from a Taylor expansion of $\bar{ \mathbb{M}}_{NT}(\hat{\theta})$ around $\theta _{0}$, together with Assumption 1(ii). It is worth stressing that the proof of (ref) requires the chaining argument which is commonly used in the theory of empirical processes. In particular, using Hoeffding's inequality and the fact that $ |\rho _{\tau }(u)-\rho _{\tau }(v)|\leq 2|u-v|$, it can be shown that, for any given $\theta _{a},\theta _{b}\in \Theta ^{r}$,
for some constant $K$. Then, along the lines of Theorem 2.2.4 of vanweak, it follows that the left-hand side of (ref) is bounded (up to a positive constant) by $\int_{0}^{\delta }\sqrt{\log D(\epsilon ,d,\Theta ^{r}(\delta ))} d\epsilon /\sqrt{NT}$. Finally, we can prove that $\int_{0}^{\delta }\sqrt{ \log D(\epsilon ,d,\Theta ^{r}(\delta ))}d\epsilon \lesssim \delta \sqrt{M}$, from which inequality (ref) follows.
Remark 1.2: Compared to \textcolor{blue}{BN (2002)}, recall that, in exchange for Assumption 1(iii) we do not require any moment of $ u_{it}$ to be finite. Thus, for the canonical AFM (e.g., Example 1) where the idiosyncratic errors have median equal to zero and satisfy Assumption 1(iii), our estimator for the case $\tau =0.5$ can be interpreted as a least absolute deviation (LAD) estimator which is robust to heavy tails and outliers. In relation to this issue, it is important to point out that the LAD estimator is related to robust PCA in the machine learning literature that aims to recover a low rank matrix from a large panel of observables. For example, the Principal Components Pursuit method proposed by candes2011robust features a combination of the $L_1$ norm (as in LAD) and a nuclear norm on the low rank matrix (see Chapter 3 of vidal2016generalized and bai2019rank for other robust PCA methods). In section 5 below, we will illustrate the robustness of the LAD estimator relative to the PCA estimator by Monte Carlo simulations.
Remark 1.3: If the true parameters do not satisfy the normalizations (ref), they can still be in the space $\Theta ^{r}$ after some normalizations. Let $H_{NT}$ be a $r\times r$ invertible matrix and define $\bar{f}_{0t}=H_{NT}^{\prime }f_{0t}$, $\bar{ \lambda}_{0i}=(H_{NT})^{-1}\lambda _{0i}$. Note that $\lambda _{0i}^{\prime }f_{0t}=\bar{\lambda}_{0i}^{\prime }\bar{f}_{0t}$. For $\{\bar{f}_{0t}\}$ and $\{\bar{\lambda}_{0i}\}$ to satisfy the normalizations (ref), we require:
where $\Sigma _{T,F}=T^{-1}\sum_{t=1}^{T}f_{0t}f_{0t}^{\prime }$, $\Sigma _{N,\Lambda }=\frac{1}{N}\sum_{i=1}^{N}\lambda _{0i}\lambda _{0i}^{\prime }$ , and $\mathbb{D}_{N}$ is a diagonal matrix with non-increasing diagonal elements. The above equalities imply that:
Thus, the rotation matrix $H_{NT}$ can be chosen as $\Sigma _{T,F}^{-1/2}\Gamma _{NT}$, where $\Gamma _{NT}$ is the matrix of eigenvectors of $\Sigma _{T,F}^{1/2}\Sigma _{N,\Lambda }\Sigma _{T,F}^{1/2}$. Note that when the eigenvalues of $\Sigma _{T,F}^{1/2}\Sigma _{N,\Lambda }\Sigma _{T,F}^{1/2}$ are distinct, its eigenvectors are unique up to signs, i.e., each eigenvector can be replaced by the negative of itself. As a result, Theorem 1 can be stated as follows:
where $\mathsf{S}$ is the diagonal matrix defined above. Notice that the rotation matrix $H_{NT}$ is slightly different from the rotation matrix of bai2003inferential. Moreover, because both $\lambda_{0i}$ and $f_{0t}$ are $\tau$-dependent, $H_{NT}$ also varies across quantiles, although we did not make it explicitly quantile dependent in the previous discussion to simplify notation.
Remark 1.4: Compared to \textcolor{blue}{BN (2002)}, our Assumption 1(iii) is admittedly strong. However, note that this assumption is made conditional on $\{f_{0t}\}$, so cross-sectional and temporal dependence of $u_{it}$ due to the common factors are still allowed for. Moreover, the independence assumption is only used to establish the sub-Gaussian inequality (ref). Thus, Assumption 1(iii) can be relaxed as long as the sub-Gaussian inequality holds.\footnote{ See van2002hoeffding for the properties of Hoeffding inequalities for martingales.}
In the previous section, we assumed the number of quantile-dependent factors $r(\tau )$ to be known at each $\tau $. In this subsection we propose two different procedures to select the correct number of factors at each quantile with probability approaching one. The first one selects the number of factors by rank minimization while the second one uses information criteria (IC). As before, the dependence of the quantile-dependent objects on $\tau $, including $r(\tau )$ , is suppressed for notational simplicity.
Let $k$ be a positive integer larger than $r$, and $\mathcal{A}^{k}$ and $ \mathcal{F}^{k}$ be compact subsets of $\mathbb{R}^{k}$. In particular, let us assume that $[ \lambda_{0i}^{\prime }\quad \mathbf{0}_{1\times (k -r ) } ]^{\prime }\in \mathcal{A}^{k} $ for all $i$.
Let $\lambda _{i}^{k},f_{t}^{k}\in \mathbb{R}^{k}$ for all $i,t$ and write $ \theta ^{k}=(\lambda _{1}^{k^{\prime }},\ldots ,\lambda _{N}^{k^{\prime }},f_{1}^{k^{\prime }},\ldots ,f_{T}^{k^{\prime }})^{\prime }$, $\Lambda ^{k}=(\lambda _{1}^{k},\ldots ,\lambda _{N}^{k})^{\prime }$, $ F^{k}=(f_{1}^{k},\ldots ,f_{T}^{k})^{\prime }$. Consider the following normalizations:
Define $\Theta ^{k}=\{\theta ^{k}:\lambda _{i}^{k}\in \mathcal{A} ^{k},f_{t}^{k}\in \mathcal{F}^{k},\text{ and }\lambda _{i}^{k},f_{t}^{k} \text{ satisfy \eqref{IC_norm}}\}$, and
Moreover, define $\hat{\Lambda}^{k}=(\hat{\lambda}_{1}^{k},\ldots ,\hat{ \lambda}_{N}^{k})^{\prime }$ and write
The first estimator of the number of factors $r$ is defined as:
where $P_{NT}$ is a sequence that goes to 0 as $N,T\rightarrow \infty $. In other words, $\hat{r}_{\text{rank}}$ is equal to the number of diagonal elements of $(\hat{\Lambda}^{k})^{\prime }\hat{\Lambda}^{k}/N$ that are larger than the threshold $P_{NT}$. We call $\hat{r}_{\text{rank}}$ the rank-minimization estimator because, as discussed below in Remark 2.1, it can be interpreted as a rank estimator of $(\hat{\Lambda}^{k})^{\prime }\hat{\Lambda }^{k}/N$.
It can then be shown that:
Remark 2.1: In the proof of Theorem 2, we show that for $k>r$, it holds that (up to sign)
where $\hat{F}^{k,r}$ is the first $r$ columns of $\hat{F}^{k}$ and $\Lambda _{0}^{\ast }=[\Lambda _{0},\mathbf{0}_{N\times (k-r)}]$. It then follows from Assumption 1 that $\hat{\sigma}_{N,j}^{k}\overset{p}{\rightarrow } \sigma _{j}>0$ for $j=1,\ldots ,r$ and $\hat{\sigma}_{N,j}^{k}=N^{-1} \sum_{i=1}^{N}\left( \hat{\lambda}_{i,j}^{k}\right) ^{2}=O_{P}(1/L_{NT}^{2})$ for $j=r+1,\ldots,k$. Thus, the first $r$ diagonal components of $(\hat{\Lambda}^{k})^{\prime } \hat{\Lambda}^{k}/N$ converge in probability to positive constants while the remaining diagonal components are all $O_{P}(1/L_{NT}^{2})$. In other words, $(\hat{\Lambda}^{k})^{\prime }\hat{\Lambda}^{k}/N$ converges to a matrix with rank $r$, and $P_{NT}$ can be viewed as a cutoff value to choose the asymptotic rank of $(\hat{\Lambda}^{k})^{\prime }\hat{\Lambda}^{k}/N$.
The second estimator of $r$ is similar to the IC-based estimator of \textcolor{blue}{BN (2002)}. Let $l$ denote a positive integer smaller or equal to $k$, and $\mathcal{A}^{l}$ and $\mathcal{F}^{l}$ be compact subsets of $\mathbb{R}^{l}$. In particular, for $l>r$, assume that $ [ \lambda_{0i}^{\prime }\quad \mathbf{0}_{1\times (l -r ) } ]^{\prime }\in \mathcal{A}^{l} $ for all $i$. Moreover, we can define $\Theta ^{l}, \hat{ \theta}^{l} , \hat{f}^l_{t},\hat{\lambda}_i^l, \hat{F}^l$ and $\hat{\Lambda} ^l$ in a similar fashion.
Define the IC-based estimator of $r$ as follows:
We can show that:
Remark 3.1: \textcolor{blue}{AB (2020)} obtain a similar result, but the difference with ours is that we only need the density function of the idiosyncratic errors to be uniformly bounded above and below, while \textcolor{blue}{AB (2020)} requires all the moments of the errors to be bounded. The reason why we can obtain the same result here with less restrictions is that our proof is based on the innovative argument discussed in Remark 1.1 and on the average convergence rate of the estimators, while the proof of \textcolor{blue}{AB (2020)} depends on the uniform convergence rate of the estimators.
Remark 3.2: Let $X$ denote the $T\times N$ matrix of observed variables, and let $\check{F} ^{l},\check{\Lambda}^{l}$ denote the matrices of PCA estimators of \textcolor{blue}{BN (2002)} when the number of factors is specified as $l$. Then \textcolor{blue}{BN (2002)}'s estimator of $r$ can be written as:
$k>r$, and $P_{NT}$ is defined as in Theorem 2 above. It can be shown that IC-based estimator $\hat{r}$ is equivalent to the number of diagonal elements in $\check{\Lambda}^{k^{\prime }}\check{ \Lambda}^{k}/N$ that are larger than $P_{NT}$. Thus, the two seemly different estimators of the number of factors are equivalent in AFM. However, due to the differences of the object functions, such equivalence does not exist in QFM.
Remark 3.3: The choice of $P_{NT}$ for $ \hat{r}_{\text{rank}}$ and $\hat{r}_{\text{IC}}$ can be different in practice. In particular, it can differ from those penalties used by \textcolor{blue}{BN (2002)}. \textcolor{blue}{AB (2020)} choose
for $\hat{r}_{\text{IC}}$, similar to $IC_{p1}$ of \textcolor{blue}{BN (2002)} . However, as shown in \textcolor{blue}{AB's (2020)} simulation results, this choice does not perform very well even for $N,T$ as large as 300.
Remark 3.4: Even though $\hat{r}_{\text{ rank}}$ and $\hat{r}_{\text{IC}}$ are both consistent estimators of $r$, the computational cost of $\hat{r}_{\text{rank}}$ is much lower than that of $ \hat{r}_{\text{IC}}$, because for $\hat{r}_{\text{rank}}$ we only estimate the model once, while for $\hat{r}_{\text{IC}}$ we need to estimate the model $k$ times. Thus, in the simulations and empirical applications we will focus on $\hat{r}_{\text{ rank}}$, and we refer to \textcolor{blue}{AB (2020)} for the corresponding simulation results of $\hat{r}_{ \text{IC}}$. In particular, we find that the choice
for $\hat{r}_{\text{rank}}$ works fairly well as long as $\min \{N,T\}$ is 100. This is also the value used in all of our simulations and applications.
The derivation of the asymptotic distribution of the QFA estimator $\hat{\theta}$ becomes a difficult task due to the non-smoothness of the check function and the problem of incidental parameters. As in the asymptotic analysis of conventional QR, one can expand the expected score function (which is smooth and continuously differentiable) and obtain a stochastic expansion for $\hat{ \lambda}_{i}-\mathsf{S}\lambda _{0i}$; yet the following term appears in the expansion:
The next step would be to show that the above expression is a higher-order term (i.e. $o_P(T^{-0.5})$) thus it does not affect the asymptotic distribution of $\hat{\lambda}_{i}$. However, due to the presence of the indicator functions in (ref), this is not an easy task. To see this, let's consider a similar problem for the PCA estimators of AFM. Let $ \check{\lambda}_{i}$ and $\check{f}_{t}$ be the PCA estimators. In the stochastic expansion of $\check{\lambda}_{i}-\lambda _{0i}$, the analogous term to (ref) happens to be:
where $\epsilon _{it}$ is the idiosyncratic error in the AFM. Note that, based on the result $ T^{-1}\sum_{t=1}^{T}\Vert \check{f}_{t}-f_{0t}\Vert =O_{P}(L_{NT}^{-1})$, one can only show that: \[ \left\| \frac{1}{T}\sum_{t=1}^{T}\epsilon _{it}(\check{f}_{t}-f_{0t}) \right\| \leq \sqrt{ \frac{1}{T}\sum_{t=1}^{T}\epsilon_{it}^2 }\cdot \sqrt{ \frac{1}{T}\sum_{t=1}^{T}\| \check{f}_{t}-f_{0t} \|^2} =O_{P}(L_{NT}^{-1}) .\] Instead, one has to use the stochastic expansion of $\check{f}_{t}- f_{0t}$ to show that $ T^{-1}\sum_{t=1}^{T}\epsilon _{it}(\check{f}_{t}-f_{0t})=O_P(L_{NT}^{-2})$ (see the proof of Lemma B.1 of bai2003inferential). Likewise, to show that (ref) is $o_{P}(T^{-0.5})$, establishing the convergence rate of $\hat{f}_{t}-\mathsf{S} f_{0t}$ is not enough, and the stochastic expansion of $\hat{f}_{t}-\mathsf{S} f_{0t}$ is required. However, due the non-smoothness of the indicator functions, it is not clear how to explore this stochastic expansion in (ref).
To overcome this problem, we proceed to define a new estimator of $\theta _{0}$, denoted as $\tilde{\theta}$, relying on the following smoothed quantile regressions (SQR):
where
$K(z)=1-\int_{-1}^{z}k(z)dz$, $k(z)$ is a continuous function with support $ [-1,1]$, and $h$ is a bandwidth parameter that goes to 0 as $N,T$ diverge.
Define
for all $i,t$. We impose the following assumptions:
The above conditions are standard in SQR, with the exception of $(v)$. Note that, like galvao2016smoothed, we need $k(z)$ to be a higher-order kernel function to control the higher-order terms in the stochastic expansions of the estimators. However, galvao2016smoothed assume that $m^{-1}<c<1/3$ (or $m\geq 4$), while we need $m^{-1}<c<1/6$ (or $m\geq 8$). The difference is due to the fact that the incidental parameters ($\lambda_{i}$ and $f_t$) in QFM enter the model interactively, while in the panel quantile models considered by these authors there are no interactive fixed-effects.
Then, we can show that:
Remark 4.1: Similar to the proof of Theorem 1, we can show that
where the extra $O_P(h^{m/2})$ term is due the approximation bias of the smoothed check function. However, Assumption 2$(v)$ implies that $ 1/L_{NT}\gg h^{m/2}$, and then it follows that average convergence rates of $ \tilde{\Lambda}$ and $\tilde{F}$ are both $L_{NT}$.
Remark 4.2: Similar to Theorems 1 and 2 of bai2003inferential, we show that the new estimator is free of incidental-parameter biases. That is, the asymptotic distribution of $\tilde{ \lambda}_{i}$ is the same as if we would observe $\{f_{0t}\}$, and likewise the asymptotic distribution of $\tilde{f}_{t}$ is the same as if $\{\lambda _{0i}\}$ were observed. The proof of this result is not trivial. To see why this is the case, first define $\varrho (u)=[\tau -K(u/h)]u$ and $\mathbb{S} _{i,T}(\lambda ,F)=T^{-1}\sum_{t=1}^{T}\varrho (X_{it}-\lambda ^{\prime }f_{t})$, then we can write $\tilde{\lambda}_{i}=\operatorname*{arg\min}_{\lambda \in \mathcal{A}}\mathbb{S}_{i,T}(\lambda ,\tilde{F})$. Expanding $\partial \mathbb{S}_{i,T}(\tilde{\lambda}_{i},\tilde{F})/\partial \lambda $ around $ (\lambda _{0i},F_{0})$ yields
where $\varrho ^{(j)}(u)=(\partial /\partial u)^{j}\varrho (u)$. The key step is to show that the last two terms on the right-hand side of the above equation are both $o_{P}(1/\sqrt{T})$. This is relatively easier for the PCA estimator of bai2003inferential, since $(\tilde{f}_{t}-\mathsf{S} f_{0t})$ has an analytical form (like e.g. in equation A.1 of bai2003inferential). In our case, we would also need a stochastic expansion for $(\tilde{f}_{t}-\mathsf{S} f_{0t})$, which in turn depends on the stochastic expansion of $( \tilde{\lambda}_{i}-\mathsf{S} \lambda _{0i})$ due to the nature of factor models. As in chen2018panel, this problem can be partly solved by showing that the expected Hessian matrix is asymptotically block-diagonal (see Lemma 11 in the Online Appendix). However, the proof of chen2018panel is only applicable to a special infeasible normalization, namely $\sum_{i=1}^{N}\lambda _{0i}\lambda _{i}=\sum_{t=1}^{T}f_{0t}f_{t}^{\prime }$, while our proof of Lemma 11 allows for normalization (ref) and can be generalized to any of the other normalizations considered by bai2013principal that uniquely pin down the rotation matrix.
Remark 4.3: As discussed in Remark 1.3, if the true parameters do not satisfy the normalizations (ref), the results of Theorem 3 can be stated as
where $\mathsf{S}$ and $H_{NT}$ are defined in Remark 1.3, $\Sigma_F= \lim_{T\rightarrow\infty}\Sigma_{T,F}$, $\Sigma_{\Lambda}= \lim_{T\rightarrow\infty}\Sigma_{N,\Lambda}$, $H = \Sigma_F^{-1/2}\Gamma$, and $\Gamma$ is the matrix of eigenvectors of $\Sigma_F^{1/2}\Sigma_\Lambda \Sigma_F^{1/2}$.
\noindentRemark 4.4: Let $l(z)$ be a continuous kernel function with support $[-1,1]$ where $l^{(j)}(z) =\partial^j l(z)/\partial z^j$ exists and $|l^{(j)}(z)|$ is bounded for $j=1,2$. Let $b$ a bandwidth. Estimators for the asymptotic variance matrices of $\tilde{\lambda}_i$ and $\tilde{f}_t$ can be simply constructed as \[ \tilde{V}_{\lambda_i} =\tau(1-\tau)\tilde{\Phi}_i^{-2} \text{ where } \tilde{\Phi}_i = \frac{1}{Tb}\sum_{t=1}^{T} l(\tilde{u}_{it}/b) \cdot \tilde{f}_t\tilde{f}_t',\] and \[ \tilde{V}_{f_t} =\tau(1-\tau)\tilde{\Psi}_t^{-1}\tilde{\Sigma}_{\Lambda}\tilde{\Psi}_t^{-1} \text{ where } \tilde{\Psi}_t = \frac{1}{Nb}\sum_{i=1}^{N} l(\tilde{u}_{it}/b) \cdot \tilde{\lambda}_i \tilde{\lambda}_i', \quad \tilde{\Sigma}_{\Lambda} = \tilde{\Lambda}' \tilde{\Lambda}/N,\] with $\tilde{u}_{it} = X_{it} -\tilde{\lambda}_i'\tilde{f}_t$. In Section A.5 of the Online Appendix we show that under Assumptions 1 and 2, the above estimators of the asymptotic covariance matrices are consistent if $b\rightarrow 0$ and $Nb^3\rightarrow\infty$. Note that this is different from the usual condition $Nb^2\rightarrow\infty$ in standard quantile regressions (see e.g. powell1984least and angrist2006quantile). Moreover, the above estimators are also consistent for the asymptotic covariance matrices discussed in Remark 4.3.
\noindentRemark 4.5: A restrictive DGP within class (1) would be a QFM where the PCA factors coincide with the quantile factors and only the factor loadings are quantile dependent. The representation for such restricted subset of QFM is as follows:
As a result, the main objects of interest are the common factors and the quantile-varying loadings. Notice that, if the factors $f_{t}$ were to be observed, using standard QR of $X_{it}$ on $f_{t}$ would lead to consistent and asymptotically normally distributed estimators of $\lambda _{i}(\tau )$ for each $i$ and $\tau \in (0,1)$. However, since $f_{t}$ are not observable, a feasible two-stage approach is to first estimate the factors by PCA, denoted as $\hat{f}_{PCA,t}$, and next run QR of $X_{it}$ on $\hat{f}_{PCA,t}$ to obtain estimates of $\lambda _{i}(\tau )$ as follows:
As explained in chen2017qfm, unlike the QFA estimators (see Remark 1.2), this two-stage procedure requires moments of the idiosyncratic term $u_{it}$ to be bounded in order to apply PCA in the first stage. However, an interesting result (see chen2017qfm, Theorem 2) is that the standard conditions on the relative asymptotics of $N$ and $T$ allowing for the estimated factors to be treated as known do not hold when applying this two-stage estimation approach. In effect, while these conditions are $T^{1/2}/N\rightarrow 0$ for linear factor-augmented regressions (see bai2006confidence) and $T^{5/8}/N\rightarrow 0$ for nonlinear factor-augmented regressions (bai2008extremum), lack of smoothness in the object (check) function at the second stage requires the stronger condition $T^{5/4}/N\rightarrow 0$. Moreover, Theorem 3 in chen2017qfm shows how to run inference on the quantile-varying loadings (e.g., testing the null that they are constant across all quantiles or a subset of them).
In this section we report the results from several Monte Carlo simulations regarding the performance of our proposed QFM methodology in finite samples. In particular, we focus on four relevant issues: (i) how well does our preferred estimator of the number of factors perform relative to other methods when the distribution of the idiosyncratic errors in an AFM exhibits heavy tails or outliers, (ii) how well do PCA and QFA estimate the true factors under the previous circumstances, (iii) how robust is the QFA estimation procedure when the errors terms are serially and cross-sectionally correlated, instead of being independent, and (iv) how good are the normal approximations given in Theorem 4 for the QFA estimators based on SQR.
As pointed out in Remark 1.2, since the consistency of our QFA estimator does not require the moments of the idiosyncratic errors to exist, at $\tau =0.5$ it can be viewed as a robust QR alternative to the PCA estimators commonly used in practice. For the same token, our estimator of the number of factors should also be more robust to outliers and heavy tails than the IC-based method of \textcolor{blue}{BN (2002)}. In this subsection we confirm the above claims by means of Monte Carlo simulations.
We consider the following DGP:
where $f_{1t}=0.8f_{1,t-1}+\epsilon _{1t}$, $f_{2t}=0.5f_{2,t-1}+\epsilon _{2t}$, $f_{3t}=0.2f_{3,t-1}+\epsilon _{3t}$, $\lambda _{ji},\epsilon _{jt}$ are all independent draws from $\mathcal{N}(0,1)$, and $u_{it}\sim i.i.d \text{ }B_{it}\cdot \mathcal{N}(0,1) + (1-B_{it})\cdot\text{Cauchy}(0,1)$, where $B_{it}$ are i.i.d Bernoulli random variables with means equal to $0.98$ and $\text{Cauchy}(0,1)$ denotes the standard Cauchy distribution. In this way, approximately $2\%$ of the idiosyncratic errors are generated as outliers.
We consider four estimators of the number of factors $r$: two estimators based on $PC_{p1}$, $ IC_{p1}$ of \textcolor{blue}{BN (2002)}, the Eigenvalue Ratio (ER) estimator proposed by ahn2013eigenvalue and our rank-minimization estimator discussed in subsection 3.2, having chosen
We set $k=8$ for all four estimators, and consider $N,T\in \{50,100,200,500\}$.
Table 1 reports the following fractions:
for each estimator having run 1000 replications.
It becomes evident from the results in Table 1 that $PC_{p1}$ and $IC_{p1}$ almost always overestimate the number factors, while the ER estimator tends to underestimate them, though to a lesser extent than what $PC_{p1}$ and $IC_{p1}$ overestimate them. By contrast, our rank-minimization estimator chooses accurately the right number of factors as long as $\min \{N,T\}\geq 100$.
Next, to compare the PCA and QFA estimators of the common factors in the previous DGP, we assume that $r=3$ is known. We first get the PCA estimator (denoted as $\hat{F}_{PCA}$), and then obtain the QFA estimator at $\tau=0.5$ (denoted $\hat{F}_{QFA}^{0.5}$) using the IQR algorithm. Next, we regress each of the true factors on $\hat{F}_{PCA}$ and $\hat{F}_{QFA}^{0.5}$ separately, and report the average $R^{2}$ from 1000 replications in Table 2 as an indicator of how well the space of the true factors is spanned by the estimated factors.\footnote{All the $R^2$ we use in this section and the next section are adjusted $R^2$.} As shown in the first three columns of Table 2, while the PCA estimators are not very successful in capturing the true common factors, the QFA estimators approximate them very well, even when $N,T$ are not too large.
As discussed earlier, the overall findings reported in Tables 1 and 2 are in line with our theoretical results. In effect, while the standard PCA estimators of \textcolor{blue}{BN (2002)} fail to capture the true factors because they require the eighth moments of all the idiosyncratic errors to be bounded (unlike the DGP above), our QFA estimators succeed to do so since they only need the density function to exist and be continuously differentiable, like in the previous DGP. Thus, this simulation exercise provides strong evidence about the substantial gains that can be achieved by using QFA rather than PCA in those cases where the idiosyncratic error terms in AFM exhibit heavy tails and outliers.
In this subsection we consider the following DGP:
where $f_{1t}=0.8f_{1,t-1}+\epsilon _{1t}$, $f_{2t}=0.5f_{2,t-1}+\epsilon _{2t}$, $f_{3t}=|g_{t}|$, $\lambda _{1i},\lambda _{2i},\epsilon _{1t},\epsilon _{2t},g_{t}$ are all independent draws from $\mathcal{N}(0,1)$ , and $\lambda _{3i}$ are independent draws from $U[1,2]$. Following \textcolor{blue}{BN (2002)}, the following specification for $e_{it}$ is used:
where $v_{it}$ are independent draws from $\mathcal{N}(0,1)$ except in the second case below. The autoregressive coefficient $\beta $ captures the serial correlation of $e_{it}$, while the parameters $\rho $ and $J$ capture the cross-sectional correlations of $ e_{it}$. We consider four cases:
For each of the previous cases and each $\tau \in \{0.25,0.5,0.75\}$, we first estimate $\hat{r}$ using our rank-minimization estimator, having set $k$ and $P_{NT}$ as described in the previous subsection. Second, we estimate $\hat{r }$ factors by means of the QFA estimation approach, which we denote $\hat{F} _{QFA}^{\tau}$. Finally, we regress each of the true factors on $\hat{F} _{QFA}^{\tau}$ and calculate the $R^{2}$s. This procedure is repeated 1000 times where, for each $\tau$, we report the averages of $\hat{r}$ and the $R^{2}$s in these 1000 replications.
The results for Case 1 and Case 2 (where the heavy tails are captured this time by a Student(3) rather than by a Cauchy distribution) are reported in Table 3 and Table 4, respectively, for $N,T\in \{50,100,200\}$. Notice that for $\tau =0.25, 0.75$, we have $r(\tau )=3$ while, for $\tau =0.5$, we get $r(\tau )=2$, since the factor $f_{3t}$ does not affect the median of $X_{it}$. It can be observed that both our selection criterion and the QFA estimators perform very well in choosing the number of QFA factors and in estimating them. It should be noticed that at $\tau =0.25,0.75$ the estimation of the scale factor $f_{3t}$ is not as good as the mean factors $f_{1t},f_{2t}$ for small $N$ and $T$. However, such differences vanish as $N$ and $T$ increase.
The results for Case 3 and Case 4 are in turn reported in Table 5 and Table 6, respectively. It can be inspected that the QFA estimators still perform well, even though the independence assumption is violated in these DGPs. Thus, despite adopting independence in Assumption 1 (iii) for tractability in the proofs (see Remark 1.4), it seems that QFA estimation still works properly when the error terms are allowed to exhibit mild serial and cross-sectional correlations.
To evaluate the normal approximations of Theorem 4 for the estimators based on SQR, we consider the following DGP: \[X_{it} = \lambda_i f_t + f_t \epsilon_{it},\] where $f_t \sim i.i.d \text{ }\mathcal{U}(1,2)$ and they are normalized such that $F'F/T=1$, $\lambda_{i} \sim i.i.d \text{ }\mathcal{N}(0,1)$ and $\epsilon_{it} \sim i.i.d \text{ }\mathcal{N}(0,1)$. Note that since our Theorem 4 is conditional on the factors and the loadings, $f_t$ and $\lambda_i$ are fixed in the simulations. To smooth the indicator function, we use the following eighth-order kernel function (see muller1984smooth): \[ k(z) = \mathbf{1}\{|z|\leq 1\} \cdot \frac{3465}{8192} \left( 7 - 105z^2 + 462 z^4 -858 z^6 + 715 z^8 -221 z^{10}\right), \] while the Epanechnikov kernel $l(z) = 0.75(1-z^2)\cdot \mathbf{1}\{|z|\leq 1\} $ is applied to estimate the variance.
Figure 1 and Figure 2 plot the histograms of the standardized estimators of the factors: $\hat{V}_{f_t}^{-1}\sqrt{N}(\tilde{f}_t -f_{0t})$ at $\tau=0.25$, $t=T/2$ from 1000 replications\footnote{We choose the signs of $\tilde{f}_t$ such that $\mathsf{S}=1$.}, where $\hat{V}_{f_t}$ is estimated using the formula in Remark 4.4. To check how the bandwidths affect the finite distributions of the estimators, we display results for different choices of $h$ and $b$.
From the reported histograms of the standardized estimators and the superimposed density function of the standard normal distribution, it becomes clear that the asymptotic distributions given in Theorem 4 provide reasonably good approximations for the finite sample distributions of the estimators based on SQR, even for $N=T=50$, and that such approximations are not very sensitive to the choice of bandwidths.
In this section we argue that QFA could provide a useful tool for causal analysis, predictive exercises and the economic interpretation of factors. In particular, we focus on applying our proposed methodology to three different datasets related to climate, macro aggregates and stock returns.
In our first empirical application we investigate how $\text{CO}_2$ emissions affect temperatures, which is a long-standing issue in climate change science and economics (see e.g. hansen1981climate, and hsiang2018economist). The dataset (coined Climate for short) we use consists of the annual changes of temperature from 441 stations from 1917 to 2018 ($N=441,T=102$), drawn from the Climate Research Unit at the University of East Anglia, where information about global temperatures across different stations in the Northern and Southern Hemisphere is collected. The annual global $\text{CO}_2$ emissions data is downloaded from The Global Change Data Lab.
Table 7 (column labeled Climate) reports the estimated number of factors using $PC_{p1}$ of \textcolor{blue}{BN (2002)}, the ER estimator and the rank-minimization estimator for a grid of quantiles ranging from 0.01 to 0.99.\footnote{In all the applications, before estimating the factors and the number of factors, each variable is standardized to have zero mean and variance equal to one.} The maximum number of factors $k$ is set to $8$ for all estimators. $PC_{p1}$ selects the maximum number of factors (8), while the ER estimator selects only one. Thus, the tendencies to overestimate (resp. underestimate) the number of PCA factors by the former (resp. latter) criteria mirror our simulation results in Section 5.1. By contrast, it can be observed that the numbers of factors estimated by the rank-minimization vary across quantiles. In particular, the number of QFA factors decreases as we move away from the median.
To compare the QFA factors (denoted as $\hat{F}_{QFA}^{\tau}$) and the PCA factors (denoted as $\hat{F}_{PCA}$), we regress each element of $\hat{F}_{QFA}^{\tau}$ on the 8 PCA factors selected by $PC_{p1}$ and compute the $R^2$ in these regressions.\footnote{We choose the number of PCA factors estimated by $PC_{p1}$ in these regressions to play conservative.} The results are shown in the upper panel of Table 8, where it becomes clear that the median factors ($\hat{F}_{QFA}^{0.5}$) are highly correlated with the PCA factors, with all the $R^2$s above $0.95$. By contrast, the QFA factors at the upper and lower quantiles ($\tau=0.01,0.05,0.95,0.99$) exhibit much lower correlations with the PCA factors, with $R^2$s around $0.60$. Thus, there seems to be room for using QFA in this application.
Next, to analyze the impact of $\text{CO}_2$ emissions on climate change, bivariate Granger non-causality tests are implemented. We regress the QFA factors at each relevant quantile on their own lags and the lagged growth rates of $\text{CO}_2$ emissions, labeled $\Delta \log ( \text{CO}_2)$, where the lag length is chosen according to BIC. Table 9 reports the p-values of these tests. The results of this novel approach to analyze quantile causality indicate that the growth rate of $\text{CO}_2$ emissions strongly Granger causes the QFA factors at the lower quantiles ($\tau=0.01,0.05$), with p-values below 0.01, as well as some of the median factors, albait to a lessser extent (p-values below 0.04). Moreover, the null of Granger non-causality is not rejected for QFA factors at the upper quantiles. Given that $\text{CO}_2$ emissions lead to global warming, the results for the lower quantiles of temperatures are in line with the evidence reported by rivas2020trends. Using a similar climate dataset but different quantile techniques to ours, these authors find that global warming over the last century seems to be mainly due to a different behaviour in the lower tail than in the central and upper tails of the distribution of global temperatures.
In the second application, we extend the diffusion-index forecasting exercise popularized by stock2002forecasting to explore the predictive power of the QFA factors. The main goal is to extract a few common factors (by both PCA and QFA) from a large panel of macroeconomic variables, and then use these factors to forecast e.g. real GDP growth and the inflation rate.
The FRED-QD dataset (coined Macro here) is used to estimate PCA and QFA factors. This is a quarterly panel consists of 211 US macroeconomic variables from 1960Q1 to 2019Q2 ($N=211, T=238$). It emulates the popular dataset used by stock2002forecasting, but also contains several additional time series. The variables in this dataset are updated in a timely manner and can be downloaded for free.\footnote{Link to the dataset: http://research.stlouisfed.org/econ/mccracken/. We refer to mccracken2016fred for the details of a very similar dataset that contains monthly macroeconomic variables.} Before estimation, each series is transformed to be stationary using Matlab codes that are also available on the FRED-QD data website.
As with the climate data, while $PC_{p1}$ selects 8 PCA factors, ER selects only one. The estimated numbers of QFA factors are reported in Table 7 (column labeled Macro). As can be seen, the number of QFA factors varies significantly across different quantiles, pointing to the existence of a nonstandard factor structure for this dataset. Moreover, the middle panel of Table 8 reports the $R^2$s of regressing each of the QFA factors on the 8 PCA factors. It becomes clear that the QFA factors at $\tau$ close to $0.5$ are all well explained by the PCA factors. However, the first QFA factor at $\tau=0.9$ (denoted $\hat{F}_{QFA}^{0.9}$) and those at $\tau=0.95,0.99$ (denoted as $\hat{F}_{QFA}^{0.95}$ and $\hat{F}_{QFA}^{0.99}$) contain some extra information that could be potentially helpful for forecasting macroeconomic variables. Since $\hat{F}_{QFA}^{0.95}$ exhibits a very high correlation with $\hat{F}_{QFA}^{0.9}$ and $\hat{F}_{QFA}^{0.99}$, we exclusively focus on the predictive power of $\hat{F}_{QFA}^{0.9}$ and $\hat{F}_{QFA}^{0.99}$ in the subsequent analysis.
Let $y_{t+1}$ denote the realized value of real GDP growth/inflation at period $t+1$. The forecasting equation we consider is as follows: \[ y_{t+1} = \alpha + \sum_{j=0}^{p_{max}} \beta_j y_{t-j} + \gamma' F_t + \epsilon_{t+1},\] where $F_t$ is vector containing several unobserved common factors extracted from the large macroeconomic dataset. The predicted value of $y_{t+1}$, based on a vector of estimated factors $\hat{F}_t$, is simply constructed as $\hat{y}_{t+1} = \hat{\alpha}+ \sum_{j=0}^{\hat{p}} \hat{\beta}_{j} y_{t-j} + \hat{\gamma}' \hat{F}_t $, where $ \hat{\alpha},\hat{\beta}_{j},$ $\hat{\gamma}$ are OLS estimates of the coefficients and $\hat{p}$ is the optimal lag length according to BIC. We compare five different specifications for $F_t$: (i) $F_t=0$, which is the benchmark AR model, (ii) AR plus $\hat{F}_t$ only including $\hat{F}_{PCA}$, (iii) AR plus $\hat{F}_t$ including $\hat{F}_{PCA}$ and $\hat{F}_{QFA}^{0.9}$, (iv) AR plus $\hat{F}_t$ including $\hat{F}_{PCA}$ and $\hat{F}_{QFA}^{0.99}$, and (v) AR plus $\hat{F}_t$ including $\hat{F}_{PCA}$, $\hat{F}_{QFA}^{0.9}$ and $\hat{F}_{QFA}^{0.99}$. Following chudik2018one, the initial estimation period is 1960Q1 to 1989Q4 (120 periods), and the forecast evaluation period is split into pre-crisis (1990Q1 to 2007Q2) and crisis and recovery (2007Q3 to 2019Q2) sub-periods. A rolling window of 120 periods is used both to estimate the coefficients and generate the rolling forecasts. In particular, following chudik2018one, the number of mean factors is estimated using $PC_{p1}$ at each rolling window, where the maximum number of factors is set equal to 5.
The mean squared error (MSE) of these procedures, and their relative MSE (R-MSE) to the benchmark AR model are reported in Table 10 for the whole evaluation period and each relevant sub-sample. As can be observed, in regards to real GDP growth, adding the upper tail QFA factors ranks better in terms of R-MSE than the AR and AR+$\hat{F}_{PCA}$ models for the three considered periods. The gains are not sizeable but yet they are relevant. As for the inflation rate, the results are weaker, though there are some gains for the crisis and recovery period.
A well-known shortcoming of point forecasts is that their uncertainty is generally unknown, hence it is difficult to quantify their precision at any given period of time. To address this problem, it has became customary among central banks to report density forecasts for important macroeconomic variables. In this respect, adrian2019vulnerable argue that a simple way of producing such densities is via QR. Following these authors, we next evaluate the predictive power of the QFA factors for forecasting the densities of real GDP growth and inflation. In particular, we first predict the conditional quantiles of the target variable $y_{t+h}$ by $\hat{q}_{\tau,t+h} = \hat{\alpha}_{\tau}+ \sum_{j=0}^{p} \hat{\beta}_{\tau,j} y_{t-j} +\hat{\gamma}_{\tau}' \hat{F}_{\tau,t}$ for $\tau\in\{0.05,0.25,0.75,0.95\}$, where $ \hat{\alpha}_{\tau},\hat{\beta}_{\tau,j},\hat{\gamma}_{\tau}$ are estimated coefficients by running QR of $y_{t+h}$ on $[1,y_t,\ldots, y_{t-p}, \hat{F}_{\tau,t}]$, and $\hat{F}_{\tau,t}$ is a vector of estimated quantile factors using the IQR algorithm.\footnote{Only 1 QFA factor is estimated at $\tau=0.05, 0.95$, whereas 5 QFA factors are estimated at $\tau=0.25, 0.75$.} Next, given the predicted quantiles: $[\hat{q}_{0.05,t+h}, \hat{q}_{0.25,t+h},\hat{q}_{0.75,t+h},\hat{q}_{0.95,t+h} ]$, the predicted density of $y_{t+h}$ is constructed as the density of a skewed $t$-distribution by matching the predicted quantiles.\footnote{We refer to adrian2019vulnerable for the details and to azzalini2003distributions for the definition and properties of the skewed $t$-distribution.} Finally, the accuracy of the density forecast is measured by the predictive score, which is the predicted density evaluated at the realized value of $y_{t+h}$. Higher predictive scores indicate more accurate predictions. The out-of-sample density forecasts are constructed using rolling windows with the most recent 120 observations, and the evaluation period is 1990Q1 to 2019Q2. Moreover, we set $p=3$, and the benchmark model is the one where $\hat{F}_{\tau,t}=0$, i.e. the quantiles of $y_{t+h}$ are predicted only using its own lags. Figure 3 displays the predictive scores of the one-quarter-ahead ($h=1$) and one-year-ahead ($h=4$) density forecasts for both variables. It can be seen that in both instances the predictive scores of the“AR + Quantile Factors” procedure is frequently above that of the “AR benchmark” model, sometimes by a large margin, indicating that the QFA factors could indeed be very informative for density forecasting of highly relevant macroeconomic variables.
Our last application concerns the interpretation of the quantile factors extracted from financial asset returns. The dataset (Finance in short) contains the monthly returns of 429 stocks from 1980M01 to 2014M12 ($N=429, T=420$), obtained from The Center of Research for Security Prices (CRSP).\footnote{The panel is balanced by only keeping stocks that have no missing observations during this time period. }
Except at $\tau=0.5$, the estimated number of QFA factors reported in Table 7 (column labeled Finance) are all equal to 1, which agrees with the choice of PCA factors by the ER estimator but again is less than the 4 factors selected by $PC_{p1}$.\footnote{$PC_{p3}$ and $IC_{p3}$ of \textcolor{blue}{BN (2002)} chose 8 factors while all the other 6 information criteria choose 4 factors.}
The lower panel of Table 8 reports the $R^2$s of regressing each of the QFA factors on the 8 PCA factors. As can be inspected, most of these factors are well explained by the PCA factors, with the exception of those at $\tau=0.01, 0.99$, where the $R^2$s are below $60\%$. Interestingly, as discussed in the Introduction this evidence is seemingly consistent with the findings of the financial literature on the existence of tail factors in the distribution of asset returns, as reported e.g. by RePEc:aah:create:2018-02. Thus, it is interesting to check whether the extra quantile factors at the lower tail and upper tail of the returns distribution could yield some confirmation of that hypothesis.
First, as shown in the upper panel of Figure 4, the QFA factors at $\tau=0.01$ and 0.99 (both with variance $=0.12$) are much less volatile than those at $\tau=0.5$ (both with variance $=1$), meaning that the tails of the distributions of returns are more stable than the median. Second, we find that the interquantile range (defined as the difference between the quantile factors at $\tau=0.99$ and $0.01$) provides a good measure of uncertainty for financial markets.\footnote{The results with the interpercentile range and interquartile range turn out to be similar.} Finally, as shown in the lower panel of Figure 4, the interquantile range is highly correlated with the volatility factor (with a correlation of 0.87) constructed by applying PCA-SQ to the squared residuals of an AFM.\footnote{ Following renault2016apt, we first project out the 8 PCA factors from the returns, and the volatility factor is obtained as the cross-sectional average of the squared residuals.} On the contrary, the correlations between the two median factors and the volatility factor only reach 0.08 and -0.05, respectively. Thus, this evidence seems supportive of the the presence of extra common factors affecting the tails and the volatility of the asset returns, with our results providing a link between them.
Approximate Factor Models (AFM) have become a leading methodology for the joint modelling of large number of economic time series with the big improvements in data collection and information technologies. This first generation of AFM was designed to reduce the dimensionality of big datasets through finding those common components (mean factors) which, by shifting the means of the observed variables with different intensities, are able to capture a large fraction of their co-movements. However, one could envisage the existence of other common factors that do not (or not only) shift the means but also affect other distributional characteristics (volatility, higher moments, extreme values, etc.). This calls for a second generation of factor models.
Inspired by the generalization of linear regressions to quantile regressions (QR), this paper proposes Quantile Factor Models (QFM) as a new class of factor models. In QFM, both factors and loadings are allowed to be quantile-dependent objects. These extra factors could be useful for identification purposes, for instance mean factors vs. volatility/skewness/kurtosis factors, as well as for forecasting purposes in factor-augmented regressions and FAVAR setups.
Using tools in the interface of QR, Principal Component Analysis (PCA) and the theory of empirical processes, we propose an estimation procedure of the quantile-dependent objects in QFM, labelled Quantile Factor Analysis (QFA), which yields consistent and asymptotically normal estimators of factors and loadings at each quantile. An important advantage of QFA is that it is able to extract simultaneously all mean and extra (non-mean) factors determining the factor structure of QFM, in contrast to PCA which can only extract mean factors. In addition, we propose novel selection criteria to estimate consistently the number of factors at each quantile. Finally, another relevant result is that QFA estimators remain valid (under some restrictive assumption on the idiosyncratic error terms -- see Assumption1 (iii) -- which are adopted to simplify the proofs) when the idiosyncratic error terms in AFM exhibit heavy tails and outliers, a case where PCA is rendered invalid.
The previous theoretical findings receive support in finite samples from a range of Monte Carlo simulations. Furthermore, it is shown in these simulations that QFA estimation performs well when we depart from some of simplifying assumptions used in the theory section for tractability, like, e.g., independence of the idiosyncratic errors. Lastly, our empirical applications to three large panel datasets of financial, macro and climate variables provide evidence that some of these extra factors may be highly relevant in practice for causality analysis, forecasting, and economic interpretation purposes.
Any time a novel methodology is proposed, new research issues emerge for future investigation. Among the ones which have been left out of this paper (some are part of our current research agenda), four topics stand out as important: