2.1 Background on DPH

Consider a Markov chain $M=\{M(n)\}_{n\in\mathbb{N}}$ evolving in a finite state-space $\mathcal{E}=\{1, 2, \dots, d\}$ according to some initial probability vector $\boldsymbol{\alpha} = (\alpha_1, \alpha_2, \dots, \alpha_d)$ and subprobability matrix $\boldsymbol{\boldsymbol{P}}=\{p_{ij}\}_{i,j \in \mathcal{E}}$ . Here, all states in $\mathcal{E}$ are assumed to be transient, and thus, the process M will be (almost surely) terminating. The transience property is mathematically equivalent to the condition that the eigenvalues of $\boldsymbol{\boldsymbol{P}}$ are strictly contained in the complex unit circle (see e.g., Corollary 1.2.63 in Bladt and Nielsen, Reference Bladt and Nielsen2017).

An alternative setup is to consider the augmented state-space $\mathcal{E} \cup \Delta$ , where $\Delta$ is an absorbing cemetery state used to indicate the termination of M. In this setup, M evolves according to the probability vector $(\alpha_1, \alpha_2, \dots, \alpha_d, 0)$ and a block-partitioned transition matrix

\begin{align*}\left(\begin{array}{c@{\quad}c}\boldsymbol{\boldsymbol{P}} & \boldsymbol{1} - \boldsymbol{\boldsymbol{P}}\boldsymbol{1} \\ \boldsymbol{0} & 1\end{array}\right),\end{align*}

where $\boldsymbol{1}$ denotes a column vector of ones of appropriate dimension (d in this case), $\boldsymbol{P}\boldsymbol{1}$ corresponds to the sum of probabilities from each state to all others in $\mathcal{E}$ , and $\boldsymbol{0}$ denotes a zero matrix of appropriate dimension (a zero row vector of dimension d in this case). The final “1” is a scalar representing the probability of remaining in the state $\Delta$ once entered. Both setups – the terminating and the absorbing ones – are interchangeable for our purposes, and we switch back and forth depending on the situation.

Define $\tau$ to be the step in which M gets terminated, that is, let

\begin{align*}\tau \, :\! =\inf\{n\in\mathbb{N}_+\, : \,M(n)\notin\mathcal{E}\}.\end{align*}

Then $\tau$ follows what is known as a DPH distribution, which we denote by $\operatorname{DPH}(\boldsymbol{\alpha},\boldsymbol{P})$ . The class of DPH distributions is particularly tractable (see, e.g., Chapter 1.2.6 in Bladt and Nielsen, Reference Bladt and Nielsen2017). For instance, its pmf over $\mathbb{N}_+$ is explicitly given by

\begin{align*}f_\tau(\ell) & \, :\! = \mathbb{P}(\tau=\ell) = \boldsymbol{\alpha}\boldsymbol{P}^{\ell - 1} \left(\boldsymbol{1}-\boldsymbol{P}\boldsymbol{1}\right).\end{align*}

Likewise, the pgf takes the tractable form

\begin{align*}\mathbb{E}(x^{\tau}) = x \boldsymbol{\alpha}(\boldsymbol{I}-x\boldsymbol{P})^{-1} \left(\boldsymbol{1}-\boldsymbol{P}\boldsymbol{1}\right),\end{align*}

where $\boldsymbol{I}$ denotes the identity matrix. Note that the above formula is valid for all $x\in\mathbb{R}$ on the real line such that the eigenvalues of $x\boldsymbol{P}$ lie within the complex unit circle (the set that contains the interval $[\! -1,1]$ ).

2.2 Definition, joint probability, and joint moments of CDPH

As mentioned in the introduction, literature on the construction of multivariate DPH distributions is rather limited. Our paper is inspired by the approach in Bladt and Yslas (Reference Bladt and Yslas2023b), adding the flexibility that the Markov chains evolve synchronously for some period of time (rather than just starting in the same initial state) before proceeding independently. The exact details behind our novel construction are as follows. Consider two jointly evolving terminating processes, $ M_1 = \{ M_1(n) \}_{n \in \mathbb{N}} $ and $ M_2 = \{ M_2(n) \}_{n \in \mathbb{N}} $ , each taking values in the state space $ \mathcal{E} \cup \mathcal{S} $ , where $ \mathcal{E} \cap \mathcal{S} = \emptyset $ and $ |\mathcal{E}| \wedge |\mathcal{S}| \geq 1 $ . We construct the joint process $ (M_1, M_2) $ such that each marginal process is a Markov chain, with a specific joint dependence structure: we let $ M_1(n) = M_2(n) $ for all n up to the first time they exit the set $ \mathcal{E} $ . After this point, they evolve independently within $ \mathcal{S} $ and eventually terminate, possibly at different times.

More explicitly, we let $ M_1 $ be a terminating Markov chain that evolves according to the block-partitioned initial probability vector $ (\boldsymbol{\alpha}, \boldsymbol{0}) $ and subprobability matrix

(2.1) \begin{equation}\boldsymbol{P}_1 \, :\! =\left(\begin{array}{c@{\quad}c}\boldsymbol{P} & \boldsymbol{U} \\ \boldsymbol{0} & \boldsymbol{Q}_1\end{array}\right),\end{equation}

where $ \boldsymbol{P} $ is a subprobability matrix over $ \mathcal{E} $ (like before), $ \boldsymbol{U} $ denotes the transition probabilities from $ \mathcal{E} $ to $ \mathcal{S} $ , and $ \boldsymbol{Q}_1 $ is a subprobability matrix over $ \mathcal{S} $ . It is further assumed that each row of the submatrix $(\boldsymbol{P}\,\,\,\boldsymbol{U})$ sums to one, so that the Markov chain $M_1$ (and also $M_2$ below) must enter $\mathcal{S}$ before termination. For $ M_2 $ , we define $ M_2(n) = M_1(n) $ for all $ n \leq \tau_{1,2} $ , where

\begin{align*}\tau_{1,2} \, :\! = \inf\{ n \in \mathbb{N}_+\, : \,M_1(n) \notin \mathcal{E} \}.\end{align*}

Obviously, $\tau_{1,2}$ follows a $\operatorname{DPH}(\boldsymbol{\alpha},\boldsymbol{P})$ distribution. After time $ \tau_{1,2} $ , we let $ M_2 $ evolve independently as a terminating Markov chain within $ \mathcal{S} $ according to the subprobability matrix $ \boldsymbol{Q}_2 $ . This construction yields $ M_2 $ as a Markov chain in $ \mathcal{E} \cup \mathcal{S} $ that evolves according to the block-partitioned initial probability vector $ (\boldsymbol{\alpha}, \boldsymbol{0}) $ and subprobability matrix

\begin{align*}\boldsymbol{P}_2 \, :\! =\left(\begin{array}{c@{\quad}c}\boldsymbol{P} & \boldsymbol{U} \\ \boldsymbol{0} & \boldsymbol{Q}_2\end{array}\right).\end{align*}

We now define the bivariate random variable $ (\tau_1, \tau_2) $ via the termination times of $ M_1 $ and $ M_2 $ so that, for $ k \in \{1, 2\} $ ,

\begin{align*}\tau_k \, :\! = \inf\{ n \in \mathbb{N}_+ \, : \, M_k(n) \notin \mathcal{E} \cup \mathcal{S} \}.\end{align*}

For each $ k \in \{1, 2\} $ , the marginal distribution of $\tau_k$ is clearly $\operatorname{DPH}((\boldsymbol{\alpha}, \boldsymbol{0}),\boldsymbol{P}_k)$ . Moreover, $ \tau_1 $ and $ \tau_2 $ are greater than $ \tau_{1,2} $ . Specifically, by writing $\tau_1=\tau_{1,2} + (\tau_1-\tau_{1,2})$ and $\tau_2=\tau_{1,2} + (\tau_2-\tau_{1,2})$ , one observes that $\tau_1$ and $\tau_2$ are dependent for two reasons:

• • Both $\tau_1$ and $\tau_2$ share the term $\tau_{1,2}$ , the time until which the Markov chains $M_1$ and $M_2$ evolve identically. The random variable $\tau_{1,2}$ can be regarded as the common shock component.

• • Although the Markov chains $M_1$ and $M_2$ evolve independently after time $\tau_{1,2}$ , the remaining times until they get absorbed (namely $\tau_1-\tau_{1,2}$ and $\tau_2-\tau_{1,2}$ ) are only conditionally independent given the state $M_1(\tau_{1,2})$ of the Markov chain $M_1$ at the time $\tau_{1,2}$ . Unconditionally, the variables $\tau_1-\tau_{1,2}$ and $\tau_2-\tau_{1,2}$ are dependent.

The joint distribution of $ (\tau_1, \tau_2) $ is provided in the following proposition.

Proposition 2.1. For $ n_1, n_2 \geq 2 $ , $ m, z_1, z_2 \geq 1 $ , and $ i \in \mathcal{E} $ , define the unconditional and conditional joint pmf’s of $ (\tau_1, \tau_2) $ by

