Extreme value theory

Extreme value theory studies limiting probability distributions of normalized sample maxima and related extremes. The standard limits for maxima include the Fréchet distribution, negative Weibull distribution, and Gumbel distribution.

The question is often practical: if a batch finishes only when its slowest task finishes, its completion time is a sample maximum. A typical individual delay or a sample mean does not answer that question. The central limit theorem describes suitably normalized sums and sample means under its assumptions; it does not give a limit law for sample maxima. Extreme value theory instead describes how the largest observation changes as the number of opportunities for a large observation grows.

For independent and identically distributed random variables $X_1,X_2,\dots$ with distribution function $F$, define the sample maximumThe event $M_n\le x$ means that every observation is at most $x$, so independence givesExtreme value theory asks whether there are constants $a_n>0$ and $c_n\in \mathbb{R}$ for which $(M_n-c_n)/a_n$ has convergence in distribution to a non-degenerate probability distribution. Equivalently, $F(c_n+a_{n}z{)}^n$ must tend to a limiting distribution function at its continuity points. The extremal types theorem says that any such limit belongs, up to a change of location and positive scale, to the Gumbel distribution, Fréchet distribution, or negative Weibull distribution. It does not guarantee that a limit exists for every $F$. ^[ref].

The maximum of independent exponential random variables gives a particularly transparent worked example: its center grows as $(\log n)/\lambda$, and its fluctuations have a Gumbel distribution limit on the scale $1/\lambda$. Here and below $\log$ denotes the natural logarithm. This normalization depends on the exponential distribution. Other survival functions can lead to the Fréchet distribution or negative Weibull distribution, with different constants $a_n,c_n$; the Gumbel distribution and the exponential example's normalization are not universal.

Independence is essential to the product calculation. If instead $X_1=X_2=\cdots =X$ for one rate-$\lambda$ exponential random variable, then $M_n=X$ for every $n$. Its distribution function stays $1-e^{-\lambda x}$ for $x\ge 0$ and zero below zero, and its median stays $(\log 2)/\lambda$. Perfect dependence creates no logarithmic growth: repeating the same random variable supplies no new chance of a larger value.

A historical starting point is R. A. Fisher and L. H. C. Tippett's 1928 paper on limiting laws of the largest or smallest sample member, published in volume 24, issue 2, pages 180–190 of the Proceedings of the Cambridge Philosophical Society. Its abstract studies the possible limiting forms and emphasizes that approach to a limit can be slow. ^[ref].

Let $n\ge 1$ be an integer and let $X_1,\dots ,X_n$ be independent and identically distributed random variables with an exponential distribution of rate $\lambda >0$. Their sample maximum $M_n$ is nonnegative. Substituting the exponential distribution's cumulative distribution function into the distribution functions of independent minima and maxima gives the exact lawFor example, at $x=(\log 2)/\lambda$ each observation is below the threshold with probability $1/2$, but all $n$ observations are below it with probability $2^{-n}$. Thus an individual observation's median becomes an increasingly unlikely upper bound for the whole batch. The general order statistic density in Siegrist's treatment, with $k=n$, also gives this law on integration. ^[ref].

Solving $(1-e^{-\lambda x}{)}^n=p$ gives the exact quantile function, for $0<p<1$,In a toy model with $\lambda =1{\mathrm{s}}^{-1}$, the median maximum is $0.693147$ seconds for $n=1$, $2.703555$ seconds for $n=10$, and $4.975147$ seconds for $n=100$. All three use the exact distribution function; the increase comes from taking a maximum over more observations. The $n=2$ case connects to the existing maximum of two exponentials as a scaled sum.

The Gumbel distribution explains the logarithmic shift visible here: after subtracting $(\log n)/\lambda$ and measuring deviations in units of $1/\lambda$, the distribution functions approach one fixed shape. The next derivation makes that convergence in distribution explicit.

For $1\le k<n$ and a positive threshold $t=X_{(n-k)}$, the Hill estimator averages the logarithmic ratios of the top $k$ order statistics to $t$. For a Fréchet distribution domain with shape $\alpha$, it estimates the reciprocal tail index $1/\alpha$. At an untied threshold it equals the empirical distribution function plug-in tail integral ${\int }_t^{\infty }(1-F_n(x))/(1-F_n(t))dx/x$. At a tied threshold the plug-in denominator counts strictly larger observations and the identity needs adjustment. Fixed $k$ does not by itself yield a consistent tail-index estimate.

