2.1 Quantities of Interest

We consider a portfolio of $n$ independent policies with data ${{\bf({D}}_i},{{\bf{Z}}_i},{N_i},{L_i},{e_i})$ , $\forall $ $i = 1, \ldots, n$ , where:

• ${\bf{D}}$ denotes a vector of unprotected policyholder information (e.g. vehicle age, area, etc.), which can be used to rate on,

• ${\bf{Z}}$ denotes protected policyholder variables (e.g. gender, ethnicity, etc.), which may not be used for pricing,

• $N$ and $L,$ respectively, denote the number of claims and the aggregated claims amount for each policy over its exposure $e$ .

In what follows, we will in turn use either the unprotected policyholder data matrix ${\bf{X}} = {\bf{D}}$ , or the combined $n \times p$ data matrix ${\bf{X}} = \left[ {{\bf{D}},{\bf{Z}}} \right]$ comprising of both ${\bf{D}}$ and ${\bf{Z}}$ . The nature of protected variables varies by jurisdiction and laws applicable in the country. In this paper, we consider fairness with respect to gender, a typical variable that is protected in Ireland. Factors such as gender, age and marital status are considered protected characteristics under EU anti-discrimination law. However, gender has a unique status among these protected characteristics; regardless of its commercial or underwriting relevance to insurance risks, it cannot be used in setting insurance prices (European Parliamentary Research Service, 2017).

Using this information, we are interested in evaluating the intrinsic policyholder risk

$${P_T} = \mathbb{E}(L|{\bf{D}},{\bf{Z}})$$

and a best feasible estimate of the expected loss for a policy

$${P_A} = \mathbb{E}(L|{\bf{D}})$$

Quantities ${P_T}$ and ${P_A}$ can be seen, respectively, as the “best estimate price” and “unaware price” (Lindholm et al., Reference Lindholm, Richman, Tsanakas and Wuthrich2022a), to which the company would apply a number of commercial margins and adjustments to account for expenses, levies, etc. Although ${P_T}$ and ${P_A}$ as defined above do not include any such additional margins and adjustments, they provide a basis for intrinsic evaluation of risk/loss. For the sake of simplicity, and as described in other works (Lindholm et al., Reference Lindholm, Richman, Tsanakas and Wüthrich2024b) in the rest of the paper, we refer to ${P_T}$ and ${P_A},$ respectively, as technical premium (TP) and actual premium (AP). In this work, we consider two distinct modelling strategies to estimate ${P_T}$ and ${P_A}$ , namely, the Poisson-gamma (product) model, and the Tweedie model.

2.2 Poisson-gamma (Product) Model

The first modelling strategy considered is the most prevalent, Poisson-gamma (product) model. This strategy consists of estimating frequency and severity as two separate steps. Claims frequencies are typically modelled by a discrete count distribution. The Poisson distribution being a common choice (Henckaerts et al., Reference Henckaerts, Côté, Antonio and Verbelen2021; Noll et al., Reference Noll, Salzmann and Wuthrich2020) and the distribution used hereafter, such that for policyholder $i \in \left\{ {1, \ldots, n} \right\}$ ,

$${N_i}\sim Poi\left( {{\mu _i}{e_i}} \right),{\rm{\;\;\;\;}}{\mu _i} \gt 0,\;\;{e_i} \gt 0$$

for some rate ${\mu _i}$ . Then the claims frequency $F$ over a unit period of exposure (say, one-year exposure), i.e. the number of claims $N$ proportional to exposure $e$ , is such that

$$\mathbb{E}\left( F \right) = \mathbb{E}\left( {{N \over e}\,\left|\, e \right \gt 0} \right)$$

or, for a given policy, evaluated as

$${F_i} = {{{N_i}} \over {{e_i}}}$$

such that $\mathbb{E}\left( {{F_i}} \right) = {\mu _i}$ and $Var\left( {{F_i}} \right) = {{{\mu _i}} \over {{e_i}}}$ , i.e. ${F_i}$ is pseudo-Poisson random variable (either over or under-dispersed depending on the sign of ${e_i} - 1$ ; usually ${e_i} \le 1$ ). In practice, we can estimate ${\mu _i}$ as the maximum likelihood estimator (MLE)

$${\hat \mu _i} = {F_i}$$

A conventional approach consists of fitting a GLM to the policyholder information $\left( {{\bf{D}},{\bf{Z}},N,e} \right)$ to explain claims frequencies for individual profiles as

(1) $${\rm{log}}\mathbb{E}(F|{\bf{X}}) = {g_F}\left( {\bf{X}} \right) = {\bf{X}}{\beta _F}$$

with real-valued regression coefficients ${\beta _F} \in {\mathbb{R}^p}$ . Model (1) may be fit to the data by minimising the Poisson deviance between claims frequency data $\{ {y_i}\} _{i = 1}^n$ (e.g. where $Y = F$ ) and a set of model predictions $\{ {g_F}\left( {{{\bf{x}}_i}} \right)\} _{i = 1}^n$ obtained from policyholder data ${\bf{X}}$ , defined as in (Henckaerts et al., Reference Henckaerts, Côté, Antonio and Verbelen2021; Venables & Ripley, Reference Venables and Ripley2013) by

(2) $${{\cal D}_F}\left( {y,{g_F}\left( {\bf{X}} \right)} \right) = 2\mathop \sum \limits_{i = 1}^n \left( {{y_i}log\left( {{{{y_i}} \over {{g_F}\left( {{{\bf{x}}_i}} \right)}}} \right) - \left( {{y_i} - {g_F}\left( {{{\bf{x}}_i}} \right)} \right)} \right)$$

Let us then define claim severity $S$ as the aggregated loss (or claim amount) $L$ divided by the number of claims from a subset of polices that have at least one or more claims reported and paid out, such that

$$\mathbb{E}\left( S \right) = \mathbb{E}\left( {{L \over N}\,\left|\, N \right\gt 0} \right)$$

or, for a given policy $i$ , evaluated as

$${S_i} = {{{L_i}} \over {{N_i}}}$$

Hereafter, following a conventional approach, $S$ is modelled by a gamma distribution to allow for skewness arising naturally in such data. Similar to equations (1) and (2), claim severities may be modelled using a distinct GLM, that we may denote as ${g_S}\left( X \right) = {\bf{X}}{\beta _S}$ with ${\beta _S} \in {\mathbb{R}^p}$ . This model may be achieved by minimising the gamma deviance between the severity observations $\{ {y_i}\} _{i = 1}^n$ (e.g. where $Y = S$ ) and a set of model predictions $\{ {g_S}\left( {{{\bf{x}}_i}} \right)\} _{i = 1}^n$ defined as in (Henckaerts et al., Reference Henckaerts, Côté, Antonio and Verbelen2021) by

