Heavy-tailed phase-type distributions: a unified approach.

A phase-type distribution is the distribution of the time until absorption in a finite state-space time-homogeneous Markov jump process, with one absorbing state and the rest being transient. These distributions are mathematically tractable and conceptually attractive to model physical phenomena due to their interpretation in terms of a hidden Markov structure. Three recent extensions of regular phase-type distributions give rise to models which allow for heavy tails: discrete- or continuous-scaling; fractional-time semi-Markov extensions; and inhomogeneous time-change of the underlying Markov process. In this paper, we present a unifying theory for heavy-tailed phase-type distributions for which all three approaches are particular cases. Our main objective is to provide useful models for heavy-tailed phase-type distributions, but any other tail behavior is also captured by our specification. We provide relevant new examples and also show how existing approaches are naturally embedded. Subsequently, two multivariate extensions are presented, inspired by the univariate construction which can be considered as a matrix version of a frailty model. We provide fully explicit EM-algorithms for all models and illustrate them using synthetic and real-life data.

Introduction

Phase-type (PH) distributions have been employed extensively in applied probability since they often provide exact and explicit solutions to complex stochastic problems. Another attractive property of PH distributions is that they form a dense class in the set of distributions in the positive half-line in the sense of weak convergence (see [Section 3.2.1] Bladt and Nielsen (2017)). However, and despite their denseness, PH distributions are always light-tailed, which may be a problem when heavy tails are present.

At least three approaches to remedy this problem have been introduced in the literature. The first one, originally introduced in Bladt et al. (2015) and called the NPH class of distributions, consists of considering PH distributions scaled by nonnegative discrete random variables, N. This construction principle has the advantage that the resulting distribution maintains the interpretation as being the absorption time of a homogeneous Markov jump process but in an infinite-dimensional state-space. This, indeed, allows for genuinely heavy tails for the resulting distribution. For instance, in Rojas-Nandayapa and Xie (2018), the authors showed that if the scaling component is unbounded (but otherwise arbitrary), then the resulting distribution is always heavy-tailed in terms of non-existent moment generating functions (see also Su and Chen (2006) for more general results). However, their different functionals are in terms of infinite-dimensional matrices, which in practice, can only be computed up to a finite number of terms. More recently, in Albrecher et al. (2021a), the authors considered continuous scaling and showed that closed-form expressions for different functionals of the resulting distributions can be obtained. They denoted this class by CPH. Another advantage of continuous scaling is that it maintains the (finite) dimension of the underlying PH.

A second approach was introduced in Albrecher et al. (2020a) by considering a time-fractional version of the underlying stochastic process dynamics, effectively moving into the semi-Markov domain. Together with subsequent multivariate extensions based on rewards (cf. Albrecher et al. (2021, 2020b)), these models were shown to be feasible models for applications such as non-life insurance modeling. More recently, Bladt (2021) showed that these models are relevant in describing lifetimes and performing the corresponding life-insurance calculations.

The third approach, introduced in Albrecher and Bladt (2019), consists of allowing the Markov jump process to be time-inhomogeneous in the construction principle of PH distributions leading to the class of inhomogeneous phase-type (IPH) distributions. An advantage of this approach is that one gains substantial flexibility on the tails: not only are heavy tails possible but also, e.g., lighter tails than exponential-decay can be obtained. Further extensions to covariate-dependent distributions can be found in Albrecher et al. (2021b), which is particularly well-suited for survival analysis applications.

Estimation of PH distributions was initially developed to calibrate such stochastic models to real-life data, and it is a well-developed topic in the literature. It is typically done via an expectation-maximization (EM) algorithm (Asmussen et al. (1996)), although other methods such as an MCMC approach have been introduced (Bladt et al. (2003)). More recent trends have moved towards considering PH-based models purely as flexible models for statistical fitting, irrespectively of their explicit and closed-form formulas. This data-driven approach is particularly attractive compared to other classical alternatives (for instance, kernel smoothing) since there is the implicit interpretation of an underlying process traversing through different states before it terminates, which is easy to justify in many application areas. Algorithms for discretely-scaled PH distributions, IPH models, and continuously-scaled PH distributions can be found, respectively, in Bladt and Rojas-Nandayapa (2018), Albrecher et al. (2020), and Albrecher et al. (2021a). To the best of the authors’ knowledge, an EM-based estimation procedure for fractional phase-type distributions (also called matrix Mittag-Leffler distributions) has not been considered before the present work, with Albrecher et al. (2020a) performing a purely numerical multi-dimensional maximum-likelihood estimation.

