Econometrics and Statistics 9 (2019) 156–170

Contents lists available at ScienceDirect

Econometrics and Statistics
journal homepage: www.elsevier.com/locate/ecosta

Nonparametric estimation and testing on discontinuity ofpositive supported densities: a kernel truncation approach
Benedikt Funke a, Masayuki Hirukawa b,∗
ab

Department of Mathematics, Technical University of Dortmund, Vogelpothsweg 87, Dortmund 44227, Germany
Faculty of Economics, Setsunan University, 17-8 Ikeda Nakamachi, Neyagawa, Osaka 572–8508, Japana r t i c l e

i n f oArticle history:
Received 19 October 2016
Revised 27 July 2017
Accepted 27 July 2017
Available online 3 August 2017
JEL Classiﬁcation:
C12
C13
C14
MSC:
62G07
62G10
62G20
Keywords:
Density estimation
Discontinuous probability density
Gamma kernel
Incomplete gamma functions
Nonparametric kernel testing
Regression discontinuity designa b s t r a c t
Discontinuity in density functions is of economic importance and interest. For instance,in studies on regression discontinuity designs, discontinuity in the density of a runningvariable suggests violation of the no-manipulation assumption. In line with this notion,
estimation and testing procedures on discontinuity in densities with positive support aredeveloped. The proposed approach is built on splitting the asymmetric, gamma kernel into
two parts at a prespeciﬁed cutoff that is suspected to be a discontinuity point and constructing two truncated kernels. The jump-size magnitude of the density at the cutoff
can be estimated nonparametrically by two kernels and a multiplicative bias correction
method. The estimator is easy to implement, and its convergence properties are deliveredby various approximation techniques on incomplete gamma functions. Based on the jumpsize estimator, two versions of test statistics for the null of continuity at a given cutoff arealso proposed. Moreover, estimation theory of the entire density in the presence of a discontinuity point is explored. Monte Carlo simulations conﬁrm nice ﬁnite-sample propertiesof the jump-size estimator and the test statistics.
© 2017 EcoSta Econometrics and Statistics. Published by Elsevier B.V. All rights reserved.

1. Introduction
The objective of this paper is to develop new estimation and testing procedures of discontinuity in density functions
with support on R+ . Inference on possibly discontinuous densities has been explored in nonparametric statistics: examples
include Liebscher (1990), Cline and Hart (1991), and Chu and Cheng (1996), to name a few. Discontinuity in densities is alsoof economic importance and interest. Local randomization of a continuous running variable is a key requirement for thevalidity of regression discontinuity designs (‘RDD’); if the value of the running variable falls into the left and right of thecutoff strategically, then treatment effects are no longer point identiﬁed due to self-selection. Therefore, detection of discontinuity in the density of the running variable at the cutoff suggests evidence of such strategic behavior or manipulationin RDD. Nonetheless, estimation and inference on discontinuity in densities have not attracted interest in econometrics up∗

Corresponding author.
E-mail addresses: benedikt.funke@mathematik.tu-dortmund.de (B. Funke), hirukawa@econ.setsunan.ac.jp (M. Hirukawa).

http://dx.doi.org/10.1016/j.ecosta.2017.07.006
2452-3062/© 2017 EcoSta Econometrics and Statistics. Published by Elsevier B.V. All rights reserved.

B. Funke, M. Hirukawa / Econometrics and Statistics 9 (2019) 156–170

157

until recently. McCrary (2008) applies a bin-based local linear regression method to estimate the jump-size magnitude of adensity at a prespeciﬁed cutoff. Subsequently, Otsu et al. (2013) propose two versions of empirical likelihood-based inference
procedures grounded on binning and local likelihood methods. While our proposal can be viewed as an extension of thesearticles, it has a unique feature. In our approach, jump sizes are estimated by means of density estimation techniques using
the kernels obtained through truncating asymmetric kernels at a given location that is suspected to be a discontinuity point,
unlike nonparametric regression or local likelihood approaches using standard symmetric kernels.
Before proceeding, it is worth explaining why we specialize in asymmetric kernel smoothing. Empirical studies on discontinuity in densities frequently pay attention to the distributions of economic variables such as (taxable or relative) incomes (Saez, 2010; Bertrand et al., 2015), wages (DiNardo et al., 1996), school enrollment counts (Angrist and Lavy, 1999)
and proportion of votes for proposed bills (McCrary, 2008). The distributions, if they are free of discontinuity points, canbe empirically characterized by two stylized facts, namely, (i) existence of a lower bound in support (most possibly at theorigin) and (ii) concentration of observations near the boundary and a long tail with sparse data. When estimating suchdensities nonparametrically using symmetric kernels, we must rely either on a boundary correction method and an adaptive
smoothing technique (e.g., variable bandwidth methods) simultaneously, or on back-transforming the density estimator fromthe log-transformed data to the original scale. The former is apparently cumbersome, and density estimates by the latteroften behave poorly (e.g., Cowell et al., 1998) although the method is popularly applied in empirical works.
Asymmetric kernels with support on R+ have emerged as a viable alternative that can accommodate the stylized facts.
Our motivation to employ asymmetric kernels can be discussed in the context of RDD. Existence of a boundary in the support of the density of the running variable may not be a serious issue unless the cutoff is near the boundary. Rather, we use
the adaptive smoothing property of asymmetric kernels as best as possible for our purposes. This capacity originates fromthe fact that their shapes vary across design points where smoothing is made. The property is particularly advantageous
for capturing and inspecting the shape of the entire density of the positively-valued running variable that is assumed to beright-skewed. Furthermore, this study focuses exclusively on the gamma kernel by Chen (20 0 0)

KG(x,b) (u ) =ux/b exp (−u/b)1(u ≥ 0 ),bx/b+1 (x/b + 1 )where x(≥ 0) and b(> 0) are the design point and smoothing parameter, respectively. Not only is the gamma kernel mostpopular among all asymmetric kernels, but also there is rich literature on approximations to the gamma and incomplete
gamma functions, which are key ingredients for our asymptotic analyses.
When the density has a discontinuity point, the jump-size magnitude at the point can be deﬁned as the difference
between left and right limits of the density at the point. While nonparametric regression (McCrary, 2008) and empirical
likelihood (Otsu et al., 2013) methods have been applied to estimate the jump size, we attempt to have our jump-size estimator preserve appealing properties of the gamma kernel. Accordingly, we split the gamma kernel into two parts at a given
cutoff where the density is suspected to be discontinuous, and make each part a legitimate kernel by re-normalization. The
left and right limits of the density can be estimated by two truncated kernels. Although the estimators are consistent andtheir variance convergences are usual O n−1 b−1/2 where n is the sample size, their bias convergences are O(b1/2 ), not theusual O(b) rate. Then, we apply the multiplicative bias correction technique by Terrell and Scott (1980) to eliminate the undesirable O(b1/2 ) biases without inﬂating the order of magnitude in variance. Moreover, we take particular care of choosing
the smoothing parameter. Speciﬁcally, the method of power-optimality smoothing parameter selection by Kulasekera andWang (1998) is tailored to inference problems on discontinuous densities.
Some readers may wonder why we avoid estimating the left and right limits of the density at the prespeciﬁed cutoff
using the beta (Chen, 1999) and gamma kernels, respectively. The reason is bias and variance convergences of the beta andgamma estimates when the cutoff is chosen as the boundary. In the vicinity
 of the cutoff, their bias convergences are usual
O(b). On the other hand, the variance convergences slow down to O n−1 b−1 . Although the inferior rate is likely to adversely

