2.1 Mathematical formulation of IPH distributions

Let $ ( J_t )_{t \geq 0}$ be a time-inhomogeneous Markov pure-jump process on the finite state space $\{1, \dots, p, p+1\}$ , where states $1,\dots,p$ are transient and $p+1$ is absorbing. In our setting, such a process can be regarded as an insurance claim evolving through different states during its lifetime until the absorbing state is reached, which determines the total size of the claim. In its most general form, the transition probabilities of the jump process:

\begin{align*}p_{ij}(s,t)={\mathbb P}(J_t=j|J_s=i),\quad i,j\in\{1,\dots,p+1\},\end{align*}

may be written in succinct matrix form in terms of the product integral as follows:

\begin{align*}\textbf{P}(s,t)=\prod_{s}^{t}(\boldsymbol{I}+\boldsymbol{\Lambda}(u) d u)\;:=\;\boldsymbol{I}+\sum_{k=1}^{\infty} \int_{s}^{t}\! \int_{s}^{u_{k}} \!\cdots \int_{s}^{u_{2}}\! {\boldsymbol\Lambda}\!\left(u_{1}\right) \cdots {\boldsymbol\Lambda}\!\left(u_{k}\right) d u_{1} \cdots \textrm{d} u_{k},\end{align*}

for $s<t,$ where $\boldsymbol{\Lambda}(t)$ is a matrix with negative diagonal elements, nonnegative off-diagonal elements, and rows sums equal to 0, commonly referred to as an intensity matrix. Since the process $ ( J_t )_{t \geq 0}$ has a decomposition in terms of transient and absorbing states, we may write

\begin{align*} \boldsymbol{\Lambda}(t)= \left( \begin{array}{c@{\quad}c} \textbf{T}(t) & \textbf{t}(t) \\ \textbf{0} & 0 \end{array} \right)\in\mathbb{R}^{(p+1)\times(p+1)}\,, \quad t\geq0\,,\end{align*}

where $\textbf{T}(t) $ is a $p \times p$ sub-intensity matrix and $\textbf{t}(t)$ is a p-dimensional column vector providing the exit rates from each state directly to $p+1$ . The intensity matrix has this form, since there are no positive rates from the absorbing state to the transient ones, and since the rows should sum to 0, we have that for $t\geq0$ , $\textbf{t} (t)=- \textbf{T}(t) \, \textbf{e}$ , where $\textbf{e} $ is a p-dimensional column vector of ones. Thus, $\textbf{T}(t)$ is sufficient to describe the dynamics of process after time zero. We write $\textbf{e}_k$ in the sequel for the k-th canonical basis vector of $\mathbb{R}^p$ , so in particular $\textbf{e}=\sum_{k=1}^p\textbf{e}_k$ .

If the matrices $\textbf{T}(s)$ and $\textbf{T}(t)$ commute for any $s<t$ , we may write

\begin{align*}\textbf{P}(s,t)=\begin{pmatrix} \exp\left(\int_s^t\textbf{T}(u)du\right)\;\;\;\; & \textbf{e}-\exp\left(\int_s^t\textbf{T}(u)du\right) \textbf{e} \\ \textbf{0}\;\;\;\; & 1 \end{pmatrix},\quad s<t.\end{align*}

The matrices $(\textbf{T}(t))_{t\ge0}$ together with the initial probabilities of the Markov process, $ \pi_{k} = {\mathbb P}(J_0 = k)$ , $k = 1,\dots, p$ , or in vector notation:

\begin{align*}\boldsymbol{\pi} = (\pi_1 ,\dots,\pi_p )^{\textsf{T}},\end{align*}

fully parametrize IPH distributions, and with the additional assumption that ${\mathbb P}(J_0 = p + 1) = 0$ we can guarantee that absorption happens almost surely at a positive time. We thus define the positive and finite random variable given as the absorption time as follows:

\begin{align*} Y = \inf \{ t > 0 : J_t = p+1 \},\end{align*}

and say that it follows an IPH distribution with representation $(\boldsymbol{\pi},(\textbf{T}(t))_{t\ge0})$ . For practical applications, however, we focus on a narrower class which allows for effective statistical analysis. We make the following assumption regarding the inhomogeneity of the sub-intensity matrix:

\begin{align*}\textbf{T}(t) = \lambda(t)\,\textbf{T},\end{align*}

with $\lambda(t)$ some parametric and positive function such that the map:

\begin{align*}y\mapsto\int_0^y \lambda(s)ds\in(0,\infty),\quad\forall y>0,\end{align*}

converges to infinity as $y\to\infty$ , and $\textbf{T}$ is a sub-intensity matrix that does not depend on time t. When dealing with particular entries of the sub-intensity matrix, we introduce the notation:

\begin{align*}\textbf{T}=(t_{kl})_{k,l=1,\dots,p},\quad \textbf{t}=-\textbf{T}\textbf{e}=(t_1,\dots,t_p)^{{\textsf{T}}},\quad \alpha_i=-t_{ii}>0,\:\:\: i=1,\dots,p.\end{align*}

We then write

\begin{align*}Y \sim \mbox{IPH}(\boldsymbol{\pi} , \textbf{T} , \lambda ),\end{align*}

where IPH stands for inhomogeneous phase type. Similarly, we write PH in place of phase type in the sequel, that is, for the case when $\lambda$ is a constant.

IPH distributions are the building blocks to model claim severity for two reasons. Firstly, they have the interpretation of being absorption times of an underlying unobserved multistate process evolving through time, which can be very natural if the losses are linked with an evolving process, such as a lifetime or a legal case; and secondly, PH distributions are the archetype generalization of the exponential distribution in that many closed-form formulas are often available in terms of matrix functions, and thus, when armed with matrix calculus and a good linear algebra software, a suite of statistical tools can be developed Please provide the expansions for BFGS, AIC, BIC, GLM, and IBNR if necessary.explicitly.

It is worth mentioning that PH distributions possess a universality property, since they are known to be dense on the set of positive-valued random variables in the sense of weak convergence. Consequently, they possess great flexibility for atypical histogram shapes that the usual probabilistic models struggle to capture. Do note, however, that the tail can be grossly misspecified if no inhomogeneity function is applied, since weak convergence does not guarantee tail equivalence between the approximating sequence and the limit.

Another feature of IPH distributions is that they can have a different tail behavior than their homogeneous counterparts, allowing for non-exponential tails. To see this, first note that if $Y \sim \mbox{IPH}(\boldsymbol{\pi} , \textbf{T} , \lambda )$ , then we may write

(2.1) \begin{equation}Y \sim g(Z) \,,\end{equation}