(3) $${{{\cal D}_S}\left\{ {y,{g_S}\left( {\bf{X}} \right)} \right\} = 2\mathop \sum \limits_{i = 1}^n {w_i}\left( {{{{y_i} - {g_S}\left( {{{\bf{x}}_i}} \right)} \over {{g_S}\left( {{{\bf{x}}_i}} \right)}} - {\rm{log}}\left( {{{{y_i}} \over {{g_S}\left( {{{\bf{x}}_i}} \right)}}} \right)} \right)}$$

where ${w_i}$ is the weight for each observation. Hereafter we use ${w_i} = {N_i}$ .

The overall policy premium $P$ may then be defined as the expected loss $L$ per unit of exposure, such that

$$P = \mathbb{E}\left( {{L \over e}} \right)$$

under the assumptions of a non-zero chance that a claim be made, and that claim frequencies and severities are independent of each other. A policy premium estimate $\hat P$ over a unit exposure may then be obtained as the product of estimates for frequency $\hat F$ and severity $\hat S$ (Denuit et al., Reference Denuit, Maréchal, Pitrebois and Walhin2007; Frees et al., Reference Frees, Meyers and Cummings2014; Parodi, Reference Parodi2014), i.e.

$$\hat P = \hat F \times \hat S$$

such that (Henckaerts et al., Reference Henckaerts, Côté, Antonio and Verbelen2021)

$$\mathbb{E}\left( {\hat P} \right) = P = \mathbb{E}\left( {F \times S} \right) = \mathbb{E}\left( {{N \over e}} \right) \times \mathbb{E}\left( {{L \over N}\,\left|\, N \right\gt 0} \right).$$

2.3 Tweedie Model

A relatively common alternative to the Poisson-gamma product model consists of modelling the quantity $\mathbb{E}\left( {\hat P} \right)$ directly as a compound variable under the Tweedie distribution (Quijano Xacur & Garrido, Reference Quijano Xacur and Garrido2015; Shi, Reference Shi2016). This member of the exponential dispersion family is characterised by a variance function of the form ${\rm{Var}}\left( Y \right) = \phi {\gamma ^\alpha }$ , where $\phi $ is the dispersion parameter, $\gamma $ is the mean, and $\alpha $ is the power parameter such that $\alpha \ne 0$ and $\alpha \ne 1$ . Given a set of $n$ independent observations $\left\{ {{y_1}, \ldots, {y_n}} \right\}$ of a Tweedie random variable $Y,$ an estimate of the premium $P$ is derived by minimising the Tweedie deviance between the observations $\{ {y_i}\} _{i = 1}^n$ and a set of model predictions $\{ {g_T}\left( {{{\bf{x}}_i}} \right)\} _{i = 1}^n$ , defined as

(4) $${{\cal D}\left( {y,{g_T}\left( {\bf{X}} \right)} \right) = \mathop \sum \limits_{i = 1}^n {2 \over \phi }\left( {{{y_i^{2 - \alpha }} \over {\left( {1 - \alpha } \right)\left( {2 - \alpha } \right)}} - {{{y_i}{g_T}{{({{\bf{x}}_i})}^{1 - \alpha }}} \over {1 - \alpha }} + {{{g_T}{{({{\bf{x}}_i})}^{2 - \alpha }}} \over {2 - \alpha }}} \right)}$$

with respect to $\left( {\phi, \gamma, \alpha } \right)$ and the model parameters for ${g_T}$ .

2.4 Premium Modelling

A number of models $g\left( {\bf{X}} \right)$ using the Poisson, gamma or Tweedie deviance ${\cal D}$ defined, respectively, by equations (2), (3) or (4), as cost function were compared. They include the conventional GLM, with a structure similar to equation (1), allowing for interaction terms (McCullagh & Nelder, Reference McCullagh and Nelder1989; Ohlsson & Johansson, Reference Ohlsson and Johansson2010) in a forward stepwise selection of the final model ESL; an elastic net and a LASSO (Tibshirani, Reference Tibshirani1996; Zou & Hastie, Reference Zou and Hastie2005), using $\delta = 0.5$ and $\delta = 1$ , respectively, and a regularisation parameter $\lambda \gt 0$ in

$${{\cal D}_\lambda }\left( {\hat \beta } \right) = \mathop \sum \limits_{i = 1}^n \,{\cal D}{({y_i},{{\bf{x}}_i}{{\hat \beta} _j})^2} + \lambda \mathop \sum \limits_{j = 1}^p \,\left( {{1 \over 2}\left( {1 - \delta } \right)\beta _j^2 + \delta \left| {{\beta _{\rm{j}}}} \right|} \right)$$

Furthermore, they also include a pruned regression tree (Breiman et al., Reference Breiman, Friedman, Stone and Olshen1984; Friedman, Reference Friedman2001) partitioning the dataset into $J$ leafs $R = \left\{ {{R_1}, \ldots, {R_J}} \right\}$ and with complexity penalisation $c \gt 0$ as described in (Henckaerts et al., Reference Henckaerts, Côté, Antonio and Verbelen2021) by optimising

where ${\hat y_R}$ denotes the average value of $Y$ within a leaf or set of leaves $R$ as defined above; a RF (Breiman, Reference Breiman2001) using an ensemble of $M$ trees; a gradient boosting model (GBM) (Friedman, Reference Friedman2001) using a penalised loss gradient at each step $m = 2, \ldots, M$ , the latter being defined by

$$ - {g_m}\left( {\bf{x}} \right) = - {\left[ {{{\partial {\cal D}\left( {y,g\left( {\bf{x}} \right)} \right)} \over {\partial g\left( {\bf{x}} \right)}}} \right]_{g\left( {\bf{x}} \right) = {g_{m - 1}}\left( {\bf{x}} \right)}}$$

as well as extreme gradient boosting (XGB) (Chen & Guestrin, Reference Chen and Guestrin2016), penalising (a second-order linearisation of) the deviance by

$$\mathop \sum \limits_{i = 1}^n \,{\cal D}\left( {{y_i},g\left( {{{\bf{x}}_i}} \right)} \right) + \xi J + {\tau \over 2}\left| {\left| \omega \right|} \right|$$

where $J$ is the number of terminal nodes or leaves in a tree, $\xi \gt 0$ a user defined penalty encouraging pruning, and for some regularisation $\tau \gt 0$ of the final leaf weights $\omega $ determined for the ensemble model by stochastic gradient descent. XGB in particular has gained a lot of popularity in the actuarial field in recent years, some research reporting it to outperform other models for insurance claims and frequency prediction (Colella & Jones, Reference Colella and Jones2023; Fauzan & Murfi, Reference Fauzan and Murfi2018; Ferrario & Hämmerli, Reference Ferrario and Hämmerli2019; Henckaerts et al., Reference Henckaerts, Côté, Antonio and Verbelen2021; Kuo & Lupton, Reference Kuo and Lupton2020).

