2.1 Pure premium computation

An insurance company offers to bear part of this risk $g(X)\leq X$ in exchange for a premium which should compensate the average cost of claim given by $p = \mathbb{E}\left[g(X)\right],$ referred to as the pure premium. We consider in this work a function g defined as

\begin{align*}g(x) = \min(\!\max(r\cdot x-d,\,0),\, l),\end{align*}

where $r\in(0,1]$ is the coverage rate, $d\gt0$ is the deductible and $l\gt0$ is the limit. We illustrate the impact of the parameters of the insurance coverage in Example 1.

Example 1. Let us consider a scenario where the risk has a Poisson-lognormal distribution $X\sim\mathsf{Poisson}(\lambda = 3)-\mathsf{LogNorm}(\mu = 0, \sigma = 1)$ and that $n = 100$ insurance coverages are proposed. These are characterized by a rate, a deductible and a limit, set randomly as

\begin{align*}r_i\sim\mathsf{Unif}([0.5, 1])\text{, }d_i\sim\mathsf{Unif}([0.5, 6])\text{, and }l_i = \infty\text{, for }i = 1, \ldots, 100.\end{align*}

Figure 1 shows the pure premiums

(2.2) \begin{equation}p_i = \mathbb{E}(g_i(X)) = \mathbb{E}(\!\min(\!\max(r_i\cdot x-d_i,\,0),\, l_i)),\text{ }i = 1,\ldots, 100.\end{equation}

as a function of the rates and deductibles.

The pure premium are increasing in the rates and decreasing in the deductible. Note that the pure premiums were estimated via simulations to overcome the lack of explicit formula for the distribution functions of X.

In practice, the rate offered to policyholders includes a loading to compensate for the variability of the risks and to cover the management costs. We describe this loading in the next section.

Figure 1.

Pure premiums as a function of the rate of coverage (r) and the deductible (d) for a $\mathsf{Poisson}(\lambda = 3)-\mathsf{LogNorm}(\mu = 0, \sigma = 1)$ risk.

2.2 From pure premiums to commercial premiums

Let $f\,:\,\mathbb{R}_+\mapsto\mathbb{R}_+$ be a nondecreasing function, such that

\begin{align*}\tilde{p} = f(p)\geq p.\end{align*}

The function f is referred to as the loading function. As the commercial premium is a function of the pure premium, then we are applying the expectation premium principle. Other premium principles are also possible like the standard deviation principle discussed in Appendix B. A simple loading function is linear in the pure premium as

\begin{align*}f(x) = (1+\eta)x,\end{align*}

where $\eta\gt0$. The loading functions used by insurance companies are unknown to us and vary from one insurance company to the other. We follow up on Example 1 in Example 2 where we randomize the linear link between pure and commercial premium.

Example 2. Take the pure premiums of Example 1 and apply the following linear loadings

\begin{align*}\eta_i \sim\mathsf{Unif}([0.5, 2]),\text{ for }i=1,\ldots, n.\end{align*}

The commercial premium then relates to the pure premium as

(2.3) \begin{equation}\tilde{p}_i = (1+\eta_i)p_i,\text{ for }i=1,\ldots, n.\end{equation}

Figure 2 displays the commercial premium as a function of the pure premium.

Figure 2.

Pure premium as a function of the commercial premium offered by various insurance companies.

We only observe the commercial premium $\tilde{p}_1,\ldots, \tilde{p}_n$ and we would like to learn from them about the risk X and the loading functions $f_1,\ldots, f_n$. We formulate our problem and describe our solution in the next section.

3.1 Optimization problem based on the pure premium

Assume that we hold a collection of pure premiums $p_{1:n}$ and consider the following optimization problem:

Problem 1. Find $\theta\in \Theta\subset\mathbb{R}^d$ to minimize $d\left(p_{1:n}, p_{1:n}^\theta\right)$, where

\begin{align*}p_i^\theta =\mathbb{E}_\theta\left[g_{i}(X)\right],\text{ for }i = 1,\ldots, n,\end{align*}

are the pure premium associated to the risk X parametrized by $\theta$ and $d(\cdot, \cdot)$ denotes a distance function over the observation space.

We measure the discrepancy between observed and model-generated pure premiums using the root mean square error (RMSE) defined as

(3.1) \begin{equation}\text{RMSE}\left(p_{1:n}, p^\theta_{1:n}\right) = \sqrt{\sum_{i=1}^nw^{\text{RMSE}}_i \left(p_i - p_i^\theta\right)^2},\end{equation}

for a candidate risk parameter $\theta$. The weights $w^{\text{RMSE}}_i \gt 0$ for $i = 1, \ldots, n$ allow us to place greater emphasis on specific data points. The statistical framework is that of minimum distance estimation. We do not have access to the full shape of the data distribution. We must base our inference on specific moments, just as in the generalized method of moments, a popular method among econometricians (see Hansen, Reference Hansen1982). The model is identifiable if there exists a unique estimator $\theta^\ast$ such that

