Inhomogeneous phase--type distributions and heavy tails

TL;DR

This paper introduces inhomogeneous phase-type distributions with heavy tails, providing a flexible framework for modeling heavy-tail phenomena and demonstrating its effectiveness on real-world insurance data.

Contribution

It extends phase-type distributions to inhomogeneous cases, enabling better modeling of heavy tails with a unified probabilistic approach.

Findings

The new distributions can model heavy tails effectively.

Application to fire insurance data shows practical utility.

The approach retains the fitting properties of traditional PH distributions.

Abstract

We extend the construction principle of phase-type (PH) distributions to allow for inhomogeneous transition rates and show that this naturally leads to direct probabilistic descriptions of certain transformations of PH distributions. In particular, the resulting matrix distributions enable to carry over fitting properties of PH distributions to distributions with heavy tails, providing a general modelling framework for heavy-tail phenomena. We also illustrate the versatility and parsimony of the proposed approach for the modelling of a real-world heavy-tailed fire insurance dataset.

Equations

$$\mathbf{\Lambda }(t)=\left(\begin{array}{cc}\mathbf{T}(t)& \mathbf{t}(t)\\ 0& 0\end{array}\right).$$

$$\tau =\inf \{t\ge 0:X_t=p+1\}$$

$$p_{ij}(s,t)=\text{mathds}P(X_t=j|X_s=i),$$

$$\mathbf{P}(s,t)=\prod _s^t(\mathbf{I}+\mathbf{\Lambda }(u)du),$$

$$\prod _s^t(\mathbf{I}+\mathbf{\Lambda }(u)du)=\mathbf{I}+\sum _{k=1}^{\infty }{\int }_s^t{\int }_s^{u_k}\cdots {\int }_s^{u_2}\mathbf{\Lambda }\left(u_1\right)\cdots \mathbf{\Lambda }\left(u_k\right)d{u}_1\cdots \text{d}{u}_k.$$

$$\mathbf{P}(s,t)=\left(\begin{array}{cc}\prod _s^t(\mathbf{I}+\mathbf{T}(u)du)& \mathbf{e}-\prod _s^t(\mathbf{I}+\mathbf{T}(u)du)\mathbf{e}\\ 0& 1\end{array}\right),$$

$$f(x)$$

$$F(x)$$

$$\text{mathds}P(X_t=j)=\mathbf{\pi }\prod _0^t(\mathbf{I}+\mathbf{T}(u)du){\mathbf{e}}_j.$$

$$f(x)dx$$

$$f(x)=\mathbf{\pi }\prod _0^x(\mathbf{I}+\mathbf{T}(u)du)\mathbf{t}(x).$$

$$1-F(x)$$

$$\text{mathds}P(\tau >s+t|\tau >s)=\frac{\mathbf{\pi }{\prod }_0^s(1+\mathbf{T}(u)du)}{\mathbf{\pi }{\prod }_0^s(1+\mathbf{T}(u)du)\mathbf{e}}\prod _s^t(1+\mathbf{T}(u)du)\mathbf{e}$$

$$\tau -s|\{\tau >s\}\sim \text{mbox}IPH\left(\mathbf{\alpha },\mathbf{S}(\cdot )\right),$$

$$\mathbf{\alpha }=\frac{\mathbf{\pi }{\prod }_0^s(1+\mathbf{T}(u)du)}{\mathbf{\pi }{\prod }_0^s(1+\mathbf{T}(u)du)\mathbf{e}}.$$

$$f(x)$$

$$F(x)$$

$${\int }_s^t{\int }_s^{u_k}\cdots {\int }_s^{u_2}\mathbf{T}\left(u_k\right)\cdots \mathbf{T}\left(u_1\right)d{u}_1\cdots d{u}_k$$

$$\prod _s^t(\mathbf{I}+\mathbf{T}(u)du)=\exp \left({\int }_s^t\mathbf{T}(u)du\right)$$

$$\left(\begin{array}{cc}\mathbf{S}(u)& 0\\ 0& \mathbf{T}(u)\end{array}\right)$$

$$f(x)$$

$$F(x)$$

$$\mathbf{t}(t)=-\mathbf{T}(t)\mathbf{e}=-\lambda (t)\mathbf{T}\mathbf{e}=\lambda (t)\mathbf{t}.$$

$$\frac{d}{dt}{\left({\int }_s^t\mathbf{T}(u)du\right)}^k=\frac{1}{k}{\left({\int }_s^t\mathbf{T}(u)du\right)}^{k-1}\mathbf{T}(t),$$

$${\bar{F}}_Y(y)$$

$$f_Y(y)=\frac{d}{dy}{F}_Y(y)=\mathbf{\pi }{e}^{\mathbf{T}{g}^{-1}(y)}\mathbf{t}\frac{1}{g^{\prime }(g^{-1}(y))}.$$

$$g^{-1}(x)={\int }_0^x\lambda (t)dt.$$

$$\tau \stackrel{d}{=}g(X),$$

$$\text{mathds}P(Y>y)$$

$$h(\mathbf{A})=\frac{1}{2\pi \text{mbox}i}{\oint }_{\gamma }h(z)(z\mathbf{I}-\mathbf{A}{)}^{-1}dz,$$

$$g^{-1}(x)={\int }_0^x\lambda (t)dt.$$

$$\text{mathds}E(\tau )=\mathbf{\pi }{L}_g(-\mathbf{T})\mathbf{t}.$$

Full text

1. Introduction

