Fractional order statistic approximation for nonparametric conditional quantile inference

Supplemental appendix: Fractional order statistic approximation for nonparametric conditional quantile inference

\doublespacing

\iftoggle{SUPPLEMENTAL}

Introduction

Quantiles contain information about a distribution's shape. Complementing the mean, they capture heterogeneity, inequality, and other measures of economic interest. Nonparametric conditional quantile models further allow arbitrary heterogeneity across regressor values. This paper concerns nonparametric inference on quantiles and conditional quantiles. In particular, we characterize the high-order accuracy of both \citeposs{Hutson1999} $L$-statistic-based confidence intervals (CIs) and our new conditional quantile CIs.

Conditional quantiles appear across diverse topics because they are fundamental statistical objects. Such topics include wages Hogg1975,Chamberlain1994,Buchinsky1994, infant birthweight Abrevaya2001, demand for alcohol ManningEtAl1995, and Engel curves AlanEtAl2005,Deaton1997, which we examine in our empirical application.

We formally derive the coverage probability error (CPE) of the CIs from Hutson1999, as well as asymptotic power of the corresponding hypothesis tests. Hutson1999 had proposed CIs for quantiles using $L$-statistics (interpolating between order statistics) as endpoints and found they performed well, but formal proofs were lacking. Using the analytic $n^{-1}$ term we derive in the CPE, we provide a new calibration to achieve $O\bigl(n^{-3/2}[\log(n)]^3\bigr)$ CPE, analogous to the HoLee2005a analytic calibration of the CIs in BeranHall1993.

The theoretical results we develop contribute to the fractional order statistic literature and provide the basis for inference on other objects of interest explored in GoldmanKaplan2014b and Kaplan2014. In particular, Theorem (ref) tightly links the distributions of $L$-statistics from the observed and `ideal' (unobserved) fractional order statistic processes. Additionally, Lemma (ref) provides Dirichlet PDF and PDF derivative approximations.

High-order accuracy is important for small samples (e.g.,\ for experiments) as well as nonparametric analysis with small local sample sizes. For example, if $n=1024$ and there are five binary regressors, then the smallest local sample size cannot exceed $1024/2^5=32$.

For nonparametric conditional quantile inference, we apply the unconditional method to a local sample (similar to local constant kernel regression), smoothing over continuous covariates and also allowing discrete covariates. CPE is minimized by balancing the CPE of our unconditional method and the CPE from bias due to smoothing. We derive the optimal CPE and bandwidth rates, as well as a plug-in bandwidth when there is a single continuous covariate.

Our $L$-statistic method has theoretical and computational advantages over methods based on normality or an unsmoothed bootstrap. The theoretical bottleneck for our approach is the need to use a uniform kernel. Nonetheless, even if normality or bootstrap methods assume an infinitely differentiable conditional quantile function (and hypothetically fit an infinite-degree local polynomial), our CPE is still of smaller order with one or two continuous covariates. Our method also computes more quickly than existing methods (of reasonable accuracy), handling even more challenging tasks in 10--15 seconds instead of minutes.

Recent complementary work of FanLiu2016 also concerns a “direct method” of nonparametric inference on conditional quantiles. They use a limiting Gaussian process to derive first-order accuracy in a general setting, whereas we use the finite-sample Dirichlet process to achieve high-order accuracy in an iid setting. FanLiu2016 also provide uniform (over $X$) confidence bands. We suggest a confidence band from interpolating a growing number of joint CIs (as in HorowitzLee2012), although it will take additional work to rigorously justify. A different, ad hoc confidence band described in Section (ref) generally outperformed others in our simulations.

If applied to a local constant estimator with a uniform kernel and the same bandwidth, the FanLiu2016 approach is less accurate than ours due to the normal (instead of beta) reference distribution and integer (instead of interpolated) order statistics in their CI in equation (6). However, with other estimators like local polynomials or that in DonaldEtAl2012, the FanLiu2016 method is not necessarily less accurate. One limitation of our approach is that it cannot incorporate these other estimators, whereas Assumption GI(iii) in FanLiu2016 includes any estimator that weakly converges (over a range of quantiles) to a Gaussian process with a particular structure. We compare further in our simulations. One open question is whether using our beta reference and interpolation can improve accuracy for the general FanLiu2016 method beyond the local constant estimator with a uniform kernel; our Lemma (ref) shows this at least retains first-order accuracy.

The order statistic approach to quantile inference uses the idea of the probability integral transform, which dates back to R.\ A.\ Fisher1932, Karl Pearson1933, and Neyman1937. For continuous $X_i\stackrel{iid}{\sim}F(\cdot)$, $F(X_i)\stackrel{iid}{\sim}\textrm{Unif}(0,1)$. Each order statistic from such an iid uniform sample has a known beta distribution for any sample size $n$. We show that the $L$-statistic linearly interpolating consecutive order statistics also follows an approximate beta distribution, with only $O(n^{-1})$ error in CDF. Although $O(n^{-1})$ is an asymptotic claim, the CPE of the CI using the $L$-statistic endpoint is bounded between the CPEs of the CIs using the two order statistics comprising the $L$-statistic, where one such CPE is too small and one is too big, for any sample size. This is an advantage over methods more sensitive to asymptotic approximation error.

Many other approaches to one-sample quantile inference have been explored. With Edgeworth expansions, HallSheather1988 and Kaplan2015 obtain two-sided $O(n^{-2/3})$ CPE. With bootstrap, smoothing is necessary for high-order accuracy. This increases the computational burden and requires good bandwidth selection in practice.\footnote{For example, while achieving the impressive two-sided CPE of $O(n^{-3/2})$, PolanskySchucany1997 admit, “If this method is to be of any practical value, a better bandwidth estimation technique will certainly be required.”} See HoLee2005b for a review of bootstrap methods. Smoothed empirical likelihood ChenHall1993 also achieves nice theoretical properties, but with the same caveats.

Other order statistic-based CIs dating back to Thompson1936 are surveyed in DavidNagaraja2003. Most closely related to Hutson1999 is BeranHall1993. Like Hutson1999, BeranHall1993 linearly interpolate order statistics for CI endpoints, but with an interpolation weight based on the binomial distribution. Although their proofs use expansions of the Renyi1953 representation instead of fractional order statistic theory, their $n^{-1}$ CPE term is identical to that for Hutson1999 other than the different weight. Prior work Bickel1967,Shorack1972 has established asymptotic normality of $L$-statistics and convergence of the sample quantile process to a Gaussian limit process, but without such high-order accuracy.

The most apparent difference between the two-sided CIs of BeranHall1993 and Hutson1999 is that the former are symmetric in the order statistic index, whereas the latter are equal-tailed. This allows Hutson1999 to be computed further into the tails. Additionally, our framework can be extended to CIs for interquantile ranges and two-sample quantile differences GoldmanKaplan2014b, which has not been done in the R\'enyi representation framework.

For nonparametric conditional quantile inference, in addition to the aforementioned FanLiu2016 approach, Chaudhuri1991 derives the pointwise asymptotic normal distribution of a local polynomial estimator. QuYoon2015 propose modified local linear estimators of the conditional quantile process that converge weakly to a Gaussian process, and they suggest using a type of bias correction that strictly enlarges a CI to deal with the first-order effect of asymptotic bias when using the MSE-optimal bandwidth rate.

Section (ref) contains our theoretical results on fractional order statistic approximation, which are applied to unconditional quantile inference in Section (ref). Section (ref) concerns our new conditional quantile inference method. An empirical application and simulation results are in Sections (ref) and (ref), respectively. Proof sketches are collected in Appendix (ref), while the supplemental appendix contains full proofs. The supplemental appendix also contains details of the plug-in bandwidth calculations, as well as additional empirical and simulation results.

Notationally, $\phi(\cdot)$ and $\Phi(\cdot)$ are respectively the standard normal PDF and CDF, $\doteq$ should be read as “is equal to, up to smaller-order terms”, $\asymp$ as “has exact (asymptotic) rate/order of”, and $A_n=O(B_n)$ as usual. Acronyms used are those for cumulative distribution function (CDF), confidence interval (CI), coverage probability (CP), coverage probability error (CPE), and probability density function (PDF).

Fractional order statistic theory

In this section, we introduce notation and present our core theoretical results linking unobserved `ideal' fractional $L$-statistics with their observed counterparts.