where $Z \sim \mbox{PH}(\boldsymbol{\pi} , \textbf{T} )$ and g is defined in terms of $\lambda$ by:

\begin{equation*}g^{-1}(y) = \int_0^y \lambda (s)ds,\quad y\ge0, \end{equation*}

which is well defined by the positivity of $\lambda$ , or equivalently (see for instance Theorem 2.9 in Albrecher and Bladt, Reference Albrecher and Bladt2019):

\begin{equation*}\lambda (s) = \frac{d}{ds}g^{-1}(s) \,. \end{equation*}

Since the asymptotic behavior of PH distributions can be deduced from its Jordan decomposition, the following result in the inhomogeneous case is available. It is a generalization of the theory developed in Albrecher et al. (Reference Albrecher, Bladt, Bladt and Yslas2021), and the proof is hence omitted.

Proposition 2.1 Let $ Y \sim\mbox{IPH}(\boldsymbol{\pi} , \textbf{T} , \lambda )$ . Then the survival function $S_Y=1-F_Y(y)$ , density $f_Y$ , hazard function $h_Y$ , and cumulative hazard function $H_Y$ of Y satisfy, respectively, as $t\to \infty$ ,

\begin{align*} S_Y(y)&= {\boldsymbol{\pi}^{\textsf{T}}}\exp \left( \int_0^y \lambda (s)ds\ \textbf{T} \right)\textbf{e} \sim c_1 [g^{-1}(y)]^{n -1} e^{-\chi [g^{-1}(y)]},\\f_Y(y) &=\lambda (y)\, {\boldsymbol{\pi}^{\textsf{T}}}\exp \left( \int_0^y \lambda (s)ds\ \textbf{T} \right)\textbf{t}\sim c_2 [g^{-1}(y)]^{n -1} e^{-\chi [g^{-1}(y)]}\lambda(y),\\ h_Y(y) &\sim c \lambda(y) ,\\ H_Y(y)& \sim k g^{-1}(y), \end{align*}

where $c_1$ , $c_2$ , c, and k are positive constants, $-\chi$ is the largest real eigenvalue of $\textbf{T}$ , and n is the dimension of the Jordan block associated with $\chi$ .

In particular, we observe that the tail behavior is determined almost entirely by function g. We also note that any model based on the intensity of the process (as will be the case for our regression model below) is asymptotically a model on the hazard of the distribution.

If g is the identity, we get Erlang tails, but for other choices, we move away from the exponential domain. The choice of g is typically determined according to the desired tail behavior required for applications. For heavy-tailed claim severities in the Fréchet max-domain of attraction, the exponential transformation $g(z)=\eta \left( \exp(z)-1\right),\:\:\eta>0$ leading to Pareto tails is the simplest choice, although the log-logistic can also have such a behavior. When data are moderately heavy-tailed, as understood by being in the Gumbel max-domain of attraction, heavy-tailed Weibull $g(z)=z^{1/\eta},\:\:\eta>0$ , or Lognormal $g(z)=\exp(z^{1/\gamma})-1,\:\:\gamma>1$ transforms can provide useful representations. In life insurance applications, the Matrix-Gompertz distribution with $g(z)=\log( \eta z + 1 ) / \eta,\:\: \eta>0$ can be particularly useful to model human lifetimes.

One tool which is particularly useful to deal with expressions involving IPH distributions is matrix calculus. Here, the main tool is Cauchy’s formula for matrices, which is as follows. If u is an analytic function and $\textbf{A}$ is a matrix, we define

\begin{align*} u( \textbf{A})=\dfrac{1}{2 \pi i} \oint_{\Gamma}u(w) (w \textbf{I} -\textbf{A} )^{-1}dw ,\quad \textbf{A}\in \mathbb{R}^{p\times p}.\end{align*}

where $\Gamma$ is a simple path enclosing the eigenvalues of $\textbf{A}$ , and $\textbf{I}$ is the p-dimensional identity matrix. This allows for functions of matrices to be well defined, which in turn gives explicit formulas for IPH distributions, simplifying calculations and numerical implementation alike.

Example 2.2 Consider the following underlying Markov process parameters:

\begin{align*}\boldsymbol{\pi}=(1,\,0,\,0)^{{\textsf{T}}},\quad\textbf{T}=\begin{pmatrix}-100\;\;\;\; & 50\;\;\;\; & 0 \\0\;\;\;\; & -1\;\;\;\; & 1/2\\0\;\;\;\; & 0\;\;\;\; & -1/100 \\\end{pmatrix},\end{align*}

and a Weibull inhomogeneity parameter $\eta=8$ , which was chosen significantly light-tailed for visualization purposes. Compared to a conventional Weibull distribution, with $\boldsymbol{\pi}=1$ , $\textbf{T}=-1/10$ , and $\eta=8$ , we see in Figure 1 that just a few additional parameters in the augmented matrix version of a common distribution can result in very different behavior in the body of the distribution (the tail behavior is guaranteed to be the same). For this structure, the case $p=3$ shows that the Markov process can traverse three states before absorption, resulting in three different modes of the density. For other matrix structures or parameters, there is in general no one-to-one correspondence between p and the number of modes of the resulting density.

Figure 1. Densities corresponding to a scalar Weibull distribution and a Matrix-Weibull distribution.

2.2 The PH regression model

We now introduce a regression model based on IPH distributions and use it to provide a segmentation model for claim severities. Classically, for pure premium calculation, the mean is the main focus. Presently, we do not aim to improve on other data-driven procedures for pure mean estimation but instead focus on the entire distribution of the loss severities. That said, a consequence of a well-fitting global model is a reasonable mean estimate in the small sample case. However, when considering large sample behavior, any scoring rule in terms of means is not a proper scoring rule for target distributions that are not fully characterized by their first moment (see Gneiting and Raftery, Reference Gneiting and Raftery2007), and thus model choice is less important in this case.

Let $\textbf{X} = (X_1, \dots, X_d)^{\textsf{T}}$ be a d-dimensional vector of rating factors associated with the loss severity random variable Y. We keep the convention that $\textbf{X}$ does not contain the entry of 1 associated with an intercept in order to obtain an identifiable regression specification (the intercept is included in the matrix $\textbf{T}$ ). Let $\boldsymbol{\beta}$ be a d-dimensional vector of regression coefficients so that $\textbf{X}^{\textsf{T}}\boldsymbol{\beta}$ is a scalar. Then, we define the conditional distribution of $Y|\textbf{X}$ as having density function, see Proposition 2.1,