The primary purpose of this paper is to present a unified theory that englobes the above approaches to produce heavy-tailed phase-type distributions. The construction principle of the proposed models is simple to conceptualize and can be seen as a matrix extension of the frailty model in survival analysis. However, the flexibility of the underlying Markov structure allows for very different objects to be constructed as special cases. More precisely, we study IPH distributions with intensity matrices scaled by any nonnegative random variable. In other words, we impose both a random and a deterministic component which modify the speed at which the finite state-space is traversed by the Markov process, such that absorption times can possess any desired tail and body behavior, in particular obtaining heavy-tailed distributions. Inhomogeneous generalizations of Albrecher et al. (2021a); Rojas-Nandayapa and Xie (2018), the matrix Mittag-Leffler models of Albrecher et al. (2020a), and randomly scaled generalizations of Albrecher and Bladt (2019); Albrecher et al. (2021b) (with the possibility of missing covariates) are all comprised in this rich class.

In terms of physical interpretation, the latent variables play different roles. The underlying Markov dynamics aim to model heterogeneity by assuming that unobserved traversing of states has occurred. In contrast, the interpretation of the scaling component is closely related to the statistical concept of frailty. Recall that frailty models (see, e.g., Wienke (2010) for a comprehensive account of such models) specify a multiplicative random effect on the hazard rate of a distribution, effectively accounting for unobserved covariates in a Cox proportional hazards model. In contrast, we specify a multiplicative random effect on the intensity function of a Markov jump process. Nonetheless, since for IPH distributions, the hazard rate and intensity function are asymptotically equivalent (cf. Albrecher et al. (2021b)), the scaling variable can also be interpreted as accounting for heterogeneity or missing covariates in an asymptotically proportional hazards model.

The secondary aim of the paper is to present multivariate models based on this construction, which can be interpreted as generalizations of the shared and correlated frailty models (cf. Wienke (2010)). We derive EM algorithms for maximum-likelihood estimation of all the proposed models, which can be implemented either in full generality or by simplifying some assumptions and tailoring the methods for the specific application. For pedagogical reasons, we build up the multivariate case from the univariate one, although a top-bottom approach is also possible.

The rest of the paper is organized as follows. In Section 2, we present an overview of the class of IPH distributions and some important properties for our present purposes. In Section 3, we introduce our main univariate model, which we call scaled inhomogeneous phase-type, derive its main properties, give several parametric examples relevant for real-life applications, and propose a generalized EM algorithm for its maximum-likelihood estimation. In Section 4, we present a multivariate extension inspired by the shared frailty model and show how estimation of the proposed models can be performed via EM algorithms. In Section 5, we present a different multivariate extension, now based on the construction principle of correlated frailty models, and derive an EM algorithm for maximum-likelihood estimation. In Section 6, we present several numerical illustrations. Finally, Section 7 concludes.

Preliminaries

This section presents the relevant preliminaries on time-inhomogeneous Markov jump-processes and their absorption times. The distributions of the latter times are the building blocks for the scaled models introduced in Section 3. For distributional equality between two random variables X, Y, we use the notation , while the notation for F a distribution function, density, or acronym is understood as X following the distribution uniquely associated with F. Unless stated otherwise, equalities between random objects hold almost surely. For two real-valued functions, g, h the terminology , as is defined as . If a is not explicitly mentioned, it is assumed to be .

Let denote a time-inhomogeneous Markov jump process on the state-space , where states are transient and state is absorbing. In this way, has an intensity matrix of the form

Since is an intensity matrix, the sum of its rows is zero for any time , and so the identity holds, where is the p–dimensional column vector of ones. Moreover, the probability transition matrix of , whereis given in terms of the product integral (see Albrecher and Bladt (2019))

To avoid degeneracies, we assume that the process starts almost surely in a non-absorbing state with probabilities given by , . In vector notation, we write . In the sequel, we follow the convention that greek boldface lowercase letters are row-vectors, while roman boldface lowercase letters are column-vectors. Thus .

The main quantity of interest of such a process for our present purposes is the time taken to reach the absorbing state, denoted bywhich has an inhomogeneous phase-type distribution (cf. Albrecher and Bladt (2019)) with representation , and we write . Application of such random variables to statistical modeling is often treated for the special case , with some known nonnegative real function, known as the intensity function, and a fixed sub-intensity matrix. We adopt this approach in the present text. Thus we may simply write . The interested reader is referred to Bladt and Nielsen (2017) for a comprehensive account of the case and Albrecher and Bladt (2019) for further reading on general IPH distributions.

The restricted class of IPH distributions is nonetheless quite versatile. Whenever , then there exists a function h such thatwhere . More specifically, the relationship between h and is given byor in terms of derivatives

To make sure that Y is positive, unbounded, and almost surely finite, we require that

