2.1. Generalized eigenvectors and Jordan normal form

Let $\mathbf{S}$ be a sub-generator matrix of dimension p. Let g denote the number of distinct eigenvalues, and let $\{-\lambda_j\}_{j \in \{1,\ldots, g\}}$ be their values. Specifically, let $-\lambda_1$ be the eigenvalue with the largest real part, which means that $-\lambda_1$ is real. Let $A_j$ be the set of indices of Jordan blocks associated with the eigenvalue $-\lambda_j$ , let $b_i$ be the order of Jordan block i, and let $B_j = \max_{i \in A_j} b_i$ be the maximal order of a Jordan block associated with the eigenvalue $-\lambda_j$ . Lastly, let $\mathbf{r}_{i,\nu}$ and $\mathbf{l}_{i,\nu}$ be respectively the right and left generalized eigenvectors of order $\nu$ of Jordan block i, so that $\mathbf{r}_{i,1}$ and $\mathbf{l}_{i,1}$ are the right and left eigenvectors. If the order of the Jordan block is $b_i = 1$ , we will suppress the last index and simply let $\mathbf{r}_i$ and $\mathbf{l}_i$ be the pair of eigenvectors. The generalized eigenvectors will be collected in the transformation matrix $\mathbf{T}$ such that

\begin{align*} \mathbf{T} &= \begin{pmatrix} \begin{pmatrix} \mathbf{r}_{1,1} & \quad \cdots & \quad \mathbf{r}_{1,b_1}\end{pmatrix} & \quad \cdots & \quad \begin{pmatrix}\mathbf{r}_{h,1} & \quad \cdots & \quad \mathbf{r}_{h,b_h} \end{pmatrix}\end{pmatrix}, \\ \mathbf{T}^{-1} &= \begin{pmatrix} \begin{pmatrix} \mathbf{l}_{1,b_1} \\ \vdots \\ \mathbf{l}_{1,1}\end{pmatrix} \\ \vdots \\ \begin{pmatrix} \mathbf{l}_{h,b_h} \\ \vdots \\ \mathbf{l}_{h,1}\end{pmatrix}\end{pmatrix},\end{align*}

where h is the number of Jordan blocks. Using this notation we can write $\mathbf{S}$ in Jordan normal form; see e.g. [Reference Gantmacher10, p 152]:

(3) \begin{align} \mathbf{S} &= \mathbf{T}\left(\sum_{j = 1}^g \sum_{i \in A_j}\left( ({-}\lambda_j) \mathbf{E}_i + \mathbf{N}_i \right)\right) \mathbf{T}^{-1}, \end{align}

where $\mathbf{E}_i$ and $\mathbf{N}_i$ are idempotent and nilpotent matrices giving the shape of each Jordan block. The matrices $\mathbf{N}_i$ have ones in their super-diagonal corresponding to the placement of the Jordan block. Example 1 illustrates the shapes of $\mathbf{E}_i$ and $\mathbf{N}_i$ .

A matrix function can be evaluated on $\mathbf{S}$ by applying it separately to all Jordan blocks; see e.g. [Reference Gantmacher10, p. 100]. Therefore we can write the matrix exponential $\exp\!\left(\mathbf{S} t \right)$ as

(4) \begin{align}\begin{split} \exp\!\left(\mathbf{S}t \right) &= \mathbf{T} \left(\sum_{j = 1}^g \sum_{i \in A_j} \sum_{\nu =1}^{b_i} \left(\dfrac{t^{\nu-1}}{(\nu-1)!} \right) \exp\!({-}\lambda_j t) \mathbf{E}_i \mathbf{N}_i^{\nu-1} \mathbf{E}_i \right)\mathbf{T}^{-1} \\ &= \sum_{j=1}^g \sum_{i \in A_j} \sum_{\nu = 1}^{B_j}\left(\sum_{\rho = 1}^{b_i -(\nu-1)} \mathbf{r}_{i,\rho} \mathbf{l}_{i,b_i + 2 -\nu - \rho}\right) \dfrac{t^{\nu-1}}{(\nu-1)!} \exp({-}\lambda_j t), \end{split} \end{align}

where the last expression is obtained by multiplying the matrices in the first expression.

Example 1. Consider the matrix $\mathbf{S}$

\begin{align*} \mathbf{S} = \begin{pmatrix} -4 & \quad 1 & \quad 2 & \quad 0 \\[1pt] 0 & \quad -2 & \quad 0 & \quad 1 \\[1pt] 0 & \quad 0 & \quad -3 & \quad 3 \\[1pt] 0 & \quad 0 & \quad 0 & \quad -3 \end{pmatrix}, \end{align*}

which has 3 distinct eigenvalues, $-\lambda_1 = -2$ , $-\lambda_2 = -3$ , and $-\lambda_3 = -4$ . The eigenvalues $-\lambda_1$ and $-\lambda_3$ have multiplicity 1, and $-\lambda_2$ has algebraic multiplicity 2 and geometric multiplicity 1. Let $A_1 =\{1\}$ , $A_2 = \{2\}$ , and $A_3 = \{3\}$ . The three right eigenvectors are then

\begin{align*} \mathbf{r}_{1,1} = \begin{pmatrix} 1/2 & \quad 1 & \quad 0 & \quad 0 \end{pmatrix}^{\prime}, \quad \mathbf{r}_{2,1} = \begin{pmatrix} 2 & \quad 0 & \quad 1 & \quad 0 \end{pmatrix}^{\prime}, \quad \mathbf{r}_{3,1} = \begin{pmatrix} 1 & \quad 0 & \quad 0 & \quad 0 \end{pmatrix}^{\prime}. \end{align*}

A generalized eigenvector $\mathbf{r}_{2,2}$ is found to be

\begin{align*} \mathbf{r}_{2,2} = \begin{pmatrix} -7/3 & \quad -1/3 & \quad 0 & \quad 1/3 \end{pmatrix}^{\prime}, \end{align*}

and thus

\begin{align*} \mathbf{T} = \begin{pmatrix} 1/2 & \quad 2 & \quad -7/3 & \quad 1 \\[1pt] 1 & \quad 0 & \quad -1/3 & \quad 0 \\[1pt] 0 & \quad 1 & \quad 0 & \quad 0 \\[1pt] 0 & \quad 0 & \quad 1/3 & \quad 0 \\ \end{pmatrix}. \end{align*}

The matrices $\mathbf{E}_i$ and $\mathbf{N}_i$ are given as

\begin{align*} \mathbf{E}_1 &= \begin{pmatrix} 1 & \quad 0 & \quad 0 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad 0 \end{pmatrix}, \quad \mathbf{E}_2 = \begin{pmatrix} 0 & \quad 0 & \quad 0 & \quad 0 \\ 0 & \quad 1 & \quad 0 & \quad 0 \\ 0 & \quad 0 & \quad 1 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad 0 \end{pmatrix}, \quad \mathbf{E}_3 \begin{pmatrix} 0 & \quad 0 & \quad 0 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad 1 \end{pmatrix}, \\ \mathbf{N}_1 & = \begin{pmatrix} 0 & \quad 0 & \quad 0 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad 0 \end{pmatrix}, \quad \mathbf{N}_2 = \begin{pmatrix} 0 & \quad 0 & \quad 0 & \quad 0 \\ 0 & \quad 0 & \quad 1 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad 0 \end{pmatrix}, \quad \mathbf{N}_3 \begin{pmatrix} 0 & \quad 0 & \quad 0 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad 0 \end{pmatrix}, \end{align*}

and using these and $\mathbf{T}$ we can write $\mathbf{S}$ in the form of Equation (4). The matrix $\mathbf{J}$ becomes

\begin{align*} \mathbf{J} = \begin{pmatrix} -\lambda_1 & \quad 0 & \quad 0 & \quad 0 \\ 0 & \quad -\lambda_2 & \quad 1 & \quad 0 \\ 0 & \quad 0 & \quad -\lambda_2 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad -\lambda_3 \end{pmatrix} = \begin{pmatrix}-2 & \quad 0 & \quad 0 & \quad 0 \\ 0 & \quad -3 & \quad 1 & \quad 0 \\ 0 & \quad 0 & \quad -3 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad -4 \end{pmatrix}. \end{align*}

2.2. A characterization of exponential phase-type representations

Theorem 2 was formulated in [Reference Bean and Green4], and it gives the necessary and sufficient conditions for a phase-type representation with an irreducible generator matrix to be that of an exponential distribution.

Theorem 2. (Bean and Green.) Let $Y \sim \text{PH}\left(\boldsymbol{\alpha},\mathbf{S}\right)$ , where $\mathbf{S}$ is an irreducible matrix. Then Y is exponential with rate $\lambda_1$ if and only if for all $\{j,\nu\} \in \{2,3,\ldots, g\} \times \{1,2,\ldots, B_j\} $

(5) \begin{align} \boldsymbol{\alpha} \mathbf{T} \left( \sum_{i \in A_j} \mathbf{E}_i \mathbf{N}_i^{\nu-1}\mathbf{E}_i\right) \mathbf{T}^{-1} \mathbf{e} = 0. \end{align}

If each Jordan block has a distinct eigenvalue, then the condition simplifies to

(6) \begin{align} \boldsymbol{\alpha} \mathbf{r}_{j,\nu} = 0 \quad \ or \ \quad \mathbf{l}_{j,B_j + 1 -\nu} \mathbf{e} = 0, \end{align}

for all $\{j,\nu\} \in \{2,3,\ldots, g\} \times \{1,2,\ldots, B_j\}$ .

In Theorem 3 we state the conditions under which an irreducible phase-type representation is that of an exponential distribution, and Theorem 2 then becomes a special case of Theorem 3.

Theorem 3. Let $Y \sim \text{PH}\left(\boldsymbol{\alpha},\mathbf{S}\right)$ , and assume that $\boldsymbol{\alpha}$ and $\mathbf{S}$ are chosen so that the representation $(\boldsymbol{\alpha},\mathbf{S})$ is irreducible. Then Y is exponentially distributed if and only if $B_1 = 1$ and

(7) \begin{align} \boldsymbol{\alpha} \mathbf{T} \left( \sum_{i \in A_j} \mathbf{E}_i \mathbf{N}_i^{\nu-1} \mathbf{E}_i\right) \mathbf{T}^{-1} \mathbf{e} = 0, \end{align}

for $\{j,\nu\} \in \{2,\ldots, g\} \times \{1,\ldots, B_j\}$ . The exponential distribution has rate $\lambda_1$ . If each Jordan block has a distinct eigenvalue, then the condition simplifies to

(8) \begin{align} \boldsymbol{\alpha} \mathbf{r}_{j,\nu} = 0 \quad \ or \ \quad \mathbf{l}_{j,B_j+1-\nu} \mathbf{e} = 0, \end{align}

for all $\{j,\nu\} \in \{2,3,\ldots, g\} \times \{1,2,\ldots, B_j\}$ .

