Computational Statistics and Data Analysis 75 (2014) 112–123

Contents lists available at ScienceDirect

Computational Statistics and Data Analysis
journal homepage: www.elsevier.com/locate/csda

Nonnegative bias reduction methods for density estimation
using asymmetric kernels
Masayuki Hirukawa a,∗ , Mari Sakudo b,c,1
aFaculty of Economics, Setsunan University, 17-8 Ikeda Nakamachi, Neyagawa, Osaka, 572-8508, Japanb

Research Institute of Capital Formation, Development Bank of Japan, 9-7, Otemachi 1-chome, Chiyoda-ku, Tokyo 100-8178, Japanc

Institute for Research in Contemporary Political and Economic Affairs, Waseda University, 1-6-1 Nishiwaseda, Shinjuku-ku, Tokyo
169-8050, Japan

articleinfo

Article history:
Received 2 July 2013
Received in revised form 10 December 2013
Accepted 7 January 2014
Available online 31 January 2014
Keywords:
Asymmetric kernel
Bias reduction
Boundary effectHigher-order bias kernel

abstract
Two classes of multiplicative bias correction (‘‘MBC’’) methods are applied to density estimation with support on [0, ∞). It is demonstrated that under sufficient smoothness ofthe true density, each MBC technique reduces the order of magnitude in bias, whereas
the order of magnitude in variance remains unchanged. Accordingly, the mean integratedsquared error of each MBC estimator achieves a faster convergence rate of O n−8/9 whenbest implemented, where n is the sample size. Furthermore, MBC estimators always generate nonnegative estimates by construction. Plug-in smoothing parameter choice rules forthe estimators are proposed, and their finite sample performance is examined via Monte
Carlo simulations.
© 2014 Elsevier B.V. All rights reserved.

1. Introduction
Hirukawa (2010) applied two classes of fully nonparametric multiplicative bias correction (‘‘MBC’’) methods originallyproposed for density estimation using symmetric kernels to estimate the density with support on [0, 1] via nonstandard
smoothing by the Beta kernel (Chen, 1999). This paper extends the analysis to estimating the density with support on [0, ∞)by asymmetric kernels (Chen, 2000; Jin and Kawczak, 2003; Scaillet, 2004). Let Kj(x,b) (·) be the asymmetric kernel indexed
by j that depends on a design point x > 0 and a smoothing parameter b > 0. Given a random sample {Xi }ni=1 drawn froma univariate distribution with density f that has support on [0, ∞), the density estimator using asymmetric kernel j can be
expressed asnfˆj,b (x) =1
n i =1Kj(x,b) (Xi ) .

(1)

Throughout, the kernel j refers to one of the Gamma (‘‘G’’), Modified Gamma (‘‘MG’’), Inverse Gaussian (‘‘IG’’), Reciprocal
Inverse Gaussian (‘‘RIG’’), Log-Normal (‘‘LN’’),2 and Birnbaum–Saunders (‘‘BS’’) kernels. Functional forms of these kernels arepresented in Table 1. Asymmetric kernels have originally emerged as an alternative to boundary correction methods; see,for instance, Karunamuni and Alberts (2005) for a brief review of the methods. Indeed, because all kernels in Table 1 have

∗ Corresponding author. Tel.: +81 72 839 8095; fax: +81 72 839 8138.
E-mail addresses: hirukawa@econ.setsunan.ac.jp (M. Hirukawa), marisakudo@gmail.com (M. Sakudo).
1 Tel.: +81 3 3244 1241; fax: +81 3 3270 7084.
2 Our definition of the LN kernel slightly differs from the original one in Jin and Kawczak (2003). This definition ensures that the leading variance of the √ 
density estimator (1) becomes n−1 b−1/2 f (x) / 2 π x for a design point x > 0 so that x/b → ∞ as b → 0.
0167-9473/$ – see front matter © 2014 Elsevier B.V. All rights reserved.
http://dx.doi.org/10.1016/j.csda.2014.01.012

M. Hirukawa, M. Sakudo / Computational Statistics and Data Analysis 75 (2014) 112–123

113

Table 1
Functional forms of asymmetric kernels.
Kernel (j)Functional form (u ≥ 0)

G (Chen, 2000)

KG(x/b+1,b) (u) = ux/b exp (−u/b) / bx/b+1 Γ (x/b + 1) .

 ρ (x)MG (Chen, 2000)

KMG(ρb (x),b) (u) = uρb (x)−1 exp (−u/b) / b b

IG (Scaillet, 2004)

x/b
for x ≥ 2b.
(1/4) (x/b)2 + 1 for x ∈ [0, 2b)1ux
KIG(x,1/b) (u) = √ 1 3 exp − 2bx
−
2
+
.
x
u

RIG (Scaillet, 2004)

KRIG(1/(x−b),1/b) (u) = √ 1−bexp − x2b

LN (Jin and Kawczak, 2003)

KLN (log x,b) (u) =BS (Jin and Kawczak, 2003)

KBS (b1/2 ,x) (u) =where ρb (x) =Γ {ρb (x)} ,



2π bu



2π bu

√1
exp
u 2π b√12x 2π b

−
 

 u
x−b

(log u−log x)2

x 1/2
u

2b

+





.

.

 x 3/2 u− 2 + x−u b1
exp − 2bux

− 2 + ux.

support on [0, ∞), they are free of boundary bias near the origin. Besides, the kernels have many other appealing properties,
including locally adaptive smoothing via changing their shapes and ‘shrinking variance’ with the position of x.3
Below we formally define two MBC estimators built on the density estimator (1). Throughout (1) is called the biasuncorrected (‘‘BU’’) estimator to distinguish it from MBC estimators. In the spirit of Terrell and Scott (1980, abbreviatedas ‘‘TS’’ hereafter), the first class of MBC techniques constructs a multiplicative combination of two density estimators
employing the same kernel but different smoothing parameters. Let fˆj,b/c (x) be the density estimator using asymmetrickernel j and smoothing parameter b/c, where c ∈ (0, 1) is some predetermined constant that does not depend on thedesign point x. Then, the TS-MBC asymmetric kernel density estimator can be defined asf˜TS ,j (x) = fˆj,b (x) 1−1 c 

− 1−c cfˆj,b/c (x).

(2)

On the other hand, the second class of MBC techniques due to Jones et al. (1995,abbreviated as ‘‘JLN’’ hereafter) utilizes asingle smoothing parameter b. In light of the identity f (x) = fˆj,b (x) f (x) /fˆj,b (x) , the JLN-MBC asymmetric kernel density
estimator can be defined asf˜JLN ,j (x) = fˆj,b (x)1  Kj(x,b) (Xi )nn i =1.

fˆj,b (Xi )(3)

Recognize that the term inside the bracket is a natural nonparametric estimator of the bias-correction term f (x) /fˆj,b (x).

Also observe that both f˜TS ,j (x) and f˜JLN ,j (x) are free of boundary bias and always generate nonnegative density estimates
everywhere by construction.
Following the convention, this paper refers to the position of x as ‘‘interior x’’ if x/b → ∞, and ‘‘boundary x’’ if x/b → κ
for some constant κ > 0, as b → 0. As demonstrated
shortly, under sufficient differentiability of f , bias convergence of each 
MBC estimator is accelerated from O (b) to O b2 , whereas the order of magnitude in variance remains unchanged from theone for (1), i.e. it is still O



