2.1. Model

Let Z be a positive, real-valued and heavy-tailed random variable. We assume that its conditional distribution given $\mathbf{X}$ satisfies

\begin{equation*}P\!\left( \left. Z \gt z\right\vert \mathbf{X}=\mathbf{x}\right) = z^{-\alpha\left( \mathbf{x}\right) }L\!\left( z;\;\mathbf{x}\right) ,\quad z \gt 0,\mathbf{x}\in \mathcal{X},\end{equation*}

where

(2.1) \begin{equation}\alpha \!\left( \mathbf{x}\right) = \sum_{i=1}^{I}\alpha _{i}\mathbb{I}_{\{\mathbf{x}\in \mathcal{A}_{i}\}} \gt 0,\quad \mathbf{x}\in \mathcal{X},\end{equation}

is the (unknown) tail index function that characterizes the dependence of the tail behavior of Z on $\mathbf{X}$ , and $L\!\left( z;\;\mathbf{x}\right) $ are slowly varying functions in the sense that $\lim_{z\rightarrow \infty}L\!\left( tz;\;\mathbf{x}\right) /L\!\left( z;\;\mathbf{x}\right) = 1$ for any $t \gt 0$ and $\mathbf{x}\in \mathcal{X}$ (see e.g., p. 37 and Appendix A3.1 in Embrechts et al., Reference Embrechts, Kluppelberg and Mikosch1997 for properties of slowly varying functions). Wang and Tsai (Reference Wang and Tsai2009) considered the same type of model but assumed that the tail index is related to the covariates through a log-linear link function (see also Li et al., Reference Li, Leng and You2020).

The tail index function takes a small number of values $\left( \alpha_{i}\right) _{i\in I}$ over a partition of the covariate space. Such an assumption may appear as practically irrelevant. It should first be kept in mind that a small number of values does not mean that the differences between these values cannot be large. This model is able to take into account strong heterogeneity. Second, since the tail index function is very important for risk management but at the same time difficult to estimate, it is interesting for the insurers to have a clear view on this function and to assume that the number of possible values is limited such they may have a good understanding of the combinations of covariates generating the smallest and the largest tail indexes. Equation (2.1) provides a parsimonious model of the tail index for an insurer wishing to use a model with low model risk, but high goodness-of-fit thanks to a suitable partition. The estimation of such a model is, however, more complex than a classical model with a fine grid partition of the covariate space for which a limited number of values would be included as a penalty, or for which the heterogeneity of estimated values would be limited (like the ridge regression-like penalty). Indeed, we impose here a spatial proximity of the areas for which the values are identical. We now present the estimation approach.

We seek to estimate the Hill index function $\xi \!\left( \mathbf{x}\right)=\alpha ^{-1}\!\left( \mathbf{x}\right) $ from a set of independent random variables $\mathcal{D}_{n}=\{\left( Z_{i},\mathbf{X}_{i}\right) _{i=1,\ldots,n}\}$ distributed as the independent pair $\left( Z,\mathbf{X}\right) $ . To do this, we introduce a family of positive threshold functions $t_{u}\!\left(\mathbf{\cdot }\right) = ut\!\left( \mathbf{\cdot }\right) $ with $u \gt 0$ and where $t\;:\;\mathcal{X}\rightarrow \mathbb{R}_{\mathbb{+}}$ satisfies $\inf_{\mathbf{x}\in \mathcal{X}}t\!\left( \mathbf{x}\right) \gt0$ . We only keep the observations $\left( Z_{i},\mathbf{X}_{i}\right) $ for which $Z_{i} \gt t_{u}\!\left( \mathbf{X}_{i}\right) $ for a large u. Let us define $Y^{(u)}=\ln \!\left( Z/t_{u}\!\left( \mathbf{X}\right) \right) $ given $Z \gt t_{u}\!\left( \mathbf{X}\right) $ and note that the distribution of $Y^{(u)}$ given $\mathbf{X}=\mathbf{x}$ may be approximated by an Exponential distribution with mean $\xi \!\left(\mathbf{x}\right) $ since

\begin{equation*}\lim_{u\rightarrow \infty }P(Y^{(u)} \gt y|\mathbf{X}=\mathbf{x)}=e^{-\alpha\left( \mathbf{x}\right) y},\quad y \gt 0\text{.}\end{equation*}

We will therefore work with the set of observations $\mathcal{D}_{n}^{(u)}=\{(Z_{i},\mathbf{X}_{i})\in \mathcal{D}_{n}\;:\;Z_{i} \gt t_{u}\!\left( \mathbf{X}_{i}\right) \}$ in order to build an estimate $\xi_{n}^{(u)}\;:\;[0,1]^{p}\rightarrow \mathbb{R}_{+}$ of the Hill index function $\xi $ , and we will use an appropriate loss function adapted to Exponential distributions. We denote by $n^{(u)}$ the cardinal number of $\mathcal{D}_{n}^{(u)}$ .

For this purpose, we consider the Gamma deviance. A deviance D is a bivariate function that satisfies the following conditions: $D\!\left(y,y\right) = 0$ and $D\!\left( y,\xi \right) \gt0$ for $y\neq \xi $ . In statistics, the deviance is used to build goodness-of-fit statistics for a statistical model. It plays an important role in Exponential dispersion models and generalized linear models. Let $f\!\left( y;\;\xi \right) $ be the density function of an Exponential random variable with mean $\xi $ . The Gamma deviance function is defined by

\begin{equation*}D\!\left( y,\xi \right) = 2\!\left( \ln f\!\left( y;\;y\right) -\ln f\!\left( y;\;\xi\right) \right) = 2\left[ \frac{y-\xi }{\xi }-\ln \!\left( \frac{y}{\xi }\right) \right] ,\quad y,\xi \gt0.\end{equation*}

