High Dimensional Latent Panel Quantile Regression with an Application to Asset Pricing

Supplementary Material for “High Dimensional Latent Panel Quantile Regression with an Application to Asset Pricing”

Introduction

A central question in asset pricing is to explain why certain assets pay higher returns than others. The Arbitrage Pricing Theory (ross1976arbitrage) and the Fama-French three factors model (fama1993common) explains the asset return variations by a linear combination of common risk factors. Assets with similar exposure to a common factor shall rise and fall together (cochrane2009asset). However, empirical evidence appears to indicate that the firm characteristics, rather than common factors, can also explain the variations in stock returns (daniel1997evidence), which suggests a characteristic-based model.

We generalize both modeling approaches and propose a “Characteristics + Latent Factors” quantile asset pricing framework. By incorporating the “Characteristics", we improve the economic interpretability and explanatory power of the model. On the other hand, the finance literature has documented a zoo of new characteristics, and the proliferation of characteristics in this “variable zoo” leads to a concern about which characteristics really provide independent information about returns (cochrane2011presidential). Our model addresses this issue by imposing a sparse structure, meaning although a large set of characteristics is available, only a much smaller subset of them might have predictive power. We also incorporate “Latent Factors” to capture the common variations in asset returns. One additional benefit of having this part is that it might help alleviate the “omitted variable bias” problem (giglio2018asset). As in the literature, typically, these latent factors are estimated via principal component analysis, which means all possible latent explanatory variables might be important for prediction although their individual contribution might be small, we term this as the dense part. \footnote{More about sparse modeling and dense modeling can be found in giannone2017economic. See also chernozhukov2017lava.} Hence, our framework allows for “Sparse + Dense" modeling with large scale panel data that consist of a large number of asset returns that are allowed to be weakly correlated across time. In addition, we focus on understanding the quantiles (hence the entire distribution) of returns rather than just the mean, in line with the recent interest in quantile factor models (e.g. ando2020quantile, chen2018quantile, ma2019estimation, feng2019nuclear, and sagner2019three). Our quantile asset pricing framework also inherits micro-foundation from the seminal quantile preference framework(manski1988ordinal,rostek2010quantile,giovannetti2013asset) and, in particular, the dynamic quantile preference framework of de2019dynamic.

Specifically, with $Y_{i,t}$ as the excess return of asset $i$ in period $t$, $X_{i,t}$ as a $p$-dimensional vector of observable characteristics such as return volatility and trading volume, we study the following high dimensional latent panel quantile regression model:

where $\theta(\tau) \in \mathbb{R}^p$ is the vector of coefficients, $g_t(\tau)$ is an $r_{\tau}$-dimensional vector of unobservable factor returns, $\lambda_i(\tau)$ represents the factor loadings which captures the sensitivity of asset $i$ on the $r_{\tau}$ factors, $\tau \in [0,1]$ is the quantile index. We allow for the possibility of quantile dependence of sensitivity to risk factors as such evidence has been reported in the literature (e.g., ando2020quantile). For notation simplicity, we often denote $\Pi_{i,t} (\tau) = \lambda_i(\tau)^{\prime} g_t(\tau)$, then $\Pi(\tau)$ is a low-rank matrix with unknown rank $r_\tau$. Thus, with $ F_{ Y_{i,t}| X_{i,t}; \theta(\tau), \lambda_i(\tau), g_{t}(\tau)}$ is the cumulative distribution function of $Y_{i,t}$ conditioning on $X_{i,t}$, $\theta(\tau)$ and $\lambda_i(\tau), g_t(\tau)$, we model the quantiles of returns (instead of expected returns) as a linear combination of the characteristics and latent factors. Our framework allows for the possibility of lagged dependent data. Here, we allow for the number of characteristics $p$, and the time horizon $T$, to grow to infinity as $n$ grows. Throughout, we focus on the case where $p$ is large, possibly much larger than $n T$, but for the true model $\theta(\tau)$ is sparse and has only $ s_{\tau} \ll p$ non-zero components.

Our framework is flexible enough that allows us to jointly answer the following three questions in asset pricing: (i) Which characteristics are important to explain the time series and cross-section of stock returns, after controlling for the factors? (ii) How much would the latent factors explain stock returns after controlling for firm characteristics? (iii) Does the relationship of stock returns and firm characteristics change across quantiles? The first question is related to the recent literature on variable selection in asset pricing using machine learning ( kozak2019shrinking, feng2019taming, han2018firm). The second question is related to an classical literature starting from 1980s on statistical factor models of stock returns (chamberlain1983arbitrage, connor1988risk and recently lettau2018estimating). The third question extends the literature in late 1990s on stock return and firm characteristics (daniel1997evidence, daniel1998characteristics) and further asks whether the relationship is heterogenous across quantiles.

There are several key features of considering prediction problem at the panel quantile model in this setting. First, stock returns are known to be asymmetric and exhibit heavy tail, thus modeling different quantiles of return provides extra information in addition to models of first and second moments. Second, quantile regression provides a richer characterization of the data, allowing heterogeneous relationship between stock returns and firm characteristics across the entire return distribution. Third, the latent factors might also be different at different quantiles of stock returns. Finally, quantile regression is more robust to the presence of outliers relative to other widely used mean-based approaches. Using a robust method is crucial when estimating low-rank structures (see e.g. she2017robust). As our framework is based on modeling the quantiles of the response variable, we do not put assumptions directly on the moments of the dependent variable.

Our main goal is to consistently estimate both the sparse part and the low-rank matrix. Recovery of a low-rank matrix, when there are additional high dimensional covariates, in a nonlinear model can be very challenging. The rank constraint will result in the optimization problem NP-hard. In addition, estimation in high dimensional regression is known to be a challenging task, which in our frameworks becomes even more difficult due to the additional latent structure. We address the former challenge via nuclear norm regularization which is similar to candes2009exact in the matrix completion setting. Without covariates, the estimation can be done via solving a convex problem, and similarly there are strong statistical guarantees of recovery of the underlying low-rank structure. We address the latter challenge by imposing $\ell_1$ regularization on the vector of coefficients of the control variables, similarly to belloni2011l1 which mainly focused on the cross-sectional data setting. Note that with regards to sparsity, we must be cautious, specially when considering predictive models (she2017robust). Furthermore, we explore the performance of our procedure under settings where the vector of coefficients can be dense (due to the low-rank matrix).

We view our work as complementary to the low dimensional quantile regression with interactive fixed effects framework as of the very recent work of feng2019nuclear, and the mean estimation setting in moon2018nuclear. However, unlike moon2018nuclear and feng2019nuclear, we allow the number of covariates $p$ to be large, perhaps $p \gg nT$. This comes with different significant challenges. On the computational side, it requires us to develop novel estimation algorithms, which turns out can also be used for the contexts in moon2018nuclear and feng2019nuclear. On the theoretical side, allowing $p \gg nT$ requires a systematically different analysis as compared to feng2019nuclear, as it is known that ordinary quantile regression is inconsistent in high dimensional settings ($p \gg nT$), see belloni2011l1.

\paragraph*{Related Literature} Our work contributes to the recent growing literature on panel quantile model. abrevaya2008effects, graham2018quantile, arellano2017quantile, considered the fixed $T$ asymptotic case. kato2012asymptotics formally derived the asymptotic properties of the fixed effect quantile regression estimator under large $T$ asymptotics, and galvao2016smoothed further proposed fixed effects smoothed quantile regression estimator. galvao2011quantile works on dynamic panel. koenker2004quantile proposed a penalized estimation method where the individual effects are treated as pure location shift parameters common to all quantiles, for other related literature see lamarche2010robust, galvao2010penalized. We refer to Chapter 19 of koenker2017handbook for a review. Furthermore, our framework can be viewed as a generalization of the model in ando2020quantile which considered panel quantile model with independent errors and low dimensional covariates.

