Return level plots for GEV-distribution

Question

I was reading An Introduction to Statistical Modeling of Extreme Values by Stuart Coles, and I ran into a problem whilst trying to replicate a basic return level graph in R. For context, I first present the definition of a return level plot (p.49 of the book):

Estimates of extreme quantiles of the annual maximum distribution are then obtained by inverting Eq. (3.2): $$z_p = \begin{cases} \mu - \frac{\sigma}{\xi}\left[1-\{-\log(1-p)\}^{-\xi}\right], & \text{for } \ \xi \neq 0, \\ \mu-\sigma \log\{ -\log(1-p)\}, & \text{for } \ \xi =0,\end{cases} \tag{3.4}$$ where $G(z_p) = 1-p$. In common terminology, $z_p$ is the return level associated with the return period $1/p$, since to a reasonable degree of accuracy, the level $z_p$ is expected to be exceeded by the annual maximum in any particular year with probability $p$.
Since quantiles enable probability models to be expressed on the scale of data, the relationship of the GEV model to its parameters is most easily interpreted in terms of the quantile expressions (3.4). In particular, defining $y_p = -\log(1-p),$ so that $$z_p = \begin{cases} \mu - \frac{\sigma}{\xi}\left[1-y_p^{-\xi}\right], & \text{for } \ \xi \neq 0, \\ \mu-\sigma \log (y_p), & \text{for } \ \xi =0,\end{cases} $$ it follows that, if $z_p$ is plotted against $y_p$ on a logarithmic scale - or equivalently, if $z_p$ is plotted against $\log(y_p)$ - the plot is linear in the case $\xi = 0$.If $\xi < 0$ the plot is convex with asymptotic limit as $p \to 0$ at $\mu - \sigma/\xi$; if $\xi> 0$ the plot is concave and has no finite bound. This graph is a return level plot. Because of the simplicity of interpretation, and because the choice of scale compresses the tail of the distribution so that the effect of extrapolation is highlighted, return level plots are particularly convenient for both model presentation and validation. Fig 3.1 shows return level plots for a range of shape parameters.

As a reminder, $\mu,\xi \in \mathbb{R}$ and $\sigma > 0$ are the parameters of the GEV distribution. Now, I wanted to replicate Figure 3.1. To do this, I did the following:

mu   <- 0; sigma <- 1
p    <- ppoints(1000)
yp   <- -log(1-p)
zp_1 <- mu - (sigma/0.2)*(1-yp^(-0.2))   #xi =  0.2
zp_2 <- mu - sigma*log(yp)               #xi =   0
zp_3 <- mu - (sigma/-0.2)*(1-yp^(0.2))   #xi = -0.2

plot(log(yp), zp_1, type = "l")
lines(log(yp), zp_2)
lines(log(yp), zp_3)

which produced the following plot:

Evidently, this plot is not the same as Figure 3.1. However, if instead of graphing $\{(\log(y_p), z_p), 0 <p<1\}$ we graph $\{(-\log(y_p), z_p), 0 <p<1\}$, we obtain the following:

plot(-log(yp), zp_1, type = "l")
lines(-log(yp), zp_2)
lines(-log(yp), zp_3)

which seems to be exactly Figure 3.1. I know the difference is just a sign, but is there something I am missing, theoretically? Did the author make a mistake? Any help is appreciated.

Answer 1

In the rest of the book, several return level plots are made by plotting the time period $T=1/p$ on the x-axis with a logarithmic scale. (like the graph on the bottom of this post)

But for this one graph, instead of plotting $\log(T) = \log(1/p) = - \log(p)$ the author chose to plot with $y_p = -\log(1-p)$.

It is not so clear why it was done because that type of return level plot does not return anywhere else in the book. Possibly it was done to show how to linearize the case of $\xi = 0$. This $y_p$ does occur a few times in the book and is explained as making it easier to evaluate the different models, because of the differences in linearity, concavity and convexity for respectively $\xi = 0$, $\xi > 0$ and $\xi < 0$.

Because the plot with $\log(y_p)$ is less typical, and because $y_p$ has no clear intuitive meaning (it is just a technical aid to linearized the graph), a little mistake with a minus sign might have slipped into to text.

The equation in the text, $z_p = \mu - \sigma \log y_p$, shows clearly that you get a line with negative slope. It may likely have been inverted to make the plot appear more like the other plots with $\log(T)$ on the x-axis.

mu   <- 0; sigma <- 1
p    <- ppoints(1000)
yp   <- -log(1-p)
zp_1 <- mu - (sigma/0.2)*(1-yp^(-0.2))   #xi =  0.2
zp_2 <- mu - sigma*log(yp)               #xi =   0
zp_3 <- mu - (sigma/-0.2)*(1-yp^(0.2))   #xi = -0.2

plot((1/p), zp_1, type = "l", log = "x",
     xlab = "time period", ylab = "return level", xaxt = "n")
axis(1,at = 10^c(0:4))
lines((1/p), zp_2)
lines((1/p), zp_3)
title("return level versus time period (1/p)")

Edit:

As mentioned in the discussion in the comments. Several plots in the book are actually also using the $1/y_p$ on the x-axis, even though the axis reads 'return period' (in the limit you do get that $T \sim 1/y_p$).

Still, some of the graphs do use the return period $T$. The difference can be seen based on the labels. When the label reads 'return period (years)', with the additional specification of units in years, then the graph relates to $T$ instead of $1/y_p$. See also the original S-PLUS files, which contains different functions for plotting (also seen in the ismev port to R here https://github.com/cran/ismev/blob/master/R/gpd.R).

See for instance Figure 9.3 which has an inflection point.