Given an iid sample $\{X_i\}_{i=1}^n$ of draws from a continuous CDF denoted\footnote{$F$ will often be used with a random variable subscript to denote the CDF of that particular random variable. If no subscript is present, then $F(\cdot)$ refers to the CDF of $X$. Similarly for the PDF $f(\cdot)$.} $F(\cdot)$, interest is in $Q(p) \equiv F^{-1}(p)$ for some $p\in(0,1)$, where $Q(\cdot)$ is the quantile function. For $u\in(0,1)$, the sample $L$-statistic commonly associated with $Q(u)$ is

where $\lfloor\cdot\rfloor$ is the floor function, $\epsilon$ is the interpolation weight, and $X_{n:k}$ denotes the $k$th order statistic (i.e.,\ $k$th smallest sample value). While $Q(u)$ is latent and nonrandom, $\hat Q^L_X(u)$ is a random variable, and $\hat Q^L_X(\cdot)$ is a stochastic process, observed for arguments in $[1/(n+1),n/(n+1)]$.

Let $\Xi_n\equiv\mathopen{}\mathclose\bgroup\originalleft\{k/(n+1)\aftergroup\egroup\originalright\}_{k=1}^n$ denote the set of quantiles corresponding to the observed order statistics. If $u\in\Xi_n$, then no interpolation is necessary and $\hat Q^L_X(u)=X_{n:k}$. As detailed in Section (ref), application of the probability integral transform yields exact coverage probability of a CI endpoint $X_{n:k}$ for $Q(p)$: $ P\mathopen{}\mathclose\bgroup\originalleft(X_{n:k}<F^{-1}(p)\aftergroup\egroup\originalright) = P\mathopen{}\mathclose\bgroup\originalleft( U_{n:k}<p\aftergroup\egroup\originalright) $, where $U_{n:k}\equiv F(X_{n:k})\sim\beta(k,n+1-k)$ is equal in distribution to the $k$th order statistic from $U_i\stackrel{iid}{\sim}\textrm{Unif}(0,1)$, $i=1,\ldots,n$ Wilks1962. However, we also care about $u\notin\Xi_n$, in which case $k$ is fractional. To better handle such fractional order statistics, we will present a tight link between the marginal distributions of the stochastic process $\hat Q^L_X(\cdot)$ and those of the analogous `ideal' (I) process

where $\tilde Q^I_U(\cdot)$ is the ideal (I) uniform (U) fractional order “statistic” process. We use a tilde in $\tilde Q^I_X(\cdot)$ and $\tilde Q^I_U(\cdot)$ instead of the hat like in $\hat Q^L_X(\cdot)$ to emphasize that the former are unobserved (hence not true statistics), whereas the latter is computable from the sample data.

This $\tilde Q^I_U(\cdot)$ in (ref) is a Dirichlet process Ferguson1973,Stigler1977 on the unit interval with index measure $\nu\mathopen{}\mathclose\bgroup\originalleft([0,t]\aftergroup\egroup\originalright)=(n+1)t$. Its univariate marginals are

The marginal distribution of $\mathopen{}\mathclose\bgroup\originalleft(\tilde Q^I_U(u_1),\tilde Q^I_U(u_2)-\tilde Q^I_U(u_1),\ldots,\tilde Q^I_U(u_k)-\tilde Q^I_U(u_{k-1})\aftergroup\egroup\originalright)$ for $u_1<\cdots<u_k$ is Dirichlet with parameters $\mathopen{}\mathclose\bgroup\originalleft(u_1(n+1),(u_2-u_1)(n+1),\ldots,(u_k-u_{k-1})(n+1)\aftergroup\egroup\originalright)$.

For all $u\in\Xi_n$, $\tilde Q^I_X(u)$ coincides with $\hat Q^L_X(u)$; they differ only in their interpolation between these points. Proposition (ref) shows $\tilde Q^I_X(\cdot)$ and $\hat Q^L_X(\cdot)$ to be closely linked in probability.

Although Proposition (ref) motivates approximating the distribution of $\hat Q^L_X(u)$ by that of $\tilde Q^I_X(u)$, it is not relevant to high-order accuracy. In fact, its result is achieved by any interpolation between $X_{n:k}$ and $X_{n:k+1}$, not just $\hat Q^L_X(u)$; in contrast, the high-order accuracy we establish in Theorem (ref) is only possible with precise interpolations like $\hat Q^L_X(u)$.

Next, we consider marginal distributions of fixed dimension $J$. We also consider the Gaussian approximation to the sampling distribution of fractional order statistics. It is well known that the centered and scaled empirical process for standard uniform random variables converges to a Brownian bridge. For standard Brownian bridge process $B(\cdot)$, we index by $u\in(0,1)$ the additional stochastic processes

The vector $\tilde Q^I_U(\mathbf u)$ has an ordered Dirichlet distribution (i.e.,\ the spacings between consecutive $\tilde Q^I_U(u_j)$ follow a joint Dirichlet distribution), while $\tilde Q^B_U(\mathbf u)$ is multivariate Gaussian. Lemma (ref) in the appendix shows the close relationship between multivariate Dirichlet and Gaussian PDFs and PDF derivatives.

Theorem (ref) shows the close distributional link among linear combinations of ideal, interpolated, and Gaussian-approximated fractional order statistics. Specifically, for arbitrary weight vector $\boldsymbol{\psi} \in \mathbb{R}^J$, we (distributionally) approximate

Our assumptions for this section are now presented, followed by the main theoretical result. Assumption (ref) ensures that the first three derivatives of the quantile function are uniformly bounded in neighborhoods of the quantiles, $u_j$, which helps bound remainder terms in the proofs. We use bold for vectors and underline for matrices.

Quantile inference: unconditional

For inference on $Q(p)$, we continue to maintain (ref) and (ref). For $p\in(0,1)$ and confidence level $1-\alpha$, define $u^h(\alpha)$ and $u^l(\alpha)$ to solve

with $\tilde Q^I_U(u)\sim\beta\bigl((n+1)u,(n+1)(1-u)\bigr)$ from (ref), parallel to (7) and (8) in Hutson1999.

One-sided CI endpoints for $Q(p)$ are $\hat Q^L_X(u^h)$ or $\hat Q^L_X(u^l)$. Two-sided CIs replace $\alpha$ with $\alpha/2$ in (ref) and use both endpoints. This use of $\alpha/2$ yields the equal-tailed property; more generally, $t\alpha$ and $(1-t)\alpha$ can be used for $t\in(0,1)$.

Figure (ref) visualizes an example. The beta distribution's mean is $u^h$ (or $u^l$). Decreasing $u^h$ increases the probability mass in the shaded region below $u$, while increasing $u^h$ decreases the shaded region, and vice-versa for $u^l$. Solving (ref) is a simple numerical search problem.

Lemma (ref) shows the CI endpoint indices converge to $p$ at a $n^{-1/2}$ rate and may be approximated using quantiles of a normal distribution.

For the lower one-sided CI, using (ref), the $1-\alpha$ CI from Hutson1999 is

Coverage probability is

where $\phi(\cdot)$ is the standard normal PDF and the $n^{-1}$ term is non-negative. Similar to the HoLee2005a calibration, we can remove the analytic $n^{-1}$ term with the calibrated CI

which has CPE of order $O\mathopen{}\mathclose\bgroup\originalleft(n^{-3/2}[\log(n)]^3\aftergroup\egroup\originalright)$. We follow convention and define $\textrm{CPE}\equiv\textrm{CP}-(1-\alpha)$, where CP is the actual coverage probability and $1-\alpha$ the desired confidence level.

By parallel argument, \citeposs{Hutson1999} uncalibrated upper one-sided and two-sided CIs also have $O(n^{-1})$ CPE, or $O\mathopen{}\mathclose\bgroup\originalleft(n^{-3/2}[\log(n)]^3\aftergroup\egroup\originalright)$ with calibration. For the upper one-sided case, again using (ref), the $1-\alpha$ Hutson CI and our calibrated CI are respectively given by

and for equal-tailed two-sided CIs,

Without calibration, in all cases the $n^{-1}$ CPE term is non-negative (indicating over-coverage).

For relatively extreme quantiles $p$ (given $n$), the $L$-statistic method cannot be computed because the $(n+1)$th (or zeroth) order statistic is needed. In such cases, our code uses the Edgeworth expansion-based CI in Kaplan2015. Alternatively, if bounds on $X$ are known a priori, they may be used in place of these “missing” order statistics to generate conservative CIs. Regardless, as $n\to\infty$, the range of computable quantiles approaches $(0,1)$.

The hypothesis tests corresponding to all the foregoing CIs achieve optimal asymptotic power against local alternatives. The sample quantile is a semiparametric efficient estimator, so it suffices to show that power is asymptotically first-order equivalent to that of the test based on asymptotic normality. Theorem (ref) collects all of our results on coverage and power.

The equal-tailed property of our two-sided CIs is a type of median-unbiasedness. If $(\hat{L},\hat{H})$ is a CI for scalar $\theta$, then an equal-tailed CI is “unbiased” under loss function $L(\theta,\hat L,\hat H)=\max\{0,\theta-\hat H,\hat L-\theta\}$, as defined in (5) of Lehmann1951. This median-unbiased property may be desirable AndrewsGuggenberger2014, although it is different than the usual “unbiasedness” where a CI is the inversion of an unbiased test. More generally, in (ref), we could replace $u^l(\alpha/2)$ and $u^h(\alpha/2)$ by $u^l(t\alpha)$ and $u^h((1-t)\alpha)$ for $t\in[0,1]$. Different $t$ may achieve different optimal properties, which we leave to future work.

Quantile inference: conditional

Setup and bias

Let $Q_{Y|X}(u;x)$ be the conditional $u$-quantile function of scalar outcome $Y$ given conditioning vector $X\in\mathcal X\subset\mathbb R^d$, evaluated at $X=x$. The object of interest is $Q_{Y|X}(p;x_0)$, for $p\in(0,1)$ and interior point $x_0$. The sample $\{Y_i,X_i\}_{i=1}^n$ is drawn iid. Without loss of generality, let $x_0=0$.

If $X$ is discrete so that $P(X=0)>0$, we can take the subsample with $X_i=0$ and compute a CI from the corresponding $Y_i$ values, using the method in Section (ref). Even with dependence like strong mixing among the $X_i$, CPE is the same $O(n^{-1})$ from Theorem (ref) as long as the subsample's $Y_i$ are independent draws from the same $Q_{Y|X}(\cdot;0)$ and $N_n\stackrel{a.s.}{\asymp}n$.

If $X$ is continuous, then $P(X_i=0)=0$, so observations with $X_i\ne0$ must be included. If $X$ contains mixed continuous and discrete components, then we can apply our method for continuous $X$ to each subsample corresponding to each unique value of the discrete subvector of $X$. The asymptotic rates are unaffected by the presence of discrete variables (although the finite-sample consequences may deserve more attention), so we focus on the case where all components of $X$ are continuous.

We now present definitions and assumptions, continuing the normalization $x_0=0$.

Definition (ref) refers to a window whose size depends on $h$: $C_h=[-h,h]$ if $d=1$, or more generally a hypercube as in Chaudhuri1991: letting $\|\cdot\|_\infty$ denote the $L_\infty$-norm,

Given fixed values of $n$ and $h$, Assumption (ref) implies that the $Y_i$ in the local sample are independent and identically distributed,\footnote{This may be the case asymptotically even with substantial dependence, although we do not explore this point. For example, PolonikYao2002 write, “Only the observations with $X_t$ in a small neighbourhood of $x$ are effectively used\ldots [which] are not necessarily close with each other in the time space. Indeed, they could be regarded as asymptotically independent under appropriate conditions such as strong mixing\ldots.”} which is needed to apply Theorem (ref). However, they do not have the quantile function of interest, $Q_{Y|X}(\cdot;0)$, but rather the biased $Q_{Y|X}(\cdot;C_h)$. This is like drawing a global (any $X_i$) iid sample of wages, $Y_i$, and restricting it to observations in Japan ($X\in C_h$) when our interest is only in Tokyo ($X=0$): our restricted $Y_i$ constitute an iid sample from Japan, but the $p$-quantile wage in Japan may differ from that in Tokyo. Assumptions (ref)--(ref)(i) and (ref) are necessary for the calculation of this bias, $Q_{Y|X}(p;C_h)-Q_{Y|X}(p;0)$, in Lemma (ref). Assumptions (ref)(ii) and (ref) (and (ref)) ensure $N_n\stackrel{a.s.}{\to}\infty$. Assumptions (ref) and (ref) are conditional versions of Assumptions (ref)(i) and (ref)(ii), respectively. Their uniformity ensures uniformity of the remainder term in Theorem (ref), accounting for the fact that the local sample's distribution, $F_{Y|X}(\cdot;C_h)$, changes with $n$ (through $h$ and $C_h$).

From (ref)(i), asymptotically $C_h$ is entirely contained within the neighborhoods implicit in (ref), (ref), and (ref). This in turn allows us to examine only a local neighborhood around $p$ (e.g.,\ as in (ref)) since the CI endpoints converge to the true value at a $N_n^{-1/2}$ rate.

The $X_i$ being iid helps guarantee that $N_n$ is almost surely of order $nh^d$. The $h^d$ comes from the volume of $C_h$. Larger $h$ lowers CPE via $N_n$ but raises CPE via bias. This tradeoff determines the optimal rate at which $h\to0$ as $n\to\infty$. Using Theorem (ref) and additional results on CPE from bias below, we determine the optimal value of $h$.

The bias characterized in Lemma (ref) is the difference between these two population conditional quantiles.

Equation (ref) is the same as in BhattacharyaGangopadhyay1990, who derive it using different arguments.

Optimal CPE order

The CPE-optimal bandwidth minimizes the sum of the two dominant high-order CPE terms. It must be small enough to control the $O(h^b+N_nh^{2b})$ (two-sided) CPE from bias, but large enough to control the $O(N_n^{-1})$ CPE from applying the unconditional $L$-statistic method. The following theorem summarizes optimal bandwidth and CPE results.

As detailed in the supplemental appendix, Theorem (ref) implies that for the most common values of dimension $d$ and most plausible values of smoothness $s_Q$, even our uncalibrated method is more accurate than inference based on asymptotic normality with a local polynomial estimator. The same comparisons apply to basic bootstraps, which claim no refinement over asymptotic normality; in this (quantile) case, even Studentization does not improve theoretical CPE without the added complications of smoothed or $m$-out-of-$n$ bootstraps.

The only opportunity for normality to yield smaller CPE is to greatly reduce bias by using a very large local polynomial if $s_Q$ is large; our approach implicitly uses a uniform kernel, so bias reduction beyond $O(h^2)$ is impossible. Nonetheless, our method has smaller CPE when $d=1$ or $d=2$ even if $s_Q=\infty$, and in other cases the necessary local polynomial degree may be prohibitively large given common sample sizes.

Figure (ref) (left panel) shows that if $s_Q=2$ and $s_X=1$, then the optimal CPE from asymptotic normality is always larger (worse) than our method's CPE. As shown in the supplement, CPE with normality is nearly $O\bigl(n^{-2/(4+2d)}\bigr)$. With $d=1$, this is $O\bigl(n^{-1/3}\bigr)$, much larger than our two-sided $O\bigl(n^{-2/3}\bigr)$. With $d=2$, $O\bigl(n^{-1/4}\bigr)$ is larger than our $O\bigl(n^{-1/2}\bigr)$. It remains larger for all $d$ since the bias is the same for both methods while the unconditional $L$-statistic inference is more accurate than normality.

Figure (ref) (right panel) shows the required amount of smoothness and local polynomial degree for asymptotic normality to match our method's CPE. For the most common cases of $d=1$ and $d=2$, two-sided CPE with normality is larger even with infinite smoothness and a hypothetical infinite-degree polynomial. With $d=3$, to match our CPE, normality needs $s_Q\ge12$ and a local polynomial of degree $k_Q\ge11$. Since interaction terms are required, an $11$th-degree polynomial has $\sum_{T=d-1}^{k_Q+d-1} \binom{T}{d-1} = 364$ terms, which requires a large $N_n$ (and yet larger $n$). As $d\to\infty$, the required number of terms in the local polynomial only grows larger and may be prohibitive in realistic finite samples.

Plug-in bandwidth

We propose a feasible bandwidth value with the CPE-optimal rate. To avoid recursive dependence on $\epsilon$ (the interpolation weight), we fix its value. This does not achieve the theoretical optimum, but it remains close even in small samples and seems to work well in practice. The CPE-optimal bandwidth value derivation is shown for $d=1$ in the supplemental appendix; a plug-in version is implemented in our code. For reference, the plug-in bandwidth expressions are collected here. The $\alpha$-quantile of $N(0,1)$ is again denoted $z_\alpha$. We let $\hat B_h$ denote the estimator of bias term $B_h$; $\hat f_X$ the estimator of $f_X(x_0)$; $\hat f_X'$ the estimator of $f_X'(x_0)$; $\hat F_{Y|X}^{(0,1)}$ the estimator of $F_{Y|X}^{(0,1)}(\xi_p;x_0)$; and $\hat F_{Y|X}^{(0,2)}$ the estimator of $F_{Y|X}^{(0,2)}(\xi_p;x_0)$, with notation from Lemma (ref).

When $d=1$, the following are our CPE-optimal plug-in bandwidths.

While we suggest the CPE-optimal bandwidths for moderate $n$, we suggest shifting toward a larger bandwidth as $n\to\infty$. Once CPE is small over a range of bandwidths, a larger bandwidth in that range is preferable since it yields shorter CIs. As an initial suggestion, we use a coefficient of $\max\{1,n/1000\}^{5/60}$ that keeps the CPE-optimal bandwidth for $n\le1000$ and then moves toward a $n^{-1/20}$ under-smoothing of the MSE-optimal bandwidth rate, as in FanLiu2016.

Empirical application

We present an application of our $L$-statistic inference to Engel1857 curves. Code is available from the latter author's website, and the data are publicly available.

BanksEtAl1997 argue that a linear Engel curve is sufficient for certain categories of expenditure, while adding a quadratic term suffices for others. Their Figure 1 shows nonparametrically estimated mean Engel curves (budget share $W$ against log total expenditure $\ln(X)$) with $95\%$ pointwise CIs at the deciles of the total expenditure distribution, using a subsample of 1980--1982 U.K.\ Family Expenditure Survey (FES) data.

We present a similar examination, but for quantile Engel curves in the 2001--2012 U.K.\ Living Costs and Food Surveys UKLCFS2012, which is a successor to the FES. We examine the same four categories as in the original analysis: food; fuel, light, and power (“fuel”); clothing and footwear (“clothing”); and alcohol. We use the subsample of households with one adult male and one adult female (and possibly children) living in London or the South East, leaving $8{,}528$ observations. Expenditure amounts are adjusted to 2012 nominal values using annual CPI data.{\interfootnotelinepenalty10000\footnote{\url{http://www.ons.gov.uk/ons/datasets-and-tables/data-selector.html?cdid=D7BT&dataset=mm23&table-id=1.1}}}

Table (ref) shows unconditional $L$-statistic CIs for various quantiles of the budget share distributions for the four expenditure categories. (Due to the large sample size, calibrated CIs are identical at the precision shown.) These capture some population features, but the conditional quantiles are of more interest.

Figure (ref) is comparable to Figure 1 of BanksEtAl1997 but with $90\%$ joint (over the nine expenditure levels) CIs instead of $95\%$ pointwise CIs, alongside quadratic quantile regression estimates. (To get joint CIs, we simply use the Bonferroni adjustment and compute $1-\alpha/9$ pointwise CIs.) Joint CIs are more intuitive for assessing the shape of a function since they jointly cover all corresponding points on the true curve with $90\%$ probability, rather than any given single point. The CIs are interpolated only for visual convenience. Although some of the joint CI shapes do not look quadratic at first glance, the only cases where the quadratic fit lies outside one of the intervals are for alcohol at the conditional median and clothing at the conditional upper quartile, and neither is a radical departure. With a $90\%$ confidence level and 12 confidence sets, we would not be surprised if one or two did not cover the true quantile Engel curve completely. Importantly, the CIs are relatively precise, too; the linear fit is rejected in 8 of 12 cases. Altogether, this evidence suggests that the benefits of a quadratic (but not linear) approximation may outweigh the cost of approximation error.

The supplemental appendix includes a similar figure but with a nonparametric (instead of quadratic) conditional quantile estimate along with joint CIs from FanLiu2016.

Simulation study

Code for our methods and simulations is available on the latter author's website.

Unconditional simulations

We compare two-sided unconditional CIs from the following methods: “L-stat” from Section (ref), originally in Hutson1999; “BH” from BeranHall1993; “Norm” using the sample quantile's asymptotic normality and kernel-estimated variance; “K15” from Kaplan2015; and “BStsym,” a symmetric Studentized bootstrap (99 draws) with bootstrapped variance (100 draws).\footnote{Other bootstraps were consistently worse in terms of coverage: (asymmetric) Studentized bootstrap, and percentile bootstrap with and without symmetry.}

Overall, L-stat and BH have the most accurate coverage probability (CP), avoiding under-coverage while maintaining shorter length than other methods achieving at least $95\%$ CP. Near the median, L-stat and BH are nearly identical. Away from the median, L-stat is closer to equal-tailed and often shorter than BH. Farther into the tails, L-stat can be computed where BH cannot.

Table (ref) shows nearly exact CP for both L-stat and BH when $n=25$ and $p=0.5$. “Norm” can be slightly shorter, but it under-covers. The bootstrap has only slight under-coverage, and K15 none, but their CIs are longer than L-stat's. Additional results are in the supplemental appendix, but the qualitative points are the same.

Table (ref) shows a case in the lower tail with $n=99$ where BH cannot be computed (because it needs the zeroth order statistic). Even then, L-stat's CP remains almost exact, and it is closest to equal-tailed. “Norm” under-covers for two $F$ (severely for Cauchy) and is almost twice as long as L-stat for the third. BStsym has less under-coverage, and K15 none, but both are generally longer than L-stat. Again, additional results are in the supplemental appendix, with similar patterns.

The supplemental appendix contains additional simulation results for $p\ne0.5$ but where BH is still computable. L-stat and BH both attain $95\%$ CP, but L-stat is much closer to equal-tailed and is shorter. The supplemental appendix also has results illustrating the effect of calibration.

Table (ref) isolates the effects of using the beta distribution rather than the normal approximation, as well as the effects of interpolation. Method “Normal” uses the normal approximation to determine $u^h$ and $u^l$ but still interpolates, while “Norm/floor” uses the normal approximation with no interpolation as in equations (5) and (6) of FanLiu2016.

Table (ref) shows several advantages of L-stat. First, for $p=0.15$, Normal and Norm/floor cannot even be computed (hence “NA”) because they require the zeroth order statistic, which does not exist, whereas L-stat is computable and has nearly exact CP ($0.905$). Second, with $p=0.25$ and $p=0.5$, the normal approximation (Normal) makes the CI needlessly longer than L-stat's CI. Third, additionally not interpolating (Norm/floor) makes the CI even longer for $p=0.25$ but leads to under-coverage for $p=0.5$. Fourth, whereas the L-stat CIs are almost exactly equal-tailed, the normal-based CIs are far from equal-tailed at $p=0.25$, where Norm/floor is essentially a one-sided CI.

Conditional simulations

For conditional quantile inference, we compare our $L$-statistic method (“L-stat”) with a variety of others. Implementation details may be seen in the supplemental appendix and available code. The first other method (“rqss”) is from the popular quantreg package in R R.quantreg. The second (“boot”) is a local cubic method following Chaudhuri1991 but with bootstrapped standard errors; the bandwidth is L-stat's multiplied by $n^{1/12}$ to get the local cubic CPE-optimal rate. The third (“QYg”) uses the asymptotic normality of a local linear estimator with a Gaussian kernel, using results and ideas from QuYoon2015, although they are more concerned with uniform (over quantiles) inference; they suggest using the MSE-optimal bandwidth (Corollary 1) and a particular type of bias correction (Remark 7). The fourth (“FLb”) is from Section 3.1 in FanLiu2016, based on a symmetrized $k$-NN estimator using a bisquare kernel; we use the code from their simulations.\footnote{Graciously provided to us. The code differs somewhat from the description in their text, most notably by an additional factor of $0.4$ in the bandwidth.} Interestingly, although in principle they are just slightly undersmoothing the MSE-optimal bandwidth, their bandwidth is very close to the CPE-optimal bandwidth for the sample sizes considered.

We now write $x_0$ as the point of interest, instead of $x_0=0$; we also take $d=1$, $b=2$, and focus on two-sided inference, both pointwise (single $x_0$) and joint (over multiple $x_0$). Joint CIs for all methods are computed using the Bonferroni approach. Uniform bands are also examined, with L-stat, QYg, and boot relying on the adjusted critical value from the Hotelling1939 tube computations in plot.rqss. Each simulation has $1{,}000$ replications unless otherwise noted.

Figure (ref) uses Model 1 from FanLiu2016: $Y_i=2.5+\sin(2X_i)+2\exp\bigl(-16X_i^2\bigr)+0.5\epsilon_i$, $X_i\stackrel{iid}{\sim}N(0,1)$, $\epsilon_i\stackrel{iid}{\sim}N(0,1)$, $X_i\protect\mathpalette{\protect\independenT}{\perp}\epsilon_i$, $n=500$, $p=0.5$. The “Direct” method in their Table 1 is our FLb. All methods have good pointwise CP (top left). L-stat has the best pointwise power (top right).

Figure (ref) (bottom left) shows power curves of the hypothesis tests corresponding to the joint (over $x_0\in\{0,0.75,1.5\}$) CIs, varying $H_0$ while maintaining the same DGP. The deviations of $Q_{Y|X}(p;x_0)$ shown on the horizontal axis are the same at each $x_0$; zero deviation implies $H_0$ is true, in which case the rejection probability is the type I error rate. All methods have good type I error rates: L-stat's is 6.2%, and other methods' are below the nominal 5%. L-stat has significantly better power, an advantage of 20--40% at the larger deviations. The bottom right graph in Figure (ref) is similar, but based on uniform confidence bands evaluated at $231$ different $x_0$. Only L-stat has nearly exact type I error rate and good power.

Next, we use the simulation setup of the rqss vignette in R.quantreg, which in turn came in part from RuppertEtAl2003. Here, $n=400$, $p=0.5$, $d=1$, $\alpha=0.05$, and

where the $U_i$ are iid $N(0,1)$, $t_3$, Cauchy, or centered $\chi^2_3$, and $\sigma(X)=0.2$ or $\sigma(X)=0.2(1+X)$. The conditional median function is graphed in the supplemental appendix. Although the function as a whole is not a common shape in economics (with multiple local maxima and minima), it provides insight into different types of functions at different points. For pointwise and joint CIs, we consider 47 equispaced points, $x_0=0.04,0.06,\ldots,0.96$; uniform confidence bands are evaluated at 231 equispaced values of $x_0$.

Figure (ref)'s first two columns show that across all eight DGPs (four error distributions, homoskedastic or heteroskedastic), L-stat has consistently accurate pointwise CP. At the most challenging points (smallest $x_0$), L-stat can under-cover by around five percentage points. Otherwise, CP is near $1-\alpha$ for all $x_0$ in all DGPs.

In contrast, with the exception of boot, the other methods can have significant under-coverage. As seen in the first two columns of Figure (ref), rqss has under-coverage (as low as 50--60% CP) for $x_0$ closer to zero. QYg has under-coverage with the $\chi^2_3$ and (especially) Cauchy. FLb has good CP except with the Cauchy, where CP can dip below 70%.

Figure (ref)'s third column shows the joint power curves. The horizontal axis of the graphs indicates the deviation of $H_0$ from the true values. For example, letting $\xi_{p,j}$ be the true conditional quantiles at the $j=1,\ldots,47$ values of $x_0$ (say, $x_j$), $-0.1$ deviation refers to $H_0:\{Q_{Y|X}(p;x_j)=\xi_{p,j}-0.1\textrm{ for }j=1,\ldots,47\}$ (which is false), and zero deviation means $H_0$ is true. Our method's type I error rate is close to $\alpha$ under all four $U_i$ distributions (5.7%, 5.8%, 7.3%, 6.3%). In contrast, other methods show size distortion under Cauchy and/or $\chi^2_3$ $U_i$; among them, boot is closest but still has $10.3\%$ type I error rate with the $\chi^2_3$. Next-best is rqss; size distortion for FLb and QYg is more serious. L-stat also has the steepest joint power curves among all methods. Beyond steepness, they are also the most robust to the underlying distribution. L-stat's type I error rate is near 5% for all four distributions. In contrast, boot ranges from only 1.2% for the Cauchy, leading to worse power, up to 10.3% for the $\chi^2_3$.

The supplemental appendix shows a comparison of hypothesis tests based on uniform confidence bands. The results are similar to the joint power curves, but with slightly higher rejection rates all around.

Figure (ref) shows pointwise power. Specifically, for a given $x_0$, this is the proportion of simulation draws in which $Q_{Y|X}(p;x_0)-0.1$ is excluded from the CI, averaged with the corresponding proportion for $Q_{Y|X}(p;x_0)+0.1$. L-stat generally has the best power among methods with correct CP (per first column of Figure (ref)).

The supplemental appendix contains results for $p=0.25$, where L-stat continues to perform well. One additional advantage is that L-stat's joint test is nearly unbiased, whereas the other joint tests are all biased.

The supplemental appendix also shows the computational advantage of our method. For example, with $n=10^5$ and $100$ different $x_0$, L-stat takes only $10$ seconds, whereas the local cubic bootstrap takes $141$ seconds; rqss is even slower.

Overall, the simulation results show the new L-stat method to be fast and accurate. Besides L-stat, the only method to avoid serious under-coverage is the local cubic with bootstrapped standard errors, perhaps due to its reliance on our newly proposed CPE-optimal bandwidth. However, L-stat consistently has better power, greater robustness across different conditional distributions, and less bias of its joint hypothesis tests.

Conclusion

We derive a uniform $O(n^{-1})$ difference between the linearly interpolated and ideal fractional order statistic distributions. We generalize this to $L$-statistics to help justify quantile inference procedures. In particular, this translates to $O(n^{-1})$ CPE for the quantile CIs proposed by Hutson1999, which we improve to $O\mathopen{}\mathclose\bgroup\originalleft(n^{-3/2}[\log(n)]^3\aftergroup\egroup\originalright)$ via calibration. We extend these results to a nonparametric conditional quantile model, with both theoretical and Monte Carlo success. The derivation of an optimal bandwidth value (not just rate) and a fast approximation thereof are important practical advantages.

Our results can be extended to other objects of interest, such as interquantile ranges and two-sample quantile differences GoldmanKaplan2014b, quantile marginal effects Kaplan2014, and entire distributions GoldmanKaplan2015c.

In ongoing work, we consider the connection with Bayesian bootstrap quantile inference, which may be a way to “relax” the iid assumption. Other future work may improve finite-sample performance, e.g.,\ by smoothing over discrete covariates LiRacine2007.

\singlespacing \iftoggle{SUPPLEMENTAL}{

}

\doublespacing \onehalfspacing \singlespacing

\iftoggle{SUPPLEMENTAL}{\setcounter{page}{2}}