affect power properties of our test, it is hard (or even impossible) to improve the rate to usual O n−1 b−1/2 . In contrast, den1/2
sity estimates −1 −1/2 using the truncated kernels have the inferior bias convergence of O(b ) with the usual variance convergenceof O n b
maintained. While the slow bias convergence potentially has a negative impact on size properties of thetest, it can be improved to usual O(b) with no additional conditions (e.g., extra smoothness in the density), by means of amultiplicative bias correction method.
Our proposal has three contributions to the literature. First, unlike the binned local linear (‘BLL’) estimation by McCrary(2008), our kernel truncation approach always generates nonnegative density estimates and is free from choosing bin widths.
Our jump-size estimator is also easy to implement. Since it has a closed form, nonlinear optimization as in Otsu et al.
(2013) is unnecessary. While incomplete gamma functions are key components of our estimation and testing procedures,standard statistical packages including GAUSS, Matlab and R prepare a command that can return values of the functions
either directly or in the form of gamma cumulative distribution functions.
Second, in delivering convergence results of asymmetric kernel estimators, we utilize mathematical tools and proof
strategies that are totally different from those for nonparametric estimators smoothed by symmetric kernels. Asymptotic
results throughout this paper are built upon a few different approximation techniques on incomplete gamma functions;such proof strategies are taken for the ﬁrst time in the econometric literature, to the best of our knowledge.

158

B. Funke, M. Hirukawa / Econometrics and Statistics 9 (2019) 156–170

Third, we also present estimation theory of the entire density in the presence of a discontinuity point. Indeed, Imbens
and Lemieux (2008) argue importance of graphical analyses in empirical studies on RDD, including inspections of densities
of running variables. It is demonstrated that density estimators smoothed by the truncated gamma kernels admit the same
bias and variance approximations as the gamma kernel density estimator does. Furthermore, the truncated gamma-kernel
density estimator is shown to be consistent even when the true density is unbounded at the origin.
The remainder of this paper is organized as follows. Section 2 presents estimation and testing procedures of the density
at a prespeciﬁed location that is suspected to be a discontinuity point. As an important practical problem, a smoothing parameter selection method is also developed. Our particular focus is on the choice method for power optimality. In Section 3,
we discuss how to estimate the entire density when the density has a discontinuity point. Convergence properties of density estimates are also explored. Section 4 conducts Monte Carlo simulations to evaluate ﬁnite-sample properties of the
proposed jump-size estimator and test statistic. An empirical application on the validity of RDD is presented in Section 5.
Section 6 summarizes the main results of the paper. All technical proofs and additional simulation results are made available
as a supplement, which is published online at the journal webpage.
∞This paper adopts the following notational conventions: for a > 0, (a ) = 0 t a−1 exp (−t )dt is the gamma function; for z a−1
 ∞ a−1a, z > 0, γ (a, z ) = 0 t
exp (−t )dt and (a, z ) = z t
exp (−t )dt = (a ) − γ (a, z ) denote the lower and upper incomplete gamma functions, respectively; 1{·} signiﬁes an indicator function; and · denotes the integer part.
2. Estimation and inference for discontinuity in the density
2.1. Setup
Suppose that we suspect discontinuity of the probability density function (‘pdf’) f(x) at a given point x = c (> 0 ), whichis assumed to be interior throughout. Also letf− (c ) := lim f (x ) and f+ (c ) := lim f (x )
x↑c

x↓c

be the lower and upper limits of the pdf at x = c, respectively. Our parameter of interest is the jump-size magnitude of thedensity at c

J (c ) := f+ (c ) − f− (c ).
To check whether f is (dis)continuous at c, we ﬁrst estimate J(c) nonparametrically and then proceed to a hypothesis testingfor the null of continuity of f at c, i.e., H0 : J (c ) = 0, against the two-sided alternative.
2.2. An issue in estimating two limits of the density
To develop a consistent estimator of J(c), we start our analysis from estimating two limits of the density at c. Let {Xi }ni=1
be a univariate random sample drawn from a distribution that has the pdf 
f. Whenat c, a rea f is indeed discontinuous
 
sonable method would be to estimate f− (c ) and f+ (c ) using sub-samples Xi− := {Xi : Xi < c} and Xi+ := {Xi : Xi ≥ c},respectively. Instead of relying on nonparametric regression or local likelihood methods, we split the gamma kernel into
two parts at c, namely,

KG(x,b) (u ) := KGL (x,b;c ) (u ) + KGU(x,b;c ) (u ),where

KGL (x,b;c ) (u ) =ux/b exp (−u/b)
ux/b exp (−u/b)1(0 ≤ u < c ) and KGU(x,b;c ) (u ) = x/b+1
1(u ≥ c ).
x/b+1b
(x/b + 1 )b
(x/b + 1 )However, neither KGL (x,b;c ) (u ) nor KGU(x,b;c ) (u ) is a legitimate kernel function in the sense that ∞0

KGL (x,b;c ) (u )du =γ (x/b + 1, c/b)and
(x/b + 1 ) ∞0

KGU(x,b;c ) (u )du =(x/b + 1, c/b).
(x/b + 1 )Therefore, we make scale-adjustments to obtain the re-normalized truncated kernels as

KG−(x,b;c ) (u ) =(x/b + 1 )(x/b + 1 )KLKU(u ) and KG+(x,b;c ) (u ) =(u ).
γ (x/b + 1, c/b) G(x,b;c )(x/b + 1, c/b) G(x,b;c )These kernels yield estimators of f− (c ) and f+ (c ) as1fˆ− (c ) =
n

n

i=1

KG−(x,b;c ) (Xi )

x=c

=

1 −
1 +
1 +
KG(c,b;c ) (Xi ) and fˆ+ (c ) =KG(x,b;c ) (Xi )=KG(c,b;c ) (Xi ).
n
n
n
x=c
n

n

n

i=1

i=1

i=1

To explore asymptotic properties of these estimators, we make the following assumptions. For notational conciseness,
expressions such as ‘f ± (c)’ are used throughout, whenever no confusions may occur.

B. Funke, M. Hirukawa / Econometrics and Statistics 9 (2019) 156–170

159

Assumption 1. {Xi }ni=1 are i.i.d. random variables drawn from a univariate distribution with a pdf f having support on R+ .
Assumption 2. There is a neighborhood N around c such that the second-order derivative of the pdf f is Hölder-continuousjjof order ς ∈ (0, 1] on N \{c}. Also let f−( ) (c ) := limx↑c d j f (x )/dx j and f+( ) (c ) := limx↓c d j f (x )/dx j for j = 1, 2. Then, f ± (c) > 0
and f±( ) (c ) < ∞.
2Assumption 3. The smoothing parameter b(= bn > 0 ) satisﬁes b + nb3/2

−1→ 0 as n → ∞.

An equivalent to Assumption 2 can be found in McCrary (2008) and Otsu
 et al. (2013). In particular, (local) Höldercontinuity of the second-order density derivative f(2) (·) implies that for N = N , N with 0 < N < c < N < ∞, there is aconstant L ∈ (0, ∞) such that

 

 

f (2 ) (s ) − f (2 ) (t ) ≤ L|s − t |ς , ∀s, t ∈ (N , c ) and f (2 ) s − f (2 ) t≤ L s −tς, ∀s , t ∈ c, N .

−1

The condition ‘ nb3/2→ 0’ in Assumption 3 establishes Proposition 2 and Theorem 1. When approximating the bias
terms of bias-corrected estimators in these statements, we require the condition to control the order of magnitude in theremainder terms. A similar rate requirement for b can be found, for example, in Hirukawa and Sakudo (2014, Assumption 2).



−1

−1

It is worth remarking that weaker conditions ‘ nb1/2→ 0’ and ‘(nb) → 0’ suﬃce for Proposition 1 and Theorem 2, respectively.
The proposition below refers to bias and variance approximations of fˆ± (c ). It is worth emphasizing that all convergences
results in this paper are built upon a few different approximation techniques on incomplete gamma functions; such proof
strategies are taken for the ﬁrst time in the econometric literature, to the best of our knowledge. Moreover, for the purposeof our subsequent analysis, the bias expansion is derived up to the second-order term.
Proposition 1. Under Assumptions 1–3, as n → ∞,Bias fˆ± (c ) = ±V ar fˆ± (c ) =2πc1/2 f±( ) (c )b1/2 +
1



