Can we fit extreme value distribution by build-in package?

Question

I try to find a package in R to fit Gumbel distribution by Block Maxima Approach using maximal likelihood function (see here) $$ G(x; \mu , \sigma)=\exp[-e^{-\frac{x-\mu}{\sigma}}]. $$

The block maxima approach is that we have $n$ iid random samples $X_1, X_2, \dots, X_n$ from $exp(1)$ (the limiting distribution of extreme value is Gumbel), and choose $m$ blocks with same size for $X_i$ and choose the maximal value of each block denoted by $Z_1, Z_2,\dots, Z_m$. Since $Z_i$ is the maximal value in $i-th$ block, we can approximate its distribution by Gumbel distribution.

In this website https://www.stat.purdue.edu/~huang251/1018.html, he used the package fevd to fit the generalized extreme value distribution with three parameters $(\mu, \sigma, \xi)$. Can we just fit two parameters? Is there any such package to fit Gumbel distribution?

Answer 1

Your simulation appears to be working as expected.

Here's a verification using optim() to directly maximise the log-likelihood function you've written yourself in the question (noting that the number of blocks, $m$, should actually be 1000 and not 100). The results are equivalent to what fevd() is estimating (within a margin of error).

# Data generating process
n <- 100   # observations within each block
m <- 1000  # number of blocks
x <- rexp(n*m, 1)
x_agg <- aggregate(x, list(gl(m, n)), max)

# Optional: Rescale to force mu=0 and sigma=1
# x_agg$x <- x_agg$x - log(n)

# Log-likelihood function: Eq. 3.9 Coles (2001) pg. 55
ll_fn <- function(par, z, m){
  mu <- par[1]
  sigma <- par[2]
  -m * log(sigma) - sum((z - mu)/sigma) - sum(exp(-((z - mu)/sigma)))
}

# Fit using optim ----
# Note fnscale set to maximise ll_fn
# hessian = TRUE for Fisher Info if needed
fit1 <- optim(par = c(1, 1), 
              fn = ll_fn, z = x_agg$x, m = m,
              hessian = TRUE, control = list(fnscale = -1))
# $par
# [1] 4.661471 1.006802
#
# $value
# [1] -1588.832

# Fit using extRemes ----
fit2 <- fevd(x_agg$x, type = "Gumbel")
# Negative Log-Likelihood Value:  1588.832 
#
# Estimated parameters:
# location    scale 
# 4.661528 1.006741 

In fact, without fitting a model at all you can verify that if the block maxima in this simulation are Gumbel distributed, this distribution can't be centred on $\mu = 0$.

Here's what the raw x_agg$x looks like just plotting the data

# Simple histogram and density overlay
hist(x_agg$x, probability = TRUE, breaks = 38, xlab = "Block Maxima - rexp(1)")
lines(density(x_agg$x), lty = 2)

The proof you've reproduced above leading you to expect the parameter estimates to approach the standard Gumbel ($\mu = 0$, $\sigma = 1$) is actually based on rescaling the block maxima by the log number of observations within each block, i.e., $\text{max}\{X_1, \dots, X_m\} - \log(n)$. This rescale does indeed result in the standard Gumbel, which you can verify using the code above, but note that when working with the raw block maxima, the parameter estimates will differ.

A more comprehensive discussion on this is provided in Chapter 3 of Coles, S. (2001). An introduction to statistical modeling of extreme values. London: Springer.