gpd: Generalized Pareto distribution modelling

In texmex: Threshold exceedences and multivariate extremes

Description

Likelihood based modelling and inference for the generalized Pareto distribution, possibly with explanatory variables.

Usage

gpd(y, data, ...)

## Default S3 method:
gpd(y, data, th, qu, phi = ~1, xi = ~1, penalty = "gaussian",
    prior = "gaussian", method = "optimize", cov="observed", start = NULL, 
    priorParameters = NULL, maxit = 10000, trace = NULL,
    iter = 10500, burn = 500, thin = 1, jump.cov, jump.const, verbose = TRUE,...)

## S3 method for class 'gpd'
print(x, digits=max(3, getOption("digits") - 3), ...)
## S3 method for class 'gpd'
summary(object, nsim=1000, alpha=0.05, ...)
## S3 method for class 'gpd'
show(x, digits=max(3, getOption("digits") - 3), ...)

## S3 method for class 'gpd'
plot(x, main=rep(NULL, 4), xlab=rep(NULL, 4), nsim=1000, alpha=0.05, ...)

## S3 method for class 'gpd'
AIC(object, ..., k=2)

## S3 method for class 'bgpd'
print(x, print.seed=FALSE, ...)
## S3 method for class 'bgpd'
plot(x, which.plots=1:3, density.adjust=2, print.seed=FALSE, ...)

Arguments

y	Either a numeric vector or the name of a variable in data.
data	A data frame containing y and any covariates.
th	The threshold for y, exceedances above which will be used to fit the GPD upper tail model.
qu	An alternative to th, a probability defined such that quantile(y,qu) equals th.
phi	Formula for the log of the scale parameter. Defaults to phi = ~ 1 - i.e. no covariates.
xi	Formula for the shape parameter. Defaults to xi = ~ 1 - i.e. no covariates.
penalty	How to penalize the likelhood. Currently, either “none”, “gaussian” or “lasso” are the only allowed arguments. If penalty is "gaussian" or "lasso" then the parameters for the penalization are specified through the prior argument, see below.
prior	If method = "optimize", just an alternative way of specifying the pentalty, and only one or neither of penalty and prior should be given. If method = "simulate", prior must be “gaussian” because no other prior distributions have been implemented.
method	Should be either “optimize” (the default) or “simulate”. The first letter or various abbreviations will do. If “optimize” is used, the (penalized) likelihood is directly optimized using optim and point estimates (either ML or MAP estimates) are returned with other information. If “simulate”, a Metropolis algorithm is used to simulate from the joint posterior distribution of the parameters.
cov	How to compute the covariance matrix of the parameters. Defaults to cov = "observed" in which case the observed information matrix is used, as given in Appendix A of Davison and Hinkley. The only other option is cov = "numeric" in which case a numerical approximation of the Hessian is used (see the help for optim). In some cases, particularly with small samples, the numerical approximation can be quite different from the closed form (cov="observed") result, and the value derived from the observed information should be preferred. However, in either case, since the underlying log-likelihood may be far from quadratic for small samples, the resulting estimates of standard errors are liable to approximate poorly the true standard errors. Also see the comments in the Details section, below.
start	Starting values for the parameters, to be passed to optim. If not provided, an exponential distribution (shape = 0) is assumed as the starting point.
priorParameters	A list with two components. The first should be a vector of means, the second should be a covariance matrix if the penalty/prior is "gaussian" or "quadratic" and a diagonal precision matrix if the penalty/prior is "lasso", "L1" or "Laplace". If method = "simulate" then these represent the parameters in the Gaussian prior distribution. If method = "optimize" then these represent the parameters in the penalty function. If not supplied: all default prior means are zero; all default prior variances are 10^4; all covariances are zero.
maxit	The number of iterations allowed in optim.
trace	Whether or not to print progress to screen. If method = "optimize", the argument is passed into optim – see the help for that function. If method = "simulate", the argument determines at how many steps of the Markov chain the function should tell the user, and in this case it defaults to trace = 1000.
iter	Number of simulations to generate under method = "simulate".
burn	The number of initial steps to be discarded.
thin	The degree of thinning of the resulting Markov chains. Defaults to 1 (no thinning). A value of 0.5 (for example) would result in every other step being discarded.
jump.cov	Covariance matrix for proposal distribution of Metropolis algorithm. This is scaled by jump.const.
jump.const	Control parameter for the Metropolis algorithm.
verbose	Whether or not to print progress to screen. Defaults to verbose=TRUE.
x, object	Object of class gpd, bgpd, summary.gpd or summary.bgpd returned by gpd or summary.gpd.
digits	Number of digits for printing.
main	In plot method for class gpd, titles for diagnostic plots. Should be a vector of length 4, with values corresponding to the character strings to appear on the titles of the pp- qq- return level and density estimate plots respectively.
xlab	As for main but labels for x-axes rather than titles.
nsim	In plot and summary methods for class gpd. The number of replicates to be simulated to produce the simulated tolerance intervals. Defaults to nsim = 1000

