GevModelling: Generalized Extreme Value Modelling

Description

A collection and description functions to compute the generalized extreme value distribution and to estimate it parameters. The functions compute density, distribution function, quantile function and generate random deviates for the GEV, for the Frechet, Gumbel, and Weibull distributions. To model the GEV three types of approaches for parameter estimation are provided: Maximum likelihood estimation, probability weighted moment method, and estimation by the MDA approach. MDA includes functions for the Pickands, Einmal-Decker-deHaan, and Hill estimators together with several plot variants. The GEV distribution functions are: ll{ dgev density of the GEV Distribution, pgev probability function of the GEV Distribution, qgev quantile function of the GEV Distribution, rgev random variates from the GEV Distribution. [dpqr]frechet Frechet Distribution, [dpqr]gumbel Gumbel Distribution, [dpqr]weibull Weibull Distribution, [dpqr]evd an alternative call for the GEV Distribution. } The GEV modelling functions are: ll{ gevSim generates data from the GEV, gevFit fits empirical or simulated data to the distribution, print print method for a fitted GEV object, plot plot method for a fitted GEV object, summary summary method for a fitted GEV object, gevrlevelPlot k-block return level with confidence intervals. } Maximum Domain of Attraction Estimators: ll{ hillPlot shape parameter and Hill estimate of the tail index, shaparmPlot variation of shape parameter with tail depth. }

Usage

dgev(x, xi = 1, mu = 0, sigma = 1, log = FALSE)
pgev(q, xi = 1, mu = 0, sigma = 1, lower.tail = TRUE)
qgev(p, xi = 1, mu = 0, sigma = 1, lower.tail = TRUE)
rgev(n, xi = 1, mu = 0, sigma = 1)
devd(x, loc = 0, scale = 1, shape = 0, log = FALSE)
pevd(q, loc = 0, scale = 1, shape = 0, lower.tail = TRUE)
qevd(p, loc = 0, scale = 1, shape = 0, lower.tail = TRUE)
revd(n, loc = 0, scale = 1, shape = 0)

gevSim(model = list(shape = 0.25, location = 0, scale = 1), n = 1000)
gevFit(x, block = NA, type = c("mle", "pwm"), gumbel = FALSE, ...)
## S3 method for class 'gevFit':
print(x, \dots)
## S3 method for class 'gevFit':
plot(x, which = "all", \dots)
## S3 method for class 'gevFit':
summary(object, doplot = TRUE, which = "all", \dots)
gevrlevelPlot(object, k.blocks = 20, add = FALSE, ...)

hillPlot(x, option = c("alpha", "xi", "quantile"), start = 15,
    end = NA, reverse = FALSE, p = NA, ci = 0.95, autoscale = TRUE, 
    labels = TRUE, ...) 
shaparmPlot(x, revert = FALSE, standardize = FALSE, tails = 0.01*(1:10), 
    doplot = c(FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE), 
    which = c(TRUE, TRUE, TRUE), doprint = TRUE, both.tails = TRUE, 
    xi.range = c(-0.5, 1.5), alpha.range = c(0, 10))

Arguments

add

[gevrlevelPlot] - whether the return level should be added graphically to a time series plot; if FALSE a graph of the profile likelihood curve showing the return level and its confidence interval is produced.

alpha.range, xi.range

[saparmPlot] - plotting ranges.

autoscale

[hillPlot] - whether or not plot should be automatically scaled; if not, xlim and ylim graphical parameters may be entered.

block

[gevFit] - the block size. Only used if method="mle" is selected. A numeric value is interpreted as the number of data values in each successive block. All the data is used, so the last block may not contain

both.tails

[shaparmPlot] - a logical, decides whether or not both tails should be investigated. By default TRUE. If FALSE only the lower tail will be investigated.

[hillPlot] - probability for asymptotic confidence band; for no confidence band set ci to zero.

doplot

a logical. Should the results be plotted? [shaparmPlot] - a vector of logicals of the same lengths as tails defining for wich tail depths plots should be created, by default plots will be generated for a tail depth of 5

doprint

[shaparmPlot] - a logical, decides whether or not for all tail depths the result for the shape parameter 1/alpha should be printed.