(2.2) \begin{eqnarray} f(y) &=& m(\textbf{X}^{\textsf{T}}\boldsymbol{\beta}) \lambda (y;\; \theta )\, {\boldsymbol{\pi}^{\textsf{T}}}\exp \left( m\left(\textbf{X}^{\textsf{T}}\boldsymbol{\beta}\right)\int_0^y \lambda (s;\; \theta )ds\ \textbf{T} \right)\textbf{t} \,, \end{eqnarray}

where $m:\mathbb{R}\to \mathbb{R}_+$ is any positive function (monotonicity is not required). We call this the phase-type regression model. In particular, we recognize an IPH distribution $ \mbox{IPH}(\boldsymbol{\pi} , \textbf{T} , \lambda(\cdot \mid \textbf{X} , \boldsymbol{\beta}, \theta ) )$ , where

(2.3) \begin{align} \lambda(t \mid \textbf{X} , \boldsymbol{\beta},\theta )\;:=\;\lambda (t;\; \theta ) m(\textbf{X}^{\textsf{T}}\boldsymbol{\beta}) \,.\end{align}

One possible interpretation of this model is that covariates act multiplicatively on the underlying Markov intensity, thus creating proportional intensities among different policies. In fact, since the intensity of an IPH distribution is asymptotically equivalent to its hazard function, we have that covariates satisfy a generalized proportional hazards specification in the far-right tail.

The conditional mean can be written in the form:

(2.4) \begin{align}\mu(Y|\textbf{X})=\int_0^\infty \boldsymbol{\pi}^{\textsf{T}}\exp\left(m(\textbf{X}^{\textsf{T}}\boldsymbol{\beta})\int_0^y \lambda(s;\; \theta )ds\,\textbf{T} \right)\textbf{e} \,dy,\end{align}

which is simply obtained by integrating the survival function over $\mathbb{R}_+$ , see Proposition 2.1. A simple special case is obtained by the following choices, giving a Gamma GLM with canonical link: take $\textbf{T}=-1$ , and $\lambda\equiv 1$ to receive

\begin{align*}\mu(Y|\textbf{X})=\int_0^\infty \exp(-m(\textbf{X}^{\textsf{T}}\boldsymbol{\beta}) y ) \,dy=\frac{1}{ m(\textbf{X}^{\textsf{T}}\boldsymbol{\beta})}.\end{align*}

Another slightly more complex special case is that of regression for Matrix-Weibull distributions, which contains the pure PH specification (when $\lambda\equiv1$ ). In this setting, it is not hard to see that

(2.5) \begin{align}\mu(Y|\textbf{X})=\int_0^\infty \boldsymbol{\pi}^{\textsf{T}}\exp\left(m(\textbf{X}^{\textsf{T}}\boldsymbol{\beta})\textbf{T} y^{ \theta }\right)\textbf{e} \,dy =\frac{\Gamma(1+ \theta ^{-1}) \boldsymbol{\pi}^{\textsf{T}}\textbf{T}^{- \theta ^{-1}}\textbf{e}}{ m(\textbf{X}^{\textsf{T}}\boldsymbol{\beta})^{ \theta ^{-1}}}.\end{align}

Figure 2 shows (2.4) as a function of a unidimensional $\textbf{X}$ for various inhomogeneity functions $\lambda$ , and fixing the underlying Markov parameters at $\boldsymbol{\pi}=(1/2,1/4,1/4)^{{\textsf{T}}}$ and

\begin{align*}\textbf{T}=\begin{pmatrix}-1\;\;\;\; & 1/2\;\;\;\; & 1/4 \\1/4\;\;\;\; & -1\;\;\;\; & 1/4\\1/4\;\;\;\; & 1/4\;\;\;\; & -1 \\\end{pmatrix}.\end{align*}

Figure 2. Mean functions of the PH regression model, as a function of a univariate regressor and for different specifications of the inhomogeneity function $\lambda$ .

2.3 Estimation

Claim severity fitting can be done by MLE in much the same way whether there are rating factors or not, and hence we exclusively present the more general case. For this purpose, we adopt a generalization of the EM algorithm, which we outline below. Such an approach is much faster to converge than naive multivariate optimization, provided a fast matrix exponential evaluation is available (see Moler and Van Loan, Reference Moler and Van Loan1978).

A full account of the EM algorithm for ordinary PH distributions can be found in Asmussen et al. (Reference Asmussen, Nerman and Olsson1996) which led to the subsequent developments and extensions to censored and parameter-dependent transformations in Olsson (Reference Olsson1996) and Albrecher et al. (Reference Albrecher, Bladt and Yslas2020b). Further details are provided in Appendix 6, since the formulas therein are the building blocks of the EM algorithm for the PH regression model.

By defining $g(\cdot;\; \theta)$ and $g(\cdot|\textbf{X}, \boldsymbol{\beta}, \theta)$ as follows in terms of their inverse functions $g^{-1}(y;\; \theta )=\int_{0}^{y}\lambda(s;\; \theta )ds$ and

\begin{align*}g^{-1}(y \,|\, \textbf{X}, \boldsymbol{\beta}, \theta )=\int_{0}^{y}\lambda(s \,|\, \textbf{X}, \theta )ds = g^{-1}(y;\; \theta ) \exp(\textbf{X}^{\textsf{T}}\boldsymbol{\beta} ),\end{align*}

we note the simple identity:

(2.6) \begin{align} g(y \,|\, \textbf{X}, \boldsymbol{\beta},\theta )=g(y \exp(-\!\textbf{X}^{\textsf{T}}\boldsymbol{\beta} ) \,;\, \theta ) \,,\end{align}

which yields that $g^{-1}( Y \,|\, \textbf{X}, \boldsymbol{\beta}, \theta ) \sim \mbox{PH}(\boldsymbol{\pi} , \textbf{T} )$ . In other words, there is a parametric transformation, depending both on the covariates and regression coefficients, which brings the PH regression model into a conventional PH distribution. The generalized EM algorithm hinges on this fact and optimizes the transformed PH distribution and the transformation itself in alternating steps.

Algorithm 1 Generalized EM algorithm for phase-type regression

Input:: positive data points $\textbf{y}=(y_1,y_2,\dots,y_N)^{{\textsf{T}}}$ , covariates $\textbf{x}_1,\dots,\textbf{x}_N$ , and initial parameters $( \boldsymbol{\pi} , \textbf{T} ,\theta)$

1. (1) Transformation: transform the data into $z_i=g^{-1}(y_i;\; \theta ) \exp(\textbf{x}_i^{\textsf{T}} \boldsymbol{\beta})$ , $i=1,\dots,N$ .

