Quantile index regression

Consider a response $Y$ and a $p$-dimensional vector of covariates $\boldsymbol{X}=(X_{1},...,X_{p})^{\prime}$. We then rewrite the quantile function of $Y$ conditional on $\boldsymbol{X}$ at (1.1) with an explicit form of $\boldsymbol{\theta}(\boldsymbol{X},\boldsymbol{\beta})$ below,

where $\mathcal{I}\subset[0,1]$ is an interval or the union of multiple disjoint intervals, $\boldsymbol{\theta}(\boldsymbol{X},\boldsymbol{\beta})=(\theta_{1}(\boldsymbol{X},\boldsymbol{\beta}),\cdots,\theta_{d}(\boldsymbol{X},\boldsymbol{\beta}))^{\prime}$, $\boldsymbol{\beta}=(\boldsymbol{\beta}_{1}^{\prime},...,\boldsymbol{\beta}_{d}^{\prime})^{\prime}$, $\boldsymbol{\beta}_{j}=(\beta_{j1},\cdots,\beta_{jp})^{\prime}$, $\theta_{j}(\boldsymbol{X},\boldsymbol{\beta})=g_{j}(\boldsymbol{X}^{\prime}\boldsymbol{\beta}_{j})$, the link functions $g_{j}^{-1}(\cdot)$s are monotonic for $1\leq j\leq d$, and the intercept term can be included by letting $X_{1}=1$. We call model (2.1) the quantile index regression (QIR) for simplicity, and the two examples of $Q(\cdot,\cdot)$ below are first introduced to illustrate the new model.

In the literature, there exist many quantile functions, which have explicit forms, such as the generalized lambda and Burr XII distributions (Gilchrist,, 2000; Fournier et al.,, 2007). For example, the generalized lambda distribution has the form of

where the indices $\theta_{3}$ and $\theta_{4}$ control the right and left tails, respectively. Note that it reduces to the Tukey lambda distribution when $\theta_{3}=\theta_{4}$. This indicates that the generalized lambda distribution can be considered if we focus on the quantiles with the full range, i.e. $\mathcal{I}=[0,1]$. On the other hand, the Tukey lambda may be a better choice if our interest is on the quantiles at one side only, i.e. $\mathcal{I}\subset[0,0.5)$ or $\mathcal{I}\subset(0.5,1]$.

Denote the observed data by $\{(Y_{i},\boldsymbol{X}_{i}^{\prime})^{\prime},i=1,...,n\}$, and they are independent and identically distributed ($i.i.d.$) samples of random vector $(Y,\boldsymbol{X})$, where $Y_{i}$ is the response, $\boldsymbol{X}_{i}=(X_{i1},...,X_{ip})^{\prime}$ contains $p$ covariates, and $n$ is the number of observations.

Let $\tau_{1}\leq\tau_{2}\leq\ldots\leq\tau_{K}$ be $K$ fixed quantile levels, where $\tau_{k}\in\mathcal{I}$ for all $1\leq k\leq K$. To achieve higher efficiency, we consider the composite quantile regression (CQR) estimator below.

where $\rho_{\tau}(x)=x\{\tau-I(x<0)\}$ is the quantile check function; see Zou and Yuan, (2008) and Kai et al., (2010). Note that $Q(\tau,\boldsymbol{\theta}(\mathbf{X},\widehat{\boldsymbol{\beta}}_{n}))$ has the non-crossing property with respect to $\tau$ since, for each $\boldsymbol{\theta}$, $Q(\tau,\boldsymbol{\theta})$ is a non-decreasing quantile function.