(3.2) \begin{equation}\theta^\ast = \underset{\theta\in \Theta}{\mathrm{arg\,min}}\,\text{RMSE}\left(p_{1:n}, p^\theta_{1:n}\right)\!.\end{equation}

Existence stems from the fact that the parameter space $\Theta$ is compact and the map $\theta\mapsto\text{RMSE}(p_{1:n}, p_{1:n}^\theta)$ is continuous. Uniqueness is more difficult to verify as it depends on the functional $g_i$’s. Given the model for X and the insurance coverages, the pure premium does not have an analytical expression making it difficult to show the convexity of (3.1). A simple necessary condition is that the number of parameters must be smaller than n, the number of moments considered. The shape of our insurance coverage function suggests that we may successfully identify the correct parameters, as demonstrated in Example 3.

Example 3 We consider the same model as in Example 1. Recall that the risk has a Poisson-lognormal distribution $X\sim\mathsf{Poisson}(\lambda = 3)-\mathsf{LogNorm}(\mu = 0, \sigma = 1)$ and that $n = 100$ insurance coverages are proposed. The rates, deductibles, and limit are set randomly as

\begin{align*}r_i\sim\mathsf{Unif}([0.5, 1])\text{, }d_i\sim\mathsf{Unif}([0.5, 6])\text{, and }l = \infty\text{, for }i = 1, \ldots, 100.\end{align*}

Let us assume that $\mu =0$ and compute the pure premium over a grid of values for $\lambda$ and $\sigma$. Figure 3 shows the $\text{RMSE}[p_{1:n}, p_{1:n}^\theta]$ depending on the value of $\lambda$ and $\sigma$ for $(\lambda,\sigma)\in [0,5]\times [0,2]$. This contour plot shows minimal RMSE values around the true parameter values $(\lambda_0, \sigma_0) = (3, 1)$.

Figure 3.

Contour plot of $\text{RMSE}(p_{1:n}, p_{1:n}^\theta)$ for $\mu = 0$ and $(\lambda,\sigma)\in [0,5]\times [0,2]$.

If uniqueness of the solution cannot be verified, a workarround consists in adding regularization terms to the discrepancy measure. We explore this direction in Section 3.2 where the actual problem is treated. The issue is further addressed by using a particle-based optimization methods to search the parameter space and provide a set of admissible candidate parameters. Such an optimization method is presented in Section 3.3. This identifiability issue from an empirical point of view in the online Supplementary Material.Footnote ¹

3.2 Optimization problem based on the commercial premiums

Assume that we hold a set of commercial rates $\tilde{p}_{1:n} = \{\tilde{p}_{1},\ldots, \tilde{p}_{n}\}$ defined as

\begin{align*}\tilde{p}_i = f_i(p_i) = f_i\left\{\mathbb{E}\left[g_{i}(X)\right]\right\},\text{ }i = 1,\ldots, n.\end{align*}

and consider the following optrimization problem.

Problem 2. Find $\theta\in \Theta\subset\mathbb{R}^d$ and $f\,:\,\mathbb{R}_+\mapsto \mathbb{R}_+$ to minimize $d\left[\tilde{p}_{1:n}, f\left(p_{1:n}^\theta\right)\right]$, where the function f is applied elementwise on $p_{1:n}^\theta$ and $d(\cdot, \cdot)$ denotes a distance function over the observation space, subject to

(3.3) \begin{equation}\tilde{p}_i \geq p_i^\theta,\text{ and } f(p_i^\theta) \geq p_i^\theta, \text{for }i = 1,\ldots, n.\end{equation}

Our first task is to find a generic function f to represent the safety loading functions $f_i$’s used by the competitors. For this, we use isotonic regression. It is a statistical technique used for fitting a nondecreasing function to a set of data points. The idea is that if two pure premiums satisfy $p_i\leq p_j$, then the commercial premiums should also verify $\tilde{p}_i\leq \tilde{p}_j$. Consider a collection of candidate pure premiums $p_{1:n}^\theta$, associated to a candidate estimate of the risk parameter $\theta$. Our data points are therefore pairs of pure and commercial premiums $(p_i^\theta,\tilde{p}_i)_{i=1,\ldots, n}$. Suppose the pure premiums have been ordered such that $p_i^\theta\leq p_j^\theta$ for $i\leq j$, isotonic regression seeks a least square fit $\tilde{p}^\theta_i$ for the $\tilde{p}_i$’s such that $\tilde{p}^\theta_i\leq \tilde{p}^\theta_j$ for $p_i^\theta\leq p_j^\theta$. It reduces to finding $\tilde{p}^\theta_1,\ldots, \tilde{p}^\theta_n$ that minimizes

