Multivariate Generalized Pareto Fitting

• Multivariate Generalized Pareto Fitting is a statistical technique for modeling joint tail behavior through limiting distributions and refined dependence structures.
• It employs componentwise threshold selection, marginal transformation, and innovative likelihood inference to capture extremal dependence in multiple dimensions.
• Recent advancements include threshold-free extensions and discrete adaptations that enhance simulation, stability, and scalability in high-dimensional settings.

A multivariate generalized Pareto fitting procedure seeks to model the joint tail behavior of a random vector via a family of limiting distributions—multivariate generalized Pareto distributions (MGPDs)—analogous to the univariate peaks-over-threshold (POT) paradigm but extended to capture extremal dependence structures in multiple dimensions. The theory is anchored in the convergence of the excess distribution over high thresholds, yielding a nondegenerate and canonical limiting law, whose structure is determined by marginal tail parameters and an extremal dependence mechanism (such as a D-norm, generator or spectral measure). The intricacies of parameterization, thresholding, inference, and diagnostics in the multivariate tail domain differentiate multivariate GP fitting from its univariate counterpart, and dictate a specialized suite of theoretical and computational tools.

1. Foundations: The Multivariate Peaks-Over-Threshold Limit and Generalized Pareto Laws

Let $X = (X_1, \ldots, X_d)$ be a random vector, $\mathbf{u} \in \mathbb{R}^d$ a high threshold vector, and consider the POT excess $Y = X - \mathbf{u}$, conditioned on $X \not\le \mathbf{u}$ (i.e., at least one component exceeds its threshold). Under the assumption that $F$ is in the domain of attraction of a max-stable law, the excess law

$$\mathcal L(Y \mid X \not\le \mathbf{u}) \to \operatorname{MGP}(\mathbf{u}, \boldsymbol{\sigma}, \boldsymbol{\xi}; S)$$

converges to a multivariate generalized Pareto (MGP) distribution, supported on the "L-shaped" region $\{\max_j y_j > 0\}$ (Naveau et al., 2024).

The canonical representation for MGPs involves:

• Marginal GPD parameters (scale $\sigma_j > 0$, shape $\xi_j \in \mathbb{R}$)
• An extremal dependence structure encoded by a "spectral" or "angular" component, typically described via an angular measure $H$ or a generator variable $\mathbf{u} \in \mathbb{R}^d$0.

A key integral representation for the multivariate GP distribution function is

$\mathbf{u} \in \mathbb{R}^d$1

with $\mathbf{u} \in \mathbb{R}^d$2 and $\mathbf{u} \in \mathbb{R}^d$3 the law of the generator (Naveau et al., 2024). Each margin, when conditioned appropriately, is univariate GPD.

This limiting law is intimately tied to the underlying extreme-value theory via point process arguments: the empirical process of threshold exceedances converges (after normalization) to a Poisson process with an intensity (exponent) measure $\mathbf{u} \in \mathbb{R}^d$4 that can be decomposed into radial and angular components (Naveau et al., 2024).

2. Parametric Constructions: Copula, Spectral, and Generator Approaches

Several representations of MGPDs or their copulas have become standard:

• Generalized Pareto Copulas (GPC): The upper-tail behavior of any copula $\mathbf{u} \in \mathbb{R}^d$5 in the domain of attraction of a multivariate EVD is approximated via

$\mathbf{u} \in \mathbb{R}^d$6

for a D-norm $\mathbf{u} \in \mathbb{R}^d$7 defined as $\mathbf{u} \in \mathbb{R}^d$8, $\mathbf{u} \in \mathbb{R}^d$9, $Y = X - \mathbf{u}$0 (Falk et al., 2018). This norm parametrizes multivariate extremal dependence structures.

