2.5 Extreme values

The largest value in a year, a century or a sample follows its own laws, just as averages follow the normal curve. Extreme value distributions describe those maxima, and let us estimate how high the 100-year flood is from far fewer than 100 years of data, with honest, and large, uncertainty.

import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm
from matplotlib.axes import Axes
from scipy import stats

from statly import style

style.use()
rng = np.random.default_rng(1)

elnino = (
    sm.datasets.elnino.load_pandas().data
)  # monthly sea surface temperature, Pacific, 1950 to 2010
annual_max = elnino.iloc[:, 1:].max(axis=1).to_numpy()  # the warmest month of each year
years = elnino.YEAR.to_numpy().astype(int)

2.5.1 The largest value has its own distribution

Most of statistics is about typical values. Engineers designing a dam, insurers pricing storm damage and hospitals planning for the worst week of a flu season care about the largest value instead. The maximum of n values behaves very differently from a single value: it moves up as n grows, and its spread changes.

fig, ax = plt.subplots(figsize=(8, 3.2))
grid = np.linspace(-4, 5.5, 400)
ax.plot(grid, stats.norm.pdf(grid), color=style.GREY, lw=2, label="one value")
for n, color in ((10, style.BLUE), (100, style.GOLD), (1000, style.CORAL)):
    maxima = rng.normal(size=(20_000, n)).max(axis=1)
    ax.hist(maxima, bins=80, density=True, color=color, alpha=0.5, label=f"largest of {n:,}")
ax.legend(fontsize=8)
ax.set(xlabel="value", ylabel="how common", yticks=[])
plt.show()

Figure 2.5.1: The distribution of the largest of n values drawn from a standard normal distribution (grey: a single value), for n = 10, 100 and 1,000, each from 20,000 simulations. The maximum moves right, slowly (the largest of 1,000 normal values is typically around 3.2), and its distribution is skewed to the right, unlike the normal curve it came from.

2.5.2 The generalised extreme value distribution

The central limit theorem says that sums of many values tend to a normal distribution, whatever the values came from. The extremal types theorem says something similar for maxima: suitably shifted and scaled, the maximum of many values tends to one family, the generalised extreme value (GEV) distribution, with a location μ, a scale σ and a shape ξ:

$$P(\text{max} \le x) = \exp\Big\{-\Big[1 + \xi\,\frac{x - \mu}{\sigma}\Big]^{-1/\xi}\Big\}$$

The shape decides the tail. ξ = 0 (the Gumbel case, the limit as ξ → 0) arises from light-tailed values such as normal ones; ξ > 0 (Fréchet) from heavy tails such as insurance losses, with no upper limit and very large extremes; ξ < 0 (reversed Weibull) when there is a hard upper limit.

grid = np.linspace(-3, 10, 600)
SHAPES = [
    (-0.3, style.BLUE, "bounded, ξ = −0.3"),
    (0.0, style.GREY, "Gumbel, ξ = 0"),
    (0.3, style.CORAL, "heavy tail, ξ = 0.3"),
]
fig, axes = style.subplots(2, height=3.1, width=10)
for xi, color, label in SHAPES:
    dist = stats.genextreme(c=-xi)  # scipy's shape parameter has the opposite sign
    axes[0].plot(grid, dist.pdf(grid), color=color, lw=2, label=label)
    axes[1].plot(grid, dist.sf(grid), color=color, lw=2)
axes[0].legend(fontsize=8)
axes[0].set(xlabel="maximum", ylabel="how common", yticks=[])
axes[1].set(
    yscale="log", ylim=(1e-5, 1.2), xlabel="maximum", ylabel="chance of exceeding (log)"
)
fig.tight_layout()
plt.show()

Figure 2.5.2: The three shapes of the generalised extreme value distribution, all with location 0 and scale 1. Left: densities. Right: the probability of exceeding each value, on a log scale. The heavy-tailed Fréchet case (ξ = 0.3) keeps a real chance of very large values long after the Gumbel case has none to speak of; the bounded case (ξ = −0.3) stops at an upper limit.