A phase–type distribution (PH) is the distribution of the time until absorption in a Markov jump process with finitely many states of which one is absorbing and the remaining being transient. PH distributions have a long history in Applied Probability (see Bladt and Nielsen (2017) and references therein), especially in areas such as insurance risk and queueing theory, where under PH assumptions exact solutions and explicit formulae can often be derived even in rather complex models (see e.g. Asmussen (2003); Asmussen and Albrecher (2010)). Though the class of phase–type distributions is dense in the class of distributions on the positive real line (in the sense of weak convergence), by construction PH tails are always light (exponential decay), and the latter property is a main concern regarding applications where heavy tails are present (such as large claims in insurance).

In Bladt et al. (2015) the class of PH distributions was extended to allow for infinite-dimensional matrices which led to a class of distributions with phase–type like properties and a genuinely heavy tail, and the estimation of such distributions was considered in Bladt and Rojas-Nandayapa (2018).

In this paper we take another approach to define a tractable and dense class of heavy-tailed distributions. The idea is to transform PH distributions into heavy-tailed ones, but rather than transforming the PH random variable directly, we transform the time scales of each state of the underlying Markov process, leading to time–dependent jump intensities. This time-scaling allows to carry over some of the computational tools and advantages of the PH class outside of the latter and leads to a probabilistic description of classes of heavy-tailed distributions akin to the one for PH distributions. Compared to the approach in Bladt et al. (2015), the present approach has the advantage of preserving the finite dimensionality, and offers flexibility and transparency in the choice of tail type. A distinctive disadvantage is that renewal theory and, consequently, ladder height techniques break down due to the inhomogeneity of time. However, we show that for certain transformations calculations are still explicit, and the flexibility of PH distributions can in this way be more efficiently transferred to heavy-tailed distributions than with previous approaches. Moreover, the resulting new class of matrix distributions is again dense in the class of all distributions on the positive half-line, which paves the way for using it in model fitting.

From a modelling perspective, the approach proposed in this paper has some attractive features. Let us recall that a particular strength of the class of PH distributions is that it extends the favourable properties of an exponential random variable (probabilistic ones such as lack of memory as well as analytic ones such as the simplicity of calculations due to the elementary and paramount role of the exponential function in real and complex analysis) to cases where the exponential random variable itself is not a good model. It does so by concatenating these properties in a Markov jump process framework, leading to matrix expressions, but still maintaining those favourable properties. The denseness of the PH class then in principle allows to approximate any random variable and its role in the respective model with these techniques. On the other hand, Pareto random variables with their heavy tails can be motivated in various ways for particular modelling applications. One may for instance interpret them in terms of scaling phenomena in nature (like with considerations of Zipf’s law). One can also simply view them as model candidates with power tails, which decay slower than exponential. Indeed, in various applications such a tail behaviour can be observed in data. Alternatively, one may also just see them as (possibly rescaled and shifted) exponentials of exponential random variables, which inherits some attractive mathematical properties, like a scaled version of lack of memory (and the latter is one reason for the particular attractiveness of Pareto random variables as modelling tools in certain applications, like reinsurance). At the same time, a pure Pareto random variable is often not a good fit to data (particularly when the entire range is considered). One way to deal with this in practice is to use spliced models (also sometimes called composite models), where one type of model is used for smaller values and a Pareto model is used for the tail (see e.g. Pigeon and Denuit (2011); Reynkens et al. (2017); Scollnik (2007)).

A number of known distributions can be obtained as the distribution of a transformed exponential random variable. This is for example indeed the case for the Pareto, the Weibull and the generalized extreme value distribution (GEV). The obvious idea is to propose a matrix version of these distributions by transforming a phase–type distributed random variable with the same transformations. We show that the resulting distributions have similar forms and properties as their original counterparts, and that they may be seen as inhomogeneous phase–type distributions. This provides a unifying framework for seeing transformed phase–type distributed random variables as absorption times.

Transformed phase–type random variables have e.g. been treated in Ahn et al. (2012), where the logarithmic transformation is used. This results in a Pareto type of tail behaviour. For other transformations we may obtain tail behaviours of other types like e.g. Weibull or GEV. The class of transformed distributions of a certain type is again dense in the class of distributions on the positive reals. In fact we will provide different sub–classes generated by mixtures of transformed Erlang distributions which are dense and can be written in an explicit and simple form. In principle a classical phase–type distribution may approximate any heavy–tailed positive distribution arbitrarily well, but since they are light–tailed, the approximating distribution will always have a distinct tail which can become an important issue if applied to situations where the tail behaviour matters (like e.g. ruin probabilities). If one instead fits a transformed phase–type distribution to data with the ”correct” tail, then not only will one be able to capture the proper tail behaviour but it will also lead to more parsimonious models (with fewer components) than fitting with the traditional PH class. On the other hand the transformed phase–type distributions may have computational advantages (and possibly methodological advantages in terms of interpretability) to other competitors like random variables with more general regularly varying tails.

From a statistical point of view one can apply standard techniques from phase–type fitting to the transformed random variables. This is for instance the case for the log–transformed phase–type distributions where one simply fits a phase–type distribution to the logarithm of the data points. In Ahn et al. (2012), a log–phase–type distribution was fitted to Danish fire insurance data. We provide an example using Dutch insurance data, where the main purpose is a comparison of our approach to a previously obtained model in Albrecher et al. (2017) which was based on splicing Erlang and Pareto distributions. Since the splicing components are in some sense implicit in the log–Phase–type setup, it is interesting to see whether we shall be able to retrieve parts of that model (and in a more automated way).

