VALUE EXTREME THEORY
Richard L. Smith Department of Statistics and Operations Research University of North Carolina Chapel Hill, NC 27599-3260 [e-mail] AMS Committee on Probability and Statistics Short Course on Statistics of Extreme Events Phoenix, January 11, 2009
1
(From a presentation by Myles Allen)
2
3
4
OUTLINE OF TALK I. Extreme value theory Probability Models Estimation Diagnostics II. Example: North Atlantic Storms III. Example: European Heatwave IV. Example: Trends in Extreme Rainfall Events
5
I. EXTREME VALUE THEORY
6
EXTREME VALUE DISTRIBUTIONS
Suppose X1, X2, ..., are independent random variables with the same probability distribution, and let Mn = max(X1, ..., Xn). Under certain circumstances, it can be shown that there exist normalizing constants an > 0, bn such that Pr Mn - b n x an = F (anx + bn)n H(x).
The Three Types Theorem (Fisher-Tippett, Gnedenko) asserts that if nondegenerate H exists, it must be one of three types: H(x) = exp(-e-x), all x (Gumbel) 0 x<0 H(x) = (Frchet) e exp(-x-) x > 0 exp(-|x|) x < 0 H(x) = (Weibull) 1 x>0 In Frchet and Weibull, > 0. e
7
The three types may be combined into a single generalized extreme value (GEV) distribution: H(x) = exp - 1 +
x- +
-1/
,
(y+ = max(y, 0)) where is a location parameter, > 0 is a scale parameter and is a shape parameter. 0 corresponds to the Gumbel distribution, > 0 to the Frchet distribution with = 1/, < 0 e to the Weibull distribution with = -1/. > 0: "long-tailed" case, 1 - F (x) x-1/ , = 0: "exponential tail" < 0: "short-tailed" case, finite endpoint at - /
8
EXCEEDANCES OVER THRESHOLDS
Consider the distribution of X conditionally on exceeding some high threshold u: Fu(y) = F (u + y) - F (u) . 1 - F (u)
As u F = sup{x : F (x) < 1}, often find a limit Fu(y) G(y; u, ) where G is generalized Pareto distribution (GPD) y -1/ . G(y; , ) = 1 - 1 + +
9
The Generalized Pareto Distribution
y -1/ . G(y; , ) = 1 - 1 + + > 0: long-tailed (equivalent to usual Pareto distribution), tail like x-1/ , = 0: take limit as 0 to get y G(y; , 0) = 1 - exp - , i.e. exponential distribution with mean , < 0: finite upper endpoint at -/.
10
The Poisson-GPD model combines the GPD for the excesses over the threshold with a Poisson distribtion for the number of exceedances. Usually the mean of the Poisson distribution is taken to be per unit time.
11
POINT PROCESS APPROACH Homogeneous case: Exceedance y > u at time t has probability 1 y- 1+ +
-1/-1
exp - 1 +
u- +
-1/
dydt
12
Illustration of point process model.
13
Inhomogeneous case: Time-dependent threshold ut and parameters t, t, t Exceedance y > ut at time t has probability 1 y - t 1 + t t t +
-1/t -1
ut - t exp - 1 + t t +
-1/t
dydt
Estimation by maximum likelihood
14
ESTIMATION
GEV log likelihood: Yi - log 1 + - i i provided 1 + (Yi - )/ > 0 for each i. 1 +1 = -N log - Poisson-GPD model: 1 N Yi = N log - T - N log - 1 + log 1 + i=1 provided 1 + Yi/ > 0 for all i. The method of maximum likelihood states that we choose the parameters (, , ) or (, , ) to maximize . These can be calculated numerically on the computer.
15
Yi - 1+
-1/
DIAGNOSTICS
Gumbel plots QQ plots of residuals Mean excess plot Z and W plots
16
Gumbel plots Used as a diagnostic for Gumbel distribution with annual maxima data. Order data as Y1:N ... YN :N , then plot Yi:N against reduced value xi:N , xi:N = - log(- log pi:N ),
1 pi:N being the i'th plotting position, usually taken to be (i- 2 )/N .
A straight line is ideal. Curvature may indicate Frchet or Weibull e form. Also look for outliers.
17
Gumbel plots. (a) Annual maxima for River Nidd flow series. (b) Annual maximum temperatures in Ivigtut, Iceland.
18
QQ plots of residuals A second type of probability plot is drawn after fitting the model. Suppose Y1, ..., YN are IID observations whose common distribution function is G(y; ) depending on parameter vector . Sup^ pose has been estimated by , and let G-1(p; ) denote the inverse distribution function of G, written as a function of . A QQ (quantile-quantile) plot consists of first ordering the observations Y1:N ... YN :N , and then plotting Yi:N against the reduced value ^ xi:N = G-1(pi:N ; ),
1 where pi:N may be taken as (i - 2 )/N . If the model is a good fit, the plot should be roughly a straight line of unit slope through the origin.
Examples...
19
QQ plots for GPD, Nidd data. (a) u = 70. (b) u = 100.
20
Mean excess plot Idea: for a sequence of values of w, plot the mean excess over w against w itself. If the GPD is a good fit, the plot should be approximately a straight line. In practice, the actual plot is very jagged and therefore its "straightness" is difficult to assess. However, a Monte Carlo technique, assuming the GPD is valid throughout the range of the plot, can be used to assess this. Examples...
21
Mean excess over threshold plots for Nidd data, with Monte Carlo confidence bands, relative to threshold 70 (a) and 100 (b).
22
Z- and W-statistic plots Consider nonstationary model with t, t, t dependent on t. Z statistic based on intervals between exceedances Tk : Zk =
Tk Tk-1
u(s)ds,
u(s) = {1 + s(u - s)/s)}-1/s . W statistic based on excess values: if Yk is excess over threshold at time Tk , Wk = Tk Yk 1 log 1 + . Tk Tk + Tk (u - Tk )
Idea: if the model is exact, both Zk and Wk and i.i.d. exponential with mean 1. Can test this with various plots.
23
Diagnostic plots based on Z and W statistics for oil company insurance data (u = 5)
24
II. NORTH ATLANTIC CYCLONES
25
Data from HURDAT Maximum windspeeds in all North Atlantic Cyclones from 1851 2007
26
TROPICAL CYCLONES FOR THE NORTH ATLANTIC
160 Max Windspeed 40 1850 60 80 100 120 140
1900 Year
1950
2000
27
POT MODELS 19002007, u=102.5
Model Gumbel GEV GEV, lin GEV, quad GEV, cubic GEV, lin , lin log GEV, quad , lin log GEV, lin , quad log p 2 3 4 5 6 5 6 6 NLLH 847.8 843.8 834.7 833.4 829.8 828.0 826.8 827.2 NLLH+p 849.8 846.8 838.7 838.4 835.8 833.0 832.8 833.2
Fitted model: = 0 + 1t, log = 2 + 3t, const 0 102.5 2.4 1 0.0158 0.049 2 2.284 0.476 3 0.0075 0.0021 0.302 0.066
28
Estimate S.E.
TROPICAL CYCLONES FOR THE NORTH ATLANTIC
160 0.99 0.95 0.90 0.75 0.50 100 40 1900 60 80
Max Windspeed
120
140
1920
1940
1960 Year
1980
2000
2020
29
Diagnostic Plots for Atlantic Cyclones
(a) (b)
0.15 Correlation 5 z value 4 3 2 1 0 0 20 40 60 Time 80 100 Observed values 6 6 5 4 3 2 1 0 0 1 2 3 4 5 6 0.10 0.05 0.00
(c)
-0.05 -0.10 -0.15 2 4 6 Lag for z 8 10
Expected values for z
(d)
6 5 w value 4 3 2 1 0 0 20 40 60 Time 80 100 Observed values 6 5
(e)
0.15 0.10 Correlation 0.05 0.00
(f)
4 3 2 1 0 0 1 2 3 4 5 6
-0.05 -0.10 -0.15 2 4 6 Lag for w 8 10
Expected values for w
30
III. EUROPEAN HEATWAVE
31
Data: 5 model runs from CCSM 18712100, including anthropogenic forcing 2 model runs from UKMO 18612000, including anthropogenic forcing 1 model runs from UKMO 20012100, including anthropogenic forcing 2 control runs from CCSM, 230+500 years 2 control runs from UKMO, 341+81 years All model data have been calculated for the grid box from 3050o N, 10o W40o E, annual average temperatures over JuneAugust Expressed anomalies an from 19611990, similar to Stott, Stone and Allen (2004)
32
CLIMATE MODEL RUNS: ANOMALIES FROM 1961-1990
6
Temperature
CCSM ANTHRO. UKMO ANTHRO. CCSM CONTROL UKMO CONTROL
-2
0
2
4
1900
1950 Year
2000
2050
2100
33
Method: Fit POT models with various trend terms to the anthropogenic model runs, 18612010 Also fit trend-free model to control runs ( = 0.176, log = -1.068, = -0.068)
34
POT MODELS 18612010, u=1
Model Gumbel GEV GEV, lin GEV, quad GEV, cubic GEV, quart GEV, quad , lin log GEV, quad , quad log p 2 3 4 5 6 7 6 7 NLLH 349.6 348.6 315.5 288.1 287.7 285.1 287.9 287.0 NLLH+p 351.6 351.6 319.5 293.1 293.7 292.1 293.9 294.9
Fitted model: = 0 + 1t + 2t2, , const 0 0.187 0.335 1 0.030 0.0054 2 0.000215 0.00003 log 0.047 0.212 0.212 0.067
35
Estimate S.E.
CLIMATE MODEL RUNS: ANOMALIES FROM 1961-1990
6
Temperature
CCSM ANTHRO. UKMO ANTHRO. CCSM CONTROL UKMO CONTROL 99 % 90 % 99 % 90 %
-2
0
2
4
1900
1950 Year
2000
36
Diagnostic Plots for Temperatures (Control)
(a)
5 4 z value 3 2 1 0 200 600 Time 1000 Observed values 5 Correlation 4 3 2 1 0 0 1 2 3 4 5 -0.2 2 4 6 Lag for z 8 10
(b)
0.2 0.1 0.0 -0.1
(c)
Expected values for z
(d)
(e)
0.2 Correlation 0.1 0.0 -0.1 -0.2 0 1 2 3 4 5 2 4
(f)
w value
3 2 1 0 200 600 Time 1000
Observed values
4
4 3 2 1 0
6 Lag for w
8
10
Expected values for w
37
Diagnostic Plots for Temperatures (Anthropogenic)
(a) (b)
0.2 Observed values 4 z value 3 2 1 0 0 200 400 600 800 Time 4 3 2 1 0 0 1 2 3 4 5 Correlation 0.1 0.0 -0.1 -0.2 2 4 6 Lag for z 8 10
(c)
Expected values for z
(d)
6 5 w value 4 3 2 1 0 0 200 400 600 800 Time Observed values 6 5
(e)
0.2 Correlation 0.1 0.0 -0.1 -0.2 0 1 2 3 4 5 2 4
(f)
4 3 2 1 0
6 Lag for w
8
10
Expected values for w
38
We now estimate the probabilities of crossing various thresholds in 2003. Express answer as N=1/(exceedance probability) Threshold 2.3: N=3024 (control), N=29.1 (anthropogenic) Threshold 2.6: N=14759 (control), N=83.2 (anthropogenic)
39
IV. TREND IN PRECIPITATION EXTREMES
(joint work with Amy Grady and Gabi Hegerl) During the past decade, there has been extensive research by climatologists documenting increases in the levels of extreme precipitation, but in observational and model-generated data. With a few exceptions (papers by Katz, Zwiers and co-authors) this literature have not made use of the extreme value distributions and related constructs There are however a few papers by statisticians that have explored the possibility of using more advanced extreme value methods (e.g. Cooley, Naveau and Nychka, to appear JASA; Sang and Gelfand, submitted) This discussion uses extreme value methodology to look for trends
40
DATA SOURCES
NCDC Rain Gauge Data (Groisman 2000) Daily precipitation from 5873 stations Select 19701999 as period of study 90% data coverage provision -- 4939 stations meet that NCAR-CCSM climate model runs 20 41 grid cells of side 1.4o 19701999 and 20702099 (A1B scenario) PRISM data 1405 621 grid, side 4km Elevations Mean annual precipitation 19701997
41
EXTREME VALUES METHODOLOGY
The essential idea is to fit a probability model to the exceedances over a high threshold at each of 5000 data sites, and then to combine data across sites using spatial statistics. The model at each site is based on the generalized extreme value distribution, interpreted as an approximate tail probability in the right hand tail of the distribution. y- Pr{Y y} t 1 + for large y, + Here x+ = max(x, 0), t is a time increment (here 1 day based on a time unit of 1 year) and the parameters , , represent the location, scale and shape of the distribution. In particular, when > 0 the marginal distributions have a Pareto (power-law) tail with power -1/.
42
-1/
TEMPORAL AND SPATIAL DEPENDENCE
Here, we make two extensions of the basic methodology. First, the parameters , , are allowed to be time-dependent through covariates. This allows a very flexible approach to seasonality, and we can also introduce linear trend terms to examine changes in the extreme value distribution over the time period of the study. The second extension is spatial smoothing: after estimating the 25-year return value at each site, we smooth the results across sites by a technique similar to kriging. We allow for spatial nonstationarity by dividing the US into 19 overlapping boxes, and interpolating across the boundaries.
43
Continental USA divided into 19 regions
44
Map of 25-year return values (cm.) for the years 19701999
45
Root mean square prediction errors for map of 25-year return values for 19701999
46
Ratios of return values in 1999 to those in 1970, using a statistical model that assumes a linear trend in the GEV model parameters
47
A B C D E F G H I J
Change 0.01 0.07 0.11 0.05 0.13 0.00 0.01 0.08 0.07 0.05
RMSPE .03 .03 .01 .01 .02 .02 .02 .01 .01 .01
K L M N O P Q R S
Change 0.08 0.07 0.07 0.02 0.01 0.07 0.07 0.15 0.14
RMSPE .01 .02 .02 .03 .02 .01 .01 .02 .02
For each grid box, we show the mean change in log 25-year return value (1970 to 1999) and the corresponding standard error (RMSPE) Stars indicate significance at 5%, 1%, 0.1%. 14 of 19 regions are statistically significant increasing: the remaining five are all in western states
48
We can use the same statistical methods to project future changes by using data from climate models. Here we use data from CCSM, the climate model run at NCAR.
49
Return value map for CCSM data (cm.): 19701999
50
Return value map for CCSM data (cm.): 20702099
51
Estimated ratios of 25-year return values for 20702099 to those of 19701999, based on CCSM data, A1B scenario
52
The climate model data show clear evidence of an increase in 25year return values over the next 100 years, as much as doubling in some places.
53
A caveat... Although the overall increase in observed precipitation extremes is similar to that stated by other authors, the spatial pattern is completely different. There are various possible explanations, including different methods of spatial aggregation and different treatments of seasonal effects. Even when the same methods are applied to CCSM data over 19701999, the results are different.
54
Extreme value model with trend: ratio of 25-year return value in 1999 to 25-year return value in 1970, based on CCSM data
55
CONCLUSIONS 1. Focus on N -year return values -- strong historical tradition for this measure of extremes (we took N = 25 here) 2. Seasonal variation of extreme value parameters is a critical feature of this analysis 3. Overall significant increase over 19701999 except for parts of western states -- average increase across continental US is 7% 4. Projections to 20702099 show further strong increases but note caveat based on point 5 5. But... based on CCSM data there is a completely different spatial pattern and no overall increase -- still leaves some doubt as to overall interpretation.
56
FURTHER READING Finkenstadt, B. and Rootzn, H. (editors) (2003), Extreme Vale ues in Finance, Telecommunications and the Environment. Chapman and Hall/CRC Press, London. (See http://www.stat.unc.edu/postscript/rs/semstatrls.pdf) Coles, S.G. (2001), An Introduction to Statistical Modeling of Extreme Values. Springer Verlag, New York.
57
THANK YOU FOR YOUR ATTENTION!
58