Every non-degenerate limit in convergence in distribution of positively scaled and centered sample maxima is a max-stable distribution. Up to positive affine transformations its law is the Gumbel distribution, a Fréchet distribution with positive shape parameter, or a negative Weibull distribution with positive shape parameter. This classifies possible limits, without asserting that every distribution function admits a limiting normalization.

A non-degenerate distribution function is max-stable when every sample maximum of $n$ independent and identically distributed random variables can be positively scaled and centered to have the original law. Thus for every integer $n\ge 1$ there are $a_n>0,b_n\in \mathbb{R}$ with $F(a_{n}x+b_n{)}^n=F(x)$. The extremal types theorem classifies these laws, up to affine changes, as the Gumbel distribution, Fréchet distribution, and negative Weibull distribution.

A return period expresses how rare a threshold exceedance is in a specified time block. Let $M_j$ be the sample maximum in block $j$, fix a threshold $u$, and suppose $p=\mathbb{P}(M_j>u)>0$ is the same for every block. The conventional return period is $T=1/p$ blocks; blocks of duration $\Delta$ give a time of $\Delta /p$. Thus an annual exceedance probability $p=0.01$ means a 100-year return period. This names a frequency scale, not a schedule or a guarantee of one exceedance in each century. The annual convention and the misleading “once every 100 years” interpretation are discussed in USGS: the 100-year flood.

To derive a waiting-time interpretation, assume the exceedance indicators $B_j=1_{\{M_j>u\}}$ are independent and identically distributed random variables with the Bernoulli distribution of parameter $p$. Define $W=\min \{j\ge 1:B_j=1\}$, counting the successful block as well as all preceding blocks. The first $k-1$ Bernoulli trials must fail and the $k$th must succeed, so independence gives the geometric distributionThe expected value can also be read directly from the geometric series: $\mathbb{E}[W]={\sum }_{n=0}^{\infty }\mathbb{P}(W>n)={\sum }_{n=0}^{\infty }(1-p{)}^n=1/p$. This counts blocks from the starting boundary, not the exact time of an exceedance within its block. Counting only failed blocks instead would give expected value $T-1$. The waiting law and its memorylessness of the geometric distribution are developed in Kyle Siegrist, Random Services: the geometric distribution.

For a horizon of $n$ complete years, at least one exceedance is the complement of no exceedances. Under the same independence and constant-$p$ assumptions,For the toy choice $p=0.01$, this gives $9.6\%$ in 10 years, $26.0\%$ in 30 years, $39.5\%$ in 50 years, and $63.4\%$ in 100 years. Even after 100 years there is a $36.6\%$ chance of no exceedance. These numbers are exact model calculations rounded to one decimal place, not measured frequencies. The median waiting count is 69 years, since $1-0.99^{68}<1/2\le 1-0.99^{69}$; the expected value is 100 years. A mean and a median answer different questions.

The exceedance count $N_n={\sum }_{j=1}^n{B}_j$ has the binomial distribution under independence, with expected value $np$. Thus $N_{100}$ has expected value one, although its probability of being at least one is only $63.4\%$: some realizations contribute several exceedances. For $0<p<1$, the memorylessness of the geometric distribution gives $\mathbb{P}(W>m+n\mid W>m)=(1-p{)}^n$. A long quiet spell does not make the next independent year's exceedance overdue, and an exceedance this year does not force a quiet interval. Two specified consecutive years both exceed with probability $p^2=0.0001$ in this toy model.

A reproducible tabletop experiment makes the assumptions tangible. Put 100 identical slips numbered 1 to 100 in a bag; after thoroughly mixing, draw one and record whether it is 100. Replace the slip and mix again before every draw. Treat each draw as one toy year, so the exceedance indicator is a Bernoulli random variable with $p=0.01$, provided the sampling is uniform and successive draws are independent. Repeat complete 30-draw runs and record the fraction containing at least one 100; separately repeat runs stopped at the first 100 and record their waiting counts. The first fraction targets $1-0.99^{30}$; the waiting-count histogram targets $0.99^{k-1}(0.01)$. Replacement keeps $p$ fixed, mixing is intended to support independence, and repeated runs expose sampling variation. This is a proposed demonstration of the toy model; no experimental results are asserted here.

