2.1. Parametric quantile regression ratemaking (PQRR)

Baione and Biancalana (Reference Baione and Biancalana2021) apply in the second stage the PQR to obtain a parsimonious model for the individual premium calculation. With PQR, the regression coefficients are modeled as parametric functions of the order of the quantile, depending on a set of parameters $\gamma$

(2.5) \begin{equation}\beta_{l}(\theta,\gamma_{l})=\gamma_{l0} + \gamma_{l1} b_{1}(\theta) + \dots + \gamma_{ld} b_{d}(\theta).\end{equation}

By convention, $b_{0}(\theta)=1$ and $l = 0,\dots,s$ . The quantity s represents the total count of numerical variables and all possible levels of the categorical risk factors, excluding the corresponding baseline level for each factor. On the other hand, d denotes the number of functions utilized for performing the PQR analysis. To choose the functions $b_1(\theta),\dots,b_d(\theta)$ , the only assumption to take into account is that their computation must lead to a $\theta$ -monotonically increasing quantile function. Baione and Biancalana (Reference Baione and Biancalana2021) suggest to use the shifted Legendre polynomials (SLPs). A polynomial of degree d is defined as $\widetilde{P}_{d}(x) = P_{d}(2x -1)$ with $\widetilde{P}_{d}$ orthogonal on [0,1]. The R package iqr, used to estimate the PQR, is optimized for this kind of orthogonal polynomials. The choice of the polynomial degree d shall be made in accordance with the dataset. In fact, QR can be affected by quantile crossing problems. This problem appears when quantile curves cross each other and this leads to a nonmonotonic behavior. In the context of individual premium estimation, this issue can be overcome by lowering the degree of the SLP. In fact, even if using higher degrees could appear better since leads to a more flexible model, a quantile crossing issue can arise due to the oscillations of polynomials. For a complete description of this method, we refer to Frumento and Bottai (Reference Frumento and Bottai2016) and Baione and Biancalana (Reference Baione and Biancalana2021). The advantage of this approach is that by calibrating only one regression we are able to model the conditional claims amount at the proper probability level in Equation (2.3) for each policyholder in the portfolio. To compute the insurance premium, Heras et al. (Reference Heras, Moreno and Vilar-Zanón2018) and Baione and Biancalana (Reference Baione and Biancalana2021) use a quantile premium principle that is a linear convex combination of the expected value of $Y_{i}$ and its $\theta$ -quantile

(2.6) \begin{equation} P^{PQRR}_{i} = (1-\alpha) \mathbb{E}[Y_{i}] + \alpha Q_{\theta}[Y_{i}],\end{equation}

for $\alpha\in [0,1]$ , where we recall that $Y_i$ is a function of the risk factors $\mathbf{x}_i$ . The rationale of this premium principle is to charge the fair value of a contract (the expected loss) with a safety loading component, which makes the premium more expensive than the mean. The unconditional expectation $\mathbb{E}[Y_{i}]$ can be computed using the law of total expectation:

(2.7) \begin{equation} \begin{split} \mathbb{E}[Y_{i}] & = \mathbb{E}\left[\mathbb{E}[Y_{i}|N_{i}]\right] \\ & = \mathbb{P}\! \left( N_{i}=0 \right) \mathbb{E}[Y_{i}|N_{i}=0] + \mathbb{P}\! \left( N_{i}>0 \right) \mathbb{E}[Y_{i}|N_{i}>0] \\ & = (1 - p_{i}) \mathbb{E}[Y_i\mid Y_i>0], \end{split} \end{equation}

and the individual premium $P^{PQRR}_{i}$ is then found with the PQR and the premium principle in Equation (2.6). The corresponding individual premium is denoted with PQRR. For its computation, it is necessary to calibrate the proper level for the loading factor $\alpha$ . The idea is that the sum of all premiums has to equate to the aggregate premium $\pi(Y)$ , where Y denotes the random aggregate loss. Therefore, we need to find the value of $\alpha$ that satisfies the following balance equation:

(2.8) \begin{equation}\pi(Y) = \sum_{i=1}^{n}[\alpha Q_{\theta}[Y_{i}] + (1 - \alpha) \mathbb{E}[Y_{i}]],\end{equation}

which yields

(2.9) \begin{equation}\alpha = \frac{\pi(Y) - \mathbb{E}[Y]}{Q_{\theta}[Y] - \mathbb{E}[Y]}.\end{equation}

Following Baione and Biancalana (Reference Baione and Biancalana2019, Reference Baione and Biancalana2021), the aggregate premium $\pi(Y)$ has to be computed as the sum of the expected aggregate loss plus a risk margin. The risk margin is selected using the solvency conditions:

(2.10) \begin{equation}\pi(Y) = (1 + 3\sigma )\mathbb{E}[Y].\end{equation}

The volatility factor $\sigma$ is set at the level of the 10% from Solvency II Standard Formula for Premium Risk in the Motor Third Party Liability insurance (EIOPA, 2009). In this way, we have indirectly found the coefficient $\alpha$ in a way that is compatible with the Solvency II regulation. We remark that $\alpha$ can be specified directly by the analyst in Equation (2.6). However, Baione and Biancalana (Reference Baione and Biancalana2021) consider this choice too important to be left at a discretionary judgment.