1−



1
f ± (c )√ 1/2 + o n−1 b−1/2 .
1
/
2
nb
πc

4
3π



f±( ) (c ) +1c (2 )f (c ) b + o(b), and2 ±ˆ Proposition
 1 implies that f± (c ) are consistent for f ± (c), and that their variance convergence has a usual rate ofO n−1 b−1/2 . Nevertheless, the bias convergence is O(b1/2 ), which is slower than the usual O(b) rate. This is an outcomeof one-sided smoothing. If f were continuous at c and smoothing were made on both sides of the design point c using the
gamma kernel, the nearly symmetric shape of the kernel would cancel out the O(b1/2 ) bias.1 In reality, because data points
used for estimating f ± (c) lie only on either the left or right side of c, the O(b1/2 ) bias never vanishes.
1/2
1/2
ˆ
A further inspection also revealsthat coeﬃcients of the leading O(b ) bias terms of f± (c ) are proportional to c ,whereas those of the leading O n−1 b−1/2 variance terms are inversely proportional to c1/2 . As the design point c moves awayfrom the origin, the biases of fˆ± (c ) may become large (in absolute value). However, the large biases are compensated bytheir shrinking variances. Similar results can be found in the bias and variance approximations to the density estimator using
the second gamma kernel by Chen (20 0 0) and those using the inverse Gaussian and reciprocal inverse Gaussian kernels by
Scaillet (2004).
It follows that when J(c) is estimated by Jˆ(c ) := fˆ+ (c) − fˆ− (c ), it also has an inferior O(b1/2 ) bias. Therefore, our goal isto propose an estimator of J(c) with a O(b) bias and a O n−1 b−1/2 variance.
2.3. Bias-corrected estimation and inference
To improve the bias convergence in estimators of f ± (c) from O(b1/2 ) to O(b) while the order of magnitude in variance remains unchanged, we propose to employ a multiplicative bias correction (‘MBC’) technique. As in Hirukawa (2010), Hirukawa
and Sakudo (2014), Hirukawa and Sakudo (2015), and Funke and Kawka (2015), the MBC method proposed by Terrell andScott (1980) is adopted.2 The method eliminates the leading bias term by constructing a multiplicative combination of two
This can be also seen by combining two estimators fˆ± (c ) as a weighted sum.
Aforementioned articles also apply another MBC method proposed by Jones et al. (1995). However, it appears that the method fails to eliminate theO(b1/2 ) bias. Their MBC estimator of f − (c ), for example, can be written as12f˘− (c ) := fˆ− (c )ᾰ− (c ) := fˆ− (c )−n1  KG(c,b;c ) (Xi ),nfˆ− (Xi )i=1

where ᾰ− (c ) serves as the ‘bias correction’ term. However, fˆ− (x ) (x < c ) has a O(b) bias, as stated in Theorem 2, so does ᾰ− (c ). Therefore, the O(b1/2 ) bias
in fˆ− (c ) never vanishes, and thus we do not pursue this type of MBC.

160

B. Funke, M. Hirukawa / Econometrics and Statistics 9 (2019) 156–170

density estimators with different smoothing parameters. In our context, for some constant δ ∈ (0, 1), the MBC estimators off ± (c) can be deﬁned asf˜± (c ) =fˆ±,b (c )(1 / 1 −δ 1 / 2)

fˆ±,b/δ (c )(−δ 1 / 2 / 1 −δ 1 / 2),

where fˆ•,b (x ) and fˆ•,b/δ (x ) signify the density estimators using smoothing parameters b and b/δ , respectively. Not only aref˜± (c ) nonnegative by construction, but also their bias and variance convergences are usual O(b) and O n−1 b−1/2 rates,
respectively, as documented in the next proposition. The proof is similar to the one for Theorem 1 of Hirukawa and Sakudo(2014), and thus it is omitted.
Proposition 2. Under Assumptions 1–3, as n → ∞,   1 2 
  1  cf±( ) (c )
4
c
1
2
(
)
(
)
Bias f˜± (c ) =−1−f (c ) + f± (c ) b + o(b), andf ± (c )3π ±2
δ 1/2 πV ar f˜± (c ) =where
λ(δ ) :=
1f (c )
λ(δ ) √ ± 1/2 + o n−1 b−1/2 ,
nb1/2
πc


1/2

1 + δ 3/2 ( 1 + δ )√− 2 2δ
2(1 + δ )1/2 1 − δ 1/2

is monotonously increasing in δ ∈ (0, 1) withlim λ(δ ) = 1 and lim λ(δ ) =
δ↓0

δ↑1

11
.
4

Proposition 2 suggests that as δ↓0 (δ↑1) or in case of oversmoothing (undersmoothing), the bias increases (decreases)
and the variance decreases (increases). It is a common practice in nonparametric kernel testing that the bias is made asymptotically negligible via undersmoothing, and thus what matters for inference is the size of λ(δ ). Because of no minimum inλ(δ ), the issue of choosing δ will be addressed in Monte Carlo simulations.
It also follows that J(c) can be consistently estimated as J˜(c ) := f˜+ (c ) − f˜− (c ). The next theorem refers to the limiting
distribution of J˜(c ).3
Theorem 1. Under Assumptions 1–3, as n → ∞,



where





d
nb1/2 J˜(c ) − J (c ) − B(c )b → N (0, V (c ) ),

(1)



 

 1  c  f (1) (c )2  f (1) (c )2 c  (2 )
4  (1 )
+
−
(1 )
(2 )
B (c ) =−−1−
f + (c ) − f − (c ) +f (c ) − f − (c ),f + (c )f − (c )3π
2 +δ 1/2 π
f + (c ) + f − (c )V (c ) = λ (δ )
,
√ 1/2
πc

and λ(δ ) is deﬁned in Proposition 2. In addition, if nb5/2 → 0 as n → ∞, then (1) reduces to







d
nb1/2 J˜(c ) − J (c ) → N (0, V (c ) ).



As indicated in Proposition 2, J˜(c ) has a O(b) bias and a O n−1 b−1/2 variance. Observe that for a given δ , the variancecoeﬃcient decreases as c increases, i.e., as the discontinuity point moves away from the origin. We can also ﬁnd that theleading bias term B(c)b cancels out if f has a continuous second-order derivative at c.
Theorem 1 also implies that given a smoothing parameter b = Bn−q for some constants B ∈ (0,∞) and q ∈ (2/5, 2/3) (which
fulﬁlls three rate requirements for b, namely, b → 0, nb3/2 → ∞ and nb5/2 → 0 as n → ∞) and V˜ (c ), a consistent estimate ofV(c), the test statistic is

T (c ) :=√nb1/2 J˜(c )V˜ (c )d→ N (0, 1 ) under H0 : J (c ) = 0.

As with other kernel-smoothed tests, the proposed test is shown to be consistent. The next proposition documents this
aspect.
3
It is possible to use different constants δ− and δ+ and/or different smoothing parameters b− and b+ for f˜− (c ) and f˜+ (c ), as long as b− and b+ shrink
to zero at the same rate. For convenience, however, we choose to employ the same δ and b.

B. Funke, M. Hirukawa / Econometrics and Statistics 9 (2019) 156–170

161

