Generalized Linear Models for Extremes: Estimation and Inference in High Dimensions

Abstract: We propose a regression model for the extreme tail of a response variable, in which covariates rescale the tail without changing its shape. A single covariate-dependent function then characterizes the entire conditional tail, in contrast to extreme quantile regression, which targets a quantile at a pre-specified level. The tail shape itself is left unrestricted: heavy-, light- and short-tailed responses are covered by the same framework. We specify the function through a link function and a linear combination of the covariates, which is in the spirit of a generalized linear model. In estimation, we match the parametric specification to the underlying tail function under a Bregman divergence, over a region localized at the largest observations. The resulting loss is convex, and an $\ell_1$-penalty allows the number of covariates to exceed the effective sample size. The tail localization makes the asymptotic theory deviate from that for classical penalized generalized linear models. Only the tail observations selected by a random threshold are used in the statistical analysis, making them dependent. We derive the convergence rate of the penalized estimator and propose a debiased estimator that is asymptotically normal, yielding confidence intervals for individual coefficients. Its asymptotic variance is determined by the covariance of the score, which under tail localization differs from the Hessian and must be estimated separately. We apply the method to automobile insurance claims data.

• The paper develops a Bregman-divergence loss and ℓ1-penalized estimator for modeling covariate effects on extreme tails, achieving the rate √(s log p/k), where k is the effective number of tail observations.
• The paper introduces a debiasing method that accounts for the failure of the information identity under random tail thresholds, delivering asymptotically valid coordinate-wise confidence intervals in high dimensions.
• Simulations show 94.6–96.7% coverage for nominal 95% intervals, while an insurance application improves exceedance prediction and identifies significant factors associated with large claims.

Model and motivation

The paper studies regression for the extreme tail of a response variable through the heteroscedastic extremes framework. The authors assume the existence of a non-negative continuous function $c(x)$ such that

$$\lim_{y\to y^+}\frac{1-F_Y(y|X=x)}{1-F_Y(y)} = c(x),$$

where $y^+$ is the upper endpoint of the marginal support of $Y$. Under this assumption, covariates rescale the tail without changing its shape: the conditional and marginal distributions share the same upper endpoint and the same extreme value index, which may be positive (heavy-tailed), zero (light-tailed), or negative (short-tailed). This restriction is adopted deliberately; the authors note that allowing a covariate-dependent extreme value index complicates inference and cite recent work cautioning against it.

The scedasis function $c$ is specified parametrically as $c(x) = g(x^\top\beta)$, where $g$ is a pre-specified positive link function, in analogy with generalized linear models (GLMs). Taking expectations yields the identifying normalization $E\{c_\beta(X)\}=1$. Compared with extreme quantile regression, which targets a single conditional quantile at a pre-specified level, this specification characterizes the entire conditional tail through one function. Prior work on the scedasis function—nonparametric estimation with scalar indices, trend detection, space-time extensions—handles only one-dimensional covariates and offers no coefficient-level inference. The present paper claims to be the first treatment of the scedasis function with high-dimensional covariates and the first to deliver coefficient-level inference in that model.

Penalized estimation via Bregman divergence

Estimation proceeds by matching $c_\beta$ to the true conditional tail ratio under a Bregman divergence generated by $\phi$, restricted to generators satisfying $$\lim_{y\to y^+}\frac{1-F_Y(y|X=x)}{1-F_Y(y)} = c(x),$$0. This choice produces two simplifications: $$\lim_{y\to y^+}\frac{1-F_Y(y|X=x)}{1-F_Y(y)} = c(x),$$1, and $$\lim_{y\to y^+}\frac{1-F_Y(y|X=x)}{1-F_Y(y)} = c(x),$$2 where $$\lim_{y\to y^+}\frac{1-F_Y(y|X=x)}{1-F_Y(y)} = c(x),$$3. The resulting population loss is

$$\lim_{y\to y^+}\frac{1-F_Y(y|X=x)}{1-F_Y(y)} = c(x),$$4

