Fractional order statistic approximation for
nonparametric conditional quantile inference
arXiv:1609.09035v1 [] 28 Sep 2016

Matt Goldman∗

David M. Kaplan∗

September 29, 2016

Abstract
Using and extending fractional order statistic theory, we characterize the O(n−1 )
coverage probability error of the previously proposed confidence intervals for population quantiles using L-statistics as endpoints in Hutson (1999). We derive an analytic
expression for the n−1 term,
which may be used to calibrate the nominal coverage level

−3/2
3
to get O n
[log(n)] coverage error. Asymptotic power is shown to be optimal.
Using kernel smoothing, we propose a related method for nonparametric inference on
conditional quantiles. This new method compares favorably with asymptotic normality and bootstrap methods in theory and in simulations. Code is provided for both
unconditional and conditional inference.
JEL classification: C21
Keywords: Dirichlet, high-order accuracy, inference-optimal bandwidth, kernel
smoothing.

c 2016 by the authors. This manuscript version is made available under the
CC-BY-NC-ND 4.0 license: http://creativecommons.org/licenses/by-nc-nd/4.0/

∗

Goldman: Microsoft, mattgold@microsoft.com. Kaplan (corresponding author): Department of Economics, University of Missouri, kaplandm@missouri.edu. We thank the co-editor (Oliver Linton), associate
editor, referees, and Yixiao Sun for helpful comments and references, and Patrik Guggenberger and Andres
Santos for feedback improving the clarity of presentation. Thanks also to Brendan Beare, Karen Messer, and
active audience members at seminars and conferences. Thanks to Ruixuan Liu for providing and discussing
code from Fan and Liu (2016). This paper was previously circulated as parts of “IDEAL quantile inference
via interpolated duals of exact analytic L-statistics” and “IDEAL inference on conditional quantiles.”

1

1

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 Hutson’s (1999) 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 (Buchinsky, 1994; Chamberlain, 1994; Hogg, 1975),
infant birthweight (Abrevaya, 2001), demand for alcohol (Manning, Blumberg, and Moulton,
1995), and Engel curves (Alan, Crossley, Grootendorst, and Veall, 2005; Deaton, 1997, pp.
81–82), which we examine in our empirical application.
We formally derive the coverage probability error (CPE) of the CIs from Hutson (1999), as
well as asymptotic power of the corresponding hypothesis tests. Hutson (1999) 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 n−3/2 [log(n)]3 CPE, analogous
to the Ho and Lee (2005a) analytic calibration of the CIs in Beran and Hall (1993).
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 Goldman and
Kaplan (2016b) and Kaplan (2014). In particular, Theorem 2 tightly links the distributions
of L-statistics from the observed and ‘ideal’ (unobserved) fractional order statistic processes.
Additionally, Lemma 7 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/25 = 32.
For nonparametric conditional quantile inference, we apply the unconditional method
2

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 Fan and Liu (2016) 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. Fan and Liu (2016) also provide
uniform (over X) confidence bands. We suggest a confidence band from interpolating a
growing number of joint CIs (as in Horowitz and Lee (2012)), although it will take additional
work to rigorously justify. A different, ad hoc confidence band described in Section 6 generally
outperformed others in our simulations.
If applied to a local constant estimator with a uniform kernel and the same bandwidth,
the Fan and Liu (2016) 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 Donald,
Hsu, and Barrett (2012), the Fan and Liu (2016) method is not necessarily less accurate.
One limitation of our approach is that it cannot incorporate these other estimators, whereas
Assumption GI(iii) in Fan and Liu (2016) includes any estimator that weakly converges (over

3

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 Fan and Liu (2016) method beyond the local constant
estimator with a uniform kernel; our Lemma 3 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. Fisher (1932), Karl Pearson (1933), and Neyman
iid

iid

(1937). For continuous Xi ∼ F (·), F (Xi ) ∼ 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, Hall and Sheather (1988) and Kaplan (2015) 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.1 See Ho and Lee
(2005b, §1) for a review of bootstrap methods. Smoothed empirical likelihood (Chen and
Hall, 1993) also achieves nice theoretical properties, but with the same caveats.
Other order statistic-based CIs dating back to Thompson (1936) are surveyed in David
and Nagaraja (2003, §7.1). Most closely related to Hutson (1999) is Beran and Hall (1993).
Like Hutson (1999), Beran and Hall (1993) linearly interpolate order statistics for CI endpoints, but with an interpolation weight based on the binomial distribution. Although their
proofs use expansions of the Rényi (1953) representation instead of fractional order statistic
For example, while achieving the impressive two-sided CPE of O(n−3/2 ), Polansky and Schucany (1997,
p. 833) admit, “If this method is to be of any practical value, a better bandwidth estimation technique will
certainly be required.”
1

4

theory, their n−1 CPE term is identical to that for Hutson (1999) other than the different
weight. Prior work (e.g., Bickel, 1967; Shorack, 1972) 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 Beran and Hall (1993) and
Hutson (1999) is that the former are symmetric in the order statistic index, whereas the
latter are equal-tailed. This allows Hutson (1999) to be computed further into the tails.
Additionally, our framework can be extended to CIs for interquantile ranges and two-sample
quantile differences (Goldman and Kaplan, 2016b), which has not been done in the Rényi
representation framework.
For nonparametric conditional quantile inference, in addition to the aforementioned Fan
and Liu (2016) approach, Chaudhuri (1991) derives the pointwise asymptotic normal distribution of a local polynomial estimator. Qu and Yoon (2015) 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 2 contains our theoretical results on fractional order statistic approximation,
which are applied to unconditional quantile inference in Section 3. Section 4 concerns our
new conditional quantile inference method. An empirical application and simulation results
are in Sections 5 and 6, respectively. Proof sketches are collected in Appendix A, 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, φ(·) and Φ(·) are respectively the standard normal PDF and CDF, = should
be read as “is equal to, up to smaller-order terms”,  as “has exact (asymptotic) rate/order
of”, and An = O(Bn ) as usual. Acronyms used are those for cumulative distribution function
(CDF), confidence interval (CI), coverage probability (CP), coverage probability error (CPE),

5

and probability density function (PDF).

2

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 {Xi }ni=1 of draws from a continuous CDF denoted2 F (·), interest is
in Q(p) ≡ F −1 (p) for some p ∈ (0, 1), where Q(·) is the quantile function. For u ∈ (0, 1), the
sample L-statistic commonly associated with Q(u) is
Q̂LX (u) ≡ (1 − )Xn:k + Xn:k+1 ,

k ≡ bu(n + 1)c,

 ≡ u(n + 1) − k,

(1)

where b·c is the floor function,  is the interpolation weight, and Xn:k denotes the kth order
statistic (i.e., kth smallest sample value). While Q(u) is latent and nonrandom, Q̂LX (u) is
a random variable, and Q̂LX (·) is a stochastic process, observed for arguments in [1/(n +
1), n/(n + 1)].
Let Ξn ≡ {k/(n + 1)}nk=1 denote the set of quantiles corresponding to the observed order
statistics. If u ∈ Ξn , then no interpolation is necessary and Q̂LX (u) = Xn:k . As detailed in
Section 3, application of the probability integral transform yields exact coverage probability
of a CI endpoint Xn:k for Q(p): P (Xn:k < F −1 (p)) = P (Un:k < p), where Un:k ≡ F (Xn:k ) ∼
iid

β(k, n + 1 − k) is equal in distribution to the kth order statistic from Ui ∼ Unif(0, 1),
i = 1, . . . , n (Wilks, 1962, 8.7.4). However, we also care about u ∈
/ Ξ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 Q̂LX (·) and those of the analogous
‘ideal’ (I) process
Q̃IX (·)

≡F

−1

2





Q̃IU (·)

,

(2)

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 (·) refers to the CDF of X. Similarly for the PDF f (·).

6

where Q̃IU (·) is the ideal (I) uniform (U) fractional order “statistic” process. We use a tilde in
Q̃IX (·) and Q̃IU (·) instead of the hat like in Q̂LX (·) to emphasize that the former are unobserved
(hence not true statistics), whereas the latter is computable from the sample data.
This Q̃IU (·) in (2) is a Dirichlet process (Ferguson, 1973; Stigler, 1977) on the unit interval
with index measure ν([0, t]) = (n + 1)t. Its univariate marginals are

Q̃IU (u) = Un:(n+1)u ∼ β (n + 1)u, (n + 1)(1 − u) .

(3)



The marginal distribution of Q̃IU (u1 ), Q̃IU (u2 ) − Q̃IU (u1 ), . . . , Q̃IU (uk ) − Q̃IU (uk−1 ) for u1 <
· · · < uk is Dirichlet with parameters (u1 (n + 1), (u2 − u1 )(n + 1), . . . , (uk − uk−1 )(n + 1)).
For all u ∈ Ξn , Q̃IX (u) coincides with Q̂LX (u); they differ only in their interpolation
between these points. Proposition 1 shows Q̃IX (·) and Q̂LX (·) to be closely linked in probability.
Proposition 1. For any fixed δ > 0 and m > 0, define U δ ≡ {u ∈ (0, 1) | ∀ t ∈ (u − m, u +
n
1
, n+1
]; then, sup Q̃IX (u) − Q̂LX (u) = Op (n−1 log(n)).
m), f (F −1 (t)) ≥ δ} and Unδ ≡ U δ ∩ [ n+1
u∈Unδ

Although Proposition 1 motivates approximating the distribution of Q̂LX (u) by that of
Q̃IX (u), it is not relevant to high-order accuracy. In fact, its result is achieved by any
interpolation between Xn:k and Xn:k+1 , not just Q̂LX (u); in contrast, the high-order accuracy
we establish in Theorem 4 is only possible with precise interpolations like Q̂LX (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(·), we index by
u ∈ (0, 1) the additional stochastic processes
−1/2
Q̃B
B(u)
U (u) ≡ u + n

and


−1
Q̃B
Q̃B
X (u) ≡ F
U (u) .

The vector Q̃IU (u) has an ordered Dirichlet distribution (i.e., the spacings between consecutive Q̃IU (uj ) follow a joint Dirichlet distribution), while Q̃B
U (u) is multivariate Gaussian.
Lemma 7 in the appendix shows the close relationship between multivariate Dirichlet and
7

Gaussian PDFs and PDF derivatives.
Theorem 2 shows the close distributional link among linear combinations of ideal, interpolated, and Gaussian-approximated fractional order statistics. Specifically, for arbitrary
weight vector ψ ∈ RJ , we (distributionally) approximate
L

L ≡

J
X

ψj Q̂LX (uj )

J
X

I

by L ≡

j=1

ψj Q̃IX (uj ),

(4)

j=1
J
X

or alternatively by LB ≡

ψj Q̃B
X (uj ).

j=1

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

Assumption A1. Sampling is iid: Xi ∼ F , i = 1, . . . , n.
Assumption A2. For each quantile uj , the PDF f (·) (corresponding to CDF F (·) in A1)
satisfies (i) f (F −1 (uj )) > 0; (ii) f 00 (·) is continuous in some neighborhood of F −1 (uj ), i.e.,
f ∈ C 2 (Uδ (F −1 (uj ))) with Uδ (x) denoting some δ-neighborhood of point x ∈ R.
Theorem 2. Define V as the J × J matrix with row i, column j entries V i,j = min{ui , uj } −
ui uj . Let A be the J × J matrix with main diagonal entries Aj,j = f (F −1 (uj )) and zeros
elsewhere, and let

Vψ ≡ ψ 0 A−1 V A−1 ψ,

X0 ≡

J
X

ψj F −1 (uj ).

j=1

Let Assumption A1 hold, and let A2 hold at ū. Given the definitions in (1), (2), and (4),
the following results hold uniformly over u = ū + o(1).
(i) For a given constant K,


L

P L < X0 + n

−1/2


K



I
−1/2
− P L < X0 + n
K
8

" J
!#

K exp{−K 2 /(2Vψ )} X ψj2 j (1 − j )
−1
−3/2
3
q
,
n
+
O
n
[log(n)]
=
f [F −1 (uj )]2
2πVψ3
j=1
where the remainder is uniform over all K.
(ii) Uniformly over K,
 



L
−1/2
I
−1/2
sup P L < X0 + n
K − P L < X0 + n
K
K∈R
" J
!#

e−1/2 X ψj2 j (1 − j )
=q
n−1 + O n−3/2 [log(n)]3 ,
2
−1
2πVψ2 j=1 f [F (uj )]





L
−1/2
B
−1/2
sup P L < X0 + n
K − P L < X0 + n
K = O n−1/2 [log(n)]3 .
K∈R

3

Quantile inference: unconditional

For inference on Q(p), we continue to maintain A1 and A2. For p ∈ (0, 1) and confidence
level 1 − α, define uh (α) and ul (α) to solve



α = P Q̃IU uh (α) < p ,




α = P Q̃IU ul (α) > p ,

(5)


with Q̃IU (u) ∼ β (n + 1)u, (n + 1)(1 − u) from (3), parallel to (7) and (8) in Hutson (1999).
One-sided CI endpoints for Q(p) are Q̂LX (uh ) or Q̂LX (ul ). Two-sided CIs replace α with
α/2 in (5) and use both endpoints. This use of α/2 yields the equal-tailed property; more
generally, tα and (1 − t)α can be used for t ∈ (0, 1).
Figure 1 visualizes an example. The beta distribution’s mean is uh (or ul ). Decreasing
uh increases the probability mass in the shaded region below u, while increasing uh decreases
the shaded region, and vice-versa for ul . Solving (5) is a simple numerical search problem.
Lemma 3 shows the CI endpoint indices converge to p at a n−1/2 rate and may be
approximated using quantiles of a normal distribution.
Lemma 3. Let z1−α denote the (1 − α)-quantile of a standard normal distribution, z1−α ≡

9

5

Determination of Upper Endpoint
n=11, u=0.65, α=0.1

5

Determination of Lower Endpoint
n=11, u=0.65, α=0.1

4
3

Density

2
1

2

3

ul = 0.42

1

Density

4

uh = 0.84

α = 10%

0

0

α = 10%

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

u

0.4

0.6

0.8

1.0

u

Figure 1: Example of one-sided CI endpoint determination,
n = 11, p = 0.65, α = 0.1.

l
I
l
Left: u makes the shaded region’s area P Q̃U (u ) > p = α. Right: similarly, uh solves

P Q̃IU (uh ) < p = α.
Φ−1 (1 − α). From the definitions in (5), the values ul (α) and uh (α) can be approximated as
p
2p − 1 2
p(1 − p) −
(z1−α + 2) + O(n−3/2 ),
6n
p
2p − 1 2
(z1−α + 2) + O(n−3/2 ).
uh (α) = p + n−1/2 z1−α p(1 − p) −
6n
ul (α) = p − n−1/2 z1−α

For the lower one-sided CI, using (5), the 1 − α CI from Hutson (1999) is


−∞, Q̂LX uh (α)



.

(6)

Coverage probability is
n



o

P Q(p) ∈ −∞, Q̂LX uh (α)
= P Q̂LX uh (α) > Q(p)

  (1 −  )z
2


Thm 2
h
h 1−α exp{−z1−α /2} −1
I
h
√
= P Q̃X u (α) > Q(p) +
n + O n−3/2 [log(n)]3
2πuh (α)(1 − uh (α))

h (1 − h )z1−α φ(z1−α ) −1
=1−α+
n + O n−3/2 [log(n)]3 ,
p(1 − p)
where φ(·) is the standard normal PDF and the n−1 term is non-negative. Similar to the Ho
and Lee (2005a) calibration, we can remove the analytic n−1 term with the calibrated CI


 
h (1 − h )z1−α φ(z1−α ) −1
L
h
−∞, Q̂X u α +
n
,
p(1 − p)
10

(7)


which has CPE of order O n−3/2 [log(n)]3 . We follow convention and define CPE ≡ CP −
(1 − α), where CP is the actual coverage probability and 1 − α the desired confidence level.
By parallel argument, Hutson’s (1999) uncalibrated upper one-sided and two-sided CIs

also have O(n−1 ) CPE, or O n−3/2 [log(n)]3 with calibration. For the upper one-sided case,
again using (5), the 1 − α Hutson CI and our calibrated CI are respectively given by


Q̂LX

 
u (α) , ∞ ,
l

  
 
` (1 − ` )z1−α φ(z1−α ) −1
L
l
n
,∞ ,
Q̂X u α +
p(1 − p)

(8)

and for equal-tailed two-sided CIs,





Q̂LX ul (α/2) , Q̂LX uh (α/2)
and

 
` (1 − ` )z1−α/2 φ(z1−α/2 ) −1
L
l α
+
n
Q̂X u
,
2
p(1 − p)
 !
 