Our work also contributes to the literature on nuclear norm penalisation, which has been widely studied in the machine learning and statistical learning literature, fazel2002matrix, recht2010guaranteed, koltchinskii2011nuclear, rohde2011estimation, negahban2011estimation, brahma2017reinforced. Recently, in the econometrics literature athey2018matrix proposes a framework of matrix completion for estimating causal effects, bai2017principal for estimating approximate factor model, chernozhukov2018inference considered the heterogeneous coefficients version of the linear panel data interactive fixed model where the main coefficients has a latent low-rank structure, bai2019robust for robust principal component analysis, and bai2019matrix for imputing counterfactual outcome.

Finally, our results contribute to a growing literature on high dimensional quantile regression. wang2012quantile considered quantile regression with concave penalties for ultra-high dimensional data; zheng2015globally proposed an adaptively weighted $\ell_1$-penalty for globally concerned quantile regression. Screening procedures based on moment conditions motivated by the quantile models have been proposed and analyzed in he2013quantile and wu2015conditional in the high-dimensional regression setting. We refer to koenker2017handbook for a review.

To sum-up, our paper makes the following contributions. First, we propose a new class of models that consist of both high dimensional regressors and latent factor structures. We provide a scalable estimation procedure, and show that the resulting estimator is consistent under suitable regularity conditions. Second, the high dimensional and non-smooth objective function require innovative strategies to derive all the above-mentioned results. In particular, our paper allows for serial dependence and this flexibility is important for the panel data case. It is well known that dealing with data dependence is a non-trivial problem, and there are additional challenges for high dimensional models even for those without incorporating the latent factors: the extension of lasso-based methods with relaxing the i.i.d. assumption and other restrictive assumptions (for instance, Gaussianity), is just beginning to occur, e.g. for time series data some recent development along this line can be found in wong2020lasso. Those lead to the novel use in our proofs of some techniques from the high dimensional statistics and econometrics literature, such as the localization argument from belloni2011l1, and the new loss function introduced in padilla2020adaptive; from spectral theory, namely, properties of nuclear norm studied in elsener2018robust, and concentrations results by chatterjee2015matrix; and from empirical process theory yu1994rates,wellner2013weak. We also generalize the sampling and smoothness assumption of belloni2011l1 by considering panel data with weak correlation across time. In particular, we refer readers to yu1994rates for thorough discussions on $\beta$-mixing.

On the theoretical side, we also present multiple results that entirely differ from those in belloni2011l1. In addition to allowing time dependence and a latent factor structure, we can consistently estimate the conditional quantiles without requiring a minimum eigenvalue condition on the behavior of the design matrix, which is typically required in the literature (e.g. belloni2011l1). Relative to approaches that incorporate latent factor structure but relies on the squared loss, the proposed estimators inherit from quantile regression certain robustness properties to the presence of outliers and heavy-tailed distributions. Finally, we apply our proposed model and estimator to a large-dimensional panel of financial data in the US stock market and find that different return quantiles have different selected firm characteristics and that the number of latent factors can be also be different.

\paragraph*{Outline} The rest of the paper is organized as follows. Section (ref) introduces the high dimensional latent quantile regression model, and provides an overview of the main theoretical results. Section (ref) presents the estimator and our proposed ADMM algorithm. Section (ref) discusses the statistical properties of the proposed estimator. Section (ref) provides simulation results. Section (ref) consists of the empirical results of our model applied to a real data set. The proofs of the main results are in the Supplementary Material.

\paragraph*{Notation} For $m \in \mathbb{N}$, we write $[m] \,=\, \{1,\ldots,m\}$. For a vector $ v \in \mathbb{R}^p$ we define its $\ell_0$ norm as $\| v\|_0 = \sum_{j =1}^{p} 1\{ v_j \neq 0 \} $, where $1\{ \cdot \}$ takes value $1$ if the statement inside $\{\}$ is true, and zero otherwise; its $\ell_1$ norm as $\Vert v \Vert_1 = \sum^p_{j=1}\vert v_j\vert$. We denote $\Vert v \Vert_{1,n,T} = \sum_j^{p}\hat{\sigma}_j\vert v_j\vert$ the $\ell_1$-norm weighted by $\hat{\sigma}_j$'s (defined in eq((ref))). The Euclidean norm is denoted by $\|\cdot \|$, thus $\|v\| = \sqrt{ \sum_{j=1}^p v_j^2 }$. If $A \in \mathbb{R}^{n \times T}$ is a matrix, its Frobenius norm is denoted by $\|A\|_F \,=\,\sqrt{ \sum_{ i= 1}^n \sum_{t=1}^T A_{i,t}^2 }$, its spectral norm by $\| A\|_2 \,=\, \sup_{x \,:\, \|x\| = 1 }\sqrt{ x^{\prime} A^{\prime} A x} $, its infinity norm by $\|A\|_{\infty} \,= \, \max\{ \vert A_{i,j }\vert \,:\, i \in [n],\,\, j \in [T] \} $ , its rank by $\text{rank}(A)$, and its nuclear norm by $\|A\|_* = \text{trace}(\sqrt{A^{\prime} A})$ where $A^{\prime}$ is the transpose of $A$. The $j$th column $A$ is denoted by $A_{\cdot,j}$. Furthermore, the multiplication of a tensor $X \in \mathbb{R}^{I_1 \times \ldots \times I_m }$ with a vector $\theta \in \mathbb{R}^{I_m}$ is denoted by $Z := X \theta \in \mathbb{R}^{ I_1 \times \ldots \times I_{m-1} }$, and, explicitly, $Z_{i_1,\ldots i_{m-1} } = \sum_{j = 1}^{I_m} X_{i_1,\ldots, i_{m-1},j } \, \theta_{j} $. We also use the notation $a \vee b= \max\{a,b\}$, $a \land b = \min\{a, b\}$, $(a)_{-} = \max\{-a,0\}$. For a sequence of random variables $\{z_j\}_{j=1}^{\infty}$ we denote by $\sigma(z_1,z_2,\ldots)$ the sigma algebra generated by $\{z_j\}_{j=1}^{\infty}$. Finally, for sequences $\{a_n\}_{n=1}^{\infty}$ and $\{b_n\}_{n=1}^{\infty}$ we write $a_n \asymp b_n$ if there exists positive constants $c_1$ and $c_2$ such that $c_1 b_n \leq a_n \leq c_2 b_n$ for sufficiently large $n$.

The Estimator and Overview of Rate Results

Basic Setting

The setting of interest corresponds to a high dimension latent panel quantile regression model, where $Y \in \mathbb{R}^{n\times T} $, and $X \in \mathbb{R}^{n \times T \times p}$ satisfying

where $i$ denotes subjects, $t$ denotes time, $\theta(\tau) \in \mathbb{R}^p$ is the vector of coefficients, $\Pi(\tau) \in \mathbb{R}^{n \times T}$ is a low-rank matrix with unknown rank $r_\tau \ll \min \{n, T\}$, $\tau \in [0,1]$ is the quantile index, and $ F_{ Y_{i,t}| X_{i,t}; \theta(\tau), \Pi_{i,t}(\tau) }$ is the cumulative distribution function of $Y_{i,t}$ conditioning on $X_{i,t}$, $\theta(\tau)$ and $\Pi_{i,t}(\tau)$. Thus, we model the quantile function at level $\tau$ as a linear combination of the predictors plus a low-rank matrix. Here, we allow for the number of covariates $p$, and the time horizon $T$, to grow to infinity as $n$ grows. Throughout the paper the quantile index $\tau \in (0,1)$ is fixed. We mainly focus on the case where $p$ is large, possibly much larger than $nT $, but for the true model $\theta(\tau)$ is sparse and has only $ s_{\tau} \ll p$ non-zero components. Mathematically, $s_\tau := \Vert \theta(\tau) \Vert_0$.

