extremeloss

Extreme-value tail estimation for large claims: peaks-over-threshold with the generalized Pareto distribution (GPD), tail analytics, and large-claim loading. Matplotlib is optional via the plot extra, and lossmodels integration via the splice extra.

The package covers threshold selection, GPD estimation, and the analytics built on a fitted tail — return levels, exceedance probabilities, and excess-layer charges — plus integration helpers for splicing an empirical body to a fitted tail. It composes with lossmodels and the pooling primitives in actuarialpy.

Empirical tail measures follow the ecosystem VaR/TVaR estimators, and simulation-backed quantities accept the shared rng argument — see Conventions.

See the API reference below for the full surface; each object’s docstring carries its own usage.

Threshold selection

Scan candidate thresholds before committing to one. threshold_diagnostic_table returns a ThresholdScan with the exceedance counts, GPD parameter estimates, and mean excesses at each candidate; mean_excess returns the same diagnostics for plotting (matplotlib via the plot extra):

import numpy as np
import extremeloss as xl

scan = xl.threshold_diagnostic_table(
    losses, np.quantile(losses, [0.90, 0.925, 0.95, 0.975])
)
scan.thresholds, scan.n_exceedances, scan.xi, scan.beta, scan.mean_excess

Take the lowest threshold above which the shape estimate \(\hat\xi\) is stable and the mean-excess function is roughly linear — the classic peaks-over-threshold trade between bias (threshold too low) and variance (too few exceedances).

Fitting the tail

fit_pot(data, threshold) selects the exceedances of the threshold and fits the GPD to their excesses; fit_gpd(excesses) is the same fit when the excesses are already in hand:

fit = xl.fit_pot(losses, threshold=u)
fit.xi, fit.beta, fit.exceedance_fraction

# unconditional, ground-up tail metrics — the fit carries the exceedance rate
xl.gpd_var(0.995, fit.threshold, fit.xi, fit.beta, fit.exceedance_fraction)
xl.gpd_tvar(0.995, fit.threshold, fit.xi, fit.beta, fit.exceedance_fraction)
xl.return_level(200, fit)        # the 1-in-200 claim == gpd_var(0.995, ...)

tail = xl.GPDTail.from_fit(fit)  # conditional law on [u, inf), for splicing

GPDFit quotes unconditional quantities; GPDTail is the conditional excess law a splice consumes. The distinction, and the \((\xi, \beta)\) parameterization, are pinned in Conventions. Sampling accepts the shared rng argument.

Uncertainty and return levels

fit_gpd and fit_pot now populate GPDFit.covariance with the observed-information covariance of (xi, beta) (None when the information matrix is not positive definite — the MLE is irregular for xi <= -1/2, and a covariance that means nothing is worse than none). That uncertainty flows into the two places it matters most.

Threshold selection with error bands. The scan reports xi and the modified scale beta* = beta - xi*u — the quantity that is actually constant above a valid threshold (raw beta drifts linearly in u even under a perfect GPD) — with standard errors, so the question becomes “where do these flatten within their bands”:

scan = el.threshold_diagnostic_table(losses, thresholds)
# ... xi | xi_se | modified_scale | modified_scale_se | n_exceedances

Return levels with confidence intervals. The loss exceeded once per T periods, with delta-method intervals over (zeta_u, xi, beta) including the binomial variance of the exceedance rate (Coles §4.3.3):

el.gpd_return_level(fit, [10, 50, 200], observations_per_period=1200)
# return_period | return_level | se | ci_low | ci_high

Point estimates agree exactly with the scalar fit.return_level(n), which remains the shorthand. A fitted tail also exposes sf and a closed-form mean_excess, so it plugs directly into ratingmodels.pooling_charge_from_severity.

The same treatment covers block maxima: fit_gev populates GEVFit.covariance, and gev_return_level carries delta-method bands (periods in blocks). Switching between the POT and block-maxima workflows no longer changes how much uncertainty you get to see.

Did the fit succeed? qq_points / pp_points give the diagnostic point sets headlessly, and parametric_bootstrap_gof tests the fit with p-values that are honest about estimation – every bootstrap replicate is refit before scoring, because comparing an estimated-parameter statistic against the known-parameter null is anti-conservative:

el.parametric_bootstrap_gof(fit, losses, n_boot=500)
# {"ks": ..., "ks_pvalue": ..., "ad": ..., "ad_pvalue": ...}

With the plot extra, plot_gpd_diagnostics(fit, losses) and plot_gev_diagnostics(fit, maxima) draw the Coles four-panel exhibit – probability, quantile, return level with its band and empirical points, density.

Splicing onto a fitted body

The handoff runs one direction: extremeloss fits the tail and returns a lossmodels.SplicedSeverity, so every downstream consumer holds the same severity class whether or not it carries an EVT tail:

sev = xl.splice_gpd_tail(body, fit)               # weight from the fit
sev = xl.fit_spliced_gpd(body, losses, threshold=u)  # weight from the data

Splicing is mass-matching, never density-continuity; pass weight= to override the empirical body mass.

From simulated portfolios

A risksim.SimulationResult feeds straight back into the tail toolkit:

res = portfolio.simulate(100_000, contract=treaty, rng=7)   # risksim
xl.tail_summary_from_risksim(res)
# n / mean / std / min / max and VaR/TVaR at (0.95, 0.99, 0.995)

API reference

class BootstrapResult(estimate: float, bootstrap_estimates: ndarray, method: str = 'percentile', ci: tuple[float, float] | None = None, stderr: float | None = None, alpha: float | None = None)

