- Email: [email protected]
A note of caution when interpreting parameters of the distribution of excesses
Advances in Water Resources 34 (2011) 1215–1221
Contents lists available at ScienceDirect
Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres
A note of caution when interpreting parameters of the distribution of excesses Pierre Ribereau a,⇑, Philippe Naveau b, Armelle Guillou c a
Université de Montpellier II, I3M, CC051, Place Eugéne Bataillon, 34090 Montpellier, France CNRS, IPSL, Laboratoire des Sciences du Climat et de l’Environnement, Orme des Merisiers/Bat. 701, C.E. Saclay, 91191 Gif-sur-Yvette, France c Université de Strabourg, CNRS, IRMA UMR 7501, 7 rue André Descartes, 67084 Strabourg Cedex, France b
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 15 November 2010 Received in revised form 10 May 2011 Accepted 10 May 2011 Available online 12 June 2011 Keywords: Generalized Pareto Distribution Reparametrization Probability Weighted Moment Maximum Likelihood Extreme Value Theory Generalized Extreme Value Distribution
In climatology and hydrology, univariate Extreme Value Theory has become a powerful tool to model the distribution of extreme events. The Generalized Pareto Distribution (GPD) is routinely applied to model excesses in space or time by letting the two GPD parameters depend on appropriate covariates. Two possible pitfalls of this strategy are the modeling and the interpretation of the scale and shape GPD parameters estimates which are often and incorrectly viewed as independent variables. In this note we first recall a statistical technique that makes the GPD estimates less correlated within a Maximum Likelihood (ML) estimation approach. In a second step we propose novel reparametrizations for two method-ofmoments particularly popular in hydrology: the Probability Weighted Moment (PWM) method and its generalized version (GPWM). Finally these three inference methods (ML, PWM and GPWM) are compared and discussed with respect to the issue of correlations. Ó 2011 Elsevier Ltd. All rights reserved.
1. Classical definition of the parameters of the distribution of excesses During the last decade, there has been a burst of interest in analyzing climatological and hydrological extremes. Such events, by definition, are rare and consequently a probabilistic framework may represent an appropriate foundation to study them. In this context, Extreme Value Theory (EVT) offers a solid mathematical blueprint (e.g. [2,4,11,15]). This theory provides a well defined class of probability density functions to model the distribution of maxima or large values. In the ideal case, the practitioner assumes that a sample of independent and identically distributed (i.i.d.) random variables is available and her/his classical goal is to infer the distributional characteristics of excesses, i.e. values above a very high threshold. According to [8], the survival distribution of such excesses can be approximated, under mild conditions, by the following Generalized Pareto Distribution (GPD) tail
( Gc;r ðxÞ ¼
1=c 1 þ rc x x exp r
r > 0; if c ¼ 0; r > 0;
there are essential in terms of inference and interpretation. Before illustrating this point, we first need to recall a few elements concerning the estimation of the GPD parameters. At least two inferential approaches, Maximum Likelihood (ML) estimation and methods-of-moments, have been frequently applied in hydrology and climatology, (e.g. [18]). Other techniques exist like Bayesian inferences (e.g. [6]) but they are not yet as popular in geosciences and they would not be discussed here. The ML estimation method applied to the GPD function has been investigated by Smith [22] and it consists in maximizing the log-likelihood of the GP distributed sample
‘ðc; rÞ ¼ n log r
⇑ Corresponding author. Tel.: +33 4 67 14 35 78; fax: +33 4 67 14 35 58. E-mail address: [email protected] (P. Ribereau). 0309-1708/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2011.05.003
c
cX i ; log 1 þ
r
i¼1
where (X1, . . . , Xn) corresponds to n i.i.d. excesses. If c > 0.5, then the resulting estimators are asymptotically normally distributed with covariance matrix given by
if c – 0;
where x P 0 when c P 0 and x 2 [0, r/c[ when c < 0. The parameters r and c are called the scale and shape parameters respectively. The names choice and the definitions of these parameters seem rather trivial and unimportant. But, as we will see in this article,
n 1þc X
RML ¼
ð1 þ cÞ2
rð1 þ cÞ
rð1 þ cÞ 2r2 ð1 þ cÞ
! :
It is important to notice that this matrix is not diagonal. This implies that the ML estimators are not independent. This contrasts with the Normal distribution for which the ML estimators of the mean and variance, the sample mean and the sample variance, are statistically independent. This decoupling is fundamental because it allows modeling flexibility and simple parameter interpretations. But the GPD differs from the Gaussian case and this limitation for the ML
1216
P. Ribereau et al. / Advances in Water Resources 34 (2011) 1215–1221
GP estimates can have an important impact. For a different shape parameter, the left panels of Fig. 2 displays the ML estimates of r (x-axis) and c (y-axis) for 500 GP samples of size 100 with a unit scale parameter r = 1 and c = 0.4, 0 and 0.4. Another way to pres^jr ^. ^ and ½r ^ jc ent this issue is to look at the conditional variables ½c ^ Basic probability theory allows us to write that the ML estimate c ^ is approximatively equal in distribution to given r
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
^ jr ^ ’ c ½c
r^ r ð1 þ cÞð1 þ 2cÞ Gaussianð0; 1Þ þ 2k 2r
RPWM ¼
rffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ^ 1 þ 2c ^ ’ r ^ jc Gaussianð0; 1Þ: ½r ðc cÞ þ r 1þc k
parameters can be written as follows
EðXÞ 2EðXGc;r ðXÞÞ
and
r¼
2EðXÞEðXGc;r ðXÞÞ EðXÞ 2EðXGc;r ðXÞÞ
:
By replacing EðXÞ and EðXGc;r ðXÞÞ by their corresponding empirical ^PWM and versions, it can be shown that, if c < 0.5, the estimators c r^ PWM are asymptotically unbiased and normally distributed with covariance matrix
and
1 ð1 2cÞð3 2cÞ 2 2 ð1 cÞð2 cÞ ð1 c2þ 2c 3Þ rðc 2Þð2 6c þ 7c 2c Þ
rðc 2Þð2 6c þ 7c2 2c3 Þ : r2 ð7 18c þ 11c2 2c3 Þ
Ref. [16] also claimed that PWMs estimators perform better than a ML estimation for small samples but the PWM approach has been criticized by Coles and Dixon [5] for imposing the constraint c < 0.5. To relax this restriction in order to model heavier tails, Diebolt et al. [10] replaced EðXÞ by a so-called Generalized Probability Weighted Moment (GPWM) E½XG1:5 c;r ðXÞ. This leads to the following system
c^GPWM ¼
^ 1 ð2:5Þ2 l ^ 1:5 4l ^ 1 2:5l ^ 1:5 2l
and
^ l ^ 2:5l r^ GPWM ¼ ^ 1:5 ^1 ; 2l1 2:5l1:5
ð1Þ
^ s corresponds to the empirical estimator of where l h i ls ¼ E XGsc;r ðXÞ . The asymptotic normality with its associated
0.4 0.6 0.8
correlation
0.2
0.0
covariance matrix is shown in [10] and it holds if c < 1.5. This constraint on the shape parameter covers the vast majority of cases encountered in climatology and hydrology. The immediate advantage of the GPWM over the ML approach resides in the computational simplicity of the GPWM estimators defined by (1). Having explicit expressions, there is no need of an optimization (maximization) scheme to find estimators, such an optimization can produce unreasonable estimates in few cases. Computing moments is also extremely fast and this computational advantage can be decisive when analyzing very large data sets from ensembles of Global Climate Model outputs [17]. In terms of correlations between the shape and the scale parameters, Fig. 1 compares the ML, PWM and GPWM theoretical asymptotic correlations (y-axis) for values of c ranging from 0.5 to 1.5 (x-axis) with r set to 1. These three methods present a strong negative correlation that roughly varies from 0.9 to 0.2. This confirms that this issue cannot be considered as minor. This graph also shows the limitation of the PWM approach, its theoretical correlation being not defined for c P 0.5.
1.0
The first equality tells us that, given a scale estimate that largely differs from r, the shape parameter is then tainted by a conditional bias. Is this issue of correlations between MLE estimates primordial for applications? To answer this question, we recall that one of the most pressing societal question within the climate community is to determine if the observed global temperature warming could lead to higher extreme events in terms of frequency and intensity. In this paper, we would not address this very complex inquiry but this leads to the statistical question on how to model non-stationarity behaviors while modeling excesses with a GPD. A possible answer within the statistical climatology community resides into letting one or two GPD parameters vary in time or space in order to capture temporal or regional changes in excesses intensity. Typically the scale parameter r or its logarithmic is supposed to be a function of some relevant covariates (temporal or spatial) and the shape parameter c remains fixed (at least in time) but may change in space. Changing the c implies that the underlying process in geophysical applications has changed. For example, a shape parameter for temperature extremes recorded over sea can be different the one obtained from extremes measured in mountainous regions because processes related to oceans can be different from the ones that drive temperatures at high altitudes. But at a specific location, it is difficult to find physical mechanisms that could explain a large change in the shape parameter. In addition the shape parameter is difficult to estimate and small variations may be difficult to infer. The issue of correlated MLE estimates exists also when analyzing maxima. For example, Maraun and his co-authors [20,21] analyzed spatio-temporal structures by plotting the MLE scale and shape estimates of a Generalized Extreme Value distribution fitted to UK precipitation maxima. Implicitly, they assumed that these GEV parameters estimates could be interpreted independently. As for the GPD estimates, this assumption can be challenged. As illustrated by Fig. 2, the negative correlation between the two ML GPD parameters is far from negligible, even for positive values of c. This issue is not restricted to the ML framework. The other popular inferential approach used in hydrology is simply based on the method-of-moments and it has the same drawback. To see this, we first need to recall its principle and its elementary properties. Basically, the explicit derivation of two well-chosen GPD moments provides a system of two equations in which the two GPD parameters are considered unknown. Solving these equations is equivalent to find estimators of r and c as function of the two moments. Obviously the choice of the two moments leads to different results and they have to be carefully chosen. For example, the Probability Weighted Moments (PWM) method introduced by Greenwood et al. [14] and Landwehr et al. [19] belongs to this class of method. Hosking and Wallis [16] studied the following two GPD r . In this case the GPD moments EðXÞ ¼ 1r c and E XGc;r ðXÞ ¼ 2ð2 cÞ
EðXÞ
c¼2
0.5
0.0
0.5
1.0
1.5
Fig. 1. Theoretical asymptotic correlations (y-axis) with c varying from 0.5 to 1.5 (x-axis) and r set to 1. The solid, dashed and dotted lines correspond to the ML, PWM and GPWM approaches, respectively.
1217
●
−0.5
shape
●
●● ●● ●● ● ● ●● ● ●● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ●● ● ● ● ●●●● ● ●●● ●● ● ● ● ● ● ● ● ●●● ●● ● ● ●●● ● ●● ●● ●● ● ●● ● ●●● ●● ● ● ●
●
−1.0
−0.5
●● ● ● ●● ● ● ●● ●● ●●●● ● ● ●● ●●●● ●● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ●● ● ● ●●●● ● ● ● ● ● ●● ● ●● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ● ● ●●●● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ●●● ● ● ● ● ●●● ●● ● ● ●● ●●● ● ● ● ● ●● ●●● ● ● ●
−1.5
−1.5
−1.0
shape
0.0
0.0
P. Ribereau et al. / Advances in Water Resources 34 (2011) 1215–1221
0.5
1.0
1.5
0.0
0.5
0.4
●
0.2 0.0 ●
−0.4
●
● ●
−0.8
−0.8
−0.6
−0.6
●
● ● ● ●● ● ●● ● ● ● ● ● ●● ● ●● ● ● ●● ● ● ●● ●● ● ● ●●● ● ● ●● ●● ● ●● ● ● ●● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ●● ●● ●● ● ● ●● ●● ● ●● ● ● ● ●● ● ● ● ● ● ● ●●●● ●●●● ●● ●● ● ● ●● ● ● ● ● ● ● ● ●● ● ●●● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ●● ●● ● ● ●●● ● ● ●● ●●●● ●● ● ● ● ● ● ●● ● ●● ●● ● ●● ● ●● ● ● ● ● ●●● ● ●● ● ●● ● ●● ● ● ● ● ●● ●● ●● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ●● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ●● ●● ● ●● ● ● ●● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ●● ● ●● ● ● ●● ●●●● ● ● ● ●●● ● ●● ● ● ●● ● ●● ● ●●● ● ●● ● ●●● ● ● ●●● ●● ● ● ●● ● ● ●● ● ● ●● ●● ●●● ● ●●●● ●●●●● ●● ● ●● ● ●● ●● ● ●● ● ●●●● ●● ● ● ●● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ●●● ● ● ● ●●● ● ● ●● ● ● ●● ●● ●●●● ●● ● ● ●● ● ● ● ●● ●● ● ●● ● ●● ● ● ●● ● ● ●● ● ● ●
−0.2
shape
0.2 0.0 −0.2
● ●
●
−0.4
shape
●
● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ●●●● ● ● ●● ●●● ●● ●● ● ● ●●●● ● ● ● ●● ● ●●● ●● ● ● ● ● ● ● ●●● ●● ● ●● ● ●●●●● ●●●●● ● ● ●● ● ● ●● ●● ●● ● ●●● ●●● ● ● ● ● ●● ● ● ●●●● ● ● ●● ● ● ●● ● ●● ●● ●● ●● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ●● ●●●● ● ● ● ● ● ● ●● ●●● ● ● ● ●●● ●● ● ● ● ● ●●● ● ●● ●● ● ● ● ● ●● ● ●● ● ● ●● ●●● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ●●●● ● ● ●● ●● ● ●● ●●●● ● ● ●● ●●● ● ● ● ● ● ● ● ● ● ● ● ●●●● ●● ● ●●● ●● ●● ● ● ●● ● ● ●●●● ● ● ●● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ●● ●● ● ● ● ●●● ● ●● ● ● ● ●● ●● ●● ● ● ● ●●● ●● ● ● ●●●● ●● ●● ● ● ●● ● ●● ● ●●● ● ● ●● ● ● ● ● ●● ●● ●● ● ● ● ● ● ●● ● ● ● ●● ●● ● ● ● ●●● ●● ● ●● ● ●● ●● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ●
0.6
0.8
1.0
1.2
1.4
1.6
0.6
0.8
1.0
1.2
1.4
1.6
scale
1.0
1.0
scale
0.8 0.6 0.4
shape
0.6 0.4 0.2
●
●
●
●
●
●
● ●
●
● ●
−0.2
● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ●●● ●●● ●● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ●●● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ●● ● ● ● ● ●● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ●● ● ● ● ●● ●●● ● ● ● ●● ● ● ● ● ● ● ●●●●● ● ●● ●● ● ● ● ● ● ● ● ● ●● ● ●●● ● ● ●● ●●● ●● ● ●●● ●● ● ● ● ● ●● ● ● ●● ●●● ●● ●● ● ● ● ● ● ● ● ●● ●● ● ● ● ●● ●● ●●● ● ●● ● ●●●● ●● ● ● ●● ●●●● ● ●● ● ● ●● ●●● ●●● ● ●●●●● ●● ● ●● ● ●● ● ●●● ● ● ●● ● ● ●● ●● ●● ● ● ● ● ● ● ● ●●● ●●●● ●● ●●● ●● ●● ● ● ● ● ● ●● ● ● ● ●●●●● ●●●● ●● ●● ● ● ●● ●● ●● ●● ●● ● ●● ● ●● ●●● ●●● ●● ●● ● ● ● ● ●● ●● ● ●● ●●● ● ● ●● ● ●●● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ●●● ●● ● ● ●● ●●●●●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ●●● ●●● ●● ● ● ● ● ●● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ●
−0.4
−0.4
−0.2
0.0
●
0.0
0.8
●
● ● ● ● ● ●● ●● ● ● ● ●●● ●●●● ● ● ● ●●● ●● ● ● ● ● ● ●● ●●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ●●● ● ●● ● ● ● ● ● ●● ●●●● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ●● ● ●●●●● ● ● ●● ●● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ●●● ● ● ● ●●● ● ●● ● ● ●● ●● ●● ●● ● ● ● ● ●● ●● ● ● ● ● ●●●● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●●● ●● ● ●●●● ●● ● ● ● ● ●●●● ●● ●● ● ●●● ●● ● ●●● ● ● ●● ● ●● ● ● ●●●● ●● ●●● ●● ● ● ● ● ●● ● ●●● ●● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●● ● ● ●● ●● ●●● ● ●● ●● ●●● ● ● ●● ● ● ● ●● ● ● ● ●●● ●●●● ● ● ●● ●● ● ● ●●● ● ●●● ● ●●● ● ● ● ●● ● ● ● ● ●● ●●● ●● ●● ●●●● ●●●●● ● ● ●● ●● ● ● ●●●● ●● ●● ●● ● ● ●● ● ● ●● ●● ● ● ●● ● ● ● ● ●● ●●● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ●● ●●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●
0.2
●
● ●
shape
1.0
scale
0.4
scale
0.4
0.6
0.8
1.0
1.2
1.4
1.6
scale
0.8
1.0
1.2
1.4
1.6
1.8
2.0
scale
Fig. 2. ML estimates for 500 samples of size 100 from a GPD (c, r = 1) with c = 0.4 (upper), 0 (middle) and 0.4 (lower). The left panels correspond to the old ML parametrization and the right panel to the new one defined by (2). The intersection of the two solid lines displays the true value of both parameters.
The left panels of Fig. 3, the analogue of the ones of Fig. 2 but for GPWMs, corroborates the fact that GPWM estimates are also dependent.
2. Alternative definitions of the GPD parameters To solve the ML estimators issue raised in the previous section, Chavez-Demoulin and Davison [3] proposed an orthogonalization
1218
P. Ribereau et al. / Advances in Water Resources 34 (2011) 1215–1221
●
−0.5
shape
0.0
● ● ●● ● ● ● ●●●● ● ● ● ●● ● ●● ●● ●●● ● ●● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ●● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ●●● ● ● ● ● ● ●● ●● ● ● ● ● ●● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●●● ●●● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ●● ●●● ● ● ●● ● ●● ●● ● ● ●● ● ●● ● ● ● ●● ● ● ●●● ● ● ● ●●● ●●● ●● ● ● ● ● ● ●● ●●
−1.0
0.0 −0.5 −1.0
shape
●
● ● ● ●●● ● ●● ● ● ● ●● ● ● ● ●● ● ●●● ●● ● ●●● ● ● ● ● ● ●● ●●● ●● ● ● ● ● ●●●● ● ● ● ●● ●● ●● ●●●● ● ●●● ●● ● ● ● ● ● ● ●● ●● ● ● ●● ● ● ●● ● ● ●● ● ● ● ●●● ● ● ● ● ● ● ●● ● ● ● ● ●● ●● ● ● ● ●● ● ● ●● ● ● ● ●● ● ●● ● ●●●● ● ●● ●● ● ● ●● ●●●● ● ●● ● ● ●●●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ●●● ● ● ● ● ● ●● ● ● ● ●● ●●●● ● ● ● ● ● ● ●●● ●● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ●●● ● ● ●● ●● ● ● ● ● ●● ●● ●● ●● ● ● ● ● ●●● ● ● ● ● ●● ● ● ●● ● ●●● ● ● ● ● ●● ● ●●● ● ● ● ● ●●● ● ● ● ● ●●● ● ● ● ● ● ●●● ●● ● ●●●● ●● ● ●● ● ● ● ● ● ●●●● ●● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ●●● ● ●● ● ● ● ●●● ● ● ●●● ● ● ●●● ● ● ●● ● ● ●●● ● ● ● ●● ●● ● ●● ● ●●● ● ●● ● ● ●●● ● ● ● ● ● ● ● ●●● ●●● ● ●● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●
● ●
●
● ●
●●
−1.5
−1.5
●
0.5
1.0
1.5
0.0
0.5
scale
−0.6
●
● ● ● ●
0.2 0.0 −0.2
shape ●
●
● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ●● ●● ● ●● ● ● ●● ●●●●●● ● ● ● ●● ● ● ● ● ●●●● ●●●● ● ● ● ● ●● ●● ● ● ● ●● ● ●●● ●●● ●● ●● ● ● ●● ● ● ●● ● ● ● ● ● ●● ●● ●● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ●●●● ●● ●●● ●●● ● ● ● ●●● ●●● ●●●● ●● ● ● ●● ● ●●● ● ● ● ●● ●● ●● ●●● ● ● ●● ● ● ●●● ● ●● ● ●● ● ●● ●● ● ●● ●● ● ●● ● ●● ● ● ●● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●●● ●●● ● ● ● ● ●● ●● ● ●●● ● ● ● ● ● ● ●●● ● ●● ●● ●● ● ● ● ● ●●● ● ●● ● ● ● ● ● ● ● ●●●● ●● ● ● ● ●● ●● ●● ●●●● ● ● ●● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ●●● ●● ● ● ●● ●● ●●●●●● ● ●●● ● ● ●● ● ●●● ● ●●● ●● ●●● ● ● ● ● ● ●●●●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ●●●● ● ●● ● ●● ● ● ● ●● ● ●● ● ● ● ●●●● ●●●● ● ●
0.4
●
−0.4
−0.4
−0.2
0.0
0.2
●
shape
● ●● ● ● ● ● ●● ● ● ● ● ●● ● ● ●● ● ●● ● ● ● ● ●●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ●●● ● ●● ●● ●●● ●● ● ● ● ● ● ● ● ●●● ● ●● ● ●● ● ● ● ●● ●● ●● ●●●● ● ●●● ●●● ●● ● ● ●●● ● ●● ●● ●● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ●●●● ●● ● ●●●● ● ● ● ●● ●● ● ● ●●● ●● ● ● ● ●● ● ● ●●● ● ●● ● ● ● ●●●● ● ●● ●● ● ●● ● ●● ●● ●●● ● ●● ● ● ● ● ● ●●●●● ●●● ● ● ●●● ●●● ●● ●● ●●● ●● ●● ● ● ● ● ●●● ● ●● ●● ●● ● ● ●●● ● ● ●●● ● ●●● ● ●●● ● ● ●● ●● ● ●● ●● ●● ● ● ● ● ●●●● ●●● ● ● ● ●● ● ● ● ● ● ● ●● ●● ● ● ● ● ●●● ● ● ● ●● ● ● ●● ● ● ● ●●● ● ● ● ● ● ● ● ●●● ● ● ●● ●●● ●● ● ● ● ● ●● ● ● ●● ● ● ● ● ●● ● ●● ● ● ● ●● ● ● ● ●● ● ●● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ●● ● ●● ●● ● ●● ● ● ●●● ● ●● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ●
●
● ●●●
●
−0.6
0.4
●
●
●
●
●
−0.8
−0.8
●
0.6
0.8
1.0
1.2
1.4
1.6
0.0
0.2
0.4
0.6
0.8
1.0
scale
1.0
1.0
scale
−0.2
0.8
●
●●● ● ●● ● ● ● ● ● ●● ●● ●●●● ● ● ● ●● ● ● ●● ● ●● ● ●● ● ●●●● ● ● ● ●● ● ● ● ● ● ●● ● ●● ● ●● ● ● ● ● ● ● ● ●●● ●●● ●●●● ●●● ●● ● ● ● ● ● ●● ● ●●●●●● ●● ●● ●● ● ● ● ● ●● ● ●● ● ● ● ●● ●● ● ●● ●●● ●●● ● ●●●● ●● ●● ●● ● ●● ● ●● ● ●● ● ●● ●●● ●● ● ●● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ●● ●● ●● ● ●● ● ● ● ●● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●●● ●● ● ●● ●● ● ● ● ● ● ●●● ●● ●● ● ●●●●●●●● ● ●● ● ● ● ● ● ● ● ● ●●● ● ● ●● ●● ● ● ●●●● ● ●● ●● ●● ●●● ● ●● ●● ●● ●●● ● ●● ● ● ● ●● ●● ●●● ● ● ● ●●● ●● ● ● ●● ●● ● ● ● ● ●●● ●● ●●●● ● ● ●● ●● ● ●● ● ●● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●● ●●● ●● ●● ●● ● ● ● ● ●● ●● ●● ● ● ●● ●● ●● ● ●●● ● ●●● ● ●●● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ●
0.4
shape ●
●
−0.4
−0.4
●
● ●
● ●
●
0.2
●● ● ● ●● ● ● ● ● ● ●● ●● ● ●●● ● ● ● ● ● ● ● ● ●●● ● ● ●● ● ●● ● ●● ● ● ●● ●● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ●●●● ●●● ● ● ●● ● ● ● ● ● ●● ●● ● ● ●● ●●● ●● ●●● ●● ● ● ●●●● ●● ●● ● ● ● ● ● ● ● ●● ● ●● ●●● ● ●● ● ● ● ● ●● ●● ● ● ● ●● ● ●● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●●●● ● ●● ●●● ●● ●● ●●●●●●● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ●● ● ● ●● ●● ● ● ●● ●●● ● ● ●●●●●● ● ● ● ● ● ●● ●● ● ● ● ● ●●●● ● ● ●●●●● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ●●●● ● ●● ●● ● ● ●●● ● ● ● ● ●● ● ●● ● ●●● ●●● ● ●●●● ●● ● ●●● ●● ● ●●● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ●●●● ● ● ●● ●● ● ● ●● ●●● ● ●● ● ● ● ●● ● ●● ● ● ● ● ● ●●● ● ● ● ● ● ●● ● ● ●● ●● ● ●● ● ●● ● ● ● ● ●● ●● ● ● ●● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ●●●● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ●
0.6
●
0.0
0.0
0.2
0.4
0.6
●
● ●
●
−0.2
0.8
●
shape
1.0
scale
0.4
0.6
0.8
1.0
1.2
1.4
1.6
0.2
scale
0.4
0.6
0.8
1.0
1.2
1.4
scale
Fig. 3. GPWM estimates for 500 samples of size 100 from a GPD (c, r = 1) with c = 0.4 (upper), 0 (middle) and 0.4 (lower). The left panels correspond to the old GPWM parametrization and the right panel to the new one defined by (4). The intersection of the two solid lines displays the true value of both parameters.
technique, see, e.g. [7], to obtain a diagonal Fisher information matrix. This leads to the following reparametrization
cI;ML :¼ c; rI;ML :¼ rð1 þ cÞ:
ð2Þ
The two advantages of this parametrization are that the shape parameter c remains unchanged and that the scale parameter r is only modified by a multiplicative factor (1 + c) which does not contain r. This latter remark implies that the GPD scaling property remains valid with the new parametrization, i.e. if X follows a GPD
1219
P. Ribereau et al. / Advances in Water Resources 34 (2011) 1215–1221
with scale parameter rI,ML then, for any positive scalar a, the variable aX also follows a GPD with scale parameter a rI,ML. Hence the interpretation of (2) is not fully identical to the old one but it is very comparable. The important bonus is that system (2) allows for individual interpretation and modeling of each GPD parameter. Fig. 2 compares the old parametrization (left panels) and the new one (right panels). It is clear that the dependence between ML estimates has been reduced for the three values of c = 0.4, 0 and 0.4, see the upper, middle and lower panels, respectively. But even with the new parametrization, one drawback of the ML estimation is still visible because of the constraint c > 0.5. Near 0.5, the ML procedure is not ideal and the upper right panel clearly displays a type of non-linear dependence for c = 0.4, see also the first column of Table 1. One could argue that negative shape parameters are unlikely to occur in climatology but recent work by Friederichs [13] hints that heavy wind gusts in Germany can reach negative shape parameters around 0.3, see Fig. 7 of [13]. Also for UK precipitation, negative shape parameters have been found along the west coast (see [20,21]). Such bounded tails can also appear while modeling relative humidity extremes (see [12]). With respect to the method-of-moments, there are currently not adequate reparametrizations for neither the PWM nor the GPWM approach. The next two paragraphs cover both cases. The matrix diagonalization detailed in Appendix A allows us to find the following new parametrization for the PWM case
cI;PWM :¼ c; rI;PWM :¼ raðcÞ;
ð3Þ
where the function a(c) is defined as follows
1 11 2 logð1 cÞ logð1 c þ 2c2 Þ þ logð2 cÞ 2 28 7 pffiffiffi 7 1 þ 4c pffiffiffi þ arctan 14 7
log aðcÞ ¼
for c < 1. Although not trivial, this expression is explicit. As for the ML parametrization captured by (2), the shape parameter is unchanged and the GPD scaling property is still valid. To make another link with (2), we can notice that the function a(c) via a Taylor expansion in the neighborhood of zero can be approximated by a0 (1 + c) with a0 = 1.138544. Concerning the GPWM approach, another matrix diagonalization that can also be found in Appendix B provides the following system
cI;GPWM :¼ c; rI;GPWM :¼ rbðcÞ;
ð4Þ
where the function b(c) is defined as follows
16 8 logð5 2cÞ logð13 10c þ 8c2 Þ 19 19 pffiffiffiffiffiffi 64 pffiffiffiffiffiffi 1 79 arctan ð10 þ 16cÞ 79 þ 1501 158
log bðcÞ ¼ logð2 cÞ þ
for c < 2. Again this expression is explicit. It preserves the shape parameter and the GPD scaling property. A Taylor expansion of
Table 2 Coverage probability of the confidence interval from 1000 GPD samples of size 100 for the different parameters and for different values of the shape parameter c and r = 1. Parameter/ Value of c
0.4
0.2
0
0.2
0.4
0.8
1.2
c with ML r with ML rI,ML from (2)
90.6 94.7 98.0
82.3 91.7 92.4
86.8 93.5 96.8
90.7 93.2 97.8
90.4 92.7 98.2
93.3 94.6 98.7
93.0a 94.1a 97.9a
c with GPWM r with GPWM rI,GPWM with (4)
96.2 96.0 94.1
96.2 96.4 95.0
96.0 95.8 93.3
96.6 96.2 94.6
94.9 94.9 93.9
97.2 98.7 97.1
95.1 94.7 95.9
a After removing all the samples in which the optimization algorithm did not converge or leaded to aberrant estimations.
order 2 of the function b(c) leads to the following approximation: 2 b0 þ b1 c þ b2 c2 with b0 = 0.5422619, b1 = 0.4046108 and b2 = 0.3265764. Hence the system (4) represents a greater departure in terms of parametrization from the ML case than the GPWM one. Fig. 3 can be compared to Fig. 2. Both ML and GPWM appears to work well in terms of decorrelation with a slight advantage for the GPWM, see also Table 1. Table 2 represents the coverage probabilities of the confidence intervals obtained for c and r (before and after the reparametrization) with level 95% for different values of the shape parameters. These results are based on the asymptotic covariance matrices given in Appendix C. As observed in this table, all the coverage probabilities are close to the nominal coverage probability, except for the ML method and c = 0.4 which is close to the value 0.5 required for the asymptotic normality. In the case where c < 0.5, the method is super-efficient (see [22]) and it is still possible to construct confidence interval by bootstrap. However, for all the values of c the GPWM method leads to better results than the ML one since the intervals in that case are conservative whatever the value of c is. 3. Conclusion This short article discusses correlation issues between the two GPD parameters estimators and the risk of interpreting such estimators. Such a problem has rarely been a point of interest within hydrological and climatological circles but, despite its simplicity, it could become more prevalent as researchers have the tendency to model spatial and temporal structures at the GPD parameters level. To decorrelate these parameters, classical matrix diagonalization techniques have been applied here and comparisons have been implemented among three inferences techniques (ML, PWM and GPWM). Overall our new GPWM parametrization appears to be the most straightforward and very fast approach in terms of its implementation (no optimization needed to maximize a likelihood). It offers an adequate decorrelation scheme and can be applied for a wide range of shape parameters. If the practitioner only focuses on the estimation of return levels in the i.i.d. case, then the proposed parametrization may not be needed because return levels combine uncertainties from the shape and scale parameters estimates. But in terms of modeling with covariates, we
Table 1 Estimated correlations obtained from 1000 GPD samples of size 100 for different values of the shape parameter c and r = 1.
a
Inference on c=
0.4
0.2
0
0.2
0.4
0.8
1.2
Classical ML ML with (2)
0.789 0.039
0.822 0.147
0.738 0.098
0.675 0.029
0.644 0.025
0.521 0.003
0.255a 0.154a
Classical GPWM GPWM with (4)
0.759 0.019
0.861 0.031
0.803 0.054
0.759 0.009
0.691 0.013
0.532 0.027
0.423 0.226
After removing all the samples in which the optimization algorithm did not converge or leaded to aberrant estimations.
1220
P. Ribereau et al. / Advances in Water Resources 34 (2011) 1215–1221
believe that decoupling the GPD parameters could facilitate the inference and the interpretation in many applications. In this paper we did not treat the choice of the threshold that define excesses. The influence of such a threshold in terms on parameters estimation and correlation must be non-negligible. This is a delicate statistical problem because the threshold is not viewed within the EVT paradigm as a parameter to be estimated but rather as a constant that has to be chosen a priori (see, e.g. [1]). Finally, the correlation issue investigated in this work is not circumvented to the GPD or to the three estimation methods studied here. For example, a powerful alternative for visualizing standard errors of the GPD estimates resides in computing log likelihood profiles for a given GPD parametrization, see Section 2.6.6 of Coles book. This approach has the advantage of producing non-symmetrical confidence intervals. However it is not clear to us what is the optimal parametrization for this method. Also, the posterior of a classical GPD Bayesian analysis should suffer from the same problem. In terms of distribution, strong correlations also exist for the Generalized Extreme Value (GEV) distribution parameters used to model maxima. Therefore it would be of interest to find new parametrizations for the GEV.
Appendix B. Generalized Probability Weighted Moment estimators
Acknowledgments
This justifies the approximatively independent reparametrization:
Part of this work has been supported by the EU-FP7 ACQWA Project (www.acqwa.ch) under Contract No. 212250, by the PEPER-GIS project, by the ANR-MOPERA project, by the ANR-McSim project and by the MIRACCLE-GICC project. The authors also would like to thanks the referees for valuable comments that have greatly improved the paper.
The GPWM method is explained in [10] for the GPD and in [9] for the GEV distribution. Using the expression of the asymptotic variance given in [10], we get:
0 5 b ðcÞ 5 13 c bðcÞ c 5c þ 4c2 ð2 cÞ 2 bðcÞ 2 2 97 þ 39c 23c2 þ 4c3 ; 4
rð2 cÞ
or similarly
97 þ 39c 23c2 þ 4c3 b0 ðcÞ 4 : ¼ bðcÞ ðc 2Þ 52 c 13 5c þ 4c2 2 Again, by integration, we derive
16 8 logð5 2cÞ logð13 10c þ 8c2 Þ 19 19 pffiffiffiffiffiffi 64 pffiffiffiffiffiffi 1 79 arctan ð10 þ 16cÞ 79 : þ 1501 158
log bðcÞ ¼ logð2 cÞ þ
8 < cI;GPWM ¼ c;
2 :r ~ I;GPWM ¼ r b0 þ b1 c þ b2 c2
ðB:1Þ
with b0 = 0.5422619, b1 = 0.4046108 and b2 = 0.3265764. Appendix C. Asymptotic covariance matrix
Appendix A. Probability Weighted Moments estimators
C.1. Maximum Likelihood
Using the asymptotic covariance matrix RPWM, we want to find the function a() such that the following matrix is diagonal:
Using the reparametrization (2), the estimators obtained by the Maximum Likelihood method are asymptotically normally distributed with covariance matrix given by:
1
0
ra0 ðcÞ aðcÞ
RPWM
ra0 ðcÞ ; 0 aðcÞ
1
!
ð1 þ cÞ2
0
0
2r2 ð1 þ cÞ2 ð1 þ 2cÞ
:
where a0 ¼ @@ac. This leads to the following equation
C.2. Probability Weighted Moments
a0 ðcÞ ð1 cÞðc 2Þð1 c þ 2c2 Þ aðcÞ þ2 6c þ 7c2 2c3 ¼ 0:
rðc 2ÞaðcÞ
Using the reparametrization (3), the corresponding estimators are asymptotically normally distributed with covariance matrix given by:
0
This equation can be rewritten as
@
a0 ðcÞ 2 6c þ 7c2 2c3 ¼ : aðcÞ ð1 cÞð2 cÞð1 c þ 2c2 Þ
This latter leads to the following approximatively independent reparametrization
with a0 = 1.138544.
0
r2 a2 ðcÞ ð1cÞð1cþ2c2 Þ
A:
With the classical GPD parametrization, the asymptotic covariance matrix of the estimators obtained by the GPWM method is the following:
1 11 2 logð1 cÞ logð1 c þ 2c2 Þ þ logð2 cÞ 2 28 7 pffiffiffi 7 1 þ 4c pffiffiffi : þ arctan 14 7
log aðcÞ ¼
cI;PWM ¼ c; r~ I;PWM :¼ a0 ð1 þ cÞ
0
C.3. Generalized Probability Weighted Moments
Finally, by integration we deduce that
1
ð1cÞð2cÞ2 ð1cþ2c2 Þ ð12cÞð32cÞ
ðA:1Þ
Cð1;1Þ rCð1;2Þ rCð1;2Þ r2 Cð2;2Þ
!
with the explicit expression of C(i,j) given in [10]. Now using the parametrization (4), the new asymptotic matrix is given by:
0 @
Cð1;1Þ 0
2
0
1 2 A: ð1;2Þ
ðcÞ r2 Cb ð1;1Þ Cð1;1Þ Cð2;2Þ C
P. Ribereau et al. / Advances in Water Resources 34 (2011) 1215–1221
References [1] Beirlant J, Dierckx G, Goegebeur Y, Matthys G. Tail index estimation and an exponential regression model. Extremes 1999;2:177–200. [2] Beirlant J, Goegebeur Y, Segers J, Teugels J. Statistics of extremes: theory and applications. Wiley series in probability and statistics, 2004. [3] Chavez-Demoulin V, Davison AC. Generalized additive modelling of sample extremes. Appl Stat 2005;54:207–22. [4] Coles SG. An introduction to statistical modeling of extreme values. Springer; 2001. [5] Coles SG, Dixon MJ. Likelihood-based inference of extreme value models. Extremes 1999;2:5–23. [6] Cooley D, Nychka D, Naveau P. Bayesian spatial modeling of extreme precipitation return levels. J Am Stat Assoc 2007;102:824–40. [7] Cox DR, Reid N. Parameter orthogonality and approximate conditional inference (with discussion). J Roy Stat Soc B 1987;49:1–39. [8] Davison A, Smith R. Models for exceedances over high thresholds (with discussion). J Roy Stat Soc B 1990;52:393–442. [9] Diebolt J, Guillou A, Naveau P, Ribereau P. Improving probability-weighted moment methods for the generalized extreme value distribution. REVSTAT Stat J 2008;6:33–50. [10] Diebolt J, Guillou A, Rached I. Approximation of the distribution of excesses through a generalized probability-weighted moments method. J Statist Plann Inf 2007;137:841–57. [11] Embrechts P, Klüppelberg C, Mikosch T. Modelling extremal events for insurance and finance. Applications of Mathematics, vol. 33. Berlin: SpringerVerlag; 1997.
1221
[12] Flecher C, Allard D, Naveau P. Truncated skew-normal distributions: moments, estimation by weighted moments and application to climatic data. METRON – Int J Stat 2010;LXVIII(3):265–79. [13] Friederichs P, Gober M, Bentzien S, Lenz A, Krampit R. A probabilistic analysis of wind gusts using extreme value statistics. Meteorol Z 2009;18: 615–29. [14] Greenwood JA, Landwehr JM, Matalas NC, Wallis JR. Probability-weighted moments: definition and relation to parameters of several distributions expressable in inverse form. Water Resour Res 1979;15:1049–54. [15] de Haan L, Ferreira A. Extreme value theory: an introduction. Springer series in operations research, 2006. [16] Hosking JRM, Wallis JR. Parameter and quantile estimation for the generalized Pareto distribution. Technometrics 1987;29:339–49. [17] IPCC Climate Change: Synthesis Report. http://www.ipcc.ch/ipccreports/, 2007. [18] Katz R, Parlange M, Naveau P. Extremes in hydrology. Adv Water Resour 2002;25:1287–304. [19] Landwehr J, Matalas N, Wallis JR. Probability weighted moments compared with some traditional techniques in estimating Gumbel parameters and quantiles. Water Resour Res 1979;15:1055–64. [20] Maraun D, Rust HW, Osborn TJ. The annual cycle of heavy precipitation accross the UK: a model based on extreme value statistics. Int J Climatol 2009;29: 1731–44. [21] Maraun D, Rust HW, Osborn TJ. The influence of synoptic airflow on UK daily precipitation extremes. Part I: observed spatio-temporal relations. Clim Dyn 2010;36:261–75. [22] Smith RL. Maximum likelihood estimation in a class of nonregular cases. Biometrika 1985;72:69–90.
Contents lists available at ScienceDirect
Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres
A note of caution when interpreting parameters of the distribution of excesses Pierre Ribereau a,⇑, Philippe Naveau b, Armelle Guillou c a
Université de Montpellier II, I3M, CC051, Place Eugéne Bataillon, 34090 Montpellier, France CNRS, IPSL, Laboratoire des Sciences du Climat et de l’Environnement, Orme des Merisiers/Bat. 701, C.E. Saclay, 91191 Gif-sur-Yvette, France c Université de Strabourg, CNRS, IRMA UMR 7501, 7 rue André Descartes, 67084 Strabourg Cedex, France b
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 15 November 2010 Received in revised form 10 May 2011 Accepted 10 May 2011 Available online 12 June 2011 Keywords: Generalized Pareto Distribution Reparametrization Probability Weighted Moment Maximum Likelihood Extreme Value Theory Generalized Extreme Value Distribution
In climatology and hydrology, univariate Extreme Value Theory has become a powerful tool to model the distribution of extreme events. The Generalized Pareto Distribution (GPD) is routinely applied to model excesses in space or time by letting the two GPD parameters depend on appropriate covariates. Two possible pitfalls of this strategy are the modeling and the interpretation of the scale and shape GPD parameters estimates which are often and incorrectly viewed as independent variables. In this note we first recall a statistical technique that makes the GPD estimates less correlated within a Maximum Likelihood (ML) estimation approach. In a second step we propose novel reparametrizations for two method-ofmoments particularly popular in hydrology: the Probability Weighted Moment (PWM) method and its generalized version (GPWM). Finally these three inference methods (ML, PWM and GPWM) are compared and discussed with respect to the issue of correlations. Ó 2011 Elsevier Ltd. All rights reserved.
1. Classical definition of the parameters of the distribution of excesses During the last decade, there has been a burst of interest in analyzing climatological and hydrological extremes. Such events, by definition, are rare and consequently a probabilistic framework may represent an appropriate foundation to study them. In this context, Extreme Value Theory (EVT) offers a solid mathematical blueprint (e.g. [2,4,11,15]). This theory provides a well defined class of probability density functions to model the distribution of maxima or large values. In the ideal case, the practitioner assumes that a sample of independent and identically distributed (i.i.d.) random variables is available and her/his classical goal is to infer the distributional characteristics of excesses, i.e. values above a very high threshold. According to [8], the survival distribution of such excesses can be approximated, under mild conditions, by the following Generalized Pareto Distribution (GPD) tail
( Gc;r ðxÞ ¼
1=c 1 þ rc x x exp r
r > 0; if c ¼ 0; r > 0;
there are essential in terms of inference and interpretation. Before illustrating this point, we first need to recall a few elements concerning the estimation of the GPD parameters. At least two inferential approaches, Maximum Likelihood (ML) estimation and methods-of-moments, have been frequently applied in hydrology and climatology, (e.g. [18]). Other techniques exist like Bayesian inferences (e.g. [6]) but they are not yet as popular in geosciences and they would not be discussed here. The ML estimation method applied to the GPD function has been investigated by Smith [22] and it consists in maximizing the log-likelihood of the GP distributed sample
‘ðc; rÞ ¼ n log r
⇑ Corresponding author. Tel.: +33 4 67 14 35 78; fax: +33 4 67 14 35 58. E-mail address: [email protected] (P. Ribereau). 0309-1708/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2011.05.003
c
cX i ; log 1 þ
r
i¼1
where (X1, . . . , Xn) corresponds to n i.i.d. excesses. If c > 0.5, then the resulting estimators are asymptotically normally distributed with covariance matrix given by
if c – 0;
where x P 0 when c P 0 and x 2 [0, r/c[ when c < 0. The parameters r and c are called the scale and shape parameters respectively. The names choice and the definitions of these parameters seem rather trivial and unimportant. But, as we will see in this article,
n 1þc X
RML ¼
ð1 þ cÞ2
rð1 þ cÞ
rð1 þ cÞ 2r2 ð1 þ cÞ
! :
It is important to notice that this matrix is not diagonal. This implies that the ML estimators are not independent. This contrasts with the Normal distribution for which the ML estimators of the mean and variance, the sample mean and the sample variance, are statistically independent. This decoupling is fundamental because it allows modeling flexibility and simple parameter interpretations. But the GPD differs from the Gaussian case and this limitation for the ML
1216
P. Ribereau et al. / Advances in Water Resources 34 (2011) 1215–1221
GP estimates can have an important impact. For a different shape parameter, the left panels of Fig. 2 displays the ML estimates of r (x-axis) and c (y-axis) for 500 GP samples of size 100 with a unit scale parameter r = 1 and c = 0.4, 0 and 0.4. Another way to pres^jr ^. ^ and ½r ^ jc ent this issue is to look at the conditional variables ½c ^ Basic probability theory allows us to write that the ML estimate c ^ is approximatively equal in distribution to given r
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
^ jr ^ ’ c ½c
r^ r ð1 þ cÞð1 þ 2cÞ Gaussianð0; 1Þ þ 2k 2r
RPWM ¼
rffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ^ 1 þ 2c ^ ’ r ^ jc Gaussianð0; 1Þ: ½r ðc cÞ þ r 1þc k
parameters can be written as follows
EðXÞ 2EðXGc;r ðXÞÞ
and
r¼
2EðXÞEðXGc;r ðXÞÞ EðXÞ 2EðXGc;r ðXÞÞ
:
By replacing EðXÞ and EðXGc;r ðXÞÞ by their corresponding empirical ^PWM and versions, it can be shown that, if c < 0.5, the estimators c r^ PWM are asymptotically unbiased and normally distributed with covariance matrix
and
1 ð1 2cÞð3 2cÞ 2 2 ð1 cÞð2 cÞ ð1 c2þ 2c 3Þ rðc 2Þð2 6c þ 7c 2c Þ
rðc 2Þð2 6c þ 7c2 2c3 Þ : r2 ð7 18c þ 11c2 2c3 Þ
Ref. [16] also claimed that PWMs estimators perform better than a ML estimation for small samples but the PWM approach has been criticized by Coles and Dixon [5] for imposing the constraint c < 0.5. To relax this restriction in order to model heavier tails, Diebolt et al. [10] replaced EðXÞ by a so-called Generalized Probability Weighted Moment (GPWM) E½XG1:5 c;r ðXÞ. This leads to the following system
c^GPWM ¼
^ 1 ð2:5Þ2 l ^ 1:5 4l ^ 1 2:5l ^ 1:5 2l
and
^ l ^ 2:5l r^ GPWM ¼ ^ 1:5 ^1 ; 2l1 2:5l1:5
ð1Þ
^ s corresponds to the empirical estimator of where l h i ls ¼ E XGsc;r ðXÞ . The asymptotic normality with its associated
0.4 0.6 0.8
correlation
0.2
0.0
covariance matrix is shown in [10] and it holds if c < 1.5. This constraint on the shape parameter covers the vast majority of cases encountered in climatology and hydrology. The immediate advantage of the GPWM over the ML approach resides in the computational simplicity of the GPWM estimators defined by (1). Having explicit expressions, there is no need of an optimization (maximization) scheme to find estimators, such an optimization can produce unreasonable estimates in few cases. Computing moments is also extremely fast and this computational advantage can be decisive when analyzing very large data sets from ensembles of Global Climate Model outputs [17]. In terms of correlations between the shape and the scale parameters, Fig. 1 compares the ML, PWM and GPWM theoretical asymptotic correlations (y-axis) for values of c ranging from 0.5 to 1.5 (x-axis) with r set to 1. These three methods present a strong negative correlation that roughly varies from 0.9 to 0.2. This confirms that this issue cannot be considered as minor. This graph also shows the limitation of the PWM approach, its theoretical correlation being not defined for c P 0.5.
1.0
The first equality tells us that, given a scale estimate that largely differs from r, the shape parameter is then tainted by a conditional bias. Is this issue of correlations between MLE estimates primordial for applications? To answer this question, we recall that one of the most pressing societal question within the climate community is to determine if the observed global temperature warming could lead to higher extreme events in terms of frequency and intensity. In this paper, we would not address this very complex inquiry but this leads to the statistical question on how to model non-stationarity behaviors while modeling excesses with a GPD. A possible answer within the statistical climatology community resides into letting one or two GPD parameters vary in time or space in order to capture temporal or regional changes in excesses intensity. Typically the scale parameter r or its logarithmic is supposed to be a function of some relevant covariates (temporal or spatial) and the shape parameter c remains fixed (at least in time) but may change in space. Changing the c implies that the underlying process in geophysical applications has changed. For example, a shape parameter for temperature extremes recorded over sea can be different the one obtained from extremes measured in mountainous regions because processes related to oceans can be different from the ones that drive temperatures at high altitudes. But at a specific location, it is difficult to find physical mechanisms that could explain a large change in the shape parameter. In addition the shape parameter is difficult to estimate and small variations may be difficult to infer. The issue of correlated MLE estimates exists also when analyzing maxima. For example, Maraun and his co-authors [20,21] analyzed spatio-temporal structures by plotting the MLE scale and shape estimates of a Generalized Extreme Value distribution fitted to UK precipitation maxima. Implicitly, they assumed that these GEV parameters estimates could be interpreted independently. As for the GPD estimates, this assumption can be challenged. As illustrated by Fig. 2, the negative correlation between the two ML GPD parameters is far from negligible, even for positive values of c. This issue is not restricted to the ML framework. The other popular inferential approach used in hydrology is simply based on the method-of-moments and it has the same drawback. To see this, we first need to recall its principle and its elementary properties. Basically, the explicit derivation of two well-chosen GPD moments provides a system of two equations in which the two GPD parameters are considered unknown. Solving these equations is equivalent to find estimators of r and c as function of the two moments. Obviously the choice of the two moments leads to different results and they have to be carefully chosen. For example, the Probability Weighted Moments (PWM) method introduced by Greenwood et al. [14] and Landwehr et al. [19] belongs to this class of method. Hosking and Wallis [16] studied the following two GPD r . In this case the GPD moments EðXÞ ¼ 1r c and E XGc;r ðXÞ ¼ 2ð2 cÞ
EðXÞ
c¼2
0.5
0.0
0.5
1.0
1.5
Fig. 1. Theoretical asymptotic correlations (y-axis) with c varying from 0.5 to 1.5 (x-axis) and r set to 1. The solid, dashed and dotted lines correspond to the ML, PWM and GPWM approaches, respectively.
1217
●
−0.5
shape
●
●● ●● ●● ● ● ●● ● ●● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ●● ● ● ● ●●●● ● ●●● ●● ● ● ● ● ● ● ● ●●● ●● ● ● ●●● ● ●● ●● ●● ● ●● ● ●●● ●● ● ● ●
●
−1.0
−0.5
●● ● ● ●● ● ● ●● ●● ●●●● ● ● ●● ●●●● ●● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ●● ● ● ●●●● ● ● ● ● ● ●● ● ●● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ● ● ●●●● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ●●● ● ● ● ● ●●● ●● ● ● ●● ●●● ● ● ● ● ●● ●●● ● ● ●
−1.5
−1.5
−1.0
shape
0.0
0.0
P. Ribereau et al. / Advances in Water Resources 34 (2011) 1215–1221
0.5
1.0
1.5
0.0
0.5
0.4
●
0.2 0.0 ●
−0.4
●
● ●
−0.8
−0.8
−0.6
−0.6
●
● ● ● ●● ● ●● ● ● ● ● ● ●● ● ●● ● ● ●● ● ● ●● ●● ● ● ●●● ● ● ●● ●● ● ●● ● ● ●● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ●● ●● ●● ● ● ●● ●● ● ●● ● ● ● ●● ● ● ● ● ● ● ●●●● ●●●● ●● ●● ● ● ●● ● ● ● ● ● ● ● ●● ● ●●● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ●● ●● ● ● ●●● ● ● ●● ●●●● ●● ● ● ● ● ● ●● ● ●● ●● ● ●● ● ●● ● ● ● ● ●●● ● ●● ● ●● ● ●● ● ● ● ● ●● ●● ●● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ●● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ●● ●● ● ●● ● ● ●● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ●● ● ●● ● ● ●● ●●●● ● ● ● ●●● ● ●● ● ● ●● ● ●● ● ●●● ● ●● ● ●●● ● ● ●●● ●● ● ● ●● ● ● ●● ● ● ●● ●● ●●● ● ●●●● ●●●●● ●● ● ●● ● ●● ●● ● ●● ● ●●●● ●● ● ● ●● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ●●● ● ● ● ●●● ● ● ●● ● ● ●● ●● ●●●● ●● ● ● ●● ● ● ● ●● ●● ● ●● ● ●● ● ● ●● ● ● ●● ● ● ●
−0.2
shape
0.2 0.0 −0.2
● ●
●
−0.4
shape
●
● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ●●●● ● ● ●● ●●● ●● ●● ● ● ●●●● ● ● ● ●● ● ●●● ●● ● ● ● ● ● ● ●●● ●● ● ●● ● ●●●●● ●●●●● ● ● ●● ● ● ●● ●● ●● ● ●●● ●●● ● ● ● ● ●● ● ● ●●●● ● ● ●● ● ● ●● ● ●● ●● ●● ●● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ●● ●●●● ● ● ● ● ● ● ●● ●●● ● ● ● ●●● ●● ● ● ● ● ●●● ● ●● ●● ● ● ● ● ●● ● ●● ● ● ●● ●●● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ●●●● ● ● ●● ●● ● ●● ●●●● ● ● ●● ●●● ● ● ● ● ● ● ● ● ● ● ● ●●●● ●● ● ●●● ●● ●● ● ● ●● ● ● ●●●● ● ● ●● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ●● ●● ● ● ● ●●● ● ●● ● ● ● ●● ●● ●● ● ● ● ●●● ●● ● ● ●●●● ●● ●● ● ● ●● ● ●● ● ●●● ● ● ●● ● ● ● ● ●● ●● ●● ● ● ● ● ● ●● ● ● ● ●● ●● ● ● ● ●●● ●● ● ●● ● ●● ●● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ●
0.6
0.8
1.0
1.2
1.4
1.6
0.6
0.8
1.0
1.2
1.4
1.6
scale
1.0
1.0
scale
0.8 0.6 0.4
shape
0.6 0.4 0.2
●
●
●
●
●
●
● ●
●
● ●
−0.2
● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ●●● ●●● ●● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ●●● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ●● ● ● ● ● ●● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ●● ● ● ● ●● ●●● ● ● ● ●● ● ● ● ● ● ● ●●●●● ● ●● ●● ● ● ● ● ● ● ● ● ●● ● ●●● ● ● ●● ●●● ●● ● ●●● ●● ● ● ● ● ●● ● ● ●● ●●● ●● ●● ● ● ● ● ● ● ● ●● ●● ● ● ● ●● ●● ●●● ● ●● ● ●●●● ●● ● ● ●● ●●●● ● ●● ● ● ●● ●●● ●●● ● ●●●●● ●● ● ●● ● ●● ● ●●● ● ● ●● ● ● ●● ●● ●● ● ● ● ● ● ● ● ●●● ●●●● ●● ●●● ●● ●● ● ● ● ● ● ●● ● ● ● ●●●●● ●●●● ●● ●● ● ● ●● ●● ●● ●● ●● ● ●● ● ●● ●●● ●●● ●● ●● ● ● ● ● ●● ●● ● ●● ●●● ● ● ●● ● ●●● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ●●● ●● ● ● ●● ●●●●●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ●●● ●●● ●● ● ● ● ● ●● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ●
−0.4
−0.4
−0.2
0.0
●
0.0
0.8
●
● ● ● ● ● ●● ●● ● ● ● ●●● ●●●● ● ● ● ●●● ●● ● ● ● ● ● ●● ●●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ●●● ● ●● ● ● ● ● ● ●● ●●●● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ●● ● ●●●●● ● ● ●● ●● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ●●● ● ● ● ●●● ● ●● ● ● ●● ●● ●● ●● ● ● ● ● ●● ●● ● ● ● ● ●●●● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●●● ●● ● ●●●● ●● ● ● ● ● ●●●● ●● ●● ● ●●● ●● ● ●●● ● ● ●● ● ●● ● ● ●●●● ●● ●●● ●● ● ● ● ● ●● ● ●●● ●● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●● ● ● ●● ●● ●●● ● ●● ●● ●●● ● ● ●● ● ● ● ●● ● ● ● ●●● ●●●● ● ● ●● ●● ● ● ●●● ● ●●● ● ●●● ● ● ● ●● ● ● ● ● ●● ●●● ●● ●● ●●●● ●●●●● ● ● ●● ●● ● ● ●●●● ●● ●● ●● ● ● ●● ● ● ●● ●● ● ● ●● ● ● ● ● ●● ●●● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ●● ●●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●
0.2
●
● ●
shape
1.0
scale
0.4
scale
0.4
0.6
0.8
1.0
1.2
1.4
1.6
scale
0.8
1.0
1.2
1.4
1.6
1.8
2.0
scale
Fig. 2. ML estimates for 500 samples of size 100 from a GPD (c, r = 1) with c = 0.4 (upper), 0 (middle) and 0.4 (lower). The left panels correspond to the old ML parametrization and the right panel to the new one defined by (2). The intersection of the two solid lines displays the true value of both parameters.
The left panels of Fig. 3, the analogue of the ones of Fig. 2 but for GPWMs, corroborates the fact that GPWM estimates are also dependent.
2. Alternative definitions of the GPD parameters To solve the ML estimators issue raised in the previous section, Chavez-Demoulin and Davison [3] proposed an orthogonalization
1218
P. Ribereau et al. / Advances in Water Resources 34 (2011) 1215–1221
●
−0.5
shape
0.0
● ● ●● ● ● ● ●●●● ● ● ● ●● ● ●● ●● ●●● ● ●● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ●● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ●●● ● ● ● ● ● ●● ●● ● ● ● ● ●● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●●● ●●● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ●● ●●● ● ● ●● ● ●● ●● ● ● ●● ● ●● ● ● ● ●● ● ● ●●● ● ● ● ●●● ●●● ●● ● ● ● ● ● ●● ●●
−1.0
0.0 −0.5 −1.0
shape
●
● ● ● ●●● ● ●● ● ● ● ●● ● ● ● ●● ● ●●● ●● ● ●●● ● ● ● ● ● ●● ●●● ●● ● ● ● ● ●●●● ● ● ● ●● ●● ●● ●●●● ● ●●● ●● ● ● ● ● ● ● ●● ●● ● ● ●● ● ● ●● ● ● ●● ● ● ● ●●● ● ● ● ● ● ● ●● ● ● ● ● ●● ●● ● ● ● ●● ● ● ●● ● ● ● ●● ● ●● ● ●●●● ● ●● ●● ● ● ●● ●●●● ● ●● ● ● ●●●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ●●● ● ● ● ● ● ●● ● ● ● ●● ●●●● ● ● ● ● ● ● ●●● ●● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ●●● ● ● ●● ●● ● ● ● ● ●● ●● ●● ●● ● ● ● ● ●●● ● ● ● ● ●● ● ● ●● ● ●●● ● ● ● ● ●● ● ●●● ● ● ● ● ●●● ● ● ● ● ●●● ● ● ● ● ● ●●● ●● ● ●●●● ●● ● ●● ● ● ● ● ● ●●●● ●● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ●●● ● ●● ● ● ● ●●● ● ● ●●● ● ● ●●● ● ● ●● ● ● ●●● ● ● ● ●● ●● ● ●● ● ●●● ● ●● ● ● ●●● ● ● ● ● ● ● ● ●●● ●●● ● ●● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●
● ●
●
● ●
●●
−1.5
−1.5
●
0.5
1.0
1.5
0.0
0.5
scale
−0.6
●
● ● ● ●
0.2 0.0 −0.2
shape ●
●
● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ●● ●● ● ●● ● ● ●● ●●●●●● ● ● ● ●● ● ● ● ● ●●●● ●●●● ● ● ● ● ●● ●● ● ● ● ●● ● ●●● ●●● ●● ●● ● ● ●● ● ● ●● ● ● ● ● ● ●● ●● ●● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ●●●● ●● ●●● ●●● ● ● ● ●●● ●●● ●●●● ●● ● ● ●● ● ●●● ● ● ● ●● ●● ●● ●●● ● ● ●● ● ● ●●● ● ●● ● ●● ● ●● ●● ● ●● ●● ● ●● ● ●● ● ● ●● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●●● ●●● ● ● ● ● ●● ●● ● ●●● ● ● ● ● ● ● ●●● ● ●● ●● ●● ● ● ● ● ●●● ● ●● ● ● ● ● ● ● ● ●●●● ●● ● ● ● ●● ●● ●● ●●●● ● ● ●● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ●●● ●● ● ● ●● ●● ●●●●●● ● ●●● ● ● ●● ● ●●● ● ●●● ●● ●●● ● ● ● ● ● ●●●●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ●●●● ● ●● ● ●● ● ● ● ●● ● ●● ● ● ● ●●●● ●●●● ● ●
0.4
●
−0.4
−0.4
−0.2
0.0
0.2
●
shape
● ●● ● ● ● ● ●● ● ● ● ● ●● ● ● ●● ● ●● ● ● ● ● ●●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ●●● ● ●● ●● ●●● ●● ● ● ● ● ● ● ● ●●● ● ●● ● ●● ● ● ● ●● ●● ●● ●●●● ● ●●● ●●● ●● ● ● ●●● ● ●● ●● ●● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ●●●● ●● ● ●●●● ● ● ● ●● ●● ● ● ●●● ●● ● ● ● ●● ● ● ●●● ● ●● ● ● ● ●●●● ● ●● ●● ● ●● ● ●● ●● ●●● ● ●● ● ● ● ● ● ●●●●● ●●● ● ● ●●● ●●● ●● ●● ●●● ●● ●● ● ● ● ● ●●● ● ●● ●● ●● ● ● ●●● ● ● ●●● ● ●●● ● ●●● ● ● ●● ●● ● ●● ●● ●● ● ● ● ● ●●●● ●●● ● ● ● ●● ● ● ● ● ● ● ●● ●● ● ● ● ● ●●● ● ● ● ●● ● ● ●● ● ● ● ●●● ● ● ● ● ● ● ● ●●● ● ● ●● ●●● ●● ● ● ● ● ●● ● ● ●● ● ● ● ● ●● ● ●● ● ● ● ●● ● ● ● ●● ● ●● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ●● ● ●● ●● ● ●● ● ● ●●● ● ●● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ●
●
● ●●●
●
−0.6
0.4
●
●
●
●
●
−0.8
−0.8
●
0.6
0.8
1.0
1.2
1.4
1.6
0.0
0.2
0.4
0.6
0.8
1.0
scale
1.0
1.0
scale
−0.2
0.8
●
●●● ● ●● ● ● ● ● ● ●● ●● ●●●● ● ● ● ●● ● ● ●● ● ●● ● ●● ● ●●●● ● ● ● ●● ● ● ● ● ● ●● ● ●● ● ●● ● ● ● ● ● ● ● ●●● ●●● ●●●● ●●● ●● ● ● ● ● ● ●● ● ●●●●●● ●● ●● ●● ● ● ● ● ●● ● ●● ● ● ● ●● ●● ● ●● ●●● ●●● ● ●●●● ●● ●● ●● ● ●● ● ●● ● ●● ● ●● ●●● ●● ● ●● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ●● ●● ●● ● ●● ● ● ● ●● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●●● ●● ● ●● ●● ● ● ● ● ● ●●● ●● ●● ● ●●●●●●●● ● ●● ● ● ● ● ● ● ● ● ●●● ● ● ●● ●● ● ● ●●●● ● ●● ●● ●● ●●● ● ●● ●● ●● ●●● ● ●● ● ● ● ●● ●● ●●● ● ● ● ●●● ●● ● ● ●● ●● ● ● ● ● ●●● ●● ●●●● ● ● ●● ●● ● ●● ● ●● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●● ●●● ●● ●● ●● ● ● ● ● ●● ●● ●● ● ● ●● ●● ●● ● ●●● ● ●●● ● ●●● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ●
0.4
shape ●
●
−0.4
−0.4
●
● ●
● ●
●
0.2
●● ● ● ●● ● ● ● ● ● ●● ●● ● ●●● ● ● ● ● ● ● ● ● ●●● ● ● ●● ● ●● ● ●● ● ● ●● ●● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ●●●● ●●● ● ● ●● ● ● ● ● ● ●● ●● ● ● ●● ●●● ●● ●●● ●● ● ● ●●●● ●● ●● ● ● ● ● ● ● ● ●● ● ●● ●●● ● ●● ● ● ● ● ●● ●● ● ● ● ●● ● ●● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●●●● ● ●● ●●● ●● ●● ●●●●●●● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ●● ● ● ●● ●● ● ● ●● ●●● ● ● ●●●●●● ● ● ● ● ● ●● ●● ● ● ● ● ●●●● ● ● ●●●●● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ●●●● ● ●● ●● ● ● ●●● ● ● ● ● ●● ● ●● ● ●●● ●●● ● ●●●● ●● ● ●●● ●● ● ●●● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ●●●● ● ● ●● ●● ● ● ●● ●●● ● ●● ● ● ● ●● ● ●● ● ● ● ● ● ●●● ● ● ● ● ● ●● ● ● ●● ●● ● ●● ● ●● ● ● ● ● ●● ●● ● ● ●● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ●●●● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ●
0.6
●
0.0
0.0
0.2
0.4
0.6
●
● ●
●
−0.2
0.8
●
shape
1.0
scale
0.4
0.6
0.8
1.0
1.2
1.4
1.6
0.2
scale
0.4
0.6
0.8
1.0
1.2
1.4
scale
Fig. 3. GPWM estimates for 500 samples of size 100 from a GPD (c, r = 1) with c = 0.4 (upper), 0 (middle) and 0.4 (lower). The left panels correspond to the old GPWM parametrization and the right panel to the new one defined by (4). The intersection of the two solid lines displays the true value of both parameters.
technique, see, e.g. [7], to obtain a diagonal Fisher information matrix. This leads to the following reparametrization
cI;ML :¼ c; rI;ML :¼ rð1 þ cÞ:
ð2Þ
The two advantages of this parametrization are that the shape parameter c remains unchanged and that the scale parameter r is only modified by a multiplicative factor (1 + c) which does not contain r. This latter remark implies that the GPD scaling property remains valid with the new parametrization, i.e. if X follows a GPD
1219
P. Ribereau et al. / Advances in Water Resources 34 (2011) 1215–1221
with scale parameter rI,ML then, for any positive scalar a, the variable aX also follows a GPD with scale parameter a rI,ML. Hence the interpretation of (2) is not fully identical to the old one but it is very comparable. The important bonus is that system (2) allows for individual interpretation and modeling of each GPD parameter. Fig. 2 compares the old parametrization (left panels) and the new one (right panels). It is clear that the dependence between ML estimates has been reduced for the three values of c = 0.4, 0 and 0.4, see the upper, middle and lower panels, respectively. But even with the new parametrization, one drawback of the ML estimation is still visible because of the constraint c > 0.5. Near 0.5, the ML procedure is not ideal and the upper right panel clearly displays a type of non-linear dependence for c = 0.4, see also the first column of Table 1. One could argue that negative shape parameters are unlikely to occur in climatology but recent work by Friederichs [13] hints that heavy wind gusts in Germany can reach negative shape parameters around 0.3, see Fig. 7 of [13]. Also for UK precipitation, negative shape parameters have been found along the west coast (see [20,21]). Such bounded tails can also appear while modeling relative humidity extremes (see [12]). With respect to the method-of-moments, there are currently not adequate reparametrizations for neither the PWM nor the GPWM approach. The next two paragraphs cover both cases. The matrix diagonalization detailed in Appendix A allows us to find the following new parametrization for the PWM case
cI;PWM :¼ c; rI;PWM :¼ raðcÞ;
ð3Þ
where the function a(c) is defined as follows
1 11 2 logð1 cÞ logð1 c þ 2c2 Þ þ logð2 cÞ 2 28 7 pffiffiffi 7 1 þ 4c pffiffiffi þ arctan 14 7
log aðcÞ ¼
for c < 1. Although not trivial, this expression is explicit. As for the ML parametrization captured by (2), the shape parameter is unchanged and the GPD scaling property is still valid. To make another link with (2), we can notice that the function a(c) via a Taylor expansion in the neighborhood of zero can be approximated by a0 (1 + c) with a0 = 1.138544. Concerning the GPWM approach, another matrix diagonalization that can also be found in Appendix B provides the following system
cI;GPWM :¼ c; rI;GPWM :¼ rbðcÞ;
ð4Þ
where the function b(c) is defined as follows
16 8 logð5 2cÞ logð13 10c þ 8c2 Þ 19 19 pffiffiffiffiffiffi 64 pffiffiffiffiffiffi 1 79 arctan ð10 þ 16cÞ 79 þ 1501 158
log bðcÞ ¼ logð2 cÞ þ
for c < 2. Again this expression is explicit. It preserves the shape parameter and the GPD scaling property. A Taylor expansion of
Table 2 Coverage probability of the confidence interval from 1000 GPD samples of size 100 for the different parameters and for different values of the shape parameter c and r = 1. Parameter/ Value of c
0.4
0.2
0
0.2
0.4
0.8
1.2
c with ML r with ML rI,ML from (2)
90.6 94.7 98.0
82.3 91.7 92.4
86.8 93.5 96.8
90.7 93.2 97.8
90.4 92.7 98.2
93.3 94.6 98.7
93.0a 94.1a 97.9a
c with GPWM r with GPWM rI,GPWM with (4)
96.2 96.0 94.1
96.2 96.4 95.0
96.0 95.8 93.3
96.6 96.2 94.6
94.9 94.9 93.9
97.2 98.7 97.1
95.1 94.7 95.9
a After removing all the samples in which the optimization algorithm did not converge or leaded to aberrant estimations.
order 2 of the function b(c) leads to the following approximation: 2 b0 þ b1 c þ b2 c2 with b0 = 0.5422619, b1 = 0.4046108 and b2 = 0.3265764. Hence the system (4) represents a greater departure in terms of parametrization from the ML case than the GPWM one. Fig. 3 can be compared to Fig. 2. Both ML and GPWM appears to work well in terms of decorrelation with a slight advantage for the GPWM, see also Table 1. Table 2 represents the coverage probabilities of the confidence intervals obtained for c and r (before and after the reparametrization) with level 95% for different values of the shape parameters. These results are based on the asymptotic covariance matrices given in Appendix C. As observed in this table, all the coverage probabilities are close to the nominal coverage probability, except for the ML method and c = 0.4 which is close to the value 0.5 required for the asymptotic normality. In the case where c < 0.5, the method is super-efficient (see [22]) and it is still possible to construct confidence interval by bootstrap. However, for all the values of c the GPWM method leads to better results than the ML one since the intervals in that case are conservative whatever the value of c is. 3. Conclusion This short article discusses correlation issues between the two GPD parameters estimators and the risk of interpreting such estimators. Such a problem has rarely been a point of interest within hydrological and climatological circles but, despite its simplicity, it could become more prevalent as researchers have the tendency to model spatial and temporal structures at the GPD parameters level. To decorrelate these parameters, classical matrix diagonalization techniques have been applied here and comparisons have been implemented among three inferences techniques (ML, PWM and GPWM). Overall our new GPWM parametrization appears to be the most straightforward and very fast approach in terms of its implementation (no optimization needed to maximize a likelihood). It offers an adequate decorrelation scheme and can be applied for a wide range of shape parameters. If the practitioner only focuses on the estimation of return levels in the i.i.d. case, then the proposed parametrization may not be needed because return levels combine uncertainties from the shape and scale parameters estimates. But in terms of modeling with covariates, we
Table 1 Estimated correlations obtained from 1000 GPD samples of size 100 for different values of the shape parameter c and r = 1.
a
Inference on c=
0.4
0.2
0
0.2
0.4
0.8
1.2
Classical ML ML with (2)
0.789 0.039
0.822 0.147
0.738 0.098
0.675 0.029
0.644 0.025
0.521 0.003
0.255a 0.154a
Classical GPWM GPWM with (4)
0.759 0.019
0.861 0.031
0.803 0.054
0.759 0.009
0.691 0.013
0.532 0.027
0.423 0.226
After removing all the samples in which the optimization algorithm did not converge or leaded to aberrant estimations.
1220
P. Ribereau et al. / Advances in Water Resources 34 (2011) 1215–1221
believe that decoupling the GPD parameters could facilitate the inference and the interpretation in many applications. In this paper we did not treat the choice of the threshold that define excesses. The influence of such a threshold in terms on parameters estimation and correlation must be non-negligible. This is a delicate statistical problem because the threshold is not viewed within the EVT paradigm as a parameter to be estimated but rather as a constant that has to be chosen a priori (see, e.g. [1]). Finally, the correlation issue investigated in this work is not circumvented to the GPD or to the three estimation methods studied here. For example, a powerful alternative for visualizing standard errors of the GPD estimates resides in computing log likelihood profiles for a given GPD parametrization, see Section 2.6.6 of Coles book. This approach has the advantage of producing non-symmetrical confidence intervals. However it is not clear to us what is the optimal parametrization for this method. Also, the posterior of a classical GPD Bayesian analysis should suffer from the same problem. In terms of distribution, strong correlations also exist for the Generalized Extreme Value (GEV) distribution parameters used to model maxima. Therefore it would be of interest to find new parametrizations for the GEV.
Appendix B. Generalized Probability Weighted Moment estimators
Acknowledgments
This justifies the approximatively independent reparametrization:
Part of this work has been supported by the EU-FP7 ACQWA Project (www.acqwa.ch) under Contract No. 212250, by the PEPER-GIS project, by the ANR-MOPERA project, by the ANR-McSim project and by the MIRACCLE-GICC project. The authors also would like to thanks the referees for valuable comments that have greatly improved the paper.
The GPWM method is explained in [10] for the GPD and in [9] for the GEV distribution. Using the expression of the asymptotic variance given in [10], we get:
0 5 b ðcÞ 5 13 c bðcÞ c 5c þ 4c2 ð2 cÞ 2 bðcÞ 2 2 97 þ 39c 23c2 þ 4c3 ; 4
rð2 cÞ
or similarly
97 þ 39c 23c2 þ 4c3 b0 ðcÞ 4 : ¼ bðcÞ ðc 2Þ 52 c 13 5c þ 4c2 2 Again, by integration, we derive
16 8 logð5 2cÞ logð13 10c þ 8c2 Þ 19 19 pffiffiffiffiffiffi 64 pffiffiffiffiffiffi 1 79 arctan ð10 þ 16cÞ 79 : þ 1501 158
log bðcÞ ¼ logð2 cÞ þ
8 < cI;GPWM ¼ c;
2 :r ~ I;GPWM ¼ r b0 þ b1 c þ b2 c2
ðB:1Þ
with b0 = 0.5422619, b1 = 0.4046108 and b2 = 0.3265764. Appendix C. Asymptotic covariance matrix
Appendix A. Probability Weighted Moments estimators
C.1. Maximum Likelihood
Using the asymptotic covariance matrix RPWM, we want to find the function a() such that the following matrix is diagonal:
Using the reparametrization (2), the estimators obtained by the Maximum Likelihood method are asymptotically normally distributed with covariance matrix given by:
1
0
ra0 ðcÞ aðcÞ
RPWM
ra0 ðcÞ ; 0 aðcÞ
1
!
ð1 þ cÞ2
0
0
2r2 ð1 þ cÞ2 ð1 þ 2cÞ
:
where a0 ¼ @@ac. This leads to the following equation
C.2. Probability Weighted Moments
a0 ðcÞ ð1 cÞðc 2Þð1 c þ 2c2 Þ aðcÞ þ2 6c þ 7c2 2c3 ¼ 0:
rðc 2ÞaðcÞ
Using the reparametrization (3), the corresponding estimators are asymptotically normally distributed with covariance matrix given by:
0
This equation can be rewritten as
@
a0 ðcÞ 2 6c þ 7c2 2c3 ¼ : aðcÞ ð1 cÞð2 cÞð1 c þ 2c2 Þ
This latter leads to the following approximatively independent reparametrization
with a0 = 1.138544.
0
r2 a2 ðcÞ ð1cÞð1cþ2c2 Þ
A:
With the classical GPD parametrization, the asymptotic covariance matrix of the estimators obtained by the GPWM method is the following:
1 11 2 logð1 cÞ logð1 c þ 2c2 Þ þ logð2 cÞ 2 28 7 pffiffiffi 7 1 þ 4c pffiffiffi : þ arctan 14 7
log aðcÞ ¼
cI;PWM ¼ c; r~ I;PWM :¼ a0 ð1 þ cÞ
0
C.3. Generalized Probability Weighted Moments
Finally, by integration we deduce that
1
ð1cÞð2cÞ2 ð1cþ2c2 Þ ð12cÞð32cÞ
ðA:1Þ
Cð1;1Þ rCð1;2Þ rCð1;2Þ r2 Cð2;2Þ
!
with the explicit expression of C(i,j) given in [10]. Now using the parametrization (4), the new asymptotic matrix is given by:
0 @
Cð1;1Þ 0
2
0
1 2 A: ð1;2Þ
ðcÞ r2 Cb ð1;1Þ Cð1;1Þ Cð2;2Þ C
P. Ribereau et al. / Advances in Water Resources 34 (2011) 1215–1221
References [1] Beirlant J, Dierckx G, Goegebeur Y, Matthys G. Tail index estimation and an exponential regression model. Extremes 1999;2:177–200. [2] Beirlant J, Goegebeur Y, Segers J, Teugels J. Statistics of extremes: theory and applications. Wiley series in probability and statistics, 2004. [3] Chavez-Demoulin V, Davison AC. Generalized additive modelling of sample extremes. Appl Stat 2005;54:207–22. [4] Coles SG. An introduction to statistical modeling of extreme values. Springer; 2001. [5] Coles SG, Dixon MJ. Likelihood-based inference of extreme value models. Extremes 1999;2:5–23. [6] Cooley D, Nychka D, Naveau P. Bayesian spatial modeling of extreme precipitation return levels. J Am Stat Assoc 2007;102:824–40. [7] Cox DR, Reid N. Parameter orthogonality and approximate conditional inference (with discussion). J Roy Stat Soc B 1987;49:1–39. [8] Davison A, Smith R. Models for exceedances over high thresholds (with discussion). J Roy Stat Soc B 1990;52:393–442. [9] Diebolt J, Guillou A, Naveau P, Ribereau P. Improving probability-weighted moment methods for the generalized extreme value distribution. REVSTAT Stat J 2008;6:33–50. [10] Diebolt J, Guillou A, Rached I. Approximation of the distribution of excesses through a generalized probability-weighted moments method. J Statist Plann Inf 2007;137:841–57. [11] Embrechts P, Klüppelberg C, Mikosch T. Modelling extremal events for insurance and finance. Applications of Mathematics, vol. 33. Berlin: SpringerVerlag; 1997.
1221
[12] Flecher C, Allard D, Naveau P. Truncated skew-normal distributions: moments, estimation by weighted moments and application to climatic data. METRON – Int J Stat 2010;LXVIII(3):265–79. [13] Friederichs P, Gober M, Bentzien S, Lenz A, Krampit R. A probabilistic analysis of wind gusts using extreme value statistics. Meteorol Z 2009;18: 615–29. [14] Greenwood JA, Landwehr JM, Matalas NC, Wallis JR. Probability-weighted moments: definition and relation to parameters of several distributions expressable in inverse form. Water Resour Res 1979;15:1049–54. [15] de Haan L, Ferreira A. Extreme value theory: an introduction. Springer series in operations research, 2006. [16] Hosking JRM, Wallis JR. Parameter and quantile estimation for the generalized Pareto distribution. Technometrics 1987;29:339–49. [17] IPCC Climate Change: Synthesis Report. http://www.ipcc.ch/ipccreports/, 2007. [18] Katz R, Parlange M, Naveau P. Extremes in hydrology. Adv Water Resour 2002;25:1287–304. [19] Landwehr J, Matalas N, Wallis JR. Probability weighted moments compared with some traditional techniques in estimating Gumbel parameters and quantiles. Water Resour Res 1979;15:1055–64. [20] Maraun D, Rust HW, Osborn TJ. The annual cycle of heavy precipitation accross the UK: a model based on extreme value statistics. Int J Climatol 2009;29: 1731–44. [21] Maraun D, Rust HW, Osborn TJ. The influence of synoptic airflow on UK daily precipitation extremes. Part I: observed spatio-temporal relations. Clim Dyn 2010;36:261–75. [22] Smith RL. Maximum likelihood estimation in a class of nonregular cases. Biometrika 1985;72:69–90.