\begin{align*}\sum_{i=1}^nw_i^{\scriptscriptstyle{\textrm{ISO}}}(\tilde{p}^\theta_i- \tilde{p}_i)^2\text{, subject to }\tilde{p}_i^\theta\leq \tilde{p}^\theta_j\text{ whenever }p_i^\theta\leq p_j^\theta,\end{align*}

where $\left(w_i^{\scriptscriptstyle{\textrm{ISO}}}\right)_{i = 1, \ldots, n}$ denotes the weights associated to each pair $(p_i^\theta,\tilde{p}_i)_{i=1,\ldots, n}$. The weights allow us to emphasize on specific data points in the same manner as the weights defined for the RMSE in (3.1). Since the $p_i^\theta$’s fall in a totally ordered space, a simple iterative procedure called the pool adjacent violators algorithm (PAVA) can be used. The algorithm that describes the procedure is provided below:

1. 1. Initialize the sequence of values to be the same as the data points $\tilde{p}_i^\ast =\tilde{p}_i$.

2. 2. Iterate through the sequence and identify “violations,” which occur when the current value is greater than the next value, that is

   \begin{align*}\tilde{p}_i^\ast\gt\tilde{p}_{i+1}^\ast\text{ for some }i=1,\ldots, n.\end{align*}
   When a violation is found, adjust the values in the associated segment of the sequence to be the average of the values,
   \begin{align*}\tilde{p}_i^\ast \leftarrow (\tilde{p}_i^\ast + \tilde{p}_{i+1}^\ast)/2,\end{align*}
   ensuring monotonicity.

3. 3. Repeat Step 2 until no violations are left.

We use the isoreg function from R to get the fitted values $\tilde{p}^\theta_i,\text{ }i=1,\ldots, n$. To complete the isotonic regression task, we shall find a function f such that $f(p_i^\theta)=\tilde{p}^\theta_i$. A common choice is a piece-wise constant function that interpolates the $\tilde{p}^\theta_i$’s. An illustration is provided in Example 4, and we outline the advantages and limitations associated to this choice in Remark 3.1.

Example 4. The isotonic fit of the data of Examples 1 and 2 is provided in Figure 4.

Figure 4.

Isotonic link between the pure and commercial premiums.

We now turn to the definition of a discrepancy measure to compare the model-generated and observed market commercial premiums. Our starting point is the RMSE defined as

(3.4) \begin{equation}\text{RMSE}\left[\tilde{p}_{1:n}, f\left(p^\theta_{1:n}\right)\right] = \sqrt{\sum_{i=1}^nw^{\text{RMSE}}_i\left[\tilde{p}_i - f\left(p_i^\theta\right)\right]^2},\end{equation}

for candidate risk parameter $\theta$ and the isotonic fit f. We seek

\begin{align*}\theta^\ast = \underset{\theta\in \Theta}{\mathrm{arg\,min}}\,\text{RMSE}\left[\tilde{p}_{1:n}, f\left(p^\theta_{1:n}\right)\right].\end{align*}

The existence of such $\theta^\ast$ is guaranteed because $\theta\mapsto \text{RMSE}\left[\tilde{p}_{1:n}, f\left(p^\theta_{1:n}\right)\right]$ only takes a finite number of values. Indeed, to each $\theta \in \Theta$ is associated a unique permutation $s^\theta\in S_n$, where $S_n$ denotes the set of all the permutations of $\{1,\ldots, n\}$, such that

\begin{align*}p^\theta_{s^\theta(1)}\leq\ldots\leq p^\theta_{s^\theta(n)}.\end{align*}

This permutation $s^\theta$ defines a unique isotonic fit f based on

\begin{align*}\tilde{p}^\theta_{s^\theta(1)}\leq\ldots\leq \tilde{p}^\theta_{s^\theta(n)},\end{align*}

leading to a given RMSE value $\text{RMSE}\left[\tilde{p}_{1:n}, f\left(p^\theta_{1:n}\right)\right]$. Concretely, for $\theta_1,\theta_2\in \Theta$, if it holds that $s^{\theta_1} = s^{\theta_2}$, then $\text{RMSE}\left[\tilde{p}_{1:n}, f(p^{\theta_1}_{1:n})\right] = \text{RMSE}\left[\tilde{p}_{1:n}, f(p^{\theta_2}_{1:n})\right]$. The application $\theta\mapsto s_n^\theta$ is surjective since $S_n^\Theta = \{s_n^\theta\text{ ; }\theta\in \Theta\}$ is finite. The fact that $\Theta$ is a continuous space implies that $\theta^\ast$ cannot be unique. The application of isotonic regression offers a straightforward interpretation of our objective, as detailed in Remark 3.2.

