Matrix Mittag--Leffler distributions and modeling heavy-tailed risks

TL;DR

This paper introduces matrix Mittag-Leffler distributions, a flexible family for modeling heavy-tailed risks, which can be interpreted through phase-type distributions or semi-Markov processes, and demonstrates their effectiveness on simulated and real insurance data.

Contribution

It defines matrix Mittag-Leffler distributions and shows their advantages over existing methods in modeling heavy-tailed risks, addressing threshold selection issues.

Findings

The class can be interpreted as inhomogeneous phase-type distributions with random scaling.

It effectively models heavy-tailed risks in insurance data.

Demonstrates advantages over traditional extreme value approaches.

Abstract

In this paper we define the class of matrix Mittag-Leffler distributions and study some of its properties. We show that it can be interpreted as a particular case of an inhomogeneous phase-type distribution with random scaling factor, and alternatively also as the absorption time of a semi-Markov process with Mittag-Leffler distributed interarrival times. We then identify this class and its power transforms as a remarkably parsimonious and versatile family for the modelling of heavy-tailed risks, which overcomes some disadvantages of other approaches like the problem of threshold selection in extreme value theory. We illustrate this point both on simulated data as well as on a set of real-life MTPL insurance data that were modeled differently in the past.

Figures

Equations

$$E_{\alpha ,\beta }(z)=\sum _{k=0}^{\infty }\frac{z^k}{\Gamma (\alpha k+\beta )},z\in \mathbb{C},$$

$$\frac{{\mathrm{d}}^m}{\mathrm{d}{z}^m}\left[z^{\beta -1}{E}_{a,\beta }\left(z^a\right)\right]=z^{\beta -m-1}{E}_{\alpha ,\beta -m}\left(z^{\alpha }\right).$$

$$E_{\alpha ,\beta }^{(k)}(z)=\frac{{\mathrm{d}}^k}{\mathrm{d}{z}^k}{E}_{\alpha ,\beta }(z)=\frac{1}{{\alpha }^k{z}^k}\sum _{j=0}^k{c}_j^{(k)}{E}_{\alpha ,\beta -j}(z)$$