Proposition 3. Under
Assumptions
1–3, as n → ∞, Pr {|T (c )| > Bn } → 1 under H1 : J(c) = 0 for any non-stochastic sequence Bnsatisfying Bn = o n1/2 b1/4 .
It is straightforward to see that the power approaches one for ﬁxed alternatives H1 : J (c ) = C1 for some constant|C1 | ∈ (0,∞). Consistency of the test for ﬁxed alternatives will be conﬁrmed numerically in Monte Carlo simulations.
Proposition 3 even refers to consistency of the test for local alternatives√of the form H1 : Jn (c ) = 0 + C2 rn for some constant|C2 | ∈ (0,∞) and a non-stochastic sequence rn satisfying rn + 1/ nb1/2 rn → 0 as n → ∞.
Our remaining tasks are to present examples of V˜ (c ) and to propose a choice method of b. The latter is discussed in
the next section, whereas there are a few candidates of V˜ (c ). Replacing f ± (c) in V(c) with their consistent estimates f˜± (c )immediately yields

V˜1 (c ) := λ(δ )f˜+ (c ) + f˜− (c )
.
√ 1/2
πc

Alternatively, it is possible to compute the gamma kernel density estimator at c

1fˆ(c ) :=
n

n


1
KG(c,b) (Xi ).
nn

KG(x,b) (Xi )

i=1

x=c

=

i=1

It follows from the proof of Proposition 1 that  γ (c/b + 1, c/b) 1
(c/b + 1, c/b) 1= + O b1/2 and= + O b1/2 .
(c/b + 1 )2
(c/b + 1 )2Then,γ (c/b + 1, c/b) ˆ(c/b + 1, c/b) ˆ
p f + (c ) + f − (c )fˆ(c ) =f (c ) +f (c ) →.
(c/b + 1 ) −(c/b + 1 ) +2As a consequence, we can obtain another estimator of V(c) as

2 fˆ(c )V˜2 (c ) := λ(δ ) √ 1/2 .
πc2.4. Smoothing parameter selection
How to choose the value of the smoothing parameter b is an important practical problem. McCrary (2008) proposes
the choice method which closely follows the literature on the BLL smoothing. Moreover, in the literature on RDD, Imbens
and Kalyanaraman (2012) and Porter and Yu (2015, Section 5.4)) discuss methods of choosing the smoothing parameter.
All these proposals rely on either a cross-validation criterion or a plug-in approach, and thus they stand on the idea of
estimation-optimality. However, once our priority is given to testing for continuity of the pdf f at a given cutoff c, suchapproaches cannot be justiﬁed in theory or practice, because estimation-optimal values may not be equally optimal for
testing purposes. Here, we have a preference for test-optimality and thus adopt the power-optimality criterion by Kulasekera
and Wang (1998), whose idea is also applied in Hirukawa and Sakudo (2016).
Below Procedure 1 of Kulasekera and Wang (1998) is tailored to our context. The procedure is a version of sub-sampling.
Let n− and n+ be the numbers of observations in sub-samples {Xi− } and {Xi+ }, respectively, where n ≡ n− + n+ . Also assume



n−



n+



n− 

n+ 

that Xi−
and Xi+are ordered samples. Then, the entire sample {Xi }ni=1 = Xi−, Xi+
can be split into M
i=1
i=1
i=1
i=1
sub-samples, where M = Mn is a non-stochastic sequence that satisﬁes 1/M + M/n → 0 as n → ∞. Given such M, (k− , k+ ) :=(n− /M, n+ /M) and k := k− + k+ , the mth sub-sample is deﬁned as{Xm,i }ki=1 := −Xm+(i−1 )Mk−  +
i=1

, Xm+(i−1 )M

k+

i=1



The test statistic using the mth sub-sample Xm,iTm (c ) :=√kb1/2 J˜m (c )V˜m (c ), m = 1, . . . , M.

ki=1

becomes, m = 1, . . . , M,

where J˜m (c ) and V˜m (c ) (which is either V˜1,m (c ) or V˜2,m (c )) are the sub-sample
analogues of J˜(c ) and V˜ (c ), respectively. Also
denote the set of admissible values for b = bn as Hn := Bn−q , Bn−q for some prespeciﬁed exponent q ∈ (2/5, 2/3) and two
constants 0 < B < B < ∞. Moreover, letπˆ M (bk ) :=1 
1{Tm (c ) > cm (α )},
M
M

m=1

where cm (α ) is the critical value for the size α test using the mth sub-sample. We pick the power-maximized bˆ k = Bˆk−q =arg maxbk ∈Hk πˆ M (bk ), and the smoothing parameter value bˆ n := Bˆn−q follows.

162

B. Funke, M. Hirukawa / Econometrics and Statistics 9 (2019) 156–170

We conclude this section by stating how to obtain bˆ n in practice. Step 1 reﬂects that M should be divergent but smaller
than both n− and n+ in ﬁnite samples. Step 3 follows from the implementation methods in Kulasekera and Wang (1998).
Finally, Step 4 corresponds to the case for more than
 one
 maximizer
 of πˆ M (bk ).
p
p
Step 1: Choose some p ∈ (0, 1) and specify M = min n− , n+ .
Step 2: Make M sub-samples of sizes (k− , k+ ) = (n− /M, n+ /M ) and put k = k− + k+ .
Step 3: Pick two constants 0 < H < H < 1 and deﬁne Hk = H , H .
Step 4: Set cm (α ) ≡ z(i.e., Pr N (0, 1 ) > zα /2 = α /2) and ﬁnd bˆ = inf arg max
πˆ M (b ) by a grid search.
α /2

bk ∈Hkk

Step 5: Recover Bˆ by Bˆ = bˆ k kq for some q ∈ (2/5, 2/3) and obtain bˆ n = Bˆn−q .

k

3. Estimation of the entire density in the presence of a discontinuity point3.1. Density estimation by truncated kernels
We are typically interested in how the shape of the pdf looks like, as well as whether it has a discontinuity point. Imbens
and Lemieux (2008) strongly recommend graphical analyses in empirical studies on RDD, including inspections of densities
of running variables. If the test in the previous section fails to reject the null of continuity of the pdf f at the cutoff c, theentire density may be re-estimated by the gamma kernel, for example. How should we estimate the entire density if thetest rejects the null?
The answer to this question is simple. It suﬃces to compute fˆ− (x ) or fˆ+ (x ) as an estimate of f(x), depending on theposition of the design point x. To put it in another way, fˆ− (x ) ( fˆ+ (x )) can be employed whenever x < c (x > c), provided
that c is the only point of discontinuity in f, as documented in Theorem 2 below. Because the theorem refers to the density
estimation at design points other than c, we impose a global smoothness condition on f as a new assumption and then state
the theorem.
Assumption 4. The second-order derivative of the pdf f is Hölder-continuous of order ς ∈ (0, 1] on R+ \{c}.
Theorem 2. Suppose that Assumptions 1, 3 and 4 hold. Then, for x > c , as n → ∞,Bias fˆ+ (x ) =f (1 ) (x ) +
x (2 )1f (x )f (x ) b + o(b), and V ar fˆ+ (x ) =√ 1/2 + o n−1 b−1/2 .
1/
2
2
nb 2 π x

On the other hand, for x < c , as n → ∞,Bias fˆ− (x ) =x( ) + f (2) b + o(b), and V ar fˆ− (x ) =2

f (1 ) x

 1



f x
√( )
+ o n−1 b−1/2
nb1/2 2 π x1/2

(2κ +1 )1
f x + o n−1 b−1
nb 22κ +1  2 (κ +1 ) ( )i f x/b → ∞i f x/b → κ ∈ (0, ∞ ).