Our problem is an ill-posed inverse problem. Ill-posedness is usually dealt with by adding a regularization to the objective function that one wants to minimize. The ratio of $p / \tilde{p}$ corresponds to what practitioners would call the expected loss ratio (LR). Our solution is based on targeting a given loss ratio. The loss ratio is a standard measure to assess the profitability of insurance lines of business. An insurance company that enters a new market is likely to have insights on the loss ratio relative to this market, for example, by having informal discussions with reinsurers, brokers, or competitors. These feedbacks may translate into the definition of a lower and upper bound denoted by $\text{LR}_{\text{low}}$ and $\text{LR}_{\text{high}}$, respectively. We can then assume that the loss ratios $\text{LR}_i = p_i / \tilde{p}_i$, for $i = 1, \ldots, n$, should fall in the range $[\text{LR}_{\text{low}}, \text{LR}_{\text{high}}]$, which we refer to as the loss ratio corridor. Assuming that $\text{LR}_{\text{high}} \lt 1$, we ensure both constraints in (3.3) by adding to our distance (3.4) two regularization terms defined as

\begin{equation*}\text{Reg}_{\text{low}}\left(\tilde{p}_{1:n}, p_{1:n}^\theta\right) = \sqrt{\sum_{i=1}^nw_i^{\text{RMSE}}\left(\tilde{p}_i - p_i^\theta \cdot \text{LR}_{\text{low}}^{-1}\right)_+^2},\end{equation*}

and

\begin{equation*}\text{Reg}_{\text{high}}\left(\tilde{p}_{1:n}, p_{1:n}^\theta\right) = \sqrt{\sum_{i=1}^nw_i^{\text{RMSE}}\left(p_i^\theta \cdot \text{LR}_{\text{high}}^{-1} - \tilde{p}_i\right)_+^2},\end{equation*}

where $(x)_+ = \max(x,0)$ denotes the positive part of x. The distance we consider within Problem 2 is now given by

\begin{align*}d\left[\tilde{p}_{1:n}, f\left(p_{1:n}^\theta\right)\right] = \text{RMSE}\left[\tilde{p}_{1:n}, f\left(p_{1:n}^\theta\right)\right] + \text{Reg}_{\text{low}}\left(\tilde{p}_{1:n}, p_{1:n}^\theta\right) + \text{Reg}_{\text{high}}\left(\tilde{p}_{1:n}, p_{1:n}^\theta\right).\end{align*}

We illustrate the impact of adding the regularization terms in Example 5.

Example 5. We consider the commercial rates of Example 2. Recall that the commercial premiums are given by

(3.5) \begin{equation}\tilde{p}_i = (1+\eta_i)p_i, \quad \text{for} \quad i=1, \ldots, n,\end{equation}

where the pure premiums are those of Example 3 and

\begin{align*}\eta_i \sim \mathsf{Unif}([0.5, 2]), \quad \text{for} \quad i=1, \ldots, n.\end{align*}

Figure 5 displays the contour plot of the discrepancy between observed and model-generated commercial rates. When comparing Figure 5(a) and (b), we note how beneficial including regularization terms is to identify the true parameter values.

Figure 5.

Contour plot of $\text{RMSE}\left[\tilde{p}_{1:n}, f\left(p^\theta_{1:n}\right)\right]$ and $d\left[\tilde{p}_{1:n}, f\left(p_{1:n}^\theta\right)\right]$ for $\mu = 0$ and $(\lambda, \sigma) \in [0, 5] \times [0, 2]$.

Regularization brings us closer to the scenario described in Section 3.1, where the pure premium is known. This corresponds to the case where $\text{LR}_{\text{low}} = \text{LR}_{\text{high}} = 1$. Regularization enables us to exclude large portions of the parameter space associated with nonsensical pure premiums given the commercial premiums. However, careful consideration is required when setting the loss ratio corridor. A narrow loss ratio corridor results in a highly precise estimation, but this estimation may be biased if the loss ratio corridor is misspecified.

The issue of identifiability persists, as multiple parameter values can still solve the optimization problem. This highlights the necessity of using particle-based optimization techniques, as described in Section 3.3, to explore the parameter space. Such techniques return a set of admissible candidate parameters and are less sensitive to initialization. For a more detailed discussion of the identifiability issue, we refer the reader to the online Supplementary Material.Footnote ²

3.3 Population Monte Carlo approximate Bayesian computation algorithm

Our solution alternates between proposing parameter values for the risk to compute the pure premiums and approximating the $f_i$’s using isotonic regression. We must accommodate the lack of tractable expressions for the pure premium

\begin{align*}p_i^\theta = \mathbb{E}_\theta\left[g_i(X)\right], \quad \text{for} \quad i = 1, \ldots, n.\end{align*}

The use of numerical methods makes a grid search procedure prohibitive from a computing time point of view. It also prevents us from using gradient-based optimization procedures. In such cases, one can turn toward particle swarm optimization algorithms or genetic algorithms to search the parameter space. Since we have decided to take a Crude Monte Carlo estimator for the pure premiums, the accuracy depends on the number of replications R of X being used. We adopt a Bayesian strategy in order to reflect the uncertainty around the pure premium calculation onto the parameters’ final estimates. Our algorithm is similar to ABC algorithms, and we simply refine the procedure laid out in the introduction to get an approximation of the posterior distribution $\pi(\theta|\mathcal{D})$.