Bases: object

Bootstrap uncertainty summary for a scalar statistic.

class GEVFit(xi: float, loc: float, scale: float, n_blocks: int, block_size: int | None = None, fit_method: str = 'mle', covariance: ndarray | None = None)

Bases: object

Fitted generalized extreme value distribution for block maxima.

property se

Standard errors of (xi, loc, scale), or None.

class GPDFit(threshold: float, xi: float, beta: float, exceedance_fraction: float, n_exceedances: int, fit_method: str = 'mle', covariance: ndarray | None = None)

Bases: object

Fitted generalized Pareto distribution above a threshold.

mean_excess(d: float) → float

Mean excess E[X - d | X > d] above the threshold, closed form.

For the GPD, \(e(d) = (\beta + \xi\,(d - u)) / (1 - \xi)\) for \(d \ge u\) and \(\xi < 1\); infinite for \(\xi \ge 1\) (the tail has no mean). Below the threshold the GPD says nothing – a ValueError rather than a silent extrapolation.

property se

Standard errors of (xi, beta), or None without covariance.

sf(x: float) → float

Unconditional survival P(X > x); alias of tail_probability.

Half of the small tail protocol (with mean_excess()) consumed by e.g. ratingmodels.pooling_charge_from_severity.

class GPDTail(threshold: float, xi: float, beta: float)

Bases: object

Conditional generalized-Pareto tail on [threshold, inf), for splicing.

A thin distribution object wrapping scipy.stats.genpareto with shape xi, scale beta and location threshold. Unlike GPDFit (which carries the exceedance rate and reports unconditional VaR/TVaR), this represents the tail conditional on X > threshold: cdf(threshold) == 0 and the density integrates to one over [threshold, inf). It exposes the cdf / pdf / sample / quantile / mean / variance interface that lossmodels.SplicedSeverity consumes as a tail. Moments raise when they do not exist (xi >= 1 for the mean, xi >= 1/2 for the variance), mirroring heavy-tailed severities.

classmethod from_fit(fit: GPDFit) → GPDTail

Build a conditional tail from a fitted GPDFit.

class TailEstimateResult(estimate: float, method: str, stderr: float | None = None, ci: tuple[float, float] | None = None, n: int | None = None, effective_n: float | None = None, threshold: float | None = None, quantile: float | None = None, diagnostics: dict[str, ~typing.Any] = <factory>)

Bases: object

Container for tail-estimation results.

class ThresholdScan(thresholds: ndarray, mean_excess: ndarray, xi: ndarray, beta: ndarray, n_exceedances: ndarray, xi_se: ndarray | None = None, modified_scale: ndarray | None = None, modified_scale_se: ndarray | None = None)

Bases: object

Threshold-diagnostic results across a grid of thresholds.

modified_scale is \(\beta^* = \beta - \xi u\): above any valid GPD threshold both \(\xi\) and \(\beta^*\) are constant in u, so the threshold-selection question is “where do these flatten within their confidence bands” – which is why the standard-error columns exist. Raw \(\beta\) drifts linearly in u even under a perfect GPD and diagnoses nothing.

block_return_level(period: float, fit: GEVFit) → float

The period-block GEV return level (the quantile at 1 - 1/period).

Thin wrapper over GEVFit.return_level() with the domain check that period exceeds one block. Periods are in blocks: with annual maxima, period=100 is the 100-year level. For confidence intervals use gev_return_level().

bootstrap_statistic(data, statistic, *, n_resamples: int = 1000, alpha: float = 0.05, random_state: int | Generator | None = None) → BootstrapResult

Nonparametric bootstrap of an arbitrary statistic, with a percentile CI.

Resamples data with replacement n_resamples times, evaluates statistic on each resample, and summarizes the bootstrap distribution with a point estimate (the statistic on the original sample), a percentile confidence interval, and a bootstrap standard error. The interval is the equal-tailed percentile interval with endpoints at the alpha/2 and 1 - alpha/2 quantiles of the resampled statistics.

Parameters:

• data (array-like) – The sample to resample.

• statistic (callable) – A function mapping a 1-D array to a float, e.g. np.mean or a tail measure. It is called once on the original data for the point estimate and once per resample.

• n_resamples (int, optional) – Number of bootstrap resamples (default 1000). Must be positive.

• alpha (float, optional) – Significance level for the two-sided percentile interval (default 0.05, i.e. a 95% interval). Must lie in (0, 1).

• random_state (int, numpy.random.Generator, or None, optional) – Seed or generator controlling the resampling. An existing Generator is used as given; anything else is passed to numpy.random.default_rng(). Fix it for reproducible intervals.

Returns:

Carries estimate (statistic on the original sample), bootstrap_estimates (the array of resampled values), ci (the (low, high) percentile bounds), stderr (standard deviation of the resampled values, ddof=1), method ("percentile"), and alpha.

Return type:

BootstrapResult

Raises:

ValueError – If n_resamples is not positive or alpha is not in (0, 1).

Notes

The percentile interval is simple and distribution-free but can under-cover for strongly skewed sampling distributions – the same caveat that applies to Wald tail intervals elsewhere in the package.

See also

bootstrap_var

Bootstrap CI for value-at-risk.

bootstrap_tvar

Bootstrap CI for tail value-at-risk.

bootstrap_tail_probability

Bootstrap CI for an exceedance probability.

bootstrap_tail_probability(losses, threshold: float, **kwargs) → BootstrapResult

Bootstrap confidence interval for an exceedance probability.