The usual way to use it is block maxima: take the largest value in each year (or month), and fit a GEV to those maxima by maximum likelihood (Chapter 2.3).

shape_c, loc, scale = stats.genextreme.fit(annual_max)
xi = -shape_c
gev = stats.genextreme(shape_c, loc, scale)

fig, axes = style.subplots(2, height=3.1, width=10)
axes[0].plot(years, annual_max, "o-", color=style.BLUE, ms=3, lw=1)
axes[0].set(xlabel="year", ylabel="warmest month (°C)")
width = style.hist(axes[1], annual_max, bins=15)
grid = np.linspace(annual_max.min() - 0.5, annual_max.max() + 1, 300)
axes[1].plot(grid, gev.pdf(grid) * len(annual_max) * width, color=style.CORAL, lw=2)
style.note(
    axes[1], f"ξ = {xi:.2f}, μ = {loc:.2f}, σ = {scale:.2f}", x=0.4, y=0.95, color=style.INK
)
axes[1].set(xlabel="warmest month (°C)", ylabel="years")
fig.tight_layout()
plt.show()

Figure 2.5.3: The warmest monthly sea surface temperature in the El Niño region of the Pacific in each year, 1950 to 2010. Left: the 61 annual maxima; the spikes are El Niño years, 1997 the largest. Right: their histogram with the fitted GEV density. The fitted shape ξ is close to 0 and slightly positive: close to the Gumbel case, with a tail a little heavier.

2.5.3 Return periods and return levels

The T-year return level is the value exceeded with probability 1/T in any one year: the “100-year flood” is the level with a 1% chance each year. It is not a level reached once every 100 years like clockwork. Over a period of n years, the chance of at least one exceedance is

$$P(\text{at least one in } n \text{ years}) = 1 - \Big(1 - \frac{1}{T}\Big)^n$$

so a house built to withstand the 100-year flood has a 26% chance of meeting a larger one during a 30-year mortgage.

periods = np.geomspace(1.1, 200, 80)
levels = gev.isf(1 / periods)


def bootstrap_levels(repeats: int = 1000) -> np.ndarray:
    """Return levels refitted to resampled years."""
    local = np.random.default_rng(0)
    out = []
    for _ in range(repeats):
        sample = local.choice(annual_max, len(annual_max))
        c, m, s = stats.genextreme.fit(sample)
        out.append(stats.genextreme(c, m, s).isf(1 / periods))
    return np.array(out)


band = np.percentile(bootstrap_levels(), [2.5, 97.5], axis=0)
ranks = np.arange(1, len(annual_max) + 1)
empirical_period = (len(annual_max) + 1) / ranks  # largest first

fig, axes = style.subplots(2, height=3.2, width=10.5)
axes[0].fill_between(periods, band[0], band[1], color=style.CORAL, alpha=0.2)
axes[0].plot(periods, levels, color=style.CORAL, lw=2)
axes[0].scatter(empirical_period, np.sort(annual_max)[::-1], s=12, color=style.BLUE, zorder=3)
axes[0].set(
    xscale="log", xlabel="return period T (years, log scale)", ylabel="return level (°C)"
)
spans = np.arange(1, 101)
for T, color in ((10, style.BLUE), (100, style.CORAL), (1000, style.GREY)):
    axes[1].plot(spans, 1 - (1 - 1 / T) ** spans, color=color, lw=2, label=f"{T}-year level")
axes[1].yaxis.set_major_formatter(style.PERCENT)
axes[1].legend(fontsize=8)
axes[1].set(xlabel="years of exposure", ylabel="chance of at least one exceedance")
fig.tight_layout()
plt.show()