Proof. First we prove that if $B_1 = 1$ and the condition in Equation (7) is true, then $Y \sim \text{exp}(\lambda_1)$ . If $ B_1 = 1$ we can write the survival function of Y as

\begin{align*} \mathbb{P}(Y > t) &= \boldsymbol{\alpha}\mathbf{T} \left(\sum_{j = 1}^g \sum_{i \in A_j} \sum_{\nu =1}^{B_j} \left(\dfrac{t^{\nu-1}}{(\nu-1)!} \right) \exp({-}\lambda_j t) \mathbf{E}_i \mathbf{N}_i^{\nu-1} \mathbf{E}_i \right)\mathbf{T}^{-1} \mathbf{e} \\ &= c_1 \exp({-}\lambda_1 t) + \sum_{j=2}^g \sum_{\nu = 1}^{B_j} c_{j,\nu} \left(\dfrac{t^{\nu-1}}{(\nu-1)!} \exp({-}\lambda_j t) \right), \end{align*}

with

(9) \begin{align} c_{1} &= \boldsymbol{\alpha} \left(\sum_{i \in A_1}\mathbf{r}_1 \mathbf{l}_1\right) \mathbf{e}, \nonumber\\ c_{j,\nu} &=\boldsymbol{\alpha} \mathbf{T}\left(\sum_{i \in A_j} \mathbf{E}_i \mathbf{N}_i^{\nu-1} \mathbf{E}_i \right)\mathbf{T}^{-1} \mathbf{e} . \end{align}

Clearly, the condition in Equation (7) implies that $c_{j,\nu} = 0$ for $\{j,\nu\} \in \{2,\ldots, g\} \times \{1,\ldots, B_j\}$ , and consequently the survival function becomes that of an exponential distribution with rate $\lambda_1$ :

\begin{equation*} \mathbb{P}\left(Y>t\right) = c_1 \exp({-}\lambda_1 t). \end{equation*}

This proves that Equation (7) is a sufficient condition.

Secondly, we show that if Y is exponentially distributed, then $B_1 = 1$ and the condition in Equation (7) holds. First we show two related statements:

• • If $B_1 > 1$ , then the distribution of Y cannot be exponential.

• • If $B_1 = 1$ , then the distribution of Y can only be exponential with rate $\lambda_1$ .

These statements imply that if Y is exponentially distributed, then $B_1 = 1$ and the rate of the exponential distribution is $\lambda_1$ . It then only remains to show that if Y is exponential with rate $\lambda_1$ , then the condition in Equation (7) holds.

The first two statements are proven simultaneously by a two-step procedure. Firstly we show that there is a term of the form $c_1 t^{B_1-1} \exp({-}\lambda_1 t)$ , $c_1 > 0$ , in the transition probability matrix $\mathbf{P}(t) = \exp\left(\mathbf{S}t \right)$ . Secondly we show that an irreducible representation implies that there is also a term of the form $k_1 t^{B_1-1} \exp({-}\lambda_1 t)$ , $k_1 > 0$ in the survival function of Y. For $B_1 > 1$ this implies that Y cannot be exponentially distributed, and for $B_1 = 1$ it implies that Y can only be exponentially distributed with rate $\lambda_1$ .

We can write the transition probability matrix $\mathbf{P}(t) = \exp\left(\mathbf{S} t \right)$ , using the Jordan normal form of $\mathbf{S}$ (see also Equation (4)), as

(10) \begin{align} \mathbf{P}(t) &= \exp\!\left(\mathbf{S}t \right) \nonumber \\ &= \mathbf{T} \left(\sum_{j=1}^g \sum_{i \in A_j} \sum_{\nu = 1}^{B_j} \dfrac{t^{\nu-1}}{(\nu-1)!}\exp({-}\lambda_j) \mathbf{E}_i \mathbf{N}_i^{\nu-1} \mathbf{E}_i \right) \mathbf{T}^{-1} \nonumber\\ &= \sum_{j=1}^g \sum_{i \in A_j} \sum_{\nu=1}^{B_j}\left( \sum_{\rho = 1}^{b_i -(\nu-1)} \mathbf{r}_{i,\rho} \mathbf{l}_{i,b_i + 2 -\nu - \rho} \right) \dfrac{t^{\nu-1}}{(\nu-1)!} \exp({-}\lambda_j t) \nonumber\\ &= \sum_{j=1}^g \sum_{\nu=1}^{B_j} \sum_{i \in A_j} \mathbf{M}_{i,\nu} \dfrac{t^{\nu-1}}{(\nu-1)!} \exp({-}\lambda_j t), \end{align}

where $\mathbf{M}_{i,\nu} = \sum_{\rho = 1}^{b_i -(\nu-1)} \mathbf{r}_{i,\rho} \mathbf{l}_{i,b_i + 2 -\nu - \rho}$ . For $\nu > b_i$ we have $\mathbf{M}_{i,\nu} = \mathbf{0}$ .

Consider the (l,k)’th element of the transition probability matrix

\begin{align*} \mathbf{P}(t)_{l,k} &= \exp\left(\mathbf{S} t\right)_{l,k} = \sum_{j=1}^g \sum_{\nu = 1}^{B_j} \sum_{i \in A_j}\left(\mathbf{M}_{i,\nu}\right)_{l,k} \dfrac{t^{\nu-1}}{(\nu-1)!} \exp({-}\lambda_j t) \\&= \sum_{j=1}^g \sum_{\nu=1}^{B_j} c_{j,\nu}^{l,k} \dfrac{t^{\nu-1}}{(\nu-1)!} \exp({-}\lambda_j t), \end{align*}

where

$$c_{j,\nu}^{l,k} = \sum_{\substack{i \in A_j \\\nu \leq b_i }} \left(\mathbf{M}_{i,\nu}\right)_{l,k}.$$

The function $t^{B_1-1} \exp({-}\lambda_1 t)$ dominates $t^{\nu-1} \exp({-}\lambda_j t)$ for all $\nu$ and $-\lambda_1 > \Re\left({-}\lambda_j\right)$ , and it dominates $t^{\nu-1} \exp\left({-}\lambda_1 t \right)$ for $\nu < B_1$ . Thus, if $c_{1,B_1}^{l,k} \neq 0$ , there exists a $t_0$ such that for $t> t_0$

(11) \begin{align}\begin{split}&\sum_{\nu=1}^{B_1-1} \lvert c_{1,\nu}^{l,k} \dfrac{t^{\nu-1}}{(\nu-1)!} \exp({-}\lambda_1 t) \rvert+ \sum_{j=2}^g \sum_{\nu = 1}^{B_j}\lvert c_{j,\nu}^{l,k} \dfrac{t^{\nu-1}}{(\nu-1)!} \exp({-}\lambda_j t) \rvert \\<& \lvert c_{1,B_1}^{l,k} \dfrac{t^{B_1-1}}{(B_1-1)!} \exp({-}\lambda_1 t) \rvert.\end{split}\end{align}

As the elements of the transition probability matrix $\mathbf{P}(t)$ are non-negative, Equation (11) implies that $c_{1,B_1}^{l,k} \geq 0$ for $l,k \in \{1,\ldots,p\} \times \{1,\ldots, p\}$ . Thus the matrix $\sum_{i \in A_1} \mathbf{M}_{i,B_1}$ can only have non-negative elements. Furthermore, it must have at least one element which is different from 0 and therefore positive. To see this, consider one Jordan block, $i \in A_1$ , with order $B_1$ . Then $\mathbf{M}_{i,B_1} = \mathbf{r}_{i,1} \mathbf{l}_{i,1}$ , where $\mathbf{r}_{i,1}$ is the right eigenvector associated with the block, and $\mathbf{l}_{i,1}$ is the left eigenvector (see Equation (10)). Both these vectors are different from the zero-vector and thus $\mathbf{M}_{i,B_1}$ contains at least one element different from 0. If there is only one Jordan block in $A_1$ of order $B_1$ , then $\mathbf{M}_{i,B_1} = \{c_{1,B_1}^{l,k}\}$ , and the non-zero element(s) of $\mathbf{M}_{i,B_1}$ must therefore be positive. If there are multiple Jordan blocks in $A_1$ of order $B_1$ , we need to argue that the non-zero elements do not cancel out. We consider the sum of these matrices:

\begin{align*} \sum_{\substack{i \in A_1\\ b_i = B_1}} \mathbf{M}_{i,B_1} &= \sum_{\substack{i \in A_1\\ b_i = B_1}} \mathbf{r}_{i,1} \mathbf{l}_{i,1} \\ &= \sum_{\substack{i \in A_1 \\ b_i = B_1}} \begin{pmatrix}\mathbf{r}_{i,1} \left(\mathbf{l}_{i,1}\right)_{1} &\mathbf{r}_{i,1} \left(\mathbf{l}_{i,1}\right)_{2}& \quad \cdots & \quad \mathbf{r}_{i,1} \left(\mathbf{l}_{i,1}\right)_{p} \end{pmatrix} \\ &= \begin{pmatrix} \sum_{\substack{i \in A_1 \\ b_i = B_1}} \mathbf{r}_{i,1} \left(\mathbf{l}_{i,1}\right)_1 & \sum_{\substack{i \in A_1 \\ b_i = B_1}} \mathbf{r}_{i,1} \left(\mathbf{l}_{i,1}\right)_2 & \quad \cdots & \quad \sum_{\substack{i \in A_1 \\ b_i = B_1}} \mathbf{r}_{i,1} \left(\mathbf{l}_{i,1}\right)_p \end{pmatrix}. \end{align*}

Thus each column of $ \sum_{\substack{i \in A_1\\ b_i = B_1}} \mathbf{M}_{i,B_1}$ can be expressed as a linear combination of the right eigenvectors of the Jordan blocks of order $B_1$ associated with $-\lambda_1$ . As the right eigenvectors are linearly independent, such a sum only equals the zero-vector if all the coefficients are zero. For column k the coefficients are $\{\left(\mathbf{l}_{i,1}\right)_k\}$ , and as $\mathbf{l}_{i,1}$ is a left eigenvector, it is different from the zero-vector. Thus there must exist at least one column with $\{\left(\mathbf{l}_{i,1}\right)_k\} \neq 0$ , and thus $\sum_{\substack{i \in A_1\\ b_i = B_1}} \mathbf{M}_{i,B_1}$ has at least one non-zero element, which by our previous conclusion must be positive. Denote the position of this element by $\{l^*,k^*\}$ . We can then write element $\{l^*,k^*\}$ of the transition probability matrix $\mathbf{P}(t)$ as

\begin{align*} \left(\mathbf{P}(t)\right)_{l^*,k^*} &= c_{1,B_1}^{l^*,k^*} \dfrac{t^{B_1-1}}{(B_1-1)!} \exp({-}\lambda_1 t) + \sum_{\nu=1}^{B_1-1} c_{1,\nu}^{l^*,k^*} \dfrac{t^{\nu-1}}{(\nu-1)!} \exp({-}\lambda_1 t) +\\ &\sum_{j=2}^g \sum_{\nu=1}^{B_j} c_{j,\nu}^{l^*,k^*} \dfrac{t^{\nu-1}}{(\nu-1)!} \exp({-}\lambda_j t), \end{align*}

