VGAM (version 0.8-1)

gpd: Generalized Pareto Distribution Family Function

Description

Maximum likelihood estimation of the 2-parameter generalized Pareto distribution (GPD).

Usage

gpd(threshold = 0, lscale = "loge", lshape = "logoff", escale = list(),
    eshape = if (lshape == "logoff") list(offset = 0.5) else
             if (lshape == "elogit") list(min = -0.5, max = 0.5) else NULL,
    percentiles = c(90, 95), iscale = NULL, ishape = NULL,
    tolshape0 = 0.001, giveWarning = TRUE, method.init = 1, zero = 2)

Arguments

threshold

Numeric, values are recycled if necessary. The threshold value(s), called $\mu$ below.

lscale

Parameter link function for the scale parameter $\sigma$. See Links for more choices.

lshape

Parameter link function for the shape parameter $\xi$. See Links for more choices. The default constrains the parameter to be greater than $-0.5$ because if $\xi \leq -0.5$ then Fisher scoring does no

escale, eshape

Extra argument for the lscale and lshape arguments. See earg in Links for general information. For the shape parameter, if the

percentiles

Numeric vector of percentiles used for the fitted values. Values should be between 0 and 100. See the example below for illustration. However, if percentiles = NULL then the mean $\mu + \sigma / (1-\xi)$ is returned; this is only de

iscale, ishape

Numeric. Optional initial values for $\sigma$ and $\xi$. The default is to use method.init and compute a value internally for each parameter. Values of ishape should be between $-0.5$ and $1$. Values of iscale

tolshape0, giveWarning

Passed into dgpd when computing the log-likelihood.

method.init

Method of initialization, either 1 or 2. The first is the method of moments, and the second is a variant of this. If neither work, try assigning values to arguments ishape and/or iscale.

zero

An integer-valued vector specifying which linear/additive predictors are modelled as intercepts only. The value must be from the set {1,2} corresponding respectively to $\sigma$ and $\xi$. It is often a good idea for the $\sigma$ parameter only

Value

• An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm and vgam. However, for this VGAM family function, vglm is probably preferred over vgam when there is smoothing.

Warning

Fitting the GPD by maximum likelihood estimation can be numerically fraught. If $1 + \xi (y-\mu)/ \sigma \leq 0$ then some crude evasive action is taken but the estimation process can still fail. This is particularly the case if vgam with s is used. Then smoothing is best done with vglm with regression splines (bs or ns) because vglm implements half-stepsizing whereas vgam doesn't. Half-stepsizing helps handle the problem of straying outside the parameter space.

Details

The distribution function of the GPD can be written $$G(y) = 1 - [1 + \xi (y-\mu) / \sigma ]_{+}^{- 1/ \xi}$$ where $\mu$ is the location parameter (known, with value threshold), $\sigma > 0$ is the scale parameter, $\xi$ is the shape parameter, and $h_+ = \max(h,0)$. The function $1-G$ is known as the survivor function. The limit $\xi \rightarrow 0$ gives the shifted exponential as a special case: $$G(y) = 1 - \exp[-(y-\mu)/ \sigma].$$ The support is $y>\mu$ for $\xi>0$, and $\mu < y <\mu-\sigma \xi$="" for="" $\xi<0$.<="" p="">

Smith (1985) showed that if $\xi <= -0.5$="" then="" this="" is="" known="" as="" the="" nonregular="" case="" and="" problems="" difficulties="" can="" arise="" both="" theoretically="" numerically.="" for="" (regular)="" $\xi=""> -0.5$ the classical asymptotic theory of maximum likelihood estimators is applicable; this is the default.

Although for $\xi < -0.5$ the usual asymptotic properties do not apply, the maximum likelihood estimator generally exists and is superefficient for $-1 < \xi < -0.5$, so it is ``better'' than normal. When $\xi < -1$ the maximum likelihood estimator generally does not exist as it effectively becomes a two parameter problem.

The mean of $Y$ does not exist unless $\xi < 1$, and the variance does not exist unless $\xi < 0.5$. So if you want to fit a model with finite variance use lshape = "elogit".

References

Yee, T. W. and Stephenson, A. G. (2007) Vector generalized linear and additive extreme value models. Extremes, 10, 1--19.

Coles, S. (2001) An Introduction to Statistical Modeling of Extreme Values. London: Springer-Verlag.

Smith, R. L. (1985) Maximum likelihood estimation in a class of nonregular cases. Biometrika, 72, 67--90.

See Also

rgpd, meplot, gev, pareto1, vglm, vgam, s.

Examples

# Simulated data from an exponential distribution (xi=0)
threshold = 0.5
gdata = data.frame(y = threshold + rexp(n=3000, rate=2))
fit = vglm(y ~ 1, gpd(threshold=threshold), gdata, trace=TRUE)
head(fitted(fit))
coef(fit, matrix=TRUE)   # xi should be close to 0
Coef(fit)
summary(fit)

fit@extra$threshold  # Note the threshold is stored here

# Check the 90 percentile
ii = fit@y < fitted(fit)[1,"90%"]
100*table(ii)/sum(table(ii))   # Should be 90%

# Check the 95 percentile
ii = fit@y < fitted(fit)[1,"95%"]
100*table(ii)/sum(table(ii))   # Should be 95%

plot(fit@y, col="blue", las=1, main="Fitted 90% and 95% quantiles")
matlines(1:length(fit@y), fitted(fit), lty=2:3, lwd=2)


# Another example
threshold = 0
gdata = data.frame(x = runif(nn <- 2000))
xi = exp(-0.8)-0.5
gdata = transform(gdata, y = rgpd(nn, scale=exp(1+0.1*x), shape=xi))
fit = vglm(y ~ x, gpd(threshold), gdata, trace=TRUE)
coef(fit, matrix=TRUE)


# Nonparametric fits
gdata = transform(gdata, yy = y + rnorm(nn, sd=0.1))
fit1 = vgam(yy ~ s(x), gpd(threshold), gdata, trace=TRUE) # Not so recommended
par(mfrow=c(2,1))
plotvgam(fit1, se=TRUE, scol="blue")
fit2 = vglm(yy ~ bs(x), gpd(threshold), gdata, trace=TRUE) # More recommended
plotvgam(fit2, se=TRUE, scol="blue")