Convenience wrapper over bootstrap_statistic() with the statistic fixed to the empirical probability \(P(X > u)\) of exceeding threshold.

Parameters:

• losses (array-like) – Loss sample to resample.

• threshold (float) – The level \(u\) whose exceedance probability is bootstrapped.

• **kwargs – Forwarded to bootstrap_statistic() (n_resamples, alpha, random_state).

Returns:

As returned by bootstrap_statistic(), for the exceedance probability.

Return type:

BootstrapResult

See also

bootstrap_statistic

The general routine this delegates to.

bootstrap_tvar(losses, q: float, **kwargs) → BootstrapResult

Bootstrap confidence interval for tail value-at-risk.

Convenience wrapper over bootstrap_statistic() with the statistic fixed to the empirical tail value-at-risk at level q – the mean loss in the worst 1 - q fraction of outcomes (also called CVaR or expected shortfall).

Parameters:

• losses (array-like) – Loss sample to resample.

• q (float) – TVaR level in (0, 1), e.g. 0.99 for the 99% TVaR.

• **kwargs – Forwarded to bootstrap_statistic() (n_resamples, alpha, random_state).

Returns:

As returned by bootstrap_statistic(), for the TVaR.

Return type:

BootstrapResult

See also

bootstrap_var

Bootstrap CI for the value-at-risk at the same level.

bootstrap_statistic

The general routine this delegates to.

bootstrap_var(losses, q: float, **kwargs) → BootstrapResult

Bootstrap confidence interval for value-at-risk (a loss quantile).

Convenience wrapper over bootstrap_statistic() with the statistic fixed to the empirical value-at-risk at level q – the q-quantile of the loss distribution.

Parameters:

• losses (array-like) – Loss sample to resample.

• q (float) – VaR level in (0, 1), e.g. 0.99 for the 99% VaR.

• **kwargs – Forwarded to bootstrap_statistic() (n_resamples, alpha, random_state).

Returns:

As returned by bootstrap_statistic(), for the VaR.

Return type:

BootstrapResult

See also

bootstrap_tvar

Bootstrap CI for the tail value-at-risk at the same level.

bootstrap_statistic

The general routine this delegates to.

component_tail_metrics(result, *, q: float = 0.99, threshold: float | None = None) → dict[str, dict[str, float]]

Per-component empirical VaR/TVaR from a risksim result’s component losses.

Reads the (n_sims, n_components) array result.component_losses (names from component_names when present) and returns, per component, the empirical VaR and TVaR at level q, plus the exceedance probability of threshold when one is given. These are marginal, per-component metrics – diversified totals belong to the portfolio-level summary.

Returns:

{name: {"var", "tvar"[, "exceedance_probability"]}}.

Return type:

dict of str -> dict of str -> float

effective_sample_size(weights) → float

Kish effective sample size of a set of importance weights.

\(\mathrm{ESS} = 1 / \sum_i \bar w_i^2\) for normalized weights \(\bar w_i\): equal to n when all weights are equal, approaching one as a single weight dominates. The headline diagnostic for weight degeneracy – the importance-sampling standard errors in this module scale by \(\sqrt{\mathrm{ESS}}\), not \(\sqrt{n}\).

empirical_tvar(losses, q)

Empirical TVaR via the average-quantile (Acerbi-Tasche) plug-in.

Implements TVaR_q = (1/(1-q)) * integral_q^1 VaR_u du exactly on the empirical distribution: with sorted losses and k = ceil(n*q),

TVaR_q = [ sum_{i>k} x_(i) + x_(k) * (k - n*q) ] / (n * (1 - q)).

Correct in the presence of ties/atoms and always >= empirical_var. Ecosystem-standard estimator shared with risksim and lossmodels.

q may be scalar (returns float) or array-like (returns array).

empirical_var(losses, q)

Empirical VaR: the order statistic x_(ceil(n*q)).

Implements VaR_q = inf{x : F(x) >= q} on the empirical distribution (identical to np.quantile(..., method="inverted_cdf")). This is the ecosystem-standard estimator shared with risksim and lossmodels; it serves as the crude-Monte-Carlo baseline that the variance-reduced estimators in this subpackage are compared against.

q may be scalar (returns float) or array-like (returns array).

estimate_exceedance_curve_is(losses, weights, thresholds) → dict[str, ndarray]

Weighted exceedance curve: self-normalized P(X > u) over a threshold grid.

The importance-sampling analogue of the empirical exceedance curve – each probability is the normalized-weight mass strictly above the threshold. Returns {"thresholds", "probabilities"} as arrays. Point estimates only; for a standard error at a single threshold use estimate_tail_probability_is().

estimate_mean_is(values, weights, *, alpha: float = 0.05) → TailEstimateResult

Self-normalized importance-sampling estimate of a mean, with a Wald CI.

\(\hat\mu = \sum_i \bar w_i x_i\) for normalized weights \(\bar w_i\). The standard error uses the weighted second moment about \(\hat\mu\) scaled by the effective sample size (Kish ESS), not n – degenerate weights widen the interval as they should.

Parameters:

• values (array-like) – Sampled values and their importance weights (normalization not required); equal lengths.

• weights (array-like) – Sampled values and their importance weights (normalization not required); equal lengths.

• alpha (float, optional) – Two-sided Wald level (default 0.05, a 95% interval).

Returns:

estimate, stderr, ci, n, effective_n, and the weight diagnostics.

Return type:

TailEstimateResult