\begin{align*}f_{\tau_1,\tau_2}(n_1, n_2) & \, :\! = \mathbb{P}(\tau_1 = n_1, \tau_2 = n_2), \\f_{\tau_1,\tau_2|M_1,M_2}(n_1, n_2|i) & \, :\! = \mathbb{P}(\tau_1 = n_1, \tau_2 = n_2 \mid M_1(0) = M_2(0) = i),\end{align*}

and define the auxiliary functions

\begin{align*}\psi(m, z_1, z_2) & \, :\! = \mathbb{P}(\tau_{1,2} = m, \tau_1 - \tau_{1,2} = z_1, \tau_2 - \tau_{1,2} = z_2), \\ \psi_i(m, z_1, z_2) & \, :\! = \mathbb{P}(\tau_{1,2} = m, \tau_1 - \tau_{1,2} = z_1, \tau_2 - \tau_{1,2} = z_2 \mid M_1(0) = M_2(0) = i).\end{align*}

Then, we have

(2.2) \begin{align}f_{\tau_1,\tau_2}(n_1, n_2) = \sum_{m=1}^{\min\{n_1, n_2\} - 1} \psi(m, n_1 - m, n_2 - m),\end{align}

(2.3) \begin{align}f_{\tau_1,\tau_2|M_1,M_2}(n_1, n_2|i) = \sum_{m=1}^{\min\{n_1, n_2\} - 1} \psi_i(m, n_1 - m, n_2 - m), \end{align}

and

(2.4) \begin{align}\psi(m, z_1, z_2) = \boldsymbol{\alpha} \boldsymbol{P}^{m - 1} \boldsymbol{U} \left( (\boldsymbol{Q}_1^{z_1 - 1} \boldsymbol{q}_1) \odot (\boldsymbol{Q}_2^{z_2 - 1} \boldsymbol{q}_2) \right),\end{align}

(2.5) \begin{align}\psi_i(m, z_1, z_2) = \boldsymbol{e}_i^\intercal \boldsymbol{P}^{m - 1} \boldsymbol{U} \left( (\boldsymbol{Q}_1^{z_1 - 1} \boldsymbol{q}_1) \odot (\boldsymbol{Q}_2^{z_2 - 1} \boldsymbol{q}_2) \right), \end{align}

where $ \boldsymbol{q}_1 = \boldsymbol{1} - \boldsymbol{Q}_1 \boldsymbol{1} $ and $ \boldsymbol{q}_2 = \boldsymbol{1} - \boldsymbol{Q}_2 \boldsymbol{1} $ denote the one-step probabilities of termination from the states in $ \mathcal{S} $ for $ M_1 $ and $ M_2 $ , respectively; $ \boldsymbol{e}_i $ denotes the i -th canonical column vector, and $ \odot $ denotes the Hadamard (entrywise) product of two vectors.

Using Proposition 2.1, one can readily show the following corollary.

Corollary 2.2. For $ \zeta_0 , \zeta_1, \zeta_2\in [0,1]$ , the joint pgf of $(\tau_{1,2}, \tau_1 - \tau_{1,2}, \tau_2 - \tau_{1,2})$ is given by

(2.6) \begin{align} \mathbb{E}&\left(\zeta_0^{\tau_{1,2}} \zeta_1^{\tau_1-\tau_{1,2}} \zeta_2^{\tau_2-\tau_{1,2}} \right) = \sum_{j\in\mathcal{S}}(\boldsymbol{\alpha}(\boldsymbol{I}\zeta_0^{-1} -\boldsymbol{P})^{-1}\boldsymbol{U})_j((\boldsymbol{I}\zeta_1^{-1} -\boldsymbol{Q}_1)^{-1}\boldsymbol{q}_1)_j((\boldsymbol{I}\zeta_2^{-1} -\boldsymbol{Q}_2)^{-1}\boldsymbol{q}_2)_j.\end{align}

Moreover, the joint pgf of $(\tau_1, \tau_2)$ is given by

(2.7) \begin{align} \mathbb{E}\left( \zeta_1^{\tau_1} \zeta_2^{\tau_2} \right) = \sum_{j\in\mathcal{S}}(\boldsymbol{\alpha}(\boldsymbol{I}(\zeta_1\zeta_2)^{-1} -\boldsymbol{P})^{-1}\boldsymbol{U})_j((\boldsymbol{I}\zeta_1^{-1} -\boldsymbol{Q}_1)^{-1}\boldsymbol{q}_1)_j((\boldsymbol{I}\zeta_2^{-1} -\boldsymbol{Q}_2)^{-1}\boldsymbol{q}_2)_j.\end{align}

Employing Corollary 2.2, the following result regarding the corresponding joint factorial moments can be obtained, which can in turn be used to compute the joint moments of $\tau_1$ and $\tau_2$ (see Remark 2.1) Here, for $n\in\mathbb{N}$ , we let $x^{[n]}=x(x-1)(x-2)\cdots (x-n+1)$ denote the falling factorial of x.

Corollary 2.3. For $n_0, n_1,n_2\in\mathbb{N}$ , the joint factorial moments of $(\tau_{1,2},\tau_1-\tau_{1,2},\tau_2-\tau_{1,2})$ are given by

(2.8) \begin{align} \mathbb{E}\bigl(\tau_{1,2}^{[n_0]}\,(\tau_1-\tau_{1,2})^{[n_1]}\,(\tau_2-\tau_{1,2})^{[n_2]}\bigr)= \sum_{j\in\mathcal{S}}(\boldsymbol{\alpha}\widetilde{\boldsymbol{P}}(n_0)\boldsymbol{U})_j(\widetilde{\boldsymbol{Q}}_1(n_1)\boldsymbol{q}_1)_j(\widetilde{\boldsymbol{Q}}_2(n_2)\boldsymbol{q}_2)_j,\end{align}

where

\begin{align*}\widetilde{\boldsymbol{P}}(n_0)& \, :\! =\left. \frac{\partial^{n_0}}{\partial \zeta_1^{n_0}}(\boldsymbol{I}\zeta_0^{-1}-\boldsymbol{P})^{-1}\right|_{\zeta_0=1} = n_0!\boldsymbol{P}^{n_0-1}(\boldsymbol{I}-\boldsymbol{P})^{-n_0-1},\qquad n_0\in\mathbb{N}_+,\\ \widetilde{\boldsymbol{Q}}_1(n_1)& \, :\! =\left. \frac{\partial^{n_1}}{\partial \zeta_1^{n_1}}(\boldsymbol{I}\zeta_1^{-1}-\boldsymbol{Q}_1)^{-1}\right|_{\zeta_1=1} = n_1!\boldsymbol{Q}_1^{n_1-1}(\boldsymbol{I}-\boldsymbol{Q}_1)^{-n_1-1},\qquad n_1\in\mathbb{N}_+,\\ \widetilde{\boldsymbol{Q}}_2(n_2)& \, :\! =\left. \frac{\partial^{n_2}}{\partial \zeta_2^{n_2}}(\boldsymbol{I}\zeta_2^{-1}-\boldsymbol{Q}_2)^{-1}\right|_{\zeta_2=1} = n_2!\boldsymbol{Q}_2^{n_2-1}(\boldsymbol{I}-\boldsymbol{Q}_2)^{-n_2-1}, \qquad n_2\in\mathbb{N}_+,\end{align*}

with $\widetilde{\boldsymbol{P}}(0)=(\boldsymbol{I}-\boldsymbol{P})^{-1}$ , $\widetilde{\boldsymbol{Q}}_1(0)=(\boldsymbol{I}-\boldsymbol{Q}_1)^{-1}$ and $\widetilde{\boldsymbol{Q}}_2(0)=(\boldsymbol{I}-\boldsymbol{Q}_2)^{-1}$ .

Remark 2.1. The joint moments of $\tau_1$ and $\tau_2$ can indeed be obtained from Corollary 2.3 as follows. For $n_1,n_2\in\mathbb{N}$ , we start by writing $\tau_1=\tau_{1,2}+(\tau_1-\tau_{1,2})$ and $\tau_2=\tau_{1,2}+(\tau_2-\tau_{1,2})$ followed by applying binomial expansions of their power to see that

(2.9) \begin{equation}\mathbb{E}\bigl(\tau_{1}^{n_1}\,\tau_{2}^{n_2}\bigr) = \sum_{j_1=0}^{n_1} \binom{n_1}{j_1} \sum_{j_2=0}^{n_2} \binom{n_2} {j_2} \mathbb{E}\bigl(\tau_{1,2}^{j_1+j_2}(\tau_1-\tau_{1,2})^{n_1-j_1} (\tau_2-\tau_{1,2})^{n_2-j_2}\bigr).\end{equation}

Further using the identity

\begin{equation*}x^n = \sum_{k=0}^n \left\{\genfrac{}{}{0pt}{}{n}{k}\right\} x^{[k]}\end{equation*}

for $n\in\mathbb{N}$ , where

\begin{equation*}\left\{\genfrac{}{}{0pt}{}{n}{k}\right\} \, :\! = \frac{1}{k!} \sum_{j=0}^k (\! -1)^{k-j} \binom{k}{j} j^n\end{equation*}

