# Peaks-over-Threshold Method Functions for fitting, analysing and risk measures according to POT/GPD ```r 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 ```r 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) ```

POT function