2. (2) E-step: compute the statistics (see Appendix 6 for precise definitions of the random variables $B_k$ , $Z_k$ , $N_{ks}$ , and $N_k$ ):

   \begin{align*} \mathbb{E}(B_k\mid \textbf{Z}=\textbf{z})=\sum_{i=1}^{N} \frac{\pi_k {\textbf{e}_k}^{{\textsf{T}}}\exp( \textbf{T} z_i) \textbf{t} }{ \boldsymbol{\pi}^{\textsf{T}} \exp( \textbf{T} x_i) \textbf{t} }\end{align*}
   \begin{align*} \mathbb{E}(Z_k\mid \textbf{Z}=\textbf{z})=\sum_{i=1}^{N} \frac{\int_{0}^{z_i}{\textbf{e}_k}^{{\textsf{T}}}\exp( \textbf{T} (z_i-u)) \textbf{t} \boldsymbol{\pi}^{\textsf{T}} \exp( \textbf{T} u)\textbf{e}_kdu}{ \boldsymbol{\pi}^{\textsf{T}} \exp( \textbf{T} z_i) \textbf{t} }\end{align*}
   \begin{align*}\mathbb{E}(N_{ks}\mid \textbf{Z}=\textbf{z})=\sum_{i=1}^{N}t_{ks} \frac{\int_{0}^{z_i}{\textbf{e}_s}^{{\textsf{T}}}\exp( \textbf{T} (z_i-u)) \textbf{t} \, \boldsymbol{\pi}^{\textsf{T}} \exp( \textbf{T} u)\textbf{e}_kdu}{ \boldsymbol{\pi}^{\textsf{T}} \exp( \textbf{T} z_i) \textbf{t} }\end{align*}
   \begin{align*}\mathbb{E}(N_k\mid \textbf{Z}=\textbf{z})=\sum_{i=1}^{N} t_k\frac{ \boldsymbol{\pi}^{\textsf{T}} \exp( \textbf{T} z_i){\textbf{e}}y_k}{ \boldsymbol{\pi}^{\textsf{T}} \exp( \textbf{T} z_i) \textbf{t} }.\end{align*}

3. (3) M-step: let

   \begin{align*}\hat \pi_k=\frac{\mathbb{E}(B_k\mid \textbf{Z}=\textbf{z})}{N}, \quad \hat t_{ks}=\frac{\mathbb{E}(N_{ks}\mid \textbf{Z}=\textbf{z})}{\mathbb{E}(Z_{k}\mid \textbf{Z}=\textbf{z})}\end{align*}
   \begin{align*}\hat t_{k}=\frac{\mathbb{E}(N_{k}\mid \textbf{Z}=\textbf{z})}{\mathbb{E}(Z_{k}\mid \textbf{Z}=\textbf{z})},\quad \hat t_{kk}=-\sum_{s\neq k} \hat t_{ks}-\hat t_k.\end{align*}

4. (4) Maximize

   \begin{align*} (\hat{\theta},\hat{\boldsymbol{\beta}}) & = argmax_{(\theta ,\boldsymbol{\beta})} \sum_{i=1}^N \log (f_{Y}(y_i;\; \hat{\boldsymbol{\pi}}, \hat{\textbf{T}}, \theta ,\boldsymbol{\beta} )) \\ & = argmax_{(\theta ,\boldsymbol{\beta})} \sum_{i=1}^N \log\left(\!m(\textbf{x}_i \boldsymbol{\beta}) \lambda (y;\; \theta )\, {\boldsymbol{\pi}^{\textsf{T}}}\exp \!\left(\!m(\textbf{x}_i \boldsymbol{\beta}\!)\int_0^y \lambda (s;\; \theta )ds\ \textbf{T}\!\right)\!\textbf{t}\!\right) \end{align*}

5. (5) Update the current parameters to $({\boldsymbol{\pi}},\textbf{T}, \theta ,\boldsymbol{\beta}) =(\hat{{\boldsymbol{\pi}}},\hat{\textbf{T}}, \hat{ \theta },\hat{\boldsymbol{\beta}})$ . Return to step 1 unless a stopping rule is satisfied.

Output: fitted representation $( \boldsymbol{\pi} , \textbf{T} ,\theta, \boldsymbol{\beta})$ .

A standard fact of the above procedure, which is not hard to verify directly, is the following result.

In Asmussen et al. (Reference Asmussen, Nerman and Olsson1996), matrix exponentiation was done for PH estimation by converting the problem into a system of ordinary differential equations (ODEs) and then subsequently solved by the Runge–Kutta method of order 4. The C implementation, called EMpht, is still available online, cf. Olsson (Reference Olsson1998). Using the same ODE approach, Albrecher et al. (Reference Albrecher, Bladt, Bladt and Yslas2021) implemented a C++ based algorithm available from the matrixdist package in R (cf. Bladt and Yslas, Reference Bladt and Yslas2021). However, when extended to IPH distributions, this approach requires reorderings and lacks the necessary speed to estimate datasets of the magnitude of those found in insurance.

For this reason, we presently choose to make use of the uniformization method, which expresses the dynamics of a continuous-time Markov jump process in terms of a discrete-time chain where jumps occur at a Poisson intensity. Please check the usage of the term ‘consists on’ in the sentence ‘For this reason, we..’ and correct if necessary.More precisely, taking $\phi= \max_{k=1,\dots,p} (\!-\!t_{kk})$ and $\textbf{Q} = {\phi}^{-1}\left( \phi \textbf{I} + \textbf{T} \right)$ , which is a transition matrix, we have that

\begin{align*} \exp( \textbf{T} y)= \sum_{n=0}^{\infty} \frac{\textbf{Q} ^{n} (\phi y)^{n}}{n !} e^{-\phi y }. \end{align*}

and a truncated approximation of this series has error smaller than

\begin{align*} \sum_{n=M+1}^{\infty} \frac{(\phi y)^{n}}{n !} e^{-\phi y } = \mathbb{P}(N_{\phi y}>M) \,, \end{align*}

with $N_{\phi y} \sim \mbox{Poisson}(\phi y) $ . In other words, with the introduction of the regularization parameter $\phi$ , we are able to tame the error incurred in the truncation approximation by considering M large enough according to a Poisson law. Further improvements can be achieved by artifacts as the following one: since trivially $e^{\textbf{T} y}=(e^{\textbf{T} y/2^{m}})^{2^m}$ , then the Poisson law with mean $\phi y /2^{m} <1 $ implies good approximations with small M for $e^{\textbf{T} y/2^{m}}$ and then $e^{\textbf{T} y}$ is obtained by simple sequential squaring.

Integrals of matrix exponentials can be computed using an ingenious result of Van Loan (Reference Van Loan1978), which states that