(1
−

)z
φ(z
)
α
h
h 1−α/2
1−α/2
+
n−1
.
Q̂LX uh
2
p(1 − p)

(9)

(10)

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 Kaplan (2015). 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 → ∞, 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 4 collects all of our results on coverage and power.
Theorem 4. Let zα denote the α-quantile of the standard normal distribution, and let h =
(n + 1)uh (α) − b(n + 1)uh (α)c and ` = (n + 1)ul (α) − b(n + 1)ul (α)c. Let Assumption A1
hold, and let A2 hold at p. Then, we have the following.

11

(i) The one-sided lower and upper CIs in (6) and (8) have coverage probability
1−α+


(1 − )z1−α φ(z1−α ) −1
n + O n−3/2 [log(n)]3 ,
p(1 − p)

with  = h for the former and  = ` for the latter.
(ii) The equal-tailed, two-sided CI in (9) has coverage probability
1−α+


[h (1 − h ) + ` (1 − ` )]z1−α/2 φ(z1−α/2 ) −1
n + O n−3/2 [log(n)]3 .
p(1 − p)

(iii) The calibrated one-sided lower, one-sided upper, and two-sided equal-tailed CIs given

in (7), (8), and (10), respectively, have O n−3/2 [log(n)]3 CPE.
(iv) The asymptotic probabilities of excluding Dn = Q(p) + κn−1/2 from lower one-sided (l),
upper one-sided (u), and equal-tailed two-sided (t) CIs (i.e., asymptotic power of the
corresponding hypothesis tests) are


Pnl (Dn ) → Φ(zα + S), Pnu (Dn ) → Φ(zα − S), Pnt (Dn ) → Φ zα/2 + S + Φ zα/2 − S ,
p
where S ≡ κf (F −1 (p))/ p(1 − p).
The equal-tailed property of our two-sided CIs is a type of median-unbiasedness. If (L̂, Ĥ)
is a CI for scalar θ, then an equal-tailed CI is “unbiased” under loss function L(θ, L̂, Ĥ) =
max{0, θ − Ĥ, L̂ − θ}, as defined in (5) of Lehmann (1951). This median-unbiased property
may be desirable (e.g., Andrews and Guggenberger, 2014, footnote 11), although it is different
than the usual “unbiasedness” where a CI is the inversion of an unbiased test. More generally,
in (9), we could replace ul (α/2) and uh (α/2) by ul (tα) and uh ((1−t)α) for t ∈ [0, 1]. Different
t may achieve different optimal properties, which we leave to future work.

12

4

Quantile inference: conditional

4.1

Setup and bias