estimate_tail_probability(data, threshold: float, *, size: int | None = None, alpha: float = 0.05) → TailEstimateResult

Estimate P(X > threshold) from simulated or observed losses.

estimate_tail_probability_cmc(conditional_probabilities, *, threshold: float | None = None, alpha: float = 0.05) → TailEstimateResult

Estimate an exceedance probability from conditional exceedance probabilities.

Parameters:

• conditional_probabilities – Samples of P(X > u | Y) from a conditioning variable Y.

• threshold – Optional tail threshold associated with the conditional probabilities.

estimate_tail_probability_is(losses, weights, threshold: float, *, alpha: float = 0.05) → TailEstimateResult

Self-normalized importance-sampling estimate of \(P(X > u)\), with a Wald CI.

The weighted exceedance-indicator mean \(\sum_i \bar w_i \, 1\{x_i > u\}\), with the standard error scaled by the effective sample size exactly as in estimate_mean_is(). This is the payoff case for importance sampling: a proposal tilted into the tail makes rare exceedances common, so far fewer simulations are needed for a given precision. The raw exceedance count rides in diagnostics.

Parameters:

• losses (array-like) – Sampled losses and their importance weights; equal lengths.

• weights (array-like) – Sampled losses and their importance weights; equal lengths.

• threshold (float) – The level \(u\) whose exceedance probability is estimated.

• alpha (float, optional) – Two-sided Wald level (default 0.05).

Returns:

estimate, stderr, ci, n, effective_n, threshold, and the weight diagnostics.

Return type:

TailEstimateResult

estimate_tvar(data, q: float, *, size: int | None = None, alpha: float = 0.05) → TailEstimateResult

Empirical tail value-at-risk with a normal-approximation confidence interval.

The ecosystem empirical TVaR at level q wrapped in a TailEstimateResult. The standard error treats the observations at or above the empirical VaR as a sample and uses their standard deviation over \(\sqrt{n_{\text{tail}}}\) – a first-order approximation that ignores the variability of the VaR threshold itself, so read the interval as indicative when the tail sample is small. data may be an array or a lossmodels-style model to sample (pass size); the tail sample size rides in diagnostics.

estimate_tvar_cmc(conditional_tail_expectations, *, q: float, threshold: float | None = None, alpha: float = 0.05) → TailEstimateResult

Estimate TVaR from conditional expectations of tail losses.

conditional_tail_expectations should contain draws of E[X | X >= VaR_q, Y] or another conditionally unbiased TVaR contribution.

estimate_tvar_is(losses, weights, q: float) → TailEstimateResult

Weighted TVaR under the ecosystem average-quantile convention.

Implements the weighted Acerbi-Tasche plug-in for TVaR_q = (1/(1-q)) * integral_q^1 VaR_u du: with values sorted, weights normalized, cumulative weights W_i and k the weighted-VaR index,

TVaR_q = [ sum_{i>k} w_i x_i + x_k (W_k - q) ] / (1 - q).

The atom at VaR contributes only the weight mass above level q, which keeps the estimator coherent with ties and makes it reduce exactly to empirical_tvar when all weights are equal.

estimate_var(data, q: float, *, size: int | None = None, alpha: float = 0.05) → TailEstimateResult

Empirical value-at-risk of simulated or observed losses, as a result object.

The ecosystem lower-quantile VaR wrapped in a TailEstimateResult with n and the quantile recorded. data may be an array of losses or a lossmodels-style model to sample (pass size). No standard error is attached for the plain empirical quantile; alpha is accepted for signature symmetry with estimate_tvar().

estimate_var_is(losses, weights, q: float) → TailEstimateResult

Weighted value-at-risk: the q-quantile under the importance weights.

The weighted analogue of the ecosystem VaR convention \(\inf\{x : F(x) \ge q\}\): with values sorted and weights normalized, the first value whose cumulative weight reaches q. Reduces exactly to the empirical VaR under equal weights. No standard error is attached (quantile standard errors under importance sampling require a density estimate); the effective sample size and weight diagnostics ride along.

estimate_var_tvar_is(losses, weights, q: float) → dict[str, TailEstimateResult]

Weighted VaR and TVaR at level q in one call.

Returns {"var": estimate_var_is(...), "tvar": estimate_tvar_is(...)} – the pair share the sorting and weighting conventions, so the TVaR always sits at or above the VaR.

exceedance_probability(losses, threshold: float) → float

Empirical exceedance probability: the fraction of losses strictly above threshold.

The plug-in estimator mean(losses > threshold) for P(X > threshold). Strict inequality, matching the exceedance convention used across the package.

extract_exceedances(data, threshold: float) → ndarray

Excesses of the data over a threshold (peaks-over-threshold data).

Returns the positive excesses \(X_i - u\) for every observation above the threshold \(u\) – the exceedance data on which a generalized Pareto tail is fitted. Note the return is excesses (measured from the threshold), not the raw exceeding values, matching the convention of the GPD fitters, which take excesses with location fixed at zero.

Parameters:

• data (array-like) – Observations to threshold.

• threshold (float) – The threshold \(u\). Only strictly-exceeding observations (x > threshold) contribute.

Returns:

The excesses x - threshold for all x > threshold, in the original data order.

Return type:

numpy.ndarray

Raises:

ValueError – If threshold is invalid, or if no observation exceeds it (an empty exceedance set cannot support a tail fit).

See also

fit_gpd

Fit a GPD to a set of excesses.

extreme_loss_summary(losses, *, thresholds=None, quantiles=(0.95, 0.99, 0.995)) → dict[str, object]

