Journal of Nonparametric Statistics

ISSN: (Print) (Online) Journal homepage: https://www.tandfonline.com/loi/gnst20

Density derivative estimation using asymmetrickernels
Benedikt Funke & Masayuki Hirukawa
To cite this article: Benedikt Funke & Masayuki Hirukawa (11 Dec 2023): Density
derivative estimation using asymmetric kernels, Journal of Nonparametric Statistics, DOI:
10.1080/10485252.2023.2291430
To link to this article: https://doi.org/10.1080/10485252.2023.2291430

Published online: 11 Dec 2023.

Submit your article to this journal

Article views: 56

View related articles

View Crossmark data

Full Terms & Conditions of access and use can be found athttps://www.tandfonline.com/action/journalInformation?journalCode=gnst20

JOURNAL OF NONPARAMETRIC STATISTICS
https://doi.org/10.1080/10485252.2023.2291430

Density derivative estimation using asymmetric kernels
Benedikt Funkea and Masayuki Hirukawaba Institute for Insurance Studies, TH Köln – University of Applied Sciences, Köln, Germany; b Faculty of

Economics, Ryukoku University, Kyoto, Japan
ABSTRACTARTICLE HISTORYThis paper studies the problem of estimating the first-order derivative of an unknown density with support on R+ or [0, 1]. Nonparametric density derivative estimators smoothed by the asymmetric,
gamma and beta kernels are defined, and their convergence properties are explored. It is demonstrated that these estimators canattain the optimal convergence rate of the mean integrated squared
error n−4/7 when the underlying density has third-order smoothness. Superior finite-sample properties of the proposed estimators
are confirmed in Monte Carlo simulations, and usefulness of theestimators is illustrated in two real data examples.

Received 5 October 2022
Accepted 30 November 2023
KEYWORDS

Beta kernel; density
derivative estimation;
digamma function; gamma
kernel; trigamma function
JEL CLASSIFICATION
CODES

C13; C14MSCS 2010

62G08; 62G20

1. Introduction
The aim of this paper is to explore statistical properties of density derivative estimators
using asymmetric kernels. Throughout we focus exclusively on univariate random variables supports of which have at least one boundary (e.g. R+ or [0, 1]). Indeed, there
are many economic and financial variables values of which are either nonnegative or
limited within a certain interval by construction. For example, wages, incomes, consumption expenditures, and insurance claims (or financial losses) are nonnegative. Variables
expressed in the forms of shares or proportions, including expenditure and budget shares,
unemployment rates, and default and recovery rates, are bounded from above and below.
In addition, distributions of these nonnegative variables tend to be right-skewed andunimodal.
Asymmetric kernels are a viable device that can capture shapes of density curves withaforementioned features effectively. Because an asymmetric kernel is based on a probability density function (‘pdf’) having the same support as that of the curves, it is free of theboundary bias by construction. Furthermore, shapes of an asymmetric kernel vary acrossdesign points at which smoothing is made, signifying an adaptive adjustment in the degreeof smoothing. In this context, asymmetric kernel smoothing exhibits similarities to, andcan be considered a variation of, variable kernel (or bandwidth) methods. This adaptive
CONTACT Benedikt Funke
Institute for Insurance Studies, TH Köln – University of
benedikt.funke@th-koeln.de
Applied Sciences, Gustav-Heinemann-Ufer 54, 50968 Köln, Germany
This article has been corrected with minor changes. These changes do not impact the academic content of the article.
© 2023 American Statistical Association and Taylor & Francis

2

B. FUNKE AND M. HIRUKAWA

smoothing attribute proves particularly advantageous when the underlying distribution
has both dense and sparse regions.
To make our discussion more concrete, let the random variable X be drawn from a
univariate distribution having the pdf f (·). Because the support of X has a boundary, we
employ an asymmetric kernel to estimate the pdf. Let KJ (x,b) (·) be an asymmetric kernel
indexed by J that depends on the design point x and the smoothing parameter b. Givenn i.i.d. observations {Xi }ni=1 , the kernel density estimator (‘KDE’) of f (x) smoothed by anasymmetric kernel KJ is defined as
1
KJ (x,b) (Xi ) .
n i=1nf̂J (x) =(1)

This class of density estimators, introduced in the pioneered work by Chen (1999, 2000),is now popular in various practical applications.
Often practitioners’ focus is not only on the density itself but also on its derivatives. From theoretical and practical standpoints, it is of particular importance andinterest to estimate the first-order density derivative f (1) (x) = ∂f (x)/∂x among density derivatives of all orders. A natural first-order derivative estimator implied by the
KDE (1) is∂1f̂J(1) (x) :=f̂J (x) =
∂x
n

n

∂
i=1

∂x

KJ (x,b) (Xi ) .

(2)

The kernel density derivative estimator of this class may be abbreviated as ‘KDDE’ hereinafter, whenever no confusions occur. While we deal exclusively with the first-orderderivative estimation, extending the analysis in this paper to estimation of higher-orderdensity derivatives is straightforward; actually, this task is more complicated, and thesederivatives appear to be less relevant for practical issues.
Among all asymmetric kernels proposed so far, our main focuses are on the beta andgamma kernels proposed by Chen (1999, 2000), respectively. Our preference is based on(i) empirical relevance of the kernels in economics and finance and (ii) analytical tractability. For (i), Hirukawa (2018, Table 1.1) presents a list of applications of these kernels toestimation and testing problems. For (ii), unlike standard symmetric kernels, exploring
statistical properties of nonparametric estimators smoothed by asymmetric kernels relies
on kernel-specific and thus diversified approaches. The density derivative estimators are
not an exception. The gamma function is a building block for the beta and gamma kernels.
There is rich literature on approximation techniques to the gamma and related functions
as they are actively studied.
Contributions of this paper can be summarised in two respects. First, this paper studies
convergence properties of asymmetric KDDEs and presents conclusive bias and varianceexpansions for the first time in the literature, to the best of our knowledge. It is demonstrated that each of the gamma and beta KDDEs has a usual O(b) bias, whereas the varianceconvergence of each estimator is O(n−1 b−3/2 ) for the interior region and O(n−1 b−3 ) for theboundary region. As in the case of derivative estimators smoothed by symmetric kernels
and thus as expected, the rates of the variance are slower than those of the corresponding KDE (1); invoke that the variance of a KDE is O(n−1 b−1/2 ) for the interior region
and O(n−1 b−1 ) for the boundary region. Second, a few new techniques are employed forJOURNAL OF NONPARAMETRIC STATISTICS

3