and any minimizer automatically satisfies the normalization constraint. With i.i.d. data $$\lim_{y\to y^+}\frac{1-F_Y(y|X=x)}{1-F_Y(y)} = c(x),$$5 and an intermediate sequence $$\lim_{y\to y^+}\frac{1-F_Y(y|X=x)}{1-F_Y(y)} = c(x),$$6, $$\lim_{y\to y^+}\frac{1-F_Y(y|X=x)}{1-F_Y(y)} = c(x),$$7, the sample loss replaces the second term by the covariate average among observations exceeding the random threshold $$\lim_{y\to y^+}\frac{1-F_Y(y|X=x)}{1-F_Y(y)} = c(x),$$8:

$$\lim_{y\to y^+}\frac{1-F_Y(y|X=x)}{1-F_Y(y)} = c(x),$$9

The structure is notable: the response enters only through the $y^+$0 upper-order statistics, while all $y^+$1 covariate vectors contribute to the convex normalization term. An $y^+$2-penalized estimator $y^+$3 accommodates $y^+$4 larger than the effective sample size.

Two features distinguish the asymptotic theory from classical penalized GLMs. First, the effective sample size is $y^+$5, not $y^+$6: although $y^+$7 observations are available, tail information resides in the largest $y^+$8. Second, the selected tail observations are dependent even under i.i.d. sampling, because they are chosen by a common random order-statistic threshold, which blocks direct application of standard GLM arguments.

Convergence rates in fixed and high dimensions

In fixed dimension, under a second-order condition on the tail approximation error (Assumption M, controlling the bias rate $y^+$9), smoothness of $Y$0 at zero, and moment conditions on the link, the estimator satisfies an argmin theorem with rate $Y$1, converging to a Gaussian limit perturbed by the penalty term $Y$2. With vanishing penalty ($Y$3), the limit is normal with sandwich covariance $Y$4, where $Y$5 and $Y$6.

In high dimensions, under bounded sub-Gaussian designs and mild restrictions on the link (satisfied by $Y$7 and softplus), with $Y$8 and the scaling condition $Y$9, the estimator achieves

$c$0

This mirrors the classical $c$1 rate of penalized GLMs with $c$2 replacing $c$3, confirming that the effective sample size governs estimation difficulty. The condition permits $c$4 to exceed $c$5; the proof verifies restricted strong convexity from the full-sample curvature term while bounding the gradient's infinity norm using Bernstein inequalities adapted to the thresholded subsample.

Debiasing when the information identity fails

The regularized estimator carries non-negligible bias relative to variance, so the authors construct a debiased estimator via sample splitting: $c$6 is computed on one split, and the score, projection direction, and correction on an independent split. For each coordinate $c$7,

$c$8

where $c$9 minimizes the estimated score covariance $c(x) = g(x^\top\beta)$0 subject to a Hessian-based constraint $c(x) = g(x^\top\beta)$1 and a technical safeguard $c(x) = g(x^\top\beta)$2.

The central methodological point is that the information matrix identity breaks down under tail localization: because the score contains the indicator $c(x) = g(x^\top\beta)$3 while the Hessian does not, the covariance of the score and the Hessian are no longer asymptotically equivalent. The debiasing procedure therefore minimizes the true score covariance in its objective while retaining the Hessian in its constraints—a structural departure from standard GLM debiasing. A remark quantifies when this matters: under the canonical exponential link ($c(x) = g(x^\top\beta)$4), the sandwich collapses to $c(x) = g(x^\top\beta)$5, so ignoring the failure would still yield valid intervals for slope coefficients but overstates intercept variance by exactly one; under any non-canonical link, the two studentizers differ at first order.

Under additional smoothness conditions on the conditional density near the tail boundary and sparsity requirements $c(x) = g(x^\top\beta)$6, the debiased estimator is asymptotically normal:

$c(x) = g(x^\top\beta)$7

with $c(x) = g(x^\top\beta)$8, yielding coordinate-wise confidence intervals. The proof handles the dependence induced by the random threshold through a conditional empirical process argument, combining tightness of a tail-indexed process, Vervaat's lemma for the normalized order statistic, and a non-degeneracy argument showing $c(x) = g(x^\top\beta)$9 stays bounded away from zero.

