Extreme Value Theory and Models

Lecture

2023-11-06

Reading

This lecture follows Chapter 3 of Coles (2001) fairly closely

Introduction

Today

  1. Introduction

  2. Inference

  3. Wrapup

Gumbel, Fréchet, and Weibull distributions

Code 
dists = [GeneralizedExtremeValue(0, 2.5, ξ) for ξ in [-0.5, 0, 0.5]]
names = ["Weibull", "Gumbel", "Fréchet"]
colors = [:blue, :red, :green]

p = plot()
for (dist, name, c) in zip(dists, names, colors)
    plot!(p, dist; label=name, color=c, xlims=(-10, 25))
    if dist.ξ > 0
        lb = dist.μ − dist.σ / dist.ξ
        scatter!(p, [lb], [0]; label=nothing, color=c)
    elseif dist.ξ < 0
        ub = dist.μ − dist.σ / dist.ξ
        scatter!(p, [ub], [0]; label=nothing, color=c)
    end
end
p

Generalized Extreme Value distribution

  • Instead of heuristically picking which of the 3 distributions to use, we can use the Generalized Extreme Value distribution (GEV)
  • Account for uncertainty in choice of extreme value distribution

Quantiles

Also known as return levels. Following Coles (2001) notation, the level \(z_p\) is exceeded with probability \(p\)

\[ z_p = \begin{cases} \mu - \frac{\sigma}{\xi} \left[ 1 - y_p^{-\xi} \right], & \quad \text{for} \quad \xi \neq 0 \\ \mu - \sigma \log y_p, & \quad \text{for} \quad \xi = 0 \\ \end{cases} \] where \(y_p = -\log(1-p)\) for concise notation

Code 
y_plot = 0:0.1:7
p = plot(; xlabel=L"$y_p$", ylabel=L"$z_p$", legend=:topleft)
for (dist, name) in zip(dists, names)
    plot!(p, y_plot, quantile(dist, 1 .- exp.(-y_plot)); label=name)
end
p

Extremal Types Theorem

  • See Coles (2001) section 3.1.4 for a proof
  • Justifies use of GEV for modeling the distribution of maxima of long sequences
  • Numerical examples here are based on analytic examples in 3.1.5

Minima

Statistical consultants hate this one simple trick:

Define \(Y_i = -X_i\) and use the GEV distribution to model the maxima of \(Y_i\)

Inference

Today

  1. Introduction

  2. Inference

  3. Wrapup

Bias-variance trade-off

  • Block size is a choice!
  • Small blocks: the limit model, leading to bias in estimation and extrapolation
  • Large blocks: few block maxima, leading to large estimation variance
  • Default choice in environmental applications is typically 1 year
    • Key assumption: each block from the same distribution (IID)
    • If you choose 3 months, need to account for differing distirbutions of, eg, temperature in summer and winter

Overview of inference methods

  1. Graphical techniques
  2. Moment techniques
  3. Order statistics based methods
  4. Likelihood approaches

Example data

data = Extremes.dataset("portpirie")
plot(
    data[!, :Year], data[!, :SeaLevel]; xlabel="Year", ylabel="Sea level (m)", legend=false
)

Maximum likelihood estimation I

Conceptually straightforward: \[ \begin{aligned} \ell(\mu, \sigma, \xi) = -m \log \sigma - (1 + 1/\xi) &\sum_{i=1}^{m} \log \left[ 1 + \xi \left( \frac{z_i - \mu}{\sigma} \right) \right]\\ &- \sum_{i=1}^{m} \left[ 1 + \xi \left( \frac{z_i - \mu}{\sigma} \right) \right]^{-1/\xi}, \end{aligned} \] provided that \[ 1 + \xi \left( \frac{z_i - \mu}{\sigma} \right) > 0 \] when \(\xi=0\) the equation is a bit different \[ \ell(\mu, \sigma) = -m \log \sigma - \sum_{i=1}^{m} \left( \frac{z_i - \mu}{\sigma} \right) - \sum_{i=1}^{m} \exp \left\{ - \left( \frac{z_i - \mu}{\sigma} \right) \right\}. \]

Tip

I don’t need you to know these equations

Maximum likelihood estimation II

Numerical implementations need to consider

  • Use the right function around \(\xi=0\)
  • Return a value of \(-\inf\) if bounds constraint is violated
  • Possible to derive approximations of uncertainty
    • Can approximate uncertainty in return levels

Maximum likelihood estimation III

gev_mle = gevfit(data.SeaLevel)
MaximumLikelihoodAbstractExtremeValueModel
model :
    BlockMaxima{GeneralizedExtremeValue}
    data :      Vector{Float64}[65]
    location :  μ ~ 1
    logscale :  ϕ ~ 1
    shape :     ξ ~ 1

θ̂  :   [3.874750223091266, -1.6192723640210762, -0.05010719929448139]

Bayesian inference

