qrmtools (version 0.0-10)

mean_excess: Mean Excess

Description

Sample mean excess function, mean excess function of a GPD and sample mean excess plot.

Usage

mean_excess_np(x, omit = 3)
mean_excess_plot(x, omit = 3,
                 xlab = "Threshold", ylab = "Mean excess over threshold", ...)
mean_excess_GPD(x, shape, scale)

Arguments

x

mean_excess_GPD()

numeric vector of evaluation points of the mean excess function of the GPD.

otherwise

numeric vector of data.

omit

number \(\ge 1\) of unique last observations to be omitted from the sorted data (as mean excess plot becomes unreliable for these observations as thresholds).

xlab

x-axis label.

ylab

y-axis label.

…

additional arguments passed to the underlying plot().

shape

GPD shape parameter \(\xi\).

scale

GPD scale parameter \(\beta\).

Value

mean_excess_np() returns a two-column matrix giving the sorted data without the omit-largest unique values (first column) and the corresponding values of the sample mean excess function (second column). It is mainly used in mean_excess_plot().

mean_excess_plot() returns invisible().

mean_excess_GPD() returns the mean excess function of a generalized Pareto distribution evaluated at x.

Details

Mean excess plots can be used in the peaks-over-threshold method for choosing a threshold. To this end, one chooses the smallest threshold above which the mean excess plot is roughly linear.

Examples

# NOT RUN {
## (Sample) mean excess function
data(fire)
ME <- mean_excess_np(fire)
stopifnot(dim(ME) == c(2164, 2),
          all.equal(ME[nrow(ME),], c(65.707491, 121.066231),
                    check.attributes = FALSE))

## A 'manual' (sample) mean excess plot
plot(ME, xlab = "Threshold", ylab = "Mean excess over threshold")

## (Sample) mean excess plot
mean_excess_plot(fire)
## => Any value in [10, 20] seems reasonable here as threshold choice
##    (one would probably go with 10 to benefit from a larger sample size).

## With mean excess functions of two fitted GPDs overlaid
u <- c(10, 20) # thresholds
fit <- lapply(u, function(u.) fit_GPD_MLE(fire[fire > u.] - u.))
q <- lapply(u, function(u.) seq(u., ME[nrow(ME),"x"], length.out = 129))
MEF.GPD <- lapply(1:2, function(k)
    mean_excess_GPD(q[[k]]-u[k], shape = fit[[k]]$par[["shape"]],
                    scale = fit[[k]]$par[["scale"]]))
mean_excess_plot(fire, ylim = range(ME, unlist(MEF.GPD)))
col <- c("royalblue3", "maroon3")
for(k in 1:2) lines(q[[k]], MEF.GPD[[k]], col = col[k])
legend("bottomright", col = rev(col), lty = rep(1, length(u)), bty = "n",
       legend = as.expression(sapply(rev(seq_along(u)),
       function(k) substitute("Threshold choice"~~u==u., list(u. = u[k])))))
# }