gumbel

[gevFit] - a logical, by default FALSE. To fit a Gumbel model with fixed shape=0 set gumbel=TRUE.

k.blocks

[gevrlevelPlot] - specifies the particular return level to be estimated; default set arbitrarily to 20.

labels

[hillPlot] - whether or not axes should be labelled.

loc, scale, shape

loc is the location parameter, scale the scale parameter, and shape is the shape parameter. The default values are loc=0, scale=1, and shape=0

log

a logical, if TRUE, the log density is returned.

lower.tail

a logical, if TRUE, the default, then probabilities are P[X <= x]<="" code="">, otherwise, P[X > x].

model

[gevSim] - a list with components shape, location and scale giving the parameters of the GEV distribution. By default the shape parameter has the value 0.25, the location is zero and the

[gevSim] - number of generated data points, an integer value. [rgev][revd] - the number of observations.

object

[summary][grlevelPlot] - a fitted object of class "gevFit".

option

[hillPlot] - whether alpha, xi (1/alpha) or quantile (a quantile estimate) should be plotted.

[qgev][qevs] - a numeric vector of probabilities. [hillPlot] - probability required when option quantile is chosen.

[pgev][pevs] - a numeric vector of quantiles.

reverse

[hillPlot] - whether plot is to be by increasing threshold, TRUE, or increasing number of order statistics FALSE.

revert

[shaparmPlot] - a logical value, by default FALSE, if set to TRUE the sign of the vector will be reverted: x = -x.

start, end

[hillPlot] - lowest and highest number of order statistics at which to plot a point.

standardize

[shaparmPlot] - a logical value, by default FALSE, if set to TRUE the vector x will be standardized: x = (x-mean(x))/sqrt(var(x)).

tails

[shaparmPlot] - a numeric vector of tail depths to be considered; by default ten values ranging from 0.1 to 1.0 in steps of 0.1 corresponding to values ranging from 1 to 10 percent.

type

a character string denoting the type of parameter estimation, either by maximum likelihood estimation "mle", the default value, or by the probability weighted moment menthod "pwm".

which

[shaparmPlot] - a vector of 3 logicals indicating which plots from the three methods will be created. The first entry decides for the Pickands estimator, the second for the Hill estimator, and the last for the Deckers-Einmahl

[dgev][devd] - a numeric vector of quantiles. [gevFit] - data vector. In the case of method="mle" the interpretation depends on the value of block: if no block size is specified then data are interpreted as blo

xi, mu, sigma

xi is the shape parameter, mu the location parameter, and sigma is the scale parameter. The default values are xi=1, mu=0, and sigma=1.

...

[gevFit] - control parameters optionally passed to the optimization function. Parameters for the optimization function are passed to components of the control argument of optim. [hillPlot]

Value

d* returns the density, p* returns the probability, q* returns the quantiles, and r* generates random variates. All values are numeric vectors. gevSim returns a vector of data points from the simulated series. gevFit returns an object of class gev describing the fit. print.summary prints a report of the parameter fit. summary performs diagnostic analysis. The method provides two different residual plots for assessing the fitted GEV model. gevrlevelPlot returns a vector containing the lower 95% bound of the confidence interval, the estimated return level and the upper 95% bound. hillPlot displays a plot. shaparmPlot returns a list with one or two entries, depending on the selection of the input variable both.tails. The two entries upper and lower determine the position of the tail. Each of the two variables is again a list with entries pickands, hill, and dehaan. If one of the three methods will be discarded the printout will display zeroes.

Details

