stability: Generalized Pareto parameter estimate stability

Description Usage Arguments Details Value References See Also Examples

View source: R/parameter_stability.R

Description

Uses maximum likelihood estimation to fit a Generalized Pareto (GP) model to threshold excesses over a range of thresholds. The threshold excesses are treated as independent and identically distributed (i.i.d.) observations. The resulting estimates and confidence intervals can be plotted, using plot.stability, to produce a crude graphical diagnostic for threshold choice.

Usage

1
2
stability(data, u_vec, prof = FALSE, conf = 95, mult = 1:2,
  plot_prof = FALSE, ...)

Arguments

data

A numeric vector of observations.

u_vec

A numeric vector of thresholds to be applied to the data. Any duplicated values will be removed. These could be set at sample quantiles of data using quantile.

prof

A logical scalar. Whether to calculate confidence intervals for the GP shape parameter ξ based on the profile-likelihood for ξ or using the MLE plus or minus a multiple of the estimated standard error (SE) of the MLE. The intervals produced by the former may be better but they take longer to calculate. Default: FALSE.

conf

A numeric scalar in (0, 100). Confidence level for the confidence intervals. Default: 95%.

mult

A numeric vector of length 2. The range of values over which the profile log-likelihood for ξ is calculated is (MLE - mult[1] c SE, MLE + mult[2] c SE), where MLE and SE are the MLE and estimated standard error of ξ and c is the constant for which this interval gives an approximate 100conf% level confidence interval for ξ when mult = c(1, 1). The default, mult = c(1, 2), works well in most cases. If the routine fails because the range of ξ is not sufficiently wide then the relevant components of mult should be increased.

plot_prof

A logical scalar. Only relevant if prof = TRUE. If plot_prof = TRUE then the profile log-likelihood is plotted for each threshold. If FALSE then nothing is plotted.

...

Further (optional) arguments to be passed to the optim function for the optimizations on which the profile-likelihood for xi is based.

Details

For each threshold in u_vec a GP model is fitted by maximum likelihood estimation to the threshold excesses, i.e. the amounts by which the data exceed that threshold. The MLEs of the GP shape parameter $ξ$ and approximate conf% confidence intervals for ξ are stored for plotting (by plot.stability) to produce a simple graphical diagnostic to inform threshold selection. This plot is used to choose a threshold above which the underlying GP shape parameter may be approximately constant. See Chapter 4 of Coles (2001). See also the vignette "Introducing threshr".

Value

An object (list) of class "stability" with components:

ests

MLEs of the GP shape parameter ξ.

ses

Estimated SEs of the MLEs of ξ.

lower

Lower limit of 100conf% confidence intervals for ξ.

upper

Upper limit of 100conf% confidence intervals for ξ.

nexc

The number of threshold excesses.

u_vec

The thresholds supplied by the user.

u_ps

The approximate sample quantiles to which the thresholds in u_vec correspond.

data

The input data.

conf

The input conf.

Each of these components is a numeric vector of length length(u_vec).

References

Coles, S. G. (2001) An Introduction to Statistical Modeling of Extreme Values, Springer-Verlag, London. http://dx.doi.org/10.1007/978-1-4471-3675-0_3

See Also

ithresh for threshold selection in the i.i.d. case based on leave-one-out cross-validation.

plot.stability for the S3 plot method for objects of class stability.

quantile.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
# Set a vector of thresholds
u_vec_gom <- quantile(gom, probs = seq(0, 0.95, by = 0.05))

# Symmetric confidence intervals
gom_stab <- stability(data = gom, u_vec = u_vec_gom)
plot(gom_stab)

# Profile-likelihood-based confidence intervals
gom_stab <- stability(data = gom, u_vec = u_vec_gom, prof = TRUE)
plot(gom_stab)

threshr documentation built on Sept. 4, 2017, 9:03 a.m.