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

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)
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).

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.

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

• How can you get the shape parameter? Jun 8, 2021 at 15:42
• logsigma is an estimate of the reciprocal of the shape parameter $\xi$ of a Fréchet distribution if I'm right. At any rate there's surely not a missing parameter estimate in the output. Jun 8, 2021 at 15:55
• The Intercept is the location; the logsigma is the scale; and $1/logsigma$ is the shape parameter? I thought there were three distinct parameters in the Frechet. I'm not following. Jun 8, 2021 at 16:05
• Is there an easy way to "confirm" that this is so by calculating the 51.4 value in the quote in the question? Jun 8, 2021 at 16:07
• (Intercept) is the maximum-likelihood estimate of the location parameter ($\sigma$ in the parameterization I've used) for the log-Fréchet distribution & therefore of the logarithm of the scale parameter for the Fréchet distribution. logsigma is the MLE of the scale parameter for the log-Fréchet distribution & therefore of the reciprocal of the shape parameter for the Fréchet distribution ($\xi$). There's no location parameter for this Fréchet distribution because the lower bound is fixed at nought (a natural enough lower bound for wind speeds, & for many other kinds of observations). Jun 8, 2021 at 18:58