2.5 Performance Evaluation

2.5.2 Post-Modelling Calibration

All model predictions (premium estimates) were calibrated to adjust for model bias before they were used to carry out analysis. The calibration factor $C$ was calculated based on the model fits obtained on the training dataset and is evaluated as

$$C = {{\mathop \sum \nolimits_{i = 1}^n \,{L_i}} \over {\mathop \sum \nolimits_{i = 1}^n \,{{\hat P}_i}}}$$

2.5.3 Premium Estimation Error and Bias

All models were first compared on the test set with respect to bootstrapped Poisson, gamma and Tweedie deviances as defined in equations (2), (3) or (4), to evaluate overall model error. Model bias was also evaluated on the test data as the difference between the total sum of observed losses and the premium estimates obtained from each model, defined as

$$\mathop \sum \limits_{i = 1}^n \,\left( {{L_i} - {{\hat P}_{T,i}}} \right)$$

2.5.4 Risk Stratification

Lorenz curves were used to compare premiums by analysing the distribution of premiums versus losses, where both quantities were ordered by increasing relativities

$${r_i} = {{{{\hat Q}_i}} \over {{{\hat P}_i}}},{\rm{\;}}i = 1, \ldots, n$$

i.e. by the ratio of the premium prediction obtained from the competing model ( $\hat Q$ ) to that of the benchmark model ( $\hat P$ ) (Henckaerts et al., Reference Henckaerts, Côté, Antonio and Verbelen2021), whose empirical c.d.f. is denoted hereafter by ${F_n}$ . Sorting policies in terms of their increasing relativity, the ordered Lorenz curve then has coordinates $\left( {{{\hat F}_L},{{\hat F}_{\hat P}}} \right)$ defined by

$${\hat F_{\hat P}}\left( s \right) = {{\mathop \sum \nolimits_{i = 1}^n \,\hat P\left( {{{\bf{x}}_i}} \right)I\left( {{F_n}\left( {{r_i}} \right) \le s} \right)} \over {\mathop \sum \nolimits_{i = 1}^n \,\hat P\left( {{{\bf{x}}_i}} \right)}},{\rm{\;\;}}s \in \left[ {0,1} \right]$$

and

$${\hat F_L}\left( s \right) = {{\mathop \sum \nolimits_{i = 1}^n \,{y_i}I\left( {{F_n}\left( {{r_i}} \right) \le s} \right)} \over {\mathop \sum \nolimits_{i = 1}^n \,{y_i}}},{\rm{\;\;}}s \in \left[ {0,1} \right]$$

where $I\left( . \right)$ is the indicator function. The curve $\left( {{{\hat F}_L},{{\hat F}_{\hat P}}} \right)$ follows the identity line $\left( {y = x} \right)$ in the case of benchmark estimation. The Gini index is calculated as twice the area between the ordered Lorenz curve and the identity line (Frees et al., Reference Frees, Meyers and Cummings2014; Henckaerts et al., Reference Henckaerts, Côté, Antonio and Verbelen2021). Larger index indicates greater risk differentiation can be obtained by using the competing premium $\hat Q$ compared to the reference premium $\hat P$ , and hence the competing pricing model is more profitable. This case corresponds to a concave Lorenz curve. The profitability point is as per the theory (Frees et al., Reference Frees, Meyers and Cummings2014) that the profitability is the area under the Lorenz curve. However, if Company B has a better estimate of elasticity and conversion/retention rates then they could operate more profitably than Company A, even if A has the more profitable pricing model. In this analysis, we assume that the price elasticity modelling capabilities between companies are similar and they do not play a role.

2.5.6 Sensitivity Analysis

2.5.6.1 Variable Importance

Variable importance was assessed to understand which variables drive premium prediction in each of the models considered. Any discrepancies in the variable importance values may indicate specific model sensitivity to the information collected on policyholders. Variable importance in regression models (GLM and its regularised alternatives) was measured directly from the fitted coefficients by normalising their magnitude (in absolute value) by the sum of all model coefficient magnitudes (Fryda, Reference Fryda2023). Variable importance in a tree was derived by taking the sum of improvements in the loss function at each split (Breiman et al., Reference Breiman, Friedman, Stone and Olshen1984)

where $\eta $ denotes the number of internal nodes (or splits) of the tree, ${v_\ell } \in {\bf{X}}$ is the variable used to perform the ${\ell ^{{\rm{th}}}}$ split, and ${({\rm{\Delta }}{\cal D})_\ell }$ measures the magnitude of the difference in deviance ${\cal D}$ (defined by equations (2) or (3) or (4)) before and after split $\ell $ . For the other tree-based models, namely RF, GBM and XGB, variable importance was calculated as the average of importances measured for the variable over all trees in the ensemble (Breiman et al., Reference Breiman, Friedman, Stone and Olshen1984; Hastie et al., Reference Hastie, Tibshirani and Friedman2009). Importance analysis was also used to assess the impact of protected variables on actual premiums derived from the various models.

2.5.6.2 Model Agnostic Evaluation

The impact of individual features on model predictions was further evaluated in model-independent ways, using SHapley Additive exPlanations (SHAP) analyses and partial dependence plots. For a prediction function $g$ , the SHAP value ${\phi _j}$ for feature ${X_{\left( j \right)}}$ is a feature importance value computed as the average marginal contribution of feature ${X_{\left( j \right)}}$ across all possible subsets of features (Lundberg & Lee, Reference Lundberg and Lee2017) as

where ${\bf{X}}$ is the whole feature set, $S$ is a subset of features not containing ${X_{\left( j \right)}}$ .

and $g\left( S \right)$ denote the premium prediction function evaluations obtained when only the features in the respective subset

and subset $S$ are known. In practice, $g\left( S \right)$ can be evaluated by integrating over all possible values of ${X_{\left( j \right)}}$ , i.e. $g\left( S \right) \approx E[g\left( {\bf{X}} \right)|S]$ , typically by sampling the missing features conditionally on the known ones ( $S$ ) (Molnar, Reference Molnar2022).

A partial dependence plot (PDP) is a popular model-agnostic visualisation tool for model interpretability. It depicts the marginal effect of one or more input features on the predicted outcome, averaging over the joint distribution of the remaining features (Friedman, Reference Friedman2001; Molnar, Reference Molnar2022). Formally, for a premium prediction function $g\left( x \right)$ and a feature ${X_{\left( j \right)}} \in {\bf{X}}$ , the partial dependence function is defined as

This marginalisation estimates the expected prediction when is fixed and all other features vary according to the empirical data distribution. For continuous features, we discretise the feature range (using deciles) into a grid of values and compute the partial dependence at each grid point.

2.5.7 Defining Actual and Technical Premiums

The ratio of actual premium to technical premium (APTP)