While our estimation method is based on maximum likelihood which treats each data point with the same weight, it may be interesting in future research to extend the respective statistical techniques to incorporate extreme value analysis considerations for the fitting procedure of this type of models, with more weight being given to larger data points.

The paper is organized as follows. In Section 2 we derive various theoretical properties of inhomogeneous PH distributions. On the basis of this representation, Section 3 gives various properties for the class of matrix-Pareto distributions and provides some details for several special cases. Subsequently, Section 4 considers other transformations of PH distributions and studies some properties for the resulting distribution classes. In Section 5 we then discuss in more detail the modelling dimension of matrix-Pareto distributions and illustrate its use for a set of Dutch fire insurance data that was recently studied by other statistical techniques in Albrecher et al. (2017). Finally, Section 6 concludes.

2. Inhomogeneous Phase–type distributions

Consider a time–inhomogeneous Markov jump process (Goodman and Johansen (1973)) $\{X_{t}\}_{t\geq 0}$ on the finite state–space $E=\{1,2,...,p,p+1\}$, where states $1,2,...,p$ are transient and $p+1$ absorbing. Thus the intensity matrix of $\{X_{t}\}_{t\geq 0}$ is of the form

[TABLE]

Here, for any time $t\geq 0$, $\bm{t}(t)=-\bm{T}(t)\bm{e}$, where $\bm{e}$ denotes the column vector $\bm{e}=(1,1,...,1)^{\prime}$. Assume that $\mathds{P}(X_{0}=p+1)=0$ and define $\bm{\pi}=(\pi_{1},...,\pi_{p})$.

Definition 2.1**.**

Let

[TABLE]

denote the time until absorption of $\{X_{t}\}$. The distribution of $\tau$ is then said to be an inhomogeneous phase–type distribution with representation $(\bm{\pi},\bm{T}(t))$ and we write $\tau\sim\mbox{IPH}(\bm{\pi},\bm{T}(t))$.

The transition matrix $\bm{P}(s,t)=\{p_{ij}(s,t)\}$, where

[TABLE]

is related to the intensity matrix in terms of the product integral (see Gill and Johansen (1990); Slavík (2007))

[TABLE]

which is defined by

[TABLE]

It is then straightforward to see that

[TABLE]

where $\bm{e}=(1,1,...,1)^{\prime}$.

Theorem 2.2**.**

Assume that $\tau\sim\mbox{IPH}(\bm{\pi},\bm{T}(t))$. Then the density $f$ and the distribution function $F$ of $\tau$ are given by

[TABLE]

Proof.

Since $(\bm{\pi},0)\bm{P}(0,t)$ is the distribution of $X_{t}$, we get from (1) that

[TABLE]

The density $f$ for $\tau$ is now readily derived as

[TABLE]

so that

[TABLE]

The distribution function $F$ for $\tau$ is obtained by

[TABLE]

∎

By an entirely similar argument we obtain the following.

Theorem 2.3**.**

If $\tau\sim\mbox{IPH}(\bm{\pi},\bm{T}(t))$ then

[TABLE]

so that

[TABLE]

where $\bm{S}(u)=\bm{T}(s+u)$ and

[TABLE]

Hence the overshoot over a certain level is again phase–type distributed. Methods for numerical evaluation of transition probabilities given by the product integral is treated in e.g. Helton and Stuckwisch (1976) or ller (1992). We shall, however, concentrate on the cases where the matrices $\bm{T}(u)$ commute, which will simplify expressions and numerical methods considerably.

Corollary 2.4**.**

If $\tau\sim\mbox{IPH}(\bm{\pi},\bm{T}(t))$ and $\bm{T}(t_{1})$ and $\bm{T}(t_{2})$ commute for all $t_{1},t_{2}\geq 0$, then the density $f$ and the distribution function $F$ of $\tau$ are given by

[TABLE]

Proof.

If $\bm{T}(t_{1})$ and $\bm{T}(t_{2})$ commute for all $t_{1},t_{2}\geq 0$, then

[TABLE]

so

[TABLE]

from which the result follows. ∎

An important special case where all matrices $\bm{T}(\cdot)$ commute, are the class of sub–intensity matrices on the form $\bm{T}(t)=\lambda(t)\bm{T}$, which corresponds to the case where all transition rates of the Markov chain are scaled in the same way. More generally, block–diagonal matrices of the form

[TABLE]

commute whenever $\{\bm{S}(u)\}$ commute and $\{\bm{T}(u)\}$ commute. This includes the case where one of the two families of matrices are constant, e.g. $\bm{S}(u)=\bm{S}$ for all $u\geq 0$. In the latter case, the resulting distribution will be a mixture of an inhomogeneous and a homogeneous phase–type distribution.

Definition 2.5**.**

If $\bm{T}(t)=\lambda(t)\bm{T}$, where $\lambda(t)$ is some known non–negative real function and $\bm{T}$ is a sub–intensity matrix, then we shall write $\tau\sim\mbox{IPH}(\bm{\pi},\bm{T},\lambda)$ instead of $\mbox{IPH}(\bm{\pi},\bm{T}(t))$.

Corollary 2.6**.**

If $\tau\sim\mbox{IPH}(\bm{\pi},\bm{T},\lambda)$, then the density $f$ and the distribution function $F$ of $\tau$ are given by

[TABLE]

Proof.

Follows immediately from

[TABLE]

∎

Remark 2.7**.**

The commutativity condition of the matrices $\bm{T}(t)$ may be slightly relaxed by assuming only quasi–commutativity of the family $\{\bm{T}(t)\}$, which is characterized by the property

[TABLE]