where $c_{1,B_1}^{l^*,k^*} > 0$ . We have hence proven the existence of a term of the form $c_{1,B_1} t^{B_1}/(B_1-1)! \exp({-}\lambda_1t)$ in the transition probability matrix $\mathbf{P}\left(t\right)$ , which was the first of the two steps.

We now move on to show that the survival function of Y, $\mathbb{P}(Y>t)$ , contains a term of the form $k_{1,B_1} \dfrac{t^{B_1-1}}{(B_1-1)!} \exp({-}\lambda_1 t)$ , $k_{1,B_1} > 0$ . Consider a fixed point in time $u > 0$ , and let $\mathbf{q}(u) = \{q_i(u)\}$ be the distribution across the states at time u. By Equation (1) we have

\begin{align*} \mathbf{q}(u) = \boldsymbol{\alpha} \exp\left(\mathbf{S} u\right). \end{align*}

As the initial vector $\boldsymbol{\alpha}$ is chosen so that the representation $(\boldsymbol{\alpha},\mathbf{S})$ is irreducible, all the elements of $\mathbf{q}(u)$ are strictly positive. Consider now the survival function of Y at time $u+t$ ,

\begin{align*} \mathbb{P}\left(Y > u + t \right) &= \boldsymbol{\alpha} \exp\left(\mathbf{S}(u+t) \right) \mathbf{e} \\ &= \boldsymbol{\alpha} \exp\left(\mathbf{S} u \right) \exp\left(\mathbf{S} t \right) \mathbf{e} \\ &= \mathbf{q}(u) \exp\left(\mathbf{S}t \right) \mathbf{e} \\ &= \sum_{l= 1}^p \sum_{k=1}^p \sum_{j =1}^g \sum_{\nu=1}^{B_j} q_{l}(u) c_{j,\nu}^{l,k} \dfrac{t^{\nu-1}}{(\nu-1)!} \exp\left({-}\lambda_j t\right)\\ &= \sum_{j=1}^g \sum_{\nu = 1}^{B_j} \tilde{c}_{j,\nu} \dfrac{t^{\nu-1}}{(\nu-1)!} \exp\left({-}\lambda_j t\right). \end{align*}

Specifically, we have

$$\tilde{c}_{1,B_1} = \sum_{l =1}^p \sum_{k=1}^p q_l(u) c_{1,B_1}^{l,k}> 0,$$

which implies that the survival function of Y has a term of the form

$$k_{1,B_1} \dfrac{t^{B_1-1}}{(B_1-1)!} \exp({-}\lambda_1 t).$$

This was the second of the two steps.

It now follows that Y cannot follow an exponential distribution if $B_1 > 1$ , and if $B_1 = 1$ then Y can only be exponentially distributed if the rate is $\lambda_1$ . By the linear independence of the exponential functions $t^{\nu}/(\nu-1)! \exp({-}\lambda_j t)$ for different values of $\{j,\nu\}$ , this happens specifically when $c_{j,\nu} = 0$ for $\{j,\nu\} \in \{2,\ldots, p\} \times \{1,\ldots, B_j\}$ . This is equivalent to the condition in Equation (7) by Equation (9). The simplification from the condition in Equation (7) to the condition in Equation (8) follows by the same arguments as used in [Reference Bean and Green4].

In [Reference Bean and Green4] it was shown that the condition in Equation (5) simplifies significantly if $\mathbf{S}$ is diagonalizable:

(12) \begin{align} \boldsymbol{\alpha} \sum_{i \in A_j} \mathbf{r}_i \mathbf{l}_i \mathbf{e} = 0, \quad j \in \{2,\ldots, g\}. \end{align}

Furthermore, if the eigenvalues are all distinct, the condition (6) simplifies to

(13) \begin{align} \boldsymbol{\alpha} \mathbf{r}_j = 0 \quad \text{ or } \quad \mathbf{l}_j \mathbf{e} = 0, \quad j \in \{2,\ldots, p\}. \end{align}

The arguments used to show the simplification do not depend on $\mathbf{S}$ being irreducible, and therefore the conditions in Equation (7) and Equation (8) can also be simplified to the conditions in, respectively, Equation (12) and Equation (13).

The generalization from an irreducible generator matrix in Theorem 2 to an irreducible representation in Theorem 3 is important. As already touched upon, many phase-type representations do not have an irreducible generator matrix, e.g. a Coxian generator or a hyper-exponential generator. By Equation (2), an irreducible representation is a representation in which there is a positive probability of visiting all the states. All phase-type representations can be made irreducible without altering the distribution, as one can simply remove the states that cannot be reached. Theorem 3 therefore covers all phase-type representations of practical use. Reducible representations sometimes occur as a result of certain operations on phase-type distributions (e.g. conditioning on a certain exit state). This is a mere curiosity, as Theorem 3 in that case applies to the irreducible part of the representation.

3.1. Eigenvectors

Theorem 1 is a corollary of Theorem 3 (and if $\mathbf{S}$ is irreducible it is a corollary of Theorem 2). If $\boldsymbol{\alpha} \mathbf{S} = -\lambda_1 \boldsymbol{\alpha}$ , then $\boldsymbol{\alpha} \mathbf{r}_{j,\nu} = 0$ for $\{j,\nu\} \in \{2,\ldots, g\} \times \{1,\ldots, B_j\}$ , because $\boldsymbol{\alpha} = \mathbf{l}_i$ for $i \in A_1$ . Similarly, if $\mathbf{S}\mathbf{e} = -\lambda_1 \mathbf{e}$ , then $\mathbf{l}_{j,\nu} \mathbf{e} = 0$ for $\{j,\nu\} \in \{2,\ldots, g\} \times \{1,\ldots, B_j\}$ . In both cases Equation (8) is immediately seen to be satisfied. The case of $\boldsymbol{\alpha}$ or $\mathbf{e}$ being an eigenvector is therefore a special case, in which the same condition in Equation (8) is satisfied for all choices of $\{j,\nu\} \in \{2,\ldots, g\} \times \{1,\ldots, B_j\}$ .

For a phase-type representation with $\mathbf{S} \mathbf{e} = -\lambda \mathbf{e}$ , the exit rate is constant, and this implies that the hazard rate is constant. A constant hazard rate is a characterizing property of exponential distributions. Let $\lambda(t)$ be the hazard rate at time t. We then have $\lambda(t) = f(t)/S(t) = \boldsymbol{\alpha}\exp(\mathbf{S}t) \mathbf{s}/\left(\boldsymbol{\alpha} \exp(\mathbf{S}t) \mathbf{e}\right) = -\lambda$ . For $\boldsymbol{\alpha} \mathbf{S} = -\lambda \boldsymbol{\alpha}$ the distribution on the transient states does not change its shape as time progresses, which also implies that the hazard-rate is constant. Let $\mathbf{q}(t) = \{q_i(t)\}_{i \in \{1,2, \ldots, p\}}$ be the distribution on the transient states at time t, that is, $\mathbb{P}\left(X_t = i\right) = q_i(t)$ . By Equation (1) we have $q_i(t) =\boldsymbol{\alpha} \exp\left(\mathbf{S} t \right) \mathbf{e}_i = \alpha_i\exp({-}\lambda t)$ , so

\begin{align*} \mathbf{q}(t) =\boldsymbol{\alpha} \exp({-}\lambda t) .\end{align*}

Such a distribution is also sometimes called quasi-stationary (see e.g. [Reference Bean3]), and it has constant hazard rate

\begin{align*} \lambda(t) = \dfrac{\boldsymbol{\alpha}\exp\left(\mathbf{S} t \right) \mathbf{s}}{\boldsymbol{\alpha}\exp\left(\mathbf{S} t \right)\mathbf{e}} = \dfrac{\sum_{i=1}^p\alpha_i \exp({-}\lambda t) s_i}{\boldsymbol{\alpha} \exp({-}\lambda t) \mathbf{e}}= \sum_{i=1}^p \alpha_i s_i.\end{align*}

Interestingly, the two cases $\boldsymbol{\alpha}\mathbf{S} = -\lambda \boldsymbol{\alpha}$ and $\mathbf{S} \mathbf{e} = -\lambda \mathbf{e}$ can be linked through the concept of time-reversal, which we will explore in the next subsection.

3.2. Time-reversal

A representation of a phase-type distribution $(\boldsymbol{\alpha},\mathbf{S})$ can be time-reversed into the representation denoted by $(\tilde{\boldsymbol{\alpha}},\tilde{\mathbf{S}})$ ; see [Reference Ramaswami21]. The distribution remains the same, but conceptually the states are visited in reverse order, i.e. the state from which absorption happens becomes the initial state and vice versa. There is a relationship between the (generalized) left (resp. right) eigenvectors of the original representation and the (generalized) right (resp. left) eigenvectors of the time-reversed representation. From this relationship we find that the key operation to link phase-type representations with $\boldsymbol{\alpha} \mathbf{S} = -\lambda \boldsymbol{\alpha}$ and $\mathbf{S}\mathbf{e} = -\lambda \mathbf{e}$ is time-reversal.

Definition 2. (Time-reversed representation.) Consider a representation $(\boldsymbol{\alpha},\mathbf{S})$ for a phase-type distribution. The time-reversed representation $(\tilde{\boldsymbol{\alpha}},\tilde{\mathbf{S}})$ is then

\begin{align*} \tilde{\boldsymbol{\alpha}} &= \dfrac{\left( \Delta(\boldsymbol{\pi}) \mathbf{s}\right)^{\prime} }{\left( \Delta(\boldsymbol{\pi}) \mathbf{s}\right)^{\prime} \mathbf{e}} = \dfrac{\left(\Delta\left( \boldsymbol{\pi}\right) \mathbf{s}\right)^{\prime}}{\boldsymbol{\pi} \mathbf{s}},\\ \tilde{\mathbf{S}} &= \Delta(\boldsymbol{\pi})^{-1} \mathbf{S}^{\prime} \Delta(\boldsymbol{\pi}),\end{align*}

where $\boldsymbol{\pi}$ is a stationary measure of the renewal process with $Y_i \sim \text{PH}\left(\boldsymbol{\alpha},\mathbf{S}\right)$ as inter-epoch times,

\begin{equation*} \boldsymbol{\pi} = \boldsymbol{\alpha}({-}\mathbf{S})^{-1},\end{equation*}

and $\Delta(\cdot)$ is the operator that creates a diagonal matrix with the input vector in the diagonal.

Theorem 4. Let $\left(\boldsymbol{\alpha},\mathbf{S}\right)$ be a phase-type representation for Y, and let $\left(\tilde{\boldsymbol{\alpha}},\tilde{\mathbf{S}} \right)$ be the time-reversed representation. If $\mathbf{l}_{j,\nu}$ is a left generalized eigenvector of order $\nu$ with eigenvalue $-\lambda_j$ for the matrix $\mathbf{S}$ , then the vector

