Maximum likelihood estimators based on the block maxima method
Abstract
The extreme value index is a fundamental parameter in univariate Extreme Value Theory (EVT). It captures the tail behavior of a distribution and is central in the extrapolation beyond observed data. Among other semiparametric methods (such as the popular Hill’s estimator), the Block Maxima (BM) and PeaksOverThreshold (POT) methods are widely used for assessing the extreme value index and related normalizing constants. We provide asymptotic theory for the maximum likelihood estimators (MLE) based on the BM method. Our main result is the asymptotic normality of the MLE with a nontrivial bias depending on the extreme value index and on the socalled second order parameter. Our approach combines asymptotic expansions of the likelihood process and of the empirical quantile process of block maxima. The results permit to complete the comparison of most common semiparametric estimators in EVT (MLE and probability weighted moment estimators based on the POT or BM methods) through their asymptotic variances, biases and optimal mean square errors.
Key words: asymptotic normality, block maxima method, extreme value index, maximum likelihood estimator, peaksoverthreshold method, probability weighted moment estimator.
AMS Subject classification: 62G32, 62G20, 62G30.
1 Introduction
The Block Maxima (BM) method, also known as Annual Maxima, after Gumbel [11] is a fundamental method in Extreme Value Theory and has been widely used. The method is justified under the Maximum Domain of Attraction (MDA) condition: for an independent and identically distributed (i.i.d.) sample with distribution function , if the linearly normalized partial maxima converges in distribution, then the limit must be a Generalized Extreme Value (GEV) distribution. In practice, one rarely exactly knows but the MDA condition holds for most common continuous distributions.
In the BM method, the initial sample is divided into blocks of the same size and the MDA condition ensures that the block maxima are approximately GEV distributed. The method is commonly used in hydrology and other environmental applications or in insurance and finance when analysing extremes  see e.g. the monographs by Embrechts et al. [9], Coles [5], Beirlant et al. [2], de Haan and Ferreira [6] and references therein.
The GEV is a three parameter distribution, with the usual location and scale parameters, and the extreme value index being the main parameter as it characterizes the heaviness of the tail. Several estimation methods have been proposed, including the classical maximum likelihood (ML) and probability weighted moments (PWM) estimators (Hosking et al. [13]). The asymptotic study of these estimators has been established for a sample from the GEV distribution and asymptotic normality holds with null bias and explicit variance (Prescott and Walden [15], Hosking et al. [13], Smith [17], Bücher and Segers [4]). The theory is made quite difficult and technical by the fact that the support of the GEV is varying with respect to its parameters. Regularity in quadratic mean of the GEV model has been proven only recently by Bücher and Segers [4] and we provide here a different and somewhat simpler proof (cf. Proposition 4.1).
However, in applications, the sample block maxima are only approximately GEV so that the classical parametric theory suffers from model misspecification. In this paper, we intend to fill this gap for ML estimators (MLE), by showing asymptotic normality under a flexible second order condition (a refinement of the MDA condition). Depending on the asymptotic block size, a non trivial bias may appear in the limit for which we provide an exact expression. Recently Ferreira and de Haan [10] showed asymptotic normality of the PWM estimators under the same conditions. They derived a uniform expansion for the empirical quantile of block maxima that is a crucial tool in our approach as well. Indeed, the MLE can be seen as a maximizer of the socalled likelihood process. Expressing the likelihood process in terms of this empirical quantile process, we are able to derive an expansion of the likelihood process that implies the asymptotic normality of the MLE. This derivation is again made quite technical by the fact that the support of the GEV is varying. Note that the asymptotic normality for the MLE of a Fréchet distribution based on the block maxima of a stationary heavytailed time series has been obtained by Bücher and Segers [3]. There the issue of parameter dependent supports is avoided but time dependence has to be dealt with. Besides, the ideas underlying their proof are quite different.
The asymptotic normality result in the present paper brings novel results to the theoretical comparison of the main semiparametric estimation procedures in EVT. On the one hand it permits to compare BM and PeaksoverThreshold (POT) methods (see e.g. Balkema and de Haan [1], Pickands [14]), the latter being another fundamental method in EVT and concurrent with BM. We discuss and compare the four different approaches – MLE/PWM estimators in the BM/POT approaches – based on exact theoretical formulas for asymptotic variances, biases and optimal mean square errors depending on the extreme value index and the second order parameter. It turns out that MLE under BM has minimal asymptotic variance among all combinations MLE/PWM and BM/POT but, on the other hand it has some significant asymptotic bias. When analysing the asymptotic optimal mean square error that balances variance and bias, the most efficient combination turns out to be MLE under POT (e.g. Drees, Ferreira and de Haan 2004). It turns out that the optimal sample size is larger for POTMLE than for BMMLE, giving a theoretical justification to the heuristic that POT allows for a better use of the data than BM.
The outline of the paper is as follows: In Section 2 we present the main theoretical conditions and results including Theorem 2.2 giving the asymptotic normality of the MLE. In Section 3 we present a comparative study of asymptotic variances and biases, optimal asymptotic mean square errors and optimal samples sizes among all combinations MLE/PWM and BM/POT. In Section 4 we state additional theoretical statements, including the local asymptotic normality of MLE under the fully parametric GEV model, and provide all the proofs. Finally, Appendix A gathers some formulas for the information matrix and for the bias of BMMLE and Appendix B provides useful bounds for the derivatives of the likelihood function that are necessary for the main proofs.
2 Asymptotic behaviour of MLE
2.1 Framework and notations
The GEV distribution with index is defined by
and the corresponding loglikelihood by
(2.1) 
For , the formula is interpreted as . The three parameter model with index , location and scale is defined by the loglikelihood
(2.2) 
A distribution is said to belong to the maxdomain of attraction of the extreme value distribution , denoted by , if there exist normalizing sequences and such that
The main aim of the BM method is to estimate the extreme value index as well as the normalizing constants and . The setup is the following. Consider independent and identically distributed (i.i.d.) random variables with common distribution function . Divide the sequence into blocks of length and define the th block maximum by
(2.3) 
For each , the variables are i.i.d. with distribution function and by the maxdomain of attraction condition
(2.4) 
This suggests that the distribution of is approximately a GEV distribution with parameters . The method consists in pretending that the sample follows exactly the GEV distribution and in maximizing the GEV loglikelihood so as to compute the MLE. A particular feature of the method is that the model is clearly misspecified since the GEV distribution appears as the limit distribution of the block maxima as the block size tends to while in practice we have to use a finite block size. As seen afterwards, we quantify the misspecification thanks to the socalled second order condition that implies an asymptotic expansion of the empirical quantile process with a non trivial bias term. When plugging this expansion in the ML equations, we obtain a bias term for the likelihood process as well as for the MLE.
The (misspecified) loglikelihood of the sample is
(2.5) 
We say that an estimator is a MLE if it solves the score equations
which we write shortly in vectorial notation
(2.6) 
A main purpose of this paper is to study the existence and asymptotic normality of the MLE under the following conditions:

Second order condition: for some positive function and some positive or negative function with ,
(2.7) with . Note that necessarily and is regularly varying with index .

Asymptotic growth for the number of blocks and the block size :
(2.8)
2.2 Main results
Before considering the MLE, we focus on the asymptotic properties of the likelihood and score processes. For the purpose of asymptotic we introduce the local parameter :
(2.9) 
Set . The local loglikelihood process at is given by
(2.10)  
and, the local score process by
(2.11)  
Clearly, the score equation (2.6) rewrites in this new variable as .
In the following, denote the quantile function of the extreme value distribution , i.e.
(2.12) 
Proposition 2.1.
Assume conditions (2.7) and (2.8). Let be a sequence of positive numbers verifying, as ,
(2.13) 
Let be the ball of center and radius . Then, uniformly for ,
(2.14) 
with the Fisher information matrix
(2.15) 
As a consequence, the local loglikelihood process is strictly concave on with high probability.
Remark 2.1.
Our main result is the following Theorem establishing the asymptotic behavior of the local likelihood process and from which the existence and asymptotic normality of MLE will be deduced.
Theorem 2.1.
Assume conditions (2.7) and (2.8). Then, the local likelihood process satisfies uniformly for in compact sets
(2.16)  
(2.17) 
where
(2.18) 
i.e., is asymptotically Gaussian with variance equal to the information matrix and mean depending on the second order condition (2.7) through
(2.19) 
and on the asymptotic block size through from (2.8).
Remark 2.2.
Remark 2.3.
Equation (2.8) requires that both the number of blocks and the block size go to infinity with a relative rate measured by the second order scaling function and a parameter . When , the bias term disappears in (2.18); this corresponds to the situation where grows to infinity very quickly with respect to so that the block size is large enough and the GEV approximation (2.4) is very good.
Existence and asymptotic normality of the MLE can be deduced from Theorem 2.1, mainly by the argmax theorem. The concavity property stated in Proposition 2.1 plays an important role in the proof of existence and uniqueness.
Theorem 2.2.
Remark 2.4.
An interesting byproduct of the strict concavity stated in Proposition 2.1 is the convergence of numerical procedures for the computation of the MLE that are implemented in software. The NewtonRaphson algorithm is commonly used to solve numerically the score equation (2.6). Strict concavity of the objective function on a large neighbourhood of the solution ensures convergence of the algorithm with high probability as soon as the initial value belongs to this neighbourhood.
3 Theoretical comparisons: BM vs POT and MLE vs PWM
The POT method uses observations above some high threshold or top order statistic and the underlying approximate model is the Generalized Pareto distribution (Balkema and de Haan [1], Pickands [14]). Estimators of the shape parameter , as well as location and scale parameters have been proposed and widely studied, including MLE and PWM (Hosking and Wallis [12]). For their asymptotic properties  under basically the same conditions as under BM in Theorem 2.2  we refer to de Haan and Ferreira [6]. Asymptotic normality of PWM estimators under BM has been established only recently by Ferreira and de Haan [10] and a comparison of PWM estimators under BM and POT has been carried out. The aim of the present section is to include our new asymptotic results for MLE estimators under BM, completing the picture in the comparison of the four different cases BM/POT and MLE/PWM.
Recall that asymptotic normality of MLE (resp. PWM estimator) holds for (resp. ). The number of selected observations corresponds to the number of blocks in BM and of selected top order statistics in POT. Similarly as in Ferreira and de Haan [10], our comparative study is restricted to the range where second order conditions for BM and POT are comparable (cf. Drees et al. [8] or Ferreira and de Haan [10]). In the following we compare MLE/PWM under BM/POT methods through their: (i) asymptotic variances (VAR), (ii) asymptotic biases (BIAS), (iii) optimal asymptotic mean square errors (AMSE) and optimal number of observations minimizing AMSE ().
(i) Asymptotic variances
The asymptotic variance depends on only and is plotted in Figure 1 where straight lines stand for MLE and dashed lines for PWM estimators. Among all four different cases, BMMLE is the one with the smallest variance within its range. Moreover, for both estimators, BM has the lowest variance indicating that BM is preferable to POT when variance is concerned.
(ii) Asymptotic biases
The asymptotic biases depend on and and are shown in Figures 2–3: POTMLE is the one with the smallest bias also in absolute value when compared to BMMLE, contrary to what was observed for variance. This is in agreement with what has been observed when comparing BMPWM and POTPWM, also shown in Figures 2–3 already analysed in Ferreira and de Haan [10]. There is again the indication that POT method is favourable to BM when concerning bias.
(iii) Optimal asymptotic MSEs and optimal number of observations
Another way to compare the estimators that combines both variance and bias information is through mean square error. One can compare these for the optimal number of observations i.e., that value for which the asymptotic mean square error (AMSE) is minimal. Similarly as in Ferreira and de Haan [10], under the conditions of Theorem 2.2, we have
with a decreasing and regularly varying function such that . It follows in particular that the optimal is different but of the same order for both estimators and methods. As for the optimal AMSE, we have
When considering ratios of optimal AMSE (or optimal number of selected observations), the regularly varying function cancels out and the asymptotic ratio does not depend on but only on and .
Figure 4 shows the contour plots of the ratio for MLE and PWM estimators. It is surprising to see a reverse behaviour in both cases: in the range of parameters considered, POT is preferable when MLE are considered, while BM is mostly preferable for PWM estimators.
In Figure 5 are shown for comparing optimal AMSE among all combinations. The green surface corresponds to MLEPOT that has always the minimal optimal AMSE in the range of parameters considered. Finally, Figure 6 reports for MLE the asymptotic ratio of optimal numbers of selected observations, that is . We can see that the optimal number of observations is larger for POT, which is in agreement with the PWM case considered in previous studies.
4 Main proofs
We start by introducing some material that will be useful for the proofs. More technical material is still postponed to Appendices.
4.1 Local asymptotic normality of the GEV model
If the observations are exactly distributed, then the choice of constants
(4.1) 
ensures that the normalised block maxima are i.i.d. with distribution . The issue of model misspecification is irrelevant in that particular case.
In this simple i.i.d. setting, a key property in the theory of ML estimation is differentiability in quadratic mean (see e.g. van der Vaart [18, Chapter 7]). A statistical model defined by the family of densities is called differentiable in quadratic mean at the point if there exists a measurable function called the score function such that
The following Proposition corresponds to Proposition 3.2 in Bücher and Segers [4]. We provide a slightly different proof in the case .
Proposition 4.1.
The three parameter GEV model with loglikelihood defined in Equation (2.2) is differentiable in quadratic mean at if and only if . The score function is then given by .
Proof of Proposition 4.1.
The density of the 3parameter GEV model is given by
if and otherwise. In the case , the function
is continuously differentiable for every and the information matrix is well defined and continuous (see Appendix A or Beirlant et al. [2] page 169). Differentiability in quadratic mean of the GEV model follows by a straightforward application of Lemma 7.6 in Van der Vaart [18].
In the case , the function is not differentiable at points such that . Going back to the definition of differentiability in quadratic mean, we need to show that
(4.2) 
This is credible because for all , the relation
entails
For further reference, we note also that, for ,
(4.3) 
A rigorous proof of (4.2) is given below. Since , we have for in a neighbourhood of so that the density vanishes outside with the right endpoint of the distribution . We also introduce . For all , the function is twice continuously differentiable, whence Taylor formula entails
for some . Together with Equation (4.3), the formula and Proposition B.1, this yields the upper bound
for all and small enough. This entails
(4.4) 
It remains to estimate the contribution of the integral between and . Recall that vanishes for . We have
The first and second integral converge to as because and . The third integral converges also to because and so that the score is square integrable (its covariance matrix is ). We deduce
(4.5) 
The fact that differentiability in quadratic mean does not hold when is proved in Bücher and Segers [4] Appendix C. They observe that for ,
which rules out differentiability in quadratic mean. We omit further details here. ∎
Differentiability in quadratic mean implies that the score function is centered with finite variance equal to the information matrix, i.e.
(4.6) 
Another important consequence of differentiability in quadratic mean is the local asymptotic normality property of the local score process. The following Corollary follows from Proposition 4.1 by a direct application of Theorem 7.2 in Van der Vaart [18].
Corollary 4.1.
4.2 The empirical quantile process associated to BM
The starting point of the proof of Proposition 2.1 and Theorem 2.1 is to rewrite the local loglikelihood process (2.10) in terms of the (normalized) empirical quantile process
(4.7) 
where are the order statistics of the block maxima sample defined by (2.3) and denotes the smallest integer larger than or equal to . The local loglikelihood process (2.10) can be rewritten as
(4.8) 
Convergence (2.4) ensures the convergence of the empirical quantile process to the “true” quantile function defined in (2.12). The following expansion of the empirical quantile process is taken from Ferreira and de Haan [10], Theorem 2.1.
Proposition 4.2.
Remark 4.1.
For Proposition 4.2, the auxiliary functions and have to be specially chosen for establishing uniform second order regular variation bounds refining (2.7), see Lemma 4.2 in [10]. However, this choice is useful for the proofs only and is irrelevant for the statements of the main results in Section 2.2.
The following Proposition provides useful technical bounds for the proof of the main results.
Proposition 4.3.
For the proof of Proposition 4.3, we need the following Lemma.
Lemma 4.1.
Let be the order statistics of i.i.d random variables with standard Fréchet distribution. Then,
where the term is uniform for .
Proof of Lemma 4.1.
Proof of Proposition 4.3.
Let be a unit Fréchet random variable, i.e. with distribution function , , and be the order statistics from the associated i.i.d. sample of size , . Note that . From Lemma 4.2 in [10],
is bounded (above and below) by
for each provided and are large enough. Hence,
is bounded (above and below) by,