Theorem 2 indicates no adversity when f(x) for x = c is estimated by fˆ± (x ); although only the bias-variance trade-off isprovided, asymptotic normality of the estimators can be established similarly to this theorem. Observe that fˆ± (x ) admit thesame bias and variance expansions as the gamma kernel density estimator fˆ(x ) does. A rationale is that as the design point
x moves away from the truncation point c, data points tend to lie on both sides of x and each truncated kernel is likely to
behave like the gamma kernel. We can also see that the variance coeﬃcient decreases as x increases. The shrinking variancecoeﬃcient as the design point x moves away from the origin reﬂects that more data points can be pooled to smooth inareas with fewer observations. This property is particularly advantageous to estimating the distributions that have a longtail with sparse data, such as those of the economic variables mentioned in Section 1.
3.2. Convergence properties of fˆ− (x ) when the density is unbounded at the origin
Clusterings of observations near the boundary are frequently observed in the distributions with positive supports. In the
study of RDD, Figure 1 of Bertrand et al. (2015) suggests that the distribution of wives’ relative income within households
has a clustering of observations near the origin, as well as a sharp drop at the point of 1/2 (i.e., the point at which wives’
income shares exceed their husbands’). Similarly, in Figure 5 of McCrary (2008), the distribution of proportion of votes for
proposed bills in the US House of Representatives appears to be unbounded at the boundary of 100%, as well as a sharp
discontinuity at the point of 50%.4 The following two theorems document weak consistency and the relative convergence offˆ− (x ) when f(x) is unbounded at x = 0.



Theorem 3. If f(x) is unbounded at x = 0, Assumptions 1 holds and b + nb2

4

−1

p
→ 0 as n → ∞, then fˆ− (0 ) → ∞.

The arguments in this section are still valid for this case, if we transform the original data X to X := 1 − X and apply them to the transformed data X .

B. Funke, M. Hirukawa / Econometrics and Statistics 9 (2019) 156–170

163

Theorem 4. Suppose that f(x) is unbounded at x = 0 and continuously differentiable in the neighborhood of the origin. In addi-−1tion, if Assumption 1 holds and b + nb2 f (x )→ 0 as n → ∞ and x → 0, thenfˆ− (x ) − f (x ) p→0f (x )
as x → 0.
It has been demonstrated by Bouezmarni and Scaillet (2005) and Hirukawa and Sakudo (2015) that the weak consistency
and relative convergence for densities unbounded at the origin are peculiar to the density estimators smoothed by the
gamma and generalized gamma kernels. The theorems ensure that fˆ− (x ) is also a proper estimate for unbounded densities.
We can deduce from Theorems 2–4 that all in all, appealing properties of the gamma kernel density estimator are inherited
to fˆ± (x ).
4. Finite-sample performance
It is widely recognized that asymptotic results on kernel-smoothed tests are not well transmitted to their ﬁnite-sample
distributions, which reﬂects that omitted terms in the ﬁrst-order asymptotics on the test statistics are highly sensitive
to their smoothing parameter values in ﬁnite samples. On the other hand, there is growing literature that reports niceﬁnite-sample properties of the estimators and test statistics smoothed by asymmetric kernels. Examples include Kristensen(2010) and Gospodinov and Hirukawa (2012) for estimation and Fernandes and Grammig (2005), Fernandes et al. (2015), andHirukawa and Sakudo (2016) for testing. To see which perspective dominates, this section investigates ﬁnite-sample performance of the estimator of the jump-size magnitude and the test statistic for discontinuity of the density via Monte Carlosimulations.
4.1. Jump-size estimation
First, we focus on the estimator of the jump-size magnitude J(c). As true densities, those of the following two asymmetric
distributions are considered:

1. Gamma: f (x ) = xα −1 exp (−x/β )1(x ≥ 0 )/{β α (α )}, (α , β ) = (2.75, 1 ).

α −1
α
2. Weibull: f (x ) = (α /β )(x/β )exp −(x/β ) 1(x ≥ 0 ), (α , β ) = (1.75, 3.5 ).
Shapes of these densities can be found in Fig. 1. For each distribution we choose two suspected discontinuity points c,
namely, 30% quantile (‘30%’) and median (‘ Med’); see Table 1 for exact values of the points. Because the gamma and Weibull
densities have modes at 1.7500 and 2.1567, respectively, the two points for each density are located on the left- and right-

Fig. 1. Shapes of true densities for Monte Carlo simulations.

164

B. Funke, M. Hirukawa / Econometrics and Statistics 9 (2019) 156–170
Table 1
Biases, standard deviations and RMSEs of estimators of J(c).
DistributioncnEstimatorJ˜(c ) with δJˆM (c )Gamma

1.7057(30%)

500

10 0 0

20 0 02.4248(Med)500

10 0 0

20 0 0Weibull

1.9419(30%)

500

10 0 0

20 0 02.8386(Med)500

10 0 0

20 0 0

0.49

0.64

0.81

Bias
StdDev
RMSE
Bias
StdDev
RMSE
Bias
StdDev
RMSE
Bias
StdDev
RMSE
Bias
StdDev
RMSE
Bias
StdDev
RMSE−0.03810.04610.0598−0.0335
0.03310.0471−0.0283
0.02500.0377−0.04070.04800.0629−0.03230.03530.0479−0.02400.02710.03620.00190.07860.07860.00190.05880.05880.0 0 06
0.04300.0430−0.00110.0626
0.0626−0.00120.04630.0463−0.0 0 04
0.03510.03510.00110.08120.08120.00100.06070.0607
0.0 0 020.04450.0445−0.00110.06480.0648−0.00120.04790.0479−0.0 0 04
0.03630.03630.0 0 06
0.08370.0837
0.0 0 06
0.0626
0.0626−0.0 0 0 00.04580.0458−0.0011
0.06680.0668−0.00130.04940.0494−0.0 0 04
0.03740.0374Bias
StdDev
RMSE
Bias
StdDev
RMSE
Bias
StdDev
RMSE
Bias
StdDev
RMSE
Bias
StdDev
RMSE
Bias
StdDev
RMSE−0.0235
0.04160.0478−0.0187
0.03020.0355−0.01440.02250.0267−0.02460.04050.0474−0.0195
0.02900.0350−0.0149
0.02180.02640.00240.06650.06650.00350.05000.05020.0017
0.0367
0.0367
0.0 0 06
0.05340.05340.0 0 020.03940.03940.0 0 070.02990.02990.00120.06840.06840.00130.05090.05090.0 0 03
0.03720.03720.0 0 04
0.05520.0552−0.0 0 01
0.04080.04080.0 0 05
0.03090.03090.0 0 05
0.07040.07040.0 0 05
0.05230.05230.0 0 01
0.03830.03830.0 0 02
0.05690.0569−0.0 0 03
0.0421
0.0421
0.0 0 04
0.03190.0319hand sides of the mode. The sample size is n ∈ {50 0, 10 0 0, 20 0 0}, and 50 0 0 replications are drawn for each combination of
the sample size n and the distribution.
The simulation study compares ﬁnite-sample performance of our jump-size estimator J˜(c ) with McCrary’s (2008) BLL
estimator JˆM (c ). The latter employs the triangular kernel K (u ) = (1 − |u| )1(|u| ≤ 1 ), and the bandwidth is chosen by themethod described on p.705 of McCrary (2008). For the former,
 the smoothing parameter b is selected by the power√
optimality criterion for two test statistics Ti (c ) := nb1/2 J˜(c )/ V˜i (c ) for i = 1, 2, where the deﬁnition of V˜i (c ) is given in

Section 2.3. Implementation details are as follows: (i) all critical values in πˆ M (bk ) are set equal to z0.025 = 1.96; (ii) (p, q)
are predetermined by ( p, q ) = (1/2, 4/9 ); (iii) the interval for bk is Hk = [0.05, 0.50]; and (iv) three different values of themixing exponent δ are considered, namely, δ ∈ {0.49, 0.64, 0.81}, so that the exponents on fˆ±,b (c ) and fˆ±,b/δ (c ) to generatef˜± (c ) are (10/3, −7/3 ), (5, −4 ) and (10, −9 ), respectively.
Table 1 presents as performance measures the bias (‘ Bias’), standard deviation (‘ StdDev’) and root-mean squared error(‘ RMSE’) of each estimator over 50 0 0 Monte Carlo samples. Since the densities are continuous at c actually, the performancemeasures are calculated on the basis of J (c ) = 0. Moreover, only the performance measures with the smoothing parameterb selected for T2 (c) are reported, because there is no substantial difference between values of b chosen for T1 (c) and T2 (c).
It can be immediately found that the RMSE shrinks with the sample size, which indicates consistency of each estimator.
Although the Bias of JˆM (c ) is larger than that of J˜(c ), the StdDev of the former is smaller, and as a consequence it tends
to yield a smaller RMSE. We can also see that J˜(c ) has extremely small biases for all cases, which conﬁrms that the MBC
technique leads to huge bias reduction. The bias-variance trade-off in terms of δ within J˜(c ) (as Proposition 2 suggests) canbe also observed.