When $\Pi_{i,t}(\tau) \,=\, \lambda_i(\tau)^{\prime} g_t(\tau)$, with $\lambda_i(\tau), \, g_t(\tau) \in \mathbb{R}^{r_{\tau}}$, this immediately leads to the following setting

where we model the quantile function at level $\tau$ as a linear combination of the covariates (as predictors) plus a latent factor structure. This is directly related to the panel data models with interactive fixed effects literature in econometrics, e.g. linear panel data model (bai2009panel), nonlinear panel data models (chen2014estimation, chen2014nonlinear).

Note, for eq ((ref)), additional identification restrictions are needed for estimating $\lambda_i(\tau)$ and $g_t(\tau)$ (see bai2013principal). In addition, in nonlinear panel data models, this creates additional difficult in estimation, as the latent factors and their loadings part induce a nonconvex quantile regression problem. However, we deal with this in the following subsection via using a nuclear norm constraint. \footnote{ Different identification conditions might result in different estimation procedures for $\lambda$ and $f$, see bai2012statistical and chen2014estimation.}

Estimator

In this subsection, we describe the high dimensional latent quantile estimator. With the sparsity and low-rank constraints in mind, a natural formulation for the estimation of $(\theta(\tau),\Pi(\tau))$ is

where $\rho_{\tau}(t) \,=\, (\tau - 1\{ t\leq 0 \})t $ is the quantile loss function as in Koenker2005, $s_{\tau}$ is a parameter that directly controls the sparsity of $\tilde{\theta}$, and $r_{\tau}$ controls the rank of the estimated latent matrix.

While the formulation in ((ref)) seems appealing, as it enforces variable selection and low-rank matrix estimation simultaneously, ((ref)) is a non-convex problem due to the constraints posed by the $\|\cdot\|_0$ and $\mathrm{rank}(\cdot)$ functions. We propose a convex relaxation of ((ref)). Inspired by the seminal works of tibshirani1996regression and candes2009exact, we formulate the problem as the following

where $ \nu_{1} >0$ and $\nu_2 >0$ are tuning parameters, and $w_1,\ldots,w_p$ are user specified weights (more on this in Section (ref) ). Notice that $\|\cdot\|_*$ is the nuclear norm defined on Page 4. The nuclear norm regularization works on the singular value of a matrix, the intuition is that via penalization with the nuclear norm, the resulting problem will be convex. Just as $\ell_1$-minimization is the tightest convex relaxation of the combinatorial $\ell_0$-minimization problem, nuclear-norm minimization is the tightest convex relaxation of the NP-hard rank minimization problem, see candes2010matrix.

In principle, one can use any convex solver software to solve ((ref)), since this is a convex optimization problem. However, for large scale problems a more careful implementation might be needed. Section (ref) presents a scheme for solving ((ref)) that is based on the ADMM algorithm (boyd2011distributed).

Summary of results

We now summarize our main results. For the model defined in ((ref)):

Although our theoretical analysis builds on the work by belloni2011l1, there are multiple challenges that we must face in order to prove the consistency of our estimator. First, the construction of the restricted set now involves the nuclear norm penalty. This requires us to define a new restricted set that captures the contributions of the low-rank matrix. Second, when bounding the empirical processes that naturally arise in our proof, we have to simultaneously deal with the sparse and dense components. Furthermore, throughout our proofs, we have to carefully handle the weak dependence assumption that can be found in Section (ref).

High Dimensional Latent Panel Quantile Regression

In this subsection, we describe the main steps of our proposed ADMM algorithm, details can be found in Section (ref). We start by introducing slack variables to the original problem ((ref)). As a result, a problem equivalent to ((ref)) is

To solve ((ref)), we propose a scaled version of the ADMM algorithm which relies on the following Augmented Lagrangian

where $\eta >0$ is a penalty parameter.

Notice that in ((ref)), we have followed the usual construction of ADMM via introducing the scaled dual variables corresponding to the constraints in ((ref)) -- those are $U_V$, $U_W$, $U_{\Pi}$, and $U_{\theta}$. Next, recall that ADMM proceeds by iteratively minimizing the Augmented Lagrangian in blocks with respected to the original variables, in our case $(V, \tilde{ \theta},\tilde{ \Pi})$ and $(W,Z_{\theta},Z_{\Pi})$, and then updating the scaled dual variables (see Equations 3.5--3.7 in boyd2011distributed). The explicit updates can be found in the Supplementary Material. Here, we highlight the updates for $Z_{\theta}$, $\tilde{\Pi}$, and $V$. For updating $Z_{\theta}$ at iteration $k+1$, we solve the problem \[ \displaystyle Z_{ \theta }^{(k+1)} \,\leftarrow \,\underset{Z_{\theta} \in \mathbb{R}^{ p} }{\arg \min}\left\{ \frac{1}{2}\| Z_{\theta} -\tilde{\theta}^{(k+1)} + U_{\theta}^{(k)} \|_F^2 + \frac{ \nu_{1}}{\eta} \sum_{j =1 }^{p} w_j \vert (Z_{\theta})_j \vert \right\}. \] This can be solved in closed form exploiting the well known thresholding operator, see the details in Section (ref). As for updating $\tilde{\Pi}$, we solve

via the singular value shrinkage operator, see Theorem 2.1 in cai2010singular.

Furthermore, we update $V$, at iteration $k+1$, via

which can be found in closed formula by Lemma 5.1 from ali2016multiple.

Finally, it is immediate to modify the proposed ADMM to the case when there are no covariates ($\beta(\tau) =0$), see Section (ref) in the Supplementary Material. Hence, our proposed estimation procedure can be applied to settings (i) with low dimensional covariates, or (ii) without covariates. However, in what follows, we focus on the high dimensional covariates setting.

Theory

The purpose of this section is to provide statistical guarantees for the estimator developed in the previous section. We focus on estimating the quantile function, allowing for the high dimensional scenario where $p$ and $T$ can grow as $n$ grows.

Estimating the quantiles

We show that our proposed estimator is consistent for estimating the conditional quantiles.

Throughout, we treat $\Pi(\tau)$ as fixed parameters. As for the data generation process, our next condition requires that the observations are independent across $i$, and weakly dependent across time.

Note that Assumption (ref) is related to the sampling and smoothness assumption of belloni2011l1. Furthermore, we highlight that similar to belloni2011l1, our framework is rich enough that avoids imposing Gaussian modeling constraints. However, unlike belloni2011l1, we consider panel data with weak correlation across time. In particular, we refer readers to yu1994rates for thorough discussions on $\beta$-mixing. Another difference with belloni2011l1 is that we do not require any smoothness assumption on the condition density function of the response given the covariates.

It is worth mentioning that the parameter $\mu$ in Assumption (ref) controls the strength of the time dependence in the data while parameter $\mu’$ will relate to the tuning parameters and will impact the rates. As we decrease the value of $\mu’$ condition ((ref)) is easier to satisfy but the tuning parameters increase and slow down the rates of convergence. Furthermore, in the case that $\{(Y_{i,t},X_{i,t})\}_{ i \in [n], t \in [T] }$ are independent our theoretical results will hold without imposing ((ref)).

Next, we require that along each dimension the second moment of the covariates is one. We also assume that the second moments can be reasonably well estimated by their empirical counterparts.

Assumption (ref) appeared as Condition D.3 in belloni2011l1. It is met by general models on the covariates, see for instance Design 2 in belloni2011l1.

Using the empirical second order moments $\{\hat{\sigma}_j^2 \}_{j=1}^p$, we analyze the performance of the constrained estimator