(14) \begin{align} \tilde{\mathbf{r}}_{j,\nu} = \left( \mathbf{l}_{j,\nu}\Delta\left(\boldsymbol{\pi}\right)^{-1}\right)^{\prime} \end{align}

is a right generalized eigenvector of order $\nu$ with eigenvalue $-\lambda_j$ for the matrix $\tilde{\mathbf{S}}$ . Similarly, if $\mathbf{r}_{j,\nu}$ is a right generalized eigenvector of order $\nu$ with eigenvalue $-\lambda_j$ for the matrix $\mathbf{S}$ , then the vector

(15) \begin{align} \tilde{\mathbf{l}}_{j,\nu} = \left(\Delta\left(\boldsymbol{\pi}\right) \mathbf{r}_{j,\nu}\right)^{\prime} \end{align}

is a left generalized eigenvector of order $\nu$ with eigenvalue $-\lambda_j$ for the matrix $\tilde{\mathbf{S}}$ .

Proof. We prove the claim by induction up to $\nu = b_j$ , where $b_j$ is the order of the relevant Jordan block. Consider first $\nu =1$ . Then $\mathbf{l}_{j,1}$ and $\mathbf{r}_{j,1}$ are the left and right eigenvectors of $\mathbf{S}$ . We show that $\tilde{\mathbf{r}}_{j,1}$ given in Equation (14) is a right eigenvector of $\tilde{\mathbf{S}}$ by direct calculation:

\begin{align*} \tilde{\mathbf{S}} \tilde{\mathbf{r}}_{j,1} &= \Delta(\boldsymbol{\pi})^{-1} \mathbf{S}^{\prime} \Delta(\boldsymbol{\pi}) \left(\mathbf{l}_{j,1} \Delta(\boldsymbol{\pi})^{-1}\right)^{\prime} = \Delta(\boldsymbol{\pi})^{-1} \mathbf{S}^{\prime} \mathbf{l}_{j,1}^{\prime}= -\lambda_j \Delta(\boldsymbol{\pi})^{-1} \mathbf{l}^{\prime}_{j,1} = - \lambda_j \tilde{\mathbf{r}}_{j,1}. \end{align*}

Similarly, the vector $\tilde{\mathbf{l}}_{j,1}$ given in Equation (15) is a left eigenvector of $\tilde{\mathbf{S}}$ :

\begin{align*} \tilde{\mathbf{l}}_{j,1} \tilde{\mathbf{S}} &= \left( \Delta(\boldsymbol{\pi}) \mathbf{r}_{j,1}\right)^{\prime}\Delta(\boldsymbol{\pi})^{-1} \mathbf{S}^{\prime} \Delta(\boldsymbol{\pi}) = \mathbf{r}^{\prime}_{j,1} \mathbf{S}^{\prime} \Delta(\boldsymbol{\pi}) = -\lambda_j \mathbf{r}^{\prime}_{j,1} \Delta(\boldsymbol{\pi}) = -\lambda_j \tilde{\mathbf{l}}_{j,1}. \end{align*}

Assume now that $\tilde{\mathbf{r}}_{j,\nu}$ given in Equation (14) is a generalized right eigenvector of order $\nu \leq b_j-1$ with eigenvalue $-\lambda_j$ for the matrix $\tilde{\mathbf{S}}$ . In the induction step we show that $\tilde{\mathbf{r}}_{j,\nu+1}$ is a generalized right eigenvector of order $\nu +1$ :

\begin{align*} \tilde{\mathbf{S}} \mathbf{r}_{j,\nu+1} &= \Delta\left(\boldsymbol{\pi}\right)^{-1} \mathbf{S}^{\prime} \Delta\left(\boldsymbol{\pi}\right) \left(\mathbf{l}_{j,\nu +1} \Delta(\boldsymbol{\pi})^{-1}\right)^{\prime} = \Delta\left(\boldsymbol{\pi}\right)^{-1} \mathbf{S}^{\prime} \mathbf{l}^{\prime}_{j,\nu +1 } = \Delta\left(\boldsymbol{\pi}\right)^{-1} \left({-} \lambda_j\mathbf{l}^{\prime}_{j,\nu+1} + \mathbf{l}^{\prime}_{j,\nu} \right) \\ &= - \lambda_j\tilde{\mathbf{r}}_{j,\nu +1} + \tilde{\mathbf{r}}_{j,\nu}. \end{align*}

Finally, assume that $\tilde{\mathbf{l}}_{j,\nu}$ given in Equation (15) is a generalized left eigenvector of order $\nu \leq b_j -1$ with eigenvalue $-\lambda_j$ for the matrix $\tilde{\mathbf{S}}$ . We then have

\begin{align*} \tilde{\mathbf{l}}_{j,\nu+1} \tilde{\mathbf{S}} &= \left( \Delta(\boldsymbol{\pi}) \mathbf{r}_{j,\nu+1}\right)^{\prime}\Delta(\boldsymbol{\pi})^{-1} \mathbf{S}^{\prime} \Delta(\boldsymbol{\pi}) = \mathbf{r}^{\prime}_{j,\nu+1} \mathbf{S}^{\prime} \Delta(\boldsymbol{\pi}) = \left({-}\lambda_j \mathbf{r}^{\prime}_{j,\nu +1 } + \mathbf{r}_{j,\nu}\right) \Delta(\boldsymbol{\pi}) \\ &= -\lambda_j \tilde{\mathbf{l}}_{j,\nu +1 } + \tilde{\mathbf{l}}_{j,\nu}, \end{align*}

which finishes the induction step and therefore the proof.

From Theorem 4 we obtain the following corollary.

Proof. Assume first that $\boldsymbol{\alpha} \mathbf{S} = - \lambda \boldsymbol{\alpha}$ . By Theorem 4 the vector $\left(\boldsymbol{\alpha} \Delta(\boldsymbol{\pi})^{-1}\right)^{\prime}$ is a right eigenvector of $\tilde{\mathbf{S}}$ . Let $\pi_i$ be the elements of $\boldsymbol{\pi} = \boldsymbol{\alpha}\left({-}\mathbf{S}\right)^{-1}$ . As $\boldsymbol{\alpha} \mathbf{S} = -\lambda \boldsymbol{\alpha}$ , we have $\pi_i = -\alpha_i/\lambda$ , where $\alpha_i$ are the elements of $\boldsymbol{\alpha}$ . Consequently $\left(\boldsymbol{\alpha} \Delta(\boldsymbol{\pi})^{-1}\right)^{\prime} = -\lambda\mathbf{e}$ , which is a scaled version of $\mathbf{e}$ . Thus $\mathbf{e}$ is also an eigenvector of $\tilde{\mathbf{S}}$ .

To show the second direction, let $ \left(\tilde{\tilde{\boldsymbol{\alpha}}},\tilde{\tilde{\mathbf{S}}} \right)$ be the time-reversed representation of $\left(\tilde{\boldsymbol{\alpha}},\tilde{\mathbf{S}}\right)$ . It then holds that $(\tilde{\tilde{\boldsymbol{\alpha}}},\tilde{\tilde{\mathbf{S}}}) = \left(\boldsymbol{\alpha},\mathbf{S}\right)$ , and therefore we only need to show $\tilde{\mathbf{S}}\mathbf{e} = -\lambda \mathbf{e} \Rightarrow\tilde{\tilde{\boldsymbol{\alpha}}} \tilde{\tilde{\mathbf{S}}} = -\lambda \tilde{\tilde{\boldsymbol{\alpha}}}$ . Let $\tilde{\boldsymbol{\pi}} = \tilde{\boldsymbol{\alpha}} ({-}\tilde{\mathbf{S}})^{-1}$ . Assume now $\tilde{\mathbf{S}} \mathbf{e} = -\lambda \mathbf{e}$ . By Theorem 4 the vector $\left(\Delta(\tilde{\boldsymbol{\pi}}) \mathbf{e}\right)^{\prime} = \tilde{\boldsymbol{\pi}}$ is a left eigenvector of $\tilde{\tilde{\mathbf{S}}}$ , and by Definition 2 we have

\begin{align*} \tilde{\tilde{\boldsymbol{\alpha}}} = \dfrac{\big(\Delta(\tilde{\boldsymbol{\pi}})\big({-}\tilde{\mathbf{S}}\mathbf{e}\big)\big)^{\prime}}{\big(\Delta(\tilde{\boldsymbol{\pi}})\big({-}\tilde{\mathbf{S}}\mathbf{e}\big)\big)^{\prime} \mathbf{e}} = \dfrac{\left(\Delta\left(\tilde{\boldsymbol{\pi}} \right) \left( \lambda \mathbf{e}\right) \right)^{\prime}}{c} = \dfrac{\lambda}{c} \left(\Delta(\tilde{\boldsymbol{\pi}}) \mathbf{e}\right)^{\prime} = \dfrac{\lambda}{c} \tilde{\boldsymbol{\pi}}. \end{align*}

The vector $\tilde{\tilde{\boldsymbol{\alpha}}}$ is then a normalized version of the left eigenvector $\left(\Delta(\tilde{\boldsymbol{\pi}}) \mathbf{e}\right)^{\prime} = \tilde{\boldsymbol{\pi}}$ and therefore is itself a left eigenvector.

Conceptually the content of Theorem 4 makes sense, as the left eigenvectors of a matrix $\mathbf{M}$ are the right eigenvectors of $\mathbf{M}^{\prime}$ . By Definition 2, $\tilde{\boldsymbol{S}}$ involves the transpose of $\mathbf{S}$ and some kind of normalization through $\Delta(\boldsymbol{\pi})$ . Corollary 1 lets us understand a Markov jump process with $\boldsymbol{\alpha}\mathbf{S} = -\lambda \boldsymbol{\alpha}$ as the time-reversal of a process with constant exit rate. For $\mathbf{S} \mathbf{e} = -\lambda \mathbf{e}$ , we have that $\boldsymbol{\pi}$ is the distribution of the states from which absorption happens, while it is the initial distribution for $\boldsymbol{\alpha} \mathbf{S} = -\lambda \boldsymbol{\alpha}$ .

In [Reference Bladt and Nielsen5,Reference Buchholz7,Reference He, Zhang and Vera12], time-reversed representations of a PH-distribution with $\mathbf{S}\mathbf{e} = -\lambda \mathbf{e}$ are studied, but the connection in Corollary 1 is not noted or proven there.

3.3. PH-simplicity

In [Reference O’Cinneide17, Definition 1], the concept of PH-simplicity was defined in the following way.

PH-simplicity is linked to the left eigenvectors of $\mathbf{S}$ , as stated in [Reference O’Cinneide17, Theorem 2].