B. Funke, M. Hirukawa / Econometrics and Statistics 9 (2019) 156–170

165

Table 2Finite-sample size properties of test statistics for discontinuity.
(%)
DistributioncnNominalT1 (c) with δ
0.49

0.64

0.81

0.49

0.64

0.81

Gamma

1.7057
(30%)

500

5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%

3.2
7.5
3.9
8.4
3.5
8.1
3.3
7.9
3.7
8.0
4.7
8.8

3.1
7.5
3.9
8.2
3.6
8.2
3.6
7.8
3.8
8.2
4.7
8.9

3.1
7.8
3.9
8.2
3.7
8.4
3.6
7.7
3.9
8.0
4.8
9.0

4.4
8.8
6.1
10.7
4.2
8.8
3.8
8.7
4.1
8.6
4.9
9.4

4.0
8.8
4.6
9.2
3.9
8.5
3.9
8.6
4.2
8.6
5.0
9.4

3.7
8.6
4.4
8.9
3.9
8.7
4.0
8.5
4.3
8.6
5.1
9.5

5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%

3.2
7.7
4.0
8.2
3.8
8.3
3.6
7.8
3.7
8.1
4.7
8.9

3.2
7.8
4.2
8.3
3.7
8.4
3.6
7.7
3.8
8.1
4.7
9.0

3.3
7.9
4.2
8.4
3.8
8.3
3.5
7.7
3.8
8.2
4.8
9.2

6.2
10.7
10.2
14.7
7.7
12.4
3.9
8.5
4.0
8.7
4.9
9.4

4.9
9.4
6.4
10.7
4.4
9.0
4.0
8.5
4.2
8.4
5.0
9.4

4.1
9.0
5.2
9.4
4.0
8.5
3.9
8.3
4.2
8.6
5.0
9.6

10 0 0
20 0 02.4248(Med)500
10 0 0
20 0 0Weibull

1.9419(30%)

500
10 0 0
20 0 02.8386(Med)500
10 0 0
20 0 0T2 (c) with δ

4.2. Testing for discontinuity
Second, size and power properties of the test statistic T(c) are investigated. In what follows, T1 (c) and T2 (c) are comparedwith McCrary’s (2008) test statistic based on the difference between logarithms of two density estimates, denoted as TM (c).
Implementation details of each test statistic are the same as described above. The Monte Carlo design in this section is
inspired by Otsu et al. (2013). Let X be drawn with probability γ from the truncated gamma or Weibull distribution withsupport on [0, c) and with probability 1 − γ from the one with support on (c, ∞). Unless γ = Pr (X ≤ c ), the gamma or
Weibull pdf is discontinuous at c. Also denote the measure of discontinuity as d := Pr (X ≤ c ) − γ , where d ∈ {0.00, 0.02, 0.04,
0.06, 0.08, 0.10} and d > 0 (⇔J(c) > 0) suggests a jump of the pdf at c. For each statistic, the empirical rejection frequenciesof the null H0 : J (c ) = 0 for d = 0 and d > 0 indicate its size and power properties, respectively.
Table 2 presents size properties of T1 (c) and T2 (c). Each test statistic exhibits mild under-rejection of the null except a
few cases, and the rejection frequencies of T2 (c) are closer to the nominal ones. The rejection frequencies tend to decrease
with δ , and substantial over-rejection of the null is not observed for δ = 0.81. Considering that δ = 0.81 also yields nearly
unbiased estimates of J(c), we set δ equal to this value for power comparisons.
Table 3 reports power properties of T1 (c) and T2 (c), in comparison with TM (c). Panel (A) refers to the results from the
gamma distribution. It can be observed that the rejection frequency of each test statistic for a given d > 0 approaches to
one with the sample size n, which indicates consistency of each test. Both T1 (c) and T2 (c) exhibit good power properties
without inﬂating their sizes, and T2 (c) appears to be more powerful than T1 (c). In contrast, TM (c) exhibits considerable size
distortions, and nonetheless its power properties look inferior to those of T1 (c) and T2 (c). It may be argued that the gamma
distribution is too advantageous to T1 (c) and T2 (c) in that both rely on the gamma kernel. Hence, the simulation study basedon the Weibull distribution could be fair, and the results are reported in Panel (B). Indeed, the size properties of TM (c) are
dramatically improved. However, it is still outperformed in terms of power properties by T1 (c) and T2 (c). Again in this case,
it appears that T2 (c) has better power properties than T1 (c). A possible rationale is that because V˜2 (c ) tends to be smaller
than V˜1 (c ), as suggested in Proposition 2, T2 (c) is likely to have a large value (i.e., tends to reject the null more often) thanT1 (c) under the alternative.
In sum, Monte Carlo results conﬁrm the following two respects. First, the MBC technique achieves huge bias reduction,and the jump-size estimator J˜(c ) yields nearly unbiased estimates. Second, the test statistics T1 (c) and T2 (c) exhibit nicepower properties without sacriﬁcing their size properties, whereas the latter appears to be more powerful than the former.
It is also worth emphasizing that the superior performance is based simply on ﬁrst-order asymptotic results. Therefore,
assistance of size-adjusting devices such as bootstrapping appears to be unnecessary, unlike most of the smoothed tests
employing conventional symmetric kernels.

166

B. Funke, M. Hirukawa / Econometrics and Statistics 9 (2019) 156–170
Table 3Finite-sample power properties of test statistics for discontinuity.
(A) Gamma distributionc1.7057
(30%)

n

500

Nominal

d
0.00

0.02

0.04

0.06

0.08

0.10

TM (c)

5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%

4.6
10.1
3.1
7.8
3.7
8.6
6.9
13.6
3.9
8.2
4.4
8.9
9.9
19.2
3.7
8.4
3.9
8.7

1.4
4.1
4.5
9.4
14.2
19.0
1.5
4.1
6.6
12.4
13.9
19.2
1.6
4.5
8.4
15.2
25.1
30.4

2.2
6.1
10.7
17.0
44.7
53.6
4.4
9.5
12.8
21.4
50.3
54.9
11.7
21.1
36.8
44.1
90.2
94.7

9.2
17.9
17.9
26.8
60.7
63.9
22.0
35.3
37.1
46.0
90.8
92.7
52.1
66.6
98.8
99.2
99.5
99.9

28.9
43.5
85.9
88.5
93.8
97.2
58.7
72.9
98.7
99.0
99.5
99.9
90.1
95.6
100.0
100.0
99.9
99.9

55.1
68.7
98.1
99.2
98.8
99.7
87.5
93.1
100.0
100.0
100.0
100.0
99.5
99.8
100.0
100.0
100.0
100.0

5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%

9.3
17.4
3.6
7.7
4.0
8.5
11.5
20.1
3.9
8.0
4.3
8.6
12.0
20.7
4.8
9.0
5.1
9.5

3.5
8.8
4.4
9.1
4.8
9.8
4.0
8.9
5.0
10.4
5.6
11.1
3.6
8.0
7.7
13.9
8.2
14.6

2.6
6.0
7.1
13.4
7.7
14.5
3.5
7.5
10.4
18.3
11.3
19.2
7.1
13.8
18.1
28.3
18.9
29.2