The density and survival function of are explicit in terms of matrix exponential formulas, and given by

The tail behavior of IPH distributions is driven by the asymptotic behavior of the function. Table 1 presents an overview of some commonly used intensities and transforms for generating parametric IPH distributions (see Bladt and Yslas (2021)). Applications and estimation can be found, for instance, in Albrecher and Bladt (2019); Albrecher et al. (2020, 2021b). Their names are inspired by the case, e.g., a matrix-Weibull distribution reduces to the regular Weibull distribution when is a matrix. In general, the additional parameters allow for more flexible modeling in the body of the distribution while preserving the same tail behavior as the scalar case.

Table 1

Some IPH distributions with their respective intensities and transforms

Distribution	λ(t)	h(z)	Parameters Domain
Matrix-Pareto	(t+η)^-1	ηexp(z)-1	η>0
Matrix-Weibull	ηt^η-1	z^1/η	η>0
Matrix-Lognormal	γ(log(t+1))^γ-1/(t+1)	exp(z^1/γ)-1	γ>1
Matrix-Loglogistic	ηt^η-1/(t^η+γ^η)	γ(exp(z)-1)^1/η	γ,η>0
Matrix-Gompertz	exp(ηt)	log(ηz+1)/η	η>0

Scaled inhomogeneous phase-type distributions

In this section, we introduce the main general specification of the paper and then derive some special cases together with a detailed analysis of their specific tail asymptotics. The central assumption underpinning our model is that an individual’s intensity function depends on an unobservable nonnegative random variable . More specifically, we focus on the case where acts multiplicatively on the intensity function, that iswhere is the baseline intensity function. If we denote by Y a random variable with intensity (2), then we have that

For the representation of these distributions, we make use of functional calculus. More specifically, if g is an analytic function and is a matrix, we can express by Cauchy’s formulawhere is the simple closed path in which encloses the eigenvalues of (cf. [Section 3.4.] Bladt and Nielsen (2017) for details).

The following result characterizes the density and survival functions of Y. In particular, observe that the asymptotic behavior of the tail of Y depends on both the shape of , the Laplace transform of , and on . In subsection 3.1, we give an in-depth asymptotic analysis of the new parametric models presented in this paper.

Proposition 3.1

Let Y be given by (3). Then we have that, for , where .

,

,

Proof

Property (1) follows fromwhere we have used functional calculus to define the Laplace transform evaluated at a matrix. Taking derivatives in the expression above yieldsfrom which (2) follows.

The following lemma shows that Y has the same distribution as the transformation of a scaled PH distribution. Such a representation is useful for simulation and for estimation, as is apparent in later sections.

Lemma 3.1

Let Y be given in terms of (3). Then, , where , independent of , and .

We now make the following formal definition of a random variable Y satisfying (3).

Definition 3.1

A random variable Y is said to have scaled inhomogeneous phase-type distribution (SIPH) with representation and scaling distribution if its survival function is given by

We write .

Remark 3.1

(Existing special cases of heavy-tailed PH models).

For and , almost surely, we obtain the class of NPH distributions introduced in Bladt et al. (2015), while for and , almost surely, we recover the CPH class in Albrecher et al. (2021a); Rojas-Nandayapa and Xie (2018).

Consider a Matrix Mittag Leffler (fractional phase-type) random variable as introduced in Albrecher et al. (2020a). Then, it can be shown thatwhere and is an independent (positive stable) random variable with Laplace transform given by , . Hence, we have that Y is SIPH distributed with and . This class of distributions is the time-fractional counterpart of PH distributions and can be interpreted as absorption times of a stochastic process that traverses through a finite number of states. The holding times of the latter are Mittag-Leffler distributed, which are regularly varying, and thus can possess abnormally large holding times compared to a Markov framework. However, the boundary case corresponds to the usual exponential holding times, and thus there is a regime-shift with respect to tail behavior.

When the scaling component degenerates to a point , we recover the class of IPH distributions. This also implies that the class of SIPH distributions, with a given and fixed intensity, is dense in the class of distributions on the positive real line. The argument is omitted, but it is a simple application of convergence through the diagonal of an array, for instance, by choosing a sequence of scalings with constant mean k and variances shrinking to zero.

Remark 3.2

Recall that for a continuous and positive random variable Y, the hazard function is given by

Sometimes, it is convenient to deal with the cumulative or integrated hazard function , which is given by

The classical frailty model in survival analysis assumes that the hazard function of an individual depends on an unobservable random variable . More specifically, it assumes that acts multiplicatively on a baseline hazard function , that is

Here, the random variable is known as the frailty. If we denote by Y the random variable with the above hazard, then the survival function of is given by

