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.
123,562 characters · 27 sections · 109 citation commands
Conditional Quantile Processes based on Series or Many Regressors
Quantile regression (QR) is a principal method for analyzing the impact of covariates on outcomes, particularly when the impact may be heterogeneous. This impact is characterized by the quantile function of the conditional distribution of the outcome given covariates and its functionals (Arias, Hallock and Sosa-Escudero AHS-E2001, Buchinsky Buchinsky1994 and Koenker K2005). For example, we can model the log of the individual demand for some good, $Y$, as a function of the price of the good, the income of the individual, and other observed individual characteristics, $X$, and an unobserved preference for consuming the good, $U$, as $$ Y = Q(U,X), $$ where the function $Q$ is strictly increasing in the unobservable $U$. With the normalization that $U\sim {\text{Uniform}}(0,1)$ and the assumption that $U$ and $X$ are independent, the function $Q(u,x)$ is the $u$-th quantile of the conditional distribution of $Y$ given $X = x$, i.e. $Q(u,x) = Q_{Y | X}(u | x)$. This function can be used for policy analysis. For example, we can determine how changes in taxes for the good could impact demand heterogeneously across individuals.
In this paper we develop the nonparametric QR-series framework for performing inference on the entire conditional quantile function $Q(u,x)$ and its linear functionals. In this framework, we approximate $Q(u,x)$ by a linear combination of series terms, $Z(x)'\beta(u)$. The vector $Z(x)$ includes transformations of $x$ that have good approximation properties such as powers, trigonometrics, local polynomials, splines, and/or wavelets. The function $u \mapsto \beta(u)$ contains quantile-specific coefficients that can be estimated from the data using the QR estimator of Koenker and Bassett KB78. As the number of series terms grows, the approximation error $Q(u,x) - Z(x)'\beta(u)$ decreases, approaching zero in the limit. By controlling the growth of the number of terms, we can obtain consistent estimators and perform inference on the entire conditional quantile function and its linear functionals. The QR-series framework also covers as a special case the so called many regressors model, which is motivated by many new types of data that emerge in the new information age, such as scanner and online shopping data.
We describe now the main results in more detail. Let $u \mapsto \widehat\beta(u)$ denote the QR estimator of $u \mapsto \beta(u)$. The first set of results provides large-sample theory for the normalized QR-series coefficient process of increasing dimension $u \mapsto \sqrt{n}(\widehat \beta(u) - \beta(u))$ that can be used to perform inference on the function $u \mapsto \beta(u)$. We note that inference on the function $u \mapsto \beta(u)$, in particular simultaneous inference on the parameters $\beta(u)$ that holds uniformly over all $u\in\mathcal U$, where $\mathcal U$ is a set of quantile indices of interest, is difficult because the standard asymptotic theory (van der Vaart and Wellner vdV-W) based on limit distributions does not help here as the process $u \mapsto \sqrt n(\widehat\beta(u) - \beta(u))$ in general does not have a limit distribution, even after an appropriate normalization. Instead, we develop high-quality inferential procedures based on the idea of coupling. The coupling is a construction of two processes on the same probability space that are uniformly close to each other with high probability. Typically, one of the processes is the process of interest and the other one is a process whose distribution is known up-to a relatively small number of parameters that can be consistently estimated from the data. Thus, being able to construct an appropriate coupling means that we are able to approximate the distribution of the process of interest by simulating the distribution of the coupling process from the data.
In this paper, we develop two couplings: pivotal and Gaussian, that is, for each sample size $n$, we construct a pivotal process and a Gaussian process on the same probability space as the data that are uniformly close to the process $u \mapsto \sqrt n(\widehat\beta(u) - \beta(u))$ with high probability. In other words, these pivotal and Gaussian processes strongly approximate the process $u \mapsto \sqrt n(\widehat\beta(u) - \beta(u))$. In addition, we develop four resampling methods (pivotal, gradient bootstrap, Gaussian, and weighted bootstrap) that allow us to approximately simulate the distribution of the pivotal (first two methods) and the Gaussian (last two methods) processes. These results provide an inference theory for the function $u \mapsto \beta(u)$ and also help to develop a unified feasible inference theory for linear functionals of the conditional quantile functions $x \mapsto Q(u,x)$, $u\in\mathcal U$, using any of the four proposed resampling methods. To the best of our knowledge, all of the above results are new.
The existence of the pivotal coupling emerges from the special nature of QR, where a (sub) gradient of the sample objective function evaluated at the true values of $\beta(u)$ is pivotal conditional on the regressors, up-to an approximation error. This coupling allows us to perform high-quality inference based on pivotal and gradient bootstrap methods without even resorting to Gaussian approximations. We also show that the gradient bootstrap method, originally introduced by Parzen, Wei and Ying ParzenWeiYing1994 in the parametric context, is effectively a means of carrying out the conditionally pivotal approximation without explicitly estimating Jacobian matrices, which may be difficult in the quantile regression context. The conditions for validity of the pivotal and gradient bootstrap methods require only a mild restriction on the growth of the number of series terms in relation to the sample size. To obtain the Gaussian coupling, we use chaining arguments and Yurinskii's construction. This coupling implies that one can use the Gaussian method to perform inference on the function $u \mapsto \beta(u)$, and we also use this coupling to show that the weighted bootstrap method works to approximate the distribution of the whole QR-series coefficient process for the same reason as the Gaussian method. The conditions for validity of the Gaussian and weighted bootstrap methods, however, may be stronger than those for the pivotal and gradient bootstrap methods.
As a corollary of our results on the process $u \mapsto \sqrt n(\widehat\beta(u) - \beta(u))$, we also demonstrate that the QR-series estimator $x\mapsto Z(x)'\widehat\beta(u)$ of the function $x\mapsto Q(u,x)$ based on either polynomials or splines has the fastest possible rate of convergence in the $L^2$ norm and that the QR-series estimator based on splines has the fastest possible rate of convergence in the sup norm in H\"{o}lder smoothness classes. These findings on optimality of series estimators in the quantile regression setting complement results in the literature on optimality of series estimators in the mean regression setting; see Newey W97, Huang H98, Cattaneo and Farrell CF13, Belloni et al BelloniChenChernozhukov2009, and Chen and Christensen CC13. In particular, our result on optimality in the sup norm is a major extension of Huang's work H03 on mean regression, as it requires us to establish some fine properties of the QR-series approximation; see the next section for details.
The second set of results provides estimation and inference methods for linear functionals of the conditional quantile functions, including
where $\mu$ is a given measure and $x_k$ is the $k$-th component of $x$. Specifically, we derive the pointwise rate of convergence and asymptotic normality of the QR-series estimators of the linear functionals. In addition, using our results on the QR-series coefficient process $u \mapsto \sqrt n(\widehat\beta(u) - \beta(u))$, we derive uniform rate of convergence and large sample inference procedures based on the pivotal, gradient bootstrap, Gaussian, and weighted bootstrap methods for the QR-series estimators of the linear functionals. These results provide solutions to a wide range of inference problems. For illustration purposes, we demonstrate how to use these results to construct uniform confidence bands for function-valued linear functionals and how to test shape constraints on the conditional quantile function $x\mapsto Q(u,x)$. It is noteworthy that all of the above results apply to function-valued parameters, holding uniformly in both the quantile index $u$ and the covariate value $x$. We also emphasize that although we do not treat non-linear/non-smooth functionals in this paper, our results on couplings and resampling methods for the process $u \mapsto \sqrt n(\widehat\beta(u) - \beta(u))$ are useful for treatment of such functionals.
The paper contributes and builds on the existing important literature on conditional quantile estimation. First and foremost, we build on the work of He and Shao HS00 that studied the many regressors model and gave pointwise limit theorems for the QR estimator in the case where only one quantile index $u$ is of interest. We go beyond the many regressors model to the series model and develop large sample estimation and inference results for the entire QR process. We also develop analogous estimation and inference results for the conditional quantile function and its linear functionals, such as derivatives, average derivatives, conditional average derivatives, and others. None of these results were available in the previous work. We also build on Lee Lee2003 that studied QR estimation of partially linear models in the series framework for a single quantile index $u$, and on Horowitz and Lee HL2005 that studied nonparametric QR estimation of additive quantile models for a single quantile index $u$ in a series framework. Our framework covers these partially linear models and additive models as important special cases, and allows us to perform inference on a considerably richer set of functionals, uniformly across covariate values and a continuum of quantile indices. After the first version of this paper appeared in BCCF11, Chao, Volgushev and Cheng Chao-17 developed related results for QR-series coefficient processes based on weak approximations.\footnote{We refer the reader to Chao-17 for a more detailed comparison with our results.} In a more applied side, Koenker and Schorfheide KS1994 used a smoothing splines QR method to estimate the conditional quantile functions of global temperature over the last century. Other important work includes Stute S86, Chaudhuri Chaudhuri-1991, Chaudhuri, Doksum and Samarov Chaudhuri-Doksum-Samarov-1997, H\"{a}rdle, Ritov, and Song hardle-ritov-song2009, Horderlein and Mammen HoderleinMammen2009, Cattaneo, Crump, and Jansson cattaneo-crump-jansson2010, Kong, Linton, and Xia kong-linton-xia-2010, Qu and Yoon QuYoon2011, and Guerre and Sabbah GS12, among others, but these papers focused on local, non-series, methods.
Our work also relies on the series literature, at least in a motivational and conceptual sense. In particular, we rely on the work of Stone Stone1982, Andrews Andews1991, Newey W97, Chen and Shen Chen-Shen1998, Chen Chen2006 and others that rigorously motivated the series framework as an approximation scheme and gave pointwise normality results for least squares estimators, and on Chen Chen2006 and van de Geer Geer2002 that gave (non-uniform) consistency and rate results for general series estimators, including quantile regression for the case of a single quantile index $u$. White White1992 established non-uniform consistency of nonparametric estimators of the conditional quantile function based on a nonlinear series approximation using artificial neural networks. In contrast to the previous results, our rate results are uniform in covariate values and quantile indices, and cover both the quantile function and its functionals. Moreover, we not only provide estimation rate results, but also derive a full set of results on feasible inference based on the couplings for the process $u \mapsto \sqrt n(\widehat\beta(u) - \beta(u))$.
While relying on previous work for motivation, our results require to develop both new proof techniques and new approaches to inference. In particular, our proof techniques rely on new maximal inequalities for function classes with growing moments and uniform entropy. In addition, as explained above, and in contrast to previous papers, our inference results heavily rely on the idea of the couplings. Yurinskii's coupling was previously used by Chernozhukov, Lee and Rosen CLR2009 to obtain a Gaussian coupling for the least squares series estimator, but the use of this technique in our context is new and much more involved. Thus, we approximate an entire QR-series coefficient process of an increasing dimension, instead of a vector of increasing dimension, by a Gaussian process. Finally, it is noteworthy that our uniform inference results on functionals, where uniformity is over covariate values, have not had analogs even in the least squares series literature until recently (the extension of our results to least squares has been recently published in Belloni et al BelloniChenChernozhukov2009).
The results developed in this paper for series (global) estimation are of interest even though some of the results have analogs in the literature on kernel (local) estimation, because series estimators have several attractive features that are not shared by kernel estimators. First, series methods represent the estimate of the whole conditional quantile function $x\mapsto Q(u,x)$ and its linear functionals via a relatively small set of parameters, that is, the estimates $\widehat \beta(u)$ of the QR-series coefficients $\beta(u)$. As long as these parameters are reported, any researcher will be able to calculate the value of the estimate, for example, of $Q(u,x)$ for any $x\in\mathcal X$, which is very convenient. Second, series methods allow to easily impose shape constraints. For example, if it is known that the function $x\mapsto Q(u,x)$ is concave, one can simply impose this constraint on the QR-series optimization problem to obtain concave estimates (imposing shape constraints is especially convenient when B-splines are used; see DeVore and Lorentz DL93). Third, as we demonstrate in the empirical example in Section (ref), when the vector $X$ contains many controls in addition to one main covariate of interest, series methods are particularly convenient to impose a partial linear form on the functions $x\mapsto Q(u,x)$, which helps curb the curse of dimensionality.
This paper does not deal with sparse models, where there are some key series terms and many “non-key” series terms which ideally should be omitted from estimation. In these settings, the goal is to find and indeed remove most of the “non-key” series terms before proceeding with estimation. Belloni and Chernozhukov BC-SparseQR obtained rate results for quantile regression estimators in this case, but did not provide inference results. Even though our paper does not explicitly deal with inference in sparse models after model selection, the methods and bounds provided herein are useful for analyzing this problem. Investigating this matter rigorously is a challenging issue, since it needs to take into account the model selection mistakes in estimation, and is beyond the scope of the present paper; however, it is a subject of our ongoing research, see Belloni, Chernozhukov and Kato belloni2013valid.
{\bf Plan of the paper.} The rest of the paper is organized as follows. In Section (ref), we describe the nonparametric QR-series model and estimators. In Section (ref), we derive asymptotic theory for the QR-series processes. In Section (ref), we give estimation and inference theory for linear functionals of the conditional quantile function. In Section (ref), we present an empirical application to the demand of gasoline and a computational experiment calibrated to the application. The computational algorithms to implement our inference methods are collected in the Appendix. The Supplemental Material BCCF16 contains some further results, and the proofs of the main results.
{\bf Notation.} In what follows, for all $x = (x_1,\dots,x_m)' \in\mathbb R^m$, we use $\|x\|$ to denote the Euclidean norm of $x$, that is, $\|x\| = (x_1^2 + \dots + x_m^2)^{1/2}$, and we use $\|x\|_{\infty}$ to denote the sup norm of $x$, that is, $\|x\|_{\infty} = \max_{1\leq j\leq m}|x_j|$. Also, we use $S^{m-1}$ to denote the unit sphere in $ {\Bbb{R}}^m$, that is, $S^{m-1} = \{x\in \mathbb R^m\colon \|x\| = 1\}$, and for $r>0$, we use $B_m(0,r)$ to denote the ball in $\mathbb R^m$ with center at $0$ and radius $r$, that is, $B_m(0,r) = \{x\in\mathbb R^m\colon \|x\|\leq r\}$. For all $m\times m$-dimensional matrices, we use $\|A\|$ to denote the operator norm of $A$ (also known as the spectral norm), that is, $\|A\| = \sup_{\alpha\in S^{m-1}}\|A\alpha\|$. For a set $I$, ${\rm diam}(I)=\sup_{v,\bar v\in I}\|v-\bar v\|$ denotes the diameter of $I$, and $\text{int}(I)$ denotes the interior of $I$. For any two real numbers $a$ and $b$, $a\vee b = \max\{a,b\}$ and $a\wedge b = \min\{a,b\}$. The relation $a_n \lesssim b_n $ means that the inequality $a_n \leq C b_n$ holds for all $n$ with a constant $C$ that is independent of $n$. We denote by $P^*$ the probability measure induced by conditioning on realization of the data $\mathcal{D}_n = (Y_i,X_i)_{i=1}^n$. We say that a random variable $\Delta_n= o_{P^*}(1)$ in $P$-probability if for any $\epsilon>0$, we have $P^*(|\Delta_n|>\epsilon) = o_P(1)$. We typically shall omit the qualifier “in $P$-probability." The operator $E$ denotes the expectation with respect to the probability measure $P$, $\mathbb{E}_n$ denotes the expectation with respect to the empirical measure, and $\mathbb{G}_n$ denotes $\sqrt{n}(\mathbb{E}_n - E)$. Finally, we use $\varepsilon'$ to denote a constant that depends only on the constant $\varepsilon$ but is such that its value can change at each appearance.
We consider the sequence of models indexed by the sample size $n$:
where $Y_{i,n}$ is a response variable, $X_{i,n}$ is a $d_n$-dimensional vector of covariates (elementary regressors) with the support $\mathcal{X}_n\subset {\Bbb{R}}^{d_n}$, $U_{i,n}$ is an unobservable individual ranking that is distributed uniformly on $(0,1)$ and is independent of $X_{i,n}$, and $Q_n\colon [0,1]\times\mathcal{X}_n\to {\Bbb{R}}$ is an unknown function such that for all $x\in\mathcal{X}_n$, the function $u\mapsto Q_n(u,x)$ is strictly increasing. Since $U_{i,n} | X_{i,n} \sim {\text{Uniform}} (0,1)$, for each $u\in(0,1)$, $x \mapsto Q_n(u,x)$ is the $u$th quantile function of $Y_{i,n}$ conditional on $X_{i,n}$, which we refer to as the conditional $u$-quantile function. For a compact set of quantile indices $\mathcal U\subset(0,1)$, we are interested in estimating the functions $x\mapsto Q_n(u,x)$ and their functionals. For given $n$, we assume that $(X_{i,n},U_{i,n},Y_{i,n})_{i=1}^n$ is a random sample from the distribution of the triple $(X^n,U^n,Y^n)$, where $n$ denotes the sample size.
The model (ref) covers two specifications of major interest:
Both models are of interest in econometrics. The NP model is important because it is very flexible as it does not impose any parametric structure on the functions $x\mapsto Q(u,x)$ and also does not require the function $(u,x)\mapsto Q(u,x)$ to be separately additive in $u$; see Matzkin M03 for extensive evidence on importance of this flexibility in economics. The MR model is also important, and different versions of this model have recently attracted much attention in the literature due to emergence of datasets with information on many variables; see Cattaneo, Jansson, and Newey CJN15 for some recent advances and also Mammen M93 for some classical results on the mean regression version of this model. As we demonstrate in this paper, both models can be treated in a unifying quantile regression series framework.
For brevity of notation, we shall omit the index $n$ whenever it does not lead to confusion, that is, we write $Y$, $X$, $U$, $Q$, $d$, and $\mathcal X$ instead of $Y^n$, $X^n$, $U^n$, $Q_n$, $d_n$, and $\mathcal{X}_n$, respectively, even though we implicitly assume that all these quantities are allowed to depend on $n$. Also, we write $(X_i,U_i,Y_i)_{i=1}^n$ instead of $(X_{i,n},U_{i,n},Y_{i,n})_{i=1}^n$.
Next, we introduce the QR-series approximation to the function $x\mapsto Q(u,x)$. We start with preparing some notation. Fix $u\in\mathcal U$ and let $x\mapsto Z(x) = (Z_1(x),\dots,Z_m(x))'$ be a vector of series approximating functions of dimension $m=m_n$, where each function $x\mapsto Z_j(x)$ maps $\mathcal{X}$ into $ {\Bbb{R}}$. Define the vector of coefficients $\beta(u) = (\beta_1(u),\dots,\beta_m(u))'$ as a solution to the QR-series approximation problem:
where $\rho_{u} (z) = (u - 1\{z<0\})z$ is the check function (Koenker K2005).\footnote{The optimization problem (ref) has a finite solution if $E[|Q(u,X)|]$ is finite. In addition, since the function $z\mapsto \rho_u(z)$ is strictly convex, the solution is unique if the matrix $E[Z(X) Z(X)']$ is non-singular, which is assumed in Condition S below. The term $\rho_{u} (Y - Q(u,X))$ does not affect the optimization problem but guarantees the existence of the solution when $E[|Y|]$ is not finite.} For the MR model, we assume that $Z(x) = x$ for all $x\in\mathcal{X}$, so that $m=d$ and the vector $\beta(u)$ defined in (ref) coincides with the vector $\beta_n(u)$ in the definition of the model, $Q_n(u,x) = x'\beta_n(u)$. For the NP model, we assume that the vector $Z$ consists of series functions with good approximation properties such as indicators, B-splines (or regression splines), polynomials, Fourier series, and/or compactly supported wavelets.\footnote{Interestingly, in the case of B-splines and compactly supported wavelets, the entire collection of series terms is dependent upon the sample size $n$.}$^,$\footnote{It is possible to combine these sets of approximating functions. For example, when we model gasoline consumption, we can simultaneously use Fourier series to capture seasonal effects and polynomials to capture long term growth.} We refer the reader to Newey W97 and Chen Chen2006 for a careful and detailed description of these series functions; see also Belloni et al BelloniChenChernozhukov2009 for an overview of recent advances on series approximating functions.
We define the QR-series approximating function $x\mapsto Z(x)'\beta(u)$ mapping $\mathcal{X}$ into $ {\Bbb{R}}$, and, for all $x\in\mathcal{X}$, the QR-series approximation error $$ R(u,x) := Q(u,x) - Z(x)'\beta(u). $$ We will assume that the QR-series approximation error asymptotically vanishes, that is, $\sup_{x\in\mathcal{X}, u \in \mathcal U}|R(u,x)|\to 0$ as $n\to\infty$. For the MR model, this assumption always holds because $R(u,x) = 0$ for all $x\in\mathcal{X}$. For the NP model, we will demonstrate that this assumption holds under appropriate conditions as long as $m = m_n \to \infty$ as $n\to\infty$. In turn, given that the QR-series approximation error asymptotically vanishes, it follows that the QR-series approximating function $x\mapsto Z(x)'\beta(u)$ approximates well the true conditional $u$-quantile function $x\mapsto Q(u,x)$.
The QR-series approximation motivates the {\em QR-series estimator} of the function $x\mapsto Q(u,x)$:
where $\widehat\beta(u)$ is the Koenker and Bassett KB78 estimator of $\beta(u)$ that solves the empirical analog of the population problem ((ref)):
where we denote $Z_i = Z(X_i)$ for all $i=1,\dots,n$. As $n$ gets large, both the estimation error $\widehat Q(u,x) - Z(x)' \beta(u)$ and the approximation error $R(u,x)$ asymptotically vanish.
Since we are interested in estimating the functions $x\mapsto Q(u,x)$ for a set of quantile indices $\mathcal U$, we solve the problem (ref) for all $u\in\mathcal U$ to obtain the QR-series coefficient process $$ \widehat \beta (\cdot) = \{ \widehat \beta(u) \colon u \in \mathcal{U}\} $$ and the QR-series estimator (ref) for all $u\in\mathcal U$. We note that obtaining this estimator is computationally easy even if $\mathcal U$ contains many quantile indices and the dimension $m$ of the vectors $Z_i$ is large. In particular, one can use the results of Portnoy and Koenker PortnoyKoenker97, who developed interior points methods with preprocessing for the problem (ref) that are very efficient and give the solution for multiple quantile indices simultaneously.
Let $\kappa\in(0,\infty]$ be some constant that is independent of $n$. Also, for $x\in\mathcal{X}$, let $\mathcal Y_x$ denote the support of the conditional distribution of $Y$ given $X = x$. Moreover, let $\bar {\mathcal U}$ denote the convex hull of $\mathcal U$. Throughout the paper, we will use the following regularity condition:
Condition S.1 requires that the data is i.i.d. but it can be extended to standard time series models at the expense of more technicalities. Condition S.2 imposes mild smoothness assumptions on the conditional density function $f_{Y|X}(y|x)$. Since it follows from simple algebra, see for example (ref), that $$ \frac{1}{f_{Y|X}(y|x)} = \frac{\partial Q(Q^{-1}(y,x),x)}{\partial u},\quad y\in\mathcal Y_x, \ x\in\mathcal X, $$ where $y\mapsto Q^{-1}(y,x)$ denotes the inverse of $u\mapsto Q(u,x)$, it is easy to provide a set of conditions in terms of the function $Q(u,x)$ that imply Condition S.2. Indeed, Condition S.2 follows if (i) $\partial Q(u,x)/\partial u$ is bounded away from zero uniformly over $u\in[0,1]$, $x\in\mathcal X$, and $n$; (ii) $\partial Q(u,x)/\partial u$ is bounded from above uniformly over $u\in \bar{\mathcal U}$, $x\in\mathcal{X}$ and $n$; (iii) $\partial^2 Q(u,x)/\partial u^2$ is bounded in absolute value from above uniformly over $u\in [0,1]$, $x\in\mathcal{X}$, and $n$.\footnote{Note that we assume that the conditional density $f_{Y|X}(y|x)$ is bounded away from zero only for $y = Q(u,x)$, where $u\in\bar{\mathcal U}$ and $x\in\mathcal X$. This allows us to avoid the stronger condition that assumes that $f_{Y|X}(y|x)$ is bounded away from zero for all $y\in\mathcal Y_x$ and $x\in\mathcal X$. The latter condition can simplify some arguments (see the proof of Lemma (ref)) but it requires the conditional density of $Y$ given $X$ to have bounded support, thus excluding some important distributions such as the Gaussian.}
For the MR model, Condition S.3 implies that there is no perfect multicollinearity among covariates, and Condition S.4 is satisfied with $\kappa = \infty$ since $R(u,x) = 0$ for all $u\in\mathcal U$ and $x\in\mathcal{X}$.
For the MR model, all conditions can be regarded as primitive. For the NP model, Conditions S.1 and S.2 are also primitive but Conditions S.3 and S.4 depend on the vector of approximating series functions $x \mapsto Z(x)$ used for the estimation. Therefore, below we provide some discussion of these conditions in the NP model. Suppose that $X$ is absolutely continuous with respect to the Lebesgue measure on $\mathcal X$ and let $f_X\colon \mathcal{X}\to {\Bbb{R}}$ denote its pdf. Then it is well-known that Condition S.3 holds if $f_X(x)$ is bounded from above and away from zero uniformly over $x\in\mathcal{X}$, and the eigenvalues of the matrix $$ \int_{x\in\mathcal{X}}Z(x)Z(x)'d x $$ are bounded from above and away from zero uniformly over $n$; see, for example, Proposition 2.1 in Belloni et al BelloniChenChernozhukov2009. In turn, the latter condition holds if, for example, the vector $Z$ consists of functions that are orthonormal on $\mathcal{X}$. In the case that the former condition is violated in the sense that the density $f_X(x)$ is not bounded away from zero uniformly over all $x\in \mathcal{X}$, one can consider a subset $\widetilde \mathcal{X}$ of $\mathcal{X}$ such that $f_X(x)$ is bounded away from zero uniformly over $x\in\widetilde \mathcal{X}$ and consider the estimation problem based on the subset of observations $i$ satisfying $X_i\in\widetilde \mathcal{X}$. This will give the estimate of $Q(u,x)$ for all $x\in\widetilde \mathcal{X}$ and $u\in\mathcal U$. As the sample size gets larger, one can increase the set $\widetilde \mathcal{X}$ to extend the estimate of $Q(u,x)$ to a larger set of points. Developing a method how this truncation should be performed in practice, however, is beyond the scope of this paper.
To provide some primitive conditions for Condition S.4 in the NP model, we need to prepare some notation. For a $d$-tuple $\alpha = (\alpha_1,\dots,\alpha_d)$ of nonnegative integers, let $D^\alpha = \partial^{\alpha_1}_{x_1}\cdots \partial^{\alpha_d}_{x_d}$. Also, for $s>0$, let $[s]$ denote the largest integer strictly smaller than $s$. For the constant $C>0$, define the H\"{o}lder ball $\Omega(s,C,\mathcal{X})$ as the set of all functions $f\colon \mathcal{X}\to {\Bbb{R}}$ such that
for all $x = (x_1,\dots,x_d)'$ and $\widetilde x = (\widetilde x_1,\dots,\widetilde x_d)'$ in $\mathcal X$ and all $d$-tuples $\alpha = (\alpha_1,\dots,\alpha_d)$ and $\beta = (\beta_1,\dots,\beta_d)$ of nonnegative integers satisfying $\alpha_1 +\dots+ \alpha_d = [s]$ and $\beta_1 + \dots + \beta_d \leq [s]$ (where the left-hand sides of the inequalities in (ref) are set to be infinity if the derivatives do not exist). For example, any $s$-times continuously differentiable function belongs to the H\"{o}lder ball $\Omega(s,C,\mathcal{X})$ for some $C>0$ as long as $\mathcal X$ is compact. Also, we say that the vector of approximating functions $Z$ consists of tensor products of polynomials if $m = J^d$ for some integer $J>0$ and $Z$ consists of all functions of the form $x = (x_1,\dots,x_d)\mapsto \prod_{j=1}^d x_j^{\alpha_j}$ for some $d$-tuple $\alpha = (\alpha_1,\dots,\alpha_d)$ of nonnegative integers such that $\alpha_j\leq J-1$ for all $j=1,\dots,d$. Finally, we say that the vector of approximating functions $Z$ consists of tensor products of B-splines of order $s_0$ if $m = J^d$ for some integer $J>0$ and $Z$ consists of all functions of the form $x = (x_1,\dots,x_d)\mapsto \prod_{j=1}^d b_{\alpha_j}(x_j)$ for some $d$-tuple $\alpha = (\alpha_1,\dots,\alpha_d)$ of nonnegative integers such that $\alpha_j\leq J-1$ for all $j=1,\dots,d$ where $b_0,\dots,b_{J-1}$ is a sequence of $J$ B-splines of order $s_0$ on the interval $[0,1]$ with uniform knot sequence; see Chen Chen2006 for more explanations about H\"{o}lder balls and B-splines. The next lemma provides a set of primitive conditions for Condition S.4 in the NP model.
The properties of the QR-series coefficient process and of the QR-series estimator depend on the choice of the approximating functions and the dimension of $Z$. Like in the analysis of series estimators of conditional mean functions (see Newey W97), the following quantity will play a crucial role in our analysis: $$ \zeta_m = \sup_{x\in\mathcal X}\|Z(x)\|. $$ Assuming that $\mathcal X = [0,1]^d$, it is well known that $\zeta_m \lesssim m$ if the vector $Z$ consists of tensor products of polynomials and $\zeta_m \lesssim m^{1/2}$ if the vector $Z$ consists of tensor products of B-splines.
As in the analysis of the parametric quantile regression, the following Jacobian matrix will also play a crucial role in the analysis:
Implementing some of our inference methods will require an estimator of $J(u)$. For the purposes of this paper, we will use Powell's Powell1984 estimator defined by
where $h$ is some bandwidth value satisfying $h = h_n \to 0$. We will also use the estimator of $\Sigma$ defined by
The properties of $\widehat J(u)$ and $\widehat \Sigma$ in our high-dimensional setting are established in Lemma (ref) in Appendix (ref) of the Supplemental Material.
In this section, we study properties of the normalized QR-series coefficient process $$ \sqrt n(\widehat \beta(\cdot) - \beta(\cdot)) = \Big\{\sqrt n(\widehat\beta(u) - \beta(u))\colon u\in\mathcal U\Big\}. $$ Specifically, we derive the rate of convergence and construct two couplings for this process. The couplings give two processes that are uniformly close to $\sqrt n(\widehat\beta(\cdot) - \beta(\cdot))$ with high probability, but are such that their distribution can be simulated. In particular, we develop four resampling methods (pivotal, gradient bootstrap, Gaussian, and weighted bootstrap) to simulate the distribution of these processes. In the next section, these results allow us to develop a unified feasible inference theory for all the functionals of interest. We also derive rates of convergence for the QR-series estimator process $\widehat Q(\cdot,\cdot) = \{\widehat Q(u,x)\colon u\in\mathcal U, x\in\mathcal X\}$ in the NP model. In particular, we show that the QR-series estimator based on either polynomials or B-splines has the fastest possible rate of convergence in the $L^2$ norm and that the QR-series estimator based on B-splines has the fastest possible rate of convergence in the sup norm.
As explained in the previous section, given an i.i.d. sample $(X_i,Y_i)_{i=1}^n$ from the distribution of the pair $(X,Y)$, we estimate the coefficient function $\beta(\cdot) = \{\beta(u)\colon u\in\mathcal U\}$ using the QR-series coefficient process $\widehat\beta(\cdot) = \{\widehat\beta(u)\colon u\in\mathcal U\}$, namely, for each $u \in \mathcal{U}$, we define $\widehat \beta(u)$ as the Koenker and Bassett KB78 estimator that solves the empirical analog (ref) of the population problem ((ref)). Our first main result is a uniform-in-$u$ rate of convergence for the QR-series coefficient process.
Theorem (ref) establishes a rate of convergence of the estimator $\widehat\beta(u)$ in our high-dimensional setting that holds uniformly over $u\in\mathcal U$. The theorem complements the rate of convergence results in the literature. Indeed, Koenker and Portnoy KoenkerPortnoy1987 established rate of convergence results that hold uniformly over $\mathcal{U}$ in the fixed-dimensional setting, and He and Shao HS00 established rates in the high-dimensional setting for the case when $\mathcal U$ is a singleton (pointwise-in-$u$ rate of convergence). Importantly, the uniform-in-$u$ rate of convergence in Theorem (ref) is the same as the pointwise-in-$u$ rate of convergence. The proof of this theorem relies on new concentration inequalities that control the behavior of the eigenvalues of the design matrix $\widehat\Sigma$. Note also that our condition $m\zeta_m^2\log^2 n = o(n)$ is similar to the analogous condition in HS00.
Theorem (ref) has an implication for the uniform-in-$u$ rate of convergence in the $L^2$ norm of the QR-series estimator in the NP model. Indeed, define $$ \|h\|_{L^2(X)} = \Big(E[|h(X)|^2]\Big)^{1/2},\quad \text{for }h\colon \mathcal{X}\to \mathbb R. $$ We then have the following corollary of Theorem (ref), which is the second main result together with Corollary (ref) below on the uniform-in-$u$ rate of convergence in the sup norm of the QR-series estimator in the NP model.
Here we derive two couplings yielding strong approximations to the process $\sqrt n(\widehat\beta(\cdot) - \beta(\cdot))$ in the form of either a pivotal or a Gaussian process, and develop four resampling methods to approximate the distribution of these processes and thus approximate also the distribution of the original process $\sqrt n(\widehat\beta(\cdot) - \beta(\cdot))$. We provide algorithms to implement the resampling methods in Appendix (ref).
Let
Note that the process $\mathbb U(\cdot) = \{\mathbb U(u)\colon u\in\mathcal U\}$ is (conditionally) pivotal since conditional on $(Z_i)_{i=1}^n$, the sequence $(U_i)_{i=1}^n$ consists of i.i.d. Uniform$(0,1)$ random variables. The following theorem, which is the third main result together with the Gaussian coupling in Theorem (ref) below, shows that the process $\sqrt n(\widehat\beta(\cdot) - \beta(\cdot))$ is strongly approximated by the (conditionally) pivotal process $J^{-1}(\cdot)\mathbb U(\cdot) = \{J^{-1}(u)\mathbb U(u)\colon u\in\mathcal U\}$.
This theorem is important because it has many useful implications. One of the implications is the following result for the uniform-in-$u$ rate of convergence in the sup norm of the QR-series estimator in the NP model.
We also note that although the uniform convergence rate based on polynomials is not optimal, the rate derived in Corollary (ref) is faster than the (trivial) uniform rate implied by the $L_2$ rate and the relation between the $L_2$-norm and sup-norm.
Another implication of Theorem (ref) is that it suggests the following high-quality method to approximate the distribution of the process $\sqrt n(\widehat\beta(\cdot) - \beta(\cdot))$, which we refer to as the pivotal method. First, simulate an i.i.d. sequence $(U_i^*)_{i=1}^n$ of Uniform$(0,1)$ random variables that are independent of the data and define
so that conditional on $(Z_i)_{i=1}^n$, the process $\mathbb U^*(\cdot) = \{\mathbb U^*(u)\colon u\in\mathcal U\}$ is a copy of the process $\mathbb U(\cdot)$, and $J^{-1}(\cdot)\mathbb U^*(\cdot) = \{J^{-1}(u)\mathbb U^*(u)\colon u\in\mathcal U\}$ is a copy of $J^{-1}(\cdot)\mathbb U(\cdot)$. Second, calculate the estimators $\widehat J(u)$ of the matrices $J(u)$ for all $u\in\mathcal U$ as in (ref) of Section (ref) (recall that $h$ in the estimators $\widehat J(u)$ is some bandwidth value satisfying $h = h_n \to 0$). Then, as shown in the next theorem, one can use the conditional distribution of the process $\widehat J^{-1}(\cdot)\mathbb U^*(\cdot) = \{\widehat J^{-1}(u)\mathbb U^*(u)\colon u\in\mathcal U\}$ given the data, which can be simulated, to approximate the distribution of the process $J^{-1}(\cdot)\mathbb U^*(\cdot)$, and, via Theorem (ref), also of the process $\sqrt n(\widehat\beta(\cdot) - \beta(\cdot))$.
This theorem is the fourth main result together with Theorems (ref), (ref), and (ref) below on gradient bootstrap, Gaussian, and weighted bootstrap methods. The pivotal method is closely related to another approach to inference, which we refer to as the gradient bootstrap method. This approach was previously introduced by Parzen, Wei and Ying ParzenWeiYing1994 for parametric models with fixed dimension. We extend it to the considerably more general series framework studied in this paper. The main idea is to generate for all $u\in\mathcal U$ the gradient bootstrap estimator $\widehat \beta^*(u)$ as the solution to the perturbed QR problem
where $\mathbb{U}^*(u)$ is defined in ((ref)). Then, as shown in the next theorem, one can use the conditional distribution of the process $\sqrt n(\widehat \beta^*(\cdot) - \widehat\beta(\cdot)) = \{\sqrt n(\widehat \beta^*(u) - \widehat\beta(u))\colon u\in\mathcal U\}$ given the data, which can be simulated, to approximate the distribution of the process $J^{-1}(\cdot)\mathbb U^*(\cdot)$, and, via Theorem (ref), also of the process $\sqrt n(\widehat\beta(\cdot) - \beta(\cdot))$.
Next, we turn to a strong approximation based on a sequence of Gaussian processes. The following theorem shows that for each $n$, one can construct a Gaussian process $G(\cdot) = G_n(\cdot) = \{G_n(u)\colon u\in\mathcal U\}$ such that the process $J^{-1}(\cdot)G(\cdot) = \{J^{-1}(u)G(u)\colon u\in\mathcal U\}$ is with high probability uniformly close to the process $\sqrt n(\widehat\beta(\cdot) - \beta(\cdot))$.
Although Theorem (ref) requires strong conditions, it is important because it suggests that one can approximate the distribution of the process $\sqrt n(\widehat\beta(\cdot) - \beta(\cdot))$ using Gaussian and weighted bootstrap methods, which are wide-spread in the literature in other contexts and which we now describe.
Let us start with the Gaussian method. Let $\widehat\Sigma^{1/2}$ denote the square root of the matrix $\widehat\Sigma$. Note that the covariance function of the process $G(\cdot)$ conditional on $(Z_i)_{i=1}^n$, given in (ref), is equal to that of the process $\widehat\Sigma^{1/2} B_m(\cdot)$, where $B_m(\cdot) = \{B_m(u)\colon u\in\mathcal U\}$ is a standard $m$-dimensional Brownian bridge, that is, a vector consisting of $m$ independent scalar Brownian bridges. Since the sample path of the Brownian bridge is continuous a.s., it follows that the process $\widehat\Sigma^{1/2} B_m(\cdot)$ is a copy of the process $G(\cdot)$, conditional on $(Z_i)_{i=1}^n$. Hence, one can simulate a standard $m$-dimensional Brownian bridge $B_m^*(\cdot) = \{B_m^*(u)\colon u\in\mathcal U\}$ that is independent of the data and define
so that conditional on $(Z_i)_{i=1}^n$, the process $G^*(\cdot) = \{G^*(u)\colon u\in\mathcal U\}$ is a copy of the process $G(\cdot)$, and $J^{-1}(\cdot)G^*(\cdot) = \{J^{-1}(u)G^*(u)\colon u\in\mathcal U\}$ is a copy of $J^{-1}(\cdot)G(\cdot)$. Let $\widehat J(u)$ be the estimators of the matrices $J(u)$ for all $u\in\mathcal U$ in (ref). Then, as shown in the next theorem, one can use the conditional distribution of the process $\widehat J^{-1}(\cdot)G^*(\cdot) = \{\widehat J^{-1}(u)G^*(u)\colon u\in\mathcal U\}$ given the data, which can be simulated, to approximate the distribution of the process $J^{-1}(\cdot)G^*(\cdot)$, and via Theorem (ref) also of the process $\sqrt n(\widehat \beta(\cdot) - \beta(\cdot))$.
Another related inference method is the weighted bootstrap method. Pr{\ae}stgaard and Wellner Praestgaard-Wellner-93, Hahn H97, Chamberlain and Imbens Chamberlain-Imbens-03, and Chen and Pouzo CP09 previously used this method in the point-wise case, where the set $\mathcal U$ is a singleton. We extend this method to obtain the distributional approximation for the process $\sqrt n(\widehat\beta(\cdot) - \beta(\cdot)) = \{\sqrt n(\widehat\beta(u) - \beta(u))\colon u\in\mathcal U\}$ where $\mathcal U$ is not a singleton and in fact can be a continuum of quantile indices. To describe the method, consider a set of weights $\pi_1,...,\pi_n$ that are i.i.d. draws from the distribution of a non-negative random variable $\pi$ with $E[\pi] =1$ and $E[\pi^2] =2$, such as the standard exponential distribution, and that are independent of the data. For all $u\in\mathcal U$, define the weighted bootstrap estimator $\widehat\beta^b(u)$ as the solution to the weighted QR problem
Then, as shown in the next theorem, one can use the conditional distribution of the process $\sqrt n(\widehat\beta^b(\cdot) - \widehat\beta(\cdot)) = \{\sqrt n(\widehat\beta^b(u) - \widehat\beta(u))\colon u\in\mathcal U\}$ given the data, which can be simulated, to approximate the distribution of the process $J^{-1}(\cdot)G^*(\cdot)$, and via Theorem (ref) also of the process $\sqrt n(\widehat\beta(\cdot) - \beta(\cdot))$.
In addition to the quantile functions $x \mapsto Q(u,x)$, $u\in\mathcal U$, we are also interested in various linear functionals of these functions. If $x$ is decomposed as $(w,v)$ and $x_k$ denotes the $k$-th component of $x$, examples of particularly useful linear functionals include
The measures $\mu$ entering the definitions above are assumed to be known but our results can also be extended to include estimated measures. To cover all examples, we denote the linear functional of interest by $\theta(u,w)$, where $w\in\mathcal W\subset \mathbb R^{d_w}$.
Let $I\subset \mathcal U \times \mathcal W$ denote the set of values of $(u,w)$ of interest. For example, if we are interested in
By the linearity of the series approximations, the function $\theta(u,w)$ can be seen as a linear functional of the quantile regression coefficients $\beta(u)$ up to an approximation error, that is,
where $\ell(w)' \beta(u)$ is the QR-series approximation, with $\ell(w)$ denoting the $m$-dimensional vector of loadings on the coefficients, and $r(u,w)$ is the remainder term, which corresponds to the QR-approximation error. Indeed, this decomposition arises from the application of different linear operators $\mathcal{A}$ to the decomposition $Q(u,\cdot) = Z(\cdot)'\beta(u) + R(u,\cdot)$ and evaluating the resulting functions at $w$:
In the three examples above the operator $\mathcal{A}$ is given by, respectively,
For notational convenience, we use the formulation ((ref)) in the analysis, instead of the motivational formulation ((ref)).
Given $\widehat Q(u,x) = Z(x)'\widehat\beta(u)$, we use the plug-in estimator $$ \widehat \theta(u,w) = \ell(w)' \widehat \beta(u),\quad (u,w)\in I, $$ to estimate $\theta(u,w)$. In cases where $\theta(u,w)$ is known to be monotone with respect to either $w$ or $u$, we show in Appendix (ref) of the Supplemental Material how to impose this restriction after estimation to improve finite sample properties of $\widehat \theta(u,w)$.
In the rest of this section, we provide rates of convergence for $\widehat\theta(u,w)$ as well as the inference tools that will be valid for inference on the QR-series approximation $$ \ell(w)' \beta(u), \quad (u,w) \in I, $$ and, provided that the QR-approximation error $r(u,w)$ is small enough relative to the estimation noise, will also be valid for inference on the functional of interest: $$ \theta(u,w), \ \ (u,w) \in I. $$ Thus, the QR-series approximation $\ell(w)'\beta(u)$ is an important penultimate target, whereas the functional $\theta(u,w)$ is the ultimate target.
We start with the rate of convergence of the estimator $\widehat\theta(u,w)$ at a particular quantile index value $u$ and a particular covariate value $w$ (pointwise rate of convergence). In principle, the point $(u,w)$ can depend on $n$, but we suppress the dependence for simplicity of notation. We use the following assumption:
Condition P. The QR-series decomposition $\theta(u,w) = \ell(w)' \beta(u) + r(u,w)$ satisfies $$ \frac{\sqrt{n}|r(u,w)|}{\|\ell(w)\|} = o(1). $$ Condition P can be understood as an undersmoothing condition. Although undersmoothing conditions are widely spread in the literature, as Belloni et al BelloniChenChernozhukov2009 pointed out, there is no theoretically justified procedure in the literature that would lead to a desired level of undersmoothing for the estimators of the linear functionals even for least squares estimators. For example, under conditions of Lemma (ref), when $w = x$, $\theta(u,w) = Q(u,x)$, so that $\ell(w) = Z(x)$ and $r(u,w) = R(u,x)$, and the vector $Z$ consists of a tensor product of B-splines of order $s_0$, Condition P holds as long as $n / m^{1+2(s\wedge s_0)/d} = o(1)$.
Based on Condition P, we derive the following theorem for the pointwise rate of convergence of $\widehat\theta(u,w)$, which is the fifth main result together with Theorem (ref) below on pointwise asymptotic normality of $\widehat\theta(u,w)$.
In order to perform inference, we consider the t-statistic $$ t(u,w) = \frac{ \widehat \theta(u,w) - \theta(u,w) }{ \widehat \sigma(u,w)}, $$ where
is a consistent estimator of
the asymptotic variance of $\widehat \theta(u,w)$, obtained by the delta method. We can carry out standard inference based on this t-statistic because $t(u,w) \to_d N(0,1)$, as we establish below.
Next, we derive the rate of convergence of the estimator $\widehat\theta(u,w)$ that holds uniformly over $(u,w)\in I$. We use the following assumption:
Condition U.
}
Condition U.1 on the dimension and the diameter of the set $I$ is mild and can be further relaxed at the expense of additional technicalities. As in the pointwise case, Condition U.2 can be understood as an undersmoothing condition. Condition U.3 requires that the vector of loadings $\ell(w)$ is bounded uniformly over $w\in\mathcal W$ in the Euclidean norm by $\zeta_{m,\theta}$ and the function $w\mapsto \ell(w)/\|\ell(w)\|$ is Lipschitz-continuous in the Euclidean norm with the Lipschitz constant $\zeta_{m,\theta}^L$. We note that the last condition is rather weak because the only requirement on the Lipschitz constant that we impose is that $\log\zeta_{m,\theta}^L \lesssim \log n$. We discuss some bounds on the constant $\zeta_{m,\theta}$ in a separate comment below.
The following theorem establishes the uniform rate of convergence of the QR-series estimator $\widehat \theta(u,w)$, which is the sixth main result.
The uniform rate of Theorem (ref) is the same as the pointwise rate of Theorem (ref) up to a small logarithmic factor. As in the pointwise rate result, the norm of the vector of loadings play a role which is controlled by $\zeta_{m,\theta}$ in Condition U.
Next, we consider inference on the function $(u,w)\mapsto \theta(u,w)$. We base inference on the t-statistic process $t(\cdot,\cdot) = \{t(u,w)\colon (u,w)\in I\}$ defined as follows:
where $\widehat \sigma^2(u,w)$, defined in ((ref)), is an estimator of the asymptotic variance $\sigma^2(u,w)$ of $\widehat\theta(u,w)$ in ((ref)). Using the results in the previous section, we construct pivotal and Gaussian couplings for this process in the following theorem, which is the seventh main result together with Theorems (ref), (ref), and (ref) below on couplings and resampling methods for the t-statistic process.
The Gaussian coupling is derived in this theorem under rather strong condition $m^7\zeta_m^6 = o(n^{1 - \varepsilon})$. It turns out that it is possible to construct the same coupling under a different set of conditions:
As in Section (ref), we can use four resampling methods to approximately simulate the distribution of the pivotal and Gaussian processes. Specifically, define the processes $\mathbb U^*(\cdot)$ and $G^*(\cdot)$ as in (ref) and (ref), respectively. Recall that conditional on $(Z_i)_{i=1}^n$, these processes are copies of the processes $\mathbb U(\cdot)$ and $G(\cdot)$, respectively, and so the processes $$ \left\{\frac{\ell(w)'J^{-1}(u)\mathbb U^*(u)/\sqrt n}{\sigma(u,w)}\colon (u,w)\in I\right\} \ \text{ and } \ \left\{\frac{\ell(w)'J^{-1}(u)\mathbb G^*(u)/\sqrt n}{\sigma(u,w)}\colon (u,w)\in I\right\} $$ are copies of the the pivotal and Gaussian processes $$ \left\{\frac{\ell(w)'J^{-1}(u)\mathbb U(u)/\sqrt n}{\sigma(u,w)}\colon (u,w)\in I\right\} \ \text{ and } \ \left\{\frac{\ell(w)'J^{-1}(u)\mathbb G(u)/\sqrt n}{\sigma(u,w)}\colon (u,w)\in I\right\}, $$ respectively. Also, define the t-statistic bootstrap process $t^*(\cdot,\cdot) = \{t^*(u,w)\colon (u,w)\in I\}$ for each method as
The following theorem shows that the conditional distribution of the t-statistic bootstrap process $t^*(\cdot,\cdot)$ given the data, which can be simulated, approximates the distribution of the pivotal (in the case of pivotal and gradient bootstrap methods) and Gaussian (in the case of Gaussian and weighted bootstrap methods) processes, and via Theorems (ref) and (ref) also of the original t-statistic process $t(\cdot,\cdot)$.
Note that in the case of weighted bootstrap method, the theorem above imposes the rather strong condition $m^7\zeta_m^6 = o(n^{1 - \varepsilon})$. Like in the case of Theorem (ref), it turns out that it is possible to obtain the same approximation as in this theorem but under a different set of conditions:
With the help of Theorems (ref) -- (ref), we can solve a wide range of inference problems. For example, we can construct uniform confidence bands for linear functionals $(u,w) \mapsto \theta(u,w)$ on $I$, and test shape constraints for the conditional quantile functions $x\mapsto Q(u,x)$. For the former problem, let $$ V = \sup_{(u,w)\in I} | t(u,w)| $$ be the maximal t-statistic. Also, let $k(1 - \alpha)$ denote the $(1-\alpha)$ quantile of the distribution of $V$. If $k(1 - \alpha)$ were known, we would have the confidence band
covering the whole function $\{\theta(u,w)\colon (u,w)\in I\}$ with probability $1 - \alpha$ exactly. However, $k(1-\alpha)$ is typically unknown, and the confidence band (ref) is infeasible. Instead, we approximate $k(1 - \alpha)$ using the resampling methods developed in this paper. Specifically, let $$ V^* = \sup_{(u,w) \in I} | t^*(u,w)| $$ be the bootstrap maximal t-statistic, and let $k^*(1-\alpha)$ be the $(1 - \alpha)$ quantile of the conditional distribution of $V^*$ given the data. This quantity can be computed numerically by Monte Carlo methods, as we illustrate in the next section via empirical examples and give precise algorithms in Appendix (ref). We then form a two-sided $(1-\alpha)$ uniform confidence band as $$ \Big\{[\dot{\iota} (u,w), \ddot{\iota}(u,w)] = [ \widehat \theta(u,w) - k^*(1-\alpha) \widehat \sigma(u,w), \ \widehat \theta(u,w) + k^*(1-\alpha) \widehat \sigma(u,w)]\colon (u,w) \in I\Big\}. $$
The following theorem establishes that this confidence band covers the whole function $\{\theta(u,w)\colon (u,w)\in I\}$ with probability $(1-\alpha)$ in large samples.
In addition to the validity of the uniform confidence band, Theorem (ref) establishes that the width of the uniform confidence band is of the same order as the uniform rate of convergence of the estimator $\widehat\theta(u,w)$.
We consider the problem of testing shape constraints for the conditional quantile functions $x\mapsto Q(u,x)$. Let $x_k$ denote the $k$-th component of $x$. We assume that the functions $x\mapsto Q(u,x)$ are twice continuously differentiable and consider three types of shape constraints:
where in the third example $\partial_x^2 Q(u,x)$ denotes the $d\times d$-dimensional matrix whose $(j,k)$-th element is given by $\partial_{x_j}\partial_{x_k}Q(u,x)$ for all $j,k = 1,\dots,d$.\footnote{Note that the twice continuously differentiable function $f\colon \mathcal X\to\mathbb R$ is concave if and only if $\alpha'\partial_x^2 f(x)\alpha\leq 0$ for all $x\in\mathcal X$ and $\alpha\in S^{d-1}$. To prove this claim, note that $f$ is concave if and only if the function $t\mapsto f(x + t \alpha)$ mapping $\{t\in \mathbb R\colon x+t\alpha\in \mathcal X\}$ to $\mathbb R$ is concave for all $x\in\mathcal X$ and $\alpha\in S^{d-1}$, which in turns holds if and only if $\alpha'\partial_x^2 f(x)\alpha\leq 0$ for all $x\in\mathcal X$ and $\alpha\in S^{d-1}$.} Note that in the first example we focus on the case where $x\mapsto Q(u,x)$ is decreasing with respect to $x_k$ but we can also consider the case where $x\mapsto Q(u,x)$ is increasing simply by replacing $Y$ by $-Y$ and $u$ by $1 - u$. Similarly, we can consider the case of convexity in the second and third examples.
Now, observe that all three shape constraints discussed above can be expressed using the same notation: $$ \theta(u,w)\leq 0,\quad\text{for all }(u,w)\in I, $$ where $\theta(u,w)$ is a linear functional with the vector of loadings being $\ell(w) = \partial_{x_k}Z(x)$ with $w = x$ in the first example, $\ell(w) = \partial^2_{x_k}Z(x)$ with $w = x$ in the second example, and $\ell(w) = (\ell_1(w),\dots,\ell_m(w))'$ where $\ell_j(w) = \alpha'\partial_x^2 Z_j(x)\alpha$ for all $j=1,\dots,m$ with $w = (x,\alpha)$ in the third example. Hence, we are interested in testing $$ H_0\colon \sup_{(u,w)\in I}\theta(u,w)\leq 0 \ \text{ against } \ H_1\colon \sup_{(u,w)\in I}\theta(u,w) > 0. $$ To test $H_0$ against $H_1$, we consider the one-sided Kolmogorov-Smirnov statistic $$ T = \sup_{(u,w)\in I} \frac{\widehat\theta(u,w)}{\widehat\sigma(u,w)}. $$ Then under $H_0$, $$ T = \sup_{(u,w)\in I} \frac{\widehat\theta(u,w)}{\widehat\sigma(u,w)} \leq \sup_{(u,w)\in I}\frac{\widehat\theta(u,w) - \theta(u,w)}{\widehat\sigma(u,w)} = \sup_{(u,w)\in I} t(u,w). $$ Note also that too large values of $T$ suggest that $H_0$ is violated. Hence, letting $\tilde k(1 - \alpha)$ denote the $(1 - \alpha)$ quantile of $\sup_{(u,w)\in I}t(u,w)$, we would like to reject $H_0$ if $T > \tilde k(1 - \alpha)$. However, such a test is not feasible because $\tilde k(1 - \alpha)$ is unknown. Instead, we approximate $\tilde k(1 - \alpha)$ using the resampling methods developed in this paper. Specifically, let $$ T^* = \sup_{(u,w)\in I}t^*(u,w) $$ be the bootstrap statistic, and let $\tilde k^*(1 - \alpha)$ be the $(1 - \alpha)$ quantile of the conditional distribution of $T^*$ given the data. This quantity can be computed numerically by Monte Carlo methods. Then we reject $H_0$ in favor of $H_1$ if $T > \tilde k^*(1 - \alpha)$. The following theorem shows that this test controls size in large samples.
This section illustrates the finite sample performance of the estimation and inference methods with two examples. All the calculations were carried out with the software \verb"R" (R08), using the package \verb"quantreg" for quantile regression (Koenker Koenker08). We refer to Appendix (ref) for implementation algorithms and to the companion computational paper RJ06 for software and additional examples.
To illustrate our methods with real data, we consider an empirical application to nonparametric estimation of the demand for gasoline. Blundell, Horowitz and Parey BHP2012, Hausman and Newey HN1995, Schmalensee and Stoker SS1999, and Yatchew and No YN2001 estimated nonparametrically the average demand function. We estimate nonparametrically the quantile demand and price elasticity functions and apply our inference methods to construct confidence bands for the average quantile price elasticity function and to test the Slutsky condition of consumer demand. We use the same data set as in Yatchew and No YN2001, which comes from the National Private Vehicle Use Survey, conducted by Statistics Canada between October 1994 and September 1996.\footnote{The data set can be downloaded from Adonis Yatchew's web site at www.economics.utoronto.ca/yatchew/.} The main advantage of this data set, relative to similar data sets for the U.S., is that it is based on fuel purchase diaries and contains detailed household level information on prices, fuel consumption patterns, vehicles and demographic characteristics. (See Yatchew and No YN2001 for a more detailed description of the data.) Our sample selection and variable construction also follow Yatchew and No YN2001. We select into the sample households with non-zero licensed drivers, vehicles, and distance driven. We focus on regular grade gasoline consumption. This selection results in a sample of 5,001 households. Fuel consumption and expenditure are recorded by the households at the purchase level.
We consider a partially linear specification for the demand function:\footnote{This partially linear specification of the demand function arises from household preferences characterized by the indirect utility function $V(w,v,u) = v^{1-\beta(u)}/[1-\beta(u)] - G(u,w)$, where $w$ is real gasoline price, $v$ is real income, and $g(u,w) = \partial_w G(u,w)$ (see Lewbel Lewbel1987, Th. 1). We thank Arthur Lewbel for pointing this out.} $$ Y= Q(U,X), \ \ \ Q(U,X)= g(U,W) + V'\beta(U), \ \ \ X = (W,V), $$ where $Y$ is the log of total gasoline consumption in liters per month; $W$ is the log of price in Canadian dollars per liter; $U$ is the unobservable preference of the household to consume gasoline; and $V$ is a vector of 28 covariates. Following Yatchew and No YN2001, the covariate vector includes the log of age, a dummy for the top coded value of age, the log of income, a set of dummies for household size, a dummy for urban dwellers, a dummy for young-single (age less than 36 and household size of one), the number of drivers, a dummy for more than 4 drivers, 5 province dummies, and 12 monthly dummies. To estimate the function $w \mapsto g(w, u)$ at each $u$, we consider three different vectors of series approximating functions $w\mapsto Z(w)$: linear, a power orthogonal polynomial of degree 6, and a cubic B-spline with 5 knots at the $\{0, 1/4, 1/2, 3/4, 1\}$ quantiles of the observed values of $W$. The series approximation to the function $(u,x) \mapsto Q(u,x)$ takes the following form: $$ Q(u,x) = Z(w)'\delta(u) + v'\gamma(u) = Z(x)'\beta(u), \quad Z(x) = (Z(w),v), \quad \beta(u) = (\delta(u), \gamma(u)). $$ The number of series terms in the power and B-spline specifications is selected by undersmoothing over the specifications chosen by applying cross validation to the corresponding least squares estimators.\footnote{There is potentially a large set of methods that can be used to choose the number of series terms (cross-validation, penalization, the method of Lepski, among others). Indeed, the problem of selecting the number of series terms is a special case of the problem of model selection, and there are several textbooks/monographs in the literature on model selection in abstract settings; for example, Massart M07 and Koltchinskii K11. However, to the best of our knowledge, there are no papers in the literature that apply to the problem of selecting the number of series terms in the nonparametric quantile regression problem studied here. Hence, we have opted to use an ad hoc method that consists of performing cross-validation as if we were to estimate the conditional mean function $x\mapsto E[Y | X = x]$, which is estimated by the series least squares method. Under the implicit assumption that the smoothness of the functions $x\mapsto Q(u,x)$ is similar to that of the function $x\mapsto E[Y | X = x]$, such a cross-validation would yield the number of series terms that approximately equalize variance and bias terms in estimating the functions $x\mapsto Q(u,x)$. We then slightly increase the number of series terms so that the bias term is of smaller order relative to the variance term (that is, to achieve undersmoothing, as stated in Condition U.2), so that valid inference can be performed.} In the next section, we analyze the size of the specification error of these series approximations in a numerical experiment calibrated to mimic this example.
The empirical results for the B-spline specification are reported in Figures (ref) and (ref).\footnote{The results for the linear and power specifications are not reported for the sake of brevity. They are similar to the results for the B-spline specification.} The first two panels of fig. (ref) plot the initial and monotonized estimates of the quantile demand surface for gasoline as a function of price and the quantile index, that is $$ (u,\exp(w)) \mapsto \theta(u,w) = \exp(g(w, u) + v'\beta(u)), $$ where the value of $v$ is fixed at the sample median values of the ordinal variables and one for the dummies corresponding to the sample modal values of the rest of the variables.\footnote{The median values of the ordinal covariates are $\$40K$ for income, $46$ for age, and $2$ for the number of drivers. The modal values for the rest of the covariates are $0$ for the top-coding of age, $2$ for household size, $1$ for urban dwellers, $0$ for young-single, $0$ for the dummy of more than 4 drivers, $4$ (Prairie) for province, and $11$ (November) for month.} The monotonized estimates are obtained using the average rearrangement over both the price and quantile dimensions proposed in Chernozhukov, Fern\'andez-Val, and Galichon CFG2010-Biometrika; see Appendix (ref). The demand surface show most noticeably non-monotone areas with respect to price at high quantiles, which are removed by the rearrangement. The last panel of fig. (ref) shows the estimate of the quantile price elasticity surface as a function of price and the quantile index, that is: $$ (u,\exp(w)) \mapsto \theta(u,w) = \partial_w g(u,w). $$ The estimates show substantial heterogeneity of the elasticity across quantiles and prices, with individuals at the upper quantiles being less sensitive to high prices.\footnote{These estimates are smoothed by local weighted polynomial regression across the price dimension (Cleveland C1979), because the unsmoothed elasticity estimates display very erratic behavior.}
Fig. (ref) shows 90% uniform confidence bands for the average quantile price elasticity function $$ u \mapsto \theta(u) = \int \partial_{w} \ g(u,w) d\mu(w), $$ over the quantile indices $\mathcal{I} = [0.1,0.9]$, where $\mu$ is the empirical distribution of $W$. The panels of the figure correspond to the pivotal, gradient bootstrap, Gaussian and weighted bootstrap methods. For the pivotal and Gaussian methods the distribution of the maximal t-statistic is obtained by 1,000 simulations. The gradient bootstrap uses 199 repetitions. The weighted bootstrap uses standard exponential weights and 199 repetitions. The confidence bands show that the evidence of heterogeneity in the elasticities across quantiles is not statistically significant, because we can trace a horizontal line within the bands. They also show that there is significant evidence of negative price sensitivity at most quantiles as the bands are bounded away from zero for most quantiles.
The Slutsky condition of consumer demand states that the compensated price elasticity is negative for all the households. Dette, Hoderlein, and Neumeyer DHN2011 showed that this condition has testable implications for the quantile demand function and its derivatives in heterogeneous demand systems with multiple goods and possible infinite dimensional unobservables. The one good version of their test is:
where $S(u,x)$ is the compensated quantile price elasticity that in our logarithmic specification takes the form $$ S(u,x) = \exp(\ell) \partial_w Q(u,x) + \exp(Q(u,x) + w) \partial_\ell Q(u,x), \ x = (w,\ell,c), $$ which is a smooth function of the quantile demand and derivatives. Here we have partitioned the covariate vector $X$ into $(W,L,C),$ where $W$ is log of price, $L$ is the log of income, and $C$ includes the rest of the covariates.
To test the functional hypothesis ($\ref{eq: slutsky}$) we use the one-sided Kolmogorov-Smirnov statistic: $$ K = \max_{(u,x) \in I} \frac{\widehat S(u,x)}{\widehat \sigma_S(u,x)}, $$ where $$ \widehat S(u,x) = \exp(\ell) \partial_w \widehat Q(u,x) + \exp(\widehat Q(u,x) + w) \partial_\ell \widehat Q(u,x), $$ is the plug-in series estimator of $S(u,x)$, $\widehat Q(u,x)$ is the series estimator of $Q(u,x)$, $\widehat \sigma_S(u,x)$ is a delta method estimator of the asymptotic standard deviation of $\widehat S(u,x)$, and $I \subseteq \mathcal{U}\times \mathcal{X} $ denotes the set of values of interest. We reject $H_0$ if the p-value of $K$ under $H_0$ is less than $\alpha$, i.e. $\sup_{P \in H_0} P(K > k) < \alpha$ where $k$ is the realized value of $K$. Dette, Hoderlein, and Neumeyer DHN2011 proposed an alternative test based on kernel estimators of the quantile function and its derivatives.
We estimate the distribution of $K$ under $H_0$ by weighted bootstrap with moment selection to reduce the asymptotic non-similarity on the boundary of composite one sided functional tests (Linton, Song, and Whang LSW2010). The weighted bootstrap version of $K$ is $$ K^*(c_n) = \max_{(u,x) \in I} \frac{\widehat S^*(u,x) - \widehat S(u,x)}{\widehat \sigma_S(u,x)} 1[|\widehat S(u,x) | < c_n \widehat \sigma_S(u,x) ], $$ where $$ \widehat S^*(u,x) = \exp(\ell) \partial_w \widehat Q^*(u,x) + \exp(\widehat Q^*(u,x) + w) \partial_\ell \widehat Q^*(u,x), $$ is the bootstrap version of $\widehat S(u,x)$, $\widehat Q^*(u,x)$ is the series estimator of $Q(u,x)$ in the weighted sample, $1[|\widehat S(u,x) | < c_n \widehat \sigma_S(u,x) ]$ is the moment selector (Chernozhukov, Hong, and Tamer CHT2007, and Andrews and Soares AS2010), and $c_n$ is a sequence of thresholds that can grow with $n$. The centering by $ \widehat S(u,x)$ imposes the least favorable null hypothesis $S(u,x) = 0$ at all the points $(u,x) \in I$ in the bootstrap to control the size of the test, whereas the moment selector discards points that are far from this hypothesis with very high probability to increase power. We consider three sequences for $c_n$: no moment selection with $c_n = 0,$ BIC moment selection with $c_n^2 = \log n,$ and LIL selection with $c_n^2 = 2 \log \log n.$ The estimator of the p-value for a realization of the statistic $k$ is the probability that $K^*(c_n)$ is greater than $k$ conditional on the data.
Table (ref) presents the results of the test of the Slutsky condition in our data set. We set the region $I$ to the product of $\{0.01, 0.02, ..., 0.99\}$ and the observed support of $X$ in the data. We obtain the p-values by weighted bootstraps with standard exponential weights and 199 replications. Here, we do not find sufficient evidence to reject the Slutsky condition at standard significance levels in any of the specifications with or without the moment selection.
To evaluate the performance of our estimation and inference methods in finite samples, we conduct a Monte Carlo experiment designed to mimic the previous empirical example. We consider the following design for the data generating process:
where $g(w) = \alpha_0 + \alpha_1 w + \alpha_2 \sin(2\pi w) + \alpha_3 \cos(2\pi w) + \alpha_4 \sin(4\pi w) + \alpha_5 \cos(4\pi w),$ $V$ is the same covariate vector as in the empirical example, $U \sim U(0, 1),$ and $\Phi^{-1}$ denotes the inverse of the CDF of the standard normal distribution. The parameters of $g(w)$ and $\beta$ are calibrated by applying least squares to the data set in the empirical example and $\sigma$ is calibrated to the least squares residual standard deviation. We consider linear, power and B-spline series methods to approximate $g(w)$, with the same number of series terms and other tuning parameters as in the empirical example. In practice, we recommend to conduct this type of Monte Carlo experiment with a data generating process that mimics the application at hand to verify the plausibility of the regularity conditions of the method.
Figures (ref) and (ref) in the Supplemental Material examine the quality of the series approximations in the population. They compare the true quantile function $$(u,\exp(w)) \mapsto \theta(u,w) = g(w) + v'\beta + \sigma \Phi^{-1}(u),$$ and the quantile price elasticity function $$(u,\exp(w)) \mapsto \theta(u,w) = \partial_w g(w),$$ to the estimands of the series approximations. In the quantile demand function the value of $v$ is fixed at the sample median values of the ordinal variables and at one for the dummies corresponding to the sample modal values of the rest of the variables (see footnote (ref)). The estimands are obtained numerically from a mega-sample (a proxy for infinite population) of $100 \times 5,001$ observations with the values of $(W,V)$ as in the data set (repeated 100 times) and with $Y$ generated from the DGP ((ref)). Although the derivative function does not depend on $u$ in our design, we do not impose this restriction on the estimands. Both figures show that the power and B-spline estimands are close to the true target functions, whereas the more parsimonious linear approximation misses important curvature features of the target functions, especially in the elasticity function.
To analyze the properties of the inference methods in finite samples, we draw 500 samples from the DGP ((ref)) with 4 sample sizes, $n$: $10,002$, $5,001$, $1,000,$ and $500$ observations. For $n=10,002$ we fix $W$ to the values in the data set repeated twice, for $n=5,001$ we fix $W$ to the values in the data set, whereas for the smaller sample sizes we draw $W$ with replacement from the values in the data set and keep them fixed across samples. To speed up computation, we drop the vector $V$ by fixing it at the sample median values of the ordinal components and at one for the dummies corresponding to the sample modal values for all the individuals. We focus on the average quantile price elasticity function $$ u \mapsto \theta(u) = \int \partial_w g(w) d\mu(w), $$ over the region $I = [0.1, 0.9]$. We estimate this function using linear, power and B-spline quantile regression with the same number of terms and other tuning parameters as in the empirical example. Although $\theta(u)$ does not change with $u$ in our design, again we do not impose this restriction on the estimators. For inference, we compare the performance of 90% confidence bands for the entire elasticity function. These bands are constructed using the pivotal, gradient bootstrap, Gaussian, and weighted bootstrap methods, all implemented in the same fashion as in the empirical example. The interval $I$ is approximated by a finite grid of 91 quantiles $\tilde{I} = \{0.10, 0.11, ..., 0.90\}$.
Table 2 reports estimation and inference results averaged across 500 simulations. The true value of the elasticity function is $\theta(u) = -0.74$ for all $u \in \tilde I$. Bias and RMSE are the absolute bias and root mean squared error integrated over $\tilde{I}$. SE/SD reports the ratios of empirical average standard errors to empirical standard deviations. SE/SD uses the analytical standard errors from expression ((ref)). The bandwidth for $\widehat J(u)$ is chosen using the Hall-Sheather option of the \verb"quantreg" \verb"R" package (Hall and Sheather Hall-Sheather1988). Length gives the empirical average of the length of the confidence band. SE/SD and length are integrated over the grid of quantiles $\tilde{I}$. Cover reports empirical coverage of the confidence bands with nominal level of 90%. Stat is the empirical average of the 90% quantile of the maximal t-statistic used to construct the bands. Table 2 shows that the linear estimator has higher absolute bias than the more flexible power and B-spline estimators, but displays lower rmse, especially for small sample sizes. The analytical standard errors provide good approximations to the standard deviations of the estimators. The confidence bands have empirical coverage close to the nominal level of 90% for all the estimators and sample sizes considered; and both bootstrap bands tend to have larger average length than the pivotal and Gaussian bands. The source of this difference in coverage might be that the bootstrap methods resample the distribution of the covariates, whereas the pivotal and Gaussian methods condition on the distribution in the sample.
All in all, these results strongly confirm the practical value of the theoretical results and methods developed in the paper. They also support the empirical example by verifying that our estimation and inference methods work quite nicely in a very similar setting.