\begin{align*} \exp \left( \left( \begin{array}{c@{\quad}c} \textbf{T} & \textbf{t} \, \boldsymbol{\pi} \\ \textbf{0} & \textbf{T} \end{array} \right) y \right)= \left( \begin{array}{c@{\quad}c} \textrm{e}^{\textbf{T} y } & \int_{0}^{y} \textrm{e}^{ \textbf{T} (y-u)} \textbf{t} \boldsymbol{\pi} \textrm{e}^{ \textbf{T} u}du \\ \textbf{0} & \textrm{e}^{\textbf{T} y } \end{array} \right)\,. \end{align*}

Hence, integrals of the form $\int_{0}^{y} \textrm{e}^{ \textbf{T} (y-u)} \textbf{t} \boldsymbol{\pi} \textrm{e}^{ \textbf{T} u}du$ can easily be computed at a stroke by one single matrix exponential evaluation of the left-hand-side expression and extracting the corresponding upper-right sub-matrix.

Bundled together with a new RcppArmadillo implementation (which is now also part of the latest Github version of matrixdist), the increase in speed was of about two to three orders of magnitude. Consequently, larger datasets with more covariates can be estimated in reasonable time. Naturally, the methods are still slower than those for simple regression models such as GLMs.

2.4 Inference and goodness of fit

Inference and goodness of fit can always be done via parametric bootstrap methods. However, re-fitting a PH regression can be too costly. Consequently, we will take a more classical approach and derive certain quantities of interest which will allow us to perform variable selection and assess the quality of a fitted model.

2.4.1 Inference for phase-type regression models

For the following calculations, we focus only in the identifiable parameters of the regression, which are the relevant ones in terms of segmentation. Furthermore, we now make the assumption that m is differentiable and write $\ell_{\textbf{y}}(\boldsymbol{\beta},\theta)$ for the joint log-likelihood function of the observed severities $\textbf{y}=(y_1,\dots,y_N)$ with corresponding rating factors $\overline{\textbf{x}}=(\textbf{x}_1,\dots,\textbf{x}_N)$ .

From (2.2), we obtain for $j = 1,\dots,d$ and $h(y;\; \theta )\;:=\;\int_0^y \lambda (s;\; \theta )ds$ ,

(2.7) \begin{align}\frac{d\ell_{\textbf{y}}(\boldsymbol{\beta},\theta)}{d\beta_j}&= \sum_{i=1}^N G_1(i,j|\boldsymbol{\pi},\textbf{T},\boldsymbol{\beta},\theta,\textbf{y},\overline{\textbf{x}}),\end{align}

(2.8) \begin{align}\frac{d\ell_{\textbf{y}}(\boldsymbol{\beta},\theta)}{d \theta }&= \sum_{i=1}^N G_2(i|\boldsymbol{\pi},\textbf{T},\boldsymbol{\beta},\theta,\textbf{y},\overline{\textbf{x}}),\end{align}

where

\begin{align*}& G_1(i,j|\boldsymbol{\pi},\textbf{T},\boldsymbol{\beta},\theta,\textbf{y},\overline{\textbf{x}})=\\[5pt] & \quad x_{ij}m'\left(\textbf{x}_i^{\textsf{T}}\boldsymbol{\beta} \right) \left(\frac{1}{m\left(\textbf{x}_i^{\textsf{T}}\boldsymbol{\beta} \right) }+\frac{{\boldsymbol{\pi}^{\textsf{T}}}\exp \left( m\left(\textbf{x}_i^{\textsf{T}}\boldsymbol{\beta} \right)h(y;\; \theta )\textbf{T} \right)h(y_i;\; \theta )\textbf{T}\textbf{t}}{{\boldsymbol{\pi}^{\textsf{T}}}\exp \left( m\left(\textbf{x}_i^{\textsf{T}}\boldsymbol{\beta} \right)h(y_i;\; \theta )\ \textbf{T} \right)\textbf{t}}\right),\\[5pt]& G_2(i|\boldsymbol{\pi},\textbf{T},\boldsymbol{\beta},\theta,\textbf{y},\overline{\textbf{x}})=\\[5pt] & \quad \frac{\frac{d}{d \theta }\lambda(y_i;\; \theta )}{\lambda(y_i;\; \theta )}+\frac{{\boldsymbol{\pi}^{\textsf{T}}}\exp \left( m\left(\textbf{x}_i^{\textsf{T}}\boldsymbol{\beta} \right)h(y_i, \theta ) \textbf{T} \right)m\left(\textbf{x}_i^{\textsf{T}}\boldsymbol{\beta})\frac{d}{d \theta }h(y_i, \theta \right)\textbf{T}\textbf{t}}{{\boldsymbol{\pi}^{\textsf{T}}}\exp \left( m\left(\textbf{x}_i^{\textsf{T}}\boldsymbol{\beta} \right)h(y_i;\; \theta ) \textbf{T} \right)\textbf{t}}\end{align*}

Such expressions should be equal to 0 whenever convergence is achieved in step 2 of Algorithm 1.

In general, the score functions (2.7) and (2.8) will not have an explicit solution when equated to 0. Concerning the $(d+1)\times(d+1)$ -dimensional Fisher information matrix, we may write it as:

(2.9) \begin{align}[\mathcal{I}]_{jk}&= \begin{cases}\! \sum_{i=1}^N G_1(i,j|\boldsymbol{\pi},\textbf{T},\boldsymbol{\beta},\theta,\textbf{y},\overline{\textbf{x}}) G_1(i,k|\boldsymbol{\pi},\textbf{T},\boldsymbol{\beta},\theta,\textbf{y},\overline{\textbf{x}}) & \!\!\! 1\le j,k \le d,\\ \sum_{i=1}^N G_1(i,j|\boldsymbol{\pi},\textbf{T},\boldsymbol{\beta},\theta,\textbf{y},\overline{\textbf{x}}) G_2(i|\boldsymbol{\pi},\textbf{T},\boldsymbol{\beta},\theta,\textbf{y},\overline{\textbf{x}}) & \!\!\! 1\le j \le d,\;k=d+1,\\ \sum_{i=1}^N G_2^2(i|\boldsymbol{\pi},\textbf{T},\boldsymbol{\beta},\theta,\textbf{y},\overline{\textbf{x}}) & \!\!\! j=k=d+1.\end{cases}\end{align}

Then, we may use the asymptotic MLE approximation, as $N\to\infty$ ,

\begin{align*}\widehat{(\boldsymbol{\beta},\theta)}\stackrel{d}{\approx}\mathcal{N}((\boldsymbol{\beta},\theta),\mathcal{I}^{-1}),\end{align*}

which can be used to compute standard errors and Neyman-type confidence intervals. For instance,

