Univariate Kernel Density Estimation

Univariate kernel density estimationBen JannETH Zürich, Switzerland, jann@soz.gess.ethz.chAugust 3, 2007Abstract. The methods and formulas used by the kdens package for Stata 9.2 arediscussed in this document. To install kdens, type ssc install kdens in Stata.Note that the moremata package is required (type ssc install moremata).Preliminary note: In addition to the methods presented below, kdens also supportsconfidence interval and variance estimation using bootstrap and jackknife techniques.The implementation of the bootstrap and jackknife for density estimation is straightforward and is therefore not discussed here.1Standard kernel density estimationLet X1 , . . . , Xn be a sample from X, where X has the probability density functionf (x). Furthermore, let w1 , . . . , wn be associated weights (set wi 1 i if there are noweights). The density of X can then be estimated as: n1 X wix XibfK (x; h) K(1)W i 1 hhPnwhere W i 1 wi and K(z) is a kernel function (see Section 9). h is the smoothingparameter (the kernel halfwidth or “bandwidth”). Formula (1) is also used, for example,by official Stata’s kdensity (see [R] kdensity).2Adaptive kernel density estimationThe adaptive kernel density estimator is defined as n1 X wix XiafbK(x; h) KW i 1 hλihλi(2)where λ, the local bandwidth factors, are based on a preliminary fixed bandwidth densityestimate. The factors are estimated asv uu G fbK (X; h)bi t, i 1, . . . , n(3)λfbK (Xi ; h)where G() stands for the geometric mean over all i. Note that G(λ) 1 and thusG(hλ) h. The estimator is based on Abramson (1982). Also see, for example,

2Univariate kernel density estimationSilverman (1986, 100–110), Fox (1990, 100–103), Salgado-Ugarte et al. (1993), SalgadoUgarte and Pérez-Hernández (2003), or Van Kerm (2003).Technical note: fbK (Xi ; h) is determined by linear interpolation if Xi falls betweenthe points at which the preliminary density estimate has been evaluated.3Approximate variance estimationAn approximate estimator for the variance of fbK (x; h) at point x is given as11Ve {fbK (x; h)} R(K)fbK (x; h) fbK (x; h)2nhn(4)whereZ K(z)2 dzR(K) (for theoretical background see, e.g., Scott 1992, 130). A simple extension of (4) to theadaptive kernel density method isaVe {fbK(x; h)} 1 a1aR(K)fbK(x; h) fbK(x; h)2nhλ(x)n(5)where λ(x) denotes the bandwidth factor at point x. However, note that (5) understatesthe true variance since the local bandwidth factors are assumed fixed.Probability weights can be taken into account by adding a penalty for the amountof variability in the weights distribution. In particular,PnVe {fbK (x; h)} 2i 1 wiW2 1R(K)fbK (x; h) fbK (x; h)2h (6)where W is the sum of weights as defined above (see Van Kerm 2003 and Burkhauseret al. 1999 for a similar approach). The assumption behind this formula, however, isthat the weights w are essentially independent from X, an assertion that may not beappropriate.1Note that (4) differs from the standard variance formula often found in the literature.The standard formula only contains the first term (see, e.g., Van Kerm 2003, Silverman1986, 40, Härdle et al. 2004) and, although the second term asymptotically disappears,has a quite substantial bias in finite samples. Estimator (4) is more accurate than thestandard formula in finite samples.1. A point could also be made that weights should be omitted from estimation entirely if the assertionof independence is, in fact, true.