see (Breuer, 2012, Section 1.4). In this case the above expressions for the density and the distribution function remain valid.

If $\tau\sim\mbox{IPH}(\bm{\pi},\bm{T},\lambda)$ and $\lambda(t)\equiv 1$, then $\tau$ has a conventional phase–type distribution (see Bladt and Nielsen (2017)), in which case we write $\tau\sim\mbox{PH}(\bm{\pi},\bm{T})$. Now consider $X\sim\mbox{PH}(\bm{\pi},\bm{T})$, $g:\mathds{R}_{+}\rightarrow\mathds{R}_{+}$ an increasing and differentiable function and $Y:=g(X)$. Then

[TABLE]

and consequently, using that $\bm{t}=-\bm{T}\bm{e}$,

[TABLE]

This suggests that a random variable $\tau\sim\mbox{IPH}(\bm{\pi},\bm{T},\lambda)$ can be obtained as a transformation of a phase–type distributed random variable $X\sim\mbox{PH}(\bm{\pi},\bm{T})$ as follows.

Theorem 2.8**.**

Let $\tau\sim\mbox{IPH}(\bm{\pi},\bm{T},\lambda)$ and let $g$ be defined in terms of its inverse function through

[TABLE]

Then

[TABLE]

where $X\sim\mbox{PH}(\bm{\pi},\bm{T})$.

Proof.

First we notice that $g$ is well-defined since $\lambda(t)>0$ for all $t$ and correspondingly $x\rightarrow\int_{0}^{x}\lambda(t)dt$ is strictly increasing. Then consider $Y=g(X)$. We have that

[TABLE]

Thus $Y\sim\tau$. ∎

If $h$ is an analytic function and $\bm{A}$ is a matrix, we define

[TABLE]

where $\gamma$ is a simple path enclosing the eigenvalues of $\bm{A}$. We refer to Section 3.4. of Bladt and Nielsen (2017) for further details.

Theorem 2.9**.**

Let $\tau\sim\mbox{IPH}(\bm{\pi},\bm{T},\lambda)$, with $g$ defined in terms of its inverse function through

[TABLE]

Assume that the Laplace transform for $g$, $L_{g}(s)$, exists for all $s>max_{i}\mbox{Re}(\lambda_{i})$, where $\lambda_{i}$ denote the eigenvalues for $\bm{T}$. Then the mean of $\tau$ is given by

[TABLE]

More generally, if $g^{\alpha}$ exists in the same domain as above, where $\alpha>0$, then

[TABLE]

where $L_{g^{\alpha}}$ denotes the Laplace transform of $g^{\alpha}$.

Proof.

Laplace transforms are analytic functions in the domain where they exist, so

[TABLE]

Similarly, if $L_{g^{\alpha}}$ exists, then it is analytic and

[TABLE]

∎

We notice the following trivial but important result.

Theorem 2.10**.**

Suppose that

[TABLE]

where $X_{1},...,X_{N}$ are independent and $X_{i}\sim\mbox{PH}(\bm{\pi},\bm{T})$. Let $X$ be distributed as the mixture of $X_{1},...,X_{N}$, with density

[TABLE]

where $\alpha_{1}+...+\alpha_{N}=1$. Then the mixture of $\tau_{1},...,\tau_{N}$ with probabilities $\alpha_{1},...,\alpha_{N}$ is distributed as $g(X)$.

3. Matrix–Pareto distributions

We start by proving the following result.

Theorem 3.1**.**

Let $Z=\exp(X)-1$ where $X\sim\mbox{PH}(\bm{\pi},\bm{T})$. Let $F_{Z},\bar{F}_{Z}$ and $f_{Z}$ denote the distribution function, survival function and density of $Z$, respectively. Then

[TABLE]

If the real part of the eigenvalue for $\bm{T}$ with maximum real part is smaller than $-\alpha$, then

[TABLE]

The Laplace transform of $Z$ is given by

[TABLE]

where

[TABLE]

is the upper incomplete $\Gamma$–function.

Proof.

Assuming that $z>0$, $x\rightarrow x\log(1+z)$ is analytic in the whole complex plane and so is then its exponential $x\rightarrow(1+z)^{x}$. Therefore $(1+z)^{\bm{T}}$ can be written in terms of the Cauchy–integral formula

[TABLE]

where $\gamma$ is a simple path in $\mathds{C}$ enclosing the eigenvalues for $\bm{T}$. The results then follow from Theorem 2.8 with $\lambda(y)={1}/(1+y)$ (which indeed results in $g(x)=e^{x}-1$) and Corollary 2.6.

The mean is straight-forward, as the mean of $e^{X}$ is simply the moment generating function of $X$ evaluated at 1. The $\alpha$–moment of $Y+1$ is obtained by noting that the Laplace transform of $e^{\alpha x}$ is $s\rightarrow 1/(s-\alpha)$ and using (5). Concerning the Laplace transform,

[TABLE]

since $\alpha\rightarrow\Gamma(\alpha,s)$ is an entire function for $s\neq 0$ (see Olver et al. (2010), 8.2) ∎

Remark 3.2**.**

In Ahn et al. (2012) similar results are proved for the distribution of $\exp(X)$ (not subtracting 1) apart from the Laplace transform. Their results are equivalent, though the presentation differs in notation and our use of functional calculus in the proof.

Remark 3.3**.**

If the real part of the eigenvalue for $\bm{T}$ with maximum real part is smaller than $-1$, then the mean of $Z=e^{X}-1$ exists and is given by

[TABLE]

Remark 3.4**.**

The matrix $\Gamma(\bm{T},s)$ is given by