is the Stirling number of the second kind (see, e.g., Chapter 24.1.4 in Abramowitz and Stegun, Reference Abramowitz and Stegun1972), (2.9) can be rewritten as

\begin{align*}\mathbb{E}\bigl(\tau_{1}^{n_1}\,\tau_{2}^{n_2}\bigr) = &\sum_{j_1=0}^{n_1} \binom{n_1}{j_1} \sum_{j_2=0}^{n_2} \binom{n_2}{j_2} \sum_{i=0}^{j_1+j_2}\left\{\genfrac{}{}{0pt}{}{j_1+j_2}{i}\right\} \sum_{\ell_1=0}^{n_1-j_1}\left\{\genfrac{}{}{0pt}{}{n_1-j_1}{\ell_1}\right\} \sum_{\ell_2=0}^{n_2-j_2}\left\{\genfrac{}{}{0pt}{}{n_2-j_2}{\ell_2}\right\} \\&\times\mathbb{E}\bigl(\tau_{1,2}^{[i]}(\tau_1-\tau_{1,2})^{[\ell_1]} (\tau_2-\tau_{1,2})^{[\ell_2]}\bigr).\end{align*}

That is, the joint moments of $\tau_1$ and $\tau_2$ can be represented as a linear combination of the joint factorial moments of $(\tau_{1,2},\tau_1-\tau_{1,2},\tau_2-\tau_{1,2})$ available in (2.8).

Note that the support of $(\tau_1, \tau_2)$ is $(\mathbb{N} + 2) \times (\mathbb{N} + 2)$ . For modeling purposes, it is often useful to consider discrete random variables that have a lattice different from $\mathbb{N} + 2$ . Thus, for $k_1, k_2 \in \mathbb{R}$ and $c_1, c_2 \gt 0$ , one can define the vector $(N_1, N_2)$ where

(2.10) \begin{align}N_1 \, :\! = c_1 (\tau_1-2) + k_1, \qquad N_2 \, :\! = c_2 (\tau_2 -2) + k_2.\end{align}

We then say in a broad sense (see Remark 2.3) that the bivariate random variable $(N_1, N_2)$ follows a bivariate CDPH distribution with support on $(c_1 \mathbb{N} + k_1) \times (c_2 \mathbb{N} + k_2)$ and parameters $(\boldsymbol{\alpha}, \boldsymbol{P}, \boldsymbol{U}, \boldsymbol{Q}_1, \boldsymbol{Q}_2)$ . Employing Proposition 2.1, we obtain the following immediate result.

Corollary 2.4. Let $(N_1, N_2)$ follow a bivariate $\operatorname{CDPH}$ distribution with support $(c_1 \mathbb{N} + k_1) \times (c_2 \mathbb{N} + k_2)$ and parameters $(\boldsymbol{\alpha}, \boldsymbol{P}, \boldsymbol{U}, \boldsymbol{Q}_1, \boldsymbol{Q}_2)$ . Then, the joint pmf of $(N_1, N_2)$ for $x_1 \in c_1 \mathbb{N} + k_1$ and $x_2 \in c_2 \mathbb{N} + k_2$ , namely $f_{N_1, N_2}(x_1, x_2) \, :\! =\mathbb{P}(N_1 = x_1,\, N_2 = x_2)$ , is given by

\begin{align*}f_{N_1, N_2}(x_1, x_2) = f_{\tau_1,\tau_2}\left( \frac{x_1 - k_1}{c_1}+ 2,\; \frac{x_2 - k_2}{c_2}+ 2 \right),\end{align*}

where $f_{\tau_1,\tau_2}(n_1, n_2) = \mathbb{P}(\tau_1 = n_1,\, \tau_2 = n_2)$ .

2.3 Relationships with other bivariate DPH distributions

Let us first note that the (bivariate version of) the so-called $\operatorname{mDPH}$ class considered in Bladt and Yslas (Reference Bladt and Yslas2023b) is a subclass of our proposed $\operatorname{CDPH}$ . Indeed, if we let $\boldsymbol{\beta}$ be a probability row vector, then choosing the parameters $\boldsymbol{\alpha} = (1)$ , $\boldsymbol{P} = (0)$ and $\boldsymbol{U} = \boldsymbol{\beta}$ for a $\operatorname{CDPH}$ over the lattice $(\mathbb{N} + 1)\times(\mathbb{N} + 1)$ (i.e., $c_1=c_2=1$ and $k_1=k_2=-1$ ) yields an $\operatorname{mDPH}$ distribution. This distribution arises from starting two terminating processes in the same state according to the initial probability vector $\boldsymbol{\beta}$ then and evolving them independently from there onward according to the subprobability matrices $\boldsymbol{Q}_1$ and $\boldsymbol{Q}_2$ , respectively.

In turn, the $\operatorname{CDPH}$ with support in $(\mathbb{N}+2) \times (\mathbb{N}+2)$ forms a subclass of the $\operatorname{MDPH}^*$ class introduced in Chapter 5.2 of Navarro (Reference Navarro2019). These distributions accumulate integer-valued rewards across each occupation epoch within a shared Markov chain, rather than bifurcating into independent processes after a common shock component concludes. Specifically, let us consider a partitioned state space $\mathcal{E} \cup (\mathcal{S} \times \mathcal{S}) \cup (\{\partial_1\} \times \mathcal{S}) \cup (\mathcal{S} \times \{\partial_2\})$ and a Markov chain $\gamma = \{\gamma(n)\}_{n \in \mathbb{N}}$ defined over it. The chain evolves within $\mathcal{E}$ while the common shock is ongoing, then jumps to $\mathcal{S} \times \mathcal{S}$ (here we take the lexicographic order of $\mathcal{S} \times \mathcal{S}$ ) while the independent components are active, and finally either terminates, jumps to $\{\partial_1\} \times \mathcal{S}$ , or jumps to $\mathcal{S} \times \{\partial_2\}$ . These transitions occur depending on whether the independent components terminate at the same time, the first coordinate terminates before the second, or the second coordinate terminates before the first. Thus, we define the bivariate random vector $(\tau_1,\tau_2)$ as

(2.11) \begin{align}\tau_1 \, :\! =\sum_{\ell=0}^{\sigma-1}\unicode{x1D7D9}_{\gamma(\ell)\in \mathcal{E}\cup(\mathcal{S}\times\mathcal{S})\cup(\mathcal{S}\times\{\partial_2\})},\qquad\tau_2 \, :\! =\sum_{\ell=0}^{\sigma-1}\unicode{x1D7D9}_{\gamma(\ell)\in \mathcal{E}\cup(\mathcal{S}\times\mathcal{S})\cup(\{\partial_1\} \times \mathcal{S})},\end{align}

where $\sigma$ corresponds to the termination time of $\gamma$ so that $\sigma \, :\! =\inf\{ n \in \mathbb{N}_+ : \gamma(n) \notin \mathcal{E} \cup (\mathcal{S} \times \mathcal{S}) \cup (\{\partial_1\} \times \mathcal{S}) \cup (\mathcal{S} \times \{\partial_2\}) \}$ . The above representation of the bivariate vector $(\tau_1,\tau_2)$ is already in the form of Equation (5.5) of Navarro (Reference Navarro2019). It remains to specify the transition probability of the terminating Markov chain $\gamma$ .

To mimic the distributional properties of the original processes $M_1$ and $M_2$ embedded in the Markov chain $\gamma$ , we take the block-partitioned parameters $(\boldsymbol{\alpha},\boldsymbol{0},\boldsymbol{0},\boldsymbol{0})$ and the subprobability matrix

(2.12) \begin{align}\boldsymbol{P}_{\mathrm{max}} \, :\! = \left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}\boldsymbol{P} & \boldsymbol{U}^* & \boldsymbol{0} & \boldsymbol{0}\\ \boldsymbol{0} & \boldsymbol{Q}_1 \otimes \boldsymbol{Q}_2 & \boldsymbol{q}_1 \otimes \boldsymbol{Q}_2 & \boldsymbol{Q}_1 \otimes \boldsymbol{q}_2 \\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{Q}_2 & \boldsymbol{0}\\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{Q}_1\end{array}\right),\end{align}

where $\boldsymbol{U}^* = \{u_{i,(j_1,j_2)}^*\}_{i \in \mathcal{E}, (j_1,j_2) \in \mathcal{S} \times \mathcal{S}}$ is defined by $u_{i,(j_1,j_2)}^* \, :\! = u_{ij_1} \delta_{j_1,j_2}$ , with $\delta_{j_1,j_2}$ denoting the Kronecker delta and $u_{ij}$ the (i,j)-th element of $\boldsymbol{U}$ , and $\otimes$ represents the Kronecker product between matrices. Then, $\sigma$ is $\operatorname{DPH}((\boldsymbol{\alpha},\boldsymbol{0},\boldsymbol{0},\boldsymbol{0}),\boldsymbol{P}_{\mathrm{max}})$ . Note that the resulting parameters are similar to those used to compute the distribution of the maximum between two discrete phase-type distributions in Theorem 1.2.67 of Bladt and Nielsen (Reference Bladt and Nielsen2017). The derivation of these parameters for the process $\gamma$ is as follows:

