welcome to part 4 of our short introduction to extreme value analysis using the
extRemes package in R.
In the recent posts about the block maxima method and the threshold excess method, we have simply assumed that all assumptions for extreme value analysis are fulfilled. However, this is most likely not the case when dealing with environmental variables. Especially the assumption of stationarity may be violated in many cases. Against the background of global climate change, it is likely that there is a a considerable trend in the time series of meteorological or other environmental variables. Of course, this trend has to be incorporated into the analysis, as the resulting return levels change over time.
The following code shows a short practical example of fitting a generalized pareto distribution to a time series of precipitation data using the
extRemes package in R. The sample data set features precipitation data of Bärnkopf (Lower Austria) from 1971 through 2013 and is provided by the hydrographic services of Austria via via eHYD.
# load packages library(extRemes) library(ismev) library(RCurl) library(ggplot2) library(Kendall) # get data from eHYD link <- getURL("http://ehyd.gv.at/eHYD/MessstellenExtraData/nlv?id=107540&file=2") df <- read.csv(text = link, sep=";", dec=",", skip=23, 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))) # derive AMS for maximum temperature ams <- data.frame() ams <- do.call("rbind", by(df, format(df$Date, "%Y"), function(x) x[which.max(x$Precip), ])) # check stationarity within AMS MannKendall(ams$Precip) # Mann-Kendall trend test summary(lm(Precip~Date, data=ams)) p <- ggplot(ams, aes(x = Date, y = Precip)) + geom_point() + geom_line() p + stat_smooth(method = "lm", formula = y ~ x, size = 1)
Both the results of the fitted linear model and the graphical impression given by the plot indicate an upward trend of annual maximum precipitation amounts. The Mann-Kendall trend test returns a very small p-value, confirming this trend. Hence, trend correction has to be applied in order to account for changing return levels over time.
# maximum likelihood estimation mle.trend <- fevd(ams$Precip, ams, location.fun=~Date, method = "MLE", type="GEV") rltrend <- return.level(mle.trend, conf = 0.05, return.period= c(2,5,10,20,50,100)) # return level plot plot(mle.trend, type="rl", ylim=c(0,200), main="Return Level Plot for Bärnkopf w/ MLE")
Reference: Coles, Stuart (2001). An Introduction to Statistical Modeling of Extreme Values. Springer Series in Statistics. London: Springer.