[TABLE]

where $\gamma$ is a simple path enclosing the eigenvalues for $\bm{T}$. If all eigenvalues are simple (i.e. have multiplicity one), then $\Gamma(\bm{T},s)$ can be readily calculated by the residual theorem, as the terms only involve $\Gamma(\lambda_{i},s)$ divided by polynomials. On the other hand, if an eigenvalue $\lambda$ has multiplicity $m>1$, then also the $(m-1)$-order derivative of a polynomial times $\Gamma(z,s)$ w.r.t. $z$ will have to be evaluated at $\lambda$, which in turn involves derivatives of $z\rightarrow\Gamma(z,s)$ of orders $1,...,m-1$.

While the distribution of $\exp(X)$ can naturally be called a log–phase–type distribution, we shall refer to distributions of $\exp(X)-1$ as matrix–Pareto distributions, since the distribution of $Z$ from Theorem 3.1 may be seen as a Pareto distribution with a matrix parameter. A previous example conforming with this idea relates to the class of matrix–exponential distributions (containing the phase–type distributions), which bears its name from the fact that its elements can be seen as exponential distributions with a matrix parameter.

To gain generality, and in particular in order to draw on the analogy with the generalized Pareto distribution $G_{\xi,\beta}(x)=1-(1+\xi x/\beta)^{-1/\xi}$ for $\xi>0$ (cf. Embrechts et al. (1997)), we introduce the following additional scaling in the transformation (for reasons that become clear in Theorem 3.6, it is useful to express the scaling in units of $\mu=\mathds{E}(X)$):

Definition 3.5** (Matrix-Pareto distribution).**

Let $X\sim\mbox{PH}(\bm{\pi},\bm{T})$ and define

[TABLE]

where $\beta\in\mathbb{R}$. If we parametrize with

[TABLE]

we say that the distribution of $Y$ follows a matrix–Pareto distribution, and write

[TABLE]

If $Y\sim\mbox{M--Pareto}(\beta,\bm{\pi},\bm{T})$, then

[TABLE]

where $Z=e^{X}-1$. Then from Theorem 3.1 we get that

[TABLE]

and

[TABLE]

Moreover, by a simple transformation of the formula in Theorem 3.1, we get that $Y$ has Laplace transform

[TABLE]

Next consider $\mathds{P}(Y>x+y|Y>x)$. From Theorem 2.3 we then get that

[TABLE]

from which we see that conditionally on $Y>x$, $Y$ has a matrix–Pareto distribution given by

[TABLE]

We collect the previous results in the following theorem.

Theorem 3.6**.**

Let $Y\sim\mbox{M--Pareto}(\beta,\bm{\pi},\bm{T})$. Denote by $\bar{F}(y)=\bar{F}_{\beta,\bm{\pi},\bm{T}}(y)=\mathds{P}(Y>y)$ the survival function of $Y$. Then

(a) $\bar{F}_{\beta,\bm{\pi},\bm{T}}(y)=\bm{\pi}\left(1+\mu\frac{y}{\beta}\right)^{\bm{T}}\bm{e}$ .

(b) $f_{Y}(y)=\frac{\mu}{\beta}\bm{\pi}(1+\mu\frac{y}{\beta})^{\bm{T}-\bm{I}}\bm{t}$.

(c) $\mathds{P}(Y>x+y|Y>x)=\bar{F}_{\beta+\mu x,\tilde{\pi},\bm{T}}(y)$, where

[TABLE]

(d) $L_{Y}(s)=e^{\beta s/\mu}\bm{\pi}\left(\frac{\mu}{s\beta}\right)^{\bm{T}}\Gamma(\bm{T},s\beta/\mu)\bm{t}$.

(e) If the real part of the eigenvalue for $\bm{T}$ with maximum real part is less than $-\alpha$, then

[TABLE]

Remark 3.7**.**

Note that these results can indeed be seen as generalizations of results from Theorem 3.4.13 of Embrechts et al. (1997), which stated such properties for the generalized Pareto distribution (which for $\xi>0$ refers to the scalar case $p=1$ in our setting). In particular, property (c) above (cf. also (Ahn et al., 2012, Eqn.(7))) is a considerable extension of (Embrechts et al., 1997, Equ.3.53), which can be quite useful when using models based on matrix-Pareto distributions. For instance, one reason for the popularity of Pareto distributions among heavy-tailed distributions in reinsurance modelling are the simple pricing formulas for XL reinsurance contracts based on (Embrechts et al., 1997, Equ.3.53), where $x$ will typically refer to the retention level. The present extension shows that these properties can be carried over to the much more general class of matrix-Pareto distributions. ∎

Example 3.8**.**

Consider the distribution of

[TABLE]

where $X$ has an Erlang distribution $\mbox{Er}_{n}(\lambda)$, i.e. $X\sim X_{1}+...+X_{n}$ for i.i.d. $X_{i}\sim$ Exp$(\lambda)$. Then $Y\sim\mbox{M-Pareto}(\beta,\bm{\pi},\bm{T})$ with $\beta=\mu=\mathds{E}(Y)$ (which is a particular type of a log-Gamma distribution). A phase–type representation of the $\mbox{Er}_{n}(\lambda)$ distribution is $\bm{\pi}=(1,0,...,0)$ and

[TABLE]

where the vector is $n$–dimensional and $\bm{T}$ is an $n\times n$ matrix. Now

[TABLE]

so

[TABLE]

is a matrix with $(i,j)$th element, $j\geq i$, equal to

[TABLE]