A non-simple generator matrix $\mathbf{S}$ can be thought of as over-parameterized, because multiple choices of $\boldsymbol{\alpha}$ lead to the same distribution, and by Theorem 5 this happens if and only if at least one left eigenvector is orthogonal to $\mathbf{e}$ . If $\mathbf{S}$ is diagonalizable, this means that the terms related to the eigenvalue associated with the orthogonal eigenvector are removed from the distribution, and there are thus more states in the phase-type representation than there are linearly independent exponential functions in the survival function, which leaves some room for obtaining the same distribution in multiple ways.

For irreducible generator matrices the concept of simplicity is closely tied to the dimension of the space of initial distributions, which makes the corresponding phase-type distributions exponential. For illustration, we consider the case of p distinct eigenvalues, which also implies that $\mathbf{S}$ is diagonalizable. If $\mathbf{l}_{k} \mathbf{e} = 0$ for some $k \in \{2,\ldots, p \}$ , then Y is exponential by Equation (13) if

\begin{align*} \boldsymbol{\alpha}\mathbf{r}_j = 0, \quad \text{for } j \in \{2,\ldots, p\}, \, j\neq k.\end{align*}

This is satisfied for $\boldsymbol{\alpha} = c_1\mathbf{l}_1 + c_k \mathbf{l}_k$ , since left and right eigenvectors of distinct eigenvalues are orthogonal to each other. If $\mathbf{l}_k$ only contains real elements, we can always choose $c_1 > 0$ and $c_k \neq 0$ such that $\boldsymbol{\alpha}$ only contains positive and real elements. This is because $\mathbf{S}$ is irreducible, so $\mathbf{l}_1$ only contains strictly positive and real elements. The space where we can choose $\boldsymbol{\alpha}$ such that the phase-type distribution is exponential is thus a (restricted) span of $(\mathbf{l}_1,\mathbf{l}_k)$ . If $\mathbf{l}_k$ has complex elements, then $\bar{\mathbf{l}}_k$ is also a left eigenvector such that $\bar{\mathbf{l}}_k\mathbf{e} = 0$ . We can therefore choose $\boldsymbol{\alpha} = c_1\mathbf{l}_1 + c_k \mathbf{l}_k + c_k \bar{\mathbf{l}}_k $ such that $\boldsymbol{\alpha}$ is a probability vector, and we can choose $\boldsymbol{\alpha}$ from a restricted span of $(\mathbf{l}_1,\mathbf{l}_k + \bar{\mathbf{l}}_k)$ .

If more left eigenvectors are orthogonal on $\mathbf{e}$ , then the space becomes a span of more vectors. Notice that since the left eigenvectors are linearly independent, the dimension of the span is equal to the number of vectors in it. The most extreme case is $\mathbf{S} \mathbf{e} = -\lambda_1 \mathbf{e}$ . Here $\mathbf{S}$ is not simple, as $\mathbf{l}_{j,\nu} \mathbf{e} = 0$ for $\{j,\nu\} \in \{2, \ldots, g\} \times \{1,\ldots, b_j\}$ . In this case all terms related to all other eigenvalues are removed and every choice of $\boldsymbol{\alpha}$ makes Y exponentially distributed.

For $\mathbf{S}$ reducible the situation is not as clear, because $\mathbf{l}_1$ may have elements of different signs and zero-elements, or it may not make the representation irreducible. Thus $c_1$ and $c_k$ cannot always be chosen so that the linear combination $\boldsymbol{\alpha} = c_1\mathbf{l}_1 + c_k \mathbf{l}_k$ is a non-negative vector with at least one positive element.

3.4. The algebraic degree of phase-type distributions

The Laplace transform of a phase-type distribution can be written as a rational function, and in [Reference O’Cinneide18] the algebraic degree of phase-type distributions was defined from the Laplace transform.

Definition 4. (O’Cinneide [Reference O’Cinneide18].) Let $L(\theta)$ be the Laplace transform of a phase-type distribution with representation $\left(\boldsymbol{\alpha},\mathbf{S}\right)$ . Then we have

\begin{align*} L(\theta) = \boldsymbol{\alpha} \left(\theta \mathbf{I} -\mathbf{S}\right)^{-1} \mathbf{s} = \dfrac{p(\theta)}{q(\theta)},\end{align*}

where $p(\theta)$ and $q(\theta)$ are coprime polynomials. The algebraic degree of the phase-type distribution with representation $\left(\boldsymbol{\alpha},\mathbf{S}\right)$ is then the degree of the polynomial $q(\theta)$ .

The algebraic degree of a phase-type distribution does not always match the number of states, p, in the representation $\left(\boldsymbol{\alpha},\mathbf{S}\right)$ . In [Reference Commault and Chemla8], the following relationship between PH-simplicity and the algebraic degree of a phase-type representation was established.

By [Reference O’Cinneide17, Theorem 1], the condition in Equation (16) implies PH-simplicity. We have already noted the link between PH-simplicity and exponential distributions. It turns out that there is also a link between the rank of the two matrices in Equations (16) and (17) and the exponential distribution.

Proof. Consider the series expansion of the survival function of $Y \sim \text{PH}\left(\boldsymbol{\alpha},\mathbf{S}\right)$ ,

(18) \begin{align}\begin{split} \mathbb{P}(Y_{\text{PH}}>y) &= \boldsymbol{\alpha} \exp(\mathbf{S} y) \mathbf{e} \\ &= \boldsymbol{\alpha}\sum_{i=0}^{\infty} \dfrac{1}{i!} \mathbf{S}^{i} y^{i} \mathbf{e}, \end{split}\end{align}

and the survival function of an exponential distribution with rate $\lambda$ ,

(19) \begin{align}\begin{split} \mathbb{P}(Y_{\exp}>y) &= \exp({-}\lambda y) \\ &= \sum_{i=0}^{\infty} \dfrac{1}{i!}({-} \lambda)^{i} y^{i}. \end{split} \end{align}

By the uniqueness of the Taylor expansion, Y is only exponentially distributed if $({-}\lambda)^{i} = \boldsymbol{\alpha}\mathbf{S}^{i} \mathbf{e}$ . The last statement has to hold for every i, and it therefore gives us an infinite set of equations. We consider the first $p+1$ equations and obtain

\begin{align*} &\boldsymbol{\alpha} \mathbf{S} \mathbf{e} = ({-}\lambda), \\ &\boldsymbol{\alpha} \mathbf{S}^2 \mathbf{e} = ({-}\lambda)^{2}, \\ & \vdots \\ &\boldsymbol{\alpha} \mathbf{S}^{p} \mathbf{e} = ({-}\lambda)^{p}, \\ &\boldsymbol{\alpha} \mathbf{S}^{p+1} \mathbf{e} = ({-}\lambda)^{p+1},\end{align*}

which can be reformulated as

(20) \begin{align} \boldsymbol{\alpha}(\mathbf{S}^{i+1}-({-}\lambda) \mathbf{S}^{i}) \mathbf{e} = 0, \quad i \in \{1,\ldots p\}. \end{align}

To see the first claim we modify the system into

\begin{align*} \boldsymbol{\alpha}(\mathbf{S}-({-}\lambda)\mathbf{I})\mathbf{S}^{i} \mathbf{e} = 0, \quad i \in \{1,\ldots p\}.\end{align*}

Now assume $\text{rank}(\mathbf{A}) = p$ . This means that $\mathbf{S}^{i} \mathbf{e}$ , $i \in \{1,\ldots ,p\}$ , is a basis for $\mathbb{R}^{p}$ . The vector $\boldsymbol{\alpha}(\mathbf{S}-({-}\lambda)\mathbf{I})$ has dimension p, but the only p-dimensional vector which is orthogonal on a basis for $\mathbb{R}^{p}$ is the zero-vector. Thus $\boldsymbol{\alpha}(\mathbf{S}-({-}\lambda)\mathbf{I}) = 0$ . But this equation is only satisfied by an eigenvector corresponding to the eigenvalue $-\lambda$ , and therefore $\boldsymbol{\alpha}$ is a left eigenvector for $\mathbf{S}$ . The argument for the second claim runs in the same way. We modify the system in (20) into

\begin{align*} \boldsymbol{\alpha}\mathbf{S}^{i}(\mathbf{S}-({-}\lambda)\mathbf{I}) \mathbf{e} = 0, \quad i \in \{1,\ldots p\},\end{align*}

and carry out the same arguments on $\mathbf{B}$ .

The two matrices which can be used to check whether the algebraic degree of a phase-type representation is p thus also reveal the ways in which an exponential distribution can arise from the representation. If Equation (16) holds, then a phase-type distribution can only be exponential if $\boldsymbol{\alpha}$ is a left eigenvector, and if Equation (17) holds, then it can only be exponential if $\mathbf{e}$ is a right eigenvector. Consider now the case where Equation (16) and Equation (17) both hold. Then both $\boldsymbol{\alpha}$ and $\mathbf{e}$ have to be eigenvectors. But if Equation (16) and Equation (17) are satisfied, this is only the case for $p = 1$ .

As for the second condition in Theorem 6, if $\mathbf{S}$ is simple, then there is only one way to construct an exponential distribution. By [Reference O’Cinneide17, Theorem 1], PH-simplicity implies $\text{rank}(\mathbf{A}) = p$ , and thus this happens if $\boldsymbol{\alpha}$ is a left eigenvector. If we time-reverse the exponential representation, then by Corollary 1, $\mathbf{e}$ will be a right eigenvector, and for any choice of initial vector the generator $\tilde{\mathbf{S}}$ will govern an exponential distribution.

4.1. Application 1: an exponential distribution outside the eigenvector structure

Perhaps the most surprising conclusion of Theorems 2 and 3 is that one can create a phase-type representation for an exponential distribution without either $\boldsymbol{\alpha}$ or $\mathbf{e}$ being a left or right eigenvector, respectively. Consider a phase-type-distributed random variable Y with representation

\begin{align*} \left(\boldsymbol{\alpha},\mathbf{S} \right) = \left(\begin{pmatrix} 1 & \quad 0 & \quad 0 \end{pmatrix}, \begin{pmatrix} - 3 & \quad 7/8 & \quad 1/8 \\[2pt] 0 & \quad -2 & \quad 2/7 \\[2pt] 0 & \quad 0 & \quad -4 \end{pmatrix} \right).\end{align*}

We have $\mathbb{P}(Y>y) = \boldsymbol{\alpha} \exp\left(\mathbf{S} y \right) \mathbf{e} = \exp({-}2 y)$ , so $Y \sim \exp(2)$ . The eigenvalues are $\lambda_1 = -2$ , $\lambda_2 = -3$ , and $\lambda_3 = -4$ . The left eigenvectors are

\begin{align*} \mathbf{l}_1 = \begin{pmatrix} 0 & \quad 1 & \quad 1/7 \end{pmatrix}, \quad \mathbf{l}_2 = \begin{pmatrix} 1 & \quad -7/8 & \quad -1/8 \end{pmatrix}, \quad \mathbf{l}_3 = \begin{pmatrix} 0 & \quad 0 & \quad 1 \end{pmatrix},\end{align*}

and the right eigenvectors are