Simulation evidence

Simulations use $g$0 with autoregressive Gaussian covariates, five active coefficients, $g$1, $g$2, $g$3, and three unit-tail-index error distributions (Pareto, Fréchet, absolute Student-$g$4); a supplementary design with finite-endpoint uniform responses confirms the framework covers short tails. Key findings:

• Selection–estimation trade-off: with $g$5, variable selection accuracy exceeds 95% in most configurations but inflates shrinkage bias on active coefficients (e.g., bias of $g$6 to $g$7 for $g$8 under the exponential link); with $g$9, accuracy drops to roughly 80% but bias roughly halves.
• Debiased coverage: nominal 95% intervals achieve coverage between 0.946 and 0.967 across configurations, with QQ plots closely tracking the diagonal.
• Cost of the naive Hessian studentizer: replacing the score covariance with the Hessian under the softplus link produces coverage between 0.850 and 0.876—substantial undercoverage—confirming empirically that the information identity failure has first-order consequences for non-canonical links.

Bias is consistently smaller under Pareto errors than Fréchet or Student-$E\{c_\beta(X)\}=1$0 errors, since the Pareto case involves no model approximation error; the other distributions add approximation bias on top of shrinkage bias.

Application to automobile insurance claims

The method is applied to a Kaggle auto claims dataset ($E\{c_\beta(X)\}=1$1, $E\{c_\beta(X)\}=1$2 after removing a collinear variable), with $E\{c_\beta(X)\}=1$3, so $E\{c_\beta(X)\}=1$4 is comparable to the effective sample size—a genuinely high-dimensional regime. On prediction of exceedance above the holdout 98% quantile, both proposed links outperform a no-covariate baseline (PE 9.87) and cross-validated penalized logistic regression (PE 9.80): PE 9.35 under the exponential link and 9.38 under softplus, averaged over 500 random splits.

Debiased inference across 50 splits, aggregated with median estimates and split-adjusted variances, identifies the same significant factors under both links: log vehicle value and its square increase tail intensity (exponential-link estimates 0.180 and 0.161), while city residence, managerial occupation, and minivan ownership decrease it (estimates $E\{c_\beta(X)\}=1$5, $E\{c_\beta(X)\}=1$6, $E\{c_\beta(X)\}=1$7). All adjusted confidence intervals exclude zero. These results separate policyholders who generate the largest claims from those who do not, information relevant to pricing and large-loss exposure management.

Limitations and open questions

The paper concedes several boundaries explicitly. There is no formal test of the core model assumption that covariates rescale rather than reshape the tail; the authors suggest a specification test based on constancy of conditional tail indices across the covariate space, analogous to existing tests, but leave its construction open. Second, the framework requires a pre-specified link function, and the theory relies on a second-order bias condition whose rate must be fast enough relative to $E\{c_\beta(X)\}=1$8; misspecification of either is not addressed. Third, inference on functionals beyond individual coefficients—conditional exceedance probabilities or extrapolated extreme quantiles via Weissman-type methods combined with $E\{c_\beta(X)\}=1$9—would require a delta-method argument involving the full score covariance, not merely single-coordinate variances; the authors note this extension is not routine and defer it. Finally, the sample-splitting construction, while simplifying the analysis, sacrifices efficiency relative to full-sample procedures, and the choice of the intermediate sequence $c_\beta$0 is not formally optimized.

Conclusion

The paper embeds the heteroscedastic extremes model into high-dimensional GLM methodology: a Bregman-divergence loss localized at the upper order statistics, an $c_\beta$1-penalized estimator with rate $c_\beta$2, and a debiasing scheme that correctly accounts for the breakdown of the information identity under tail localization. The theoretical distinction between score covariance and Hessian is shown to be consequential in practice, and the insurance application demonstrates that the procedure yields interpretable, statistically validated risk factors for tail behavior when the number of covariates rivals the effective tail sample size.