We start by setting a prior distribution $\pi(\theta)$ over the parameter space that we sequentially improve through intermediate distributions characterized by a sequence of tolerance levels $(\epsilon_g)_{g\geq0}$ that decrease gradually as $\infty = \epsilon_0 \gt \epsilon_1 \gt \epsilon_2 \gt \ldots \gt 0$. Each intermediate distribution (called a generation and denoted by g) is represented by a cloud of weighted particles$\left(\theta_j, w_j^g\right)_{j=1, \ldots, J}$. We approximate each intermediate posterior distribution using a multivariate kernel density estimator (kde) denoted by $\pi_{\epsilon_g}(\theta|\mathcal{D})$. The parameters of the algorithm are the number of generations G, the population size J (the number of particles in the cloud), and the number of Monte Carlo replications R of X.

The algorithm is initialized by setting $\epsilon_0 = \infty$ and $\pi_{\epsilon_0}(\theta|\mathcal{D}) = \pi(\theta)$. For generation $g \geq 1$, we hold an intermediate distribution $\pi_{\epsilon_{g-1}}(\theta|\mathcal{D})$ from which we can sample particles $\theta^\ast \sim \pi_{\epsilon_{g-1}}(\theta|\mathcal{D})$. We compute the associated pure premium

\begin{align*}p_i^{\theta^{\ast}} = \mathbb{E}_{\theta^\ast}\left[g_i(X)\right], \quad \text{for} \quad i = 1, \ldots, n.\end{align*}

The pure premiums are computed via Monte Carlo simulations. The accuracy depends on the number R of copies of X involved in the Monte Carlo estimations. We then fit the isotonic regression model

\begin{align*}\tilde{p}_i = f\left(p_i^{\theta^{\ast}}\right) + e_i, \quad \text{for} \quad i = 1, \ldots, n,\end{align*}

where $e_i$ is an error term that captures the mismatch between the true value of the pure premium and its empirical counterpart estimated by the competitor insurance company using its historical data and the company-specific loading function. We further compare the observed commercial premiums to the model-generated ones via the distance defined in Section 3.2 with

\begin{align*}d\left[\tilde{p}_{1:n}, f\left(p_{1:n}^{\theta^\ast}\right)\right] = \text{RMSE}\left[\tilde{p}_{1:n}, f\left(p_{1:n}^{\theta^\ast}\right)\right] + \text{Reg}_1\left(\tilde{p}_{1:n}, p_{1:n}^{\theta^\ast}\right) + \text{Reg}_2\left(\tilde{p}_{1:n}, p_{1:n}^{\theta^\ast}\right).\\[8pt]\end{align*}

If the distance satisfies $d\left[\tilde{p}_{1:n}, f\left(p_{1:n}^{\theta^\ast}\right)\right] \lt \epsilon_{g-1}$, then we keep the associated particle $\theta^\ast$. New particles are proposed until we reach J accepted particles denoted by $\theta_1^g, \ldots, \theta^g_J$. We also store the distances $d_1^g, \ldots, d_J^g$. We need to set the next tolerance threshold $\epsilon_g$, which is used to calculate the particle weights

\begin{align*}w_j^g \propto \frac{\pi(\theta^g_j)}{\pi_{\epsilon_{g-1}}(\theta)}\mathbb{I}_{d_j^g \lt \epsilon_{g-1}}, \quad j = 1, \ldots, J.\end{align*}

The tolerance threshold is chosen so as to maintain a specified effective sample size (ess) of $J/2$ as in Del Moral et al. (Reference Del Moral, Doucet and Jasra2012a). Following Kong et al. (Reference Kong, Liu and Wong1994), the ess is estimated by $1/\sum_{j=1}^J(w_j^g)^2$. This weighted sample then allows us to update the intermediate distribution as

\begin{align*}\pi_{\epsilon_g}(\theta|\mathcal{D}) = \sum_{j=1}^J w_j^g K_H(\theta - \theta^g_j),\end{align*}

where $K_H$ is a multivariate kde with smoothing matrix H. A common choice for the kde is the multivariate Gaussian kernel with a smoothing matrix set to twice the empirical covariance matrix of the cloud of particles $\{\theta_j^g, w_j^g\}$ as in Beaumont et al. (Reference Beaumont, Cornuet, Marin and Robert2009). The procedure is summarized in Algorithm 1.

Algorithm 1

Population Monte Carlo Approximate Bayesian Computation

