Transcript Extreme Value Analysis FISH 558 Decision Analysis in Natural Resource Management 12/4/2013

Extreme Value Analysis
FISH 558 Decision Analysis in Natural
Resource Management
12/4/2013
Noble Hendrix
QEDA Consulting LLC
Affiliate Faculty UW SAFS
Lecture Overview
• Motivating examples of extreme events
• Generalized Extreme Value
– Statistical Development
– Case Study: the white cliffs of Dover
• Generalized Pareto Distribution
– Statistical Development
– Case Study: whale strikes in SE Alaska
• Additional resources
2
Why should we care about extreme
events?
• They are rare by definition, so why spend
much time thinking about them?
• Often the consequences of the event have
significant impacts to the system – mortality,
colonization, episodic recruitment
• We tend to focus on averages, but extremes
may be more important in some situations.
• We may also be interested in estimating
extremes beyond what has been observed
3
0.010
0.005
0.000
density
0.015
0.020
Distribution of outcomes
|||||||||||||||
|||||||||||||||||||||| |||| | |
0
50
||
100
x
150
4
0.010
0.005
0.000
density
0.015
0.020
Distribution of outcomes
|||||||||||||||
|||||||||||||||||||||| |||| | |
0
50
||
100
x
150
5
0.010
0.005
0.000
density
0.015
0.020
Distribution of outcomes
|||||||||||||||
|||||||||||||||||||||| |||| | |
0
50
||
100
x
150
6
0.010
0.005
0.000
density
0.015
0.020
Distribution of outcomes
|||||||||||||||
|||||||||||||||||||||| |||| | |
0
50
||
100
x
150
7
0.010
0.005
0.000
density
0.015
0.020
Distribution of outcomes
|||||||||||||||
|||||||||||||||||||||| |||| | |
0
50
||
100
x
150
8
Motivation
100 year floodplain
9
Motivation
Surpassing the 100 year floodplain
• Road and home
construction
based on flood
frequency and
intensity i.e.,
100 year
floodplain
10
Motivation
Hurricanes
11
Financial Markets
12
Central Limit Theorem
Consider sequence of iid random variables, X1,
… Xn
We know that sum Sn = X1 + … + Xn, when
normalized lead to the CLT:
Statistical Foundations
13
Generalized Extreme Value
Fisher-Tippet Asymptotic Theorem
Define maxima of sequence of random
variables Mn = max(X1, …, Xn)
For normalized maxima, there is also a nondegenerate distribution H(x), which is a GEV
distribution
14
Generalized Extreme Value
Cumulative Density Function
u – location
s – scale
v - shape
15
Generalized Extreme Value
Variants of the GEV
Shape parameter v defines several
distributions:
Gumbel: v = 0
Weibull: v < 0
Fréchet: v > 0
16
Generalized Extreme Value
Shapes of GEV
0.1
0.6
shape = 0.0
shape = 0.5
shape = -0.5
0.0
Fréchet
0.4
0.2
Probability
Gumbel
0.2
0.3
0.8
0.4
Weibull
0.0
Density
Distribution function
1.0
Density function
0
1
2
3
x
4
5
0
1
2
3
4
5
x
17
Generalized Extreme Value
Applicability
Almost all common continuous distributions
converge on H(x) for some value of v
• Weibull – beta
• Gumbel – normal, lognormal, hyperbolic,
gamma, chi-squared
• Fréchet – Pareto, inverse gamma, Student t,
loggamma
18
Generalized Extreme Value
Minima
What about minima?
min(X1, …, Xn) = - max(-X1, … ,-Xn)
If H(x) is the limiting distribution for maxima,
then 1 – H(-x) is the limiting distribution for
minima, so can also be handled
19
Generalized Extreme Value
Estimation
Obtain data from an unknown distribution F
• Let’s assume that there is an extreme value
distribution Hv for some value of v
• The true distribution of the n-block
maximum Mn can be approximated for large
enough n with a GEV distribution H(x)
• Fit model to repeated observations of an nblock maximum, thus m blocks of size n
20
Generalized Extreme Value
Example - Data
Annual sea level height at Dover, Britain
between 1912 and 1992
4.0
3.6
3.2
Sea Level (m)
4.4
Dover, Britain
1920
1940
1960
Year
1980
21
Generalized Extreme Value
Example - Data
Annual sea level height at Dover, Britain
between 1912 and 1992
1.0
0.5
0.0
Density
1.5
Dover, Britain
3.2
3.4
3.6
3.8
4.0
4.2
4.4
4.6
Sea Level (m)
22
Generalized Extreme Value
R package evd
>
>
>
>
require(evd)
data(sealevel)
sl.no<-na.omit(sealevel[,1])
fgev(sl.no)
Call: fgev(x = sl.no)
Deviance: -5.022368
Estimates
loc
3.59252
scale
0.20195
Standard Errors
loc
scale
0.02642 0.01874
shape
-0.02107
shape
0.07730
23
Generalized Extreme Value
Diagnostics
Empirical
0.4
0.0
Model
3.5 4.0 4.5 5.0
Quantile Plot
0.8
Probability Plot
0.0
0.2
0.4
0.6
Empirical
0.8
1.0
3.4
3.6
3.8
4.0
4.2
4.4
Model
24
Return Level Plot
Return level – “how long
to wait on average until
see another event equal
to or more extreme”
Return Level Plot
3.5 4.0 4.5 5.0
Return Level
Generalized Extreme Value
0.2
1.0
5.0
20.0
Return Period
100.0
If H is the distribution of
the n-block maximum,
the k return level is the 1
– 1/k quantile of H
25
Generalized Extreme Value
4
0
-4
-2
log likelihood
2
0
-2
log likelihood
0
-2
-6
-4
-4
-8
-6
-6
log likelihood
2
2
4
4
Profile likelihood of parameters
3.50 3.55 3.60 3.65
location (u)
0.16
0.20
0.24
scale (s)
0.28
-0.2
0.0
0.2
shape (v)
26
Generalized Extreme Value
Limitations
• Limitations of the GEV:
– Used for block maxima, e.g., annual
precipitation, annual flow,
– Only 1 exceedance per block
– May ignore some important observations,
– Some go so far as to say it is a wasteful method!
(McNeil et al. 2005 Quantitative Risk Management, Princeton)
27
Generalized Pareto Distribution
GEV has largely been surpassed by another
method for extremes over a threshold
Pickands (1975) developed a model for
excesses y over threshold a
Pickands 1975 Annals of Stats 3:119
28
Generalized Pareto Distribution
a – threshold
b – scale
v - shape
29
Generalized Pareto Distribution
Shapes of GPD
Distribution function
0.8
0.6
0.4
shape = 0.0
shape = 0.5
shape = -0.5
0.2
Probability
0.6
0.4
0.2
Positive shape =
limitless loss
0.0
Density
0.8
1.0
Density function
0
1
2
x
3
4
0
1
2
3
4
x
30
Generalized Pareto Distribution
Applicability
For any continuous distributions that
converge on H(x) for some value of v, which
was most of the continuous distributions of
interest
The same distributions will converge on G(x)
as an excess distribution as the threshold a is
raised
31
Generalized Pareto Distribution
Estimation
Obtain data from an unknown distribution F
Calculate Yj = Xj – a for Na that exceed
threshold a
maximize log-likelihood:
32
Generalized Pareto Distribution
Threshold Estimation
Have an interesting problem:
• Need a value of threshold a that must be
high enough to satisfy the theoretical
assumptions
• Need enough data above the threshold a so
that the parameters are well estimated
• Use a sample mean residual life plot to help
identify a reasonable threshold value a
33
Generalized Pareto Distribution
Sample Mean Residual Life Plot
Let Y = X – a0. At threshold a0, if Y is GPD with parameters b and v then
E(Y) = b/(1 – v), v < 1
This is true for all thresholds ai > a0, but the scale parameter bi must be
appropriate to the threshold ai
E(X-ai| X > ai) = (bi + v*ai)/(1-v),
Thus E(X - a| X > a) is a linear function of a where GPD appropriate, so can
plot E(x-ai) (where x are our observed data) versus ai.
This is the sample mean residual life plot, and confidence intervals added
by assuming E(x-a) are approximately normally distributed
34
Generalized Pareto Distribution
Example - Data
Quantifying strike rates of whales in southeast
Alaska
35
Generalized Pareto Distribution
Distances to Whales
80
60
40
20
0
Whale distance metric
100
Minimum distances (i.e., D < 0) are where losses
occur, so transform distance D into a positive loss
metric, where value of 100 equates to D = 0
0
200
400
Observation number
600
800
36
Generalized Pareto Distribution
0.010
0.005
0.000
Density
0.015
Whale Distance Metric
0
20
40
60
80
100
Whale distance metric
37
Generalized Pareto Distribution
Threshold determination
30
20
10
0
Mean Excess
40
Mean Residual Life Plot
50
60
70
80
90
• Looking for
discontinuities in the
mean excess, E(x-ai), at
different threshold
values ai
• Identified value of 70
as the threshold
(equates to a distance
of 300m between
whales and ships)
Threshold
38
Generalized Pareto Distribution
100
60
library(POT)
mrlplot(w.metric, xlim =
c(50,90) )
tcplot(w.metric, u.range =
c(50, 90) )
20
Modified Scale
Threshold determination
50
60
70
80
90
-0.6
Discontinuity in scale and
shape estimates when
threshold a > 70
-1.0
Shape
-0.2
Threshold
Mean residual life plot
(previous slide) indicates a =
70
50
60
70
80
Threshold
90
39
Generalized Pareto Distribution
Estimation
> fitgpd(w.metric, thresh
= 70, est = "mle")
Estimator: MLE
Deviance: 974.4418
AIC: 978.4418
Standard Error Type:
observed
Standard Errors
scale
shape
1.53542 0.07452
Varying Threshold: FALSE
Threshold Call: 70
Number Above: 151
Proportion Above: 0.1946
Asymptotic Variance
Covariance
scale
shape
scale
2.357530 -0.106864
shape -0.106864
0.005553
Estimates
scale
shape
14.8380 -0.4706
Optimization Information
Convergence: successful
40
Generalized Pareto Distribution
Diagnostics
0.0
0.2
0.4
0.6
Empirical
0.8
1.0
90
80
70
------------------------------------ -------------- ------------------------ ------------------------------------- ---------------------------------- ----------------------------------------------------------- ------------------------------
100
QQ-plot
Empirical
0.4
0.0
Model
0.8
Probability plot
-- ------ - - - ------------ ---------------------------- --------------------------- -------------------- ---------------------------- ----------------------70
75
80
85
90
95
100
Model
41
Generalized Pareto Distribution
-978
-986
-982
log likelihood
-978
-982
-986
log likelihood
-974
-974
Likelihood profiles
10
12
14
16
scale (b)
18
20
-0.6
-0.5
-0.4
-0.3
-0.2
shape (v)
42
Generalized Pareto Distribution
0
-2
-6
-4
a = 60
a = 70
a = 80
-10
-8
relative log likelihood
-4
-6
-8
-10
-12
relative log likelihood
-2
0
Likelihood profiles with different thresholds
10
15
20
scale (b)
25
30
-0.8
-0.6
-0.4
-0.2
0.0
0.2
shape (v)
relative log likelihood - likelihood relative to maximum for that threshold value
43
Generalized Pareto Distribution
Empirical and Estimated
0.02
0.04
Empirical
GPD
0.00
Density
0.06
Density Plot
70
80
90
100 110 120
Comparison of empirical
(no observed strikes)
and GPD model
estimates for a = 70
• Since 2000, 2
confirmed strikes
• GPD provides better
characterization of
risk
Quantile
44
Generalized Pareto Distribution
Return Level
95 100
90
Return level – how many
encounters where whales
are less than 300m until a
strike?
75
80
85
Conditional return level of
approx. 500
70
Whale distance metric
Return Level Plot
1
5
50
500
Absolute return level of
approx. 2500 (1 in 5
encounters has an
encounter < 300m)
Return period (observations)
45
Summary: GEV and EVT
• Generalized Extreme Value (GEV) distribution
– Used for block maxima, e.g., maximum sea-level per
year
– Data loss due to only block maxima
• Generalized Pareto Distribution (GPD)
– Used for points over a threshold
– All exceedances above some limit are used
– Question about how to deal with selecting a
threshold value
46
Additional Resources
Books and Papers
Coles, S. 2001. An Introduction to Statistical Modelling of Extreme Values.
Springer Series in Statistics. London.
McNeil, A. J., Frey, R., & Embrechts, P. 2005. Quantitative risk management:
concepts, techniques, and tools. Princeton University Press.
Embrechts, P. 1997. Modelling extremal events: for insurance and finance
(Vol. 33). Springer.
Bayesian GPD Modeling
Coles, S. and L. Pericchi. 2003. Anticipating catastrophes through extreme
value modeling. Applied Statistics 52(4): 405–416.
Jagger. T. H. and J. B. Elsne 2004. Climatology models for extreme
hurricane winds near the United States. Journal of Climate 19: 3220-3236.
47
Additional Resources
Fitting models in R and BUGS
A few R packages
• Points over Threshold (POT)
• Extreme Value Distributions (evd)
• extRemes
• Quantitative Risk Management (QRM)
• evdbayes
BUGS
• OpenBUGS – GEV and GPD
• WinBUGS/JAGS – GPD with 1’s trick
48