## Generate losses to work with
set.seed(271)
X <- rt(1000, df = 3.5) # in MDA(H_{1/df}); see MFE (2015, Section 16.1.1)
## (Sample) mean excess plot and threshold choice
mean_excess_plot(X[X > 0]) # we only use positive values here to see 'more'
## => Any value in [0.8, 2] seems reasonable as threshold at first sight
## but 0.8 to 1 turns out to be too small for the degrees of
## freedom implied by the GPD estimator to be close to the true value 3.5.
## => We go with threshold 1.5 here.
u <- 1.5 # thresholds
## An alternative way
ME <- mean_excess_np(X[X > 0])
plot(ME, xlab = "Threshold", ylab = "Mean excess over threshold")
## Mean excess plot with mean excess function of the fitted GPD
fit <- fit_GPD_MLE(X[X > u] - u)
q <- seq(u, ME[nrow(ME),"x"], length.out = 129)
MEF.GPD <- mean_excess_GPD(q-u, shape = fit$par[["shape"]], scale = fit$par[["scale"]])
mean_excess_plot(X[X > 0]) # mean excess plot for positive losses...
lines(q, MEF.GPD, col = "royalblue", lwd = 1.4) # ... with mean excess function of the fitted GPD
Run the code above in your browser using DataLab