Thus, the unconditional survival function of Y is given by

Furthermore, model (4) can incorporate covariates in a similar way to the Cox’s proportional hazards model viawhere is a q-dimensional parameter row vector. Note that when the frailty degenerates to , one recovers the proportional hazards model, meaning that the frailty model generalizes the proportional hazards model. Commonly employed distributions as frailties include the Gamma and the positive stable distributions, among others.

In Albrecher et al. (2021b), it was shown that the intensity function of an IPH distribution is asymptotically equivalent to its hazard function. More specifically, we have that as with a positive constant. In particular, when , the previous asymptotic result becomes equality. It follows that the frailty model is a special case of our more general matrix specification of SIPH distributions, when .

Remark 3.3

(Incorporating regressors). As in the frailty model, we can introduce covariates into (2) via

In this case, we write to denote a random variable with above intensity. Note that the proportional intensities model introduced in Albrecher et al. (2021b) is retrieved if the scaling distribution degenerates to for all individuals. Consequently, the SIPH model is a generalization of the proportional intensities model.

In what follows, we mostly restrict ourselves to the model (2) without covariates, the extension being straightforward but somewhat distracting to the current train of thought. Moreover, we assume that is a continuous random variable unless stated otherwise.

Novel examples

Next, we present a suite of new examples that arise naturally as matrix extensions of some well-known frailty models, providing along the way some insight into the precise asymptotic behavior of the proposed models. In Appendix 1, the definitions of the different classes of heavy-tailed distributions are provided.

Example 3.4

(Gamma Scaling). Consider , , with Laplace transform

Then, the survival function of Y is given by

As for the matrix-Pareto type II laws introduced in Albrecher et al. (2021a), taking more general , , results in the same class of distributions. For this reason, we work only with . Consider now the particular case , , then

We call this the Matrix-Burr distribution.

Regarding the asymptotic behavior, we have thatwhere C is a positive constant, which follows from an eigenvalue decomposition of . The first-order precise asymptotics for the different intensities from Table 1 are provided in Table 2, where D, b, and c denote positive real-valued constants, which may change between intensities, but we write the same symbol for display purposes. Throughout the rest of this section, we use the same notational convention.

Table 2

Asymptotics for Gamma scaling

Intensity	Precise asymptotics	Class
Pareto	D(log(by))^-α	Slowly varying
Weibull	Dy^-αη	Regularly varying
Lognormal	D(log(y))^-αη	Slowly varying
Loglogistic	D(log(by))^-α	Slowly varying
Gompertz	Dexp(-by)	Exponential

(Positive stable scaling). Consider positive stable with stability parameter . Then

As a particular case, take , . Then

It was noted in Albrecher et al. (2021a) that is a PH representation. Thus, some simple calculations show that these distributions span the same class as the matrix-Weibull laws introduced in Albrecher and Bladt (2019). This is in contrast to the class of CPH distributions with stable mixing in Albrecher et al. (2021a), which only span the matrix-Weibull laws with .

Regarding their asymptotic behavior, we have

Table 3 gives the precise asymptotics for the different intensities of Table 1.

Table 3

Asymptotics for positive stable scaling

Intensity	Precise asymptotics	Class
Pareto	Dexp(-b(log(cy))^α)	Slowly varying
Weibull	Dexp(-by^-αη)	Weibull-type
Lognormal	Dexp(-b(log(y))^αη)	Slowly varying for αη<1 Regularly varying for αη=1 Lognormal-type for αη>1
Loglogistic	Dexp(-b(log(cy))^α)	Slowly varying
Gompertz	Dexp(-bexp(cy))	Gumbel

Example 3.6

(Inverse Gaussian scaling). Consider inverse Gaussian scaling with parameters and and density

Then, the corresponding Laplace transform of is given by

We take the particular case and . In this way

Thus,

Regarding the asymptotic behavior, we have that

Table 4 gives the precise asymptotics for the different intensities of Table 1.

Table 4

Asymptotics for inverse Gaussian scaling

Intensity	Precise asymptotics	Class
Pareto	Dexp(-b(log(cy))^1/2)	Slowly varying
Weibull	Dexp(-by^η/2)	Weibull-type
Lognormal	Dexp(-b(log(y))^η/2)	Slowly varying for η<2 Regularly varying for η=2 Lognormal-type for η>2
Loglogistic	Dexp(-b(log(cy))^1/2)	Slowly varying
Gompertz	Dexp(-bexp(cy))	Gumbel

Example 3.7

(PVF scaling). Consider the family of power variance function (PVF) distributions with Laplace transformwhere , and This family includes the Gamma, inverse Gaussian and the positive stable distributions as particular cases. Here we assume that , which results in