(5) $$\rho = {{{P_A}} \over {{P_T}}}$$

may be used to assess how the removal of protected information from the model may impact the unaware price ${P_A}$ relative to the best estimate risk ${P_T}$ . An APTP ratio greater than 1 indicates a higher premium than required to cover the expected costs of the policy for that policyholder profile. On the other hand, an APTP ratio lower than 1 indicates the customer paid less than the expected cost of insuring the risk. Univariate analyses of APTP ratios were carried out to illustrate how some typical protected variables affected the pricing methodologies listed in Section 2.4. Since our quantitations of ${P_A}$ and ${P_T}$ are raw model output (see Section 2.1), $\rho $ may be seen as a proxy for an APTP ratio.

2.5.8 Statistical Distance Between Populations

In order to compare premium distributions across fairness groups, we used the Wasserstein distance (Pichler, Reference Pichler2014), with order 1. Given two premium estimators $P$ and $P{\rm{'}}$ , with probability distributions ${F_P}$ and ${F_{P{\rm{'}}}}$ , this distance is defined as

(6) $$W\left( {P,P{\rm{'}}} \right) = \mathop \smallint \nolimits_{ - \infty }^\infty \,\left| {{F_P}\left( u \right) - {F_{P{\rm{'}}}}\left( u \right)} \right|du$$

In practice, this quantity can be evaluated by using the empirical distribution functions for sets of premium estimates $\{ {P_i}\} _{i = 1}^n$ and $\{ {P^{\prime}_{i{\rm{}}}}\} _{i = 1}^n$ , i.e. ${F_Q}\left( u \right) = {1 \over n}\sum\nolimits_{i = 1}^n {\mkern 1mu} I\left( {{Q_i} \le u} \right)$ for either $Q = P$ or $Q = P{\rm{'}}$ , in a numerical approximation of the integral term. Wasserstein distance as defined in equation (6) is an actual metric that measures the area between the marginal distributions ${F_P}$ and ${F_{P{\rm{'}}}}$ (Chhachhi & Teng, Reference Chhachhi and Teng2023; De Angelis & Gray, Reference De Angelis and Gray2021). The Wasserstein distance yields a positive value that increases with the distance between the two distributions that is bounded upwards. In this analysis, we use the first order of the Wasserstein distance to mitigate sensitivity to the natural skewness in the distributions of premiums (Chhachhi & Teng, Reference Chhachhi and Teng2023).

2.5.9 Characterisation of Fairness

Fairness is commonly assessed by analysing actual premium distributions with respect to protected variables. The nature of the departure from a fair actual premium may be investigated in more detail, especially to assess whether it affects a specific protected group. Building on prior studies (Dolman & Semenovich, Reference Dolman and Semenovich2018; Lindholm et al., Reference Lindholm, Richman, Tsanakas and Wuthrich2022a; Mosley & Wenman, Reference Mosley and Wenman2022; Xin & Huang, Reference Xin and Huang2024), a framework for the evaluation model fairness can be extended to assess compliance of the actual premiums ${P_A}$ with the following three axiom.

Axiom 1 (Demographic Parity). The actual premium ${P_A}$ must not depend on the protected features. In particular, where the feature is a categorical variable with $L$ levels, $Z \in \left\{ {{l_1}, \ldots, {l_L}} \right\}$ (e.g. gender), then it should be the case that

$$\mathbb{P}({P_A}|Z = {l_i}) = \mathbb{P}({P_A}|Z = {l_j}){\rm{\;}}\forall i \ne j$$

Axiom 2 (Actuarial Group Fairness). The distribution of actual premium ${P_A}$ conditional on technical premium ${P_T}$ must not depend on the protected information $Z$ , i.e.

$$\mathbb{P}({P_A}|{P_T},Z = {l_i}) = \mathbb{P}({P_A}|{P_T},Z = {l_j}){\rm{\;}}\forall i \ne j$$

Axiom 3 (Calibration). The distribution of actual loss $Y$ conditional on actual premium ${P_A}$ must be conditionally independent of any protected feature $Z$ , i.e.

$$\mathbb{P}(Y|{P_A},Z = {l_i}) = \mathbb{P}(Y|{P_A},Z = {l_j}){\rm{\;}}\forall i \ne j$$

Without loss of generality, $Z$ could be a numerical variable and the conditional distributions may be evaluated at discretised $Z$ -value bins.

Alignment with Axiom 1 may be assessed quantitatively by measuring the magnitude of disparities between probability distributions $\mathbb{P}({\hat P_A}|Z = {l_i})$ and $\mathbb{P}({\hat P_A}|Z = {l_j})$ for any two groups ${l_i}$ and ${l_j}$ of a given protected variable $Z$ . In our framework, we propose using the Wasserstein distance between the distributions of two subgroups of policyholders. As an example, for a binary categorical variable $Z \in \left\{ {{l_1},{l_2}} \right\}$ over the ${n_{test}}$ data points:

(7) $${d_Z} = W\left( {{{\hat P}_A}\left| {Z = {l_1},{{\hat P}_A}} \right|Z = {l_2}} \right)$$

Similarly, adherence to Axiom 2 may be measured on the basis of the Wasserstein distance in the test dataset for a given category $\left\{ {Z = l} \right\}$ , and particular TP percentiles of interest. As an example, subsetting the dataset into bands defined by the TP quintiles $\left\{ {{{\tilde P}_{T,0}},{{\tilde P}_{T,1}}, \ldots, {{\tilde P}_{T,5}}} \right\}$ such that

(8) $${{\cal P}_{\cal T}} \in \left\{ {\left( {{{\tilde P}_{{T_{k - 1}}}},{{\tilde P}_{{T_k}}}} \right],k \in \left\{ {1,2,3,4,5} \right\}} \right\},$$

to measure adherence for each band we calculate the Wasserstein distance as:

(9) $${d_{l,{{\cal P}_{\cal T}}}} = W\left( {{{\hat P}_A}\left| {\left( {{{\hat P}_T} \in {{\cal P}_{\cal T}},Z = {l_1}} \right),{{\hat P}_A}} \right|\left( {{{\hat P}_T} \in {{\cal P}_{\cal T}},Z = {l_2}} \right)} \right)$$

Adherence to Axiom 3 may be measured on the basis of the Wasserstein distance of observed losses in the test dataset between protected groups in a given actual premium quintile band, as follows:

(10) $${d_{L,Z}} = W\left( {L\left| {{{\hat P}_A},Z = {l_1},L} \right|{{\hat P}_A},Z = {l_2}} \right)$$

In each case, a lower Wasserstein distance from equations (7), (9) or (10) indicates better conformity to the axioms. In addition, two-sided Mann–Whitney tests were performed at the 5% significance level to assess whether any of the conditional distributions compared in equations (7), (9) and (10) were statistically different. These quantitations are reported on in Section 4.1.