4.4
8.7
12.3
21.1
13.6
22.4
9.2
16.0
20.7
32.1
22.2
33.1
23.9
36.0
37.8
50.2
38.8
51.4

9.1
16.5
20.6
31.0
21.9
32.3
23.1
36.1
35.5
49.0
36.8
50.3
55.6
68.3
60.7
72.6
61.7
73.5

18.7
29.7
30.1
42.5
32.7
44.2
46.4
61.7
53.2
65.8
55.0
67.0
83.9
91.6
80.2
87.9
85.5
90.7

T1 (c)T2 (c)
10 0 0TM (c)T1 (c)T2 (c)

20 0 0TM (c)T1 (c)T2 (c)

2.4248(Med)

500

TM (c)T1 (c)T2 (c)

10 0 0TM (c)T1 (c)T2 (c)

20 0 0(%)

Test

TM (c)T1 (c)T2 (c )

(continued on next page)

5. Empirical illustration
This section applies our estimation and testing procedures of discontinuity in densities to real data. We employ the data
sets on fourth and ﬁfth graders of Israeli elementary schools used by Angrist and Lavy (1999). The data sets are madepublic on the Angrist Data Archive web page, and they are often utilized in empirical application parts of the closely related
literature (e.g., Otsu et al., 2013; Feir et al., 2016).
Following Maimonides’ rule, Israeli public schools make each class size no greater than 40. As a result of strategic behavior on schools’ and/or parents’ sides, the density of school enrollment counts for each grade may be discontinuous atmultiples of 40. Then, setting the cutoff c = 40, 80, 120, 160 for enrollment densities of fourth and ﬁfth graders, we estimate
the jump size and conduct the test for the null of continuity at each cutoff. Speciﬁcally, the results from our truncatedgamma-kernel approach are compared with those from McCrary’s (2008) BLL method. T2 (c) with δ = 0.81 is chosen as ourtest statistic because of its better ﬁnite-sample properties. The smoothing parameter for our approach and the bandwidth
for the BLL method are chosen in the same manners as in Section 4.
Table 4 presents estimation and testing results on discontinuity in enrollment densities, where fˆ−M (c ) and fˆ+M (c ) are
BLL estimates of left and right limits of the density at the cutoff c, respectively. The table shows remarkable differences
between estimation results from McCrary’s (2008) and our procedures. The former ﬁnds upward jump estimates only at
c = 40 for each grade. On the other hand, the latter yields upward jump estimates at c = 40, 80 and downward jump (ordrop) estimates at c = 120, 160 for each grade. Testing results also differ. While McCrary’s (2008) test rejects the null ofcontinuity at the cutoff only for three cases (i.e., c = 40, 120 for fourth graders and c = 40 for ﬁfth graders), rejections of
the null by our test include additional two cases (i.e., c = 160 for fourth graders and c = 120 for ﬁfth graders) as well as thethree cases. This appears to reﬂect better ﬁnite-sample power properties of T2 (c) reported in Section 4.

B. Funke, M. Hirukawa / Econometrics and Statistics 9 (2019) 156–170

167

Table 3 (continued)(B) Weibull distributioncn1.9419(30%)

500

Nominal

d
0.00

0.02

0.04

0.06

0.08

0.10

TM (c)

5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%
5%
10%

3.4
7.6
3.3
7.9
4.1
9.0
4.2
8.8
4.2
8.4
5.2
9.4
4.5
9.5
3.8
8.3
4.0
8.5
4.6
9.4
3.5
7.7
3.9
8.3
5.7
11.2
3.8
8.2
4.2
8.6
6.6
12.7
4.8
9.2
5.0
9.6

1.9
5.0
4.9
9.5
16.7
21.8
2.2
5.2
6.5
12.4
17.5
23.2
3.1
7.1
8.3
14.8
33.2
39.1
2.5
6.1
4.5
9.2
4.8
9.7
2.5
6.2
5.1
10.7
5.7
11.3
3.1
6.7
7.8
14.1
8.1
14.6

4.0
9.0
14.4
19.8
42.4
53.0
7.5
15.4
12.9
21.2
51.0
56.5
18.2
30.1
53.7
58.0
87.4
93.0
3.0
6.6
7.2
13.7
7.6
14.4
4.8
9.1
10.2
18.1
11.3
19.0
10.0
17.7
18.0
28.1
18.6
28.7

11.8
21.2
19.1
26.9
57.5
61.0
26.3
39.3
42.1
49.1
88.5
91.1
53.7
66.6
98.8
99.5
98.8
99.6
6.0
11.8
12.3
21.2
14.6
22.9
12.7
22.5
20.2
31.5
21.3
32.5
31.8
46.0
36.2
49.0
37.8
49.9

28.5
41.2
85.0
88.8
90.6
95.5
55.0
68.2
98.4
99.1
98.9
99.7
84.9
91.5
100.0
100.0
99.9
100.0
12.6
22.2
19.9
30.3
22.1
32.2
31.0
45.0
34.3
47.3
36.7
48.9
63.9
76.1
58.6
70.8
64.7
74.6

48.8
61.8
97.1
98.6
97.8
99.2
79.3
87.4
99.9
100.0
99.9
100.0
97.5
99.0
100.0
100.0
100.0
100.0
24.4
36.8
29.1
41.1
43.6
51.9
55.2
68.7
50.9
63.8
61.7
70.6
87.6
93.3
79.8
87.1
99.2
99.5

T1 (c)T2 (c)
10 0 0TM (c)T1 (c)T2 (c)

20 0 0TM (c)T1 (c)T2 (c)

2.8386(Med)

500

TM (c)T1 (c)T2 (c)

10 0 0TM (c)T1 (c)T2 (c)

20 0 0(%)

Test

TM (c)T1 (c)T2 (c)

Note. The value of δ for each of T1 (c) and T2 (c) is set equal to 0.81.

Table 4
Estimation and testing for the discontinuity of densities of school enrollmentsncBinned local linear methodfˆ−M (c )fˆ+M (c )Truncated kernel methodJˆM (c )TM (c)f˜− (c )f˜+ (c )J˜(c )T2 (c)

0.00980.00900.00440.0 0 05

0.00640.0 0 03
−0.0020−0.0 0 08

5.76
0.24
−3.55
−2.88

0.01160.01030.00430.00100.00740.0017−0.0014−0.0 0 04

6.28
1.25
−2.84
−1.28

2059

40
80
120
160

0.00460.01030.00610.00110.00960.00970.00390.0 0 090.0050−0.0 0 06−0.0022-0.0 0 03(a) Fourth graders:
5.61
0.0034−0.62
0.0086−3.35
0.0063−0.84
0.0013

2029

40
80
120
160

0.00550.01070.00540.00140.01140.00980.00450.00110.0059−0.0 0 09−0.0 0 09−0.0 0 03(b) Fifth graders:
6.29
0.0042−0.98
0.0087−1.20
0.0057−0.80
0.0014

Note: The value of δ for T2 (c) is set equal to 0.81. Values of test statistics in bold faces indicate signiﬁcance at the 5% level.

For convenience, density estimates with (possible) discontinuity points at c = 40, 120 are plotted in Fig. 2, where smoothing parameter values chosen for the corresponding discontinuity tests are applied for the entire density estimates. The ﬁgure
illustrates that the truncated gamma density estimators tend to capture peaks and troughs more clearly.
Furthermore, it follows from Table 4 that enrollment densities for fourth and ﬁfth graders are suspected to be discontinuous at c = 40, 120, 160 and c = 40, 120, respectively. Then, assuming that the densities are in fact discontinuous at all thesecutoffs, we draw Fig. 3 by extending the truncated kernel method to the cases of multiple cutoffs. Smoothing parameter

168

B. Funke, M. Hirukawa / Econometrics and Statistics 9 (2019) 156–170