2.2. Quantile regression premium calculation (QRPC)

In the actuarial literature, Baione and Biancalana (Reference Baione and Biancalana2019) propose another method to estimate individual premiums using the QR.

In the second stage of the Two-Part model, a unique probability level is fixed for the conditional total claim amount ${Y}_i$ given that $Y_i > 0$ , in order to have $\theta^{*}_{i} \equiv \theta^{*}$ . As a consequence, the probability level for the total claim amount is no longer unique, but it is different for each policyholder, depending on $p_{i}$ and $\theta^{*}$ . So, it is possible to run a unique QR using the probability level $\theta^{*}$ for ${Y}_i$ whenever it is positive. With this intuition, they propose the following premium principle:

(2.11) \begin{equation}P^{QRPC}_{i}=(1-p_{i}) Q_{\theta^{*}}[{Y}_i\mid Y_i > 0].\end{equation}

To calibrate the unique probability level $\theta^*$ for all policyholders, Baione and Biancalana (Reference Baione and Biancalana2019) suggest finding the value of $\theta$ such that all individual premiums sum up to the aggregate premium $\pi(Y)$

(2.12) \begin{equation}\pi(Y)=\sum_{i=1}^{n}P^{QRPC}_{i}=\sum_{i=1}^n (1-p_{i}) Q_{\theta}[{Y}_i\mid Y_i > 0],\end{equation}

where $\pi(Y)$ is the aggregate premium calculated via Equation (2.10). In particular, $\theta^*$ can be found by optimizing the value of $\theta$ in Equation (2.12).

Using the QRPC and PQRR models, we find two distinct sets of individual insurance premiums since we use different premium principles and different methods to assess the probability level of interest. In our application in Section 5 we show that, despite differences in premiums at the individual level, at the global level the importance risk factors that influence premiums are the same.

3.1. Shapley value

The Shapley value (Shapley, Reference Shapley, Kuhn and Tucker1953) is an attribution method originated from game theory. It is used to allocate the total value generated by a coalition of players into the contribution of each individual player in a coalition.

Denote with $\phi_j$ the Shapley value of the j-th player. Consider a set of players $K = \{1,2, \dots,k\}$ that participate in a game. We denote with $u\subseteq K$ any possible subcoalition of players and with $\text{val}\,:\,2^{K} \to \mathbb{R}$ the value function of the game, with $\text{val}(\emptyset) = 0$ . The value generated by the players in the subset u is $\text{val}(u)$ and the total value of the game to be allocated is $\text{val}(K)$ .

Consider the following desirable properties for an attribution method (Winter, Reference Winter, Aumann and Hart2002):

• • Efficiency: The sum of Shapley values of all players equates to the total value of the game $\sum_{j=1}^{k}\phi_{j} = \text{val}(K)$ . This means that the total payoff of the game is allocated among all players.

• • Symmetry: If players $j_1$ and $j_2$ make the same marginal contribution to any coalition, $\text{val}(u \cup \{j_1\}) = \text{val}(u \cup \{j_2\})$ for all $u \subseteq K$ with $j_1,j_2 \notin K$ , then $\phi_{j_1}=\phi_{j_2}$ . Thanks to this axiom, the order in which players are considered does not influence the value of the coalition; players with the same marginal contribution receive the same payoff.

• • Dummy: If j is a dummy player, $\text{val}(u \cup \{j\}) = \text{val}(u)$ for all $u \subseteq K$ , then $\phi_{j}=0$ . In other words, if the player’s contribution to the coalition is zero, the payoff of the player is zero (axiom of the null-player).

• • Additivity: If val and $\text{val}'$ have Shapley values $\phi$ and $\phi'$ , respectively, then the game with value $\text{val}(u) + \text{val}'(u)$ has Shapley value $\phi_{j}+ \phi'_{j}$ for all $j \in K$ . This means that if two games are combined, the payoff assigned to the j-th player is equal to the sum of the two payoffs he would get playing the two games separately.

For the generic player j, Shapley (Reference Shapley, Kuhn and Tucker1953) proves that the unique value satisfying these four properties can be expressed as:

(3.5) \begin{equation}\phi_{j} = \sum_{u \subseteq K\setminus\{j\}}\dfrac{(k - |u| -1)!|u|!}{k!}(\text{val}(u \cup \{j\}) - \text{val}(u))\end{equation}

This index is based on the marginal increase $\text{val}(u \cup \{j\}) - \text{val}(u)$ generated when the j-th player joins the coalition u. The marginal contribution is then averaged over all possible coalitions.

The Shapley value is attracting increasing attention in the sensitivity analysis literature as a measure of variable importance (Owen, Reference Owen2014; Song et al., Reference Song, Nelson and Staum2016). Owen (Reference Owen2014) suggests to use as value function the Sobol’ index

(3.6) \begin{equation}\text{val}(u)=\underline{\tau}^{2}_{u}.\end{equation}

The Shapley values obtained with this value function are called the Shapley effects (Song et al., Reference Song, Nelson and Staum2016). Shapley effects enjoy appealing properties:

1. 1. $\phi_j\geq 0$ for all $j=1,2,\dots,k$ . This holds by Equation (3.5) and because $\underline{\tau}^{2}_{u\cup \{j\}}-\underline{\tau}^{2}_{u}$ is positive as the variance explained increases when a new risk factor is considered for all $u\subseteq K$ .

2. 2. $\sum_{j=1}^k \phi_j=\textrm{Var}(W)$ . This follows from the efficiency property, meaning that we allocate the variance of the model output.

Remarkably, Properties 1 and 2 hold regardless the dependence structure and marginal distributions of the risk factors. This feature contributes significantly to the widespread adoption of Shapley effects as sensitivity measures in the interpretation of complex models.Footnote 1

The computation of the Shapley effects is not a trivial issue. Song et al. (Reference Song, Nelson and Staum2016) propose a double loop Monte Carlo for the estimation of the value function to be inserted in the Shapley representation in Equation (3.5). Plischke et al. (Reference Plischke, Rabitti and Borgonovo2021) propose an algorithm based on the concept of Möbius inverses, which reduce the number of value functions to be evaluated with respect to the approach of Song et al. (Reference Song, Nelson and Staum2016). With this algorithm, the number of value functions to be evaluated is remarkably reduced from $k!\cdot k$ to $2^k-1$ . Furthermore, it should be noted that the Plischke algorithm offers a way to estimate Shapley values for interactions. However, since it goes beyond the scope of this work, a thorough exploration of this topic is not included here. For more detailed insights, interested readers can refer to the works of Plischke et al. (Reference Plischke, Rabitti and Borgonovo2021) and Lundberg et al. (Reference Lundberg, Erion, Chen, DeGrave, Prutkin, Nair, Katz, Himmelfarb, Bansal and Lee2020).

Denote with $\text{mob}(u)$ the Möbius inverses of the value function $\text{val}(u)$ with:

(3.7) \begin{equation} \text{val}(u) = \sum_{v \subseteq u}\text{mob}(v) \quad \text{and} \quad \text{mob}(u) = \sum_{v \subseteq u}({-}1)^{|u|+|v|} \text{val}(v).\end{equation}

Note that, in the case of independent risk factors, we find $\textrm{mob}(u)=\sigma_u^2$ . Equivalently to the representation in Equation (3.5), Shapley values can be written as

(3.8) \begin{equation}\phi_{j}=\sum_{u:j \in u}\dfrac{\text{mob}(u)}{|u|} .\end{equation}

Following Plischke et al. (Reference Plischke, Rabitti and Borgonovo2021), let all $2^{k}-1$ non-vanishing subsets be indexed by $u_{m}$ with $m=1,\dots,2^{k}-1$ . The value functions are included in a vector H with coordinates $H_{m}=\text{val}(u_{m})$ . To obtain subset inclusion, a matrix Z is proposed:

(3.9) \begin{equation}Z_{mf}=\begin{cases} 1 & \text{if $u_{m} \subseteq u_{f}$} \\[5pt] 0 & \text{otherwise} .\end{cases}\end{equation}

Then, Möbius inverses are obtained with $M=HZ^{-1}$ and Shapley effects are given by

(3.10) \begin{equation}\phi_{j}=\sum_{m:j\in u_{m}}\frac{M_{m}}{|u_{m}|} .\end{equation}

Now, the Shapley effects computation requires an estimate of the value function in Equation (3.6) to compute the Möbius inverses. Differently from the computer experiment literature, we need a method to estimate the Sobol’ indices from given data. A solution to this issue comes from the recent work of Mase et al. (Reference Mase, Owen and Seiler2019), who propose the estimation of $\text{Var}(\mathbb{E}[W|\mathbf{X}_{u}])$ using a cohort approach, as presented in the next section.

3.2. The cohort approach

Assume we have available at hand a dataset of observations ${\{(w_i,\mathbf{x}_i)\}}_{i=1}^n$ , where $w_i$ denotes the quantity of interest. The cohorts are subsets of observations sharing the values of risk factors similar to the ones of a target observation denoted by t. Let $x_{ij}$ be the value of the j-th risk factor for the i-th observation (with $j=1, \dots,k$ and $i=1, \dots, n$ ) and let $x_{tj}$ be the target value for the selected risk factor j. It is necessary to construct a similarity function $z_{tj}$ for each risk factor j. In particular, $z_{tj}:X_{j} \to \{0,1\}$ . The similarity function is equal to 1 when the subject i is considered similar to the target subject t for the risk factor j. In the case of real-valued covariates, it can be appropriate to define similarity using a measure of relative distance. A possible similarity function is

(3.11) \begin{equation} z_{tj}(x_{ij})= \begin{cases} 1 & |x_{ij} - x_{tj}| \le \delta_{j}, \\[5pt] 0 & \text{otherwise}, \end{cases}\end{equation}

where $\delta_{j}$ is the specified distance. If the variables are continuous, Mase et al. (Reference Mase, Owen and Seiler2019) consider $\delta_j=0.1\cdot\text{range}(x)$ , while $\delta_{j}=0$ if they are discrete. Now it is possible to define the cohort of t with respect to risk factors in u as

