Learn R Programming
QRM (version 0.4-31)

POT: Peaks-over-Threshold Method

Description

Functions for fitting, analysing and risk measures according to POT/GPD

Usage

fit.GPD(data, threshold = NA, nextremes = NA, type = c("ml", "pwm"),
        information = c("observed", "expected"),
        optfunc = c("optim", "nlminb"), verbose = TRUE, ...)
plotTail(object, ppoints.gpd = ppoints(256), main = "Estimated tail probabilities",
         xlab = "Exceedances x", ylab = expression(1-hat(F)[n](x)), ...)
showRM(object, alpha, RM = c("VaR", "ES"),
       like.num = 64, ppoints.gpd = ppoints(256),
       xlab = "Exceedances x", ylab = expression(1-hat(F)[n](x)),
       legend.pos = "topright", pre.0.4.9=FALSE, ...)
findthreshold(data, ne)
MEplot(data, omit = 3., main = "Mean-Excess Plot", xlab = "Threshold",
       ylab = "Mean Excess", ...)
xiplot(data, models = 30., start = 15., end = 500., reverse = TRUE,
       ci = 0.95, auto.scale = TRUE, labels = TRUE, table = FALSE, ...)
hill(data, k, tail.index = TRUE)
hillPlot(data, option = c("alpha", "xi", "quantile"), start = 15,
         end = NA, reverse = FALSE, p = NA, ci = 0.95,
         auto.scale = TRUE, labels = TRUE, ...)
plotFittedGPDvsEmpiricalExcesses(data, threshold = NA, nextremes = NA)
RiskMeasures(out, p)

Arguments

alpha

numeric, probability level(s).

auto.scale

logical, whether plot should be automatically scaled.

ci

numeric, probability for asymptotic confidence bands.

data

numeric, data vector or timesSeries.

end

integer, maximum number of exceedances to be considered.

information

character, whether standard errors should be calculated with “observed” or “expected” information. This only applies to maximum likelihood type; for “pwm” type “expected” information is used if possible.

k

number (greater than or equal to 2) of order statistics used to compute the Hill plot.

labels

logical, whether axes shall be labelled.

legend.pos

if not NULL, position of legend().

pre.0.4.9

logical, whether behavior previous to version 0.4-9 applies (returning the risk measure estimate and confidence intervals instead of just invisible()).

like.num

integer, count of evaluations of profile likelihood.

main,xlab,ylab

title, x axis and y axis labels.

models

integer, count of consecutive gpd models to be fitted; i.e., the count of different thresholds at which to re-estimate \(\xi\); this many \(\xi\) estimates will be plotted.

ne

integer, count of excesses above the threshold.

nextremes

integer, count of upper extremes to be used.

object

list, returned value from fitting GPD

omit

integer, count of upper plotting points to be omitted.

optfunc

character, function used for ML-optimization.

verbose

logical indicating whether warnings are given; currently only applies in the case where type="pwm" and xi > 0.5.

option

logical, whether "alpha", "xi" (1 / alpha) or "quantile" (a quantile estimate) should be plotted.

out

list, returned value from fitting GPD.

p

vector, probability levels for risk measures.

ppoints.gpd

points in (0,1) for evaluating the GPD tail estimate.

reverse

logical, plot ordered by increasing threshold or number of extremes.

RM

character, risk measure, either "VaR" or "ES"

start

integer, lowest number of exceedances to be considered.

table

logical, printing of a result table.

tail.index

logical indicating whether the Hill estimator of \(alpha\) (the default) or \(1/alpha\) is computed.

threshold

numeric, threshold value.

type

character, estimation by either ML- or PWM type.

...

ellpsis, arguments are passed down to either plot() or optim() or nlminb().

Details

hillplot(): This plot is usually calculated from the alpha perspective. For a generalized Pareto analysis of heavy-tailed data using the gpd function, it helps to plot the Hill estimates for xi. See pages 286--289 in QRM. Especially note that Example 7.28 suggests the best estimates occur when the threshold is very small, perhaps 0.1 of the sample size (10--50 order statistics in a sample of size 1000). Hence one should NOT be using a 95 percent threshold for Hill estimates. MEplot(): An upward trend in plot shows heavy-tailed behaviour. In particular, a straight line with positive gradient above some threshold is a sign of Pareto behaviour in tail. A downward trend shows thin-tailed behaviour whereas a line with zero gradient shows an exponential tail. Because upper plotting points are the average of a handful of extreme excesses, these may be omitted for a prettier plot. plotFittedGPDvsEmpiricalExcesses(): Build a graph which plots the GPD fit of excesses over a threshold u and the corresponding empirical distribution function for observed excesses. RiskMeasures(): Calculates risk measures (VaR or ES) based on a generalized Pareto model fitted to losses over a high threshold. xiplot(): Creates a plot showing how the estimate of shape varies with threshold or number of extremes.

See Also

GEV

Examples

Run this code
# NOT RUN {
data(danish)
plot(danish)

MEplot(danish)

xiplot(danish)

hillPlot(danish, option = "alpha", start = 5, end = 250, p = 0.99)
hillPlot(danish, option = "alpha", start = 5, end = 60, p = 0.99)

plotFittedGPDvsEmpiricalExcesses(danish, nextremes = 109)
u <- quantile(danish, probs=0.9, names=FALSE)
plotFittedGPDvsEmpiricalExcesses(danish, threshold = u)

findthreshold(danish, 50)
mod1 <- fit.GPD(danish, threshold = u)

RiskMeasures(mod1, c(0.95, 0.99))
plotTail(mod1)

showRM(mod1, alpha = 0.99, RM = "VaR", method = "BFGS")
showRM(mod1, alpha = 0.99, RM = "ES", method = "BFGS")

mod2 <- fit.GPD(danish, threshold = u, type = "pwm")
mod3 <- fit.GPD(danish, threshold = u, optfunc = "nlminb")

## Hill plot manually constructed based on hill()

## generate data
set.seed(1)
n <- 1000 # sample size
U <- runif(n)
X1 <- 1/(1-U) # ~ F_1(x) = 1-x^{-1}, x >= 1 => Par(1)
F2 <- function(x) 1-(x*log(x))^(-1) # Par(1) with distorted SV function
X2 <- vapply(U, function(u) uniroot(function(x) 1-(x*log(x))^(-1)-u,
                                    lower=1.75, upper=1e10)$root, NA_real_)

## compute Hill estimators for various k
k <- 10:800
y1 <- hill(X1, k=k)
y2 <- hill(X2, k=k)

## Hill plot
plot(k, y1, type="l", ylim=range(y1, y2, 1),
     xlab=expression("Number"~~italic(k)~~"of upper order statistics"),
     ylab=expression("Hill estimator for"~~alpha),
     main="Hill plot") # Hill plot, good natured case (based on X1)
lines(k, y2, col="firebrick") # Hill "horror" plot (based on X2)
lines(x=c(10, 800), y=c(1, 1), col="royalblue3") # correct value alpha=1
legend("topleft", inset=0.01, lty=c(1, 1, 1), bty="n",
       col=c("black", "firebrick", "royalblue3"),
       legend=as.expression(c("Hill estimator based on"~~
                               italic(F)(x)==1-1/x,
                              "Hill estimator based on"~~
                               italic(F)(x)==1-1/(x~log~x),
                              "Correct value"~~alpha==1)))

## via hillPlot()
hillPlot(X1, option="alpha", start=10, end=800)
hillPlot(X2, option="alpha", start=10, end=800)
# }

Run the code above in your browser using DataLab