3
$\begingroup$

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.

$\endgroup$
0

1 Answer 1

2
+100
$\begingroup$

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
$\endgroup$
18
  • $\begingroup$ How can you get the shape parameter? $\endgroup$ Jun 8, 2021 at 15:42
  • $\begingroup$ 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. $\endgroup$ Jun 8, 2021 at 15:55
  • $\begingroup$ 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. $\endgroup$ Jun 8, 2021 at 16:05
  • $\begingroup$ Is there an easy way to "confirm" that this is so by calculating the 51.4 value in the quote in the question? $\endgroup$ Jun 8, 2021 at 16:07
  • 1
    $\begingroup$ (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). $\endgroup$ Jun 8, 2021 at 18:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.