$$c_{j}^{(k)}=\left\{\begin{array}[]{ll}{(1-\beta-\alpha(k-1))c_{0}^{(k-1)}},&{j=0},\\ c_{j-1}^{(k-1)}+(1-\beta-\alpha(k-1)+j)c_{j}^{(k-1)},&j=1,\ldots,k-1,\\ {1},&{j=k}.\end{array}\right.$$

$$E_{\alpha ,\beta }(\mathbf{A})=\sum _{k=0}^{\infty }\frac{{\mathbf{A}}^k}{\Gamma (\alpha k+\beta )}.$$

$$E_{\alpha ,\beta }(\mathbf{A})=\frac{1}{2\pi \text{mbox}i}{\int }_{\gamma }{E}_{\alpha ,\beta }(z)(z\mathbf{I}-\mathbf{A}{)}^{-1}\mathrm{d}z,$$

$${\mathbf{J}}_i=\left(\begin{array}{ccccc}{\lambda }_i& 1& 0& \cdots & 0\\ 0& {\lambda }_i& 1& \cdots & 0\\ 0& 0& {\lambda }_i& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮⋮⋮& ⋮\\ 0& 0& 0& \cdots & {\lambda }_i\end{array}\right),$$

$$E_{\alpha ,\beta }(\mathbf{A})=\mathbf{P}\text{mbox}diag(E_{\alpha ,\beta }({\mathbf{J}}_1),...,E_{\alpha ,\beta }({\mathbf{J}}_r)){\mathbf{P}}^{-1},$$

$$E_{\alpha ,\beta }({\mathbf{J}}_i)=\left(\begin{array}{ccccc}{E}_{\alpha ,\beta }({\lambda }_i)& E_{\alpha ,\beta }^{(1)}({\lambda }_i)& \frac{E_{\alpha ,\beta }^{(2)}({\lambda }_i)}{2!}& \cdots & \frac{E_{\alpha ,\beta }^{(m_i-1)}({\lambda }_i)}{(m_i-1)!}\\ 0& E_{\alpha ,\beta }({\lambda }_i)& E_{\alpha ,\beta }^{(1)}({\lambda }_i)& \cdots & \frac{E_{\alpha ,\beta }^{(m_i-2)}({\lambda }_i)}{(m_i-2)!}\\ 0& 0& E_{\alpha ,\beta }({\lambda }_i)& \cdots & \frac{E_{\alpha ,\beta }^{(m_i-3)}({\lambda }_i)}{(m_i-3)!}\\ ⋮& ⋮& \cdots & ⋮⋮⋮& ⋮\\ 0& 0& 0& \cdots & E_{\alpha ,\beta }({\lambda }_i)\end{array}\right),$$

$$F_{\delta ,\alpha }(x)$$

$$f_{\delta ,\alpha }(x)$$

$$\frac{1}{1+(\delta u{)}^{\alpha }}.$$

$$X\stackrel{d}{=}\delta Z{R}^{1/\alpha },$$

$$F_R(x)=\frac{2}{\pi \alpha }\left[\arctan \left(\frac{x}{\sin (\alpha \pi /2)}+\cot (\alpha \pi /2)\right)-\frac{\pi }{2}\right]+1.$$

$$\mathbb{E}(e^{-uX})=(1+u^{\alpha }{)}^{-\beta }.$$

$$F_{\alpha ,\beta }(x)$$

$$X\stackrel{d}{=}{W}^{1/\alpha }{S}_{\alpha },$$

$$\mathbb{E}(e^{-u{S}_{\alpha }})=\exp (-u^{\alpha }).$$

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

$$\varphi (u)=\mathbf{\pi }(u^{\alpha }\mathbf{I}-\mathbf{T}{)}^{-1}\mathbf{t},0<\alpha \le 1,$$

$$f(x)=\mathbf{\pi }(x\mathbf{I}-\mathbf{T}{)}^{-1}\mathbf{t}.$$

$$f^{(n)}(g(u))=(-1{)}^{n}n!\mathbf{\pi }{\left(u^{\alpha }\mathbf{I}-\mathbf{T}\right)}^{-n-1}\mathbf{t},$$

$$g^{(j)}(u)=\alpha (\alpha -1)\cdots (\alpha -j+1)u^{\alpha -j}$$

$$\frac{{\mathrm{d}}^n}{\mathrm{d}{x}^n}f(g(x))=\sum \frac{n!}{m_1!m_2!\cdots m_n!}\cdot f^{\left(m_1+\cdots +m_n\right)}(g(x))\cdot \prod _{j=1}^n{\left(\frac{g^{(j)}(x)}{j!}\right)}^{m_j},$$

$$1\cdot m_1+2\cdot m_2+3\cdot m_3+\cdots +n\cdot m_n=n,$$

$$\text{mbox}sign\left(f^{(m_1+m_2+\cdots +m_n)}(g(u))\prod _{i=1}^n{g}^{(i)}(u{)}^{m_i}\right)$$

$$\varphi (u)=\mathbf{\pi }(g(u)\mathbf{I}-\mathbf{T}{)}^{-1}\mathbf{t}$$

$$f(x)=x^{\alpha -1}\mathbf{\pi }{E}_{\alpha ,\alpha }\left(\mathbf{T}{x}^{\alpha }\right)\mathbf{t}.$$

$${\int }_0^{\infty }{e}^{-ux}f(x)\mathrm{d}x$$

$$\mathbb{E}(e^{-uX})=\mathbf{\pi }(u^{\alpha }\mathbf{I}-\mathbf{T}{)}^{-1}\mathbf{t}.$$

$$f(x)=\mathbf{\pi }{e}^{\mathbf{T}x}\mathbf{t}$$

Full text

1. Introduction

The modeling of heavy-tailed risks is a classical topic in probability, statistics and its applications to understand and interpret data, see e.g. Embrechts et al. [13], Beirlant et al. [5] and Klugman et al. [20]. The folklore heavy-tailed distributions like Pareto, (heavy-tailed) Weibull and lognormal distributions can often serve as very useful benchmark models, particularly when only a few data points are available. In addition, the Pareto distribution is simple to work with and has an intuitive justification in terms of limit properties of extremes. However, in situations with more (but not an abundance of) data points, one often empirically observes that a simple Pareto distribution does not serve as a good model across the entire range of the distribution. This is also the case for more general parametric families like Burr or Benktander distributions. Traditionally, and according to a main paradigm of the extreme value statistics approach, this is handled by only using the largest available data points to estimate the tail behavior, and model the bulk of the distribution separately by another distribution, finally splicing together the respective parts (see e.g. Albrecher et al. [1, Ch.4] for details). In insurance practice one often refers to this separate modeling as the modeling of attritional and large claims, and the resulting models are frequently referred to as composite models [26]. A natural problem in this context is how to choose the threshold between the separate regions, often boiling down to the compromise of not leaving too few data points for the tail modeling. At the same time, the consequences of that choice can be considerable, for instance for the determination of solvency capital requirements in insurance (cf. [1, Ch.6] for illustrations). A considerable effort has therefore been made to develop techniques and criteria for an appropriate choice of such thresholds, see e.g. [5] for an overview and [7] for a recent contribution in that direction.

If the confidence in the relevance of available data points for the description of the (future) risk is sufficiently high, another possible approach is to use a tractable, but much larger family of distributions and identify a good fit. A particularly popular candidate for such an approach is the class of phase-type (PH) distributions, see e.g. Asmussen et al. [4]. The class of PH distributions is dense (in the sense of weak convergence) in the class of distributions on the positive real line, meaning that they can approximate any positive distribution arbitrarily well. They are, however, light-tailed, which may be a problem in applications which require a heavier tail and where the quantities of interest heavily depend on the tail behavior (as e.g. for ruin probabilities, cf. [3]). Fitting heavy-tailed distributions with a PH distribution can then lead to requiring many phases (rendering its use computationally cumbersome), and the resulting model will still not capture the tail behavior in a satisfactory manner. Two approaches to remedy this problem are Bladt & Rojas-Nandayapa [10] and Bladt et al. [9]. In Albrecher & Bladt [2] recently another direction was suggested, namely to transform time in the construction of PH distributions (as absorption times of Markov jump processes), leading to inhomogeneous phase-type (IPH) distributions. For suitable transformations, this approach allows to transport the versatility of PH distributions into the domain of heavy-tailed distributions, by introducing dense classes of genuinely heavy-tailed distributions. As a by-product, it was shown in [2] that a class of matrix-Pareto distributions can be identified, where the scalar parameter of a classical Pareto distribution is replaced by a matrix, providing an intuitive and somewhat natural extension of the Pareto distribution, much as the matrix-exponential distribution, which is a powerful extension of the classical exponential distribution, see e.g. Bladt & Nielsen [8].

In this paper we establish another matrix version of a distribution, namely the Mittag-Leffler distribution (first studied by Pillai [27]). While the identification of matrix versions of distributions is of mathematical interest in its own right, we will show that the resulting matrix Mittag-Leffler distributions (and its power transforms) have favorable properties for the modeling of heavy tails, and it can outperform some other modeling approaches in a remarkable way. Furthermore, we will identify this class of distributions as a particular extension of the IPH class, where the scaling is random. In addition, we will establish the matrix Mittag-Leffler distribution as the absorption time of a semi-Markov process with (scalar) Mittag-Leffler distributed inter-arrival times, extending the role of the exponential distribution for the inter-arrival in continuous-time Markov chains.

The Mittag-Leffler function was first introduced in [24] and over the years turned out to be a crucial object in fractional calculus. It can be seen as playing the same role for fractional differential equations as the exponential function does for ordinary differential equations, see e.g. Gorenflo et al. [17] for a recent overview. A recent application of fractional calculus for a particular risk model in insurance can be found in Constantinescu et al. [12]. Mittag-Leffler functions with matrix argument were first introduced in Chikrii and Eidel’man [11] and play a prominent role for identifying solutions of systems of fractional differential equations, see e.g. Garrappa and Popolizio [16]. Here we will use them to define the class of matrix Mittag-Leffler (MML) distributions, which enjoy some attractive mathematical properties and are heavy-tailed with a regularly varying tail with index $\alpha<1$. While such extremely heavy tails with resulting infinite mean can be relevant in the modeling of operational risk [25] and possibly insurance losses due to natural catastrophes [1], in most applications of interest the tails are slightly less heavy. We therefore enlarge the class of matrix Mittag-Leffler distributions by also including its power transforms and estimate the corresponding power together with the other parameters from the data in the fitting procedure. The index of regular variation of this larger class of distributions can now be any positive number. Whereas the number of needed phases for a PH or IPH fit can be very large also due to multi–modality or other irregularities in the shape of the main body of the distribution, we will see that the class of matrix Mittag-Leffler distribution and its power transforms (which we call power matrix Mittag-Leffler (PMML) distributions) offers a significant reduction in the number of phases needed to obtaining adequate fits. For this reason, it can even be worthwhile to scale light-tailed data points to heavy-tailed ones first, then apply a matrix Mittag-Leffler fit to the latter and transform the fit back to the original light-tailed scaling. This procedure is to some extent the reverse direction of the philosophy that underlied the PH fitting of heavy-tailed distributions. We will illustrate the potential advantage of this alternative approach in the numerical section at the end of the paper.

The remainder of the paper is organized as follows. Section 2 recollects some useful definitions and properties of Mittag-Leffler functions, Mittag-Leffler distributions and PH distributions. Section 3 then defines matrix Mittag-Leffler distributions and derives a number of its properties. We also give three explicit examples. Section 4 establishes MML distributions as IPH distributions under a particular random scaling, which allows to intuitively understand the additional flexibility gained from using MML distributions for the modeling of heavy-tailed risks. We then also establish MML distributions as the absorption times of a semi-Markov process with ML distributed interarrival times, which is yet another perspective on the potential of the MML class as a modeling tool. Finally, Section 5 is devoted to the modeling of data using (power-transformed) MML distributions. We first illustrate the convincing performance of the numerical fitting procedure to something as involved as tri–modal data. Secondly, we consider an MTPL data set taken from [1] and already studied by various other means in the literature. We show that a plain maximum-likelihood fit to this data set gives a convincing fit to the entire range of the data with remarkably few parameters, and even identifies the tail index with striking accuracy when compared to recent extreme value techniques as in [7], without having to choose a threshold for the tail modeling at all. We then also provide an example where transforming light-tailed data into heavy-tailed ones, fitting with a PMML distribution and transforming back can lead to a much better fit for the same number of parameters than a classical phase-type distribution. We then also discuss the signature of MML distributions in the tail in terms of the behavior of the Hill plot, which allows to develop an intuitive guess as to when MML distributions are particularly adequate for a fitting procedure of heavy tails. Finally, Section 6 concludes.

2. Some relevant background

2.1. Mittag–Leffler functions

The Mittag–Leffler (ML) function is defined by

[TABLE]

where $\beta\in\mathbb{R}$ and $\alpha>0$. The ML function is an entire function if $\beta>0$, and it satisfies (see e.g. Erdelyi et al.  [14, p.210])

[TABLE]

This implies that (see [16, Prop.2])

[TABLE]

where

[TABLE]

For a matrix $\bm{A}$, we may define its ML function as

[TABLE]

If $\beta>0$, one can then express the entire ML function of a matrix $\bm{A}$ by Cauchy’s formula

[TABLE]

where $\gamma$ is a simple path enclosing the eigenvalues of $\bm{A}$. If $\bm{A}$ has a Jordan normal form $\bm{A}=\bm{P}\,\mbox{diag}(\bm{J}_{1},...,\bm{J}_{r})\,\bm{P}^{-1}$ with

[TABLE]

then we may equivalently express $E_{\alpha,\beta}(\bm{A})$ by

[TABLE]

where

[TABLE]

and $m_{i}$ is the dimension of $\bm{J}_{i}$.

In either case, we shall need to evaluate the derivatives of ML functions at the eigenvalues, which by (1) implies the evaluation of ML functions with possibly negative indices $\beta-j$. This is not a problem, but it is important that initially $\beta>0$ to ensure that the ML function is entire and thereby the existence of the Cauchy integral formula is guaranteed.

For further properties on the ML function we refer e.g. to [14], [16], [22] and [18].

2.2. Mittag–Leffler distributions

A random variable having Mittag-Leffler (ML) distribution was defined in Pillai [27] through the cumulative distribution function and consequently density given by

[TABLE]

with Laplace transform

[TABLE]

A convenient representation, due to Kozubowski [21], for a ML random variable $X$ is

[TABLE]

where $Z$ is standard exponential and $R$ has cumulative distribution function

[TABLE]

The tail behaviour of $R$ is equivalent to that of a Cauchy random variable, and hence $X$ is regularly varying with parameter $\alpha$ in the tail (see e.g. [23, Prop.1.3.9]).

The following extension (for the case $\delta=1$) will also play a role in the sequel: a random variable $X$ is said to follow a generalized Mittag-Leffler (GML) distribution with parameters $\alpha$ $(0<\alpha\leq 1)$ and $\beta>0$, if its Laplace transform is given by

[TABLE]

The corresponding cumulative distribution function then is

[TABLE]

where $\mathfrak{B}(x,y)$ is the Beta function (see Jose et al. [19]). The analogous representation for a GML variable $X$ is

[TABLE]

where $W$ is Gamma with scale parameter 1 and shape parameter $\beta$, and $S_{\alpha}$ is a random variable with Laplace transform given by

[TABLE]

2.3. Phase–type distributions

A random variable $\tau$ is said to be phase–type distributed with generator (or representation) $(\bm{\pi},\bm{T})$, and we write $\tau\sim\mbox{PH}(\bm{\pi},\bm{T})$, if it is the time until absorption of a (time–homogeneous) Markov jump process $\{X_{t}\}_{t\geq 0}$ with state–space $E=\{1,2,...,p,p+1\}$ where states $1,...,p$ are transient and state $p+1$ is absorbing. The row vector $\bm{\pi}=(\pi_{1},...,\pi_{p})$ is the initial distribution, $\pi_{i}=\mathop{}\!\mathbb{P}(X_{0}=i)$, and $\bm{T}=\{t_{ij}\}_{i,j=1,...,p}$ where $t_{ij}$ denotes the transition rates of jumps between transient states $i$ and $j$. We assume that $\pi_{1}+\cdots\pi_{p}+1$, i.e. $X_{0}$ cannot start in the absorbing state which would have caused an atom at zero. The intensity matrix for $\{X_{t}\}_{t\geq 0}$ can be written as

[TABLE]

where $\bm{t}=(t_{1},...,t_{p})^{\prime}$ is a column vector of exit rates, i.e. $t_{i}$ is the rate of transition from state $i$ to the absorbing state $p+1$. We notice that $-\bm{T}\bm{e}=\bm{t}$, where $\bm{e}=(1,1,...,1)^{\prime}$ is the $p$–dimensional column vector of ones, since the row sums in $\bm{\Lambda}$ must all be zero. Hence $(\bm{\pi},\bm{T})$ fully parametrises the Markov process.

The class of Phase–type distributions is dense in the class of distributions on the positive reals, meaning that they may approximate any positive distribution arbitrarily well. On the other hand, they constitute a class of probabilistically tractable distributions, which often allows for exact solutions to complex stochastic problems, and frequently in a closed form. The theory is well developed with numerous applications in insurance risk and queueing theory (see e.g. [8] and reference therein). Phase–type distributions are light-tailed (i.e., their tail has an exponential decay), which makes them inadequate for modelling certain phenomena like insurance risks with heavy-tailed claims. Recently, Albrecher & Bladt [2] proposed an extension of the PH construction principle to time–inhomogeneous Markov processes, in which case the absorption times can also be heavy-tailed with a wide spectrum of possible tail shapes.

If $\tau\sim\mbox{PH}(\bm{\pi},\bm{T})$, then its density function is given by $f_{\tau}(x)=\bm{\pi}e^{\bm{T}x}\bm{t}$, its distribution function by $F_{\tau}(x)=1-\bm{\pi}e^{\bm{T}x}\bm{e}$ and its Laplace transform by $L_{\tau}(s)=\bm{\pi}(s\bm{I}-\bm{T})^{-1}\bm{t}$, where $\bm{I}$ denotes the identity matrix. The (fractional) moments are $\mathds{E}(\tau^{\alpha})=\Gamma(\alpha+1)\bm{\pi}(-\bm{T})^{-\alpha}\bm{e}$. For further details on Phase–type distributions we refer to [8].

3. Matrix Mittag–Leffler distributions

Let us now derive a matrix version of the Mittag-Leffler distribution by defining its Laplace transform and identifying the distribution associated to it. To this end, in view of (2) consider the function

[TABLE]

where $(\bm{\pi},\bm{T})$ is a PH generator.

Theorem 3.1**.**

$\phi(u)$* is the Laplace transform of a probability distribution.*

Proof.

Let $g(u)=u^{\alpha}$ and

[TABLE]

Then $\phi(u)=f(g(u))$. Now $\bm{T}-u^{\alpha}\bm{I}$ is a sub–intensity matrix for all $u\geq 0$ and therefore $(u^{\alpha}\bm{I}-\bm{T})^{-1}$ is a non–negative matrix (Green matrix, see [8, p.134]). Thus

[TABLE]

which has sign $(-1)^{n}$. Concerning $g$,

[TABLE]

has sign $(-1)^{j+1}$. We shall employ Faá di Bruno’s formula,

[TABLE]

where the summation is over $n$–tuples for which

[TABLE]

to determine the sign of $\phi^{(n)}(u)$. Notice that

[TABLE]

Hence all terms in the summation have the same sign $(-1)^{n}$, and we conclude that also the sum itself has sign $(-1)^{n}$, i.e. $\mbox{sign}(\phi^{(n)}(u))=(-1)^{n}$. Since $\phi(0)=\bm{\pi}(-\bm{T})\bm{t}=\bm{\pi}(-\bm{T})(-\bm{T}\bm{e})=\bm{\pi}\bm{e}=1$, the result then follows with Bernstein’s theorem (see [15, p.439]). ∎

Remark 3.2**.**

Note that the proof does not rely on the special form of $u^{\alpha}$, and it follows that if $g$ is any function with $g(0)=0$ and $-g$ completely monotone, then

[TABLE]

is the Laplace transform of a probability distribution, which may also be useful in other contexts.

Theorem 3.3**.**

Let $X$ be a random variable with Laplace transform (4). Then the density function of $X$ is given by

[TABLE]

Proof.

We show that $f$ has the required Laplace transform. For $u$ sufficiently large one has

[TABLE]

The result for all $u\in[0,\infty)$ then follows by analytic continuation. ∎

Definition 3.4**.**

Let $(\bm{\pi},\bm{T})$ be a PH generator and let $0<\alpha\leq 1$. A random variable $X$ is said to have a matrix Mittag–Leffler distribution, if its Laplace transform is given by

[TABLE]

In this case we write $X\sim\mbox{MML}(\alpha,\bm{\pi},\bm{T})$.

Remark 3.5**.**

The proof Theorem 3.3 shows that for any triplet $(\bm{\pi},\bm{T},\bm{t})$ for which

[TABLE]

is a density function, the function $\phi(u)=\bm{\pi}(u^{\alpha}\bm{I}-\bm{T})^{-1}\bm{t}$ is indeed a Laplace transform of a probability distribution. Distributions with density (5) are referred to as matrix–exponential distributions, and contain the class of phase–type distributions as a strict subset. We have stated the above definition in terms of a phase–type generator, but it is hence clear that the construction also applies to any matrix–exponential distribution. While most results in the following could be stated in terms of a matrix–exponential triplet, for the sake of simplicity and notation we restrict the parameters to phase–type generators only.

Corollary 3.6**.**

Let $X\sim\mbox{MML}(\alpha,\bm{\pi},\bm{T})$. Then the cumulative distribution function for $X$ is given by

[TABLE]

Proof.

The derivative of this function is

[TABLE]

which indeed coincides with the density. It remains to note that $F$ satisfies the boundary condition $F(0)=0$. ∎

One can now realize that the matrix ML distribution provides an extension of representation (3) in the following way:

Theorem 3.7**.**

Let $X\sim\mbox{MML}(\alpha,\bm{\pi},\bm{T})$. Then

[TABLE]

where $W\sim\mbox{PH}(\bm{\pi},\>\bm{T})$, and $S_{\alpha}$ is an independent (positive stable) random variable with Laplace transform given by $\exp(-u^{\alpha})$.

Proof.

Simply note that

[TABLE]

∎

Corollary 3.8**.**

Let $X\sim\mbox{MML}(\alpha,\bm{\pi},\bm{T})$. The fractional moments of order $\rho<\alpha\leq 1$ are given by

[TABLE]

Proof.

It is known [28] that the fractional moments of a random variable $S_{\alpha}$ with Laplace transform $\exp(-u^{\alpha})$ are given by

[TABLE]

From [8] we know that the $\nu$th fractional moment of a random variable $W$ with PH($\bm{\pi},\>\bm{T}$) distribution is given by

[TABLE]

By Theorem 3.7, and setting $\nu=\rho/\alpha$, $X$ will have the $\rho$th fractional moment given by

[TABLE]

∎

Remark 3.9**.**

Representation (6) is not only useful for establshing closed-form formulas, but it also suggests a simple and efficient simulation technique for $X$. Simulation algorithms for PH and stable distributions are for instance available in the statistical software R via the packages actuar and stabledist, respectively.

Example 3.1**.**

Consider $X\sim\mbox{MML}(\alpha,\bm{\pi},\bm{T})$ where $(\bm{\pi},\bm{T})$ is the PH representation of an Erlang distribution with $p$ ($p\in{\mathbb{N}}$) stages and intensity $\lambda$, i.e. $\bm{\pi}=(1,0,...,0)$ and

[TABLE]

In this case

[TABLE]

so

[TABLE]

where $\gamma$ is a simple path enclosing $-\lambda x^{\alpha}$, and where in the last step we used the residue theorem. Note that this corresponds to the GML random variable given already for general shape parameter $p=\beta\in{\mathbb{R}}_{+}$ in (3) (although this explicit form of the density was not given in [19]). Figure 1 depicts the density for several choices of parameters. The parameter $\alpha$ controls the heaviness of the tail (and more generally the deviation from exponentiality of the Mittag-Leffler function), the parameter $p$ determines the shape of the body of the distribution (since larger $p$ implies a more pronounced Erlang component), and $\lambda$ is a scaling parameter.

∎

Example 3.2**.**

Mixture of Erlang distributions form the simplest sub–class of PH distributions which are dense in the class of distributions on the positive real line (in the sense of weak convergence). Let $h$ be the density

[TABLE]

where $f_{{Erl}}(x;{p_{i}},\lambda_{i})=\lambda_{i}^{p_{i}}x^{p_{i}-1}\exp(-\lambda_{i}x)/(p_{i}-1)!$ denotes the density of an Erlang distribution and $\theta_{i}\geq 0$ are weights with $\sum\theta_{i}=1$. Then it follows immediately from Example 3.1 that the corresponding MML distribution has density

[TABLE]

In Figure 2, a trimodal distribution is considered, corresponding to $X\sim\mbox{MML}(\alpha,\bm{\pi},\bm{T})$, for a mixture of three Erlang PH components. The densities of $\log(X)$ and $X$ are both depicted. ∎

Example 3.3**.**

A Coxian phase–type distribution has a representation of the form

[TABLE]

where all $\lambda_{i}$, $i=1,\dots,p,$ are distinct. In this case

[TABLE]

so

[TABLE]

In Figure 3 four such densities are plotted.

∎

4. Sample path respresentations

We now provide two different representations of the MML distribution as sample path properties of a stochastic process. The first one will be as an absorption time of a randomly scaled time-inhomogeneous Markov jump process, and the second one as an absorption time of a particular semi-Markov process (where random time scaling is not needed).

4.1. Random time-inhomogeneous phase–type distributions

We recall from Section 2.3 that a random variable is PH distributed if it is the time until absorption of a time–homogeneous Markov jump process on a finite state–space, where one state is absorbing and the remaining states are transient. In this section we show that a MML distribution can be interpreted as a time–inhomogeneous phase–type distribution with random intensity matrix.

Let us define a random time-inhomogeneous Markov jump process $X_{t}$ as a jump process with state space $E=\{1,2,\dots,p,p+1\}$, where $p+1$ is an absorbing state, the remaining states being transient, and the intensity matrix given by

[TABLE]

where $\bm{t}(t)=-\bm{T}(t)\bm{e}$, $\mathop{}\!\bm{e}=(1,1,\dots,1)^{T}$, $\bm{0}=(0,0,\dots,0)$, $\mathop{}\!\bm{\pi}=(\pi_{1},\dots,\pi_{p})$,

[TABLE]

and the independent, positive random variable $Y$ is a random scaling factor. For the resulting absorption time

[TABLE]

we write $\tau\sim\mbox{RIPH}(Y,\mathop{}\!\bm{\pi},\mathop{}\!\bm{T}(t))$. Note that the special case $Y\equiv 1$ corresponds to the IPH class in [2].

If furthermore we can write $\mathop{}\!\bm{T}(t)=\lambda(t)\mathop{}\!\bm{T}$, the intensity matrix (7) takes the form

[TABLE]

in which case we write $\tau\sim\text{RIPH}(Y,\mathop{}\!\bm{\pi},\mathop{}\!\bm{T},\lambda)$. The following result is then immediate.

Theorem 4.1**.**

Let $\tau\sim\text{RIPH}(Y,\mathop{}\!\bm{\pi},\mathop{}\!\bm{T},\lambda)$ for a random variable $Y$ with density $g$. Then the density $f$ and distribution function $F$ of $\tau$ are given by

[TABLE]

Furthermore, if $\lambda(t)$ is a strictly positive function and we define $h$ by

[TABLE]

then

[TABLE]

where $\tau_{0}\sim PH(\mathop{}\!\bm{\pi},\mathop{}\!\bm{T})$.

Combining Theorem 3.7 and Theorem 4.1, the matrix Mittag-Leffler random variable $X\sim\mbox{MML}(\alpha,\bm{\pi},\bm{T})$ can hence be interpreted as a particular random scaling of a time-inhomogeneous phase-type distribution $\tau_{0}\sim PH(\mathop{}\!\bm{\pi},\mathop{}\!\bm{T})$ with $h(x)=x^{1/\alpha}$ (translating into $\lambda(t)=\alpha t^{\alpha-1}$) and heavy-tailed random scaling factor $Y=S_{\alpha}$. As we will illustrate in Section 5, this represents a particularly versatile yet simple class of random variables for fitting real data.

Remark 4.2**.**

For any random variable $W$ with cumulative distribution function $F_{W}$, it is possible to write

[TABLE]

where $E$ is a unit mean exponential random variable. For modeling purposes, the rationale in [2] can be interpreted as approximating the transformation function $h_{W}$ before-hand by some function $h$ (in absence of the knowledge of $W$) and then replacing $E$ with a general PH distribution, providing flexibility for the fit with often explicit formulas for the resulting random variable. The matrix-Pareto distributions defined in [2] are then the special case $Y\equiv 1$ and (up to a constant) $h(x)=e^{x}-1$ in (8), which entails $\lambda(t)=1/(1+t)$. The fitting in that case is particularly parsimonious for distributions ’close’ to a Pareto distribution (where the distance concept here is then inherited from the distance in the PH domain after the log-transform). The general representation (8), in contrast, allows to introduce a potential heavy-tail behavior also through the random scaling factor $Y$, providing more flexibility for the shape of the function $h(x)=x^{1/\alpha}$ (through the choice of $\alpha$) in a fitting procedure while keeping the resulting expressions tractable.

Remark 4.3**.**

The form (8) may also suggest to consider – for modelling purposes – the somewhat simpler case of a PH variable multiplied by a standard Pareto variable with tail index $\beta>0$, that is $h(x)=x$ and $f_{Y}(y)=\beta x^{-\beta-1}$, $x\geq 1$. For general $\tau_{0}\sim PH(\mathop{}\!\bm{\pi},\mathop{}\!\bm{T})$ it is straightforward to see that then

[TABLE]

where $m_{\beta}$ is the $\beta$-th moment of $\tau_{0}$, and $F_{m_{\beta}}$ is its $\beta$-th moment distribution.

For $\tau_{0}\sim\text{Exp}(\lambda)$ this simplifies to

[TABLE]

where

[TABLE]

is the Whittaker M function and $(x)_{n}$ is the Pochhammer symbol. Inserting a matrix into the third argument of this function, one may now proceed again with the matrix version of Cauchy’s formula and by Jordan decomposition, potentially giving rise to a theory similar to the one for Mittag-Leffler distributions. However, this direction is not the focus of the present paper. In addition, a key difference between the above product construction and the MML distribution will be discussed in Section 4.2 below in the context of a non-random path representation, for which the fine properties of the ML distribution play a crucial role.

4.2. Semi–Markov framework

Let $E^{\ast}=\{1,2,...,p\}$ be a state space and let $\bm{Q}=\{q_{ij}\}_{i,j\in E}$ denote a transition matrix for some Markov chain $\{Y_{n}\}_{n\in\mathbb{N}}$ defined on $E^{\ast}$. We assume that $q_{ii}=0$ for all $i$. $\{Y_{n}\}_{n\in\mathbb{N}}$ will serve as an embedded Markov chain in a Markov renewal process.

Let $\alpha\in(0,1]$ and $\lambda_{i}>0$. For $i=1,...,p,$ let $T^{i}_{n}$ be i.i.d. random variables with a Mittag–Leffler distribution $\mbox{ML}(\alpha,\lambda_{i})$, where we use the parametrisation such that the density of a generic $T^{i}$ is given by

[TABLE]

An alternative common parametrisation is obtained in terms of the parameter $\rho_{i}$ that satisfies $\rho_{i}^{-\alpha}=\lambda_{i}$.

We now construct a semi–Markov process $\{X_{t}\}_{t\geq 0}$ as follows. Let $S_{0}=0$ and

[TABLE]

Define

[TABLE]

Then $\{X_{t}\}_{t\geq 0}$ changes states according to the Markov chain $Y_{n}$, $S_{n}$ denotes the time of the $n$’th jump, and the sojourn times in states $i$ are Mittag–Leffler distributed with paramters $(\alpha,\lambda_{i})$. The construction is illustrated in Figure 4.

Define the intensity matrix $\bm{\Lambda}=\{\lambda_{ij}\}_{i=1,...,p}$ by

[TABLE]

and let

[TABLE]

Theorem 4.4**.**

We have

[TABLE]

Proof.

Conditioning on the time of the first jump, we get that

[TABLE]

Taking Laplace transforms, and using that

[TABLE]

together with $\lambda_{ik}=\lambda_{i}q_{ik}$, we get that

[TABLE]

or

[TABLE]

Now using $\lambda_{ii}=-\lambda_{i}$, we get

[TABLE]

In matrix form this amounts to

[TABLE]

which has the solution

[TABLE]

The right-hand side is the Laplace transform of ${E}_{\alpha,1}(\bm{\Lambda}t^{\alpha})$, establishing the result. ∎

Next we consider the case where $E=\{1,2,...,p,p+1\}$ and where the states $1,...,p$ are transient and state $p+1$ is absorbing (with respect to the Markov chain $\{Y_{n}\}_{n\in\mathbb{N}}$). This means that $\{Y_{n}\}_{n\in\mathbb{N}}$ has a transition matrix of the form

[TABLE]

and regarding the intensities we set $\lambda_{p+1}=0$. The matrix $\bm{\Lambda}$ then is of the form

[TABLE]

We notice the following useful result.

Lemma 4.5**.**

[TABLE]

where $\bm{t}=-\bm{T}\bm{e}$ (rows sum to zero).

Proof.

By definition,

[TABLE]

∎

Thus the restriction of $E_{\alpha,1}(\bm{\Lambda}x^{\alpha})$ to the transient states $1,...,p$ equals $E_{\alpha,1}(\bm{T}x^{\alpha})$ and is hence the sub–transition matrix between the transient states.

Theorem 4.6**.**

Let $\{X_{t}\}_{t\geq 0}$ be a semi-Markov process, where the matrix $\bm{\Lambda}$ has the form

[TABLE]

Let $\tau=\inf\{t\geq 0:\,X_{t}=p+1\}$ denote the time until absorption. Then $\tau$ has a MML($\alpha,\bm{\pi},\bm{T}$) distribution, with cumulative distribution function given by

[TABLE]

Proof.

For $E^{\ast}=\{1,2,\ldots,p\}$, the events $\{\tau>u\}$ and $\left\{X_{u}\in E^{\ast}\right\}$ coincide. Thus,

[TABLE]

∎

Remark 4.7**.**

The proofs above heavily depend on the form of the Laplace transform of the ML distribution, and its similarity with the exponential function. Note that the above construction naturally extends the definition of PH distributions as absorption times of continuous–time Markov chains, the latter being the limit case $\alpha\to 1$. In general, such a semi-Markov representation will hence not be available for other product distributions.

5. Statistical modeling using MML distributions

In this section we present some examples of MML distribution fitting to data. Let us start with an illustration of the maximum likelihood fitting performance to simulated data.

Example 5.1**.**

We simulate $300$ observations from $X\sim\mbox{MML}(\alpha,\bm{\pi},\bm{T})$, with a mixture of Erlang PH component, with parameters chosen in such a way that the log-data is trimodal. The corresponding maximum likelihood fit is depicted in Figure 5 (for visualization purposes the scale of the $x$-axis is logarithmic). The true parameters are $m=3$, $\alpha=0.9$, $p_{1}=p_{2}=p_{3}=3$, $\lambda_{1}=10,\>\lambda_{2}=1,\>\lambda_{3}=0.1$, $\theta_{1}=0.3,\>\theta_{2}=0.3\>\theta_{3}=0.4$, whereas the maximum likelihood estimator is found to be

[TABLE]

The negative log-likelihood at the fitted parameters was $1020.102$, compared to $1023.963$ at the true parameters (i.e., in the likelihood sense, the fitted model even outperforms the true model for the simulated data points). Note that the tri-modal shape of the underlying density is nicely identified here (which in pure PH fitting would typically not work as smoothly). We also provide a Hill plot of the simulated data, where the potentially heavy-tailed behavior can be clearly appreciated. $\Box$

The fact that MML distributions behave in a Pareto manner with parameter $\alpha\in(0,1]$ in the tail can be seen from (6). Since such a tail will be too heavy for most applications, we introduce a simple power transformation to gain flexibility for the tail behavior, which particularly allows lighter tails as well.

Definition 5.1**.**

Let $X\sim\mbox{MML}(\alpha,\bm{\pi},\bm{T})$. For $\nu>0$, we define

[TABLE]

*and refer to it as the class of Power-MML (PMML) distributions. *

The density function of a $\mbox{PMML}(\alpha,\bm{\pi},\bm{T},\nu)$ distribution is given by

[TABLE]

which will be needed for the maximum likelihood procedure below.

Remark 5.2**.**

The introduction of the PMML class allows for an adaptive transformation of the data during the fitting procedure. The interpretation of $\nu$ is then as the power to which the data should be taken in order for the latter to be most adequately fit by a pure MML distribution. As the number of MML components grows, the product $\alpha\nu$ is expected to estimate the tail index. However, this estimate might be far off when the matrix $\bm{T}$ is not large enough in order for the global fit to be adequate. In those cases, the power transform will tend to improve the fit of the body of the distribution, rather than the tail. When compared to the approach taken in [2] (fitting a PH density to log-transformed heavy tailed variables) one can consider the present procedure as adaptive selection of the transformation function, as opposed to fixing it to be the logarithm.

Example 5.2**.**

We consider a real-life motor third party liability (MTPL) insurance data set which was thoroughly studied in [1], mainly from a heavy-tailed perspective (referred to as ”Company A” there). The data set originally consists of $837$ observations, having the interpretation of claim sizes reported to the company during the time frame 1995-2010. The data are right-censored, and were analyzed recently in [6] using perturbed likelihood with censoring techniques. For the present purpose, we solely focus on the ultimates, which consists of imputing an expert prediction of the final claim amount for all claims which are still open, i.e. right-censored. We restrict our analysis here to the largest $800$ observations, since the inclusion of the 37 smallest claims lead to a sub-optimal fit, but are somehow irrelevant for modeling purposes. For convenience, we divided the claim sizes by $100,000$.

The heavy-tailed nature of the data suggests that using MML distributions to model the claim sizes is appropriate. Recently, in [7], a tail index of $\alpha^{-1}=0.48$ was suggested through an automated threshold selection procedure, using a novel trimming approach for the Hill estimator. Since this (or also other much rougher pre-analysis techniques like Pareto QQ-plots) suggests a finite mean, we employ the PMML distributions for the present purpose. This is in fact advised as a general procedure, since in situations where a pure MML fit is appropriate, the fitting procedure will suggest a value for $\nu$ close to 1 anyway. The maximum likelihood procedure identifies here a surprisingly simple PMML distribution as adequate, namely with a PH component just being a simple exponential random variable:

[TABLE]

More complex PH components turn out to indeed numerically degenerate into this simple model again. The resulting model density is hence given by

[TABLE]

The adequate fit can be appreciated in Figure 6. Observe that the maximum likelihood approach is concerned not only with the tail behaviour but also with adequately fitting the body of the distribution. Nonetheless, the tail index of the PMML fit is given by $(\hat{\alpha}\cdot\hat{\nu})^{-1}=0.4761427$, which is strikingly(!) close to the $0.48$ suggested in [7].

A previous approach to describe the entire data set by one model was given in [1, p.99], where a splicing point was suggested for this data set at around the $20$th largest order statistic, based purely on expert opinion. A semi-automated approach in [7] suggested splicing at the $14$th largest data point. Notice that not only does our model fit the data well and is much more parsimonious, but it also circumvents threshold or splicing point selection completely. $\Box$

Phase-type distributions are weakly dense in the set of all probability distributions on the positive real line. However, often a very large dimension of the PH distribution is needed to get a decent fit to data. Here, we show an example of how the class of PMML distributions can be used to reduce the dimension of a PH fit, thanks to the increased flexibility that the randomization with an $\alpha$-stable distribution and the power function $(\cdot)^{1/\nu}$ provide.

Example 5.3**.**

We consider $n=500$ simulated data points $X_{1},\dots,X_{n}$ following a mixture of two Erlang(40) distributions. The PH representation is of dimension $80$, with parameters

[TABLE]

Since the data is light-tailed, and PMML are heavy-tailed, we consider the transformed observations

[TABLE]

which are Pareto in the tail. We then proceed to fit a PMML distribution to the transformed data, but with a much lower matrix dimension. Concretely, we consider a mixture of two Erlang distributions of three phases each for the PH component of the PMML representation. In this way we are led to the maximum likelihood estimates

[TABLE]

The (back-transformed) fitted density is plotted in Figure 7, along with a histogram of the original PH data points. We also include a fitted density using a pure PH distribution of the same dimension and kind: a mixture of two Erlangs of three phases each. We observe how transforming the data into the heavy-tail domain, fitting a PMML and then back-transforming adds only two extra parameters and improves the estimation dramatically. Additionally, a Hill plot of the transformed data is provided, which shows that the tail index empirically could correspond to $\alpha^{-1}\approx 0.1$, such that it is necessary to use the PMML class, as opposed to only the MML. For reference, the resulting tail is $(\hat{\alpha}\cdot\hat{\nu})^{-1}=0.223427$, but here the quantification of the tail behaviour is not the main focus of the estimation, and used only qualitatively. In fact, a quick calculation shows that the true index is $\alpha^{-1}=1/50=0.02$, and it is well-known that it is a hard task to estimate tail indices in the transition area between Fréchet and Gumbel domains of attraction. The additional Hill plots of simulated paths of the estimated model in Figure 7 show how the hump at the middle of the Hill plot is not a random fluctuation, but rather a systematic feature.

Notice that we took the exponential transformation (13) of the simulated data because in this case we know that their tail is exponential and hence the transforms will be regularly varying, in accordance with the PMML distribution. In general, the exponentiation of light-tailed data does not imply that the resulting underlying distribution is regularly varying in the tail. Hence, for real data, a preliminary assessment of the tail behavior and the one of their exponential transforms is recommended.

Remark 5.3**.**

When fitting a MML or PMML distribution to data, one has to decide upon the dimension of the underlying phase–type representation. This problem arises similarly when fitting phase–type distributions to light-tailed data, and there are no generally accepted and well established methods for model selection, since PH distributions may be well overparametrised so that penalized methods such AIC or BIC indices will not work in general. The order of the PMML (or phase–type distribution) is therefore usually chosen by fitting a range of models of different dimension and then comparing the fit and likelihood values (which, as opposed to information indices, are comparable).

6. Conclusion

In this paper we define the class of matrix Mittag-Leffler distributions and derive some of its properties. We identify this class as a particular case of inhomogeneous phase-type distributions under random scaling with a stable law, which together with its power transforms is surprisingly versatile for modeling purposes. In addition, the class is shown to correspond to absorption times of semi–Markov processes with Mittag–Leffler distributed interarrival times, providing a natural extension of the phase–type construction. We illustrate with several examples that this class can simultaneously fit the main body and the tail of a distribution with remarkable accuracy in a parsimonious manner. It turns out that the flexibility of this heavy-tailed class of distributions can even make it worthwhile to transform data into the heavy-tailed domain, fitting the resulting data points and then transforming them back. It will be an interesting direction for future research to further explore the potential of such fitting procedures, both from a theoretical and practical perspective.

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