Skip to contents

Theory

The generalized extreme value (GEV) distribution arises as the non-degenerate limiting distribution of maxima. A key property of the generalized extreme value distribution is max-stability, meaning that the distribution of maxima of TT \in \mathbb{N} is the same up to a location-scale transform: if F𝖦𝖤𝖵F \sim \mathsf{GEV}, then there exist location and scale constants bb \in \mathbb{R} and a>0a>0 such that FT(ax+b)=F(x)F^T(ax+b) = F(x). The parameters of the new distribution are easily derived: if Xi𝖦𝖤𝖵(μ,σ,ξ)X_i \stackrel{\cdot}{\sim} \mathsf{GEV}(\mu,\sigma, \xi), then max{X1,,XT}𝖦𝖤𝖵(μT,σT,ξ)\max\{X_1, \ldots, X_T\} \sim \mathsf{GEV}(\mu_T, \sigma_T, \xi) with μT=μσ(1Tξ)/ξ\mu_T = \mu - \sigma(1-T^\xi)/\xi and σT=σTξ\sigma_T = \sigma T^\xi for ξeq0\xi \mathbb{N}eq 0, or μT=μ+σlog(T)\mu_T = \mu +\sigma \log(T) and σT=σ\sigma_T = \sigma if ξ=0\xi=0.

In practice, the GEV is fitted to data that are partitioned into kk blocks of sizes mm, assumed eequal for simplicity, so that n=kmn = km. The quality of the generalized extreme value approximation increases with the block size mm. For fixed block sizes, the estimated shape parameter generally differs from its asymptotic counterpart and this discrepancy usually introduces bias if one extrapolates far beyond the range of the observed data.

Let F(x)F(x) denote a thrice-differentiable distribution function with endpoint x*x^* and density f(x)f(x) and define s(x)=F(x)log{F(x)}/f(x)s(x)=-F(x)\log\{F(x)\}/f(x). The existence of the limit ξ=limns{bn}\xi_{\infty} = \lim_{n \to \infty} s'\{b_n\} is necessary and sufficient for convergence to an extreme value distribution GG, limnFn(anx+bn)=exp{(1+ξx)1/ξ}=G(x).\begin{align*} \lim_{n \to \infty} F^n(a_nx+b_n) = \exp\left\{-(1+\xi_\infty x)^{-1/\xi_\infty}\right\}=G(x). \end{align*}

Smith (1987) shows that, for any x{y:1+ξy>0}x \in \{y:1+\xi_\infty y >0\}, there exists zz such that log[F{v+xs(v)}]log{F(v)}={1+xs(z)}1/s(z),v<z<v+xs(v).\begin{align*} \frac{-\log[F\{v+xs(v)\}]}{-\log\{F(v)\}} = \left\{1+xs'(z)\right\}^{-1/s'(z)}, \qquad v < z < v+xs(v). \end{align*} For each n1n \geq 1, setting v=bnv=b_n and an=s(bn)a_n=s(b_n) yields Fn(anx+bn)=exp[{1+s(z)x}1/s(z)]+O(n1)\begin{align*} F^n(a_nx+b_n)=\exp\left[-\left\{1+s'(z)x\right\}^{-1/s'(z)}\right] + \mathrm{O}(n^{-1}) \end{align*} for z[min(anx+bn,bn),max(anx+bn,bn)]z \in [\min(a_nx+b_n, b_n), \max(a_nx+b_n, b_n)], depending on the support of FF. The ultimate approximation replaces s(z)s'(z) by s(x*)=ξs'(x^*)=\xi_{\infty}, whereas Smith (1987) suggests instead suggests taking s(bn)s'(b_n), which is closer to s(anx+bn)s'(a_nx+b_n) than is s(x*)s'(x^*).

The maximum likelihood estimate of the shape should change as the threshold or the block size increases. Consider for simplicity yearly maxima arising from blocking m=nym=n_y observations per year. We are interested in the distribution of the maximum of NN years and thus the target has roughly shape ξNny\xi_{Nn_y}, but our estimates will instead give approximately ξny\xi_{n_y}; an extrapolation error arises from this mismatch between the shape values, which would be constant if the observations were truly max-stable. The curvature of ss' determines how stable the estimates of ξt\xi_t are when the extrapolation window increases.

We illustrate the previous discussion by looking at block maximum from a log-normal distribution.

library(mev)
set.seed(123)
x  <- seq(1, 30, length = 200)
fitted <- t(
  replicate(
    n = 200,
    fit.gev(
      xdat = apply(
        X = matrix(rlnorm(30 * 100), 
               ncol = 30), 
        MARGIN = 1, 
        FUN = max),
      method = "BFGS")$est))
penult.bm.lnorm <-
  mev::smith.penult(model = "bm",
                    m = (m <- 30),
                    family = "lnorm")

par(mfrow = c(1, 2), mar = c(5, 5, 1, 1))
hist(
  fitted[, 3],
  probability = TRUE,
  breaks = 20,
  xlab = "estimated shape",
  main = ""
)
segments(
  y0 = -0.5,
  y1 = 0,
  x0 = penult.bm.lnorm$shape,
  col = "red",
  lwd = 2
)
segments(
  y0 = -0.5,
  y1 = 0,
  x0 = 0,
  col = "black",
  lwd = 2
)

