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
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)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(QRM)
Loading required package: gsl
Loading required package: Matrix
Loading required package: mvtnorm
Loading required package: numDeriv
Loading required package: timeSeries
Loading required package: timeDate
Attaching package: 'QRM'
The following object is masked from 'package:base':
lbeta
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/QRM/POT.Rd_%03d_medium.png", width=480, height=480)
> ### Name: POT
> ### Title: Peaks-over-Threshold Method
> ### Aliases: POT fit.GPD findthreshold plotTail MEplot showRM xiplot
> ### RiskMeasures hill hillPlot plotFittedGPDvsEmpiricalExcesses
> ### Keywords: distributions
>
> ### ** Examples
>
> 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)
[1] 17.06847
> mod1 <- fit.GPD(danish, threshold = u)
>
> RiskMeasures(mod1, c(0.95, 0.99))
p quantile sfall
[1,] 0.95 9.401016 25.64404
[2,] 0.99 27.452048 69.01994
> 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")
Warning message:
In fit.GPD(danish, threshold = u, type = "pwm") :
Asymptotic standard errors not available for PWM Type when xi > 0.5.
> 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)
>
>
>
>
>
> dev.off()
null device
1
>