\begin{align*}(\hat\beta_1-1.96 \sqrt{[\mathcal{I}^{-1}]_{11}},\:\hat\beta_1+1.96 \sqrt{[\mathcal{I}^{-1}]_{11}}),\end{align*}

constitutes an approximate 95% confidence interval for the parameter $\beta_1$ . We obtain p-values in much the same manner.

Observe that we have deliberately omitted inference for the PH parameters $\boldsymbol{\pi}$ and $\textbf{T}$ . The reason being, this is a particularly difficult task due to the non-identifiability issue of PH distributions, namely that several representations can result in the same model. Papers such as Bladt et al. (Reference Bladt, Esparza and Nielsen2011) and Zhang et al., Reference Zhang, Zheng, Okamura and Dohi2021 provide methods of recovering the information matrix for these parameters, but the validity of the conclusions drawn from such approach are not fully understood. Consequently, and similarly to other regression models, we will only perform estimation on the regression coefficients themselves and use the distributional parameters merely as a vehicle to obtain $\boldsymbol{\beta}$ .

Remark 2.2 We performed estimation using numerical optimization, with and without using the gradient of the log-likelihood for step 2 in Algorithm 1. Not using the gradient was always faster, since when dealing with matrix exponentials the evaluation of the derivatives is equally costly as the objective function itself, and thus best avoided. In practice, the Fisher information matrix may be near-singular when large numbers of covariates are considered. The negative of the the Hessian matrix obtained by numerical methods such as the Broyden-Fletcher-Goldfarb-Shanno algorithm is a more reliable alternative in this case, whose inverse may also be used to approximate the asymptotic covariance matrix.

2.4.2 Goodness of fit for phase-type regression models

Below is a brief description of a goodness-of-fit diagnostic tool for the above models, which, when ordered and plotted against theoretical uniform order statistics, constitute a PP plot. The idea is to transform the covariate-specific data into a parameter-free scale using the probability integral transformation (PIT).

Indeed, by applying the PIT transform to the data:

(2.10) \begin{align}r_i=\left({\boldsymbol{\pi}^{\textsf{T}}}\exp \left(m\left(\textbf{x}_i^{\textsf{T}} \boldsymbol{\beta} \right) \int_0^{y_i} \lambda (s;\; \theta )ds\ \textbf{T} \right)\textbf{e}\right),\quad i=1,\dots N.\end{align}

The null hypothesis under the PH regression model is

\begin{align*}\mathcal{H}_0:Y_i\sim \mbox{IPH}\left({\boldsymbol{\pi}},\textbf{T},m\left(\textbf{x}_i^{\textsf{T}}\boldsymbol{\beta}\right)\lambda\right),\quad i=1,\dots,N,\end{align*}

and it standard see that the dataset:

\begin{align*}\mathcal{D}=\{r_1,\:r_2,\dots, r_N\}\end{align*}

constitutes a sample of uniform variables.

Thus, it only remains to assess the goodness of fit of $\mathcal{D}$ against a standard uniform distribution, for which a suite of statistical tests and visual tools can be applied.

2.5 On the choice of matrix dimension and structure

Often a general sub-intensity matrix structure is too general to be either practical or parsimonious, and special structures can do the job almost as well with fewer parameters (see, e.g., Bladt and Nielsen, Reference Bladt and Nielsen2017, Section 3.1.5, and also Section 8.3.2). This is because IPH distributions are not identifiable: two matrices of different dimension and/or structure may yield exactly the same distribution.

A zero in the off-diagonal of the sub-intensity matrix $\textbf{T}$ indicates absence of jumps between two states, whereas a zero in the diagonal is not possible. Below we describe the most commonly used special structures for $\textbf{T}$ (and its corresponding $\boldsymbol{\pi}$ ). Note that we borrow the names of the structures from the homogeneous case, but the resulting distribution will in general not be linked to such name. For instance, a Matrix-Pareto of exponential structure will yield a Pareto distribution, and not an exponential one.

Exponential structure. This is the simplest structure one can think of and applies only for $p=1$ , that is, the scalar case. We have $\boldsymbol{\pi}^{\textsf{T}}=1$ , $\textbf{T}=-\alpha<0$ , and so $\textbf{t}=\alpha$ . The corresponding stochastic process is schematically depicted in the top-left panel of Figure 3.

Figure 3. Underlying Markov structures. Names are borrowed from the corresponding PH representations but apply to our inhomogeneous setup as well. The state 0 is added for schematic reasons but is not part of the actual state space of the chain. The (F) general case has the intensities $t_{ij}$ and $t_{ji}$ between each pair of states $i,j\in\{0,\dots,p\}$ omitted for display purposes.

Erlang structure. This structure traverses p identical states consecutively from beginning to end, or in matrix notation:

\[ \boldsymbol{\pi}^{\textsf{T}} =(1, \:0, \dots, \:0 ),\]

\[ \textbf{T} = \begin{pmatrix} -\alpha\;\;\;\; & \alpha\;\;\;\; & \;\;\;\; & 0 \\ \;\;\;\; & \ddots\;\;\;\; &\ddots\;\;\;\; & \vdots\\ \;\;\;\; & \;\;\;\; & -\alpha\;\;\;\; & \alpha\\ 0\;\;\;\; & \;\;\;\; & \;\;\;\; &-\alpha \end{pmatrix}, \, \quad \textbf{t}= \begin{pmatrix}0 \\ \vdots \\ 0 \\\alpha \end{pmatrix},\quad \alpha>0. \]

The corresponding stochastic process is schematically depicted in the top-right panel of Figure 3.

Hyperexponential PH distribution. This structure will always give rise to a mixture of scalar distributions, since it starts in one of the p states but has no transitions between them thereafter. The matrix representation is the following:

\[ \boldsymbol{\pi}^{\textsf{T}} =(\pi_1, \: \pi_2, \dots, \: \pi_p),\]

\[ \textbf{T} = \begin{pmatrix} -\alpha_1\;\;\;\; & \;\;\;\; & \;\;\;\; & 0 \\ \;\;\;\; & \;\;\;\; & \ddots\;\;\;\; & \\ 0\;\;\;\; & \;\;\;\; & \;\;\;\; & -\alpha_p \end{pmatrix}, \, \quad \textbf{t} = \begin{pmatrix}\alpha_1 \\ \vdots \\\alpha_p \end{pmatrix},\quad \alpha_i>0, \:\:\: i=1,\dots,p. \]

The corresponding stochastic process is schematically depicted in the center-left panel of Figure 3.

Coxian structure. This structure is a generalization of the Erlang one, where different parameters are allowed for the diagonal entries, and also absorption can happen spontaneously from any intermediate state. The parametrization is as follows:

\[ \boldsymbol{\pi}^{\textsf{T}} =(1,\: 0, \dots,\: 0),\]