Regarding the asymptotic behavior, we have thatwhich results in the same asymptotics of Table 3 for the positive stable case, but with replaced by .

Example 3.8

(Compound Poisson scaling). Consider a compound model with i.i.d. random variables independent of N. In general, the Laplace transform of is given by

In particular, for and , we obtain

Thus,

Note that this distribution has an atom at infinity with probability , corresponding to the probability of . In survival analysis terms, this means that an individual may never experience the event of interest with such probability. Considering instead of N removes such an atom.

Example 3.9

(Discrete scaling). Assume that is a discrete random variable taking values in with corresponding probabilities , that is, , . Then,

Define the linear transformation on given by

Then, we can rewrite the survival function of Y aswhere denotes the Kronecker product, and is a column vector of ones of appropriate dimension. This can be thought of as an infinite-dimensional IPH distribution. The case recovers the class of NPH distributions introduced in Bladt et al. (2015).

Note that another approach to study the asymptotic behavior, and that is particularly convenient in the discrete scaling case, is to use the representation , so thatand employ the asymptotics of . For instance, taking , we have that is regularly varying with index (see Albrecher et al. (2021a) for details). This leads to the same asymptotic results in Table 2 for the different choices of intensities . For the discrete scaling, we could take, for instance, with Zeta distribution leading to the same asymptotic results.

As a second case, take with Weibull-type tail so that VZ has Weibull-type tail with shape parameter in (0, 1) (see Rojas-Nandayapa and Xie (2018)). Thus, the asymptotic behavior for the different intensities resemble those in Table 3.

Example 3.10

(Missing covariates in the proportional intensities model). Consider the proportional intensities model (also known as PH regression) introduced in Albrecher et al. (2021b) with vectors of observed and unobserved covariates and , respectively. Namely, the intensity is of the form

Given that the vector is unknown, the model cannot be employed in practice. However, we can assume thatis an unobserved random variable independent of . In this way, the scaled intensity model can be employed to account for the effect of omitted covariates by considering a parametric model for . Such additional random component can thus help account for additional variability observed in data that cannot be explained by a simpler model.

Parameter estimation

In order to derive an EM algorithm for SIPH distribution, we first recall the corresponding algorithm for CPH distributions in Albrecher et al. (2021a) (see Bladt and Rojas-Nandayapa (2018) for the discrete scaling case). Consider an i.i.d. sample from a CPH distributed random variable Y, which we will also denote by . Here, we assume that the scaling component belongs to a parametric family depending on the parameter vector and denote by its corresponding density. We now make the following definitions. Let be the number of times the underlying Markov jump process of Y starts in state k, the total number of transitions from state k to l until absorption, the number of times that k was the last state to be visited before absorption, and finally, let be the cumulated time that the Markov jump process spent in state k. The detailed routine for estimation of CPH distributions is given in Algorithm 1.

We now derive a generalized EM algorithm for maximum-likelihood estimation of SIPH distributions. Assume that is a nonnegative parametric function depending on the vector . Let , thenwhere . In particular, this implies that , meaning that is scaled PH distributed. Consider now an i.i.d. sample from this Y, then the EM algorithm for parameter estimation is the following.

Proposition 3.2

Algorithm 2 increases the likelihood function at each iteration. Since for fixed p, the likelihood of SIPH distributions is bounded, convergence towards a (possibly local) maximum is guaranteed.

Proof

By the change of variable theorem, we have that

Consider parameter values after the i-th iteration. Then the data log-likelihood after the i-th iteration is given by

In the -th iteration, we first obtain in 1. so that

Finally, by 2.

Remark 3.4

The optimization problemof Algorithm 2 is computationally heavy. However, observe that fewer iterations of any optimization routine are sufficient for the proof and conclusion of Proposition 3.2 to hold, and full convergence of (5) is not necessary. For instance, one step of the routine can already provide good results.

Remark 3.5

(Incorporating right-censoring). Algorithm 2 can be modified to work with censored data. We illustrate the changes by considering only the case of right-censoring since it is the most common scenario in survival analysis applications. However, left-censoring and interval-censoring can be treated by similar means. In such a case, we no longer observe but instead only that . By monotonicity of h, we have that , which can be interpreted as a censored observation of a scaled PH distributed random variable. Moreover, in Albrecher et al. (2021a) (and Bladt and Rojas-Nandayapa (2018)), a modified EM algorithm for the estimation of scaled PH distributions is presented for the case of censored observations. This means that the main change in Algorithm 2 is in step 2, where we must now compute

Estimation for fractional PH distributions