Note that $D\!\left( y,\xi \right) \sim \!\left( y-\xi \right) ^{2}/\xi^{2}$ as $y\rightarrow \xi $ , and therefore, the Gamma deviance and the $L^{2}$ distance are equivalent when y is close to $\xi $ . The Gamma deviance is not only more appropriate than the $L^{2}$ distance because the observations $Y_{i}^{(u)}$ are asymptotically distributed as Exponential random variables, but also because it prevents the estimates from taking negative values.

Thereafter, we will consider families of estimators satisfying

\begin{equation*}\xi _{n}^{\left( u\right) }\left( \mathbf{\cdot }\right) = \arg \min_{f\in\mathcal{F}_{n}^{(u)}}\sum_{i\in \mathcal{I}_{n}^{(u)}}D(Y_{i}^{(u)},f\!\left(\mathbf{X}_{i}\right)\! )\end{equation*}

where $\mathcal{I}_{n}^{(u)}=\{i\;:\;(Z_{i},\mathbf{X}_{i})\in \mathcal{D}_{n}^{(u)}\}$ and $\mathcal{F}_{n}^{(u)}=\mathcal{F}_{n}(\mathcal{D}_{n}^{(u)})$ are a class of functions $f\;:\;\mathcal{X}\rightarrow \mathbb{R}_{+} $ that may depend on the data $\mathcal{D}_{n}^{(u)}$ .

2.2. Tree-based tail index estimators

Tree-based algorithms, such as Classification and Regression Trees (CART) (Breiman et al., Reference Breiman, Friedman, Olshen and Stone1984), random forests (Breiman, Reference Breiman2001) and boosted regression trees (Friedman, Reference Friedman2001), are popularly used in all kinds of data science problems because they are considered to be among the best supervised learning methods. They constitute a class of predictive models with high accuracy and capability of interpretation (see e.g., Chapters 9 and 10 in Hastie et al., Reference Hastie, Tibshirani and Friedman2009 for a general introduction).

The CART algorithm is the cornerstone of this class of algorithms. It makes a tree by splitting the sample into two or more homogeneous subsets based on most significant splitter/differentiator in explanatory variables. Both random forests and boosted regression trees create tree ensembles by using randomization during the tree creations. However, a random forest builds the trees in parallel and average them on the prediction, whereas a boosted regression tree creates a series of trees, and the prediction receives incremental improvement by each tree in the series.

2.2.1. Tail regression trees

A decision tree is derived from a rule-based method that generates a binary tree through a recursive partitioning algorithm that splits subsets (called nodes) of the data set into two subsets (called subnodes) according to the minimization of a split/heterogeneity criterion computed on the resulting sub-nodes. The root of the tree is $\mathcal{X}$ itself. Each split involves a single variable. While some variables may be used several times, others may not be used at all.

For tail index regression, the split criterion will be based on the Gamma deviance. To properly define it, we let A be a generic node and $n^{(u)}\!\left( A\right) $ be the number of data points of $\mathcal{I}_{n}^{(u)}$ such that $\mathbf{x}$ belongs to A. A cut in A is a pair (j, x), where j is an element of $\{1,...,p\}$ , and x is the position of the cut along coordinate j, within the limits of A. Let $\mathcal{C}_{A}$ be the set of all such possible cuts in A. Then, for any $(j,x)\in \mathcal{C}_{A}$ , the Gamma deviance split criterion takes the form

\begin{eqnarray*}L_{n}\!\left( j,x\right) &=&\frac{1}{n^{(u)}\!\left( A\right) }\sum_{i\in\mathcal{I}_{n}^{(u)}}D\!\left( Y_{i}^{(u)},\bar{Y}(A)\right) \mathbb{I}_{\{\mathbf{X}_{i}\in A\}}-\frac{1}{n^{(u)}\!\left( A\right) }\sum_{i\in \mathcal{I}_{n}^{(u)}}D\!\left( Y_{i}^{(u)},\bar{Y}(A_{L})\right) \mathbb{I}_{\{\mathbf{X}_{i}\in A,X_{i}^{(j)}\leq x\}} \\[5pt] &&-\frac{1}{n^{(u)}\!\left( A\right) }\sum_{i\in \mathcal{I}_{n}^{(u)}}D\!\left(Y_{i}^{(u)},\bar{Y}(A_{U})\right) \mathbb{I}_{\{\mathbf{X}_{i}\in A,X_{i}^{(j)} \gt x\}},\end{eqnarray*}

where $A_{L}=\{\mathbf{x}\in A\;:\;x^{(j)}\leq x\}$ , $A_{U}=\{\mathbf{x}\in A\;:\;x^{(j)} \gt x\}$ and $\bar{Y}(A)$ (resp., $\bar{Y}(A_{L})$ , $\bar{Y}(A_{U})$ ) is the average of the $Y_{i}^{(u)}$ ’s such that $\mathbf{X}_i$ belongs to A (resp., $A_{L}$ , $A_{U}$ ). Note that $L_{n}\!\left( j,x\right) $ simplifies in the following way

\begin{equation*}L_{n}\!\left( j,x\right) = \ln \!\left( \bar{Y}(A)\right) -\frac{n^{(u)}\!\left(A_{L}\right) }{n^{(u)}\!\left( A\right) }\ln \!\left( \bar{Y}(A_{L})\right) -\frac{n^{(u)}\!\left( A_{U}\right) }{n^{(u)}\!\left( A\right) }\ln \!\left( \bar{Y}(A_{U})\right) .\end{equation*}

At each node A, the best cut $\left( j_{n}^{\ast },x_{n}^{\ast }\right) $ is finally selected by maximizing $L_{n}\!\left( j,x\right) $ over $\mathcal{C}_{A}$ . The best cut is always performed along the best cut direction $j_{n}^{\ast }$ , at the middle of two consecutive data points.