Let QY |X (u; x) be the conditional u-quantile function of scalar outcome Y given conditioning
vector X ∈ X ⊂ Rd , evaluated at X = x. The object of interest is QY |X (p; x0 ), for p ∈ (0, 1)
and interior point x0 . The sample {Yi , Xi }ni=1 is drawn iid. Without loss of generality, let
x0 = 0.
If X is discrete so that P (X = 0) > 0, we can take the subsample with Xi = 0 and
compute a CI from the corresponding Yi values, using the method in Section 3. Even with
dependence like strong mixing among the Xi , CPE is the same O(n−1 ) from Theorem 4 as
a.s.

long as the subsample’s Yi are independent draws from the same QY |X (·; 0) and Nn  n.
If X is continuous, then P (Xi = 0) = 0, so observations with Xi 6= 0 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 x0 = 0.
Definition 1 (local smoothness). Following Chaudhuri (1991, pp. 762–3): if, in a neighborhood of the origin, function g(·) is continuously differentiable through order k, and its kth
derivatives are uniformly Hölder continuous with exponent γ ∈ (0, 1], then g(·) has “local
smoothness” of degree s = k + γ.
Assumption A3. Sampling of (Yi , Xi0 )0 is iid, for continuous scalar Yi and continuous vector
Xi ∈ X ⊆ Rd . The point of interest X = 0 is in the interior of X , and the quantile of interest
is p ∈ (0, 1).
Assumption A4. The marginal density of X, denoted fX (·), satisfies 0 < fX (0) < ∞ and
has local smoothness sX = kX + γX > 0.
13

Assumption A5. For all u in a neighborhood of p, QY |X (u; ·) (as a function of the second
argument) has local smoothness3 sQ = kQ + γQ > 0.
√
Assumption A6. As n → ∞, the bandwidth satisfies (i) h → 0, (i’) hb+d/2 n → 0 with
b ≡ min{sQ , sX + 1, 2}, (ii) nhd /[log(n)]2 → ∞.
Assumption A7. For all u in a neighborhood of p and all x in a neighborhood of the origin,

fY |X QY |X (u; x); x is uniformly bounded away from zero.
Assumption A8. For all y in a neighborhood of QY |X (p; 0) and all x in a neighborhood
of the origin, fY |X (y; x) has a second derivative in its first argument (y) that is uniformly
bounded and continuous in y, having local smoothness sY = kY + γY > 2.
Definition 2 refers to a window whose size depends on h: Ch = [−h, h] if d = 1, or more
generally a hypercube as in Chaudhuri (1991, pp. 763): letting k · k∞ denote the L∞ -norm,
Ch ≡ {x : x ∈ Rd , kxk∞ ≤ h},


Nn ≡ # {Yi : Xi ∈ Ch , 1 ≤ i ≤ n} .

(11)

Definition 2 (local sample). Using Ch and Nn defined in (11), the “local sample” consists
of Yi values from observations with Xi ∈ Ch ⊂ Rd , and the “local sample size” is Nn .
Additionally, let the local quantile function QY |X (p; Ch ) be the p-quantile of Y given X ∈ Ch ,

satisfying p = P Y < QY |X (p; Ch ) | X ∈ Ch ; similarly define the local CDF FY |X (·; Ch ),
local PDF fY |X (·; Ch ), and derivatives thereof.
Given fixed values of n and h, Assumption A3 implies that the Yi in the local sample
are independent and identically distributed,4 which is needed to apply Theorem 4. However, they do not have the quantile function of interest, QY |X (·; 0), but rather the biased
QY |X (·; Ch ). This is like drawing a global (any Xi ) iid sample of wages, Yi , and restricting
3