To study the asymptotic properties of $\widehat{\boldsymbol{\beta}}_{n}$, we first consider $\boldsymbol{\beta}_{0}=(\boldsymbol{\beta}_{01}^{\prime},...,\boldsymbol{\beta}_{0d}^{\prime})^{\prime}=\arg\min_{\boldsymbol{\beta}}\bar{L}(\boldsymbol{\beta})$, where $\bar{L}(\boldsymbol{\beta})=E[\sum_{k=1}^{K}\rho_{\tau_{k}}\{Y-Q(\tau_{k},\boldsymbol{\theta}(\mathbf{X},\boldsymbol{\beta}))\}]$ is the population loss function, and it may not be unique, i.e. the CQR estimation at (2.2) may suffer from the identification problem. For the sake of illustration, let us consider the case without including covariates $\mathbf{X}$, i.e. we estimate $\boldsymbol{\theta}$ with a sequence of quantile levels $\tau_{k}\in\mathcal{I}$ and $1\leq k\leq K$. To this end, it requires that two different values of $\boldsymbol{\theta}$ can not yield the same quantile function $Q(\tau_{k},\boldsymbol{\theta})$ across all $K$ levels. In other words, if there exists $\boldsymbol{\theta}\neq\boldsymbol{\theta}^{*}$ that yield $Q(\tau_{k},\boldsymbol{\theta})=Q(\tau_{k},\boldsymbol{\theta}^{*})$ for all $K$ quantiles, then $\boldsymbol{\theta}$ and $\boldsymbol{\theta}^{*}$ are not identifiable. In sum, to guarantee that $\boldsymbol{\beta}_{0}$ is the unique minimizer of the population loss, we make the following assumption on the quantile function $Q(\tau,\boldsymbol{\theta})$.

Intuitively, for any quantile function $Q(\tau,\boldsymbol{\theta})$, one can always increase the number of quantile levels $K$ to make Assumption 1 hold. However, it may also depend on the structures of quantile functions and the number of indices. For the sake of illustration, we state the sufficient and necessary condition for the Tukey lambda distribution that satisfies Assumption 1.

Assumption 1, together with an additional assumption on covariates $\mathbf{X}$, allows us to show that $\boldsymbol{\beta}_{0}$ is the unique minimizer of $\bar{L}(\boldsymbol{\beta})$; see the following theorem, which is critical to establish the asymptotic properties of $\widehat{\boldsymbol{\beta}}_{n}$.

To demonstrate the consistency of $\widehat{\boldsymbol{\beta}}_{n}$ given below, we assume that the parameter space $\Theta\subset\mathbb{R}^{dp}$ is compact, and the true parameter vector $\boldsymbol{\beta}_{0}$ is an interior point of $\Theta$.

Denote

and

Note that the objective function $L_{n}(\boldsymbol{\beta})$ is non-convex and non-differentiable, and this makes it challenging to establish the asymptotic normality of $\widehat{\boldsymbol{\beta}}_{n}$. We overcome the difficulty by making use of the bracketing method in Pollard, (1985). Moreover, to estimate the asymptotic variance in Theorem 3, we first apply the nonparametric method in Hendricks and Koenker, (1991) to estimate the quantities of $f_{Y}\left\{Q(\tau_{k},\boldsymbol{\theta}(\mathbf{X},\boldsymbol{\beta}_{0}))|\mathbf{X}\right\}$ with $1\leq k\leq K$, and then matrices $\Omega_{0}$ and $\Omega_{1}$ can be approximated by the sample averaging.

In addition, based on the estimator $\widehat{\boldsymbol{\beta}}_{n}$, one can use $Q(\tau^{*},\boldsymbol{\theta}(\mathbf{X},\widehat{\boldsymbol{\beta}}_{n}))$ to further predict the conditional quantile at any level $\tau^{*}\in\mathcal{I}$, and the corresponding theoretical justification can be established by directly applying the delta-method (van der Vaart,, 1998, Chapter 3).

This subsection considers the case with high-dimensional covariates, i.e., $p\gg n$, and the true parameter vector $\boldsymbol{\beta}_{0}$ is assumed to be $s$-sparse, i.e. the number of nonzero elements in $\boldsymbol{\beta}_{0}$ is no more than $s>0$. A regularized CQR estimation can then be introduced,

where $\Theta$ is given in Theorem 4, $p_{\lambda}$ is a penalty function, and it depends on a tuning (regularization) parameter $\lambda\in\mathbb{R}^{+}$ with $\mathbb{R}^{+}=(0,\infty)$.