the bias and variance expansions of the gamma and beta KDDEs. These expansions involve
evaluations of moments of log-transformed gamma and beta random variables. We address
this issue by exploring a moment-generating function-based approach inspired by Ye andChen (2017). Finally, combining this approach with double-inequalities on digamma andtrigamma functions by Gordon (1994) helps us determine expressions of the leading bias
and variance terms of the estimators.
There are potentially many applications of the KDDEs. First, these can be used to evaluate peaks, valleys and inflection points of the underlying density. A direct application arises
from the field of actuarial science. Understanding the mode of an actuarial loss distribution
is of particular significance for insurance companies, as the mode, representing the mostprobable payout amount of an insurance obligation, plays an essential role in the calculation of an appropriate premium. Second, the KDDEs can be employed to estimate density
scores, which are a key ingredient for (indirect) average derivative estimation of nonparametric regression curves. Third, there are several important economic applications. The
density derivative is a key component of the Engel curve, which is essential for empiricallyverifying the law of demand. It also plays a significant role in the estimation of elasticities; see Hildenbrand and Hildenbrand (1986) for a comprehensive understanding of theseconcepts, their motivations, and the economic conclusions drawn from them. Interested
readers may also consult Härdle et al. (1990) for further applications.
The remainder of this paper is organised as follows. In Section 2, the gamma and beta
KDDEs are defined and their convergence properties are explored under some regularity
conditions. Section 3 conducts Monte Carlo simulations to investigate finite-sample properties of the KDDEs in comparison with several competing density derivative estimators.
In Section 4, we apply the gamma and beta KDDEs to two real datasets to illustrate theirusefulness. Section 5 concludes. All proofs are provided in the Appendix.
This paper adopts the following notational conventions: for a > 0, (a) = ∞ a−1
1exp(−t) dt is the gamma function; for p, q > 0, B(p, q) = 0 yp−1 (1 − y)q−1 dy
0 t
denotes the beta function; (x) = ∂ log (x)/∂x =  (1) (x)/ (x) and 1 (x) =  (1) (x) =∂(x)/∂x are the digamma and trigamma functions, respectively; 1{·} signifies anindicator function; and (r)m = (r + m)/ (r) = r(r + 1) · · · (r + m − 1) represents thedPochhammer symbol. The expression ‘X = Y’ reads ‘A random variable X obeys the distribution Y’. The expression ‘Xn ∼ Yn ’ is used whenever Xn /Yn → 1 as n → ∞. Lastly, inorder to describe different asymptotic properties of an asymmetric kernel estimator acrosspositions of the design point x ∈ R+ (x ∈ [0, 1]), we denote by ‘interior x’ and ‘boundary
x’ a design point x that satisfies x/b → ∞ (x/b → ∞ and (1 − x)/b → ∞) and x/b → κ(x/b or (1 − x)/b → κ) for some κ ∈ (0, ∞) as n → ∞, respectively.

2. Estimators and their large-Sample properties
2.1. The estimators
Due to their popularity in a wide range of empirical studies and analytical tractability, hereinafter we confine our discussion on kernels chosen for the KDDE (2) to the gamma andbeta ones. Extending the scope of asymmetric kernels to other distribution-specific ones
is, in principle, feasible but falls outside the primary focus of this paper.

4

B. FUNKE AND M. HIRUKAWA

The first estimator employs the gamma kernel (J = G)
KG(x,b) (u) =ux/b exp (−u/b)
1 {u ≥ 0}bx/b+1  (x/b + 1)

for x ∈ R+ by Chen (2000). The gamma KDDE is defined asf̂G(1) (x) =

n
n
1 ∂
1
KG(x,b) (Xi ) :=LG(x,b) (Xi ) KG(x,b) (Xi ) ,n i=1 ∂x
n i=1

where

x

1
log u − log b − +1 .
bbThe second estimator is built on the beta kernel (J = B)
LG(x,b) (u) =KB(x,b) (u) =ux/b (1 − u)(1−x)/b
1 {u ∈ [0, 1]}
B {x/b + 1, (1 − x) /b + 1}

for x ∈ [0, 1] by Chen (1999). The beta KDDE is given byf̂B(1) (x) =1 ∂1
KB(x,b) (Xi ) :=LB(x,b) (Xi ) KB(x,b) (Xi ) ,
n
∂x
n
n

n

i=1

i=1

where


x

u
1−x
1
log
−
+1 +
+1 .
LB(x,b) (u) =
b
1−u
b
b
Before proceeding, we inspect shapes of two weight functions ∂KG(x,b) (u)/∂x and∂KB(x,b) (u)/∂x used for density derivative estimation. Figure 1 presents plots of the functions evaluated at five different design points; more specifically, ∂KG(x,b) (u)/∂x for x =Figure 1. Shapes of Weight Functions with a Fixed Smoothing Parameter (b = 0.1) at Various Design
Points.

JOURNAL OF NONPARAMETRIC STATISTICS

5

1, 2, 3, 4, 5 and ∂KB(x,b) (u)/∂x for x = 0.1, 0.3, 0.5, 0.7, 0.9 are presented. It is worthemphasising that for all plotted functions, the smoothing parameter value is fixed at
b = 0.1.
A few remarks are in order. First, each plot produces an asymmetric S-shape, and thesign changes from negative to positive around (but not exactly at) the design point, ingeneral. An exception is ∂KB(x,b) (u)/∂x for x = 1/2; in this case the sign change occurs
exactly at 1/2. To make this argument more precise, let u∗J be the interior solution to∂KJ (x,b) (u)/∂x = 0 for a given interior design point x. Then,⎧ ⎨ x + b/2 + O b2
  for J = G and x ∈ (0, ∞)x + (x − 1/2) b + O b2 for J = B and x ∈ (0, 1/2) ∪ (1/2, 1) .
u∗J =
(3)
⎩
x
for J = B and x = 1/2
The derivation of (3) is documented in the Appendix. Moreover, the weight functions implied by the gamma and beta kernels√have a sharp contrast to the one impliedby the normal kernel
KN (u) = exp(−u2 /2)/ 2π. The weight function ∂KN (u)/∂u =√−u exp(−u2 /2)/ 2π generates a reverse S-shape and is symmetric with respect to theorigin.
Second, shapes of the weight functions vary across the design points at which smoothingis made, and they become flatter as the design point moves away from the boundary forthe gamma case and from both boundaries for the beta case. Hence, it can be observedthat when smoothing is made at a boundary, the weighting functions put their maximum
weights (in magnitude) on the boundary. An eminent feature of an asymmetric kernel isthat the kernel shape changes according to the location of the design point. The weightfunctions do inherit the property of altering the amount of smoothing in a locally adaptive
manner by using a single smoothing parameter. This property makes them much more
appealing in empirical work, especially when data is unequally spaced.
2.2. Regularity conditions
To deliver convergence properties of each KDDE, we impose the following common
regularity conditions.
Assumption 2.1: The i.i.d. random sample {Xi }ni=1 is drawn from a univariate distribution
having support either on R+ (for the gamma estimator) or on [0, 1] (for the beta estimator).
Assumption 2.2: The third-order derivative of f (·) is uniformly bounded and Höldercontinuous of order ς ∈ (0, 1] over the entire support.
Assumption 2.3: The smoothing parameter b(= bn > 0) satisfies b + (nb3 )−1 → 0 as
n → ∞.
Random sampling in Assumption 2.1 makes it easier to derive the dominant terms inbias and variance of the estimators. Observe that supports of the kernel and underlying
distribution match. Therefore, no boundary bias occurs.
Boundedness and Hölder-continuity of the third-order density derivative f (3) (·) in
Assumption 2.2 are key requirements for the bias expansion of the estimators. The former

6

B. FUNKE AND M. HIRUKAWA

helps determine the leading bias term, whereas the latter controls the order of magnitudein the remainder term in the bias expansion. Boundedness and (Hölder or some other classof) continuity of the second-order derivative are typically assumed for density estimation.
Now that our target quantity is the first-order density derivative, these regularity conditions
are consistently transferred to the third-order derivative. Notice that Hölder-continuity off (3) (·) implies that there is a constant C ∈ (0, ∞) such that



 (3)
(4)
f (u) − f (3) (v) ≤ C |u − v|ς , ∀ u, v.
The condition (nb3 )−1 → 0 in Assumption 2.3 is stronger than (nb)−1 → 0, which istypically imposed for the variance expansion of KDEs smoothed by asymmetric kernels;see Chen (1999, 2000) for more details. It will be shown shortly that the MSE-optimalsmoothing parameter for each KDDE becomes b∗ = O(n−2/7 ) for interior x and b† =O(n−1/5 ) for boundary x; these convergence rates are indeed within the required range.
2.3. Convergence properties of the estimators
2.3.1. Local property
Bias-Variance Tradeoff. Asymptotic properties of the gamma and beta KDDEs are explored
below. We start from delivering their bias-variance tradeoffs in the following theorem.
(1)

Theorem 2.1: If Assumptions 2.1–2.3 hold, then f̂J (x) for J ∈ {G, B} admits the following
expansions as n → ∞.
(i) Bias: E f̂J(1) (x) := f (1) (x) + BJ x, f b + o (b) ,

where

⎧