.

alpha	In plot and summary methods for class gpd. A (1 - alpha)% simulation envelope is produced. Defaults to alpha = 0.05
k	Constant used in calculation of AIC=-2*loglik + k*p, defaults to k=2.
print.seed	Whether or not to print the seed used in the simulations, or to annotate the plots with it. Defaults to print.seed=FALSE.
which.plots	In plot method for class bgpd. Which plots to produce. Option 1 gives kernel density estimates, 2 gives traces of the Markov chains with superimposed cumulative means, 3 gives autocorrelation functions. Defaults to which.plots=1:3.
density.adjust	In plot method for class bgpd. Passed into density. Controls the amount of smoothing of the kernel density estimate. Defaults to density.adjust=2.
...	Further arguments to be passed to methods.

Details

We use the following parameterisation of the GPD:

P(X ≤ x) = 1 - (1 + ξ x / σ)^(-1/ξ)

for x ≥ 0 and 1 + ξ x / σ ≥ 0. The scale parameter is sigma (σ) and the shape parameter is xi (ξ).

Working with the log of the scale parameter improves the stability of computations, making a quadratic penalty more appropriate and enabling the inclusion of covariates in the model for the scale parameter, which must remain positive. We therefore work with φ=log(σ). All specification of priors or penalty functions refer to φ rather than σ. A quadratic penalty can be thought of as a Gaussian prior distribution, whence the terminology of the function.

Parameters of the GPD fitted to excesses above threshold th are estimated by using penalized maximum likelihood (method="optimize"), or by simulating from the posterior distribution of the model parameters using a Metropolis algorithm (method="simulate"). In the latter case, start is used as a starting value for the Metropolis algorithm; in its absence, the maximum penalized likelhood point estimates are computed and used.

When a summary or plot is performed, a pointwise (1 - alpha)% tolerance envelope is simulated, based on quantiles of the fitted model. Since the ordered observations will be correlated, if any observation is outside the envelope, it is likely that a chain of observations will be outside the envelope. Therefore, if the number outside the envelope is a little more than alpha%, that does not immediately imply a serious shortcoming of the fitted model.

When method = "optimize", the plot function produces diagnostic plots for the fitted generalized Pareto model. These differ depending on whether or not there are covariates in the model. If there are no covariates then the diagnostic plots are PP- and QQ-plots, a return level plot (produced by plotrl.gpd) and a histogram of the data with superimposed generalized Pareto density estimate. These are all calculated using the data on the original scale. If there are covariates in the model then the diagnostics consist of PP- and QQ- plots calculated by using the model residuals (which will be standard exponential devaiates under the GPD model), and plots of residuals versus fitted model parameters.

The PP- and QQ-plots show simulated pointwise tolerance intervals. The region is a (1 - alpha)% region based on nsim simulated samples.

