arXiv:1611.08049v2 [math.ST] 10 Dec 2016

A new kernel estimator of hazard ratio and its
asymptotic mean squared error
Taku MORIYAMA
Graduate School of Mathematics
Yoshihiko MAESONO
Faculty of Mathematics

Kyushu University, Motooka 744, Fukuoka 819-0395, Japan
Abstract
The hazard function is a ratio of a density and survival function, and
it is a basic tool of the survival analysis. In this paper we propose a
kernel estimator of the hazard ratio function, which are based on a modification of Ćwik and Mielniczuk [2]’s method. We study nonparametric
estimators of the hazard function and compare those estimators by means
of asymptotic mean squared error (AM SE). We obtain asymptotic bias
and variance of the new estimator, and compare them with a naive estimator. The asymptotic variance of the new estimator is always smaller
than the naive estimator’s, so we also discuss an improvement of AM SE
using Terrell and Scott [10]’s bias reduction method. The new modified
estimator ensures the non-negativity, and we demonstrate the numerical
improvement.

Keywords: Kernel estimator; hazard ratio; nonparametric estimator; mean
squared error

1

Introduction

Rosenblatt [8] has proposed a kernel estimator of a probability density function
f (·). After that, many researchers have discussed various kernel estimators for
distribution, regression, hazard functions etc. Since most of the kernel type
estimators are biased, many papers have studied bias reductions. Those are
based on a higher order kernel, transformations of estimators etc. Terrell and
Scott [10] have proposed efficient transformation. We will apply their method to
the kernel estimator in Section 4. On the other hand, Ćwik and Mielniczuk [2]
have discussed a new kernel estimator of a density ratio which they called ‘direct’
estimator. Asymptotic mean squared error (AM SE) of the direct estimator is

1

different from a naive estimator, and the variance is always smaller. Here let us
introduce their direct estimator.
Let X1 , X2 , · · · , Xn be independently and identically distributed (i.i.d.) random variables with a distribution function F (·), and Y1 , Y2 , · · · , Yn be i.i.d. random variables with a distribution function G(·). f (·) and g(·) are the density
functions of F (·) and G(·), and we assume g(x0 ) 6= 0 (x0 ∈ R). A naive estimator of f (x0 )/g(x0 ) is defined as replacing f (x0 ) and g(x0 ) by the kernel density
estimators fb(x0 ) and gb(x0 ) where
n

1 X
fb(x0 ) =
K
nh i=1



x0 − Xi
h



n

and

1 X
gb(x0 ) =
K
nh i=1



x0 − Yi
h



.

K(·) is a kernel function and h is a bandwidth which satisfies h → 0 and
nh → ∞ (n → ∞). We call fb(x0 )/b
g(x0 ) ‘indirect’ estimator. The direct type
estimator has been proposed by Ćwik and Mielniczuk [2], which is given by


Z
fb
1 ∞
Gn (x0 ) − Gn (y)
dFn (y),
(x0 ) =
K
g
h −∞
h

where Fn (·) and Gn (·) are the empirical distribution functions of X1 , · · · , Xn
and Y1 , · · · , Yn , respectively. Chen, Hsu, and Liaw [1] have obtained an explicit
form of its AM SE.
Extending their method, we propose a new ‘direct’ estimator of the hazard
ratio function in Section 2 and investigate its AM SE. We will compare the
naive and direct kernel estimators in Section 3. It is demonstrated that our
direct estimator performs asymptotically better, especially when the underlying distribution is exponential or gamma. Though the bias term of the direct
estimator is large in some cases, the asymptotic variance is always smaller. As
mentioned before, there are many bias reduction methods and then we discuss
the bias reduction in Section 4. Applying the Terrell and Scott [10]’s method,
we will show that the modified direct estimator gives us good performance both
theoretically and numerically. Proofs of them are given in Appendices.

2
2.1

Hazard ratio estimators and asymptotic properties
Hazard ratio estimators

Let us assume that the density function f (·) of Xi satisfies f (x0 ) 6= 0 (x0 ∈ R).
The hazard ratio function is defined as
H(x) =

f (x)
.
1 − F (x)

The meaning of H(x)dx is conditional probability of ‘death’ in [x, x+dx] given a
survival to x. This is a fundamental measure in survival analysis, and the range
2