Let us denote by $\mathcal{T}_{n}^{(u)}=\{A_{i,n}^{(u)}\;:\;i=1,\ldots,I_{n}^{(u)}\}$ the partition of $\mathcal{X}$ obtained where $I_{n}^{(u)}$ is the total number of leaf nodes in the tree. Then, the estimate of $\xi $ takes the form of a piecewise step function

\begin{equation*}\xi _{n}^{\left( u\right) }(\mathbf{x};\;\mathcal{T}_{n}^{(u)})=\sum_{i=1}^{I_{n}^{(u)}}\bar{Y}(A_{i,n}^{(u)})\mathbb{I}_{\left\{\mathbf{x}\in A_{i,n}^{(u)}\right\}},\end{equation*}

where $\mathbb{I}_{\{\mathbf{x}\in A_{i,n}^{(u)}\}}$ is the indicator function that $\mathbf{x}$ is in leaf node $A_{i,n}^{(u)}$ of the tree partition.

Decision tree output is very easy to understand and does not require any statistical knowledge to read and interpret it. Decision tree is one of the fastest way to identify most significant variables and relationships between two or more variables. One of the most practical issues is overfitting. A first way to solve it is to set constraints on model parameters: the users may for example fix the minimum number of observations which are required in a node to be considered for splitting, or the maximum depth of the tree (the number of edges from the root node to the leaf nodes of the tree), or else the maximum number of terminal nodes or leaves in the tree. A second way is to use pruning that consists in reducing the size of the decision tree by removing sections of the tree that are non-critical. By reducing the complexity, pruning improves predictive accuracy and mitigates overfitting.

2.2.2. Tail random forest

A random forest is a predictor consisting of a collection of several randomized regression trees which are built on different subsets of covariates and observations (see e.g., Scornet et al., Reference Scornet, Biau and Vert2015 for a formal presentation with the precise description of the algorithm). Let $\Theta $ be a random variable independent of $\mathcal{D}_{n}^{(u)}$ that will characterize the set of covariates among the components of $\mathbf{X}=(X^{(1)},...,X^{(p)})$ and the set of observations among $\mathcal{D}_{n}^{(u)}$ that will be used to build a tail regression tree. Let $\Theta _{1},\ldots ,\Theta _{M}$ be independent random variables, distributed as $\Theta $ and independent of $\mathcal{D}_{n}^{(u)}$ . For the j-th tail regression tree partition $\mathcal{T}_{n}^{(u)}\!\left( \Theta_{j}\right) $ obtained from the subset of covariates and observations characterized by $\Theta _{j}$ , we denote by $\xi _{n}^{\left( u\right) }(\mathbf{\cdot};\Theta _{j},\mathcal{D}_{n}^{(u)})$ the estimate of $\xi $ .

For our tail index regression problem, the trees will be combined through a harmonic mean to form the forest estimate (see Maillart and Robert, Reference Maillart and Robert2021)

(2.2) \begin{equation}\frac{1}{\xi _{M,n}^{\left( u\right) }(\mathbf{x};\;\Theta _{1},\ldots ,\Theta_{M})}=\frac{1}{M}\sum_{j=1}^{M}\frac{1}{\xi _{n}^{\left( u\right) }\left(\mathbf{x};\;\Theta _{j},\mathcal{D}_{n}^{(u)}\right)},\end{equation}

or equivalently

\begin{equation*}\alpha _{M,n}^{\left( u\right) }\left(\mathbf{x};\;\Theta _{1},\ldots ,\Theta _{M},\mathcal{D}_{n}^{(u)}\right)=\frac{1}{M}\sum_{j=1}^{M}\alpha _{n}^{\left( u\right)}\left(\mathbf{x};\;\Theta _{j},\mathcal{D}_{n}^{(u)}\right),\end{equation*}

where $\alpha _{M,n}^{\left( u\right) }$ and $\alpha _{n}^{\left( u\right) }$ denote the respective tail index estimates. Note that such an aggregation is different from the one done for the usual random forest and ensures that the Gamma deviance loss of $\xi _{M,n}^{\left( u\right) }$ will be smaller than the average of the individual Gamma deviance losses of the $\xi _{n}^{\left(u\right) }$ ’s (see Maillart and Robert, Reference Maillart and Robert2021). Let us denote by $\mathcal{R}_{n}^{(u)}=\{B_{j,n}^{(u)}\;:\;j=1,\ldots ,J_{n}^{(u)}\}$ the partition of rectangles of $\mathcal{X}$ obtained from crossing the partitions of the regression trees $\mathcal{T}_{n}^{(u)}\!\left( \Theta _{1}\right) ,\ldots ,\mathcal{T}_{n}^{(u)}\!\left( \Theta _{M}\right) $ . The tail random forest estimate of $\xi $ also takes the form of a piecewise step function

\begin{equation*}\xi _{M,n}^{\left( u\right) }(\mathbf{x};\;(\Theta_{m}))=\sum_{j=1}^{J_{n}^{(u)}}b_{j,n}\mathbb{I}_{\{\mathbf{x}\in B_{j,n}^{(u)}\}},\end{equation*}

where

\begin{equation*}\frac{1}{b_{j,n}}=\frac{1}{M}\sum_{j=1}^{M}\frac{1}{\xi _{n}^{\left(u\right) }(\mathbf{x};\;\Theta _{j},\mathcal{D}_{n}^{(u)})}\mathbb{I}_{\{\mathbf{x}\in B_{j,n}^{(u)}\}}\text{.}\end{equation*}

As for the usual random forests, the algorithm needs two additional parameters: the number of pre-selected covariates for splitting, the number of sampled data points in a tree. Tail random forests can be used to rank the importance of variables in a natural way as described in Breiman’s original paper (Breiman, Reference Breiman2001). Tail random forests achieve higher accuracy than a single tail regression tree and suffer less from the overfitting issue, but they sacrifice the intrinsic interpretability present in decision trees and may appear as a black-box approach for statistical modelers.

2.2.3. Tail gradient tree boosting

