Extreme Value Theory
Advanced Statistical Methods
Modeling the Tails That Matter Most
Extreme value theory provides the mathematical framework for modeling rare, high-impact events using the GEV distribution and peaks-over-threshold methods. It quantifies return periods for events far beyond ordinary observations.
- Insurance — Estimate catastrophic loss probabilities for flood, earthquake, and hurricane risk pricing
- Structural engineering — Design buildings and bridges to withstand rare wind and seismic loads
- Finance — Model value-at-risk and expected shortfall for extreme market crashes
EVT gives you the tools to prepare for events that standard distributions underestimate.
Block Maxima Method
Peaks-Over-Threshold (POT)
Threshold Selection
import numpy as np
from scipy import stats, optimize
class ExtremeValueAnalysis:
def __init__(self, data):
self.data = np.asarray(data)
def fit_gev_mle(self):
shape, loc, scale = stats.genextreme.fit(-self.data)
return -shape, loc, scale
def fit_gpd_pot(self, threshold):
exceedances = self.data[self.data > threshold] - threshold
if len(exceedances) < 10:
raise ValueError("Too few exceedances")
shape, loc, scale = stats.genpareto.fit(exceedances)
return shape, scale, len(exceedances)
def hill_estimator(self, k=None):
sorted_data = np.sort(self.data)
n = len(sorted_data)
if k is None:
k = int(n ** 0.5)
log_ratios = np.log(sorted_data[n-k:n] / sorted_data[n-k-1])
alpha_hat = 1.0 / np.mean(log_ratios)
alpha_std = alpha_hat / np.sqrt(k)
return alpha_hat, alpha_std
def mean_residual_life(self, thresholds):
mrl = []
counts = []
for u in thresholds:
exceed = self.data[self.data > u]
if len(exceed) > 0:
mrl.append(np.mean(exceed - u))
counts.append(len(exceed))
else:
mrl.append(np.nan)
counts.append(0)
return np.array(mrl), np.array(counts)
def return_level(self, T, u, sigma, xi, n_total, n_exceed):
rate = n_exceed / n_total
z_T = u + (sigma / xi) * ((rate * T) ** xi - 1)
return z_T
def gev_quantile(self, p, shape, loc, scale):
return stats.genextreme.ppf(1-p, -shape, loc=loc, scale=scale)
def block_maxima(self, block_size):
n_blocks = len(self.data) // block_size
maxima = [np.max(self.data[i*block_size:(i+1)*block_size])
for i in range(n_blocks)]
return np.array(maxima)