where $\nu_1,\nu_2 \,>\, 0$ are tuning parameters and $\|\tilde{\theta}\|_{1,n,T} \,:=\, \sum_{j =1 }^{p} \hat{\sigma}_j \vert \tilde{\theta}_j\vert$, \[ \displaystyle \hat{Q}_{\tau}(\tilde{\theta},\tilde{\Pi}) \,=\, \frac{1}{nT}\sum_{i=1}^n \sum_{t=1}^T \rho_{\tau}( Y_{i,t} - X_{i,t}^{\prime}\tilde{\theta} - \tilde{\Pi}_{i,t}), \] with $\rho_{\tau}$ as defined in Section (ref).

Our main result of this subsection is provided next.

Theorem (ref) shows that we can consistently estimate the conditional quantiles at a rate that combines the sparse and dense signal magnitudes. This result differs from belloni2011l1 in several ways. First, belloni2011l1 work with cross sectional data and focus on the estimation of the vector of parameter coefficients instead of the conditional quantiles. Therefore the assumptions for Theorem (ref) and those in belloni2011l1 are very different. For example, belloni2011l1 required a minimum eigenvalue condition on the behavior of the design matrix, whereas we can avoid this. Furthermore, unlike belloni2011l1, Theorem (ref) allows for panel data with dependence across time. This translates into a technical challenge to arrive at Theorem (ref), since it is not possible to use the analysis from belloni2011l1 which relies heavily on independence. In fact, the minimum eigenvalue condition in belloni2011l1 (Condition D.4) can be difficult to be verified in practice without the independence data assumption, and hence Theorem (ref) has a different setting than those results in belloni2011l1. Also, our setting and estimator involve the latent factors which make both the estimator and theory more complex than the framework in belloni2011l1. Lastly, while some of the ideas in the proof of Theorem (ref) come from padilla2020adaptive, all of the technical steps involved in the proof of Theorem (ref) are novel.

On another note, the rate of $T$ in Theorem (ref) must be such that $T \rightarrow \infty$ and $T$ satisfies ((ref)). The parameter $\mu'$ in ((ref)) controls the terms $c_T=\lceil T^{1/(1+\mu')}\rceil$ and $d_T= \lfloor T/(2c_T)\rfloor$ both of which determine the tuning parameters and imply specific rates of convergence. Cearly, $c_T \cdot d_T \asymp T$. Moreover, the parameter $c_T$ measures the strength of the dependence in the data, whereas $d_T$ can be interpreted as the effective number of independent samples across time. In particular, when $c_T \asymp 1$, we have that $\{\epsilon_{i,t}\}$ are basically independent across time, and so the final rate depends on $T$. As a different example, since we would like $d_T \gg c_T$, with $\mu'=2$ we have $c_T \sim T^{1/3}$ and $d_T \sim T^{2/3}$ so that $\nu_1 \sim (T^{-1/6}+ \sqrt{T^{1/3}/n})\sqrt{\log(npT)}$ and $\nu_2 \sim \frac{1}{\sqrt{n}T^{2/3}}+\frac{1}{n T^{1/3}}$. In this case condition ((ref)) requires the parameter $\mu$ that governs the mixing speed to be such that $npT \ll T^{\mu/3}$.

Estimating the coefficients and latent factors

We show that our proposed estimator is consistent for estimating the vector of coefficients and latent factors separately in a broad range of models, and in some cases attains minimax rates, as in candes2009tight.

In order to obtain our main result of this subsection, we first provide some additional notation and assumptions. For a fixed $\tau >0$, we assume that ((ref)) holds. We also let $T_{\tau}$ be the support of $\theta(\tau)$, thus \[ T_{\tau } \,=\, \left\{j \in [p] \,:\, \theta_j(\tau) \neq 0 \right\}, \] and we write $s_{\tau} \,=\, \vert T_{\tau} \vert $, and $r_{\tau} \,=\, \text{rank}(\Pi(\tau))$.

As it can been seen in Lemma (ref) from the Supplementary Material, the error of our estimator defined in ((ref)), $(\hat{\theta}(\tau) - \theta(\tau) , \hat{\Pi}(\tau) - \Pi (\tau) ) $, belongs to a restricted set, which in our framework is defined as

for an appropriate positive constant $C_0$.

Similar in spirit to other high dimensional settings such as those in candes2007dantzig, bickel2009simultaneous, belloni2011l1 and dalalyan2017prediction, we impose an identifiability condition involving the restricted set which is expressed next and will be used in order to attain our main results. This is the key difference with the analysis in Section (ref). Before arriving at our next condition, we introduce some notation.

For $m\geq 0$, we denote by $\overline{T}_{\tau }(\delta,m) \subset \{ 1,\ldots,p\} \backslash T_{\tau }$ the support ot the $m$ largest components, excluding entries in $T_{\tau }$, of the vector $(\vert\delta_1\vert, \ldots, \vert \delta_p \vert)^T$. We also use the convention $\overline{T}_{\tau }(\delta,0) \,=\,\emptyset$.

Few comments are in order. First, if $\Delta =0$ then ((ref)) and ((ref)) become the restricted identifiability and nonlinearity conditions as of belloni2011l1. Second, the denominator of ((ref)) contains the term $\|\Delta\|_F /(\sqrt{nT})$. To see why this is reasonable, consider the case where $\mathbb{E}(X_{i,t}) =0$, and $X_{i,t}$ are i.i.d.. Then \[ J_{\tau }(\delta,\Delta) \,=\, \underline{f} \, \mathbb{E}( (\delta^{\prime} X_{i,t})^2 ) + \frac{ \underline{f} }{nT} \|\Delta\|_F^2. \] Hence, $\|\Delta\|_F /(\sqrt{nT})$ appears also in the numerator of ((ref)) and it is not restrictive its presence in the denominator of ((ref)).

We now state our result for estimating $\theta(\tau)$ and $\Pi(\tau)$.

Theorem (ref) gives an upper bound on the performance of $(\hat{\theta}(\tau),\hat{\Pi}(\tau))$ for estimating the vector of coefficients $\theta(\tau)$ and the latent matrix $\Pi(\tau)$. For simplicity, consider the case of i.i.d data. Then the convergence rate of our estimation of $\theta(\tau)$, under the Euclidean norm, is in the order of $(\sqrt{s_{\tau}} + \sqrt{ r_{\tau} })/\min\{ \sqrt{n}, \sqrt{d_T} \} $, if we ignore all the other factors. Hence, we can consistently estimate $\theta(\tau)$ provided that $ \max\{s_{\tau},r_{\tau} \} << \min\{ n,T \} $. This is similar to the low-rank condition in negahban2011estimation. In the low dimensional case $s_{\tau} = O(1)$, the rate $\sqrt{ r_{\tau} }/\min\{ \sqrt{n}, \sqrt{d_T} \} $ matches that of Theorem 1 in moon2018nuclear. However, unlike moon2018nuclear, our estimator is based on a loss function that is robust to outliers, and our assumptions also allow for weak dependence across time, making our framework potentially more general. Furthermore, the same applies to our rate on the mean squared error for estimating $\Pi(\tau)$, which also matches that in Theorem 1 of moon2018nuclear.

With regards to the novelty of Theorem (ref), we highlight that, while its proof is similar in spirit to that of Theorem 2 in belloni2011l1, there are some significant differences. First, the construction of the restricted set used in Theorem (ref) involves two different penalties which makes challenging to disentangle the behavior of $\hat{\beta}(\tau)$ and $\hat{\Pi}$, whereas belloni2011l1 only had to dealt with one penalty. Second, the empirical processes in belloni2011l1 all involved independent data whereas as the proof of Theorem (ref) handles the time dependence of our model and the latent factors.