Here $\gamma$ is a simple path about the multiple eigenvalue $-\lambda$. For $j<i$ the elements are clearly zero. Defining

[TABLE]

we get by the residue theorem

[TABLE]

where $\phi^{(k)}(s)$ denotes the $k$th order derivative of $\phi$. But

[TABLE]

so we conclude that

[TABLE]

Then the tail of $Y$ is seen to be the sum of the terms in the first row of $(1+y)^{\bm{T}}$,

[TABLE]

The density is given by

[TABLE]

since $\bm{t}=(0,0,...,0,\lambda)$, so the density is essentially given by $\lambda$ times the $(1,n)$th element of $(1+x)^{\bm{T}}$. ∎

An immediate consequence of Example 3.8 is that if $X$ is a mixture of Erlang distributions $\mbox{Er}_{n_{i}}(\lambda_{i})$, $i=1,...,N$, i.e.

[TABLE]

where $\alpha_{1}+...+\alpha_{N}=1$, then $Y$ is a mixture of distributions on the form (6), i.e.

[TABLE]

cf. also Theorem 2.10. Mixtures of Erlang distributions constitute the smallest and simplest sub–class of phase–type distributions which is dense in the class of distributions on the positive real line, and has often been employed in applications (see e.g. Willmot and Woo (2007)). Since $f:x\rightarrow\exp(x)-1$ and $f^{-1}:x\rightarrow\log(x+1)$ are continuous functions, by the continuous mapping theorem it is clear that also the class of mixtures of log–Erlang distributions is dense in the class of distributions on the positive half-line (this denseness was also already noted in (Ahn et al., 2012, Th.1) for general log-PH distributions). Due to its potential for modelling purposes, we formulate the above result as a theorem.

Theorem 3.9**.**

The class of matrix–Pareto distributions generated by mixtures of Erlang distributions

[TABLE]

contains densities of the form

[TABLE]

The class is dense (in the sense of weak convergence) in the class of distributions on the positive real line.

This means that any distribution may be approximated arbitrarily closely by either a mixture of Erlang distributions or a mixture of matrix–Pareto distributions generated by Erlangs. Whether one chooses one over the other class when approximating a distribution in practice should depend on the tail of the distribution to be approximated. If the tail is of Pareto or log–Gamma type, then one will obtain a much better approximation using the matrix–Pareto, as for the mixtures of Erlangs the number of phases needed would be large and still the resulting approximation will not capture the tail behaviour well. We illustrate this important point in Section 5 by fitting a general matrix–Pareto distribution to a heavy-tailed dataset, and it will turn out that an excellent fit is in fact obtained by a matrix–Pareto which is close to being generated by a mixture based on Erlangs.

Example 3.10**.**

Consider the distribution of

[TABLE]

where $X$ has a generalized Erlang distribution $\mbox{Er}_{n}(\lambda_{1},...,\lambda_{n})$, i.e. $X\sim X_{1}+...+X_{n}$ for i.i.d. $X_{i}\sim\mbox{Exp}(\lambda_{i})$, $i=1,...,n$. Assume for simplicity that all $\lambda_{i}$ are different. Then a similar calculation as above reveals that the density for $Y$ is given by

[TABLE]

∎

Example 3.11**.**

While the derivation of Theorem 3.6(c) basically relies on a probabilistic argument using transition probabilities, most distributional properties in our construction rely only on the matrix–exponential forms of the phase–type distributions. To emphasize this point, we here consider a distribution which can be written in a matrix–exponential form but which is not phase–type. For such distributions the matrix–exponentials are no longer transition matrices. Distributions whose density can be written in terms of a matrix–exponential, including the phase–type distributions, are known as matrix–exponential distributions (see Bladt and Nielsen (2017) for further details) and consist exactly of the distributions with rational Laplace transform. Concretely,

[TABLE]

is an example of a density that does not belong to a random variable of phase–type but has rational Laplace transform

[TABLE]

and can be written in matrix–exponential form as

[TABLE]

Let $\bm{\pi}=(101,0,0)$, $\bm{t}=(0,0,1)^{\prime}$ and let $\bm{T}$ denote the matrix inside the exponential. Define the resolvent by

[TABLE]

and let

[TABLE]

Then the density for $Y$ is given by

[TABLE]

and since only the first and last elements of $\bm{\pi}$ and $\bm{t}$ are different from zero, respectively, only the $(1,3)-$element in $\bm{R}$ is of importance. Since

[TABLE]

the density for $Y$ is readily calculated (using the residue theorem) to be

[TABLE]

which simplifies to

[TABLE]

In Figure 1, both $f$ and $f_{Y}$ are plotted. As might be expected, one observes that the shape of the distribution is somewhat preserved while the tail is stretched out.

∎

4. Analogous generalisations of other classes of distributions

We now list a number of distributions obtained through a transformation of a PH random variable other than the exponential function. All the these transformations, as opposed to the M–Pareto distribution, are parameter–dependent, and as such fitting procedures to data are less straightforward.

4.1. Matrix-Weibull distributions

Let $X\sim\mbox{PH}_{p}(\bm{\pi},\bm{T})$. Inspired by the construction of a Weibull random variable as a power of an exponentially distributed random variable, one can define

[TABLE]

for some $\beta>0$. The distribution function of $Y$ is immediately seen to be

[TABLE]

with corresponding density

[TABLE]

$Y$ can be called a matrix-Weibull distribution, since the scale parameter of the usual Weibull distribution with cumulative distribution function $F_{Wei}(y)=1-\exp(-cy^{\beta})$ is now replaced by a matrix. Alternatively, one may refer to $Y$ as a power-PH distribution. Its mean is given by