nb1/2

−1 

for interior x. Accordingly, the mean integrated squared error (‘‘MISE’’) of each MBC

estimator for interior x takes the form of O b4 + n−1 b−1/2 . Therefore, when best implemented, each estimator can achieve −8/9 the convergence rate of O n
in MISE, which is faster than O n−4/5 , the MISE-optimal convergence rate within the classof nonnegative kernel estimators (Stone, 1980). Moreover, to implement MBC estimators employing the G and MG kernels,
this paper proposes plug-in methods of choosing the smoothing parameter b with gamma density used as a reference.
A few articles other than Hirukawa (2010) have investigated bias reduction methods for density estimation via
nonstandard smoothing when the support has a boundary. Hagmann and Scaillet (2007) and Gustafsson et al. (2009) study
semi-parametric MBC methods for density estimation with support on [0, ∞). Each method employs asymmetric kernels
at the bias correction step after initial parametric density estimation. Unlike MBC methods in this paper, their approaches
do not improve the bias convergence in order of magnitude. Moreover, Leblanc (2010) explores a bias reduction method for
estimating the density with support on [0, 1] using Bernstein polynomials, and establishes acceleration in bias convergence.
However, he adopts an additive bias correction, and thus the bias-corrected estimator does not always generate nonnegative
estimates unlike the one in Hirukawa (2010).
The remainder of this paper is organized as follows. Section 2 presents asymptotic properties of two MBC estimators. Section 3 proposes plug-in methods of choosing the smoothing parameter b for MBC estimators using the G and MG kernels, andconducts Monte Carlo simulations to check finite sample properties of the estimators. Section 4 applies two MBC techniques

3 It is an open question whether the asymmetric kernels studied here may fit with nonparametric analysis of functional or infinite-dimensional data by
Ferraty and Vieu (2006). While they consider the asymmetric kernels that take the form of K ((X − x) /b) for a data point X , design point x, and smoothing
parameter b, none of the kernels in Table 1 can be expressed in this form.

114

M. Hirukawa, M. Sakudo / Computational Statistics and Data Analysis 75 (2014) 112–123

to real data to illustrate how they work. Section 5 summarizes the main results of the paper and briefly mentions extensionsof MBC to joint density estimation. Outline proofs and formulae of plug-insmoothing parameters are given in the Appendix.
∞This paper adopts the following notational conventions. Γ (α) = 0 yα−1 exp (−y) dy (α > 0) denotes the gammadfunction. The expression ‘X = Y ’ reads ‘‘A random variable X obeys the distribution Y ’’. Lastly, the expression ‘Xn ∼ Yn ’ isused whenever Xn /Yn → 1 as n → ∞.
2. Main results
2.1. Asymptotic properties of MBC estimators
To develop convergence properties of MBC estimators, we make the following assumptions:
Assumption 1. f has four continuous and bounded derivatives, and f (x) > 0 for a given design point x > 0.
Assumption 2. The smoothing parameter b = bn (> 0) satisfies b → 0 and nbrj +5/2 → ∞ as n → ∞, where
1/2
1
3/2


rj =

for j = G, MG, RIG
for j = LN , BSfor j = IG.

The smoothness condition on f in Assumption 1 is standard for consistency of density estimators using fourth-orderkernels, whereas the positivity
for MBC. Assumption 2 implies that the convergence rate of the smoothing of f (x) is requiredparameter b is slower than O n−1/(rj +5/2) . We require this condition to control the order of magnitude in remainder terms

when approximating the bias of each MBC estimator. It will be shown shortly that the MSE-optimal smoothing parameterfor each estimator becomes b∗ = O n−2/9 for interior x and bĎ = O n−1/(rj +9/2) for boundary x; these convergence rates

are indeed within the required range.
Delivering bias and variance approximations for MBC estimators (2) and (3) requires (i) a second-order expansion of thebias termand (ii) an approximation to the leading variance of the BU estimator (1). For (i), Assumptions 1 and 2 imply thatE fˆj,b (x) = f (x) + a1,j (x, f ) b + a2,j (x, f ) b2 + o b2 , where a1,j (x, f ) and a2,j (x, f ) are kernel-specific functions that
depend on the design point x and derivatives of f . Using properties of the random variable corresponding to each kernel, we
can specify
forms of a1,j (x, f ) and a2,j (x, f ) as in Table 2. Furthermore, for (ii), the varianceapproximation
is given explicit



 


 −1

 −1

−1  √ r −1
f (x) := nb1/2
vj (x) f (x) for interior x and Var fˆj,b (x) = O nbrj +1/2
∼ nb1/2
2 πx jby Var fˆj,b (x)for boundary x.
Below we present two theorems on the approximations to bias and variance terms of two MBC estimators.
Theorem 1. If Assumptions 1 and 2 hold, then for a1,j (x, f ) and a2,j (x, f ), the bias of the TS-MBC estimator using kernel j canbe approximated byBias f˜TS ,j (x) ∼

1

 
1

a21,j (x, f )c (1 − c )

2

f (x)






− a2,j (x, f ) b2 :=1
c (1 − c )pj (x) b2 .

For vj (x) and rj , the variance of the TS-MBC estimator can be approximated by −1 
 n−1 b−1/2 λ (c ) vj (x) f (x) + o nb1/2
Var f˜TS ,j (x) =O nbrj +1/2 −1


where

λ (c ) =for boundary x,√

1 + c 5/2 (1 + c )1/2 − 2 2c 3/2



for interior x(1 + c )1/2 (1 − c )2.

Theorem 2. If Assumptions 1 and 2 hold, then the bias of the JLN-MBC estimator using kernel j can be approximated byBias f˜JLN ,j (x) ∼ −f (x) a1,j (x, h) b2 := qj (x) b2 ,where a1,j (x, h) is obtained by replacing f = f (x) in a1,j (x, f ) withh = hj (x, f ) :=a1,j (x, f )
f (x)

.

M. Hirukawa, M. Sakudo / Computational Statistics and Data Analysis 75 (2014) 112–123

115

Table 2
Explicit forms of a1,j (x, f ) and a2,j (x, f ).
Kernel (j)a1,j (x, f )Gf (x ) +a2,j (x, f )

x ′′
f
2

′


 x f ′′ (x)for x ≥ 2b

 ξ (x ) f ′ (x )for x ∈ [0, 2b) ,

2

MG

2

(x).

b

f ′′ (x) + 56 xf ′′′ (x) + x8 f ′′′′ (x).

x
x2


 f ′′′ (x) + f ′′′′ (x)381  2x  ′′ξb (x) + ξb (x) +f (x )
2

where ξb (x) = ρb (x) − bx = O (1).
IG

x3 ′′
f
2

RIG
LN
BS

for x ≥ 2b
for x ∈ [0, 2b) .

b

(x ).

x5 ′′′
f
2

x ′′
f
2

(x ).

1 ′′
f
2

(x) + 2x f ′′′ (x) + x8 f ′′′′ (x).

x ′
f
2

(x) + x2 f ′′ (x).

x ′
f
8

(x) + 78 x2 f ′′ (x) + 34 x3 f ′′′ (x) + x8 f ′′′′ (x).

x ′
f
2

x2

2

(x) + 2 f ′′ (x).

6(x) + x8 f ′′′′ (x).
23 2 ′′
x f
4

4

x4