One-call tail summary of a loss sample: moments, VaR/TVaR table, exceedances.

Returns n / mean / std / min / max plus a var_tvar row per quantile – empirical VaR, TVaR, and their tail ratio, under the ecosystem estimators; pass thresholds to add the empirical exceedance_curve. This is the summary tail_summary_from_risksim() produces for a simulation result.

Parameters:

• losses (array-like) – The loss sample.

• thresholds (array-like, optional) – Grid for the exceedance curve; omitted when None.

• quantiles (tuple of float, optional) – Levels for the VaR/TVaR rows (default (0.95, 0.99, 0.995)).

Returns:

Summary statistics as above.

Return type:

dict

fit_block_maxima(data, block_size: int, method: str = 'mle', *, drop_last: bool = True) → GEVFit

Block the data and fit the GEV in one call.

Equivalent to fit_gev(make_blocks(data, block_size, drop_last=drop_last), block_size=block_size); arguments as for those two functions.

fit_gev(block_maxima, method: str = 'mle', *, block_size: int | None = None) → GEVFit

Fit a generalized extreme value distribution to block maxima.

Maximum-likelihood GEV fit in the package’s parameterization – xi positive for heavy tails (SciPy’s genextreme shape is c = -xi). The observed-information covariance of (xi, loc, scale) is attached when the information matrix is positive definite (None otherwise; the MLE is irregular for \(\xi \le -1/2\)).

Parameters:

• block_maxima (array-like) – The per-block maxima (see make_blocks()); at least two.

• method (str, optional) – Only "mle" is currently supported.

• block_size (int, optional) – Recorded on the result for bookkeeping; not used by the fit.

Returns:

Fitted (xi, loc, scale) with n_blocks and, when available, the parameter covariance.

Return type:

GEVFit

See also

fit_block_maxima

Blocking and fitting in one call.

gev_return_level

Return levels with delta-method intervals.

fit_gpd(excesses, threshold: float = 0.0, method: str = 'mle') → GPDFit

Fit a generalized Pareto distribution to excess losses.

fit_pot(data, threshold: float, method: str = 'mle') → GPDFit

Peaks-over-threshold fit: a GPD for the excesses of data over threshold.

Selects the observations strictly above the threshold \(u\), fits a generalized Pareto distribution to their excesses \(x - u\) by maximum likelihood (location fixed at zero), and returns a GPDFit carrying the empirical exceedance rate n_exceedances / n – so the fit quotes unconditional, ground-up tail quantities (sf, VaR/TVaR, return levels) out of the box. The observed-information covariance of (xi, beta) is attached when the information matrix is positive definite (None otherwise; the MLE is irregular for \(\xi \le -1/2\)).

Parameters:

• data (array-like) – Ground-up loss observations (the full sample, not pre-extracted excesses – for excesses already in hand use fit_gpd()).

• threshold (float) – The POT threshold \(u\). Choose it with threshold_diagnostic_table().

• method (str, optional) – Only "mle" is currently supported.

Returns:

Fitted (xi, beta) above threshold, with exceedance_fraction, n_exceedances, and (when available) the parameter covariance.

Return type:

GPDFit

Raises:

ValueError – If no observation exceeds the threshold, or method is not "mle".

See also

fit_gpd

The same fit when the excesses are already extracted.

threshold_diagnostic_table

Scan candidate thresholds before committing.

fit_pot_from_lossmodel(model, *, size: int, threshold: float)

Sample a lossmodels-style severity and POT-fit its tail in one call.

Draws size losses from model (anything sample_lossmodel() accepts) and returns fit_pot(losses, threshold=threshold) – how a fitted severity’s tail reads through the peaks-over-threshold lens.

fit_spliced_gpd(body, data, *, threshold: float, weight: float | None = None)

Fit a GPD tail above threshold (peaks-over-threshold) and splice it onto body, returning a lossmodels.SplicedSeverity.

Parameters:

• body (severity model) – Any fitted body severity (e.g. a lossmodels Lognormal).

• data (array-like) – Loss sample used to fit the tail and, by default, to set the body mass.

• threshold (float) – Peaks-over-threshold cutoff u.

• weight (float, optional) – Body mass P(X <= u). Defaults to 1 - exceedance_fraction from the POT fit (the empirical fraction at or below the threshold), consistent with the fitted exceedance rate.

• package. (Requires the lossmodels)

splice_gpd_tail(body, fit, *, weight: float | None = None)

Splice an already-fitted GPD tail (a GPDFit) onto body.

Returns a lossmodels.SplicedSeverity whose body is body (any fitted body severity) and whose tail is the conditional GPD of fit above its threshold. The mixing weight defaults to the body mass implied by the fit, 1 - fit.exceedance_fraction (i.e. P(X <= threshold)).

gev_return_level(fit, return_periods, confidence_level: float = 0.95)

Block-maxima return levels with confidence intervals.

The T-block return level is the GEV quantile at 1 - 1/T – identically what GEVFit.return_level() computes for a single period. Confidence intervals are by the delta method over (xi, loc, scale) (Coles, 2001, section 3.3.3); they require the fit to carry a parameter covariance (populated by fit_gev() when the information matrix is positive definite). Periods are in blocks: with annual maxima, T = 100 is the 100-year level.

Returns:

return_period, return_level, se, ci_low, ci_high (se/bounds are nan without a covariance).

Return type:

dict of str -> numpy.ndarray

qq_points(fit, data) → dict[str, ndarray]

Quantile-quantile point set for a fitted tail.

