Residuals in Generalized Pareto Distribution

Question

I'm learning generalized Pareto distribution for fitting extreme value data. I came across an R package evir that is able to plot residuals. Residuals from a GPD would also follow an exponential distribution. GPD pdf for a random variable $y$ is given as.

$y = f(y|u,\xi,\beta) = \frac{1}{\beta}(1+\xi\frac{y-u}{\beta})^{{}-1-\frac{1}{\xi}}$

where $u$ is the threshold, $\xi$ is the shape parameter and $\beta$ is scale parameter, and $\xi \ne 0$ and $\beta >0$.

I'm not able to follow how the residuals are calculated for GPD.

My question is how to calculate residuals from a fitted GPD?

Here is a reproducible example:

## load the evir package
library(evir)

## Data danish for illustration
data(danish)


## Estimate parameters xi and beta, u = 10
fitmodel <- gpd(data = danish,threshold = 10)


## Paremters xi and beta
##$par.ests
##xi      beta 
##0.4968062 6.9745523 


## Plot residuals and QQ plot for residuals
plot(fitmodel)


##Output
#Make a plot selection (or 0 to exit): 
  
#1: plot: Excess Distribution
#2: plot: Tail of Underlying Distribution
#3: plot: Scatterplot of Residuals
#4: plot: QQplot of Residuals

#Selection: 3

#Selection: 4

Here are the plots of residuals (selection:3) and also QQ plot (Selection:4).

Answer 1

As suggested by the comment of @whuber, since the package documentation does not tell precisely what is plotted, the answer is to be found in the code of the plot method for the S3 class "gpd" which is the class of fitmodel object in the code chunk. The answer is on the line #270 of the R code. The shown points are $[i, \, e_i]$ with $$ e_i := \widehat{\xi}^{-1} \log\{ 1 + \widehat{\xi} z_i / \widehat{\beta} \} = -\log S_{\text{GPD}}(z_i;\, \widehat{\beta},\, \widehat{\xi}) $$ where $z_i:= y_i -u$ is the $i$-th the excess over the threshold and $S_{\text{GPD}}(z;\,\beta,\,\xi)$ denotes the survival function of the GPD with scale $\beta$ and shape $\xi$. Of course, the GPD location is $\mu = 0$ since the concern is the distribution of the excesses.

The distribution of the $z_i$ should be nearly a standard exponential because for a r.v. $Z$ the transformed r.v. $S_Z(Z)$ is standard uniform hence $-\log S_Z(Z)$ is standard exponential. The $z_i$ can be called "generalised residuals". Note that the transformation used is increasing hence a large value of the observation (excess) corresponds to a large value of the generalised residual. Inasmuch the standard exponential is well-known, we have an idea of what an "abnormally large" residual is. For example $\Pr\{Z > 3\} \approx 0.05$. The choice of the exponential distribution here is quite arbitrary: a uniform distribution could have been used as well. But the standard exponential is a good choice when the data can be regarded as coming from a GPD with a small shape $\xi \approx 0$, which is quite common in some applications: we then use a "mild" transformation of the excesses. Another choice would be the standard normal, but the strong deformation of the small excesses $\approx 0$ would be unnatural.

Generalised residuals are useful for "non-stationary models" as described in Chap 6 of Stuart Coles' book, with parameters depending on covariates, typically depending a time variable or functions of time. Then a pattern in the generalised residuals plotted against the time (or another covariate) would indicate a misspecification. I think that this motivated adding a lowess smooth on the plot produced by evir, because the index $i$ of the excess $z_i$ above is often the rank of the exceedance in time order.