Estimating return levels from maxima of non-stationary random sequences using the Generalized PWM method

Since the pioneering work of Landwehr et al. (1979), Hosking et al.(1985) and their collaborators, the Probability Weighted Moments (PWM) method has been very popular, simple and efficient to estimate the parameters of the Generalized Extreme Value (GEV) distribution when modeling the distribution of maxima (e.g., annual maxima of precipitations) in the Identically and Independently Distributed (IID) context. When the IID assumption is not satisfied, a flexible alternative, the Maximum Likelihood Estimation (MLE) approach offers an elegant way to handle nonstationarities by letting the GEV parameters to be time dependent. Despite its qualities, the MLE applied to the GEV distribution does not always provide accurate return level estimates, especially for small sample sizes or heavy tails. These drawbacks are particularly true in some non-stationary situations. To reduce these negative effects, we propose to extend the PWM method to a more general framework that enables us to model temporal covariates and provide accurate GEV-based return levels. Theoretical properties of our estimators are discussed. Small and moderate sample sizes simulations in a non-stationary context are analyzed and two brief applications to annual maxima of CO 2 and seasonal maxima of cumulated daily precipitations are presented.


Introduction
Extreme value theory provides a solid mathematical foundation (e.g.Embrechts et al., 1997;Beirlant et al., 2004;de Haan and Ferreira, 2006) for studying maxima in fields like hydrology or climatology (e.g.Katz et al., 2002).In the Identically and Independently Distributed (IID) setup, this theory states that maxima should follow the Generalized Extreme Correspondence to: P. Ribereau (pribere@math.univ-montp2.fr)Value (GEV) distribution whenever re-normalized maxima of a random sample converges to a non-degenerate random variable (e.g.Coles, 2001).From a statistical point of view, this means that the cumulative probability distribution function of maxima from IID samples is very likely to be correctly fitted by the following GEV distribution where σ >0, γ =0 and µ ∈ R are called the GEV scale, shape and location parameters, respectively, and with the constraint 1+γ x−µ σ >0.If γ →0, then Eq. ( 1) corresponds to the Gumbel case and is equal to exp − exp{− x−µ σ } with x ∈ R. To estimate the three parameters of a GEV distribution, several methods have been developed, studied and compared during the last twenty years.In 1979, Greenwood et al. (1979) and Landwehr et al. (1979) introduced the socalled Probability Weighted Moments (PWM).This methodof-moments approach (see Appendix A) has been very popular in hydrology (Hosking et al., 1985) and climatology because of its conceptual simplicity, its easy implementation and its good performance for most distributions encountered in geosciences.In 1985, Smith studied and implemented the Maximum Likelihood (ML) method for the GEV density (see Appendix C).According to Hosking et al. (1985), the PWM approach is superior to the MLE for small GEV distributed samples.Coles and Dixon (1999) argued that the PWM method assumes a priori on the shape parameter which is equivalent to assume a finite mean for the studied distribution.To integrate this condition in the ML approach, these authors proposed a penalized MLE scheme with the constraint γ <1.If this condition is satisfied, then the ML approach is as competitive as the PWM one, even for small samples.Still, the debate over the advantages and drawbacks of both estimation methods is not closed.For example, the classical and penalized ML approaches impose a restriction on the lower values of γ , i.e. we need γ >−0.Diebolt et al. (2008) introduced the concept of the General-ized Probability Weighted Moments (GPWM) for the GEV.It broadens the domain of validity of the PWM approach, allowing heavier tails to be fitted.Despite this advantage, the ML method has kept a strong advantage over a PWM approach: its inherent flexibility in a non-stationary context.
When studying climatological and hydrological data, it is not always possible to assume that the distribution of the maxima remains unchanged in time.For example, trends can be present in extreme values of different hydroclimatological series (e.g., Kharin et al., 2007;IPCC Report, 2007).
The MLE can easily integrate temporal covariates within the GEV parameters (e.g., Katz et al., 2002;Coles, 2001;El Adlouni et al., 2007) and conceptually, the MLE procedure remains the same if the GEV parameters vary in time (see El Adlouni and Ouarda (2008) for a comparison study of different methods for non-stationary GEV models).In practice, numerical problems quickly arise and estimated ML return levels can be misleading in some non-stationary situations (see Tables 1 and 2), especially for strong heavy tails (γ ≥0.4).With these limitations in mind, our aim is to propose and to study a novel GPWM procedure that can handle temporal covariates.Our main motivation is to keep the interesting GPWM properties identified in the IID case while adding the needed flexibility to handle non-stationarities that are often present in real case studies.This will provide a valuable alternative to the MLE for non-stationary GEV distributed data.
Before closing this introduction, we would like to emphasize that, beyond the three aforementioned estimation methods (MLE, PWM and GPWM), there exists other variants and extensions.For example, Hosking (1990) proposed and studied the L-moments, Zhang (2007) recently derived a new and interesting method-of-moments, and Bayesian estimation has also generated a lot of interest in extreme value analysis (e.g.Lye et al., 1993;Coles, 2001;Cooley et al., 2007).But we will neither compare nor discuss these alternatives with respect to our proposed method for the following reasons.The objective of this work is not to write a review paper that compares all existing estimation methods.Instead our aim is to extend the applicability of a wellknown approach (PWM) by working with a larger class of estimators (GPWM) and by proposing an algorithm that encompasses temporal dependences.We focus on comparing our approach with the classical MLE because the later may be the most popular in statistical climatology and hydrology for non-stationary time series analysis.

