Transcript Extreme Value Distributions

Methodological summary of flood
frequency analysis
A.Zempléni
(Eötvös Loránd University, Budapest)
13.04.2004
Analysis of extreme values
• Classical methods: based on annual maxima
• Peaks-over-threshold methods: utilize all floods
higher than a given (high) threshold.
• Multivariate modelling
– Bayesian approach (dependence among
parameters)
– Joint behaviour of extremes
Extreme-value distributions
Let X1, X2,…,Xn be independent, identically distributed
random variables. If we can find norming constants an, bn
such that
[max(X1, X2,…, Xn)-an]/ bn
has a nondegenerate limit, then this limit is necessarily a
max-stable or so-called extreme value distribution.
The conditions are related to the smoothness of the density
of the sample elements, are fulfilled by all of the important
parametric families.
Characterisation of extreme-value
distributions
• Limit distributions of normalised maxima:

F ( x)  exp(  x )
Frechet:
(x>0)
 is a positive parameter.

F
(
x
)

exp(

(

x
)
) (x<0)
Weibull:

F ( x)  exp(  exp(  x))
Gumbel:
(Location and scale parameters can be
incorporated.)
Another parametrisation
The distribution function of the generalised
extreme-value (GEV) distribution:
1/ 
 


z


z




 0.
G , , ( z )  exp  1   
if 1  

 








: location, : scale, : shape parameters;
>0 corresponds to Frechet,
=0 to Gumbel
<0 to Weibull distribution
Examples
for GEVdensities
Check the conditions
• Are the observations (annual maxima)
– independent? It can be accepted for most of the
stations.
– identically distributed? Check by
• comparing different parts of the sample. For details,
see the next talk.
• fitting models, where time is a covariate.
– follow the GEV distribution?
Tests for GEV distributions
• Motivation: limit distribution of the maximum of
normalised iid random variables is GEV, but
– the conditions are not always fulfilled
– in our finite world the asymptotics is not always
realistic
• Usual goodness-of-fit tests:
– Kolmogorov-Smirnov
– χ2
Not sensitive for the tails
Alternatives

•
2
(
F
(
x
)

F
(
x
))
dF ( x)
Anderson-Darling test: A2   n
F ( x)(1  F ( x))

Computation:
n
A2  n   (2i  1)(log zi  log( 1  z n 1i )) / n
i 1
where zi=F(Xi). Sensitive in both tails.

2
• Modification:
(
F
(
x
)

F
(
x
))
2
n
B 


1  F ( x)
dF ( x)
(for maximum; upper tails). Its computation:
n
n
i 1
i 1
B 2  n / 2   (2i  1) log( 1  z n 1i ) / n   zi / n
Further alternatives
• Another test can be based on the stability property
of the GEV distributions: for any m N there exist
am, bm such that F(x)=Fm(amx+bm) (x R)
The test statistics:
2
h(a, b)  sup n Fn ( x)  Fn (a2 x  b2 )
x
Alternatives for estimation:
• To find a,b which minimize h(a,b) (computerintensive algorithm needed).
• To estimate the GEV parameters by maximum
likelihood and plug these in to the stability
property.
Limit distributions
• Distribution-free for the case of known
parameters. For example:
lim sup n Fn ( x)  Fn (a2 x  b2 )  sup B( y )  y B
2
n 
x
 y
y
where B denotes the Brownian Bridge over [0,1].
• As the limits are functionals of the normal
distribution, the effect of parameter estimation by
maximum likelihood can be taken into account by
transforming the covariance structure.
• In practice: simulated critical values can also be
used (advantage: small-sample cases).
Power studies
• For typical alternatives, the test A-D seems to
outperform B. The power of h very much depends
on the shape of the underlying distribution.
• The probability of correct decision (p=0.05):
n
Test
Distr.
100 200 400 100
200
200
400
NB
exp
Normal
B
0.02 0.27 0.49
0.17
0.58
0.05
0.08
A-D
0.31 0.62 0.96
0.72
0.97
0.21
0.34
h
0.67 0.87 0.99
0.75
0.91
0.10
0.14
Applications
• For specific cases, where the upper tails play the
important role (e.g. modified maximal values of
real flood data), B is the most sensitive.
• When applying the above tests for the flood data
(annual maxima; windows of size 50), there were
only a couple of cases when the GEV hypothesis
had to be rejected at the level of 95%.
• Possible reasons: changes in river bed properties
(shape, vegetation etc).
An example for rejection:
Szolnok water level, 1931-80
Q-Q plot
400
0.0000
0.0005
500
0.0010
0.0015
Observed values
cm
600
700
0.0020
800
0.0025
900
0.0030
Histogram
300
500
700
900
400
500 600
700 800
Model
cm
900
Estimation methods
• Maximum likelihood, based on the unified
parametrisation (GEV) is the most widely
used, with optimal asymptotic properties, if
ξ>-0.5 (it is superefficient for -0.5>ξ>-1).
We have applied it, with good results.
• Probability-weighted moments (PWM)
• Method of L-moments
Robustness of maximum
likelihood estimators
• The effect of small observations is limited:
in our case (negative shape parameters)
halving the smallest 3 values, the difference
in return level estimators was not more than
5-8%.
• However, for positive shape parameters the
effect of smaller values seem to be larger.
Further investigations
• Confidence bounds should be calculated, possible
methods
– based on asymptotic properties of maximum likelihood
estimator
– profile likelihood
– resampling methods (bootstrap, jackknife)
– Bayesian approach
• Estimates for return levels, including confidence
bounds
Confidence intervals
• For maximum likelihood:
– By asymptotic normality of the estimator:
ˆ ~ N ( , i ,i )
where  i,i is the (i,i)th element of the inverse of
the information matrix
– By profile likelihood
• For other nonparametric methods by
bootstrap.
Profile likelihood
• One part of the parameter vector is fixed, the
maximization is with respect the other
components: l ( )  max l ( , )
p
i
 i