Consider the loss function $L_{n}(\boldsymbol{\beta})=\sum_{k=1}^{K}\sum_{i=1}^{n}\rho_{\tau_{k}}\{Y_{i}-Q(\tau_{k},\boldsymbol{\theta}(\mathbf{X}_{i},\boldsymbol{\beta}))\}$ defined in (2.2), and $Q(\tau_{k},\boldsymbol{\theta}(\mathbf{X}_{i},\boldsymbol{\beta}))$ usually is nonconvex with respect to $\boldsymbol{\beta}$. As a result, $L_{n}(\boldsymbol{\beta})$ will be nonconvex although the check loss $\rho_{\tau}(\cdot)$ is convex, and there is no more harm to use nonconvex penalty functions. Specifically, we consider the component-wise penalization,

where $p_{\lambda}(\cdot)$ is possibly nonconvex and satisfies the following assumption.

The above is the $\mu$-amenable assumption given in Loh and Wainwright, (2015) and Loh, (2017), and the penalty function is required not too far from the convexity. Note that the popular penalty functions, including SCAD (Fan and Li.,, 2001) and MCP (Zhang,, 2010), satisfy the above properties.

In the literature of nonconvex penalized quantile regression, Jiang et al., (2012) studied nonlinear quantile regressions with SCAD regularizer from the asymptotic viewpoint, while it can only handle the case with $p=o(n^{1/3})$. Wang et al., 2012b and Sherwood et al., (2016) considered the case that $p$ grows exponentially with $n$, and their proving techniques heavily depend on the condition that the loss function should be represented as a difference of the two convex functions. However, $L_{n}(\boldsymbol{\beta})$ does not meet this requirement since quantile function $Q(\tau,\boldsymbol{\theta})$ can be nonconvex.

On the other hand, non-asymptotic properties recently have attracted considerable attention in the theories of high-dimensional analysis; see, e.g., Belloni and Chernozhukov, (2011); Sivakumar and Banerjee, (2017); Pan and Zhou, (2021). This subsection attempts to study them for our proposed quantile estimators, while it is a nontrival task since existing results only focused on linear quantile regression. Loh and Wainwright, (2015) and Loh, (2017) studied the non-asymptotic properties for M-estimators with both nonconvex loss and regularizers, while they required the loss function to be twice differentiable. The technical proofs in the Appendix follow the framework in Loh and Wainwright, (2015) and Loh, (2017), and some new techniques are developed to tackle the nondifferentiability of the quantile check function.

Let $\boldsymbol{\theta}(\boldsymbol{\gamma})=(g_{1}(\gamma_{1}),\ldots,g_{d}(\gamma_{d}))$ with $g_{j}^{-1}(\cdot)$s being link functions and $\boldsymbol{\gamma}=(\gamma_{1},\ldots,\gamma_{d})$, and we can then denote $Q(\tau,\boldsymbol{\gamma}):=Q(\tau,\boldsymbol{\theta}(\boldsymbol{\gamma}))$. Moreover, by letting $\gamma_{j}(\mathbf{X},\boldsymbol{\beta})=\mathbf{X}\boldsymbol{\beta}_{j}$ for $1\leq j\leq d$ and $\boldsymbol{\gamma}(\mathbf{X},\boldsymbol{\beta})=(\gamma_{1}(\mathbf{X},\boldsymbol{\beta}),\ldots,\gamma_{d}(\mathbf{X},\boldsymbol{\beta}))$, we can further denote $Q(\tau,\boldsymbol{\gamma}(\mathbf{X},\boldsymbol{\beta})):=Q(\tau,\boldsymbol{\theta}(\mathbf{X},\boldsymbol{\beta}))$.

The differentiable assumption of quantile functions allows us to use the Lipschitz property and multivariate contraction theorem. The boundedness of covariates is to assure that the bounded difference inequality can be used, and it can be relaxed with more complicated and lengthy technical details (Wang and He,, 2021).