Generalized Extreme Value Distribution: Computes density, distribution function, quantile function and generates random variates for the Generalized Extreme Value Distribution, GEV, for the Frechet, Gumbel, and Weibull distributions. Parameter Estimation: gevFit estimates the parameters either by the probability weighted moment method, method="pwm" or by maximum log likelihood estimation method="mle". As a limiting case the Gumbel distribution can be selected. The summary method produces diagnostic plots for fitted GEV or Gumbel models. Methods: print.gev, plot.gev and summary.gev are print, plot, and summary methods for a fitted object of class gev. Concerning the summary method, the data are converted to unit exponentially distributed residuals under null hypothesis that GEV fits. Two diagnostics for iid exponential data are offered. The plot method provides two different residual plots for assessing the fitted GEV model. Two diagnostics for iid exponential data are offered. Return Level Plot: gevrlevelPlot calculates and plots the k-block return level and 95% confidence interval based on a GEV model for block maxima, where k is specified by the user. The k-block return level is that level exceeded once every k blocks, on average. The GEV likelihood is reparameterized in terms of the unknown return level and profile likelihood arguments are used to construct a confidence interval. Hill Plot: The function hillPlot investigates the shape parameter and plots the Hill estimate of the tail index of heavy-tailed data, or of an associated quantile estimate. This plot is usually calculated from the alpha perspective. For a generalized Pareto analysis of heavy-tailed data using the gpdFit function, it helps to plot the Hill estimates for xi. Shape Parameter Plot: The function shaparmPlot investigates the shape parameter and plots for the upper and lower tails the shape parameter as a function of the taildepth. Three approaches are considered, the Pickands estimator, the Hill estimator, and the Decker-Einmal-deHaan estimator.

References

Coles S. (2001); Introduction to Statistical Modelling of Extreme Values, Springer. Embrechts, P., Klueppelberg, C., Mikosch, T. (1997); Modelling Extremal Events, Springer.

Examples

Run this code

## SOURCE("fExtremes.52A-GevModelling")

## *gev  -
   xmpExtremes("Start: GEV Frechet >")
   # Create and plot 1000 GEV/Frechet distributed rdv:
   par(mfrow = c(3, 3))
   r = rgev(n = 1000, xi = 1)
   plot(r, type = "l", main = "GEV/Frechet Series")
   ## Plot empirical density and compare with true density:
   ## Omit values greater than 500 from plot
   hist(r[r<10], n = 25, probability = TRUE, xlab = "r", 
     xlim = c(-5, 5), ylim = c(0, 1.1), main = "Density")
   x = seq(-5, 5, by=0.01)
   lines(x, dgev(x, xi = 1), col = 2)
   ## Plot df and compare with true df:
   plot(sort(r), (1:length(r)/length(r)), 
     xlim = c(-3, 6), ylim = c(0, 1.1),
     cex = 0.5, ylab = "p", xlab = "q", main = "Probability")
   q = seq(-5,5, by=0.1)
   lines(q, pgev(q, xi=1), col=2)
## Compute quantiles, a test:
   qgev(pgev(seq(-5, 5, 0.25), xi = 1), xi = 1) 
   
## *gev -
   xmpExtremes("Next: GEV Gumbel >")
   # Create and plot 1000 Gumbel distributed rdv:
   ##> r = rgev(n = 1000, xi = 0)
   ##> plot(r, type = "l", main = "Gumbel Series")
   ## Plot empirical density and compare with true density:
   ##>hist(r[abs(r)<10], nclass = 25, freq = FALSE, xlab = "r", 
   ##>   xlim = c(-5,5), ylim = c(0,1.1), main = "Density")
   ##>x = seq(-5, 5, by = 0.01)
   ##>lines(x, dgev(x, xi = 0), col=2)
   ## Plot df and compare with true df:
   ##>plot(sort(r), (1:length(r)/length(r)), 
   ##>  xlim = c(-3, 6), ylim = c(0, 1.1),
   ##>   cex=0.5, ylab = "p", xlab="q", main="Probability")
   ##>q = seq(-5, 5, by = 0.1)
   ##>lines(q, pgev(q, xi = 0), col = 2)
   ## Compute quantiles, a test:
   ##>qgev(pgev(seq(-5, 5, 0.25), xi = 0), xi = 0)   

## *gev -
   xmpExtremes("Next: GEV Weibull >")
   # Create and plot 1000 Weibull distributed rdv:
   r = rgev(n = 1000, xi = -1)
   plot(r, type = "l", main = "Weibull Series")
   ## Plot empirical density and compare with true density:
   hist(r[abs(r)<10], nclass = 25, freq = FALSE, xlab = "r", 
     xlim=c(-5,5), ylim=c(0,1.1), main="Density")
   x = seq(-5, 5, by=0.01)
   lines(x, dgev(x, xi = -1), col = 2)
   ## Plot df and compare with true df:
   plot(sort(r), (1:length(r)/length(r)), 
     xlim = c(-3, 6), ylim = c(0, 1.1),
     cex = 0.5, ylab = "p", xlab = "q", main = "Probability")
   q=seq(-5, 5, by = 0.1)
   lines(q, pgev(q, xi = -1), col = 2)
   ## Compute quantiles, a test:
   qgev(pgev(seq(-5, 5, 0.25), xi = -1), xi = -1)   