Our sQ corresponds to variable p in Chaudhuri (1991); Bhattacharya and Gangopadhyay (1990) use
sQ = 2 and d = 1.
4
This may be the case asymptotically even with substantial dependence, although we do not explore
this point. For example, Polonik and Yao (2002, p. 237) write, “Only the observations with Xt in a small
neighbourhood of x are effectively used. . . [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. . . .”

14

it to observations in Japan (X ∈ Ch ) when our interest is only in Tokyo (X = 0): our
restricted Yi constitute an iid sample from Japan, but the p-quantile wage in Japan may
differ from that in Tokyo. Assumptions A4–A6(i) and A8 are necessary for the calculation
of this bias, QY |X (p; Ch ) − QY |X (p; 0), in Lemma 5. Assumptions A6(ii) and A7 (and A3)
a.s.

ensure Nn → ∞. Assumptions A7 and A8 are conditional versions of Assumptions A2(i) and
A2(ii), respectively. Their uniformity ensures uniformity of the remainder term in Theorem
4, accounting for the fact that the local sample’s distribution, FY |X (·; Ch ), changes with n
(through h and Ch ).
From A6(i), asymptotically Ch is entirely contained within the neighborhoods implicit
in A4, A5, and A8. This in turn allows us to examine only a local neighborhood around p
−1/2

(e.g., as in A5) since the CI endpoints converge to the true value at a Nn

rate.

The Xi being iid helps guarantee that Nn is almost surely of order nhd . The hd comes
from the volume of Ch . Larger h lowers CPE via Nn but raises CPE via bias. This tradeoff
determines the optimal rate at which h → 0 as n → ∞. Using Theorem 4 and additional
results on CPE from bias below, we determine the optimal value of h.
Definition 3 (steps to compute CI for QY |X (p; 0)). First, Ch and Nn are calculated as
in Definition 2. Second, using the Yi from observations with Xi ∈ Ch , a p-quantile CI is
constructed as in Hutson (1999). If additional discrete conditioning variables exist, then
repeat separately for each combination of discrete conditioning values. This procedure may
be repeated for any number of x0 . For the bandwidth, we recommend the formulas in Section
4.3.
The bias characterized in Lemma 5 is the difference between these two population conditional quantiles.
Lemma 5. Define b as in A6 and let Bh ≡ QY |X (p; Ch ) − QY |X (p; 0). If Assumptions A4,
A5, A6(i), A7, and A8 hold, then the bias is of order |Bh | = O(hb ). Defining
ξp ≡ QY |X (p; 0),

(0,1)

FY |X (ξp ; 0) ≡

∂
FY |X (ξp ; x)
∂x
15

,
x=0

(0,2)

FY |X (ξp ; 0) ≡

∂2
FY |X (ξp ; x)
∂x2

,
x=0

with d = 1, kX ≥ 1, and kQ ≥ 2, the bias is
(0,2)
f (0)FY |X (ξp ; 0)
2 X

Bh = −h

(0,1)

+ 2fX0 (0)FY |X (ξp ; 0)

6fX (0)fY |X (ξp ; 0)

+ o(h2 ).

(12)

Equation (12) is the same as in Bhattacharya and Gangopadhyay (1990), who derive it
using different arguments.

4.2

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(hb + Nn h2b ) (two-sided) CPE from bias, but large
enough to control the O(Nn−1 ) CPE from applying the unconditional L-statistic method.
The following theorem summarizes optimal bandwidth and CPE results.
Theorem 6. Let Assumptions A3–A8 hold. The following results are for the method in
Definition 3. For a one-sided CI, the bandwidth h∗ minimizing CPE has rate h∗  n−3/(2b+3d) ,
corresponding to CPE of order O(n−2b/(2b+3d) ). For a two-sided CI, the optimal bandwidth
rate is h∗  n−1/(b+d) , and the optimal CPE is O(n−b/(b+d) ). Using the calibration in Section
3, if p = 1/2, then the nearly (up to log(n)) CPE-optimal two-sided bandwidth rate is h∗ 

n−5/(4b+5d) , yielding CPE of order O n−6b/(4b+5d) [log(n)]3 ; if p 6= 1/2, then h∗  n−3/(b+3d)

and CPE is O n−3b/(2b+6d) [log(n)]3 . The nearly CPE-optimal calibrated one-sided bandwidth

rate is h∗  n−2/(b+2d) , yielding CPE of order O n−3b/(2b+4d) [log(n)]3 .
As detailed in the supplemental appendix, Theorem 6 implies that for the most common
values of dimension d and most plausible values of smoothness sQ , 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
16

using a very large local polynomial if sQ is large; our approach implicitly uses a uniform
kernel, so bias reduction beyond O(h2 ) is impossible. Nonetheless, our method has smaller
CPE when d = 1 or d = 2 even if sQ = ∞, and in other cases the necessary local polynomial
degree may be prohibitively large given common sample sizes.
sQ = 2

800

10

600

sQ

400

−0.4
1

2

3

4

sQ
terms

0

−1.0

L−stat
L−stat (calib)
C91
5

1

d

2

3

4

5

6

0

200

5

−0.6
−0.8

CPE exponent

−0.2

L−stat
always
better

Number of polynomial terms

15

0.0

sQ needed to match L−stat CPE

7

d

Figure 2: Two-sided CPE comparison between new (“L-stat”) method and the local polynomial asymptotic normality method based on Chaudhuri (1991). Left: with sQ = 2 and
sX = 1, writing CPE as nκ , comparison of κ for different methods and different values of d.
Right: required smoothness sQ for the local polynomial normality-based CPE to match that
of L-stat, as well as the corresponding number of terms in the local polynomial, for different
d.
Figure 2 (left panel) shows that if sQ = 2 and sX = 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 n−2/(4+2d) . With d = 1, this is O n−1/3 , much



larger than our two-sided O n−2/3 . With d = 2, O n−1/4 is larger than our O n−1/2 . 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 2 (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
17

sQ ≥ 12 and a local polynomial of degree kQ ≥ 11. Since interaction terms are required,
Pk +d−1 T 
an 11th-degree polynomial has TQ=d−1 d−1
= 364 terms, which requires a large Nn (and
yet larger n). As d → ∞, the required number of terms in the local polynomial only grows
larger and may be prohibitive in realistic finite samples.

4.3

Plug-in bandwidth

We propose a feasible bandwidth value with the CPE-optimal rate. To avoid recursive dependence on  (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 α-quantile of N (0, 1) is again denoted zα . We let B̂h denote the
(0,1)
estimator of bias term Bh ; fˆX the estimator of fX (x0 ); fˆX0 the estimator of fX0 (x0 ); F̂Y |X
(0,1)

(0,2)

(0,2)

the estimator of FY |X (ξp ; x0 ); and F̂Y |X the estimator of FY |X (ξp ; x0 ), with notation from
Lemma 5.
When d = 1, the following are our CPE-optimal plug-in bandwidths.
• For one-sided inference, let
2/7



z1−α


ĥ+− = n−3/7  h
i1/2 h
i
(0,2)
(0,1)
3 p(1 − p)fˆX
fˆX F̂Y |X + 2fˆX0 F̂Y |X

,

ĥ++ = −0.770ĥ+− .

(13)

(14)

For lower one-sided inference, ĥ+− should be used if B̂h < 0, and ĥ++ otherwise. For
upper one-sided inference, ĥ++ should be used if B̂h < 0, and ĥ+− otherwise.
• For two-sided inference with general p ∈ (0, 1),
1/3
p
2
(B̂h /|B̂h |)(1 − 2p) + (1 − 2p) + 4 
ĥ = n−1/3 
,
(0,2)
(0,1)
0
ˆ
ˆ
2 fX F̂
+ 2f F̂


Y |X

18

X

Y |X

(15)

(0,2)
(0,1)
which simplifies to ĥ = n−1/3 fˆX F̂Y |X + 2fˆX0 F̂Y |X

−1/3

with p = 0.5.

While we suggest the CPE-optimal bandwidths for moderate n, we suggest shifting toward a larger bandwidth as n → ∞. 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 ≤ 1000
and then moves toward a n−1/20 under-smoothing of the MSE-optimal bandwidth rate, as in
Fan and Liu (2016, p. 205).

5

Empirical application

We present an application of our L-statistic inference to Engel (1857) curves. Code is
available from the latter author’s website, and the data are publicly available.
Banks, Blundell, and Lewbel (1997) 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 (Office for National Statistics and Department for Environment, Food and Rural Affairs, 2012), 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.5
Table 1 shows unconditional L-statistic CIs for various quantiles of the budget share
5

http://www.ons.gov.uk/ons/datasets-and-tables/data-selector.html?cdid=D7BT&dataset=
mm23&table-id=1.1

19

Table 1: L-statistic 99% CIs for various unconditional quantiles (p) of the budget share
distribution, for different categories of expenditure described in the text.
Category
p = 0.5
p = 0.75
p = 0.9
food
(0.1532,0.1580) (0.2095,0.2170) (0.2724,0.2818)
fuel
(0.0275,0.0289) (0.0447,0.0470) (0.0692,0.0741)
clothing (0.0135,0.0152) (0.0362,0.0397) (0.0697,0.0761)
alcohol
(0.0194,0.0226) (0.0548,0.0603) (0.1012,0.1111)
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 3 is comparable to Figure 1 of Banks et al. (1997) 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 − α/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 Fan and Liu (2016).

6

Simulation study

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

Joint: Food
0.50−quantile
0.75−quantile
0.90−quantile

●

●

●

●
●

0.10

Joint: Clothing/footwear

0.12

6.5

●

●

●
●

●

0.00

●

6.5

●

●

●

6.0

●
●

●
●

7.0

●
●

●
●

●
●

6.0

log expenditure

●
●

7.0

0.75−quantile
0.90−quantile

●
●

●
●

6.5

0.05

budget share
●

●
●

●
●

log
● expenditure
0.50−quantile

0.20

0.10
0.08
0.06

●
●

0.06

6.0

0.75−quantile
0.90−quantile

●

●

Joint: Alcohol

7.0

0.04
0.02
0.00

●

log
● expenditure
0.50−quantile

●

●

●
●
●

●

●
●

●
●

●

●

0.02

●

●

0.15

●
●

●

0.10

●

0.50−quantile
0.75−quantile
0.90−quantile

0.08

0.10

●

0.04

budget share

0.30
0.25
0.20

●

6.0

budget share

●

●

0.15

budget share

0.35

●

Joint: Fuel, light, and power

●

●

●
●

●

●

●

●

●
●

●

6.5

7.0

log expenditure

Figure 3: Joint (over the nine expenditure levels) 90% confidence intervals for quantile Engel
curves: food (top left), fuel (top right), clothing (bottom left), and alcohol (bottom right).

21

6.1

Unconditional simulations

We compare two-sided unconditional CIs from the following methods: “L-stat” from Section
3, originally in Hutson (1999); “BH” from Beran and Hall (1993); “Norm” using the sample
quantile’s asymptotic normality and kernel-estimated variance; “K15” from Kaplan (2015);
and “BStsym,” a symmetric Studentized bootstrap (99 draws) with bootstrapped variance
(100 draws).6
Overall, L-stat and BH have the most accurate coverage probability (CP), avoiding undercoverage 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 2: CP and median CI length, 1 − α = 0.95; n, p, and distributions of Xi (F ) shown
in table; 10,000 replications. “Too high” is the proportion of simulation draws in which the
lower endpoint was above the true F −1 (p), and “too low” is the proportion when the upper
endpoint was below F −1 (p).
n
25
25
25
25
25
25
25
25
25
25
25
25
25
25
25

p
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5

F
Normal
Normal
Normal
Normal
Normal
Uniform
Uniform
Uniform
Uniform
Uniform
Exponential
Exponential
Exponential
Exponential
Exponential

Method
L-stat
BH
Norm
K15
BStsym
L-stat
BH
Norm
K15
BStsym
L-stat
BH
Norm
K15
BStsym

CP Too low
0.953
0.022
0.955
0.021
0.942
0.028
0.971
0.014
0.942
0.028
0.953
0.022
0.954
0.021
0.908
0.046
0.963
0.018
0.937
0.031
0.953
0.024
0.954
0.024
0.924
0.056
0.968
0.022
0.941
0.039

Too high
0.025
0.024
0.030
0.015
0.030
0.025
0.025
0.046
0.020
0.032
0.023
0.022
0.020
0.010
0.020

Length
0.99
1.00
1.02
1.19
1.13
0.37
0.37
0.35
0.44
0.45
0.79
0.80
0.75
0.96
0.93

Table 2 shows nearly exact CP for both L-stat and BH when n = 25 and p = 0.5.
6

Other bootstraps were consistently worse in terms of coverage: (asymmetric) Studentized bootstrap,
and percentile bootstrap with and without symmetry.

22

“Norm” can be slightly shorter, but it under-covers. The bootstrap has only slight undercoverage, 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 3: CP and median CI length, as in Table 2.
n
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99

p
0.037
0.037
0.037
0.037
0.037
0.037
0.037
0.037
0.037
0.037
0.037
0.037
0.037
0.037
0.037

F
Method CP Too low
Normal
L-stat 0.951
0.023
Normal
BH
NA
NA
Normal
Norm 0.925
0.016
Normal
K15
0.970
0.009
Normal BStsym 0.950
0.020
Cauchy
L-stat 0.950
0.022
Cauchy
BH
NA
NA
Cauchy
Norm 0.784
0.082
Cauchy
K15
0.957
0.002
Cauchy BStsym 0.961
0.002
Uniform L-stat 0.951
0.024
Uniform
BH
NA
NA
Uniform Norm 0.990
0.000
Uniform
K15
0.963
0.028
Uniform BStsym 0.924
0.053

Too high
0.026
NA
0.059
0.021
0.030
0.028
NA
0.134
0.041
0.037
0.026
NA
0.010
0.009
0.022

Length
1.02
NA
0.83
1.55
1.20
39.37
NA
18.90
36.55
48.77
0.07
NA
0.12
0.11
0.08

Table 3 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 6= 0.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 4 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 approx-

23

imation to determine uh and ul but still interpolates, while “Norm/floor” uses the normal
approximation with no interpolation as in equations (5) and (6) of Fan and Liu (2016, Ex.
2.1).
iid

Table 4: CP and median CI length, n = 19, Yi ∼ N (0, 1), 1 − α = 0.90, 1,000 replications,
various p. In parentheses below CP are probabilities of being too low or too high, as in Table
2. Methods are described in the text.

Method
L-stat
Normal
Norm/floor

Two-sided CP
(Too low, Too high)
p = 0.15
p = 0.25
p = 0.5
0.905
0.901
0.898
(0.048,0.047) (0.050,0.049) (0.052,0.050)
NA
0.926
0.912
(NA,NA)
(0.062,0.012) (0.045,0.043)
NA
0.913
0.876
(NA,NA)
(0.083,0.004) (0.087,0.037)

Median length
p = 0.15 p = 0.25 p = 0.5
1.20
1.03
0.93
NA

1.22

1.00

NA

1.47

0.91

Table 4 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.

6.2

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
(Koenker, 2012). The second (“boot”) is a local cubic method following Chaudhuri (1991)
but with bootstrapped standard errors; the bandwidth is L-stat’s multiplied by n1/12 to get
24

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 Qu and Yoon (2015),
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 Fan and Liu (2016), based on a symmetrized kNN estimator using a bisquare kernel; we use the code from their simulations.7 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 x0 as the point of interest, instead of x0 = 0; we also take d = 1, b = 2,
and focus on two-sided inference, both pointwise (single x0 ) and joint (over multiple x0 ).
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
Hotelling (1939) tube computations in plot.rqss. Each simulation has 1,000 replications
unless otherwise noted.

Figure 4 uses Model 1 from Fan and Liu (2016, p. 205): Yi = 2.5+sin(2Xi )+2 exp −16Xi2 +
iid

iid

0.5i , Xi ∼ N (0, 1), i ∼ N (0, 1), Xi ⊥
⊥ 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 4 (bottom left) shows power curves of the hypothesis tests corresponding to the
joint (over x0 ∈ {0, 0.75, 1.5}) CIs, varying H0 while maintaining the same DGP. The deviations of QY |X (p; x0 ) shown on the horizontal axis are the same at each x0 ; zero deviation
implies H0 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 4 is similar, but based on uniform confidence bands evaluated
at 231 different x0 . Only L-stat has nearly exact type I error rate and good power.
7

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.

25

0.5

1.0

FLb

40

60

80

boot
QYg

0

50
0.0

boot
QYg
FLb

L−stat
rqss

20

80
70

1−α
L−stat
rqss

60

CP (%)

90

Rejection Probability (%)

100

Pointwise Power

100

Pointwise Coverage Probability

1.5

0.0

Joint Power
X Curves

0.5

1.0

1.5

Uniform Power
Curves
X
QYg
FLb

80

100

rqss
boot

60

60
40

0

20

α
L−stat

40

Rejection Probability (%)

QYg
FLb

20

rqss
boot

80

α
L−stat

0

Rejection Probability (%)

100

Deviation: −0.098387 or 0.098387

−0.2

−0.1

0.0

0.1

0.2

−0.2

Deviation of Null from Truth

−0.1

0.0

0.1

0.2

Deviation of Null from Truth

Figure 4: Results from DGP in Model 1 of Fan and Liu (2016), n = 500, p = 0.5. Top left:
pointwise CP at x0 ∈ {0, 0.75, 1.5}, interpolated for visual ease. Top right: pointwise power
at the same x0 against deviations of ±0.1. Bottom left: joint power curves. Bottom right:
uniform power curves.

26

Next, we use the simulation setup of the rqss vignette in Koenker (2012), which in turn
came in part from Ruppert, Wand, and Carroll (2003, §17.5.1). Here, n = 400, p = 0.5,
d = 1, α = 0.05, and
iid

Xi ∼ Unif(0, 1),

Yi =

p

Xi (1 − Xi ) sin 2π(1 + 2−7/5 )/(Xi + 2−7/5 ) + σ(Xi )Ui ,

(16)

where the Ui are iid N (0, 1), t3 , Cauchy, or centered χ23 , and σ(X) = 0.2 or σ(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, x0 = 0.04, 0.06, . . . , 0.96; uniform confidence
bands are evaluated at 231 equispaced values of x0 .
Figure 5’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 x0 ), L-stat can under-cover by around five percentage
points. Otherwise, CP is near 1 − α for all x0 in all DGPs.
In contrast, with the exception of boot, the other methods can have significant undercoverage. As seen in the first two columns of Figure 5, rqss has under-coverage (as low as
50–60% CP) for x0 closer to zero. QYg has under-coverage with the χ23 and (especially)
Cauchy. FLb has good CP except with the Cauchy, where CP can dip below 70%.
Figure 5’s third column shows the joint power curves. The horizontal axis of the graphs
indicates the deviation of H0 from the true values. For example, letting ξp,j be the true
conditional quantiles at the j = 1, . . . , 47 values of x0 (say, xj ), −0.1 deviation refers to
H0 : {QY |X (p; xj ) = ξp,j − 0.1 for j = 1, . . . , 47} (which is false), and zero deviation means
H0 is true. Our method’s type I error rate is close to α under all four Ui distributions (5.7%,
5.8%, 7.3%, 6.3%). In contrast, other methods show size distortion under Cauchy and/or χ23
Ui ; among them, boot is closest but still has 10.3% type I error rate with the χ23 . Next-best
is rqss; size distortion for FLb and QYg is more serious. L-stat also has the steepest joint

27

Pointwise Coverage Probability

X

Pointwise Coverage Probability

X

X

Pointwise Coverage Probability

X

X

0.6

0.4

0.6

X

0.8

80
70

0.2

0.4

0.6

X

boot
QYg
FLb
0.8

80
60
40

0.3

20
80
60
40
20

boot
QYg
FLb

−0.2

Joint Power Curves
−0.1

0.0

0.1

0.2

0.3

Deviation of Null from Truth

80
60
40

α
L−stat
rqss

20

Rejection Probability (%)

0.8

1−α
L−stat
rqss

60

CP (%)
0.2

boot
QYg
FLb

50

90
80
70
50

60

1−α
L−stat
rqss

−0.3

−0.3

boot
QYg
FLb

−0.2

Joint Power Curves
−0.1

0.0

0.1

0.2

0.3

Deviation of Null from Truth

100

0.4

100

0.2

90

100

Pointwise
Coverage
Probability
0.2
0.4
0.6
0.8

0.2

0

boot
QYg
FLb

0.1

80

CP (%)

80

1−α
L−stat
rqss

60
50

boot
QYg
FLb

Rejection Probability (%)

0.8

70

90
80
70
50

60

1−α
L−stat
rqss

0.0

α
L−stat
rqss

100

0.6

−0.1

60

100

0.4

Joint Power Curves

40

0.2

0.8

90

100

0.6

−0.2

Deviation of Null from Truth

0

boot
QYg
FLb

Pointwise Coverage Probability
0.4

−0.3

α
L−stat
rqss

20

CP (%)

80

1−α
L−stat
rqss

60

boot
QYg
FLb

50

50

70

90
80
70

1−α
L−stat
rqss

boot
QYg
FLb

0

0.8

100

0.6

Rejection Probability (%)

100

0.4

90

100

0.2

0.8

Rejection Probability (%)

100

Pointwise Coverage Probability

0.6

α
L−stat
rqss

boot
QYg
FLb

0

CP (%)

70

boot
QYg
FLb

X

0.4

60

CP (%)

1−α
L−stat
rqss

60
50

boot
QYg
FLb

Pointwise Coverage Probability

0.2

CP (%)

80

90

100
90
80

CP (%)

70
50

60

1−α
L−stat
rqss
0.2

CP (%)

Joint Power Curves

100

Pointwise Coverage Probability

−0.3

−0.2

−0.1

0.0

0.1

0.2

0.3

Deviation of Null from Truth

Figure 5: Pointwise CP (first two columns) and joint power curves (third column), 1 − α =
0.95, n = 400, p = 0.5, DGP in (16). Distributions of Ui are, top row to bottom row: N (0, 1),
t3 , Cauchy, and centered χ23 . Columns 1 & 3: σ(x) = 0.2; Column 2: σ(x) = (0.2)(1 + x).

28

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 χ23 .
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.

0.4

0.6

100
80
60

0.8

0.2

Pointwise
X Power

0.6

0.8

100

FLb

60

80

boot
QYg

40

60
40
20

0

0

L−stat
rqss

20

FLb

Rejection Probability (%)

boot
QYg

80

L−stat
rqss

0.4

X Power
Pointwise
Deviation: −0.1 or 0.1

Deviation: −0.1 or 0.1

100

FLb

0

0.2

Rejection Probability (%)

boot
QYg

40

60
40
20

L−stat
rqss

20

FLb

Rejection Probability (%)

boot
QYg

80

L−stat
rqss

Pointwise Power

0

Rejection Probability (%)

100

Pointwise Power

0.2

0.4

0.6

0.8

0.2

0.4

0.6

X

X

Deviation: −0.1 or 0.1

Deviation: −0.1 or 0.1

0.8

Figure 6: Pointwise power (described in text), 1 − α = 0.95, n = 400, p = 0.5, DGP from
(16), σ(x) = 0.2. The Ui are N (0, 1) (top left), t3 (top right), Cauchy (bottom left), and
centered χ23 (bottom right).
Figure 6 shows pointwise power. Specifically, for a given x0 , this is the proportion of

29

simulation draws in which QY |X (p; x0 ) − 0.1 is excluded from the CI, averaged with the
corresponding proportion for QY |X (p; x0 ) + 0.1. L-stat generally has the best power among
methods with correct CP (per first column of Figure 5).
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 = 105 and 100 different x0 , 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 CPEoptimal bandwidth. However, L-stat consistently has better power, greater robustness across
different conditional distributions, and less bias of its joint hypothesis tests.

7

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 Hutson (1999), which we improve to O n−3/2 [log(n)]3 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 (Goldman and Kaplan, 2016b), quantile marginal effects
(Kaplan, 2014), and entire distributions (Goldman and Kaplan, 2016a).

30

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 finitesample performance, e.g., by smoothing over discrete covariates (Li and Racine, 2007).

References
Abrevaya, J. (2001). The effects of demographics and maternal behavior on the distribution
of birth outcomes. Empirical Economics 26 (1), 247–257.
Alan, S., T. F. Crossley, P. Grootendorst, and M. R. Veall (2005). Distributional effects
of ‘general population’ prescription drug programs in Canada. Canadian Journal of Economics 38 (1), 128–148.
Andrews, D. W. K. and P. Guggenberger (2014). A conditional-heteroskedasticity-robust
confidence interval for the autoregressive parameter. Review of Economics and Statistics 96 (2), 376–381.
Banks, J., R. Blundell, and A. Lewbel (1997). Quadratic Engel curves and consumer demand.
Review of Economics and Statistics 79 (4), 527–539.
Beran, R. and P. Hall (1993). Interpolated nonparametric prediction intervals and confidence
intervals. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 55 (3),
643–652.
Bhattacharya, P. K. and A. K. Gangopadhyay (1990). Kernel and nearest-neighbor estimation of a conditional quantile. Annals of Statistics 18 (3), 1400–1415.
Bickel, P. J. (1967). Some contributions to the theory of order statistics. In Proceedings
of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1:
Statistics. The Regents of the University of California.
Buchinsky, M. (1994). Changes in the U.S. wage structure 1963–1987: Application of quantile
regression. Econometrica 62 (2), 405–458.
Chamberlain, G. (1994). Quantile regression, censoring, and the structure of wages. In
Advances in Econometrics: Sixth World Congress, Volume 2, pp. 171–209.
Chaudhuri, P. (1991). Nonparametric estimates of regression quantiles and their local Bahadur representation. Annals of Statistics 19 (2), 760–777.
Chen, S. X. and P. Hall (1993). Smoothed empirical likelihood confidence intervals for
quantiles. Annals of Statistics 21 (3), 1166–1181.
DasGupta, A. (2000). Best constants in Chebyshev inequalities with various applications.
Metrika 51 (3), 185–200.
David, H. A. and H. N. Nagaraja (2003). Order Statistics (3rd ed.). New York: Wiley.
Deaton, A. (1997). The analysis of household surveys: a microeconometric approach to
development policy. Baltimore: The Johns Hopkins University Press.
Donald, S. G., Y.-C. Hsu, and G. F. Barrett (2012). Incorporating covariates in the measurement of welfare and inequality: methods and applications. The Econometrics Journal 15 (1), C1–C30.
Engel, E. (1857). Die productions- und consumtionsverhältnisse des königreichs sachsen. Zeitschrift des Statistischen Bureaus des Königlich Sächsischen, Ministerium des
Inneren 8–9, 1–54.
31

Fan, X., I. Grama, and Q. Liu (2012). Hoeffding’s inequality for supermartingales. Stochastic
Processes and their Applications 122 (10), 3545–3559.
Fan, Y. and R. Liu (2016). A direct approach to inference in nonparametric and semiparametric quantile models. Journal of Econometrics 191 (1), 196–216.
Ferguson, T. S. (1973). A Bayesian analysis of some nonparametric problems. Annals of
Statistics 1 (2), 209–230.
Fisher, R. A. (1932). Statistical Methods for Research Workers (4th ed.). Edinburg: Oliver
and Boyd.
Goldman, M. and D. M. Kaplan (2016a). Evenly sensitive KS-type inference on distributions.
Working paper, available at http://faculty.missouri.edu/~kaplandm.
Goldman, M. and D. M. Kaplan (2016b). Nonparametric inference on conditional quantile
differences, linear combinations, and vectors, using L-statistics. Working paper, available
at http://faculty.missouri.edu/~kaplandm.
Hall, P. and S. J. Sheather (1988). On the distribution of a Studentized quantile. Journal
of the Royal Statistical Society: Series B (Statistical Methodology) 50 (3), 381–391.
Ho, Y. H. S. and S. M. S. Lee (2005a). Calibrated interpolated confidence intervals for
population quantiles. Biometrika 92 (1), 234–241.
Ho, Y. H. S. and S. M. S. Lee (2005b). Iterated smoothed bootstrap confidence intervals for
population quantiles. Annals of Statistics 33 (1), 437–462.
Hogg, R. (1975). Estimates of percentile regression lines using salary data. Journal of the
American Statistical Association 70 (349), 56–59.
Horowitz, J. L. and S. Lee (2012). Uniform confidence bands for functions estimated nonparametrically with instrumental variables. Journal of Econometrics 168 (2), 175–188.
Hotelling, H. (1939). Tubes and spheres in n-space and a class of statistical problems.
American Journal of Mathematics 61, 440–460.
Hutson, A. D. (1999). Calculating nonparametric confidence intervals for quantiles using
fractional order statistics. Journal of Applied Statistics 26 (3), 343–353.
Jones, M. C. (2002). On fractional uniform order statistics. Statistics & Probability Letters 58 (1), 93–96.
Kaplan, D. M. (2014). Nonparametric inference on quantile marginal effects. Working paper,
available at http://faculty.missouri.edu/~kaplandm.
Kaplan, D. M. (2015). Improved quantile inference via fixed-smoothing asymptotics and
Edgeworth expansion. Journal of Econometrics 185 (1), 20–32.
Kaplan, D. M. and Y. Sun (2016). Smoothed estimating equations for instrumental variables
quantile regression. Econometric Theory XX (XX), XX–XX. Forthcoming.
Koenker, R. (2012). quantreg: Quantile Regression. R package version 4.81.
Kumaraswamy, P. (1980). A generalized probability density function for double-bounded
random processes. Journal of Hydrology 46 (1–2), 79–88.
Lehmann, E. L. (1951). A general concept of unbiasedness. Annals of Mathematical Statistics 22 (4), 587–592.
Li, Q. and J. S. Racine (2007). Nonparametric econometrics: Theory and practice. Princeton
University Press.
Manning, W., L. Blumberg, and L. Moulton (1995). The demand for alcohol: the differential
response to price. Journal of Health Economics 14 (2), 123–148.
Muir, T. (1960). A Treatise on the Theory of Determinants. Dover Publications.
32

Neyman, J. (1937). »Smooth test» for goodness of fit. Skandinavisk Aktuarietidskrift 20 (3–
4), 149–199.
Office for National Statistics and Department for Environment, Food and Rural Affairs
(2012). Living Costs and Food Survey. 2nd Edition. Colchester, Essex: UK Data Archive.
http://dx.doi.org/10.5255/UKDA-SN-7472-2.
Pearson, K. (1933). On a method of determining whether a sample of size n supposed to have
been drawn from a parent population having a known probability integral has probably
been drawn at random. Biometrika 25, 379–410.
Peizer, D. B. and J. W. Pratt (1968). A normal approximation for binomial, F , beta, and
other common, related tail probabilities, I. Journal of the American Statistical Association 63 (324), 1416–1456.
Polansky, A. M. and W. R. Schucany (1997). Kernel smoothing to improve bootstrap confidence intervals. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 59 (4), 821–838.
Polonik, W. and Q. Yao (2002). Set-indexed conditional empirical and quantile processes
based on dependent data. Journal of Multivariate Analysis 80 (2), 234–255.
Pratt, J. W. (1968). A normal approximation for binomial, F , beta, and other common,
related tail probabilities, II. Journal of the American Statistical Association 63 (324),
1457–1483.
Qu, Z. and J. Yoon (2015). Nonparametric estimation and inference on conditional quantile
processes. Journal of Econometrics 185 (1), 1–19.
Rényi, A. (1953). On the theory of order statistics. Acta Mathematica Hungarica 4 (3),
191–231.
Robbins, H. (1955). A remark on Stirling’s formula. The American Mathematical
Monthly 62 (1), 26–29.
Ruppert, D., M. P. Wand, and R. J. Carroll (2003). Semiparametric Regression. Cambridge
Series in Statistical and Probabilistic Mathematics. Cambridge University Press.
Shorack, G. R. (1972). Convergence of quantile and spacings processes with applications.
Annals of Mathematical Statistics 43 (5), 1400–1411.
Shorack, G. R. and J. A. Wellner (1986). Empirical Processes with Applications to Statistics.
New York: John Wiley & Sons.
Stigler, S. M. (1977). Fractional order statistics, with applications. Journal of the American
Statistical Association 72 (359), 544–550.
Thompson, W. R. (1936). On confidence ranges for the median and other expectation distributions for populations of unknown distribution form. Annals of Mathematical Statistics 7 (3), 122–128.
Wilks, S. S. (1962). Mathematical Statistics. New York: Wiley.

A

Proof sketches and additional lemmas

The following are only sketches of proofs. The full proofs, with additional intermediate steps
and explanations, may be found in the supplemental appendix.

33

Sketch of proof of Proposition 1
For any u, let k = b(n +1)uc and  = (n + 1)u − k ∈ [0, 1). If  = 0, then the objects Q̃IX (u),
Q̂LX (u), and F −1 Q̃LU (u) are identical and equal to Xn:k . Otherwise, each lies in between Xn:k
and Xn:k+1 due to monotonicity of the quantile function and k/(n + 1) ≤ u < (k + 1)/(n + 1):
Xn:k = Q̃IX (k/(n + 1)) ≤ Q̃IX (u) ≤ Q̃IX ((k + 1)/(n + 1)) = Xn:k+1 ,
Xn:k ≤ Q̂LX (u) = (1 − )Xn:k + Xn:k+1 ≤ Xn:k+1 ,






Xn:k = F −1 Q̃LU k/(n + 1) ≤ F −1 Q̃LU (u) ≤ F −1 Q̃LU (k + 1)/(n + 1) = Xn:k+1 .
Thus, differences between the processes can be bounded by the maximum (over k) spacing
Xn:k+1 − Xn:k . Using the assumption that the density is uniformly bounded away from
zero over the interval of interest (and applying a maximal inequality from Bickel (1967,
eqn. (3.7))), this in turn can be bounded by a maximum of uniform order statistic spacings
Un:k+1 − Un:k . The marginal distribution (Un:k+1 − Un:k ) ∼ β(1, n) can then be used to bound
the probability as needed.

Lemma for PDF approximation
P
Lemma 7. Let ∆k be a positive (J + 1)-vector of natural numbers such that J+1
j=1 ∆kj =
Pj
n + 1, minj {∆kj } → ∞, and minj {n − ∆kj } → ∞, and define kj ≡ i=1 ∆ki and k ≡
(k1 , . . . , kJ )0 . Let X ≡ (X1 , . . . , XJ )0 be the random J-vector such that
∆X ≡ (X1 , X2 − X1 , . . . , 1 − XJ )0 ∼ Dirichlet(∆k).
Take any sequence an that satisfies conditions a) an → ∞, b) an n−1 [max{∆kj }]1/2 → 0, and
c) a3n [min{∆kj }]−1/2 → 0. Define Condition ?(an ) as satisfied by vector x if and only if
n
o
−1/2
max n∆kj |∆xj − ∆kj /n| ≤ an .
Condition ?(an )
j