Denote by $\mathcal{B}_{R}(\boldsymbol{\beta}_{0})=\{\boldsymbol{\beta}\in\mathbb{R}^{dp}:\|\boldsymbol{\beta}-\boldsymbol{\beta}_{0}\|\leq R\}$ the Euclidean ball centered at $\boldsymbol{\beta}_{0}$ with radius $R>0$, and let $\lambda_{\min}(\boldsymbol{\beta})$ be the smallest eigenvalue of matrix

The above assumption guarantees that the population loss $\bar{L}(\boldsymbol{\beta})=E[n^{-1}L_{n}(\boldsymbol{\beta})]$ is strongly convex around the true parameter vector $\boldsymbol{\beta}_{0}$. Specifically, let $\bar{\mathcal{E}}(\Delta)=\bar{L}(\boldsymbol{\beta}_{0}+\Delta)-\bar{L}(\boldsymbol{\beta}_{0})-\Delta^{\prime}{\partial\bar{L}(\boldsymbol{\beta}_{0})}/{\partial\boldsymbol{\beta}}$ be the first-order Taylor expansion. Then, by Assumption 6, we have that $\bar{\mathcal{E}}(\Delta)\geq 0.5f_{\min}\lambda_{\min}^{0}\|\Delta\|^{2}$ for all $\Delta$ such that $\|\Delta\|\leq R$; see Lemma A.5 in the Appendix for details. We next obtain the non-asymptotic estimation bound of $\widetilde{\boldsymbol{\beta}}_{n}$.

In practice, we can choose $\lambda\asymp\sqrt{\log p/n}$, and it then holds that $\|\widetilde{\boldsymbol{\beta}}_{n}-\boldsymbol{\beta}_{0}\|\lesssim\sqrt{s\log p/n}$, which has the standard rate of error bounds; see, e.g., (Loh,, 2017). Moreover, for the predicted conditional quantile of $Q(\tau^{*},\boldsymbol{\theta}(\mathbf{X},\widetilde{\boldsymbol{\beta}}_{n}))$ at any level $\tau^{*}\in\mathcal{I}$, it can be readily verified that $|Q(\tau^{*},\boldsymbol{\theta}(\mathbf{X},\widetilde{\boldsymbol{\beta}}_{n}))-Q(\tau^{*},\boldsymbol{\theta}(\mathbf{X},{\boldsymbol{\beta}}_{0}))|$ has the same convergence rate as $\|\widetilde{\boldsymbol{\beta}}_{n}-\boldsymbol{\beta}_{0}\|$. Finally, the above theorem requires the minimization (2.3) to be conducted in $\Theta=\mathcal{B}_{R}(\boldsymbol{\beta}_{0})$, which is unknown but fixed. This enables us to solve the problem by conducting a random initialization in optimizing algorithms.

This subsection provides algorithms to search for the CQR estimator at (2.2) and regularized estimator at (2.3).

For the CQR estimation without penalty at (2.2), we employ the commonly used gradient descent algorithm to search for estimators, and the $(r+1)$th update is given by

where $\widehat{\boldsymbol{\beta}}_{n}^{(r)}$ is from the $r$th iteration, and $\eta^{(r)}$ is the step size. Note that the quantile check loss is nondifferentiable at zero, and $\nabla L_{n}({\boldsymbol{\beta}}^{(r)})$ in the above refers to the subgradient (Moon et al.,, 2021) instead. In practice, too small step size will cause the algorithm to converge slowly, while too large step size may cause the algorithm to diverge. We choose the step size by the backtracking line search (BLS) method, which is shown to be simple and effective; see Bertsekas, (2016). Specifically, the algorithm starts with a large step size and, at $(r+1)$th update, it is reduced by keeping multiplying a fraction of $b$ until $L_{n}(\boldsymbol{\beta}^{(r+1)})-L_{n}(\boldsymbol{\beta}^{(r)})<-a\eta^{(r)}\|\nabla L_{n}(\boldsymbol{\beta}^{(r)})\|_{2}^{2}$, where $a$ is another hyper-parameter. The simulation experiments in Section 4 work well with the setting of $(a,b)=(0.3,0.5)$.