Tail gradient tree boosting combines weak tree “learners” (small depth tail regression trees) into a single strong learner in the following iterative way:

1. – Initialize the model with a small depth tail regression tree $\xi_{0,n}^{\left( u\right),g }$ . Let M be an integer.

2. – For $m=1,\ldots ,M$ , compute the pseudo-residuals

   \begin{equation*}r_{i,m}^{(u)}=-\left. \frac{\partial D\!\left( y_{i},\xi \right) }{\partial\xi }\right\vert _{\xi = \xi _{m-1,n}^{\left( u\right),g }(\mathbf{x}_{i})},\quad i\in \mathcal{I}_{n}^{(u)}\text{.}\end{equation*}

3. – At the m-th step, the algorithm fits a regression tree $h_{m,n}^{(u)}\!\left( \mathbf{x}\right) $ to pseudo-residuals $(r_{i,m}^{(u)})_{i\in \mathcal{I}_{n}^{(u)}}$ such that

   \begin{equation*}h_{m,n}^{(u)}\!\left( \mathbf{x}\right) = \sum_{k=1}^{K_{m,n}^{(u)}}c_{k,m}\mathbb{I}_{\{\mathbf{x}\in C_{k,m,n}^{(u)}\}},\end{equation*}

where $\mathcal{G}_{m,n}^{(u)}=\{C_{k,m,n}^{(u)}\;:\;k=1,\ldots ,K_{m,n}^{(u)}\}$ is the associated partition of the tree and $\left( c_{k,m}\right)_{k=1,\ldots ,K_{m,n}^{(u)}}$ are the predicted values in each region.

1. – The model update rule is then defined as

   \begin{equation*}\xi _{m,n}^{\left( u\right),g }(\mathbf{x})=\xi _{m-1,n}^{\left( u\right),g }(\mathbf{x})+\sum_{k=1}^{K_{m,n}^{(u)}}\gamma _{k,m}\mathbb{I}_{\{\mathbf{x}\in C_{k,m}\}}\end{equation*}
   where
   \begin{equation*}\gamma _{k,m}=\arg \min_{\gamma }\sum_{\mathbf{x}_{i}\in C_{k,m}}D\!\left(y_{i},\xi _{m-1,n}^{\left( u\right),g }(\mathbf{x}_{i})+\gamma \right) .\end{equation*}

One of the outputs of this algorithm is a partition of rectangles of the covariate space: $\mathcal{G}_{n}^{(u)}=\{C_{j,n}^{(u)}\;:\;j=1,\ldots,K_{n}^{(u)}\}$ such that the gradient tree boosting estimate of $\xi $ takes the form of a piecewise step function

\begin{equation*}\xi _{M,n}^{\left( u\right),g }(\mathbf{x})=\sum_{j=1}^{K_{n}^{(u)}}c_{j}\mathbb{I}_{\{\mathbf{x}\in C_{j,n}^{(u)}\}}\end{equation*}

where the superscript g is used to refer to the gradient approach.

The number of gradient boosting iterations M appears as a regularization parameter. Increasing M reduces the error on training set, but setting it too high may lead to overfitting. An optimal value of M is often selected by monitoring prediction error on a separate validation dataset. At each iteration of the algorithm, it is also possible to only consider a subsample of the training set (drawn at random without replacement) to fit the weak tree learner. Friedman (Reference Friedman2001) observed a substantial improvement accuracy with this modification in the case of the usual gradient boosting. Subsampling may introduce randomness into the algorithm and help prevent overfitting, acting also as a kind of regularization. Another useful regularization techniques for gradient-boosted trees is to penalize model complexity of the learned model. The model complexity is in general defined as the number of leaves in the learned trees.

2.3. Hierarchical clustering with spatial constraints

2.3.2. ClustGeo algorithm by-products

ClustGeo algorithm of Chavent et al. (Reference Chavent, Kuentz-Simonet, Labenne and Saracco2018) provides a partition of the covariate space: $\mathcal{P}_{n^{(u)}}=\{\mathcal{A}_{i,n^{(u)}}\;:\;i=1,\ldots ,I_{n^{(u)}}\}$ . For each subset $\mathcal{A}_{i,n^{(u)}}$ of the partition, it can be assumed that the distribution of $Y_{j}^{(u)}=\ln \!\left( Z_{j}/t_{u}\!\left( \mathbf{X}_{j}\right) \right) $ given $Z_{j}>t_{u}\!\left( \mathbf{X}_{j}\right) $ is approximately an Exponential distribution with parameter $\hat{\alpha}_{i,n^{(u)}}$ (obtained as the reciprocal of the average of the Hill index predictions over the subset $\mathcal{A}_{i,n^{(u)}}$ ). We can then derive several interesting by-products.

First, we can propose from the Central Limit Theorem a 95% confidence interval for $\hat{\alpha}_{i,n^{(u)}}$

\begin{equation*}\left( \hat{\alpha}_{i,n^{(u)}}-1.96\frac{\hat{\alpha}_{i,n^{(u)}}}{\sqrt{n_{i}^{(u)}}},\hat{\alpha}_{i,n^{(u)}}+1.96\frac{\hat{\alpha}_{i,n^{(u)}}}{\sqrt{n_{i}^{(u)}}}\right)\end{equation*}

where $n_{i}^{(u)}$ is the number of the $\mathbf{X}_{j}$ ’s belonging to $\mathcal{A}_{i,n^{(u)}}$ such that $Z_{j}>t_{u}\!\left( \mathbf{X}_{j}\right) $ . Indeed, $\hat{\alpha}_{i,n^{(u)}}$ is very close to the reciprocal of the averages of the $Y_{j}^{(u)}$ ’s.

Second, it is possible to build tail-quantile estimates. For $z>t_{u}\!\left( \mathbf{x}\right) $ , the probability of exceedance is given by