(x) + 34 x3 f ′′′ (x) + 8 f ′′′′ (x).

For vj (x) and rj , the variance of the JLN-MBC estimator can be approximated by



−1 
 n−1 b−1/2 vj (x) f (x) + o nb1/2
Var f˜JLN ,j (x) =O nbrj +1/2 −1for interior x
for boundary x.

2.2. Discussions
2.2.1. Local property
Leading bias and variance terms. Because the support of asymmetric kernels matches that of the true density f , both TS- andJLN-MBC estimators are free of boundary bias. Morethese estimators reduce the orderof magnitude in bias
 importantly,from O (b) to O b2 , while their variances are O

 



nb1/2

−1

for interior x and O

nbrj +1/2−1

for boundary x. Observethat orders of variances remain unchanged from those for the corresponding BU estimator (1). We can also compare pj (x)and qj (x) in leading bias terms with the corresponding ones for nonnegative symmetric kernels. As stated in Jones andSignorini (1997, Sections 3.2–3.3), when the symmetric kernels are employed, the term corresponding to pj (x) for TS-MBC

2/f (x), and the term corresponding to qj (x) for JLN-MBC is proportional tof (x) f (x) /f (x) . The reason why pj (x) and qj (x) take more complicated forms is that while odd-order moments of 
symmetric kernels are exactly zero, those of asymmetric kernels around the design point x are often O (b) or O b2 . As aresult, extra density derivatives are included in pj (x) and qj (x).
is a linear combination of f ′′′′ (x) and f ′′ (x)



′′

 ′′

The variance of JLN-MBC estimators is first-order asymptotically equivalent to that of the corresponding BU estimator(1) for interior x. While the semi-parametric MBC density estimator by Hagmann and Scaillet (2007) also yields the sameleading variance term, this technique does not reduce the bias in order of magnitude. In contrast, when JLN-MBC is applied
for the density estimation using nonnegative symmetric kernels, the leading variance term tends to be larger (although not
inflated in order of magnitude) because the multiplier in the variance term involves the roughness (or squared integral) of
the ‘twiced’ kernel (Stuetzle and Mittal, 1979). Besides, since the multiplier λ (c ) in the variance for TS-MBC estimators isincreasing in c ∈ (0, 1), ranging from 1 to 27/16, the variance of these estimators tends to be larger than that of the BU
estimator (1) for interior x. Lastly (but not least importantly), asymptotic variances of both TS- and JLN-MBC estimators forinterior x are proportional to x−rj , even after MBC is made. The property of ‘shrinking variance’ with the position of x is equivalent to the strategy of using longer bandwidths over the tail region where the data are sparse. Therefore, MBC estimators
share the appealing variance property with the corresponding BU estimator.
Mean squared error (‘‘MSE’’). For interior x, the MSEs of f˜TS ,j (x) and f˜JLN ,j (x) can be approximated by

p2j (x)

MSE f˜TS ,j (x) =c 2 (1 − c )b4 + n−1 b−1/2 λ (c ) vj (x) f (x) + o b4 + n−1 b−1/2 ,



2





MSE f˜JLN ,j (x) = q2j (x) b4 + n−1 b−1/2 vj (x) f (x) + o b4 + n−1 b−1/2 .

The MSE-optimal smoothing parameters are2/9bTS ,j = c (1 − c ) λ (c )
∗

 2

2



vj (x) f (x)8p2j (x)2/9
n−2/9 ,

and

and

116

M. Hirukawa, M. Sakudo / Computational Statistics and Data Analysis 75 (2014) 112–123

b∗JLN ,j =vj (x) f (x)8q2j (x)2/9
n−2/9 ,

which yield the optimal MSEs
8/9 −8/992/9n,MSE ∗ f˜TS ,j (x) ∼ 8/9 γ (c ) pj (x) vj (x) f (x)88/9 −8/99 2/9MSE ∗ f˜JLN ,j (x) ∼ 8/9 qj (x) vj (x) f (x)
n
,
8
where

√

1 + c 5/2 (1 + c )1/2 − 2 2c 3/2γ (c ) =



and

8/9
.

c 1/4 (1 + c )1/2 (1 − c )9/4

Observe that the MSE-optimal smoothing parameters are O n−2/9 = O h∗2 , where h∗ is the MSE-optimal bandwidth for
fourth-order kernel estimators, or TS- or JLN-MBC estimators using nonnegative symmetric kernels. As a result, the optimal
MSEs of f˜TS ,j (x) and f˜JLN ,j (x) for interior x become O n−8/9 , as with MBC estimation using the second-order kernels. The









convergence rate is faster than O n−4/5 , the optimal convergence rate in the MSE of the corresponding BU estimator (1) forinterior x. On the other hand, for boundary x, the MSEs of f˜TS ,j (x) and f˜JLN ,j (x) are O b4 + n−1 b−(rj +1/2) , which yields the

MSE-optimal smoothing parameter bĎ = O n−1/(rj +9/2) and the optimal MSE of O n−4/(rj +9/2) . The optimal convergence









rate of MSEs is indeed faster than O n−2/(rj +5/2) , that of the BU estimator for boundary x.



2.2.2. Global property
The undesirable convergence rates over boundary regions do not affect the global properties of the MBC estimators. By
the trimming argument in Chen (2000), the MISEs of the MBC estimators are





 ∞

b4



MISE f˜TS ,j (x) =c2

MISE f˜JLN ,j (x) = b4

(1 − c ) ∞2p2j (x) dx +0

q2j (x) dx +0λ (c )

 ∞

nb1/2

1

 ∞

nb1/2

0



vj (x) f (x) dx + o b4 + n−1 b−1/2 ,

and

0



vj (x) f (x) dx + o b4 + n−1 b−1/2 ,provided that p2j (x) , q2j (x), and vj (x) are integrable.4 The MISE-optimal smoothing parameters are then given by2b∗∗
TS ,j = c (1 − c ) λ (c )

 2

2/9

 ∞
∗∗

bJLN ,j =vj (x) f (x) dx
∞ 2
8 0 qj (x) dx

 ∞vj (x) f (x) dx
∞
8 0 p2j (x) dx

2/9
n−2/9 ,

0

and

2/9
n−2/9 .

0

Therefore, the optimal MISEs become

MISE ∗∗ f˜TS ,j (x) ∼MISE ∗∗ f˜JLN ,j (x) ∼9γ (c )
8/9

 ∞

8

9

 ∞

88/9

0

2/9  ∞p2j (x) dx

0

8/9
vj (x) f (x) dx
n−8/9 ,and0

q2j (x) dx

2/9  ∞vj (x) f (x) dx

8/9

n−8/9 .

0

Furthermore, the multiplier γ (c ) in the optimal MISE for the TS-MBC estimator is minimized at c ∗ ≈ 0.2636; this value is
exclusively considered in subsequent analyses.
3. Finite sample performance
3.1. Monte Carlo setup
We evaluate the finite sample performance of two classes of MBC estimators via Monte Carlo simulations. From now
on we concentrate on the G and MG kernels due to their popularity in the literature. This simulation study compares thefollowing three classes of estimators: (i) BU estimators (1) [BU-G, BU-MG]; (ii) TS-MBC estimators (2) [TS-G, TS-MG]; and(iii) JLN-MBC estimators (3) [JLN-G, JLN-MG]. The value of the constant c in each TS-MBC estimator is set equal to the4 Whenever the integrated squared bias is considered, pMG(x) and qMG (x) refer to those for interior x (i.e. x ≥ 2b).