For the regularized estimation at (2.3), we adopt the composite gradient descent algorithm (Loh and Wainwright,, 2015), which is designed for a nonconvex problem and fits our objective functions well. Consider the SCAD penalty, which satisfies Assumption 4 with $L=1$ and $\mu=1/(\alpha-1)$. We then can rewrite the optimization problem at (2.3) into

where, from Assumption 4, $g(\boldsymbol{\beta})=\{\sum_{j=1}^{d}p_{\lambda}(\boldsymbol{\beta}_{j})+\mu\|\boldsymbol{\beta}\|_{2}^{2}/2\}/\lambda$ is convex. As a result, similar to the composite gradient descent algorithm in Loh and Wainwright, (2015), the $(r+1)$th update can be calculated by

which has a closed-form solution of

with $z=(\boldsymbol{\beta}^{(r)}-\eta\nabla\widetilde{L}_{n}(\boldsymbol{\beta}^{(r)}))/(1+\mu\eta)$ and $\nu=\eta/(1+\mu\eta)$, where the step size $\eta$ is chosen by the BLS method.

There are two types of hyper-parameters in the penalized estimation at (2.3): the tuning parameter $\lambda$ and quantile levels of $\tau_{k}$ with $1\leq k\leq K$. We can employ validation methods to select the tuning parameter $\lambda$ such that the composite quantile check loss is minimized.

The selection of $\tau_{k}$’s is another important task since it will affect the efficiency of resulting estimators. Suppose that we are interested in some high quantiles of $\tau_{m}^{*}$ with $1\leq m\leq M$, and then the QIR model can be assumed to the interval of $\mathcal{I}=[\tau_{0},1]$, which contains all $\tau_{m}^{*}$’s. We may further choose a suitable interval of $[\tau_{L},\tau_{U}]\subset\mathcal{I}$ such that $\tau_{k}$’s can be equally spaced on it, i.e. $\tau_{k}=\tau_{L}+k(\tau_{U}-\tau_{L})/(K+1)$ for $1\leq k\leq K$, where it can be set to $\tau_{0}=\tau_{L}$. As a result, the selection of $\tau_{k}$’s is equivalent to that of $[\tau_{L},\tau_{U}]$.

We may choose $\tau_{U}$ such that it is close to $\tau_{m}^{*}$’s, while a reliable estimation can be afforded at this level. The selection of $\tau_{L}$ is a trade-off between estimation efficiency and model misspecification; see Wang et al., 2012a ; Wang and Tsai, (2009). On one hand, to improve estimation efficiency, we may choose $\tau_{L}$ close to 0.5 since the richest observations will appear at the middle for most real data. On the other hand, we have to assume the parametric structure over the whole interval of $[\tau_{L},1]$, i.e. more limitations will be added to the real example. The criterion of prediction errors (PEs) is hence introduced,

where we will choose $\tau_{L}$ with the minimum value of PEs; see also Wang et al., 2012a .

In practice, the cross validation method can be used to select $\lambda$ and $\tau_{L}$ simultaneously. Specifically, the composite quantile check loss and PEs are both evaluated at validation sets. For each candidate interval of $[\tau_{L},\tau_{U}]$, the tuning parameter $\lambda$ is selected according to the composite quantile check loss, and the corresponding value of PE is also recorded. We then will choose the value of $\tau_{L}$, which corresponds to the minimum value of PEs among all candidate intervals.

This subsection conducts simulation experiments to evaluate the finite-sample performance of the low-dimensional composite quantile regression (CQR) estimation at (2.2).

The Tukey Lambda distribution in Example 2 is used to generate the $i.i.d.$ sample,