## gevSim -
## gevFit -
   # Simulate GEV Data:
   xmpExtremes("Start: Simulte GEV Sample >")
   # Use default length n=1000
   x = gevSim(model = list(shape = 0.25, location =0 , scale = 1))
   # Fit GEV Data by Probability Weighted Moments:
   fit = gevFit(x, type = "pwm") 
   print(fit)
   # Summarize Results:
   par(mfcol = c(3, 2))
   summary(fit)
   
## gevFit -
   # Fit GEV Data by Max Log Likelihood Method:
   xmpExtremes("Next: Estimate Parameters >")
   fit = gevFit(x, type = "mle") 
   print(fit)
   # Summarize Results:
   summary(fit) 
    
## gevSim -
## gevFit -
   # Simulate Gumbel Data:
   xmpExtremes("Next: Simulte Gumbel Sample >")
   # Use default length n=1000
   ##> x = gevSim(model = list(shape = 0, location = 0, scale = 1))
   # Fit Gumbel Data by Probability Weighted Moments:
   ##> fit = gevFit(x, type = "pwm", gumbel = TRUE) 
   ##> print(fit)
   # Summarize Results:
   ##> par(mfcol = c(3, 2))
   ##> summary(fit)  
   
## Fit Gumbel Data by Max Log Likelihood Method:
   xmpExtremes("Next: Estimate Parameters >")
   ##> fit = gevFit(x, type = "mle", gumbel = TRUE) 
   ##> print(fit)
   # Summarize Results:
   ##> summary(fit)     
   ##> xmpExtremes("Press any key to continue >") 
   
## Return levels based on GEV Fit:
   # BMW Stock Data:
   xmpExtremes("Next: Compute BMW Return Levels >")
   par(mfrow = c(2, 1))
   data(bmw)
   # Fit GEV to monthly Block Maxima:
   fit = gevFit(-bmw, block = "month") 
   # Calculate the 40 month return level
   gevrlevelPlot(fit, k.block = 40, main = "BMW: Return Levels")   
   
## Return levels based on GEV Fit:
   xmpExtremes("Next: Compute SIEMENS Return Levels >")
   # Siemens Stock Data:
   data(siemens)
   # Fit GEV to monthly Block Maxima:
   fit = gevFit(-siemens, block = "month") 
   # Calculate the 40 month return level
   gevrlevelPlot(fit, k.block = 40, main = "SIEMENS: Return Levels")
   
## Interactive Plot:
   ##> par(mfrow = c(1, 1), ask = TRUE)
   ##> plot(fit)
   
## hillPlot -
   xmpExtremes("Start: Hill Estimator >")
   # Hill plot of heavy-tailed Danish fire insurance data 
   # and BMW stock data for estimated 0.999 quantile
   par(mfrow = c(2, 2))
   data(bmw)
   hillPlot(bmw)
   hillPlot(bmw, option = "quantile", end = 500, p = 0.999)
   data(danish)
   hillPlot(danish)
   hillPlot(danish, option = "quantile", end = 500, p = 0.999)
   
## shaparmPlot -
   xmpExtremes("Next: Shape Parameter Plots >")
   par(mfcol = c(3, 2), cex = 0.6)
   data(bmw)
   shaparmPlot(bmw)

## shaparmPlot -
   xmpExtremes("Next: Simulated Frechet Data >")
   par(mfcol = c(3, 2), cex = 0.6)
   set.seed(4711)
   x = rgev(10000, xi = 1/4)
   shaparmPlot(x, revert = TRUE, both.tails = FALSE)
   lines(c(0.01, 0.1), c(4, 4), col = "steelblue3") # True Value

Run the code above in your browser using DataLab