Stationarity supports a time-invariant annual probability distribution, but does not imply independence. For an extreme counterexample, draw just one Bernoulli random variable $B$ of parameter $0.01$ and set $B_j=B$ for every year. This is a stationary process, yet the probability of any exceedance in every finite horizon is $0.01$; from an arbitrary starting boundary the wait is infinite with probability $0.99$. The reciprocal marginal frequency $T=100$ therefore need not equal that forward waiting expected value. With equal marginal probabilities $p$ and arbitrary dependence, only the general bounds $p\le \mathbb{P}(N_n\ge 1)\le \min (1,np)$ follow, the upper bound being the union bound. Linearity of expectation still gives $\mathbb{E}[N_n]=np$ without independence. Mere zero correlation does not restore the product formula; even pairwise independence does not imply joint independence. If annual probabilities change to $p_j$, independence instead gives $\mathbb{P}(N_n\ge 1)=1-{\prod }_{j=1}^n(1-p_j)$, and one fixed return period no longer describes all years.

A return level answers the inverse question to a return period: which threshold has the specified rarity? Let $F(z)=F_{\mathrm{annual}}(z)=\mathbb{P}(M_j\le z)$ be the common cumulative distribution function of annual sample maxima in a model with stationarity. For a continuous, strictly increasing $F$ and $T>1$, the $T$-year return level isIt has the same units as $M_j$, whereas $T$ is a count of annual blocks. Independence between years is unnecessary for this one-year quantile definition; it is needed for the geometric waiting and multi-year formulas under return period. For a discrete probability distribution, use the generalized quantile function $F^{-1}(q)=\inf \{z:F(z)\ge q\}$. Then $\mathbb{P}(M_j>z_T)\le 1/T$, with equality not guaranteed at jumps; an exact requested exceedance probability may be unattainable.

For a concrete standard model, choose the maxima convention of the Gumbel distribution, with location $\mu$ and scale $b>0$:The second formula follows by taking logarithms twice and solving $F(z_T)=1-1/T$. The maxima cumulative distribution function and its inverse are given in NIST: extreme value type I distribution, maximum case. With the deliberately chosen toy values $\mu =10$ and $b=2$ in arbitrary units, $z_{10}=14.501$, $z_{100}=19.200$, and $z_{1000}=23.815$, rounded to three decimals. In particular, $19.200$ is approximately the annual quantile at cumulative probability $0.99$; it is not a prediction of the largest observation in the next 100 years. Under independence, the probability that all those annual sample maxima stay below the exact $z_{100}$ is $0.99^{100}\approx 0.366$.

The time block matters. If a toy year contains $m$ independent daily blocks, each with exceedance probability $q$ at one fixed threshold, its annual exceedance probability is $p=1-(1-q{)}^m$, giving annual return period $1/[1-(1-q{)}^m]$. The reciprocal $1/q$ counts daily blocks. Dependence among daily exceedances requires their joint law instead; a threshold crossed on several successive days still produces only one exceedance indicator for that year's sample maximum. The peaks-over-threshold method and declustering of extremes address a different choice of extreme-event sampling.

A return level plot graphs $T$ horizontally and $z_T$ vertically, often using a logarithmic scale for $T$. To read it, choose a return period, move up to the curve, and read the corresponding return level. Conversely, a threshold's horizontal coordinate gives the reciprocal of its annual exceedance probability. The usual plot of fitted quantiles against return periods, and optional pointwise confidence intervals, are documented in Mathieu Ribatet, POT package: return level plot.

The near-straight shape has a mathematical explanation: since $-\log (1-1/T)\sim 1/T$ as $T\to \infty$, the Gumbel distribution gives $z_T=\mu +b\log T+o(1)$. Multiplying a large return period by ten therefore increases the return level by approximately $b\log 10$, about $4.605$ units for this toy curve. A curve indexed by rarity puts several thresholds on one common scale while retaining their measurement units.