Figure 2.5.4: Left: return levels of the annual maximum sea temperature from the fitted GEV (line), with a 95% bootstrap band (shaded, from refitting to 1,000 resamples of the 61 years) and the observed maxima placed at their empirical return periods (dots). Beyond the length of the record, the band widens quickly: the 100-year level is far less certain than the 10-year level. Right: the chance of at least one exceedance of the T-year level within a given number of years.

2.5.4 Peaks over a threshold

Using only one value per year throws away the other large values. The peaks-over-threshold approach keeps every value above a high threshold u, and models the excess over u. For a high enough threshold, excesses follow a generalised Pareto distribution with the same shape ξ as the GEV:

$$P(X - u > y \mid X > u) = \Big(1 + \xi\,\frac{y}{\sigma_u}\Big)^{-1/\xi}$$

The threshold is a trade-off: too low and the approximation fails; too high and too few values are left. A mean excess plot helps: above a suitable threshold, the average excess grows roughly in a straight line with the threshold.

daily = stats.genpareto(0.2, scale=4).rvs(size=365 * 30, random_state=rng) * (
    rng.uniform(size=365 * 30) < 0.4
)
thresholds = np.linspace(2, 45, 60)
mean_excess, half_width = [], []
for u in thresholds:
    over = daily[daily > u] - u
    mean_excess.append(over.mean())
    half_width.append(1.96 * over.std(ddof=1) / np.sqrt(len(over)))
mean_excess, half_width = np.array(mean_excess), np.array(half_width)
THRESHOLD = 20
excesses = daily[daily > THRESHOLD] - THRESHOLD
gp_shape, _, gp_scale = stats.genpareto.fit(excesses, floc=0)

fig, axes = style.subplots(2, height=3.1, width=10)
axes[0].fill_between(
    thresholds, mean_excess - half_width, mean_excess + half_width, color=style.BLUE, alpha=0.2
)
axes[0].plot(thresholds, mean_excess, color=style.BLUE, lw=2)
axes[0].axvline(THRESHOLD, color=style.CORAL, ls="--", lw=1)
axes[0].set(xlabel="threshold u (mm)", ylabel="mean excess over u (mm)")
width = style.hist(axes[1], excesses, bins=30)
grid = np.linspace(0, excesses.max(), 300)
axes[1].plot(
    grid,
    stats.genpareto(gp_shape, scale=gp_scale).pdf(grid) * len(excesses) * width,
    color=style.CORAL,
    lw=2,
)
style.note(
    axes[1],
    f"{len(excesses)} days above {THRESHOLD} mm\nξ = {gp_shape:.2f}",
    x=0.45,
    y=0.95,
    color=style.INK,
)
axes[1].set(xlabel=f"excess over {THRESHOLD} mm", ylabel="days")
fig.tight_layout()
plt.show()

Figure 2.5.5: Thirty years of simulated daily rainfall with a heavy tail. Left: the mean excess over each threshold (with 95% bands); it rises roughly in a straight line, as a generalised Pareto tail predicts, and gets noisy at high thresholds where few days remain. Right: the excesses over 20 mm (histogram) with the fitted generalised Pareto density; the shape estimate is close to the true tail used in the simulation, ξ = 0.2.

2.5.5 Cautions

• Extreme value theory assumes the maxima come from a stable process. A changing climate, a new dam or a new policy breaks that, and the model must include the trend.
• Extremes in consecutive days cluster (one storm, several wet days); thresholds should count clusters, not days.
• Extrapolating to return periods much longer than the record is possible only through the model, and the uncertainty grows fast; report the interval, not just the level.

2.5.6 In one picture

• Maxima follow their own distribution; it moves up and changes shape with n.
• Block maxima → GEV with location, scale and shape ξ (Gumbel ξ = 0, heavy tail ξ > 0, bounded ξ < 0).
• The T-year return level has chance 1/T each year; over n years, $1 - (1 - 1/T)^n$.
• Peaks over threshold → generalised Pareto; choose u with the mean excess plot.
• Long return periods are extrapolations: show the widening uncertainty.

Next: simulation, the tool that lets a computer answer probability questions by trying.