\[ \textbf{T} = \begin{pmatrix} -\alpha_1\;\;\;\; & q_1\alpha_1\;\;\;\; & \;\;\;\; & 0 \\ \;\;\;\; & -\alpha_2\;\;\;\; & q_2\alpha_2\;\;\;\; & \\ \;\;\;\; & \;\;\;\; & \ddots \ddots\;\;\;\; & \\\;\;\;\; & \;\;\;\; & -\alpha_{p-1}\;\;\;\; & q_{p-1}\alpha_{p-1} \\ 0 & & & -\alpha_p \end{pmatrix}, \quad \textbf{t} = \begin{pmatrix}(1-q_1)\alpha_1 \\ (1-q_{2})\alpha_{2} \\ \vdots \\ (1-q_{p-1})\alpha_{p-1} \\\alpha_p \end{pmatrix},\]

with $\alpha_i>0,\:q_i\in[0,1], \:i=1,\dots,p$ . The corresponding stochastic process is schematically depicted in the center-right panel of Figure 3.

Generalized Coxian structure. This structure has the same $\textbf{T}$ and $\textbf{t}$ parameters as in the Coxian case, the only difference being that the initial vector is now distributed among all p states, that is,

\[ \boldsymbol{\pi}^{\textsf{T}} =(\pi_1, \: \pi_2,\dots, \:\pi_p ).\]

The corresponding stochastic process is schematically depicted in the bottom-left panel of Figure 3.

Hence, when selecting between fitted dimensions and structures, it is advised to keep track of the negative log-likelihood (rather than AIC or BIC) and its decrease when increasing the number of phases or changing the structure. Theoretically, larger matrices will always yield better likelihood, but in practice, it is often the case that the likelihood stops increasing (effectively getting stuck in a good local optimum), or increases are not so significant. Usually, EM algorithms for PH distributions for marginal distributions estimate effectively and in reasonable time up to about 20–30 phases, while for a PH regression the number is closer to 5–10. The additional inhomogeneity function of IPH distributions, however, makes models more parsimonious than the PH counterparts, and usually less than 10 phases are needed.

Once the dimension and structure is chosen, one can perform model selection with respect to the regression coefficients. These fall into the usual framework, and the AIC or BIC criteria can then be used to compare and select between various models.

4.1 Marginal estimation

In the first step, we estimate the marginal behavior of the data, without rating factors. The model we choose to estimate is an IPH distribution with the exponential transform, namely, a Matrix-Pareto distribution. This choice is motivated by a preliminary analysis of the data by which a Hill plot indicated the presence of regularly varying tails, indicating heavy-tailedness. Since the data on the log-scale have a pronounced hump in the middle, a dimension larger than 1 or 2 is needed, so we chose 5. Finally, since there are no particular specificities in the right or left tail, a Coxian sub-matrix structure can do an equally good job as a general structure, and thus we select the former in the name of parsimony.

We also estimate a 20-dimensional Matrix-Weibull distribution on the data to illustrate the denseness property of IPH distributions. The Matrix-Weibull distribution generates modes more easily than the Matrix-Pareto, and we choose such a high dimension in order to illustrate the denseness of IPH distributions. However, this model is over-parametrized and thus mainly of academic interest.

As the first reference, we also fit a two-parameter Gamma distribution to the data. We also consider a splicing (or composite) model, cf. Albrecher et al. (Reference Albrecher, Beirlant and Teugels2017) and Reynkens et al. (Reference Reynkens, Verbelen, Beirlant and Antonio2017) (see also Miljkovic and Grün, Reference Miljkovic and Grün2016), defined as follows. Let $f_1$ and $\:f_2$ be two densities supported on the positive real line. A spliced density based on the the latter components is then given by:

\begin{align*}f_S(x)=v \frac{f_1(x)}{F_1(t)} 1_{(0,t]}(x) + (1-v) \frac{f_2(x)}{1- F_2(t)} 1_{(t,\infty)}(x), \quad t>0,\:\: v\in[0,1],\end{align*}

where t is the splicing location. We consider a popular implementation found in the ReIns package in R, Reynkens and Verbelen (Reference Reynkens and Verbelen2020), where

\begin{align*}f_1(x)=\sum_{j=1}^{M} \omega_{j} \frac{x^{r_{j}-1} \exp (-x / \zeta)}{\zeta^{r_{j}}\left(r_{j}-1\right) !},\:\:x>0,\quad f_2(x)=\frac{1}{\nu t}\left(\frac{x}{t}\right)^{-\frac{1}{\nu}-1},\:\:x>t,\end{align*}

are, respectively, a ME densities with common scale parameter and a Pareto density, which we call the ME-Pareto specification.

The parameter estimates are the following for the Gamma distribution, found by MLE,

\begin{align*}\text{Gamma shape}= 0.75, \quad \text{Gamma scale} = 2925.67,\end{align*}

while for the splicing model we obtain, selecting $\nu$ through EVT techniques (Hill estimator) and the remaining parameters through MLE (using an EM algorithm), and choosing M from 1 to 10 according to the best AIC:

\begin{align*}&M=5, \quad \zeta=363.7,\quad \textbf{r}=(1,\, 2 ,\, 4,\, 11,\, 21),\\ &\boldsymbol{\omega}=(0.26,\, 0.07,\, 0.53,\, 0.09,\, 0.02),\quad \nu =0.56,\quad v=0.96,\end{align*}

and $t\approx10,000$ was selected by visual inspection according to when the Hill plot stabilizes for the sample. Finally, for the IPH Matrix-Pareto distribution, we obtain the following estimates by applying Algorithm 1 without covariates (which is implemented in the matrixdist package in R, Bladt and Yslas, Reference Bladt and Yslas2021):

\begin{align*}\boldsymbol{\pi}=\textbf{e}_1, \quad \textbf{T}= \begin{pmatrix}-12.61\;\;\;\; & 12.48\;\;\;\; & 0\;\;\;\; & 0\;\;\;\; & 0\\0\;\;\;\; & -12.61\;\;\;\; & 10.33\;\;\;\; & 0\;\;\;\; & 0\\0\;\;\;\; & 0\;\;\;\; & -1.99\;\;\;\; & 1.99\;\;\;\; & 0\\0\;\;\;\; & 0\;\;\;\; & 0\;\;\;\; & -7.34\;\;\;\; & 7.34\\0\;\;\;\; & 0\;\;\;\; & 0\;\;\;\; & 0\;\;\;\; & -7.34\\\end{pmatrix},\quad \eta = 1149.57,\end{align*}