Ben Jann43Exact variance estimationThe variance of fbK (x; h) can be written as12E [Kh (x X) E{Kh (x X)}]n 1 E{Kh (x X)2 } E{Kh (x X)}2nV {fbK (x; h)} (7)where Kh (z) 1/h K(z/h) andZ Kh (x y)f (y) dyE{Kh (x X)} E{fbK (x; h)} A natural estimators for (7) is1Vb {fbK (x; h)} n( 2nx Xi1 X wiK fbK (x; h)2W i 1 h2h)or, in the case of the adaptive method,() 2n1x Xi1 X wiaa2bbbV {fK (x; h)} K fK (x; h)n W i 1 (hλi )2hλi(8)(9)where w is assumed to represent frequency weights or analytic weights (for similarformulas see, e.g., Hall 1992 and Fiorio 2004). Again note that (9) will be downwardbiased because the local bandwidth factors are assumed fixed.If the weights are sampling weights, an estimator for (7) may be derived as 2nx Xi1 X 2 1bbbwKV {fK (x; h)} 2 fK (x; h)W i 1 i hh(10)In many cases the approximate estimator is quite good and using the exact formulais not worth the extra computational effort. However, the exact formula should be usedif the data contain sampling weights. Furthermore, note that both the exact and theapproximate variance formulas assume h fixed. Data dependent choice of h, however,may result in additional variability of the density estimate, especially in regions withhigh curvature.5Confidence intervalsPointwise confidence intervals for f (x) are constructed asqfbK (x; h) z1 α/2 Vb {fbK (x; h)}(11)

4Univariate kernel density estimationwhere z1 α/2 is the (1 α/2) quantile of the standard normal distribution.Confidence intervals such as (11) may have bad coverage due to the well-known biasin fbK . One suggestion to improve coverage is to use an undersmoothed density estimateto construct the confidence intervals (Hall 1992, Fiorio 2004). The implementation ofthis approach in kdens usesn1/5hus h τnwhere τ is the undersmoothing parameter. τ should be larger than 1/5; τ 1/4 is thedefault choice in kdens once undersmoothing is requested.6Density estimation for bounded variablesTwo simple methods to estimate the density of variables with bounded domain are therenormalization method and the reflection method. Both methods produce consistentestimates, but have a bias of order h near the boundaries (compared to the usual biasof order h2 in the interior). Various more advanced correction techniques with biasof order h2 at the boundaries have been proposed (see, e.g., Jones and Foster 1996or Karunamuni and Alberts 2005). One of these methods is the linear combinationtechnique discussed in Jones (1993).6.1RenormalizationThe most natural way to deal with the boundary problem is to use a standard kerneldensity estimate and then locally rescale it relative to the amount of local kernel massthat lies within the support of X (see, e.g., Jones 1993, 137). Let L be the lowerboundary of the support of X (e.g. L 0) and U be the upper boundary (e.g. U 1).Furthermore, letZua0 (l, u) K(y) dylThe renormalization version of the standard estimator then is1n fbK (x; h)fbK(x; h, L, U ) U xa0 L xh , h(12)for x [L, U ]. Furthermore, the boundary renormalization analogue to the adaptiveestimator in (2) can be written as n1 X1wix Xin,ab fK (x; h, L, U ) K(13)W i 1 a L x , U x hλihλi0 hλihλinThe approximate variance estimator for fbKis L x U x1 b h , h1 nnn2Ve {fbK(x; h, L, U )} fb (x; h, L, U ) fbK (x; h, L, U )nh a0 L x , U x 2 Knhh(14)

Ben Jann5whereZb(l, u) uK(y)2 dyl(see Jones 1993). The exact variance estimator can be written asnVb {fbK(x; h, L, U )} 1b {fbK (x; h)} VU x 2,a0 L xhh(15)in the most simple case. The renormalization variance estimators for the adaptive kernelmethod and for data containing sampling weights can also be easily derived, but theresulting formulas are more complicated.6.2ReflectionThe reflection estimator approaches the boundary problem by “reflecting” the dataat the boundaries (see, e.g., Silverman 1986, 30, Ćwik and Mielniczuk 1993). In thestandard case the reflection estimator is given asn1 X wi rrK (x; Xi , h, L, U )fbK(x; h, L, U ) W i 1 h(16)for x [L, U ], whereK r (x; X, h, L, U ) K x Xh Kx X 2Lh Kx X 2Uh The reflection technique, that is, replacing K with K r , can be applied analogouslyto the adaptive estimator and also the exact variance estimators are easily derived.Unfortunately, however, the reflection solution for the approximate variance estimatoris more complex and not supported by kdens.6.3Linear CombinationLet z (x X)/h, l (L x)/h, and u (U x)/h. The linear combination techniquethen replaces K(z) byK lc (z; l, u) whereZa2 (l, u) a1 ( u, l)zK(z)a2 (l, u)a0 (l, u) (a1 ( u, l))2 la1 ( u, l) ZyK(y) dy, ua2 (l, u) (17)uy 2 K(y) dylSimilar to the reflection technique, exact variance estimates can be obtained by simply plugging K lc into the standard formulas. Approximate variance estimation is notsupported.