\begin{align*} \mathbf{r}^{\prime}_1 = \begin{pmatrix} 7/8 & \quad 1 & \quad 0 \end{pmatrix}, \quad \mathbf{r}^{\prime}_2 = \begin{pmatrix} 1 & \quad 0 & \quad 0 \end{pmatrix}, \quad \mathbf{r}^{\prime}_3 = \begin{pmatrix} 0 & \quad -1/7 & \quad 1 \end{pmatrix}.\end{align*}

Notice that the right and left eigenvectors are scaled so that

\begin{align*} \mathbf{T} = \begin{pmatrix} \mathbf{r}_1 & \quad \mathbf{r}_2 & \quad \mathbf{r}_3 \end{pmatrix},\quad \mathbf{T}^{-1} = \begin{pmatrix} \mathbf{l}_1 \\[3pt] \mathbf{l}_2 \\[3pt] \mathbf{l}_3 \end{pmatrix}.\end{align*}

Clearly, $\boldsymbol{\alpha}$ is not a left eigenvector and $\mathbf{e}$ is not a right eigenvector. However, we have

\begin{align*} \boldsymbol{\alpha} \mathbf{r}_3 &= 0, \\ \mathbf{l}_2 \mathbf{e} &=0,\end{align*}

and, by Theorem 3, Y is therefore exponential. In particular, consider the sub-generator of states 2 and 3,

\begin{align*} \mathbf{S}^{*} = \begin{pmatrix} -2 & \quad 2/7 \\ 0 & \quad -4 \end{pmatrix}.\end{align*}

This matrix has eigenvalues $\lambda_1^{*} = -2$ and $\lambda_2^{*} = -4$ . The corresponding left eigenvectors are $\mathbf{l}_1^* = \begin{pmatrix} 7/8 & \quad 1/8 \end{pmatrix}$ and $\mathbf{l}_2^* = \begin{pmatrix} 0 & \quad 1\end{pmatrix}$ . Now the vector of rates from state 1 into states 2 and 3, $\mathbf{S}_{1;(2,3)}$ , equals $\mathbf{l}_1^{*}$ , so if the Markov jump process enters state 2 or 3, the time in these states is exponential with rate 2. The exit rate from state 1 is also 2, and this gives an intuitive understanding of why the distribution is exponential.

4.2. Application 2: a bivariate exponential distribution

One motivation for studying new generators for exponential phase-type distributions is to use these representations to create multivariate exponential distributions. MPH $^*$ distributions, introduced in [Reference Kulkarni14], are a generalization of phase-type distributions to multiple dimensions, and we can use phase-type representations for exponential distributions to create MPH $^*$ representations for multivariate exponential distributions. MPH $^*$ distributions introduce a $(p \times m)$ -dimensional reward matrix, $\mathbf{R} = \{r_{ij}\}$ , with all entries non-negative to create a multivariate phase-type distribution with m variables. The construction uses one instance of the underlying Markov jump process, $\{X_t\}_{t \geq 0}$ , but scales the time spent in the transient states differently according to the reward matrix. Let $Z_i$ be the time that $\{X_t\}_{t \geq 0}$ spends in state i before absorption:

\begin{align*} Z_i = \int_{0}^Y \mathbf{1}\left(X_t = i\right) dt.\end{align*}

For a univariate phase-type-distributed variable Y, we then have

\begin{align*} Y = \sum_{i=1}^p Z_i,\end{align*}

and for the MPH $^*$ construction we have

\begin{align*} Y_j = \sum_{i=1}^p r_{ij} Z_i.\end{align*}

An appealing property of MPH $^*$ distributions is that the marginal distribution of e.g. $Y_j$ is of phase type. A phase-type representation for the marginal distribution of $Y_j$ can be found through transformation by rewards (see e.g. [Reference Bladt and Nielsen6, Theorem 8.1.6]). Let $E_j^{+}$ and $E_j^{0}$ be the sets of states in the state space $E = \{1,2,\ldots, p\}$ that assign a positive and a zero reward to $Y_j$ , respectively:

\begin{align*} i \in E_j^{+} \quad & \text{if } r_{ij} > 0, \\ i \in E_j^{0} \quad & \text{if } r_{ij} = 0 .\end{align*}

New let $\boldsymbol{\alpha}_j^{+} = \{\alpha_i\}_{i \in E_j^{+}}$ and $\boldsymbol{\alpha}_j^{0} = \{\alpha_i\}_{i \in E_j^{0}}$ be such that $\boldsymbol{\alpha}_j^{+}$ contains the initial probabilities for the states in $E_j^{+}$ , and $\boldsymbol{\alpha}_j^{0}$ contains the initial probabilities for the states in $E_j^{0}$ . Similarly, let $\mathbf{S}_j^{++} = \{s_{ik}\}_{i,k \in E_j^{+}}$ , $\mathbf{S}_j^{+0} = \{s_{ik}\}_{i \in E_j^{+},k \in E_j^{0}}$ , $\mathbf{S}_j^{0+} = \{s_{ik}\}_{i \in E_j^{0},k \in E_j^{+}}$ , and $\mathbf{S}_j^{00} = \{s_{ik}\}_{i,k \in E_j^{0}}$ . A phase-type representation for $Y_j$ , $(\boldsymbol{\alpha}_j,\mathbf{S}_j)$ , is now given by

(21) \begin{align}\begin{split} \boldsymbol{\alpha}_j &= \boldsymbol{\alpha}_j^{+} + \boldsymbol{\alpha}_j^{0} \Big({-}\mathbf{S}_j^{00} \Big)^{-1} \mathbf{S}_j^{0+}, \\ \mathbf{S}_j &= \Delta(\mathbf{r}_j^+)^{-1} \Big(\mathbf{S}_j^{++} + \mathbf{S}_j^{+0} \Big({-}\mathbf{S}_j^{00} \Big)^{-1} \mathbf{S}_j^{0+}\Big), \end{split}\end{align}

where $\mathbf{r}_j^+ = \{r_{ij}\}_{j \in E^+}$ .

Some known bivariate exponential distributions have an MPH $^*$ representation. This is the case for Kibble’s bivariate exponential distribution (introduced in [Reference Kibble13]), Marshall and Olkin’s bivariate exponential distribution (introduced in [Reference Marshall and Olkin15]), Raftery’s bivariate exponential distribution (introduced in [Reference Raftery19, Reference Raftery20]), and the Farlie–Gumbel–Morgenstern distribution with exponential marginals; see [Reference Baker2, Reference Bladt and Nielsen5]. The MPH $^*$ representations for these distributions are given in [Reference Bladt and Nielsen6, Chapter 8], and if the formula in Equation (21) is applied, the marginal distributions are all seen to satisfy one of the conditions in Theorem 1.

In this example we construct an MPH $^*$ representation for a bivariate exponential distribution such that the marginal representations do not satisfy the conditions of Theorem 1, and we explore the impact this has on the correlation structure. A simple exponential phase-type distribution, which does not have $\boldsymbol{\alpha} \mathbf{S} = -\lambda \mathbf{S}$ or $\mathbf{S} \mathbf{e} = - \lambda \mathbf{e}$ , arises if we make two ‘versions’ of one state and make them behave differently. In a multivariate setting, this allows for some interesting models, as one version may be shared across variables and another not. As an example, consider a representation for a 2-state Coxian distribution,

(22) \begin{align}\begin{split} \left(\boldsymbol{\alpha},\mathbf{S} \right) = \left(\begin{pmatrix} 1 & \quad 0 \end{pmatrix}, \begin{pmatrix} - \mu_1 & \quad \mu_1 q \\ 0 & \quad -\mu_2 \end{pmatrix} \right). \end{split} \end{align}

If $\mu_1 (1-q) = \mu_2$ , then the exit rate is constant, and $Y \sim \text{PH}\left(\boldsymbol{\alpha},\mathbf{S}\right)$ is exponential with rate $\mu_2 = \lambda$ . One way to construct a bivariate exponential distribution from this representation is to let the first state be shared. An MPH $^*$ representation is then

\begin{align*} \left( \boldsymbol{\alpha},\mathbf{S},\mathbf{R}\right) = \left(\begin{pmatrix} 1 & \quad 0 & \quad 0 \end{pmatrix}, \begin{pmatrix} - 1 & \quad q & \quad q (1-q) \\[3pt] 0 & \quad -1 & \quad q \\[3pt] 0 & \quad 0 & \quad -1 \end{pmatrix}, \begin{pmatrix} \dfrac{1-q}{\lambda^{(1)}} & \quad \dfrac{1-q}{\lambda^{(2)}} \\[8pt] \dfrac{1}{\lambda^{(1)}} & \quad 0 \\[8pt] 0 & \quad \dfrac{1}{\lambda^{(2)}} \end{pmatrix} \right),\end{align*}

where the rate of $Y_1$ is $\lambda^{(1)}$ and the rate of $Y_2$ is $\lambda^{(2)}$ . This distribution has a singularity, because there is a positive probability $(1-q)^2$ that $\lambda^{(1)} Y_1 = \lambda^{(2)} Y_2$ , and the distribution resembles that of Marshall and Olkin, which also has such a singularity; see [Reference Marshall and Olkin15]. This feature can be problematic when real data are considered. To remedy this, we make two (or really three) versions of the first state. One version is shared, and the other two only count for one of the variables. If the shared state is entered, the process always transits to the two individual states with rates $\lambda^{(1)}$ and $\lambda^{(2)}$ . If the individual states are entered, the process absorbs. This gives the MPH $^*$ representation

(23) \begin{align}\begin{split} \boldsymbol{\alpha} & \quad =\begin{pmatrix} q & \quad (1-q) & \quad 0 & \quad 0 & \quad 0 \end{pmatrix}, \\ \mathbf{S}& \quad = \begin{pmatrix} -1 & \quad 0 & \quad 0 & \quad 1 & \quad 0 \\ 0 & \quad -1 & \quad 1 & \quad 0 & \quad 0 \\ 0 & \quad 0 & \quad -1& \quad 0 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad -1 & \quad 1\\ 0 & \quad 0 & \quad 0 & \quad 0 & \quad -1 \end{pmatrix}, \\ \mathbf{R}& \quad = \begin{pmatrix} \dfrac{1-q}{\lambda^{(1)}} & \quad \dfrac{1-q}{\lambda^{(2)}} \\[12pt] \dfrac{1-q}{\lambda^{(1)}}& \quad 0 \\[8pt] 0 & \quad \dfrac{1-q}{\lambda^{(2)}} \\[8pt] \dfrac{1}{\lambda^{(1)}} & \quad 0 \\[8pt] 0 & \quad \dfrac{1}{\lambda^{(2)}} \end{pmatrix}. \end{split} \end{align}

The two marginal representations, $(\boldsymbol{\alpha}_1,\mathbf{S}_1)$ and $(\boldsymbol{\alpha}_2,\mathbf{S}_2)$ , are found using transformation by rewards (see Equation (21)):