where $\textbf{e}_1=(1,\,0,\,0,\,0,\,0)^{\textsf{T}}.$ The 404 fitted parameters of the Matrix-Weibull, most of the zero, are omitted, since they do not provide further intuition on IPH fitting.

Table 2 provides numerical results on the three fitted models, and in particular it shows that an IPH Matrix-Pareto distribution is well aligned with the state-of-the-art models for severity modeling in insurance such as composite methods. The AIC and BIC are not good measures for IPH distributions, since they can over-penalize them when the underlying sub-intensity matrix is not of minimal order. However, we see that even with this hinderance, the performance in terms of information criteria is positive.

Table 2. Summary of fitted marginal models to severities from the freMPL dataset.

The smallest AIC and BIC values are given in boldface.

The quality of the fit is also assessed in Figure 6. In particular, we observe that the data have a big spike at the log-severity around the value $7.3$ , which poses problems to the first three distributions. This disturbance manifests itself in the PP plot as a slight “S” shape. Not surprisingly, the visual quality of the fit of the splicing and IPH models is, as expected, far superior to its Gamma counterpart, while the difference between the two Pareto-tailed models is not very noticeable. Their implied tail indices agree, as given respectively by the formulas $\hat\xi=\nu=0.56,$ and

\begin{align*}\hat\xi = -1/\max\{\Re\text{Eigen}(\hat{\textbf{T}})\} = 0.50,\end{align*}

where $\Re\text{Eigen}(\textbf{A})$ denotes the set of real eigenvalues of a matrix $\textbf{A}$ . For reference, recall that in terms of the Pareto tail parameter, the following relationship holds $\alpha = 1/\xi$ . As a sanity check, we see that both estimates for the tail index fall perfectly within the Hill estimator’s bounds. In particular, the estimated models have finite mean (their tail index is less than 1; alternatively, their Pareto index is larger than 1). We would like to remark that fitting a conventional Pareto distribution as a global model will result in a particularly bad fit, such that in the IPH case we have managed use matrix parameters to our advantage, in order to automatically engineer an all-inclusive model, featuring no threshold selection.

Figure 6. Fitted distributions to claim severities from the MPL dataset. The log-transform is used exclusively for visual purposes, the estimation having been carried out in the usual scale. KS refers to the Kolmogorv–Smirnov statistic for testing uniformity.

4.2 Incorporating rating factors

Having already analyzed the marginal behavior of the French MPL data, we proceed to incorporate rating factors through the PH regression model that we have presented above. In Table 3, we show summarizing statistics for claim severity, while two further tables of summarizing statistics for continuous rating factors and categorical rating factors, respectively, can be found in Supplementary Material.

Table 3. Summary statistics of the French MPL dataset: claim severity.

To reduce the large number of categories and variables, we perform a preprocessing step where we individually, for each categorical variable, prune a regression tree according to the best predictive performance on the mean of the log-severities (since the variables exhibit heavy-tails). We then merge categories together according to the resulting optimal tree. For all the continuous variables, we perform a simple shift-and-scale transformation.

The Matrix-Pareto is the only distribution which according to the Hill estimator will have the correct tail behavior, and thus we consider it as a good target model for PH regression. However, we also consider a model which does not have the correct tail behavior, namely the Matrix-Weibull distribution. The latter is convenient since it has a closed-form formula for mean prediction as per Equation (2.5).

The result of the fit using all the processed rating factors is given in Table 4 and Figure 7, with all the coefficients of the five-dimensional IPH models multiplied by $-1$ , to be comparable with the log-link Gamma GLM. The significance is obtained from a normal approximation using the implied Fisher matrix based on the Hessian matrix from the numerical optimization, as detailed in the goodness-of-fit subsection above. When considering less covariates, Equation (2.9) can be equivalently used, yielding very close results. The AIC and BIC are well below what the corresponding GLM can achieve, even if this is not a generous metric for IPH models.

Table 4. Summary for GLM and PH regression models for the freMPL dataset.

The smallest AIC and BIC values are given in boldface.

Figure 7. Coefficients and p-values of IPH and GLM regression. For display: IPH coefficients multiplied by $-1$ and intercept of GLM omitted.

We observe that the loss ratio (where we understand as losses the claim severities, i.e., disregarding frequencies) is kept much better for the Weibull variant than for the Pareto, despite the latter being a better model for large claims. In practice, this inconvenience can be amended by regressing the losses with respect to the premium. Since estimation is done via MLE, it is not uncommon for a model with misspecified tail behavior to globally perform better. The fact that the Weibull PH regression performs almost as well as the GLM when predicting the mean is remarkable, since the former model does not specifically target averages, and the latter does.

In Figure 8, we illustrate the behavior of the aggregation of the implied mean predictions across all policies. We consider prediction from all three models, for selected rating factors and their categories, as computed by the general formula (2.4) in the Pareto case, and the explicit formula (2.5) in the Weibull case. Each of the aggregate predictions are normalized to sum to 1 for display purposes. We again see that despite the Weibull PH regression not being specifically designed for such task, it comes very close to the performance of the GLM.

Figure 8. Left panels: aggregate observed losses versus aggregate implied model premia (expected value), normalized to sum to 1, for the Gamma GLM, and Pareto and Weibull PH regressions; right panels: number of claims within each category (right).

In Figure 9, we illustrate the prediction of quantiles, as implied by the fitted conditional distribution functions of each model. We average the quantile estimates across policies corresponding to each of the rating factors and their categories. We observe that matrix distributions perform much better in this respect, and the correct specification of the Pareto tail can be appreciated in the bottom right panel, corresponding to the quantile at level $0.9$ . Figure 10 further supports the global estimation improvement in terms of the PP plot goodness-of-fit diagnostic.

Figure 9. Empirical quantiles by coverage category versus mean (accoss all other covariates) quantiles implied by the Gamma GLM, and Pareto and Weibull PH regressions.

Figure 10. Ordered PITs from Equation (2.10) versus uniform order statistics for the French MPL dataset. KS refers to the Kolmogorv–Smirnov statistic for testing uniformity.

In terms of out-of-sample sample predictive performance, the mean square error (MSE) when holding out the first 20% of the observations and training on the remaining 80% is given by ^{Footnote 1} $3.66$ , $3.69$ , and $3.70$ , for the Gamma GLM, Pareto, and Weibull PH regression models, respectively. In general, if the risk manager is solely interested in mean prediction, our experiments suggest that other data-driven methods may outperform PH regression models, although the advantage is usually small. However, the tables turn when considering the entire distribution of loss severities. Thus, a key takeaway here is that PH regression models should be used when the aim is to capture the entire distributional properties of claim severities correctly. In contrast, models specifically designed for mean prediction may be more appropriate if only the average is of interest.