\begin{eqnarray*} P\!\left( Z\gt z|\mathbf{X}=\mathbf{x}\right) &=&P\!\left( Z\gt t_{u}\!\left( \mathbf{x}\right) |\mathbf{X}=\mathbf{x}\right)P(Z \gt z|\mathbf{X}=\mathbf{x},Z \gt t_{u}\!\left( \mathbf{x}\right)\!\mathbf{)} \\[5pt] &=&P\!\left( Z \gt t_{u}\!\left( \mathbf{x}\right) |\mathbf{X}=\mathbf{x}\right)P(\!\ln \!\left( Z/t_{u}\!\left( \mathbf{x}\right) \right) \gt\ln \!\left(z/t_{u}\!\left( \mathbf{x}\right) \right) |\mathbf{X}=\mathbf{x},Z \gt t_{u}\!\left(\mathbf{x}\right)\! \mathbf{)}.\end{eqnarray*}

Since, for large u,

\begin{equation*}P(\!\ln \!\left( Z/t_{u}\!\left( \mathbf{x}\right) \right) \gt y |\mathbf{X}=\mathbf{x},Z>t_{u}\!\left( \mathbf{x}\right) \mathbf{)}\simeq e^{-\alpha _{i}y},\quad y>0,\mathbf{x}\in \mathcal{A}_{i},\end{equation*}

we deduce that, for $\mathbf{x}\in \mathcal{A}_{i}$ , this probability is approximated for large u by

\begin{equation*}P\!\left( Z \gt t_{u}\!\left( \mathbf{x}\right) |\mathbf{X}=\mathbf{x}\right) \left(\frac{z}{t_{u}\!\left( \mathbf{x}\right) }\right) ^{-\alpha _{i}}.\end{equation*}

The probability $P\!\left( Z>t_{u}\!\left( \mathbf{x}\right) |\mathbf{X}=\mathbf{x}\right) $ may be estimated by an empirical probability based on observations for which the covariate $\mathbf{X}$ belongs to a neighborhood of $\mathbf{x}$ . Then, for a large probability level (close to one), an estimate of its associated quantile is then easily derived by inverting the previous formula with respect to z and replacing $P\!\left( Z>t_{u}\!\left( \mathbf{x}\right) |\mathbf{X}=\mathbf{x}\right) $ by its empirical counterpart and $\alpha _{i}$ by its estimate $\hat{\alpha}_{i,n^{(u)}}$ .

2.4. Equivalent tree: a tail tree surrogate

Our procedure described in the previous subsections provides a small-size partition of the covariate space and ultimately makes accurate predictions by modeling the complex structure of the partition. However, it is not able to explain easily how the partition subsets are made from the covariate vector $\mathbf{X}$ , and then, it can be viewed as a black-box model producing predictions. Therefore, we attach to our procedure a global tree surrogate model that approximates the partition-based rules while providing an explainable model using the initial covariates. This surrogate model combines a binary encoder representation of the additive tree ensemble with a regression tree. It is called the equivalent tree model.

Let us denote by $R_{g}$ a rectangle of the partition of the covariate space obtained from the additive tree ensemble. We are interested in a binary representation of $R_{g}$ . We assume that the additive tree ensemble has L internal nodes/splits of the form $X_{i_{l}}>z_{l}$ for $l=1,\ldots ,L$ . The rectangle $R_{g}$ may then be represented by the binary vector $\boldsymbol{\eta }_{g}\in \{0,1\}^{L}$ such that $\mathbf{x}\in R_{g}$ if $x_{i_{l}}>z_{l}$ for $l=1,\ldots ,L$ . Figure 1 illustrates this representation on a simple example where $p=2$ and $L=4$ .

Figure 1. A simple example of the binary representation of the rectangles of the covariate space obtained from an additive tree ensemble.

We now propose an alternative partition based on a regression tree and proceed as follows:

1. 1. we associate to each rectangle $R_{g}$ its predicted value (as a new label), its binary representation (as a new feature vector) and a weight corresponding to the proportion of the observations $(Y_{i}^{(u)},\mathbf{X}_{i})$ such that $\mathbf{X}_{i}$ belongs to it;

2. 2. we build a maximal regression tree from these new data.

By choosing an appropriate depth, we obtain an approximation of the large-size partition of the additive tree ensemble with a reasonable explanation. Note that the fidelity measured by the $R^{2}$ measure (i.e., the percentage of variance that our surrogate model is able to capture from the additive tree ensemble) increases at each split. Since the regression tree method is applied to binary predictors that are linked to the splits made by the tree-based model, it divides the covariate space from the same rectangles as the tree-based model. The regression tree method thus provides a nested set of trees that approximates faithfully the large-size partition.

It should be noted that additive tree ensemble models natively provide local model interpretations for analyzing predictions. Actually, the vector of covariates associated to a prediction belongs to a unique rectangle of the large-size partition whose edges are intervals belonging to the support of the covariate components. Such local explanations can be used to answer questions like: why did the model make this specific prediction? or what effect did this specific feature value have on the prediction?. But additive tree ensemble models need additional surrogate models for good global explanations.

3.1. Simulation studies

This section investigates how our approach performs on simulated data. We consider two toy models where $p=2$ , $\mathcal{X}=[0,1]\times \lbrack 0,1]$ , $X_{1}$ and $X_{2}$ are two independent random variables uniformly distributed over [0, 1],

\begin{equation*}P\!\left( \left. Z>z\right\vert \mathbf{X}=\mathbf{x}\right) = z^{-\alpha\left( \mathbf{x}\right) },\quad z>1,\mathbf{x}\in \mathcal{X},\end{equation*}

and the partition subsets are squares for Model 1 and parts of nested discs for Model 2. More specifically, the respective partitions $\mathcal{P}^{\left( 1\right) }$ and $\mathcal{P}^{\left( 2\right) }$ of $\mathcal{X}$ are of size 4 and are depicted in Figure 2.