6Univariate kernel density estimation7Binned estimation7.1DensityKernel density estimators such as (1) involve lots of computations. One solution toreduce processing time is to estimate the density based on binned data. The binnedkernel density estimator is defined as m1 X c gj g efK (gj ; h) K,Whhj 1, . . . , m(18) 1where g1 , . . . , gm is a grid of m Pequally spaced evaluation points and c1 , . . . , cm arethe associated grid counts withc W . Given the equal spacing of grid points,equation (18) has a discrete convolution structure and can be calculated using fastFourier transform (see Wand and Jones 1995, 182–188, for details). This makes theestimator very fast. The grid counts are computed using linear binning as follows.Let g and g be the two nearest grid points below and above observation Xi . Thenwi (g Xi )/(g g ) is added to the grid count at g and wi (Xi g )/(g g )is added to the count at g . Note that the results from (18) are usually quite accurateeven for relatively small m. The rule-of-thumb given by Hall and Wand (1996, 182) isthat ,,grid sizes of about 400–500 are adequate for a wide range of practical situations“.The binned version of the adaptive kernel density estimator is m1 X c gj g aefK (gj ; h) K,Whλ hλ j 1, . . . , m(19) 1where λ denotes the local bandwidth factors. Unfortunately, the computational shortcutused for (18) is not applicable to (19).7.2VarianceVariance estimation is straightforward with binned data. For example,11Ve {feK (gj ; h)} R(K)feK (gj ; h) feK (gj ; h)2nhn(20)If sampling weights are applied, a reasonable variance formula for the binned estimatoris 2m1 X1gj g Vb {feK (gj ; h)} 2K feK (gj ; h)(21)c (w2 )Whh 1with c(w2 ) representing linearly binned squared weights.

Ben Jann7.37Density derivatives and density functionalsAdvanced automatic bandwidth selection involves estimating density functionals of theformZZR(r) R(f (r) ) f (r) (z)2 dz ( 1)r f (2r) (z)f (z) dzwhere f (r) denotes the rth derivative of f . A binned approximation estimator employingthe gaussian kernel can be written asmX(2r)e(r) (m, h) ( 1)r 1cj feφ (gj ; h)RφW j 1where m1 X c (r) gj g (r)eφfφ (gj ; h) ,Whr 1hj 1, . . . , m(22)(23) 1and φ(r) denotes the rth derivative of φ, the standard normal density. Equation (23)can be solved as the convolution of fast Fourier transforms.7.4Bounded variablesBinned versions of the estimators for bounded variables are usually simple to derive.For example, the binned renormalization density estimator in the standard case isnfeK(gj ; h, L, U ) 1IK L gj U gjh ,h feK (gj ; h),j 1, . . . , m(24)for gj [L, U ].An exception is the estimation of density derivatives and density functionals wherethe renormalization and the linear combination methods have no easy solutions. Fortunately, the reflection technique is a valuable alternative in this situation.8Data-dependent bandwidth selectionIt can be shown from asymptotic theory that the bandwidth 15R(K)hopt 2 }2 R(f 00 )n{σKwhere2σK Zz 2 K(z) dz,ZR(K) (25){K(z)}2 dzis “optimal” in the sense that it minimizes the asymptotic mean integrated squared2error (AMISE). Note that σK, the kernel variance, and R(K), the kernel “roughness”,are known properties of the chosen kernel function. However, R(f 00 ), where f 00 denotesthe second derivative of f , is unknown.