Other metrics may be considered to complement indicators from equations (7), (9) and (10). Quantitatively, in what follows we also evaluate the distance between distribution medians, to provide additional context around the value of the Wasserstein distance. Complementary assessment may also be performed qualitatively. For the qualitative assessment, we consider evaluating the statistical significance of a nonparametric Wilcoxon test applied to the samples of premium estimates across protected groups.

2.6 Implementation

All analyses were carried out using R version 4.2.0 (R Core Team, 2021) and its H2O interface (Fryda, Reference Fryda2023). Some of our implementation builds upon source code developed by (Henckaerts et al., Reference Henckaerts, Côté, Antonio and Verbelen2021), available at https://github.com/henckr/distRforest. The RF and regression tree models rely on the distRforest (Henckaerts et al., Reference Henckaerts, Côté, Antonio and Verbelen2021) and rpart (Therneau & Atkinson, Reference Therneau and Atkinson2019) R packages. Dedicated packages were also used for model glmnet (Friedman et al., Reference Friedman, Hastie and Tibshirani2010; Zou & Hastie, Reference Zou and Hastie2005), GBM (Ridgeway, Reference Ridgeway2014; Southworth, Reference Southworth2015) SHAP analyses and PDPs were implemented using the iml package in R (Molnar, Reference Molnar2022), which supports both individual and grouped feature effects. To calculate the Wasserstein distance, the R package Transport (Dominic et al., Reference Dominic, Björn and Nicolas2023) was used. All source code for our analyses is available in open access on Github at https://github.com/tisrani/Pricing.

3.1 Datasets

For this study we used two open-access datasets which comprised of both frequency and severity. These datasets yielded the opportunity to create a benchmark on two distinct types of cover namely: third-party liability and fully comprehensive cover. Furthermore, they allowed exploring the impact of different sets of covariates such as driver age. Additionally, Dataset 2 contained information on protected variables, such as gender, that allowed for a realistic benchmark on fairness.

Dataset 1 comprised of two parts, namely freMTPL2freq and freMTPL2sev from the R package CASdatasets (Charpentier, Reference Charpentier2014; Dutang & Charpentier, Reference Dutang and Charpentier2015), containing data from a French motor third-party liability book of business over one calendar year. The merged dataset contained the numbers of claims and aggregated amounts incurred by each of 678,013 policies, of which 24,944 had non-zero claims. Dataset 1 also included information on the policy (bonus-malus), main driver’s age and location (area, density and region), and vehicle (age, brand, power, fuel type). We refer the reader to (Noll et al., Reference Noll, Salzmann and Wuthrich2020) for a more detailed description of Dataset 1. Dataset 1 has also been used in a number of other works (Ciatto et al., Reference Ciatto, Verelst, Trufin and Denuit2023; Delcaillau et al., Reference Delcaillau, Ly, Papp and Vermet2022; Denuit et al., Reference Denuit, Charpentier and Trufin2021; Ferrario et al., Reference Ferrario, Noll and Wuthrich2020; Krasniqi et al., Reference Krasniqi, Bardet and Rynkiewicz2022; Lorentzen & Mayer, Reference Lorentzen and Mayer2020; Schelldorfer & Wuthrich, Reference Schelldorfer and Wuthrich2019; Su & Bai, Reference Su and Bai2020).

Dataset 2 also had two parts, namely, pg17trainpol and pg17trainclaim (Dutang & Charpentier, Reference Dutang and Charpentier2015) from the R package CASdatasets. Dataset 2 contains observations of 64,515 fully comprehensive cover policies with one full year of exposure. Out of the 64,515 policies, there were claims in 8,648 of them. Since there was no information on the type of claims, we could not identify how much of the losses were related to third-party claims. Hence, an additional 35,485 third-party-only liability policies were removed from this second dataset. This yielded a more statistically balanced dataset, and also allowed for a more effective benchmark design, because now Datasets 1 and 2 thus contained a single and distinct type of policy cover. Dataset 2 also included variables on the policy (bonus-malus, coverage, payment frequency, pay-as-you-drive), main driver (INSEE codeFootnote ¹ , age, gender, number of years with full licence), and vehicle (type, make, model, age, number of cylinders, DINFootnote ² , fuel type, maximum speed vehicle can achieve, weight and value). Information on a second driver and other secondary variables were ignored to simplify the analysis. Furthermore, as the vehicle age was already known, the vehicle sales dates were ignored. Dataset 2 has also been used in other works (Bove et al., Reference Bove, Aigrain, Lesot, Tijus and Detyniecki2022; Brouste et al., Reference Brouste, Dutang and Rohmer2024; Havrylenko & Heger, Reference Havrylenko and Heger2022; Simjanoska, Reference Simjanoska2022).

Claims with extremely low and high amounts outside the range of (€50, €10,000) were removed in both datasets. Following this, Datasets 1 and 2 contained $n = 24,123$ and $n = 8,500$ policies with non-zero claims, respectively. The base level for each factor (e.g. vehicle gas) was set at the level with the highest number of policies.

Following a conventional actuarial representation of policyholder information, categorical variables (for Dataset 1, region, vehicle brand and area; for Dataset 2, vehicle make and vehicle model) were recategorised into broader, well-defined groups, to select key differentiating features and use interpretable modelling approaches to enhance the significance and understanding of predictions for each group. Here, a tree-based binning strategy (Henckaerts et al., Reference Henckaerts, Côté, Antonio and Verbelen2021) was applied to create categorical risk factors with levels being separated based on the observed claims values at each level. Tree-based binning strategy is typically used to improve the generalisability of estimated premiums and to simplify the model in line with the parsimony principle. The distributions of these recategorised variables are provided in Appendix A. Following this step, information on vehicle type in Dataset 2 was removed due to its high correlation with other vehicle features based on Spearman’s or Cramer’s correlation, as appropriate.

The overall distributions of claims numbers ( $N$ ) and amounts ( $L$ ) in the datasets after data cleaning and pre-processing are shown in Figure 1 for both datasets. These were found to be globally similar in shape, although the distribution of claims amount in Dataset 1 has a relatively fatter tail.

Figure 1. Distributions of claims amounts (top) and number of claims (bottom) for Dataset 1 (left) and Dataset 2 (right).

3.2 Which Factors Drive the Premiums?

3.2.1 Variable Importance

Analyses of variable importance showed both similarities and differences between models, influenced by the nature of modelling strategy: the contrast between linear regression and tree-based modelling strategy. These similarities and discrepancies can be observed in Table 1. The figure also indicates that only a small number of variables effectively drove any of the models, particularly for the estimation of frequency in both datasets, and of severity in Dataset 2.