[TABLE]

where the last line follows from Theorem 3.4.4 of Bladt and Nielsen (2017). More generally, the $\theta$th moment ($\theta>0$) can be deduced similarly to be

[TABLE]

The moment generating function is found by expansion,

[TABLE]

which is valid for $\beta>1$.

If $(\bm{\pi},\bm{T})$ is the representation of an $\mbox{Er}_{n}(\lambda)$–distribution (see Example 3.8), then

[TABLE]

Therefore, as for the matrix–Pareto distribution, a dense class of distributions with a Weibull tail can then be defined through mixtures of the form

[TABLE]

It is evident that this class will be a very appropriate choice for approximating distributions which have tails (suspected to be) close to a Weibull.

4.2. Exponential–PH distribution

Definition 4.1**.**

Let $X\sim\mbox{PH}_{p}(\bm{\pi},\bm{T})$. Then define an exponential-PH distribution with parameters $\mu$, $\sigma$, $\bm{\pi}$ and $\bm{T}$ as the distribution of the random variable

[TABLE]

The distribution function of $Y$ is given by

[TABLE]

with support on $(-\infty,\infty)$. The density is obtained through differentiation as

[TABLE]

One can interpret (7) as a matrix extension of the usual Gumbel distribution with cumulutative distribution function $F_{Gu}(y)=\exp(-\exp(-(y-\mu)/\sigma))$, but since here the matrix $\bm{T}$ replaces the constant $-1$ rather than a parameter, it is less obvious to call it a matrix-Gumbel distribution, and one may prefer the name exponential-PH distribution. Concerning the mean we proceed as follows. Since

[TABLE]

(use (4) or see e.g. Bladt and Nielsen (2017) p. 180) we get that

[TABLE]

where $\gamma$ is Euler’s constant. In particular, if $X\sim\mbox{Exp}(1)$ then $\bm{T}=-1$, the term $\log(-\bm{T})=0$ and $\mathds{E}(Y)=\mu+\gamma\sigma$.

Concerning the Laplace transform for $Y$,

[TABLE]

which follows from

[TABLE]

and where the last step follows from Theorem 3.4.6 of Bladt and Nielsen (2017) (p. 175). For the special (scalar) case of $X\sim\mbox{Exp}(1)$, we recuperate the known formula $e^{-\mu s}\Gamma(1+\sigma s)$ of the Laplace transform of the Gumbel distribution.

If $(\bm{\pi},\bm{T})$ is the representation of an $\mbox{Er}_{n}(\lambda)$–distribution (see Example 3.8), then

[TABLE]

so a dense class of distributions with a Gumbel type of tail can be defined in terms of distributions of the form

[TABLE]

4.3. Shifted-power–PH distribution

More generally, inspired by the Generalized Extreme Value distribution $\mbox{GEV}(\mu,\sigma,\xi)$ with cumulative distribution function

[TABLE]

we see that for $X\sim\mbox{Exp}(1)$ one has

[TABLE]

That is, for $\xi\neq 0$ one can construct $Y\sim\mbox{GEV}(\mu,\sigma,\xi)$ from an exponentially distributed random variable $X\sim\mbox{Exp}(1)$ by letting

[TABLE]

Taking the limit $\xi\rightarrow 0$ then leads back to

[TABLE]

which is the Gumbel case treated in Section 4.2, so we focus on $\xi\neq 0$ in the sequel.

Letting now $X\sim\mbox{PH}(\bm{\pi},\bm{T})$ replace the exponentially distributed random variable above, we apply the transformation (9) for $\xi\neq 0$. The cumulative distribution function for $Y$ is then readily obtained as

[TABLE]

and the density is obtained through differentiation

[TABLE]

where

[TABLE]

Again, the constant -1 in (8) is replaced by the matrix $\bm{T}$ in (10) now, so that one could view the latter as a certain matrix-GEV distribution, but the name shifted-power–PH distribution may be considered more appropriate. The mean of $Y$ is

[TABLE]

which follows from Theorem 3.4.6 of Bladt and Nielsen (2017) whenever $\xi<1$. Higher order integer moments can be calculated recursively by expanding the power of $Y$ in terms of $X$ and using Theorem 3.4.6 of Bladt and Nielsen (2017) again.

If $(\bm{\pi},\bm{T})$ is the representation of an $\mbox{Er}_{n}(\lambda)$–distribution, then

[TABLE]

so a dense class of distributions with GEV–type tails can be defined in terms of mixtures

[TABLE]

5. Modelling with a matrix–Pareto distribution

As mentioned earlier, the denseness of PH distributions in the class of all distributions on the positive half-line makes them an attractive modelling tool. However, traditionally one disadvantage of this class is that by its exponentially decaying tails it may require many states to be a reasonable fit to distributions with a heavy tail (even on finite time intervals). From the denseness of matrix–Pareto distributions in the class of all distributions on the positive half-line (Theorem 3.9 or (Ahn et al., 2012, Th.1)), it becomes clear that if one reckons a heavy-tailed distribution to describe a given set of data-points well, then a fitting procedure with matrix–Pareto distributions may lead to a much more parsimonious (and natural) model (quantitatively expressed by a low number of states and a sparse intensity matrix $\bm{T}$) of the underlying distribution yet keeping the advantageous properties of phase-type distributions. In Ahn et al. (2012), a Danish fire insurance dataset was investigated in this direction. Here, in the context of the present paper, we would like to illustrate the potential of the inhomogeneous PH approach proposed in this paper for a set of 1282 Dutch fire insurance claims from the period 2000-2014 (all in excess of 1 mio.  EUR), which were already studied in detail and modelled by other means in Albrecher et al. (2017). Indeed, after testing various models, in (Albrecher et al., 2017, p.107) a splicing model with a mixed Erlang component for the bulk, a Pareto distribution for the range to the right of that, and another Pareto distribution for the tail was identified to be an appropriate choice for that dataset. Apart from seeking an adequate fit by a matrix–Pareto distribution, we shall pay special attention to the aforementioned (possibly mixed) Erlang structure and compare it to the splicing method.

