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.
84,289 characters · 16 sections · 84 citation commands
Regularized Quantile Regression with Interactive Fixed Effects
Panel data models are widely applied in economics and finance. Allowing for rich heterogeneity, interactive fixed effects are important components in such models in many applications. In applications such as asset pricing, it could be desirable to explain or forecast an outcome variable at certain quantile levels. However, the well studied mean regression with interactive fixed effects (e.g. pesaran2006estimation, bai2009panel and moon2015linear) misses such distributional heterogeneity.
In this paper, we consider a panel data model where the conditional quantile of an outcome variable is linear in the covariates and in the product of time- and individual-fixed effects. These fixed effects are unobservables that may be correlated with the covariates. The number of covariates is allowed to grow slowly to infinity with $N$ and $T$. Meanwhile, we allow the coefficients, the set of the effective fixed effects and the realization of each of them to all be quantile level dependent, generating large modeling flexibility.
To estimate the model, this paper proposes a nuclear norm penalized estimator. By deriving its theoretical error bound, we show that the estimator can consistently estimate the coefficients and the (realizations of) the interactive fixed effects uniformly in quantile level. The estimator solves a convex problem and computes fast in practice even in large panel data sets based on our proposed augmented Lagrangian multiplier algorithm. Implementing the estimator does not require pre-estimation of the number of fixed effects or the fixed effects themselves.
To illustrate the estimator, let us consider a simple example where only the coefficients are quantile level dependent: For a panel data set $\{(Y_{it},X_{it}):i=1,...,N;t=1,...,T\}$, suppose the $u$-th conditional quantile of $Y_{it}$ given $p$ covariates $X_{it}$ and $r$ time- and individual-fixed effects $(F_{t},\Lambda_{i})$ is $q_{Y_{it}|X_{it},F_{t},\Lambda_{i}}(u)\coloneqq X_{it}'\beta_{0}(u)+F_{t}'\Lambda_{i}$ where both $F_{t}$ and $\Lambda_{i}$ are $r\times 1$ vectors. The fixed effects form an $N\times T$ matrix $L_{0}\coloneqq (\Lambda_{1},...,\Lambda_{N})'(F_{1},...,F_{T})$ whose rank is at most $r$. Thus, the dense matrix $L_{0}$ is low-rank when $r$ is small relative to $N$ and $T$. Exploiting such low-rankness, our estimator, inspired by the seminal work by candes2009exact, jointly estimates $(\beta_{0}(u),L_{0})$ for a given quantile level $u\in (0,1)$ by solving
where $\rho_{u}$ is the standard check function in the quantile regression literature, $\lambda$ is the positive penalty coefficient, $||L||_{*}$ is the nuclear norm of the $N\times T$ matrix $L$, and $\mathcal{B}$ and $\mathcal{L}$ are convex parameter spaces about which we will be specific later.
The key component of the estimator is the convex nuclear norm penalty. Summing the singular values of a matrix, the nuclear norm is to the rank, counting the nonzero singular values, what the convex $\ell_{1}$-norm is to the nonconvex $\ell_{0}$-norm of the vector of the singular values. Hence, the nuclear norm penalty can be viewed as the matrix counterpart of the LASSO penalty in regression with high-dimensional regressors. Being a convex surrogate of the rank functional, we show that this penalty is effective to deliver consistent estimates under low-rankness of $L_{0}$.
The main benefits of setting up the minimization problem as in (ref) are that the objective function is convex in $(\beta,L)$ and that one does not need to know $r$ before implementation. To highlight these benefits, let us consider a natural alternative estimator
ando2019quantile study a similar estimator where the coefficients are $i$-specific. The objective function in (ref) is nonconvex in the parameters $(\beta,\Lambda,F)$. Due to nonsmoothness of the check function $\rho_{u}$, there does not exist a closed form solution to the minimization problem (ref). Implementation needs to be carried out iteratively. Then, nonconvexity leads to two potential issues. First, one may obtain a local minimum which can be arbitrarily far from the global one. Second, solving the minimization problem may be computationally intensive, especially in large panel data sets. In the simulation experiments in the paper, we will see that the penalized estimator we propose indeed outperforms this alternative estimator in the computation aspect as it saves computation time by much. Meanwhile, to make the alternative estimator (ref) feasible, $r$ needs to be known or pre-estimated. This step results in additional computation burden and a mis-specified $r$ may lead to inconsistent estimates.
Although as we show in this paper, the nuclear norm penalty slows down the rate of convergence of the coefficient estimator when the number of regressors is small, we find from the simulation experiments that its finite sample bias is actually comparable with the alternative estimator (ref) even though we use the true $r$ for the latter. Moreover, we propose a consistent estimator of $r$ based on our penalized estimator. With this rank estimator, treating our consistent penalized estimator as an initial value for the iterative estimator (ref) may potentially avoid the local-minimum problem, reduce computation time, and remove the bias in the penalized estimator. Finally, although the coefficient estimator may have a slower rate of convergence, we show that the error bound on the estimator of $L_{0}$ can be nearly optimal in squared Frobenius norm, not affected by penalization.
With the dense latent component and the nonsmooth objective function involved, deriving the estimator's uniform error bound is challenging. We prove new results on random matrices for this purpose. Moreover, we develop novel theoretical arguments which relax some usually made assumptions or replace some high-level technical conditions in the panel data quantile regression literature with primitive ones that are easier to interpret. These results will be introduced later in related sections and may be of independent interest.
This paper adds to the literature of panel data quantile regression. Since koenker2004quantile, panel data quantile regression began to draw increasing attention. abrevaya2008effects, lamarche2010robust, canay2011simple, kato2012asymptotics, galvao2013estimation and galvao2016smoothed study quantile regression with one-way or two-way fixed effects. harding2014estimating consider interactive fixed effects with endogenous regressors. They require the factors to be pre-estimated or known. chen2019fixed considers quantile regression with interactive fixed effects. They need to first estimate the time fixed effects, or, the factors, that are assumed to be quantile-level-independent. They then estimate the coefficients and the individual fixed effects via smoothed quantile regression. chen2019quantile propose a quantile factor model without regressors. They estimate the factors and the factor loadings via nonconvex minimization similar to (ref). Pre-estimation of the number of factors is needed. ando2019quantile consider quantile regression with heterogeneous coefficients and a factor structure. They propose both a frequentist and a Bayesian estimation procedure. The number of factors also needs to be estimated first, and the minimization problem is noncovex. Both chen2019quantile and ando2019quantile establish consistency pointwise in quantile level, while we focus on uniform consistency. On the technical side, both impose stronger assumptions on the conditonal density of the outcome variable than our paper\footnote{We will discuss these differences in detail in Appendix (ref).}. In our simulation study in Section (ref), we compare our estimator with ando2019quantile and find the estimates are comparable while our procedure is computationally more efficient.
Another literature this paper speaks to is on nuclear norm penalized estimation. This literature was initially motivated by low-rank matrix completion or recovery problems in computer science and statistics (e.g. candes2009exact, ganesh2010dense, zhou2010stable, candes2011robust, hsu2011robust, negahban2011estimation, agarwal2012noisy and negahban2009unified among others). In this literature, the outcome matrix is usually modeled as the sum of a low-rank matrix and some other matrices that are for instance, sparse or Gaussian. The primary goal is to estimate the low-rank or the sparse matrix. This setup is different from our paper. Nuclear norm penalized estimation and matrix completion related topics have also gained interest in econometrics recently. athey2018matrix, moon2019nuclear\footnote{moon2019nuclear also briefly discuss nuclear norm penalized quantile regression with a single regressor as an extension. Using a different approach than this paper, they focus on pointwise (in quantile level) convergence rate of the coefficient estimator. In this paper, we obtain uniform rates for both the coefficients and the low-rank component. Also, the number of covariates can be more than one and growing to infinity slowly.} and chernozhukov2018inference investigate nuclear norm penalized mean regression with interactive fixed effects. beyhum2019square also consider mean regression with interactive fixed effects but they use a square-root nuclear norm penalty. bai2019robust propose a nuclear norm regularized median regression for robust principal component anaylsis for fat tailed data. bai2019matrix consider imputation of missing data and counterfactuals. bai2019rank study penalized estimation for approximate factor models with singular values thresholding. chao2015factorisable consider penalized multi-task quantile regression where there are multiple outcome variables and the coefficient matrix is low-rank. ma2020detecting apply nuclear norm penalized logistic regression to the study of an undirected network formation model.
A recent paper by belloni2019high studies quantile regression with both interactive fixed effects and high-dimensional regressors. This work was done in parallel and our paper is independent of it. In that paper, they use a nuclear norm constraint for the low-rank matrix and an additional $\ell_{1}$-norm constraint on the coefficients to deal with the high-dimensional regressors. In contrast, we focus on low-dimensional regressors, although we do allow the number of regressors to slowly grow to infinity. On the other hand, unlike our paper that derives a uniform error bound, they focus on convergence rate pointwise in quantile level. As aforementioned, achieving uniformity is challenging and requires the new results on random matrices developed in this paper. Moreover, some of our assumptions are weaker or more primitive. We defer a detailed discussion on this aspect until Appendix (ref). On the computation side, our algorithm differs from theirs and works fast in our simulation experiments. We view these two papers as complementary.
The rest of the paper is organized as follows. Section (ref) introduces the model and the estimator. Section (ref) previews the main results and provides a proof sketch to highlight the challenges and the key theoretical arguments we develop. Section (ref) discusses the main results and their consequences. Section (ref) demonstrates a Monte Carlo simulation study by comparing our estimator with two alternative approaches. Section (ref) concludes. The algorithm and implementation details are in Appendix (ref). An alternative approach to proving consistency and a comparison of the assumptions in this paper with the most related literature are in Appendix (ref). Appendix (ref) collects all the proofs. Appendix (ref) presents some additional simulation results.
Besides the nuclear norm $||\cdot||_{*}$, four additional matrix norms are used in the paper: Let $||\cdot||$, $||\cdot||_{F}$, $||\cdot||_{1}$, and $||\cdot||_{\infty}$ denote the spectral norm, the Frobenius norm, the $\ell_{1}$-norm and the maximum norm. When applied to a vector, the Frobenius norm is equal to the Euclidean norm. For two generic $N\times T$ matrices $A$ and $B$, $\left\langle A, B\right\rangle\coloneqq \sum_{i,t}A_{it}B_{it}$ denotes the inner product of $A$ and $B$. For two generic real numbers, $a\lor b$ and $a\land b$ return the maximum and the minimum of $a$ and $b$, respectively.
We consider a panel data set $\{(Y_{it},X_{it}):i=1,...,N;t=1,...,T\}$ where $Y_{it}$ is a scalar outcome and $X_{it}$ is a $p\times 1$ vector of covariates. Let $Y=(Y_{it})_{i,t}$ and $X_{j}=(X_{j,it})_{i,t}$ ($j=1,...,p$) be $N\times T$ matrices of the outcome and the $j$-th covariate. Let $\mathcal{U}$ be a compact subset of $(0,1)$. For any $u\in\mathcal{U}$, there are $\bar{r}$ possibly $u$-dependent time- and individual-fixed effects. For $k=1,...,\bar{r}$, let $F_{k}(u)=(F_{1k}(u),...,F_{Tk}(u))'$ be the $k$-th time fixed effect. Let $\Lambda_{k}(u)=(\Lambda_{1k}(u),...,\Lambda_{Nk}(u))'$ be the $k$-th individual fixed effects. Let $W=(X_{1},...,X_{p},\{F_{k}(u)\}_{k=1,...,\bar{r},u\in\mathcal{U}},\{\Lambda_{k}(u)\}_{k=1,...,\bar{r},u\in\mathcal{U}})$. Assume for all $u\in\mathcal{U}$, the conditional quantile of outcome $Y_{it}$ in matrix notation satisfies the following model with probability one:
where $\mathbbm{1}_{k}(u)\in\{0,1\}$ determines whether the $k$-th fixed effect $F_{k}(u)$ or $\Lambda_{k}(u)$ affects the $u$-th conditional quantile of $Y$ at all. Model (ref) allows both the set of the effective fixed effects and the realizations of them to depend on $u$. Throughout, we allow the fixed effects to be either random or deterministic. The covariates can be correlated with them when they are random. Similar setups of the fixed effects or factor structures in panel data quantile regression can be found in ando2019quantile and chen2019quantile.
The fixed effects in equation (ref) form an $N\times T$ matrix $L_{0}(u)\coloneqq \sum_{k=1}^{\bar{r}}\mathbbm{1}_{k}(u)\Lambda_{k}(u)F_{k}(u)'$. The rank of the matrix $L_{0}(u)$ is at most $r(u)\coloneqq\sum_{k=1}^{\bar{r}}\mathbbm{1}_{k}(u)\leq \bar{r}$ by construction. In this paper, we assume $\bar{r}$ is fixed. Thus, $L_{0}(u)$ is low-rank when $N$ and $T$ are large.
Let $\beta_{0}(u)=(\beta_{0,j}(u))_{j=1,...,p}$. This paper focuses on consistently estimating $(\beta_{0}(u),L_{0}(u))$ uniformly in $u\in\mathcal{U}$. When $L_{0}(u)$ is random, consistency is in terms of its realization.
Now let us present a few models which admit the conditional quantile function (ref).
Now we introduce our estimator of $(\beta_{0}(u),L_{0}(u))$. For a generic $N\times T$ matrix $Z$, define $\bm{\rho}_{u}(Z)\coloneqq \sum_{i,t}\rho_{u}(Z_{it})\equiv \sum_{i,t} Z_{it}(u-\mathbbm{1}(Z_{it}\leq 0))$. By exploiting the linearity of the conditional quantile function (ref) in $(\beta_{0}(u),L_{0}(u))$ and the low-rankness of $L_{0}(u)$, this paper proposes the following nuclear norm penalized quantile regression estimator to jointly estimate $\beta_{0}(u)$ and $L_{0}(u)$ for any $u\in\mathcal{U}$:
where $\lambda>0$. The parameter space $\mathcal{L}\coloneqq\{L\in\mathbb{R}^{N\times T}:||L||_{\infty}\leq \alpha_{NT}\}$ is convex and compact and $\alpha_{NT}\geq 1$ can be $(N,T)$-dependent. In particular, we allow $\alpha_{NT}$ to grow to infinity with $N$ and $T$. We need $\alpha_{NT}$ for technical reasons to be discussed in Section (ref). In Appendix (ref), we show that we can drop $\alpha_{NT}$ to make $\mathcal{L}=\mathbb{R}^{N\times T}$ under a different set of assumptions.
Two remarks on the estimator are in order .
In this section, we first briefly summarize the main theoretical results of this paper, then outline the proof strategy to highlight the challenges arising from the nonsmooth objective function and the dense low-rank common component $L_{0}(u)$. We also introduce and discuss the assumptions we make motivated by these challenges.
Under the assumptions to be introduced in this section, this paper shows that for some universal constant $C_{error}>0$, the estimator $(\hat{\beta}(u),\hat{L}(u))$ defined in equation (ref) satisfies the following inequality with probability approaching one (w.p.a.1):
where
The error bound implies uniform consistency of the estimator given a fixed $\bar{r}$ and a fixed or slowly growing $p$ and $\alpha_{NT}$. Based on the order of the error bound, we also propose a consistent estimator of the number of the effective fixed effects $r(u)$ for each $u\in\mathcal{U}$. We will discuss these result with more details in Section (ref), but now let us first sketch how the error bound is derived.
To derive the error bound, it is helpful to first exploit some simple implications from the definition of the estimator to sharpen the space where the estimation errors $(\hat{\Delta}_{\beta}(u),\hat{\Delta}_{L}(u))\coloneqq (\hat{\beta}(u)-\beta_{0}(u),\hat{L}(u)-L_{0}(u))$ lie so that the analysis can be conducted in this smaller space instead of $\mathbb{R}^{p}\times \mathbb{R}^{N\times T}$. For this purpose, we make the following assumption.
The independence requirement in part i) is for simplicity so that some inequalities for random matrices can be easily applied in the proof. The same assumption when $\mathcal{U}$ is a singleton so that $W$ only contains the covariates and the fixed effects at $u$ can be found in ando2019quantile and chen2019quantile as well. Moderate serial correlation in $V_{it}(u)$ can be allowed at a cost of more technical conditions. On the other hand, serial correlation in the covariates is indeed allowed, and $p$ is allowed to be growing in $N$ and $T$; part ii) in Assumption (ref) holds as long as for instance, for every $j\leq p$, Chebyshev's inequality holds for $\sum_{i,t}X_{j,it}^{2}/NT$ and its variance multiplied by $p$ is $o(1)$.
It turns out that under Assumption (ref) alone, the estimation error $(\hat{\Delta}_{\beta}(u),\hat{\Delta}_{L}(u))$ lies in a cone uniformly in $u\in\mathcal{U}$ w.p.a.1. The cone has nice properties for deriving the error bound. To characterize the cone, let us introduce some notation. Let $R(u)\Sigma(u) S(u)'$ be a singular value decomposition of $L_{0}(u)$. Following candes2009exact, let $\Phi(u)\coloneqq\{M\in \mathbb{R}^{N\times T}:\exists A\in\mathbb{R}^{r(u)\times T}\ \text{and}\ B\in\mathbb{R}^{N\times r(u)}\ s.t.\ M=R(u)A+BS(u)'\}$. Denote the orthogonal projection of a generic $N\times T$ matrix $W$ onto this space by $\mathcal{P}_{\Phi(u)}W$, then
We then have the following lemma.
By Lemma (ref), define cone $\mathcal{R}_{u}$ by
We thus have $(\hat{\Delta}_{\beta}(u),\hat{\Delta}_{L}(u))\in\mathcal{R}_{u}$ uniformly in $u\in\mathcal{U}$ w.p.a.1 under Assumption (ref) and under $||L_{0}(u)||_{\infty}\leq \alpha_{NT}$ for all $u\in\mathcal{U}$.
The key property of $\mathcal{R}_{u}$ is that for any element $(\Delta_{\beta},\Delta_{L})\in\mathcal{R}_{u}$ and for any $u\in\mathcal{U}$, $||\Delta_{L}||_{*}$ and $||\Delta_{L}||_{F}$ can be of the same order, a property that low-rank matrices also share\footnote{Note that this property is non-trivial because in general the nuclear norm $||\Delta_{L}||_{*}$ can be as large as $\sqrt{N\land T}||\Delta_{L}||_{F}$. }. To see why this is true, note that by equation (ref), the rank of $\mathcal{P}_{\Phi(u)}\Delta_{L}$ is at most $3r(u)$ for all $u$. Hence, for all $u\in\mathcal{U}$, \[||\mathcal{P}_{\Phi(u)}\Delta_{L}||_{*}\leq \sqrt{3r(u)}||\mathcal{P}_{\Phi(u)}\Delta_{L}||_{F}\leq \sqrt{3r(u)} ||\Delta_{L}||_{F}\leq \sqrt{3r(u)}||\Delta_{L}||_{*}, \] where the first and the last inequalities are by the relationship between the nuclear norm and the Frobenius norm. The second inequality is due to $\left\langle\mathcal{P}_{\Phi(u)}\Delta_{L},\Delta_{L}-\mathcal{P}_{\Phi(u)}\Delta_{L}\right\rangle=0$ and by the Pythagoras formula. As a consequence, elements in $\mathcal{R}_{u}$ satisfy
which implies that if $C_{Cone}\sqrt{p(N\land T)\log(pNT)}||\Delta_{\beta}||_{F}/||\Delta_{L}||_{*}=o(1)$, $||\Delta_{L}||_{*}$ and $||\Delta_{L}||_{F}$ are of the same order. Note that since the estimation error $(\hat{\Delta}_{\beta}(u),\hat{\Delta}_{L}(u))$ is in $\mathcal{R}_{u}$ w.p.a.1, $||\hat{\Delta}_{L}(u)||_{F}$ and $||\hat{\Delta}_{L}(u)||_{*}$ are also of the same order w.p.a.1. This implies that, although it is unclear whether the nuclear norm penalty makes some of the singular values of $\hat{\Delta}_{L}(u)$ to be exact zero so that $\hat{\Delta}_{L}(u)$, and in turn, $\hat{L}(u)$, is low-rank, at least the singular values of $\hat{\Delta}_{L}(u)$ must not be of the same order w.p.a.1.
In the sense that the nuclear norm and the Frobenius norm of the matrix elements in $\mathcal{R}_{u}$ can be of the same order, the cone $\mathcal{R}_{u}$ in Lemma (ref) matches those obtained in the broad literature of nuclear norm penalized estimation under different objective functions (see e.g. agarwal2012noisy, negahban2012restricted, athey2018matrix, chernozhukov2018inference). Different from the mentioned literature, to establish uniformity in Lemma (ref) under the nonsmooth objective function, we derive new uniform bounds on some norms of random matrices whose entries are jump processes (see Lemma (ref) in Appendix (ref) for details). These results may be of independent interest.
Under Lemma (ref), we can conduct all the subsequent analysis conditional on the event that $(\hat{\Delta}_{\beta}(u),\hat{\Delta}_{L}(u))\in\mathcal{R}_{u}$ for all $u\in\mathcal{U}$ to exploit property (ref) of $\mathcal{R}_{u}$. Now let us sketch the derivation of the error bound (ref) to motivate the other assumptions we make, introduce the theoretical challenges, and discuss how we overcome them.
Let $\mathcal{D}\coloneqq \mathbb{R}^{p}\times \{\Delta_{L}\in\mathbb{R}^{N\times T}:||\Delta_{L}||_{\infty}\leq 2\alpha_{NT} \}$. Recall that $\mathcal{L}=\{L\in\mathbb{R}^{N\times T}:||L||_{\infty}\leq\alpha_{NT}\}$. If $||L_{0}(u)||_{\infty}\leq \alpha_{NT}$ for all $u\in\mathcal{U}$ w.p.a.1, then by $\hat{L}(u)\in\mathcal{L}$ and by Lemma (ref), we have $(\hat{\Delta}_{\beta}(u),\hat{\Delta}_{L}(u))\in \mathcal{R}_{u}\cap\mathcal{D}$ for all $u\in\mathcal{U}$ w.p.a.1. For a generic sample $(Z_{it})_{i,t}$ and a function $f$, let $\mathbb{G}_{u}(f(Z_{it}))\coloneqq \sum_{i,t}[f(Z_{it})-\mathbb{E}(f(Z_{it})|W)]/\sqrt{NT}$ where recall $W\coloneqq (X_{1},...,X_{p},\{F_{k}(u)\}_{k=1,...,\bar{r},u\in\mathcal{U}},\{\Lambda_{k}(u)\}_{k=1,...,\bar{r},u\in\mathcal{U}})$. Since $\mathcal{R}_{u}$ is a cone and $\mathcal{D}$ is convex and contains zero, by convexity of the objective function in equation (ref), it is sufficient to show that it is zero probability to have a $u\in\mathcal{U}$ such that:
where $(\Delta_{\beta},\Delta_{L})$ is a generic element in $\mathcal{R}_{u}\cap\mathcal{D}$ and recall that $V(u)\coloneqq Y-\sum_{j=1}^{p}X_{j}\beta_{0,j}(u)-L_{0}(u)$.
We prove inequality (ref) can not hold for any $u\in\mathcal{U}$ in two steps:
After we find these bounds, we conclude that inequality (ref) is impossible uniformly in $u\in \mathcal{U}$ if the lower bound in (ref) is always greater in than the upper bound in (ref).
In (ref), we are to lower bound the conditional expectation by $(||\Delta_{\beta}||_{F}^{2}+||\Delta_{L}||_{F}^{2}/NT)$ up to some multiplicative factor. The major theoretical difficulty arises from high-dimensionality of $\Delta_{L}$. Note that (ref) is independent of the nuclear norm penalty and a similar difficulty exists even if we maintain the factor structure (e.g. ando2019quantile and chen2019quantile) instead of treating the interactive fixed effects as a single low-rank matrix. Therefore, the theoretical arguments developed in the paper to overcome the difficulty can be applied to other estimators for panel data quantile regression with interactive fixed effects or factor structures.
To illustrate the difficulty, consider a simple case without covariates. The conditional expectation under consideration can then be simplified as $\mathbb{E}\left[\bm{\rho}_{u} \left(V(u)-\Delta_{L}\right)-\bm{\rho}_{u}(V(u))|W\right]$. By Knight's identity knight1998limiting and by the definition of $V_{it}(u)$, it can be rewritten as
where $F_{V_{it}(u)|W}$ is the conditional cumulative distribution function of $V_{it}(u)$. The key problem is as follows. The magnitude of some $\Delta_{L,it}$s can be arbitrarily large or grow to infinity with $N$ and $T$ even when $||\Delta_{L}||_{F}^{2}/NT\to 0$. Hence, if one adopts the standard argument in quantile regression to first-order Taylor expand $F_{V_{it}(u)|W}(s)$ around $0$, it is insufficient to obtain a strictly positive lower bound on quantity (ref) if one only assumes that the conditional density of $V_{it}(u)$ at zero, $f_{V_{it}(u)|W}(0)$, is continuous and positive uniformly in $i$, $t$ and $u$ almost surely.
To overcome the difficulty, the literature on panel data quantile regression with factor structures or interactive fixed effects often assumes that elements in $L_{0}(u)$ (or the fixed effects) lie in a fixed compact space and the conditional density $f_{V_{it}(u)|W}(s)$ is bounded away from 0 on any compact interval of $s$ (e.g. ando2019quantile and chen2019quantile). Alternatively, belloni2019high adopt higher order Taylor expansion by assuming differentiability of the conditional density with bounded derivatives. They also impose a high-level condition on the estimation error matrix which is not straightforward to interpret\footnote{A more detailed comparison between these approaches and ours is in Appendix (ref)}.
In this paper, we develop two novel sets of theoretical arguments, or, approaches, to overcome the difficulty to achieve a positive quadratic lower bound. Under both approaches, $f_{V_{it}(u)|W}$ bounded away from zero in an arbitrarily small neighborhood of $0$ is sufficient. Meanwhile, we do not need the conditional density function to be differentiable. Formally, we make the following assumption:
Note that $\delta$ can be arbitrarily small. A sufficient condition for Assumption (ref) to hold is that $f_{V_{it}(u)|W}(0)>0$ uniformly in $i,t$ and $u\in\mathcal{U}$ and the functions $\{f_{V_{it}(u)|W}\}_{i,t,u}$ are equicontinuous at $0$ for all realizations of $W$. This assumption can be shown to be weaker than the assumptions on the conditional density in belloni2019high, ando2019quantile and chen2019quantile. See Appendix (ref) for a detailed discussion.
Under Assumption (ref), our first approach requires a compact parameter space $\mathcal{L}$ for the matrix component as in equation (ref), but unlike ando2019quantile and chen2019quantile, the boundary $\alpha_{NT}$ is not fixed and is allowed to grow to infinity. Our second approach relaxes $\mathcal{L}$ to be $\mathbb{R}^{N\times T}$ at the cost of more technical conditions.
We present the main idea behind our first approach now and discuss our second approach in Appendix (ref). Since we are considering a case without the covariates, let us redefine $\mathcal{D}=\{\Delta_{L}\in\mathbb{R}^{N\times T}:||\Delta_{L}||_{\infty}\leq 2\alpha_{NT}\}$. We show in Appendix (ref) that each integral in the summation (ref) is decreasing in the absolute upper limit. So each integral satisfies
where $\delta$ is the constant in Assumption (ref). For the integral on the right side, even if $\Delta_{L,it}$ is diverging, $|\Delta_{L,it}|/2\alpha_{NT}\leq 1$ so $(1\land \delta)\Delta_{L,it}/2\alpha_{NT}$ must lie in the region where the conditional density is positive. Then first-order Taylor expanding $F_{V_{it}(u)|W}(s)$ around $0$ yields a strictly positive quadratic lower bound. The issue is thus resolved.
When there are covariates, the upper limits in the integrals in quantity (ref) become $\Delta_{L,it}+X_{it}'\Delta_{\beta}$. To make this approach still valid at least in the ball $||\Delta_{\beta}||_{F}^{2}+||\Delta_{L}||_{F}^{2}/NT\leq \gamma^{2}$, we make the following assumption so that $X_{it}'\Delta_{\beta}$ is also bounded by some function of $\alpha_{NT}$ w.p.a.1.
Assumption (ref) is reasonably mild because we allow $\alpha_{NT}$ to be $(N,T)$-dependent and to grow to infinity. For instance, if $\alpha_{NT}$ has order $\log(NT)$, part i) is satisfied in all the models in Examples (ref) to (ref) if all the individual- and time-fixed effects are sub-Gaussian with $\bar{r}$ fixed. Note that this does not rule out the case where $Y_{it}$ itself has a heavy tail. The restriction on the covariates' tails is even milder. For instance, if $\alpha_{NT}=\log(NT)$ and $p$ is fixed, given the order of $\gamma$ defined in equation (ref), some heavy-tailed distributions are allowed for the covariates.
With the covariates, recall that $\mathcal{D}\coloneqq \mathbb{R}^{p}\times \{\Delta_{L}\in\mathbb{R}^{N\times T}:||\Delta_{L}||_{\infty}\leq 2\alpha_{NT} \}$. We have the following lemma.
Finally, in order to obtain the error bounds on $\hat{\beta}(u)$ and $\hat{L}(u)$ separately, we impose the following identification assumption.
Part i) in Assumption (ref) guarantees that the individual coefficients $\beta_{0,j}(u)$s can be separately consistently estimated. Part ii) is a version of the widely adopted restricted strong convexity condition in the literature on low-rank matrix recovery or nuclear norm penalized estimation (e.g. agarwal2012noisy, negahban2011estimation,negahban2012restricted, negahban2009unified, belloni2019high, chernozhukov2018inference, etc.). It says that any linear combinations of the covariate matrices $X_{j}$s must lie sufficiently far from the matrix elements in the cone $\mathcal{R}_{u}$. The constant $C_{RSC}$ is determined by the joint distribution of the covariates and $L_{0}(u)$. When the covariates are more correlated with $L_{0}(u)$, this constant tends to become smaller. Consequently, the error bound on the estimation error would be larger. We can observe this pattern in the Monte Carlo simulations in Section (ref) and Appendix (ref) later. Restricted strong convexity can be interpreted as an identification condition. To see this, note that $(\Delta_{\beta},L_{0}(u))\in\mathcal{R}_{u}$ for all $\Delta_{\beta}\in\mathbb{R}^{p}$ because $\mathcal{P}_{\Phi(u)}L_{0}(u)=L_{0}(u)$ by construction. Thus by letting $\Delta_{L}=L_{0}(u)$, condition (ref) implies that $L_{0}(u)$ is not equal to any linear combination of the $X_{j}$s, a necessary condition in order to identify $L_{0}(u)$.
Under Assumption (ref), the bound obtained in Lemma (ref) can be shown to be further lower bounded by $(C_{\sigma}\sigma_{min}^{2}\land 1)C_{min}C_{RSC}(||\Delta_{\beta}||_{F}^{2}+||\Delta_{L}||_{F}^{2}/NT)/\alpha_{NT}^{2}$ for some $C_{\sigma}>0$ w.p.a.1 (see the proof of Theorem (ref) for details). (ref) is thus completed.
The key property we use to obtain a tight enough upper bound on the absolute empirical process and the penalty difference in equation (ref) is equation (ref) implied by Lemma (ref).
For illustrative purpose, let us again assume there are no covariates, so equation (ref) implies that $||\Delta_{L}||_{*}$ has the same order as $||\Delta_{L}||_{F}$, which, in the ball $||\Delta_{L}||_{F}^{2}/NT\leq\gamma^{2}$ we consider, is at most $\sqrt{NT}\gamma$. For the absolute penalty difference in equation (ref), by triangle inequality, $|||L_{0}(u)+\Delta_{L}||_{*}-||L_{0}(u)||_{*}|$ is upper bounded by $||\Delta_{L}||_{*}\leq 4\sqrt{3r(u)}||\Delta_{L}||_{F}\leq 4\sqrt{3r(u)}\sqrt{NT}\gamma$ by property (ref). Note that without the property, the upper bound would be $\sqrt{N\land T}\sqrt{NT}\gamma$ instead, much greater than the current one as $r(u)$ is fixed.
Property (ref) also helps to obtain a tight upper bound on the absolute process $|\mathbb{G}_{u}|$ in equation (ref). Let $A$ denote a generic $N\times T$ matrix. When deriving an upper bound on $|\mathbb{G}_{u}|$, we often need to upper bound terms in the form of $|\sum_{i,t}A_{it}\Delta_{L,it}|$, or equivalently, $|\langle A,\Delta_{L}\rangle|$. There are at least two possible upper bounds for this inner product:
where inequality (ref) is by Cauchy-Schwartz and inequality (ref) is by Lemma 3.2 in candes2009exact. In general, these two upper bounds can be of the same order; although $||A||$ can be as small as $||A||_{F}/\sqrt{N\land T}$, $||\Delta_{L}||_{*}$ can be as large as $\sqrt{N\land T}||\Delta_{L}||_{F}$ without any restrictions. However, now in the ball $||\Delta_{L}||_{F}\leq \sqrt{NT}\gamma$, under equation (ref) without covariates, $||\Delta_{L}||_{*}$ can be bounded by $\sqrt{4\bar{r}}\sqrt{NT}\gamma$. As $\bar{r}$ is fixed, the order of the right side of inequality (ref) can then be only $1/(\sqrt{N\land T})$ of that of inequality (ref), providing a tight enough upper bound on the inner product.
Finally, in order to bound the absolute process and the penalty difference uniformly in $u\in\mathcal{U}$, we make the following assumption on smoothness in our estimands. This assumption is trivially satisfied if $\mathcal{U}$ is a singleton.
The smoothness assumption for $\beta_{0}(\cdot)$, equation (ref), is the same as in belloni2011l1. The smoothness condition on $L_{0}(\cdot)$, equation (ref), is a matrix counterpart. One can verify that $L_{0}(\cdot)$ in Examples (ref) satisfies condition (ref) if the quantile function of $\epsilon$, $q_{\epsilon}(\cdot)$, is Lipschitz continuous. For Example (ref), $L_{0}(\cdot)$ satisfies condition (ref) if $q_{\epsilon}(\cdot)$ is Lipschitz continuous and if $\sum_{i,t}(\sum_{m=1}^{\bar{r}_{2}}F_{tk}^{b}\Lambda^{b}_{ik})/NT$ converges in probability to a constant. Note that condition (ref) rules out some cases where the set of the effective fixed effects changes on $\mathcal{U}$. To see this, suppose in Example (ref), there exists a jumping point $u_{0}\in \mathcal{U}$ such that $r(u)<r(u')$ for any $u<u_{0}\leq u'$. Suppose the individual- and the time-fixed effects do not depend on $u$ and the first $r(u)$ of them at $u'$ are the same as those at $u$. Then \[\frac{1}{\sqrt{NT}}||L_{0}(u')-L_{0}(u)||_{F}=\sqrt{\frac{1}{NT}\sum_{i,t}\left(\sum_{k=r(u)+1}^{r(u')}F_{tk}\Lambda_{ik}\right)^{2}} \] which may converge in probability to a positive constant if the law of large number holds for it. Nevertheless, this assumption is not restrictive even in this situation when there are only a finite number of such jumping points in $(0,1)$: Let $\{u_{0,k}:k=1,2,...,K\}$ $(K<\infty)$ be the set of such jumping points with $u_{0,k}<u_{0,k+1}$ for all $k=1,...,K-1$. If Assumption (ref) holds for compact interval $\mathcal{U}_{k}\subset (u_{0,k},u_{0,k+1})$ for each $k$, we can then establish uniform error bound over each $\mathcal{U}_{k}$ and the uniform bound over $\bigcup_{k=1}^{K}\mathcal{U}_{k}$ is immediately obtained.
We have the following lemma bounding the absolute process in equation (ref). The bound on the penalty difference is straightforward and is thus omitted here.
In this section, we formally present the main result on the error bound on $(\hat{\beta}(u),\hat{L}(u))$. We also derive a consistent estimator of the number of fixed effects $r(u)$, i.e., the rank of $L_{0}(u)$.
Under the assumptions proposed in Section (ref), we have the following theorem.
Theorem (ref) implies uniform consistency of our estimator of $\beta_{0}(u)$ and $L_{0}(u)$ over $\mathcal{U}$ for a fixed $\bar{r}$ and slowly growing $\alpha_{NT}$ and $p$. The key determinants of the rate of convergenge are $p\log(pNT)/NT$ and $1/(N\land T)$. This part matches the nuclear norm penalized mean regression literature (athey2018matrix, moon2019nuclear and chernozhukov2018inference). We will discuss this part with more details in this section. Other determinants of the error bound include the number of fixed effects ($\bar{r}$), the quality of lower-bounding the conditional expectation in (ref) ($\alpha_{NT}^{2}$ and $\underline{f}$), and the strength of identification ($\sigma_{min}^{2}$ and $C_{RSC}$). The way that these parameters affect the error bound is expected: A larger $\bar{r}$ results in a relatively higher-rank common component, making the estimation problem more difficult. A greater quadratic lower bound (a greater $\underline{f}$ or a smaller $\alpha_{NT}$) implies that the objective function tends to be more sensitive to perturbations around the global minimum. Finally, stronger identification (a larger $C_{RSC}$ or a larger $\sigma_{min}^{2}$) makes it easier to separate the estimation errors $\hat{\Delta}_{\beta,j}(u)$s and $\hat{\Delta}_{L}(u)$.
agarwal2012noisy provide a minimax result that implies the (near) optimality of the error bound on $\hat{L}(u)$ obtained in Theorem (ref) when $\alpha_{NT}=\log(NT)$ and $p\log(pNT)/NT=O(\bar{r}/(N\land T))$. In agarwal2012noisy, they study a model where an observable $N\times T$ matrix $Y$ is the sum of a low-rank matrix, a sparse matrix and a noise matrix of i.i.d. Gaussian entries. To apply their result to our model, consider a special case where $u=0.5$ and $Y=L_{0}+V$ where $V$ is an $N\times T$ matrix of i.i.d. $N(0,v^{2})$ entries. This model both satisfies the conditional quantile model (ref) studied in this paper at $u=0.5$ with $\beta(0.5)=0$ and $q_{Y|L_{0}}(0.5)=L_{0}$, and also satisfies their setup with the sparse component being exactly zero. Their Theorem 2 (p.1195) shows that the lower bound on the minimax risk in the squared Frobenius norm over the family $\{L_{0}:\text{rank}(L_{0})\leq \bar{r},||L_{0}||_{\infty}\leq\alpha_{NT}\}$ has the order $\bar{r}/(N\land T)$. Comparing the order of the lower bound and our upper bound on the estimation error under the aforementioned order of $p$, they are equal up to a factor of $\alpha_{NT}^{4}\log(NT)$.
From the error bound on $\hat{L}(u)$, we can obtain the order of the singular values of the estimation error $\hat{\Delta}_{L}(u)$ by Weyl's theorem. Let $\sigma_{1}(u)\geq\cdots\geq \sigma_{r(u)}(u)>0$ be the nonzero singular values of $L_{0}(u)$, and $\hat{\sigma}_{1}(u)\geq\cdots\geq \hat{\sigma}_{N\land T}(u)$ be the singular values of $\hat{L}(u)$. We have the following corollary.
Note that the nonzero singular values of $L_{0}(u)$ has order $\sqrt{NT}$ if $L_{0}(u)$ is formed by strong factors and factor loadings or if elements in $L_{0}(u)$ are $O(1)$. Although Theorem (ref) and Corollary (ref) are silent about whether the estimated low-rank component $\hat{L}(u)$ is low-rank or not, Corollary (ref) says that for large enough $N$ and $T$, as long as $p$ and $\alpha_{NT}$ are such that $\gamma=o(1)$, there does exist an arbitrarily large gap between the largest $r(u)$ and the remaining $((N\land T)-r(u))$ singular values of $\hat{L}(u)$. Specifically, the first $r(u)$ singular values of $\hat{L}(u)$ are of the order of $\sqrt{NT}$ while the other singular values are of the order of $\sqrt{NT}\gamma$. This confirms the intuition in Remark (ref).
This implication naturally leads to an estimator of $r(u)$. Let $\hat{r}(u)=\sum_{k}\mathbbm{1}(\hat{\sigma}_{k}(u)\geq C_{r})$ for an $(N,T)$-dependent $C_{r}$ such that $\sqrt{NT}\gamma=o(C_{r})$ and $C_{r}=o(\sqrt{NT})$. The following corollary establishes consistency of this estimator.
Theorem (ref) implies uniform consistency of $\hat{\beta}(u)$ if $\alpha_{NT}$ and $p$ are not too large such that $\gamma=o(1)$. When $\bar{r}/(N\land T)=O(p\log(pNT)/NT)$, the rate of convergence of $\hat{\beta}(u)$ is nearly optimal under $\alpha_{NT}=O(\log(NT))$. However, the rate of convergence of $\hat{\beta}(u)$ is slower than optimal when the number of covariates is small such that $p\log(pNT)/(NT)=o(\bar{r}/(N\land T))$. This is because of penalization and because the estimation error of the coefficient estimator and that of the low-rank matrix estimator are not orthogonal under the check function-involved objective function. In this case, the penalized estimator $\hat{\beta}(u)$ can be used as a first-step estimator; using it as an initialization with the rank estimator proposed in Section (ref), one can adopt the iterative estimator (ref) described in the Introduction without penalization. Although this is back to a nonconvex problem, the initial value is already lying in a small neighborhood of the true parameter by uniform consistency. Some rounds of iterations even before convergence is reached may correct the penalization bias and achieve the optimal rate (see moon2019nuclear and chernozhukov2018inference for mean regressions). The benefits of adopting such a two-step procedure instead of a fully iterative approach are three-fold. First, it provides a consistent initial guess of $\beta_{0}(u)$, potentially avoiding the issue of convergence to a local minimum. Second, as a by-product, the penalized estimation step also provides a rank estimator which is needed for the iterative approach. Third, from the Monte Carlos in Section (ref) and Appendix (ref), we can see that the penalized estimator saves a significant amount of computation time than the fully iterative one.
In this section, we illustrate the finite sample performance of our estimator using Monte Carlo simulations.
We consider the following data generating process, which is a special case of Example (ref) and is adapted from ando2019quantile:
where the $U_{it}$s are independently drawn from $\text{Unif}[0,1]$. The coefficients satisfy
The indicator functions $\mathbbm{1}_{k}(\cdot):(0,1)\mapsto\{0,1\},k=1,2,3,$ satisfy
We draw the time fixed effects $F_{1t},F_{2t}$ and $F_{3t}$ independently from $\text{Unif}[0,2]$. We generate the individual fixed effects as $\Lambda_{ki}(U_{it})=\chi_{ki}+0.1U_{it}$ where the $\chi_{ki}$s are independently drawn from $\text{Unif}[0,1]$ for $k=1,2,3$. The covariates are generated by:
where $\eta_{j,it}$ are independent $\text{Unif}[0,2]$ for all $j=1,2,3$. The parameter $\phi\in\{0.1,0.2,0.3\}$ governs the correlation between the covariates and the fixed effects. Finally, $\varepsilon_{it}$ is generated by $G^{-1}(U_{it})$ where $G$ is the cumulative distribution function of either the standard normal distribution or student's $t$-distribution with degree of freedom 2.
Since $X_{j,it}$ and $F_{kt}$ are positive for all $i,t,j$ and $k$ almost surely and $\beta_{j}(\cdot)$, $\Lambda_{ki}(\cdot)$ and $G^{-1}(\cdot)$ are strictly increasing (almost surely) for all $i,j$ and $k$, $Y_{it}$ is strictly increasing in $U_{it}$ almost surely. Let $\textbf{1}_{N\times T}$ be an $N\times T$ matrix of all ones. The $u$-th conditional quantile of $Y$ is thus $q_{Y|W}(u)=\sum_{j=1}^{3}X_{j}\beta_{j}(u)+L_{0}(u)$ where
From the model, one can see that the rank of $L_{0}(u)$ and the set of effective fixed effects vary in $u$. Also, the model allows the covariates to be correlated with $L_{0}(u)$. Higher correlation (greater $\phi$) would make it more difficult to separately estimate $\beta(u)$s and the common component $L_{0}(u)$ since it tends to yield a smaller restricted strong convexity constant $C_{RSC}$. Finally, $Y_{it}$ is allowed to have heavy tails because $\varepsilon_{it}$ can be student's $t$-distributed.
To illustrate the performance of the estimator, we conduct Monte Carlo simulations with various sample sizes at $\phi\in\{0.1,0.2,0.3\}$ and $u\in\{0.2,0.5,0.8\}$ with $G$ equal to the cumulative distribution function of either the standard normal or student's t-distribution with degree of freedom 2. As for the sample size, $(N,T)\in \{(200,200),(300,300),(400,400),(500,500),\newline(200,300),(300,200), (200,400),(400,200),(200,500),(500,200)\}$. The first four sample sizes with $N=T$ allow us to see the convergence of the estimator. The other six sample sizes with $N\neq T$ and $N\land T=200$ allow us to test the theory which suggests that the rate of convergence does not depend on $N\lor T$ under a fixed $p$.
We use multiple measures to evaluate the estimator's performance. For the three coefficients $\beta_{j}(u)$s, we compute their average squared bias $\text{Bias}_{\beta}^{2}$ and variance $\text{Var}_{\beta}$ over 100 simulation replications\footnote{Specifically, $\text{Bias}_{\beta}^{2}\coloneqq\sum_{j=1}^{3}\left(\sum_{b=1}^{100}(\hat{\beta}_{j,b}(u)-\beta_{j}(u))/100\right)^{2}/3$ where $\hat{\beta}_{j,b}(u)$ is the estimator of $\beta_{j}(u)$ in the $b$-th simulation. The variance $\text{Var}_{\beta}\coloneqq\sum_{j=1}^{3}\left[\sum_{b=1}^{100}\hat{\beta}_{j,b}(u)^{2}/100-\left(\sum_{b=1}^{100}\hat{\beta}_{j,b}(u)/100\right)^{2}\right]/3$.}. For the low-rank component, we compute the average squared Frobenius norm of the estimation error:
where $\hat{L}_{b}(u)$ is the estimated $L_{0}(u)$ in the $b$-th simulation ($b=1,...,100$). Besides these measures, we also look at the mean squared error of the estimated conditional quantile function $\hat{q}_{Y|W}(u)$, defined as follows ando2019quantile:
Note that $\text{MSE}_{q}$ can be small even if $\text{MSE}_{L}$ is large; the latter is also affected by $C_{RSC}$, which in turn, is affected by the correlation between the covariates and $L_{0}(u)$ (determined by parameter $\phi$). Hence, by comparing $\text{MSE}_{q}$ with $\text{MSE}_{L}$ across different values of $\phi$, we can see how restricted strong convexity affects the results.
Finally, we also record the average computation time of the estimator.
To make comparison between the nuclear norm penalized estimator proposed in this paper and alternative estimators, we compute the following iterative estimator adapted from ando2019quantile as introduced in the Introduction, and the pooled estimator:
The finite sample bias of these two estimators provide us with two benchmarks for bias analysis; the bias of any reasonably estimator should be close to the iterative estimator and at least smaller than the pooled estimator. As discussed in the Introduction, $(\hat{\beta}^{It}(u),\hat{\Lambda}(u),\hat{F})$ is obtained by a nonconvex problem which needs the number of fixed effects $r(u)$ to be known. Nonconvexity may result in convergence to local minima, but when the global minimum is indeed achieved and $r(u)$ is correctly set, its finite sample bias is presumably small as there is no penalization. On the other hand, the pooled estimator $\hat{\beta}^{Po}$ completely ignores the fixed effects. Since the covariates are correlated with the fixed effects by construction, $\hat{\beta}^{Po}$ is inconsistent. By comparing the bias of our penalized estimator with these two estimators, we can see what regularization buys and costs.
For our penalized estimator, we present the details about our algorithm in Appendix (ref). The algorithm is adapt from the Augmented Lagrangian Multiplier method proposed in lin2010augmented, candes2011robust and yuan2013sparse. The method was originally designed for $u=0.5$ with no covariates. We extend it to accommodate any $u\in(0,1)$ with covariates. From the simulation results, the new algorithm works fast and well. For $\lambda$, here we simply set it to be $\lambda=\log(NT)\sqrt{N\lor T}/(3.6NT)$. This choice of $\lambda$ slightly slows down the theoretical rate of convergence because of the additional $\log(NT)$, but it works no matter what $C_{Cone}$ is. Alternatively, one may use cross-validation or adopt the BIC criterion proposed in belloni2019high to select $\lambda$.
For the iterative estimator, we compute $(\hat{\beta}^{It}(u),\hat{\Lambda}(u),\hat{F})$ by iteratively running quantile regression to update $\beta(u),\Lambda(u)$ and $F$ using the true number of the fixed effects $r(u)$. This algorithm is adapted from ando2019quantile which was for $i$-specific $\beta(u)$. Specicically, for each $i$, obtain $\Lambda_{i}(u)$ by quantile regression using the time series variation by fixing $\beta$ and $F$. For each $t$, update $F_{t}$ by quantile regression using the cross-sectional variation by fixing $\Lambda(u)$ and $\beta$. Then fixing $\Lambda(u)$ and $F$, update $\beta$ by pooled quantile regression. Iterate until converged. For initialization, for each $u$, use $\hat{\beta}^{Po}(u)$ as the initial value of $\beta(u)$, obtain the residual matrix $R$, and initialze the $F$ matrix by the eigenvector matrix of $R'R$ multiplied by $\sqrt{T}$. The termination criterion for the iteration, adopted from ando2019quantile, is set to be the same as that for the nuclear norm penalized estimator (see Appendix (ref) for details) to make the results, especially the computation time, comparable.
Tables (ref) and (ref) present the results for $\phi=0.2$ with standard normal and student $t$-distributed $\varepsilon_{it}$, respectively. The results for $\phi=0.1$ and $0.3$ are provided in Appendix (ref). The Columns Nu, It and Po report results of our penalized estimator, the iterative estimator (ref) and the pooled estimator (ref). The column Time reports on average how many seconds each estimator takes. All experiments were performed in MATLAB 2019b under Windows 10 on a desktop computer with 10-core 2.8GHz Intel i9 processor and 16GB RAM using parallel computing. From the results, we have the following key observations.
First, for the penalized estimator (Columns Nu), we can see that in all specifications, the bias and variance of $\hat{\beta}(u)$ shrink as $N$ and $T$ both increase for all $u$ considered (see the $\text{Bias}^{2}_{\beta}\times 100$ and $\text{Var}_{\beta}\times 10^{4}$ columns). The same holds for the estimation errors of the estimated low-rank common component $\hat{L}(u)$ and the conditional quantile function $\hat{q}_{Y|W}(u)$. In particular, the results are robust to heavy-tailed data (Table (ref)). On the other hand, when $N\land T$ is fixed, increasing $N\lor T$ does not lead to big improvement in the penalized estimator's performance. This echoes Theorem (ref) which says the rate of convergence is dominated by $1/(N\land T)$ for a fixed $p$. Meanwhile, we can see that the average squared estimation error of the conditional quantile function, $\text{MSE}_{q}$, is smaller than that of the low-rank common component, $\text{MSE}_{L}$. This is due to the separation issue and restricted strong convexity, arisen from the correlation between the covariate matrices and the low-rank component. Together with the results in Appendix (ref), we can see that $\text{MSE}_{L}$ becomes smaller and closer to $\text{MSE}_{q}$ when $\phi$, and thus the correlation between the covariates and the low-rank component, gets smaller. This is because a smaller $\phi$ tends to yield a larger restricted strong convexity constant $C_{RSC}$.
Next, let us compare the performance of the different estimators we consider. Comparing the pooled estimator (Column Po) and the penalized estimator, we can see that even though the latter is biased due to regularization, the bias is still much smaller than that of the pooled estimator, whose bias does not shrink at all as the sample size increases due to endogeneity. Now let us compare our penalized estimator with the iterative estimator (Columns It), which seems not to suffer from the local minima issue in these experiments. The penalized estimator has similar performance to the iterative one with relatively larger bias and smaller variance
in many cases (e.g. $u=0.5$ in all cases). Admittedly, the penalized estimates of the coefficients tend to have a larger mean squared error since the theoretical convergence rate is slow. However, this is not always the case especially when the sample size is relatively small. For instance, when $u=0.8$, $(N,T)=(200,200)$ for standard normal error, and $(N,T)=(200,200),(300,300)$, etc. for student $t$-distributed error, the penalized estimator
has smaller bias and smaller mean squared error in the coefficients than the iterative estimator.
In terms of computation speed, the penalized estimator in these experiments is much faster than the iterative one. The former rarely takes more than half a minute and in most cases only a few seconds. But the iterative estimator takes about 20-40 times longer, and when the sample size is large ($N=T=500$), it can take almost 15 minutes. Note that when we compute the iterative estimator, we treat the number of fixed effects as known. In practice, this additional parameter also needs to be estimated, so the actual computation time can be even longer.
In this paper, we study a conditional quantile panel data model with interactive fixed effects. By exploiting the low-rankness of the matrix formed by the interactive fixed effects, we propose a nuclear norm penalized estimator. The estimator jointly estimates the coefficients and the low-rank matrix by solving a convex problem. We derive a uniform error bound on the estimator and in turn, establish uniform consistency. Based on the error bound, we also construct a consistent estimator of the number of fixed effects at any given quantile level. From the Monte Carlo simulations, our estimator performs well and computation is much more efficient than a related iterative estimator.
We conjecture that after the penalized estimator and the rank estimator are obtained, by using them as the initial value, a few rounds of iterations based on the iterative estimator's minimization problem would remove the bias and restore the rate of convergence of the coefficient estimator that is slowed down by penalization. Inference may also be available based on such post-penalization procedures. We leave these for future work.