A GEV regression model
One of the most simples, most frequently used and most studied models in statistics is the classical regression where Y=(Y 1 , . . ., Y n ) t represents an observational vector of length n, X is a n×p known matrix of explanatory variables and β a vector of unknown parameters of length p that characterizes the relationship between observations and explanatory variables.For example, in a time series context where a cycle of length T can be present, a model like Y i =β 0 +β 1 cos(2π i/T )+ i corresponds to a matrix A regression like Eq. ( 2) can also model linear and even polynomial trends and consequently, it can handle various types of non-stationarities in time.Classically, the vector in Eq. ( 2) is assumed to be a zero-mean Gaussian vector, but this hypothesis is not reasonable whenever the observations can be considered as maxima.As already mentioned in the Introduction, Extreme Value Theory tells us that a more adequate fit for maxima should be the GEV distribution defined by Eq. ( 1).For this reason, we assume in this paper that the vector in Eq. ( 2) consists of IID random variables from a GEV distribution defined by Eq. ( 1) and with parameters (0, σ, γ ).Note that µ is set to zero here because β 0 in the vector β usually plays the role of the location parameter.Overall, this is equivalent to state that the sequence of maximum Y i represents a sequence of independent GEV distributed random variables with the same scale parameter σ , the same shape parameter γ but with a varying location parameter µ i that depends on the covariate X, i.e. µ i =(Xβ) i .Classically, most estimation methods for estimating β in a regression model assume that the "noise" is zero-mean.In our case, we have imposed that i follows a GEV distribution with µ=0.But this does not imply that the mean i is null because the GEV density is not always symmetric around µ. To be in accordance with the zero-mean constraint, we just have to re-parametrize our model where and E( 1 ) corresponds to the mean value of 1 .
To estimate return levels, we propose to implement the following three-steps algorithm: Step a: Implement an existing regression method on Eq. ( 4) to find regression estimates that we call ( β 1 , . . ., β p−1 ).Step b: Apply the GPWM approach described in Appendix B to estimate the GEV parameters of the "pseudo-residuals " Step c: From the GEV estimates obtained from Step b, compute the desired return levels.
Although each of these three steps appears to be simple, they deserve careful examinations.In order to apply Step b, the pseudo-residuals ( ε 1 , . . ., ε n ) should be IID and GEV distributed.This is not true because β j has to be estimated in Step a (it would be true if we knew β j ).This issue is less relevant if the sample size n gets larger.For example, if the classical least squares regression 1 is used for Step a, our estimation of β j asymptotically becomes better and our IID GEV assumption for our residuals truer.More precisely, a the variance of all i is finite (i.e.γ < 1 2 ) the pseudo-residuals ε i can be asymptotically viewed as independent random variables distributed according to a GEV(β 0 , σ, γ ) distribution.This result is general because it does not impose a specific type of distribution for the in Eq. ( 2) but only a finite variance condition.Still extreme values like maxima do not necessary obey the condition γ < 1 2 , it may be more appropriate to implement a robust/resistant regression in Step a.
Such a remark was confirmed by our simulation study.In Sect. 3 we take advantage of the Least Trimmed Squares (LTS) regression (Venables and Ripley, 2002, p. 156-163).Basically, the influence of very high values that could mislead the estimation of the β j 's is trimmed by minimizing an error/cost function that is only based on the core data.This LTS method usually gives accurate estimators even in the presence of large values but it needs a long computation time.This is not an important issue for maxima because of the rarity of the observations.Hence, for Step a, we advice to apply such a resistant/robust method for heavy tail distributions.
Concerning Step c, the return level z t for the fixed time period t can be easily defined in the IID case.More precisely, the return z t corresponds to the 1−1/t quantile, i.e. it is the level that, in average, is crossed one time during the time period t.For the GEV(µ, σ, γ ) distribution, this means that In a non-stationary context, no unique definition exists for z t because extreme values recorded during the time 1, . . ., n have a different distribution than the extremes occurring during any translated time 1+C, . . ., n+C for any C =0.With respect to our regression model Eq. ( 4), the IID return level defined by Eq. ( 5) can be adapted in the following way This means the trend at time t plus an IID return level based on the residual of the regression.
To estimate beyond the range of observations, e.g. to compute the centennial return level with only 50 years of data, one needs to extrapolate the trends in Eq. ( 9) captured by (Xβ) t = µ t .As it is well know in statistics, it is very dangerous to extrapolate a linear regression outside the sample support.Hence, extra caution is required when interpreting z * t in real applications because it means that the trend will remain the same until time t.
3 Analysis of non-stationary maxima
To estimate all parameters and return levels with our approach the same procedure is followed.Step a is implemented by regressing with a LTS method.Then Step b is applied.Finally, Step c is used to estimate the return level z * t defined by Eq. ( 9) for t=10×n.
To obtain the return level estimate from the ML approach, we directly obtain estimation of β 0 , β 1 , γ and σ by maximizing the log-likelihood function (see Appendix C).Abnormal ML estimates have been eliminated by removing all samples that provide ML estimates of γ greater than 2.52 .
For both estimation methods, the bias between true and estimated return values is displayed in Table 1 and the standard deviation in Table 2.These tables clearly showed that, for this specific example, the MLE does not provide adequate estimates whenever the sample size is small (n<50) and the shape parameter is large, especially in terms of standard deviation.In contrast, our GPWM regression method does a better job at handling small samples and large γ (except for γ =1 where both methods fail), especially in terms of standard deviation.Other simulations have been done, in particular a comparison between the two methods in a linear dependence case.Since the conclusions are the same, we do not include them.

