fit_GPD function

Parameter Estimators of the Generalized Pareto Distribution

Method-of-moments estimator, probability weighted moments estimator, log-likelihood and maximum-likelihood estimator for the parameters of the generalized Pareto distribution (GPD).

fit_GPD_MOM(x)
fit_GPD_PWM(x)

logLik_GPD(param, x)
fit_GPD_MLE(x, init = c("PWM", "MOM", "shape0"),
            estimate.cov = TRUE, control = list(), ...)

Arguments

• x: numeric vector of data. In the peaks-over-threshold method, these are the excesses (exceedances minus threshold).

• param: numeric(2) containing the value of the shape $xi$ (a real) and scale $beta$

  (positive real) parameters of the GPD in this order.

• init: character string specifying the method for computing initial values. Can also be numeric(2)

  for directly providing $xi$ and $beta$.

• estimate.cov: logical indicating whether the asymptotic covariance matrix of the parameter estimators is to be estimated (inverse of observed Fisher information (negative Hessian of log-likelihood evaluated at MLE)) and standard errors for the estimators of $xi$ and $beta$

  returned, too.

• control: list; passed to the underlying optim().

• ...: additional arguments passed to the underlying optim().

Returns

fit_GEV_MOM() and fit_GEV_PWM() return a numeric(3) giving the parameter estimates for the GPD.

logLik_GPD() computes the log-likelihood of the GPD (-Inf if not admissible).

fit_GPD_MLE() returns the return object of optim()

and, appended, the estimated asymptotic covariance matrix and standard errors of the parameter estimators, if estimate.cov.

Details

fit_GPD_MOM() computes the method-of-moments (MOM) estimator.

fit_GPD_PWM() computes the probability weighted moments (PWM) estimator of Hosking and Wallis (1987); see also Landwehr et al. (1979).

fit_GPD_MLE() uses, as default, fit_GPD_PWM() for computing initial values. The former requires the data x

to be non-negative and adjusts $beta$ if $xi$ is negative, so that the log-likelihood at the initial value should be finite.

Author(s)

Marius Hofert

References

McNeil, A. J., Frey, R. and Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques, Tools. Princeton University Press.

Hosking, J. R. M. and Wallis, J. R. (1987). Parameter and Quantile Estimation for the Generalized Pareto Distribution. Technometrics 29 (3), 339--349.

Landwehr, J. M., Matalas, N. C. and Wallis, J. R. (1979). Estimation of Parameters and Quantiles of Wakeby Distributions. Water Resourches Research 15 (6), 1361--1379.

Examples

## Simulate some data
xi <- 0.5
beta <- 3
n <- 1000
set.seed(271)
X <- rGPD(n, shape = xi, scale = beta)

## Fitting via matching moments
(fit.MOM <- fit_GPD_MOM(X))
stopifnot(all.equal(fit.MOM[["shape"]], xi,   tol = 0.52),
          all.equal(fit.MOM[["scale"]], beta, tol = 0.24))

## Fitting via PWMs
(fit.PWM <- fit_GPD_PWM(X))
stopifnot(all.equal(fit.PWM[["shape"]], xi,   tol = 0.2),
          all.equal(fit.PWM[["scale"]], beta, tol = 0.12))

## Fitting via MLE
(fit.MLE <- fit_GPD_MLE(X))
(est <- fit.MLE$par) # estimates of xi, mu, sigma
stopifnot(all.equal(est[["shape"]], xi,   tol = 0.12),
          all.equal(est[["scale"]], beta, tol = 0.11))
fit.MLE$SE # estimated asymp. variances of MLEs = std. errors of MLEs

## Plot the log-likelihood in the shape parameter xi for fixed
## scale beta (fixed as generated)
xi. <- seq(-0.1, 0.8, length.out = 65)
logLik <- sapply(xi., function(xi..) logLik_GPD(c(xi.., beta), x = X))
plot(xi., logLik, type = "l", xlab = expression(xi),
     ylab = expression("GPD log-likelihood for fixed"~beta))

## Plot the profile likelihood for these xi's
## (with an initial interval for the nuisance parameter beta such that
##  logLik_GPD() is finite)
pLL <- sapply(xi., function(xi..) {
    ## Choose beta interval for optimize()
    int <- if(xi.. >= 0) {
               ## Method-of-Moment estimator
               mu.hat <- mean(X)
               sig2.hat <- var(X)
               shape.hat <- (1-mu.hat^2/sig2.hat)/2
               scale.hat <- mu.hat*(1-shape.hat)
               ## log-likelihood always fine for xi.. >= 0 for all beta
               c(1e-8, 2 * scale.hat)
           } else { # xi.. < 0
               ## Make sure logLik_GPD() is finite at endpoints of int
               mx <- max(X)
               -xi.. * mx * c(1.01, 100) # -xi * max(X) * scaling
               ## Note: for shapes xi.. closer to 0, the upper scaling factor
               ##       needs to be chosen sufficiently large in order
               ##       for optimize() to find an optimum (not just the
               ##       upper end point). Try it with '2' instead of '100'.
           }
    ## Optimization
    optimize(function(nuis) logLik_GPD(c(xi.., nuis), x = X),
             interval = int, maximum = TRUE)$maximum
})
plot(xi., pLL, type = "l", xlab = expression(xi),
     ylab = "GPD profile log-likelihood")