Let kvk∞ ≡ maxj∈{1,...,k} |vj | denote the maximum norm of vector v = (v1 , . . . , vk )0 .
(i) Condition ?(an ) implies
n
o
−1/2
max n∆kj |∆xj − ∆kj /(n + 1)| = O(an ),
j
n
o
−1/2
max n∆kj |∆xj − (∆kj − 1)/(n − J)| = O(an ),
j

(17)
(18)

where ∆kj /(n + 1) and (∆kj − 1)/(n − J) are respectively the mean and mode of ∆Xj .
(ii) At any point of evaluation ∆x satisfying Condition ?(an ), the log Dirichlet PDF of

34

∆X may be uniformly approximated as

log f∆X (∆x) = D −

J+1
2 X

(n − J)
2

j=1



∆xj −

∆kj −1
n−J

2

∆kj − 1

+ Rn ,



J+1
1X
n
J
log
,
D ≡ log(n/2π) +
2
2 j=1
∆kj − 1
and Rn = O(a3n k∆k−1/2 k∞ ) uniformly (over ∆x). We also have the uniform (over
∆x) approximations


∂ log[f∆X (∆x)]
(n − J)2
∆kj − 1
= (n − J) −
∆xj −
+ O(a2n nk∆k−1 k∞ ),
∂∆xj
∆kj − 1
n−J

