Hey there,

welcome to part 2 of our short introduction to extreme value analysis using the `extRemes`

package in R.

As some of you might know, there are two common approaches for practical extreme value analysis. Today, we will focus on the first of these two approaches, called the block maxima method.

This approach for modelling extremes of a (time) series of observations is based on the utilization of maximum or minimum values of these observations within a certain sequence of constant length. For a sufficiently large number *n* of established blocks, the resulting peak values of these *n* blocks of equal length can be used for fitting a suitable distribution to these data. While the block size is basically freely selectable, a trade-off has to be made between bias (small blocks) and variance (large blocks). Usually, the length of the sequences is often chosen to correspond to a certain familiar time period, in most cases a year. The resulting vector of annual maxima (or minima, respectively) is calles “Annual Maxima (Minima) Series” or simply AMS.

According to the Fisher–Tippett–Gnedenko theorem, the distribution of block maxima can be approximated by a generalized extreme value distribution.

The following code shows a short practical example of fitting a generalized extreme value distribution to a time series of precipitation data using the `extRemes`

package in R. The sample data set features precipitation data of Oberwang (Salzburg, Austria) from 1981 through 2013 and is is provided by the hydrographic services of Austria via eHYD.

# load packages library(extRemes) library(RCurl) # get data from eHYD link <- getURL("http://ehyd.gv.at/eHYD/MessstellenExtraData/nlv?id=105700&file=2") df <- read.csv(text = link, sep=";", dec=",", skip=21, header=FALSE) df[df == "Lücke"] = NA colnames(df) <- c("Date","Precip") df$Date <- as.Date(df$Date, format = "%d.%m.%Y") df$Precip <- as.numeric(gsub(",", ".", as.character(df$Precip))) head(df) # derive AMS for maximum precipitation ams <- data.frame() ams <- do.call("rbind", by(df, format(df$Date, "%Y"), function(x) x[which.max(x$Precip), ])) amsval <- ams$Precip # generalized maximum-likelihood fitting of the GEV distribution mlefit <- fevd(amsval, method = "GMLE", type="GEV") plot(mlefit) # diagnostic plots # return levels: rlmle <- return.level(mlefit, conf = 0.05, return.period= c(2,5,10,20,50,100)) # fitting of GEV distribution based on L-moments estimation lmomfit <- fevd(amsval, method = "Lmoments", type="GEV") plot(lmomfit) # diagnostic plots # return levels: rlmom <- return.level(lmomfit, conf = 0.05, return.period= c(2,5,10,20,50,100)) # return level plots par(mfcol=c(1,2)) # return level plot w/ MLE plot(mlefit, type="rl", main="Return Level Plot for Oberwang w/ MLE", ylim=c(0,200), pch=16) loc <- as.numeric(return.level(mlefit, conf = 0.05,return.period=100)) segments(100, 0, 100, loc, col= 'midnightblue',lty=6) segments(0.01,loc,100, loc, col='midnightblue', lty=6) # return level plot w/ LMOM plot(lmomfit, type="rl", main="Return Level Plot for Oberwang w/ L-Moments", ylim=c(0,200)) loc <- as.numeric(return.level(lmomfit, conf = 0.05,return.period=100)) segments(100, 0, 100, loc, col= 'midnightblue',lty=6) segments(0.01,loc,100, loc, col='midnightblue', lty=6) # comparison of return levels out <- t(data.frame(mle=as.numeric(rlmle), lmom=as.numeric(rlmom))) colnames(out) <- c(2,5,10,20,50,100) round(out,1)

In this case, L-Moments estimation reacts more sensitive to the influence of the outlier, which is quite unusual. In most cases, L-moments estimation is more robust than maximum likelihood estimation, but our example apparently is the exception that proves the rule. In addition, I have made the observation that maximum likelihood estimation works more reliable in other R packages in some cases (e.g. `fExtremes`

, `ismev`

, `evd`

, `evir`

).

*Reference: Coles, Stuart (2001). An Introduction to Statistical Modeling of Extreme Values. Springer Series in Statistics. London: Springer.*

## 20 Comments

You can post comments in this post.

i want a help of you

Muhammad Arif 1 year ago

Sure. Can you specify your request?

Matthias 1 year ago

plz send me the r-coding for fitting generalized pareto distribution,estimating the parameters estimation by mle method

Muhammad Arif 1 year ago

There is another post about fitting a generalized pareto distribution to a partial duration series in R using the extRemes package. Both maximum likelihood estimation and L-moments estimaiton are covered there:

Introduction to Extreme Value Analysis in R – Part 3: Peak over Threshold Approach

Matthias 1 year ago

my research topic is “non stationary frequency analysis of extreme precipitation”.but i did not know my track. how i will proceed my job?

please help me.

thanks for your kindness!

Muhammad Arif 1 year ago

I also wrote a post about dealing with non-stationarity, you can find further information there:

Introduction to Extreme Value Analysis in R – Part 4: Dealing with trends

The

`fevd()`

function from the package`extRemes`

is quite powerful and well documented, I recommend to have a look at its help page.Matthias 1 year ago

Hi I am currently doing a thesis on Garch-EVT approach to estimate Value at Risk and Expected Shortfall. Can you send me the R procedure for this approach? Thanks

Gaby 7 months ago

Hi Gaby,

I have not yet worked with GARCH-EVT-copulas, if that’s what you intend to do. I am pretty sure that there is no specific R package targeted at this topic. So you probably have to split your procedure into several steps and deal with each step individually. Stack overflow is most likely a very good place to get some information and see the workflow of other people working with this topic. For a start, also take a look at

`fGarch`

(which is part of the`Rmetrics`

-suite, just as the excellent package`fExtremes`

) or`rugarch`

for GARCH modelling.In addition, you might want to have a look at the package

`GEVStableGarch`

which seemingly employs maximum likelihood extimation of GARCH models with a GEV distribution, maybe this is sufficient for your needs.Regards,

Matthias

Matthias 7 months ago

Thanks for your reply.

I’ll just be working with garch and evt. The procedure consists of applying a garch model to the return series and then extract the residuals from the model. Then we apply evt on the residuals extracted from the garch model. However I am having difficulties selecting the threshold from the residuals so as to apply EVT-POT or block maxima method to it. Can you help me?

Regards

Gaby

Gaby 7 months ago

I’m shocked that I found this info so easily.

Denver 4 months ago

Hi

Your website is very useful and make this topic easy to follow

How did you change your points on your return level plots to green?

misha 3 months ago

Hi,

unfortunately, the plotting methods for objects of type

`fevd`

are – albeit comprehensive – hardly customizable. Since I didn’t like the default plotting, I modified the plotting functions`plot.fevd.mle`

and`plot.fevd.lmoments`

to support better visualization options within the`extRemes`

framework.Regards,

Matthias

Matthias 3 months ago

Matthias,

i am still learning phase of R.

Would you be able to provide me with some guidance on how to get this done. Thanks

misha 3 months ago

Hi Misha,

I have written the following function for creating a simple return level plot (base graphics) from

`fevd`

objects that have been fitted of type GEV.You can control both the color of points and the lines of the fitted function as well as the number of ticks on your y-axis.

Please note that this is a very simple approach, I'm mainly re-using code from

`plot.fevd.mle()`

.In order to create a similar plot for the threshold excess approach, you would have to make slight adjustments to parameters of the GP and to the plotting positions of the extreme values:

Using

`rlplot_gev(mlefit, pointcolor = "darkblue")`

on the code of my post should yield return level plots similar to the default`plot(mlefit, type = "rl", xlim = c(1, 200), ylim = c(25, 225))`

Regards,

Matthias

Matthias 3 months ago

Hey!

I am currently writing my Bachelor Thesis on evds.

Would it be okay for you if I used your code as an example?

Lukas 3 months ago

Hi Lukas,

yes, sure. However, if you have enough time it might be worth to have a deeper look at the packages presented above.

This post is really meant as an introduction to the topic.

Regards,

Matthias

Matthias 3 months ago

When modelling the effects of extreme events (say of variable X) on volatility of prices of a commodity (say Y), is possible to fit GEV to X and then insert it as an exogenous variable in mean or volatility equation? If not, how can I approach this issue?

Waiguru 4 weeks ago

Hi,

I am not really sure if I understand what you are intending to do.

What is the question you are trying to answer?

Modelling the influence of X on Y sounds more like a standard regression model in the first place.

Matthias 4 weeks ago

Hello!

I have a question:

How can I use R and the Gumbel distribution to predict discharge on specific return periods?

And how do I know the return period?

Thank you!

Abeer Haddad 3 weeks ago

Hi,

basically, you can stick to my example code above.

If you specifically want to fit a Gumbel distribution, simply set the

`type`

argument in the`fevd()`

function to`type="Gumbel"`

.Return levels can be obtained by using the function

`return.level()`

.Regards,

Matthias

Matthias 1 week ago

## Post A Reply