Table 1. Variable importance metrics for frequency (left), severity (centre) and Tweedie (right) models trained on Dataset 1 (top) and Dataset 2 (bottom), as defined in Section 2.5.6. The greyed-out cells indicate variables that were removed during the stepwise selection in the GLM model building process. A value of 0% in the table corresponds to an importance below 0.5%

Considering the product modelling approach in Dataset 1, overall the variable importance values and patterns were similar across models for frequency estimation, with only minor discrepancies across model types. This was also largely the case for severity estimation, although the XGB model tended to leverage more variables than the other models. A similar conclusion could also be drawn for the Tweedie modelling approach in Dataset 1. Under the Tweedie approach, variable importance values from XGB tended to be more spread out across the top 6 variables in the table, and regression models in this framework used more variables. For the regression model, increased importance was observed for variables area and vehicle gas.

In Dataset 2, different patterns of variable importance were observed between regression-based and tree-based techniques. However, similar patterns of variable importance were observed within each modelling family. In the product modelling approach, tree-based models placed more emphasis on vehicle value, age and weight for frequency estimation, while vehicle make was more important in regression models. For severity estimation, tree-based techniques placed higher importance to vehicle weight and engine cylinder, whilst bonus-malus was a large contributor in regression settings. Under the Tweedie modelling approach, vehicle value had a larger contribution in tree-based models; on the other hand, variables vehicle make, pay-as-you-drive and vehicle cylinder variables were more important in the regression models.

3.2.2 SHAP Analysis

The SHAP analysis presented in Table 2 offers additional insights into feature importance, and allows for a direct comparison of the product and Tweedie modelling frameworks.

Table 2. SHAP value metrics for Product (left) and Tweedie (right) models trained on Dataset 1 (top) and Dataset 2 (bottom), as defined in Section 2.5.6. The greyed-out cells indicate variables that were removed during the stepwise selection in the GLM model-building process. A value of 0% in the table corresponds to an SHAP value below 0.5%

For the product framework (Table 2, left), distinct patterns emerge. For Dataset 1 (top-left panel), the two regularised models have identical SHAP profiles. In contrast, GLM displays a different pattern, with more evenly distributed SHAP values. This differentiation was not seen in the model-based variable importance analysis of Section 3.2.1, and suggests that regularisation may enhance strong predictors for this data, in a way that is not impacted by collinearity (since the patterns for LASSO and Elastic Net are almost identical). RF and GBM yield comparable SHAP profiles, suggesting a shared interpretation of the feature space, but different from that of XGB, which again differs from a common trend seen in Section 3.2.1 across all tree-based models. Bonus-malus is relevant in the regularised models (16%) and in XGB (28%), along with at least one other predictor (density, vehicle power and vehicle gas for the linear models; vehicle and driver ages for XGB). However, for Dataset 2 (bottom-left panel), the regularised linear models yield different SHAP profiles, as do RF and GBM. In most models, the distributions of SHAP values are more spread out for Dataset 2. Density and vehicle speed are influential predictors for a number of models, but there are no clear strong predictors across all models. Noticeably, gender contributes around 10% of the risk evaluation in elastic net and GBM, which would potentially raise ethical and fairness implications.

For the Tweedie modelling framework (Table 2, right), we observe different SHAP-based variable rankings, and more consistency across all ML models (for both datasets), whereas SHAP values for GLM are spread more uniformly across the range of features. The fact that SHAP profiles differ so much between the product and Tweedie frameworks emphasises how the choice of loss function alters the learning dynamics, leading to a redistribution of predictive importance across features. For Dataset 1 (top-right panel), bonus-malus in particular yields high importance across all ML models. For Dataset 2 (bottom-right panel), vehicle age and vehicle value arise as important predictors. Density and bonus-malus appear with moderate importance across most models. As driver age and gender appear meaningful in a number of models, this may raise issues around pricing fairness.

3.2.3 Partial Dependence Analysis

Figure 2 presents a selection of partial dependence patterns that emerge when comparing the product and Tweedie framework predictions across models for Dataset 1. Both product and Tweedie models demonstrate a sharp increase in average prediction with increasing bonus-malus (top left), with slightly stronger gradients for Tweedie XGB and RF at higher premiums. This greater sensitivity is consistent with the SHAP-based findings. For both pricing frameworks, all models yielded comparable variations in average premiums with respect to a number of variables, such as region (bottom left). However, more variability was observed across models and/or frameworks for certain variables, such as vehicle age (top right) and density (bottom right). For instance, for vehicle age, trends differed across models, but those differences were globally consistent across product and Tweedie frameworks. For density, however, variations in average prediction appeared consistent across all models in the product framework, but different in the Tweedie framework. As these variables appeared important (see Section 3.2.1), the subset of PDP graphs shown in Figure 2 is sufficient to highlight the importance of an independent sensitivity analysis of any given model with respect to the information available for pricing. In summary, the interpretation of several variables in terms of the pricing outcome is inconsistent across models.

Figure 2. Partial dependence plot by feature level for categorical and discretised (by deciles) numerical variables for models trained on Dataset 1. Each subplot shows the average predicted premium across levels of a single feature, under the product (solid line) and Tweedie (dashed line) modelling strategies.

For Dataset 2 (Figure 3), trends were comparable across most variables, and in particular, for all numerical variables (left and top-right panels). The magnitude of the change in average premiums with respect to changes in levels for a few categorical variables showed more significant differences, although trends always had the same direction. This was the case, for example, for payment frequency (bottom-right panel).

Figure 3. Partial dependence plot by feature level for categorical and discretised (by deciles) numerical variables for models trained on Dataset 2. Each subplot shows the average predicted premium across levels of a single feature, under the product (solid line) and Tweedie (dashed line) modelling strategies.

3.3 Deviance and Bias

Figure 4 shows the bootstrap distributions of Poisson, gamma and Tweedie deviances, obtained from the test subsets for Datasets 1 and 2. Two-sample one-sided Mann–Whitney tests were carried out at the 5% significance level on the bootstrapped test deviances to assess statistical significance of relative differences in performance between any two models.

Figure 4. Bootstrap distributions of test Poisson (top), gamma (middle) and Tweedie (bottom) deviances, for Datasets 1 (left) and 2 (right). The solid blue line indicates the deviance measure from a null model (average value observed from training data).

For estimation of claims frequency, XGB was found to yield the lowest Poisson deviance overall, outperforming all other models for Dataset 1 ( $p \lt {10^{ - 4}}$ ). For Dataset 2, both regularised models and the RF model performed at par with GLM (all $p \gt 0.05$ ). However, the boosting models performed significantly better than the GLM (GBM: $p \lt {10^{ - 4}}$ , XGB: $p = 0.0356$ ).

For estimation of claims severity, GLM performed at par with all other models ( $p \gt 0.05$ ) for Dataset 1. For Dataset 2, regularised models performed comparably with GLM (LASSO: $p = 0.2871$ , E.Net: $p = 0.2829$ ). However, the boosting models and RF performed significantly better than GLM (all $p \lt {10^{ - 4}}$ ).