For a fitted curve, a pointwise confidence interval concerns uncertainty in the estimated quantile at a specified $T$; it is distinct from a prediction interval for a future annual sample maximum, and pointwise intervals need not form a simultaneous confidence band. A far-right curve value can also be an extrapolation beyond the observation record: extending the line adds no observations to support the chosen tail model. The diagram above uses fixed toy parameters, so it shows neither estimation uncertainty nor an observational record. Its interpretation throughout assumes the stated annual block and stationarity.

For scale $\sigma >0$ and shape $\xi$, the generalized Pareto distribution on nonnegative excesses has cumulative distribution functionThe case $\xi =0$ is the continuous limit $1-e^{-y/\sigma }$. Positive $\xi$ gives a power-law tail, zero gives an exponential distribution, and negative $\xi$ gives a finite upper endpoint $-\sigma /\xi$. Its role in modelling high-threshold excesses is explained by the Pickands-Balkema-de Haan theorem.

For the positive-shape generalized Pareto distribution with shape $\psi >0$ and scale $\sigma >0$, the expected Fisher information matrix for one observation isThe orthogonalization equation for two statistical parameters gives ${\partial }_{\psi }\sigma =-\sigma /(1+\psi )$ at fixed new nuisance coordinate. Thus $\eta =\sigma (1+\psi )$ is a global interest-respecting reparametrization. In coordinates $(\psi ,\eta )$, the information is diagonal with entries $(1+\psi {)}^{-2}$ and $[{\eta }^2(1+2\psi ){]}^{-1}$.

The peaks-over-threshold method fits a generalized Pareto distribution to exceedances of a high threshold and also estimates their frequency. Together these determine tail probabilities and return levels. For dependent observations, threshold exceedances may occur in clusters; declustering of extremes can separate extreme events before fitting a model for independent events.

Declustering groups dependent threshold exceedances into extreme-event clusters and retains a representative, commonly the cluster maximum. It aims to make the events used in a peaks-over-threshold method approximately independent. The event rate must be estimated consistently with the retained clusters when computing return levels.

For a distribution in an extreme-value domain of attraction, the conditional distribution of the excess above a high threshold is asymptotically approximated by a generalized Pareto distribution with a threshold-dependent scale and a limiting shape parameter. This theorem justifies the peaks-over-threshold method. Threshold selection must balance approximation error against the reduced number of excess observations.

The standard Gumbel distribution for maxima has distribution function $\Lambda (z)=e^{-e^{-z}}$ for every real $z$. Its probability density function is the derivative of its distribution function:This article uses the maxima convention throughout. Reflecting a standard Gumbel distribution random variable $Z$ to $-Z$ gives the minima convention, with distribution function $1-\exp (-e^z)$. The two conventions have opposite directions of asymmetry, so their formulas must not be interchanged. ^[ref].

For a location $\mu \in \mathbb{R}$ and positive scale $b>0$, let $Y=\mu +bZ$. This location-scale family has cumulative distribution function and probability density functionSolving $\Lambda (z)=p$ gives the standard quantile function $q(p)=-\log (-\log p)$, $0<p<1$, and hence the standard medianThe median of $Y$ is therefore $\mu -b\log (\log 2)$. Location $\mu$ is not the median: at $x=\mu$ the cumulative distribution function is $e^{-1}\approx 0.367879$.

To derive the Gumbel distribution as an extreme-value limit, take the maximum of independent exponential random variables with rate $\lambda >0$ and setThe exact distribution function isThe finite-$n$ lower boundary is required because $M_n\ge 0$. For each fixed real $z$, eventually $z>-\log n$, and the natural logarithm of the first branch satisfiesExponentiating proves convergence in distribution to the standard Gumbel distribution. It is the familiar limit $(1-u/n{)}^n\to e^{-u}$ with $u=e^{-z}$: a simple calculation reveals a whole limiting shape. The name honors Emil Gumbel; Siegrist gives this exponential-maxima derivation and the standard Gumbel distribution formulas. ^[ref].

