mixle.analysis.extreme module

Extreme-value analysis: tails, threshold exceedances, and support endpoints.

The mean of a sample tells you nothing reliable about its rare extremes – floods, market crashes, record temperatures live in the tail, where ordinary fitting has almost no data. Extreme-value theory gives the limiting shapes that govern tails and provides estimators built only from the largest order statistics:

• gpd_fit() / peaks_over_threshold() – the Generalized Pareto Distribution for threshold exceedances (the POT method), by maximum likelihood or probability-weighted moments, with return_level() for “the once-in-m-observations level”.

• hill_estimator() / moment_estimator() – tail-index estimators from the top-k order statistics (Hill for heavy tails; the Dekkers–Einmahl–de Haan moment estimator for any tail).

• mean_residual_life() – the mean-excess plot for choosing the POT threshold (linear in the threshold once the GPD regime is reached).

• endpoint_estimator() – the finite right endpoint of a bounded support (Hall-type / GPD), generic to frontier analysis, reliability limits, and image edges.

• record_times() / n_records() – running-maximum records and their count.

class GPDFit(shape, scale, threshold, n_exceedances, n_total, method)

Bases: object

Fitted Generalized Pareto Distribution for exceedances over a threshold.

Parameters:

• shape (float)

• scale (float)

• threshold (float)

• n_exceedances (int)

• n_total (int)

• method (str)

shape

tail index xi (> 0 heavy/Pareto tail, = 0 exponential, < 0 bounded).

Type:

float

scale

scale beta (> 0).

Type:

float

threshold

the threshold u the exceedances were measured over.

Type:

float

n_exceedances / n_total

exceedance and full-sample sizes (for return levels).

method

"mle" or "pwm".

Type:

str

shape: float

scale: float

threshold: float

n_exceedances: int

n_total: int

method: str

property endpoint: float

Right endpoint of the support (finite when shape < 0, else inf).

gpd_fit(exceedances, *, threshold=0.0, method='mle', n_total=None)

Fit a Generalized Pareto Distribution to threshold exceedances.

Parameters:

• exceedances (ndarray) – the excesses x - u for observations above the threshold (all positive). If you have raw data, use peaks_over_threshold() instead.

• threshold (float) – the threshold u (stored for return-level computation).

• method (str) – "mle" (Newton/Nelder-Mead on the GPD likelihood) or "pwm" (probability-weighted moments, closed form, robust for xi < 0.5).

• n_total (int | None) – full sample size before thresholding (defaults to the number of exceedances).

Returns:

A GPDFit.

Return type:

GPDFit

peaks_over_threshold(data, threshold, *, method='mle')

Peaks-over-threshold: select exceedances above threshold and fit a GPD to the excesses.

Parameters:

• data (ndarray)

• threshold (float)

• method (str)

Return type:

GPDFit

return_level(fit, period)

POT return level: the level exceeded on average once per period observations.

x_m = u + (beta/xi) [ (m zeta_u)^xi - 1 ] with zeta_u = n_exceed/n_total the exceedance rate (period = m). For xi = 0 it reduces to u + beta log(m zeta_u).

Parameters:

• fit (GPDFit)

• period (float)

Return type:

float

hill_estimator(data, k)

Hill estimator of the tail index xi = 1/alpha from the top k order statistics.

xi_hat = (1/k) sum_{i=1}^{k} log(X_(n-i+1) / X_(n-k)) – consistent for heavy (Pareto-type, xi > 0) tails. For a Pareto tail with exponent alpha this estimates 1/alpha.

Parameters:

• data (ndarray)

• k (int)

Return type:

float

moment_estimator(data, k)

Dekkers–Einmahl–de Haan moment estimator of the extreme-value index (any tail sign).

Generalises Hill to xi of either sign by combining the first two log-moments of the top k exceedances; works for heavy, light, and bounded (xi < 0) tails.

Parameters:

• data (ndarray)

• k (int)

Return type:

float

mean_residual_life(data, thresholds)

Mean-excess (mean-residual-life) function for POT threshold selection.

e(u) = mean(X - u | X > u). Over a range where the GPD fits, e(u) is approximately linear in u (slope xi/(1-xi)); the lowest threshold from which the plot is linear is the choice.

Returns:

{'threshold', 'mean_excess', 'n_exceed'}.

Parameters:

• data (ndarray)

• thresholds (ndarray)

Return type:

dict[str, ndarray]

endpoint_estimator(data, k, *, method='gpd')

Estimate the finite right endpoint of a bounded support (frontier / boundary estimation).

Fits a GPD to the top k exceedances over X_(n-k); when the tail index xi is negative the support is bounded and the endpoint is x+ = X_(n-k) - beta/xi (which, by the GPD support constraint, necessarily exceeds the observed maximum). Generic to econometric frontier analysis, reliability limits, and image-edge localisation. Returns inf if the estimated tail is unbounded (xi >= 0).

Parameters:

• data (ndarray) – the sample.

• k (int) – number of upper order statistics (exceedances) used.

• method (str) – "gpd" – GPD-MLE endpoint.

Returns:

The estimated right endpoint (inf if unbounded).

Return type:

float

record_times(data)

Indices at which a new running maximum (upper record) occurs, including the first observation.

Parameters:

data (ndarray)

Return type:

ndarray

n_records(data)

Number of upper records. For an i.i.d. sequence of length n the expectation is H_n.

Parameters:

data (ndarray)

Return type:

int