How to extract the shape parameter of a Fréchet fitted model using the R SPREDA package?

Question

I'm trying to follow this post, which fits a Frechet distribution to some wind measurements as follows:

require('SPREDA')
require(e1071)
require(extRemes)
require(survival)

#The R functions for performing the extreme value analysis (for maxima)
#can be downloaded from: 
source("http://www.datall-analyse.nl/R/eva_max.R")

#data from Castillo et al., Table 1.1, p. 9-10
#yearly maximum wind speed (in miles/hour)
evobs <- scan("http://www.datall-analyse.nl/blog_data/extremes_Table1-1.txt")

#explore data visually
options(repr.plot.width=7, repr.plot.height=7)
hist(evobs)

#Fréchet
frechetmod <- Lifedata.MLE(Surv(evobs) ~ 1, dist="frechet")
frechetmod

At some point, the author uses the output or the function Lifedata.MLE() to calculate the following return level:

#note: for the wind data, a return period of 20 means that once every #20 years the wind speed is (on average) expected to be larger than #muG+qlev(1-1/20)*sigmaG=49.4 mph (in case of a Gumbel distribution), #or exp(muF+qlev(1-1/20)*sigmaF)=51.4 mph (in case of a Fréchet distribution) #(this expected wind speed is also called the return level)

which I suppose requires the xi shape parameter to calculate the 51.4 mph, using the quantile formula:

$$X_T=\mu +\frac{\sigma\left( 1 - \left( -\log\left(1-1/T \right)\right)^\xi\right)}{\xi}$$

where $T$ is the return period ($1/(1-F)$, where $F$ is the distribution function).

REFERENCE: Overeem et al. (2010), Water Resources Research, 46, "Extreme value modeling of areal rainfall from weather radar"

However, this parameter is not included in the output of frechetmod:

Call:
model.frame(formula = Surv(evobs) ~ 1)

Coefficients:
(Intercept)    logsigma 
  3.3579367   0.1957559 

Loglikelihod:
-173.4035 (df=2)

How can I get this shape parameter estimate?

---

The $51.4$ value in the quoted paragraph is obtained in the post as:

frechetmod <- Lifedata.MLE(Surv(evobs) ~ 1, dist="frechet")
frechetmod

muF <- coef(frechetmod)[1]
sigmaF <- coef(frechetmod)[2]

exp(muF+qlev(1-1/20)*sigmaF)

using the Package ‘SPREDA’ and the qlev() call, which gives the quantile of the Standard Largest Extreme Value Distribution, but I don't know how to reconcile with the closed equation above.

Answer 1

I believe you're estimating the location ((Intercept)) & scale (logsigma) parameters of a log-Fréchet distribution (corresponding to a Fréchet distribution with a fixed lower bound of nought–a natural enough lower bound for wind speeds, & for many other kinds of observations). So if the density function for a Fréchet random variable $Y$, parametrized with scale $\beta$ & shape $\alpha$, is

$$f_Y(y) = \alpha\cdot\left(\frac{y}{\beta}\right)^{-(\alpha+1)} \cdot \exp\left[-\left(\frac{y}{\beta}\right)^{-\alpha}\right]$$

& given $Z = \log(Y)$

$$\begin{align} f_Z(z) &= \alpha\cdot\left(\frac{\exp z}{\beta}\right)^{-(\alpha+1)} \cdot \exp\left[-\left(\frac{\exp z}{\beta}\right)^{-\alpha}\right] \cdot \left| \frac{\operatorname{d}\exp z}{\operatorname{d}z} \right| \\ &= \alpha\cdot \exp \{-\alpha (z-\log\beta) - \exp[-\alpha (z-\log\beta)]\} \end{align} $$

the location parameter of the log-Fréchet distribution will therefore be $\log\beta$ & the scale parameter $\frac{1}{\alpha}$. (In a regression model the former would be a function, perhaps a linear function, of the predictors—here of course you've fitted an intercept-only model.) And of course maximum-likelihood estimates are equivariant to parameter transformations.

Check the help pages and/or source code for the Lifedata.MLE function to confirm. (Note that plugging $y=51.4$, $\beta=\exp(3.3579)$ & $\alpha=\frac{1}{0.19576}$ into the Fréchet distribution function

$$ F_Y(y) = \exp\left[-\left(\frac{y}{\beta}\right)^{-\alpha}\right] $$

gives $0.95$, which is reassuring.)

---

The above parametrization of the Fréchet can be extended to include a location parameter $\tau$, with the distribution function

$$ F_Y(y) = \exp\left[-\left(\frac{y - \tau}{\beta}\right)^{-\alpha}\right] $$ It's a common one, with $\tau$ conveniently interpretable as the lower bound on $Y$. The parametrization you're using corresponds to the distribution function $$ F_Y(y) = \exp\left[-\left(\frac{y - (\mu + \sigma/\xi)}{-\sigma/\xi}\right)^{1/\xi}\right] $$ & I suppose is more suited to describing Fréchet distributions within the broader Fisher–Tippett family.

From these it's easy to see the relations:

$$ \begin{align} \alpha &= \frac{-1}{\xi}\\ \beta &= \frac{-\sigma}{\xi}\\ \tau &= \mu + \frac{\sigma}{\xi} \end{align} $$

A numerical check with R confirms 51.4 mph as the return level:

alpha <- 1/0.1957559
beta <- exp(3.3579367)
tau <- 0

zeta <- -1/alpha
sigma <- -beta*zeta
mu <- tau - sigma/zeta

T <- 20
F <- (T-1)/T

y <- mu+sigma*(1-(-log(F))^zeta)/zeta