• • The process evolves in $\mathcal{E}$ according to $\boldsymbol{P}$ , mimicking the common shock component $\tau_{1,2}$ , eventually jumping from a state $i \in \mathcal{E}$ to a state $(j,j) \in \mathcal{S} \times \mathcal{S}$ with probability $u_{ij}$ .

• • It then evolves in $\mathcal{S} \times \mathcal{S}$ according to $\boldsymbol{Q}_1 \otimes \boldsymbol{Q}_2$ such that a transition from $(i_1, i_2) \in \mathcal{S} \times \mathcal{S}$ to $(j_1, j_2) \in \mathcal{S} \times \mathcal{S}$ occurs with probability $q_{1,i_1j_1}\,q_{2,i_2j_2}$ (where, for $k=1,2$ , we use $q_{k,ij}$ to denote the (i,j)-th element of $\boldsymbol{Q}_k$ ). Then, either both components terminate simultaneously according to the exit column vector $\boldsymbol{q}_1 \otimes \boldsymbol{q}_2$ ; the first terminates while the second remains active, moving to $\{\partial_1\} \times \mathcal{S}$ ; or the second terminates while the first remains active, moving to $\mathcal{S} \times \{\partial_2\}$ .

• • If the process jumps to $\{\partial_1\} \times \mathcal{S}$ or $\mathcal{S} \times \{\partial_2\}$ , it continues evolving according to the matrix $\boldsymbol{Q}_2$ or $\boldsymbol{Q}_1$ , respectively.

The above results are summarized in the following proposition.

Proposition 2.5. Concerning bivariate $\operatorname{DPH}$ distributions, our proposed $\operatorname{CDPH}$ class satisfies the relationships

\begin{align*}&{\text{$\operatorname{mDPH}$ $\subset$ $\operatorname{CDPH}$ $\operatorname{(with}$ $\operatorname{support}$ $(\mathbb{N}+1) \times (\mathbb{N}+1))$}},\\&{\text{$\operatorname{CDPH}$ $\operatorname{(with}$ $\operatorname{support}$ $(\mathbb{N}+2) \times (\mathbb{N}+2))$ $\subset$ $\operatorname{MDPH}^*$}},\end{align*}

where $\operatorname{mDPH}$ is the class introduced by Bladt and Yslas Reference Bladt and Yslas(2023b) and $\operatorname{MDPH}^*$ is the class introduced by Navarro Reference Navarro(2019).

It is important to emphasize the practical advantages of our CDPH class relative to the more general MDPH^* framework from Navarro (Reference Navarro2019). While CDPH is a subclass of MDPH^*, its tailored construction for common shocks yields an explicit and closed-form joint pmf (see Proposition 2.1). This analytical tractability, which is not a feature of the general MDPH^* class, is the key that enables direct likelihood-based inference and allows for the data-driven estimation of the underlying common shock distribution, as implemented via the EM algorithm in Section 4.

2.4 Further distributional and closure properties

Below, we explore further properties of the $\operatorname{CDPH}$ class. In particular, we investigate certain operations on the vector $(\tau_1, \tau_2)$ that belong to the class of $\operatorname{DPH}$ or $\operatorname{CDPH}$ distributions. Apart from theoretical significance, such closure properties can also be of practical importance. For example, closure of the CDPH class under mixtures and sums means that the resulting bivariate random vectors still follow $\operatorname{CDPH}$ distributions but with updated parameters, and the closed-form joint pmf of the $\operatorname{CDPH}$ class is still conveniently applicable for computation purposes.

Minimum and maximum between $\tau_1$ and $\tau_2$

Consider again the Markov chain $\gamma$ driven by (2.12), with initial vector $(\boldsymbol{\alpha}, \boldsymbol{0}, \boldsymbol{0}, \boldsymbol{0})$ , state-space $\mathcal{E} \cup (\mathcal{S} \times \mathcal{S}) \cup (\{\partial_1\} \times \mathcal{S}) \cup (\mathcal{S} \times \{\partial_2\})$ , and termination time $\sigma$ . One can easily see that $\sigma$ indeed corresponds to the maximum $\tau_{\max} \, :\! =\max\{\tau_1, \tau_2\}$ of $\tau_1$ and $\tau_2$ . Moreover, the minimum $\tau_{\min} \, :\! =\min\{\tau_1, \tau_2\}$ is simply the exit time of the Markov chain $\gamma$ from the set $\mathcal{E} \cup (\mathcal{S} \times \mathcal{S})$ . Defining the matrix

(2.13) \begin{align}\boldsymbol{P}_{\mathrm{min}} \, :\! = \left(\begin{array}{c@{\quad}c}\boldsymbol{P} & \boldsymbol{U}^* \\ \boldsymbol{0} & \boldsymbol{Q}_1 \otimes \boldsymbol{Q}_2 \\ \end{array}\right),\\[-12pt]\nonumber\end{align}

the next proposition follows as a direct consequence.

Proposition 2.6. Given a $\operatorname{CDPH}(\boldsymbol{\alpha}, \boldsymbol{P}, \boldsymbol{U}, \boldsymbol{Q}_1, \boldsymbol{Q}_2)$ random vector $(\tau_1,\tau_2)$ , the maximum $\tau_{\max}$ follows a $\operatorname{DPH}((\boldsymbol{\alpha},\boldsymbol{0},\boldsymbol{0},\boldsymbol{0}),\boldsymbol{P}_{\mathrm{max}})$ distribution and the minimum $\tau_{\min}$ follows a $\operatorname{DPH}((\boldsymbol{\alpha},\boldsymbol{0}),\boldsymbol{P}_{\mathrm{min}})$ distribution, where $\boldsymbol{P}_{\mathrm{max}}$ and $\boldsymbol{P}_{\mathrm{min}}$ are given by (2.12) and (2.13), respectively.

Aggregation of coordinates $\tau_1 + \tau_2$

When considering the sum of the coordinates of $(\tau_1,\tau_2)$ , it is helpful to represent it as $\tau_1 + \tau_2 = 2\min\{\tau_1, \tau_2\} + (\max\{\tau_1, \tau_2\} - \min\{\tau_1, \tau_2\})$ . That is, $\tau_1 + \tau_2$ counts twice the steps in which both $\tau_1$ and $\tau_2$ have not terminated and then adds the steps during which the larger value exceeds the smaller one. Consequently, we have

(2.14) \begin{align}\tau_1+\tau_2 = \sum_{\ell=0}^{\sigma-1}\left(2 \cdot \unicode{x1D7D9}_{\gamma(\ell) \in \mathcal{E} \cup (\mathcal{S} \times \mathcal{S})} + \unicode{x1D7D9}_{\gamma(\ell) \in (\{\partial_1\} \times \mathcal{S}) \cup (\mathcal{S} \times \{\partial_2\})}\right).\end{align}

We have the following proposition.

Proposition 2.7. Given a $\operatorname{CDPH}(\boldsymbol{\alpha}, \boldsymbol{P}, \boldsymbol{U}, \boldsymbol{Q}_1, \boldsymbol{Q}_2)$ random vector $(\tau_1,\tau_2)$ , the sum $\tau_1+\tau_2$ follows a $\operatorname{DPH}$ distribution.

Mixture of CDPH random vectors

For each $i\in \{1, 2, \dots, m\}$ , suppose that the random vector $(\tau_1(i),\tau_2(i))$ is CDPH distributed. Let a mixture be formed where the component $(\tau_1(i),\tau_2(i))$ is chosen with probability $w_i \gt 0$ . The following proposition shows that the class of CDPH distributions is closed under finite mixture.

Proposition 2.8. Suppose that a mixture is formed such that a bivariate random vector $(\tau_{1,\mathrm{mix}},\tau_{2,\mathrm{mix}})$ is chosen to be $(\tau_1(i),\tau_2(i))$ with probability $w_i \gt 0$ , such that $(\tau_1(i),\tau_2(i))$ is $\operatorname{CDPH}(\boldsymbol{\alpha}(i), \boldsymbol{P}(i), \boldsymbol{U}(i), \boldsymbol{Q}_1(i), \boldsymbol{Q}_2(i))$ distributed and $\sum_{i=1}^m w_i = 1$ . Then, $(\tau_{1,\mathrm{mix}},\tau_{2,\mathrm{mix}})$ follows a $\operatorname{CDPH}(\boldsymbol{\alpha}_{\mathrm{mix}}, \boldsymbol{P}_{\mathrm{mix}}, \boldsymbol{U}_{\mathrm{mix}}, \boldsymbol{Q}_{1,\mathrm{mix}}, \boldsymbol{Q}_{2,\mathrm{mix}})$ distribution, where the parameters are given by