Annual maxima of CO 2 concentrations
We consider the annual concentrations of CO 2 at Amsterdam in the Indian Ocean from 1981 to 2001.These observations are available on the website of RAMCES (http://www.ipsl.jussieu.fr/services/Observations/fr/RAMCES.htm)where an extensive bibliography can be found on this subject, for example in Gros et al. (1999).Each of the 21 measures is assumed to follow a GEV(β 0 +β 1 ×year, σ, γ ) distribution.With the GPWM and ML approaches, we have estimated γ , σ , β 0 and β 1 and also the return level for a period of 200 years.The results are listed in Table 3.We observe that the two methods lead to very similar results for all the parameters estimated.This corroborates our simulations analysis.For light tailed random variables such as CO 2 concentration maxima, the GEV parameters estimates are fairly similar for values of γ near zero, see Table 1.

Seasonal maximum of cumulated daily precipitations
We treat here seasonal maxima of cumulated daily precipitations recorded in the Orgeval basin (France) during 31 years.
We suppose that each of the 124 measures (Y i,j ) i=1,...,31;j =1,..,4 follows a GEV β 0 +β 1 × cos π 2 ×j +β 2 ×i, σ, γ distribution.The indice i corresponds to the year while the indice j corresponds to the season.The parameter β 2 indicates the trend of the series.Table 4 provides the estimates of γ , σ , β 0 , β 1 and β 2 and the 100-year return level obtained from the GPWM and ML approaches.
In contrast to our previous example summarized by Table 3, the MLE and GPWM procedures give very different estimates of the 100-year return level estimates in Table 4.To better understand this difference, we re-estimate the GEV parameters but with only the first 20 years.Table 5 summarizes our findings.As observed in Tables 1 and 2, the properties of a GPWM estimated shape parameter around 0-0.2 does not vary greatly when the sample size changes.This is not the case for the ML approach where the estimates of γ can double by changing the sample size, see Table 5.This also impacts the estimates of the return levels.

Conclusions
In this paper, we propose an extension of the PWM method which can be viewed as an alternative to the MLE method that enable us to model temporal covariates and provide accurate return levels.We illustrate our approach by a simulation study and by applying our approach to two time series of maxima.Compared to the ML method, the GPWM method performs better and is computationally easy, but we cannot obtain confidence intervals and, up to now, we can only have a non-stationary location parameter, whereas with the MLE method the three parameters of the GEV distribution can be non-stationary.Concerning the estimation of return levels in such non-stationary cases, we would like to conclude with a word of caution.Extrapolating the trend beyond the range of observations is always a delicate and sometimes dangerous operation (since we assume that the trend will remain the same in the future), especially when dealing with extremes.Hence high return levels in a non-stationary context must be interpreted with extreme care.

Table 1 .
5 to have 1034 P. Ribereau et al.: Return levels for non-stationary maxima Bias of the estimated return levels obtained by the GPWM and ML method (for a return period t=10×n).

Table 2 .
Standard deviation of the estimated return levels obtained by the GPWM and ML method (for a return period t=10×n).