(24) \begin{align}\begin{split} \left( \boldsymbol{\alpha}_1,\mathbf{S}_1\right) & \quad = \left( \begin{pmatrix} q & \quad (1-q) & \quad 0 \end{pmatrix}, \begin{pmatrix} -\dfrac{\lambda^{(1)}}{1-q} & \quad 0 & \quad \dfrac{\lambda^{(1)}}{1-q} \\[5pt] 0 & \quad -\dfrac{\lambda^{(1)}}{1-q} & \quad 0 \\[5pt] 0 & \quad 0 & \quad - \lambda^{(1)} \end{pmatrix}\right), \\ \left( \boldsymbol{\alpha}_2,\mathbf{S}_2\right) & \quad = \left( \begin{pmatrix} q & \quad (1-q) & \quad 0 \end{pmatrix}, \begin{pmatrix} -\dfrac{\lambda^{(2)}}{1-q} & \quad 0 & \quad \dfrac{\lambda^{(2)}}{1-q} \\[5pt] 0 & \quad -\dfrac{\lambda^{(2)}}{1-q} & \quad 0 \\[5pt] 0 & \quad 0 & \quad - \lambda^{(2)} \end{pmatrix}\right). \end{split} \end{align}

Direct calculation of the density function confirms that $Y_1 \sim \exp(\lambda^{(1)})$ and $Y_2 \sim \exp(\lambda^{(2)})$ .

The bivariate exponential distribution bears similarities to Raftery’s bivariate distribution; see [Reference Raftery19,Reference Raftery20]. Let $\begin{pmatrix} Y_1 & \quad Y_2 \end{pmatrix}$ be two random variables which are distributed according to Raftery’s bivariate exponential distribution. Then we can write

\begin{align*} \begin{pmatrix} Y_1 & \quad Y_2 \end{pmatrix} =\,& \mathbf{1}(J=0) \begin{pmatrix} W_1 \dfrac{1-q}{\lambda^{(1)}} & \quad W_2 \dfrac{1-q}{\lambda^{(2)}}\end{pmatrix} + \\ &\mathbf{1}(J=1) \begin{pmatrix} W_1 \dfrac{1-q}{\lambda^{(1)}} + V \dfrac{1}{\lambda^{(1)}} & \quad W_2 \dfrac{1-q}{\lambda^{(2)}} + V \dfrac{1}{\lambda^{(2)}} \end{pmatrix},\end{align*}

and for the representation in Equation (23),

\begin{align*} \begin{pmatrix} Y_1 & \quad Y_2 \end{pmatrix} =\,& \mathbf{1}(J=0) \begin{pmatrix} W_1 \dfrac{1-q}{\lambda^{(1)}} & \quad W_2 \dfrac{1-q}{\lambda^{(2)}}\end{pmatrix} + \\ &\mathbf{1}(J=1) \begin{pmatrix} W \dfrac{1-q}{\lambda^{(1)}} + V_1 \dfrac{1}{\lambda^{(1)}} & \quad W \dfrac{1-q}{\lambda^{(2)}} + V_2 \dfrac{1}{\lambda^{(2)}} \end{pmatrix}.\end{align*}

Both distributions are thus a mixture of a bivariate exponential distribution with independent marginals, and a bivariate, generalized Erlang distribution where one exponential term is shared. The difference lies in the rate of the shared exponential term of the bivariate Erlang distribution. In Raftery’s bivariate exponential distribution, the rate is $\lambda^{(i)}$ . In the distribution discussed here, the rate is $\lambda^{(i)}/(1-q)$ , which is larger; the dependence is therefore not as strong as in Raftery’s bivariate exponential distribution. Figure 1 illustrates the differences in densities for three different values of q. In Raftery’s bivariate exponential distribution the correlation is

\begin{align*} \text{Cor}(Y_1,Y_2)= 2q - q^2,\end{align*}

and for the distribution in Equation (23) the correlation is

\begin{align*} \text{Cor}(Y_1,Y_2) = q^3 - 3q^2 + 2q.\end{align*}

Figure 1.

The top three panels show the contour lines of the joint density of the new bivariate exponential distribution with the representation in Equation (23) for three different values of q. The bottom three panels show the joint density for Raftery’s bivariate exponential distribution for the same values of q. In all cases we take $\lambda^{(1)} = \lambda^{(2)} = 1$ . For $q = 0.2$ the dependence structures are similar, but for $q = 0.5$ and $q = 0.8$ they are different, and the dependence is not as strong in the new model.

The marginal representations in Equation (24) have $\boldsymbol{\alpha}_1 \mathbf{S}_1 \neq - \lambda^{(1)} \boldsymbol{\alpha}_1$ , $\mathbf{S}_1 \mathbf{e} \neq - \lambda^{(1)} \mathbf{e}$ , and similarly for $(\boldsymbol{\alpha}_2,\mathbf{S}_2)$ . The model is therefore an example of a bivariate exponential distribution that falls outside the eigenvector structure. The eigenvalues of $\mathbf{S}_1$ are $-\lambda^{(1)}$ with multiplicity 1 and $-\lambda^{(1)}/(1-q)$ with multiplicity 2. The left eigenvectors are

\begin{align*} \mathbf{l}_1 = \begin{pmatrix} 0 & \quad 0 & \quad 1 \end{pmatrix}, \quad \mathbf{l}_{2} = \begin{pmatrix} 0 & \quad 1 & \quad 0 \end{pmatrix}, \quad \mathbf{l}_{3} = \begin{pmatrix} 1 & \quad 0 & \quad \dfrac{1}{-q}\end{pmatrix}, \\\end{align*}

and the right eigenvectors are

\begin{align*} \mathbf{r}^{\prime}_1 = \begin{pmatrix} \dfrac{1}{q} & \quad 0 & \quad 1\end{pmatrix}, \quad \mathbf{r}^{\prime}_2 = \begin{pmatrix} 0 & \quad 1 & \quad 0 \end{pmatrix}, \quad \mathbf{r}^{\prime}_3 = \begin{pmatrix} 1 & \quad 0 & \quad 0 \end{pmatrix}. \quad\end{align*}

The eigenvectors $\mathbf{l}_1$ and $\mathbf{r}_1$ are associated with $-\lambda^{(1)}$ , while the eigenvectors $\mathbf{l}_2$ , $\mathbf{l}_3$ , $\mathbf{r}_2$ , and $\mathbf{r}_3$ are associated with $-\lambda^{(1)}/(1-q)$ . As $\mathbf{S}_1$ is diagonalizable, the condition in Theorem 3 corresponds to the condition in Equation (12), and we have

\begin{align*} \boldsymbol{\alpha}_1 \left( \mathbf{r}_2 \mathbf{l}_2 + \mathbf{r}_3 \mathbf{l}_3 \right) \mathbf{e} = \begin{pmatrix} q & \quad (1-q) & \quad 0 \end{pmatrix} \begin{pmatrix} 1 & \quad 0 & \quad -\dfrac{1}{q}\\ 0 & \quad 1 & \quad 0 \\ 0 & \quad 0 & \quad 0\end{pmatrix} \mathbf{e} = 0.\end{align*}

Similar calculations can be carried out for $\mathbf{S}_2$ .

4.3. Application 3: when is a Coxian distribution exponential?

Theorems 2 and 3 give the same necessary and sufficient conditions for phase-type distributions to be exponential, namely,

$$\boldsymbol{\alpha} \mathbf{T} \left( \sum_{i \in A_j} \mathbf{E}_i \mathbf{N}_i^{\nu-1} \mathbf{E}_i\right) \mathbf{T}^{-1} \mathbf{e} = 0$$

for $\{j,\nu\} \in \{2,\ldots, p\} \times \{1,\ldots, B_j\}$ . The strength of that result lies in the necessity of the condition, and here Theorem 3 is a significant improvement of Theorem 2, as it characterizes all representations of practical use. In this example we will use Theorem 3 to describe the conditions under which a Coxian distribution is exponential. This case is particularly interesting because it is studied in [Reference Buchholz7,Reference He, Zhang and Vera12], which optimize the parameters in a Coxian representation for an exponential distribution that leads to the broadest correlation band for bivariate exponential distributions. In [Reference He, Zhang and Vera12], however, the Coxian representations are subject to stronger conditions than necessary. With the necessary conditions in place, one could find the truly optimal parameters within the class of Coxian distributions; see [Reference Buchholz7]. It turns out that a Coxian distribution can only be exponential if it has constant exit rate, that is, $-\mathbf{S} \mathbf{e} = \lambda \mathbf{e}$ . First we will present a theorem which states that if a phase-type representation is such that $\mathbf{S}$ has no right eigenvectors orthogonal to $\boldsymbol{\alpha}$ , then the phase-type distribution is exponential if and only if $-\mathbf{S} \mathbf{e} = \lambda \mathbf{e}$ . We will use this to prove our statement about Coxian generators.

Proof. It follows from Theorem 1 that if $\mathbf{S} \mathbf{e} = -\lambda \mathbf{e}$ , then $Y \sim \exp(\lambda)$ . The difficult part is to show that if $\mathbf{S}$ has no right eigenvectors orthogonal to $\boldsymbol{\alpha}$ and Y is exponentially distributed, then $-\mathbf{S}\mathbf{e} = \lambda \mathbf{e}$ . Assume therefore that Y is exponentially distributed with rate $\lambda$ and that $\mathbf{S}$ has no right eigenvectors orthogonal to $\boldsymbol{\alpha}$ . By Theorem 3 a phase-type representation is exponential if and only if $B_1 = 1$ and

\begin{align*} c_{j,\nu} = \boldsymbol{\alpha} \mathbf{T} \left(\sum_{i \in A_j} \mathbf{E}_i \mathbf{N}_i^{\nu-1} \mathbf{E}_i\right) \mathbf{T}^{-1} \mathbf{e} = 0, \end{align*}

for $\{j,\nu\} \in \{2,\ldots, g\} \times \{1,\ldots,B_j\}$ . Let us begin with the case where each Jordan block has a distinct eigenvalue. In that case the set $A_j$ only contains one element, which we will take to be j. We use the same approach as in the proof of Theorem 2 (see [Reference Bean and Green4]) and write out the constants $c_{j,\nu}$ . For $B_j$ we have

\begin{align*} c_{j,B_j} = \left(\boldsymbol{\alpha} \mathbf{r}_{j,1}\right)\left(\mathbf{l}_{j,1} \mathbf{e} \right). \end{align*}

As $\left(\boldsymbol{\alpha} \mathbf{r}_{j,1}\right) \neq 0$ and $c_{j,B_1} = 0$ by assumption, we have $\left(\mathbf{l}_{j,1} \mathbf{e} \right) = 0$ . We can now write $c_{j,B_j-1}$ as

\begin{align*} c_{j,B_j-1} = \left(\boldsymbol{\alpha} \mathbf{r}_{j,1} \right) \left( \mathbf{l}_{j,2} \mathbf{e} \right) + \left(\boldsymbol{\alpha} \mathbf{r}_{j,2} \right) \left( \mathbf{l}_{j,1} \mathbf{e} \right) = \left(\boldsymbol{\alpha} \mathbf{r}_{j,1} \right) \left( \mathbf{l}_{j,2} \mathbf{e} \right),\end{align*}

