return.level: Calculation of the return levels.

View source: R/extremes.R

return.levelR Documentation

Calculation of the return levels.

Description

Calculate arbitrary return levels and their error estimates for generalized extreme value (GEV) and generalized Pareto (GP) distributions.

Usage

return.level(x, return.period = 100, error.estimation = c("none", "MC",
  "MLE", "bootstrap"), model = c("gev", "gpd"),
  monte.carlo.sample.size = 100, bootstrap.sample.size = 100,
  threshold = NULL, total.length = NULL,
  thresholded.time.series = NULL, extreme.type = c("max", "min"),
  silent = FALSE, mc.cores = NULL)

Arguments

x

Either of class numeric, climex.fit.gev, climex.fit.gpd or a list of those objects.

error.estimation

Method for calculating the standard errors of the fitted results. The errors of the GEV/GP parameters will be calculated as the square roots of the diagonal elements of the inverse of the hessian matrix. The latter will be evaluated at the maximum likelihood estimates (MLE) of the GEV/GP parameters.

For all three methods of estimating the fitting errors of the return levels underlying series of threshold exceedances or block maxima is required. In case the user supplies numerical values to specify GEV/GP parameters and not the output of either the fit.gev or fit.gpd function no error estimation for the return level can be performed.

MLE: The standard error of the return level is calculated using the Delta method and the maximum likelihood estimates of the GEV/GP parameters. Note: For positive shape parameters bigger than 0.3 this approach tends to highly overestimates the errors of the return levels.

MC: Alternative one can use a Monte Carlo method for which monte.carlo.sample.size samples of the same size as x will be drawn from a GEV/GP distribution constituted by the obtained MLE of the GEV/GP parameters of x. The standard error is then calculated via the square of the variance of all fitted GEV/GP parameters and calculated return levels. Note: In its essence this approach is not an estimation of the error involved in fitting the time series to a GEV/GP distribution. It is rather the mean error of fitting a GEV/GP-distribution with the same length and parameters as estimated ones.

bootstrap: Using this option the provided time series x will be sampled with replacement bootstrap.sample.size times and with the same length as the original time series. The standard errors of the GEV/GP parameters and return levels of all those sampled series is calculated and returned as an estimate of the fitting error. Note: Since the data is (hopefully) GEV/GP-distributed, such a sampling has to be treated with a lot of care.

Sometimes the inversion of the hessian fails (since the are some NaN in the hessian) when calculating the error estimates using the maximum likelihood approach (MLE) (which is also the reason why the ismev package occasionally does not work). In such cases the Monte Carlo (MC) method is used as a fallback.

none skips the calculation of the error. Default = "MLE".

model

Determining whether to calculate the initial parameters of the GEV or GP function. Default = "gev".

monte.carlo.sample.size

Number of samples used to obtain the Monte Carlo estimate of the standard error of the fitting. Default = 100.

bootstrap.sample.size

Number of samples with replacements to drawn from the original series x in order to determine the standard errors for the GP parameters and return levels. Default = 100.

threshold

Optional threshold for the GP model. If present it will be added to the return level to produce a value which fits to underlying time series. Default = NULL.

total.length

Total number of observations in the time series the exceedance were obtained from (before! applying the threshold). This argument is needed to calculate the standard error of the return level via the delta method of the MLE in the GP model. Default = NULL.

thresholded.time.series

Time series used with fit.gpd on which already a threshold (the one supplied here as well) was applied. Necessary to transform the return level for numerical input and the GP model from m-th observation return level to annual return level. If omitted the return level will be per observation. Default = NULL.

extreme.type

String indicating whether to calculate the quantiles at the right ("max") or left ("min") side of the PDF of the series. Default = "max".

silent

Throws an warning whenever the "gpd" model is used and the thresholded.time.series is not supplied. Since this can be annoying one can also disable it. Default = FALSE.

mc.cores

A numerical input specifying the number of cores to use for the multi core application of the function (see detectCores). This functionality is only available if the input is a list of different objects. If NULL, the function will be calculated classically with only one core. Default = NULL.

Details

Uses the rlevd function at its core (a port from the extRemes package) but also can handle the outputs of the fit.gev and fit.gpd functions, is capable of calculating numerous return levels at once, and also calculates the errors of the return levels. For the errors the ML fit is using the option hessian=TRUE (if not done already) or uses a Monte Carlo based approach. If no fitting object is provided, no errors will be calculated.

This function is also capable of working with lists of fit or numerical objects.

Value

A list containing the estimates "upper.limit" and their standard errors "error" if the input was a fitting object or a numerical input. If, on the other hand, the input was a list of such objects, the output is a list of the aforementioned list.

See Also

Other extremes: block.list, block.xts, block, decluster.list, decluster.xts, decluster, extremal.index, gev.density, gpd.density, qevd, return.level.climex.fit.gev, return.level.climex.fit.gpd, return.level.list, return.level.numeric, revd, rlevd, threshold.list, threshold.xts, threshold, upper.limit.climex.fit.gev, upper.limit.climex.fit.gpd, upper.limit.list, upper.limit.numeric, upper.limit

Examples

fit.results <- fit.gev( block( anomalies( temp.potsdam ) ) )
return.level( fit.results, return.period = c( 10, 50, 100 ),
              error.estimation = "MLE" )

theGreatWhiteShark/climex documentation built on July 13, 2022, 9:11 a.m.