The user must configure several aspects of the algorithm. The prior assumptions $\pi(\theta)$ determine the parameter space that will be searched. The loss ratio corridor $[\text{LR}_{\text{low}}, \text{LR}_{\text{high}}]$ sets up the two regularization terms, ensuring that parameters associated with unreasonable pure premiums are excluded. The prior settings and loss ratio corridor can be guided by expert opinions. The population size J drives the quality of the posterior distributions approximations through the cloud of particles. A large J also enhances the chances of finding global optimums, as more particles improve the coverage of the parameter space. A greater number R of Monte Carlo simulations ensures the accuracy of the pure premium evaluation. Both R and J contribute to the stability of the algorithm’s results over several runs. The number of generations G relates to the tolerance level $\epsilon$, which in turn drives the narrowness of the posterior distribution output by the ABC algorithm. As one would expect, the computational time for the algorithm increases with higher values of R, J, and G. Therefore, the choice of suitable values for G, J, and R can be made in consideration of a predetermined computational time budget. A practical solution to set G is to stop the algorithm whenever the difference between two consecutive tolerance levels is lower than some threshold $\Delta_\epsilon$ or if we reach a minimum tolerance level $\epsilon_{\min}$. Convergence results for ABC algorithms are readily available in the literature as discussed in Remark 3.3.

Remark 3.3. The final output of our ABC algorithm is the following expression:

(3.6) \begin{equation}\pi_{\epsilon_G}(\theta|\mathcal{D}) = \sum_{j=1}^J w_j^G K_H(\theta - \theta^G_j),\end{equation}

which represents an approximate posterior distribution of $\theta$ given an iid sample $x_{1:R}$ of X. The data $\mathcal{D}$ can be interpreted as a set of summary statistics calculated based on a sample $x_{1:R}$ of size R, which corresponds to the number of Monte Carlo replications used to calculate the pure premium. Several types of convergence can be studied:

1. 1. $R \to \infty$

2. 2. $\epsilon_G \to 0$

3. 3. $J \to \infty$

Convergence (1) determines the accuracy of $d\left[\tilde{p}_{1:n}, f\left(p_{1:n}^{\theta^\ast}\right)\right]$ when using simulation to evaluate the pure premium $p_{1:n}^\theta$. These are Monte Carlo estimators of $\mathbb{E}_\theta[g_i(X)]$ for $i = 1, \ldots, n$, and convergence occurs at a rate proportional to $1/\sqrt{R}$.

Regarding convergence (2), as $\epsilon_G$ decreases toward 0, our ABC estimator converges to the distribution of $\theta$ conditional on $\mathcal{D}$. This is equivalent to the true posterior distribution only if $\mathcal{D}$ consists of sufficient statistics for the risk model X. Rubio and Johansen (Reference Rubio and Johansen2013, Proposition 2) justify justify the use of non-sufficient statistics, under certain conditions, to make inferences about $\theta$. Their results remain valid when using a kernel density estimator such as (3.6) to represent the posterior distribution. The bias of the ABC posterior has been estimated to be $\mathcal{O}(\epsilon_G^2)$ in Barber et al. (Reference Barber, Voss and Webster2015).

Convergence (3) pertains to the empirical probability measure approaching the true probability measure. Central limit theorems for this convergence are discussed in Del Moral et al. (Reference Del Moral, Doucet and Jasra2006) and (Reference Del Moral, Doucet and Jasra2012b).

We illustrate the posterior distribution evolution along the algorithm iterations in Example 6.

Example 6. We follow up on Examples 3 and 5. We aim to fit the model

\begin{align*}X\sim\mathsf{Poisson}(\lambda)-\mathsf{LogNorm}(\mu = 0, \sigma).\end{align*}

The prior assumptions are as follows

\begin{align*}\lambda \sim \mathsf{Unif}([0,10]),\text{ and }\sigma \sim \mathsf{Unif}([0,5]).\end{align*}

The algorithm parameters are set to

\begin{align*}J = 1000\text{, }R = 1000\text{, }\Delta_\epsilon = 0.1\text{, and }LR\in [0.3, 0.66].\end{align*}

The algorithm halts at the $9{\text{th}}$ generation reaching a tolerance level of $\epsilon = 0.87$. Figure 6 shows the sequence intermediate posterior distributions for $\lambda$ and $\sigma$.

Figure 6.

Intermediate posterior distributions of $\lambda$ and $\sigma$.

After the algorithm terminates, it is customary to focus on the last generations of particles for inference. Pointwise estimators are derived from this final set of particles. Two commonly used estimators include the Mean A Posteriori (map) obtained by averaging the particles in the last cloud and the Mode A Posteriori (mode), which is the mode of the empirical distribution within the final cloud of particles. The simulation study, conducted in the following section, is designed to investigate the convergence behavior and to compare the characteristics of the map and mode estimators.

5.1 Evolution and growth of the pet insurance market

Pet insurance is a product designed to cover the costs of veterinary care for pets. It operates on a similar principle to human health insurance, providing a way for pet owners to manage the financial risks associated with unexpected medical expenses for their animals. Usually the expenses are covered in case of an accident or a disease. Pet owners can choose from different policy options based on their budget and coverage needs. Policies may vary in terms of deductibles (d), coverage limits (l), and coverage rates (r). The cost of premiums can depend on various factors, including the pet’s age, breed, health condition, and the level of coverage selected.