Figure 2. Partitions $\mathcal{P}^{\left( 1\right) }$ and $\mathcal{P}^{\left( 2\right) }$ of $\mathcal{X}$ with their associated tail index values.

Geometric figures such as squares are a priori easily identified by gradient tree boosting algorithms. However, disks do not appear to be natural partitions for this type of algorithms. These two extreme cases will allow us to observe the efficiency of our approach for sets of partitions with potentially complex shapes.

Since the slowly varying functions of the family of conditional survival distribution functions of Z are equal to 1, it is not necessary to choose a family of thresholds for these datasets or to select observations whose values are higher than the thresholds. The sizes of the datasets are chosen to be equal to $n=50,000$ . The datasets $\mathcal{D}^{(1)}$ and $\mathcal{D}^{(2)}$ are given by the sets of values $\{\left( Z_{i},\mathbf{X}_{i}\right) _{i=1,\ldots ,n}\}$ for each model. In Figure 2, each point represents an observation (the color of a point is given by the value of its associated observation).

We choose the tail gradient tree boosting method to build the large-size partitions. The gradient tree boosting combines 100 weak tree learners. We performed a grid search with a 5-fold cross-validation to find a set of good hyper-parameters. We evaluated the performance of the models with the Gamma deviance. The values of the hyper-parameters are presented in Table 1 (see https://docs.h2o.ai/h2o/latest-stable/h2o-docs/data-science/gbm.html for the definitions of the hyper-parameters). We simulated test datasets with the same size to understand how the models generalize.

Table 1. Choice of the hyper-parameters.

Figure 3 provides the large-size partitions obtained for $\mathcal{D}^{(1)}$ and $\mathcal{D}^{(2)}$ . The color inside a rectangle represents the value of the prediction of the Hill index function for this rectangle. The tail gradient tree boosting model performs well on $\mathcal{D}^{(1)}$ which is not surprising because the subsets of the partition $\mathcal{P}^{\left( 1\right) }$ are characterized by intersections of conditions that only depend on $X_{1}$ or $X_{2}$ , but not both. The learning task is more difficult and challenging for $\mathcal{D}^{(2)}$ because the subsets of the partition $\mathcal{P}^{\left( 2\right) }$ depend on an intricate way of $X_{1}$ and $X_{2}$ .

Figure 3. Tail gradient tree boosting partitions for $\mathcal{D}^{(1)}$ and $\mathcal{D}^{(2)}$ .

We then run ClustGeo to gather the rectangles of the tail gradient tree boosting partitions in order to reveal the partitions $\mathcal{P}^{\left( 1\right) }$ and $\mathcal{P}^{\left( 2\right) }$ . The matrix $D_{1}$ which gives the dissimilarities in the spatial space is defined in the following way: if two rectangles are adjacent, the distance value is set to $0.1,$ and if this is not the case, the distance value is set to a large value $d=9.10^{6}$ . Figure 4 shows the fidelity curves (which are defined as the curves of the $R^{2}$ measures between the aggregated partition and the initial partition created by the tail gradient tree boosting for different values of the weight parameter $\gamma \in \lbrack 0,1]$ ). On the basis of these curves, we choose $\gamma = 0.3$ for $\mathcal{D}^{(1)}$ and $\gamma = 0.1$ for $\mathcal{D}^{(2)}$ , while fixing the size K of the small-size partitions to 4. The fidelity is then equal to $94.91$ % for $\mathcal{D}^{(1)}$ and $91.18$ % for $\mathcal{D}^{(2)}$ .

Figure 4. Fidelity curves for $\mathcal{D}^{(1)}$ and $ \mathcal{D}^{(2)}$ .

Figure 5. Estimated small-size partitions for $\mathcal{D}^{(1)}$ and $\mathcal{D}^{(2)}$ .

Figure 6. QQ-plots comparing the distributions of the normalized observations inside each subsets of the estimated partitions with the unit Exponential distribution.

Figure 7. Fidelity curves of the equivalent tree models for $\mathcal{D}^{(1)}$ and $\mathcal{D}^{(2)}$ .

Figure 8. Approximated small-size partitions of the equivalent tree models for $\mathcal{D}^{(1)}$ and $\mathcal{D}^{(2)}$ .

Figure 9. QQ-plots comparing the distributions of the normalized observations inside each subsets of the estimated partitions with the unit Exponential distribution.

Figure 5 provides the estimated small-size partitions for $\mathcal{D}^{(1)}$ and $\mathcal{D}^{(2)}$ . For $\mathcal{D}^{(1)}$ , we observe that the hierarchical clustering algorithm has some difficulties in separating the region where the Hill index is equal to $0.2$ from the one where it is equal to $0.25$ . This reason is that these values are actually too close. Moreover, it creates a fourth small subset for the partition in the region where the Hill index is equal to $0.25$ without providing a very different predicted value. A partition with three subsets would therefore be sufficient here. Nevertheless, the shapes of the estimated subsets of the partition are close to the shapes of the true subsets. For $\mathcal{D}^{(2)}$ , we observe that the circular subsets for the tail indexes equal to $0.25$ , $0.33$ and $0.5$ are relatively well estimated, but the subset for the value $0.25$ is the union of two disjoint subsets. We conclude that ClustGeo does a good job for identifying the small-size partitions.

In Figure 6, QQ-plots compare the distributions of the normalized observations inside each subsets of the estimated partitions (they have been divided by their mean) with the unit Exponential distribution. The alignment of the points in the 95% pointwise confidence intervals (based on order statistics of the unit Exponential distribution) shows good fits and validates the assumption of Exponential distributions for $\log \!\left( Z\right) $ . The tail indexes for the orange and yellow subsets are very well estimated.