Interestingly, it is expected that the rate in Theorem (ref) is optimal. To elaborate on this point, consider the simple case where $n = T$, $\theta = 0$, $\tau = 0.5$, and $e_{i,t} := Y_{i,t} - \Pi_{i,t}(\tau)$ are mean zero i.i.d. sub-Gaussian($\sigma^2)$. The latter implies that \[ \mathbb{P}( \vert e_{1,1}\vert > z) \,\leq \, C_1 \exp\left( -\frac{z^2}{2 \sigma^2} \right),\,\,\,\, \] for a positive constant $C_1$, and for all $z >0$. Then by Theorem 2.3 in candes2009tight, we have the following lower bound for estimating $\Pi(\tau)$:

Notably, the lower bound in ((ref)) matches the rate implied by Theorem (ref), ignoring other factors depending on $s_{\tau}$, $\kappa_0$, $\kappa_m$, $p$ and $\phi_n$. However, we highlight that the upper bound ((ref)) in Theorem (ref) holds without the perhaps restrictive condition that the errors are sub-Gaussian.

We conclude this section with a result regarding the estimation of the factors and loadings of the latent matrix $\Pi(\tau)$. This is expressed in Corollary (ref) below and is immediate consequence of Theorem (ref) and Theorem 3 in yu2014useful.

A particularly interesting instance of Corollary (ref) is when $$(\sigma_1(\tau ))^2,\,(\sigma_{r_{\tau} -1}( \tau))^2 - (\sigma_{r_{\tau} }( \tau))^2\,\asymp \,n T, $$ a natural setting if the entries of $\Pi(\tau)$ are $O(1)$. Then the upper bound ((ref)) becomes \[ \underset{ O \in \mathbb{O}_{r_{\tau}} }{\min} \| \hat{g}(\tau) O - g(\tau)\|_F \,=\, O_{\mathbb{P}}\left( \frac{ \phi_n \sqrt{ c_T \,\log(p c_T \vee n ) }(\sqrt{1+s_{\tau}}+ \sqrt{\frac{r_{\tau}}{\log(p c_T \vee n ) }} ) }{\kappa_{0}^2 \underline{f}^{1/2} }\left( \frac{1}{\sqrt{n}} +\frac{1}{\sqrt{d_T}} \right ) \right), \] whereas ((ref)) is now \[ \frac{\| \hat{\lambda}(\tau) - \lambda(\tau)\|_F^2}{nT}\,=\,O_{\mathbb{P}}\left( \frac{ r_{\tau}\,\phi_{n}^2 c_T \log(p c_T \vee n )(1+s_{\tau}+ r_{\tau}/ \log(p c_T \vee n ))}{\kappa_{0}^4 \, \, \underline{f} }\left( \frac{1}{n} +\frac{1}{d_T} \right ) \right). \]

The conclusion of Corollary (ref) allows us to provide an upper bound on the estimation of factors ($g(\tau)$) and loadings ($\lambda(\tau)$) of the latent matrix $\Pi(\tau)$. Notice that we are not claiming that we provide consistent estimation of the number latent factors, as Theorem (ref) only guarantees consistent estimation of $\Pi(\tau)$. However, other authors, e.g. moon2015linear, have observed that estimation can be possible even if $r_{\tau}$ is unknown.

Nearly low rank quantiles

We conclude our theory section by studying the case where the matrix $X\theta(\tau)$ has nearley low rank. We make this formal by imposing the condition that $X\theta(\tau)$ can be perturbed into a low-rank matrix. While the result in this subsection differs in the choice of tuning parameters from those in Theorem (ref), our result here suggests that even in this low-rank setting there exists a choice of tuning parameters that allows our estimator to provide consistent estimation.

We view our setting below as an extension of the linear model in chernozhukov2018inference to the quantile framework and reduced rank regression. The specific condition is stated next.

Notice that in Assumption (ref), $\xi$ is an approximation error. In the case $\xi=0$, the condition implies that $\mathrm{rank}( X\theta(\tau) ) = O(r_{\tau})$ with probability close to one.

Next, exploiting Assumption (ref), we show that ((ref)) provides consistent estimation of the quantile function, namely, of $X \theta(\tau) + \Pi(\tau)$.

Interestingly, unlike Theorem (ref), Theorem (ref) does not show that we can estimate $\theta(\tau)$ and $\Pi(\tau)$ separately. Instead, we show that $\hat{\Pi}(\tau)$, the estimated matrix of latent factors, captures the overall contribution of both $\theta(\tau)$ and $\Pi(\tau)$. This is expected since Assumption (ref) states that, with high probability, $X \theta(\tau)$ has rank of the same order as of $\Pi(\tau)$. Notably, $\hat{\Pi}(\tau)$ is able to estimate $X \theta(\tau) + \Pi(\tau)$ via requiring that the value of $ \nu_{1}$ increases significantly with respect to the choice in Theorem (ref), while keeping $\nu_2 \asymp c_T( \sqrt{n} +\sqrt{d_T} )/(nT)$.

As for the convergence rate in Theorem (ref) for estimating $\Pi(\tau)$, this is of the order $r_{\tau} c_T(n^{-1} + d_T^{-1})$, if we ignore $\underline{f}$, $\overline{f}^{\prime}$, and $\phi_n$. When the data are independent, the rate becomes of the order $r_{\tau} ( n^{-1} + T^{-1} ) $. In such framework, our result matches the minimax rate of estimation in candes2010matrix for estimating an $n \times T$ matrix of rank $r_{\tau}$, provided that $n \asymp T$, see our discussion in Section (ref).

Simulation

In this section, we evaluate the performance of our proposed approach ($\ell_1$-NN-QR) with extensive numerical simulations focusing on the median case, namely the case when $\tau$ = 0.5. As benchmarks, we consider the $\ell_1$-penalized quantile regression studied in belloni2009computational, and similarly we refer to this procedure as $\ell_1$-QR. We also compare with the mean case, which we denote it as $\ell_1$-NN-LS as it combines the $\ell_2$-loss function with $\ell_1$ and nuclear norm regularization. We consider different generative scenarios. For each scenario we randomly generate 100 different data sets and compute the estimates of the methods for a grid of values of $ \nu_{1}$ and $ \nu_{2}$. Specifically, these tuning parameters are taken to satisfy $ \nu_{1} \in \{10^{-4}, 10^{-4.5},\ldots,10^{-8} \}$ and $ \nu_{2} \in \{10^{-3}, 10^{-4} ,\ldots,10^{-9} \} $. Given any choice of tuning parameters, we evaluate the performance of each competing method, averaging over the 100 data sets, and report values that correspond to the best performance. These are referred as optimal tuning parameters and can be thought of as oracle choices.

We also propose a modified Bayesian Information Criterion (BIC) to select the best pair of tuning parameters. Given a pair $( \nu_{1}, \nu_{2})$, our method produces a score $(\hat{\theta}(\tau),\hat{\Pi}(\tau))$. Specifically, denote $\hat{s}_{\tau} = \vert \{ j\,:\, \hat{\theta}_j (\tau) \neq 0 \}\vert$ and $\hat{r}_{\tau} = \mathrm{rank}(\hat{\Pi}(\tau))$,

where $c_1 >0$ is a constant. The intuition here is that the first term in the right hand side of ((ref)) corresponds to the fit to the data. The second term includes the factor $\log (nT) /2$ to emulate the usual penalization in BIC. The number of parameters in the model with choices $ \nu_{1}$ and $ \nu_{2}$ is estimated by $\hat{s}_{\tau}$ for the vector of coefficients, and $(1+n+T) \cdot \hat{r}_{\tau}$ for the latent matrix. The latter is reasonable since $\hat{\Pi}(\tau)$ is potentially a low rank matrix and we simply count the number of parameters in its singular value decomposition. As for the extra quantity $c_1$, we have included this term to balance the dominating contribution of the $(1+n+T) \cdot \hat{r}_{\tau} $. We find that in practice $c_1 = \log^2 (nT)$ gives reasonable performance in both simulated and real data. This is the choice that we use in our experiments. Then for each of data set under each design, we calculate the minimum value of $\mathrm{BIC}( \nu_{1}, \nu_{2})$, over the different choices of $ \nu_{1}$ and $ \nu_{2} $, and report the average over the 100 Monte Carlo simulations. We refer to this as BIC-$\ell_1$-NN-QR.