The pet insurance market has been witnessing significant growth globally, driven by increasing pet ownership (especially with so-called pandemic pets, i.e. animals adopted during 2020 lockdowns), rising veterinary costs, and the changing role that a pet plays in a families social structure. This latter factor is also influenced by changing societal views and the increased awareness of the importance of health and welfare of pets, which in turn comes with increased consideration of regular veterinary health checks. In order to offset the cost associated with such expenditures, there has begun to be a broader interest in households purchasing pet insurance.

To date, the adoption and acceptance of pet insurance still varies significantly across regions of the world. Nordic countries, such as Sweden, have historically had a very high penetration rate with around 70% of pets insured. Some Anglo-Saxon countries (UK and Germany mostly) have seen significant growth in the pet insurance market during the last decades, leading to 30% of penetration rate. Other developed countries, like France, have significantly lower market sizes, with less than 10% of pets that are insured, which suggests high growth potential. The market place for pet insurance in the USA is currently also experiencing sustained growth. According to MarketWatch nationwide survey,Footnote ³ about $44.6\%$, of pet owners stated they currently have pet insurance. The North American Pet Health Insurance Association (NAPHIA) conducted a survey for its 2022 State of the Industry Report and found that over $4.41$ million pets were insured in North America in 2021, up from $3.45$ million in 2020. The report also revealed that pet insurance premiums totaled $\$2.84$ billion in 2021, marking a $30.5\%$ increase from the previous year.

This growth continues to spur increases in the capital investments associated with such an insurance line of business:

• • in Sweden, Lassie has raised 11m euros in 2022 and 23m euros in 2023 ;

• • in the UK, ManyPets has raised $350m at a valuation higher than $2bn in 2021;

• • in France, Dalma has raised 15m euros in 2022.

• • JAB Holding Company has invested around 2 billion dollars in 2021 to create the Pinnacle Pet Group and the Independence Pet Holdings. Their purpose is to become the pet insurance leaders, respectively, in Europe and North America through multiple acquisitions of historic players.

Hence, the pet insurance market is becoming more competitive with an increasing number of insurance companies or brokers offering pet insurance policies. To capture new market share as a new agent, it is essential to offer differentiated products, such as coverage combinations with no deductibles and higher limits.

5.2 Data description

During the week of the 18 May 2024, we have collected 1080 quotes from five insurance companies. Each row of our datasets corresponds to a yearly premium collected from some insurance company website associated to a specific insurance coverage and a specific dog. We therefore find the coverage parameters that are the coverage rate r, the deductible d, and the limit l. Recall that the compensation for an annual expense of amount X is calculated as $\min\left[\max(r \cdot X - d, 0), l\right]$. We also have the rating factors that reduces for pet insurance in France to specie, breed, age, and gender. Table 1 provides a list of the variables in the datasets.

Table 1.

List of the variables of our datasets.

The first five rows are given in Table 2.

Table 2.

First five rows of our datasets.

We have collected rates associated to 12 risk classes given in Table 3.

Table 3.

The 12 risk classes under study.

It means that we have 90 quotes to study each risk class. Figure 9 provides a visual overview of the range of insurance coverage options available in the pet insurance market.

Figure 9.

Overview of the insurance coverages offered by the five insurance companies operating in the French market under study.

Remark 5.1 briefly discusses data quality.

We conduct two separate studies. In Section 5.3, we focus on a specific risk class associated to a female, 4-year-old Australian sheperd. Several claim models are compared. One is selected to look into the pricing strategies of the actors. In Section 5.4, we look into the quotes of various risk classes that we investigate using a single model.

5.3 Analysis of one risk class using several models

We study 90 quotes from five insurers operating in the pet insurance market for a specific risk class associated to a 4-year-old female Australian Shepherd. We need to make some parametric assumptions to model claim frequency and severity. Classical claim frequency distributions include the Poisson, Binomial, and Negative Binomial distributions, which allow us to accommodate equidispersion, underdispersion, and overdispersion, respectively. For claim severity, we have chosen the Gamma and Lognormal distributions. The Gamma distribution is a common choice for modeling claim severity when using GLM, but it is characterized by a light tail. The Lognormal distribution has a thicker tail, making larger claim sizes more likely to occur. We limit ourselves to two-parameter models: one parameter for claim frequency and another for claim severity. So we take the $\mathsf{Poisson}(\lambda)$ to model the claim frequencyFootnote ⁴ with the following prior setting:

(5.1) \begin{equation}\lambda \sim\mathsf{Unif}([0,10]).\end{equation}

