The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
123,562 characters
Conditional Quantile Processes based on Series or Many Regressors
\title[]{Conditional Quantile Processes based on Series or Many Regressors}
\author[]{Alexandre Belloni \and Victor Chernozhukov \and Denis Chetverikov \and \\ Iv\'an Fern\'andez-Val}
\date{\tiny{\sc first version} May, 2006, {\sc this version} of \today.
The main results of this paper, particularly the pivotal method for
inference based on the entire quantile regression process, were
first presented at the NAWM Econometric Society, New Orleans, January, 2008
and also at the Stats in the Chateau in September 2009. We
are grateful to Arun Chandraksekhar, Ye Luo, Denis Tkachenko, and Sami Stouli for
careful readings of several versions of the paper. We thank the journal editor, two anonymous referees, Gary Chamberlain, Andrew Chesher, Holger Dette,
Roger Koenker, Tatiana Komarova, Arthur Lewbel, Oliver Linton, Whitney Newey, Zhongjun Qu, and seminar participants at the Econometric Society meeting, CEME Econometrics of Demand conference, CEMMAP master-class, CIREQ High-Dimensional Problems in Econometrics conference, ERCIM conference, ISI World Statistics Congress, Oberwolfach Frontiers in Quantile Regression workshop, Stats in the Chateau, Austin, BU, CEMFI, Columbia, Duke, Harvard/MIT, NUS, Rutgers, Sciences Po, SMU, Upenn, Virginia, and Yale for many
useful suggestions. We are grateful to Adonis Yatchew for giving us permission to use the data set in the empirical application. We gratefully acknowledge research support from
the NSF. The \texttt{R} package \texttt{quantreg.nonpar} implements some of the methods of this paper \cite{RJ06}.}
\maketitle
\begin{abstract}
\footnotesize{Quantile regression (QR) is a principal regression
method for analyzing the impact of covariates on outcomes. The
impact is described by the conditional quantile function and its
functionals. In this paper we develop the nonparametric QR-series
framework, covering many regressors as a special case, for
performing inference on the entire conditional quantile function and
its linear functionals. In this framework, we approximate the
entire conditional quantile function by a linear combination of
series terms with quantile-specific coefficients and estimate the
function-valued coefficients from the data. We develop large sample
theory for the QR-series coefficient process, namely we obtain
uniform strong approximations to the QR-series coefficient
process by conditionally pivotal and Gaussian processes. Based on these two strong approximations, or couplings, we develop four resampling methods (pivotal, gradient bootstrap, Gaussian, and weighted bootstrap) that can be used for inference on the entire QR-series coefficient function.
We apply these results to obtain estimation and inference methods
for linear functionals of the conditional quantile function, such as
the conditional quantile function itself, its partial derivatives,
average partial derivatives, and conditional average partial
derivatives. Specifically, we obtain uniform rates of convergence and show how to use the four resampling methods mentioned above for inference on the functionals. All of the above results are for function-valued
parameters, holding uniformly in both the quantile index and the
covariate value, and covering the pointwise case as a by-product. We
demonstrate the practical utility of these results with an
empirical example, where we estimate the price elasticity function
and test the Slutsky condition of the individual demand for gasoline, as indexed by the individual
unobserved propensity for gasoline consumption. }
\end{abstract}
\section{Introduction}
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 \cite{AHS-E2001}, Buchinsky \cite{Buchinsky1994} and Koenker \cite{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 \cite{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 \textit{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 \cite{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
\cite{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 \cite{W97}, Huang \cite{H98}, Cattaneo and Farrell \cite{CF13}, Belloni et al \cite{BelloniChenChernozhukov2009}, and Chen and Christensen \cite{CC13}. In particular, our result on optimality in the sup norm is a major extension of Huang's work \cite{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
\begin{itemize}
\item[(i)] the conditional quantile function itself, $(u,x) \mapsto Q(u, x)$,
\item[(ii)] the partial derivative function, $(u,x) \mapsto \partial_{x_k} Q(u, x)$,
\item[(iii)] the average partial derivative function, $u \mapsto \int \partial_{x_k} Q(u, x) d \mu(x)$, and
\item[(iv)] the conditional average partial derivative, $(u,x_{k}) \mapsto \int \partial_{x_k} Q(u, x) d \mu(x|x_k)$,
\end{itemize}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 \cite{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 \cite{Lee2003} that studied QR estimation
of partially linear models in the series framework for a single quantile index $u$, and on Horowitz and Lee \cite{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 \cite{BCCF11}, Chao, Volgushev and Cheng \cite{Chao-17} developed related results for QR-series coefficient processes based on weak approximations.\footnote{We refer the reader to \cite{Chao-17} for a more detailed comparison with our results.} In a more applied side, Koenker and Schorfheide \cite{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 \cite{S86},
Chaudhuri \cite{Chaudhuri-1991}, Chaudhuri, Doksum and Samarov
\cite{Chaudhuri-Doksum-Samarov-1997}, H\"{a}rdle, Ritov, and Song
\cite{hardle-ritov-song2009}, Horderlein and Mammen \cite{HoderleinMammen2009}, Cattaneo, Crump, and Jansson
\cite{cattaneo-crump-jansson2010}, Kong, Linton, and Xia
\cite{kong-linton-xia-2010}, Qu and Yoon \cite{QuYoon2011}, and Guerre and Sabbah \cite{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 \cite{Stone1982}, Andrews
\cite{Andews1991}, Newey \cite{W97}, Chen and Shen
\cite{Chen-Shen1998}, Chen \cite{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 \cite{Chen2006} and van de Geer \cite{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 \cite{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 \cite{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 \cite{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 \cite{DL93}). Third, as we demonstrate in the empirical example in Section \ref{sub: empirical example}, 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 \cite{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 \cite{belloni2013valid}.
{\bf Plan of the paper.} The rest of the paper is organized as follows.
In Section \ref{Sec:Series}, we describe the nonparametric QR-series
model and estimators. In Section \ref{Sec:theory_coeff}, we derive
asymptotic theory for the QR-series processes. In Section
\ref{Sec:Functionals}, we give estimation and inference theory for
linear functionals of the conditional quantile function. In Section \ref{Sec:examples}, 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 \cite{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.
\section{Quantile Regression Series Framework}\label{Sec:Series}
\subsection{Model}
We consider the sequence of models indexed by the sample size $n$:
\begin{equation}\label{eq: model}
Y_{i,n} = Q_n(U_{i,n},X_{i,n}),\quad i=1,\dots,n,
\end{equation}
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 \eqref{eq: model} covers two specifications of major interest:
\begin{itemize}
\item [1.] {\bf Nonparametric model:} In this model, the data generating process does not depend on $n$, that is, for all $n\geq 1$, we have $Y^n = Y$, $X^n = X$, $U^n = U$, and $Q_n = Q$ for a response variable $Y$, a $d$-dimensional vector of covariates $X$ with the support $\mathcal{X}\subset {\Bbb{R}}^d$, an unobservable individual ranking $U$ satisfying $U|X\sim {\text{Uniform}}(0,1)$, and a function $Q\colon[0,1]\times\mathcal{X}\to {\Bbb{R}}^d$ with the property that for all $x\in\mathcal{X}$, the function $u\mapsto Q(u,x)$ is strictly increasing. We assume that the functions $x\mapsto Q(u,x)$ are smooth but we do not impose any parametric structure on them. Throughout the paper, we refer to this specification as the NP model.
\item [2.] {\bf Many regressors model:} In this model, the dimension $d_n$ is allowed to grow with $n$ but the function $Q_n$ is linear in its second argument: $Q_n(u,x) = x'\beta_n(u)$ for some vector of coefficients $\beta_n(u)\in {\Bbb{R}}^{d_n}$ and all $u\in\mathcal U$ and $x\in\mathcal{X}_n$. Throughout the paper, we refer to this specification as the MR model.
\end{itemize}
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 \cite{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 \cite{CJN15} for some recent advances and also Mammen \cite{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$.
\subsection{QR-Series Approximation}
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 \textit{QR-series approximation problem}:
\begin{equation}\label{define beta}
\min_{\beta \in \Bbb{R}^{m}}E\Big[\rho_{u} (Y - Z(X)'\beta) - \rho_u(Y - Q(u,X))\Big],
\end{equation}
where $\rho_{u} (z) = (u - 1\{z<0\})z$ is the check function (Koenker \cite{K2005}).\footnote{The optimization problem \eqref{define beta} 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 \eqref{define beta} 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 \cite{W97} and Chen \cite{Chen2006} for a careful and detailed description of these series functions; see also Belloni et al \cite{BelloniChenChernozhukov2009} for an overview of recent advances on series approximating functions.
We define the \textit{QR-series approximating function} $x\mapsto Z(x)'\beta(u)$ mapping $\mathcal{X}$ into $ {\Bbb{R}}$, and, for all $x\in\mathcal{X}$, the \textit{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)$.
\subsection{QR-Series Estimator} The QR-series approximation motivates the {\em QR-series estimator} of the function $x\mapsto Q(u,x)$:
\begin{equation}\label{eq: quantile regression estimator}
\widehat Q(u,x) = Z(x)'\widehat\beta(u),\quad x\in\mathcal{X},
\end{equation}
where $\widehat\beta(u)$ is the Koenker and Bassett \cite{KB78} estimator of $\beta(u)$ that solves the empirical analog of the population problem (\ref{define beta}):
\begin{equation}\label{eq: empirical analog problem}
\min_{\beta\in {\Bbb{R}}^m}\mathbb{E}_n[\rho_u(Y_i - Z_i'\beta)],
\end{equation}
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 \eqref{eq: empirical analog problem} for all $u\in\mathcal U$ to obtain the \textit{QR-series coefficient process}
$$
\widehat \beta (\cdot) = \{ \widehat \beta(u) \colon u \in \mathcal{U}\}
$$
and the QR-series estimator \eqref{eq: quantile regression estimator} 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 \cite{PortnoyKoenker97}, who developed interior points methods with preprocessing for the problem \eqref{eq: empirical analog problem} that are very efficient and give the solution for multiple quantile indices simultaneously.
\subsection{Main Regularity Conditions}
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:
\begin{samepage}
\noindent
\textbf{Condition S.}\text{ }\textit{\begin{itemize}
\item[S.1] The data form a triangular array of random variables so that for any given $n$, the data $\mathcal{D}_n = \{(X_i,Y_i) : 1 \leq i \leq n\}$ is an i.i.d. random sample from the distribution of the pair $(X,Y)$.
\item[S.2] (i) The conditional density $f_{Y|X}(y|x)$ is bounded from above uniformly over $y\in\mathcal Y_x$, $x \in \mathcal{X}$, and $n$; (ii) $f_{Y|X}(Q(u,x)|x)$ is bounded away from zero uniformly over $u\in\bar{\mathcal U}$, $x\in\mathcal{X}$, and $n$; and (iii) the derivative of $y\mapsto f_{Y|X}(y|x)$ is continuous and bounded in absolute value from above uniformly over $y\in\mathcal Y_x$, $x \in \mathcal{X}$, and $n$.
\item[S.3] The eigenvalues of the Gram matrix $\Sigma=E[Z(X) Z(X)'] $ are bounded
from above and away from zero uniformly over $n$.
\item[S.4] The approximation error $R(u,x)$ satisfies $\sup_{x \in \mathcal{X}, u \in \mathcal{U}} |R(u,x)| \lesssim m^{-\kappa}$.
\end{itemize}
}
\end{samepage}
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 \eqref{eq: first derivative Q}, 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{Lemma:AUX_L2sparse}) 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 \cite{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
\begin{equation}\label{eq: holder smoothness}
|D^\alpha f(x) - D^\alpha f(\widetilde x)|\leq C\Big(\textstyle{\sum_{j=1}^d} (x_j - \widetilde x_j)^2\Big)^{(s - [s])/2}\text{ and }|D^\beta f(x)|\leq C
\end{equation}
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 \eqref{eq: holder smoothness} 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 \cite{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.
\begin{lemma}[Verification of Condition S.4 in the NP model for polynomials and B-splines]\label{lem: approximation error}
Consider the NP model. Suppose that Conditions S.2 and S.3 hold. In addition, suppose that $\mathcal{X} = [0,1]^d$. Moreover, suppose that $Q(u,\cdot)\in\Omega(s,C,\mathcal{X})$ for all $u\in\mathcal U$ and some $s,C>0$. If the vector of approximating functions $Z$ consists of tensor products of polynomials and $s>d$, then
\begin{equation}\label{eq: approximation error 1}
(E[|R(u,X)|^2])^{1/2}\lesssim m^{-s/d}\text{ and }\sup_{x\in\mathcal{X}}|R(u,x)|\lesssim m^{1-s/d},
\end{equation}
uniformly over $u\in\mathcal U$. Also, if the vector of approximating functions $Z$ consists of tensor products of B-splines of order $s_0$, $s\wedge s_0 > d$, and $X$ has the pdf $f_X(x)$ bounded from above and away from zero uniformly over $x\in\mathcal{X}$, then
\begin{equation}\label{eq: approximation error 2}
(E[|R(u,X)|^2])^{1/2}\lesssim m^{-(s\wedge s_0)/d}\text{ and }\sup_{x\in\mathcal{X}}|R(u,x)|\lesssim m^{-(s\wedge s_0)/d},
\end{equation}
uniformly over $u\in\mathcal U$. Thus, under the presented conditions, Condition S.4 is satisfied with $\kappa= s/d - 1$ in the case of polynomials and with $\kappa = (s\wedge s_0)/d$ in the case of B-splines.
\end{lemma}
\begin{remark}[Importance of Lemma \ref{lem: approximation error}]
Lemma \ref{lem: approximation error} makes precise the nature of the QR-series approximation in the NP model and plays a crucial role in our derivation of the convergence rate of the QR-series estimator for the NP model in the next section. Indeed, it is well-known from the approximation theory that under the assumptions of the lemma, in the case of polynomials, for example, there exists $\beta^m(u)$ such that $\sup_{x\in\mathcal{X}}|Q(u,x) - Z(x)'\beta^m(u)|\lesssim m^{-s/d}$; see Chen \cite{Chen2006}. However, this result does not help in our analysis because the QR-series estimator $\widehat\beta_n(u)$ converges in probability to $\beta(u)$, which may or may not be equal to $\beta^m(u)$. We therefore need to derive a bound on $\sup_{x\in\mathcal{X}}|Q(u,x) - Z(x)'\beta(u)|$.
The part of the lemma concerning the B-splines case is particularly important because it allows us to prove in the next section that the QR-series estimator based on B-splines achieves the fastest possible rate of convergence in the sup norm. This part of the lemma is a major extension of a result in Huang \cite{H03}, who obtained similar inequalities with the QR-series approximation error replaced by the least-squares-series approximation error.
\qed
\end{remark}
\begin{remark}[Other series approximating functions]
Other popular choices of the series approximating functions include Fourier series and compactly supported wavelets. Although we do not provide formal results for these choices, we note that under conditions similar to those in Lemma \ref{lem: approximation error}, one can show that Condition S.4 holds with $\kappa = 1/2 - s/d$ in the case of Fourier series and with $\kappa = -(s\wedge s_0)/d$ in the case of compactly supported wavelets, where $s_0$ is the order of the wavelets.
\qed
\end{remark}
\subsection{Additional Notation}\label{lab: additional notation}
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 \cite{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:
\begin{equation}\label{eq: J matrix}
J(u) = E\Big[f_{Y|X}(Q(u,X) | X) Z(X) Z(X)'\Big], \quad u\in\mathcal U.
\end{equation}
Implementing some of our inference methods will require an estimator of $J(u)$. For the purposes of this paper, we will use Powell's \cite{Powell1984} estimator defined by
\begin{equation}\label{inf-est-hatJ}
\widehat J(u) = \frac{1}{2h} \mathbb{E}_n\Big[
1\{ |Y_i - Z_i'\widehat\beta(u)| \leq h \} \cdot Z_i Z_i'\Big],
\end{equation}
where $h$ is some bandwidth value satisfying $h = h_n \to 0$. We will also use the estimator of $\Sigma$ defined by
\begin{equation}\label{inf-est-hatSigma0}
\widehat \Sigma = \mathbb{E}_n[ Z_i Z_i'].
\end{equation}
The properties of $\widehat J(u)$ and $\widehat \Sigma$ in our high-dimensional setting are established in Lemma \ref{covariance} in Appendix \ref{App:EmpiricalProcess} of the Supplemental Material.
\section{Asymptotic Theory for QR-Series Coefficient Processes} \label{Sec:theory_coeff}
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.
\subsection{Uniform-in-$u$ Rate of Convergence}
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 \cite{KB78} estimator that solves the
empirical analog \eqref{eq: empirical analog problem} of the population problem (\ref{define beta}).
Our first main result is a uniform-in-$u$ rate of convergence for the
QR-series coefficient process.
\begin{theorem}[Uniform-in-$u$ rate of convergence for QR-Series coefficient process]\label{Thm:SeriesRates} Suppose that Condition S holds. In addition, suppose that $m \zeta_m^2\log^2 n =o(n)$ and $m^{-\kappa}\log n = o(1)$. Then
$$ \sup_{u \in
\mathcal{U}} \| \widehat \beta(u) - \beta(u) \|
\lesssim_P \sqrt{m/n}.$$
\end{theorem}
Theorem \ref{Thm:SeriesRates} 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 \cite{KoenkerPortnoy1987} established rate of convergence results that hold uniformly over $\mathcal{U}$ in the fixed-dimensional setting, and He and Shao \cite{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{Thm:SeriesRates} 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 \cite{HS00}.
Theorem \ref{Thm:SeriesRates} 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{Thm:SeriesRates}, which is the second main result together with Corollary \ref{cor: sup rate qr series estimator} below on the uniform-in-$u$ rate of convergence in the sup norm of the QR-series estimator in the NP model.
\begin{corollary}[Uniform-in-$u$ $L^2$ rate of convergence for QR-series estimator in the NP model]\label{cor: l2 rate qr series estimator}
Consider the NP model. Suppose that (i) Condition S.1-3 holds. In addition, suppose that (ii) $\mathcal{X} = [0,1]^d$ and that (iii) $Q(u,\cdot)\in\Omega(s,C,\mathcal{X})$ for all $u\in\mathcal U$ and some $s,C>0$. If the vector of approximating functions $Z$ consists of tensor products of polynomials, $m^3\log^2 n = o(n)$, and $m^{1 - s/d}\log n = o(1)$, then
\begin{equation}\label{eq: L2 rate polynomials}
\sup_{u\in\mathcal U}\|\widehat Q(u,\cdot) - Q(u,\cdot)\|_{L^2(X)} \lesssim_P \sqrt{m/n} + m^{-s/d}.
\end{equation}
Also, if the vector of approximating functions $Z$ consists of tensor products of B-splines of order $s_0$, $s\wedge s_0 > d$, $m^2 \log^2 n = o(n)$, $m^{-(s\wedge s_0)/d} \log n = o(1)$, and $X$ has the pdf $f_X(x)$ bounded from above and away from zero uniformly over $x\in\mathcal{X}$, then
\begin{equation}\label{eq: L2 rate bsplines}
\sup_{u\in\mathcal U}\|\widehat Q(u,\cdot) - Q(u,\cdot)\|_{L^2(X)} \lesssim_P \sqrt{m/n} + m^{-(s\wedge s_0)/d}.
\end{equation}
\end{corollary}
\begin{remark}[QR-series estimator achieves the fastest possible rate of convergence in the $L^2$ norm]
Consider the NP model and suppose that conditions (i)--(iii) of Corollary \ref{cor: l2 rate qr series estimator} hold. If $Z$ consists of a tensor product of polynomials and $s>d$, setting $m = C n^{d/(d+2 s)}$ for some constant $C>0$ satisfies conditions that $m^3\log^2 n = o(n)$ and $m^{1 - s/d}\log n = o(1)$, and so substituting this $m$ into the bound \eqref{eq: L2 rate polynomials} gives
\begin{equation}\label{eq: optimal L2 rate of convergence}
\sup_{u\in\mathcal U}\|\widehat Q(u,\cdot) - Q(u,\cdot)\|_{L^2(X)} \lesssim_P n^{-s/(d + 2s)}.
\end{equation}
Similarly, if $Z$ consists of a tensor product of B-splines of order $s_0$, $s_0\geq s>d$, and $X$ has the pdf $f_X(x)$ bounded from above and away from zero uniformly over $x\in\mathcal X$, setting $m = C n^{d/(d+2 s)}$ for some constant $C>0$ satisfies conditions that $m^2\log^2 n = o(n)$ and $m^{-(s\wedge s_0)/d}\log n = o(1)$, and so substituting this $m$ into the bound \eqref{eq: L2 rate bsplines} again gives \eqref{eq: optimal L2 rate of convergence}. Note that the rate in \eqref{eq: optimal L2 rate of convergence} is the optimal $L^2$ rate of convergence for the estimators of nonparametric conditional quantile functions; see Chaudhuri \cite{C91}. Thus, the QR-series estimator based on either polynomials or B-splines has the fastest possible $L^2$ rate of convergence, and as we demonstrate, this rate is actually achieved uniformly in $u\in\mathcal U$. The same results can also be shown for the QR-series estimator based on Fourier series and compactly supported wavelets. This is one of the attractive properties of the QR-series estimator.\qed
\end{remark}
\subsection{Uniform Strong Approximations (Couplings) and Resampling Methods}
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{app:algorithms}.
\subsubsection{Pivotal Coupling:} Let
\begin{equation}\label{Def:U}
\mathbb{U}(u) = \frac{1}{\sqrt{n}} \sum_{i=1}^n Z_i (u-1\{U_i \leq u\}),\quad u\in\mathcal U.
\end{equation}
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{theorem: strong} 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\}$.
\begin{theorem}[Pivotal Coupling]\label{Thm:MainULA}
Suppose that Condition S holds. In addition, suppose that $m^3\zeta_m^2 =o(n^{1-\varepsilon})$ and $m^{-\kappa+1} = o(n^{-\varepsilon})$ for some constant $\varepsilon>0$. Then
$$
\sqrt{n}\left(\widehat \beta(u) -
\beta(u)\right) = J^{-1}(u)\mathbb{U}(u) +
r(u),\quad u\in\mathcal U,
$$
where
$$
\sup_{u \in \mathcal{U}}\| r(u) \| \lesssim_P \frac{m^{3/4} \zeta_m^{1/2} \log^{1/2} n}{n^{1/4}} + \sqrt{m^{1-\kappa}\log n} = o(n^{-\varepsilon'})
$$
for some $\varepsilon'>0$.
\end{theorem}
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.
\begin{corollary}[Uniform-in-$u$ sup-rate of convergence for QR-series estimator in the NP model]\label{cor: sup rate qr series estimator}
Consider the NP model. Suppose that (i) Condition S.1-3 holds. In addition, suppose that (ii) $\mathcal{X} = [0,1]^d$ and that (iii) $Q(u,\cdot)\in\Omega(s,C,\mathcal{X})$ for all $u\in\mathcal U$ and some $s,C>0$. If the vector of approximating functions $Z$ consists of tensor products of polynomials and for some $\varepsilon > 0$, $m^5 = o(n^{1-\varepsilon})$ and $m^{2 - s/d} = o(n^{-\varepsilon})$, then
$$
\sup_{u\in\mathcal U} \sup_{x\in\mathcal X}|\widehat Q(u,x) - Q(u,x)| \lesssim_P \sqrt{m^2\log n/n} + m^{1-s/d}.
$$
Also, if the vector of approximating functions $Z$ consists of tensor products of B-splines of order $s_0$, $X$ has the pdf $f_X(x)$ bounded from above and away from zero uniformly over $x\in\mathcal{X}$, and for some $\varepsilon > 0$, $m^4 = o(n^{1 - \varepsilon})$ and $m^{1 - (s\wedge s_0)/d} = o(n^{-\varepsilon})$, then
\begin{equation}\label{eq: sup rate bsplines}
\sup_{u\in\mathcal U} \sup_{x\in\mathcal X}|\widehat Q(u,x) - Q(u,x)| \lesssim_P \sqrt{m\log n/n} + m^{-(s\wedge s_0)/d}.
\end{equation}
\end{corollary}
\begin{remark}[B-splines version of QR-series estimator achieves the fastest possible rate of convergence in the sup norm]
Consider the NP model and suppose that conditions (i)--(iii) of Corollary \ref{cor: sup rate qr series estimator} hold. In addition, suppose that $Z$ consists of a tensor product of B-splines of order $s_0$, $s_0\geq s > 3d/2$, and $X$ has the pdf $f_X(x)$ bounded from above and away from zero uniformly over $x\in\mathcal X$. Then setting $m = C(n/\log n)^{d/(d + 2s)}$ for some constant $C>0$ satisfies conditions that $m^4 = o(n^{1-\varepsilon})$ and $m^{1 - (s\wedge s_0)/d} = o(n^{-\varepsilon})$ for some $\varepsilon>0$, and so substituting this $m$ into the bound \eqref{eq: sup rate bsplines} gives
$$
\sup_{u\in\mathcal U}\sup_{x\in\mathcal X}|\widehat Q(u,x) - Q(u,x)|\lesssim_P \left(\frac{\log n}{n}\right)^{s/(d + 2s)},
$$
which is the optimal rate of convergence in the sup norm for an estimator of the nonparametric conditional quantile function; see Chaudhuri \cite{C91}. Thus, the QR-series estimator based on B-splines has the fastest possible rate of convergence in the sup norm, and as we demonstrate, this rate is actually achieved uniformly in $u\in\mathcal U$.\footnote{The same results can also be shown for the QR-series estimator based on compactly supported wavelets.} This is another attractive property of the QR-series estimator. \qed
\end{remark}
We also note that although the uniform convergence rate based on polynomials is not optimal, the rate derived in Corollary \ref{cor: sup rate qr series estimator} is faster than the (trivial) uniform rate implied by the $L_2$ rate and the relation between the $L_2$-norm and sup-norm.
\subsubsection{Resampling Methods Based on Pivotal Coupling:} Another implication of Theorem \ref{Thm:MainULA} 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
\begin{equation}\label{Def:U*}
\mathbb{U}^*(u) =
\frac{1}{\sqrt{n}} \sum_{i=1}^n Z_i (u-1\{U^*_i \leq u\}),\quad u\in\mathcal U,
\end{equation}
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 \eqref{inf-est-hatJ} of Section \ref{lab: additional notation} (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{Thm:MainULA}, also of the process $\sqrt n(\widehat\beta(\cdot) - \beta(\cdot))$.
\begin{theorem}[Pivotal Method]\label{Thm:MainULAfeasible}
Suppose that Condition S holds. In addition, suppose that $h\sqrt m = o(n^{-\varepsilon})$, $m^2 \zeta_m^2 = o(n^{1-\varepsilon}{h})$, and $m^{-\kappa+1/2} = o(n^{-\varepsilon})$ for some constant $\varepsilon > 0$. Then
$$
\widehat J^{-1}(u)\mathbb{U}^*(u) = J^{-1}(u)\mathbb{U}^*(u) +
r(u),\quad u\in\mathcal U,
$$
where
$$
\sup_{u \in \mathcal{U}}\| r(u) \| \lesssim_P \sqrt{\frac{\zeta_m^2 m^2\log n}{n{h}}} + m^{-\kappa+1/2} + {h}\sqrt{m} = o(n^{-\varepsilon'})
$$
for some $\varepsilon'>0$. The stated bound continues to hold in $P$-probability if we replace
the unconditional probability $P$ by the conditional probability
$P^*$.
\end{theorem}
This theorem is the fourth main result together with Theorems \ref{Thm:MainULAstar}, \ref{thm: gaussian method}, and \ref{Thm:MainBootstrap} 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
\cite{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
\begin{eqnarray}\label{define betastar}
\min_{\beta \in \Bbb{R}^m}\Big(\mathbb{E}_n
[\rho_{u} (Y_i - Z_i'\beta)] - \mathbb{U}^*(u)'\beta/\sqrt{n}\Big),
\end{eqnarray}
where $\mathbb{U}^*(u)$ is defined in (\ref{Def:U*}). 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{Thm:MainULA}, also of the process $\sqrt n(\widehat\beta(\cdot) - \beta(\cdot))$.
\begin{theorem}[Gradient Bootstrap Method]\label{Thm:MainULAstar}
Suppose that Condition S holds. In addition, suppose that $m^3\zeta_m^2=o(n^{1-\varepsilon})$ and
$m^{-\kappa+1/2} = o(n^{-\varepsilon})$ for some constant $\varepsilon>0$. Then
$$
\sqrt{n} \left(\widehat \beta^*(u) -
\widehat \beta(u)\right) = J^{-1}(u)\mathbb{U}^*(u) + r(u),
$$
where
$$
\sup_{u \in \mathcal{U}} \| r(u) \| \lesssim_P \frac{m^{3/4} \zeta_m^{1/2} \log^{1/2} n}{n^{1/4}} + m^{-\kappa+1/2} = o(n^{-\varepsilon'})
$$
for some $\varepsilon'>0$.
The stated bound continues to hold in $P$-probability if we replace
the unconditional probability $P$ by the conditional probability
$P^*$.
\end{theorem}
\begin{remark}[Comparison of pivotal and gradient bootstrap methods]
Both the pivotal and gradient bootstrap methods have their own advantages. Perhaps the main advantage of the gradient bootstrap method relative to the pivotal method is that it does not require estimating the matrices $J(u)$, $u\in \mathcal{U}$, which is important because estimating these matrices requires a potentially subjective choice of the bandwidth $h$. In fact, implementing the gradient bootstrap method does not require any choice of smoothing parameters, making it particularly convenient for empirical researchers. On the other hand, an advantage of the pivotal method relative to the gradient bootstrap method is that it is computationally simple as it does not require solving the quantile optimization problem for each simulation of the process $\mathbb U^*(\cdot)$. \qed
\end{remark}
\subsubsection{Gaussian Coupling:} 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))$.
\begin{theorem}[Gaussian Coupling]\label{theorem: strong}
Suppose that Condition S holds. In addition, suppose that $m^{7} \zeta_m^6 = o(n^{1-\varepsilon})$ and $m^{-\kappa + 1} = o(n^{-\varepsilon})$ for some constant $\varepsilon>0$. Then
$$
\sqrt n\Big(\widehat\beta(u) - \beta(u)\Big) = J^{-1}(u) G(u) + r(u),\quad u\in\mathcal U,
$$
where $G(\cdot) = G_n(\cdot)$ is a process on $\mathcal U$ that, conditionally on $(Z_i)_{i=1}^n$, is zero-mean Gaussian with a.s. continuous sample paths and the covariance function
\begin{equation}\label{eq: covariance matrix gaussian process thm 5}
E\Big[G(u_1)G(u_2)'\mid (Z_i)_{i=1}^n\Big] = \mathbb{E}_n[Z_i Z_i'](u_1\wedge u_2 - u_1 u_2), \ \text{ for all $u_1$ and $u_2$ in $\mathcal U$},
\end{equation}
and
$$
\sup_{u\in\mathcal U}\|r(u)\| = o_P(n^{-\varepsilon'})
$$
for some $\varepsilon'> 0$.
\end{theorem}
\begin{remark}[Conditions of Theorem \ref{theorem: strong}]
Note that the strong approximation to the process $\sqrt n(\widehat\beta(\cdot) - \beta(\cdot))$ by the Gaussian process $J^{-1}(\cdot)G(\cdot)$ constructed in Theorem \ref{theorem: strong} requires the condition that $m^7 \zeta_m^6 = o(n^{1 - \varepsilon})$, which is more restrictive than the corresponding condition in Theorem \ref{Thm:MainULA}, $m^3\zeta_m^2 = o(n^{1 - \varepsilon})$, required for the strong approximation to the process $\sqrt n(\widehat\beta(\cdot) - \beta(\cdot))$ by the pivotal process $J^{-1}(\cdot)\mathbb U(\cdot)$. We note that this restrictive condition is sufficient but we do not know whether it is necessary. This condition is a consequence of a step in the proof of Theorem \ref{theorem: strong} that relies upon Yurinskii's coupling. Therefore, improving that step
through the use of another coupling could potentially lead to significant
improvements in the conditions of the theorem; see, in particular, Theorem \ref{thm: gaussian coupling for t process} in the next section. See also \cite{K94} and \cite{CNS15}, where a Hungarian coupling is derived that may give a result similar to that in Theorem \ref{theorem: strong} but under somewhat weaker conditions if $d$ is small and the vector of approximating functions $Z$ consists of a tensor products of B-splines or wavelets. \qed
\end{remark}
\subsubsection{Resampling Methods Based on Gaussian Coupling:} Although Theorem \ref{theorem: strong} 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 \eqref{eq: covariance matrix gaussian process thm 5}, 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
\begin{equation}\label{eq: gaussian process to simulate}
G^*(u) = G_n^*(u) = \widehat \Sigma^{1/2} B_m^*(u),\quad u\in\mathcal U,
\end{equation}
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 \eqref{inf-est-hatJ}. 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{theorem: strong} also of the process $\sqrt n(\widehat \beta(\cdot) - \beta(\cdot))$.
\begin{theorem}[Gaussian Method]\label{thm: gaussian method}
Suppose that Condition S holds. In addition, suppose that $h\sqrt m = o(n^{-\varepsilon})$, $m^2\zeta_m^2 = o(n^{1 - \varepsilon} h)$, and $m^{-\kappa + 1/2} = o(n^{-\varepsilon})$ for some constant $\varepsilon > 0$. Then
$$
\widehat J^{-1}(u) G^*(u) = J^{-1}(u) G^*(u) + r(u),\quad u\in\mathcal U,
$$
where
$$
\sup_{u\in\mathcal U}\|r(u)\| \lesssim_P \sqrt{\frac{m^2 \zeta_m^2\log n}{n{h}}} + m^{-\kappa+1/2} + {h}\sqrt{m} = o(n^{-\varepsilon'})
$$
for some $\varepsilon' > 0$. The stated bound continues to hold in $P$-probability if we replace
the unconditional probability $P$ by the conditional probability
$P^*$.
\end{theorem}
Another related inference method is the weighted bootstrap method. Pr{\ae}stgaard and Wellner~\cite{Praestgaard-Wellner-93}, Hahn~\cite{H97}, Chamberlain and Imbens~\cite{Chamberlain-Imbens-03}, and Chen and Pouzo~\cite{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
\begin{equation}\label{eq: weighted bootstrap problem}
\widehat \beta^b(u) \in \arg \min_{\beta \in {\Bbb{R}}^m} \mathbb{E}_n[ \pi_i \rho_u(Y_i-Z_i'\beta)].
\end{equation}
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{theorem: strong} also of the process $\sqrt n(\widehat\beta(\cdot) - \beta(\cdot))$.
\begin{theorem}[Weighted Bootstrap Method]\label{Thm:MainBootstrap}
Suppose that Condition S holds. In addition, suppose that $m^{7} \zeta_m^6 = o(n^{1-\varepsilon})$ and $m^{-\kappa + 1} = o(n^{-\varepsilon})$ for some constant $\varepsilon>0$. Moreover, suppose that the random variable $\pi$ is non-negative and satisfies $E[\pi] = 1$, $E[\pi^2] = 2$, $E[\pi^4]\lesssim 1$. Finally, suppose that $\max_{1\leq i\leq n} \pi_i \lesssim \log n$. Then
\begin{equation}\label{eq: main implication thm 6}
\sqrt n\Big(\widehat\beta^b(u) - \widehat\beta(u)\Big) = J^{-1}(u)G^*(u) + r(u),
\end{equation}
where $G^*(\cdot) = G^*_n(\cdot)$ is a process on $\mathcal U$ that, conditionally on $(Z_i)_{i=1}^n$, is zero-mean Gaussian with a.s. continuous sample paths and the covariance function \eqref{eq: covariance matrix gaussian process thm 5}, and
$$
\sup_{u\in\mathcal U}\|r(u)\| \lesssim_P o(n^{-\varepsilon'})
$$
for some $\varepsilon' >0$. Moreover, the stated bound continues to hold in $P$-probability if we replace the
unconditional probability $P$ by the conditional probability $P^*$.
\end{theorem}
\begin{remark}[Comparison of Gaussian and weighted bootstrap methods]
The comparison of the Gaussian and weighted bootstrap methods is similar to that of the pivotal and gradient bootstrap methods. Again both methods have their own advantages. The main advantage of the weighted bootstrap method is arguably that it does not require estimating the matrices $J(u)$, $u\in \mathcal{U}$, which allows us to bypass the need to select a bandwidth $h$. An advantage of the Gaussian method is that it is computationally simple as it does not require solving the quantile optimization problem for each simulation of weights $(\pi_i)_{i=1}^n$.\qed
\end{remark}
\begin{remark}[Comparison of resampling methods based on the pivotal and Gaussian couplings]
Although it is difficult to compare the resampling methods based on the pivotal coupling (pivotal and gradient bootstrap methods) with those based on the Gaussian coupling (Gaussian and weighted bootstrap methods) from a theoretical point of view, our results suggest that the former methods might be more accurate than the latter ones. Indeed, the methods based on the pivotal coupling require weaker conditions (see, however, Theorem \ref{thm: gaussian coupling for t process} in the next section, where it is possible to substantially weaken conditions required for the Gaussian coupling in some examples) and, in addition, developing the Gaussian coupling requires a ``double approximation'': in order to construct a Gaussian process that strongly approximates the original process $\sqrt n(\widehat\beta(\cdot) - \beta(\cdot))$, we first construct a coupling of the latter process with the pivotal process, and then we construct a coupling of the Gaussian process with the pivotal process, so that the Gaussian process is coupled with the original process $\sqrt n(\widehat\beta(\cdot) - \beta(\cdot))$ via the pivotal process. In the numerical examples of Section \ref{subsec:simulations} and the companion computational paper \cite{RJ06}, however, we find that the performance of the four methods is similar in finite samples.
In addition, the Gaussian coupling is important because of the existence of well-developed extreme value theory for Gaussian processes; see, for example, Leadbetter, Lindgren and Rootzen \cite{LLR83}. In combination with the Gaussian coupling, this theory can be used to develop inferential procedures for some linear functionals without relying upon resampling methods (that is, with non-bootstrap critical values) like in Rio \cite{R94}. Moreover, the Gaussian coupling is important because of existence of anti-concentration inequalities for Gaussian processes (Lemma \ref{lemma:Anti}), which are useful to construct uniform confidence bands for linear functionals in the next section.
\qed
\end{remark}
\section{Linear Functionals of the Conditional Quantile Function}\label{Sec:Functionals}
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
\begin{itemize}
\item[1.] the derivative: \quad $\theta(u,x) = \partial_{x_k} Q(u,x)$;
\item[2.] the average derivative: \quad $\theta(u) = \int \partial_{x_k} Q(u,x) d\mu(x)$;
\item[3.] the conditional average derivative: \ $\theta(u,w) = \int \partial_{x_k} Q(u,w,v) d\mu(v|w)$.
\end{itemize}
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
\begin{itemize}
\item the function $\theta(u,w)$ at a particular point $(u,w)$, then $I = \{(u,w)\}$,
\item the function $u\mapsto \theta(u,w)$ having fixed $w$, then $I =\mathcal{U}\times\{w\}$,
\item the function $w\mapsto \theta(u,w)$ having fixed $u$, then $I = \{u\}\times\mathcal{W}$,
\item the entire function $(u,w)\mapsto \theta(u,w)$, then $I = \mathcal{U}\times\mathcal{W}$.
\end{itemize}
\subsection{QR-Series Approximation} 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,
\begin{equation}\label{eq: cool}
\theta(u,w) = \ell(w)' \beta(u) + r(u,w), \quad (u,w) \in I,
\end{equation}
where $\ell(w)' \beta(u)$ is the \textit{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
\textit{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$:
\begin{equation}\label{eq: uncool}
\left(\mathcal{A} Q(u,\cdot)\right)[w] = \left(\mathcal{A} Z(\cdot)\right)[w]'\beta(u) + \left(\mathcal{A} R(u,\cdot)\right)[w].
\end{equation}
In the three examples above the operator $\mathcal{A}$ is given by,
respectively,
\begin{itemize}
\item[1.] a differential operator: $(\mathcal{A}g) [x]= (\partial_{x_k}g) [x] $, so that
$$\ell(x) =\partial_{x_k} Z(x), \ \ \ r(u,x) = \partial_{x_k} R(u,x);$$
\item[2.] an integro-differential operator: $\mathcal{A} g=\int \partial_{x_k} g(x) d\mu(x)$, so that
$$\ell = \int \partial_{x_k} Z(x) d \mu(x), \ \ \ r(u) = \int \partial_{x_k} R(u,x) d \mu(x); $$
\item[3.] a partial integro-differential operator: $(\mathcal{A}g) [w]=\int \partial_{x_k} g(w,v) d\mu(v|w)$, so that
$$\ell(w) =\int \partial_{x_k} Z(w,v)d\mu(v|w), \ \ \ r(u, w) = \int \partial_{x_k} R(u,w,v) d\mu(v|w).$$
\end{itemize}
For notational convenience, we use the formulation (\ref{eq: cool}) in the analysis, instead of the motivational formulation (\ref{eq: uncool}).
\subsection{QR-Series Estimator} 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{subset: mono} 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.
\subsection{Pointwise Asymptotic Theory}
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:
\noindent
\textbf{Condition P.} \textit{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 \cite{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{lem: approximation error}, 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{Thm:DistributionsInferentialpointwise} below on pointwise asymptotic normality of $\widehat\theta(u,w)$.
\begin{theorem}[Pointwise Convergence Rate for Linear Functionals]
\label{theorem: pointwise rate} Suppose that the conditions of
Theorem \ref{Thm:MainULA} hold. In addition, suppose that Condition P holds. Then
$$
| \widehat \theta(u,w) - \theta(u,w)| \lesssim_P \frac{\|\ell(w)\|}{\sqrt{n}}.
$$
\end{theorem}
\begin{remark}[Rates and norm of vector of loadings] The rate of convergence of $\widehat \theta(u,w)$ depends on the functional $\theta(u,w)$ through the norm of the vector of loadings $\ell(w)$. For example, if we are interested in the coefficient $\beta_1(u)$, so that $\theta(u,w) = \beta_1(u)$, which might be a parameter of interest in the Many regressors (MR) model, then $\ell(w) = (1, 0, \ldots, 0)'$, and so $\|\ell(w) \| = 1$, yielding a $\sqrt{n}$-consistent estimator $\widehat\theta(u,w)$. See Comment \ref{comment: xi} below for additional examples of linear functionals with bounds on $\|\ell(w)\|$. \qed
\end{remark}
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
\begin{equation}\label{Def:hatsigma2n}
\widehat \sigma^2(u,w) = u(1-u) \ell(w)' \widehat J^{-1}(u) \widehat \Sigma \widehat J^{-1}(u)\ell(w)/n
\end{equation}
is a consistent estimator of
\begin{equation}\label{Def:sigma2n}
\sigma^2(u,w) = u(1-u) \ell(w)'J^{-1}(u) \Sigma J^{-1}(u)\ell(w)/n,
\end{equation}
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.
\begin{theorem}[Pointwise Inference for Linear Functionals]\label{Thm:DistributionsInferentialpointwise}
Suppose that the conditions of Theorem \ref{Thm:MainULA} hold. In addition, suppose that Condition P holds, $h = o(1)$ and $m\zeta_m^2\log^2 n = o(n h)$. Then
$$
t(u,w) \to_d N(0,1).
$$
\end{theorem}
\begin{remark}[Using resampling methods for pointwise inference]
Although it is possible to establish validity of all the resampling methods from the previous section to perform pointwise inference on linear functionals, we do not show these results here because they will follow as a special case from our results below on uniform inference for linear functionals. We provide an implementation algorithm to perform pointwise inference using the resampling methods in Appendix \ref{app:algorithms}. \qed
\end{remark}
\subsection{Uniform Asymptotic Theory} 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:
\noindent
\textbf{Condition U.}\textit{
\begin{itemize}
\item[U.1] The set $I$ is such that its dimension $d_I$ is fixed and its diameter is bounded uniformly over $n$.
\item[U.2] For some $\varepsilon > 0$, the QR-approximation error $r(u,w)$ satisfies
$$
\sqrt{n} \sup_{(u,w)\in I}\frac{ | r(u,w) |}{\|\ell(w)\|} = o(n^{-\varepsilon}).
$$
\item[U.3] The vector of loadings $\ell(w)$ satisfies
$$
\|\ell(w)\| \leq \zeta_{m,\theta}\text{ and } \ \left\|\frac{\ell(w)}{\|\ell(w)\|} - \frac{\ell(w')}{\|\ell(w')\|}\right\| \leq \zeta_{m,\theta}^L\|w - w'\|
$$
for all $w,w'\in\mathcal W$, where $\log \zeta_{m,\theta}^L\lesssim \log n$.
\end{itemize}
}
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.
\begin{remark}[Primitive bounds on $\zeta_{m,\theta}$]\label{comment: xi}
The uniform rate of convergence for the estimator $\widehat\theta(u,w)$ derived below in Theorem \ref{theorem: uniform rate} will crucially depend on the constant $\zeta_{m,\theta}$ appearing in Condition U. Here we discuss some bounds on this constant. For brevity, we only discuss the case of B-splines and refer to Newey \cite{W97} and Chen \cite{Chen2006} for other choices of approximating functions. We assume that $\mathcal X = [0,1]^d$ and that the vector of approximating functions $Z$ consists of tensor products of B-splines of order $s_0$. As discussed above, then $\zeta_m = \sup_{x\in\mathcal X} \|Z(x)\|\lesssim \sqrt m$ and it is also possible to verify that for all positive integers $\alpha \leq s_0$, $\sup_{x\in\mathcal X}\|\partial^{\alpha}_{x_k} Z(x)\| \lesssim m^{1/2 + \alpha / d}$; see for example Chen and Christensen \cite{CC13}. Then \\
\begin{tabular}{ll}
$\bullet$ For &$\theta(u,w) = Q(u,x)$, $\ell(w)=Z(x)$ and $\zeta_{m,\theta} \lesssim m^{1/2}$;\\
$\bullet$ For &$\theta(u,w) = \partial_{x_k} Q(u,x)$, $\ell(w)= \partial^{\alpha}_{x_k} Z(x)$ and $\zeta_{m,\theta} \lesssim m^{1/2 + \alpha/d}$;\\
$\bullet$ For &$\theta(u) = \int \partial_{x_k} Q(u,x) d\mu(x)$ with ${\rm supp}(\mu) \subset {\rm int}(\mathcal{X})$ and $|\partial_{x_k} \mu(x)| \lesssim 1$,\\
&$\ell = \int \partial_{x_k} Z(x) \mu(x) \ dx = -\int Z(x) \partial_{x_k} \mu(x) \ dx$ and $\zeta_{m,\theta} \lesssim 1$;\\
\end{tabular}
\noindent
see Newey \cite{W97} for more explanations on the last bound. \qed
\end{remark}
\subsubsection{Uniform-in-$u$ Rate of Convergence:} The following theorem establishes the uniform rate of convergence of the QR-series estimator $\widehat \theta(u,w)$, which is the sixth main result.
\begin{theorem}[Uniform Convergence Rate for Linear Functionals]
\label{theorem: uniform rate} Suppose that the conditions of Theorem
\ref{Thm:MainULA} hold. In addition, suppose that Condition U hold. Then
$$
\sup_{(u,w)\in I} | \widehat \theta(u,w) - \theta(u,w)| \lesssim_P \sqrt{\frac{\zeta^2_{m,\theta}\log n}{n}}.
$$
\end{theorem}
The uniform rate of Theorem \ref{theorem: uniform rate} is the same as the pointwise rate of Theorem \ref{theorem: pointwise rate} 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.
\begin{remark}[Comparison of Theorem \ref{theorem: uniform rate} and Corollary \ref{cor: sup rate qr series estimator}] When $\theta(u,w)$ is the conditional quantile function, the convergence rate of Theorem \ref{theorem: uniform rate} is asymptotically equivalent to the rate of Corollary \ref{cor: sup rate qr series estimator} under the undersmoothing condition U.2. For example, in the case of B-splines, $\zeta^2_{m,\theta}\log n/ n = m \log n/ n$ and $m^{-(s\wedge s_0)/d} = o(\sqrt{m\log n/n})$ under U.2.
\end{remark}
\subsubsection{Gaussian and Pivotal Couplings for $t$-Statistic Processes:} 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:
\begin{equation}\label{tnprocess}
t(u,w) = \frac{ \widehat \theta(u,w) - \theta(u,w) }{ \widehat \sigma(u,w)},
\end{equation}
where $\widehat \sigma^2(u,w)$, defined in (\ref{Def:hatsigma2n}), is an estimator
of the asymptotic variance $\sigma^2(u,w)$ of $\widehat\theta(u,w)$ in (\ref{Def:sigma2n}). 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{thm: gaussian coupling for t process}, \ref{thm: resampling methods t process}, and \ref{thm: weighted bootstrap alternative conditions} below on couplings and resampling methods for the t-statistic process.
\begin{theorem}[Pivotal and Gaussian Couplings for t-statistic Process]\label{thm: couplings for t process}
Suppose that Conditions S and U hold. If $h = o(n^{-\varepsilon})$, $m\zeta_m^2 = o(n^{1 - \varepsilon} h)$, $m^3\zeta_m^2 = o(n^{1 - \varepsilon})$, and $m^{-\kappa + 1} = o(n^{-\varepsilon})$ for some constant $\varepsilon > 0$, then
\begin{equation}\label{eq: pivotal coupling t process}
\sup_{(u,w)\in I}\left|t(u,w) - \frac{\ell(w)'J^{-1}(u)\mathbb U(u)/\sqrt n}{\sigma(u,w)}\right| \lesssim_P o(n^{-\varepsilon'})
\end{equation}
for the process $\mathbb U(\cdot)$ defined in \eqref{Def:U} for some $\varepsilon' > 0$. Also, if $h = o(n^{-\varepsilon})$, $m\zeta_m^2 = o(n^{1 - \varepsilon} h)$, $m^7\zeta_m^6 = o(n^{1 - \varepsilon})$, and $m^{-\kappa + 1} = o(n^{-\varepsilon})$ for some constant $\varepsilon > 0$, then
\begin{equation}\label{eq: gaussian coupling t process}
\sup_{(u,w)\in I}\left|t(u,w) - \frac{\ell(w)'J^{-1}(u) G(u)/\sqrt n}{\sigma(u,w)}\right| \lesssim_P o(n^{-\varepsilon'})
\end{equation}
for the process $G(\cdot)$ defined in Theorem \ref{theorem: strong} for some $\varepsilon' > 0$.
\end{theorem}
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:
\begin{theorem}[Gaussian Coupling for t-statistic Process under Alternative Conditions]\label{thm: gaussian coupling for t process}
Suppose that Conditions S and U hold. In addition, suppose that $h = o(n^{-\varepsilon})$, $m\zeta_m^2 = o(n^{1 - \varepsilon} h)$, $m^3\zeta_m^2 = o(n^{1 - \varepsilon})$, and $m^{-\kappa + 1} = o(n^{-\varepsilon})$ for some constant $\varepsilon > 0$. Moreover, suppose that $(1 + \zeta_{m,\theta}^L)^{2 d_I}\zeta_m^2 = o(n^{1 - \varepsilon})$. Then \eqref{eq: gaussian coupling t process} holds for the same process $G(\cdot)$ as that used in Theorem \ref{thm: couplings for t process}.
\end{theorem}
\begin{remark}[Comparison of conditions for the Gaussian coupling in Theorems \ref{thm: couplings for t process} and \ref{thm: gaussian coupling for t process}]
The conditions of Theorems \ref{thm: couplings for t process} and \ref{thm: gaussian coupling for t process} required for the Gaussian coupling are non-nested. In particular, Theorem \ref{thm: gaussian coupling for t process} requires the condition $m^3\zeta_m^2 = o(n^{1 - \varepsilon})$ that is weaker than the corresponding condition in Theorem \ref{thm: couplings for t process}, $m^7\zeta_m^6 = o(n^{1 - \varepsilon})$, but it also requires the condition $(1 + \zeta_{m,\theta}^L)^{2 d_I} \zeta_m^2 = o(n^{1 - \varepsilon})$ that is stronger than the corresponding condition in Theorem \ref{thm: couplings for t process}, $\log \zeta_{m,\theta}^L \lesssim \log n$. However, in most cases of practical importance, the conditions of Theorem \ref{thm: gaussian coupling for t process} are substantially weaker than those of Theorem \ref{thm: couplings for t process}. For example, consider the NP model and suppose that we are interested in the conditional quantile function $Q(u,x)$ itself, so that $\ell(\omega) = Z(x)$. Further, suppose that $\mathcal X = [0,1]^d$ and that the vector of approximating functions $z$ consists of tensor products of B-splines. Then $\zeta_m \lesssim \sqrt m$, as discussed above, and it is also possible to show that $\zeta_{m,\theta}^L \lesssim m^{1/d}$. Hence, in this case Theorem \ref{thm: gaussian coupling for t process} requires that $m^{4\vee (2/d + 3)} = o(n^{1 - \varepsilon})$ since $d_I = 1 + d$ whereas Theorem \ref{thm: couplings for t process} requires $m^{10} = o(n^{1 - \varepsilon})$.
\end{remark}
\subsubsection{Resampling Methods:} As in Section \ref{Sec:theory_coeff}, 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 \eqref{Def:U*} and \eqref{eq: gaussian process to simulate}, 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
\begin{equation*}
\begin{array}{llll}
\\
\mbox{pivotal method:} & \displaystyle t^*(u,w) = \frac{ \ell(w)' \widehat J^{-1}(u) \mathbb{U}^*(u)/\sqrt{n} }{ \widehat \sigma(u,w)};\\
\\
\mbox{gradient bootstrap method:} & \displaystyle t^*(u,w) = \frac{ \ell(w)' (\widehat \beta^*(u)- \widehat\beta(u)) }{ \widehat \sigma(u,w)}; \\
\\
\mbox{Gaussian method:} & \displaystyle t^*(u,w) = \frac{ \ell(w)' \widehat J^{-1}(u) G^*(u)/\sqrt{n} }{ \widehat \sigma(u,w)}; \\
\\
\mbox{weighted bootstrap method:} & \displaystyle t^*(u,w) = \frac{ \ell(w)' (\widehat \beta^b(u)- \widehat \beta(u)) }{ \widehat \sigma(u,w)}.
\end{array}
\end{equation*}
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{thm: couplings for t process} and \ref{thm: gaussian coupling for t process} also of the original t-statistic process $t(\cdot,\cdot)$.
\begin{theorem}[Validity of Resampling Methods for t-statistic Process]\label{thm: resampling methods t process}
Suppose that Conditions S and U hold. In addition, suppose that $h = o(n^{-\varepsilon})$ and $m\zeta_m^2 = o(n^{1 - \varepsilon} h)$ for some constant $\varepsilon > 0$. Moreover, suppose that (i) the conditions of Theorem \ref{Thm:MainULAfeasible} hold in the case of the pivotal method, (ii) the conditions of Theorem \ref{Thm:MainULAstar} hold in the case of gradient bootstrap method, (iii) the conditions of Theorem \ref{thm: gaussian method} hold in the case of Gaussian method, and (iv) the conditions of Theorems \ref{Thm:MainBootstrap} hold in the case of weighted bootstrap method. Then for the pivotal and gradient bootstrap methods,
$$
\sup_{(u,w)\in I}\left|t^*(u,w) - \frac{\ell(w)'J^{-1}(u)\mathbb U^*(u)/\sqrt n}{\sigma(u,w)}\right| \lesssim_P o(n^{-\varepsilon'})
$$
for some $\varepsilon' > 0$.
In addition, for the Gaussian and weighted bootstrap methods,
$$
\sup_{(u,w)\in I}\left|t^*(u,w) - \frac{\ell(w)'J^{-1}(u) G^*(u)/\sqrt n}{\sigma(u,w)}\right| \lesssim_P o(n^{-\varepsilon'})
$$
for some $\varepsilon' > 0$. Moreover, the stated bounds continue to hold in $P$-probability if we replace the unconditional probability $P$ by the conditional probability $P^*$.
\end{theorem}
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{thm: gaussian coupling for t process}, it turns out that it is possible to obtain the same approximation as in this theorem but under a different set of conditions:
\begin{theorem}[Weighted Bootstrap Method for t-statistic Process under Alternative Conditions]\label{thm: weighted bootstrap alternative conditions}
Suppose that Conditions S and U hold. In addition, suppose that $h = o(n^{-\varepsilon})$, $m\zeta_m^2 = o(n^{1 - \varepsilon} h)$, $m^3\zeta_m^2 = o(n^{1 - \varepsilon})$, and $m^{-\kappa + 1} = o(n^{-\varepsilon})$ for some constant $\varepsilon > 0$. Moreover, suppose that $(1 + \zeta_{m,\theta}^L)^{2 d_I}\zeta_m^2 = o(n^{1 - \varepsilon})$. Finally, suppose that the conditions of Theorem \ref{Thm:MainBootstrap} on the weights $\pi_i$ hold. Then for the weighted bootstrap method,
$$
\sup_{(u,w)\in I}\left|t^*(u,w) - \frac{\ell(w)'J^{-1}(u) G^*(u)/\sqrt n}{\sigma(u,w)}\right| \lesssim_P o(n^{-\varepsilon'})
$$
for some $\varepsilon' > 0$. Moreover, the stated bound continues to hold in $P$-probability if we replace the unconditional probability $P$ by the conditional probability $P^*$.
\end{theorem}
\subsection{Uniform Confidence Bands}
With the help of Theorems \ref{thm: couplings for t process} -- \ref{thm: weighted bootstrap alternative conditions}, 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
\begin{equation}\label{eq: exact band}
\Big\{[\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\}
\end{equation}
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 \eqref{eq: exact band} 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{app:algorithms}. 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.
\begin{theorem}[Uniform Confidence Bands]\label{theorem: inference using couplings}
Suppose that Conditions S and U hold. In addition, suppose that $m^3\zeta_m^2 = o(n^{1 - \varepsilon})$ and $m^{-\kappa + 1} = o(n^{-\varepsilon})$ for some constant $\varepsilon > 0$. Moreover, suppose that (i) $h\sqrt m = o(n^{-\varepsilon})$ and $m^2\zeta_m^2 = o(n^{1 - \varepsilon} h)$ in the case of pivotal and Gaussian methods and (ii) $h = o(n^{-\varepsilon})$ and $m\zeta_m^2 = o(n^{1 - \varepsilon} h)$ in the case of gradient and weighted bootstrap methods. Finally, suppose that the conditions of Theorem \ref{Thm:MainBootstrap} on the weights $\pi_i$ hold in the case of the weighted bootstrap method.
(1) Then
\begin{equation}\label{eq: conservative}
P \Big( V \leq k^*(1-\alpha) \Big) = 1-\alpha + o(1).
\end{equation}
(2) As a consequence,
\begin{equation}\label{eq: conservative coverage}
P \Big( \theta(u,w) \in [\dot{\iota}(u,w), \ddot{\iota}(u,w)], \mbox{ for all } (u,w) \in I \Big) = 1-\alpha + o(1).
\end{equation}
(3) The width of the confidence band $2k^*(1-\alpha) \widehat
\sigma(u,w)$ obeys
\begin{equation}\label{eq: nonconservative width}
2k^*(1-\alpha) \widehat \sigma(u,w) \lesssim_P \sqrt{\frac{\zeta_{m,\theta}^2\log n}{n}}
\end{equation}
uniformly over $(u,w)\in I$.
\end{theorem}
In addition to the validity of the uniform confidence band, Theorem \ref{theorem: inference using couplings} 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)$.
\begin{remark}[Related literature]
The construction of uniform confidence bands for nonparametric functions has been of large interest both in econometrics and statistics at least from the seventies. Early constructions can be traced back at least to the classic work \cite{BR1973} by Bickel and Rosenblatt. More recent contributions include Claeskens and Keilegom \cite{CK03}, Horowitz and Lee \cite{HL09}, Gin\'{e} and Nickl \cite{GN10}, and Chernozhukov, Chetverikov and Kato \cite{CCK2013}, among many others. Most of the constructions in the literature rely on a two-step strategy. First, the distribution of an estimator of the function of interest is approximated by some Gaussian process uniformly over its domain. Second, extreme value theory is employed to obtain the limit distribution of the supremum of the absolute value of the Gaussian process and its appropriate quantile is used to choose the width of the confidence band. A widely understood problem of this construction, however, is that the limit distribution on the second step may not exist and even if it does, it is often difficult to derive its explicit form. This distribution depends both on the function of interest and on the estimator considered, so that treatment of any new estimation problem requires a separate theorem, and in fact considerable efforts have been devoted to derive this distribution even in relatively simple settings, like density estimation based on projection kernels; see Gin\'{e} and Nickl \cite{GN10} and references therein. We avoid this problem: instead of deriving the limit distribution on the second step, we rely upon resampling methods developed in this paper. As a result, our construction yields asymptotically exact uniform confidence bands that work generically for all linear functionals of the conditional quantile functions. Our strategy is related to that used in Chernozhukov, Chetverikov and Kato \cite{CCK2013} for the problem of density estimation and is built on Chernozhukov, Lee and Rosen \cite{CLR2009}, who proposed a related strategy for inference on the minimum of a function. \qed
\end{remark}
\subsection{Test of Shape Constraints}
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:
\begin{itemize}
\item[(i)] Monotonicity of $x\mapsto Q(u,x)$ with respect to $x_k$: $\partial_{x_k}Q(u,x) \leq 0$ for all $x\in\mathcal X$ and $u\in \mathcal U$;
\item[(ii)] Concavity of $x\mapsto Q(u,x)$ with respect to $x_k$: $\partial^2_{x_k}Q(u,x) \leq 0$ for all $x\in\mathcal X$ and $u\in \mathcal U$;
\item[(iii)] Concavity of $x\mapsto Q(u,x)$ with respect to $x$: $\alpha'\partial_x^2Q(u,x)\alpha \leq 0$ for all $\alpha\in S^{d-1}$, $x\in\mathcal X$, and $u\in \mathcal U$,
\end{itemize}
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.
\begin{theorem}[Test of Shape Constraints]\label{thm: shape constraints}
Suppose that the conditions of Theorem \ref{theorem: inference using couplings} hold. Then under $H_0$,
$$
P\Big( T > \tilde k^*(1 - \alpha)\Big) \leq \alpha + o(1).
$$
Moreover, if $\mathcal M = \mathcal M_n$ is a set of data-generating processes that satisfy $H_0$ and is such that the conditions of Theorem \ref{theorem: inference using couplings} hold uniformly over this set, then
$$
\sup_{M \in \mathcal M}P_M\Big( T > \tilde k^*(1 - \alpha) \Big) \leq \alpha + o(1),
$$
where $P_M$ denotes probability under the data-generating process $M$.
\end{theorem}
\section{Examples}\label{Sec:examples}
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" (\cite{R08}), using the
package \verb"quantreg" for quantile regression (Koenker \cite{Koenker08}). We refer to Appendix \ref{app:algorithms} for implementation algorithms and to the companion computational paper \cite{RJ06} for software and additional examples.
\subsection{Empirical Example}\label{sub: empirical example}
To illustrate our methods with real data, we consider an empirical
application to nonparametric estimation of the demand for gasoline. Blundell, Horowitz and Parey \cite{BHP2012}, Hausman and Newey
\cite{HN1995}, Schmalensee and Stoker \cite{SS1999}, and Yatchew and No \cite{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 \cite{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 \texttt{ 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 \cite{YN2001} for a more detailed description of the data.) Our sample selection and variable construction also follow
Yatchew and No \cite{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 \cite{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 \cite{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 \cite{M07} and Koltchinskii \cite{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{fig: surfaces} and \ref{fig: average elasticity cis}.\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{fig: surfaces} 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{\label{ft:cov_values}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 \cite{CFG2010-Biometrika}; see Appendix \ref{subset: mono}. 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{fig: surfaces} 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 \cite{C1979}),
because the unsmoothed elasticity estimates display very erratic
behavior.}
Fig. \ref{fig: average elasticity cis} 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 \cite{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:
\begin{equation}\label{eq: slutsky}
H_0 : S(u,x) \leq 0, \text{ for all } (u,x) \in \mathcal{U} \times \mathcal{X}, \text{ vs } H_1: S(u,x) > 0, \text{ for some }(u,x) \in \mathcal{U}\times \mathcal{X},
\end{equation}
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 \cite{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 \cite{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 \cite{CHT2007}, and Andrews and Soares \cite{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{table: slutsky} 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.
\begin{table}[h]
\centering
\caption{Test of Slutsky Condition}\label{table: slutsky}
\begin{tabular}{ccccc}\hline\hline
\multicolumn{2}{c}{ } & \multicolumn{3}{c}{P-value$^{\dag}$ }\\
\cline{3-5}
Specification & $K$ stat & No selection & BIC selection & LIL selection \\
\hline
Linear & 0.47 & 0.95 & 0.76 & 0.58 \\
Power & 3.63 & 0.30 & 0.30 & 0.28 \\
B-spline & 2.30 & 0.96 & 0.96 & 0.96 \\
\hline\hline
\multicolumn{5}{l}{\footnotesize{$^{\dag}$P-values obtained by weighted bootstrap with standard exponential weights }}\\
\multicolumn{5}{l}{\footnotesize{ and 199 replications.}}
\end{tabular}
\end{table}
\subsection{Numerical Example}\label{subsec:simulations}
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:
\begin{equation}\label{eq: dgp}
Y = g(W) + V'\beta + \sigma \Phi^{-1}(U),
\end{equation}
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{fig: estimands demand} and \ref{fig: estimands
elasticity} 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{ft:cov_values}). 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{eq: dgp}). 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{eq:
dgp}) 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{Def:hatsigma2n}).
The bandwidth for $\widehat J(u)$ is chosen using the Hall-Sheather
option of the \verb"quantreg" \verb"R" package
(Hall and Sheather \cite{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.