As performance measure we use a scaled version (see Tables (ref)-(ref)) of the squared distance between the true vector of coefficients $\theta$ and the corresponding estimate. We also consider a different metric, the “Quantile error" (koenker1999goodness):

which measures the average squared error between the quantile functions at the samples and their respective estimates. Since our simulations consider models with symmetric mean zero error, the above metric corresponds to the mean squared error for estimating the conditional expectation.

Next, we provide a detailed description of each of the generative models that we consider in our experiments. In each model design the dimensions of the problem are given by $ n \in \{100,500\}$, $p \in \{5,30\}$ and $T \in \{100,500\}$, and we also consider the instance $n=p=T=50$. The covariates $\{X_{i,t}\}$ are i.i.d $N(0,I_p)$.

Design 1. (Location shift model) The data is generated from the model

where $\sqrt{3} \epsilon_{i,t} \overset{i.i.d.}{\sim} t(3) $, $i = 1,\ldots,n$ and $t = 1,\ldots,T$, with $t(3)$ the Student's t-distribution with 3 degrees of freedom. The scaling factor $\sqrt{3}$ simply ensures that the errors have variance $1$. In ((ref)), we take the vector $\theta \in \mathbb{R}^p$ to satisfy \[ \theta_j =

\] We also construct $\Pi \in \mathbb{R}^{n \times T}$ to be rank one, defined as $ \Pi_{i,t} = 5i\,(\mathrm{cos}( 4\pi t/T ))/n.$

Design 2. (Location-scale shift model) We consider the model

where $\epsilon_{i,t} \overset{i.i.d.}{\sim} N(0,1) $, $i = 1,\ldots,n$ and $t = 1,\ldots,T$. The parameters in $\theta$ and $\Pi$ in ((ref)) are taken to be the same as in ((ref)). The only difference now is that we have the extra parameter $\theta \in \mathbb{R}^{p}$, which we define as $\theta_j = j/(2p)$ for $j \in \{1,\ldots,p\}$.

Design 3. (Location shift model with random factors) This is the same as Design 1 with the difference that we now generate $\Pi$ as

where

Design 4. (Location-scale shift model with random factors) This is a combination of Designs 2 and 3. Specifically, we generate data as in ((ref)) but with $\Pi$ satisfying ((ref)) and ((ref)).

The results in Tables (ref)-(ref) show a clear advantage of our proposed method against the benchmarks across the four designs we consider. This is true for estimating the vector of coefficients, and under the measure of quantile error. Importantly, our approach is not only the best under the optimal choice of tuning parameters but it remains competitive with the BIC type of criteria defined with the score ((ref)). In particular, under Designs 1 and 2, the data driven version of our estimator, BIC-$\ell_1$-NN-QR, performs very closely to the ideally tuned one $\ell_1$-NN-QR. In the more challenging settings of Designs 3 and 4, we noticed that BIC-$\ell_1$-NN-QR performs reasonably well compared to $\ell_1$-NN-QR.

Empirical Performance of the “Characteristics + Latent Factor” Model in Asset Pricing

Data Description

We use data from CRSP and Compustat to construct 24 firm level characteristics that are documented to explain the cross section and time series of stock returns in the finance and accounting literature. The characteristics we choose include well-known drivers of stock returns such as beta, size, book-to-market, momentum, volatility, liquidity, investment and profitability. Table (ref) in the Supplementary Material lists details of the characteristics used and the methods to construct the data. We follow the procedures of green2017characteristics to construct the characteristics of interest. The characteristics used in our model are standardized to have zero mean and unit variance. Figure (ref) plots the histogram of monthly stock returns and $9$ standardized firm characteristics. Each of them have different distribution patterns, suggesting the potential nonlinear relationship between returns and firm characteristics, which can be potentially captured by our quantile model.

Our empirical design is closely related to the characteristics model proposed by daniel1997evidence, daniel1998characteristics. To avoid any “data snooping” issue cause by grouping, we conduct the empirical analysis at individual stock level. Specifically, we use the sample period from January 2000 to December 2018, and estimate our model using monthly returns (228 months) from 1306 firms that have non-missing values during this period.

A “Characteristic + Latent Factor” Asset Pricing Model

We apply our model to fit the cross section and time series of stock returns (lettau2018estimating). There are $n$ assets (stocks), and the return of the each asset can potentially be explained by $p$ observed asset characteristics (sparse part) and $r$ latent factors (dense part). The asset characteristics are the covariates in our model. Our model imposes a sparse structure on the $p$ characteristics so that only the characteristics having the strongest explanatory powers are selected by the model. The part that's unexplained by the firm characteristics are captured by latent factors.

Suppose we have $n$ stock returns ($R_{1}$,...,$R_{n}$), and $p$ observed firm characteristics ($X_{1}$,...,$X_{p}$) over $T$ periods. The return quantile at level $\tau$ of portfolio $i$ in time $t$ is assumed to be the following: \[

\] where $X_{i,t-1,k}$, e.g. $k=1$ or $k=p$, is the $k$-th characteristic (e.g. the book-to-market ratio) of asset $i$ in time $t-1$. The coefficient $\theta_k$ captures the extent to which assets with higher/lower characteristic $X_{i,t,k}$ delivers higher average return. The term $g_{t}$ contains the $r_{\tau}$ latent factors in period $t$ which captures systematic risks in the market, and $\lambda_{i}$ contains portfolio $i$'s loading on these factors (i.e. exposure to risk).

There is a discussion in academic research on “factor versus characteristics” in late 1990s and early 2000s. The factor/risk based view argues that an asset has higher expected returns because of its exposure to risk factors (e.g. Fama-French 3 factors) which represent some unobserved systematic risk. An asset's exposure to risk factors are measured by factor loadings. The characteristics view claims that stocks have higher expected returns simply because they have certain characteristics (e.g. higher book-to-market ratios, smaller market capitalization), which might be independent of systematic risk (daniel1997evidence,daniel1998characteristics). The formulation of our model accommodates both the factor view and the characteristics view. The sparse part is similar to daniel1997evidence,daniel1998characteristics, in which stock return are explained by firm characteristics. The dense part assumes a low-dimensional latent factor structure where the common variations in stock returns are driven by several “risk factors”.

Empirical Results

We first get the estimates $\hat{\theta}(\tau)$ and $\hat{\Pi}(\tau)$ at three different quantiles, $\tau = \{0.1, 0.5, 0.9\}$ using our proposed ADMM algorithm. We then decompose $\hat{\Pi}(\tau)$ into the products of its $\hat{r}_\tau$ principal components $\hat{g}(\tau)$ and their loadings $\hat{\lambda}(\tau)$ via eq((ref)). The $(i,k)$-th element of $\hat{\lambda}(\tau)$, denoted as $\hat{\lambda}_{i,k}(\tau)$, can be interpreted as the exposure of asset $i$ to the $k$-th latent factor (or in finance terminology, “quantity of risk”). And the $(t,k)$-th elements of $\hat{g}(\tau)$, denoted as $\hat{g}_{t,k}(\tau)$, can be interpreted as the compensation of the risk exposure to the $k$-th latent factor in time period $t$ (or in finance terminology, “price of risk”). The model are estimated with different tuning parameters $ \nu_{1}$ and $ \nu_{2}$, and use our proposed BIC to select the optimal tuning parameters. The details of the information criteria can be found in equation ((ref)).

The tuning parameter $ \nu_{1}$ governs the sparsity of the coefficient vector $\theta$. The larger $ \nu_{1}$ is, the larger the shrinkage effect on $\theta$. Figure (ref) illustrate the effect of this shrinkage. With $\nu_2$ fixed, as the value of $ \nu_{1}$ increases, more coefficients in the estimated $\theta$ vector shrink to zero. From a statistical point of view, the “effective characteristics” that can explain stock returns are those with non-zero coefficient $\theta$ at relatively large values of $ \nu_{1}$.