We consider three claim severity distributions including $\mathsf{LogNorm}(\mu =0, \sigma)$, $\mathsf{LogNorm}(\mu, \sigma = 1)$, and $\mathsf{Gamma}(\alpha, \beta = 1)$. The prior settings over the parameters of the claim size distributions are as follows:

(5.2) \begin{equation}\mu \sim \mathsf{Unif}([-10,10]), \quad \sigma \sim \mathsf{Unif}([0,10]), \quad \text{and} \quad \alpha \sim \mathsf{Unif}([0, 10^5]).\end{equation}

Combining the distributions for the claim frequency and severities results in a total of three loss models. The population size in the abc algorithm is set to $J = 1000$. The pure premiums are computed using $R = 2000$ Monte Carlo replications. The algorithm stops whenever the difference between two consecutive tolerance levels is lower than $\Delta_\epsilon = 1$. The bounds for the loss ratio corridor are set to $\text{LR}_{\text{low}} = 40\%$ and $\text{LR}_{\text{high}} = 70\%$. The posterior distributions of the parameters for each model are provided in Figure 10.

For all the models, the algorithm updates the prior distribution in an informative way. Table 4 provides the tolerance levels (ranked in increasing order) during the last iteration of the abc algorithm for the loss models.

Table 4.

Tolerance level during the last iteration of the abc algorithm of each loss model.

Figure 10.

Posterior distribution of the parameters of the loss models when fitted to the pet insurance dataset.

The final tolerance levels lie between $93.72$ and $113.15$, which is higher than the tolerance obtained in the simulation study which was around 33 for 50 data points and 50 for 200 data points. This discrepancy indicates misspecifications that stem from our assumptions about insurance companies adhering to the expectation principle for premium calculation and the models employed for claim frequency and claim amounts. Table 5 reports the estimations of the parameters of all the model using the map and the mode.

Table 5.

map and mode estimator for the parameters of the loss models.

Table 6 reports the estimations of the average total claim amounts and the average loss ratio for all the models for all models when fitted using the map and the mode.

Table 6.

map and mode estimators of the average loss ratio and average total claim amounts.

We note that the risk level, characterized here by the expected total claim amount, is similar for all the models, maybe a bit higher for the model having the $\mathsf{LogNorm}(\mu = 0,\sigma)$ as claim sizes distribution. We further look into the loading function approximated via the isotonic regression. We estimate the pure premium for each model using the map as an estimator of the model parameters, and we plot the isotonic regression function to explain the commercial premium in Figure 11.

Figure 11.

Isotonic link between pure and commercial premium for the different loss models.

The isotonic fits of the loading function across all the models are similar which means that the models all agree on a common ordering of the pure premiums of the various insurance coverages. To highlight the explanatory power of our methodology, let us focus on the $\mathsf{Poisson}(\lambda)-\mathsf{LogNorm}(\mu,\sigma=1)$ loss model. Note that the choice of the loss model is somewhat arbitrary because the information extracted from the data in Figure 11 is relatively consistent across the considered models. In Figure 12, we present a plot that illustrates the relationship between the commercial premium and the pure premium for the $\mathsf{Poisson}(\lambda)-\mathsf{LogNorm}(\mu,\sigma=1)$ model. Different insurance companies are indicated by distinct colors, providing a visual representation of each company’s respective rates.

Figure 12.

Commercial premium as a function of the pure premium for the $\mathsf{Poisson}(\lambda)-\mathsf{LogNorm}(\mu,\sigma=1)$ depending on the insurance carrier.

The accuracy of the loss model fitting enables us to condense the three-dimensional information of the rate of coverage, deductible, and limit into a single metric: the pure premium. Subsequently, isotonic regression unveils the relationship between commercial and pure premiums, providing a link between the two. The distinctions among various players in the pet insurance market come to light through the color-coded points, offering insights into the pricing strategies adopted by industry participants.

5.4 Analysis of several risk classes with one model

The $\mathsf{Poisson}(\lambda)-\mathsf{LogNorm}(\mu,\sigma=1)$ model is fitted to the data within each risk classes (90 quotes) of Table 3. The prior settings are given by

\begin{align*}\lambda \sim\mathsf{Unif}([0,10]) \text{, and }\mu\sim\mathsf{Unif}([-10,10]).\end{align*}

The algorithm’s hyperparameters are similar to that of the previous subsection with

\begin{align*}J = 1000, \text{ }R = 2 000, \text{ and }\Delta_\epsilon = 1.\end{align*}

Figure 13 shows the posterior predictive distribution of the expected total claim amounts and the averaged loss ratio within each risk class.

Figure 13.

Posterior predictive distribution of $\mathbb{E}(X)$ and average loss ratio within each risk class.

Figure 13(a) allows us to compare the different risk classes. We note that an older dog is more expensive on average and that the breeds may be ordered as Australian Sheperd, Golden-Retriever, German Sheperd, and French Bulldog in terms of riskiness. Figure 13(b) indicates that the loss ratios are arround 62–65% for all the risk classes.