Finally, we compare the estimated partitions with the outputs of the equivalent tree model. Figure 7 provides the fidelity curves ( $R^{2}$ measure) between the approximated partitions and the initial partition created by the tail gradient tree boosting. We decide to choose $K=4$ subsets for $\mathcal{D}^{(1)}$ with a high fidelity ( $98.10$ %) and $K=7$ subsets for $\mathcal{D}^{(2)}$ with a good fidelity ( $79.24$ %). Figure 8 shows the approximated small-size partitions for $\mathcal{D}^{(1)}$ and $\mathcal{D}^{(2)}$ . The equivalent tree model provides the exact partition for $\mathcal{D}^{(1)}$ , doing better than the hierarchical clustering algorithm. Although the approximated small-size partition differs from the true one for $\mathcal{D}^{(2)}$ , the explanations that could be made would be very convincing. In Figure 9, QQ-plots compare the distributions of the normalized observations inside each subsets of the approximated partitions with the unit Exponential distribution. The alignment of the points also shows good fits.

3.2. A case study for the insurance industry

We consider a NOAA tornado dataset containing 60,652 events from 1950 to 2011. The variables that we kept from this dataset are given in Table 2 (the other variables were quantitative variables with too many missing values or categorical variables). We combined the variables PROPDMGEXP and PROPDMG to a single variable PROPERTY DAMAGE containing the estimated costs (in dollars) of the damages caused by tornadoes. We only retained strictly positive amounts of PROPERTY DAMAGE and ended up with 39,036 observations. We finally extracted YEAR and MONTH from the variable DATE.

Table 2. Descriptive table of variables.

We are interested in characterizing the extremal behavior of property damages caused by tornadoes. The documentation provided with the dataset does not mention whether inflation has been taken into account to evaluate property damages, but after analyzing the amounts of damages, we concluded that it was not the case and decided to adjust these amounts for inflation (we used the historic American inflation based upon the consumer price index because it was the only index with a large enough history). We also decided to add contextual information about the population densities since the damages caused by tornadoes are correlated to the human constructions in the area they strike. We used a shapefile containing geographical frontiers of US counties that we linked with the census data to extract population densities by year and county (DENSITY) (we did not have information on the populations affected by the tornadoes, and therefore, we assumed that the affected population densities were proportional to the population densities). We thereafter consider as our variable of interest the following variable: DAMAGE BY DENSITY $=$ PROPERTY DAMAGE DENSITY.

Figure 10 illustrates the linear relationships between $\log \!\left( \text{DAMAGE BY DENSITY}\right) $ with respectively $\log \!\left( \text{LENGTH}\right) $ and $\log \!\left( \text{WIDTH}\right) $ , while Figure 11 provides a scatter plot showing that the longer and the wider a tornado track, the greater the damage will be.

Figure 10. Linear relationships between $\log \!\left( \text{DAMAGE BY DENSITY}\right) $ with respectively $\log \!\left( \text{LENGTH}\right) $ and $\log \!\left( \text{WIDTH}\right) $ .

Figure 11. Scatter plot of $\log \!\left( \text{LENGTH}\right) $ and $\log \!\left( \text{WIDTH}\right) $ .

Figure 12 displays a map of the United States with the population density per county as well as the locations of a random subsample of the tornadoes depicted with different colors depending on the amount of damages.

Figure 12. Population density in the US and tornado locations.

After a detailed study of the tails of the distributions of the variable DAMAGE BY DENSITY, we concluded that it was necessary to make a deterministic transformation of this variable so that the observations are compatible with the hypothesis of a Pareto-type distribution. The tails of the distributions are in fact of Weibull-type with a coefficient $\theta $ that can be estimated following the approach developed in Girard (Reference Girard2004). The values of the estimators of the Weibull tail-coefficient $\hat{\theta}_{n}$ based on the $k_{n}$ upper order statistics are given in the left panel of Figure 13. In practice, the choice of the parameter $k_{n}$ is the key problem to obtain a correct estimation of $\theta $ . If $k_{n}$ is too small, the variance of $\hat{\theta}_{n}$ may be high and conversely, if $k_{n}$ is too large, the bias may be important. We retained the approach proposed in Girard (Reference Girard2004) and selected the values $k_{n}=5000$ and $\hat{\theta}_{n}=4$ . Therefore, we transform the variable DAMAGE BY DENSITY in the following way

\begin{equation*}Z=\exp \!\left( \left( \text{DAMAGE BY DENSITY}\right)^{1/4}\right) .\end{equation*}

Figure 13 displays a QQ-plot comparing the distribution of $Y=\log \!\left(Z\right) $ with an Exponential distribution. It should be kept in mind that the threshold exceedance distribution of Y is expected to be a mixture of Exponential distributions with different parameters in different subsets of the covariate space. It is therefore not possible to obtain a perfect alignment of the points as in the case of a simple Exponential distribution. However, this graph allows us to think that the choice of $\hat{\theta}_{n}$ is consistent with our working hypothesis.

Figure 13. Left panel: Weibull tail-coefficient estimator $\hat{\theta}_{n}$ ; Right panel: QQ-plot comparing the distributions of the observations with the Exponential distribution.

The covariates that have been retained for building an additive tree model are MAG, LENGTH, WIDTH, LATITUDE, LONGITUDE, YEAR. We divided the observations into two subsets : a training set (80% of observations) and a validation set (20% of observations). As mentioned previously, we chose the threshold function for which the function t is uniformly constant per region where the regions were obtained from products of median class intervals of the covariates. In a first step, we fitted several tail gradient tree boosting models on all the data with not optimized hyper-parameters to choose the best threshold function. The selected threshold function was the one that retains 90% of the largest values in each region. Then, we fitted again tail gradient tree boosting models on the subset of data to optimize the values of the hyper-parameters. The choice of the hyper-parameters are given in Table 3.

The grid search for the hyper-parameters was performed with a 5-fold cross-validation. We selected the model with the smallest Gamma deviance, and we denote by GBM Gamma this model (the Gamma deviance on the train set is equal to $0.8389$ , on the validation set $0.8488$ and by cross-validation $0.8422$ ).