as $\left(\mathbf{l}_{j,1} \mathbf{e} \right) = 0$ . From this it follows that $\left(\mathbf{l}_{j,2} \mathbf{e} \right) = 0$ . It is easy to continue this procedure and realize that $\mathbf{l}_{j,\nu} \mathbf{e} = 0$ for $\{j,\nu\} \in \{2,\ldots,g\} \times \{1,\ldots, B_j\}$ . The case where there is more than one Jordan block associated with the same eigenvalue is more complicated. Assume for simplicity that there are at most two Jordan blocks associated with the same eigenvalue. This does not change the proof, but makes notation easier. Let $j_1,j_2$ be the indices of the Jordan blocks associated with the eigenvalue $-\lambda_j$ . Let them furthermore be ordered so that $b_{j_1} \geq b_{j_2}$ . Then $B_j = b_{j_1}$ . For $\nu > b_{j_2}$ the constant $c_{j,\nu}$ takes the same form as in the previous calculations, and we conclude $\mathbf{l}_{j_1,\rho} \mathbf{e} = 0$ for $\rho \leq (b_{j_1} - b_{j_2})$ . For $\nu = b_{j_2}$ we have

\begin{align*} c_{j,b_{j_2}} &= \left(\boldsymbol{\alpha} \mathbf{r}_{j_2,1}\right) \left( \mathbf{l}_{j_2,1} \mathbf{e} \right) +\left(\boldsymbol{\alpha} \mathbf{r}_{j_1,1} \right) \left(\mathbf{l}_{j_1,b_{j_1} - b_{j_2} + 1} \mathbf{e}\right) + \sum_{i=1}^{b_{j_2}-1} \left( \boldsymbol{\alpha} \mathbf{r}_{j_1, 1+i}\right) \left( \mathbf{l}_{j_1,b_{j_1}- b_{j_2} + 1 - i} \mathbf{e} \right) \\ &=\left(\boldsymbol{\alpha} \mathbf{r}_{j_2,1}\right) \left( \mathbf{l}_{j_2,1} \mathbf{e} \right) + \left(\boldsymbol{\alpha} \mathbf{r}_{j_1,1} \right) \left(\mathbf{l}_{j_1,b_{j_1}-b_{j_2} + 1} \mathbf{e}\right) \\ &= \boldsymbol{\alpha} \left(\mathbf{r}_{j_2,1} \left(\mathbf{l}_{{j_2},1} \mathbf{e}\right) +\mathbf{r}_{j_1,1} \left(\mathbf{l}_{{j_1},b_{j_1} - b_{j_2} + 1} \mathbf{e}\right) \right).\end{align*}

The second equality is true because $\left( \mathbf{l}_{j_1,\rho} \mathbf{e} \right) = 0$ for $\rho \leq b_{j_1} - b_{j_2}$ . Now $\left(\mathbf{l}_{{j_2},1} \mathbf{e}\right)$ and $\Big(\mathbf{l}_{{j_1},b_{j_1} - b_{j_2} + 1} \mathbf{e}\Big)$ are both scalar numbers. The vectors $\mathbf{r}_{j_2,1}$ and $\mathbf{r}_{j_1,1}$ are both right eigenvectors associated with the eigenvalue $-\lambda_j$ , and any linear combination of these eigenvectors is therefore also a right eigenvector associated with $-\lambda_j$ . Thus

$$\mathbf{v} = \mathbf{r}_{j_2,1}\left(\mathbf{l}_{{j_2},1} \mathbf{e}\right) +\mathbf{r}_{j_1,1} \left(\mathbf{l}_{{j_1},b_{j_1} - b_{j_2} + 1} \mathbf{e}\right) $$

is a right eigenvector, and by definition is not orthogonal to $\boldsymbol{\alpha}$ . As $c_{j,b_{j_2}} = 0$ and as $\mathbf{r}_{j_2,1}$ and $\mathbf{r}_{j_1,1}$ are linearly independent, we have $\mathbf{l}_{{j_2},b_{j_2}} \mathbf{e} = 0$ and $\mathbf{l}_{{j_1}, b_{j_1} - b_{j_2} +1} \mathbf{e} = 0$ . We can now continue the same procedure, as for $\nu = b_{j_2} -1$ we obtain

\begin{align*} c_{j,b_{j_2}-1} =& \left(\boldsymbol{\alpha} \mathbf{r}_{j_2,1}\right) \left( \mathbf{l}_{j_2,2} \mathbf{e} \right) + \left(\boldsymbol{\alpha} \mathbf{r}_{j_1,1} \right) \left(\mathbf{l}_{j_1,b_{j_1} - b_{j_2}+2} \mathbf{e}\right) + \left(\boldsymbol{\alpha} \mathbf{r}_{j_2,2}\right) \left( \mathbf{l}_{j_2,1} \mathbf{e} \right) + \\ &\left(\boldsymbol{\alpha} \mathbf{r}_{j_1,2} \right) \left(\mathbf{l}_{j_1,b_{j_1} - b_{j_2} + 1} \mathbf{e}\right) + \sum_{i=1}^{b_{j_2}} \left( \boldsymbol{\alpha} \mathbf{r}_{j_1, i+ 2}\right) \left( \mathbf{l}_{j_1,b_{j_1}-b_{j_2} + 1- i} \mathbf{e} \right) \\ =& \left(\boldsymbol{\alpha} \mathbf{r}_{j_2,1}\right) \left( \mathbf{l}_{j_2,2} \mathbf{e} \right) + \left(\boldsymbol{\alpha} \mathbf{r}_{j_1,1} \right) \left(\mathbf{l}_{j_1,b_{j_1}- b_{j_2}-2} \mathbf{e}\right),\end{align*}

and we can carry out the same argument to see that

$$\left( \mathbf{l}_{j_2,2} \mathbf{e} \right) = \left(\mathbf{l}_{j_1,b_{j_1}- b_{j_2}+ 2} \mathbf{e}\right) =0.$$

We continue this procedure to establish $\mathbf{l}_{j,\nu} \mathbf{e} = 0$ for $\{j,\nu\} \in \{2,\ldots,g\} \times \{1,\ldots, B_j\}$ . If there are more than two Jordan blocks associated with the same eigenvalue, we can use the same argument; we will just have three or more right eigenvectors in the sums.

Now if $\mathbf{l}_{j,\nu} \mathbf{e} = 0$ for $\{j,\nu\} \in \{2,\ldots,g\} \times \{1,\ldots, B_j\}$ , then $\mathbf{e}$ is a right eigenvector associated with the last eigenvalue, $-\lambda_1$ . To see this, note first that any right eigenvector associated with $-\lambda_1$ is orthogonal to all the vectors in $\{\mathbf{l}_{j,\nu}\}_{\{2,\ldots,g\} \times \{1,\ldots, B_j\}}$ , because $\mathbf{T}$ contains the (generalized) right eigenvectors of $\mathbf{S}$ and $\mathbf{T}^{-1}$ contains the (generalized) left eigenvectors of $\mathbf{S}$ . The vector $\mathbf{e}$ is also orthogonal to the subspace $\{\mathbf{l}_{j,\nu}\}_{\{2,\ldots,g\} \times \{1,\ldots, B_j\}}$ , and as $B_1 = 1$ it must therefore belong to the subspace spanned by all the right eigenvectors associated with $-\lambda_1$ . Thus $\mathbf{e}$ is itself a right eigenvector associated with $-\lambda_1$ . This shows that if $Y \sim \exp(\lambda)$ and there is no right eigenvector orthogonal to $\boldsymbol{\alpha}$ , then $\mathbf{e}$ is a right eigenvector. This finishes the proof.

It is now easy to show our claimed statement about Coxian representations.

Proof. By Theorem 8 it is sufficient to show that a Coxian distribution does not have any right eigenvectors orthogonal to $\boldsymbol{\alpha} = \mathbf{e}_1$ .

Consider therefore a Coxian phase-type representation

\begin{align*} \left(\boldsymbol{\alpha},\mathbf{S}_C \right) = \left(\begin{pmatrix} 1 & \quad 0 & \quad 0 \cdots & \quad 0 & \quad 0 \end{pmatrix}, \begin{pmatrix} -\mu_1 & \quad \mu_1 q_1 & \quad 0 & \quad \cdots & \quad 0 & \quad 0 \\ 0 & \quad -\mu_2 & \quad \mu_2 q_2 & \quad \cdots & \quad 0 & \quad 0 \\ 0 & \quad 0 & \quad -\mu_3 & \quad \cdots & \quad 0 & \quad 0 \\ \vdots & \quad \vdots & \quad \vdots & \quad \ddots & \quad \vdots & \quad \vdots \\ 0 & \quad 0 & \quad 0 & \quad \cdots & \quad -\mu_{p-1} & \quad \mu_{p-1} q_{p-1} \\ 0 & \quad 0 & \quad 0 & \quad \cdots & \quad 0 & \quad -\mu_p \end{pmatrix} \right), \end{align*}

where $q_i > 0$ for all i, as the representation is irreducible. Now assume that there exists a vector $\mathbf{v}$ , different from the zero-vector, that satisfies

1. (1) $\boldsymbol{\alpha} \mathbf{v} = 0$ ,

2. (2) $\mathbf{S}_C \mathbf{v} = -\lambda \mathbf{v}$ .

The vector $\mathbf{v}$ is different from the zero-vector, so it must have at least one non-zero element. Let $v_i$ be the ith element of $\mathbf{v}$ . Assume first that $v_1 \neq 0$ . This clearly cannot be the case, as $\boldsymbol{\alpha} \mathbf{v} = \alpha_1 v_1 = v_1$ , which contradicts (1); therefore $v_1 = 0$ . Assume next that $v_2 \neq 0$ . But by (2) this cannot be, as $0 = \lambda v_1 = \left(\mathbf{S}_C \mathbf{v}\right)_1 = - \mu_1 v_1 + \mu_1 q_1 v_2$ , and $\mu_1 q_1 > 0$ and $v_1 = 0$ . Thus $v_2 = 0$ . Assume now that $v_3 \neq 0$ . But by the same argument as for $v_2$ , this cannot be true, as $ 0=\lambda v_2 = \left(\mathbf{S}_C \mathbf{v}\right)_2 =-\mu_2 v_2 + \mu_2 q_2 v_3$ . Therefore also $v_3 = 0$ . We can continue this argument all the way to $v_p$ , and it is clear that $\mathbf{v} = \mathbf{0}$ . But this is a contradiction, and therefore no vector different from the zero-vector satisfies $\boldsymbol{\alpha} \mathbf{v} = 0$ and $\mathbf{S}_C \mathbf{v} = - \lambda \mathbf{v}$ . By Theorem 8 a Coxian distribution can therefore only be exponential if $\mathbf{S}_C \mathbf{e} = -\lambda \mathbf{e}$ .