gev_bayes = gevfitbayes(data.SeaLevel; niter=10_000, warmup=2_000)
BayesianAbstractExtremeValueModel
model :
    BlockMaxima{GeneralizedExtremeValue}
    data :      Vector{Float64}[65]
    location :  μ ~ 1
    logscale :  ϕ ~ 1
    shape :     ξ ~ 1

sim :
    MambaLite.Chains
    Iterations :        2001:10000
    Thinning interval : 1
    Chains :        1
    Samples per chain : 8000
    Value :         Array{Float64, 3}[8000,3,1]

Profile likelihood

  1. For a range of values of the shape parameter \(\xi\):
    1. Maximize the likelihood with respect to \(\mu\) and \(\sigma\)
    2. Compute the likelihood
  2. Plot the (log) likelihood as a function of \(\xi\)

Moment matching

  • Moment: mean, variance, skewness, kurtosis, etc
  • Moment matching: estimate the moments of the samples then match them to the theoretical moments of the distribution
  • Probability-weighted moments (Hosking et al., 1985): weight data points by their cumulative probability (i.e., more weight for biggest values)
  • \(L\) moments (Hosking, 1990): Linear combinations of order statistics, less sensitive to outliers

Probability-weighted moments II

gev_pwm = gevfitpwm(data.SeaLevel)
pwmAbstractExtremeValueModel
model :
    BlockMaxima{GeneralizedExtremeValue}
    data :      Vector{Float64}[65]
    location :  μ ~ 1
    logscale :  ϕ ~ 1
    shape :     ξ ~ 1

θ̂  :   [3.8731723562720766, -1.5932320395836068, -0.051477125862911276]

Uncertainty

cint(gev_mle)
3-element Vector{Vector{Float64}}:
 [3.820004234825991, 3.929496211356541]
 [-1.819669858589598, -1.4188748694525544]
 [-0.24268345866324326, 0.1424690600742805]
cint(gev_bayes)
3-element Vector{Vector{Float64}}:
 [3.817789629088061, 3.928899159987298]
 [-1.7886588596429014, -1.3888832120165522]
 [-0.2119157154228395, 0.1597075738279036]
cint(gev_pwm)
3-element Vector{Vector{Float64}}:
 [3.82219678818851, 3.9329570198054546]
 [-1.8238458137877533, -1.4129406745797473]
 [-0.23891541608078742, 0.09151391904454692]

Return levels

function weibull_plot_pos(y)
    N = length(y)
    ys = sort(y; rev=false) # sorted values of y
    nxp = xp = [r / (N + 1) for r in 1:N] # exceedance probability
    xp = 1 .- nxp
    return xp, ys
end

function plot_rl_extremes(fit, obs)
    minval = 1 + 1 / (length(obs) + 1)
    return_period = 10 .^ range(log10(minval), log10(250); length=50)
    x_ticks = [2, 5, 10, 20, 50, 100, 250]

    p = plot(;
        xscale=:log10,
        xticks=(x_ticks, string.(x_ticks)),
        xlabel="Return Period",
        ylabel="Return Level",
        legend=false,
    )

    rl = zero(return_period)
    ub = zero(return_period)
    lb = zero(return_period)

    for (i, rt) in enumerate(return_period)
        r = returnlevel(fit, rt) # special object
        rl[i] = mean(r.value)
        lb[i], ub[i] = cint(r)[1]
    end

    plot!(p, return_period, rl)
    plot!(p, return_period, lb; fillrange=ub, fillalpha=0.35, linewidth=0)

    xp, ys = weibull_plot_pos(obs)
    scatter!(p, 1 ./ xp, ys)

    return p
end

Comparison

Code 
p1 = plot_rl_extremes(gev_mle, data.SeaLevel)
title!(p1, "MLE")

p2 = plot_rl_extremes(gev_bayes, data.SeaLevel)
title!(p2, "Bayesian Inference")

p3 = plot_rl_extremes(gev_pwm, data.SeaLevel)
title!(p3, "Probability-weighted moments")

plot(p1, p2, p3; link=:all)

Wrapup

Today

  1. Introduction

  2. Inference

  3. Wrapup

Limitations

  • Nonstationarity
  • Spatial dependence and multivariate extremes
  • Small samples
  • UQ
  • Model checking / validation

References

Coles, S. (2001). An introduction to statistical modeling of extreme values. London ; Springer.
Hosking, J. R. M. (1990). L-Moments: Analysis and Estimation of Distributions Using Linear Combinations of Order Statistics. Journal of the Royal Statistical Society. Series B (Methodological), 52(1), 105–124. Retrieved from https://www.jstor.org/stable/2345653
Hosking, J. R. M., Wallis, J. R., & Wood, E. F. (1985). Estimation of the Generalized Extreme-Value Distribution by the Method of Probability-Weighted Moments. Technometrics, 27(3), 251–261. https://doi.org/10.1080/00401706.1985.10488049