Nonparametric Threshold Detection
for Cost Distributions
Benedikt Funke
TH Köln

Masayuki Hirukaway
Ryukoku University

4th April 2024

Abstract
The tail of a distribution beyond some threshold is of importance and interest in academics and industries. In this paper, we develop a nonparametricthreshold detection method designed specially for cost distributions. Our approach originates from tailoring the existing techniques for change point or jump
location detection in statistics to stylized facts of cost distributions. Considering that the threshold is located in the right-tail region, we propose to employ
the asymmetric gamma kernel when constructing the diagnostic function. It isdemonstrated that our jump location estimator is consistent with a faster convergence rate than the parametric one and asymptotically normal when suitably
implemented. Since this estimator tends to underestimate the jump location
in …nite samples, we advocate correcting its bias. Monte Carlo simulations
and real data examples illustrate attractive properties and practical relevanceof our proposal in several di¤erent use cases.
Keywords: distribution of costs; gamma kernel; heavy-tailedness; incomplete gamma function; cross validation; threshold detection.
JEL Classi…cation Codes: C13; C14; G22.
MSC 2010 Codes: 62G07; 62G20; 62P05; 91B30.

Corresponding Author ; Institute for Insurance Studies, TH Köln - University of Applied Sciences,
Gustav-Heinemann-Ufer 54, 50968 Köln, Germany; e-mail: benedikt.funke@th-koeln.de.
yFaculty of Economics, Ryukoku University, 67 Tsukamoto-cho, Fukakusa, Fushimi-ku, Kyoto
612-8577, Japan; e-mail: hirukawa@econ.ryukoku.ac.jp.

1

Introduction

The tail of a distribution beyond some threshold gathers special interest in academics
and in several industrial sectors. It is a source of large market ‡uctuations or evenvolatility clustering in …nance. The region in an income distribution corresponds tothe ‘elite income’ tier. In the …eld of non-life insurance, a handful of severe losses
within a collection of policies beyond some threshold can signi…cantly contribute to
the total claim amount. In geophysics and hydrology, tail analyses are conductedfrom the viewpoint of catastrophic disasters.
At the same time, because distributions in these examples behave quite di¤erently
below and above a threshold, it is also widely recognized that a single model cannot
capture characteristics over the while ranges of such distributions. In fact, the tailpart of an income distribution above a threshold is modelled and studied separatelyfrom the part that is below the threshold and contains most of observations (called
the bulk part hereinafter); see Cowell, Ferreira and Litch…eld (1998) and Jenkins(2017), for instance. Moreover, modelling the whole range of a loss distribution
is of particular importance and interest in actuarial science, particularly in non-life
insurance. Then, splicing (e.g., Klugman et al., 2019) or composite modelling (e.g.,
Cooray and Ananda, 2005; Scollnik and Sun, 2012) -- the practice of modelling theparts below and above a threshold di¤erently -- is frequently employed.
Several research questions arise in this context. For example, where does thetail part begin in an income distribution?

Or where should we switch to the tailmodel in a loss distribution? The aim of this paper is to propose a nonparametricmethod of …nding a threshold in these distributions. Incomes and actuarial losses,
as well as wages and consumption expenditures, are examples of cost variables. Ourestimation strategy is grounded on prior knowledge (or stylized facts) about shapes of
distributions of such variables. Distributions of cost variables have support on R+ witha natural boundary at the origin and are highly right-skewed with a concentration ofobservations in the vicinity of the origin and a long tail with sparse data. As revealed
shortly, we employ a speci…c kernel to capture these characteristics while avoiding
possible misspeci…cation by relying on a particular parametric model.
In addition, we translate the problem of threshold estimation into that of change
point (or jump location) detection, which has been actively studied in statistics.
Our translation or interpretation of the problem is motivated by the original idea of1

splicing. In this idea, two di¤erent probability density functions (pdfs) are …tted forthe regions below and above a threshold or a splicing point, and continuity of the
distribution at this point is not required.1 Accordingly, the words “threshold” and“jump”are used exchangeably hereinafter, whenever no confusion may arise.
Following the literature on change point detection (e.g., Chu and Cheng, 1996;
Couallier, 1999; Huh, 2002), we take the absolute di¤erence of two kernel density
estimates as the diagnostic function and de…ne its maximizer as the jump location
estimator.2 Nonetheless, our procedure remarkably di¤ers from previous ones. Whilethey are designed to work in the middle part (i.e., near the peak) of a distribution,a threshold that divides the bulk and tail parts is supposed to be located on theright tail and far away from the peak of the distribution. To overcome such dif…culty while exploiting prior information on shapes of cost distributions, we adopt
the asymmetric gamma kernel by Chen (2000); our contribution is the …rst work inwhich an asymmetric kernel is applied for the problem of jump detection in density
or regression curves, to the best of our knowledge. The kernel is de…ned as the pdf
of the gamma distribution G (x=b + 1; b) for a design point x (0) and a smoothing

parameter b (> 0), and it is in the form ofKG(x;b) (u) =where

ux=b exp ( u=b)
1 fu
bx=b+1 (x=b + 1)

0g ;

(1)

( ) is the gamma function, and 1 f g is an indicator function. This kernel

has an adaptive smoothing property via changing shapes automatically across design
points, and a single smoothing parameter su¢ ces to generate a variety of shapes.
Thanks to this property, the gamma kernel can calibrate right-skewed densities withsupport on R+ (e.g., those of cost distributions) well.
p

Our jump location estimator is super-consistent (i.e., its convergence rate exceeds
n, where n is the sample size) and asymptotically normal when suitably implemen-ted. While the estimator shares its asymptotic properties with the existing literature (e.g., Chu and Cheng, 1996; Couallier, 1999; Huh, 2002), our proof strategy is1

Many authors (e.g., Cooray and Ananda, 2005; Scollnik and Sun, 2012) additionally imposedi¤erentiability of the pdf at the splicing point to make the entire density smooth and reduce thenumber of parameters. This practice has computational advantage now that the splicing point canbe expressed as a function of other model parameters. However, such a restriction on parameters
results in less ‡exibility. From this viewpoint, Reynkens et al. (2017) …rst estimate the splicing
point and then compute estimates of remaining model parameters.
2
An alternative approach is proposed by Desmet et al. (2010), who transform kernel density
estimation into kernel regression estimation via prebinning and applying an existing method of
discontinuity detection for nonparametric regression curves. This approach has the disadvantagethat prebinning substantially reduces the sample size which can be used for regression estimation.

2

totally di¤erent.

The gamma kernel does not admit the location-scale transform-

ation Kh (u x) = K f(u x) =hg =h with a bandwidth h (> 0) or exchangeability
between u and x, unlike standard symmetric kernels. Instead, approximations to
the incomplete gamma, digamma and polygamma functions studied by Funke andHirukawa (2019, 2023) are tailored for the technical proofs. This proof strategy isnovel and of independent interest; see the Appendix and/or the Supplemental Material for more details. It is also demonstrated that a uniform approximation to ourdiagnostic function has a unique maximum. So far existence of the maximum in the
diagnostic function has been simply suggested, not formally proven, in the literature,
although it is a key ingredient for consistency of the jump location estimator. A
concern is that the jump location estimator tends to generate negative biases in …nite
samples. However, an elementary bias correction can instantly improve bias properties of the estimator without in‡ating its variance, and thus we advocate putting thebias-corrected version to practical use.
Acknowledging that threshold estimation is a notoriously di¢ cult problem, many
authors have proposed various threshold detection methods, mainly in the …eld of
non-life insurance. The methods are divided roughly into three categories, namely,(i) heuristic approaches, (ii) graphical diagnostics and (iii) automated procedures.
For (i), a threshold is de…ned as a …xed quantile (DuMouchel, 1983) or determined bysome formula depending on the sample size (Loretan and Phillips, 1994); see Scarrottand MacDonald (2012) for more details. Despite no theoretical justi…cation, thesemethods are used by actuaries in practical applications. Examples of (ii) include
the Hill plot and its variants (e.g., Resnick 1997), the mean excess plot (Davison andSmith, 1990), and the peaks-over-threshold method (e.g., Cunnane, 1973).

These

are easy to grasp and thus used regularly, whereas there is room for practitioners’
discretion at the stage of identifying a threshold. Recent research has shifted toward(iii) such as the minimum Kolmogorov-Smirnov (KS) distance procedure (Clauset etal., 2009; Drees et al., 2020), sequential goodness-of-…t testing (e.g., Bader et al.,2018), and the minimum quantile discrepancy and automated Eye-Balling methods(Danielsson et al., 2019). These procedures circumvent arbitrariness but rely on acertain parametric model of the tail part including the generalized Pareto distribution(GPD). Our threshold detection procedure will be able to serve as a more objectiveand ‡exible alternative to these methods.
Furthermore, there will be potentially many applications of our proposal. The
de…nition of a ‘cost’is not strict, actually. Our threshold estimation procedure is expected to work equally for non-cost variables (e.g., quantities demanded, transaction
3

volumes, etc.), as long as shapes of their distributions have similarities to those ofcosts. It can be also employed for threshold detection prerequisite for tracing out theevidence of illegal trading (James et al., 2023) and implementing extreme changes inchanges (CIC) estimation (Sasaki and Wang, 2023), for instance.
The remainder of this paper is organized as follows.

Section 2 overviews thejump location estimation using the gamma kernel and then recommends its practicalimplementation. In Section 3, convergence properties of the jump location estimator, namely, strong consistency and asymptotic normality, are explored.

Section4 conducts Monte Carlo simulations to compare …nite-sample behaviors of our jump
location estimator with those of several existing threshold or jump detection methods.
Our aim is to show numerically the advantage of our proposal over these alternatives.
In Section 5, the jump location estimation is applied to a couple of real world datasets. Section 6 concludes. Proofs of theorems and propositions are provided in the
Appendix. Proofs of lemmata are deferred to the Supplementary Material, which isavailable on the second author’s webpage.
This paper adopts the following notational conventions: ‘an

bn ’ means thatan =bn converges to 1; ‘an = o (bn )’signi…es that an =bn converges to 0; ‘an = O (bn )’means that an =bn is bounded; and we say that ‘an
bn ’ if there exist constants
0 < c1 < c2 < 1 so that c1 an
bn
c2 an . For a function h (x) and a pointc, h(c ) = limx"c h(x), h(c+ ) = limx#c h(x) and h(m) (x) = dm h (x) =dxm denote theleft and right limits, and the mth-order derivative, respectively. The abbreviationd‘a:s:’stands for “almost surely”. Finally, the expression ‘X = Y ’reads “A randomvariable X obeys the distribution Y .”2

Our Proposal: An Informal Overview

2.1

Estimation of the Jump Location

It is suspected that f (x), the pdf of a ‘cost’variable X 2 R+ , is discontinuous at t0

on a prespeci…ed closed interval I0 := t; t with 0 < t < t < 1. Throughout it isassumed that the interval I0 is located in the right tail part of the underlying cost

distribution. Prior knowledge on the interval is not at all unrealistic, because quite
often practitioners have a rough idea about the location of the threshold through,for example, preliminary threshold estimates, historical experiences and/or empiricalquantiles. Against this background, our method can complement existing practices
and objectify the graphical analysis procedures mentioned in Section 1.
The model below basically follows those of Chu and Cheng (1996) and Couallier4

(1999).

A similar local structure can be also found in jump detection problems

