Description Usage Arguments Details Value Note Author(s) References See Also Examples
Likelihood based modelling and inference for the generalized Pareto distribution, possibly with explanatory variables.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25  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, ...)

y 
Either a numeric vector or the name of a variable in 
data 
A data frame containing 
th 
The threshold for 
qu 
An alternative to 
phi 
Formula for the log of the scale parameter. Defaults to

xi 
Formula for the shape parameter. Defaults to

penalty 
How to penalize the likelhood. Currently, either
“none”, “gaussian” or “lasso” are the only allowed arguments.
If 
prior 
If 
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 
cov 
How to compute the covariance matrix of the parameters. Defaults to

start 
Starting values for the parameters, to be passed to 
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 
maxit 
The number of iterations allowed in 
trace 
Whether or not to print progress to screen. If 
iter 
Number of simulations to generate under 
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 
Control parameter for the Metropolis algorithm. 
verbose 
Whether or not to print progress to screen. Defaults to

x, object 
Object of class 
digits 
Number of digits for printing. 
main 
In 
xlab 
As for 
nsim 
In 
.
alpha 
In 
k 
Constant used in calculation of AIC=2*loglik + k*p, defaults to 
print.seed 
Whether or not to print the seed used in the simulations,
or to annotate the plots with it. Defaults to

which.plots 
In 
density.adjust 
In 
... 
Further arguments to be passed to methods. 
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 QQplots, 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 QQplots 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.
If method = "optimize"
, an object of class gpd
:
convergence 
Output from 
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 
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 loglikelihood. 
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 
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 
burn 
The number of steps of the Markov chain that are to be treated as the burnin 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 burnin and applying any thinning. 
There are summary, plot, print and coefficients methods available for these classes.
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).
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.
A. C. Davison and R. L. Smith, Models for exceedances over high thresholds, Journal of the Royal Statistical Society B, 53, 393 – 442, 1990
rl.gpd
, predict.gpd
, gpd.declustered
.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16  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)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.