\begin{align*}\boldsymbol{\alpha}_{\mathrm{mix}} \, :\! = \begin{pmatrix}w_1 \boldsymbol{\alpha}(1) & w_2 \boldsymbol{\alpha}(2) & \cdots & w_m \boldsymbol{\alpha}(m)\end{pmatrix},\\[-24pt]\nonumber\end{align*}

\begin{align*}\boldsymbol{P}_{\mathrm{mix}} \, :\! = \left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}\boldsymbol{P}(1)& \boldsymbol{0} & \cdots & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{P}(2) & \cdots & \boldsymbol{0} \\ \vdots & \vdots & \ddots & \vdots \\ \boldsymbol{0} & \boldsymbol{0} & \cdots & \boldsymbol{P}(m)\end{array}\right),\quad\boldsymbol{U}_{\mathrm{mix}} \, :\! = \left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}\boldsymbol{U}(1) & \boldsymbol{0} & \cdots & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{U}(2) & \cdots & \boldsymbol{0} \\ \vdots & \vdots & \ddots & \vdots \\ \boldsymbol{0} & \boldsymbol{0} & \cdots & \boldsymbol{U}(m)\end{array}\right),\end{align*}

\begin{align*}\boldsymbol{Q}_{1,\mathrm{mix}} \, :\! = \left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}\boldsymbol{Q}_1(1) & \boldsymbol{0} & \cdots & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{Q}_1(2) & \cdots & \boldsymbol{0} \\ \vdots & \vdots & \ddots & \vdots \\ \boldsymbol{0} & \boldsymbol{0} & \cdots & \boldsymbol{Q}_1(m)\end{array}\right),\quad\boldsymbol{Q}_{2,\mathrm{mix}} \, :\! = \left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}\boldsymbol{Q}_2(1) & \boldsymbol{0} & \cdots & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{Q}_2(2) & \cdots & \boldsymbol{0} \\ \vdots & \vdots & \ddots & \vdots \\ \boldsymbol{0} & \boldsymbol{0} & \cdots & \boldsymbol{Q}_2(m)\end{array}\right).\end{align*}

Sum of CDPH random vectors

Suppose that the random vector $(\tau_1(i),\tau_2(i))$ is CDPH distributed for each $i\in \{1, 2, \dots, m\}$ , and the random vectors $\{(\tau_1(i),\tau_2(i))\}_{i=1}^m$ are independent. Then one could be interested in the joint distribution of the sum $(\tau_{1,\mathrm{sum}},\tau_{2,\mathrm{sum}}) \, :\! =(\sum_{i=1}^m\tau_1(i),\sum_{i=1}^m\tau_2(i))$ . Indeed, this follows a CDPH distribution. In the next proposition, we consider the case $m=2$ and explicitly state the resulting parameters.

Proposition 2.9. Consider the bivariate random vector $(\tau_{1,\mathrm{sum}},\tau_{2,\mathrm{sum}})=(\tau_1(1)+\tau_1(2),\tau_2(1)+\tau_2(2))$ , where $(\tau_1(i),\tau_2(i))$ is $\operatorname{CDPH}(\boldsymbol{\alpha}(i), \boldsymbol{P}(i), \boldsymbol{U}(i), \boldsymbol{Q}_1(i), \boldsymbol{Q}_2(i))$ distributed for $i\in\{1,2\}$ , and the random vectors $(\tau_1(1),\tau_2(1))$ and $(\tau_1(2),\tau_2(2))$ are independent. Then, $(\tau_{1,\mathrm{sum}},\tau_{2,\mathrm{sum}})$ follows a $\operatorname{CDPH}(\boldsymbol{\alpha}_{\mathrm{sum}}, \boldsymbol{P}_{\mathrm{sum}}, \boldsymbol{U}_{\mathrm{sum}}, \boldsymbol{Q}_{1,\mathrm{sum}}, \boldsymbol{Q}_{2,\mathrm{sum}})$ distribution, where the parameters are given by

(2.15) \begin{align}\boldsymbol{\alpha}_{\mathrm{sum}} \, :\! = (\boldsymbol{\alpha}(1),\boldsymbol{0}),\\[-24pt]\nonumber\end{align}

(2.16) \begin{align}\boldsymbol{P}_{\mathrm{sum}} \, :\! = \left(\begin{array}{c@{\quad}c}\boldsymbol{P}(1) & \boldsymbol{U}(1) \otimes \boldsymbol{\alpha}(2) \\ \boldsymbol{0} & \boldsymbol{I} \otimes \boldsymbol{P}(2)\end{array}\right),\quad\boldsymbol{U}_{\mathrm{sum}} \, :\! = \left(\begin{array}{c@{\quad}c}\boldsymbol{0} & \boldsymbol{0} \\ \boldsymbol{I} \otimes \boldsymbol{U}(2) & \boldsymbol{0}\end{array}\right),\end{align}

(2.17) \begin{align}\boldsymbol{Q}_{1,\mathrm{sum}} \, :\! = \left(\begin{array}{c@{\quad}c}\boldsymbol{Q}_1(1) \otimes \boldsymbol{I} & \boldsymbol{q}_1(1) \otimes \boldsymbol{I} \\ \boldsymbol{0} & \boldsymbol{Q}_1(2)\end{array}\right),\quad\boldsymbol{Q}_{2,\mathrm{sum}} \, :\! = \left(\begin{array}{c@{\quad}c}\boldsymbol{I} \otimes \boldsymbol{Q}_2(1) & \boldsymbol{I} \otimes \boldsymbol{q}_2(1) \\ \boldsymbol{0} & \boldsymbol{Q}_2(2)\end{array}\right),\end{align}

with $ \boldsymbol{q}_{1}(1) = \boldsymbol{1} - \boldsymbol{Q}_{1}(1) \boldsymbol{1} $ and $ \boldsymbol{q}_{2}(1) = \boldsymbol{1} - \boldsymbol{Q}_{2}(1) \boldsymbol{1} $ . The state space corresponding to the transitions specified by $\boldsymbol{P}_{\mathrm{sum}}$ is $\mathcal{E}_{\mathrm{sum}} \, :\! = \mathcal{E}(1) \cup (\mathcal{S}(1) \times \mathcal{E}(2))$ , and that for the transitions specified by $\boldsymbol{Q}_{1,\mathrm{sum}}$ and $\boldsymbol{Q}_{2,\mathrm{sum}}$ is ${\mathcal{S}_{\mathrm{sum}} \, :\! =} (\mathcal{S}(1) \times \mathcal{S}(2)) \cup \mathcal{S}(2)$ . Here, with the obvious notation, $\mathcal{E}(i)\cup\mathcal{S}(i)$ is the state space of the underlying terminating Markov chains that define the vector $(\tau_1(i),\tau_2(i))$ .

It is clear that $(\tau_{1,\mathrm{sum}},\tau_{2,\mathrm{sum}})=(\sum_{i=1}^m\tau_1(i),\sum_{i=1}^m\tau_2(i))$ for values of m greater than two still follows a CDPH distribution by summing the CDPH bivariate random vectors one at a time and using Proposition 2.9 repeatedly.

3.1 Independent summands

Suppose that $ X_1 $ is independent of $ X_2 $ . Below, we provide information about the distribution of $ (Y_1, Y_2) $ under this assumption.

Proposition 3.1. For $\theta_1,\theta_2\ge0$ , the joint Laplace transform $\mathcal{L}_{Y_1,Y_2}(\theta_1,\theta_2) \, :\! =\mathbb{E}\left( \exp(-\theta_1 Y_1 - \theta_2 Y_2) \right)$ of the bivariate random sum in (3.1) is given by

\begin{align*}\mathcal{L}_{Y_1,Y_2}(\theta_1,\theta_2) = \boldsymbol{\alpha} \left( \boldsymbol{I} a^{-1}_{1,\theta_1} a^{-1}_{2,\theta_2} - \boldsymbol{P} \right)^{-1} \boldsymbol{U} \left( \left( \left( \boldsymbol{I} a^{-1}_{1,\theta_1} - \boldsymbol{Q}_1 \right)^{-1} \boldsymbol{q}_1 \right) \odot \left( \left( \boldsymbol{I} a^{-1}_{2,\theta_2} - \boldsymbol{Q}_2 \right)^{-1} \boldsymbol{q}_2 \right) \right),\end{align*}

where $ a_{k,\theta_k} \, :\! = \mathbb{E}\left( e^{-\theta_k X_{k}} \right) $ is the Laplace transform of $X_{k}$ for $ k \in \{1, 2\} $ .

Proposition 3.1 provides a characterization of the distribution of $ (Y_1, Y_2) $ in terms of the joint Laplace transform while leaving the distributions of $X_1$ and $X_2$ unspecified. It is easy to see that if $a_{1,\theta_1}$ and $a_{2,\theta_2}$ are rational in $\theta_1$ and $\theta_2$ , respectively, then $\mathcal{L}_{Y_1,Y_2}(\theta_1,\theta_2)$ is rational in $ \theta_1 $ and $ \theta_2 $ . Since univariate rational Laplace transform characterizes matrix-exponential distributions (see, e.g., Chapter 4 in Bladt and Nielsen, Reference Bladt and Nielsen2017), while multivariate rational Laplace transform characterizes the class of MME distributions proposed by Bladt and Nielsen (Reference Bladt and Nielsen2010), we have the following corollary as an immediate consequence.