(3.12) \begin{equation}C_{t,u}=\{i \in 1\,:\,n \quad | \quad z_{tj}(x_{ij})=1, \forall j \in u\}.\end{equation}

By definition, $C_{t,u}$ is the set of subjects that are similar to the target subject t for all risk factors $j \in u$ . By construction, $C_{t,u}$ is never empty since it contains at least the target observation $x_t$ . Using the cohorts, Mase et al. (Reference Mase, Owen and Seiler2019) found estimates of the Sobol’ indices as follows. Consider the cohort means

(3.13) \begin{equation} \frac{1}{|C_{t,u}|} \sum_{i \in C_{t,u}}w_i = \bar{w}_{t,u},\end{equation}

where $|C_{t,u}|$ is the cardinality. It holds that the cohort average $\bar{w}_{t,u}$ is an estimate of the conditional mean for the target $\mathbb{E}[W| X_{t,u}]$ . Since each observation has probability $n^{-1}$ , we obtain an estimate of the value function:

(3.14) \begin{equation}\text{val}(u)= \frac{1}{n} \sum_{t=1}^{n}{\left( \bar{w}_{t,u} - \bar{w}\right)}^{2},\end{equation}

where $\bar{w}$ is an estimate of the model output mean $\mathbb{E}[W]$ . Hence, to compute the Shapley effects, we use the value function in Equation (3.14), estimated via cohort approach, to calculate the Möbius inverses in Equation (3.7). Finally, Shapley effects are obtained via Equation (3.8), adopting the algorithm proposed by Plischke et al. (Reference Plischke, Rabitti and Borgonovo2021). A MATLAB implementation can be found at https://github.com/giovanni-rabitti/ratingsystemshapley.

4.1. Analytical example

We illustrate our procedure with an analytical example. Assume that the generic premium $P^*$ can be written as

(4.2) \begin{equation} P^*=\sum_{j=1}^{k} \lambda_{j} X_{j}\end{equation}