When method = "simulate" the plot function produces diagnostic plots for the Markov chains used to simulate from the posterior distributions for the model parameters. If the chains have converged on the posterior distributions, the trace plots should look like "fat hairy caterpillars" and their cumulative means should converge rapidly. Moreover, the autocorrelation functions should converge quickly to zero.

When method = "simulate" the print and summary functions give posterior means and standard deviations. Posterior means are also returned by the coef method. Depending on what you want to do and what the posterior distributions look like (use plot method) you might want to work with quantiles of the posterior distributions instead of relying on standard errors.

Value

If method = "optimize", an object of class gpd:

convergence	Output from optim relating to whether or not the optimizer converged.
message	A message telling the user whether or not convergence was achieved.
threshold	The threshold of the data above which the GPD model was fit.
penalty	The type of penalty function used, if any.
coefficients	The parameter estimates as computed under maximum likelihood or maximum penalized likelihood.
rate	The proportion of observations above the threshold.
call	The call to gpd that produced the object.
y	The response data above the threshold for fitting.
X.phi	The design matrix for the log of the scale parameter.
X.xi	The design matrix for the shape parameter.
priorParameters	See above.
data	The original data (above and below the threshold for fitting).
residuals	Data above the threshold for fitting after transformation to standard exponential scale by using the fitted GPD.
loglik	The value of the optimized log-likelihood.
cov	The estimated covariance of the parameters in the model.
se	The estimated standard errors of the parameters in the model.

If method = "simulate", an object of class bgpd:

call	The call to gpd that produced the object.
threshold	The threshold above which the model was fit.
map	The point estimates found by maximum penalized likelihood and which were used as the starting point for the Markov chain. This is of class gpd and methods for this class (such as resid and plot) may be useful.
burn	The number of steps of the Markov chain that are to be treated as the burn-in and not used in inferences.
thin	The degree of thinning used.
chains	The entire Markov chain generated by the Metropolis algorithm.
y	The response data above the threshold for fitting.
data	The original data (above and below the threshold for fitting).
X.phi	The design matrix for the log of the scale parameter.
X.xi	The design matrix for the log of the shape parameter.
acceptance	The proportion of proposals that were accepted by the Metropolis algorithm.
seed	The seed used by the random number generator.
param	The remainder of the chain after deleting the burn-in and applying any thinning.

There are summary, plot, print and coefficients methods available for these classes.

Note

When there are estimated values of xi <= -0.5, the regularity conditions of the likelihood break down and inference based on approximate standard errors cannot be performed. In this case, the most fruitful approach to inference appears to be by the bootstrap. It might be possible to simulate from the posterior, but finding a good proposal distribution might be difficult and you should take care to get an acceptance rate that is reasonably high (around 40% when there are no covariates, lower otherwise).

Author(s)

Janet E. Heffernan, Harry Southworth. Some of the internal code is based on the gpd.fit function in the ismev package and is due to Stuart Coles.

References

A. C. Davison and R. L. Smith, Models for exceedances over high thresholds, Journal of the Royal Statistical Society B, 53, 393 – 442, 1990

See Also

rl.gpd, predict.gpd, gpd.declustered.

Examples

  x <- rnorm(1000)
  mod <- gpd(x, qu = 0.7)
  mod
  par(mfrow=c(2, 2))
  plot(mod)
  
  x <- runif(100,-0.2,0.2)
  data <- data.frame(x=x,y=rgpd(100,sigma=exp(3 + 2*x),xi=x))
  mod <- gpd(y, data, phi = ~x, xi = ~x, th = 0)
  plot(mod)
  
# Following lines commented out to keep CRAN robots happy
#  mod <- gpd(x, qu=.7, method="sim")
#  mod
#  par(mfrow=c(3, 2))
#  plot(mod)

texmex documentation built on May 2, 2019, 4:56 p.m.