In the following subsection, we establish a related result by showing that if $ (X_1, X_2) $ belongs to the subclass of multivariate phase-type distributions ( $ \operatorname{MPH}^* $ ) proposed in Kulkarni (Reference Kulkarni1989) then so does $ (Y_1, Y_2) $ . The main advantage of $ \operatorname{MPH}^* $ distributions is that, unlike the MME class, they possess a simple probabilistic interpretation.

3.2 Dependent summands

Here we assume that $ (X_1, X_2) $ follows a bivariate version of Kulkarni’s $ \operatorname{MPH}^* $ class, which we briefly describe next. Such an $ \operatorname{MPH}^* $ distribution is constructed using a single terminating continuous-time Markov jump process $J=\{J(t)\}_{t\ge 0}$ with finite state space $\mathcal{C}$ and a system of state-dependent rewards. Specifically, let $ (\boldsymbol{\pi}, \boldsymbol{S}, \boldsymbol{R}) $ be the parameters, where

• • $ \boldsymbol{\pi} $ is the initial probability vector;

• • $ \boldsymbol{S} $ is the subintensity matrix governing the continuous-time transitions between transient states; and

• • $ \boldsymbol{R} =\{r_{kj}\}$ is the reward matrix, where $ r_{kj} $ represents the instantaneous nonnegative reward associated to the j-th state in $\mathcal{C}$ and k-th coordinate for $k\in\{1,2\}$ .

With these elements, J evolves according to $(\boldsymbol{\pi},\boldsymbol{S})$ and the coordinates $X_1$ and $X_2$ are defined via

\begin{align*}X_k \, :\! =\int_0^\delta r_{k,J(s)} \,d s,\qquad k\in\{1,2\},\end{align*}

where $\delta \, :\! =\inf\{t\ge 0\, : \,J(t)\notin\mathcal{C}\}$ is the termination time of J. The random vector $ (X_1, X_2) $ obtained through this construction is said to follow an $ \operatorname{MPH}^*(\boldsymbol{\pi}, \boldsymbol{S}, \boldsymbol{R})$ distribution.

This construction allows for capturing dependencies between $ X_1 $ and $ X_2 $ through the shared underlying Markov process J. Apart from including independence as a special case, the class of $ \operatorname{MPH}^* $ distributions is known to be dense within the class of bivariate distributions with support on the positive quadrant. Unfortunately, the $ \operatorname{MPH}^* $ class is not particularly tractable in the general case as the joint distribution function cannot be expressed in closed form. However, its joint Laplace transform is known to be

\begin{equation*}\mathbb{E} \left( e^{-\theta_1X_1-\theta_2X_2} \right) = \boldsymbol{\pi}\left( \boldsymbol{\Delta}_{\theta_1,\theta_2,\boldsymbol{R}} - \boldsymbol{S} \right)^{-1} (\! -\boldsymbol{S}\boldsymbol{1}),\end{equation*}

for $\theta_1,\theta_2\ge0$ , where

\begin{equation*}\boldsymbol{\Delta}_{\theta_1,\theta_2,\boldsymbol{R}}=\operatorname{diag}\{\theta_1r_{1j}+\theta_2r_{2j}: j\in\mathcal{C}\}.\end{equation*}

Our aim is to prove that, under this setup, the bivariate random vector $(Y_1,Y_2)$ defined in (3.1) with each $(X_{1,\ell},X_{2,\ell})$ distributed as the generic bivariate random vector $(X_1,X_2)$ follows an $\operatorname{MPH}^*$ distribution as well.

\begin{align*}{\boldsymbol{\pi}_{\mathrm{sum}} \, :\! =}\, (\boldsymbol{\alpha}, \boldsymbol{0}, \boldsymbol{0}, \boldsymbol{0}) \otimes \boldsymbol{\pi},\end{align*}

subintensity matrix

\begin{align*}{\boldsymbol{S}_{\mathrm{sum}} \, :\! =}\left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}\boldsymbol{I} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{I} & \boldsymbol{0} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{I} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{I}\end{array}\right) \otimes \boldsymbol{S} + \left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}\boldsymbol{P} & \boldsymbol{U}^* & \boldsymbol{0} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{Q}_1 \otimes \boldsymbol{Q}_2 & \boldsymbol{q}_1 \otimes \boldsymbol{Q}_2 & \boldsymbol{Q}_1 \otimes \boldsymbol{q}_2 \\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{Q}_2 & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{Q}_1\end{array}\right) \otimes (-\boldsymbol{S}\boldsymbol{1}\boldsymbol{\pi}),\end{align*}

and reward matrix

Here, for $k\in\{1,2\}$ , the quantity $\boldsymbol{r}_k$ represents the row vector $(r_{k1},r_{k2},\ldots,r_{k|\mathcal{C}|})$ such that

\begin{align*}\boldsymbol{R} = \begin{pmatrix} \boldsymbol{r}_1 \\ \boldsymbol{r}_2 \end{pmatrix}.\end{align*}

With the random vector $(Y_1,Y_2)$ fully characterized as an $\operatorname{MPH}^*$ distribution, the joint moments of the pair follow directly from Theorem 8.1.5 in Bladt and Nielsen (Reference Bladt and Nielsen2017). For completeness, we explicitly state the formulas for the first- and second-order moments below, which are essential for calculating dependence measures like covariance and correlation.

Corollary 3.4. Suppose that each bivariate random vector in the sequence $\{(X_{1,\ell}, X_{2,\ell})\}_{\ell \in \mathbb{N}_+}$ follows an $\operatorname{MPH}^*$ distribution with parameters $ (\boldsymbol{\pi}, \boldsymbol{S}, \boldsymbol{R}) $ , so that $(Y_1, Y_2)$ defined in (3.1) is $\operatorname{MPH}^*$ distributed with parameters $(\boldsymbol{\pi}_{\mathrm{sum}}, \boldsymbol{S}_{\mathrm{sum}}, \boldsymbol{R}_{\mathrm{sum}})$ given in Proposition 3.3. For $k\in \{1,2\}$ , let $\boldsymbol{D}_k = \mathrm{diag}(\boldsymbol{r}_{k, \mathrm{sum}})$ be the diagonal matrix with the elements of $\boldsymbol{r}_{k, \mathrm{sum}}$ in its diagonal, where $\boldsymbol{r}_{k, \mathrm{sum}}$ is the k-th row of the reward matrix $\boldsymbol{R}_{\mathrm{sum}}$ . Then, for each $k\in \{1,2\}$ , the first moment of the compound sum is given by

$$ \mathbb{E}(Y_k) = \boldsymbol{\pi}_{\mathrm{sum}} (\! -\boldsymbol{S}_{\mathrm{sum}}^{-1})\boldsymbol{r}_{k, \mathrm{sum}}^\intercal,$$

and the second moment is

$$ \mathbb{E}(Y_k^2) = 2 \boldsymbol{\pi}_{\mathrm{sum}} (\! -\boldsymbol{S}_{\mathrm{sum}}^{-1})\boldsymbol{D}_k (\! -\boldsymbol{S}_{\mathrm{sum}}^{-1})\boldsymbol{r}_{k, \mathrm{sum}}^\intercal. $$

Furthermore, the second-order cross-moment is given by

$$ \mathbb{E}(Y_1 Y_2) = \boldsymbol{\pi}_{\mathrm{sum}} (\! -\boldsymbol{S}_{\mathrm{sum}}^{-1})\boldsymbol{D}_1(\! -\boldsymbol{S}_{\mathrm{sum}}^{-1})\boldsymbol{r}_{2, \mathrm{sum}}^\intercal + \boldsymbol{\pi}_{\mathrm{sum}}(\! -\boldsymbol{S}_{\mathrm{sum}}^{-1})\boldsymbol{D}_2(\! -\boldsymbol{S}_{\mathrm{sum}}^{-1})\boldsymbol{r}_{1, \mathrm{sum}}^\intercal. $$

4.1 The fully observed case

Let the parameters be collected in $\boldsymbol{\Theta}: = (\boldsymbol{\alpha} , \boldsymbol{P},\boldsymbol{U},\boldsymbol{Q}_1,\boldsymbol{Q}_2)$ . The likelihood is given by

\begin{align*}\mathcal{L}^O(\boldsymbol{\Theta} \mid \{(n_1^{(m)},n_2^{(m)})\}_{m=1}^n)&=\prod_{m=1}^nf_{\tau_1,\tau_2}(n_1^{(m)},n_2^{(m)};\boldsymbol{\Theta}),\end{align*}

where we place emphasis on the dependence of the pmf on $\boldsymbol{\Theta}$ by writing $f_{\tau_1,\tau_2}(\cdot,\cdot;\boldsymbol{\Theta})$ instead of $f_{\tau_1,\tau_2}(\cdot,\cdot)$ (see Proposition 2.1).