For a GPDFit, empirical quantiles are the sorted losses above the threshold and theoretical quantiles are the fitted conditional GPD quantiles at the plotting positions (i - 0.5) / n (both in original loss units, threshold included). For a GEVFit, the block maxima against fitted GEV quantiles. A good fit puts the points on the 45-degree line – deviations in the upper corner are exactly where a tail model earns or loses its keep.

Returns:

theoretical, empirical (equal-length arrays), n.

Return type:

dict

pp_points(fit, data) → dict[str, ndarray]

Probability-probability point set: model cdf vs plotting positions.

Complements qq_points(): PP is most sensitive in the body of the fitted range, QQ in the tail. Both on the unit square.

parametric_bootstrap_gof(fit, data, n_boot: int = 500, rng=None) → dict

Goodness-of-fit test for a fitted GPD or GEV, done honestly.

Kolmogorov-Smirnov and Anderson-Darling statistics of the data against the fitted model, with p-values from a parametric bootstrap that refits inside every replicate (see the module docstring for why the refit is not optional). A-D weights the tails, K-S the body; report both, because a tail model can pass one and fail the other.

Returns:

ks, ks_pvalue, ad, ad_pvalue, n, n_boot.

Return type:

dict

gpd_return_level(fit, return_periods, observations_per_period: float = 1.0, confidence_level: float = 0.95)

Return levels with confidence intervals from a POT/GPD fit.

The T-period return level is the loss exceeded once per T periods on average: with exceedance rate \(\zeta_u\) and m observations per period, it solves \(P(X > r) = 1/(T\,m)\):

\[r_T = u + \frac{\beta}{\xi}\left[(T\,m\,\zeta_u)^{\xi} - 1\right]\]

(\(u + \beta\log(T m \zeta_u)\) as \(\xi \to 0\)). Confidence intervals are by the delta method over \((\zeta_u, \xi, \beta)\), including the binomial variance of the exceedance rate (Coles, 2001, §4.3.3); they require the fit to carry a parameter covariance (populated by fit_gpd() / fit_pot when the information matrix is positive definite).

Parameters:

• fit (GPDFit) – A fit whose exceedance_fraction reflects the full dataset (i.e. from fit_pot, not raw fit_gpd on excesses alone).

• return_periods (float or array-like) – Periods T in the same period unit as observations_per_period.

• observations_per_period (float) – Observations per period (e.g. claims per year), so T * observations_per_period * exceedance_fraction is the expected number of threshold exceedances in T periods – it must exceed 1 for the return level to sit above the threshold.

• confidence_level (float) – Wald interval level.

Returns:

return_period, return_level, se, ci_low, ci_high (se/bounds are nan without a covariance).

Return type:

dict of str -> numpy.ndarray

gpd_tail_probability(x: float, threshold: float, xi: float, beta: float, exceedance_fraction: float) → float

Unconditional GPD tail probability \(P(X > x)\) above a POT threshold.

For \(x > u\), returns \(\zeta_u \, (1 + \xi (x - u)/\beta)^{-1/\xi}\) (exponential form as \(\xi \to 0\)), where \(\zeta_u\) is the exceedance rate; for \(x \le u\) it returns \(\zeta_u\) itself – the GPD says nothing below its threshold. Zero beyond the finite upper endpoint when \(\xi < 0\).

Parameters:

• x (float) – The loss level.

• threshold (float) – POT threshold \(u\) and the fitted GPD shape and scale.

• xi (float) – POT threshold \(u\) and the fitted GPD shape and scale.

• beta (float) – POT threshold \(u\) and the fitted GPD shape and scale.

• exceedance_fraction (float) – The exceedance rate \(\zeta_u = P(X > u)\), e.g. GPDFit.exceedance_fraction.

Returns:

Ground-up (unconditional) \(P(X > x)\).

Return type:

float

gpd_tvar(p: float, threshold: float, xi: float, beta: float, exceedance_fraction: float) → float

Unconditional GPD tail value-at-risk (expected shortfall) at level p.

Closed form on top of gpd_var():

\[\mathrm{TVaR}_p = \frac{\mathrm{VaR}_p + \beta - \xi u}{1 - \xi}, \qquad \xi < 1.\]

Infinite for \(\xi \ge 1\) (the tail has no mean) – a ValueError rather than a silent inf. Arguments as for gpd_var().

gpd_var(p: float, threshold: float, xi: float, beta: float, exceedance_fraction: float) → float

Unconditional GPD value-at-risk: the p-quantile of the ground-up loss.

Inverts the POT tail: with exceedance rate \(\zeta_u\),

\[\mathrm{VaR}_p = u + \frac{\beta}{\xi} \left[\left(\frac{1 - p}{\zeta_u}\right)^{-\xi} - 1\right]\]

(\(u + \beta \log(\zeta_u / (1 - p))\) as \(\xi \to 0\)). Valid only when the quantile lands in the fitted tail, i.e. \(1 - p < \zeta_u\); otherwise a ValueError – below the threshold the GPD has nothing to say.

Parameters:

• p (float) – Quantile level in (0, 1), e.g. 0.995.

• threshold (float) – POT threshold \(u\) and the fitted GPD shape and scale.

• xi (float) – POT threshold \(u\) and the fitted GPD shape and scale.

• beta (float) – POT threshold \(u\) and the fitted GPD shape and scale.

• exceedance_fraction (float) – The exceedance rate \(\zeta_u\), e.g. GPDFit.exceedance_fraction.

Returns:

The ground-up p-quantile.

Return type:

float

See also

gpd_tvar

The matching expected shortfall.

gpd_return_level

Return levels with delta-method intervals.

hill_curve(data, k_grid=None) → dict[str, ndarray]

Hill tail-index estimate across a grid of k, for a Hill plot.

Evaluates hill_estimator() at each k in k_grid. Plotting the returned hill values against k produces the standard Hill plot; the tail index is read from a region where the estimate is roughly flat, balancing the variance of small k against the bias of large k.

Parameters:

• data (array-like) – Strictly positive observations.

• k_grid (array-like of int, optional) – Values of k at which to evaluate the estimator. When None (the default), uses 1, 2, ..., max(2, len(data) // 4) so the grid stays within the upper quarter of the sample.

Returns:

k – the grid of order-statistic counts used, and hill – the corresponding Hill estimates, aligned elementwise.

Return type:

dict of str -> numpy.ndarray

Raises:

ValueError – If data contains a non-positive value, or if any k in the grid violates 1 <= k < len(data).

See also

hill_estimator

The per-k estimator evaluated here.

hill_estimator(data, k: int) → float

Hill estimator of the tail index from the k largest observations.

For positive data with an approximately Pareto upper tail, the Hill estimator of the tail index \(\gamma = 1/\alpha\) (equivalently the shape parameter \(\xi\) of the corresponding GPD) is the mean log spacing of the k largest order statistics above the (k+1)-th:

\[\hat{\gamma}_k = \frac{1}{k} \sum_{i=1}^{k} \bigl(\log X_{(n-i+1)} - \log X_{(n-k)}\bigr)\]

where \(X_{(1)} \le \dots \le X_{(n)}\) are the sorted values. The estimate is sensitive to the choice of k: small k gives high variance, large k introduces bias by reaching into the distribution body. Sweep k with hill_curve() and read the tail index from a stable region of the resulting plot.

Parameters:

• data (array-like) – Strictly positive observations. Order is irrelevant; the values are sorted internally.

• k (int) – Number of upper order statistics used, satisfying 1 <= k < len(data).

Returns:

The Hill tail-index estimate \(\hat{\gamma}_k\).

Return type:

float

Raises:

ValueError – If data contains a non-positive value, or if k is outside 1 <= k < len(data).

See also

hill_curve

Hill estimate across a grid of k.

pickands_estimator

Alternative tail-index estimator valid for any real tail index.

importance_sampling_diagnostics(weights) → dict[str, float]

Weight-degeneracy diagnostics for an importance sample.

Normalizes the weights and returns effective_n (Kish ESS), max_weight / min_weight, coefficient_of_variation, and entropy (natural log). A max weight near one, or an ESS far below n, says the proposal is missing the target’s mass and the estimates carry more uncertainty than their nominal n suggests. These diagnostics ride along on every estimate_*_is result.

layer_tail_metrics(result, *, q: float = 0.99, threshold: float | None = None) → dict[str, dict[str, float]]

Per-layer empirical VaR/TVaR from a risksim result’s layer losses.

The layer analogue of component_tail_metrics(): reads the (n_sims, n_layers) array result.layer_losses (names from layer_names when present) and returns, per layer, the empirical VaR and TVaR at level q, plus the exceedance probability of threshold when one is given.

Returns:

{name: {"var", "tvar"[, "exceedance_probability"]}}.

Return type:

dict of str -> dict of str -> float

log_importance_weights(log_target_density, log_proposal_density, *, normalize: bool = True) → ndarray

Importance weights from log-densities, computed stably in log space.

\(w_i \propto \exp(\log f(x_i) - \log g(x_i))\) for target \(f\) and proposal \(g\). With normalize=True (default) the weights are self-normalized via logsumexp – the numerically safe route when the densities span many orders of magnitude, exactly the regime importance sampling is used for. With normalize=False the raw (exponentiated) ratios are returned.

Parameters:

• log_target_density (array-like) – Log-densities of the two distributions at the sampled points; equal lengths.

• log_proposal_density (array-like) – Log-densities of the two distributions at the sampled points; equal lengths.

• normalize (bool, optional) – Self-normalize the weights to sum to one (default True).

Returns:

The importance weights.

Return type:

numpy.ndarray

losses_from_risksim(result, *, view: str = 'losses') → ndarray

Extract a 1-D loss array from a risksim-style simulation result.

Accepts any object exposing the requested attribute – view may be "losses", "gross" / "gross_losses", "retained" / "retained_losses", or "ceded" / "ceded_losses" – and returns it as a float array. Duck-typed: nothing risksim-specific is required beyond the attribute. Raises ValueError for an unknown view or one that is None on this result, TypeError if the attribute is absent.

make_blocks(data, block_size: int, *, drop_last: bool = True) → ndarray

Block maxima: the maximum of each consecutive block_size observations.

Partitions the data, in the given order, into consecutive blocks and returns each block’s maximum – the sample a GEV fit consumes. With drop_last=True (default) a trailing partial block is discarded, so every maximum comes from a full block; drop_last=False keeps the partial block’s maximum.

Parameters:

• data (array-like) – Observations, in block order (e.g. chronological).

• block_size (int) – Observations per block.

• drop_last (bool, optional) – Whether to drop a trailing partial block (default True).

Returns:

The block maxima.

Return type:

numpy.ndarray

Raises:

ValueError – If block_size exceeds the data length, or fewer than two blocks result (a GEV fit needs at least two maxima).

See also

fit_block_maxima

Blocking and fitting in one call.

mean_excess(data, thresholds) → dict[str, ndarray]

Empirical mean-excess function \(e(u) = E[X - u \mid X > u]\) on a grid.

For each candidate threshold, the average excess of the strictly-exceeding observations, with the exceedance counts alongside (nan where nothing exceeds). Linearity of \(e(u)\) in \(u\) is the classic peaks-over-threshold diagnostic: above any valid GPD threshold the mean excess is linear with slope \(\xi / (1 - \xi)\).

Parameters:

• data (array-like) – Loss observations.

• thresholds (array-like) – Candidate thresholds \(u\).

Returns:

thresholds, mean_excess, n_exceedances. With the plot extra, feed this to plot_mean_excess.

Return type:

dict of str -> numpy.ndarray

See also

threshold_diagnostic_table

The full scan, with GPD fits and error bands.

pickands_estimator(data, k: int) → float

Pickands estimator of the tail index from ordered tail spacings.

Unlike the Hill estimator, the Pickands estimator is valid for any real tail index \(\gamma\) (light, heavy, or bounded tails) and requires no positivity of the data beyond the internal ordering. Using the sorted values \(X_{(1)} \le \dots \le X_{(n)}\), it compares spacings at the k-th, 2k-th, and 4k-th largest observations:

\[\hat{\gamma}_k = \frac{1}{\log 2}\, \log\!\left( \frac{X_{(n-k+1)} - X_{(n-2k+1)}}{X_{(n-2k+1)} - X_{(n-4k+1)}} \right)\]

As with any tail-index estimator the result depends on k; it is typically read from a stable region of a plot over k.

Parameters:

• data (array-like) – Strictly positive observations. Order is irrelevant; the values are sorted internally.

• k (int) – Order-statistic spacing parameter, satisfying 4k < len(data) + 1 so that the 4k-th largest observation exists.

Returns:

The Pickands tail-index estimate \(\hat{\gamma}_k\).

Return type:

float

Raises:

ValueError – If data contains a non-positive value, if k violates 4k < len(data) + 1, or if the ordered tail spacings are not both strictly positive (which can occur with ties or a short tail).

See also

hill_estimator

Log-spacing estimator for heavy (positive-index) tails.

return_level(period: float, fit: GPDFit) → float

The loss exceeded on average once per period observations, under a POT fit.

Thin wrapper over GPDFit.return_level() – the fit’s unconditional quantile at 1 - 1/period – with the domain check that period exceeds one. Point estimates agree exactly with gpd_return_level(), which adds period units (observations_per_period) and delta-method confidence intervals.

return_period(probability: float) → float

The return period 1 / p of an event with exceedance probability p.

The reciprocal convention: an event with per-observation exceedance probability p recurs once per 1/p observations on average. Raises ValueError outside 0 < p < 1.

sample_lossmodel(model, size: int) → ndarray

Sample losses from a lossmodels-style severity or aggregate model.

stabilize_weights(weights, *, clip_quantile: float | None = None, renormalize: bool = True) → ndarray

Clip extreme importance weights at an upper quantile, then renormalize.

With clip_quantile=c, weights are capped at their own c-quantile – the standard variance-for-bias trade when a few weights dominate the sample; renormalize=True (default) rescales to sum to one afterwards. With clip_quantile=None this only (optionally) normalizes.

Parameters:

• weights (array-like) – Nonnegative importance weights (need not be normalized).

• clip_quantile (float, optional) – Cap level in (0, 1]; None disables clipping.

• renormalize (bool, optional) – Rescale to sum to one after clipping (default True).

Returns:

The stabilized weights.

Return type:

numpy.ndarray

See also

importance_sampling_diagnostics

Check whether stabilization is needed.

tail_summary_from_risksim(result, *, view: str = 'losses', thresholds=None, quantiles=(0.95, 0.99, 0.995)) → dict[str, object]

Tail summary of a risksim simulation result: extract the view, then summarize.

extreme_loss_summary(losses_from_risksim(result, view=view), ...) in one call: n / mean / std / min / max, empirical VaR/TVaR at each quantile, and the exceedance curve when thresholds is given.

threshold_diagnostic_table(data, thresholds) → ThresholdScan

Scan candidate POT thresholds: mean excess, GPD fits, and error bands.

For each candidate \(u\) with at least five exceedances, records the exceedance count, the empirical mean excess, the fitted GPD (xi, beta) (via fit_pot()), and the modified scale \(\beta^* = \beta - \xi u\) with delta-method standard errors from the fit covariance (gradient \((-u, 1)\) over \((\xi, \beta)\)). Above any valid threshold both \(\xi\) and \(\beta^*\) are constant in \(u\) – raw \(\beta\) drifts linearly even under a perfect GPD – so the selection question is where the two flatten within their bands. Candidates with fewer than five exceedances yield nan rows.

Parameters:

• data (array-like) – Ground-up loss observations.

• thresholds (array-like) – Candidate thresholds \(u\).

Returns:

Arrays over the grid: thresholds, mean_excess, xi, beta, n_exceedances, xi_se, modified_scale, modified_scale_se.

Return type:

ThresholdScan

Fitting from an Experience

extremeloss.integrations.actuarialpy fits GPD tails straight from a claims-listing Experience (or an ExperienceSet): fit_gpd_from_experience(exp, threshold=...) – extracting the excesses is structural; selecting the threshold stays the caller’s judgment. Aggregated experience tabs are refused. The core package stays array-level; install the [actuarialpy] extra or the openactuarial meta-package.