∂ log[f∆X (x)]
∂ log[f∆X (∆x)]
=−
+ O a2n nk∆k−1 k∞ .
∂∆kj /n
∂∆xj
(iii) Uniformly over all x ∈ RJ satisfying Condition ?(an ),
1
log[fX (x)] = D − (x − k/(n + 1))0 H(x − k/(n + 1)) + O(a3n k∆k−1/2 k∞ ),
2

∂ log[fX (x)]
= −H(x − k/(n + 1)) + O a2n nk∆k−1 k∞ ,
∂x
∂ log[fX (x)]
= H(x − k/(n + 1)) + O(a2n nk∆k−1 k∞ ),
∂k/(n + 1)
where the constant D is the same as in part (ii), and the J × J matrix H has non−1
and one off the diagonal
zero elements only on the diagonal H j,j = n2 ∆kj−1 + ∆kj+1
−1
2
H j,j+1 = H j+1,j = −n ∆kj+1 . The covariance matrix for x, V/n ≡ H −1 , has row i,
column j elements
V i,j = min(ki , kj )(n + 1 − max(ki , kj ))/[n(n + 1)],

(19)

connecting the above with the conventional asymptotic normality results for sample
quantiles. That is,
i
h

fX (x) = φV/n (x − k/(n + 1)) 1 + O a3n k∆k−1/2 k∞ ,
(20)