Table (ref) reports the relationship between tuning parameter $ \nu_{2}$ and rank of estimated $\Pi$ at different quantiles. It shows that the tuning parameter $ \nu_{2}$ governs the rank of matrix $\Pi$, and that as $ \nu_{2}$ increases, we penalize more on the rank of matrix $\Pi$ through its nuclear norm.

The left panel of Table (ref) reports the estimated coefficients in the sparse part when we fix the tuning parameters at $\log_{10}( \nu_{1})=-3.5$ and $\log_{10}( \nu_{2})=-4$. The signs of some characteristics are the same across the quantiles, e.g. size (mve), book-to-market (bm), momentum (mom1m, mom12m), accurals (acc), book equity growth (egr), leverage (lev), and standardized unexpected earnings (sue). However, some characteristics have heterogenous effects on future returns at different quantiles. For example, at the 10% quantile, high beta stocks have high future returns, which is consistent with results found via the CAPM; while at $50\%$ and 90% quantile, high beta stocks have low future returns, which conforms the “low beta anomaly” phenomenon. Volatility (measured by both range and idiosyncratic volatility) is positively correlated with future returns at 90% quantile, but negatively correlated with future returns at 10% and 50% percentile. The result suggests that quantile models can capture a wider picture of the heterogenous relationship between asset returns and firm characteristics at different parts of the distribution (koenker2000galton).

Table (ref) reports the selected optimal tuning parameters $ \nu_{1}$ and $ \nu_{2}$ for different quantiles. The tuning parameters are selected via BIC based on ((ref)) as discussed in Section (ref). For every $ \nu_{1}$ and $ \nu_{2}$, we get the estimates $\widetilde{\theta}( \nu_{1}, \nu_{2})$ and $\widetilde{\Pi}( \nu_{1}, \nu_{2})$ and the number of factors $r=\text{rank}(\widetilde{\Pi}( \nu_{1}, \nu_{2}))$. The $\theta$ vector is sparse with non-zero coefficients on selected characteristics. The 10% quantile of returns has only 1 latent factor, and 3 selected characteristics. The median of returns has 7 latent factors and 2 selected characteristics. The 90% quantile of returns has 2 latent factors and 7 selected characteristics. Range is the only characteristic selected across all 3 quantiles. Idiosyncratic volatility is selected at 10% and 90% quantiles, with opposite signs. 1-month momentum is selected at 50% and 90% percentiles, with negative sign suggesting reversal in returns.

Overall, the empirical evidence suggests that both firm characteristics and latent risk factors have valuable information in explaining stock returns. In addition, we find that the selected characteristics and number of latent factors differ across the quantiles.

Interpretation of Latent Factors

Table (ref) below reports the variance in the matrix $\Pi$ explained by each Principal Component (PC) or latent factor. At upper and lower quantiles, the first PC dominates. At the median there are more latent factors accounting for the variations in $\Pi$, with second PC explaining 13.8% and third PC explaining 6.8%.

We also found the first PC captures the market returns in all three quantiles: Figure (ref) plots the first principal component against the monthly returns of S&P500 index, showing that they have strong positive correlations.

Acknowledgements

We would like to thank the Editors, the Associate Editors and two anonymous referees for their detailed reviews, which helped to improve the paper substantially. We are also grateful to Victor Chernozhukov, Iv\'{a}n Fern\'{a}ndez-Val, Bryan Graham, Hiroaki Kaido, Anna Mikusheva, Whitney Newey, Eric Renault, Jeremy Smith, and Vasilis Syrgkanis for helpful discussions.

\setcounter{section}{0} \setcounter{figure}{0} \setcounter{table}{0} \setcounter{equation}{0}

\setcounter{page}{1}

Implementation Details of the Proposed ADMM Algorithm

Denoting by $P_{+}(\cdot)$ and $P_{-}(\cdot)$ the element-wise positive and negative part operators, the ADMM proceeds doing the iterative updates

where $\eta>0$ is the penalty, see boyd2011distributed.

The update for $\tilde{\theta}$ is \[ \displaystyle \tilde{\theta}^{(k+1)} \,\leftarrow \, \left[ \sum_{i=1}^n\sum_{t=1}^T X_{i,t} X_{i,t}^{\prime} +I_p \right]^{-1} \left[ - \sum_{ i= 1}^n \sum_{t=1}^T X_{i,t}A_{i,t} +Z_{\theta}^{(k)} + U_{\theta}^{(k)} \right], \] where \[ A := W^{(k)}+ Z_{\Pi}^{(k)} + U_W^{(k)} - Y. \] The update for $\tilde{\Pi}$ is \[ \tilde{\Pi}^{(k+1)} \,\leftarrow\, P \text{ diag}\left( \max\left\{ 0, v_j - \frac{ \nu_2}{\eta} \right\}_{ 1\leq j \leq l} \right) Q^{\prime}, \] where \[ Z_{\Pi}^{(k)} + U_{ \Pi }^{(k)} \,=\, P \text{ diag}( \{ v_j \}_{1\leq j \leq l} ) Q^{\prime}. \] Furthermore, for $Z_{\theta}$, \[ Z_{\theta,j}^{(k+1)} \,\leftarrow \, \text{sign}(\tilde{\theta}_j^{(k+1)} - U_{\theta,j}^{(k)} ) \left[ \vert \tilde{\theta}_j^{(k+1)} - U_{\theta,j}^{(k)} \vert - \frac{ \nu_{1} w_j }{\eta} \right]. \] Finally, defining \[ \tilde{A} \,=\, - Y + X\tilde{\theta}^{(k+1)} + U_W^{(k)}, \,\,\,\, \tilde{B} \,=\,-V^{(k+1)}-U_V^{(k)},\,\,\,\, \tilde{C} \,=\, -\tilde{ \Pi}^{(k+1)} + U_{\Pi}^{(k)}, \] the remaining updates are \[ Z_{\Pi}^{(k+1)} \,\leftarrow\, \frac{ -\tilde{A} - 2\tilde{C} + \tilde{B} }{3}, \] and \[ W^{(k+1)} \,\leftarrow\, -\tilde{A} -\tilde{C} - 2Z_{\Pi}^{(k+1)}. \]

Estimation without Covariates

Note, when there are no covariates, our proposed ADMM can be simplified. In this case, we face the following problem

This can be thought as a convex relaxation of the estimator studied in chen2018quantile. Problem ((ref)) is also related to the setting of robust estimation of a latent low-rank matrix, e.g. elsener2018robust. However, our approach can also be used to estimate different quantile levels. As for solving ((ref)), we can proceed by doing the iterative updates

and

where $\eta> 0 $ is the penalty parameter (boyd2011distributed). The minimization in ((ref)) is similar to ((ref)), whereas ((ref)) can be done similarly as in ((ref)).

Estimation in unbalanced designs

We now explore the setting of unbalanced designs, specifically, instead of fully observing $(X,Y)$, we now assume that we only observe $\{ (X_{i,t},Y_{i,t}) \}_{(i,t) \in \mathcal{I}}$ for a set $\mathcal{I}\subset [n] \times [T]$. For instance, if $n=2$ and $T=3$ but for $i=1$ the data at time $t=2$ is not available, then $\mathcal{I}$ would be $\{ (1,1),(1,3),(2,1),(2,2),(2,3)\}$.

For this setting, inspired by ((ref)), we formulate the problem

where $ \nu_{1} >0$ and $\nu_2 >0$ are tuning parameters, and $w_1,\ldots,w_p$ are user specified weights. Thus, comparing with ((ref)), we now only apply the loss function to the indices for which there is data available. Following ((ref)), we write ((ref)) as