A key distinction of the matrix Mittag-Leffler distribution (or fractional PH), with respect to the other models introduced in Section 3.1, is that the transformation and the mixing distribution depend on the same parameter . This makes statistical estimation very challenging by ad-hoc methods, and thus embedding into the SIPH class is useful for this purpose. Note that the transformation parameters are different from the scaling component’s parameters for the previously presented models, and this last scenario is the central assumption in the derivation of Algorithm 2. Thus, special treatment must be taken for the estimation of matrix Mittag-Leffler distributions when seen as SIPH distributions. This is now solved by employing a modified EM algorithm, the details given in Algorithm 3.

By the same method of proof of Algorithm 2, one can show that Algorithm 3 increases the likelihood in each iteration, and hence we omit the details for brevity.

Shared scaling

This section presents a multivariate extension of SIPH distributions, inspired by the construction principle of the shared frailty model. The key idea is to think of an underlying random variable which is a common scaling factor to all the coordinates of an independent random vector, creating dependency and heavy-tailedness all at once through the same mechanism.

A class of multivariate CPH distributions

Before going into full generality, we consider the case where there is no deterministic time-transform component. This allows for a more transparent treatment with explicit formulas. Thus, consider the conditionally independent random variables given such that

Then, the joint survival function of is given bywhere denote the Kronecker sum and product, respectively. In particular, this yields the joint densitywhere and is the derivative of order d of , which can again be shown by the use of functional calculus through Cauchy’s formula. Moreover, marginally we get continuously scaled PH behavior:

Alternatively, it is easy to see that has representation , where are independent distributed random variables independent of , . Indeed,

These multivariate distributions were studied from another perspective in Furman et al. (2021), where the authors derived some properties in the context of risk management. We presently derive some probabilistic properties, provide an estimation method, and extend the class to allow for deterministic time transforms. In the next section we also allow for scaling of different components of the random vector by different (but correlated) scaling random variables. Since these distributions will be the building blocks of the more general time-inhomogeneous multivariate models presented in Section 4.3, a good understanding of the former facilitates the treatment of the latter.

Example 4.1

(Gamma scaling). Consider , , then the joint survival function of is given by

This distribution can be seen to be a matrix version of Mardia’s multivariate Pareto distribution (see Mardia et al. (1962)).

Parameter estimation: multivariate CPH distributions

We now present a generalized EM algorithm for maximum-likelihood estimation of the class of multivariate CPH distributions introduced previously. The complete data is the scaling component together with the conditionally independent Markov jump processes paths. We further assume that belongs to a parametric family depending on the vector and denote by its corresponding density.

Consider observations , , from a multivariate CPH distributed random vector, and let denote the whole data set. We also denote by and the sets of parameters and , respectively, and and to refer to the entries of and , . In order to write down the complete likelihood , we need the following definitions. For each , let be the number of times the underlying Markov jump process of starts in state k, the total number of transitions from state k to l until absorption, the number of times that k was the last state to be visited before absorption, and finally, let be the cumulated time that the Markov jump process spent in state k.

Then, the complete likelihood is given bywith corresponding log-likelihood (discarding the terms which do not depend on any parameters)

Regarding the E-step, which consists of computing the conditional expectation of the log-likelihood given the observed data, the calculations are somewhat similar to those of Albrecher et al. (2021a). We illustrate the procedure by computing the conditional expectation of the logarithmic term. Consider one (generic) data point () and let . Then

The formulas for all the other statistics are derived by similar calculations.

Concerning the M-step, consisting of maximizing the conditional expected log-likelihood in terms of the parameters, for the parameter of the scaling component we have in full generality

Regarding the PH component’s parameters, the entries of the sub-intensity matrix can be found by direct differentiation of the log-likelihood, while for the vector of initial probabilities, we can employ a Lagrange multiplier argument. We omit further details for brevity. We summarize the complete procedure in Algorithm 4.

A class of multivariate SIPH distributions

We now proceed to incorporate deterministic time-inhomogeneity into the shared scaling construction. Consider conditionally independent random variables given by

Thenandwhere , . Note that has representation , which can be seen as follows

Example 4.2

(Positive stable scaling). Take positive stable with stability parameter , then

For the particular case , , , we have

This joint distribution can be seen to be a matrix-parameter version of the multivariate Weibull distribution introduced in Manatunga and Oakes (1999).

Remark 4.1

Covariates can be incorporated into the model by assuming that the intensities are of the form

Remark 4.2

(Shared frailty model). In the shared frailty model, it is assumed that a group of individuals is conditionally independent given the frailty. In this way, the conditional joint survival function of , , is given bywhere are baseline cumulative hazards, . Thus, the joint survival function of is given by

Using thatthe above joint survival function can be rewritten as