for nonparametric regression (e.g., Wu and Chu, 1993a,b; Joo and Qiu, 2009) anddeconvolution (e.g., Delaigle and Gijbels, 2006). It is assumed that the pdf f (x) forx 2 I0 can be modelled locally asf (x) = g (x) + d0 1 fx < t0 g ;(2)where g (x) is a su¢ ciently smooth function, and t0 2 I0 is the discontinuity point(or jump location). The jump size is given by
f t+
0 ;

d0 = f t0
where jd0 j 2 (0; 1) is assumed throughout.

Our problem is how to estimate the jump location t0 nonparametrically. Supposethat there are n i.i.d. observations fXi gni=1 at hand. If t0 were located in the bulk
region, as in the existing literature, splitting the entire sample into two sub-samples
near t0 would cause no serious issue. Two sample sizes are roughly the same, and thus
both left and right limits of a density can be estimated equally well. This is clearly
not the case in our problem. Because t0 is located in the tail region, sample-splittingnear t0 results in imbalance in sample sizes of two sub-samples and an imprecise
density estimate from the right sub-sample. In view of this, we use the entire sample
to estimate both limits of the pdf. To do so, we introduce ‘shifted’gamma kernels
KG(x;b;i.e.,) ( ), which are de…ned as pdfs of gamma distributions G f(xKG(x;b;) (u) :=

u(x
b(x

)=b+1

) =b + 1; bg,)=bexp ( u=b)
1 fu
f(x) =b + 1g

where b (= bn > 0) is the smoothing parameter,0g ;n > 0) plays the role of a shift(=parameter, and each parameter shrinks towards zero at a certain rate.

Obviously,) ( ) collapse to Chen’s (2000) original gamma kernel (1) when

= 0. The