M. Hirukawa, M. Sakudo / Computational Statistics and Data Analysis 75 (2014) 112–123

117

Table 3
True distributions considered in Monte Carlo simulations.
Distribution

Density function f (x) , x ≥ 0

1. Gammaxα−1 exp (−x/β) / {β α Γ (α)} , (α, β) = (1.5, 1).

2. Weibull(α/β) (x/β)α−1 exp {− (x/β)α } , (α, β) = (1.5, 1.5).
(x−µ)2√2
exp − 2σ 2, (µ, σ ) = (0, 1.5).
2πσ  x−µ  
 
2
2
exp − s/ 1 + exp − x−µ, (µ, s) = (0, 1).
sslog x−µ)2(1√
exp − 2σ 2, (µ, σ ) = (0, 0.75).
x 2πσ

3. Half-Normal
4. Half-Logistic
5. Log-Normal
6. Paretoρλρ / (x + λ)ρ+1 , (λ, ρ) = (1, 2).

7. Burrαβ xα−1 / (1 + xα )β+1 , (α, β) = (1.5, 2.5).

8. Generalized Gamma

γ xα−1 exp {− (x/β)γ } /{β α Γ (α/γ )}, (α, β, γ ) = (5, 2, 2.5).





ρλρ
(log x−µ)2
p √ 1 exp − 2σ 2
+ (1 − p) (x+λ)
ρ+1 ,
x 2πσ

9. Log-Normal and(p, µ, σ , λ, ρ) = (0.7, 0, 0.5, 1, 2).





ρλρ
(log x−µ)2
p √ 1 exp − 2σ 2
+ (1 − p) (x+λ)
ρ+1 ,
x 2πσ

Pareto Mixture 1
10. Log-Normal and(p, µ, σ , λ, ρ) = (0.3, 0, 0.5, 1, 2).

Pareto Mixture 2

MISE-optimal c ∗ = 0.2636. Ten true distributions are considered, as listed in Table 3. All these distributions are popularly
chosen as models of, for example, the income distribution, the distribution of insurance claims and the baseline hazard.
For each distribution, 1000 data sets of sample size n = 100, 200 or 500 are simulated. All density estimates are evaluatedon an equally spaced grid of 500 points over the interval [0, 5]. As performance measures for each estimator f¯ , we computethe root integrated squared error (‘‘RISE’’) RISE f¯ (x) =





 
∞

2

f¯ (x) − f (x)

dx

1/2

and the integrated absolute bias (‘‘IAB’’)0 ∞ E f¯ (x) − f (x) dx. In our reports, the integrals are approximated over the 500 points, and the expected
IAB f¯ (x) =0value is replaced by the simulation average.

3.2. Choices of smoothing parameters
Choosing the smoothing parameter b is an important practical issue. To expedite computations, as in Hirukawa (2010),
we develop plug-in methods for TS- and JLN-MBC estimators that use gamma density as a reference.5 The plug-in smoothing
parameters for f˜TS ,MG (x) and f˜JLN ,G (x) (called ‘‘gamma-referenced smoothing parameters’’ hereafter) are defined as theminimizers of asymptotic weighted mean integrated squared errors (‘‘AWMISEs’’)b̂GR–TS = arg min AWMISE f˜TS ,MG (x)



b



 ∞

b4

λ (c )(x) wTS (x) dx + √2 π nb1/2

 ∞

g (x)



√ wTS (x) dx ,x0˜b̂GR–JLN = arg min AWMISE fJLN ,G (x)
b 

 ∞
 ∞
1
g ( x)= arg min b4
q̃2G (x) wJLN (x) dx + √
w
x
dx
,
√
(
)
JLN
b
2 π nb1/2 0x0= arg minbc 2 (1 − c )2

0

p̃2MG

and



where g (x) = xα−1 exp (−x/β) / {β α Γ (α)} is the density function for the gamma distribution with parameters (α, β), andp̃MG (x) and q̃G (x) can be obtained by replacing f (x) in pMG (x) and qG (x) with g (x). The weighting functions are chosen as

wTS (x) = x5 and wJLN (x) = x to ensure finiteness of integrals. The parameters (α, β) are replaced by their estimates α̂, β̂

via method of moments or maximum likelihood, where the latter is exclusively used throughout. Analytical expressions ofb̂GR–TS and b̂GR–JLN , as well as b̂GR–BU (= the gamma-referenced smoothing parameter for fˆMG (x)), are given in the Appendix.
We do not pursue the gamma-referenced smoothing parameter for f˜TS ,G (x); since extra terms are involved in pG (x), theminimizer of its AWMISE takes a much more complicated form than b̂GR–TS . Moreover, although it is possible to derive the5 Alternatively, b could be chosen via a version of cross-validation methods. However, this appears to be a hard problem; indeed, even for nonnegative
symmetric kernels, Jones and Signorini (1997) defer automatic bandwidth selection in this direction to future work.

118

M. Hirukawa, M. Sakudo / Computational Statistics and Data Analysis 75 (2014) 112–123

Table 4
Simulation results.
Panel (a): n = 100
Estimator

BU-G

RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: Ave
SD
IAB

BU-MG

TS-G

TS-MG

JLN-G

JLN-MG

Estimator

BU-G

RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: Ave
SD
IAB

BU-MG

TS-G

TS-MG

JLN-G

JLN-MG

Distribution 1

Distribution 2

Distribution 3

Distribution 4

Distribution 5

ROT

GR

ROT

GR

ROT

GR

ROT

GR

ROT

GR

0.0389(0.0120)0.01750.0386(0.0113)0.01800.0415(0.0110)0.01660.0445(0.0084)
0.02350.0388(0.0126)
0.0156
0.0388(0.0090)0.01170.0382(0.0123)
0.01250.0386(0.0119)0.01260.0394(0.0132)
0.00620.0408(0.0127)
0.00740.0365(0.0130)0.01040.0384(0.0104)
0.00770.0397(0.0139)0.01830.0389(0.0126)
0.01850.0434(0.0145)0.01970.0436(0.0110)0.02230.0386(0.0154)0.01790.0394(0.0109)
0.01170.0398(0.0138)
0.01270.0396(0.0131)0.01250.0414(0.0150)
0.00560.0421(0.0144)0.0067
0.0365(0.0149)
0.01290.0385(0.0120)0.00850.0350(0.0122)0.01350.0369(0.0144)0.00980.0289(0.0131)0.01100.0370(0.0151)0.01700.0318(0.0138)
0.01170.0307(0.0141)0.00470.0371(0.0130)0.01080.0397(0.0149)
0.00770.0381(0.0155)
0.00440.0431(0.0177)0.00580.0327(0.0143)0.00990.0325(0.0148)0.00420.0330(0.0115)0.01370.0336(0.0134)0.01020.0264(0.0119)0.01000.0362(0.0133)
0.02020.0282(0.0127)0.01030.0267(0.0129)0.00320.0345(0.0122)0.01050.0364(0.0140)0.0076
0.0345(0.0143)0.00370.0393(0.0168)
0.00580.0302(0.0135)
0.00750.0303(0.0139)0.00230.0710(0.0140)0.03280.0658(0.0158)
0.03830.0890(0.0118)
0.03880.0926(0.0114)0.05630.0777(0.0147)
0.03290.0849(0.0129)0.0406

0.0534(0.0153)0.01980.0469(0.0165)
0.02110.0509(0.0156)
0.01190.0450(0.0158)
0.01340.0566(0.0140)0.02100.0636(0.0133)
0.0253Distribution 6

Distribution 7

Distribution 8

Distribution 9

Distribution 10

ROT

GR

ROT

GR

ROT

GR

ROT

GR

ROT

GR

0.1287(0.0449)
0.04110.1176(0.0448)0.05520.1346(0.0498)
0.04330.1764(0.0418)
0.09670.1185(0.0535)
0.03750.1512(0.0523)
0.05380.0947(0.0321)0.02770.0836(0.0272)
0.03540.0797(0.0356)
0.02130.1171(0.0307)0.05100.0762(0.0260)
0.02190.1043(0.0235)
0.0331

0.0690(0.0201)
0.02300.0667(0.0188)
0.02690.0753(0.0173)
0.02410.0841(0.0157)
0.03650.0713(0.0215)
0.0211
0.0688(0.0159)
0.02040.0623(0.0206)
0.0152
0.0616(0.0201)
0.01740.0626(0.0224)0.00870.0636(0.0211)0.01080.0676(0.0210)0.0191
0.0663(0.0148)0.01790.0525(0.0139)0.04050.0497(0.0139)0.03740.0640(0.0124)0.05390.0572(0.0122)0.04610.0596(0.0108)0.04080.0510(0.0096)0.03380.0385(0.0136)0.01590.0379(0.0133)
0.01470.0405(0.0143)0.00590.0407(0.0142)0.00520.0366(0.0117)0.0186
0.0338(0.0106)
0.01630.0663(0.0233)
0.03380.0655(0.0234)
0.03930.0763(0.0231)0.04160.0795(0.0272)
0.05350.0755(0.0279)
0.03590.0700(0.0271)
0.0361

0.0549(0.0168)
0.0237
0.0560(0.0179)
0.02670.0529(0.0171)0.01780.0558(0.0199)
0.0211
0.0619(0.0160)
0.0276
0.0595(0.0173)
0.02690.0756(0.0303)
0.03280.0698(0.0275)
0.03200.0720(0.0337)
0.03220.1000(0.0319)0.06030.0666(0.0345)0.02620.0789(0.0350)
0.02520.0603(0.0191)
0.02170.0581(0.0172)
0.02240.0516(0.0190)0.01540.0711(0.0221)0.02660.0506(0.0131)0.01730.0624(0.0117)0.0204(continued on next page)gamma-referenced smoothing parameter for f˜JLN ,MG (x) in a similar way,6 our preliminary simulation results indicate thatthe formula frequently generates large values and thus we do not advocate its use. From the viewpoint of practical relevance,b̂GR–TS and b̂GR–JLN are simply employed for f˜TS ,G (x) and f˜JLN ,MG (x) in our simulations, respectively. Similarly, b̂GR–BU is chosen

as the smoothing parameter for fˆG (x).
Besides, a very simple formula is frequently applied in the literature (e.g. Gustafsson et al., 2009). In this respect,a ‘‘rule-of-thumb’’ smoothing parameter is also considered for each estimator. More precisely, we additionally employb̂ROT –BU = σ̂x n−2/5 for BU-G, BU-MG and b̂ROT –MBC = σ̂x n−2/9 for four MBC estimators, where σ̂x is the sample standard
deviation of observations.
3.3. Simulation results
Table 4 presents Monte Carlo results. For each distribution, results are qualitatively similar across three sample sizes,
and values of performance measures decrease with the sample size. As regards smoothing parameters, except average
RISEs for Distributions 3 and 4, the gamma-referenced method in general works better than the rule-of-thumb method,

6 Choosing x5 as the weighting function, we can derive the smoothing parameter asarg min b4
b

 ∞
0

q̃2MG (x) x5 dx +

√

1

2 π

nb1/2

 ∞
0

g (x ) 5√ x dxx=

4α β 9/2 Γ (α + 9/2) Γ (α)√

4 π (α − 1) (α − 2) Γ (2α)
2

2

2/9

n−2/9 .

M. Hirukawa, M. Sakudo / Computational Statistics and Data Analysis 75 (2014) 112–123

119

Table 4 (continued)
Panel (b): n = 200
Estimator

BU-G

BU-MG

TS-G

TS-MG

JLN-G

JLN-MG

RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB

Estimator

BU-G

BU-MG

TS-G

TS-MG

JLN-G

JLN-MG

RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IABDistribution 1

Distribution 2

Distribution 3

Distribution 4

Distribution 5

ROT

GR

ROT

GR

ROT

GR

ROT

GR

ROT

GR

0.0319(0.0094)
0.01400.0309(0.0086)
0.01410.0366(0.0087)
0.01440.0375(0.0057)
0.01940.0328(0.0097)
0.01290.0332(0.0064)0.00990.0306(0.0097)
0.00990.0303(0.0091)0.00980.0309(0.0102)0.00520.0317(0.0093)
0.00620.0295(0.0100)
0.00840.0315(0.0075)
0.00650.0319(0.0103)0.01460.0302(0.0093)0.01440.0372(0.0109)0.01710.0357(0.0078)0.01890.0312(0.0111)0.01450.0319(0.0073)
0.0096

0.0314(0.0100)
0.01010.0305(0.0093)0.00970.0313(0.0104)0.00480.0314(0.0098)
0.00550.0286(0.0105)
0.01020.0302(0.0081)0.00690.0278(0.0092)
0.01050.0288(0.0106)
0.00760.0225(0.0100)
0.00930.0284(0.0113)0.01340.0244(0.0104)
0.00920.0234(0.0105)
0.00380.0291(0.0096)0.00840.0309(0.0108)
0.00600.0285(0.0111)0.00350.0319(0.0126)
0.00440.0249(0.0106)
0.00780.0246(0.0108)
0.00330.0261(0.0087)
0.0106
0.0263(0.0100)
0.00790.0205(0.0092)
0.00840.0280(0.0101)
0.0161
0.0217(0.0097)
0.00800.0205(0.0097)
0.00260.0271(0.0091)0.00820.0283(0.0103)
0.00590.0258(0.0104)
0.00310.0292(0.0121)0.00430.0229(0.0101)
0.00570.0230(0.0103)0.00190.0592(0.0109)0.02720.0523(0.0123)0.03060.0820(0.0094)
0.03510.0837(0.0090)0.05000.0686(0.0114)0.02840.0766(0.0101)
0.0361

0.0424(0.0112)0.01590.0362(0.0121)0.01640.0403(0.0110)0.01020.0339(0.0111)0.0106
0.0468(0.0102)0.0172
0.0534(0.0100)
0.0214

Distribution 6

Distribution 7

Distribution 8

Distribution 9

Distribution 10

ROT

GR

ROT

GR

ROT

GR

ROT

GR

ROT

GR

0.1155(0.0391)
0.03660.1018(0.0401)
0.04820.1298(0.0445)0.04170.1706(0.0395)
0.09420.1121(0.0484)
0.03510.1453(0.0481)0.05200.0780(0.0236)
0.02270.0649(0.0194)
0.02790.0677(0.0275)
0.01840.1041(0.0237)0.0451
0.0638(0.0197)
0.01800.0905(0.0178)
0.02800.0572(0.0156)
0.01890.0539(0.0143)
0.02170.0681(0.0133)
0.02150.0725(0.0110)0.03100.0622(0.0166)
0.01800.0601(0.0110)0.01760.0498(0.0156)
0.01220.0478(0.0148)0.01360.0495(0.0163)
0.00740.0489(0.0147)
0.00890.0580(0.0161)
0.01580.0574(0.0103)
0.01540.0425(0.0104)
0.03250.0403(0.0105)
0.03020.0558(0.0091)0.04720.0502(0.0091)0.04090.0505(0.0079)
0.03430.0428(0.0068)
0.02830.0301(0.0098)
0.01250.0297(0.0097)
0.01170.0301(0.0101)
0.00480.0303(0.0100)
0.00440.0289(0.0085)0.01500.0265(0.0073)
0.01310.0571(0.0184)
0.02980.0551(0.0181)0.03380.0712(0.0191)0.03960.0728(0.0235)0.0501
0.0694(0.0238)0.0341
0.0632(0.0220)0.03380.0445(0.0124)0.01960.0442(0.0135)
0.02130.0426(0.0127)0.01580.0432(0.0144)0.01750.0533(0.0116)0.02430.0502(0.0130)
0.02330.0668(0.0250)0.02910.0598(0.0221)0.02870.0677(0.0297)
0.0311
0.0951(0.0295)
0.05970.0609(0.0308)0.02510.0743(0.0313)0.02480.0492(0.0139)0.01750.0458(0.0128)0.0181
0.0422(0.0133)
0.01390.0605(0.0156)
0.02330.0427(0.0096)0.01550.0552(0.0079)
0.0194(continued on next page)and improvement in performance measures by the former is often substantial. In particular, a smaller IAB by the gammareferenced method can be attributed to the fact that on average it yields a smaller smoothing parameter value (although theresults are unreported). Therefore, we may safely concentrate on the results from the gamma-referenced method hereafter.
Table 4 reveals that overall MBC works. For each combination of sample size and distribution, MBC estimators at least
tend to deliver smaller IABs than two BU ones. The results suggest that there is no uniformly superior MBC estimator; rather,which MBC estimator performs best depends on distributions.
It is hard to judge whether for a given MBC estimator, the MG kernel improves performance measures over the G kernel,
whereas it appears to be advantageous to employ the former when no MBC is made. In particular, TS-MG often performs mostpoorly among all six estimators. Its inferior performance can be attributed to the following two respects. First, TS-MBC estimation relies on two smoothing parameters b and b/c. Controlling both b and b/c is a cumbersome task. Because 0 < c < 1,
the density estimator using b/c tends to be oversmoothed, which is potentially a source of a large bias in every TS-MBC estimator. On the other hand, if we make b too short in order to have a reasonable length of b/c, additional variability is introduced to the other estimator using b due to undersmoothing. Second, when the MG kernel is employed, these two smoothing
parameters also play a role of determining the boundary region explicitly (e.g. [0, 2b) for the density estimator using b). Unless b is short enough, there is a relatively small interior region for the density estimator using b/c. This aspect is also thought
to worsen the performance measures of TS-MG, in particular for small to medium sample sizes. In conclusion, when TS-MBC
estimation is applied, it is desirable to use the smoothing parameter choice method that tends to provide a small value consistently. Our experiment indicates that for a given distribution and a given sample size, b̂GR–TS on average yields a smaller
value than b̂ROT –MBC does, which explains why the former is preferable in general for TS-MBC estimation over the latter.

120

M. Hirukawa, M. Sakudo / Computational Statistics and Data Analysis 75 (2014) 112–123

Table 4 (continued)
Panel (c): n = 500
Estimator

BU-G

BU-MG

TS-G

TS-MG

JLN-G

JLN-MG

RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB

Estimator

BU-G

BU-MG

TS-G

TS-MG

JLN-G

JLN-MG

RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IAB
RISE: AveSD
IABDistribution 1

Distribution 2

Distribution 3

Distribution 4

Distribution 5

ROT

GR

ROT

GR

ROT

GR

ROT

GR

ROT

GR

0.0247(0.0066)
0.01030.0230(0.0058)
0.01030.0316(0.0065)
0.01190.0307(0.0036)0.01520.0266(0.0071)
0.00980.0276(0.0040)0.00800.0231(0.0066)
0.00730.0219(0.0060)
0.00710.0230(0.0070)0.00430.0231(0.0059)
0.00490.0228(0.0069)
0.00630.0247(0.0047)0.00530.0238(0.0072)
0.01080.0215(0.0061)0.01000.0307(0.0077)0.01420.0278(0.0052)
0.01500.0239(0.0077)
0.01110.0244(0.0045)
0.00740.0229(0.0068)
0.00750.0216(0.0061)
0.00660.0219(0.0069)0.0041
0.0218(0.0062)0.00390.0210(0.0072)
0.00770.0223(0.0050)
0.00520.0202(0.0060)
0.00750.0209(0.0069)
0.00560.0160(0.0065)
0.00720.0201(0.0071)
0.00970.0171(0.0065)
0.00660.0164(0.0067)
0.00280.0212(0.0062)0.00600.0223(0.0070)0.00440.0195(0.0069)
0.00240.0218(0.0079)
0.00310.0173(0.0066)
0.00540.0172(0.0069)
0.00240.0190(0.0058)
0.00760.0191(0.0066)
0.00580.0146(0.0062)0.00660.0199(0.0065)
0.01160.0152(0.0063)0.00560.0144(0.0063)
0.00200.0197(0.0060)
0.00580.0205(0.0067)
0.00440.0177(0.0065)
0.00220.0200(0.0076)
0.00310.0159(0.0063)
0.00380.0161(0.0066)0.00140.0457(0.0079)
0.02080.0380(0.0091)
0.02220.0731(0.0070)0.03050.0727(0.0070)0.04220.0576(0.0080)0.02320.0658(0.0078)0.03030.0313(0.0080)
0.01190.0263(0.0087)
0.01190.0304(0.0078)0.00840.0234(0.0079)
0.0076
0.0366(0.0070)0.01330.0421(0.0075)
0.0168Distribution 6

Distribution 7

Distribution 8

Distribution 9

Distribution 10

ROT

GR

ROT

GR

ROT

GR

ROT

GR

ROT

GR

0.0983(0.0326)0.03080.0824(0.0335)
0.03930.1228(0.0394)
0.03910.1629(0.0368)
0.08960.1040(0.0428)0.03250.1370(0.0438)
0.0486

0.0587(0.0154)0.01700.0452(0.0115)0.0196
0.0544(0.0183)0.01470.0889(0.0164)0.03720.0499(0.0133)
0.01380.0740(0.0115)
0.02170.0443(0.0111)0.01430.0394(0.0101)
0.01570.0607(0.0097)
0.01830.0599(0.0074)
0.0246
0.0526(0.0119)0.01460.0508(0.0073)
0.01420.0369(0.0110)0.00900.0336(0.0099)0.00930.0370(0.0116)0.00600.0352(0.0093)
0.00670.0479(0.0115)0.01250.0475(0.0068)
0.01220.0313(0.0072)
0.02370.0297(0.0072)
0.02200.0459(0.0063)0.03890.0414(0.0061)0.03420.0400(0.0055)
0.02700.0332(0.0044)
0.02190.0217(0.0065)
0.00880.0214(0.0065)
0.00820.0207(0.0066)
0.00330.0208(0.0066)
0.00300.0212(0.0057)
0.01130.0190(0.0049)
0.00940.0461(0.0141)0.02430.0439(0.0137)
0.02720.0651(0.0155)
0.03650.0654(0.0188)0.0452
0.0619(0.0195)
0.03100.0560(0.0170)0.03040.0330(0.0090)0.01450.0323(0.0091)0.01550.0324(0.0088)0.01270.0313(0.0094)
0.0136
0.0430(0.0080)0.01980.0403(0.0085)0.0191

0.0540(0.0188)
0.02280.0472(0.0162)0.02320.0598(0.0239)
0.02690.0868(0.0249)
0.0541
0.0524(0.0244)
0.02190.0662(0.0249)
0.02270.0370(0.0097)
0.0126
0.0332(0.0089)
0.01330.0336(0.0090)0.01170.0499(0.0110)0.01960.0352(0.0068)
0.01330.0474(0.0053)0.0175

Note: ‘‘ROT’’ and ‘‘GR’’ in column headings denote ‘‘rule-of-thumb’’ and ‘‘gamma-referenced’’ smoothing parameter choice methods. ‘‘Ave’’ and ‘‘SD’’ aresimulation averages and standard deviations of RISEs.
Table 5
Summary statistics of real data.
Variable

Mean

SD

Min

Max

Percentiles
25%

50%

75%

Per Capita Income (in USD)
Belgian Female Wage (in BEF)
US Male Wage (in USD)

3,790
414
958

4,156
154
404

224
88
115

25,646
1,364
3,078

744
314
668

2,005
384
905

5,817
483
1,160

4. Empirical illustrations
This section applies the MBC density estimators to real data sets. Our focus is on income and wage data. They have anatural boundary at the origin, and their distributions are empirically characterized by concentration of observations near
the boundary and a long tail with sparse data. Hence, we are motivated to estimate the densities using asymmetric kernels.
The data we use include: (a) per capita income of 114 countries (in US dollars); (b) gross hourly wage rate for 579
Belgian females (in Belgian francs); and (c) monthly earnings of 935 US males (in US dollars). The first and third data are
originally used in Romer (1993) and Blackburn and Neumark (1992), and these data sets are now available under the names‘‘openness’’ and ‘‘wage2’’ as supplemental materials for Wooldridge (2013). The second is taken from the data set ‘‘bwage’’,a supplemental material for Verbeek (2012). Table 5 reports summary statistics of the data.
The density of each data is estimated by BU-G, TS-G and JLN-G. The data are first converted to the scale-adjusted ones
by dividing them by either 103 or 104 , and then the resulting density estimates are back-transformed to the ones in theM. Hirukawa, M. Sakudo / Computational Statistics and Data Analysis 75 (2014) 112–123

4.0

1.00

3.5
0.00

Per Capita Income (x 10^4 U.S. Dollars)

0.75

Density

0.25
0.25

0.50

0.75

1.00

1.25

1.50

Wage (x 10^3 Belgian Francs)

0.00

0.0

0.5

0.5
0.0
0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00

G
TS–G
JLN–G

0.50

1.5

2.0

2.5

G
TS–G
JLN–G

1.0

Density

3.0

3.0
2.5
2.0
1.5
1.0

Density

G
TS–G
JLN–G

Panel (c): U.S. Male
Wage (n = 935)

1.25

Panel (b): Belgian
Female Wage (n = 579)

3.5

Panel (a): Per Capita
Income (n = 114)

121

0.0

0.5

1.0

1.5

2.0

2.5

3.0

Wage (x 10^3 U.S. Dollars)Fig. 1. Density estimates from real data.

original scale.7 Plug-in smoothing parameters b̂GR–BU , b̂GR–TS and b̂GR–JLN are employed for BU-G, TS-G and JLN-G, respec-tively. Their values (based on the scale-adjusted data) are: b̂GR–BU , b̂GR–TS , b̂GR–JLN = (0.0434, 0.0655, 0.1752) for (a);(0.0041, 0.0054, 0.0395) for (b); and (0.0105, 0.0152, 0.0677) for (c). The plots of density estimates are presented in Fig. 1.
We can observe that while estimates from BU-G and TS-G largely look alike, the one from JLN-G is considerably different.
Hirukawa (2010) reports the tendency of JLN-MBC that when the BU estimator underestimates (overestimates) the density,
its corresponding JLN-MBC estimator corrects the estimate in an upward (downward) direction. Fig. 1 indicates that JLN-G indeed exhibits such smooth-out tendencies, and as a result it substantially smooths away the part that BU-G estimates wiggly.
5. Concluding remarks
This paper has demonstrated that two MBC techniques studied in Hirukawa (2010) for density estimation with support on[0, 1] can be extended to density estimation with support on [0, ∞) using asymmetric kernels. Under sufficient smoothness of the true density, both bias reduction methods are shown to improve the order of magnitude in bias from O (b) to O b2 ,while the order of magnitude in variance remains unchanged.
Two classes of MBC density estimators are by construction
nonnegative, and establish a faster convergence rate of O n−8/9 in MSE for the interior part when best implemented. Monte
Carlo simulations confirm bias reduction via two MBC methods.
While this paper deals exclusively with the univariate case, it appears to be possible to extend the analysis to joint density estimation. In the multivariate (‘‘MV’’) case, the MV-BU estimator of the density f with support on Rd+ using a random

nsample {Xi }ni=1 = (X1i , . . . , Xdi )⊤ i=1 is fˆj,b (x) = n−1

(Xi ), where Kj(x,b) (u) := dℓ=1 Kj(xℓ ,b) (uℓ ) is the productasymmetric kernel given a design point x = (x1 , . . ., xd )⊤ and a smoothing
parameter b1 = · · · = bd = b. If f has four cond(rj +1/2)+2
tinuous and bounded partial derivatives and b + 1/ nb→ 0, then it can be shown that each of the MV-TS-MBC
1/(1−c ) −c /(1−c )estimator f˜TS ,j (x) = fˆj,b (x)fˆj,b/c (x)for some c ∈ (0, 1) and the MV-JLN-MBC estimatorn1  Kj(x,b) (Xi )f˜JLN ,j (x) = fˆj,b (x)n i=1 fˆj,b (Xi )dadmits the expansion f (x)+ O b2 + n−1/2 ℓ=1 b−(1/2)(1/2+rj 1ℓ ) , where 1ℓ := 1 (xℓ /b → κℓ > 0) is the indicator functionni=1 Kj(x,b)

that takes 1 if xℓ lies in the boundary region. Analytical expressions of leading bias and variance terms and their finite sample
performance will be addressed in a separate paper.
Acknowledgements
We would like to thank the Associate Editor, two anonymous referees, Benedikt Funke, Matthias Hagmann, Susumu Imai,James MacKinnon and seminar participants at Queen’s University and Development Bank of Japan for insightful comments
and suggestions. The first author gratefully acknowledges financial support from Japan Society for the Promotion of Science(grant number 23530259). The views expressed herein and those of the authors do not necessarily reflect the views of theDevelopment Bank of Japan.

7 See the horizontal axis of each panel in Fig. 1 for details.

122

M. Hirukawa, M. Sakudo / Computational Statistics and Data Analysis 75 (2014) 112–123

Appendix A. Sketches of the proofs
The proof of each theorem requires kernel-specific arguments, which include Taylor expansions and properties of therandom variable corresponding to the kernel. To save space, we provide only outline proofs when the Gamma kernel is
employed. Full-length proofs are available upon request.
Outline Proof of Theorem 1. For the bias, it follows from the proof of Theorem 1 in Hirukawa (2010) and Assumption 2 that
1

 
1

a21,G (x, f )c (1 − c )

2

f (x)




E f˜TS ,G (x) = f (x) +




 
− a2,G (x, f ) b2 + o b2 .

For the variance, recognize that



1



Var f˜TS ,G (x) =(1 − c )

2













Var fˆG,b (x) − 2cCov fˆG,b (x) , fˆG,b/c (x) + c 2 Var fˆG,b/c (x)Then, approximate each term as in the proof of Theorem 1 in Hirukawa (2010).

 + O n− 1 .

Outline Proof of Theorem 2. For the bias, write h (x) = a1,G (x, f ) /f (x). It follows from the procedures in Section A.2.1 of
Hirukawa (2010, pp. 490–491) that







E f˜JLN ,G (x) = f (x) − f (x) h′ (x) +


 
 
x ′′
h (x) b2 + o b2 = f (x) − f (x) a1,G (x, h) b2 + o b2 .
2

For the variance, the procedures in Section A.2.2 of Hirukawa (2010, p. 492) imply that

1  KG(x/b+1,b) (Xi )nf˜JLN ,G (x) ∼ f (x)2−

f (Xi )n i=1

fˆG,b (Xi ).

f (Xi )

By Stirling’s approximation to the gamma function inside the integral and the trimming argument in Chen (2000), the leading
variance on the right-hand side can be approximated by n−1 Ab (x) f (x), where −1/2
b


√ 1/2

−1
b Γ (2x/b + 1)
2 πx∼
Ab (x) = 2x/b+1 2
 Γ (2κ + 1) b−12Γ (x/b + 1)
22κ+1 Γ 2 (κ + 1)
Then, the result immediately follows.

if x/b → ∞if x/b → κ.

Appendix B. Formulae for gamma-referenced smoothing parameters
The analytical expression of b̂GR–TS is2/9 4α β 9/2 Γ (α + 9/2) Γ (α) 2/9 −2/9b̂GR–TS = c (1 − c ) λ (c )n,√16 π CTS (α) Γ (2α)
 2

2

where
CTS (α) =

1
36

+
−

2





3
3
1
1
(α − 2)2 α −(α − 1)2 − (α − 2) α −(α − 1)2 α α +

1
9
1
3

2



(α − 2) α −

3

6



2



(α − 1) α α +1(α − 1) α α +2(α + 1) α +

2

1



2
3
2



(α + 1) +

2

1
4
1



(α − 1)2 α α +

(α + 2) + α α +
9

1
2



 α

4 β 5/2 Γ (α + 1/2) Γ (α)

2/9√4 π Γ (2α)

n−2/9 .

Moreover, the gamma-referenced smoothing parameter for fˆMG (x) is defined asb̂GR–BU = arg min AWMISE fˆMG (x)b 2 ∞= arg min
b

b

4

0

x2 g ′′ (x)



2

1

w (x) dx + √2 π nb1/2

 ∞g ( x)

0

x

2





3
(α + 1) α +2(α + 1) α +On the other hand, b̂GR–JLN takes a much simpler form. It is given byb̂GR–JLN =1√ w (x) dx ,

3
2





5
.
(α + 2) α +2

M. Hirukawa, M. Sakudo / Computational Statistics and Data Analysis 75 (2014) 112–123

123

where the weighting function w (x) is chosen as w (x) = x3 to ensure finiteness of integrals. It follows that b̂GR–BU can be
expressed as αb̂GR–BU =4 β 5/2 Γ (α + 5/2) Γ (α)√

8 π CBU (α) Γ (2α)

2/5

n−2/5 ,

where
CBU (α) =



1
1
2
2
2
−
2
−
1
−
1
3
α
−
4
α
α
+
−
−
2
−
1
α
+
(α
) (α
)
(α
)(
)
(α
) (α
)
4
2
2






1
1
1
3
− (α − 1) α α +. (α + 1) + α α +(α + 1) α +
1

2

4

2

2

References
Blackburn, M., Neumark, D., 1992. Unobserved ability, efficiency wages, and interindustry wage differentials. Quart. J. Econ. 107, 1421–1436.
Chen, S.X., 1999. Beta kernel estimators for density functions. Comput. Statist. Data Anal. 31, 131–145.
Chen, S.X., 2000. Probability density function estimation using gamma kernels. Ann. Inst. Statist. Math. 52, 471–480.
Ferraty, F., Vieu, P., 2006. Nonparametric Functional Data Analysis: Theory and Practice. Springer, New York, NY.
Gustafsson, J., Hagmann, M., Nielsen, J.P., Scaillet, O., 2009. Local transformation kernel density estimation of loss distributions. J. Bus. Econom. Statist. 27,
161–175.
Hagmann, M., Scaillet, O., 2007. Local multiplicative bias correction for asymmetric kernel density estimators. J. Econometrics 141, 213–249.
Hirukawa, M., 2010. Nonparametric multiplicative bias correction for kernel-type density estimation on the unit interval. Comput. Statist. Data Anal. 54,
473–495.
Jin, X., Kawczak, J., 2003. Birnbaum–Saunders and lognormal kernel estimators for modelling durations in high frequency financial data. Ann. Econ. Finance
4, 103–124.
Jones, M.C., Linton, O., Nielsen, J.P., 1995. A simple bias reduction method for density estimation. Biometrika 82, 327–338.
Jones, M.C., Signorini, D.F., 1997. A comparison of higher-order bias kernel density estimators. J. Amer. Statist. Assoc. 92, 1063–1073.
Karunamuni, R.J., Alberts, T., 2005. On boundary correction in kernel density estimation. Stat. Methodol. 2, 191–212.
Leblanc, A., 2010. A bias-reduced approach to density estimation using Bernstein polynomials. J. Nonparametr. Stat. 22, 459–475.
Romer, D., 1993. Openness and inflation: theory and evidence. Quart. J. Econ. 108, 869–903.
Scaillet, O., 2004. Density estimation using inverse and reciprocal inverse Gaussian kernels. J. Nonparametr. Stat. 16, 217–226.
Stone, C.J., 1980. Optimal rates of convergence for nonparametric estimators. Ann. Statist. 8, 1348–1360.
Stuetzle, W., Mittal, Y., 1979. Some comments on the asymptotic behavior of robust smoothers. In: Gasser, T., Rosenblatt, M. (Eds.), Smoothing Techniques
for Curve Estimation: Proceedings of a Workshop Held in Heidelberg, 2–4 April, 1979. Springer-Verlag, Berlin, pp. 191–195.
Terrell, G.R., Scott, D.W., 1980. On improving convergence rates for nonnegative kernel density estimators. Ann. Statist. 8, 1160–1163.
Verbeek, M., 2012. A Guide to Modern Econometrics, fourth ed. John Wiley & Sons, Chichester, UK.
Wooldridge, J.M., 2013. Introductory Econometrics: A Modern Approach, fifth ed. South-Western Cengage Learning, Mason, OH.