p30 <- penult.bm.lnorm
x = seq(0, 100, length = 400)
N <- 1000
N30 = N / 30
plot(
  x,
  N * exp((N - 1) * plnorm(x, log.p = TRUE)) *
    dlnorm(x),
  type = "l",
  bty = "l",
  ylim = c(0, 0.1),
  xlab = "x",
  ylab = "density"
)
# Get parameters of maximum of N GEV
maxstabp <- function(loc, scale, shape, N) {
  if (!isTRUE(all.equal(shape, 0, check.attributes = FALSE))) {
    mut <- loc - scale * (1 - N ^ shape) / shape
    sigmat = scale * N ^ shape
    return(c(
      loc = mut,
      scale = sigmat,
      shape = shape
    ))
  } else{
    mut <- loc + scale * log(N)
    return(c(
      loc = mut,
      scale = scale,
      shape = shape
    ))
  }
}
p30e <-
  maxstabp(
    loc = p30$loc,
    scale = p30$scale,
    shape = p30$shape,
    N = N30
  )
lines(
  x,
  evd::dgev(
    x,
    loc = p30e['loc'],
    scale = p30e['scale'],
    shape = p30e['shape']
  ),
  col = "red",
  lwd = 2,
  lty = 2
)
p30l <-
  maxstabp(
    loc = p30$loc,
    scale = p30$scale,
    shape = 0,
    N = N30
  )
lines(
  x,
  evd::dgev(
    x,
    loc = p30e['loc'],
    scale = p30l['scale'],
    shape = 0
  ),
  col = 1,
  lty = 3,
  lwd = 1
)
legend(
  x = "topright",
  legend = c("exact", "ultimate", "penultimate"),
  col = c(1, 1, 2),
  lwd = c(2, 2, 1),
  bty = "n",
  lty = c(1, 3, 2)
)

The left panel illustrates how maximum likelihood estimates of the parameters for repeated samples from a log-normal distribution are closer to the penultimate approximation than to the limiting tail index ξ=0\xi_{\infty}=0. Using max-stability, we could get an estimate of the distribution of the maximum of NnyNn_y observations. In this case, the penultimate shape for m=1000m=1000 is approximately ξ1000=0.215\xi_{1000}=0.215, compared to ξ300.284\xi_{30} \approx 0.284, which is far from the limiting value ξ=0\xi_{\infty}=0. The extrapolated density estimate is based on F11000/30F_1^{1000/30}, the generalized extreme value density associated to the distribution function of the penultimate approximation F11000/30F_1^{1000/30}, with F1𝖦𝖤𝖵(a30,b30,ξ30)F_1 \sim\mathsf{GEV}(a_{30}, b_{30}, \xi_{30}). The penultimate approximation F2F_2 is more accurate than the ultimate approximation 𝖦𝖤𝖵(a30,b30,ξ)\mathsf{GEV}(a_{30}, b_{30}, \xi_\infty), which is too short tailed.

The function smith.penult computes the penultimate approximation for parameter models in R for the distribution assuming a model family with arguments dfamily, qfamily and pfamily exist. It returns the penultimate parameters for a given block size or quantile in case of threshold models. These functions are particularly useful for simulation studies in which one wants to obtain an approximation for FNF^N for large NN, to get an approximation to the asymptotic distribution of a test statistic for finite NN; the 𝖦𝖤𝖵\mathsf{GEV} penultimate approximation is more numerically stable for certain models and its moments and cumulants are easily derived.

The following plot illustrates the penultimate shape parameter for the normal distribution based on threshold exceedances. In finite samples, the upper tail is bounded, but the quality of the approximation by the Gumbel distribution increases as uu \to \infty, albeit very slowly.

u <- seq(0.8, 0.9999, by = 0.0001)
penult <- smith.penult(family = "norm",
                       method = "pot",
                       u = qnorm(u))
plot(
  u,
  penult$shape,
  bty = "l",
  type = 'l',
  xlab = 'quantile',
  ylab = 'penultimate shape'
)

While the above discussion has been restricted to block maxima and the generalized extreme value distribution, similar results exists for peaks over threshold using a generalized Pareto distribution. For threshold exceedances, we start with the reciprocal hazard r(x)={1F(x)}/f(x)r(x) = \{1-F(x)\}/f(x) in place of s(x)s(x). The normalizing constants are (u,r(u))(u, r(u)) and the penultimate shape parameter for the generalized Pareto distribution is r(u)r'(u).

We can relate threshold exceedances to the distribution of maxima as follows: suppose a generalized Pareto distribution FF is fitted to threshold exceedances above a threshold uu, with parameters (ζu,σu,ξ)(\zeta_u, \sigma_u, \xi), where ζu\zeta_u denotes the unknown proportion of points above the threshold uu. If there are nyn_y observations per year on average, then the TT-year return level is z1/T=u+σu/ξ{(Tnyζu)ξ1}z_{1/T}=u+\sigma_u/\xi\{(Tn_y\zeta_u)^\xi-1\}. We take FζuTnyF^{\zeta_uTn_y} as the approximation to the TT-year maxima distribution, conditional on exceeding uu.

Exercice

  1. Sample observations from the standard Gaussian distribution 𝖭𝗈(0,1)\mathsf{No}(0,1) for m=25,50,100,1000m = 25, 50, 100, 1000 and k=1000k=1000 replications.
  2. Fit a generalized extreme value distribution to block maxima of size mm. Obtain an estimate of the distribution of m=1000000m = 1000000 observations from a standard Gaussian using max-stability (the functions fit.gev and gev.mle can be useful).
  3. Plot the density of the resulting approximations and compare them to the density of the true distribution Nϕ(x)ΦN1(x)N\phi(x)\Phi^{N-1}(x).
  4. Superpose the density of the penultimate approximation for m=1000000m = 1000000 obtained using the Smith penultimate approximation (smith.penult).

References

Smith, Richard L. 1987. “Approximations in Extreme Value Theory.” Edited by University of North Carolina Chapel Hill Center for Stochastic Processes.