Fig. 2. Density estimates of school enrollments.
Note: In each panel, solid and dashed lines are density estimates via the truncated gamma kernels and the binned local linear method, respectively. The‘ × ’ symbols indicate binned data points.

Fig. 3. Density estimates of school enrollments via the truncated kernel method in the presence of multiple cutoffs.
Note: Solid and dashed lines are density estimates for fourth and ﬁfth graders, respectively.

B. Funke, M. Hirukawa / Econometrics and Statistics 9 (2019) 156–170

169

values are averages of those chosen for the corresponding discontinuity tests at the cutoffs. The ﬁgure should be interpretedwith caution; it is prepared only for illustrative purposes and thus with no theoretical justiﬁcations.
6. Conclusion
This paper has developed estimation and testing procedures on discontinuity in densities with positive support. Our proposal is built on smoothing by the gamma kernel. To preserve its appealing properties, we split the gamma kernel into twoparts at a prespeciﬁed cutoff and construct two truncated kernels after re-normalization. The jump-size magnitude of thedensity at the cutoff can be estimated nonparametrically by two truncated kernels and the MBC technique by Terrell andScott (1980). The estimator is easy to implement, and its convergence properties are explored by means of various approximation techniques on incomplete gamma functions. Given the jump-size estimator, two versions of test statistics for thenull of continuity at the cutoff are also proposed, and a smoothing parameter selection method under the power-optimality
criterion is tailored to our testing procedure. Furthermore, estimation theory of the entire density in the presence of a discontinuity point is provided. It is demonstrated that density estimators smoothed by the truncated gamma kernels admit
the same bias and variance approximations as the gamma kernel density estimator does. Monte Carlo simulations indicate that the jump-size estimator is nearly unbiased when there is no jump in the true density, and that the test statistics
with power-optimal smoothing parameter values plugged in enjoy more power than McCrary’s (2008) BLL-based test does,
without sacriﬁcing their size properties.
We conclude this paper by noting a few research extensions. First, the assumption of a single (given) point of discontinuity may be relaxed. It is worth investigating the cases for more than one (prespeciﬁed) point of discontinuity or those foreven unknown (ﬁnite) number of discontinuity points. For the latter, locations of discontinuity points are estimated ﬁrst andthen the corresponding upper and lower limits of the density can be evaluated at each estimated location. Second, while
our focus has been exclusively on univariate densities, the discontinuity analysis may be extended to multivariate densities.
Acknowledgments
We would like to thank the Associate Editor and two anonymous referees for their constructive comments and suggestions. The ﬁrst author gratefully acknowledges ﬁnancial support from the German Science Foundation (DFG) through the
Collaborative Research Center SFB 823. The second author gratefully acknowledges ﬁnancial support from Japan Society of
the Promotion of Science (grant number 15K03405).
Supplementary material
Supplementary material associated with this article can be found, in the online version, at 10.1016/j.ecosta.2017.07.006 .
References
Angrist, J.D., Lavy, V., 1999. Using Maimonides’ rule to estimate the effect of class size on scholastic achievement. Q. J. Econ. 114, 533–575.
Bertrand, M., Kamenica, E., Pan, J., 2015. Gender identity and relative income within households. Q. J. Econ. 130, 571–614.
Bouezmarni, T., Scaillet, O., 2005. Consistency of asymmetric kernel density estimators and smoothed histograms with application to income data. Econom.
Theory 21, 390–412.
Chen, S.X., 1999. Beta kernel estimators for density functions. Comput. Stat. Data Anal. 31, 131–145.
Chen, S.X., 20 0 0. Probability density function estimation using gamma kernels. Ann. Inst. Stat. Math. 52, 471–480.
Chu, C.K., Cheng, P.E., 1996. Estimation of jump points and jump values of a density function. Stat. Sin. 6, 79–95.
Cline, D.B., Hart, J.D., 1991. Kernel estimation of densities with discontinuities or discontinuous derivatives. Stat. J. Theor. Appl. Stat. 22, 69–84.
Cowell, F.A., Ferreira, F.H.G., Litchﬁeld, J.A., 1998. Income distribution in Brazil 1981–1990: parametric and nonparametric approaches. J. Income Distrib. 8,
63–76.
DiNardo, J., Fortin, N.M., Lemieux, T., 1996. Labor market institutions and the distribution of wages, 1973–1992: a semiparametric approach. Econometrica
64, 1001–1044.
Feir, D., Lemieux, T., Marmer, V., 2016. Weak identiﬁcation in fuzzy regression discontinuity designs. J. Bus. Econ. Stat. 34, 185–196.
Fernandes, M., Grammig, J., 2005. Nonparametric speciﬁcation tests for conditional duration models. J. Econom. 127, 35–68.
Fernandes, M., Mendes, E.F., Scaillet, O., 2015. Testing for symmetry and conditional symmetry using asymmetric kernels. Ann. Inst. Stat. Math. 67, 649–671.
Funke, B., Kawka, R., 2015. Nonparametric density estimation for multivariate bounded data using two non-negative multiplicative bias correction methods.
Comput. Stat. Data Anal. 92, 148–162.
Gospodinov, N., Hirukawa, M., 2012. Nonparametric estimation of scalar diffusion models of interest rates using asymmetric kernels. J. Empir. Financ. 19,
595–609.
Hirukawa, M., 2010. Nonparametric multiplicative bias correction for kernel-type density estimation on the unit interval. Comput. Stat. Data Anal. 54,
473–495.
Hirukawa, M., Sakudo, M., 2014. Nonnegative bias reduction methods for density estimation using asymmetric kernels. Comput. Stat. Data Anal. 75, 112–123.
Hirukawa, M., Sakudo, M., 2015. Family of the generalised gamma kernels: a generator of asymmetric kernels for nonnegative data. J. Nonparametr. Stat.
27, 41–63.
Hirukawa, M., Sakudo, M., 2016. Testing symmetry of unknown densities via smoothing with the generalized gamma kernels. Econometrics 4 (Article 28).
Imbens, G.W., Kalyanaraman, K., 2012. Optimal bandwidth choice for the regression discontinuity estimator. Rev. Econ. Stud. 79, 933–959.
Imbens, G.W., Lemieux, T., 2008. Regression discontinuity designs: a guide to practice. J. Econom. 142, 615–635.
Jones, M.C., Linton, O., Nielsen, J.P., 1995. A simple bias reduction method for density estimation. Biometrika 82, 327–338.
Kristensen, D., 2010. Nonparametric ﬁltering of the realized spot volatility: a kernel-based approach. Econom. Theory 26, 60–93.
Kulasekera, K.B., Wang, J., 1998. Bandwidth selection for power optimality in a test of equality of regression curves. Stat. Probab. Lett. 37, 287–293.
Liebscher, E., 1990. Kernel estimators for probability densities with discontinuities. Stat. J. Theor. Appl. Stat. 21, 185–196.
McCrary, J., 2008. Manipulation of the running variable in the regression discontinuity design: a density test. J. Econom. 142, 698–714.

170

B. Funke, M. Hirukawa / Econometrics and Statistics 9 (2019) 156–170

Otsu, T., Xu, K.-L., Matsushita, Y., 2013. Estimation and inference of discontinuity in density. J. Bus. Econ. Stat. 31, 507–524.
Porter, J., Yu, P., 2015. Regression discontinuity designs with unknown discontinuity points: testing and estimation. J. Econom. 189, 132–147.
Saez, E., 2010. Do taxpayers bunch at kink points? Am. Econ. J. Econ. Policy 2, 180–212.
Scaillet, O., 2004. Density estimation using inverse and reciprocal inverse gaussian kernels. J. Nonparametr. Stat. 16, 217–226.
Terrell, G.R., Scott, D.W., 1980. On improving convergence rates for nonnegative kernel density estimators. Ann. Stat. 8, 1160–1163.