While direct optimization via gradient methods is possible and relatively quick for very small dimensions, the EM algorithm is more precise and efficient for moderate matrix dimensions $|\mathcal{E}|$ and $|\mathcal{S}|$ . Thus, we introduce the following latent variables:

1. 1. $A_i$ : the number of times the processes $\{M_1^{(m)}\}_{m=1}^n$ start in state $i\in\mathcal{E}$ .

2. 2. $N^{A}_{ij}$ : the total number of jumps of $\{M_1^{(m)}\}_{m=1}^n$ from state $i\in\mathcal{E}$ to $j\in\mathcal{E}$ .

3. 3. $N^{T}_{ij}$ : the total number of jumps of $\{M_1^{(m)}\}_{m=1}^n$ from state $i\in\mathcal{E}$ to $j\in\mathcal{S}$ .

4. 4. $N^{B,k}_{ij}$ : the total number of jumps of $\{M_k^{(m)}\}_{m=1}^n$ from state $i\in\mathcal{S}$ to $j\in\mathcal{S}$ , where $k\in\{1,2\}$ .

5. 5. $N^{B,k}_{i}$ : the total number of jumps of $\{M_k^{(m)}\}_{m=1}^n$ from state $i\in\mathcal{S}$ to the absorbing state, where $k\in\{1,2\}$ .

Having specified the above sufficient statistics, we are able to reconstruct the paths of each of the Markov jump processes. Recall that $\alpha_i$ is the i-th element of the row vector $\boldsymbol{\alpha}$ for $i\in \mathcal{E}$ as well as the notation $\boldsymbol{P}=\{p_{ij}\}_{i,j \in \mathcal{E}}$ , $\boldsymbol{Q}_k=\{q_{k,ij}\}_{i,j \in \mathcal{S}}$ and $\boldsymbol{U}=\{u_{ij}\}_{i\in \mathcal{E},j\in \mathcal{S}}$ . For $k\in\{1,2\}$ , we further define $q_{k,i}$ to be the i-th element of the column vector $\boldsymbol{q}_k$ for $i\in \mathcal{S}$ . Then, we can write the complete likelihood as

\begin{align*}&\mathcal{L}^C( \boldsymbol{\Theta}\mid \{M_1^{(m)}\}_{m=1}^n, \{M_2^{(m)}\}_{m=1}^n )\\&=\left(\prod_{i\in\mathcal{E}} \alpha_i^{A_i}\right) \left(\prod_{i,j\in\mathcal{E}} p_{ij}^{N^{A}_{ij}}\right)\left(\prod_{i\in\mathcal{E}}\prod_{j\in\mathcal{S}} u_{ij}^{N^{T}_{ij}}\right)\left(\prod_{k=1}^2\prod_{i,j\in\mathcal{S}} q_{k,ij}^{N^{B,k}_{ij}}\right)\left(\prod_{k=1}^2\prod_{i\in\mathcal{S}} q_{k,i}^{N^{B,k}_{i}}\right).\end{align*}

The above specification is particularly convenient since it belongs to the exponential dispersion family of distributions, and thus has fully explicit maximum likelihood estimators (with a ‘ $\hat{\, \, }$ ’ placed on top of the corresponding parameter) given by

\begin{align*}\hat \alpha_{i}&=\frac{A_i}{n},\quad\hat p_{ij}=\frac{N_{ij}^A}{\sum_{\ell\in\mathcal{E}}N_{i\ell}^{A}+\sum_{\ell\in\mathcal{S}}N_{i\ell}^{T}},\quad\hat u_{ij}=\frac{N_{ij}^T}{\sum_{\ell\in\mathcal{E}}N_{i\ell}^{A}+\sum_{\ell\in\mathcal{S}}N_{i\ell}^{T}},\\ \hat q_{k,ij}&=\frac{N_{ij}^{B,k}}{\sum_{\ell\in\mathcal{S}}N_{i\ell}^{B,k} + N_{i}^{B,k}},\quad\hat q_{k,i}=\frac{N_{i}^{B,k}}{\sum_{\ell\in\mathcal{S}}N_{i\ell}^{B,k} + N_{i}^{B,k}}.\end{align*}

4.2 The EM algorithm

Since the full-trajectory data are not observed, we employ the EM algorithm to find the maximum likelihood estimate iteratively. This implies that at each iteration, the conditional expectations of the sufficient statistics, given the absorption times, are computed, corresponding to the E-step. Subsequently, such conditional expectations are replaced into $\mathcal{L}^C$ instead of the statistics themselves, and so upon maximization, we obtain updated parameters $( \boldsymbol{\Theta} )$ , commonly referred to as the M-step.

The conditional expectations required in the E-step are given as follows:

Lemma 4.1. Let $\xi \, :\! =\{(N^{(m)}_1,N^{(m)}_2)=(n^{(m)}_1,n^{(m)}_2)\}_{m=1}^n$ be the observed data. The conditional expectations required in the E-step are given as follows:

\begin{align*} \mathbb{E}[A_i\mid \xi]=\sum_{m=1}^n \frac{\alpha_i{\kern1.5pt} f_{\tau_1,\tau_2|M_1,M_2}(n_1^{(m)},n_2^{(m)}|i)}{ f_{\tau_1,\tau_2}(n_1^{(m)},n_2^{(m)}) }\,. \end{align*}

\begin{align*} \mathbb{E}[N_{ij}^A\mid \xi]=\sum_{m=1}^n \sum_{t\ge1} \frac{(\boldsymbol{\alpha} \boldsymbol{P}^{t-1})_i\,p_{ij}{\kern1.5pt} f_{\tau_1,\tau_2|M_1,M_2}(n_1^{(m)}-t,n_2^{(m)}-t|j) }{ f_{\tau_1,\tau_2}(n_1^{(m)},n_2^{(m)}) }\,. \end{align*}

\begin{align*}\mathbb{E}[N_{ij}^{T}\mid \xi]=\sum_{m=1}^n \sum_{t\ge1} \frac{ (\boldsymbol{\alpha} \boldsymbol{P}^{t-1})_i \, u_{ij} (\boldsymbol{Q}_1^{n_1^{(m)}-t-1}\boldsymbol{q}_1)_j(\boldsymbol{Q}_2^{n_2^{(m)}-t-1}\boldsymbol{q}_2)_j }{ f_{\tau_1,\tau_2}(n_1^{(m)},n_2^{(m)}) }\,.\end{align*}

\begin{align*}\mathbb{E}[N_{ij}^{B,1}\mid \xi]=\sum_{m=1}^n \sum_{t\ge2}\sum_{\ell=1}^{t-1}\sum_{r\in\mathcal{S}} \frac{ (\boldsymbol{\alpha} \boldsymbol{P}^{\ell-1}\boldsymbol{U})_r (\boldsymbol{Q}_1^{t-\ell-1})_{ri} \, q_{1,ij} (\boldsymbol{Q}_1^{n_1^{(m)}-t-1}\boldsymbol{q}_1)_j (\boldsymbol{Q}_2^{n_2^{(m)}-\ell-1}\boldsymbol{q}_2)_r }{ f_{\tau_1,\tau_2}(n_1^{(m)},n_2^{(m)}) }\,.\end{align*}

\begin{align*}\mathbb{E}[N_{ij}^{B,2}\mid \xi]=\sum_{m=1}^n \sum_{t\ge2}\sum_{\ell=1}^{t-1}\sum_{r\in\mathcal{S}} \frac{ (\boldsymbol{\alpha} \boldsymbol{P}^{\ell-1}\boldsymbol{U})_r (\boldsymbol{Q}_2^{t-\ell-1})_{ri} \, q_{2,ij} (\boldsymbol{Q}_2^{n_2^{(m)}-t-1}\boldsymbol{q}_2)_j (\boldsymbol{Q}_1^{n_1^{(m)}-\ell-1}\boldsymbol{q}_1)_r }{ f_{\tau_1,\tau_2}(n_1^{(m)},n_2^{(m)}) }\,.\end{align*}

\begin{align*}\mathbb{E}[N_{i}^{B,1}\mid \xi]=\sum_{m=1}^n \sum_{t\ge2}\sum_{\ell=1}^{t-1}\sum_{r\in\mathcal{S}} \frac{ (\boldsymbol{\alpha} \boldsymbol{P}^{\ell-1}\boldsymbol{U})_r (\boldsymbol{Q}_1^{t-\ell-1})_{ri} \, q_{1,i} (\boldsymbol{Q}_2^{n_2^{(m)}-\ell-1}\boldsymbol{q}_2)_r }{ f_{\tau_1,\tau_2}(n_1^{(m)},n_2^{(m)}) }\,.\end{align*}

\begin{align*}\mathbb{E}[N_{i}^{B,2}\mid \xi]=\sum_{m=1}^n \sum_{t\ge2}\sum_{\ell=1}^{t-1}\sum_{r\in\mathcal{S}} \frac{ (\boldsymbol{\alpha} \boldsymbol{P}^{\ell-1}\boldsymbol{U})_r (\boldsymbol{Q}_2^{t-\ell-1})_{ri} \, q_{2,i} (\boldsymbol{Q}_1^{n_1^{(m)}-\ell-1}\boldsymbol{q}_1)_r }{ f_{\tau_1,\tau_2}(n_1^{(m)},n_2^{(m)}) }\,.\end{align*}