In order to fit the totality of the data by one single distribution, we propose alternatively a matrix–Pareto distribution (log–Phase–type), which also contains the Pareto tail as a special case. Thus we must fit a phase–type distribution to the logarithm of the data. Since the log–data have a support contained in $[8.5,20]$, we decide to subtract 8.5 from all log data which makes it possible to reduce the number of phases significantly. Hence the model we consider is

[TABLE]

where $y_{i}=\log(x_{i})-8.5\ (i=1,\ldots,N=1282)$ are i.i.d. transformed data points $x_{1},...,x_{N}$. The implied model for the original generic data $X$ is then

[TABLE]

i.e. a (scaled) log–phase–type distribution.

Utilizing an EM-algorithm (cf. Asmussen et al. (1996)), we use 20 phases to obtain an adequate fit to the data and the result is presented in Figure 2.

The choice of 20 phases is taken so as to allow a sufficiently flexible phase–type fit to the data, the concrete magnitude being supported by trial experiments. Though the starting point of the EM algorithm then was a general phase–type distribution, the resulting estimate reduces to a much simpler structure, a phenomenon which is commonly observed in phase–type fitting. Concretely, the structure of the fitted intensity matrix is as follows

[TABLE]

where $a=3.322438$, $a_{1}=3.292191$, $a_{2}=3.322884$, $a_{3}=3.292191$, $b=3.288628$, $b_{1}=2.063919$ and $c=48.21421$ . The estimator for the initial distribution is

[TABLE]

Hence we may describe this distribution in terms of 8 parameters. The distribution consists of a very fast state (the last one), which can be entered directly with the very small probability 0.00078, representing a possibility to have very small outcomes (one can interpret this to stem from some data points in that region due to the left-truncation of the data at $y_{i}=0$). Furthermore, there are two blocks of 9–dimensional convolutions of exponential distributions (i.e., Er9 distributions) with intensities $a$ and $b$, respectively. Exits to the absorbing state are possible from states 9,10,19 and 20.

The fit of the phase–type data (cf. Figure 2) looks adequate in the main body and most of the tail (QQ–plot), which is quite remarkable when compared to the much more complex model suggested in (Albrecher et al., 2017, p.107). The seven most extreme tail points start to deviate from the fit, being most pronounced for the last three points. This could be expected since maximum likelihood treats all points of the distribution equally as opposed to extreme value methods focussing on the tail fit. In any case, the identified matrix-Pareto distribution seems to be an excellent and parsimonious fit to the data over quite a large range.

From the concrete numerical values, one can also see that a simple matrix-Pareto distribution with underlying Er19 distribution (with rate around 3.3, and hence altogether only two parameters(!)) might be a (surprisingly) reasonable description of the data as well. To further pursue this point, we fitted an $\mbox{Er}_{19}(\lambda)$ (i.e. a convolution of 19 i.i.d. exponential distributions with intensities $\lambda$) to the data resulting in an estimated intensity of $\hat{\lambda}=3.340752$. The Erlang distribution indeed also fits the data well in a first approximation (cf. the dashed blue line in Figure 2), though it is clear that the full PH is a more suitable model for small values. The log–likelihoods for the full phase–type and Erlang fits are -1250.092 and -1366.735, respectively.

In conclusion, one may argue that we have obtained a parsimonious, and still somewhat similar model as the one in Albrecher et al. (2017) but using a quite different rationale, and without the need to decide about the location of splicing points. In that respect, it has the advantage to be rather easily implemented in practice.

6. Conclusion

In this paper we proposed a simple way to carry over the advantages of PH distributions as a modelling tool to distributions with a different tail behaviour. This can be achieved by introducing time-inhomogeneity in the Markov jump process underlying the construction of the PH distribution. In case the respective time scaling is the same for all states of the Markov process, this leads to a number of new simple families of distributions (including matrix-Pareto distributions, (shifted) power-PH distributions, exponential-PH distributions), that all inherit the denseness in the class of distributions on the positive real line, but can lead to much more parsimonious descriptions than a direct PH fit, when the tail of the distribution is in fact a corresponding transform of a PH tail. In view of its relevance in applications, we focused on heavy tails and in particular established various properties of matrix-Pareto distributions. We finally studied the modelling performance of the latter on a real-world dataset from fire insurance, with rather promising results.

There are several possible directions for future research. One may be to more directly merge light- and heavy-tail fitting by time-scaling only some of the states of the Markov process, and possibly design an EM algorithm that automatically decides on the basis of the dataset how many components of each type are needed for a good description of the data. It should also be feasible to extend the approach proposed in this paper to complement splicing models under censoring, see e.g.  Reynkens et al. (2017). Finally, adapting extreme value techniques to the dense distribution classes proposed in this paper may provide interesting alternatives for parsimonious modelling with emphasis on the appropriateness of the fit in the tail, but still good performance for smaller values.

Acknowledgement. H.A. acknowledges financial support from the Swiss National Science Foundation Project 200021_168993.