• Spectral (Angular) Representations: Using a spectral measure $Y = X - \mathbf{u}$1 on the unit simplex, facilitating the decomposition of the exponent measure into radial and angular parts, and yielding the Pickands dependence function and stable tail dependence function.
• Generator-Based Parameterizations: The GP law may be constructed via various "generator" approaches, including "T-generator" (Poisson intensity), "U-generator" (tilted), and "R-generator" (real-scale), leading to specific density formulas and simulation schemes (Rootzén et al., 2016, Kiriliouk et al., 2016).

Table: Principal Parametric Constructions

Construction	Main Parameter(s)	Extreme-Dependence Modeling
D-norm/GPC	Generator $Y = X - \mathbf{u}$2, D-norm $Y = X - \mathbf{u}$3	GPC, via analytic representation
Angular measure $Y = X - \mathbf{u}$4	Spectral measure on simplex $Y = X - \mathbf{u}$5	Spectral/mixture models
Generator approaches	$Y = X - \mathbf{u}$6, $Y = X - \mathbf{u}$7, $Y = X - \mathbf{u}$8 random vectors	Flexible, explicit density

The D-norm framework enables practical implementation of fitting and simulation, while spectral and generator approaches are advantageous for flexible modeling of specific dependence regimes (logistic, Hüsler–Reiss, Dirichlet, etc.) (Falk et al., 2018, Rootzén et al., 2016, Kiriliouk et al., 2016).

3. Fitting Algorithms: Threshold Choice, Marginal Transformation, and Likelihood Inference

Multivariate GP fitting comprises the following methodological structure:

1. Componentwise Threshold Selection: Choose thresholds $Y = X - \mathbf{u}$9 per margin (via mean residual life, stability, or quantile plots). Empirical diagnostics (e.g., $X \not\le \mathbf{u}$0-plots of exceedance frequency at high quantile levels) inform the multivariate threshold (Kiriliouk et al., 2016, Rootzén et al., 2016).
2. Marginal Fitting: Univariate GPD parameters $X \not\le \mathbf{u}$1, $X \not\le \mathbf{u}$2 are estimated on the threshold excesses, usually by maximum likelihood (or PWM).
3. Standardization: Transform observed excesses via fitted marginal cdfs to pseudo-uniform or exponential margins, ensuring comparability across components.
4. Dependence Modeling: Above the chosen threshold, parametric dependence structures (e.g., D-norm, spectral, generator-based copula) are fitted to the multivariate excesses.
5. Likelihood Inference:
   • Full Likelihood: Joint loglikelihood over all dimensions, utilizing exact or approximate MGP density (Naveau et al., 2024, Rootzén et al., 2016).
   • Censored Likelihood: For observations with some components below threshold, integrate out the censored coordinates. This ensures essentially unbiased estimation and avoids tail-bias (Kiriliouk et al., 2016, Rootzén et al., 2016).
   • Composite Likelihood: For high dimensions, employ pairwise or low-dimensional marginal likelihoods for computational tractability, preserving consistent inference (Naveau et al., 2024).
6. Goodness-of-Fit Diagnostics: Threshold stability for parameters; marginal QQ-plots; uniformity tests for scaled residuals; empirical tail-dependence function comparison (Falk et al., 2018, Kiriliouk et al., 2016).

The piecing-together (PT) approach joins a fitted GPD copula in the tail with an arbitrary copula in the body, ensuring continuity at the threshold cut-point. The PT-copula is given by

$X \not\le \mathbf{u}$3

where $X \not\le \mathbf{u}$4 arises from the original copula and $X \not\le \mathbf{u}$5 from the GPD copula, with gluing at threshold (Aulbach et al., 2011).

4. Recent Extensions: Threshold-Free and Discrete Modeling

Recent developments extend classical GP fitting methods:

• Threshold-Free Extended GP Models: Multivariate extended GPDs (eGPD) specify a joint distribution agreeing with GPD tails but valid for all intensities, using hierarchical constructions (e.g., radial and angular decomposition) and avoiding subjective threshold choice (Alotaibi et al., 7 Sep 2025, Gaetan et al., 2 Oct 2025). Dependence for both tails is modeled via simplex-valued random vectors controlling upper- and lower-tail dependence. Parameter inference leverages modern neural Bayes estimators and semiparametric MLE using Bernstein polynomials and penalized regression on the radial/ angular components.
• Discrete Multivariate GPDs: For count data, as in drought analysis, the MDGPD employs a generator representation adapted from the continuous domain, with likelihood-free Bayesian estimation using permutation-invariant (DeepSets) neural architectures (Aka et al., 24 Jun 2025). Bootstrap generators provide end-to-end simulation and inference for multivariate discrete thresholds.
• Mixture and Model Selection Paradigms: Mixture-based formulations, such as multi-tree mixtures of Hüsler–Reiss components for spatial extremes, provide non-stationary and high-dimensional models with interpretable graphical structures (Cisneros et al., 2023).

5. Simulation, Stability Properties, and Asymptotics

Simulation from the fitted MGPD or GPC is essential for risk analysis and uncertainty quantification. Standard schemes include:

• Generator-Based Simulation: Simulate $X \not\le \mathbf{u}$6 from the D-norm generator, $X \not\le \mathbf{u}$7, set $X \not\le \mathbf{u}$8, obtain Pareto marginals, then invert marginals by the fitted univariate cdfs (Falk et al., 2018).
• Spectral- and T-generator Based Simulation: Draw a random vector $X \not\le \mathbf{u}$9 on the corresponding space, $F$0, set $F$1, then $F$2, and transform to the desired scale (Rootzén et al., 2016, Mourahib et al., 2023).

MGPDs possess sum-stability properties: for equal shape parameter, any nonnegative weighted sum is again univariate GP. Asymptotic properties (e.g., consistency and normality for estimator sequences) follow Z-estimator theory under the regular variation regime (Kiriliouk et al., 2016, Mourahib et al., 2023).

6. Practical Applications, Diagnostics, and Open Issues

Multivariate GP fitting underpins joint risk estimation for regulatory compliance, hydrological extremes, financial portfolios, and environmental phenomena. Empirical studies use jointly estimated marginal and dependence parameters to infer joint exceedance probabilities, construct confidence intervals, and test stability or independence hypotheses (Falk et al., 2018, Aulbach et al., 2012, Cisneros et al., 2023).

Key challenges and topics include:

• Threshold Selection vs. Global Tail Fitting: Choice of thresholds affects inference bias and variance; threshold-free approaches address this at the cost of more flexible (and complex) parametricism (Alotaibi et al., 7 Sep 2025).
• High-dimensional Structures: Scalability in $F$3 exploits sparsity (e.g., vine copulas, tree-based mixtures), composite likelihood, and regularization (Cisneros et al., 2023, Naveau et al., 2024).
• Identifiability: Multiple generator parameterizations may yield the same conditional excess distribution; caution is needed in interpretation (Naveau et al., 2024).
• Diagnostics: Empirical angular and radial tail dependence, Clarke’s $F$4-plots, sum-stability, and uniformity tests of residuals are standard tools for model validation (Falk et al., 2018, Kiriliouk et al., 2016).

7. Model Selection, Uncertainty Quantification, and Outlook

Model selection proceeds by comparing alternative generator classes (e.g., logistic, Hüsler–Reiss, Dirichlet), using information criteria (AIC), likelihood-ratio tests for nested models, and by performing goodness-of-fit diagnostics on marginal and multivariate features (Kiriliouk et al., 2016, Gaetan et al., 2 Oct 2025). Bayesian inference strategies, including posterior sampling and credible interval construction, accommodate threshold uncertainty and marginal dependence.

Composite likelihood and likelihood-free (simulation-based) methods are increasingly pivotal for high-dimensional and discrete/threshold-free extensions (Alotaibi et al., 7 Sep 2025, Aka et al., 24 Jun 2025, Hu et al., 2024). The multivariate peaks-over-threshold framework, anchored in the point-process perspective and generator-based construction, remains foundational, but advances in neural inference, mixture modeling, and semiparametric methods are broadening both theoretical scope and applied reach.