Here $f_{\tau_1,\tau_2}$ and $f_{\tau_1,\tau_2|M_1,M_2}$ are given by (2.2) and (2.3) with the help of (2.4) and (2.5).

5.1 Bivariate Poisson distribution with common shocks

We begin by demonstrating that our proposed CDPH model can effectively fit a bivariate Poisson distribution. Although the Poisson distribution itself is not a special case of the DPH distribution, we aim to show that the CDPH model provides a good approximation.

To define the bivariate Poisson distribution, we let $(N_1^*,N_2^*)$ be a bivariate Poisson distributed pair such that $N_1^* \, :\! = V_1 + Z$ and $N_2^* \, :\! = V_2 + Z$ , where $V_1$ , $V_2$ , and Z are independent Poisson random variables with parameters $\lambda_{V_1}$ , $\lambda_{V_2}$ , and $\lambda_Z$ , respectively. Under this construction, the marginal distributions are again Poisson with intensities $\lambda_{N_1^*} \, :\! = \lambda_{V_1} + \lambda_Z$ and $\lambda_{N_2^*} \, :\! = \lambda_{V_2} + \lambda_Z$ . The joint pmf of $(N_1^*,N_2^*)$ is provided in Equation (2) of Holgate (Reference Holgate1964).

To evaluate the effectiveness of our model, we simulate data from $(N_1^*,N_2^*)$ and then fit $(N_1^*+2,N_2^*+2)$ using our CDPH approach.Footnote ¹ , The EM algorithm guarantees convergence to local minima, so trying different starting points is generally recommended. An interesting aspect of the study is to observe how the proposed model performs as the value of the common shock rate $\lambda_Z$ increases, while keeping $\lambda_{N_1^*}$ and $\lambda_{N_2^*}$ fixed. Specifically, we fix $\lambda_{N_1^*} = 5$ and $\lambda_{N_2^*} = 4$ and consider three different values of $\lambda_Z$ , namely 1, 2, and 3. The results for the case with $\lambda_Z=1$ are presented in Figure 1, while the other two cases are collected in the Appendix. For each value of $\lambda_Z$ , we explore three different model configurations with dimensions $(|\mathcal{E}|, |\mathcal{S}|) = (2,1)$ , $(|\mathcal{E}|, |\mathcal{S}|) = (3,2)$ , and $(|\mathcal{E}|, |\mathcal{S}|) = (4,3)$ . Since matrix representations are not unique and hence not identifiable, the use of Akaike information criterion (AIC) and Bayesian information criterion (BIC) is not appropriate in this context. Instead, we track the evolution of the likelihood as a function of the EM iteration step to assess model performance.

Figure 1.

Bivariate Poisson distribution with common shock intensity $\lambda_Z=1$ . In the top panels, we show the empirical pmf (black circles) versus the fitted pmf (by increasing matrix dimension: red triangles, blue crosses, and green diamonds), for the distribution of the common shock (left), and the two marginals (center and right). In the middle panel, we have two heatmaps, corresponding to the best-fitted bivariate pmf (left), and the absolute difference between the best fitted and empirical bivariate pmf (right). The bottom panel shows the evolution of the likelihood for the three models as the number of EM steps increases to 500.

The simulations suggest that as the common shock component becomes more pronounced, our model captures it more effectively. That is, the fitted pmf of the common shock component retrieved by our CDPH model offers a more accurate representation of the empirical pmf of the actual common shock Z. Nevertheless, in all cases, the model provides a satisfactory fit, which is quantified through heatmap plots, showing the best-fitted pmf among the three considered models (which is the most complex one with dimensions $(|\mathcal{E}|, |\mathcal{S}|) = (4,3)$ ), as well as the absolute difference between the empirical pmf and the best-fitted one.

5.2 Bivariate Poisson–Lindley distribution

In addition to the bivariate Poisson case, we now investigate the ability of our model to capture overdispersion in claim counts, where the variance is not necessarily equal to the mean. The Poisson–Lindley is a mixture distribution that serves as a natural candidate for this setting. Specifically, the bivariate Poisson–Lindley distribution is constructed by considering conditionally independent Poisson random variables, where the Poisson parameter is mixed over a Lindley distribution that depends on a parameter $\theta$ (see Definition 2 in Gómez-Déniz et al., Reference Gómez-Déniz, Sarabia and Balakrishnan2012). The joint pmf in the bivariate case is also provided in Equation (17) of Gómez-Déniz et al. (Reference Gómez-Déniz, Sarabia and Balakrishnan2012).

To examine the performance of our model, we simulate $(N_1^*, N_2^*)$ from this distribution and apply our CDPH model to the transformed data $(N_1^*+2, N_2^*+2)$ . Concretely, we simulate a Lindley random variable L with $\theta=2$ and subsequently construct $N_1^*$ and $N_2^*$ as independent Poisson with parameters 2L and 3L, respectively. This is then repeated 10,000 times. We again consider three model configurations with dimensions $(|\mathcal{E}|, |\mathcal{S}|) = (2,1)$ , $(|\mathcal{E}|, |\mathcal{S}|) = (3,2)$ , and $(|\mathcal{E}|, |\mathcal{S}|) = (4,3)$ . The evolution of the likelihood across EM iterations is tracked, as in the previous study, to assess convergence and fitting accuracy. The results are provided in Figure 2, and they further illustrate the adaptability of our CDPH approach in modeling overdispersed count data. Unlike the bivariate Poisson setting, the Poisson–Lindley mixture does not explicitly possess a common shock component. Remarkably, however, our model is still able to effectively capture and mimic dependencies between the two components, demonstrating its practical modeling flexibility. In particular, our methodology is able to detect the presence of the inherent (or hidden) common shock component and recover its distribution, which can be useful for risk management purposes.

Figure 2.

Poisson–Lindley distribution with parameter $\theta=2$ . In the top panels, we show the empirical pmf (black circles) versus the fitted pmf (by increasing matrix dimension: red triangles, blue crosses, and green diamonds) for the two marginals (center and right), and the fitted values for a suggested common shock (left). In the middle panel, we have two heatmaps, corresponding to the best-fitted bivariate pmf (left), and the absolute difference between the best fitted and empirical bivariate pmf (right). The bottom panel shows the evolution of the likelihood for the three models as the number of EM steps increases to 500.

5.3 Application to insurance data

We now apply our model to real-world insurance data, previously analyzed in Vernic (Reference Vernic2000). The dataset can be found in Table 1 of Vernic (Reference Vernic2000) and was originally modeled using the bivariate generalized Poisson distribution (BGPD). It consists of claim frequencies from two related insurance categories, with a total of 708 claims.

We fit this dataset to our proposed CDPH distribution, applying the same shift transformation as before. The plots in Figure 3 illustrate the results for the same three different model configurations: $(|\mathcal{E}|, |\mathcal{S}|) = (2,1)$ , $(|\mathcal{E}|, |\mathcal{S}|) = (3,2)$ , and $(|\mathcal{E}|, |\mathcal{S}|) = (4,3)$ . As with the previous cases, we track the likelihood evolution over EM iterations and evaluate the fit using probability plots. The CDPH model in general provides a good fit to the insurance dataset, even in the case of low dimensions. It is interesting to note that a non-negligible common shock component is estimated quite robustly by the three models, which by the previous simulation studies suggests that there could either be a true underlying common shock or that the bivariate distribution is well approximated by the CDPH class. Despite the structural differences between BGPD and CDPH, our approach is also able to capture dependencies effectively.

Figure 3.

Bivariate claim frequency insurance data. In the top panels, we show the fitted pmf (by increasing matrix dimension: red triangles, blue crosses, and green diamonds), for the distribution of the suggested common shock (left), and the two marginals (center and right). In the middle panel, we have two heatmaps, corresponding to the best-fitted bivariate pmf (left), and the absolute difference between the best fitted and empirical bivariate pmf (right). The bottom panel shows the evolution of the likelihood for the three models as the number of EM steps increases to 500.

In fact, we also compare our proposed models directly to the BGPD models of Vernic (Reference Vernic2000), using the $\chi^2$ measure of goodness-of-fit as in that reference. The benchmark BGPD models yield $\chi^2$ values of $25.94$ and $23.61$ , respectively. Among our three matrix-based models, the first two, while simpler, produce higher $\chi^2$ values of $58.13$ and $29.08$ , indicating poorer fit. However, the third and most complex matrix model outperforms both BGPD benchmarks, with a lower $\chi^2$ value of $18.93$ . We expect that models with even larger matrix dimensions (not shown here) could improve fit further, though this comes at the cost of increased parameter complexity. Addressing this trade-off – ideally through predictive performance rather than classical criteria like AIC or BIC, which are known to be inadequate for matrix distributions – is left for future research.