⎪ 1  (2)(3)
  ⎨ 3f (x) + xf (x)BJ x, f = 2
3
1⎪⎩−2f (1) (x) + (1 − 2x) f (2) (x) + x (1 − x) f (3) (x)22for J = G
for J = B.

(ii) Variance:
⎧ 1


⎪
⎨ 3/2 VJ (x) f (x) + o n−1 b−3/2
nb(1)
Var f̂J (x) :=⎪⎩ 1  (2κ + 1) f (x) + o n−1 b−3
3
2κ+2
2
nb κ2
 (κ + 1)


where



 √
−1
4 πx3/2
VJ (x) =  √−14 π {x (1 − x)}3/2

for interior x,
for boundary x

for J = G.
for J = B

Before discussing the above theorem in detail, we describe a few new strategies taken
for its proof in the Appendix. The bias and variance expansions in the theorem involve

JOURNAL OF NONPARAMETRIC STATISTICS

7

evaluations of moments of log-transformed gamma and beta random variables. We tackle
this issue by developing a moment-generating function-based approach inspired by Ye andChen (2017). Finally, combining this approach with Gordon’s (1994) double-inequalities
on digamma and trigamma functions yields the leading bias and variance terms.
Theorem 2.1 implies that the pointwise bias convergence of f̂J(1) (x) is O(b), just likethat of the corresponding KDE. Observe that while BG (x, f ) is free of f (1) (x), BB (x, f )dinvolves this quantity. This is because the moments E{LB(x,b) (θx )(θx − x)j } for θx =Beta{x/b + 1, (1 − x)/b + 1} are rational functions of (x, b), whereas the momentsdE{LG(x,b) (ξx )(ξx − x)j } for ξx = G(x/b + 1, b) are polynomials of (x, b); see Appendix for
more details. In particular, E{LB(x,b) (θx )(θx − x)} is not exactly 1 unlike E{LG(x,b) (ξx )(ξx −x)}, and as a consequence, the term including f (1) (x) is left over in BB (x, f ).
(1)
Moreover, the pointwise variance convergence of f̂J (x) slows down from that of the
−1
−3/2
−1
−3
) and O(n b ) for interior and boundary x, respeccorresponding KDE; it is O(n b(1)tively, to be more precise. It follows from the bias and variance expansions that f̂J (x) is
consistent for f (1) (x) under the set of assumptions.
Furthermore, Chen (2000) calls the quantity VJ (x) the variance coefficient. Vari√ √ance coefficientsfor the gamma and beta KDEs are known to be (2 π x)−1 and√ √{2 π x(1 − x)}−1 , respectively. Observe that each variance coefficient decreases as thedesign point x moves away from the boundaries. For the gamma and beta KDDEs, this√√feature is preserved in that their variance coefficients are (4 π x3/2 )−1 and [4 π{x(1 −x)}3/2 ]−1 , respectively. It can be inferred from this property that the gamma KDDE works
well for estimating the density derivative over sparse areas with fewer observations (e.g.
the right-tail part), as in the gamma KDE. The sparse areas are typically observed in theapplications of asymmetric kernel functions mentioned in Section 1, which emphasises thepractical relevance of our approach once again.
Mean Squared Error (‘MSE’). It follows from Theorem 2.1 that the MSE for interior x of(1)f̂J (x) is given by



 2
(1)
MSE f̂J (x) = b2 BJ x, f+

1
1
VJ (x) f (x) + o b2 + 3/2
3/2
nb
nb

.

(5)

The smoothing parameter value that minimises the two leading terms on the right-handside of (5) is2/7
3VJ (x) f (x)∗bJ =n−2/7 .
 2
4 BJ x, f
Observe that the MSE-optimal smoothing parameter b∗J = O(n−2/7 ) = O(h∗2 ), where h∗is the MSE-optimal bandwidth for KDDEs using nonnegative symmetric kernels (see, e.g.
Wand and Jones 1995). Therefore, when best implemented, the approximation to the MSE
becomes

 
4/7


 2 3/7 1
1
(1)∗
MSE f̂J (x) ∼ 7
BJ x, fVJ (x) f (x)n−4/7 .
3
4
(1)

Observe that the optimal MSE of f̂J (x) for interior x becomes O(n−4/7 ), which is also the
optimal convergence rate in the MSE of the nonnegative symmetric KDDE. On the other8

B. FUNKE AND M. HIRUKAWA

hand, for boundary x, MSE{f̂J(1) (x)} = O(b2 + n−1 b−3 ), which yields the MSE-optimalsmoothing parameter bJ = O(n−1/5 ) and the optimal MSE of O(n−2/5 ).
†2.3.2. Global property(1)
The inferior rate in the optimal MSE of f̂J (x) for boundary x does not affect its global
property. If {BJ (x, f )}2 dx and VJ (x)f (x) dx are both finite, then applying the trimming argument in Chen (1999, 2000) yields the mean integrated squared error (‘MISE’) off̂J(1) (x) as
MISE f̂J(1) (x)


 2
2
=b
BJ x, f
dx +

1
nb3/2



VJ (x) f (x) dx + o b2 +

1
nb3/2

.

(6)

The smoothing parameter value that minimises the two leading terms on the right-handside of (6) is
2/7
 
3
f
dx
V
(x)
(x)
J
b∗∗
n−2/7 .

 2
J =
4 BJ x, fdx
Therefore, when best implemented, the approximation to the MISE becomes
MISE

∗∗





f̂J(1) (x) 
3/7  
4/7
 2
1 
1
∼7
BJ x, fdx
VJ (x) f (x) dx
n−4/7 .
3
4

It is worth remarking that the rate n−4/7 coincides with the optimal convergence rate of
the MISE for a KDDE under third-order smoothness in the underlying density f that is
implied by Theorem 1 of Müller and Gasser (1979) and Theorem of Stone (1980).

3. Finite-sample performance
In the Monte Carlo study below, finite-sample performances of the asymmetric KDDEs arecompared with those of other competing density derivative estimators. For each estimator,
how to select the tuning parameter is an important practical issue. We start our analysis
from the ‘oracle’ method, where a tuning parameter value that minimises a performancemeasure is taken as the optimal one. While assessing an estimator with the oracle tuning
parameter value plugged in can be interpreted as examining its performance in the bestcase scenario, this method is infeasible in reality. To make the asymmetric KDDEs fully
operational, we also examine a cross-validation (‘CV’) method.
3.1. Oracle method
3.1.1. Case #1: gamma KDDE
Our simulation study starts from appraising the gamma KDDE. Three density functions
with support on R+ are examined. The first case is the Weibull distribution with density

JOURNAL OF NONPARAMETRIC STATISTICS

and derivative
  
x k
k  x k−1,
exp −f (x) =
λ λ
λ

9





(k, λ) = (2, 1) ⇒ f (1) (x) = −4x2 + 2 exp −x2 .

The second case is the beta prime distribution (also known as the beta distribution of thesecond kind) with density and derivativef (x) =1xα−1,
B (α, β) (1 + x)α+β(α, β) = (3, 3) ⇒ f (1) (x) =−120x2 + 60x.
(1 + x)7

The third case is the Lomax distribution (also known as the Pareto distribution of the secondkind) by Lomax (1954) with density and derivativex −(α+1)2α1+f (x) =, (α, λ) = (1, 1) ⇒ f (1) (x) = −.
λ
λ(1 + x)3
Shapes of these densities and derivatives can be found in Figure 2. Observe that the Weibull
distribution has a thin tail with an exponential decay, whereas two others have thick tails
with polynomial decays. For each distribution, 1000 data sets of sample size n ∈ {250, 500}
are simulated.
The density derivative estimators compared are: (i) the KDDE using the normal kernel[KDDE-N]; (ii) the KDDEs obtained via the second- and third-order local polynomial
density estimators using the triangular kernel KT (u) = (1 − |u|)1{|u| ≤ 1} by Cattaneo
et al. (2020) [LPDE2-T, LPDE3-T]; and (iii) the KDDE using the gamma kernel [KDDE-G].
The choice of the normal kernel is based on the recommendation by Härdle et al. (1990),
who study the CV bandwidth choice for density derivative estimation and report superior
performance of this kernel.
The idea of local polynomial density estimation by Cattaneo et al. (2020) comes fromlocal polynomial regression estimation. Local polynomial regression-based techniques

Figure 2. True Densities, True Density Derivatives and Pointwise Averages of Density Derivative Estimates (Case #1, n = 500).

10B. FUNKE AND M. HIRUKAWA

typically require prebinning of the data, which results in choosing an additional tuning
parameter (i.e. the bin width) and decreasing the sample size at the stage of density estimation. By smoothing out the empirical cumulative distribution function (‘cdf’) directly,
Cattaneo et al. (2020) preserve the original sample size and avoid choosing an extra tuning
parameter. The KDDE can be 
obtained from their estimation procedure in the following manner. Let F̂(x) := (1/n) ni=1 1{Xi ≤ x} be the empirical cdf, which is a consistentxestimator of the unknown cdf F(x) = 0 f (u) du. Then, using n data points {F̂(Xi )}ni=1
as regressands, we run a local polynomial regression of order p (≥ 1). Also let rp (u) :=(1, u, u2 , . . . , up ) and e2 := (0, 0, 1, 0, . . . , 0) be a pth-order polynomial expansion andthe third (p + 1)-dimensional unit vector, respectively. Finally, f (1) (x) can be estimated by(1)f̂LP(x) := e2 β̂(x), whereβ̂ (x) = arg minb∈Rp+1 2F̂ (Xi ) − rp (Xi − x) b KT

n 

i=1

Xi − x,h

and h is the bandwidth. We consider two orders of polynomials, namely, p = 2, 3, because(1)of the difference in bias convergences of f̂LP(x) between the polynomial orders. To be more
(1)
2
precise, f̂LP (x) for p = 2 has O(h ) and O(h) biases for interior and boundary x, respectively, whereas the bias convergence of the one for p = 3 is O(h2 ) everywhere. The varianceconvergence of the estimator for p = 2, 3 is O(n−1 h−3 ); see the supplemental appendix of
Cattaneo et al. (2020) for more details.
As performance measures of an estimator f̄ (1) (·), we adopt the root integrated squared
error (‘RISE’)

 
∞
2
RISE f̄ (1) =f̄ (1) (x) − f (1) (x) dx
0and the integrated absolute deviation (‘IAD’)   ∞ (1)
IAD f̄ (1) =f̄ (x) − f (1) (x) dx,
0where each integral is approximated by the trapezoidal rule on an equally-spaced grid of
551 points over the interval [0.0, 5.5].
The value of a tuning parameter a ∈ {h, b} is chosen via the oracle method. For eachsimulated sample the RISE is computed over a grid of a, and then a minimiser of the RISE is
obtained. The minimiser is found on an equally-spaced grid of 100 points over the interval[0.005, 0.500].
Table 1 reports averages and standard deviations of two performance measures and tuning parameter values over 1000 Monte Carlo samples. It can be immediately found that foreach distribution, KDDE-G outperforms all other estimators. Moreover, its tuning parameter values are by far the smallest among all (and this tendency continues to hold for thebeta KDDE in Case #2 below). This is not surprising, because the MISE-optimal smooth−2/7 ) tends to be smaller than theing parameter of two asymmetric KDDEs b∗∗
J = O(nMISE-optimal bandwidth for a symmetric KDDE h∗∗ = O(n−1/7 ). While KDDE-N performs much better than two LPDEs for Weibull and beta prime distributions, it performs
poorly for the Lomax distribution, failing to capture concavity of the curve. There is slightJOURNAL OF NONPARAMETRIC STATISTICS

11

Table 1. Monte Carlo Results (Oracle: Case #1).
n = 250n = 500

Estimator

RISE

IAD

h or bRISE

IADh or b

KDDE-N0.4504(0.0756)
0.7956(0.0956)
0.7782(0.1443)
0.3867(0.1061)

0.4866(0.0999)
0.9786(0.0986)
0.9569(0.1268)
0.4296(0.1206)

Weibull
0.1858(0.0264)0.3866(0.0781)0.4314(0.0815)
0.0656(0.0133)

0.4008(0.0591)
0.7618(0.0689)
0.7065(0.1015)0.3255(0.0804)0.4206(0.0776)
0.9419(0.0718)
0.8986(0.0881)0.3612(0.0894)

0.1628(0.0220)0.3344(0.0758)
0.4125(0.0846)
0.0532(0.0105)

0.5519(0.0874)
0.7843(0.1261)
0.6782(0.1987)
0.4497(0.1327)0.7395(0.1032)
0.9856(0.1440)0.8424(0.1726)0.4442(0.1057)

Beta Prime
0.1647(0.0232)
0.3430(0.0870)0.4299(0.0892)
0.0560(0.0162)

0.4881(0.0709)
0.7365(0.0945)0.5907(0.1396)
0.3979(0.1088)

0.6546(0.0873)
0.9091(0.1223)0.7303(0.1213)0.3834(0.0854)0.1401(0.0161)0.2939(0.0744)
0.4030(0.0851)0.0441(0.0109)

0.9875(0.0184)
0.5343(0.1323)0.5302(0.1862)
0.3263(0.1459)

0.8223(0.0383)
0.6678(0.1264)0.6507(0.1659)0.2762(0.0969)

Lomax
0.5000(0.0000)
0.4210(0.0846)
0.4199(0.1325)0.1602(0.0612)0.9869(0.0126)0.5043(0.1023)
0.4727(0.1560)
0.2891(0.1258)

0.8054(0.0277)
0.6069(0.0973)
0.5601(0.1388)0.2354(0.0763)0.5000(0.0000)
0.3829(0.0919)0.4230(0.1009)
0.1291(0.0441)

LPDE2-T
LPDE3-T
KDDE-G

KDDE-NLPDE2-T
LPDE3-T
KDDE-G

KDDE-NLPDE2-T
LPDE3-T
KDDE-G

Note: Simulation averages and standard deviations (in parentheses) of two performance measures (‘RISE’ and ‘IAD’) andtuning parameter values (‘h or b’) are reported.

advantage of LPDE3-T over LPDE2-T, which appears to reflect better bias properties of theformer. For this reason, we do not include LPDE2-T in Figure 2.
Figure 2 exhibits plots of pointwise averages of KDDE-N, LPDE3-T and KDDE-G
for n = 500 over 1000 Monte Carlo samples. Although no significant differences can befound when x > 3 for all distributions, there are substantial differences over x ∈ [0, 3]. The
density derivative for the Weibull distribution is initially decreasing, and KDDE-G alone
successfully traces the shape. For the beta prime case, the density derivative starts from theorigin and has a sharp mode in the vicinity of the origin. While all three estimators fail tostart from the origin, KDDE-G captures the curvature better than two others.
KDDE-N exhibits erratic behaviour for the Lomax distribution; more specifically, theplot of KDDE-N is initially decreasing and then increasing, whereas the truth is monotone increasing everywhere. As discussed
in Section 2.1, the weight function for KDDE-N√∂KN (u)/∂u = −u exp(−u2 /2)/ 2π exhibits a reverse S-shape, unlike ∂KG(x,b) (u)/∂x(or ∂KB(x,b) (u)/∂x). When a density f (x) has support on R, it dies out as x → ±∞.
KDDE-N works well in this scenario, because the corresponding density derivative takes
a reverse S-shape. However, a density with support on R+ or [0, 1] is allowed to be monotone decreasing everywhere, and the corresponding density derivative becomes monotone
increasing in this case. Moreover, the amount of smoothing by symmetric kernels with asingle bandwidth parameter is fixed everywhere. Under these conditions KDDE-N tends tomistakenly assume a reverse S-shape of the density derivative and then produces a reverse

12

B. FUNKE AND M. HIRUKAWA

S-curve with a fictitious decreasing part. Also observe that the oracle bandwidth for thisestimator is persistently chosen at the upper bound 0.5. Such an oversmoothing bandwidthreflects that KDDE-N focuses more on smoothing out the flat part on the right tail. Noticethat although KDDE-G is built on an S-shape weight function, it can successfully handle
such an irregular shape by varying shapes of the function across design points.
3.1.2. Case #2: beta KDDE
Now we assess the beta KDDE. Again we investigate three density functions with support
on [0, 1]. The first case is the beta distribution with density and derivativef (x) =

xα−1 (1 − x)β−1,
B (α, β)√(α, β) = (2.5, 2) ⇒ f (1) (x) = −21.875 x (x − 0.6) .

The second case is the truncated normal distribution with density and derivative
2√ 1 exp − (x−μ)
2
2σ
2π σ

f (x) = 
 , (μ, σ ) = (0.2, 0.4)1−μ
 σ −  −μσ2(−125x + 25) exp − (5x−1)
8,⇒ f (1) (x) =√8 2π { (2) −  (−0.5)}
where (·) is the standard normal cdf. The third case is the Kumaraswamy distribution by
Kumaraswamy (1980) and Jones (2009) with density and derivativeβ−1
f (x) = αβxα−1 1 − xα,(α, β) = (1, 4) ⇒ f (1) (x) = −12 (1 − x)2 .

Shapes of these densities and derivatives can be found in Figure 3. For each distribution,
1000 data sets of sample size n ∈ {250, 500} are simulated, as before.
The KDDE using the beta kernel, labelled as KDDE-B, is compared with KDDE-N,
LPDE2-T and LPDE3-T. To implement all these estimators, we choose their tuning parameters via the oracle method. The RISE and IAD are again considered as performancemeasures, where each integral is approximated by the trapezoidal rule on an equally-spaced
grid of 101 points over the interval [0, 1].
Table 2 reports averages and standard deviations of two performance measures andtuning parameter values over 1000 Monte Carlo samples. We can see the dominanceof KDDE-B over the other three estimators. As before, pointwise averages of KDDE-N,
LPDE3-T and KDDE-B for n = 500 over 1000 Monte Carlo samples are plotted in Figure 3.
It can be immediately found that KDDE-B captures the entire shape of each density derivative well. KDDE-N continues to produce a reverse S-curve for each of three distributions,and plots of LPDE3-T are often off target.
3.2. CV method
Before proceeding to real data examples in the next section, we should investigate a CV
method. Our construction of a CV criterion is inspired by Härdle et al. (1990). The integrated squared error (‘ISE’) of KDDE-J using the smoothing parameter b for J ∈ {G, B}

JOURNAL OF NONPARAMETRIC STATISTICS

13

Figure 3. True Densities, True Density Derivatives and Pointwise Averages of Density Derivative Estimates (Case #2, n = 500).
Table 2. Monte Carlo Results (Oracle: Case #2).
n = 250
Estimator
KDDE-NLPDE2-T
LPDE3-T
KDDE-B

KDDE-NLPDE2-T
LPDE3-T
KDDE-B

KDDE-NLPDE2-T
LPDE3-T
KDDE-Bn = 500

IAD

h or bRISE

IADh or b1.3051(0.1652)2.2846(0.2247)
2.0021(0.3645)
1.1704(0.2668)

0.9098(0.1431)
1.8768(0.1823)1.6521(0.2705)
0.9069(0.2205)

Beta
0.0969(0.0154)0.2402(0.0589)0.4164(0.1015)0.0481(0.0103)

1.1891(0.1433)
2.1695(0.1656)
1.8529(0.2560)
1.0036(0.2013)0.8086(0.1280)
1.7873(0.1279)
1.5679(0.1844)0.7716(0.1676)

0.0831(0.0128)0.1931(0.0489)
0.3887(0.1044)0.0383(0.0074)

0.7895(0.0905)
0.9275(0.1519)
0.9925(0.1818)0.7031(0.2409)0.6871(0.0759)0.8087(0.1398)
0.8720(0.1639)0.5787(0.2047)

Truncated Normal
0.2397(0.0171)
0.3667(0.0921)0.3728(0.1096)0.1420(0.0606)

0.7896(0.0644)0.8988(0.1108)
0.9213(0.1168)
0.6243(0.2020)0.6857(0.0543)0.7892(0.1092)
0.8143(0.1160)
0.5052(0.1610)0.2377(0.0125)
0.3226(0.0901)
0.3853(0.1101)
0.1067(0.0326)

5.3098(0.0366)
2.7114(0.4150)
2.6403(0.5454)1.5498(0.5570)

3.6365(0.1751)1.9286(0.2333)
1.9679(0.3184)
0.9562(0.2980)

Kumaraswamy
0.3857(0.0780)0.2483(0.0971)
0.3727(0.1164)0.0735(0.0215)

5.3133(0.0238)2.6707(0.3271)2.5344(0.4104)
1.3378(0.4523)3.6553(0.1206)
1.9093(0.1707)1.9236(0.2533)
0.8121(0.2314)0.3846(0.0621)0.2076(0.0873)0.3334(0.1210)0.0601(0.0167)

RISE

Note: Simulation averages and standard deviations (in parentheses) of two performance measures (‘RISE’ and ‘IAD’) andtuning parameter values (‘h or b’) are reported.

14

B. FUNKE AND M. HIRUKAWA

is defined as
ISEJ (b) =

 

2

f̂J(1),b (x) − f (1) (x)

dx
 

 
2
2
(1)
(1)
(1)
f (1) (x) dx,=
f̂J ,b (x) dx − 2 f̂J ,b (x) f (x) dx +where f̂J(1),b (x) signifies the dependence of KDDE-J on b, and the third term of the righthand side does not depend on b and thus may be safely omitted from the CV criterion.
We also assume that f (0) = 0 for J = G and that f (0) = f (1) = 0 for J = B, where theseassumptions are called the ‘zero-boundary condition’ hereinafter. Then, by integration byparts, it holds thatf̂J(1),b (x) f (1) (x) dx = − f̂J(2),b (x) dx.
Therefore, the CV criterion is given by
CVJ1 (b) =

 

(1)

f̂J ,b (x)22  (2)f̂(Xi ) ,
n i=1 J ,b,−i
n

dx +

(7)

where the superscript ‘1’ in CVJ1 (b) signifies the CV criterion for the first-order density
derivative estimator, andf̂J(2),b,−i (x) =

n

 
1
∂2
KJ (x,b) Xj2
n−1∂x
j=1,j=iis the second-order density derivative estimate using the sample with the ith observation
eliminated. Specifically,∂2
KJ (x,b) (u)
∂x2
⎧


1
⎪
⎪
L2
(u) − 2 1 bx + 1 KG(x,b) (u)⎪⎪⎨ G(x,b)
b
= 
 
⎪

⎪
1
x
1−x
⎪ 2
⎪
+ 1 + 1
+1
⎩ LB(x,b) (u) − 2 1bbb

for J = G.


KB(x,b) (u)for J = B

Finally, the CV smoothing parameter b̂ is defined as the minimiser of CVJ1 (b).
The same Weibull and beta distributions as in Section 3.1 are chosen to check finitesample performances of KDDE-G and KDDE-B, respectively, because each distribution
satisfies the zero-boundary condition. For each distribution, 1000 data sets of sample sizen ∈ {250, 500} are simulated. Judging from performance measures in the corresponding
oracle cases, we focus exclusively on KDDE-N as a competitor, where the CV bandwidthĥ is chosen via the method by Härdle et al. (1990). As before, b̂ and ĥ are found on an
equally-spaced grid of 100 points over the interval [0.005, 0.500].
Table 3 reports averages, standard deviations and medians of two performance measures and CV tuning parameter values over 1000 Monte Carlo samples. The reason forJOURNAL OF NONPARAMETRIC STATISTICS

15

Table 3. Monte Carlo Results (CV).
n = 250n = 500

Estimator

RISE

IAD

ĥ or b̂

RISE

IAD

ĥ or b̂

KDDE-N

4.7562(10.9551)[0.7318]
0.4767(0.1397)[0.4550]

5.1578(11.4842)[0.8514]0.5202(0.1285)[0.5167]

Weibull (Case #1)
0.1277(0.0905)[0.1200]
0.0685(0.0245)[0.0750]

3.4975(7.4160)[0.6075]
0.3815(0.0903)[0.3712]3.8651(8.0525)[0.6965]
0.4213(0.0893)[0.4181]0.1141(0.0775)[0.1150]
0.0570(0.0174)[0.0600]

9.9746(17.7338)[1.9212]1.8429(0.9501)[1.4146]

7.4058(13.1133)[1.5355]
1.4384(0.7288)[1.1383]

Beta (Case #2)
0.0680(0.0469)[0.0650]
0.0443(0.0243)[0.0500]

6.6297(11.3958)[1.6902]1.4427(0.6113)[1.1894]

4.9283(8.5476)[1.3296]
1.1237(0.4790)[0.9381]0.0606(0.0398)[0.0550]
0.0361(0.0181)[0.0400]

KDDE-G

KDDE-N
KDDE-B

Note: Simulation averages, standard deviations (in parentheses) and medians (in brackets) of two performance measures(‘RISE’ and ‘IAD’) and CV tuning parameter values (‘ĥ or b̂’) are reported.

presenting medians is that KDDE-N generates extremely large averages and standard deviations of the performance measures. This phenomenon can be attributed to the fact that its
CV bandwidths are systematically shorter than the corresponding oracle ones. Such shortbandwidths introduce additional variability to the resulting derivative estimates and makes
them highly volatile and unstable. In contrast, CV smoothing parameters for KDDE-G and
KDDE-B are very close to the corresponding oracle ones. As a consequence, these estimators do not induce extra bias or variability even after they are implemented by the CVsmoothing parameters.

4. Real data examples
The first-order density derivative serves as a tool for finding locations of peaks and valleys
in the underlying pdf. Consequently, we are motivated to see how precisely two asymmetric KDDEs can detect such locations. For this exercise, these estimators are applied to two
datasets. Strictly speaking, features of curve estimates including peaks, valleys and theirlocations vary across resolutions (i.e. smoothing parameter values). Accordingly, it is desirable to conduct a scale-space analysis like the SiZer map by Chaudhuri and Marron (1999).
However, for simplicity, we exclusively examine KDEs and KDDEs implemented by CVsmoothing parameters.
4.1. Data descriptions
The first dataset, used in Vella and Verbeek (1998), comprises a sample of 545 full-timeworking males who completed their schooling by 1980 over the period from 1980 to 1987,
and it is originally taken from the National Longitudinal Survey (Youth Sample). We utilise
hourly wages, labelled as ‘Wage’, in year 1980. Since all observations are recorded in the
log-transformed scale, they are exponentiated to restore the original scale on R+ before theanalysis. The second one is the dataset on the Old Faithful Geyser in Yellow Stone National
Park, Wyoming, USA. While the literature suggests that several different versions of the16

B. FUNKE AND M. HIRUKAWA

Table 4. Descriptive Statistics.
n

Mean

SD

545

4.59

2.20
1.14

272

3.49

SK

Min.

Q1

Q2

Q3

Max.

(a) US Hourly Wages (in US dollars)
1.06
0.33
3.21

4.25

5.71

16.81

(b) Eruption Durations (in minutes)
−0.42
1.60
2.16

4.00

4.46

5.10

Note: n = sample size; Mean = average; SD = standard deviation; SK = skewness; Min. = minimum value; Q1 = first quartile; Q2 = median (i.e. second quartile); Q3 = third quartile; and Max.
= maximum value.

data exist, we employ what is used in Wasserman (2006). Our focus is on durations of 272
consecutive eruptions. Both datasets analysed in this section are openly available; see thedata availability statement for the URLs at which these can be downloaded.
4.2. Estimation methods
We estimate densities and density derivatives of wages and eruption durations nonparametrically using the gamma and beta kernels, respectively. Table 4 presents summary
statistics of two variables. While observations of wages may be directly smoothed by the
gamma kernel, those of eruption durations are not confined within [0, 1]. Considering thatthe longest eruption duration is less than six minutes, we map observations in the original
scale, denoted as X, onto the unit interval by dividing them by six, before employing thebeta kernel. Then, the density and density derivative estimates evaluated at a design point
x in the original scale become f̂B (x/6)/6 and f̂B(1) (x/6)/36, respectively, where f̂B (·) and(1)f̂B (·) are the beta KDE and KDDE using the transformed observations U = X/6 ∈ [0, 1].
Density and density derivative estimations are implemented by CV. We may safely apply
the CV criterion for density derivative estimation (7), as Table 4 indicates that the zeroboundary condition is fulfilled for densities of both wages (in the original scale) anderuption durations (in the transformed scale). Furthermore, the CV criterion for KDE-J
is
 
n
2
2
CVJ0 (b) =f̂J ,b,−i (Xi ) ,f̂J ,b (x) dx −ni=1where f̂J ,b,−i (x) is the density estimate using the sample with the ith observation eliminated. CV smoothing parameters are found on an equally-spaced grid of 200 points overthe interval [0.0025, 0.5000].
4.3. Results
Figure 4 displays plots of density and density derivative estimates of two variables. First,
CV smoothing parameter values for the KDE-G and KDDE-G using the wage data are(coincidentally) 0.1375 and 0.1375, respectively. The density estimate looks single-peaked,
having the peak located around 3.57. The density derivative estimate supports the locationof the peak in that its sign changes from positive to negative between 3.56 and 3.57.
Second, CV smoothing parameter values for the KDE-B and KDDE-B using the eruption duration data (in the transformed scale on [0, 1]) are 0.0025 and 0.0050, respectively.

JOURNAL OF NONPARAMETRIC STATISTICS

17

Figure 4. Estimates of Densities and Density Derivatives from Real Data.

The distribution of eruption durations is empirically known to be bimodal. The density
estimate captures this shape successfully and indicates that that two peaks are locatedaround 1.89 and 4.47, whereas a valley exists around 3.06. The density derivative estimatechanges its sign from positive to negative between 1.92 and 1.93 and between 4.44 and
4.45, and from negative to positive between 2.97 and 2.98. Again, locations of two peaks
and a valley of the density estimate roughly coincide with the information provided by thedensity derivative estimate. All in all, both KDDEs work fine in practice in the sense thateach derivative estimate can well detect locations of local maximums and minimums of thecorresponding density estimate.

5. Conclusion
In this paper, we have proposed to employ asymmetric kernels to estimate first-order density derivatives. The gamma and beta KDDEs are defined, and their large-sample properties
including the novel bias and variance expansions are explored. It is demonstrated that theseestimators can attain the optimal convergence rate of the MISE n−4/7 if the underlying density has third-order smoothness. Superior finite-sample properties of these estimators havebeen confirmed in Monte Carlo simulations, and two real data examples have illustrated
their usefulness.
It is of importance and interest to see whether these estimators can be applicable to otherproblems than detection of locations of peaks and valleys in a pdf. In particular, incorporating them into (indirect) average derivative estimation of nonparametric regression curves
appears to be promising. We leave this task for our future research.

Acknowledgments
The authors would like to thank the editor and two anonymous referees for their constructive
comments and suggestions.

Disclosure statement
No potential conflict of interest was reported by the author(s).

18

B. FUNKE AND M. HIRUKAWA

Funding
This work was supported by the Japan Society for the Promotion of Science under grant numbers
19K01595 and 23K01340 (M. Hirukawa).

Data availability
The datasets used in Section 4 are openly available at the following URLs:• http://qed.econ.queensu.ca/jae/1998-v13.2/vella-verbeek/ (US hourly
wages)• https://www.stat.cmu.edu/ ∼ larry/all-of-nonpar/data.html (eruption
durations of the Old Faithful Geyser).

References
Cattaneo, M.D., Jansson, M., and Ma, X. (2020), ‘Simple Local Polynomial Density Estimators’,
Journal of the American Statistical Association, 115, 1449–1455.
Chaudhuri, P., and Marron, J.S. (1999), ‘SiZer for Exploration of Structures in Curves’, Journal of theAmerican Statistical Association, 94, 807–823.
Chen, S.X. (1999), ‘Beta Kernel Estimators for Density Functions’, Computational Statistics & Data
Analysis, 31, 131–145.
Chen, S.X. (2000), ‘Probability Density Function Estimation Using Gamma Kernels’, Annals of theInstitute of Statistical Mathematics, 52, 471–480.
Gordon, L. (1994), ‘A Stochastic Approach to the Gamma Function’, American Mathematical
Monthly, 101, 858–865.
Härdle, W., Marron, J.S., and Wand, M.P. (1990), ‘Bandwidth Choice for Density Derivatives’, Journalof the Royal Statistical Society, Series B, 52, 223–232.
Hildenbrand, K., and Hildenbrand, W. (1986), ‘On the Mean Income Effect: A Data Analysis of the
U.K. Family Expenditure Survey’, in Contributions to Mathematical Economics: in Honor of Gérard
Debreu, ed. W. Hildenbrand and A. Mas-Colell, Amsterdam: North-Holland, pp. 247–268.
Hirukawa, M. (2018), Asymmetric Kernel Smoothing: Theory and Applications in Economics andFinance, Singapore: Springer.
Jones, M.C. (2009), ‘Kumaraswamy’s Distribution: A Beta-Type Distribution with Some Tractability
Advantages’, Statistical Methodology, 6, 70–81.
Kumaraswamy, P. (1980), ‘A Generalized Probability Density Function for Double-Bounded Random Processes’, Journal of Hydrology, 46, 79–88.
Lomax, K.S. (1954), ‘Business Failures: Another Example of the Analysis of Failure Data’, Journal of
the American Statistical Association, 49, 847–852.
Müller, H.-G., and Gasser, T. (1979), ‘Optimal Convergence Properties of Kernel Estimates of
Derivatives of a Density Function’, in Smoothing Techniques for Curve Estimation: Proceedingsof a Workshop Held in Heidelberg, April 2–4, 1979, ed. T. Gasser and M. Rosenblatt, Berlin:
Springer-Verlag, pp. 144–154.
Stone, C.J. (1980), ‘Optimal Rates of Convergence for Nonparametric Estimators’, Annals of Statistics,
8, 1348–1360.
Vella, F., and Verbeek, M. (1998), ‘Whose Wages Do Unions Raise? A Dynamic Model of Unionism
and Wage Rate Determination for Young Men’, Journal of Applied Econometrics, 13, 163–183.
Wand, M.P., and Jones, M.C. (1995), Kernel Smoothing, Boca Raton, FL: Chapman & Hall/CRC.
Wasserman, L. (2006), All of Nonparametric Statistics, New York: Springer.
Ye, Z.-S., and Chen, N. (2017), ‘Closed-Form Estimators for the Gamma Distribution Derived From
Likelihood Equations’, American Statistician, 71, 177–181.

JOURNAL OF NONPARAMETRIC STATISTICS

19

Appendix
A.1. A list of useful results
Before proceeding, we summarise useful results. All results but two recursive formulae are statedas lemmata. First two lemmata refer to the moments of log-transformed gamma and beta randomvariables. The third lemma provides approximations to certain types of digamma and trigamma(1)functions. These are key ingredients for the bias and variance expansions of f̂J (x).

A.1.1. Recursive formulae for digamma and trigamma functions
For z > 0,1
Ψ (z + 1) = Ψ (z) + , andz
1Ψ1 (z + 1) = Ψ1 (z) − 2 .
z

(A1)
(A2)

A.1.2. Lemma A.1
d

Let X = G(α, β). Then,
E X m log X = β m (α)m log β + Ψ (α + m) , 2
E log2 X = log β + Ψ (α) + Ψ1 (α) .

m ∈ {0, 1, 2, . . .} ,and

A.1.3. Lemma A.2
d

Let X = Beta(γ , δ). Then,XE X m log
1−XX
E log2
1−X

=

(γ )m
{Ψ (γ + m) − Ψ (δ)} ,(γ + δ)mm ∈ {0, 1, 2, . . .} ,and= {Ψ (γ ) − Ψ (δ)}2 + Ψ1 (γ ) + Ψ1 (δ) .

A.1.4. Lemma A.3
For z > 0,


z
 
b
+ 1 = log
+ O b2 ,
+
b
b
2z
z
 b
 2
Ψ1
+1 = +O b ,
b
z
Ψ

z

and

as b → 0.

A.1.5. Proof of Lemma A.1

E(X log X) is available in the proof of Theorem 3.1 of Ye and Chen (2017). E(log X) and E(log2 X)
are also provided in the proof of Theorem 3.2 of Ye and Chen (2017). Finally, E(X m log X) for m ≥ 2
can be derived analogously to the proof of Lemma A.2 below.

A.1.6. Proof of Lemma A.2
The proof strategy largely comes from Ye and Chen (2017). Observe that the moment generatingdfunction of log{X/(1 − X)} for X = Beta(γ , δ) is given by




X
X


=E
M
(z) = E exp z log
Xlog 1−X
1−X
1−X

z

=

 (γ + z)  (δ − z).
 (γ ) (δ)

20

B. FUNKE AND M. HIRUKAWA

Therefore,


E log

E log2X
1−X=  (γ ) −  (δ) ,X
1−X= { (γ ) −  (δ)}2 + 1 (γ ) + 1 (δ) .

andTo obtain E[X m log{X/(1 − X)}] for m ≥ 1, consider that 1uγ −1 (1 − u)δ−1Xu=
um log
E X m log
du
1−X
1−u
B (γ , δ)0uγ +m−1 (1 − u)δ−1u
B (γ + m, δ) 1
log
du,
=
B (γ , δ)1−uB (γ + m, δ)0dwhere the integral corresponds to E[log{Y/(1 − Y)}] for Y = Beta(γ + m, δ). Then, the result
immediately follows.

A.1.7. Proof of Lemma A.3
By Theorems 4 and 5 of Gordon (1994),
log y −

 
1
1
1
1
<  y < log y −
−
− 
 ,
2
2y 12y2y 12 y + 1/14 2

and

  1
1
1
1
1
1
+ 
+
3 < 1 y < + 2 + 3
y 2y2
y
2y
6y
6 y + 1/14
for all y > 0. The stated results can be demonstrated by putting y = z/b + 1 and using (A1) and (A2).


A.2. Derivation of (3)(a) For J = G:
It suffices to solve
LG(x,b) (u) =


x
1
+1 =0
log u − log b − bbfor u. It follows from Lemma A.3 thatlog u = log b + Then,u∗G = x exp




 
b
+ 1 = log x ++ O b2 .
b
2x

x


 
 
 
b
b
b
+ O b2 = x 1 ++ O b2 = x + + O b2 .
2x2x2(b) For J = B:
Again, we solve
LB(x,b) (u) =


u
1
log
b
1−u

−

x
b


+1 +

1−x
+1
b

=0

for u. If x = 1/2, then this equation reduces to log{u/(1 − u)} = 0 and thus u∗B = 1/2 = x triviallyholds. If x = 1/2, then, by Lemma A.3,

x

 
1−x
x
1 − 2xu=
+1 −
+ 1 = log−
log
b + O b2 ,
1−u
b
b
1−x
2x (1 − x)JOURNAL OF NONPARAMETRIC STATISTICS

so that

 

 2
x
1 − 2x
exp −
b+O b
1−x
2x (1 − x)



 2
1 − 2x
x
1−
.
b+O b
1−x
2x (1 − x)

u
=
1−u
=
Therefore,



21



 
b + O b2


 
1−2x
1 − 2(1−x)
b + O b2





 
 
1 − 2x1 − 2x1+= x−
b + O b2
b + O b2
2 (1 − x)
2 (1 − x)

u∗B =

x−

1−2x
2(1−x)=x+ x− 
1b + O b2 .
2



A.3. Proof of Theorem 2.1
A.3.1. Proof of (i)(a) For J = G:dLet ξx = G(x/b + 1, b). Because the gamma kernel is the pdf of G(x/b + 1, b), we may write

  ∞


(1)
(A3)
E f̂G (x) =LG(x,b) (u) KG(x,b) (u) f (u) du = E LG(x,b) (ξx ) f (ξx ) .
0

Substituting
f (2) (x)f (3) (x̄)(ξx − x)2 +(ξx − x)32!
3!
with x̄ = tξx + (1 − t)x for some t ∈ [0, 1] into the right-hand side of (A3) yields
f (ξx ) = f (x) + f (1) (x) (ξx − x) +

3

 

f (j) (x) 
E LG(x,b) (ξx ) (ξx − x)jE f̂G(1) (x) =
j!
j=0



f (3) (x̄) − f (3) (x)3+E
LG(x,b) (ξx ) (ξx − x)
3!

=

3



 
AG j + RG say .

j=0

It also follows from Lemma A.1 and a property of gamma random variables thatx
  xxE LG(x,b) (ξx ) ξxm = bm−1
+1
+m+1 −
+1 .
b
b
b
m
By this and (A1),
E LG(x,b) (ξx ) = 0,
E LG(x,b) (ξx ) (ξx − x) = 1,
E LG(x,b) (ξx ) (ξx − x)2 = 3b, and
E LG(x,b) (ξx ) (ξx − x)3 = 3xb + 11b2 .

22

B. FUNKE AND M. HIRUKAWA

Therefore,
3

j=0



 
1  (2)
AG j = f (1) (x) +3f (x) + xf (3) (x) b + O b2 .
2

The remaining task is to demonstrate that RG = o(b). It follows from (4) and the Cauchy-Schwarz
inequality that
1/2
 1/2 
 C  2
Ct ζ 
|RG | ≤.
E |ξx − x|2(3+ζ )E LG(x,b) (ξx ) |ξx − x|3+ζ ≤E LG(x,b) (ξx )66Lemmata A.1 and A.3 jointly imply that


x

 
1
1
+ 1 = 2 O (b) = O b−1 .
E L2G(x,b) (ξx ) = 2 1bbb
Also, by Hölder’s inequality and E|ξx − x|8 = O(b4 ),(3+ζ )/4

= O b3+ζ .
E |ξx − x|2(3+ζ ) ≤ E |ξx − x|8Hence, |RG | ≤ O b−1/2 O b(3+ζ )/2 = O b1+ζ /2 = o (b) ,and thus the bias expansion is demonstrated.
(b) For J = B:dLet θx = Beta{x/b + 1, (1 − x)/b + 1}. Since the beta kernel is
(1)
Beta{x/b + 1, (1 − x)/b + 1}, it holds that E{f̂B (x)} := 3j=0 AB (j) + RB , where

the

pdf

of

  f (j) (x) 
E LB(x,b) (θx ) (θx − x)j ,AB j =j!

f (3) (x̄) − f (3) (x)3RB = E
LB(x,b) (θx ) (θx − x) ,3!
and x̄ = tθx + (1 − t)x for some t ∈ [0, 1] as before. It also follows from Lemma A.2 and a property
of beta random variables that

x



(x/b + 1)m   x
E LB(x,b) (θx ) θxm = b−1
+m+1 −+1 .
bb(1/b + 2)mCombining this with (A1) yields

E LB(x,b) (θx ) = 0,
E LB(x,b) (θx ) (θx − x) = 
1= 1 − 2b + O b2 ,1 + 2b
2x + 3b2x
E LB(x,b) (θx ) (θx − x)2 =−(1 + 2b) (1 + 3b) 1 + 2b 
= 3 (1 − 2x) b + O b2 , andE LB(x,b) (θx ) (θx − x)3 =3x2 + 12xb + 11b2
(1 + 2b) (1 + 3b) (1 + 4b)

3x2
3x (2x + 3b)+(1 + 2b) (1 + 3b) 1 + 2b
 
= 3x (1 − x) b + O b2 .
−It can be also shown that |RB | ≤ O(b1+ζ /2 ) = o(b), which establishes the stated result.

JOURNAL OF NONPARAMETRIC STATISTICS

23

A.3.2. Proof of (ii)(a) For J = G:
It follows from


 1 



(1)
2
Var f̂G (x) = E L2G(x,b) (X) KG(x,b)(X) + O n−1nd2that we focus on approximating E{L2G(x,b) (X)KG(x,b)(X)}. Let ηx = G(2x/b + 1, b/2). Then, as in the
proof of (i),




2
E L2G(x,b) (X) KG(x,b)(X) := BG E L2G(x,b) (ηx ) f (ηx ) ,where, by Chen (2000, p. 474),
BG =In addition,b−1  (2x/b + 1)∼
2x/b+12 2 (x/b + 1) √ √ b−1/2 / 2 π x b−1  (2κ + 1) / 22κ+1  2 (κ + 1)

if x/b → ∞.
if x/b → κE L2G(x,b) (ηx ) f (ηx ) ∼ f (x) E L2G(x,b) (ηx ) ,and it follows from Lemma A.1 that the expectation can be simplified as

2

x


2x
1
2
+1 −
+ 1 + log 2 + 1

E LG(x,b) (ηx ) = 2bbb
It follows from Lemma A.3 that
 −1

 b−1
b / (2x)+ O (1) ∼E L2G(x,b) (ηx ) =b−2 / (2κ)
2x

2x
+1
b


.

if x/b → ∞,if x/b → κwhich leads to the variance expansion.
(b) For J = B:
By

 1 



(1)
2
Var f̂B (x) = E L2B(x,b) (X) KB(x,b)(X) + O n−1 ,nwe have2
E L2B(x,b) (X) KB(x,b)(X) := BB E L2B(x,b) (ϑx ) f (ϑx ) ,dwhere ϑx = Beta{2x/b + 1, 2(1 − x)/b + 1}. Lemma of Chen (1999) implies thatB {2x/b + 1, 2 (1 − x) /b + 1}B2 {x/b + 1, (1 − x) /b + 1}
 −1/2  √ √b/ 2 π x (1 − x)∼b−1  (2κ + 1) / 22κ+1  2 (κ + 1)

BB =if x/b → ∞ and (1 − x) /b → ∞.
if x/b → κ or (1 − x) /b → κIn addition, E{L2B(x,b) (ϑx )f (ϑx )} ∼ f (x)E{L2B(x,b) (ϑx )}, where the expectation reduces to
E L2B(x,b) (ϑx )


x
1−x
2x
2 (1 − x)
1
+1 −
+1 −
+1 +
+1

= 2
b
b
b
b
b


2x
2 (1 − x)+ 1 + 1
+1
+1
b
b

2

24

B. FUNKE AND M. HIRUKAWA

by Lemma Lemma A.2. Then, Lemma A.3 implies thatb−1
E L2B(x,b) (ϑx ) =+ O (1)2x (1 − x)
 −1b / {2x (1 − x)} if x/b → ∞ and (1 − x) /b → ∞,∼if x/b → κ or (1 − x) /b → κ
b−2 / (2κ)
which completes the proof.