KG(x;b;kernels can be interpreted as those designed to smooth the data o¤ the target design
point x by a margin of . In addition, they put the maximum weight at sightly leftor right of x because they have their modes at x.

We exploit as the source of our jump detection the di¤erence between two density
estimates that is created by the shift parameter . Let the shifted density estimators
be
n
1X
^
^
f (x) := fb; (x) :=KG(x;b; ) (Xi ) :n i=1

Also de…neJ^ (x) := f^ (x)5f^+ (x) ;where a single, common value is chosen for the smoothing parameter b in both density
estimates, as in Chu and Cheng (1996), Couallier (1999) and Huh (2002). Following
these articles, we also utilize J^ (x) as the diagnostic function for threshold detection.
The estimator of the jump location t0 , denoted as t^, is de…ned as the maximizer ofJ^ (x) on x 2 I0 , i.e.,t^ := arg max J^ (x) :
x2I0

2.2

Recommended Estimation Procedure

Monte Carlo results in Section 4 indicate that t^ tends to underestimate the jump
location. However, it turns out that the negative bias can be alleviated substantially,
with no additional cost of spread, by an elementary bias correction. In practice, we
recommend the bias-corrected version of the jump location estimatort~ := t^ + b
for a suitably chosen smoothing parameter b. A theoretical foundation of the bias
correction can be found in Remark 2, and its necessity is also visualized in Figure 1.
Superior …nite-sample properties of t~ over t^ are con…rmed in Section 4.
We can compute t~ in the following steps:
1. Prespecify the interval I0 that is likely to cover the jump location t0 .
2. Put the shift parameter= b for= 0:70 and select the smoothing parameterb via the modi…ed likelihood cross-validation method given by (10)-(11). Boththe exponentin

and the choice method for b are based on our judgmentsfrom the Monte Carlo study and real data examples.
3. Find the maximizer of J^ (x) on x 2 I0 and take it as t^.
4. Obtain the bias-corrected estimate t~ = t^+ b using the value of b selected in Step
2.

3

Large-Sample Properties of the Jump Location
Estimator

In this section convergence properties of the jump location estimator t^ are documented. Our particular focus is on its consistency and asymptotic normality. In thecourse of this, we demonstrate existence of a unique maximum in a certain uniform
6

approximation to J^ (x) on x 2 I0 , which constitutes a key condition for consistencyof t^. The asymptotic distribution of t^ also hints that a simple form of its leading
bias enables us to derive the bias-corrected estimator t~.

3.1

Regularity Conditions

Convergence results below rely on the fact that J^ (x) can be approximated by thedi¤erence between two incomplete gamma functions.

To deliver the results, we

impose the following regularity conditions.
Assumption 1. fXi gni=1 2 R+ are i.i.d. random variables.
Assumption 2.
(i) The pdf f (x) is uniformly bounded on x 2 R+ .
(ii) The local structure (2) holds and g (3) (x) is Lipschitz continuous and bounded
on x 2 I0 .
Assumption 3. Tuning parameters b and
b3=4
and

+

b

+
1=2

ln n
nb3=2

for somesatisfy b;! 0,b1=2 ln n
! 0;
n 2

(3)

= O (1)(4)

2 [0; 1), as n ! 1.

All these assumptions are standard for uniform approximations to asymmetrickernel estimators. Similar conditions can be found, for example, in Hirukawa et al.
(2022). It follows from Assumption 2(ii) that f (1) (t0 ) = f (1) (t+0 ). Boundedness of thethird-order derivative of the smoothedg (3) ( ) on I0 is also an importanton component
ingredient for approximations to E J^(p) (x) for p = 0; 1; 2. This type of condition
has been often imposed in simulation studies on jump detection (e.g., Wu and Chu,
1993a,b; Chu and Cheng, 1996).
Assumption 3 controls the shrinkage rates of tuning parameters b and . The
condition (3) draws the following important conclusions: (i) b = o ( ); (ii) b1=2 =o(2=b); and (iii)=o3=b3=2 . These are frequently used to controls remainder

terms in the asymptotic expansions. It also follows from b = o ( ) and7= o b1=2

that although the shift parametershould shrink to zero more slowly than thesmoothing parameter b, the convergence rate of must not be too slow (or must be
faster than b1=2 , to be more precise). Couallier (1999), for instance, also imposesa similar rate requirement.
holds.

The condition (3) also implies that ln n= nb1=2 ! 0

The other condition (4) is an additional technical requirement for strong
uniform consistency of t^.

3.2

Consistency

Below asymptotic properties of the jump location estimator t^ are explored. Our
analysis starts from a uniform approximation to f^ (x) on I0 , which is documented
in the next proposition. To save space, we adopt the following shorthand notation
whenever no confusion may arise: Kx (u) = KG(x;b; ) (u); Kx (u) = KG(x;b) (u);a = (x) =b; and z0 = t0 =b.

Proposition 1. If Assumptions 1-3 hold, then
Z t0
n
o
^
E f (x) = g (x) + d0
Kx (u) du g (1) (x)
0
n
o
x (2)
(1)
+ g (x) + g (x) b + o (b)2for every x 2 I0 , andsup f^ (x)
x2I0

!
r
n
o
ln
n
E f^ (x) = O
a:s:;
nb1=2

(5)

(6)

as n ! 1.
A direct outcome from Proposition 1 is that
sup J^ (x)x2I0

It also follows from J^ (x)

!
r
o
n
ln n
a:s:
E J^ (x) = O
nb1=2

n
o
E J^ (x)J^ (x)on
E J^ (x) that

!
r
o
n
ln
n
J^ (x) = E J^ (x) + O
a:s:
nb1=2

uniformly on I0 .

no^This reveals that E J (x) constitutes the dominant term inJ^ (x) . In other words, the e¤ect of the jump location on the magnitude of J^ (x)8noappears only in the value of E J^ (x) in a …rst-order asymptotic sense. This also

plays a key role in the proof of Theorem
1obelow; see the Appendix for more details.
nno^By (5) and Assumption 2(ii), E J (x) can be further approximated by E J^ (x) :=jd0 j J (x) + O ( ) uniformly on I0 , where
Z t0
Z t0
Kx+ (u) duKx (u) duJ (x) =
0
0
Z t0
Z t0
Kx+ (u) duKx (u) du=00= P a + 1; z0

P a+ + 1; z0 ;and P (a; z) :=(a; z) = (a) is a normalized version of the lower incomplete gammaRzfunction (a; z) = 0 ta 1 exp ( t) dt for a; z > 0. The reasons why P (a + 1; z0 )dP (a+ + 1; z0 ) holds are that P (a + 1; z0 ) = Pr (Yz0 ) for Y = G (a + 1; 1)(i.e., P (a + 1; z0 ) are cumulative distribution functions (cdfs) of Y

evaluated atz0 ) and that the larger the shape parameter is, the ‡atter the gamma distribution
becomes.
We are about to demonstrate strong consistency of t^.

Before proceeding, it iscurious whether J (x) on x 2 I0 indeed has a unique maximum at t0 (or within ashrinking neighborhood of t0 even if it is not maximized exactly at this point). In
reality, however, it is quite cumbersome to look into the local property of J (x) analytically. Fortunately, several approximations to the incomplete gamma function areavailable, and we rely on one of them.

More speci…cally, we employ equation (1)of Pagurova (1965) to approximate the normalized lower incomplete gamma functions P (a + 1; z0 ) around the standard normal cdf. The next proposition refers toproperties of the approximation and the maximizer of the approximated function.
Proposition 2. If Assumption 3 holds, then the followings hold true.
(i)J (x) can be further approximated by J (x) = Q (x)=b1=2 + O

formly on I0 as n ! 1, whereQ (x) :=

x + t0
x3=2

x t0
p
bx

and

( ) is the pdf of N (0; 1).

(ii)

Q (x) on I0 has a unique maximum at x = t 2 (t0FIGURE 1 ABOUT HERE
9b; t0 ).

3=b3=2 uni-

Before establishing strong consistency of t^, we show by some numerical illustrationthat maximizing J^ (x) is a well-de…ned problem. A discontinuous density f (x) andthe diagnostic function J^ (x) are drawn in Panel (a) of Figure 1. Model 1-A in
Section 4 is chosen for this illustration. The density is discontinuous at t0 = 4 withthe jump size d0 = 0:15. The diagnostic function is computed from a Monte Carlosample of sample size 500 under tuning parameters b = 0:05 and = b0:70 .
The uniform approximation J^ (x)jd0 j J (x) jd0 j Q (x) =b1=2 also tempts
us to make a visual inspection of shapes of these three curves. The curves aroundthe true threshold t0 = 4 are plotted in Panel (b) of Figure 1. Notice that the panel
magni…es the area surrounding t0 to visualize preciseness of the approximations. It
can be immediately found that all three curves are single-peaked around the true
threshold, which con…rms well-de…nedness of the optimization problem. Approximating J^ (x) by jd0 j J (x) looks decent, whereas the discrepancy between the two

curves suggests that the approximation errors which are asymptotically negligible
may not be ignored in …nite samples. Furthermore, jd0 j Q (x) =b1=2 approximates
jd0 j J (x) quite well; rather, they are almost indistinguishable. For reference, Panel(b) also indicates maximizers of J^ (x) and Q (x) are t^ 3:8693 and t
3:9502,
respectively. It follows from b = 0:05 that the latter con…rms Proposition 2(ii).
The function Q (x) has the following properties. Observe that
(
)
r
2
2
2
(t
=b)
(x
t
)
x + t0
t
=b
0
0
0
exp
b 1=2 Q (x) =
2
3
2 x
2t0 x
t0
where

p

(t20 =b) = (2 x3 ) exp(t20 =b) (x;t0 )2 = (2t20 x) is the pdf of the inverse Gaus-

sian distribution IG (t0 ; t20 =b). Because the shape parameter of this distribution
t20 =b ! 1, the pdf is close to a normal one for a su¢ ciently small b > 0. In addition,
the distribution has mean t0 and variance bt0 . It follows that the pdf roughly behaves

like N (t0 ; bt0 ), and thus, heuristically, the shape of Q (x) also looks like a bell curve
centered around t0 .
The theorem below delivers strong consistency of t^ for t0 .

The proof of thistheorem closely follows that of Theorem 1 in Chu and Cheng (1996); see the Appendix
for more details.3p
In principle, weak consistency t^ ! t0 su¢ ces, and some readers may wonder why Theorem 2.1 of
Newey and McFadden (1994) is not employed. Q (x) is uniquely maximized at x = t , I0 is compact,
and Q (x) is continuous. The problem is that while b1=2 =J^ (x) is uniformly approximated by
jd0 j Q (x), Q (x) still depends on n through b. Therefore, it is not possible to rely on the theorem
directly.
3

10

Theorem 1. If Assumptions 1-3 hold, then t^

3.3

t0 = O (b) a:s:, as n ! 1.

Asymptotic Normality

The theorem below documents asymptotic normality of t^. The asymptotic distribution is derived indirectly, as in Chu and Cheng (1996, Theorem 1), Couallier (1999,
Théorème 2) and Delaigle and Gijbels (2006, Theorem 3.1). The indirect derivation comes from the fact that t^ solves the …rst-order condition J^(1) t^ = 0. Then,a mean-value expansion of the left-hand side around t^ = t0 is made, and suitableapproximations to the incomplete gamma, digamma and polygamma functions are
utilized in the expansion; see the Supplemental Material for more details. This ispossible because unlike f (x), its estimates f^ (x) are smooth functions even at t0 due
to di¤erentiability of shifted gamma kernels Kx ( ) with respect to x.
Theorem 2. If Assumptions 1-3 hold, then
r

n
b1=2

t^

t0

d

( b) ! N (0; V0 ) := N

p 1=2
3 t0
0;
4d20

f (t0 ) + f (t+0)
2!

as n ! 1.
Remark 1. While it is di¢ cult to obtain asymptotic bias and variance of t^ in lightof the indirect nature, the asymptotic distribution in Theorem 2 implies the …rst two
moments of t^. The dominant bias term of t^ is b regardless of the position of t0 .
The expression of the term is much simpler than what is obtained by Couallier (1999,
Théorème 2). The di¤erence arises from di¤erent assumptions on the local structureof the pdf f on I0 ; invoke that our Assumption 2(ii) follows the setup by Chu andCheng (1996). As discussed in Remark 2 below, simplicity of the dominant bias
term enables us to make the bias correction of t^ straightforward. Moreover, V0 , thecoe¢ cient of the dominant variance term, suggests that the larger the magnitude of
discontinuity jd0 j, the easier the estimation of t0 . It can be also recognized that thefarther t0 moves away from the origin, the less precise its estimator becomes. Finally,both bias and variance terms are free of the shift parameter. A similar result isobtained in Théorème 2 of Couallier (1999); to put it another way,
convergence properties of t^ in a …rst-order asymptotic sense.

does not a¤ectRemark 2. As will be seen in the next section, t^ tends to yield negative biases in…nite samples, i.e., it is likely to underestimate the jump location, which coincides
11with what Theorem 2 predicts.

However, the theorem also suggests that the bias
can be corrected straightforwardly by adding b to t^. This is the foundation of thebias-corrected estimator t~ = t^ + b described in Section 2.2. Indeed, the proofs ofLemma A6 and Theorem 2 jointly imply that the leading bias of t~ is O (2), whereas

its variance is still O b1=2 =n . It will be con…rmed in the Monte Carlo study shortlythat t~ is a remedy for better …nite-sample properties.
Remark 3. An approximation to the mean squared error (AMSE) of t^ is
AM SE t^ = b2 +

b1=2
b1=2
V0 = O b2 +
n
n

;

(7)

where O (b2 ) and O b1=2 =n terms are leading squared bias and variance of t^, respectively. The AMSE for the jump location parameter implied by Théorème 2 ofCouallier (1999) is in the form of O (h4 + h=n), where h is the the bandwidth forstandard symmetric kernels.

It can be found that this AMSE and (7) are of thesame order of magnitude by recognizing that b

h2 . Furthermore, it follows from

(7) that no bias-variance trade-o¤ occurs, because a smaller b makes both squared
bias and variance terms smaller.
Remark 4. Some readers may wonder how to pick b andof t^. Then, putb for some
2 (1=2; 3=4) and b

for super-consistencynfor some
2(0; 1= (21=2)). Suchand b jointly satisfy (3) and draw the following threeconclusions on the convergence rate of t^:
1. For any

2 (1=2; 3=4), we can always put
1o (n ), or t^ becomes super-consistent.

2. Again for any> 1=2. In this case, AM SE t^ =2 (1=2; 3=4), we may even choose= 2=3, which balances or-ders of magnitude in the squared bias and variance so that O (b2 ) = O b1=2 =n =O n 4=3 . As a consequence, AM SE t^ = O n 4=3 . It is also clear that
the AMSE convergence rate of t^ is determined by the exponent . AM SE t^ =O (b2 ) (i.e., the squared bias dominates) for
2=3, and AM SE t^ = O b1=2 =n(i.e., the squared bias becomes asymptotically negligible) otherwise. The lattercase corresponds to an ‘undersmoothing’scenario so that nb3=2 ! 0 holds. As

a consequence, the asymptotic normality statement in Theorem 2 reduces topdn=b1=2 t^ t0 ! N (0; V0 ).
12

3. The best possible rate is AM SE t^ = O (n 2+" ) for an arbitrarily small " > 0.
The rate can be attained by setting and slightly above 1=2 and slightly
below 2, respectively. Chu and Cheng (1996) and Couallier (1999) also reportthat their jump location estimators can attain the same convergence rate underthe best case scenario.
Furthermore, it is not hard to see that for…nd some

and b de…ned above, we can always

2 [0; 1) satisfying (4). To see this, observe that (4) holds if nb3=2!1

at a polynomial rate. The rate requirement is attained for cases 1 and 2, for example,
by setting = 2=3 and choosing any 2 (0; 1). For case 3, slightly below 2 andslightly below 1 can jointly establish a polynomial divergence of nb3=2.

Remark 5. As long as f (x) can be locally modelled as or well-approximated by (2),both t^ and t~ become super-consistent when implemented as in Remark 4. It follows
that regardless of whether a parametric (e.g., the GPD) or nonparametric model (e.g.,Markovitch and Krieger, 2000) is …tted to the tail part, our jump location estimatorcan be safely used as a threshold estimate without deteriorating the convergence rate
for the tail model. In addition, Table 1 of Clauset et al. (2009) lists examples of
non-power law distributions that behave like the GPD. Clauset et al. (2009) evenargue that …tting a power law distribution in their procedure has nothing to do witha plausible match of the distribution with the data, and they recommend a goodnessof-…t test as a post-estimation analysis. Super-consistency of our estimators does noharm to convergence rates of the test statistics, either.

4

Finite-Sample Performance

4.1

Monte Carlo Design

We consider two alternative models in the simulation study below. In the …rst case,the univariate random variable X 2 R+ is drawn from a log-normal distribution,

which is an example of a heavy-tailed distribution and known as its power law-likebehavior (Clauset et al., 2009). What di¤ers from a usual log-normal distribution isthat a quadratic term is added to the pdf on the interval [0; t0 ) = [0; 4). Speci…cally,
the pdf f (x) is
f (x) =

1
1 + (2=3) Dt0

"

x

1
p

2

exp

(

(ln x
)2
2 2

13

)

#

+ (x) ; ( ; ) =

1 3
;
5 4

;

where(x) := D 1

t0 ) =t0 g2 1 fx < t0 g and the jump size d0 = f t0

f(x

f t+
= D= f1 + (2=3) Dt0 g. The shift parameter D takes two values, and D =0
1=4; 3=22 yield d0 = 0:15; 0:10, respectively. The former and latter cases are labelled
as “Model 1-A” and “Model 1-B”. It also follows from f (1) (t0 ) = f (1) (t+
0 ) that (2)is satis…ed in the neighborhood of the jump location t0 = 4.
This design is in some sense similar to that of Chu and Cheng (1996), who consider
the density f (x) constructed by splicing left and right sides of two normal distributions
with zero mean but di¤erent variances. In their design, f (0 ) 6= f (0+ ) but f (1) (0 ) =f (1) (0+ ), and thus (2) holds in the neighborhood of the origin (aside from the factthat the jump occurs in the middle part).
In the second case, the nonnegative random variable X is generated by some
distribution spliced at t0 = 4. The pdf f (x) in this case takes the general formf (x) = fL (x) 1 fx < t0 g + (1cL ) fR (x) 1 fx

t0 g ;

where fL (x) is some density function truncated at t0 , fR (x) is another density funcRt
Rttion with support on [t0 ; 1), and cL := 0 0 f (x) dx = 0 0 fL (x) dx ensures unity of
the integral of f (x) over its entire support R+ ; in other words, fL (x) and fR (x)represent bulk and tail models, respectively. This scenario is labelled as “Model 2”.
Throughout the Weibull distribution with densityfL (x) =

x

1

exp

n

x

o

; ( ; )=

3;

11
4

:

is considered as the bulk part of Model 2. Densities of the following distributions are
examined for the tail part, and two cases are denoted as “Model 2-A” and “Model
2-B”, depending on the tail model:8no (1+1= )(x t0 )
1>>
1+
1 fx t0 g ; ( ; s) = 14 ; 4[A: GPD]
< s
n so
n
o
k
fR (x) =x t0 +1 k 1
x t0 +1 kk
exp 1`
exp
1 fx t0 1g ;
>
`
`
`
>
:
(k; `) = 41 ; 1[B: Translated Weibull]
where f t0

= fL (3) and f t+
= (1
0

;

cL ) fR (3). Drees et al. (2022, p.83) argue

that discontinuity of the density at the threshold is an easy scenario for the thresholddetection method by Clauset et al. (2009). Model 2-A is most favorable to existing
threshold detection methods because of a power law distribution for the tail modellingas well as the discontinuity. The tail part in Model 2-B, also known as a stretchedexponential distribution, re‡ects that again it seemingly behaves like a power law
distribution (Clauset et al., 2009). Because this pdf is unbounded at the boundary
14

of the support t0

1, it is truncated at t0 . Moreover, these models violate (2) byconstruction but may be more realistic, because it is hard to judge whether the local
structure (2) indeed holds in real data. Our aim is to examine how our approach
can behave under such unfavorable circumstances. We simulate 1000 Monte Carlo
replications of fXi gni=1 with sample size n 2 f250; 500g from each model.

Table 1 presents the mode of the distribution of X, the constant cL , left and right
limits of the density at the discontinuity point f t0 , and the jump size d0 = f t0
f t+
0 . In each model, more than 95% of observations concentrate on the interval[0; t0 ) (i.e., in the bulk region), and the density has a long tail with a polynomial
decay rate. These features reasonably mimic properties of cost distributions.
TABLE 1 ABOUT HERE
Our estimation procedure for t0 is implemented as follows. There are two optim-

izations required, namely, (i) the one for tuning parameters (b; ) and (ii) the otherfor jump location search. For (i), Remark 4 suggests 2 (1=2; 3=4), and thus we
restrict our attention to four values, namely, 2 f0:55; 0:60; 0:65; 0:70g. A few cross-validation (CV) methods for b are investigated, and their details are deferred to thenext section. For each CV method, candidates of b are taken from 100 equally-spaced
grids over the interval [0:005; 0:500]. For (ii), after (b; ) are determined, the jump
location is searched via a numerical optimization routine for the diagnostic functionJ^ (x) on the interval I0 = [3; 5].4 Once the jump location estimator t^ is obtained
as the maximizer of J^ (x) , the bias-corrected estimator is computed as t~ = t^ + b.
Finite-sample performances of t^ and t~ are compared with those of existing (i)kernel-smoothed jump location estimation procedures and (ii) automated thresholddetection methods.

For (i), we focus on the procedure by Chu and Cheng (1996)[CC]. The CC diagnostic function is J^CC (x) = f^1 (x) f^2 (x) , where f^j (x) =P(nh) 1 ni=1 Kj f(Xi x) =hg ; j = 1; 2, for the kernels K1 and K2 to be speci…ed
shortly and a common bandwidth h (> 0). As in our method, the maximizer ofJ^CC (x) on I0 is de…ned as the jump location estimator. The kernels K1 and K2
4

The reason why di¤erent algorithms are utilized for (i) and (ii) is as follows. While a numerical
optimization routine substantially reduces computation time for (i), it often …nds local extrema
and corner solutions because of a high degree of nonlinearity in CV criteria. A grid search cancircumvent these issues. In contrast, J^ (x) is concave on I0 , and a numerical optimization routine
helps expedite computation for (ii).

15are fourth-order polynomial ones. These are
K1 (u) = 0:4857

3:8560u + 2:8262u2 + 19:1631u3 + 11:9952u4

1 fu 2 [ 1; 0:2012]g ;and K2 (u) = K1 ( u) for all u. These kernels are also employed for jump location
detection in nonparametric regression curves by Wu and Chu (1993a,b). The CC procedure is implemented as in ours. After the bandwidth value is found via grid search
for a CV criterion in the next section, the jump location is searched via numerical
optimization for the diagnostic function J^CC (x) on I0 .
For (ii), we investigate the followings: (a) the minimum KS distance procedure
between the empirical and GPD-based distribution functions by Clauset et al. (2009)[KS]; (b) the minimum quantile discrepancy criterion for the mean absolute deviation
between empirical and GPD-based quantiles by Danielsson et al. (2019) [Q-MAD]; (c)the minimum quantile discrepancy criterion for the sup-norm between empirical and
GPD-based quantiles by Danielsson et al. (2019) [Q-SUP]; and (d) the automated
Eye-Balling method based on tail index estimates by Danielsson et al. (2019) [AEB].
All simulations are conducted on R. In particular, R-packages “poweRlaw” and“tea” are employed to implement automated threshold detection methods (a) and(b)-(d), respectively.

4.2

Smoothing Parameter Selection

Selecting the smoothing parameter b is the most important practical issue. In ourcontext, values of (b; ) must be determined before jump location search so that thediagnostic function can be …xed on I0 . However, Remark 3 does not help resolve
this issue. There is no optimal choice for b on the basis of the bias-variance tradeo¤. Theorem 2 provides no guidance forguarantee that any, either, because it does not automaticallysatisfying (3) works equally well in …nite samples. Furthermore,to the best of our knowledge, there is no decisive conclusion on selecting the tuning
parameter in the context of jump location estimation; in fact, Chu and Cheng (1996)
adopt …xed bandwidths in their Monte Carlo study.
Taking the dependence ofon b into account, we tailor Huh’s (2012) approachto construct a few CV criterion functions. Before proceeding, put = b for a given. Accordingly, f^ (x) are rewritten as f^b (x; ), which signify the dependence ofdensity estimates on (b; ). Also letf^b; i (x; ) :=

1
n

n
X

1 j=1;j6=i
16

Kx (Xj )be density estimates using the sample with the ith observation eliminated. Finally,P
denote the number of observations falling into I0 as n0 := ni=1 1 fXi 2 I0 g.
The CV criterion function by Huh (2012) is de…ned as the sum of CV criteria for
two density estimates that construct the diagnostic function for jump location estimation. We incorporate this idea into three CV criterion functions. The minimizer of
each criterion function is taken as the corresponding CV smoothing parameter. The…rst one is the least-squares cross-validation (LSCV) criterion. It is de…ned as+
CVLS (b; ) = CVLS (b; ) + CVLS(b; ) ;

where

Z n
o2
CVLS (b; ) :=f^b (x; ) dx
I0

(8)

2 X ^
f (Xi ; ) :n0 i:X 2I b; i

(9)

0

i

The remaining two criteria are likelihood-based ones. One is the simple likelihood^ 2 (h) of Marron (1985) andcross-validation (LCV) criterion, which is analogous to L
equation (2.1) of Van Es (1991). It is given by
CVL (b; ) = CVL (b; ) + CVL+ (b; ) ;
where

X

CVL (b; ) :=i:Xi 2I0onln f^b; i (Xi ; )is the negative log-likelihood. The other is the modi…ed LCV criterion, denoted as^ 5 (h) of Marron (1985) and equation (2.2) of Van Es
MLCV, which corresponds to L
(1991). It takes the form of
CVM L (b; ) = CVM L (b; ) + CVM+L (b; ) ;

(10)

where
"

CVM L (b; ) :=
=

"
n
X
i=1

X

i:Xi 2I0

X

i:Xi 2I0

P

n Z
X

n
o
ln f^b; i (Xi ; )

i=1

n
o
ln f^b; i (Xi ; )

Xi
b

+ 1;

t
b

P

I0

Xi
b

KXi (u) du

+ 1;

t
b

#

#

;

(11)

and the second term is intended to eliminate the endpoint e¤ect of the interval
I0 = t; t .

Corresponding jump location estimates are labelled “SG-LS”, “SG-L”
17and “SG-ML”, where “SG” abbreviates “shifted gamma”. For their bias-correctedversions, we put “-BC”at the end.
Finally, choosing the bandwidth h is also required to implement the CC procedure.
The LSCV analogous to (8)-(9) (i.e., …nding a minimizer of the sum of two LSCV
criteria for f^1 (x) and f^2 (x)) is adopted.
TABLE 2 ABOUT HERE

4.3

Results

Table 2 presents several performance measures of jump location estimators. These
include the bias (“Bias”), standard deviation (“SD”) and root mean squared error(“RMSE”) of each threshold estimator over 1000 Monte Carlo samples. In addition,
for CC and SG estimators, Monte Carlo averages and standard deviations (in parentheses) of CV tuning parameters are reported for reference. Furthermore, there is no
guarantee that automated threshold detection methods necessarily yield threshold estimates falling into the interval I0 . For these methods, counts of threshold estimates
inside I0 out of 1000 Monte Carlo samples (“# fI0 g”) are provided.

We start from examining the results from automated threshold detection methods. It can be found that Q-MAD generates the smallest variance among all four
automated methods, resulting in the smallest RMSEs. In addition, more than 90%
of estimates from Q-MAD are inside I0 for each model and sample size, despite no restriction on the parameter space. While Models 2-A and 2-B are thought to be morefavorable than 1-A and 1-B for these automated methods, Q-MAD is comparable to
SG methods in terms of RMSE for the latter case. There are also general tendencies
of underestimation by KS and Q-MAD and overestimation by AEB. In particular,
the degree of overestimation by AEB is considerable.
CC also looks comparable to SG methods. It consistently overestimates the jump
location. However, the bias reduces with the sample size, and jointly by the decreasein dispersion, its RMSE becomes smaller as the sample size is larger.
We can …nd a general tendency in the results from SG methods. It can be immediately found that the initial estimate t^ tends to be negatively biased, as Theorem 2
predicts. Moreover, a short b makes both the bias and variance small, as suggested
in Remark 3.

Accordingly, it can be reasonably conjectured that a CV algorithm
that can generate a small smoothing parameter value ^b will contribute to the initial
estimate t^ with good quality.

Clearly, MLCV alone ful…lls this requirement. In
contrast, LSCV and LCV consistently give a large value of ^b, which leads to a con18

siderably biased t^. This is the source of their poor performance, and the sizable bias
cannot be corrected completely even after adding a large ^b to t^.
In what follows, we look into SG-ML and SG-ML-BC more carefully. Although
the shift parameter does not enter the asymptotic normality result in Theorem 2,
the choice of(or the exponent, to be more precise) matters in …nite samples.

For each model and sample size, smoothing parameter values via MLCV do not varymuch across four values of . On the other hand, = 0:70 (i.e., the smallest shift)
always results in the smallest RMSE among four initial SG-ML estimates. Because
bias correction is made by simply adding ^b to the initial estimate t^, additional variability through bias correction (i.e., estimation error of the correction term) is not
introduced. As a consequence, superiority of the smallest shift is maintained after
the bias correction. Indeed SG-ML-BC with = 0:70 outperforms others in termsof RMSE for Models 1-A and 1-B. Its RMSE is not the smallest for Models 2-A
and 2-B; an exception is Model 2-B with n = 500, to be more precise.

For thesecases, however, the RMSE is simply larger than those from other bias corrected SG
methods, and it is still smaller than those from CC and four automated methods.
Invoke that the local structure (2) holds in the neighborhood of the true jump
location for Models 1-A and 1-B, whereas it is violated for remaining two models.
It can be con…rmed from Monte Carlo results that SG-ML-BC with = 0:70 is of
most practical relevance and importance, because its RMSE is by far the smallest
in favorable scenarios and less than those from other competing methods even inunfavorable scenarios with deviation from (2).

5Real Data Examples

5.1

Data Description

In this section our jump location estimation procedure is applied to two datasets on
cost variables that are made publicly available. The one is of actuarial losses andthe other of wages. Below each dataset is discussed in detail.
The …rst dataset is of Danish …re insurance claims. Since the seminal analysis by
McNeal (1997), the dataset has been arguably most popularly chosen for empiricalstudies on non-life insurance. We extract the one named “danish” in R-package“SMPracticals”. The dataset contains 2; 492 …re insurance losses denominated inmillions of Danish kroner from 1980 to 1990.
The second dataset is taken from Merged Outgoing Rotation Group Earnings
Data of the US Current Population Survey (CPS), also known as CPS Labor Extract,19

which is available on the webpage of the National Bureau of Economic Research. We
extract hourly wages (the variable “earnhre”) earned by males from the dataset in
1979. Before proceeding, all observations denominated originally in US cents areconverted into US dollars.

The original sample size is 54; 769.

Considering thecomputation burden of kernel smoothing, we downsize the dataset to a sub-sample ofsample size 2; 709 by random sampling, which roughly accounts for 5% of the original.
Table 3 reports summary statistics of the datasets. Each distribution is rightskewed with a long tail and thus reasonably represents stylized facts of a ‘cost’distribution. It can be also found that the distribution of the downsized dataset of US
male wages well represents that of the original.
TABLES 3-4 ABOUT HERE

5.2

Estimation Details

Based on simulation results in Section 4, we compare SG-ML(-BC) with

= 0:70with KS, Q-MAD, Q-SUP, AEB, and CC. The tuning parameter for each of SG-ML(BC) and CC is taken from 100 equally-spaced grids over the interval [0:005; 0:500],
and then the threshold estimate t^ is searched via a numerical optimization routine of
the corresponding diagnostic function.
The prespeci…ed interval I0 for each dataset roughly covers the tail part of the
distribution. In addition, Cooray and Ananda (2005) and Scollnik and Sun (2012)
estimate several di¤erent parametric composite models from the Danish …re insurance
dataset and obtain threshold estimates ranging roughly from 1 through 3. Using
the same dataset, Reynkens et al. (2017) report the threshold estimate of 17 via agraphical method. The interval I0 for this dataset also incorporates these empirical…ndings.

5.3

Results

Table 4 presents the estimation results. For the Danish …re insurance data, there isno consensus in threshold estimates among four GPD-based detection methods. The
estimation result from SG-ML is fairly close to the one from KS, so is the one fromits bias-corrected version SG-ML-BC. In contrast, the result from CC is problematicand recognized as an example of estimation failure in that the threshold estimate isa corner solution.
On the other hand, threshold estimates from four GPD-based methods for the
US male wage data are similar one another. SG-ML and SG-ML-BC yield threshold20

estimates near 10, which is considerably close to the results from Q-MAD and Q-SUP.
Once again, CC may have failed to estimate the threshold, because its estimate lies
almost on the boundary.
As the location of the threshold is unknown in each distribution, it is hard tojudge among all threshold estimators considered. Nonetheless, it is safe to say that
SG-ML(-BC) can serve as a good alternative to existing threshold detection methods.

6

Conclusion

It is widely recognized that a single model cannot describe the while range of a cost
distribution well. Accordingly, it is of growing importance and interest to …nd thelocation of a threshold at which two di¤erent models are spliced, whereas this estimation problem is known to be notoriously di¢ cult. This paper has explored a methodof estimating the threshold of a cost distribution nonparametrically. Developmentof our method starts from tailoring the existing techniques for change point or jump
location detection in statistics to stylized facts of cost distributions. The diagnostic
function is the absolute di¤erence of two kernel density estimates, and the maximizerof the function over a prespeci…ed interval is de…ned as the jump location estimator.
Because the threshold is located in the right tail region, we propose to compute two
density estimates using the entire sample and shifted gamma kernels.

Our jump

location estimator is shown to be super-consistent and asymptotically normal whensuitably implemented. The proof strategy is also new in that approximations to
the incomplete gamma, digamma and polygamma functions are utilized. Since theestimator tends to underestimate the jump location and its dominant bias term takes
a simple form, we advocate correcting its bias in a straightforward manner. It is con…rmed in the Monte Carlo study that as a result of the bias correction, the negative
bias is substantially eliminated with no additional cost of spread. Several data-based
methods of choosing the smoothing parameter are also assessed via simulations. Our
proposal is …nally applied to three kinds of real world datasets. Judging from simulation results and real data examples, we recommend SG-ML-BC for practical use.

21

A

Appendix

The Appendix provides technical proofs of theorems and propositions.

To savespace, we defer proofs of lemmata to the Supplemental Material. Before proceeding,
additional short-handed notation is introduced, and a few useful formulae related to
the gamma function are presented.

A.1

Additional Notation(x) = d ln (x) =dx = (1) (x) = (x) and (m) (x) = dm (x) =dxm signify the digamma and polygamma functions, respectively. In addition, the following notationis adopted in the proofs: R (a; z) = z a exp ( z) = (a + 1) for a; z > 0; K_ c (u) =...
33• c (u) = @ 2 Kx (u) =@x2 j ; K@Kx (u) =@xjx=c ; K
c (u) = @ Kx (u) =@x jx=c ; Hi =x=c
K_ t (Xi ) K_ t+ (Xi ); and z = (t0
) =b.
0

A.2

0

Useful Formulae on the Gamma Function

Stirling’s Formula.
(a + 1) =p

2 aa+1=2 exp ( a) 1 +1+O a 212a

as a ! 1:(A1)

Recursive Formula for the Lower Incomplete Gamma Function.
z a exp ( z) for a; z > 0:

(a + 1; z) = a (a; z)(A2)

Recursive Formula for the Polygamma Function.
(m)

A.3

(a + 1) =

(m)

(a) +

( 1)m m!
for a > 0 and m 2 f0; 1; 2; :::g :am+1

Proof of Proposition 1

The proof requires the following lemmata.
Lemma A1.

As b ! 0,
sup
x2I0

sup
x2I0

x
x
b
+1
ln
+
b
b
2x
b
b
(1) x
+1
b
x 2x2
x
sup (m)
+1
b
x2I0
22

= O b2 ;= O b3 ; and= O (bm ) for m

2:

(A3)

o
n
1=2
+ t 1=2 . Then, as n ! 1,
Lemma A2. Let C0 := 2 max t 3=2 ; t
r
2
3=2
_
sup
K (u)
b
C0 :
x

(x;u)2I0 R+

Lemma A3. (Van der Vaart and Wellner, 1996, Lemma 2.2.9) Let X1 ; : : : ; Xnbe independent random variables with bounded ranges [ M; M ] and zero means.
Then,
Pr

n
X

Xi > xi=1

for all x and v
A.3.1

P
V ar ( ni=1 Xi ).

!

2 exp

x2
2 (v + M x=3)

Proof of Proposition 1

Proof of (5). Because
Z t0
o Z 1
n
^
Kx (u) g (u) du + d0
E f (x) =Kx (u) du;0

it su¢ ces to show thatZ 1
Kx (u) g (u) du = g (x)

0

g

(1)

(x)

0

n
o
x (2)
(1)
+ g (x) + g (x) b + o (b) :2(A4)

Observe that the left-hand side of (A4) can be expressed as E fg (X )g for Xd=

G (a + 1; b). A second-order Taylor expansion of g (X ) around X = x yields
E g X= g (x) + g (1) (x) E X

1
x + g (2) (x) E X
22+ Rg ;(A5)

g (2) (x) E (X

x)

x

2

where x is between X

and x, and Rg := (1=2) g (2) (x )is the remainder term.

Then, by the property of gamma random variables andAssumption 3, E (Xx) =+ b and E (X2x) = xb +2

3 b + 2b2 =xb + o (b). Rg = o (b) can be also shown by using Lipschitz continuity of g (2) ( ).
Then, (A4) is demonstrated by substituting these results into (A5).
Proof of (6). The proof for this part closely follows the one for Theorem 2 ofp
Hirukawa et al. (2022). For ease of exposition, let an := ln n= (nb1=2 ), & in (x) :=(1=n) [Kx (Xi ) E fKx (Xi )g], Nn := n1+ b 3=2 for a su¢ ciently small > 0, andjI0 j := tt. The proof for this part takes the following two steps.

1. Split the interval I0 into Nn equally-spaced grids to create Nn sub-intervals withlength Nn 1 jI0 j, and replace the supremum with a maximization over …nite Nn
sub-intervals.

23

2. Employ Lemma A3 (Bernstein’s inequality) to bound the remainder term.
Step 1. Let Ij ; j = 1; : : : ; Nn ; be the jth sub-interval. Also let xj be the rightmost point in Ij with x0

t and xNn

t. Suppose that the design point x falls into

Ij . By the mean-value theorem Lemma A2 and Assumption 3,
Kx (u)

Kxj (u)sup(x;u)2I0 R+

K_ x (u) sup jx

xj j

x2Ij

O b 3=2 O Nn 1 = O n (1+ )O (an ) :It follows from Cr -inequality thatnX

max sup

1 j Nn x2IjnX

& in (x)i=1& in (xj ) = O (an ) :

(A6)

i=1

Step 2. Before employing Bernstein’s inequality in Lemma A3, we must determine two bounds M and v. First, it holds that for a given x > (> 0),) 1=2 f1 + o (1)gp:2

b 1=2 (xKx (u)(A7)

This can be con…rmed by recognizing that Kx (u) is maximized at u = x

andemploying (A1). Then, by (A7),
sup

b 1=2

Kx (u)(x;u)2I0 R+

r

2
;
t

and thus, by Cr -inequality,& in (x)rr
2 1 1=2
2 a2n
2
=: M:
n b
=2
t
t ln nSecond, Assumption 2(i) implies that there is some constant C 2 (0; 1) so thatsupx2R+ f (x)C. Then,
( n
)
n
X
X
V ar
& in (x) =V ar & in (x)
i=1

i=1

=

n
X

i=1
n
X

E & in (x)E

1
K (Xi )n x

Z 1

Kx (u)

i=1

C
n

24

2

0

2

2

du:

Also letA (x) :=
Observe thatZ 1

Kx (u)

2

du = A (x)b 1 (2a + 1):
22a +1 2 (a + 1)

Z 1

u2a exp f u= (b=2)g(b=2)2a +1 (2a + 1)

0

0

(A8)

du = A (x) ;because the integrand in the middle term is the pdf of G (2a + 1; b=2). Using (A1),(x) 1=2 = x 1=2 f1 + o (1)g and o (1)A (x) =

b 1=2 (xuniformly on I0 . In sum,
( n
)
X
V ar
& in (x)1 for a su¢ ciently large n gives) 1=2 f1 + o (1)g
p
2

b 1=2
p
t

C
C a2n
p n 1 b 1=2 = p=: v:tt ln ni=1

Lemma A3 establishes that for such M and v and an arbitrarily chosen K > 0,
9
8
8
9
s
>
>
n
<X
=
=
<
C
K 2 ln nqPr:& in (x) > K p an
2 exp
q
:
>
t ;
;
: 2 1 + 2 2 Kan = pC >
i=1
3
t
t

q
p
p
It follows from an = o (1) that (2=3) 2= ( t)Kan = C= tciently large n. Accordingly,
8
9
s
n
<X
=
C
K 2 ln np
Pr
2 exp& in (x) > Kan:t ;2 (1 + 1)i=1

1 holds for a su¢ -= 2n

In the end,

Pr

8
<

n
X

s

9
=

C
p an
:1 j Nn
t ;
i=1
9
8
s
Nn
n
=
<X
X
C
max Pr& in (xj ) > K p an
1 j Nn
:
t ;
i=1
i=1
8
9
s
n
<X
=
C
Nn max Pr& in (xj ) > K p an
1 j Nn:t ;max& in (xj ) > K

i=1

2= O Nn n K =4 :

25

K2
4

:

Moreover, (4) implies that(1b 3=2 = O n 1

b2
ln n

1

1

)

O n1

1

;

where the last inequality holds because b =2 = ln n is convergent. Then, picking K =p2
2 2 (1 + ) + 1= (1) yields Nn n K =4 = O n (1+ ) so that
8
9
s
n
1
1
= X
<
X
X
C
1
& in (xj ) > K p an
Pr
max
O
< 1:
1+
:1 j Nnt ; n=1nn=1
i=1Therefore, by the Borel-Cantelli lemma,max

1 j Nnn
X(A9)& in (xj ) = O (an ) a:s:

i=1

It follows from (A6) and (A9) thatno
sup f^ (x) E f^ (x)
x2I0

max

1 j Nnn
X& in (xj ) + max sup
1 j Nn x2Ij

i=1

= O (an ) a:s:n
X

& in (x)i=1nX& in (xj )i=1This completes the proof.

A.4

Proof of Proposition 2

Proof of (i). Dividing both sides of (A2) bya (a ), we haveP a + 1; z0 = P a ; z0

(a + 1) and using(a + 1) =R a ; z0 :Then, J (x) can be rewritten as
J (x) = P a ; z0p

P a+ ; z0

R a+ ; z0

R a ; z0

We start from approximating P (a ; z0 ).

Let y

26

(A10)

solve the equation z0 = a +a y . Then,
z0 a
t0 (x
y = p
= p
a
b (x

:

)
)

:

It follows from equation (1) of Pagurova (1965) that
P a ; z0 = P a ; a +
where

p

a y

=

1
p
3 a

y

(3)

y

(A11)

+ Ra ;( ) is the cdf of N (0; 1), and the remainder term Ra is (1=a ) times somelinear combination of higher-order derivatives oforder Taylor expansion of(y ) around( ) evaluated at y . Now, a third-= 0 yields

y
t0 x
p
bx
1
t0 + x+2
4x3
3t0 + x+4x5=2

t0 + x
t0 x
p
3=2
2x
b1=2
bx
2
t0 x
(1)
p
b
bx
2
3
t0 x
p
+
O
b1=2
b3=2
bx

=

where the O3=b3=2 rate is uniform on I0 . It follows that

y

y

+

=

1
b1=2
p = p
x
a

of

(3)

2

y

where the O (

=
2

(A13)

+ O 3=2 :
b1=2
b
p
(3)
(y ) = a . Notice that

1

2x

+O

2

(A14)

;

) rate is uniform on I0 . Moreover, a second-order Taylor expansion(y ) around
(3)

3

t0 x
p
bx

t0 + xx3=2

Next we work on the term involving

where the O (

(A12)

;

= 0 yields(3)

t0 + x
2x3=2

t0 x
p
bx

2

t0 x
p
bx

(4)

b1=2

+O

b

;

=b) rate is again uniform on I0 . By (3), we have
1
p
a

(3)

1
p
a+

y

(3)

y+
3

O b1=2 O= O( ) = o

b1=2

(A15)

b3=2

uniformly on I0 .
Furthermore, using the above argument leads to
3

jRa

Ra+ j

O (b) O

b1=2
27

=O

b1=2 = o

b3=2

(A16)

again uniformly on I0 . Combining (A11), (A13), (A15) and (A16) gives+

P a ; z0
where the O

P a ; z0 =3

t0 + x
x3=2

3

t0 x
p
bx

+O

b1=2

b3=2

;

(A17)

=b1=2 rate is uniform on I0 .
R (a+ ; z0 ). Using

We should also evaluate the order of magnitude in R (a ; z0 )(A1) yields
z0a exp ( z0 )(a + 1)
1 + O (1=a )
1
z0 a
p
p
=
exp a
z0
a
2
a
p
where f1 + O (1=a )g = 2 = O (1) uniformly on I0 . In addition,
R a ; z0 =

for
1

z0 a
exp a
z0 =: exp aa:= z0 =a > 0. It follows from ln +1

O (1). Then, by (A14), we have

ln

(A18)

;

(A19)

+1

0 that exp fa (ln

+1

)g

3

R a ; z0

R a+ ; z0b1=2 = o

O

(A20)

b3=2

uniformly on I0 . Part (i) is established by substituting(A17)
o and (A20) into (A10)npand then denoting Q (x) := (t0 + x) =x3=2(t0 x) = bx .
Proof of (ii). Notice that
(1)

Q

q (x)
(x) = pexp
2 2 bx7=2

(

(x

t0 )2
2bx

)

(A21)

;

where q (x) := bx (x + 3t0 ) (x t0 ) (x + t0 )2 . Heuristically, q (x)(x t0 ) (x + t0 )2
for b 0, and we can see that Q (x) is maximized at x t0 . To be more precise,the maximizer of Q (x) is given by the solution of q (x) = 0 or(x

bx (x + 3t0 ) =2t0 ) (x + t0 ) . It is not hard to see that for a su¢ ciently small b > 0, the latterhas a unique solution on R+ .
Let t be the solution. Our argument so far suggests that t

t0 for a su¢ ciently

small b > 0. Now we even conjecture that t = t0 + cb for some jcj < 1. To verify
this conjecture, consider that 0 = q (t ) =b (t0 + cb) (4t0 + cb) + c (2t0 + cb)2 .

Because b > 0, c solves (t0 + cb) (4t0 + cb) + c (2t0 + cb)2 = 0. This equation can befurther rewritten as b (c + 1) (c + 2t0 =b)2 = t0 c. Observe that the left- and righthand sides are cubic and downward sloping linear functions of c, respectively. For asu¢ ciently small b > 0,2t0 =b <1 holds, and in this case we can recognize thatthe equation has a unique solution c 2 ( 1; 0). This completes the proof.
28

A.5

Proof of Theorem 1

Because t^that jt

t^

t0

t0 j

t + jtt0 j and we have already known in Proposition 2(ii)b, it su¢ ces to show that
t^

t

(A22)

= O (b) a:s:

The proof for (A22) closely follows the one for Theorem 1 of Chu and Cheng (1996).
We keep using the same notation as in the proof of (6) in Proposition 1.
In
addition, de…ne:= fx : x 2 I0 ; jx

t j > b jI0 jg.

Also let w be an element of

En := fx0 ; x1 ; : : : ; xNn g that is closest to t , i.e.,jw

t j = min jxj

t j , w = arg min js

0 j Nn

t j:

s2En

^
= 0. However, the set
n To establish (A22), we
o need to show that Pr t 2
supx2 J^ (x)J^ (w) is larger than the set t^ 2
, and thus Pr t^ 2onJ^ (w) is the case. In what follows, we demonstrate that
Pr supx2 J^ (x)Pr sup J^ (x)J^ (w)

(A23)

= 0:

x2

Let x := arg supx2no
E J^ (w)no
E J^ (x) . Then, by (5) and Proposition 2,no
sup E J^ (x) = jd0 j fQ (w)Q (x )g

x2

+O

b1=2

3

b3=2

:

For a su¢ ciently large n, w is closer to t than x . Then, Q (w) Q (x ) > 0
is the case. In addition, Q (w) Q (x ) = Q(1) (w) (w x ) for some w between
w and x by the mean value theorem.
(1)

Q (x ) = Q(w) jw

x j.

Combining these, we may write Q (w)Now, by (A21), Q(1) (x) = O (b 1 ) uniformly on

x 2 I0 . Because Q(1) (w) > 0, there is some constant C > 0 so that b Q(1) (w)or Q(1) (w)
C=b holds. Moreover,
jw

xj

jx

tj

jw

tj

b jI0 jNn 1 jI0 j = b 1

Since n (1+ ) b1=2 is convergent, we may put n (1+ ) b1=2
n so that jw(1=2) b jI0 j. Therefore,nono
E J^ (w)
sup E J^ (x)n (1+ ) b1=2 jI0 j :1=2 for a su¢ ciently largexj

3

3L + Ox2

where

L = Ln :=1
jd0 j C jI0 j
6
29

b1=2

C,

:

b3=2

;

By this result and a straightforward calculation,
n
o
^
^
^
^
sup J (x)J (w)2 sup J (x) E J (x)x2

33L + O

b3=2

x2I0

:

However, combining the de…nition of L with (3) and (6) implies thatno3
2 sup J^ (x) E J^ (x)L <0
L + O 3=2
b
x2I0
o
n
with probability one. Therefore, Pr supx2 J^ (x) < J^ (w) = 1, and (A23) isestablished. This completes the proof.

A.6

Proof of Theorem 2

The proof requires the following lemmata.
Lemma A4. For a; z;(a + m + 1; z)=: pm P (a; z)(a + 1)mwhere pm+1 => 0 and m 2 f0; 1; 2; : : :g,(a + m + 1) pm , rm+1 =(a + m + 1) rm +( z)m+1 , and p0 = r0 = 1.

Lemma A5. As n ! 1,
P z ; z0 =12p

R z ; z0 = p

2

1
p

2

t0

b1=2
1
p
2 t0

b1=2
2t0

+O
1
2t0

b1=2

J^(1) (t0 )

2

b

!

Lemma A7. As n ! 1,)(r
nb5=2 ^(1)
3
V ar
J (t0 ) ! p 5=22
2 t0
Lemma A8. If t^

for any

; and

b1=2

Lemma A6. As n ! 1,
E

rm R (a; z) ;3

+ O max b;

r

2 d0

30

:

f (t0 ) + f (t+0)
2

t0 = op b1=2 , then, as n ! 1,rb3=2 ^(2)2 d0
p
J ( )!
3=2
t0

between t^ and t0 .

3=2

t0

b

:

:

Lemma A9. As n ! 1, E jHi j3 = O (
A.6.1

3

=b4 ).

Proof of Theorem 2

By a mean-value expansion of the …rst-order condition J^(1) t^ = 0, we have
0 = J^(1) (t0 ) + J^(2) t t^ t0no h= E J^(1) (t0 ) + J^(1) (t0 )noi
E J^(1) (t0 ) + J^(2) t

t^

(A24)

t0

for some t between t^ and t0 . Theorem 1 indicates that the condition for Lemma A8
is satis…ed.
Rearranging (A24), we obtain
893^(1) (t0 ) =
<
J
E
n 4^
5
t t0
1=2
:
;
b
J^(2) tno32r^(1) (t0 ) E J^(1) (t0 )
J
n 4
5:
=
b1=2
J^(2) t2rThen, by Lemmata A6 and A8, the leading bias term for t^ becomesno
E J^(1) (t0 )= b + op (b) :J^(2) thpJ^(1) (t0 )To demonstrate asymptotic normality of n=b1=2 b3=2 =we also check Lyapunov’s condition. This quantity can be expressed as
r
n
n
X
X
b5=2
Yi :=
fHi E (Hi )g :
n 2
i=1
i=1

(A25)

n
oi
E J^(1) (t0 ) ,Then, by Cr -inequality, Jensen’s inequality (due to the convexity of y 3 for y
Lemma A9,3

E jYi j

8

b5=2
n 2

3=2

E jHi j3 = O n 3=2 b 1=4 :It also follows from Lemma A7 that V ar (Yi ) = O (n 1 ). Therefore,Pn3
1
i=1 E jYi j
=O p
! 0;
Pn
3=2
nb1=2
f i=1 V ar (Yi )gand Lyapunov’s condition is established.

310) andNow we are allowed to employ a central limit theorem, in conjunction with Lemma
A7, to obtain
r

n

b3=2

b1=2

h

J^(1) (t0 )noidE J^(1) (t0 ) ! N

30; p 5=2
2 t0

Then, it holds that for the right-hand side of (A25),
n
o3
2
r
3=2
^(1) (t0 ) E J^(1) (t0 )
b
=
J
n 4
d
5!
N
1=2
b
(b3=2 = ) J^(2) tf (t0 ) + f (t+
0)
2

p 1=2
3 t0
0;
4d20

!

:

f (t0 ) + f (t+0)
2!

by Lemma A8. This completes the proof.

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

Declaration of Interest
The authors report that there are no competing interests to declare.

Declaration of Generative AI in Scienti…c Writing
The authors do not use generative arti…cial intelligence (AI) or AI-assisted technologies in the writing process.

Data Statement
The datasets used in Section 5 are openly available; see Section 5.1 for more details.

32References[1] Bader, B., J. Yan, and X. Zhang (2018): “Automated Threshold Selection for
Extreme Value Analysis via Ordered Goodness-of-Fit Tests with Adjustment for
False Discovery Rate,”Annals of Applied Statistics, 12, 310 - 329.
[2] Chen, S. X. (2000): “Probability Density Function Estimation Using Gamma
Kernels,”Annals of the Institute of Statistical Mathematics, 52, 471 - 480.
[3] Chu, C. K., and P. E. Cheng (1996): “Estimation of Jump Points and Jump Values of a Density Function,”Statistica Sinica, 6, 79 - 95.
[4] Clauset, A., C. R. Shalizi, and M. E. J. Newman (2009): “Power-Law Distributions in Empirical Data,”SIAM Review, 51, 661 - 703.
[5] Cooray, K., and M. M. A. Ananda (2005): “Modeling Actuarial Data with a Composite Lognormal-Pareto Model,”Scandinavian Actuarial Journal, 5, 321 - 334.
[6] Couallier, V. (1999): “Estimation Non Paramétrique d’une Discontinuité dansune Densité,” Comptes Rendus de l’Académie des Sciences - Série I, 329, 633 636.
[7] Cowell, F. A., F. H. G. Ferreira, and J. A. Litch…eld (1998): “Income Distributionin Brazil 1981-1990: Parametric and Non-Parametric Approaches,” Journal of
Income Distribution, 8, 63 - 76.
[8] Cunnane, C. (1973): “A Particular Comparison of Annual Maxima and Partial
Duration Series Methods of Flood Frequency Prediction,”Journal of Hydrology,
18, 257 - 271.
[9] Danielsson, J., L. M. Ergun, L. de Haan, and C. G. de Vries (2019): “Tail Index
Estimation: Quantile-Driven Threshold Selection,”Bank of Canada Sta¤ Working Paper 2019-28.
[10] Davison, A. C., and R. L. Smith (1990): “Models for Exceedances over High
Thresholds,”Journal of the Royal Statistical Society, Series B, 52, 393 - 442.
[11] Delaigle, A., and I. Gijbels (2006): “Estimation of Boundary and Discontinuity
Points in Deconvolution Problems,”Statistica Sinica, 16, 773 - 788.
[12] Desmet, L., I. Gijbels, and A. Lambert (2010): “Estimation of Irregular Probability Densities,” in J. Antoch, M. Hušková and P. K. Sen (eds.), Nonparametrics and Robustness in Modern Statistical Inference and Time Series Analysis: A
Festschrift in Honor of Professor Jana Jureµcková. Beachwood, OH: Institute ofMathematical Statistics, 46 - 61.
[13] Drees, H., A. Janßen, S. I. Resnick, and T. Wang (2020): “On a Minimum Distance Procedure for Threshold Selection in Tail Analysis,” SIAM Journal on
Mathematics of Data Science, 2, 75 - 102.
[14] DuMouchel, W. H. (1983): “Estimating the Stable Index in Order to Measure
Tail Thickness: A Critique,”Annals of Statistics, 11, 1019 - 1031.
[15] Funke, B., and M. Hirukawa (2019): “Nonparametric Estimation and Testing on
Discontinuity of Positive Supported Densities: A Kernel Truncation Approach,”Econometrics and Statistics, 9, 156 - 170.
33

[16] Funke, B., and M. Hirukawa (2023): “Density Derivative Estimation Using Asymmetric Kernels,”Journal of Nonparametric Statistics, forthcoming.
[17] Hirukawa, M., Murtazashvili I., and A. Prokhorov (2022): “Uniform Convergence
Rates for Nonparametric Estimators Smoothed by the Beta Kernel,” Scandinavian Journal of Statistics, 49, 1353 - 1382.
[18] Huh, J. (2002): “Nonparametric Discontinuity Point Estimation in Density or
Density Derivatives,”Journal of the Korean Statistical Society, 31, 261 - 276.
[19] Huh, J. (2012): “Bandwidth Selections Based on Cross-Validation for Estimationof a Discontinuity Point in Density”[in Korean], Journal of the Korean Data &Information Science Society, 23, 765 - 775.
[20] James, R., H. Leung, and A. Prokhorov (2023): “A Machine Learning Attack on
Illegal Trading,”Journal of Banking and Finance, 148, Article No.106735.
[21] Jenkins, S. P. (2017): “Pareto Models, Top Incomes and Recent Trends in UK
Income Inequality,”Economica, 84, 261 - 289.
[22] Joo, J.-H., and P. Qiu (2009): “Jump Detection in a Regression Curve and Its
Derivative,”Technometrics, 51, 289 - 305.
[23] Klugman, S. A., H. H. Panjer, and G. E. Willmot (2019): Loss Models : From
Data to Decisions, 5th Edition. Hoboken, NJ: John Wiley & Sons.
[24] Loretan, M., and P. C. B. Phillips (1994): “Testing the Covariance Stationarity
of Heavy-Tailed Time Series: An Overview of the Theory with Applications to
Several Financial Datasets,”Journal of Empirical Finance, 1, 211 - 248.
[25] Markovitch, N. M., and U. R. Krieger (2000): “Nonparametric Estimation ofLong-Tailed Density Functions and Its Application to the Analysis of World
Wide Web Tra¢ c,”Performance Evaluation, 42, 205 - 222.
[26] Marron, J. S. (1985): “An Asymptotically E¢ cient Solution to the Bandwidth
Problem of Kernel Density Estimation,”Annals of Statistics, 13, 1011 - 1023.
[27] McNeal, A. J. (1997): “Estimating the Tails of Loss Severity Distributions Using
Extreme Value Theory,”ASTIN Bulletin, 27, 117 - 137.
[28] Newey, W. K., and D. L. McFadden (1994): “Large Sample Estimation and Hypothesis Testing,” in R. F. Engle and D. L. McFadden (eds.), Handbook of Econometrics, Volume IV. Amsterdam: Elsevier, 2111 - 2245.
[29] Pagurova, V. I. (1965): “An Asymptotic Formula for the Incomplete Gamma
Function,” USSR Computational Mathematics and Mathematical Physics, 5,
162 - 166.
[30] Resnick, S. I. (1997): “Discussion of the Danish Data on Large Fire Insurance
Losses,”ASTIN Bulletin, 27, 139 - 151.
[31] Reynkens, T., R. Verbelen, J. Beirlant, and K. Antonio (2017): “Modelling Censored Losses Using Splicing: A Global Fit Strategy with Mixed Erlang and Extreme Value Distributions,”Insurance: Mathematics and Economics, 77, 65 - 77.
[32] Sasaki, Y., and Y. Wang (2023): “Extreme Changes in Changes,” Journal ofBusiness & Economic Statistics, forthcoming.
34[33] Scarrott, C., and A. MacDonald (2012): “A Review of Extreme Value Threshold
Estimation and Uncertainty Quanti…cation,” REVSTAT - Statistical Journal,
10, 33 - 60.
[34] Scollnik, D. P. M., and C. Sun (2012): “Modeling with Weibull-Pareto Models,”
North American Actuarial Journal, 16, 260 - 272.
[35] Van der Vaart, A. W., and J. A. Wellner (1996): Weak Convergence and Empirical
Processes: With Applications to Statistics. New York: Springer-Verlag.
[36] Van Es, B. (1991): “Likelihood Cross-Validation Bandwidth Selection for Nonparametric Kernel Density Estimators,”Journal of Nonparametric Statistics, 1,
83 - 110.
[37] Wu, J. S., and C. K. Chu (1993a): “Kernel-Type Estimators of Jump Points andJump Values of a Regression Function,”Annals of Statistics, 21, 1545 - 1566.
[38] Wu, J. S., and C. K. Chu (1993b): “Nonparametric Function Estimation and
Bandwidth Selection for Discontinuous Regression Functions,”Statistica Sinica,
3, 557 - 576.
Figure 1: The Diagnostic Function and Its Approximations

Table 1: Characteristic Numbers of Underlying Distributions
Model
1 A
B
2 A
B

Distribution
Log-Normal + Quadratic (D = 1=4)
Log-Normal + Quadratic (D = 3=22)
Splicing with Weibull & GPD
Splicing with Weibull & Translated Weibull

35

Mode
0:69590:69592:4023
2:4023cL
0:96590:9583
0:95390:9539f t0
0:17280:12790:1063
0:1063f t+
0
0:02280:02790:0115
0:0115d0
0:15000:10000:09480:0948

Table 2: Monte Carlo Resultsn = 250
RM SE

E stim a to r

B ia sSDCC

0:22600:63400:6731KS
Q -M A D
Q -S U P
AEB

0:87700:29900:10390:6615

0:5392
0:15690:40140:8305

1:0295
0:3376
0:4146
1:0618^ or ^hb# fI0 g

B ia sSDn = 500
RM SE^ or ^hb# fI0 gM o d e l 1 -A : L o g -N o rm a l + Q u a d ra tic (D = 1=4)

S G -L S

0:55
0:60
0:65
0:70

0:66800:65780:6501
0:6333

0:4375
0:44230:4444
0:4497

0:7985
0:7927
0:7874
0:7767

S G -L S -B C

0:55
0:60
0:65
0:70

0:33780:33150:32800:31920:2250
0:22730:22820:2313

0:40580:40200:3996
0:3942

S G -L

0:55
0:60
0:65
0:70

0:56280:51180:47580:45900:44620:4494
0:44680:4433

0:7182
0:6811
0:6527
0:6381

S G -L -B C

0:55
0:60
0:65
0:70

0:2904
0:26680:24970:24300:24740:2546
0:2575
0:2598

0:3814
0:36880:3586
0:3557

S G -M L

0:55
0:60
0:65
0:70

0:1375
0:04650:00350:0352

0:4524
0:3836
0:3204
0:2737

0:47280:38640:32040:2759S G -M L -B C

0:55
0:60
0:65
0:70

0:16900:08000:03160:0011

0:43780:3714
0:30980:26390:4693
0:37990:3114
0:2639CC

0:32000:6223

0:6998KS
Q -M A D
Q -S U P
AEB

0:9832
0:3227
0:07200:94180:6012
0:26880:63800:8936

1:15250:42000:6421
1:2983

0:0116(0:0092)((((

0:3302
0:3263
0:3221
0:3141))))(0:2076)(0:2052)(0:2011)(0:1966)((((

0:0315
0:03350:0351
0:0364

644
1000
976
729

(0:2182)(0:2206)(0:2220)(0:2243)((((

0:2725
0:24490:2262
0:2160

)
)
)
)

)
)
)
)

(0:0210)(0:0196)(0:0188)(0:0187)
(
(
(
(

)
)
)
)

0:0487

0:5775

0:5796

0:5234
0:32290:08000:7181

0:7046
0:1962
0:67080:57450:87780:37780:67550:9197

0:7461
0:7351
0:71850:7053

0:40790:41490:4243
0:43100:8503
0:8441
0:83450:8266

0:37450:36850:36000:35350:20780:21070:2155
0:21920:42830:42450:4196
0:41600:6063
0:5363
0:48980:46080:42690:43200:42880:42190:74150:6886
0:65100:6247

0:31490:2814
0:25870:24600:23570:24170:24370:2437

0:3933
0:37100:3554
0:3463

0:32650:1547
0:06080:0101

0:5146
0:4346
0:3497
0:2831

0:6094
0:4613
0:35500:2832

0:3514
0:1845
0:0933
0:04430:4993
0:4216
0:33900:2741

0:61050:4602
0:3516
0:2777

0:0073(0:0037)((((

0:3716
0:3666
0:3585
0:3518))))(0:1987)(0:1989)(0:1939)(0:1877)((((

0:02490:0297
0:0326
0:0342

707
997
913
708

(0:2027)(0:2066)(0:2113)(0:2147)((((

0:29140:25480:2311
0:2148

)
)
)
)

)
)
)
)

(0:0183)(0:0173)(0:0166)(0:0157)
(
(
(
(

)
)
)
)

M o d e l 1 -B : L o g -N o rm a l + Q u a d ra tic (D = 3=22)

S G -L S

0:55
0:60
0:65
0:70

0:9517
0:94720:93630:93030:20140:21090:2304
0:23980:97270:9703
0:9642
0:9607

S G -L S -B C

0:55
0:60
0:65
0:70

0:4777
0:47550:47050:46790:0973
0:10290:1122
0:11650:48750:4865
0:4837
0:4822

S G -L

0:55
0:60
0:65
0:70

0:9171
0:89580:87890:8711

0:2572
0:2857
0:3064
0:31380:95250:94030:93080:9259S G -L -B C

0:55
0:60
0:65
0:70

0:4617
0:4521
0:44430:44070:12940:14390:15730:1616

0:4794
0:4745
0:4713
0:4694

S G -M L

0:55
0:60
0:65
0:70

0:0586
0:0557
0:1146
0:14980:52520:4497
0:3854
0:3557

0:5285
0:4532
0:4021
0:3859S G -M L -B C

0:55
0:60
0:65
0:70

0:0965
0:0133
0:0696
0:10300:5064
0:43320:3706
0:34080:51550:4334
0:3771
0:35600:0139(0:0126)((((

0:47400:4717
0:46580:4625(0:1070)(0:1115)(0:1220)(0:1273)((((

0:4554
0:4437
0:43470:4304))))(0:1321)(0:1469)(0:1556)(0:1592)((((

0:03790:0423
0:04500:0468

)
)
)
)

)
)
)
)

(0:0278)(0:0279)(0:0275)(0:0275)
(
(
(
(

36

)
)
)
)

594
993
925
602

0:14090:59420:6107

0:6664
0:36380:15491:01000:6955
0:2716
0:96580:6164

0:9632
0:45400:9782
1:1832

0:9940
0:9923
0:9914
0:9856

0:0711
0:08080:0856
0:1130

0:99650:9956
0:9951
0:9921

0:4973
0:49650:4961
0:4931

0:0321
0:0367
0:03900:05450:4983
0:49790:4976
0:4961

0:98090:9652
0:9605
0:95100:12950:17370:18400:2031

0:9894
0:98070:97800:9724

0:4912
0:4833
0:4816
0:47700:06490:0878
0:09250:10290:49550:49120:4904
0:48800:2704
0:0642
0:0434
0:0943

0:5556
0:47290:3983
0:3433

0:61790:47720:4007
0:35600:3012
0:10400:00200:0455

0:53600:45490:3822
0:32800:61490:4667
0:38220:3311

0:0086(0:0056)((((

0:49670:49580:4953
0:4925(0:0392)(0:0444)(0:0468)(0:0590)((((

0:4897
0:48190:47890:4740))))(0:0660)(0:0875)(0:0939)(0:1029)((((

0:03080:03980:0454
0:0489

)
)
)
)

)
)
)
)

(0:0246)(0:0253)(0:0259)(0:0260)
(
(
(
(

)
)
)
)

660
995
837
535

Table 2: (continued)

E stim a to r

B ia sSDn = 250
RM SE

CC

0:3261

0:66790:7432KS
Q -M A D
Q -S U P
AEB

1:3273
0:2051
0:3085
4:3704

0:9856
0:7303
1:62982:9114

1:65330:7585
1:65885:2514^ or ^hb# fI0 g

B ia sSDn = 500
RM SE^ or ^hb# fI0 gM o d e l 2 -A : S p lic in g w ith W e ib u ll & G P D

S G -L S

0:55
0:60
0:65
0:70

0:62380:62680:62900:63080:1250
0:12600:12720:1282

0:6362
0:6393
0:64170:6437

S G -L S -B C

0:55
0:60
0:65
0:70

0:4736
0:47170:47080:47000:08480:0841
0:08400:08390:4811
0:4792
0:4782
0:4775

S G -L

0:55
0:60
0:65
0:70

0:6196
0:6216
0:6231
0:62390:1186
0:12330:1282
0:1346

0:63090:6337
0:6362
0:6383

S G -L -B C

0:55
0:60
0:65
0:70

0:43480:44490:4514
0:4547

0:0715
0:07590:0811
0:0881

0:4407
0:4513
0:4586
0:4632

S G -M L

0:55
0:60
0:65
0:70

0:4484
0:45180:4506
0:4334

0:3354
0:3036
0:29680:3124

0:56000:5443
0:5395
0:5342

S G -M L -B C

0:55
0:60
0:65
0:70

0:4152
0:4186
0:4171
0:39940:33130:3005
0:29400:3092

0:5311
0:51530:5102
0:5051

0:0104(0:0228)((((

0:1503
0:1550
0:15830:1608

)
)
)
)

91
975
864
81

(0:0644)(0:0675)(0:0705)(0:0728)((((

0:18480:17670:17180:1692))))(0:0668)(0:0660)(0:0648)(0:0638)((((

0:0333
0:0332
0:03350:0340))))(0:0082)(0:0074)(0:0073)(0:0075)
(
(
(
(

)
)
)
)

0:19990:6443

0:6746

0:65560:1488
1:0775
4:2648

1:8320
1:48842:9302
2:0070

1:9458
1:49583:1221
4:71350:6202
0:6237
0:62650:6285

0:0914
0:0914
0:09190:09350:62690:6304
0:6333
0:6354

0:4674
0:46590:4652
0:46490:0615
0:0601
0:0592
0:0592

0:4714
0:46980:46900:4687

0:6224
0:6242
0:62590:6275

0:0831
0:0871
0:0898
0:09150:62790:6303
0:6323
0:6342

0:4362
0:4464
0:4531
0:45740:0491
0:0522
0:0547
0:0566

0:43890:4494
0:45640:46090:42590:45680:4697
0:4733

0:3665
0:29980:26320:25080:56190:5464
0:5384
0:5356

0:39280:4237
0:43600:43900:3612
0:29640:26090:2487

0:5336
0:5171
0:5081
0:5046

0:0060(0:0023)((((

0:15280:15780:16130:1635))))(0:0470)(0:0466)(0:0459)(0:0452)((((

0:0331
0:0332
0:0337
0:0343

162
905
784
29

(0:0450)(0:0473)(0:0495)(0:0514)((((0:18620:17790:17280:1701

)
)
)
)

)
)
)
)

(0:0079)(0:0064)(0:0057)(0:0057)
(
(
(
(

)
)
)
)

M o d e l 2 -B : S p lic in g w ith W e ib u ll & Tra n sla te d W e ib u llCC

0:3302

0:67880:7548KS
Q -M A D
Q -S U P
AEB

1:47140:2656
0:1266
6:49590:9325
1:3378
1:8973
13:3971

1:7420
1:3639
1:9015
14:8889

S G -L S

0:55
0:60
0:65
0:70

0:6123
0:6154
0:61750:6193

0:12780:12830:12960:1307

0:62550:6287
0:63090:6329S G -L S -B C

0:55
0:60
0:65
0:70

0:47390:47300:4721
0:47200:08990:08880:0883
0:0881

0:4824
0:48120:48030:4801

S G -L

0:55
0:60
0:65
0:70

0:6036
0:60280:6022
0:6023

0:13280:13970:1461
0:15180:6181
0:61880:6196
0:6211

S G -L -B C

0:55
0:60
0:65
0:70

0:45340:46000:46350:4654

0:0857
0:0925
0:0993
0:1061

0:46150:4692
0:4741
0:4774

S G -M L

0:55
0:60
0:65
0:70

0:4056
0:40990:3956
0:3742

0:3411
0:3227
0:3177
0:32800:53000:5217
0:5074
0:4976

S G -M L -B C

0:55
0:60
0:65
0:70

0:3732
0:3774
0:3628
0:3412

0:33700:3194
0:3147
0:32490:50280:49440:4802
0:4711

0:0116(0:0351)((((

0:13840:14250:14530:1473))))(0:0591)(0:0619)(0:0645)(0:0665)((((

0:1502
0:14280:13860:1368))))(0:0621)(0:0614)(0:0605)(0:0596)((((

0:0325
0:0326
0:03290:0330))))(0:0087)(0:0082)(0:0078)(0:0077)
(
(
(
(

)
)
)
)

124
972
921
398

0:16090:6716

0:6906

1:0546
0:4715
2:7196
3:7522

3:0694
4:2932
16:7906
5:4395

3:2455
4:3190
17:0094
6:6081

0:5991
0:6024
0:60500:60690:09680:0974
0:0982
0:0993

0:60680:6103
0:61300:61500:4626
0:46200:4616
0:4616

0:06700:06580:0652
0:06490:46750:4667
0:4662
0:4662

0:5977
0:59690:59690:5976

0:09680:1013
0:1041
0:1056

0:6055
0:6054
0:60590:60680:45000:45650:46050:46300:0611
0:0656
0:0687
0:0706

0:4541
0:4612
0:4656
0:4683

0:37690:40480:4224
0:41180:36880:3125
0:25950:25590:5273
0:5114
0:49580:4848

0:34180:3693
0:38650:37580:36250:30790:25670:2532

0:49820:48080:46400:45320:0058((((0:13640:14050:14340:14530:14770:14040:1364
0:1346))))(0:0443)(0:0438)(0:0430)(0:0423)((((

0:03500:03550:03590:0359))))(0:0427)(0:0448)(0:0467)(0:0484)((((

Note: “Bias”, “SD” and “RMSE” are biases, standard deviations and root meansquared errors of jump location or threshold estimates over 1000 Monte Carlo samples,
respectively. For kernel-smoothed jump location estimation procedures, Monte Carloaverages and standard deviations (in parentheses) of CV tuning parameters are presen^ or ^b”. For automated threshold detection methods, counts
ted under the heading “hof threshold estimates falling into I0 are reported under the heading “# fI0 g”.
37(0:0020)))))(0:0101)(0:0090)(0:0080)(0:0078)
(
(
(
(

)
)
)
)

173
936
801
376

Table 3: Descriptive Statistics of Datasets
Data

n

M ean

SD

SK

M in.

Q1

Q2

Q3

90%

95%

99%

M ax.

(A) D anish Fire Insurance Losses (in m illions of D anish kroner)
Original

2; 492

3:063

7:975

19:884

0:313

1:157

1:634

2:646

5:080

8:406

24:614

263:250

9:980
9:954

11:000
11:000

14:793
14:64099:99040:000(B) U S M ale H ourly Wages (in U S dollars)
Original
Downsized

54; 769
2; 709

6:136
6:080

2:972
2:858

3:013
1:525

0:500
1:000

3:850
3:750

5:630
5:550

8:000
8:000

Note: n = sample size; Mean = average; SD = standard deviation; SK = skewness;
Min. = minimum value; Q1 = …rst quartile; Q2 = median (i.e., second quartile); Q3= third quartile; 90% = 90% quantile; 95% = 95% quantile; 99% = 99% quantile;and Max. = maximum value.
Table 4: Results of Threshold Estimation
Data

Estimator

Estimate of t0

Estimator

I0

Estimate of t0^ or ^bh(A) Danish Fire Insurance Losses (in millions of Danish kroner)
Original

KS
Q-MAD
Q-SUP
AEB

1:375
29:037
11:123
25:288CC
SG-ML
SG-ML-BC

0:70[1; 30]

30:0000:005[1; 30]

1:861
2:096

0:235

[5; 15]

14:9980:005[5; 15]

9:524
9:664

0:140

(B) US Male Hourly Wages (in US dollars)
Downsized

KS
Q-MAD
Q-SUP
AEB

10:400
10:000
10:000
13:500CC
SG-ML
SG-ML-BC38

0:70