where $\{U_{i}\}$ are independent and follow $\text{Uniform}(0,1)$, $\mathbf{X}_{i}=(1,X_{i1},X_{i2})^{\prime}$, $\{(X_{i1},X_{i2})^{\prime}\}$ is an $i.i.d.$ sequence with 2-dimensional standard normality. The true parameter vector is $\boldsymbol{\beta}_{0}=(\boldsymbol{\beta}_{01}^{\prime},\boldsymbol{\beta}_{02}^{\prime},\boldsymbol{\beta}_{03}^{\prime})^{\prime}$, and we set the location parameters $\boldsymbol{\beta}_{01}=(1,0.5,-1)^{\prime}$, the scale parameters $\boldsymbol{\beta}_{02}=(1,0.5,-1)^{\prime}$ and the tail parameters $\boldsymbol{\beta}_{03}=(1,-1,1)^{\prime}$. For the tail index $\theta_{3}(\mathbf{X}_{i},\boldsymbol{\beta})$, before generating the data, we first scale each covariate into the range of $[-0.5,0.5]$ such that a relatively stable sample can be generated. In addition, $g_{1}$, $g_{2}$ and $g_{3}$ are the inverse of link functions. We choose identity link for the location index and softplus-related link for the scale and tail indices, i.e., $g_{1}(x)=x,g_{2}(x)=\textrm{softplus}(x)$ and $g_{3}(x)=1-\textrm{softplus}(x)$, where $\textrm{softplus}(x)=\log(1+\exp(x))$ is a smoothed version of $x_{+}=\max\{0,x\}$ and hence the name. Note that $g_{2}(x)>0$ and $g_{3}(x)<1$. We consider three sample sizes of $n=500$, 1000 and 2000, and there are 500 replications for each sample size.

The algorithm for CQR estimation in Section 3 is applied with $K=10$ and $\tau_{k}$’s being equally spaced over $[\tau_{L},\tau_{U}]$. We consider three quantile ranges of $(\tau_{L},\tau_{U})=(0.5,0.99)$, $(0.7,0.99)$ and $(0.9,0.99)$, and the estimation efficiency is first evaluated. Figure 1 gives the boxplots of three fitted location parameters $\widehat{\boldsymbol{\beta}}_{1n}=(\widehat{\beta}_{1,1},\widehat{\beta}_{1,2},\widehat{\beta}_{1,3})^{\prime}$. It can be seen that both bias and standard deviation decrease as the sample size increase. Moreover, when $\tau_{L}$ decreases, the quantile levels with richer observations will be used for the estimation and, as expected, both bias and standard deviation will decrease. Boxplots for fitted scale and tail parameters show a similar pattern and hence are omitted to save the space.

We next evaluate the prediction performance of $Q(\tau^{*},\boldsymbol{\theta}(\mathbf{X},\widehat{\boldsymbol{\beta}}_{n}))$ at two interesting quantile levels of $\tau^{*}=0.991$ and 0.995. Consider two values of covariates, $\boldsymbol{X}=(1,0.1,-0.2)^{\prime}$ and $(1,0,0)^{\prime}$, and the corresponding tail indices are $\theta_{3}(\mathbf{X},\boldsymbol{\beta}_{0})=-0.1139$ and $-0.3132$, respectively. Note that the Tukey lambda distribution can provide a good approximation to Cauchy and normal distributions when the tail indices are $-1$ and $0.14$, respectively, and it becomes more heavy-tailed when the tail index decreases (Freimer et al.,, 1988). The prediction error in terms of squared loss (PES), $[Q(\tau^{*},\boldsymbol{\theta}(\mathbf{X},\widehat{\boldsymbol{\beta}}_{n}))-Q(\tau^{*},\boldsymbol{\theta}(\mathbf{X},{\boldsymbol{\beta}}_{0}))]^{2}$, is calculated for each replication, and the corresponding sample mean refers to the commonly used mean square error. Table 1 presents both the sample mean and standard deviation of PESs across 500 replications as in Wang et al., 2012a . A clear trend of improvement can be observed as the sample size becomes larger, and the prediction is more accurate at the 99.1-th quantile level for almost all cases.

This subsection conducts simulation experiments to evaluate the finite-sample performance of the high-dimensional regularized estimation at (2.3).

