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.
96,070 characters · 24 sections · 103 citation commands
Universal Factor Models
Factor models have a wide application in economics and finance. In these models, observable outcomes are affected by a few number of latent economic variables called factors. In asset pricing, researchers often assume that there are a small number of latent common shocks that drive asset returns. In macroeconomics, common shocks may have heterogeneous effects on cyclical variations. In synthetic control, variations in individuals' potential outcomes are assumed to be affected by latent time trends.
A frequent challenge for researchers is the possible existence of weak factors. When the outcome variable has little exposure to some factors, learning about them becomes difficult. Moreover, their presence can distort inference regarding other factors that exert a more prominent influence onatski2012asymptotics,bai2023approximate,giglio2025test. Because of this issue, popular methods, such as the principal component analysis (PCA, bai2002determining, bai2003inferential) and quantile factor analysis (QFA, chen2021quantile), often assume that all factors have strong impact on the mean or quantile of $Y$ at the quantile level of interest, usually referred to as the strong factor assumptions.
However, we argue in this paper that violating this type of strong factor assumptions does not necessarily imply that the factors are inherently weak; it could only imply the poor choice of the moment or quantile of the outcome variable used in estimation; this argument echoes a similar spirit in giglio2025test where they state that “the strength or weakness of a factor ... should not be viewed as a property of the factor itself; rather, it should be viewed as a property of the set of test assets used in estimation” (p.250). As a consequence, failure of PCA and QFA does not necessarily imply that the factors cannot be well-estimated. For illustration, consider a random coefficient model for an outcome variable $Y_{it},i=1,\ldots,N,t=1,\ldots,T$, $$Y_{it}=\sum_{j=1}^{r}\beta_{j}(U_{it})\lambda_{0,ij}^{*}f_{0,tj}^{*},$$ where $\lambda_{0,ij}^{*}$ and $f_{0,tj}^{*}$, $j=1,\ldots,r$, are $r$ latent factor loadings and factors for $i$ and $t$, respectively, and the unobserved $U_{it}$ is independent of the factors and uniformly distributed. Suppose $\beta_{j}(\cdot)$ is integrable on $(0,1)$ for all $j$ with the integral denoted by $\bar{\beta}_{j}$ and $\sum_{j=1}^{r}\beta_{j}(\cdot)\lambda_{0,ij}^{*}f_{0,tj}^{*}$ is strictly increasing on $(0,1)$, then the conditional mean and $\tau$-th quantile of $Y_{it}$ are
respectively. Let $f_{0,t}^{*}=(f_{0,t1}^{*},\ldots,f^{*}_{0,tr})'$ and $\lambda_{0,i}^{*}=(\lambda_{0,i1}^{*},\ldots,\lambda_{0,ir}^{*})'$. Suppose $\sum_{t=1}^{T}f^{*}_{0,t}f^{*'}_{0,t}/T$ and $\sum_{i=1}^{n}\lambda^{*}_{0,i}\lambda^{*'}_{0,i}/N$ are both positive definite uniformly in $T$ and $N$. Then the strength of factors in the sense of PCA and QFA respectively refers to the magnitude of the $\bar{\beta}_{j}$s (i.e., $\int_{0}^{1}\beta_{j}(\tau)d\tau$) and $\beta_{j}(\tau)$. One can see from here that a weak factor in the sense of PCA can be a strong factor in the sense of QFA, vice versa, and a weak factor at some $\tau$ in the sense of QFA can be a strong factor in the same sense just at a different $\tau'$. In fact, factor $f_{0,tj}^{*}$ is not inherently weak as long as its overall strength $\int_{\mathcal{U}}[\beta_{j}(\tau)]^{2}d\tau$ is not small for some $\mathcal{U}\subset (0,1)$.
This paper proposes a new factor analysis framework, the universal factor model (UFM), that only imposes restrictions on the overall strength of factors rather than on the strength particularly regarding the mean or quantile of $Y_{it}$ at a certain $\tau$. Our model has the following form, which nests the previous example
Let $\Lambda_{0}^{*}(\cdot)$ and $F_{0}^{*}$ be $N\times r$ and $T\times r$ matrices collecting all the loading functions and factors, respectively. We require that the factors have an overall strong impact on $Y_{it}$ in the sense that $\int_{\mathcal{U}}\Lambda^{*'}_{0}(\tau)\Lambda^{*}_{0}(\tau)d\tau/N$ is positive definite uniformly in $N$ on a fixed compact interval $\mathcal{U}\subseteq (0,1)$; since we do not require integrability to hold on the entire $(0,1)$, this restriction does not rule out $Y_{it}$ with heavy or fat tails. We refer to $F_{0}^{*}$ as universal factors.
One can view our model (ref) as a generalization of standard panel data models with fixed effects without observable regressors. For instance, by letting $f_{0,t}^{*}=(1,f_{0,t})'$ and $\lambda^{*}_{0,i}(U_{it})\coloneqq g(\lambda_{0,i},U_{it})\coloneqq (\lambda_{0,i}+\Phi^{-1}(U_{it}),1)$ for some inverse cumulative distribution function $\Phi^{-1}$, we obtain an additive fixed effect model $Y_{it}=\lambda_{0,i}+f_{0,t}+\varepsilon_{it}$ where $\varepsilon_{it}\coloneqq \Phi^{-1}(U_{it})$. Interactive fixed effects are also allowed by, for instance, redefining $\lambda^{*}_{0}(U_{it})$ as $(\lambda_{0,i}+\Phi^{-1}(U_{it}),\lambda_{0,i}+1)'$. Since fixed effects are commonly viewed as genuine latent economic variables, our $f_{0,t}^{*}$ captures latent time-specific variables while $\lambda^{*}_{it}(U_{it})$ absorbs the time-invariant unobservable individual heterogeneity $\lambda_{0,i}$ and time-varying idiosyncratic shocks $U_{it}$. One can also compare our model with the standard quantile regression model with observed interactive terms: $Y_{it}=\beta_{0}(U_{it})+X_{1i}\beta_{1}(U_{it})+X_{2t}\beta_{2}(U_{it})+X_{1i}X_{2t}\beta_{3}(U_{it})$, where $X_{1i}$ and $X_{2t}$ are observable. One can reparameterize this model as $Y_{it}=X_{i}^{*'}(U_{it})X_{t}^{*}$ by letting $X_{i}^{*}(U_{it})\coloneqq (\beta_{0}(U_{it})+X_{1i}\beta_{1}(U_{it}),\beta_{2}(U_{it})+X_{1i}\beta_{3}(U_{it}))'$ and $X_{t}^{*}=(1,X_{2t})'$. Hence, just like $X_{1i}$ and $X_{2t}$, the factors and loadings in our model, although latent, are still genuine economic variables. A factor is weak only when the corresponding coefficient (marginal effect) is weak, very much like the case of weak instruments, which are themselves not weak (small); indeed, instruments are weak when they have small effects (correlation) on the endogenous variable.
Our model (ref) induces a quantile factor model (QFM, chen2021quantile) by assuming that $U_{it}$ is independent of $F_{0}^{*}$ and every component in $\Lambda_{0}^{*}(\cdot)F_{0}^{*'}$ is strictly increasing. Out assumption on the strength of the factors is weaker. To be specific, the strong factor assumption in chen2021quantile requires all factors that have any impact on the $\tau$-th quantile of $Y_{it}$ to have a strong impact. In contrast, our assumption holds as long as each factor is strong at some potentially different $\tau\in\mathcal{U}$.
A second advantage of our UFM modeling over QFM is that the factors in QFM are $\tau$-dependent; QFM thus only captures a subset of the universal factors $F_{0}^{*}$, and we show in this paper that in many models, such $\tau$-dependence is an unavoidable artificial parameterization only for their strong factor assumption to hold. Such ad hoc parameterization leads to difficulties in interpretation since in economics and finance, factors are simply unobserved latent variables, not to be $\tau$-dependent. Moreover, confusion may arise when the researcher wishes to use the recovered factors to augment regression in a different scenario. Our UFM modeling approach together with our relaxed requirement on the strength of factors avoids such artifact while entertaining a wider class of data generating processes.
Our model also induces a mean factor model, or, approximate factor model (AFM) in the form of bai2003inferential if $\Lambda^{*}(\cdot)$ is integrable on $(0,1)$. We show that our requirement on the strength of the factors is also weaker than their strong factor assumption for PCA to obtain a $\sqrt{N}$-asymptotically normal estimator of the factors.
To estimate the universal factors, we first propose our baseline estimator, referred to as the universal factor analysis (UFA). It begins by obtaining preliminary estimators of $F^{*}_{0}$ and $\Lambda^{*}_{0}(\tau)$ at various $\tau$s by iteratively conducting standard quantile regression or the smoothed version such as fernandes2021smoothing and he2023smoothed, which delivers both nice theoretical properties and computational efficiency. Using them, we estimate $F^{*}_{0}\int_{\tau\in\mathcal{U}}\Lambda_{0}^{*'}(\tau)\Lambda_{0}^{*}(\tau)d\tau F_{0}^{*'}/NT$. Its eigenvectors multiplied by $\sqrt{T}$ are the final estimated factors under our normalization. UFA consistently estimates the space spanned by the factors at the $\sqrt{N}$-rate regardless of whether smoothed quantile regression is adopted. In this paper, we only present the smoothed estimator for the sake of space.
The asymptotic expansions of the smoothed UFA estimator involve the densities of $Y_{it}$ evaluated at the $(i,t)$-th common component $\lambda_{0,i}^{*'}(\tau)f_{0,t}^{*}$ conditional on the universal factors. These $(i,t)$-specific nuisance parameters impose theoretical challenges to obtain sharper theoretical results. To sharpen the results to achieve asymptotic normality, we propose a two-stage inverse density weighted estimator, referred to as IDW-UFA, based on a new sample splitting strategy. In the first stage, we estimate the factors, loadings and thus $(i,t)$-specific inverse conditional densities by UFA using appropriate subsamples; we directly estimate these inverse densities exploiting the factor structure of our model instead of inverting an estimated density, achieving numerical stability. We then estimate the factors and loadings again in the second stage using the full dataset by weighting the original objective function in UFA by those estimated inverse densities to eliminate the incidental nuisance parameters in the asymptotic expansions.
The new sample splitting approach differs from the standard approach (e.g. chernozhukov2018double) in two ways. First, unlike the standard approach where different subsamples are usually governed by the same nuisance parameters, our nuisance parameters here are $(i,t)$-specific, so the ones obtained using one subsample do not apply to another. Our new approach solves this problem by delicately constructing subsamples. Second, our estimated nuisance parameters have a negligible impact on the second stage where the full dataset is used. Because of this, the IDW-UFA estimator does not suffer from the loss of efficiency or instability. This sample splitting strategy may be of independent interest; it can be applied to other scenarios where one needs to estimate a latent factor structure prior to estimating the main model; it also applies to estimating QFMs under the strong factor assumption; bai2021matrix tackle missing data problems in factor models in a similar spirit.
The UFA-IDW estimator achieves asymptotic normality (up to rotation) at the $\sqrt{N}$- or $\sqrt{T}$-rate for individual factors, loadings and common components, assuming that $N$ and $T$ are of the same order. The rotation matrix matches those in bai2003inferential and bai2023approximate for PCA.
In addition to our estimator of the $\tau$-specific loadings, we also develop a $\sqrt{T}$-asymptotically normal estimator of the factor loadings in the UFM-induced AFM based on our estimated universal factors. Hence, our estimator is useful even if the researcher is only interested in recovering the mean factors and their effects. Unlike PCA, it does not require the mean factors to be strong.
Finally, we propose an eigen/singular-value-thresholding type estimator for the total number of universal factors. We also propose consistent selectors of factors that have an arbitrary tolerated level of strength on the mean or quantile of $Y$.
Our paper adds to the growing literature on weak factors. onatski2012asymptotics, bai2023approximate, jiang2023revisiting, fan2024can, and choi2025high provide theoretical and simulation evidence showing that the PCA may be inconsistent or have a slower rate of convergence when the factors in an AFM are weak. bai2019rank propose a ridge penalized estimator which selects out the relevant strong factors in an AFM in a data driven way. onatski2010determining estimate the number of factors in an AFM with weak mean factors. Other methods developed in the presence of weak factors in mean models often impose conditions requiring the eigenvalues of $\Lambda_{0}^{\mu*'}\Lambda_{0}^{\mu*}/N$ not too small or the loading matrix to be sparse, e.g. de2008forecasting,lettau2020estimating,bailey2021measurement, freyaldenhoven2022factor,uematsu2022estimation. In our paper, these eigenvalues can be arbitrarily small as long as our identification condition mentioned earlier is satisfied. Meanwhile, we do not need to assume sparsity for $\Lambda_{0}^{\mu *}$. Moreover, our method handles both mean and quantile models. The weak factor problem is also relevant in the literature of risk premium estimation. anatolyev2022factor and giglio2025test develop methods to tackle weak factors in that setting. The focus and models are different from our paper. In panel data regression with interactive fixed effects, armstrong2022robust propose a robust estimation approach to construct bias-aware confidence intervals for the slope coefficients on the regressors when the fixed effects are weak. In contrast, we focus on the factors and loadings.
The rest of the paper is organized as follows. We formally set up the model in Section (ref). We introduce our baseline estimator UFA and provide its rate of convergence in Section (ref). Section (ref) presents the inverse density weighted estimator and derives its asymptotic distribution. Section (ref) demonstrates how to estimate the mean loadings. Section (ref) provides a consistent estimator of the total number of factors and consistent selectors of the factors in a QFM or an AMF with arbitrary tolerated strength. Section (ref) examines the finite sample performance of the estimator by Monte Carlo simulations. Section (ref) concludes. The Appendix and Online Appendix collect all the proofs.
Assume that observable $Y_{it}$ ($i=1,\ldots,N;t=1,\ldots,T$) follows the model in the Introduction,
where $\lambda^{*}_{0,i}(\tau)$ is an $r\times 1$ vector of factor loadings. We treat $\lambda_{0,i}^{*}(\tau)$ as deterministic whereas $f_{0,t}^{*}$ as realizations of some underlying random variables $f_{0,t}^{0*}$. All the statements in the rest of the paper are implicitly conditional on $f_{0,t}^{0*}=f_{0,t}^{*}$.\footnote{Alternatively, we can view the loading functions as random and condition on the realization of them as well.} Let the $T\times r$ matrix $F^{*}_{0}$ collect all the factors. We call our factors, i.e., the columns in $F_{0}^{*}$, the universal factors in contrast to the mean factors in AFMs such as in bai2003inferential and quantile factors in the QFM in chen2021quantile. Those factors are subsets of the universal factors since they only include the universal factors that have strong impact on the mean and a certain quantile of $Y_{it}$, respectively.
Throughout the paper, we maintain the following assumption which admits a linear conditional quantile representation for $Y_{it}$:
Let $q_{Y_{it}|f_{0,it}^{*}}(\tau),\tau\in (0,1)$, denote the $\tau$-th conditional quantile of $Y_{it}$. Then by Assumption (ref),
Let the $N\times T$, $N\times r$ and $T\times r$ matrices $Y$, $\Lambda_{0}^{*}(\tau)$ and $F^{*}_{0}$ collect all the observables, loadings and factors, respectively. We can rewrite (ref) in matrix form as:
Model (ref) also implies a factor structure for the conditional mean of $Y_{it}$, provided the existence of the latter, or equivalently, integrability of $\lambda_{0,i}^{*}(\cdot)$ on $(0,1)$. Let $\bar{\lambda}_{0,i}^{*}\coloneqq \int_{0}^{1}\lambda_{0,i}^{*}(\tau)d\tau$. Then
Therefore, by letting $\nu_{it}\coloneqq Y_{it}-\bar{\lambda}_{0,i}^{*'}f^{*}_{0,t}$,
In the rest of this section, we first introduce our assumption on the strength of the universal factors. We then compare our model and assumption with QFM and AFM.
Let $\mathcal{B}$ be a compact subset of $\mathbb{R}$ and $\mathcal{U}$ be a compact subset of $(0,1)$.
Part (ii) of Assumption (ref) is standard in the factor model literature. Under Parts (i) and (ii), the $r$ largest eigenvalues of the matrix considered in (iii) are bounded away from 0. Part (iii) further imposes distinctiveness of eigenvalues to guarantee uniqueness (up to column signs) of the eigenvectors. The boundedness of factors and loadings in part (iv) is identical to chen2021quantile.
Part (i) of Assumption (ref), on the other hand, is new. It imposes restrictions on the strength of universal factors. It allows the singular values of the loading matrix to be small at any given $\tau$. Note that integrability of the loadings implied by this assumption does not rule out heavy- or fat-tailed $Y_{it}$ because of the compactness of $\mathcal{U}$. Now for better illustration, we compare our model and assumption with the QFM in chen2021quantile and AFM in bai2002determining, bai2003inferential, and bai2023approximate.
chen2021quantile study the following quantile factor model:
where the factors are $\tau$-dependent. Let the number of factors at $\tau$ be $r(\tau)$. They require all the $r(\tau)$ factors to be strong at every $\tau$ of interest: For sufficiently large $N$,
Similar to our Assumption (ref)-(ii), they also maintain that for sufficiently large $T$,
Under model (ref) and by (ref), the columns in quantile factors $F^{q*}_{0}(\tau)$ form a subspace of the column space of universal factors $F^{*}_{0}$. The strong factor assumption (ref) thus says that, for every $\tau$, the factors in $F^{*}_{0}$ have either sufficiently large or exactly zero impact on the $\tau$-th conditional quantile of $Y_{it}$. However, our Assumption (ref)-(i) allows any universal factor to have arbitrarily weak impact at any $\tau\in\mathcal{U}$, as long as each universal factor is “strong” in the sense of chen2021quantile in a neighborhood of some $\tau$, and these neighborhoods can be different for different factors.
Besides the more stringent requirements on the strength of factors, another drawback of the QFM is the ad hoc dependence of the factors on the quantile level $\tau$. In chen2021quantile, they provide four examples where $Y_{it}$ is modeled as a function of factors, loadings and idiosyncratic shocks. As latent economic variables, none of the factor in those models depends on $\tau$, and all those four models can be nested by our model (ref). However, once transformed into the quantile model, dependence in $\tau$ can become unavoidable to satisfy their strong factor assumption (ref). Such artificial dependence on $\tau$ creates barriers to understand the underlying model and to use the obtained factors to augment regression in other scenarios. In contrast, our quantile representation can avoid such $\tau$-dependence because of our milder restrictions on the factors' strength.
The following example, taken from chen2021quantile, illustrates these points.
Recall that our model (ref) implies the following AFM provided that $\mathbb{E}(Y_{it}|f^{*}_{0,t})$ exists:
Compared with the standard AFM under the strong factor assumption:
we can see that the strong factor assumption for AFM is equivalent to:
Therefore, when all the factors are strong in the sense of bai2002determining and bai2003inferential, our Assumption (ref)-(i) holds by the Cauchy-Schwarz inequality.
Now we consider the case when some universal factors are weak mean factors. Let $\kappa(N)$ be the order of the smallest eigenvalue of $\left(\int_{0}^{1}\Lambda_{0}^{*}(\tau)d\tau\right)'\left(\int_{0}^{1}\Lambda_{0}^{*}(\tau)d\tau\right)$. If all factors all strong, $\kappa(N)\asymp N$. onatski2012asymptotics show that PCA is inconsistent if $\kappa(N)\asymp 1$ and bai2023approximate, jiang2023revisiting, fan2024can and choi2025high relax the strong factor assumption by allowing for $\kappa(N)=o(N)$ while still requiring $\kappa(N)\to\infty$. Under different requirements on $\kappa(N)$, these works show that PCA achieves $\sqrt{\kappa(N)}$-average and pointwise-in-$t$ rate of convergence for factors. In particular, fan2024can obtains $\sqrt{\kappa(N)}$-asymptotic normality for PCA when $\kappa(N)$ can be as small as $\log(N)$.
In contrast, our Assumption (ref) puts no restrictions on $\kappa(N)$ and even allows for $\kappa(N)=O(1)$ or $o(1)$, while still achieves $\sqrt{N}$-asymptotic normality. We illustrate this in the following example.
Under our conditional quantile model (ref), one can in principle use quantile regression to estimate the unknown factors and loadings. Since $\Lambda^{*}_{0}(\cdot)$ as a function on $\mathcal{U}$ is an infinite dimensional parameter to compute, we discretize $\mathcal{U}$ for feasibility. As for quantile regression, $\sqrt{N}$-rate can be obtained when estimating the space spanned by the factors no matter whether one uses the standard or some smoothed version of quantile regression. To save space, in the paper we only present the estimator based on smoothed quantile regression following fernandes2021smoothing and he2023smoothed because similar to chen2021quantile, smoothed quantile regression also leads to uniform consistency for factors at each $t$.
Let $(\tau_{1},\ldots,\tau_{M})$ be an equally spaced grid on $\mathcal{U}$.
By this observation, by letting $M$ be a large fixed constant or $M=h(T)$ where $h$ is a known increasing function, we can focus on estimating the factors and loadings at $\tau_{m},m=1,\ldots,M$, which leads to both an implementable estimator and tractable theoretical derivations. In what follows, we let $M$ be a sufficiently large fixed constant for simplicity. It is straightforward to extend the results in the paper to the case where $M$ is slowing growing with $T$.
Next, we impose normalization to the true parameters for identification. Let
where the right-hand side is an eigendecomposition of the left-hand side; the matrix $F_{0}/\sqrt{T}$ collects all the eigenvectors of the left-hand side, and $\sum_{m=1}^{M}\Lambda^{'}_{0}(\tau_{m})\Lambda_{0}(\tau_{m})/MN$ is a diagonal matrix with diagonal entries $\sigma_{1}^{2}\geq \cdots\geq \sigma_{r}^{2}$ equal to the eigenvalues.\footnote{The diagonal matrix of the eigenvalues still has the additive structure in $m$ as on the right-hand side of equation (ref). This is because by definition of eigenvectors, there exists an $m$-independent full-rank matrix $H$ such that $F_{0}^{*}/\sqrt{T}=F_{0}H/\sqrt{T}$. So, the eigenvalue matrix is equal to $\sum_{m=1}^{M}H\Lambda^{*'}_{0}(\tau_{m})\Lambda^{*}_{0}(\tau_{m})H'/MN$, and our $\Lambda_{0}(\tau_{m})=\Lambda_{0}^{*}(\tau_{m})H'$.} Then Lemma (ref) implies that \[\frac{F_{0}'F_{0}}{T}=I_{r}\text{ and }\sigma_{1}^{2}>\cdots>\sigma_{r}^{2}>0.\] Moreover, $f_{0,t}$ and $\lambda_{0,i}(\tau_{m})$ are also uniformly bounded; with a bit abuse of notation, still assume that both are in $\mathcal{B}^{r}$ for all $i,t$ and $m$. Meanwhile, $\lambda_{0,i}(\cdot)$ is also Lipschitz uniformly on $\mathcal{U}$ under Assumption (ref). Similar to chen2021quantile, in the rest of the paper, we will treat these diagonalized $F_{0}$ and $\Lambda_{0}(\tau_{m})$s as our parameters of interest for simplicity.
To achieve uniform consistency for the factors and loadings, we adopt the smoothed quantile regression in fernandes2021smoothing and he2023smoothed. Unlike horowitz1998bootstrap, galvao2016smoothed and chen2021quantile where the indicator function in the check function is smoothed by some CDF kernel, this version of smoothed quantile regression keeps the check function but smooths the empirical distribution. fernandes2021smoothing and he2023smoothed show that it has superior theoretical and computational properties compared to the traditional smoothing techniques.
Specifically, let $k(\cdot)$ be some smooth kernel function and $h$ be some bandwidth converging to 0 as $N,T\to\infty$. Under diagonalization of the true factors and loadings, we impose the following normalization for the estimator. Let $\Lambda(\cdot)\coloneqq (\Lambda(\tau_{m}))_{m=1,\ldots,M}$. Define the following parameter spaces:
By construction, the diagonalized true loadings and factors satisfy $(\Lambda_{0}(\cdot),F_{0})\in\Xi\times\mathcal{F}$.
Our baseline estimator, universal factor analysis (UFA), is defined as follows.
where $\rho_{\tau}(\cdot)$ is the check function at $\tau$. Note that we treat $r$ as known for the moment. A consistent estimator of $r$ is introduced in Section (ref).
\sloppy Let $\hat{R}_{h,t}(f;\Lambda(\cdot))\coloneqq\sum_{m}\int \rho_{\tau_{m}}(s)\sum_{i}k\left((s-(Y_{it}-\lambda_{i}'(\tau_{m})f))/h\right)ds/(MNh)$ and $\hat{R}_{h,i,\tau}(\lambda;F)\coloneqq\int \rho_{\tau}(s)\sum_{t}k\left((s-(Y_{it}-\lambda'f_{t}))/h\right)ds/(Th)$. We propose the following algorithm.
In Algorithm (ref), steps 1.1 and 1.2 can be carried out by gradient descent type methods by the smoothness of the objective function. Note that in these steps, $\lambda$ and $f$ are simply treated as $r\times 1$ real vectors without normalization. Compared to ando2019quantile and chen2021quantile, the most distinctive feature of Algorithm (ref) is the normalization from steps 2.1 to 2.4. The rationale behind step 2.4 is that $L_{0}(\tau)F_{0}/T=\Lambda_{0}(\tau)$ for all $\tau$ by definition.
We show in this section that, under the assumption $N$ and $T$ having the same order, the estimator (ref) can estimate the space spanned by the factors at $\sqrt{N}$-rate. This is equal to the rate of the QFA estimator in chen2021quantile and PCA in bai2003inferential under strong factors. We also derive uniform rates for $\hat{\lambda}_{i}(\tau_{m})$ and $\hat{f}_{t}$ that match the preliminary pointwise rates in Lemma S.5 in chen2021quantile.
Define $\varepsilon_{it}(\tau_{m})=Y_{it}-\lambda_{0,i}'(\tau_{m})f_{0,t}$. Denote the density of $\varepsilon_{it}(\tau)$ conditional on $F_{0}$ by $\textsf{f}_{\tau,it}(\cdot)$; note that $\textsf{f}_{\tau,it}(0)$ is also equal to the density of $Y_{it}$ conditional on the factors evaluated at the true common component. We now impose the following assumptions.
Similar to ando2019quantile and chen2021quantile, the independence assumption in Assumption (ref) is made so that we can adopt some concentration inequalities from the random matrix theory. Note that only conditional independence is assumed, so serial or cross-sectional correlation among $Y_{it}$s are allowed, captured by the correlation among the factors and loadings. Assumption (ref) ensures that the estimators satisfy the first order conditions. Assumptions (ref) to (ref) are similar to galvao2016smoothed and chen2021quantile. Note that Assumptions (ref)-(i) does not rule out the case of unbounded support since the lower bound $\underline{\textsf{f}}_{C}$ is $C$ specific. Assumption (ref) implies differentiability of $\lambda_{0,i}(\cdot)$. Assumption (ref) says that the kernel function is of bounded variation on $\mathbb{R}$ and has order $\gamma$. Compared with the aforementioned literature, our requirements on $c$ and $\gamma$ are stronger, needed for asymptotic normality for the inverse density weighted estimator in the next section; in this section, we can relax them to be $\gamma \geq 4$ and $\gamma^{-1}<c<1/2$. Due to Assumption (ref)-(ii), which is made for simplicity, we will use $N$ and $T$ exchangeably when discussing rates of convergence.
Let $\zeta_{NT}\coloneqq \sqrt{1/N}+\sqrt{1/T}$. Let $\|\cdot\|_{F}$ denote the Frobenius norm of a matrix. For a real number $a$, let $\text{sgn}(a)=1$ if $a\geq 0$ and $\text{sgn}(a)=-1$ if $a< 0$. Let $\hat{F}_{j}$ and $F_{0,j}$ be the $j$-th column in the $T\times r$ matrices $\hat{F}$ and $F_{0}$, respectively. Let $H_{NT,1}\coloneqq\text{diag}(\text{sgn}(\hat{F}_{j}'F_{0,j}))$ be an $r\times r$ diagonal matrix.
Theorem (ref)-(i) shows that, even if (some of) the factors are weak to the mean or some quantiles of $Y$, the space they span can be consistently estimated at the optimal rate, in contrast to bai2023approximate and chen2021quantile. Note that similar to chen2021quantile, Theorem (ref)-(i) can also be shown by replacing the smoothed objective function with the standard check-function-based objective function; this result does not rely on smoothing.
The uniform rate in Theorem (ref)-(ii), on the other hand, is slower than the optimal rate $\sqrt{\log N}/\sqrt{N}$. However, the current rate is sufficient to establish $\sqrt{N}$-asymptotic normality and to achieve the optimal uniform rate for the inverse density estimator to be introduced in the following section.
The main challenge to derive $\sqrt{N}$- and $\sqrt{T}$-asymptotic normality for $\hat{f}_{t}$ and $\hat{\lambda}_{i}(\tau_{m})$ defined in (ref) is the heterogeneity of the conditional density $\textsf{f}_{\tau,it}(0)$ in $i$ and $t$. To see it, let $\eta_{h,\tau_{m},it}\coloneqq K\left((\lambda_{0,i}'(\tau_{m})f_{0,t}-Y_{it})/h\right)-\mathbb{E}\left[K\left((\lambda_{0,i}'(\tau_{m})f_{0,t}-Y_{it})/h\right)\right]$ where $K(c)=\int_{-\infty}^{c}k(z)dz$. For simplicity, assume $H_{NT,1}=I_{r}$. Let $Q_{F,t}\coloneqq \sum_{m=1}^{M}\sum_{i=1}^{N}\textsf{f}_{\tau_{m},it}(0)\lambda_{0,i}(\tau_{m})\lambda_{0,i}'(\tau_{m})/MN$ and $Q_{\Lambda,mi}\coloneqq \sum_{t=1}^{T}\textsf{f}_{\tau_{m},it}(0)f_{0,t}f_{0,t}'/T$. Assumptions (ref) and (ref) ensure that $Q_{F,t}$ and $Q_{\Lambda,mi}$ are invertible for sufficiently large $N$ and $T$. We can show that the Taylor expansion of the first order conditions leads to the following stochastic expansion for $\hat{f}_{t}$:
where $A_{1}=\sum_{i=1}^{N}\sum_{m=1}^{M}\sum_{s=1}^{T}\textsf{f}_{\tau_{m},it}(0)\textsf{f}_{\tau_{m},is}(0)\lambda_{0,i}(\tau_{m})\lambda_{0,i}'(\tau_{m})\hat{f}_{s}f'_{0,s}Q_{\Lambda,mi}^{-1}/MNT$. The first term on the right-hand side is $\sqrt{N}$-asymptotically normal under our assumptions. The term $Q_{F,t}^{-1}A_{1}$ serves as an $r\times r$ rotation matrix on $f_{0,t}$. However, this rotation matrix loses the simplicity and interpretability compared to the rotation matrix in PCA bai2003inferential,bai2023approximate.
More importantly, terms such as \[A_{2t}=\sum_{i=1}^{N}\sum_{m=1}^{M}\sum_{s=1}^{T}\eta_{h,\tau_{m},it}\cdot\textsf{f}_{\tau_{m},is}(0)Q_{\Lambda,mi}^{-1}f_{0,s}\lambda_{0,i}'(\tau_{m})\left(\hat{f}_{s}-f_{0,s}\right)/MNT\] appear in the remainder. It is unclear whether the order of $A_{2}$ is $o_{p}(1/\sqrt{N})$: Although $A_{2}$ contains a mean zero random variable $\eta_{h,\tau_{m},it}$, the variable depends on $t$ instead of $s$, so it is not to be averaged out in the time dimension, and is correlated with $\tilde{f}_{s}$. The conditional density $\textsf{f}_{\tau_{m},is}(0)$, on the other hand, depends on all indices to be summed over; it thus prevents us to separately handle the average of $\eta_{h,\tau_{m},it}$ over $i$ and the average of $\hat{f}_{s}-f_{0,s}$ over $s$. Hence, we can only show that $A_{2}$ is $O_{p}(\zeta_{NT})$ by Theorem (ref).\footnote{An important special case is when $\textsf{f}_{\tau_{m},is}(0)$ also has a factor structure, i.e., there exist some $a_{\tau_{m},i}$ and $b_{\tau_{m},s}$ such that $\textsf{f}_{\tau_{m},is}(0)=a_{\tau_{m},i}'b_{\tau_{m},s}$. In that case, $A_{2}$ and in fact the whole remainder term are $o_{p}(\zeta_{NT})$ and thus UFA is asymptotically normal. One such an example is a model where $\textsf{f}_{\tau_{m},is}(0)$ only depends on $i$ or $s$; the Monte Carlo design for normal approximation in Section 5.3 in chen2021quantile satisfies this condition. Another example is when $r=1$: One can verify that $\textsf{f}_{\tau_{m},is}(0)=1/(\lambda_{0,i}'(\tau_{m})f_{0,t})$ (see equation (ref)), so it has a factor structure if $r=1$.}
Now if, instead, the objective function in (ref) is weighted by $1/\textsf{f}_{\tau_{m},it}(0)$ for each $m,i$ and $t$, then $Q_{F,t}$ becomes $\Phi\coloneqq \sum_{m=1}^{M}\sum_{i=1}^{N}\lambda_{0,i}(\tau_{m})\lambda_{0,i}'(\tau_{m})/MN$ whereas $Q_{\Lambda,mi}=I_{r}$. For $A_{2}$, since $\textsf{f}_{\tau,is}(0)$ is now cancelled out, we can rewrite the summations as follows
Indeed, we can show that the whole remainder is $o_{p}(1/\sqrt{N})$. Moreover, $Q_{F,t}^{-1}A_{1}$ now becomes $\left(F_{0}'\hat{F}/T\right)'$. We thus have the following equation which admits $\sqrt{N}$-asymptotic normality of the estimator.
One nice feature of this expansion is that the rotation matrix $F_{0}'\hat{F}/T$ is exactly equal to the rotation matrix $H_{NT,2}$ for PCA in Lemma 3 in bai2023approximate under $F_{0}'F_{0}/T=I_{r}$. bai2023approximate show that this matrix is equivalent to\footnote{Equivalence is in the sense that the difference of those two rotation matrices is $o_{p}(1/\sqrt{N})$.} the rotation matrix in bai2003inferential for PCA for an AFM under strong factors. Hence, it draws a close analogy between our inverse density weighted quantile estimator and PCA. We will further explain the intuition behind after we formally introduce our estimator in this section.
The above analysis is based on that $\textsf{f}_{\tau,it}(0)$ is known which in most applications is not the case. This motivates us to consider how to estimate those densities in a way such that the estimation error does not lead to further complications.
The observation in Section (ref) motivates us to reconstruct the objective function (ref) by weighting the kernel function at each $(m,i,t)$ by a consistent estimator of $1/\textsf{f}_{\tau_{m},it}(0)$. However, the estimation error of these inverse densities can be correlated with $\eta_{h,\tau_{m},it}$, causing technical challenges. Specifically, in the asymptotic expansion (ref), the first term on the right-hand side, which drives the asymptotic distribution, is now modified as follows:
where $\widehat{1/\textsf{f}_{\tau_{m},it}(0)}$ is some uniformly consistent estimator of $1/\textsf{f}_{\tau_{m},it}(0)$. If $\widehat{1/\textsf{f}_{\tau_{m},it}(0)}$ and $\eta_{h,\tau_{m},it}$ are correlated, it is unclear whether the second term is $o_{p}(1/\sqrt{N})$.
In this subsection, we introduce a novel subsample estimation approach to estimate $1/\textsf{f}_{\tau,it}(0)$ so that the problem above is avoided; the key observation is that we only need three quarters of data to estimate the whole set of factors and loadings. This idea may be of independent interest and can be applied to other two-stage estimation procedures that involve estimating some latent factor structure via in the first stage. An idea sharing a similar spirit can be found in bai2021matrix where they focus on factor analysis with missing data.
Recall that $\textsf{f}_{\tau,it}(0)$ is equal to the conditional density of $Y_{it}$ at $L_{0,it}(\tau)\coloneqq\lambda_{0,i}'(\tau)f_{0,t}$ given $f_{0,t}$. Since only $\lambda_{0,i}(\cdot)$ depends on $\tau$, one can show that (e.g. koenker1999goodness)
where $\lambda_{0,i}^{(1)}(\tau)$ is the derivative of $\lambda_{0,i}(\cdot)$ evaluated at $\tau$. For each $\tau$, we can approximate this derivative by numerical differentiation. koenker1999goodness approximate it by the three-point central difference formula; letting $h_{d}$ be some bandwidth, the approximation error is $O(h_{d}^{2})$. We propose to use a five-point difference formula (FPDF) to reduce the approximation error to $O(h_{d}^{4})$. For instance, we can use the central-difference-FPDF to approximate $\lambda_{0,i}^{(1)}(\tau)$ as follows:\footnote{For $\tau$ that is close to 0 or 1, one can instead use forward-difference-FPDF or backward-difference-FPDF, respectively, to mitigate the performance drop near the boundaries. The forward-difference formula reads \[ \frac{-25\lambda_{0,i}(\tau)+48\lambda_{0,i}(\tau+h_{d})-36\lambda_{0,i}(\tau+2h_{d})+16\lambda_{0,i}(\tau+3h_{d})-3\lambda_{0,i}(\tau+4h_{d})}{12h_{d}}, \] whereas the backward-difference formula is obtained by flipping the sign of $h_{d}$. Both forward and backward formulae have approximation error $O(h_{d}^{4})$ as the central-difference version.}
We can thus estimate the inverse density $1/\textsf{f}_{\tau_{m},it}(0)$ by substituting proper estimators of the factors and loadings at $\tau_{m}$ into (ref).
To avoid dependence between density estimation and the second stage estimation, we develop the following subsample estimation approach. Let $\mathcal{N}_{1}\coloneqq\{1,\ldots,\lfloor N/2\rfloor\}$, $\mathcal{N}_{2}\coloneqq\{\lfloor N/2\rfloor+1,\ldots,N\}$, $\mathcal{T}_{1}\coloneqq\{1,\ldots,\lfloor T/2\rfloor\}$, $\mathcal{T}_{2}\coloneqq\{\lfloor T/2\rfloor+1,\ldots,T\}$, where $\lfloor c\rfloor$ equals the largest integer that is no greater than $c$. Let $N_{a}=|\mathcal{N}_{a}|$ and $T_{b}=|\mathcal{T}_{b}|,a,b\in\{1,2\}$. Divide the $N\times T$ data matrix $Y$ into four regions:
where formally, for instance, $Top\ Left=\{Y_{it}:i\in\mathcal{N}_{1},t\in\mathcal{T}_{1}\}$. Let the subsample $Top\coloneqq Top\ Left\cup Top\ Right$. Subsamples $Bottom$, $Left$ and $Right$ are defined similarly.
The key observation is that we can estimate the full set of factors $\{f_{0,t}\}$ using, for example, $Top$ only, and estimate the full set of loadings $\{\lambda_{0,i}(\tau_{m}\pm h_{d}),\lambda_{0,i}(\tau_{m}\pm 2h_{d})\}$ using $Left$ only. Together, they consist of only approximately $3/4$ of the full data set. Under Assumption (ref), these estimated factors and loadings are by construction independent of all $Y_{it}$ in $Bottom\ Right$ even if they share the same $t$ and $i$ indices.
One subtlety when applying this idea is the mismatch of rotation matrices. If we separately obtain the factors and loadings using $Top$ and $Left$ by UFA, respectively, their rotation matrices may not cancel out when computing the product. Hence, we propose a sequential process. For illustration, suppose we are to estimate $f_{0,t}$ and $\lambda_{0,i}(\cdot)$ for every $(i,t)\in\mathcal{N}_{2}\times \mathcal{T}_{2}$, i.e., $Y_{it}\in Bottom\ Right$. We first obtain the whole set $\{\hat{f}_{s}^{top}\}$ where $s=1,\ldots,T$ by UFA using $Top$. Then, taking out the $i$-th row in $Left$, regress these $T_{1}$ data points of $Y$ onto the first half of the estimated factors $\{\hat{f}^{top}_{1},\ldots,\hat{f}^{top}_{\lfloor T/2\rfloor}\}$ by smoothed quantile regression. Denote the obtained loading by $\hat{\lambda}^{(t,l)}_{i}(\cdot)$, where $(t,l)$ indicates that the estimate is obtained using $Y$ in $Left$ and the estimated factors using $Top$. The estimators $\hat{\lambda}^{(t,l)}_{i}(\tau)$ and $\hat{f}^{top}_{t}$ are independent of any data point in $Bottom\ Right$ even though $(i,t)$ falls into that region. Now the product $\hat{\lambda}^{(t,l)'}_{i}(\tau)\hat{f}^{top}_{t}$ can consistently estimate $\lambda_{0,i}'(\tau)f_{0,t}$ for any $\tau$ because the rotation matrices for $\hat{\lambda}^{(t,l)}_{i}(\cdot)$ and $\hat{f}^{top}_{t}$ automatically cancel out by construction when doing the multiplication.
Generally, for an arbitrary fixed $(i,t)\in \mathcal{N}_{a}\times \mathcal{T}_{b}$, we estimate the inverse density by the following estimator:
where
We now make the following assumption and derive the uniform rate of convergence of (ref).
Assumption (ref) is the subsample counterpart of Assumption (ref). Note that the matrices $\sum_{m=1}^{M}\sum_{i\in\mathcal{N}_{a}}\lambda_{0,i}(\tau_{m})\lambda_{0,i}'(\tau_{m})/MN$ and $\sum_{t\in\mathcal{T}_{b}}f_{0,t}f_{0,t}'/T$ are in general no longer diagonal. Let $\psi_{NT}\coloneqq \zeta_{NT}/(hh_{d})+h_{d}^{4}$. We have the following theorem.
Finally, let us revisit the problem raised in the beginning of this subsection. We can see that the second term in equation (ref) can be split into two parts: \[ \frac{1}{N}\sum_{i\in \mathcal{N}_{1}}\sum_{m=1}^{M}\left(\widehat{\frac{1}{\textsf{f}_{\tau_{m},it}(0)}}-\frac{1}{\textsf{f}_{\tau_{m},it}(0)}\right)\eta_{h,\tau_{m},it}\lambda_{0,i}(\tau_{m})+\frac{1}{N}\sum_{i\in \mathcal{N}_{2}}\sum_{m=1}^{M}\left(\widehat{\frac{1}{\textsf{f}_{\tau_{m},it}(0)}}-\frac{1}{\textsf{f}_{\tau_{m},it}(0)}\right)\eta_{h,\tau_{m},it}\lambda_{0,i}(\tau_{m}). \] These two parts are correlated because the estimated density functions in each part are correlated with the $\eta_{h,\tau_{m},it}$ in the other part by construction. However, within each part, all the $\widehat{1/\textsf{f}_{\tau_{m},it}(0)}$s are independent of the $\eta_{h,\tau,it}$s. Hence, conditional on the $\widehat{1/\textsf{f}_{\tau_{m},it}(0)}$s, both parts are $o_{p}(1/\sqrt{N})$ by the Hoeffding's inequality.
Once we obtain the estimated densities, we use the full sample to estimate the factors and loadings. Hence, our estimator does not lose efficiency. Specifically, we define our inverse density weighted estimator (IDW-UFA) as follows:
Define $\tilde{R}_{h,t}(f;\Lambda(\cdot))$ and $\tilde{R}_{h,i,\tau}(\lambda;F)$ similar to $\hat{R}_{h,t}(f;\Lambda(\cdot))$ and $\hat{R}_{h,i,\tau}(\lambda;F)$ in Section (ref) with the estimated inverse densities as weights. We implement our estimator by the following algorithm.
In this subsection, we show that our inverse density weighted estimator can both estimate the spaces spanned by the factors and loadings at $\sqrt{N}$ or $\sqrt{T}$ rate, and is pointwise asymptotically normal up to rotation.
One key step to achieve asymptotic normality is to recenter the estimator. Following a similar argument as in Theorem (ref), we can show that, for instance, $\|\tilde{F}-F_{0}\tilde{H}_{NT,1}\|_{F}/\sqrt{T}=O_{p}(\zeta_{NT})$ for a diagonal matrix $\tilde{H}_{NT,1}$ whose $j$-th diagonal entry is $\text{sgn}(\tilde{F}_{j}'F_{0,j})$. However, to achieve asymptotic normality for $\tilde{f}_{t}$ for each $t$, letting $H_{NT,2}\coloneqq F_{0}'\tilde{F}/T$, we show that we need to recenter $\tilde{f}_{t}$ around $H_{NT,2}'f_{0,t}$ instead of $\tilde{H}_{NT,1}'f_{0,t}$. This result echoes bai2003inferential and bai2023approximate as $H_{NT,2}$ is equivalent to all their rotation matrices under $F_{0}'F_{0}/T=I_{r}$ as mentioned in Section (ref).
Recall $\Phi\coloneqq \sum_{m=1}^{M}\sum_{i=1}^{N}\lambda_{0,i}(\tau_{m})\lambda'_{0,i}(\tau_{m})/(MN)$. For any $\tau^{*}\in\{\tau_{1},\ldots,\tau_{M}\}$, let
In some applications, the parameters of interest are the factors and loadings affecting the mean, rather than the quantiles, of the outcome variable. However, directly estimating an AFM using PCA requires strong factors. Our method provides an alternative approach to estimate the mean factor loadings.
Under (ref), we have
if $\mathbb{E}(Y_{it}|f_{0,t})$ exists. The mean factor loading $\bar{\lambda}_{0,i}$ is by construction $\int_{0}^{1}\lambda_{0,i}(\tau)d\tau$. We can estimate $\bar{\lambda}_{0,i}$ for each $i$ by solving the following least square problem:
where $\tilde{f}_{t}$ is obtained by estimator (ref). Our normalization gives a simple analytical solution:
Let the mean common component matrix be $\bar{L}_{0}\coloneqq \bar{\Lambda}_{0}F_{0}'$. The estimator $\tilde{\bar{\lambda}}_{i}$ has the following properties.
Two remarks are in order. First, the boundedness assumption on $Y_{it}$ in the theorem is for simplicity; it can be replaced by, for example, the existence of higher order moments of $Y_{it}$. Second, all the variances can be consistently estimated by plugging in the estimated factors, loadings and conditional densities, under a similar argument as Remark (ref).
So far, we have assumed that $r$ is known or can be consistently estimated. In this section, we first propose a consistent estimator of $r$ that is robust to weak factors. We achieve this by estimating the common component $L_{0}(\tau_{m})$ for each $m$ by a nuclear norm penalized estimator that does not require strong factors. We also introduce estimators of the number of factors that have any tolerated level of influence on the conditional quantile or the conditional mean of the outcome variable; these “strong” factor selectors can be useful in applications where the researcher would like to include factors that have relatively large influence.
For a constant $C>0$ and compact interval $\mathcal{B}_{L}\subset\mathbb{R}$, for each $m=1,\ldots,M$, define
where $\|\cdot\|_{*}$ is the nuclear norm of a matrix. Applying the results in feng2023nuclear without regressors, we can show that $\hat{L}^{pel}(\tau_{m})$ is consistent of $L_{0}(\tau_{m})$ in the average squared Frobenius norm with the rate equal to $O_{p}(\log(NT)\max\{1/N,1/T\})$ uniformly in $m$, regardless of the order of the singular values of $L_{0}(\tau_{m})$. We can then estimate $\sum_{m=1}^{M}L_{0}'(\tau_{m})L_{0}(\tau_{m})/(MNT)$ by $\sum_{m=1}^{M}\hat{L}^{pel'}(\tau_{m})\hat{L}^{pel}(\tau_{m})/(MNT)$, where under $F_{0}'F_{0}/T=I_{r}$, all the nonzero eigenvalues of the former have order $O(1)$ by Lemma (ref). Therefore, between the $r$-th and the $(r+1)$-th largest eigenvalues of $\sum_{m=1}^{M}\hat{L}^{pel'}(\tau_{m})\hat{L}^{pel}(\tau_{m})/(MNT)$, denoted by $\hat{\sigma}^{2}_{r}$ and $\hat{\sigma}^{2}_{r+1}$, we can show that there is a sufficiently large gap with probability approaching 1. We thus propose the following thresholding estimator for $r$:
where $C_{r}$ is any sequence of $(N,T)$ satisfying $C_{r}\to 0$ and $\sqrt{\log(NT)}/(C_{r}\sqrt{\min\{N,T\}})\to 0$. The following theorem shows consistency of $\hat{r}$.
In applications, researchers may be interested in an AFM or a QFM at a specific quantile level, and only wish to include sufficiently influential factors. In this section, we propose a method in a similar spirit of $\hat{r}$ to select factors in each model that have any tolerated strength.
We start from the AFM (ref). Let the nonzero singular values of $\bar{L}_{0}/\sqrt{NT}\coloneqq \bar{\Lambda}_{0}F'_{0}/\sqrt{NT}$ be $\bar{\sigma}_{1}\geq \ldots \geq \bar{\sigma}_{r}$. Under $F_{0}'F_{0}/T=I_{r}$, the strong factor condition in the literature of AFM (e.g. bai2003inferential) refers to the case where $\bar{\sigma}_{j}$ has order $O(1)$, away from 0. A factor is weak in the sense of bai2023approximate refers to the case where $\bar{\sigma}_{j}$ has order $O(N^{\alpha_{j}/2-1/2})$ for some $j$ and $0<\alpha_{j}<1$.
We can show that, by the first part of Theorem (ref) and by $N\asymp T$, $|\tilde{\bar{\sigma}}_{j}-\bar{\sigma}_{j}|$ is $O_{p}(N^{-1/2})=o_{p}(N^{\alpha_{j}/2-1/2})$ for any $\alpha_{j}>0$ uniformly in $j=1,\ldots,\min\{N,T\}$, where $\tilde{\bar{\sigma}}_{j}$ is the $j$-th largest singular value of $\tilde{\bar{\Lambda}}\tilde{F}'/\sqrt{NT}$. Hence, for any $\alpha\in (0,1]$, we estimate the number of factors that influence the conditional mean of $Y$ with strength at least $\alpha$ by
where $C$ is an arbitrary constant.
Similarly, let $\sigma_{j}(\tau_{m})$ be the $j$-th largest singular value of $L_{0}(\tau_{m})/\sqrt{NT}$. The strong factor assumption in chen2021quantile is the case that $\sigma_{j}(\tau_{m})$ is $O(1)$ and away from 0 for all $j=1,\ldots,r$. Now similar to bai2023approximate, consider weak factors such that $\sigma_{j}(\tau_{m})$ has order $O(N^{\alpha_{j}/2-1/2})$ for some $j$ and $0<\alpha_{j}<1$. We estimate the number of factors that influence the conditional quantile of $Y$ at $\tau_{m}$ with strength at least $\alpha$ by
where $\tilde{\sigma}_{j}(\tau_{m})$ is the $j$-th largest singular value of $\tilde{\Lambda}(\tau_{m})\tilde{F}'/\sqrt{NT}$. \
Our factor selectors provide alternative approaches to select empirically relevant factors compared to the existing methods. Unlike freyaldenhoven2022factor who essentially imposes sparsity on $\bar{\Lambda}_{0}$ and requires that the smallest $\alpha_{j}$ to be greater than $1/2$ or bai2019rank who use a ridge penalty to filter out factors that have relatively small influence on the mean of $Y$, our method can select out either mean or quantile factors with any $\alpha>0$.
In this section, we demonstrate the finite sample performance of our estimators. We first compare the performance of UFA and IDW-UFA with QFA and PCA in estimating the factor space when there is a relatively weak quantile/mean factor. We then present the quality of Gaussian approximation of IDW-UFA in finite samples.
We let $r=1$ for simplicity. We consider sample size $(N,T)\in\{(50,50),(100,100),(150,150)\}$. We first draw $N\times 1$ and $T\times 1$ vectors $F^{*}_{0}$ and $\Lambda^{*}_{0,base}$ independently from $\text{Unif}[0,2]$. Draw an $N\times T$ matrix $U$ independently from $\text{Unif}[0,1]$. Let $\beta(U_{it})\coloneqq -0.99+2U_{it}$ and $\lambda_{0,i}^{*}(U_{it})\coloneqq \beta(U_{it})\lambda_{0,base,i}^{*}$. Construct $Y$ by $Y_{it}=\lambda_{0,i}^{*}(U_{it})f_{0,t}^{*}$. By construction, $q_{Y|F_{0}^{*}}(\tau)=\Lambda^{*}(\tau)F_{0}^{*'},\tau\in (0,1)$.
We can verify that $F_{0}^{*}$ is a strong universal factor because $\int_{0}^{1}\|\Lambda_{0}^{*}(\tau)\|^{2}_{F}d\tau/N$ is well bounded away from 0. However, $F_{0}^{*}$ is a “relatively weak” quantile factor near $\tau=0.5$ because $\|\Lambda_{0}^{*}(\tau)\|^{2}_{F}/N$ is close to 0 when $\tau$ is around 0.5, and a “relatively weak” mean factor because $\|\int_{0}^{1}\Lambda_{0}^{*}(\tau)d\tau)\|^{2}_{F}/N$ is close to 0, too; “relatively weak” because these values, though close to zero, are still fixed, not diminishing as $N$ and $T$ grow to infinity. Consequently, when the sample size is sufficiently large, we should expect that QFA in chen2021quantile and PCA still work. However, as shown below, these estimators' performance in the sample size we consider is not as well as UFA and IDW-UFA.
To implement our estimators, we first estimate $r$ by $\hat{r}$ proposed in Section (ref). We set $C$ in (ref) equal to $0.2$ and $C_{r}$ in (ref) equal to $1/(12(\min(N,T))^{1/3})$. Table (ref) presents the average, maximum and minimum $\hat{r}$ in 1000 simulation repetitions. The results show that our estimator of $r$ performs very good; for sample size $(50,50)$, only in 16 repetitions $r$ is overestimated by 1. Under larger sample size, $r$ is correctly estimated in all repetitions.
Next, we estimate the factor by UFA and IDW-UFA. To satisfy Assumption (ref), we let $\gamma=14$, set $h=1/\min(N,T)^{1/13}$, and choose the following fourteenth-order Gaussian-based kernel (see wand1990gaussian and marron1992exact): $k(z)\coloneqq \left(\sum_{i=0}^{6}c_{2i}z^{2i}\right)\phi(z)$ where $\phi$ is the standard normal density and
As discussed below Assumption (ref), a fourth-order kernel is sufficient to deliver $\sqrt{N}$-consistency for factor space estimation. In fact, we find in simulations that even a second-order kernel has similar performance. To save space, we present results under the current kernel just to be coherent with the assumption. For the inverse density estimator, we set $h_{d}=0.04$. We fix $M=9$, $\tau_{m}\in \{0.1,0.2,\ldots,0.9\}$, and use $\hat{r}$ as the number of factors for computation. For the initial guess $F^{0}$ and $\Lambda^{0}(\cdot)$ in Algorithm (ref) to compute UFA, we use the ones proposed in Remark (ref), whereas we use $\hat{F}$ and $\hat{\Lambda}(\cdot)$ obtained from UFA as the initial guess for Algorithm (ref).
In addition to the two estimators proposed in this paper, we also estimate the factor by PCA and by smoothed QFA chen2021quantile at each $\tau_{m}$ to see how weak mean/quantile factor affects the performance. For these two estimators, we directly use the true $r$ for the number of factors. For the smoothed QFA, we slightly modify the approach in chen2021quantile by using our smoothed objective function.
Note that since the smoothed QFA is also obtained from a nonconvex optimization problem without an analytical solution, no matter how smoothing is conducted, the initial guess should play a crucial role. We use two sets of initial guesses for it. The first one is the same as the one for our UFA; denoted by $\text{ini}_{UFA}$. This initial guess, by utilizing the information across all $M$ quantile levels, is already consistent in the average squared Frobenius norm. In practice, however, if a researcher is to use QFA, she is taking a quantile-level-specific approach, so it is more likely that she only solves the nuclear norm penalized estimator (ref) at one $\tau$, obtains $\hat{L}^{pel}(\tau)$, sets $\hat{F}^{pel}(\tau)$ equal to the eigenvectors corresponding to the $r$ largest eigenvalues of $\hat{L}^{pel'}(\tau)\hat{L}^{pel}(\tau)/NT$ multiplied by $\sqrt{T}$ and $\hat{\Lambda}^{pel}(\tau)=\hat{L}^{pel}(\tau)\hat{F}^{pel}(\tau)/T$, and uses $(\hat{F}^{pel}(\tau),\hat{\Lambda}^{pel}(\tau))$ as the initial guess, denoted by $\text{ini}_{\tau}$. We compute QFA under either initial guess.
To evaluate the performance, we regress the true factor $F_{0}^{*}$ on the estimated factors and obtain the adjusted $R^{2}$. A higher adjusted $R^{2}$ indicates that the estimated factor space is closer to the true one. Table (ref) shows the results averaged over 1000 simulation repetitions. Since our estimators and PCA are not quantile-level-specific, they have identical results across $\tau_{m}$.
Three observations can be drawn from Table (ref). First, the performance of UFA and IDW-UFA are very similar, confirming Theorems (ref) and (ref). The adjusted $R^{2}$, already high when $(N,T)=(50,50)$, increases as the sample size increases. Second, QFA, on the other hand, has inferior performance, especially when the factor becomes weaker around $\tau=0.5$. The initial guess is indeed important, but even under the already consistent initial guess $\text{ini}_{UFA}$, QFA still only has an adjusted $R^{2}$ of $0.51$ at $\tau=0.5$ under $(N,T)=(150,150)$. Under the initial guess $\text{ini}_{\tau}$ that is more likely to be adopted in practice, QFA no longer works at $\tau=0.5$ as the adjusted $R^{2}$ is close to 0. Third, PCA does not work due to the weak strength of $F_{0}^{*}$ as a mean factor.
We now demonstrate the quality of the normal approximation in Theorem (ref) for our IDW-UFA estimator. We follow the same data generating process as in Section (ref). However, now we only draw the factors and loadings once and keep them fixed across the 1000 simulation repetitions, the same as chen2021quantile. We estimate the fixed factors and loadings using IDW-UFA with sample size $(N,T)\in\{(50,50),(75,75),(100,100),(150,150)\}$. The implementation details are identical to Section (ref) except that now we use the true number of factors $r=1$ to avoid the rare cases of $\hat{r}\neq r$ where the dimensionality of $\tilde{f}_{t}$ is different from that of $f_{0,t}$.
Before presenting the results, first note that in Theorem (ref), $\tilde{f}_{t}$ is asymptotically normal up to rotation. Same as bai2003inferential, this rotation matrix $H_{NT,2}$ is unknown as it depends on the true $F_{0}$. However, under $r=1$, Theorem (ref) shows that $H_{NT,2}=\tilde{H}_{NT,1}+o_{p}(1/\sqrt{N})$, and thus by normalizing the sign of $\tilde{f}_{t}$ and $f_{0,t}$ which leads to $\tilde{H}_{NT,1}=1$, we can drop the rotation matrix and directly look at $\tilde{f}_{t}-f_{0,t}$. To verify, Table (ref) presents $|H_{NT,2}-1|$ averaged over 1000 repetitions under different sample size. It can be seen that the difference is indeed very small, and gets closer to 0 as the sample size increases.
Now, we consider the standardized estimated factor and common component. Define $\tilde{f}^{std}_{t}\coloneqq (\tilde{f}_{t}-f_{0,t})/se_{\tilde{f}_{t}}$, where $se_{\tilde{f}}$ is the standard error, constructed by the plug-in estimator described in Remark (ref). Note that it is not the standard deviation of the estimate across simulation repetitions. Under Theorem (ref), the distribution of $\tilde{f}^{std}_{t}$ approaches to standard normal as $N,T\to\infty$ for any fixed $t$. Define $\tilde{L}_{it}^{std}(\tau)$ similarly. Let $i=\lfloor N/2\rfloor$, $t=\lfloor T/2\rfloor$, where $\lfloor a \rfloor$ is the largest integer that is no greater than $a$. Table (ref) presents the sample mean and standard deviation of $\tilde{f}^{std}_{t}$ and $\tilde{L}^{std}_{it}(\tau)$ for $\tau=0.2,0.5,0.8$ over 1000 simulation repetitions. Note that the factor is relatively weak at $\tau=0.5$ whereas $0.2$ and $0.8$ are quantile levels near the boundaries. From the results, the sample mean is in general close to 0; it is relatively large at $(N,T)=(50,50)$, but soon gets smaller when the sample size reaches $(75,75)$. The standard deviation is close to 1, especially for $N>50$ and $T>50$.
Finally, we plot the histograms (scaled to be density functions) of these standardized estimates, superimposed by the standard normal density curve. Figure (ref) shows that the normal approximation is reasonably good in all cases, and performs better when the sample size increases.
This paper proposes a new factor analysis framework. By collecting all the factors that may have impact on the outcome, namely, the universal factors, our framework can induce an AFM and a QFM, but does not require the strong factor assumptions maintained in those models.
We build two estimators for the universal factors and loadings. Both estimators achieve the optimal rate when estimating the space spanned by the factors, regardless of the strength of the factors with respect to the mean or quantile of the outcome at a given quantile level. Built on a novel sample splitting strategy, the inverse density weighted estimator further achieves $\sqrt{N}$-asymptotic normality for individual factors and loadings. Monte Carlo shows that our estimators have superior performance to QFA and PCA in the presence of weak factors.
In addition, we propose estimators concerning the number of factors. A weak-factor-robust consistent estimator for the total number of universal factors is constructed. We also develop consistent factor selectors that can select out factors having any tolerated level of influence in an AFM or a QFM.