Discrete Extremes

arXiv:1707.05033v1 [math.ST] 17 Jul 2017

Adrien Hitz, Richard Davis and Gennady Samorodnitsky

Our contribution is to widen the scope of extreme value analysis applied to discretevalued data. Extreme values of a random variable X are commonly modeled using
the generalized Pareto distribution, a method that often gives good results in practice.
When X is discrete, we propose two other methods using a discrete generalized Pareto
and a generalized Zipf distribution respectively. Both are theoretically motivated and
we show that they perform well in estimating rare events in several simulated and real
data cases such as word frequency, tornado outbreaks and multiple births.

1

Introduction

Extreme quantile estimation is an important but difficult problem in statistics, especially
when the quantile is beyond the range of the data. In the univariate case, extreme value
theory motivates the choice of a parametric family called the generalized Pareto distribution
(GPD) which is used to model the tail and estimate the probability of rare events (Pickands
[1975]). Let X be a random variable taking values in [0, xF ) with survival function F̄X ,
where xF ∈ R+ ∪ {∞} and R+ = (0, ∞). Suppose that there exists a strictly positive
sequence au such that
d

a−1
u (X − u) | X ≥ u → Z,

(1)

as u → xF , for some Z following a non-degenerate probability distribution on [0, ∞), where
d
→ denotes weak convergence. A stunning result is that this assumption is sufficient to
characterize the limiting distribution: Z follows a generalized Pareto distribution (GPD),
defined by its survival function

x −1/ξ
1{x<τ } ,
1 − FGPD (x; σ, ξ) = F̄GPD (x; σ, ξ) = 1 + ξ
σ

x ≥ 0,

with τ = ∞ if ξ ≥ 0, τ = σ/|ξ| if ξ < 0, where 1{x<τ } is 1 if x < τ and 0 otherwise, and
σ > 0. We use the convention that if ξ = 0, then (1 + ξx)1/ξ = ex . Condition (1), written
as X ∈ MDAξ , means that X is in the maximum domain of attraction of an extreme value
distribution with shape parameter ξ (see e.g. Resnick [1987]). Sometimes one says that the

Adrien Hitz, University of Oxford, 24-29 St Giles, Oxford OX1 3LB, UK.
Gennady Samorodnitsky, Cornell University, 220 Rhodes Hall, Ithaca, NY.
Richard Davis, Columbia University, 1255 Amsterdam Avenue, New York, NY.

law F of X is in MDAξ . In this case, the sequence of cumulative distribution functions of
a−1
u (X − u) | X ≥ u converges uniformly to FGPD on [0, ∞). Thus, for large u,
−1
P (X − u > x | X ≥ u) = P {a−1
u (X − u) > au x | X ≥ u}

is often approximated by
≈ F̄GPD (x; σau , ξ) ,

(2)

motivating the approximation of the distribution of exceedances of X above a large threshold u by a GPD (Davison and Smith [1990]). In the continuous case, where X is assumed
to have a continuous distribution, most common distributions belong to some maximum
domain of attraction, and the GPD approximation to the tail of the distribution can be
applied. This approximation often works well in practice, although it can be poor when u
is not large enough, for instance if X is normal and u is around the 90th percentile of its
distribution.
Two issues are apparent with applying the GPD method of approximating the distribution tail to a discrete distribution. First, a necessary condition for a discrete distribution FX to be in some maximum domain of attraction in the case xF = ∞ is that FX
is long-tailed, i.e., F̄X (u + 1)/F̄X (u) → 1 as u → ∞, see Shimura [2012] and Anderson
[1970, 1980]. Without being in a maximum domain of attraction the GPD approximation
(2) does not necessarily apply. Note that many common discrete distributions, including
geometric, Poisson and negative binomial distributions, are not long-tailed. The second
issue in approximating a discrete distribution by a GPD, a continuous distribution, is that
ties are not allowed. To overcome these limitations, we suggest two alternative methods
of modeling the tails of discrete observations, each one relying on a specific assumption
on the underlying distribution. For this work, we will only consider approximations to
distributions with infinite support.
We say that a discrete random variable X with non-negative values is in the discrete
maximum domain of attraction, which we write as X ∈ D-MDAξ , if there exists a random
variable Y ∈ MDAξ with ξ ≥ 0 such that Pr(X ≥ k) = Pr(Y ≥ k) for k = 0, 1, 2, . . . , i.e.,
d
the equality in distribution, X = bY c, holds. We call Y an extension of X and such an
extension is not unique. Shimura [2012] showed that X ∈ MDAξ for some ξ ≥ 0 if and only
if X ∈ D-MDAξ and X is long-tailed. It was also shown by Shimura [2012] that geometric,
Poisson and negative binomial distributions are in D-MDA. This set is, therefore, strictly
larger than the set of all discrete distributions in MDAξ .
d
Let X ∈ D-MDAξ and Y ∈ MDAξ be the corresponding extension satisfying X = bY c.
Then, for large integers u,
P (X − u = k | X ≥ u) = P (Y − u ≥ k | Y ≥ u) − P (Y − u ≥ k + 1 | Y ≥ u)
2

which we will approximate as in (2) by
≈ pD-GPD (k; σau , ξ),

(3)

where pD-GPD is the probability mass function of the discrete generalized Pareto distribution (D-GPD) defined by
pD-GPD (k; σ, ξ) = F̄GPD (k; σ, ξ) − F̄GPD (k + 1; σ, ξ),

(4)