In particular, this means that the survival copula of is an Archimedean copula. Note that the shared frailty model is a particular case of the class of multivariate SIPH distributions introduced here when .

We now study the dependence structure of multivariate SIPH distributions. When , the survival copula of is an Archimedean copula. To study the more general case, note that all the transformations presented in Table 1 are strictly increasing. This means that the copulas for models based on these intensities are the same as the ones of the models presented in Section 4.1, and thus it is enough to study the later case. Define the coefficient of upper tail dependence as

Proposition 4.1

Let be regularly varying with index . Thenwhere , .

Given the definition of our model, Proposition 1 of Section 2 in Engelke et al. (2019) yieldswhere are PH(), and

Moreover, is PH distributed with the same vector of initial probabilities and sub-intensity matrix , . This implies thatwhich now yields

Note that the resulting explicit expression for is in terms of the parameters of the PH components. For instance, when considering , the survival copula of the model can be different from the Clayton copula, for which . In Figure 1, we take the same value and plot the implicit copula of two multivariate CPH distributions, one with upper tail dependence coefficient smaller than and the other larger than , achieved solely by changing the parameters of the PH components.

Fig. 1

Simulation of implicit copulas of multivariate SIPH with (left), and multivariate SIPH with (right)

Parameter estimation: multivariate SIPH distributions

If we assume that is a parametric function depending on the vector , , and let . Then we can use that to formulate a generalized EM algorithm for maximum-likelihood estimation, which generalizes Algorithm 2 to the multivariate case.

Correlated scaling

We now extend the scaling of the sub-intensity matrix of SIPH distributions to the case where we condition on a random vector, the scaling factors being the components of such vector. We consider first the conditionally PH case, i.e. when no deterministic time-transform is present, and a scaling vector and such that the random variables are conditionally independent given with laws

Then, in full generality, the joint survival function of is given by

Consider the bivariate case. Then, using functional calculus, we have that that the joint survival function takes the explicit formwhere is the joint Laplace transform of , that is

Note that , where are independent distributed random variables, , independent of . Indeed,

Parameter estimation: correlated CPH distributions

The maximum-likelihood estimation of this class of multivariate distributions can be performed via a generalized EM algorithm. The derivation is done similarly to Algorithm 4 and thus omitted for brevity. Again, for estimation, we assume that belongs to a parametric family depending on the vector and denote by its corresponding joint density. The resulting detailed routine is provided in Algorithm 6.

Remark 5.1

This algorithm suffers from the curse of dimensionality. The integrals above must typically be computed numerically, given that explicit expressions are not available. Thus, the number of summands needed for the approximation increases rapidly with the dimension. It is also important to mention that correlated frailty models are typically employed only in the bivariate case. In such a case, the above algorithm is computationally feasible, thus its relevance.

Correlated SIPH distributions

We now introduce an analogous model to the correlated frailty model based on IPH distributions, effectively the most general of our models. Consider a multivariate random scaling component and , both in in , such that are conditionally independent given with conditional distribution

The joint survival function of is then given by

In the bivariate case, we have by simple calculations (using functional calculus) the explicit expression

Note that an alternative representation for is , where are independent PH distributed random variables independent of . The proof is akin to those of previous sections.

Now we consider a specific example with explicit joint density, namely the correlated Gamma case.

Example 5.1

(Correlated Gamma scaling). Inspired by Yashin et al. (1995), we consider such thatwhere , , , are independent. Then we have that

This yields

One typically sets and . In this way , , and .

Remark 5.2

(Estimation). Maximum-likelihood estimation can be performed via a modified EM algorithm, which is in the same form as Algorithm 5 with the only change in step 1, where we now employ Algorithm 4. We omit further details.

Remark 5.3

(Correlated frailty). The correlated frailty model assumes that the frailties of individuals are correlated and not necessarily shared. More specifically, in a bivariate correlated frailty model, the conditional joint density of is

In this way, the joint survival function of is given by

This is indeed a particular case of the correlated intensities model introduced in the present section when .

Numerical illustrations

In this section, we present some numerical illustrations of practical relevance. In the first example, we test the performance of Algorithm 3 for the estimation of matrix Mittag-Leffler distributions in a simulation study. In the second example, we consider the fitting of a SIPH distribution to a theoretical given distribution. In the third example, we fit a SIPH to a real-life insurance data set. As a final example, we perform a simulation study for a multivariate CPH distribution. In all cases, we ran the generalized EM algorithms until the changes in the successive log-likelihoods became negligible.

Matrix mittag-leffler distributions

We generated an i.i.d. sample of size 1, 000 from a matrix Mittag-Leffler distribution of 4 phases with parameters