We then run ClustGeo to gather subsets of the tail gradient tree boosting partition. We first choose $D_{1}$ as the distance which is set to $0.1$ if two rectangles are adjacent and to a large value if this is not the case. To help us choose the mixing parameter $\gamma $ , we plot the proportions (resp. normalized proportions) of explained pseudo-inertia of the partitions in K clusters obtained with the ClustGeo procedure for a range of $\gamma $ (the pseudo-inertia of a cluster is calculated from the dissimilarity matrix and not from the data matrix). When the proportion (resp. normalized proportion) of explained pseudo-inertia based on $D_{0}$ decreases, the proportion (resp. normalized proportion) of explained pseudo-inertia based on $D_{1}$ increases. The plots are given in Figure 14.

We also compute the fidelity curves with respect to $\gamma$ (see Figure 15). On the basis of these plots, we choose $\gamma = 0.1$ and $K=12$ ( $R^{2}$ measure is equal to $88.38$ %.). Figure 16 displays box plots of the predicted values by GBM Gamma for each subset of the estimated partition, as well as box plots of their relative differences with their averages (i.e., by taking the averages as the reference values). For almost every subsets, more than 50% of the predictions have relative differences less than 15% in absolute value.

Table 3. Selected model.

Figure 14. Proportions (resp. normalized proportions) of explained pseudo-inertias according to $\gamma $ in the left panel (resp. in the right panel).

Figure 15. Fidelity curves according to $\gamma $ .

Figure 16. Left panel: GBM Gamma predictions by subset of the partition; Right panel: Relative differences.

We finally provide in the left panel of Figure 17 the QQ-plots comparing the distributions of the normalized observations inside each subsets of the estimated partitions with the unit Exponential distribution and in the right panel the box plots of the predicted values for each subset of the estimated partitions and for the covariates: LATITUDE, LENGTH, LONGITUDE, WIDTH, YEAR. We note that the empirical distributions inside the subsets of the estimated partition are not all well approximated by Exponential distributions. We observe that the covariates LENGTH and WIDTH are the covariates that most influence the means of observations through the subsets of the estimated partition.

Figure 17. Left panel: QQ-plots comparing the distributions of the normalized observations inside each subset of the estimated partitions with the unit Exponential distribution; Right panel: Box plots of the predicted values for each subset of the estimated partitions and for the covariates: $\text{LATITUDE, LENGTH, LONGITUDE, WIDTH, YEAR}$ .

Figure 18. Proportions (resp. normalized proportions) of explained pseudo-inertias according to $\gamma $ in the left panel (in the right panel).

Figure 19. Fidelity curves according to $\gamma $ .

For $D_{1}$ , we now take the Euclidean distance between the gravity centers of the rectangles. We obtain Figures 18 and 19.

On the basis of these plots, we chose $\gamma = 0.1$ and $K=12$ (the $R^{2}$ measure is equal to $89.44$ %.). Figure 20 displays box plots of the predicted values by GBM Gamma for each subset of the estimated partition, as well as box plots of their relative differences with their means. The ranges of the box plots are slightly larger than for the non-Euclidean distance.

Figure 20. Left panel: GBM Gamma’s predictions by subset of the partition; Right panel: Relative errors.

Figure 21. Left panel: QQ-plots comparing the distributions of the normalized observations inside each subset of the estimated partitions with the unit Exponential distribution; Right panel: Box plots of the predicted values for each subset of the estimated partitions and for the covariates: $\text{LATITUDE, LENGTH, LONGITUDE, WIDTH, YEAR}$ .

Figure 22. Fidelity curve.

We also provide in the left panel of Figure 21 the QQ-plots comparing the distributions of the normalized observations inside each subsets of the estimated partition with the unit Exponential distribution and in the right panel the box plots of the predicted values for each subset of the estimated partitions and for the covariates: LATITUDE, LENGTH, LONGITUDE, WIDTH, YEAR. The assumption of an Exponential distribution for each subset of the partition is now more convincing. Moreover, the box plots provide evidence of clearer links between the covariates considered and the empirical distributions of the observations through the subsets.

We finally consider the Equivalent tree method. Figure 22 displays the fidelity curve which led us to also choose $K=12$ subsets with the $R^{2}$ measure equal to $75.89$ %.

The tree that leads to the approximated partition is depicted in Figure 23. The covariates LENGTH and WIDTH appear to be the most discriminating variables. The year 1993 also appears to be an important year because a significant change in the tail index can be observed before and after this year. A comparison with the regression tree obtained by the methodology developed in Farkas et al. (Reference Farkas, Lopez and Thomas2021) is presented in Appendix B.

Figure 24 provides box plots of the predicted values by the equivalent tree for each subset of the estimated partition, as well as box plots of their relative differences with their averages. Figure 25 shows that the approximated partition of the equivalent tree gives distributions of normalized observations inside each subset that are relatively close to the unit Exponential distribution.

Figure 23. Surrogate tree ( $R^{2} = 75$ %).

Figure 24. Left panel: Equivalent tree predictions by subset of the partition; Right panel: Relative errors.

Figure 25. Left panel: QQ-plots comparing the distributions of the normalized observations inside each subset of the estimated partitions with the unit Exponential distribution; Right panel: Box plots of the predicted values for each subset of the estimated partitions and for the covariates: $\text{LATITUDE, LENGTH, LONGITUDE, WIDTH, YEAR}$ .

From these different analyses, we can draw the following conclusions:

1. – The small-size partition obtained from the additive tree ensemble and the hierarchical clustering with spatial constraints based on the Euclidean distance has a better fidelity than the one obtained with the tree surrogate model for the same number of classes.

2. – But the empirical distributions of the observations within each subset of the partitions are even so close to Exponential distributions for the tree surrogate model while providing natural explanations for the classes in terms of covariates and a simple division of the covariate space.