for k = 0, 1, 2, . . .. The discrete generalized Pareto distribution has been used by Prieto
et al. [2014] to model road accidents, while various aspects of discrete Pareto-type distributions were studied in Krishna and Pundir [2009], Buddana and Kozubowski [2014], and
Kozubowski et al. [2015]. Notice that the scale parameter in (3) is undetermined since the
extension Y is non-unique.
We now discuss an alternative assumption on the distribution of the discrete random
variable X that will allow us to construct another approximation of its tail. Let pX be
the probability mass function of X and suppose that there exists a non-negative random
variable Y ∈ MDAξ/(1+ξ) with ξ ≥ 0 such that pX (k) = c F̄Y (k) for k = d, d + 1, d + 2, . . . ,
for some c > 0 and d ∈ N0 = {0, 1, . . .}. We denote this condition by pX ∈ D-MDAξ/(1+ξ)
and call F̄Y an extension of pX . Then for large integers u,
pX (u + k)/pX (u)
P (X − u = k | X ≥ u) = P∞
i=0 pX (u + i)/pX (u)
P (Y > u + k)/P (Y > u)
= P∞
i=0 P (Y > u + i)/P (Y > u)
which can be approximated as in (2),

≈ pGZD k; (1 + ξ)σau , ξ ,

(5)

where
1 + ξ σk

pGZD (k; σ, ξ) = P
∞

i=0

−1/ξ−1

1 + ξ σi

−1/ξ−1 ,

k = 0, 1, 2, . . . .

(6)

We call the probability mass function in (6) generalized Zipf distribution (GZD). In the
case ξ > 0, the GZD is a Zipf–Mandelbrot distribution whose probability mass function is
usually written in the form p(k) = (k + q)−s /Hs,q , for k = 0, 1, 2, . . . , s > 1 and q > 0,
where Hs,q is the Hurwitz-Zeta function (Mandelbrot [1953]). The GZD is of this form
with s = 1 + 1/ξ and q = σ/ξ. When q = 1 the distribution is a Zipf law which is
sometimes presented as the counterpart of the Pareto distribution because its probability
mass function, after a shift, can be written in a homogeneous form (Arnold [1983]). Zipftype families of discrete distributions have been fitted to various data sets such as word
3

frequencies (Booth [1967]), city sizes (Gabaix [1999]), company sizes (Axtell [2001]) and
numbers of website hits (Clauset et al. [2009]). In the case ξ = 0, the generalized Zipf
distribution is a geometric distribution with probability of success p = 1 − e−1/σ (and so
is the discrete generalized Pareto distribution).
We will show that the geometric, Poisson and negative binomial probability mass functions belong to D-MDA0 , so the tail approximation (5) makes sense in these cases. In
addition, we will see that pX ∈ D-MDAξ/(1+ξ) for ξ ≥ 0 implies X ∈ MDAξ (under an
additional condition in the case ξ = 0).
Equations (3) and (5) motivate the approximation of the tails of some discrete distributions by a D-GPD and GZD, respectively. If ξ ≥ 0, then
sup
k=0,1,2...

fGPD (k; σ, ξ)
− 1 −→ 0,
σ→∞
q(k; σ, ξ)

for either q = pD-GPD or q = pGZD . One thus expects the GPD, D-GPD and GZD approximations to give similar results when the scale parameter σ is large.
In Section 2, we justify theoretically the three approximations in the case ξ > 0 by
showing that if pX ∈ MDAξ/(1+ξ) , then for any sequence ku ⊂ N0 such that supu ku /u < ∞,
P (X = ku + u | X ≥ u)
→ 1,
q(ku ; ξu, ξ)

k = 0, 1, 2, . . . ,

where q ≡ fGPD , pD-GPD and pGZD (see Theorem 2.3 and Proposition 2.4). A similar
justification is provided in the case ξ = 0 for the D-GPD and GZD approximations (Theorem 2.5). We further present an important invariance property satisfied by the D-GPD
(Proposition 2.7).
In Section 3, we simulate data from different discrete distributions and compare the
ability of the three approximations for estimating the probability of regions far from the
origin. The D-GPD and GZD approximations outperform the GPD when there are many
tied observations, otherwise the results are similar. As opposed to the GZD, the D-GPD
has a closed-form survival and probability mass function allowing for an exact likelihood
based inference.
In Section 4, we study a data set counting the occurrence of British words in a corpus
and show that the three distributions are appropriate to describe the tails of these word
frequencies. We also use the D-GPD and GZD to model the length of French words, the
number of tornado outbreaks in the United States over 50 years, and the number of births
in the United States in the last 20 years, illustrating how they offer potential methods for
estimating the probability of rare events when the data are discrete.

4

2

Theoretical Results

In this section we describe a number of properties of the approximation procedures introduced in Section 1.
Proposition 2.1. If X ∈ D-MDAξ , then there exists a positive sequence (au , u = 1, 2, . . .)
such that
lim

sup

u∈N0 , u→∞ k=0,1,2,...

P (X = u + k | X ≥ u) − pD-GPD (k; au , ξ) = 0 ,

(7)

where pD-GPD is defined in (4).
We remark that if au → ∞, then the result is less interesting because the two terms in
(7) converge to 0.
Proof. By assumption, there exists a random variable Y ∈ MDAξ for ξ ≥ 0 and a positive
d
function (ãu , u > 0) such that X = bY c and the sequence of functions P ã−1
u (Y − u) ≥ x |

Y ≥ u , x ≥ 0, converges uniformly, as u → ∞, to the function F̄GPD (x; σ, ξ), x ≥ 0, for
some σ > 0 and ξ ≥ 0. For a positive integer u we let au = ãu σ. Then
| P (X = u + k | X ≥ u) − pD-GPD (k; au , ξ) |

sup
k=0,1,2,...

=

sup
k=0,1,2,...



−1
−1
−1
P ã−1
u (Y − u) ≥ ãu k | Y ≥ u − P ãu (Y − u) ≥ ãu (k + 1) | Y ≥ u

− F̄GPD (k; au , ξ) + F̄GPD (k + 1; au , ξ)

≤ 2 sup P ã−1
u (Y − u) ≥ x | Y ≥ u − F̄GPD (x; σ, ξ) → 0
x≥0