We then fitted a matrix Mittag-Leffler distribution with the same number of phases to the resulting sample using Algorithm 3, obtaining the following parameters:

Observe that we can somewhat retrieve the parameters by keeping in mind possible permutation of states (since their labels are not relevant). Figure 2 shows that the algorithm recovers the structure of the data. Moreover, note that , which determines the heaviness of the tail, is close to the original value . As further evidence of the quality of the fit, we have that the log-likelihood of the fitted model is , while using the original distribution parameters and structure, we obtain .

Fig. 2

Histogram of log-simulated data versus density of the fitted matrix Mittag-Leffler model (left), and corresponding QQ-plot (right)

Matrix-Weibull

Algorithm 2 can be easily modified to approximate given theoretical distributions. As in the PH case (Asmussen et al. (1996)), the idea consists of considering sequences of empirical distributions with increasing sample size. For instance, if we denote by g the theoretical given density that we want to approximate, in step 1, we have that as ,

The rest of the formulas in step 1 are adapted through the same limit. Regarding step 2, we have

As a concrete example, we consider a Matrix-Weibull distribution (as introduced in Albrecher and Bladt (2019), having no random scaling component) with density functionand parameters

Then we fitted a SIPH distribution of 3 phases with baseline intensity , , and positive stable scaling. The fitted parameters are the following

The quality of the approximation is supported by Figure 3, which shows that we recover the shape of the original distribution. Moreover, the product , which determines the heaviness of the tail, can be compared with for the given theoretical model.

Fig. 3

Density of the original matrix-Weibull versus density of the fitted SIPH (left), and corresponding QQ-plot (right)

Real-life data

The Gamma-Gompertz frailty model is commonly employed for modeling human mortality at old ages (see, e.g., Missov (2013); Vaupel et al. (1979)). In the present example, we propose using SIPH distributions with Gamma scaling and Gompertz baseline intensity for modeling this type of data.

As a concrete case of study, we consider the lifetimes of the Swedish population that die in the year 2011 between ages 50-100. This data was obtained from the Human Mortality Database (HMD). We add covariate information by considering a separation between females () and males () in the population. Then we fitted a SIPH distribution of 4 phases with general Coxian structure in the PH component. The estimated parameters are

Figure 4 shows that the fitted distribution provides a reasonable model for both groups. If an even closer fit is sought, other parameters of the model need to be regressed as well.

Fig. 4

Histogram of lifetimes of the Swedish female population that died in 2011 at ages 50 to 100 versus density of the fitted SIPH (left), and corresponding plot for the male population(right)

Multivariate example

We generated an i.i.d. sample of size 2, 500 from a bivariate CPH distribution with parametersand Gamma scaling with . Note that the upper tail dependence coefficient of the theoretical model is , while the empirical estimator of the sample is . Then we fitted a bivariate CPH model of same dimensions using Algorithm 4 obtaining the parameters

Figure 5 shows that we recover the structure of both marginals. Regarding the dependence structure, this is supported by Figure 6, where we offer some contour plots. Moreover, note that the parameter that determines the heaviness of the tails of the marginals is close to the original model and that the coefficient of upper tail dependence is close to the original (and sample) one. Finally, note that the original model’s log-likelihood is , compared with for the fitted model.

Fig. 5

Histograms of log-simulated data versus densities of the fitted distribution

Fig. 6

Contour plot of the sample (left), contour plot of original distribution (center), and contour plot of fitted distribution (right)

Conclusion

We have provided a phase-type-based model which can result in non-exponential tail behavior by introducing random and deterministic transformations. The resulting model is generally tractable in terms of matrix calculus through the Laplace transform of the random component, and thus closed-form formulas allow for statistical and probabilistic treatments, for instance, for fully explicit generalized EM algorithms. In the univariate case, the current three main ways of generating heavy-tailed phase-type distributions fall into our framework, and several new models are introduced to complement the existing suite of hidden Markov models. In the multivariate case, we obtain generalizations of well-known frailty models with fully explicit densities, contrary to other approaches of multivariate phase-type distributions in the literature (in terms of rewards or copulas). We finally show the feasibility of the statistical implementation of our models using four different examples.

Heavy-tailed phase-type distributions are statistically attractive since their interpretation in terms of an underlying evolving process is natural in many domains of application which involve processes that traverse numerous states through time, for instance, human lifetimes or legal cases. With the models and algorithms provided in this paper, we aim to provide a clearer picture of the possibilities and limitations of Markov models for practitioners that require non-standard but interpretable models. A promising further direction of research for generating uni- and multivariate scaled phase-type distributions is to consider a general stochastic process as time-change, which for certain choices may provide fully explicit functionals and estimation procedures while remaining conceptually simple.