In terms of Tweedie deviance, GLM and regularised models performed at the same level (LASSO: $p = 0.9817$ and E.Net: $p = 0.1699$ ). However, the tree-based models yielded a lower Tweedie deviance than GLM ( $p \lt {10^{ - 6}}$ ). Additionally, the boosting models outperformed the RF ( $p \lt {10^{ - 6}}$ ). For Dataset 2, GLM, RF and regularised models perform at par (all $p \gt 0.05$ ). However, the boosting models performed significantly better than GLM ( $p \lt {10^{ - 6}}$ ).

In all comparisons, Poisson, gamma and Tweedie deviances from GLM estimates were either significantly larger than, or comparable to those of other models. Globally, in any of these deviance comparisons, at least two tree-based methods out of RF, GBM and XGB tended to outperform models based on the linear regression principle. The only exception to the tree-based methods domination was under the product modelling strategy for Dataset 2 and for severity prediction in Dataset 1. Another key finding was that with respect to all deviance measures in both datasets, none of the ML models stood out as a single optimal technique.

Based on the training dataset, we calculated a calibration factor, so that the sum of the overall predicted premium equals the actual amount of claims paid for the portfolio. The derived value of the factor was then used to calibrate the test dataset. Based on 250 bootstrapped data runs of the test samples, we calculated a measure of the portfolio-level bias by comparing the sum of predicted premiums to the actual claims amount. The confidence intervals for bootstrapped bias calculated for each model and dataset are given in Table 3. All intervals include zero, indicating an absence of significant bias, which was expected given the models were calibrated. The variation in bias values is higher in Dataset 2 than in Dataset 1 due to the relatively smaller size of Dataset 2.

Table 3. Bootstrapped confidence intervals for model bias on the predicted premium values

3.4 Risk Stratification

Table 4 below provides the bootstrapped mean Gini indices (defined in Section 2.5.4) and corresponding standard errors obtained from the comparison of the machine learning (ML) methodologies with the reference GLM approach, on the test subset. Under the product modelling strategy, this comparison demonstrated that all tree-based models (RF, GBM and XGB) significantly improved policyholder risk profiling into high- and low-risk policies when compared to the GLM model ( $p \lt {10^{ - 6}}$ ).

Table 4. Bootstrapped Gini indices (and corresponding standard errors) for comparison of risk stratification from the ML pipelines against that of the reference GLM

For the product model, improvement appeared strongest from RF and XGB for Dataset 1 and from GBM for Dataset 2. This observation, similar to the results from the previous section, shows that these models performed relatively differently for the two datasets. We note that the regularised linear regression models (elastic net and LASSO) also yielded improvement over GLM for Dataset 1 ( $p \lt {10^{ - 6}}$ ), but not for Dataset 2 ( $p = 0.2725$ ).

For the Tweedie model, all tree-based models (RF, GBM and XGB) significantly improved policyholder risk profiling compared to the GLM ( $p \lt {10^{ - 6}}$ ). Improvement appeared strongest from RF and GBM for Dataset 1 and from GBM for Dataset 2, which, compared to the results from the previous section, shows that these models performed relatively similarly for the two datasets. We note that the regularised linear regression models (elastic net and LASSO) also yielded improvement over GLM for Dataset 1 ( $p \lt {10^{ - 6}}$ ), but not for Dataset 2 ( $p = 0.3164$ ). Globally, these results demonstrate that the ML methodologies led to a positive economic benefit (i.e. higher profit, assuming modelling of price elasticity capabilities are similar and ignoring customer loyalty), compared to the GLM pricing approach.

3.5 Loss Ratios and Conversion Rates in Open-Market Competition

Market competition is a key consideration in the pricing process. We build upon the work of Li et al. (Reference Li, Li and Shen2021) to develop a dynamic pricing game. The game compares pricing strategies based on LRs. Multiple insurers (models) compete to sell insurance contracts based on their predicted premiums. In this setup, we construct a zero-sum game entirely based on premiums. Figure 5 shows the bootstrap distributions of total (aggregated) actual premiums, total claims and LRs for each dataset in a direct two-player competition between the GLM, XGB and RF methodologies.

Figure 5. Bootstrap distributions of total book premiums (top), total claims ( ${2^{nd}}$ row) and loss ratios ( ${3^{rd}}$ row) comparison between two models, using, respectively, a GLM and an XGB methodology for premium pricing, on Dataset 1 (left) and Dataset 2 (right). The dashed line indicates a loss ratio of 1. The same analysis was carried out to compare RF with GLM, with loss ratio results (bottom) consistent with the comparison between XGB and GLM.

The LR charts in Figure 5 for GLM versus XGB (third row) and GLM versus RF (bottom) show a similar result when GLM is compared to these tree-based models. However, the difference in LRs varies in intensity when comparing GLM to XGB versus GLM to RF. This demonstrates that the RF and XGB models compete with GLM in different ways. The top and second row charts show that the comparisons between models yield different findings between Datasets 1 and 2.

For Dataset 1, under the product model strategy, when comparing GLM and XGB, the XGB tended to yield lower premiums for a larger book of business. The lower premium prediction resulted in XGB securing a higher amount of total premium, with the XGB median total premium sitting higher at around € 2,700k compared to the GLM median at around € 2,100k. However, the distributions of total claims payable were comparable for both models. This led to a lower ratio and more profitable business when using an XGB pricing model compared to a GLM model. This aligns with the results of the risk stratification analysis, where the XGB model was found to improve risk profiling over GLM. To make the comparison fair, models were calibrated to quote (on average) the same premium for the considered group of policies.

For Dataset 1 again, under the Tweedie model strategy, the trend was different, with GLM predicting lower individual premiums (not shown in the figure) compared to XGB, which resulted in the XGB securing a lower amount of business, i.e. a lower total book premium, as shown in Figure 5, top left. However, the business secured with XGB tended to attain a lower LR because of comparably lower overall claims payable. This result suggests that XGB yields a lower LR (risk-averse) pricing, which leads to a comparatively lower overall book premium while still retaining higher profits.

In Dataset 2, the distributions of LRs between GLM, XGB and RF are comparable. Hence, there is no clear indication that ML models outperform GLM in terms of LR. This can be driven by the relative difference in size and type of policies between Datasets 1 and 2. Here, the results from the product and Tweedie model strategies are consistent. Comparing GLM and XGB models in an open-market scenario also indicated improvement in the distribution of XGB LRs. However, there were differences in the type and amount of business the two models earned. For Dataset 2 (under product model), XGB earned a lower amount of total premium (with medians of € 830k and € 760k for GLM and XGB, respectively), but with significantly lower amounts of claims payable (with medians of € 900k and € 790k for GLM and XGB, respectively). Although XGB led to writing a lower number of policies compared to RF, these policies had a lower LR.