as u → ∞ over the integers.
The following auxiliary lemma is elementary (as the sum can be sandwiched between
two integrals).
Lemma 2.2. If ξ > 0, then, as u → ∞,
u1/ξ H1+1/ξ,u → ξ,
where Hs,q =

P∞

i=0 (q

+ i)−s is the Hurwitz-Zeta function.

Recall that a positive and measurable function f on [1, ∞) is regularly varying if there
exists a positive function ` such that
lim

u→∞

f (ux)
→ `(x),
f (u)
5

x ≥ 1.

In this case, there exists α ∈ R such that `(x) = xα and we write f ∈ RVα (see e.g.
Bingham et al. [1989]). If f ∈ RV−α for α ≥ 0, then
lim sup

u→∞ x∈[1,b]

f (ux)
− x−α → 0,
f (u)

(8)

for b = ∞ if α > 0, and for any b < ∞ if α = 0. If f ∈ RV−α for α > 0, then by Potter’s
bounds (see e.g. Resnick [1987]) for any  > 0 there is u ∈ (0, ∞) such that
e− x−α− ≤

f (ux)
≤ e x−α+ ,
f (u)

x ≥ 1,

(9)

for u ≥ u . We say that X is regularly varying if F̄X ∈ RV−α for some α > 0, a necessary
and sufficient condition for X ∈ MDA1/α .
The following result sheds some light on the approximation suggested in (5).
Theorem 2.3. If pX ∈ D-MDAξ/(1+ξ) for ξ > 0, then X ∈ MDAξ and for any sequence
of nonnegative integers (ku , u = 1, 2, . . .) such that supu ku /u < ∞,
lim

u∈N0 , u→∞

P (X = ku + u | X ≥ u)
= 1.
pGZD (ku ; ξu, ξ)

(10)

Proof. By assumption, there exists a survival function F̄ such that F̄ (k) = c pX (k) for c >
0, k large enough and F̄ ∈ RV−1/ξ−1 . The last condition is equivalent to F̄ (b·c) ∈ RV−1/ξ−1
(Shimura [2012]). Therefore,
P∞
R duxe
R∞
F̄ (byc)dy − ux F̄ (byc)dy
Pr(X ≥ ux)
i=duxe pX (i)
ux
R∞
=
= P∞
Pr(X ≥ u)
i=u pX (i)
u F̄ (byc)dy
R∞
u x F̄ (buzc)/F̄ (buc)dz − (duxe − ux)F̄ (buxc)/F̄ (buc)
R∞
=
→ x−1/ξ , x ≥ 1,
u 1 F̄ (buzc)/F̄ (buc)dz
applying (9) and dominated convergence. Thus, F̄X ∈ RV−1/ξ and X ∈ MDAξ .
For the second part of the theorem, we have
P
−1/ξ−1
F̄ (u + ku )/F̄ (u) ∞
P (X = ku + u | X ≥ u)
i=0 (1 + i/u)
P
=
.
pGZD (ku ; ξu, ξ)
(1 + ku /u)−1/ξ−1 ∞
i=0 F̄ (u + i)/F̄ (u)
By the uniform convergence (8) and the the fact that ku grows at most linearly fast, we
conclude that
F̄ (u + ku )/F̄ (u)
→ 1,
(1 + ku /u)−1/ξ
as u → ∞ over the integers. Second, Lemma 2.2 yields
u−1

∞
X

(1 + i/u)−1/ξ−1 → ξ.

i=0

6

Third, it follows from (9) that for  ∈ (0, 1/ξ), there exists u > 0 such that for u ≥ u ,

∞
∞ 
X
X
ξe
i −1−1/ξ+
→
,
u−1
F̄ (u + i)/F̄ (u) ≤ u−1 e
1+
u
1 − ξ
i=0

i=0

using again Lemma 2.2. A lower bound is found in the same manner and we let  → 0 to
conclude.
We now present a tail equivalence property between the probability mass and density
functions of the GZD, D-GPD and GPD. A direct consequence is that the denominator
pGZD in (10) can be replaced either by pD-GPD or by fGPD .
Proposition 2.4. If ξ ≥ 0, then
lim

sup

σ→∞ k=0,1,2,...

pD-GPD (k; σ, ξ)
pD-GPD (k; σ, ξ)
− 1 = lim sup
− 1 = 0.
σ→∞
pGZD (k; σ, ξ)
k=0,1,2,... fGPD (k; σ, ξ)

Proof. Suppose first that ξ > 0. Then
−1/ξ
−1/ξ
1 + ξ σk
− 1 + ξ k+1
pD-GPD (k; σ, ξ)
σ
=

1
k −1/ξ−1
fGPD (k; σ, ξ)
1
+
ξ
σ
σ
(

−1/ξ )
ξ
= 1− 1+
(σ + ξk) → 1,
σ + ξk
uniformly in k = 0, 1, 2, . . .. Furthermore,
∞
X
fGPD (k; σ, ξ)
= σ −1
(1 + ξi/σ)−1/ξ−1 → 1
k=0,1,2,... pGZD (k; σ, ξ)

sup

i=0

by Lemma 2.2. In the case ξ = 0,
pD-GPD (k; σ, 0)/fGPD (k; σ, 0) = pGZD (k; σ, 0)/fGPD (k; σ, 0) = σ(1 − e−1/σ ) → 1
as σ → ∞.
Next we extend Theorem 2.3 to the case ξ = 0. Recall that a distribution F is in MDA0
if and only if the survival function has a representation
 Z x

1
F̄ (x) = c(x) exp −
dy , −∞ < x < xF ,
(11)
0 a(y)
where c(·) is a positive function with c(x) → c > 0 as x ↑ xF , and a(·) is a positive,
differentiable function a(·) with limx↑xF a0 (x) = 0. If c(x) = c on (d, xF ) for some d < xF ,
then we say that the distribution F satisfies the von Mises condition. The function a(·)
in (11) is sometimes called the auxiliary function. Note, however, that it is only uniquely
defined (on (d, xF )) under the von Mises condition; see Embrechts et al. [2013]. Recall that
in the sequel we only consider the case of unbounded support, i.e. xF = ∞.
7

Theorem 2.5. Suppose that pX ∈ D-MDA0 and, moreover, that a distribution F such that
pX (k) = F̄ (k) has the property that an auxiliary function of F̄ satisfies limx→∞ a(x) = σ ∈
[0, ∞]. Then
lim

u∈N0 , u→∞

P (X = k + u | X ≥ u) = pGe (k; σ),

k = 0, 1, 2, . . . ,

(12)

where pGe (k; σ) = (1 − e−1/σ ) e−k/σ is the probability mass function of a geometric distribution if 0 < σ < ∞, and pGe (k; ∞) = pGe (k; 0) = 0. Furthermore, if σ ∈ [0, ∞), then
X ∈ D-MDA0 .
Proof. Note that for large integers u,
F̄ (k + u)/F̄ (u)
.
P (X = k + u | X ≥ u) = P∞
i=0 F̄ (i + u)/F̄ (u)

(13)

We have for every i = 0, 1, 2, . . .,
 Z i

c(i + u)
exp −
1/a(u + y)dy → e−i/σ
F̄ (i + u)/F̄ (u) =
c(u)
0
as u → ∞. If 0 < σ < ∞, then the dominated convergence theorem gives us
∞
X

F̄ (i + u)/F̄ (u) →

∞
X

e−i/σ = 1/(1 − e−1/σ ) ,

i=0

i=0

and (12) follows. If σ = ∞,
∞
X

F̄ (i + u)/F̄ (u) → ∞

i=0

by Fatou’s lemma, and (12), once again, follows. If σ = 0, the claim follows from the fact
that the denominator in (13) cannot be smaller than 1.
For the second part of the proposition, it follows from
 Z n

1
pX (n) = c(n) exp −
dy
0 a(y)
for all n and a(y) → σ ∈ [0, ∞) that
pX (n)
= 1 − e−1/σ ∈ (0, ∞) ,
n→∞ P (X ≥ n)
lim

which immediately implies that X ∈ D-MDA0 as well.
To summarize, for a discrete random variable X, the conditions X ∈ MDA, X ∈
D-MDA and pX ∈ D-MDA are related to each other as follows. If ξ ≥ 0, then X ∈ MDAξ
if and only if X ∈ D-MDAξ and X is long-tailed. If ξ > 0, then pX ∈ D-MDAξ/(1+ξ) implies
8

X ∈ D-MDAξ ; the same implication holds in the case ξ = 0 if the auxiliary function of the
extension of pX satisfies a(x) → σ ∈ (0, ∞) as x → ∞.
The condition pX ∈ D-MDA is satisfied, among others, by the Zipf–Mandelbrot, geometric, Poisson and negative binomial distributions, as shown in the next example.
Example 2.6. The probability mass function of a Zipf–Mandelbrot distribution with parameters s > 1 and q > 0 is in D-MDA1/s because it is regularly varying of order −s.
The probability mass function of a geometric distribution belongs to D-MDA0 as it
coincides up to a constant with the survival function of an exponential distribution. The
latter distribution clearly satisfies the von Mises condition and thus is a member of MDA0 .
The auxiliary function is, in fact, equal (eventually) to 1/λ, where λ is the rate of the
exponential distribution.
The probability mass function p of a Poisson distribution with rate λ > 0 coincides on
k = 0, 1, 2, . . . with the function
λx e−λ
,
g(x) =
Γ(x + 1)
a continuous function on R+ satisfying limx→∞ g(x) = 0. Moreover,
d
log g(x) = −ψ0 (x + 1) + log λ ,
dx
where ψ0 is the polygamma function of order 0. Since ψ0 (x) → ∞ as x → ∞, we see that
g 0 (x) < 0 for x sufficiently large. Therefore, F̄Y (x) = g(x)/g(d) is a survival function on
[d, ∞) for some d ≥ 0. Furthermore,


1
ψ1 (x + 1)
d
− 0
=−
,
dx
g (x)
(ψ0 (x + 1) − log λ)2
where ψ1 = ψ00 is is the polygamma function of order 1. Since ψ1 (x) → 0 as x → ∞,
we conclude that F satisfies the von Mises condition, with the auxiliary function a(x) =
(ψ0 (x + 1) − log λ)−1 → 0 as x → ∞. Therefore, the Poisson probability mass function is
in D-MDA0 .
Similarly, the probability mass function of the negative binomial distribution with probability of success p ∈ (0, 1) and number of successes r > 0 is also in D-MDA0 because it
coincides on {0, 1, 2, . . .} with the function
g(x) =

pr Γ(x + r)
(1 − p)x ,
Γ(r) Γ(x + 1)

a continuous function on R+ . It is simple to check that limx→∞ g(x) = 0, and g 0 (x) < 0 for
x large enough, so that F̄Y (x) = g(x)/g(d) is a survival function on [d, ∞) for some d ≥ 0.

9

Furthermore, g(x) ∼ cxr−1 (1 − p)x for large x, where c is a positive constant. Therefore,
F̄Y is of the form (11) with the auxiliary function
a(x) =

1
, x large,
− log(1 − p) − (r − 1)/x

and so it converges to −1/ log(1 − p) as x → ∞.
We conclude this section with a discussion designed to provide some intuition on how
the approximation methods suggested above differ, assuming that both apply in a given
situation. First of all, one would expect the D-GPD and GZD approximations to perform
similarly when ξ is close to zero because both approximating distributions coincide with a
geometric distribution when ξ = 0. Second, Proposition 2.4 suggests that, regardless of the
underlying justification, if one uses either pD-GPD (k; σ, ξ), fGPD (k; σ, ξ) or pGZD (k; σ, ξ) as
an approximation to P (X − u = k | X ≥ u), one should not expect major differences as
long as one chooses the scale parameter σ to be large. This would always be the case if X
is long-tailed, since the scale parameter is chosen to be proportional to the normalization
sequence au defined in (1), which grows to infinity if and only if X is long-tailed.
When using a continuous distribution, such as the generalized Pareto distribution, to
approximate the probabilities related to a discrete distribution, it is also common to use a
“continuity correction” and shift the argument in the continuous approximation by some
δ ∈ [0, 1). In our situation this amounts to replacing fGPD (k; σ, ξ) by fGPD (k + δ; σ, ξ),
some δ ∈ [0, 1). When ξ > 0, it is elementary to check that, as σ → ∞,
pD-GPD(σ,ξ) (k)
(1 + ξ)(2δ − 1)
=1+
+ O(σ −2 )
fGPD(σ,ξ) (k + δ)
2σ
for every k = 0, 1, 2, . . .. Therefore, the approximations by pD-GPD(σ,ξ) (k) and fGPD (k +
δ; σ, ξ) with large σ are most similar when δ = 1/2, a property that will be illustrated in
the empirical part. Similarly, in the case ξ = 0, as σ → ∞,
pD-GPD(σ,ξ) (k)
2δ − 1
= σeδ/σ (1 − e−1/σ ) = 1 +
+ O(σ −2 )
fGPD(σ,ξ) (k + δ)
2σ
for every k = 0, 1, 2, . . ., and the fastest convergence to the unity is again found when
δ = 1/2.
The final result of this section accomplishes two things. Its first part is related to the
well-known invariance property of the generalized Pareto distribution: its residual lifetime
is again generalized Pareto distributed. More precisely, if Y has the generalized Pareto
distribution with scale parameter σ and shape parameter ξ ≥ 0, then the exceedance Y − u
has, given Y ≥ u, the generalized Pareto distribution with scale parameter σ+ξu and shape
parameter ξ, for all u ≥ 0. This invariance property is important when approximating
10

the exceedance distribution using generalized Pareto distributions because changing the
threshold does not alter the distributional assumptions used in the approximation. It is
easily checked that the D-GPD has an analogous property. Moreover, the D-GPD also
exhibits the property that discretizing a GPD using different types of rounding does not
affect the fact that a D-GPD is obtained, and the shape parameter remains invariant.
Proposition 2.7. Let Y follow a GPD with scale parameter σ > 0 and shape parameter
ξ ≥ 0. Let 0 < h ≤ 1. If X = bλY + 1 − hc, then for any integer u ≥ 1 − h, the distribution
of X −u | X ≥ u is a D-GPD with shape parameter λσ +ξ(u+h−1) and scale parameter ξ.

3

Empirical Results

0.0

0.5

1.0

1.5

2.0

2.5

3.0

0.0 0.5 1.0 1.5 2.0 2.5 3.0

●

Fitted Quantiles

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●

●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

2.0
1.5
1.0
0.0

0.5

Fitted Quantiles

2.5

3.0

We now assess the performance of the generalized Pareto distribution (GPD), the discrete
generalized Pareto distribution (D-GPD) and the generalized Zipf distribution (GZD) approximations for estimating extreme quantiles from several simulated and real data sets,
starting with a simple example to illustrate the methods.

●

Empirical Quantiles

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

0.0

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

0.5

1.0

1.5

2.0

2.5

3.0

Empirical Quantiles

Figure 1: Quantile-quantile plots for the fit of a GPD (left) and D-GPD (right) to X − 3 | X ≥ 3,
where X is Poisson distributed with rate 1 and the sample size is 5000. A good fit occurs on the
right because the points are contained between the two lines, but not on the left.

3.1

Simulated Data

Let X be Poisson distributed with rate λ = 1 and consider an i.i.d. sample of size n = 5000.
We are interested in inferring the distribution of the exceedances X − u | X ≥ u above a
11

large threshold u, say u = 3 which is the 95th empirical percentile of the data in this
case. Since X ∈ D-MDA0 , we can approximate the distribution of X − u | X ≥ u
by a D-GPD as explained in Section 1. Moreover, pX , the probability mass function
of X, also belongs to D-MDA0 and satisfies the additional assumption in Theorem 2.5 as
shown in Example 2.6, motivating the D-GPD and GZD approximations. However, X is
not in MDA and thus the GPD approximation does not necessarily apply. In order to
compare these three approximations, we fit a GPD, D-GPD and GZD to the observations
above u = 3, estimating their two parameters σ and ξ by maximum likelihood. Figure 1
compares the fitted quantiles of the GPD (left) and D-GPD (right) relative to the empirical
quantiles. Since the data are discrete, we applied a slightly different graphical method
−1
than the standard quantile-quantile (QQ) plots: FGPD
{i/(n + 1); σ, ξ}, for i = 1, . . . , n,
were plotted against empirical quantiles, where (σ, ξ) are the estimated parameters. A
good fit occurs when the points accumulate between the two diagonal lines in Figure 1
and touch the bottom line, as would be the case if the quantiles of a continuous random
variable Y were plotted against those of bY c. One clearly sees that the D-GPD fits well the
observations compared to the GPD which delivers a poor fit. We mention that the D-GPD
and GZD approximations produce very similar estimates in this case; the QQ-plot of the
GZD (not displayed here) looks visually identical to the one of the D-GPD. The estimated
scale parameters σ for the GPD, D-GPD and GZD are 1.91, 0.71 and 0.69 respectively.
Increasing λ would increase these estimates and render the three methods indistinguishable
as expected from Proposition 2.4.
We now compare the performances of the GPD, D-GPD and GZD approximations in
estimating the probability of a rare event in the following simulated case. Let
Y ∼ IG(2, 1),

X = bY c,

(14)

where IG(α, β) denotes an inverse-gamma distribution with probability density function
f (x) = Γ(α)−1 β α x−α−1 exp(−β/x), x > 0. We repeat 500 times the experiment described
below. An i.i.d. sample of size 8000 is drawn from the distribution of X. From these
observations, the goal is to estimate the probability of the extreme region
pe = P (X ≥ bqe c),

(15)

where qe is the 99.99 percentile of Y, i.e., the value exceeded once every 10 000 times
on average. The strategy pursued is to select an integer threshold u as the 95th empirical
percentile of the sample, fit a parametric distribution to the exceedances X −u | X ≥ u, and
use it to extrapolate the tail and estimate pe . It holds X ∈ D-MDA1/2 and X ∈ MDA1/2 ,
which motivates the choice of a GPD and D-GPD as seen in Section 1. We implement two
variants of the GPD approximation: the first has no continuity correction and the second
12

250
150
100
50

Frequency

200

100 150 200 250 300

0

0

50

Frequency

0

20

40

60

80

100

100

1000

10000

Word Frequency

Figure 2: On the left: frequency plot of a sample of X − 2 | X ≥ 2 of size 701 for X defined in
(14). On the right: frequency plot of the 5588 most frequent words in a British corpus (x axis on
log-scale).

shifts the observations by −δ = − 21 . It is not immediate if the probability mass function
of X is in D-MDA and we thus apply the GZD approximation heuristically.
As a benchmark, we will also estimate pe from a sample of the continuous variable Y
(as opposed to its discretization X). In this context, the GPD approximation is motivated
by the fact that Y ∈ MDA1/2 , and we thus fit a GPD to Y − u | Y ≥ u.
A frequency plot of a sample of X − u | X ≥ u is displayed in Figure 2 on the lefthand side. For each model, we compute maximum likelihood estimators for σ and ξ by
performing a two dimensional maximization using the function optim of R (R Core Team
[2015]) with starting values (1, 1). We then compute pe and approximate 90% confidence
intervals under asymptotic normality of the estimators. Table 1 displays: the average
parameters pe , ξ and σ over the 500 experiments, the average length of the confidence
intervals, the true length and their coverage. True length is the length the intervals should
have had to contain the estimates across the 500 experiments 90% of the time. Coverage
indicates the proportion of time the truth lies in the confidence interval.
It appears that the D-GPD and GZD approximations accurately estimate pe from the
discretized data with a coverage close to the correct one of 90%, and that their performance
is good relative to the situation of full information where the continuous data are available.
On the other hand, the GPD approximation with δ = 21 is inaccurate and the GPD
with δ = 0 performs very poorly (misleadingly, the latter has a larger likelihood at the

13

pe × 103
truth
Y −u|Y ≥u
GPD
X −u|X ≥u
D-GPD
GZD
GPDδ= 21
GPDδ=0

ξ

σ

mean (cov, len)
0.10

true length

mean (cov, len)
0.5

mean (len)

0.10 (87%, 0.16)

0.16

0.49 (95%, 0.22)

1.19 (0.30)

0.10 (86%, 0.17)
0.11 (88%, 0.17)
0.04 (20%, 0.07)
7.93 (83%, 23.97)

0.16
0.17
0.08
11.25

0.49 (95%, 0.23)
0.50 (95%, 0.24)
0.37 (22%, 0.18)
8.27 (0%, 1.24)

1.19 (0.33)
1.39 (0.29)
1.43 (0.32)
0.00 (0.00)

Table 1: Performance of several methods in estimating the probability pe of a rare event defined
in (15). A generalized Pareto distribution (GPD) is fitted to the observations of Y exceeding
u = 2, and the following distributions are fitted to exceedances of X = bY c: a discrete generalized
Pareto distribution (D-GPD), a generalized Zipf distribution (GZD) and a GPD with continuity
correction δ. The table displays average maximum likelihood estimators for pe , ξ and σ accross
500 experiments. Coverage (cov) and average length (len) of 90% confidence intervals are shown
between brackets. In each experiment, about 700 observations exceeded the threshold.

maximum likelihood estimate than the former). This illustrates why the D-GPD or GZD
approximations should be preferred to the GPD approximation in situations similar to this
simulated case.
We know that the D-GPD and GZD coincide with a geometric distribution when ξ = 0
and that they are asymptotically equivalent when ξ > 0 and σ → ∞ (Proposition 2.4).
Although σ is not particularly large in the above simulated case, the D-GPD and GZD
approximations deliver very similar results. We point out that they outperform the Poisson
and negative binomial distributions which would estimate pe very poorly.
We remark few differences between the GPD fitted to exceedances of Y and the D-GPD
fitted to exceedances of its discretization X, except that the former method yields slightly
shorter confidence intervals. As a consequence of the invariance property in Proposition 2.7,
the two methods would still give similar estimates of ξ if X would be defined by other forms
of rounding than X = bY c, such as using rounding or ceiling function.
In conclusion, we showed that the D-GPD and GZD approximations were able to estimate the probability of rare events in this example. These findings are supported by
two complementary simulated cases with Y ∈ MDAξ for ξ = 0 and ξ < 0 (Hitz [2016],
Chapter 2). The D-GPD provides a more efficient method because its probability mass,
survival and quantile functions are closed-form, but both distributions are useful to describe extreme values of discrete random variables, which we will further illustrate on real
data sets.

14

600

800

●
●

0

●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

400

Fitted Quantiles

10000

●

200

50000
30000

●

●

0

Fitted Quantiles

●

0 10000

30000

50000

●●
●
●●
●●
●●
●
●
●
●●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

0

Empirical Quantiles

200

400

600

800

Empirical Quantiles

Figure 3: QQ-plots for a D-GPD fitted to the frequencies of the 5588 most frequent words in a
British corpus. On the right, only percentiles below 99% are plotted. Dashed lines denote pointwise
90% confidence intervals.

3.2

Word Frequencies and Word Lengths

We consider a data set based on the British National corpus counting written and spoken
English words (British National Corpus, [2007]) and are interested in modeling the frequency of the most popular words. The six most frequent words in the corpus are: “the”,
“of”, “and”, “a”, “in” and “to.” Let X be the frequency at which a word occurs. We select
a large threshold u and fit a GPD, a D-GPD and a GZD to X − u | X ≥ u by maximum likelihood. The GPD is implemented with continuity correction (δ = 12 ) and without
(δ = 0). Selecting an appropriate threshold is crucial when estimating high quantiles and
can be based on techniques such as mean residual plots (see e.g. Davison and Smith [1990]).
As our focus here is rather on describing the tail distribution of these word frequencies,
we choose a relatively low threshold exceeded by 5588 observations. A frequency plot of
X | X ≥ u is displayed on the right-hand side in Figure 2.
The D-GPD delivers a good fit as revealed by the QQ-plot in Figure 3; indeed, most
observations lie within the pointwise 90% confidence intervals obtained by simulating 2000
times from the fitted model and computing empirical quantiles at each simulation. Table 2
presents maximum likelihood estimators for σ and ξ with 90% confidence intervals. Apart
from the case of the GPD with δ = 0, all models provide similar results as expected from
Proposition 2.4 when σ is large. A Kolmogorov–Smirnov test for discrete data (Arnold
and Emerson [2011]) leads to the same conclusion (in order to perform the test for the

15

0

●
●
●
●
●
●
●
●
●

10
8
6
2

Fitted Quantiles

●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●

●
●
●
●
●

0

●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●

2

●

4

12
8
6
4
0

2

Fitted Quantiles

10

●
●
●
●
●
●

4

6

8

10

●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●

0

Empirical Quantiles

●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●

2

●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●

4

●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●

6

●
●
●
●
●

●
●
●
●
●

8

Empirical Quantiles

Figure 4: QQ-plots for a D-GPD (left) and GPD with δ = 12 (right) fitted to the length of the
2875 French words consisting of 15 letters or more. A good fit occurs when the lowest part of each
accumulation of points is close to the line.

GPD, which is a continuous distribution, we assumed that data were rounded realizations
of the fitted model). The analysis of frequencies of French words in a collection of books
and movie subtitles leads to analogous results (Hitz [2016], Chapter 2). Since the GZD
coincides with a Zipf–Mandelbrot distribution when ξ > 0, the above results are consistent
with a common hypothesis in linguistic that word frequencies follow a Zipf-type law (see
e.g. Booth [1967]).
We now consider the set of 150 000 words in the French lexical (New et al. [2004]) and
denote by X the length of a word. The longest French word, for instance, is “anticonstitutionnellement”, consisting of 25 letters. Contrary to the word frequencies, this data
set contains many tied observations and we want to see if this translates into a marked
difference between the methods. We thus fit the models to X − u | X ≥ u with u = 15
(the 98% empirical percentile of the data), leaving 2875 exceedances. The D-GPD and
GZD deliver a good fit and similar estimations between each other, and this time they
clearly outperform the GPD approximation with δ = 21 as shown by QQ-plots in Figure 4
and discrete Kolmogorov–Smirnov tests in Table 2 (p-values are here computed by Monte
Carlo simulation). Notice that the negative binomial also fits well the observations in this
case.
To summarize, we have illustrated the adequacy of the D-GPD and GZD in describing
the frequencies of the most common and longest words from large corpora. This data
analysis supports the conclusion drawn earlier from the simulated case: the D-GPD and

16

Model

p-value

Word Frequency
British
D-GPD
GZD
GPDδ= 12
GPDδ=0
Negative binomial

0.56
0.56
0.55
0.00
0.00

Word length
French
D-GPD
GZD
GPDδ= 12
Negative binomial

0.85
0.84
0.00
0.75

NLL

ξ

σ

27896.6
27896.6

0.88[0.84,0.93]
0.88[0.84,0.93]
0.88[0.84,0.93]
0.93[0.89,0.97]

22.38[21.42,23.34]
22.82[21.87,23.77]
22.39[21.43,23.34]
20.89[19.96,21.81]

0.02[−0.01,0.06]
0.02[−0.01,0.06]
−0.04[−0.06,−0.01]

1.36[1.30,1.43]
1.37[1.32,1.43]
1.51[1.45,1.57]

29663.2

3894.0
3894.0
3893.9

Table 2: Fit of a GPD, D-GPD and GZD to frequencies of the most frequent words in a British corpus and length of the longest French words. The table displays the p-value of discrete Kolmogorov–
Smirnov tests, negative log-likelihood (NLL) and maximum likelihood estimators with 90% confidence intervals. The sample sizes are 5588 for word frequencies and 2875 for word lengths.

GZD are preferred over the GPD to model extremes of discrete data when tied observations
are frequent.

3.3

Tornadoes

Accurately assessing the risk of environmental hazards is crucial for insurance companies in
particular, and we illustrate here how the D-GPD and GZD approximations can be useful
techniques for this purpose.
We consider the data set studied in Tippett et al. [2016] which reports the number of
extreme tornadoes per outbreak in the United States between 1965 and 2015. They defined
an “extreme outbreak” as a sequence of twelve or more tornadoes occurring close to one
another in time and that are rated F1 and greater on the Fujita scale, or EF1 and greater
on the enhanced Fujita scale [Fuhrmann et al., 2014].
Let X be the number of such tornadoes per extreme outbreak. The authors found that
observations from X − u | X ≥ u for u = 12 were well modeled by a GPD with linear
temporal trend in the scale parameter and continuity correction δ = 12 . The GPD, however,
is not a discrete distribution and the D-GPD and GZD seem appropriate choices for this
type of data. We fit these three distributions and find that they all achieve comparable
quality of fit as shown in Table 3; linear trend in σ leads to significative improvements in
17

each case.
Model
D-GPD
GPDδ= 12

NLL
1439.92
1439.93

ξ
0.27[0.16,0.37]
0.26[0.16,0.37]

σ0
4.81[3.64,5.99]
4.86[3.68,6.04]

σt
6.11[3.74,8.48]
6.13[3.75,8.50]

Table 3: Fit of a D-GPD and GPD (with linear trend in the scale parameter) to the number of
tornadoes per extreme tornado outbreak in the United States. Only the 435 outbreaks with more
than 12 tornadoes were considered. The table displays negative log-likelihood (of the discretized
model in the case of the GPD) and maximum likelihood estimates with 90% confidence intervals.

3.4

Multiple Births

We now turn our attention to a data set consisting of very small integer values to see if
the D-GPD and GZD can describe tail distributions in this special case. We examine data
counting multiple births in the United States from 1995 to 2014 (Hamilton et al. [2015]);
its frequency table reads
single
78 178 588

twin
2 500 340

triplet
117 603

quadruplet
8 108

quint. or more
1 353.

Let X be the number of children at birth. We fit a right-censored D-GPD, GZD,
negative binomial and Poisson distribution to X C − u | X C ≥ u for u = 1, where X C =
min(X, 5). As shown in the table below, which displays Bayesian information criteria (BIC)
and maximum likelihood estimates, the D-GPD and GZD outperform the Poisson and
negative binomial distributions, and are useful methods if one must estimate, for example,
the probability that an American women delivers sextuplets.

D-GPD
GZD
Negative binomial
Poisson

BIC
546 441.2
546 440.6
546 547
552 284

ξ
0.06[0.06,0.06]
0.06[0.06,0.07]

σ
0.30[0.30,0.30]
0.32[0.32,0.32]

The applicability of the D-GPD and GZD approximations for estimating tail distributions when the selected threshold u is a small integer should be more rigorously explored.
Future work could also assess the validity of the approximations in the case ξ < 0, and
compare them to a broader class of discrete distributions such as discrete compound Poisson distributions. Since the D-GPD and GZD delivered similar performances in the data
18

analysis carried out in this article, it would be interesting to further understand how they
relate to one another.

Acknowledgements The first author is grateful to the Berrow Foundation for financial support and would like to thank Robin Evans, Gesine Reinert, David Steinsaltz and
Jonathan Tawn for valuable advice and encouragement. This research was partially supported by the ARO grant W911NF-12-10385.

References
C. W. Anderson. Extreme value theory for a class of discrete distributions with applications to
some stochastic processes. Journal of Applied Probability, 7:99–113, 1970.
C. W. Anderson. Local limit theorems for the maxima of discrete random variables. Mathematical
Proceedings of the Cambridge Philosophical Society, 88(1):161–165, 1980.
B. C. Arnold. Pareto distribution. International Cooperative Publishing House, Maryland, 1983.
T. B. Arnold and J. W. Emerson. Nonparametric goodness-of-fit tests for discrete null distributions.
The R Journal, 3(2):34–39, 2011.
R. L. Axtell. Zipf distribution of US firm sizes. Science, 293:1818–1820, 2001.
N. H. Bingham, C. M. Goldie, and J. L. Teugles. Regular Variation. Cambridge University Press,
1989.
A. D. Booth. A law of occurrences for words of low frequency. Information and Control, 10(4):
386–393, 1967.
British National Corpus. Version 3 BNC XML edition, 2007.
A. Buddana and T. J. Kozubowski. Discrete Pareto distributions. Economic Quality Control, 29
(2):143–156, 2014.
A. Clauset, C. R. Shalizi, and M. E. J. Newman. Power-law distributions in empirical data. SIAM
Review, 51(4):661–703, 2009.
A. C. Davison and R. L. Smith. Models for exceedances over high thresholds. Journal of the Royal
Statistical Society. Series B. Methodological, 52(3):393–442, 1990.
P. Embrechts, C. Klüppelberg, and T. Mikosch. Modelling extremal events: for insurance and
finance. Springer, New York, 2013.
C. M. Fuhrmann, C. E. Konrad, M. M. Kovach, J. T. McLeod, W. G. Schmitz, and P. G. Dixon.
Ranking of tornado outbreaks across the United States and their climatological characteristics.
Weather and Forecasting, 29(3):684–701, 2014.

19

X. Gabaix. Zipf’s law and the growth of cities. The American Economic Review, 89(2):129–132,
1999.
B. E. Hamilton, J. A. Martin, M. J. Osterman, S. C. Curtin, and M. TJ. Births: Final data for
2014. National Vital Statistics Reports, 64(1), 2015.
A. S. Hitz. Modelling of Extremes. PhD thesis, University of Oxford, 2016.
T. J. Kozubowski, A. K. Panorska, and M. L. Forister. A discrete truncated Pareto distribution.
Statistical Methodology, 26:135–150, 2015.
H. Krishna and P. S. Pundir. Discrete Burr and discrete Pareto distributions. Statistical Methodology, 6(2):177–188, 2009.
B. Mandelbrot. Contribution à la théorie mathématique des jeux de communication. Publications
de l’Institut de statistique de l’Université de Paris, 1953.
B. New, C. Pallier, M. Brysbaert, and L. Ferrand. Lexique 2: A new French lexical database.
Behavior Research Methods, Instruments, & Computers, 36(3):516–524, 2004.
J. Pickands, III. Statistical inference using extreme order statistics. The Annals of Statistics, 3:
119–131, 1975.
F. Prieto, E. Gómez-Déniz, and J. M. Sarabia. Modelling road accident blackspots data with the
discrete generalized Pareto distribution. Accident Analysis & Prevention, 71:38–49, 2014.
R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for
Statistical Computing, 2015.
S. I. Resnick. Extreme values, regular variation, and point processes. Springer, New York, 1987.
T. Shimura. Discretization of distributions in the maximum domain of attraction. Extremes, 15(3):
299–317, 2012.
M. K. Tippett, C. Lepore, and J. E. Cohen. More tornadoes in the most extreme US tornado
outbreaks. Science, 354(6318):1419–1423, 2016.

20