For the data generating process at (4.1), we consider three dimensions of $p=50$, $100$ and $150$, and the true parameter vectors are extended from those in Section 4.1 by adding zeros, i.e. $\boldsymbol{\beta}_{01}=(1,0.5,-1,0,...,0)^{\prime}$, $\boldsymbol{\beta}_{02}=(1,0.5,-1,0,...,0)^{\prime}$ and $\boldsymbol{\beta}_{03}=(1,-1,1,0,...,0)^{\prime}$, which are vectors of length $p$ with $3$ non-zero entries. As a result, all true parameters $\boldsymbol{\beta}_{0}=(\boldsymbol{\beta}_{01}^{\prime},\boldsymbol{\beta}_{02}^{\prime},\boldsymbol{\beta}_{03}^{\prime})^{\prime}$ make a vector of length $3p$ with $s=9$ non-zero entries. The sample size is chosen such that $n=\left\lfloor cs\log p\right\rfloor$ with $c=5$, 10, 20, 30, 40 and 50, where $\lfloor x\rfloor$ refers to the largest integer smaller than or equal to $x$. All other settings are the same as in the previous subsection.

The algorithm for regularized estimation in Section 3 is used to search for the estimators, and we generate an independent validation set of size $5n$ to select tuning parameter $\lambda$ by minimizing the composite quantile check loss; see also Wang et al., 2012b . Figure 2 gives the estimation errors of $\|\widetilde{\boldsymbol{\beta}}_{n}-\boldsymbol{\beta}_{0}\|$. It can be seen that $\|\widetilde{\boldsymbol{\beta}}_{n}-\boldsymbol{\beta}_{0}\|$ is roughly proportional to the quantity of $\sqrt{s\log p/n}$, and this confirms the convergence rate in Theorem 4. Moreover, the estimation errors approach zero as the sample size $n$ increase, and we can then conclude the consistency of $\widetilde{\boldsymbol{\beta}}_{n}$. Finally, when $\tau_{L}$ increases, the quantile levels with less observations will be used in the estimating procedure, and hence larger estimation errors can be observed.

We next evaluate the prediction performance at quantile levels $\tau^{*}=0.991$ and $0.995$, and covariates $\mathbf{X}$ take values of $(1,0.1,-0.2,0,\cdots,0)^{\prime}$ and $(1,0,0,0,\cdots,0)^{\prime}$, similar to those in the previous subsection. Table 2 gives mean square errors of the predicted conditional quantiles $Q(\tau^{*},\boldsymbol{\theta}(\mathbf{X},\widetilde{\boldsymbol{\beta}}_{n}))$, as well as the sample standard deviations of prediction errors in squared loss, with $p=50$. It can be seen that larger sample size leads to much smaller mean square errors. Moreover, when $\tau_{L}$ is larger, the prediction also becomes worse, and it may be due to the lower estimation efficiency. Finally, similar to the experiments in the previous subsection, the prediction at $\tau^{*}=0.991$ is more accurate for almost all cases. The results for the cases with $p=100$ and 150 are similar and hence omitted.

Finally, we consider the following criteria to evaluate the performance of variable selection: average number of selected active coefficients (size), percentage of active and inactive coefficients both correctly selected simultaneously ($\textrm{P}_{\textrm{AI}}$), percentage of active coefficients correctly selected ($\textrm{P}_{\textrm{A}}$), percentage of inactive coefficients correctly selected ($\textrm{P}_{\textrm{I}}$), false positive rate (FP), and false negative rate (FN). Table 3 reports the selecting results with $p=50$ and $c=10$, 30 and 50. When $\tau_{L}$ is larger, both $\textrm{P}_{\textrm{AI}}$ and $\textrm{P}_{\textrm{A}}$ decrease, and it indicates the increasing of selection accuracy. In addition, performance improves when sample size gets larger. The results for $p=100$ and 150 are similar and hence omitted.

This subsection gives technical proofs of Lemma 1 and Theorems 1-3 in section 2.2. Two auxiliary lemmas are also presented at the end of this subsection: Lemma A.1 is used for the proof of Lemma 1, and Lemma A.2 is for that of Theorem 3. The proofs of these two auxiliary lemmas are given in a supplementary file.