of applications is wide. For example, actuaries call it “force of mortality” and use
it to estimate insurance money. In extreme value theory, it is called “intensity
function” and used to determine the form of an extreme value distribution (see
Gumbel [3])
The naive nonparametric estimator of H(x0 ) is given by Watson [12]
e 0) =
H(x

where
n

1 X
fb(x0 ) =
K
nh i=1



x0 − Xi
h



fb(x0 )
,
1 − Fb (x0 )

and

n

1X
W
Fb (x0 ) =
n i=1



x0 − Xi
h



.

Here K(·) is the kernel function and W (·) is an integral of K(·)
Z u
W (u) =
K(t)dt.
−∞

Using properties of the kernel density estimator, Murthy [5] proves the consise 0 ). In random censorship model, Tanner,
tency and asymptotic normality of H(x
Wong et al. [9] proves them, using the Hájek’s projection method. Patil [6] gives
its mean integrated squared error (M ISE) and discusses the optimal bandwidth
in both uncensored and censored settings. For dependent data, M SE of the indirect estimator is obtained by Quintela-del Rı́o [7]. Using Vieu [11]’s results,
he has obtained modified M ISE which avoids that the denominator equal to 0.
In this paper, we assume that the support of the kernel K(·) is [−d, d] (d > 0)
and the data do not have a censoring.
Extending the idea of Ćwik and Mielniczuk [2], we propose a new ‘direct’
estimator of the hazard ratio function as follows:

Z ∞ 
y − ηn (y) − (x0 − ηn (x0 ))
b 0) = 1
dFn (y),
H(x
K
h −∞
h
where

ηn (y) =

Z

y

Fn (u)du = n−1

−∞

X

I(Xi ≤ y)(y − Xi ).

b
It is easy to see that H(x)
is a smooth function. We will discuss its asymptotic
properties.

2.2

Asymptotic properties

For brevity, we use a notation
Ai,j =

Z

∞

K i (u)uj du.

−∞

The proofs of the following theorems are given in Appendix. For the direct
hazard estimator, we have the following AM SE.
3

Theorem 1 Let us assume that (i) f (·) is 4 times continuously differentiable
at x0 , (ii) K is symmetric and the support is [−d, d] for some positive number
b 0 ) is given
d > 0 and (iii) A1,4 and A2,0 are bounded. Then, the M SE of H(x
by



2

f
b
E H(x0 ) −
(x0 )
1−F




h4 2
A2,0
f
{(1 − F ){(1 − F )f ′′ + 4f f ′ } + 3f 3 }2
=
(x(1)
(x0 ) +
A
0)
4 1,2
(1 − F )10
nh 1 − F


1
+O h6 +
.
nh1/2

e 0 ) is
On the other hand, under some regularity conditions, M SE of H(x
given by Patil [6] as follows:

=




2
f
e
E H(x0 ) −
(x0 )
1−F




h4 2
A2,0
{(1 − F )f ′′ + f f ′ }2
f
(x0 ) +
(x0 )
A
4 1,2
(1 − F )4
nh (1 − F )2


1
.
+O h6 +
nh1/2

(2)

The asymptotic variances are second terms of (1) and (2), and then the direct
estimator has small variance because of 0 < 1−F (x0 ) < 1. Taking the derivative
b 0 ),
of the AM SE with respect to h, we have an optimal bandwidth h = h∗ of H(x
where
!1/5


A2,0
(1 − F )9 f
∗
−1/5
(x0 )
.
h =n
A21,2 {(1 − F ){(1 − F )f ′′ + 4f f ′ } + 3f 3 }2
Though h∗ depends on unknown functions, it is available by replacing them by
their estimators (plug-in method). Similarly, the following optimal bandwidth
e 0 ) is obtained
of the indirect H(x
h

∗∗

=n

−1/5

!1/5


A2,0
(1 − F )2 f
(x0 )
.
A21,2 {(1 − F )f ′′ + f f ′ }2

b
Further we can show the asymptotic normality of the direct H.

Theorem 2 Let us assume the same conditions of Theorem 1 and h = cn−d (0 <
b 0 ) holds.
c, 15 ≤ d < 21 . Then, the following asymptotic normality of H(x




√
f
d
2
b
nh H(x0 ) −
→ N (0, V1 ),
(x0 ) − h B1 −
1−F
4

where


A1,2 (1 − F ){(1 − F )f ′′ + 4f f ′ } + 3f 3
(x0 ) and
2
(1 − F )5

B1 =

V1 = A2,0




f
(x0 ).
1−F

Using the Slutsky’s theorem, the asymptotic normality of the indirect estimator
is easily obtained.
The direct estimator is superior in the sense of the asymptotic variance, so we
will consider the bias reduction in Section 4. We have the following higher-order
asymptotic bias.
Theorem 3 Let us assume that (i′ ) f is 6 times continuously differentiable at
x0 , (ii′ ) K is symmetric and the support is [−d, d] for some positive number
d > 0 and (iii′ ) A1,6 is bounded. Then, the higher-order asymptotic bias of
b 0 ) is
H(x





f
b 0) −
E H(x
(x0 ) = h2 a2 (x0 ) + h4 a4 (x0 ) + O h6 + n−1 ,
1−F
where

a2 (x0 ) =
a4 (x0 ) =



A1,2 −m2 m′′′ + 4mm′ m′′ − 3(m′ )3
(x0 ),
2
m5

A1,4 −60m2 (m′ )2 m′′′ + 15m3 m′′ m′′′ + 11m3 m′ m(4) − m4 m(5)
24
m9

210m(m′ )3 m′′ − 73m2 m′ (m′′ )2 − 105(m′ )5
(x0 )
+
m9

and m(x) = 1 − F (x).

3

Comparison of kernel hazard estimators

b 0)
In some special cases, we will investigate the AM SE values of the direct H(x
e
b
and indirect H(x0 ). We show that the proposed estimator H(x0 ) performs
asymptotically better when F (·) is the exponential or gamma distribution.
Here we consider that F (·) is the exponential, uniform, gamma, weibull
or beta distribution. The cumulative distribution function of the exponential
distribution Exp(1/λ) is F (x) = 1 − exp(−λx) and the hazard ratio is constant,
that is H(x) = λ. It is one of the most common and an important model of the
e 0)
survival analysis. When F (·) is the exponential, the asymptotic biases of H(x
b
and H(x0 ) vanish and the AM SEs are
h
i
h
i
b 0 ) = λ A2,0 < AM SE H(x
e 0 ) = λ exp(λx0 )A2,0 .
AM SE H(x
nh
nh

Thus the proposed estimator is always asymptotically better regardless of the
parameter λ and the point x0 .
5

Next, let us assume that F (·) is the uniform distribution (F (x) = x/b (0 <
x < b)), then the hazard ratio is H(x) = (b − x)−1 . In this model, the hazard
ratio increases drastically in the tail area. These AM SEs are the followings
i
h
b 0) =
AM SE H(x

h
i
e 0) =
AM SE H(x

h4 2
1
9b4
1
A2,0
A1,2
+
10
4
(b − x0 )
nh b − x0
b
1
A2,0 .
nh (b − x0 )2

e 0 ) vanishes and the variance of H(x
b 0)
We find that the asymptotic bias of H(x
is smaller. Their asymptotic performance depends on x0 and b, but AM SE of
b 0 ) is smaller when the life span b is large.
the proposed H(x
Last, we consider that F is the gamma Γ(p, 100), weibull W (q, 100) or beta
distribution 100 × B(r, s) where p, q, r and s are their shape parameters. Their
scale σ = 100 is moderate as life span. Γ(p, σ) is the distribution of a sum
of p i.i.d. random variables of Exp(σ), and so it is most important case. Its
asymptotic squared bias, variance and AM SE values for some fixed points x0
[p]
[p]
eε represent
are found in Table 1, where we omit power terms of h. γ
bε and γ
b 0 ) and H(x
e 0 ) when x0 is ε-th quantile of Γ(p, 100). The
those values of H(x
kernel is the Epanechnikov with A1,2 = 1/5, A2,0 = 3/10 and h = n−1/5 . The
coefficients n−4/5 are omitted here.
The weibull distribution W (q, σ) is also important in survival analysis because the hazard ratio is proportional to polynomial degree (q − 1), that is
H(x) = qσ q xq−1 . W (1, σ) is the exponential distribution. The beta distribution is often used to describe distribution whose support is finite, and it has
plentiful shapes. Table 2, 3 and 4 give each of the least AM SE values using
[q]
[r,s]
b 0)
h∗ or h∗∗ (in Section 2.2) where wε and βε
stand AM SE values of H(x
e 0 ) when x0 is ε-th quantile of W (q, 100) or B(r, s) × 100. We demonand H(x
b performs asymptotically better in most cases
strate the proposed estimator H
of the gamma Γ(p, 100). In the weibull cases, the asymptotic performance of
our estimator is good and comparable in the beta cases.

4

Bias reduction and simulation study

As discussed in Section 3, the direct estimator has good performance, especially
small variance. If we could reduce the bias, we can get a better estimator. There
are many bias reduction methods, but the variance reduction is quite difficult.
b 0 ). The bias term is
In this section, we discuss reduction of the bias of H(x
complicated, but if we use 4-th order kernel, we have A1,2 = 0. Thus we can
reduce the convergence order of the bias from O(h2 ) to O(h4 ). A simple way to
construct the 4-th order kernel is proposed by Jones [4]. However the estimator
takes negative value in some cases, though the hazard ratio is always nonnegative. On the other hand, Terrell and Scott [10]’s method reduces asymptotic
bias without loss of non-negativity. We propose the following modified direct

6

Table 1: Asymptotic Bias2 , V ar and AM SE when F is gamma
[1/2]

b0.05
γ
[1/2]
γ0.05
e
[1/2]
γ0.25
b
[1/2]
γ0.25
e
[1/2]
γ0.75
b
[1/2]
γ0.75
e
[1/2]
γ0.95
b
[1/2]
γ0.95
e
[10]

γ
b0.05
[10]
γ
e0.05
[10]
γ
b0.25
[10]
γ
e0.25
[10]
γ
b0.75
[10]
γ
e0.75
[10]
γ
b0.95
[10]
γ
e0.95

Bias2
7.22×10−2
6.53×10−2
1.33×10−8
7.74×10−9
5.35×10−13
6.00×10−14
1.01×10−12
6.93×10−16

V ar
4.01×10−2
4.22×10−2
9.52×10−3
1.27×10−2
4.29×10−3
1.72×10−2
3.58×10−3
7.16×10−2

AM SE
0.112
0.107
9.52×10−3
1.27×10−2
4.29×10−3
1.72×10−2
3.58×10−3
7.16×10−2

2.49×10−18
7.16×10−19
2.61×10−18
2.40×10−18
2.84×10−16
1.05×10−17
1.84×10−13
8.70×10−18

1.56×10−4
1.64×10−4
4.77×10−4
6.36×10−4
1.07×10−3
4.28×10−3
1.45×10−3
2.90×10−2

1.56×10−4
1.64×10−4
4.77×10−4
6.36×10−4
1.07×10−3
4.28×10−3
1.45×10−3
2.90×10−2

[1/2]

b0.1
γ
[1/2]
γ0.1
e
[1/2]
γ0.5
b
[1/2]
γ0.5
e
[1/2]
γ0.9
b
[1/2]
γ0.9
e
[1/2]
γ0.975
b
[1/2]
γ0.975
e
[10]

b0.1
γ
[10]
γ0.1
e
[10]
γ0.5
b
[10]
γ0.5
e
[10]
γ0.9
b
[10]
γ0.9
e
[10]
γ0.975
b
[10]
γ0.975
e

Bias2
8.24×10−5
6.72×10−5
2.36×10−11
6.83×10−12
1.69×10−15
2.93×10−15
1.47×10−11
2.33×10−16

V ar
2.11×10−2
2.33×10−2
5.65×10−3
1.13×10−2
3.76×10−3
3.76×10−2
3.46×10−3
0.139

AM SE
2.11×10−2
2.34×10−2
5.65×10−3
1.13×10−2
3.76×10−3
3.76×10−2
3.46×10−3
0.139

2.16×10−18
1.73×10−21
1.37×10−17
7.91×10−18
1.19×10−14
9.83×10−18
2.75×10−12
7.58×10−18

2.55×10−4
2.83×10−4
7.72×10−4
1.54×10−3
1.32×10−3
1.32×10−2
1.56×10−3
6.24×10−2

2.55×10−4
2.83×10−4
7.72×10−4
1.54×10−3
1.32×10−3
1.32×10−2
1.56×10−3
6.24×10−2

Table 2: AM SE values when F is gamma and h = h∗ or h∗∗
[1/2]

b0.05
γ
[1/2]
γ0.05
e
[1/2]
γ0.75
b
[1/2]
γ0.75
e
[10]

γ
b0.05
[10]
γ
e0.05
[10]
γ
b0.75
[10]
γ
e0.75

AM SE
7.44×10−2
7.60×10−2
7.40×10−5
1.45×10−4
4.48×10−7
3.64×10−7
5.39×10−6
8.46×10−6

[1/2]

b0.1
γ
[1/2]
γ0.1
e
[1/2]
γ0.9
b
[1/2]
γ0.9
e
[10]

b0.1
γ
[10]
γ0.1
e
[10]
γ0.9
b
[10]
γ0.9
e

AM SE
1.14×10−2
1.19×10−2
2.11×10−5
1.48×10−4

[1/2]

γ
b0.25
[1/2]
γ
e0.25
[1/2]
γ
b0.95
[1/2]
γ
e0.95
[10]

6.45×10−7
1.68×10−7
1.34×10−5
2.05×10−5

γ
b0.25
[10]
γ
e0.25
[10]
γ
b0.95
[10]
γ
e0.95

AM SE
1.06×10−3
1.20×10−3
7.26×10−5
1.86×10−4
1.11×10−6
1.37×10−6
2.51×10−5
3.76×10−5

[1/2]

b0.5
γ
[1/2]
γ0.5
e
[1/2]
γ
b0.975
[1/2]
γ
e0.975
[10]

b0.5
γ
[10]
γ0.5
e
[10]
γ
b0.975
[10]
γ
e0.975

AM SE
1.97×10−4
2.67×10−4
1.21×10−4
2.54×10−4
2.26×10−6
3.53×10−6
4.57×10−5
6.75×10−5

Table 3: AM SE values when F is weibull and h = h∗ or h∗∗

[1/2]

w
b0.05
[1/2]
w
e0.05
[1/2]
w
b0.75
[1/2]
w
e0.75
[10]

w
b0.05
[10]
w
e0.05
[10]
w
b0.75
[10]
w
e0.75

AM SE
4.11×10−2
4.20×10−2
9.22×10−6
1.53×10−5
1.20×10−4
1.11×10−4
1.38×10−2
1.33×10−2

[1/2]

w
b0.1
[1/2]
w
e0.1
[1/2]
w
b0.9
[1/2]
w
e0.9
[10]

w
b0.1
[10]
w
e0.1
[10]
w
b0.9
[10]
w
e0.9

AM SE
5.68×10−3
5.93×10−3
3.31×10−6
8.61×10−6

[1/2]

w
b0.25
[1/2]
w
e0.25
[1/2]
w
b0.95
[1/2]
w
e0.95
[10]

2.65×10−4
2.23×10−4
5.49×10−2
5.96×10−2

w
b0.25
[10]
w
e0.25
[10]
w
b0.95
[10]
w
e0.95

7

AM SE
3.85×10−4
4.30×10−4
3.59×10−7
7.68×10−6
9.00×10−4
4.79×10−4
0.135
0.153

[1/2]

w
b0.5
[1/2]
w
e0.5
[1/2]
w
b0.975
[1/2]
w
e0.975
[10]

w
b0.5
[10]
w
e0.5
[10]
w
b0.975
[10]
w
e0.975

AM SE
4.25×10−5
5.50×10−5
2.67×10−6
7.90×10−6
3.40×10−3
2.34×10−3
0.309
0.360

Table 4: AM SE values when F is beta and h = h∗ or h∗∗

[1/2,1/2]
βb0.05
[1/2,1/2]
e
β0.05
[1/2,1/2]
βb0.75
[1/2,1/2]
βe
0.75
[2,5]
βb0.05
[2,5]
βe0.05
[2,5]
βb0.75
[2,5]
e
β0.75
[5,2]
βb0.05
[5,2]
βe0.05
[5,2]
b
β0.75
[5,2]
e
β0.75

AM SE
7.61×10−3
7.77×10−3
2.35×10−3
1.61×10−3
1.55×10−4
1.83×10−4
3.09×10−3
1.84×10−3
7.18×10−5
6.57×10−5
9.66×10−3
4.23×10−3

hazard estimator

[1/2,1/2]
βb0.1
[1/2,1/2]
e
β0.1
[1/2,1/2]
βb0.9
[1/2,1/2]
βe
0.9
[2,5]
βb0.1
[2,5]
βe0.1
[2,5]
βb0.9
[2,5]
e
β0.9
[5,2]
βb0.1
[5,2]
βe0.1
[5,2]
b
β0.9
[5,2]
e
β0.9

AM SE
1.20×10−3
1.25×10−3
0.167
0.116

[1/2,1/2]
βb0.25
[1/2,1/2]
e
β0.25
[1/2,1/2]
βb0.95
[1/2,1/2]
βe
0.95
[2,5]
βb0.25
[2,5]
βe0.25
[2,5]
βb0.95
[2,5]
e
β0.95
[5,2]
βb0.25
[5,2]
βe0.25
[5,2]
b
β0.95
[5,2]
e
β0.95

2.70×10−4
2.18×10−4
1.01×10−2
7.22×10−3
1.38×10−4
1.16×10−4
6.68×10−2
4.01×10−2

AM SE
1.42×10−4
1.54×10−4
4.57
3.17
5.35×10−4
2.91×10−4
2.41×10−2
1.88×10−2
4.20×10−4
2.84×10−4
0.261
0.171

[1/2,1/2]
βb0.5
[1/2,1/2]
e
β0.5
[1/2,1/2]
βb0.975
[1/2,1/2]
βe
0.975
[2,5]
βb0.5
[2,5]
βe0.5
[2,5]
βb0.975
[2,5]
βe0.975
[5,2]
βb0.5
[5,2]
βe0.5
[5,2]
b
β0.975
[5,2]
βe0.975

AM SE
1.46×10−4
1.11×10−4
127
87.9
1.15×10−3
5.67×10−4
5.70×10−2
4.72×10−2
1.72×10−3
6.73×10−4
0.981
0.676

b I (x) = {H
b h (x)} 43 {H
b 2h (x)}− 31 ,
H

b h (·) and H
b 2h (·) are the direct estimator with the bandwidth h and 2h,
where H
b I (x0 )
respectively. By [10] we can get an asymptotic bias of H
2a22 (x0 ) − 4a4 (x0 )H(x0 ) 4
h ,
H(x0 )

where a2 (x0 ) and a4 (x0 ) are given in Theorem 2.3. Whereas, the order of the
variance is


4b
1b
I
b
V ar[H (x0 )] = V ar Hh (x0 ) − H2h (x0 ) + O(n−1 ) = O((nh)−1 ),
3
3

and then the order is same as the original. The order of the optimal bandwidth
is O(h−1/9 ), and the optimal convergence rate is improved to O(n−8/9 ). We
confirm this in numerical simulations.
Table 5 shows the results of 10,000 times simulation of improvement rate
(%) of
b I ]/M ISE[H]}),
b
1 − {M ISE[H

when F (·) is the beta B(r, s) and a sample size is n = 100. Their kernel functions
b is the asymptotically optimal
are the Epanechnikov. The bandwidth h† of H
b
bI
in the sense of M ISE[H] which is available from Theorem 2.1 and that of H
† 5/9
is h
modified to be optimal in the sense of the convergence rate. Table 5
b I (x) is superior. To obtain
shows that how well the bias reduction works, and H
b I (x)] is quite difficult and so
an explicit form of the asymptotic variance V ar[H
we postpone it in future work.

8

b I ] - F is B(r, s) ∗ σ , h = h∗ or h∗∗
Table 5: Improvement rate (%) of M ISE[H
Improvement rate (%)
r = 1/2, s = 1/2
r = 1, s = 1
r = 3, s = 3
r = 2, s = 5
r = 5, s = 2

σ = 4−1
0.301
0.593
6.74
12.9
3.96

σ=1
1.25
1.05
8.35
13.9
5.10

σ=4
1.52
1.24
8.27
15.0
5.11

σ = 42
1.64
1.38
8.28
14.8
5.05

σ = 43
1.78
1.30
8.35
15.0
5.05

References
[1] Chen, S.M., Hsu, Y.S., and Liaw, J.T. (2009). On kernel estimators of
density ratio. Statistics 43, 463–479.
[2] Ćwik, J., and Mielniczuk, J. (1989). Estimating density ratio with application to discriminant analysis. Communications in Statistics-Theory and
Methods, 18, 3057–3069.
[3] Gumbel, E. (1958). Statistics of extremes. Columbia Univ. press, New York.
[4] Jones, M., and Signorini, D. (1997). A comparison of higher-order bias
kernel density estimators. Journal of the American Statistical Association
92, 1063–1073.
[5] Murthy, V. (1965). Estimation of jumps, reliability and hazard rate. The
Annals of Mathematical Statistics, 36, 1032–1040.
[6] Patil, P. (1993). On the least squares cross-validation bandwidth in hazard
rate estimation. The Annals of Statistics, 21, 1792–1810.
[7] Quintela-del Rı́o, A. (2007). Plug-in bandwidth selection in kernel hazard
estimation from dependent data. Computational Statistics & Data Analysis,
51, 5800–5812.
[8] Rosenblatt, M. (1956). Remarks on some nonparametric estimates of a
density function. The Annals of Mathematical Statistics, 27, 832–837.
[9] Tanner, M.A., Wong, W.H., et al. (1983). The estimation of the hazard
function from randomly censored data by the kernel method. The Annals
of Statistics, 11, 989–993.
[10] Terrell, G.R., and Scott, D.W. (1980). On improving convergence rates for
nonnegative kernel density estimators. The Annals of Statistics, 8, 1160–
1163.
[11] Vieu, P. (1991). Quadratic errors for nonparametric estimates under dependence. Journal of Multivariate Analysis, 39, 324–347.
[12] Watson, G., and Leadbetter, M. (1964). Hazard analysis II. Sankhyā: The
Indian Journal of Statistics, Series A, 26, 101–116.
9

5

Appendices: Proofs of Theorems

Proof of Theorem 2.1
For the sake of simplicity, we use the following notations
Z z
Z z
η(z) =
F (u)du,
γ(z) =
η(u)du,
−∞

M (z) = z − η(z),

−∞

m(z) = M ′ (z) = 1 − F (z).

To begin with, we consider the following stochastic expansion of the direct
estimator
b 0)
H(x


Z
1
M (y) − M (x0 )
=
dFn (y)
K
h
h


Z
1
M (y) − M (x0 )
′
+ 2
{[η(y) − ηn (y)] − [η(x0 ) − ηn (x0 )]}dFn (y)
K
h
h


Z
M (y) − M (x0 )
1
{[η(y) − ηn (y)] − [η(x0 ) − ηn (x0 )]}2 dFn (y)
K ′′
+ 3
h
h
+···
= J1 + J2 + J3 + · · ·
(say).
Since J1 is a sum of i.i.d. random variables, the expectation is obtained
directly


 Z

1
M (y) − M (x0 )
E[J1 ] = E
dFn (y)
K
h
h


Z
1
M (y) − M (x0 )
=
f (y)dy
K
h
h


Z
f
(M −1 (M (x0 ) + hu))du
=
K(u)
1−F




f
h2
(1 − F ){(1 − F )f ′′ + 4f f ′ } + 3f 3
=
(x0 ) + A1,2
(x0 ) + O(h4 ).
1−F
2
(1 − F )5
Combining the following second moment


Z
1
M (y) − M (x0 )
2
f (y)dy
K
nh2
h


Z
f
1
(M −1 (M (x0 ) + hu))du
K 2 (u)
=
nh
1−F
1
f
=
(x0 )A2,0 + O(n−1 ),
nh 1 − F
we have the variance as follows



f
1
(x0 )A2,0 + O(n−1 ).
V [J1 ] =
nh 1 − F
10

Next, we consider the following representation of J2


n
n
1 X X ′ M (Xi ) − M (x0 )
J2 = 2 2
Q(Xi , Xj ),
K
n h i=1 j=1
h
where
Q(xi , xj ) = [η(xi ) − J(xi − xj )] − [η(x0 ) − J(x0 − xj )]
and
J(y − Xi ) = I(Xi ≤ y)(y − Xi ).
Using the conditional expectation, we get the following equation
E[J2 ] =

=
=
=

=



 
n
M (Xi ) − M (x0 )
1 X
′
Q(Xi , Xj )
E K
nh2 j=1
h





n
X
1
M (Xi ) − M (x0 )
E
Q(Xi , Xj ) Xi 
E K ′
nh2
h
j=1


 
1
M (Xi ) − M (x0 )
′
{η(Xi ) − [η(x0 ) − J(x0 − Xi )]}
E K
nh2
h
Z

1
K ′ (u) η(M −1 (M (x0 ) + hu)) − η(x0 ) + J(x0 − M −1 (M (x0 ) + hu))
nh
f
(M −1 (M (x0 ) + hu))du
×
1−F
 
Z
f
1
1
′
.
(x0 )du = O
K (u) O(hu)
nh
1−F
n

Next we have
J22

=





n
n
n
n
M (Xk ) − M (x0 )
1 X X X X ′ M (Xi ) − M (x0 )
′
K
K
n4 h4 i=1 j=1
h
h
k=1 l=1

=

1
n4 h4

×Q(Xi , Xj )Q(Xk , Xl )
n X
n X
n X
n
X
J22 (i, j, k, l)

(say).

i=1 j=1 k=1 l=1

Taking a conditional expectation, we obtain that if all of (i, j, k, l) are different,

 
E[J22 (i, j, k, l)] = E E J22 (i, j, k, l)|Xi , Xk = 0,
and

E[J22 (i, j, k, l)] = 0
E[J22 (i, j, k, l)] = 0

(if i = j and all of (i, k, l) are different),
(if i = k and all of (i, j, l) are different),

E[J22 (i, j, k, l)] = 0

(if i = l and all of (i, j, k) are different).
11

Then the terms in which ‘j = l and all of (i, j, k) are different’ are the main
term of E[J22 ]. If j = l and all of (i, j, k) are different, we have
E[J22 (i, j, k, l)]



 
n(n − 1)(n − 2)
M (Xk ) − M (x0 )
M (Xi ) − M (x0 )
′
′
K
E K
n4 h4
h
h

×Q(Xi , Xj )Q(Xk , Xj ) .

=

Using the conditional expectation of Q(Xi , Xj )Q(Xk , Xj ) given Xi and Xk , we
find
oi
h n
E E Q(Xi , Xj )Q(Xk , Xj ) Xi , Xk
h
= E η(Xi )η(x0 ) + η(Xk )η(x0 ) − η 2 (x0 ) + 2γ(x0 ) − η(Xi )η(Xk )
−(x + Xi − 2 min(x, Xi ))η(min(x, Xi )) − 2γ(min(x, Xi ))
−(x + Xk − 2 min(x, Xk ))η(min(x, Xk )) − 2γ(min(x, Xk ))

i
+(Xi + Xk − 2 min(Xi , Xk ))η(min(Xi , Xk )) + 2γ(min(Xi , Xk )) .

Therefore, the entire expectation of the last row is



 
M (Xk ) − M (x0 )
M (Xi ) − M (x0 )
K′
E K′
h
h

=

i
×(Xi + Xk − 2 min(Xi , Xk ))η(min(Xi , Xk )) + 2γ(min(Xi , Xk ))




Z "Z y
M (y) − M (x0 )
M (w) − M (x0 )
′
′
K
K
h
h
∞
+

Z

y

∞

K′



× {(−w + y)η(w) + 2γ(w)} f (y)dw



M (y) − M (x0 )
M (w) − M (x0 )
′
K
h
h
#

× {(w − y)η(y) + 2γ(y)} f (w)dw f (y)dy.

12

Finally we get
 



M (Xi ) − M (x0 )
M (Xk ) − M (x0 )
′
′
E K
K
h
h

i
×{Xi + Xk − 2 min(Xi , Xk )}η(min(Xi , Xk )) + 2γ(min(Xi , Xk ))


Z
M (y) − M (x0 )
f (y)dy
= −h2 K ′
h
"




f
M (y) − M (x0 )
{η(x0 ) + (−x0 + y)F (x0 )}
(x0 )
×
W
h
m2
#
 ′


f m − f m′
+{(−x0 + y)η(x0 ) + 2γ(x0 )}
(x0 ) +O(h)
m3
"


M (y) − M (x0 )
+
1−W
h
#!


 ′

f m − f m′
f
(x0 ) + O(h)
× η(y) 2 (x0 ) + {(x0 − y)η(y) + 2γ(y)}
m
m3




 ′
Ff2
f
f ′ m − f m′
f m − f m′
= h4
2η
(x
)
+ O(h5 ).
(x
)
+
+
2γ
0
0
m4
m3
m2
m3

After similar calculations of the other terms, if j = l and all of (i, j, k) are
different, we get
 
1
E[J22 (i, j, k, l)] = O
.
n
In addition, it is easy to see that the expectations of the other combinations of
(i, j, k, l) are o(n−1 ). As a result, we have
E[J22 ] = O(n−1 ) and V [J2 ] = O(n−1 ).
The moments of J3 is also obtained in a similar manner and we find that
b 0 ) comes from J1 . By CauchyJ3 + · · · is negligible. The main bias of H(x
Schwarz inequality, we find that the main term of the variance comes from J1 .
Now, we get the AM SE of the direct estimator and the proof of Theorem 2.1.
Proof of Theorem 2.2 and 2.3
It follows from the previous discussion that




√
f
b
nh H(x0 ) −
(x0 )
1−F




√
f
=
nh J1 −
(x0 ) + oP (1)
1−F




√
f
2
(x0 ) − h B1 + oP (1).
= B1 + nh J1 −
1−F
13

Since J1 is a sum of i.i.d. random variables and the expectation of the second
term is o(1), the asymptotic normality of Theorem 2.2 holds. Furthermore,
these arguments show that
h
i
b 0 ) = E[J1 ] + O(n−1 )
E H(x


Z
f
(M −1 (M (x0 ) + hu))du + O(n−1 ),
=
K(u)
1−F
and we get directly Theorem 2.3 by the Taylor expansion.

14