3.6 Overall Comparison

Figure 6 captures the essential metrics (namely bias, Poisson deviance, gamma deviance, Tweedie deviance, Gini index and conversion rate) in the form of radar charts to summarise the previous analyses and allow for overall model comparison on Datasets 1 and 2. For Dataset 2, additional metrics (calibration, actuarial group fairness and demographic parity) are shown on the radar chart that help measure model fairness. These metrics are described in more detail in Section 4 on group fairness.

Figure 6. Radar charts comparing the reference GLM with LASSO and the tree-based models on all key performance indicators.

In this analysis, the scales were calibrated independently for each dataset and key performance indicator. This calibration ensured that in each chart, the minimum and maximum values achieved among all models for a given metric were represented as the minimum and maximum achieved for that metric. A maximum value in any metric indicates comparably superior, i.e. preferable, model performance. Since the elastic net and LASSO yielded comparable performance, only the LASSO is represented in these summaries.

For Dataset 1, the overall footprints of the RF and XGB models appear larger than those of other models. The tree-based models led to the most competitive frequency, severity and Tweedie premium estimation performance for both datasets. Under the product model, for Dataset 1, XGB also tended to be the optimal model in terms of the Gini index and conversion rate. In contrast, GBM yielded the lowest bias on the overall portfolio level but appeared inferior for risk stratification (measured by Gini index). Under the Tweedie model, for Dataset 1, GBM was the optimal model in terms of the Gini index and conversion rate. In contrast, XGB yielded the lowest deviance on the overall portfolio-level prediction.

For Dataset 2, the GBM and XGB appeared as the two most competitive models for most metrics except for gamma deviance (product) and bias (Tweedie). Every time, GBM yielded the most effective risk stratification and conversion rate and, hence, was the most profitable model. (Recall that these relative differences were found to be statistically significant in previous sections.) Overall, GLM overall appeared weaker than all other methods for most criteria.

From the fairness aspect, GBM was the optimal model when calibration was used as the fairness measure. For demographic parity and actuarial group fairness RF (under the product model) and XGB/GBM (under the Tweedie model) were optimal. These measures are discussed in more detail in the next section.

4.1 Group Fairness, Demographic Parity and Calibration

Figures in each section below present outputs specific to the evaluation of the three axioms defined in Section 2.5.9.

4.1.1 Demographic Parity (Axiom 1)

A comparison of test set distributions of actual premium ${\hat P_A}$ with respect to protected variable gender ${Z_{gender}}$ showed that the distributions from the XGB model ( $\mathbb{P}({\hat P_{{A_{XGB}}}}|{Z_{gender}} = male)$ and $\mathbb{P}({\hat P_{{A_{XGB}}}}|{Z_{gender}} = female)$ ) were comparable for both product (Figure 7, left) and Tweedie (Figure 7, right) modelling strategies, unlike the corresponding GLM model. In addition, the Wasserstein distance ${d_{gender}}$ defined similar to equation (7) by

$${d_{gender}} = W\left( {{{\hat P}_A}\left| {{Z_{gender}} = male,{{\hat P}_A}} \right|{Z_{gender}} = female} \right)$$

calculated for the test set distribution were 20.0 and 19.1, respectively, for GLM and XGB under the product modelling, and 20.2 and 15.5 under the Tweedie modelling strategies. Table 5 gives the Wasserstein distances, the differences in distribution medians, and the p-values of two-sided Mann–Whitney (Wilcoxon) tests for all the models used in the study. Based on the test outcomes reported in Table 5, we can conclude that the actual premium distributions differed significantly between males and females.

Figure 7. Distributions of the GLM and XGB model actual premium, by gender, under the Product (left) and Tweedie (right) model strategies.

Table 5. Measure of the distance ( ${d_{gender}}$ ), difference in distribution medians and Wilcoxon test p-values of the actual premium models under the product model (left) and Tweedie model (right) strategies

4.1.2 Actuarial Fairness (Axiom 2)

A comparison of test data distributions of actual premiums ( $\mathbb{P}({\hat P_A}|{\hat P_T},{Z_{gender}})$ ) conditional on predicted TP bands (defined by TP quintiles as in equation (8)) and on gender is provided in Figure 8. The results showed greater alignment between distributions of XGB-derived ratios for male, female and overall drivers, i.e. greater alignment with Axiom 2, compared to those from GLM. Discrepancies ${d_{l,{\cal Y}}}$ defined by equation (9) were calculated for the five TP bands successively, as

$${d_{l,{P_T}}} = W\left( {{{\hat P}_A}\left| {{{\hat P}_T},{Z_{gender}} = male,{{\hat P}_A}} \right|{{\hat P}_T},{Z_{gender}} = female} \right)$$

Figure 8. Distributions of actual premium across technical premium quintile bands for the XGB model (left) and the GLM model (right) under the product (top) and Tweedie (bottom) modelling strategies.

The Wasserstein distances for the GLM and XGB models for the product and Tweedie modelling strategies are given in Tables 6 and 7, respectively. This quantitation also indicates that the distance between the XGB distributions of both groups were closer compared to the same distributions from the GLM model.

Table 6. For the Product model, Wasserstein distances, median difference and Wilcoxon test p-values comparing actual premium distributions by technical premium band for actuarial group fairness (left) and calibration (right), as well as mean APTP ratios and associated standard deviation by actual premium band (centre). Significant p-values are italicised

Table 7. For the Tweedie model, Wasserstein distances, median difference and Wilcoxon test p-values comparing actual premium distributions by technical premium band for actuarial group fairness (left) and calibration (right), as well as mean APTP ratios and corresponding standard deviation by actual premium band (centre). Significant p-values are italicised

4.1.3 Calibration (Axiom 3)

A comparison of actual claim distributions conditional on actual premium bands ( $\mathbb{P}(L|{\hat P_A},{Z_{gender}})$ ) with respect to the protected variable gender is provided in Figure 9. Figure 9 does not demonstrate overall large discrepancies across gender groups from GLM-based pricing, for the majority of the data. Quantitatively, ${d_{L,Z}}$ defined by equation (10) were evaluated for each actual premium band defined as in equation (8), but using the AP, as follows:

$${d_{L,{P_A}}} = W\left( {L\left| {{{\hat P}_A},{Z_{gender}} = male,L} \right|{{\hat P}_A},{Z_{gender}} = female} \right)$$

Figure 9. Distributions of actual claims by actual premium quintile band for the XGB and GLM models under the Product and Tweedie modelling strategies.

These discrepancies between male and female groups for XGB and GLM actual claims by actual premium bands for the product and Tweedie modelling strategies are also given in Tables 6 and 7, respectively.