i
i
l() is the log-likelihood function; =(i , -i )
Let X1,…,Xn be iid observations. Under the regularity
conditions for the maximum likelihood estimator,
asymptotically
2{l (ˆ)  l p ( i )} ~  k2
(a chi-squared distribution with k degrees of
freedom, if i is a k-dimensional vector).
Use of the profile likelihood
• Confidence interval construction for a parameter
of interest: C  { i : 2{l (ˆ)  l p ( i )}  c }
where c is the 1- quantile of the 12 distribution.
• Testing nested models:
M1() vs. M0 (the first k components of  =0).
l1( M1 ), l0 (M0 ) are the maximized log-likelihood
functions and D:=2{l1( M1 )- l0 (M0 )}.
M0 is rejected in favor of M1 if D>c
(c is the 1- quantile of the k2 distribution).
Return levels
• zp: return level, associated with the return period
1/p (the expected time for a level higher than zp to
appear is 1/p): G , , ( z p )  1  p.

(1  y p ) if   0

   log y p
if  = 0
• The quantiles of the GEV: z p   
where y p  log( 1  p)
• Remark: the probability that it actually appears
before time 1/p is more than 0.5 (approx. 0.63 if p
is small)
Return level plots
10
It can be used for diagnostics,
if the observed data points
are also plotted.
0
5
return level
on a logarithmic scale
-Linear if  = 0



-Convex, with a limit

if  < 0
-Concave, if if  > 0.
15
Continuous:  = 0.2
broken:  = -0.2
1
5
10
50
return period (years)
100
500
1000
-643
-644
-645
-646
Profile Log-likelihood
Profile likelihood
can be calculated
(the return level is
considered as one of
the parameters)
-642
Example: profile likelihood for 100year return level (Vásárosnamény)
900
920
940
Return Level
960
Investigation of the estimators
• Backtest: estimators based on data from a shorter
window. Quite often too many floods are observed
above the estimated level - simulation studies may
confirm if this is a significant deviation from the
iid case (for details see a later talk about
resampling techniques).
• Alternative model: linear trend in the location
parameter (the other parameters are supposed to
be constant).
• Centred time-scale is used: t*=(t-50.5)
Some results with time-varying
location parameter
Location, type Linear estimator
for 
Tivadar,h
509.3+0.761t*
Increment of
loglikelihood
0.75
Tivadar,q
1225 +1.774t*
0.1
Namény, h
616.0+1.205t*
4.18
Szolnok, h
644.4+1.321t*
5.08
Polgár, h
520.8 +1.190t*
5.70
Polgár, q
1709 +0.455t*
0.01
Peaks over threshold methods
If the conditions of the theorem about the GEV-limit
of the normalised maxima hold, the conditional probability
of X-u, under the condition that X>u, can be given as
1 / 
 y 
H ( y )  1  1  ~  ,
  
~ )  0 , where ~     (u   ).
if y>0 and (1  y / 
H(y) is the so called generalized Pareto distribution
(GPD).
 is the same as the shape parameter of the
corresponding GEV distribution.
0.0
0.5
1.0
1.5
2.0
2.5
3.0
Densities of GPD with =1; solid: =0.5, dotted:
=-0.1, dots-and-lines: =-0.7, broken: =-1.3
0
1
2
3
4
5
Peaks over threshold methods
• Advantages:
– More data can be used
– Estimators are not affected by the small
“floods”
• Disadvantages:
– Dependence on threshold choice
– Original daily observations are dependent;
declustering not always obvious (see FerroSegers, 2003 for a recent method).
Inference
• Similar to the annual maxima method:
– Maximum likelihood is to be preferred
– Confidence bounds can be based on profile likelihood
– Model fit can be analyzed by P-P plots and Q-Q plots
or formal tests (similar to those presented earlier)
– Return levels/upper bounds can be estimated
• Our results for the flood data: sometimes slightly lower
return level estimators (reasons have to be analyzed) .
GPD fit: Vásárosnamény, water level
Density Plot
0.005
shape=-0.51,
0.001
0.002
0.003
the upper endpoint of
its 95% conf. int.: 1085 cm
0.000
f(x)
0.004
estimated
upper endpoint=940 cm
600
650
700
750
800
x
850
900
950
800
850
return level for year
30
900
950
1000
1050
Return level estimators by parts
of the dataset: Vásárosnamény
1950
1960
1970
1980
window:
50
1990
2000
Future
Our plans: to incorporate
most recent data into
the analyzis
Plans for the future
(engineers):
• to build temporal
reservoirs
• to utilise our results in
levy construction
So we may hope to
prevent such events
to happen again.
Some references
1.
Ferro, T. A.- Segers, J. (2003): Inference for clusters of
extreme values. Journal of Royal Statistical Soc. Ser. B.
65, p. 545-556.
2.
Kotz, S. – Nadarajah, S. (2000): Extreme Value
Distributions. Imperial College Press.
3.
Zempléni, A. (1996): Inference for Generalized Extreme
Value Distributions Journal of Applied Statistical
Science 4, p. 107-122.
4.
Zempléni, A. Goodness-of-fit tests in extreme value
theory. (In preparation.)