where each variable $X_{j}$ is considered to be discrete (continuous variables are typically discretized for rating systems construction – see Henckaerts et al., for $j=1, \ldots, k$ . Every realization of $\mathbf{X}$ is $\mathbf{x,}$ which is the vector of the characteristics of a policyholder (for the sake of simplicity, we drop the index i in this analytical example).

Assume that the rating system is constructed using the selected variables $u\subseteq K$ . In this system, a generic class is indexed by r with $r=1, \dots, R$ and is uniquely identified by fixing the levels of the vector $\mathbf{X}_u=\mathbf{x}_u$ . In this class, a generic premium is given by

(4.3) \begin{equation} P_r = \sum_{j\in u} \lambda_{j} x_j + \sum_{j \notin u} \lambda_{j} X_{j}.\end{equation}

In other words, since the other variables $\mathbf{X}_{-u}$ can vary, $P_r $ is a random variable representing the individual policyholders’ premiums in the class r. In particular, when $u=\{1,\dots,k\}$ , then r is the cohort of policyholders sharing the same levels for all variables and $P_r$ becomes deterministic. In the present framework, the following result holds.

Proof. First of all, by the additivity and efficiency axioms of the Shapley values and the independence of the $X_{i}$ s, the variance of $P^*$ is

(4.4) \begin{equation} \text{Var}(P^*)=\sum_{j=1}^{k}\phi_j,\end{equation}

where $\phi_j=\lambda_{j}^2 \text{Var}(X_{j})$ for all $j=1,\dots,k$ . Analogously, the variance of the individual premiums of policyholders in the class r is

(4.5) \begin{equation} \text{Var}\left(P_r\right)=\sum_{j\not\in u}\phi_{j}\end{equation}

because of Equation (4.3). Moreover, the average class premium in r is

(4.6) \begin{equation} \mathbb{E}\left[P_r\right]=\sum_{j\in u} \lambda_{j} x_j + \sum_{j \notin u} \lambda_{j} \mathbb{E}\left[X_{j}\right].\end{equation}

This mean value represents the premium that all policyholders in class r are required to pay. Then, it holds

(4.7) \begin{equation} \begin{split} \mathbb{P}\! \left( \left|P_r-\mathbb{E}\left[P_r\right]\right| \geq \xi \right) &= \mathbb{P}\! \left( \left|\sum_{j \notin u} \lambda_{j} X_j - \sum_{j \notin u} \lambda_{j}\mathbb{E}\left[X_{j}\right]\right|\geq \xi \right)\\ & \leq\frac 1{\xi^2}{\text{Var}\left(\sum_{j \notin u} \lambda_{j} X_j\right)}=\frac 1{\xi^2}\sum_{j \notin u}\phi_j \end{split}\end{equation}

by Chebyshev’s inequality applied to the random variable $\sum_{j \notin u} \lambda_{j} X_j$ , which represents the residual variability within the system constructed by fixing $\mathbf{X}_u$ . If $\mathbf{X}_u$ are the rating factors with the highest Shapley effects, then the upper bound of $\mathbb{P}\! \left( \left|P_r-\mathbb{E}\left[P_r\right]\right| \geq \xi \right)$ is the lowest possible, because for any fixed $\xi$ the value $\sum_{j \notin u}\phi_j$ is minimized. Hence, for any class in this rating system, the probability that individual premiums deviate from the class mean is minimized.

It is worth commenting Proposition 4.1. First of all, Equation (4.5) is a within sum of squares index for the individual premiums in the class r (since we are summing up the residual variances of the individual premiums in the class r). However, differently from clustering methods, our aim is to minimize this index by changing the selected rating factors and not by aggregating factor levels. Secondly, we assumed a linear pricing model with independent and discrete inputs. This choice was necessary to obtain a closed form expression of the Shapley effects (as a matter of fact, Owen and Prieur, Reference Owen and Prieur2017 managed to provide analytical expressions of Shapley effects only for some specific bivariate cases).

4.2. Homogeneity indices for system comparison

In general, numerical methods represent the only possible approach to estimate these sensitivity indices. Consequently, it becomes necessary to include statistical indicators to assess the reliability of the our procedure. Consider the premium $P^*_{i}$ , which is the loading premium computed at the individual level for each policyholder with QRPC or PQRR presented in Equations (2.6) and (2.11). To measure the difference of the individual premiums when computed with all risk factors and with the selected risk factors with the Shapley effects, we propose the mean-squared difference (MSD) defined as

(4.8) \begin{equation} \text{MSD} = \frac{1}{n}\sum_{i=1}^n {\left( P_i^* - P_i^{\rm{red}} \right)}^2,\end{equation}

where $P_i^{\rm{red}}$ represents the class premium of the i-th policyholder computed with the reduced model. The lower the MSD, the closer the individual premium is to the class premium in the rating system. We observe that $P_i^{\rm{red}}$ coincides with the class premium to which the i-th policyholder belongs to (in every risk class all policyholders pay the same average price), so it is the empirical counterpart of Equation (4.6). To emphasize the dependency on the class r rather than on the policyholder i, we can denote the same average premium as $P_r^{\rm{red}}$ . Hence, we can equivalently rewrite the MSD as

(4.9) \begin{equation} \text{MSD} = \frac 1n \sum_{r=1}^R \sum_{{i_r}=1}^{n_r} {\left(P^*_{i_r} - P_r^{\rm{red}}\right)}^2,\end{equation}

where $n_r$ is the number of insured in the risk class r, with $\sum_{r=1}^R n_r=n$ . Equation (4.9) tells us that the MSD is the sum of the deviations of the individual premiums with respect to their class means, averaged for all policyholders in all classes. Hence, the numerator of the MSD is a sum of squared distances, which can be interpreted as a sum of within-classes variations.

We can also finally consider two measures of deviation used in the actuarial literature (Henckaerts et al., Reference Henckaerts, Antonio, Clijsters and Verbelen2018), namely the goodness of variance fit (GVF)

(4.10) \begin{equation} \text{GVF} = 1- \frac{\sum_{r=1}^R \sum_{i_r=1}^{n_r} {\left(P^*_{i_r} - P_r^{\rm{red}}\right)}^2}{\sum_{i=1}^n {\left( P_i^* - \overline{P} \right)}^2},\end{equation}

where $\overline{P} = 1/n\sum_{i=1}^{n}P_i^*$ and the tabular accuracy index (TAI)

(4.11) \begin{equation} \text{TAI} = 1- \frac{\sum_{r=1}^R \sum_{i_r=1}^{n_r} |P^*_{i_r} - P_r^{\rm{red}}|}{\sum_{i=1}^n | P_i^* - \overline{P} |}.\end{equation}

The denominator in the previous formulas measures the deviation of each individual premium from the global mean. The numerator measures the deviation of each individual premium from its class mean. The closer the GVF and TAI values are to 1, the more uniform the classes become in terms of their homogeneity. Lastly, note it holds

(4.12) \begin{equation} \text{GVF} = 1- \frac{n}{n-1}\frac{\text{MSD}}{\sum_{j=1}^k\phi_j},\end{equation}

formally connecting the GVF with the MSD and the Shapley effects. We remark that the GVF provides insights on how close the rating system is with respect to the lowest heterogeneity case.

In practice, the rating classes have different means and heterogeneous variances and they are unbalanced. Thus, the problem of their comparison in terms of internal variability intensifies. In general, when the means of classes differ significantly, comparing the absolute differences in variability (i.e. TAI) or the variances (i.e. MSD) can be misleading. For example, in Equation (4.9) the MSD does not explicitly consider the relative variations within each class. When comparing rating systems, it is important to account for this effect in order to ensure an accurate and effective comparison. This is because the absolute differences may be influenced more by the class means rather than the inherent variability within each class. However, this is not the case when considering measures of relative variability, which account for these differences in means. One of such measures is the coefficient of variation (Olivieri and Pitacco, Reference Olivieri and Pitacco2015)

(4.13) \begin{equation}\text{CV}=\dfrac{\sqrt{\sigma^{2}_{P}}}{\overline{P}}\end{equation}

where $\sigma_{P}^{2}$ is the variance of the premiums. Olivieri and Pitacco (Reference Olivieri and Pitacco2015) state that the CV is suitable to compare the riskiness of classes composed of a different number of policyholders. Mahmoudvand and Oliveira (Reference Mahmoudvand and Oliveira2018) use the CV to measure concentration risk in a loan portfolio, while in the actuarial literature Donnelly (Reference Donnelly2014) adopts it to quantify the mortality risk in a defined-benefit pension scheme.

When comparing rating systems in terms of overall homogeneity of their risk classes, we can adopt the coefficients of variation for every risk class and then average. However, a simple arithmetic average of all coefficients of variations of the risk classes of a rating system might not be a good indicator of the general homogeneity of that system. Namely, a risk class with fewer insured that might be more homogeneous will receive the same weight as a risk class with more insured, producing a distortion. Since different rating systems are based on different combinations of the selected risk factors, they are based on different risk classes. We want to aggregate the coefficient of variations of the risk classes of a rating system in order to make comparisons between rating systems. We consider a linear combination of the coefficients of variation. Precisely, an arithmetic mean would be a misleading indicator of homogeneity since classes with different sizes would receive the same weight. Conversely, adopting weighted means can give the right weight to the homogeneity of risk classes composed of a large number of insured. Hence, we propose to compare rating systems using the weighted mean of coefficient of variations (WMCV)

(4.14) \begin{equation} \text{WMCV}=\dfrac{\sum_{r=1}^{R} n_{r}\cdot \text{CV}_{r}}{\sum_{r=1}^{R}n_{r}}.\end{equation}

With this approach, we are able to firstly compute the homogeneity of each risk class (with the coefficient of variation) and, secondly, we take the weighted mean to have a measure for the entire rating system. We consider the WMCV the most appropriate index to compare rating systems in terms of overall homogeneity.

5.1. Individual premium estimation

The first step of the Two-Part model consists of the estimation of the probability of having no claims. In this case, we use a GLM regression of the binomial family with a logit link function adjusted for the exposure in Equation (2.1). The R code for the logit exposure adjusted is provided by de Jong and Heller (Reference de Jong and Heller2008). In Table A.1, we report the parameters of the logit GLM regression obtained via the maximum likelihood estimation.

In the second step, we estimate the severity component. We first perform the PQR to find the premiums with the PQRR approach (see Section 2.1). The selected function $b(\theta)$ to implement the procedure is a SLP of degree 1. In particular, the regression coefficient has the following formulation: $\beta_{l}(\theta, \gamma_{l})= \gamma_{l0} + \gamma_{l1} (2\theta - 1)$ , $l =0, \dots, 27$ .Footnote 2

Baione and Biancalana (Reference Baione and Biancalana2021) use in the same dataset a SLP of degree 3, since in their numerical application the purpose is the estimation of insurance premiums for the eight risk classes that compose the rating system they construct. In our case, we need to use a SLP of degree 1 since a degree equal to 3 produces numerical instability. In particular, using a polynomial of the second or third degree, we obtain some exceptional small/large values (outliers) for the severity component. Figure A.2 in the Appendix A shows that the goodness-of-fit of this model is satisfactory despite the low degree. Estimated parameters of the PQR are reported in Table A.4 in the Appendix A. In Table 2, the Wald test for the significance of rating factors is shown. From Table 2, we can see that the most statistically significant risk factors are $\textit{vehicle body}$ , $\textit{agecat}$ , and $\textit{area}$ . However, as previously explained, a selection based on the p-values might be misleading due to the lack of interpretability of the ranking in terms of importance of the risk factors. To solve this critical issue, we use the Shapley effects approach to be able to obtain factor prioritization. The Wald test for the $b(\theta)$ components is also provided in Table A.5 in the Appendix. Using Equation (2.6), we estimate the quantile premium principle. As previously done in Heras et al. (Reference Heras, Moreno and Vilar-Zanón2018) and Baione and Biancalana (Reference Baione and Biancalana2021) on this dataset, we use a Gamma GLM with log link function to estimate the conditional expectation of the loss distribution (regression coefficients are provided in Table A.3 in the Appendix).

Table 2. Wald test for rating factors.

Finally, we need to calibrate the loading factor $\alpha$ using Equation (2.9). The value obtained through the calibration procedure is $\alpha = 6.84\%$ . The severity component can also be estimated using QR and the individual premiums can be found by applying the QRPC approach (explained in Section 2.2): this method requires to find the probability level $\theta^{*}$ such that the balance Equation (2.12) holds. In our application, the resulting value is $\theta^{*} = 78.60\%$ . We recall that we are able to run a unique QR with this probability level. The obtained parameters estimates are shown in Table A.2 (see the Appendix). At this point, we can compute the individual insurance premium using the quantile premium principle in Equation (2.11).

5.2. Application of Shapley effects to the rating system construction

Once individual premiums are estimated, we apply the Shapley effects algorithm to the datasets $\left\{\left(P^{QRPC}_i,\mathbf{x}_i\right)\right\}_{i=1}^n $ and $\left\{(P^{PQRR}_i,\mathbf{x}_i)\right\} _{i=1}^n $ to identify the two risk factors that, in combination, explain a greater amount of variance compared to any other pair of variables. The utilization of Shapley effects is further justified by the dependence analysis depicted in Figure A.1 in the Appendix. This analysis reveals that the risk factors are not independent, strengthening the rationale for employing the Shapley effects algorithm. Given the presence of six discrete risk factors, it is important to define the notion of similarity between two observations within a subset u. Specifically, two observations are considered similar with respect to subset u if they share the same values for the coordinates indexed by u. In the case of discrete risk factors, the similarity function in Equation (3.11) is set to $\delta_{\text{j}}=0$ (as done in Mase et al., Reference Mase, Owen and Seiler2019). By applying the algorithm to the premiums obtained with QRPC and PQRR, we display the results in Figure 1. Risk factors $\textit{agecat}$ and $\textit{vehicle body}$ explain together $82.45\%$ of individual premium total variance in case of the QRPC approach and $77.77\%$ in case of the PQRR approach. These two risk factors are the greatest drivers of the individual premiums. Conversely, vehicle value is negligible and could be “frozen” to a constant value, reducing the model complexity (see, for instance, Saltelli et al., Reference Saltelli, Marco, Terry, Campolongo, Cariboni, Gatelli, Saisana and Tarantola2008 for a discussion on freezing irrelevant factors). Also, gender and vehicle age have a low impact on the variability of the individual premiums. These insights enable model interpretability since the analyst knows the main drivers of the pricing models. We note that, despite some differences in premiums at the individual level between the two pricing methods, at the global level the importance of risk factors is the same for both QRPC and PQRR premiums, as shown in Figure 1. Only the fourth and fifth most important risk factors are switched between QRPC and PQRR: vehicle age and gender are, respectively, the fourth and fifth most important risk factors for QRPC, while their roles are inverted for the PQRR. However, the magnitude of the fourth highest Shapley effect is relatively small (around 3–4% for the two pricing methods). We lastly note that discretizing the variable veh_value does not change its irrelevant importance. If it is left continuous, its Shapley effects are $0.0014$ for the QRPC and $0.0047$ for the PQRR, while as a discrete variable they are $0.0036$ and $0.0150$ , respectively.

Figure 1. Shapley effects for the individual premiums computed with QRPC and PQRR. The y-axis shows the Shapley effects for each risk factor (displayed on the x-axis). By the efficiency property, the Shapley effect normalized by the variance can be interpreted as the proportion of explained variance by a single risk factor.

5.3. Rating systems comparison

The estimation of Shapley effects provides us insights about the most important risk factors from a statistical point of view. On the other hand, we need an actuarial valuation to demonstrate that the rating system we have constructed is the one with less variability in the insurance premium. We adopt the indices presented in Subsection 4.2 as a measure to assess the variability within each risk class in a rating system and to evaluate the riskiness of the overall rating system.

With our dataset, composed of 6 risk factors, we can construct 15 different systems made up of 2 risk factors each (remember that we refer to the discretized variable $\textit{vehicle value}$ ). In the literature, two examples of rating systems constructed with this dataset are provided. Heras et al. (Reference Heras, Moreno and Vilar-Zanón2018) build a system using the risk factors $\textit{vehicle age}$ and $\textit{agecat}$ selected with the p-values resulting from the logistic regression model. However, $\textit{vehicle age}$ is not an important risk factor as Figure 1 shows. The second application is offered in Baione and Biancalana (Reference Baione and Biancalana2019) and (2021): the rating system is constructed by using $\textit{gender}$ and $\textit{vehicle age}$ to simplify the representation of their models.

For every rating system based on two risk factors, we estimate the weighted mean of the variation coefficient via Equation (4.14) and the number of non-empty risk classes. Furthermore, we compute the MSD between the individual premiums using the quantile models with all risk factors included and quantile models based on two risk factors only. Note that, for the latter case of two factors, the PQRR and the QRPC for each rating system were recalibrated.Footnote 3

The results, as presented in Table 3, demonstrate that the rating system built on the factors of $\textit{vehicle body}$ and $\textit{agecat}$ exhibits the lowest overall variability in average premiums within each risk class. Additionally, from Table 3 we note the following remarkable aspects. To begin with, all rating systems constructed using the risk factor $\textit{agecat}$ exhibit low values of WMCV. This was expected since $\textit{agecat}$ explains the largest fraction of the variance of the premiums. The MSD is also minimized for both models, when the risk factors are selected using Shapley effects. Consequently, the GVF is maximized as per Equation (4.12) note that it takes values around 0.8, which are close to 1 (recall that a GVF close 1 indicates a high degree of homogeneity). Furthermore, it is important to point out that the low WMCV and the high GVF are not a consequence of the high number of classes. In fact, the selected rating system with $\textit{vehicle body}$ and $\textit{agecat}$ has 78 risk classes, but also the rating system composed of $\textit{vehicle body}$ and $\textit{area}$ , for instance, has a high number of risk classes (76). The difference between the two is that the latter has classes that are overall not homogeneous. We highlight that by using the ranking provided by the Shapley effects, it is possible to select the two risk factors explaining together the highest fraction of the variance subject to a desirably reduced number of classes.

Table 3. Weighted mean coefficient of variation, mean-squared difference, and goodness of variance fit.

Figure 2. Boxplots showing the differences between the premiums estimated using the complete model and two different reduced models. One of the reduced models is suggested by Heras et al. (Reference Heras, Moreno and Vilar-Zanón2018), and the other is constructed using Shapley values.

Finally, we computed the TAI in Equation (4.11). However, this index is not informative for both the PQRR and the QRPC, since all values are close to 1. For instance, the minimum and the maximum values of the index for the QRPC are 0.9999854 and 0.9999936, respectively (confirming the rating system with vehicle body and agecat as the best one in any case). From now on, we consider as the best reduced model the one including the two risk factors with the highest Shapley effects. By using the indicators in Table 3, there is evidence that the premium that a policyholder is charged in the rating class in this system is “closer” to what this policyholder would have been priced at the individual level. We provide a graphical analysis to investigate this point. In Figure 2, we compare the differences of the individual premiums found with the complete regression model (including all risk factors), with the best reduced model and the model with the two risk factors selected by Heras et al. (Reference Heras, Moreno and Vilar-Zanón2018) (denoted with HMVZ).

Figure 2 illustrates that the model constructed with the two primary risk factors explaining the majority of the variance produces premiums that closely align with the individual premiums estimated using the complete model. This observation holds true irrespective of whether the PQRR or QRPC approach is employed for computing the individual premium. In contrast, the HMVZ model demonstrates a greater disparity from the estimated individual premiums. This is coherent with the findings from Table 3 and suggests that constructing a rating system using variables selected based on global sensitivity indices is more appropriate than relying only on p-values. We perform further analyses to delve deeper into the tail behavior. In Figure 3, we compute the empirical probability that the absolute value of the difference between the premiums exceeds a fixed threshold. Figure 3 shows that the tail of the best model converges to zero with a higher pace with respect to the HMVZ model. This fact confirms empirically our findings in Proposition 4.1 for a non-linear pricing model for which Chebyshev bounds cannot be obtained in a closed form. Finally, we want to demonstrate the validity of our Shapley effects approach even for the case in which the actuarial analyst wants to construct the rating system using more than two risk factors. To show this, we compute the WMCV for all possible rating systems composed of 1–6 rating factors. In Figure 4, we plot the WMCV for each rating system versus the number of selected risk factors in the system, for the insurance prices found via the PQRR approach (left panel) and QRPC approach (right panel), respectively. Let a be the total number of risk factors, $a=6$ in our case, and let b be the number of selected risk factors, $b= 1, \dots, 6$ in our case. The number of possible rating system we can construct for a selected number of risk factor b is given by $\binom{a}{b}=\frac{a!}{b!(a-b)!}$ . For this reason, the number of gray bullets in the plot increases up to $b=3$ and then decreases. The linked black bullets represent the rating system constructed with the risk factors found through the Shapley effects algorithm. As expected, these rating systems are the ones with lower WMCV for every number of selected risk factors. The WMCV shows a decreasing trend as the number of selected risk factors that compose the rating systems increases, since the premium computed within each risk class is clearly more accurate by using more risk factors. This is clearly intuitive. However, our approach shows that the reduction in heterogeneity varies depending on the considered risk factors. Specifically, the downward trend is more pronounced when considering up to three selected risk factors, after which the decrease becomes less significant, particularly in the case of PQRR. This is in accordance with the results obtained via the Shapley effects algorithm, since the three risk factors $\textit{agecat}$ , $\textit{vehicle body}$ , and $\textit{area}$ explain approximately $94 \%$ of the variance in both the PQRR model and QRPC model. Hence, considering also the fourth most important risk factors does not lead to a drastic decrease of the WMCV especially for the PQRR case. Note that if we select four risk factors, the rating system with smaller WMCV is for the PQRR case the one composed by $\textit{agecat}$ , $\textit{vehicle body}$ , $\textit{area}$ , and $\textit{gender}$ , while for the QRPC case the one composed by $\textit{agecat}$ , $\textit{vehicle body}$ , $\textit{area}$ , and $\textit{vehicle age}$ , in accordance with the ranking obtained with the Shapley effects algorithm. Finally, we remark that there is a trade-off between the number of risk classes and the overall homogeneity of the system. One could for instance directly estimate the WMCV without computing the Shapley effects. However, we suggest to estimate the Shapley effects at first and then compute the WMCV of the desired combination of the risk factors. Indeed, the WMVC needs to be estimated for all possible systems and this step has to be repeated every time the analyst wants to include more risk factors in the construction, as it happens in Figure 4.

Figure 3. Empirical probability that the individual premiums, computed with the complete model, differ from their class premiums by more than a threshold $\xi$ , computed with the Shapley model and the HMVZ model. In both panels (PQRR on the left and QRPC on the right), given the same threshold, this probability is higher for the HMVZ model. This implies that ceteris paribus in the latter system more policyholders pay a class premium differing from the individual ones the fixed difference $\xi$ .

Figure 4. WMCV for rating systems composed by a different number of risk factors. Black dots represent the WMCV of each rating system based on PQRR (left panel) and QRPC (right panel). On the x-axis, we display the number of risk factors that compose the rating systems. The dotted line links rating systems constructed using risk factors found via the Shapley effects algorithm applied to insurance prices obtained via PQRR and QRPC.