Then, as in Section (ref), we obtain the ADMM updates given by \[ V_{i,t}^{(k+1)}\,\leftarrow\,P_{+}\left(W_{i,t}^{(k)}-(U_{V})_{i,t}^{(k)}-\frac{\tau}{\vert \mathcal{I}\vert \eta}\right)+P_{-}\left(W_{i,t}^{(k)}-(U_{V})_{i,t}^{(k)}-\frac{\tau}{ \vert \mathcal{I}\vert \eta}\right)\\ \label{eqn:iter4} \] for $(i,t) \in \mathcal{I}$, and \[ V_{i,t}^{(k+1)}\,\leftarrow\, W_{i,t} - (U_{V})_{i,t}^{(k)} \] for $(i,t) \notin \mathcal{I}$, and with the rest of updates given exactly as in Section (ref).

Notice that our formulation is similar in spirit to athey2018matrix. In particular, the objective function ((ref)) resembles Equation (4.3) in athey2018matrix, with the main difference that we allow for covariates and work with the quantile loss.

Proof of Theorem (ref)

Auxiliary lemmas for proof of Theorem (ref)

Throughout, we use the notation \[ Q_{\tau}(\tilde{\theta},\tilde{\Pi}) \,=\, \mathbb{E}(\hat{Q}_{\tau}(\tilde{\theta},\tilde{\Pi}) ). \] Moreover, as in yu1994rates, we define the sequence $\{(\tilde{Y}_{i,t}, \tilde{X}_{i,t} ) \}_{i \in [n], t \in [T]}$ such that

Here, we define $ \Lambda \,:=\, \{ H_1,H_1^{\prime},\ldots,H_{d_T},H_{d_T}^{\prime},R \}$ with

We also use the symbol $\mathcal{L}(\cdot)$ to denote the distribution of a sequence of random variables.

Next, define the scores $a_{i,t} \,=\, \tau - 1\{ Y_{i,t} \leq X_{i,t}^{\prime} \theta(\tau) + \Pi_{i,t}(\tau) \} $, and $\tilde{a}_{i,t} \,=\, \tau - 1\{ \widetilde{Y}_{i,t} \leq \widetilde{X}_{i,t}^{\prime} \theta(\tau) + \Pi_{i,t}(\tau) \} $.

We start by controlling an empirical process involving the scores $\{a_{i,t}\}$. This is given next.

Next we proceed to control the complexity of the set $\{ \Delta \in \mathbb{R}^{ n \times T } \,:\, \| \Delta \|_* \leq 1\,\,\, \}$ in terms of the scores $\{a_{i,t}\}$.

We will proceed to exploit Lemmas (ref) and (ref) to show that $(\hat{\beta}(\tau) -\beta(\tau), \hat{\Pi}(\tau) -\Pi(\tau) )$ belongs to the restricted set with probability approaching one. Before that, we recall an important property relating the nuclear norm to the rank of a matrix.

Our next result shows how the function $ Q_{\tau}$ changes locally in the restricted set around $\theta(\tau), \Pi(\tau)$.

Next we study the behavior of the $\ell_2$ and nuclear norms in the restricted set $A_{\tau}$.

Using all the previous lemmas our next results provides the control of the empirical process associated with our estimator. Control of this empirical process leads to the convergence rates obtained in the paper.

Combining all the previous lemmas we prove Theorem (ref) in the next subsection.

Proof of Theorem (ref)

Proof of Theorem (ref)

Auxiliary lemmas for proof of Theorem (ref)

Proof of Theorem (ref)

The proof of Theorem (ref) proceeds by exploiting Lemmas (ref) and (ref). By Lemma (ref), we have that

with high probability. Therefore, we assume that ((ref)) holds. Hence, if $ \hat{\Delta} \notin A_{\tau}^{\prime \prime}$, then

If $ \hat{\Delta} \in A_{\tau}^{\prime \prime}$, then we proceed as in the proof of Theorem (ref) by exploiting Lemma (ref), and treating $X\theta(\tau) + \Pi(\tau)$ as the latent factors matrix, the design matrix as the matrix zero, $A_{\tau}^{\prime \prime}$ as $A_{\tau}$, and \[ \kappa_0 \,=\, \underline{f}^{1/2}, \] in Assumption (ref). This leads to

and the claim in Theorem (ref) follows combining ((ref)) and ((ref)).

Proof of Corollary (ref)

First notice that by Theorem (ref) and Theorem 3 in yu2014useful,

Furthermore, \[

\] where the third inequality follows from Weyl's inequality, and the last one from ((ref)).

Proof of Theorem (ref)

Conditioning on $\Pi$ and $\{X_{i,t}\}$, as in yu1994rates, we define the sequence $\{\tilde{\varepsilon}_{i,t} \}_{i \in [n], t \in [T]}$ such that

Here, we define $ \Lambda \,:=\, \{ H_1,H_1^{\prime},\ldots,H_{d_T},H_{d_T}^{\prime},R \}$ with

We then let \[ \tilde{Y}_{i,t}\, = \, X_{i,t}^{\prime} \theta(\tau)+\Pi_{i,t}(\tau) + \tilde{\varepsilon}_{i,t} G(X_{i,t}) \] for all $i\in [n]$ and $t\in [T]$.

Furthermore, we define the scores $a_{i,t} \,=\, \tau - 1\{ Y_{i,t} \leq X_{i,t}^{\prime} \theta(\tau) + \Pi_{i,t}(\tau) \} $, and $\tilde{a}_{i,t} \,=\, \tau - 1\{ \widetilde{Y}_{i,t} \leq X_{i,t}^{\prime} \theta(\tau) + \Pi_{i,t}(\tau) \} $.

Auxiliary lemmas for proof of Theorem (ref)

Throughout we use the notation from Section (ref) and $(\hat{\beta}(\tau), \hat{\Pi}(\tau) )$ denotes the estimator defined in ((ref)).

We begin by defining \[ \displaystyle \hat{M}( \tilde{\beta},\tilde{ \Pi} ) \,=\, \frac{1}{nT}\sum_{i=1}^n \sum_{t=1}^{T} \left[ \rho_{\tau}(Y_{i,t} - X_{i,t}^{\prime} \tilde{\beta} - \tilde{ \Pi}_{i,t} ) - \rho_{\tau}(Y_{i,t} - X_{i,t}^{\prime} \beta(\tau) - \Pi_{i,t}(\tau) ) \right] \] and \[ \displaystyle M(\tilde{\beta},\tilde{ \Pi}) \,=\, \mathbb{E}\left( \hat{M}( \tilde{\beta},\tilde{ \Pi} ) \big| \{ X_{i,t}\} \right). \] We also set

with the notation in Assumption 2.

We start by recalling a result from padilla2020adaptive involving the behavior of $M$ locally around the true quantiles.

We now proceed to construct a restricted set $K$ where the solution $(\hat{\beta}(\tau),\hat{\Pi}(\tau)) $ lies with high probability.

Next we combine the previous two results to arrive at an upper bound on the estimation error of the quantiles.

We now proceed to give an upper bound on the second term in the right hand side of ((ref)).

We now proceed to control the first term in the right hand side of ((ref)).

Proof of Theorem (ref)

The proof follows from Lemmas (ref)--(ref) by setting \[ \eta\,\asymp\, m_n \left[\frac{ \phi_{n,T} \sqrt{c_T\log p} }{n\sqrt{ d_T}} \, +\,\frac{\psi_{n,T} ( \sqrt{n} +1)}{ nd_T }\right], \] for any sequence $m_n$ satisfying $m_{n} \rightarrow \infty$.

Note: Estimated under different values of turning parameter $ \nu_{2}$, when $ \nu_{1}=10^{-5}$ is fixed. The results are reported for quantiles 10%, 50% and 90%.