For the toy delay model $\lambda =1{\mathrm{s}}^{-1}$, $n=100$, the limiting Gumbel distribution approximation to $M_{100}$ has $\mu =(\log 100)\mathrm{s}$ and $b=1\mathrm{s}$. Its predicted median iscompared with the exact median $-\log (1-2^{-1/100})\mathrm{s}\approx 4.975147\mathrm{s}$. The approximation is about $0.003464$ seconds lower in this example. It is not an exact finite-sample probability distribution: the approximating Gumbel distribution assigns positive probability to negative delays, while the exact sample maximum cannot be negative.

The maximum domain of attraction of the Gumbel distribution is characterized by an exponential tail ratio on the scale of a Gumbel auxiliary function. This allows many original probability distributions to share the limiting shape, but it does not make the centering $\log n$ appropriate to all of them.

If $Y$ has cumulative distribution function $F(y)=\exp (-e^{-y})$, put $Z=e^{-Y}$. For $z>0$, the decreasing transformation gives $\mathbb{P}(Z\le z)=1-F(-\ln z)=1-e^{-z}$, while this probability is zero for $z\le 0$. Hence $Z$ has the rate-one exponential distribution. Conversely, if $Z$ has that exponential distribution, then $-\ln Z$ has the standard Gumbel distribution.

For fixed $a>0$ and $r>1$, the factor $\exp (-\lambda e^{-x})$ is a translated Gumbel distribution cumulative function, with transition near $x=L$. Integration by parts transforms the integral to an exponentially small boundary term plus $(r-1{)}^{-1}{\int }_0^{\lambda e^{-a}}{e}^{-u}(L-\log u{)}^{1-r}du$. Expanding the latter against its finite logarithmic moments gives the displayed two terms because ${\int }_0^{\infty }{e}^{-u}\log udu=-\gamma$, with $\gamma$ Euler's constant. The smooth cutoff shifts the first correction even though replacing it by a step gives the leading tail.

For the Gamma distribution with fixed positive integer shape $m$ and rate one, the survival function is $e^{-t}{\sum }_{j=0}^{m-1}{t}^j/j!$. Its sample maximum centered by $\log n+(m-1)\log \log n-\log ((m-1)!)$ has convergence in distribution to the standard Gumbel distribution, without an additional scale factor.

A positive Gumbel auxiliary function makes $\overline{F}(t+a(t)x)/\overline{F}(t)\to e^{-x}$ as $t$ approaches the right endpoint of a distribution, for every real $x$. Under the Von Mises conditions for extreme values, the reciprocal hazard function is a valid choice.

For $\alpha >0$, the standard negative Weibull distribution has distribution function ${\Psi }_{\alpha }(x)=e^{-(-x{)}^{\alpha }}$ for $x<0$ and one for $x\ge 0$. This is the probability distribution of the negative of a standard positive Weibull distribution variable.

For $\alpha >0$, the standard Fréchet distribution has distribution function ${\Phi }_{\alpha }(x)=e^{-x^{-\alpha }}$ for $x>0$ and zero for $x\le 0$. Its maximum domain of attraction consists of distribution functions with an infinite right endpoint of a distribution and a survival function with regular variation of index $-\alpha$.

A distribution function $F$ lies in the maximum domain of attraction of $G$ if $F(a_{n}x+b_n{)}^n\to G(x)$ at every continuity point of $G$ for some $a_n>0$ and real $b_n$. This expresses convergence in distribution of normalized sample maxima.

For a distribution function with tail $1-F(x)=1/(x\log x)$, let $a_n>1$ solve $a_n\log a_n=n$. Then $n(1-F(a_{n}x))\to 1/x$ for $x>0$, so the normalized sample maximum tends to the shape-one Fréchet distribution. The inequalities $n/\log n<a_n<n/(\log n-\log \log n)$ for large $n$ show $a_n\sim n/\log n$, allowing that simpler scale too.

For an absolutely continuous distribution function with positive probability density function near its right endpoint of a distribution, write $\eta =f/(1-F)$ for its hazard function. The limits $t\eta (t)\to \alpha$ at an infinite endpoint or $(x_F-t)\eta (t)\to \alpha$ at a finite endpoint imply attraction to the Fréchet distribution or negative Weibull distribution, respectively. A continuously differentiable reciprocal hazard function with $(1/\eta {)}^{\prime }\to 0$ implies attraction to the Gumbel distribution, provided the survival function tends to zero at the endpoint.