∂fX (x)
∂
=
φV/n (x − k/(n + 1)) + O a4n nJ/2 nk∆k−1 k∞ .
∂x
∂x

(21)

(iv) For the Dirichlet-distributed ∆X, Condition
?(an ) is violated
with only exponentially


−1 −a2n /2
decaying (in n) probability: 1 − P ?(an ) = O an e
.
(v) If instead there are asymptotically fixed components of the parameter
vector, the largest

of which is ∆kj = M < ∞, then with M = 1, 1 − P ?(an ) ≤ e−an −1 . With M ≥ 2,
35

√



for any η > 0, 1 − P ?(an ) = o a−1
.
n exp −an M (1/2 − η)

Sketch of proof of Lemma 7
The proof of part (i) uses the triangle inequality and the fact that the mean and mode differ
from ∆kj /n by O(1/n).
For part (ii), since ∆X ∼ Dirichlet(∆k), for any ∆x that sums to one,
log(f∆X (∆x)) = log(Γ(n + 1)) +

J+1 
X


(∆kj − 1) log(∆xj ) − log(Γ(∆kj )) .

(22)

j=1

Applying Stirling-type bounds in Robbins (1955) to the gamma functions,
h(∆x)

D

}|
{ z
}|
{
z




J+1
J+1
X
J
1X
n
n∆xj
log
+
(∆kj − 1) log
−J
log(f∆X (∆x)) = log(n/(2π)) +
2
2 j=1
∆kj − 1
∆kj − 1
j=1
+ O(k∆k−1 k∞ ),

(23)

where D is the same constant as in the statement of the lemma.
We then expand h(·) around the Dirichlet mode, ∆x0 . The cross partials are zero, the
first derivative terms sum to zero, and the fourth derivative is smaller-order uniformly over
∆x satisfying Condition ?(an ):
≡R1n =O(n−1 )

J+1
J+1
z }| {
X
1X
hj (∆x0 )(∆xj − ∆x0j ) +
hj,j (∆x0 )(∆xj − ∆x0j )2
h(∆x) = h(∆x0 ) +
2 j=1
j=1
≡R2n =O(a3 k∆k−1/2 k∞ )

≡R3n =O(a4 k∆k−1 k∞ )

n
n
z
}|
{ z
}|
{
J+1
J+1
X
X
1
1
+
hj,j,j (∆x0 )(∆xj − ∆x0j )3 +
hj,j,j,j (∆x̃)(∆xj − ∆x0j )4 ,
6 j=1
24 j=1

(24)

where the quadratic term expands to the form in the statement of the lemma.
The derivative with respect to ∆x is computed by expanding hj (∆x) = (∆kj − 1)∆x−1
j
P
around the mode, and then simplifying with Condition ?(an ) and the fact that J+1
(∆x
j−
j=1
∆x0,j ) = 1 − 1 = 0. The derivative with respect to ∆kj /n is computed from an expansion of
h(∆x) around the mode, reusing many results from the original computation of the PDF.
For part (iii), the results are intuitive given part (ii), so we defer to the supplemental
appendix. It is helpful that the transformation from the values Xj to spacings ∆Xj is
unimodular.
For part (iv), we use Boole’s inequality along with the beta tail probability bounds from
DasGupta (2000) and our beta PDF approximation from Lemma 7(ii).
P
For part (v), since ∆kj = M < ∞ is a fixed natural number, we can write ∆Xj = M
i=1 δi ,
where each δi is a spacing between consecutive uniform order statistics. The marginal distribution of each δi is β(1, n). Using the corresponding CDF formula and Boole’s inequality
36

leads to the result. The other bound may be derived using equation (6) in Inequality 11.1.1
in Shorack and Wellner (1986, p. 440), as seen in the supplemental appendix, but for our
quantile inference application we only use the first result since it is better for M = 1.
For violations of Condition ?(an ) in the other direction, the probability is zero for large
enough n since P (∆Xj < 0) = 0 and M 1/2 − an < 0 for large enough n.

Lemma for proving Theorem 2
First, we introduce notation. From earlier, u0 ≡ 0 and uJ+1 ≡ 1. For all j, kj ≡ b(n +1)uj c,
j ≡ (n + 1)uj − kj . Let ∆k denote the (J + 1)-vector such that ∆kj = kj − kj−1 , let
ψ = (ψ1 , . . . , ψJ )0 be the fixed weight vector from (4), and
Yj ≡ Un:kj ∼ β(kj , n + 1 − kj ),
Λj ≡ Un:kj +1 − Un:kj ∼ β(1, n),
√
Zj ≡ n(Yj − uj ),
X≡

J
X

ψj F −1 (Yj ),

∆Y ≡ (Y1 , Y2 − Y1 , . . . , 1 − YJ ) ∼ Dirichlet(∆k),
Vj ≡
X0 ≡

j=1

W≡

√


√  −1
n F (Yj ) − F −1 (uj ) ,
J
X

(25)

ψj F −1 (uj ),

j=1
0

n(X − X0 ) = ψ V,

1/2

W,Λ ≡ W + n

J
X

j ψj Λj [Q0 (uj ) + Q00 (uj )(Yj − uj )],

j=1

where the preceding variables are all understood to vary with n.
Let φΣ (·) be the PDF of a mean-zero multivariate normal distribution with covariance
Σ.
Lemma 8. Let Assumption A2 hold at ū, and let each element of Y and Λ satisfy
Condition ?(an ) (as defined in Lemma 7) with an = 2 log(n). The following results
hold uniformly over any u = ū + o(1).
(i) Let C be a J-vector of random interpolation coefficients as defined in Jones (2002):
each Cj ∼ β(j , 1 − j ), and they are mutually independent and independent of all other
random variables. Then,

n1/2 (LL − X0 ) − W,Λ = O n−3/2 [log(n)]3 ,
(26)

n1/2 (LI − X0 ) − WC,Λ = O n−3/2 [log(n)]3 .
(ii) Define V as the J × J matrix with row i, column j elements V i,j = min{ui , uj }(1 −
max{ui , uj }), and define A = diag{f (F −1 (u))},
i.e., Ai,j = f (F −1 (ui )) if i = j and

zero if i 6= j. Define Vψ ≡ ψ 0 A−1 V A−1 ψ ∈ R. For any realization λ of Λ =
(Λ1 , . . . , ΛJ ) satisfying Condition ?(2 log(n)),
sup
{w:?(2 log(n)) holds}


fW,Λ |Λ (w | λ)
− 1 = O n−1/2 [log(n)]3 ,
φVψ (w)

37

sup
{w:?(2 log(n)) holds}


∂fW,Λ |Λ (w | λ) ∂φVψ (w)
−
= O n−1/2 [log(n)]3+J ,
∂w
∂w

where the notation φVψ (·) denotes the PDF of a normal random variable with mean
zero and variance Vψ . For any value ˜ ∈ [0, 1)J , uniformly over K satisfying Condition
?(an ),



∂ 2 FW,Λ |Λ (K | λ)
∂φVψ (w)
2 0
2 2
−3/2
5+J
=
nψ
Q
(u
)
λ
+
O
n
[log(n)]
.
j
j
j
∂2j
∂w
w=K
=˜


Sketch of proof of Lemma 8
For part (i), with a Taylor expansion, the object LL may be rewritten as
L

L = X0 + n

−1/2

W,Λ +

J
X


L
L
ψj j νj,1
+ νj,2
,

j=1
000

L
νj,1
≡

Q (ũj )
[Yj − uj ]2 Λj ,
2

L
νj,2
≡

Q00 (ỹj ) 2
Λj ,
2

(27)


where ∀j, ỹj ∈ (Yj , Yj + Λj ) and ũj is between uj and Yj . The remainder is O n−3/2 [log(n)]3
by applying Condition ?(an ), noting that A2 uniformly bounds the quantile function derivatives for large enough n under Condition ?(an ). The argument for LI is essentially the
same.
For part (ii), since Λ contains finite spacings, we cannot apply Lemma 7(iii). Instead,
we use the result from 8.7.5 in Wilks (1962, p. 238),
(Λ1 , . . . , ΛJ , 1 − Λ1 − · · · − ΛJ ) ∼ Dirichlet(1, . . . , 1, n + 1 − J).
Directly approximating the corresponding PDF yields
log fΛ (λ) = J log(n) − n

J
X

λj + O(n−1 log(n)).

(28)

j=1


For the joint density of {Y, Λ}, define T ≡ ∆Y1 , Λ1 , ∆Y2 − Λ1 , . . . , ΛJ , ∆YJ+1 − ΛJ =
0
T Y0 , Λ0 , where det(T ) = 1 can be shown. Now
T ∼ Dirichlet(∆k1 , 1, ∆k2 − 1, . . . , 1, ∆kJ+1 − 1).
Using the formula for the PDF of a transformed vector and plugging in the Dirichlet PDF
formula for T yields the joint log PDF of Y and Λ. Combining this with the marginal log
PDF of Λ in (28) yields the log conditional PDF.
For the PDF of W, we can use the formula for the PDF of a transformed random
vector
and then expand around the uj . The transformation √
from Z ≡ Q(Y) to V ≡
√
n[Q(Y) − Q(u)] (from (25)) is straightforward centering and n-scaling. The last transformation is from V to W = ψ 0 V, as defined in (25). For the special case J = 1, as in our
38

quantile inference application, this step is trivial since W = V. Altogether, up to this point,


fW|Λ (w | λ) = φVψ (w) 1 + O n−1/2 [log(n)]3 ,

(29)

and it remains to account for the difference between W and W,Λ .
Altogether, it can be shown that the PDF of W conditional on Condition ?(an ) and Λ is


Z
Z
w − ψ1 v1 − · · · − ψJ−1 vJ−1
fW|?(an ),Λ (w | λ) = · · · fV|?(an ),Λ v1 , . . . , vJ−1 ,
| λ dv1 · · · dvJ−1 .
ψJ
?(an )

√ P
To transition to W,Λ , define η = n Jj=1 j ψj Λj [Q0 (uj ) + Q00 (uj )(Yj − uj )], so W,Λ =
W+η. Conditional on W = w, Y1 = y1 , . . . , YJ−1 = yJ−1 , the value of YJ is fully determined.
Additionally conditioning on Λ = λ, the value of η is fully determined. Along with the
implicit function theorem, this can be used to derive the final normal approximation.
Results for the PDF derivative follow the same sequence of transformations; details are
left to the supplemental appendix.
For the last result in this part of the lemma, in addition to Condition ?(2 log(n)) and A2,
we use the law of iterated expectations for CDFs.

Sketch of proof of Theorem 2
For part (i), we start by restricting attention to cases where the largest of the J spacings
between relevant uniform order statistics, Un:b(n+1)uj c+1 − Un:b(n+1)uj c , and the largest difference between the Un:b(n+1)uj c and uj satisfy Condition ?(2 log(n)) as in Lemma 7. By Lemma
7(iv,v), the error from this restriction is smaller-order. We then use the representation of
ideal uniform fractional order statistics from Jones (2002), which is equal in distribution
to the linearly interpolated form but with random interpolation weights Cj ∼ β(j , 1 − j )
instead of fixed j , where each Cj is independent of every other random variable we have.
The leading term in the error is due to Var(Cj ), and by plugging in other calculations from
Lemma 8, we see that it is uniformly O(n−1 ) and can be calculated analytically.
For part (ii), the first result comes from the FOC
" J
!#
∂ K exp{−K 2 /(2Vψ )} X ψj2 j (1 − j )
q
n−1
0=
2
−1
∂K
[f (F (uj ))]
2πV 3
j=1
ψ

p
whose solution K = Vψ is plugged into the expression in Theorem 2(i).
The additional result for LB in part (ii) follows from the Dirichlet PDF approximation
in Lemma 8(ii).

Sketch of proof of Lemma 3
The results are based on the Cornish–Fisher-type expansion from Pratt (1968) and Peizer
and Pratt (1968), solving for the high-order constants.

39

Sketch of proof of Theorem 4

For CP, let Rn = O n−3/2 [log(n)]3 be the remainder from Theorem 2(i).
For a lower one-sided CI,
uh (α) = p + O(n−1/2 ),

J = 1,

h = (n + 1)uh (α) − b(n + 1)uh (α)c,

uh (α)(1 − uh (α))
, X0 = F −1 (uh (α)),
f (F −1 (uh (α)))2
p
 −1

z1−α uh (α)(1 − uh (α))
1/2
−1
h
+ O(n−1/2 )
K=n
F (p) − F
u (α) = −
−1
h
f (F (u (α)))
p
−1/2
= −z1−α Vψ + O(n
),

Vψ =

where the first and last lines use Lemma 3, and the last line uses Assumption A2. Then, the
rate of coverage probability error is






L
h
L
h
−1/2
P Q̂X u (α) < Q(p) = P Q̂X u (α) < X0 + n
K



K exp{−K 2 /(2Vψ )}
h (1 − h )
q
+ Rn
= P Q̂IX uh (α) < X0 + n−1/2 K + n−1
[f (F −1 (uh (α)))]2
2πVψ3
= α − n−1 z1−α

h (1 − h )
φ(z1−α ) + O(n−3/2 ) + Rn ,
p(1 − p)

(30)

where f (F −1 (uh (α))) is uniformly (for large enough n) bounded away from zero by A2 since
uh (α) = p + O(n−1/2 ) → p. The argument for the lower endpoint is similar.
Two-sided CP comes directly from the two one-sided results, replacing α with α/2. For
the yet-higher-order calibration, the results follow from plugging in the proposed α̃.
L
The results for power are derived using the normal approximation Q̃B
X (u) of Q̂X (u),
along with a first-order Taylor approximation and arguments that the remainder terms are
negligible.

Sketch of proof of Lemma 5
Since the result is similar to other kernel bias results, and since the special case of d = 1 and
b = 2 is already given in Bhattacharya and Gangopadhyay (1990), we leave the proof to the
supplemental appendix and provide only a very brief sketch here. The approach is to start
from the definitions of QY |X (p; Ch ) and QY |X (p; x),
)
Z (Z
Z
QY |X (p;Ch )

p=

QY |X (p;x)

fY |X (y; x) dy fX|Ch (x) dx,
Ch

−∞

Z (Z

QY |X (p;Ch )

0=

fY |X (y; x) dy,
−∞

)
fY |X (y; x) dy fX|Ch (x) dx.

Ch

p=

QY |X (p;x)

40

so

After a change of variables to w = x/h, an expansion around w = 0 is taken, and the bias
can be isolated. If b = 2, kQ ≥ 2, and kX ≥ 1, then a second-order expansion is justified;
otherwise, the smoothness determines the order of both the expansion and the remainder.

Sketch of proof of Theorem 6
As in Chaudhuri (1991), we consider a deterministic bandwidth sequence, leaving treatment
of a random (data-dependent) bandwidth to future work. Whereas n is a deterministic
a.s.
sequence, Nn is random, but Nn  nhd as shown in Chaudhuri (1991). Another difference
with the unconditional case is that the local sample’s distribution, FY |X (·; Ch ), changes with
n (through h). The uniformity of the remainder term in Theorem 4 relies on the properties
of the PDF in Assumption A2. In the conditional case, we show that these properties hold
uniformly over the PDFs fY |X (·; Ch ) as h → 0, for which we rely on A4, A5, A7, and A8.
In the lower one-sided case, let Q̂LY |Ch (uh ) be the Hutson (1999) upper endpoint, with
notation analogous to Section 2, with uh = uh (α). The CP of the lower one-sided CI is


L
P QY |X (p; 0) < Q̂Y |Ch (uh ) = 1 − α + CPEU + CPEBias ,
(31)
where CPEU is CPE due to the unconditional method and CPEBias comes from the bias:



CPEU ≡ P QY |X (p; Ch ) < Q̂LY |Ch (uh ) − (1 − α) = O Nn−1 ,




CPEBias ≡ P QY |X (p; 0) < Q̂LY |Ch (uh ) − P QY |X (p; Ch ) < Q̂LY |Ch (uh ) .
1/2

Using Lemmas 5 and 8, or alternatively Theorem 2, one can show CPEBias = O(Nn hb ).
Then, one can solve for the h that equates the orders of CPEU and CPEBias , i.e., so that
1/2
Nn−1  Nn hb , using Nn  nhd .
With two-sided inference, the lower and upper endpoints have opposite bias effects. For
the median, the dominant terms of these effects cancel completely. For other quantiles,
there is a partial, order-reducing cancellation. The calculations, which use Theorem 2, are
extensive and thus left to the supplemental appendix. Ultimately, it can be shown that twosided CP is 1 − α plus terms of O(Nn−1 ), O(Bh ), and O(Bh2 Nn ), in addition to smaller-order
remainders. With the new CPE terms, one can again solve for the h that sets the orders
equal.

41

