# An extended class of multivariate counting processes and its main properties

Ji Hwan Cha, Sophie Mercier

<figure><img src="image_1.png" /></figure>

## To cite this version:

Ji Hwan Cha, Sophie Mercier. An extended class of multivariate counting processes and its main properties. Probability in the Engineering and Informational Sciences, 2025, 39 (1), pp.83-100. 10.1017/S0269964824000111. hal-04787950

HAL Id: hal-04787950

https://hal.science/hal-04787950v1

Submitted on 18 Nov 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

RESEARCH ARTICLE

# An extended class of multivariate counting processes and its main properties

Ji Hwan Cha<fnref n="1" /> and Sophie Mercier<fnref n="2" />

**Corresponding author:** Ji Hwan Cha; Email: jhcha@ewha.ac.kr

**Keywords:** Multivariate counting process; Multivariate generalized Pólya process; Positive dependence; Restarting property; Stochastic properties

**MSC:** Primary 60K10; Secondary 62P30

## Abstract

In this paper, a new multivariate counting process model (called Multivariate Poisson Generalized Gamma Process) is developed and its main properties are studied. Some basic stochastic properties of the number of events in the new multivariate counting process are initially derived. It is shown that this new multivariate counting process model includes the multivariate generalized Pólya process as a special case. The dependence structure of the multivariate counting process model is discussed. Some results on multivariate stochastic comparisons are also obtained.

## 1. Introduction

Until now, a variety of univariate counting processes for modeling univariate random recurrent events have been developed and studied intensively in the literature. However, in practice, stochastically dependent series of events are also frequently observed, leading to multivariate counting processes. For instance, in queueing models, bivariate point processes arise as the input and output processes (Daley [11]). In reliability applications, the occurrence of recurrent failures in two or more parts in a system are frequently positively dependent. In finance area, a bankruptcy of a financial company in one group may also affect those in other groups (Allen and Gale [2]). In econometrics, multivariate point processes are frequently used to model multivariate market events (see Bowsher [6]). In insurance, two types of recurrent claims can be modeled by a bivariate point process (Partrat [21]). For more plenty of examples, see Cox and Lewis [10]. Although some multivariate counting process models have been developed in the literature, practically available models which can meet practical needs are very limited and there still exists a big gap between the need for proper models in various applications and available useful models.

The main contribution of this paper is to develop a new general class of multivariate counting processes which has mathematical tractability and applicability. Specifically, in the multivariate counting process model developed in this paper, for example, the distribution for the numbers of events in a time interval and the stochastic intensity of the process can be obtained explicitly. This is practically important because it allows explicit expression of the likelihood function in estimation procedure, which increases the utility of the model considerably. Furthermore, the developed counting process model has sufficient flexibility because the baseline intensity functions contained in the model have general forms. In addition, the developed model is very general in the sense that it includes an existing model as a special case.

The paper is organized as follows. In Section 2, a new class of bivariate counting processes is defined and its basic properties are derived. Furthermore, the corresponding marginal processes and the future

<footnote n="1">¹Department of Statistics, Ewha Womans University, Seoul, Republic of Korea</footnote>

<footnote n="2">²Universite de Pau et des Pays de l'Adour, E2S UPPA, CNRS, LMAP, Pau, France</footnote>
process of the developed model are also characterized. In Section 3, the stochastic intensity functions of the model are derived and it is shown that the model includes the bivariate generalized Pólya process in Cha and Giorgio [8] as a special case. In Section 4, the bivariate process is generalized to the multi-variate case, and the results from the previous section are extended. In Section 5, multivariate stochastic comparisons for the numbers of events and the arrival times of the events are studied. In addition, the dependence structure of the process is analyzed.

## 2. Bivariate Poisson generalized gamma process and its basic properties

In this section, the bivariate Poisson generalized gamma process is defined and its basic properties are derived. To define the new counting process model, we first introduce the generalized gamma distribution proposed in Agarwal and Kalla [1] and Ghitany [15]. A random variable $\Phi$ is said to follow the generalized gamma distribution (GGD) with parameters $(\nu, k, \alpha, l)$, where $\nu \ge 0, k, \alpha, l > 0$, if its probability density function (pdf) is given by

$$f(\phi) = \frac{\alpha^{k-\nu} \phi^{k-1} \exp\{-\alpha\phi\}}{\Gamma_{\nu}(k, \alpha l) (\phi + l)^{\nu}}, \quad \phi > 0, \qquad (2.1)$$

where

$$\Gamma_{\nu}(k, \beta) = \int_{0}^{\infty} \frac{y^{k-1} \exp\{-y\}}{(y + \beta)^{\nu}} dy,$$

for all $\beta > 0$, with

$$\Gamma_{\nu}(k, \alpha l) = \int_{0}^{\infty} \frac{y^{k-1} \exp\{-y\}}{(y + \alpha l)^{\nu}} dy = \int_{0}^{\infty} \frac{\alpha^{k-\nu} y^{k-1} \exp\{-\alpha y\}}{(y + l)^{\nu}} dy. \qquad (2.2)$$

The function $\Gamma_{\nu}(k, \beta)$ is called the generalized gamma function (see Kobayashi [20]) and if $\nu = 0$, then

$$\Gamma_{0}(k, \beta) = \int_{0}^{\infty} y^{k-1} \exp\{-y\} dy = \Gamma(k), \quad \forall k > 0.$$

Thus, when $\nu = 0$, it can be seen that the pdf in (2.1) becomes that of a gamma distribution with parameter $(k, \alpha)$. Hence, the GGD includes the gamma distribution as a special case.

Coming back to the general case, one can note from Gupta and Ong [17] that

$$\Gamma_{\nu}(k, \beta) = \frac{\Gamma(k)}{\beta^{\nu-k}} \varphi(k, k - \nu + 1; \beta),$$

where

$$\varphi(a, c; x) = \frac{1}{\Gamma(a)} \int_{0}^{\infty} \frac{e^{-xt} t^{a-1}}{(1+t)^{a-c+1}} dt,$$

is the confluent hypergeometric function of the second kind. This allows an easy computation of $\Gamma_{\nu}(k, \beta)$, as $\varphi(a, c; x)$ and $\Gamma(k)$ are implemented in most statistical or mathematical software. From Ghitany [15], the moment generating function of $\Phi$ is given by
$$M_{\Phi}(s) = E(e^{s\Phi}) = \left(1 - \frac{s}{\alpha}\right)^{\nu-k} \frac{\Gamma_{\nu}(k, (\alpha - s)l)}{\Gamma_{\nu}(k, \alpha l)}, s < \alpha, \quad (2.3)$$

and its r-th moment about the origin is

$$E(\Phi^r) = \alpha^{-r} \frac{\Gamma_{\nu}(k+r, \alpha l)}{\Gamma_{\nu}(k, \alpha l)}, \quad r \in \mathbb{N}^*. \quad (2.4)$$

Now, we will define the bivariate Poisson generalized gamma process using the GGD. Let {$\mathbf{N}(t)$, $t \ge 0$}, where $\mathbf{N}(t) = (N_1(t), N_2(t))$, be a bivariate process. Then, the corresponding “pooled” point process {$M(t)$, $t \ge 0$}, where $M(t) = N_1(t) + N_2(t)$, can be defined. In this paper, we will consider regular (also known as *orderly*) multivariate point processes. In a univariate point process {$N(t)$, $t \ge 0$}, regularity is intuitively the nonoccurrence of multiple events in a small interval (see e.g., Cox and Lewis [10], Finkelstein [13, 14]; see also Cha and Giorgio [8]). Note that there are two types of regularity in multivariate point processes: (i) marginal regularity and (ii) regularity. For a multivariate point process, we say the process is *marginally regular* if its marginal processes, considered as univariate point processes, are all regular. The multivariate process is said to be *regular* if the pooled process is regular. Throughout this paper, we will assume that the multivariate process {$\mathbf{N}(t)$, $t \ge 0$} of our interest is a regular process. In the following, we shall use the notation $\Phi \sim \mathcal{GG}(v, k, \alpha, l)$ to represent that the continuous random variable $\Phi$ follows the GGD with parameters $(v, k, \alpha, l)$ and {$N(t)$, $t \ge 0$} $\sim \mathcal{NHPP}(\eta(t))$ to indicate that the counting process {$N(t)$, $t \ge 0$} follows the nonhomogeneous Poisson process (NHPP) with intensity function $\eta(t)$.

**Definition 2.1. (Bivariate Poisson Generalized Gamma Process)** A bivariate counting process {$\mathbf{N}(t)$, $t \ge 0$} is called the bivariate Poisson generalized gamma process (BPGGP) with the set of parameters $(\lambda_1(t), \lambda_2(t), v, k, \alpha, l)$, $\lambda_i(t) > 0, \forall t \ge 0, i = 1, 2, v \ge 0, k, \alpha, l > 0$, if

(i) {$N_i(t)$, $t \ge 0$} | $(\Phi = \phi) \sim \mathcal{NHPP}(\phi\lambda_i(t))$, $i = 1, 2$, *independent*;

(ii) $\Phi \sim \mathcal{GG}(v, k, \alpha, l)$.

Throughout this paper, the BPGGP with the set of parameters $(\lambda_1(t), \lambda_2(t), v, k, \alpha, l)$ will be denoted by BPGGP($\lambda_1(t)$, $\lambda_2(t)$, $v, k, \alpha, l$).

Based on Definition 2.1, we now derive some basic properties of BPGGP and with that aim, let us introduce $\Lambda_i(t) \equiv \int_0^t \lambda_i(s)ds$, $i = 1, 2, t \ge 0$.

**Proposition 2.2.**

(i) {$\mathbf{M}(t) = (M_1(t), M_2(t))$, $t \ge 0$} is a BPGGP(1, 1, v, k, \alpha, l) if and only if {$\mathbf{N}(t) = (M_1(\Lambda_1(t)), M_2(\Lambda_2(t)))$, $t \ge 0$} is a BPGGP($\lambda_1(t)$, $\lambda_2(t)$, v, k, \alpha, l).

(ii) For $c > 0$, let $\tilde{\alpha} = \alpha c$, $\tilde{l} = l/c$ and $\tilde{\lambda}_i(t) = \lambda_i(t)/c$ for $i = 1, 2$. Then, a BPGGP($\tilde{\lambda}_1(t)$, $\tilde{\lambda}_2(t)$, v, k, $\tilde{\alpha}$, $\tilde{l}$) is a BPGGP($\lambda_1(t)$, $\lambda_2(t)$, v, k, \alpha, l).

The proof is similar to that of Proposition 1 from Cha and Mercier [9] in the univariate case and it is omitted. As in that paper, the second point of Proposition 2.2 shows that the BPGGP model as given in Definition 2.1 is not identifiable. An additional constraint should hence be added such as $l \equiv 1$ for instance, wherever statistical procedures are studied (which is not the case in the present paper).

We now study the distributions for the numbers of events, which are of major interest for any counting process model.
**Theorem 2.3.** For $0 \le u_{i1} < u_{i2} < \cdots < u_{im}, i = 1, 2,$

$$
P(N_i(u_{i2}) - N_i(u_{i1}) = n_i, i = 1, 2) \\
= \left[ \prod_{i=1}^{2} \frac{(\Lambda_i(u_{i2}) - \Lambda_i(u_{i1}))^{n_i}}{n_i!} \right] \frac{\alpha^{k-\nu}}{(\alpha + \sum_{i=1}^{2} (\Lambda_i(u_{i2}) - \Lambda_i(u_{i1})))^{k+n_1+n_2-\nu}} \\
\times \frac{\Gamma_\nu(k+n_1+n_2, (\alpha + \sum_{i=1}^{2} (\Lambda_i(u_{i2}) - \Lambda_i(u_{i1})))l)}{\Gamma_\nu(k, \alpha l)},
$$

and

$$
P(N_i(u_{ij}) - N_i(u_{ij-1}) = n_{ij}, i = 1, 2, j = 1, 2, \dots, m)
$$
$$
= \left[ \prod_{i=1}^{2} \prod_{j=1}^{m} \frac{(\Lambda_i(u_{ij}) - \Lambda_i(u_{ij-1}))^{n_{ij}}}{n_{ij}!} \right] \frac{\alpha^{k-\nu}}{(\alpha + \sum_{i=1}^{2} \sum_{j=1}^{m} (\Lambda_i(u_{ij}) - \Lambda_i(u_{ij-1})))^{k+\sum_{i=1}^{2} \sum_{j=1}^{m} n_{ij}-\nu}} \\
\times \frac{\Gamma_{\nu}(k + \sum_{i=1}^{2} \sum_{j=1}^{m} n_{ij}, (\alpha + \sum_{i=1}^{2} \sum_{j=1}^{m} (\Lambda_i(u_{ij}) - \Lambda_i(u_{ij-1})))l)}{\Gamma_{\nu}(k, \alpha l)}.
$$

*Proof.* From the definition of BPGGP(λ₁(t), λ₂(t), ν, k, α, l),

$$
\begin{align*}
& P(N_i(u_{i2}) - N_i(u_{i1}) = n_i, i = 1, 2) \\
&= \int_0^\infty \left[ \prod_{i=1}^2 \frac{(\phi(\Lambda_i(u_{i2}) - \Lambda_i(u_{i1}))^{n_i} \exp\{-\phi(\Lambda_i(u_{i2}) - \Lambda_i(u_{i1}))\}}{n_i!} \right] \frac{\alpha^{k-\nu}}{\Gamma_\nu(k, \alpha l)} \frac{\phi^{k-1} \exp\{-\alpha\phi\}}{(\phi+l)^\nu} d\phi \\
&= \left[ \prod_{i=1}^2 \frac{(\Lambda_i(u_{i2}) - \Lambda_i(u_{i1}))^{n_i}}{n_i!} \right] \frac{\alpha^{k-\nu}}{\Gamma_\nu(k, \alpha l)} \\
&\quad \times \int_0^\infty \frac{\phi^{k+n_1+n_2-1} \exp\{-\phi(\alpha + \sum_{i=1}^2 (\Lambda_i(u_{i2}) - \Lambda_i(u_{i1})))\}}{(\phi+l)^\nu} d\phi \\
&= \left[ \prod_{i=1}^2 \frac{(\Lambda_i(u_{i2}) - \Lambda_i(u_{i1}))^{n_i}}{n_i!} \right] \frac{\alpha^{k-\nu}}{(\alpha + \sum_{i=1}^2 (\Lambda_i(u_{i2}) - \Lambda_i(u_{i1})))^{k+n_1+n_2-\nu}} \\
&\quad \times \frac{\Gamma_\nu(k+n_1+n_2, (\alpha + \sum_{i=1}^2 (\Lambda_i(u_{i2}) - \Lambda_i(u_{i1})))l)}{\Gamma_\nu(k, \alpha l)}.
\end{align*}
$$

In a similar way,

$$
\begin{align*}
& P(N_i(u_{ij}) - N_i(u_{ij-1}) = n_{ij}, i = 1, 2, j = 1, 2, \dots, m) \\
&= \int_0^\infty \left[ \prod_{i=1}^2 \prod_{j=1}^m \frac{(\phi(\Lambda_i(u_{ij}) - \Lambda_i(u_{ij-1}))^{n_{ij}} \exp\{-\phi(\Lambda_i(u_{ij}) - \Lambda_i(u_{ij-1}))\}}{n_{ij}!} \right] \\
&\qquad \times \frac{\alpha^{k-\nu}}{\Gamma_\nu(k, \alpha l)} \frac{\phi^{k-1} \exp\{-\alpha\phi\}}{(\phi+l)^\nu} d\phi \\
&= \left[ \prod_{i=1}^2 \prod_{j=1}^m \frac{(\Lambda_i(u_{ij}) - \Lambda_i(u_{ij-1}))^{n_{ij}}}{n_{ij}!} \right] \frac{\alpha^{k-\nu}}{(\alpha + \sum_{i=1}^2 \sum_{j=1}^m (\Lambda_i(u_{ij}) - \Lambda_i(u_{ij-1})))^{k+\sum_{i=1}^2 \sum_{j=1}^m n_{ij}-\nu}} \\
&\qquad \times \frac{\Gamma_\nu(k + \sum_{i=1}^2 \sum_{j=1}^m n_{ij}, (\alpha + \sum_{i=1}^2 \sum_{j=1}^m (\Lambda_i(u_{ij}) - \Lambda_i(u_{ij-1})))l)}{\Gamma_\nu(k, \alpha l)}.
\end{align*}
$$
The joint moments of $(N_1(t), N_2(t))$ also are of practical interest for the applications. They are obtained in the following theorem.

**Theorem 2.4.** Let {$N(t), t \ge 0$} be the BPGGP with the set of parameters $(\lambda_1(t), \lambda_2(t), v, k, \alpha, l)$. Then the following properties hold.

(i) The joint moment generating function of $(N_1(t), N_2(t))$ is given by

$$M_{N_1(t_1),N_2(t_2)}(s_1,s_2) = \left(1 - \frac{(\Lambda_1(t_1)(e^{s_1}-1) + \Lambda_2(t_2)(e^{s_2}-1))}{\alpha}\right)^{v-k} \times \frac{\Gamma_v[k, (\alpha - (\Lambda_1(t_1)(e^{s_1}-1) + \Lambda_2(t_2)(e^{s_2}-1)))l]}{\Gamma_v(k, \alpha l)}$$

for all $s_1$ and $s_2$ such that $\Lambda_1(t_1)(e^{s_1}-1) + \Lambda_2(t_2)(e^{s_2}-1) < \alpha$.

(ii) We have

$$E[N_1(t_1)^{r_1} N_2(t_2)^{r_2}] = \sum_{i_1=0}^{r_1} \sum_{i_2=0}^{r_2} \frac{(\Lambda_1(t_1))^{i_1} (\Lambda_2(t_2))^{i_2}}{\alpha^{i_1+i_2}} \frac{\Gamma_v(k+i_1+i_2, \alpha l)}{\Gamma_v(k, \alpha l)} \begin{Bmatrix} r_1 \\ i_1 \end{Bmatrix} \begin{Bmatrix} r_2 \\ i_2 \end{Bmatrix}$$

for all $r_1, r_2 \in \mathbb{N}$, where the {$\text{braces}$} denote Stirling numbers of the second kind.

(iii) The covariance of $(N_1(t_1), N_2(t_2))$ is given by

$$Cov(N_1(t_1), N_2(t_2)) = \frac{\Lambda_1(t_1)\Lambda_2(t_2)}{\alpha^2} \left[ \frac{\Gamma_v(k+2, \alpha l)}{\Gamma_v(k, \alpha l)} - \left( \frac{\Gamma_v(k+1, \alpha l)}{\Gamma_v(k, \alpha l)} \right)^2 \right], \quad (2.5)$$

and the corresponding Pearson's correlation coefficient is:

$$\rho_{(N_1(t_1), N_2(t_2))} = \left( \sqrt{1 + \frac{C(k, \alpha, l, v)}{\Lambda_1(t_1)}} \sqrt{1 + \frac{C(k, \alpha, l, v)}{\Lambda_2(t_2)}} \right)^{-1},$$

where

$$C(k, \alpha, l, v) = \alpha \left( \frac{\Gamma_v(k+2, \alpha l)}{\Gamma_v(k+1, \alpha l)} - \frac{\Gamma_v(k+1, \alpha l)}{\Gamma_v(k, \alpha l)} \right)^{-1}.$$

Proof.

(i) Conditioning on $\Phi$, we can write

$$M_{N_1(t_1),N_2(t_2)}(s_1,s_2) = E \left[ E \left( e^{s_1 N_1(t_1) + s_2 N_2(t_2)} | \Phi \right) \right],$$

where $[N_i(t)|\Phi = \phi]$, $i = 1, 2$ are independent and Poisson distributed with parameter $\Lambda_i(t_i)\phi$, respectively. By using the moment generating function of a Poisson distribution (see e.g., Ross [23]), we have:
$$
\begin{align*}
M_{N_1(t_1), N_2(t_2)}(s_1, s_2) &= E[\exp\{(\Lambda_1(t_1)(e^{s_1}-1) + \Lambda_2(t_2)(e^{s_2}-1))\Phi\}] \\
&= M_{\Phi}[(\Lambda_1(t_1)(e^{s_1}-1) + \Lambda_2(t_2)(e^{s_2}-1))] \\
&= \left(1 - \frac{(\Lambda_1(t_1)(e^{s_1}-1) + \Lambda_2(t_2)(e^{s_2}-1))}{\alpha}\right)^{\nu-k} \\
&\quad \times \frac{\Gamma_{\nu}[k, (\alpha - (\Lambda_1(t_1)(e^{s_1}-1) + \Lambda_2(t_2)(e^{s_2}-1)))l]}{\Gamma_{\nu}(k, \alpha l)}
\end{align*}
$$

for all $s_1$ and $s_2$ such that $\Lambda_1(t_1)(e^{s_1}-1) + \Lambda_2(t_2)(e^{s_2}-1) < \alpha$, due to (2.3).

(ii) By a similar procedure, we have:

$$
\begin{align*}
E[N_1(t_1)^{r_1} N_2(t_2)^{r_2}] &= E[E(N_1(t_1)^{r_1} N_2(t_2)^{r_2} | \Phi)] \\
&= E[E(N_1(t_1)^{r_1} | \Phi) E(N_2(t_2)^{r_2} | \Phi)],
\end{align*}
$$

due to the conditional independence between $N_1(t_1)$ and $N_2(t_2)$ given $\Phi$. Based on Formula (3.4) in Riordan [22], which provides the moments of Poisson distributions with respect to Stirling numbers of the second kind, we now get:

$$
E[N_1(t_1)^{r_1} N_2(t_2)^{r_2}] = E \left[ \sum_{i_1=0}^{r_1} \sum_{i_2=0}^{r_2} (\Lambda_1(t_1))^{i_1} (\Lambda_2(t_2))^{i_2} \Phi^{i_1+i_2} \begin{Bmatrix} r_1 \\ i_1 \end{Bmatrix} \begin{Bmatrix} r_2 \\ i_2 \end{Bmatrix} \right].
$$

The result follows, based on (2.4).

(iii) The computation of the covariance now is a direct consequence from point (ii). The variance of each $N_i(t_i)$, $i = 1, 2$, can be derived in the same way from point (ii) (or using Theorem 2 from Cha and Mercier [9], based on the fact that $\{N_i(t), t \ge 0\}$, $i = 1, 2$, are univariate Poisson generalized gamma processes, see Proposition 2.7 later on). We obtain

$$
\begin{equation}
\begin{split}
\mathrm{Var}(N_i(t_i)) ={}& \frac{\Lambda_i(t_i) \Gamma_v(k+1, \alpha l)}{\alpha \Gamma_v(k, \alpha l)} \\
& + \left(\frac{\Lambda_i(t_i)}{\alpha}\right)^2 \left[ \frac{\Gamma_v(k+2, \alpha l)}{\Gamma_v(k, \alpha l)} - \left(\frac{\Gamma_v(k+1, \alpha l)}{\Gamma_v(k, \alpha l)}\right)^2 \right],
\end{split}
\tag{2.6}
\end{equation}
$$

for $i = 1, 2$. The result for Pearson's correlation coefficient then is a routine computation, remembering that

$$
\rho(N_1(t_1), N_2(t_2)) = \frac{\mathrm{Cov}(N_1(t_1), N_2(t_2))}{\sqrt{\mathrm{Var}(N_1(t_1))} \sqrt{\mathrm{Var}(N_2(t_2))}}.
$$

□

**Remark 2.5.** Based on (2.6) and (2.5), it is clear that $\mathrm{Cov}(N_1(t_1), N_2(t_2)) > 0$, which shows that $N_1(t_1)$ and $N_2(t_2)$ always are positively correlated, whatever the parameters of the BPGGP are, and whatever the times $t_1$ and $t_2$ are. Also, $\rho_{(N_1(t_1), N_2(t_2))}$ is non decreasing with respect to $(t_1, t_2)$, with

$$
\lim_{(t_1, t_2) \to (0,0)^+} \rho(N_1(t_1), N_2(t_2)) = 0,
$$

$$
\lim_{(t_1, t_2) \to (\infty, \infty)} \rho(N_1(t_1), N_2(t_2)) = 1.
$$

Now, in the following proposition, some important properties of the BPGGP will be stated. For this, we employ the same notations as those in Cha and Giorgio [8]. Denote by $\mathcal{H}_{Pt-} \equiv \{M(u), 0 \le u < t\}$
the history of the pooled process in $[0, t)$. Define $M(t-)$ as the total number of events in $[0, t)$ and $T_i$ as the time from 0 until the arrival of the ith event in $[0, t)$ of the pooled process $\{M(t), t \ge 0\}$. Then $\mathcal{H}_{P_{i-}}$ can equivalently be defined in terms of $M(t-)$ and the sequential arrival points of the events $0 \le T_1 \le T_2 \le \dots \le T_{M(t-)} < t$ in $[0, t)$. Similarly, define the marginal histories of the marginal processes $\mathcal{H}_{it-} \equiv \{N_i(u), 0 \le u < t\}, i = 1, 2$. Then, $\mathcal{H}_{it-} \equiv \{N_i(u), 0 \le u < t\}$ can also be defined in terms of $N_i(t-)$ and the sequential arrival points of the events $0 \le T_{i1} \le T_{i2} \le \dots \le T_{iN_i(t-)} < t$ in $[0, t), i = 1, 2$, where $N_i(t-)$ is the total number of events of type i point process in $[0, t), i = 1, 2$. In the following, we also use the definition of “p(t)-thinning” in Cha and Giorgio [8]. Also, we define univariate Poisson generalized Gamma process to define the marginal process.

**Definition 2.6. (Poisson Generalized Gamma Process)** A counting process $\{N(t), t \ge 0\}$ is called the Poisson generalized gamma process (PGGP) with the set of parameters $(\lambda(t), v, k, \alpha, l)$, $\lambda(t) > 0, \forall t \ge 0, v \ge 0, k, \alpha, l > 0$, if

(i) $\{N(t), t \ge 0\}|(\Phi = \phi) \sim \mathcal{N}\mathcal{H}\mathcal{P}\mathcal{P}(\phi\lambda(t));$

(ii) $\Phi \sim \mathcal{GG}(v, k, \alpha, l)$.

See Cha and Mercier [9] for various properties of the PGGP.

For convenience, we now introduce the following notations: $N_{ui}(t) \equiv N_i(u+t) - N_i(u)$, $\Lambda_i(t) \equiv \int_0^t \lambda_i(u)du$, $i = 1, 2$, $\lambda(t) \equiv \lambda_1(t) + \lambda_2(t)$, $\Lambda(t) = \int_0^t \lambda(u)du = \Lambda_1(t) + \Lambda_2(t)$, and $p_i(t) \equiv \lambda_i(t)/\lambda(t)$, $i = 1, 2$.

**Proposition 2.7.** Let $\{N(t), t \ge 0\}$ be the BPGGP with the set of parameters $(\lambda_1(t), \lambda_2(t), v, k, \alpha, l)$. Then

(i) The pooled process $\{M(t), t \ge 0\}$ is PGGP($\lambda(t), v, k, \alpha, l$).

(ii) The process $\{N(t), t \ge 0\}$ is constructed by thinning of $\{M(t), t \ge 0\}$ with thinning probabilities $p_i(t), i = 1, 2$: $\{(M_{p_1(t)}, t), M_{p_2(t)}(t)), t \ge 0\}$.

(iii) Given $(\mathcal{H}_{1u-}, \mathcal{H}_{2u-})$, $\{N_u(t), t \ge 0\}$, where $N_u(t) = (N_{u1}(t), N_{u2}(t))$, is BPGGP($\lambda_1(t+u), \lambda_2(t+u)$, $v, k+n_1+n_2, \alpha+\Lambda(u), l$), where $n_i$ is the realization of $N_i(t-)$, $i = 1, 2$, respectively.

(iv) For any fixed $u \ge 0$, $\{N_u(t), t \ge 0\}$ is “unconditionally” BPGGP($\lambda_1(t+u), \lambda_2(t+u), v, k, \alpha, l$).

(v) The marginal process $\{N_i(t), t \ge 0\}$ is PGGP($\lambda_i(t), v, k, \alpha, l)$, $i = 1, 2$.

*Proof.* Properties (i), (ii), (iv) and (v) obviously hold. Let $\{N(t), t \ge 0\}$ be PGGP($\lambda(t), v, k, \alpha, l$). Then, at an arbitrary time $u > 0$, given $\{N(u-) = n, T_1 = t_1, T_2 = t_2, \dots, T_n = t_n\}$, the conditional future process $\{N_u(t), t \ge 0\}$, where $N_u(t) \equiv N(u+t) - N(u)$, is a PGGP with the set of parameters $(\lambda(u+t), v, k+n, \alpha+\Lambda(u), l)$ (see Cha and Mercier [9]). Then property (iii) also obviously holds due to property (ii). $\square$

Properties (iii) and (iv) state about the conditional and unconditional restarting properties of the BPGGP. For more details on the restarting property, see Cha [7] and Cha and Giorgio [8].

### 3. Further properties of bivariate Poisson generalized gamma process

An efficient characterization for a multivariate point process can be done through the stochastic intensity approach (see Cox and Lewis [10], Cha and Giorgio [8]). As mentioned in Cha and Giorgio [8], a marginally regular bivariate process can be specified by the following *complete intensity functions*:
$$
\begin{align}
\lambda_{1t} &\equiv \lim_{\Delta t \to 0} \frac{P(N_1(t, t + \Delta t) \ge 1 | \mathcal{H}_{1t-}; \mathcal{H}_{2t-})}{\Delta t} = \lim_{\Delta t \to 0} \frac{P(N_1(t, t + \Delta t) = 1 | \mathcal{H}_{1t-}; \mathcal{H}_{2t-})}{\Delta t}, \nonumber \\
\lambda_{2t} &\equiv \lim_{\Delta t \to 0} \frac{P(N_2(t, t + \Delta t) \ge 1 | \mathcal{H}_{1t-}; \mathcal{H}_{2t-})}{\Delta t} = \lim_{\Delta t \to 0} \frac{P(N_2(t, t + \Delta t) = 1 | \mathcal{H}_{1t-}; \mathcal{H}_{2t-})}{\Delta t}, \nonumber \\
\lambda_{12t} &\equiv \lim_{\Delta t \to 0} \frac{P(N_1(t, t + \Delta t) N_2(t, t + \Delta t) \ge 1 | \mathcal{H}_{1t-}; \mathcal{H}_{2t-})}{\Delta t} \tag{3.1}
\end{align}
$$

where $N_i(t_1, t_2)$, $t_1 < t_2$, represents the number of events in $[t_1, t_2)$, $i = 1, 2$, respectively (see Cox and Lewis [10]). For a regular process, $\lambda_{12t} = 0$, and it is sufficient to specify just $\lambda_{1t}$ and $\lambda_{2t}$ in (3.1) in order to define a regular process.

**Theorem 3.1.** The complete intensity functions of the BPGGP with the set of parameters $(\lambda_1(t), \lambda_2(t), v, k, \alpha, l)$ are given by

$$ \lambda_{it} = \frac{1}{(\alpha + \Lambda_1(t) + \Lambda_2(t))} \frac{\Gamma_v(k + N_1(t-) + N_2(t-) + 1, (\alpha + \Lambda_1(t) + \Lambda_2(t))l)}{\Gamma_v(k + N_1(t-) + N_2(t-) - (\alpha + \Lambda_1(t) + \Lambda_2(t))l)} \lambda_i(t), \quad i = 1, 2. \quad (3.2) $$

*Proof.* Observe that

$$
\begin{align*}
\lambda_{1t} &= \lim_{\Delta t \to 0} \frac{P(N_1(t, t + \Delta t) = 1 | \mathcal{H}_{1t-}; \mathcal{H}_{2t-})}{\Delta t} \\
&= E_{(\Phi | \mathcal{H}_{1t-}; \mathcal{H}_{2t-})} \left[ \lim_{\Delta t \to 0} \frac{P(N_1(t, t + \Delta t) = 1 | \Phi; \mathcal{H}_{1t-}; \mathcal{H}_{2t-})}{\Delta t} \right],
\end{align*}
$$

where $E_{(\Phi | \mathcal{H}_{1t-}; \mathcal{H}_{2t-})} [\cdot]$ stands for the expectation with respect to the conditional distribution of $(\Phi | \mathcal{H}_{1t-}; \mathcal{H}_{2t-})$ and

$$ \lim_{\Delta t \to 0} \frac{P(N_1(t, t + \Delta t) = 1 | \Phi; \mathcal{H}_{1t-}; \mathcal{H}_{2t-})}{\Delta t} = \Phi \lambda_1(t). $$

Thus, $\lambda_{1t} = E_{(\Phi | \mathcal{H}_{1t-}; \mathcal{H}_{2t-})} [\Phi \lambda_1(t)]$.

Similar to the procedure described in Cha [7], the conditional distribution of $(\Phi | \mathcal{H}_{it-} = \mathbf{h}_{it-}, i = 1, 2)$, where $\mathbf{h}_{it-} = (t_{i1}, t_{i2}, \dots, t_{in_i}, n_i)$ is the realization of $\mathcal{H}_{it-}$, $i = 1, 2$, respectively, is given by

$$
\left( \phi^{n_1+n_2} \exp \left\{ -\phi \int_0^t \sum_{i=1}^2 \lambda_i(x) dx \right\} f(\phi) \right) \times \left( \int_0^\infty \phi^{n_1+n_2} \exp \left\{ -\phi \int_0^t \sum_{i=1}^2 \lambda_i(x) dx \right\} f(\phi) d\phi \right)^{-1},
$$

where we recall that $f$ stands for the probability density function of $\Phi$. Then,

$$
\begin{align*}
\lambda_{1t} &= E_{(\Phi | \mathcal{H}_{1t-}; \mathcal{H}_{2t-})} [\Phi \lambda_1(t)] \\
&= \frac{\int_0^\infty \phi^{n_1+n_2+1} \exp \left\{ -\phi \int_0^t \sum_{i=1}^2 \lambda_i(x) dx \right\} f(\phi) d\phi}{\int_0^\infty \phi^{n_1+n_2} \exp \left\{ -\phi \int_0^t \sum_{i=1}^2 \lambda_i(x) dx \right\} f(\phi) d\phi} \lambda_1(t),
\end{align*}
$$
which extends Proposition 4.1 in Grandell [16] to a bivariate and non homogeneous mixed Poisson process. Now, one can check that

$$ \int_0^\infty \phi^n \exp \left\{ -\phi \int_0^t \sum_{i=1}^2 \lambda_i(x) dx \right\} f(\phi) d\phi \\ = \frac{\alpha^{k-\nu}}{(\alpha + \Lambda_1(t) + \Lambda_2(t))^{k+n-\nu}} \frac{\Gamma_\nu(k+n, (\alpha + \Lambda_1(t) + \Lambda_2(t))l)}{\Gamma_\nu(k, \alpha l)}. $$

Therefore,

$$ \lambda_{1t} = \frac{1}{(\alpha + \Lambda_1(t) + \Lambda_2(t))} \frac{\Gamma_\nu(k + N_1(t-) + N_2(t-) + 1, (\alpha + \Lambda_1(t) + \Lambda_2(t))l)}{\Gamma_\nu(k + N_1(t-) + N_2(t-) + 1, (\alpha + \Lambda_1(t) + \Lambda_2(t))l)} \lambda_1(t). $$

The intensity function $\lambda_{2t}$ can be obtained symmetrically. □

**Proposition 3.2.** Let {$N(t), t \ge 0$} be a BPGGP with the set of parameters {$\lambda_1(t), \lambda_2(t), v = 0, k, \alpha, l$} such that $\lambda_i(t) \equiv \phi_i(t) \exp\{\alpha(\Phi_1(t) + \Phi_2(t))\}, i = 1, 2$, where $\Phi_i(t) \equiv \int_0^t \phi_i(s)ds, i = 1, 2$. Then {$N(t), t \ge 0$} is a BVGPP($\phi_1(t), \phi_2(t), 1/\alpha, k/\alpha$) (whatever $l$ is).

*Proof.* Under the given specific setting, it can be shown that the complete intensity functions in (3.2) becomes those in BVGPP given in Definition 1 of Cha and Girogio [8]. □

The result of Proposition 3.2 is also clear from the definition of the BPGGP and the characterization of the BVGPP in Theorem 2 of Cha and Giorgio [8].

It can be shown that

$$ \eta(k) \equiv \frac{\Gamma_\nu(k+1, \alpha l)}{\Gamma_\nu(k, \alpha l)}, $$

is increasing in $k > 0$ for any $v \ge 0, \alpha, l > 0$ (see the proof of Proposition 8 in Cha and Mercier [9]). Thus, the complete intensity functions in (3.2) are increasing in $N_1(t-) + N_2(t-)$ . This implies that the proneness to the future event occurrence in each marginal process is increasing with the number of events, previously occurred in the pooled process.

## 4. Multivariate Poisson generalized gamma process

In this section, we study the multivariate Poisson generalized gamma process (MPGGP), by extending the results obtained in the previous sections. As in Cha and Giorgio [8], let {$N(t), t \ge 0$}, where $N(t) = (N_1(t), N_2(t), \dots, N_m(t))$, be a multivariate process and define the corresponding “pooled” point process {$M(t), t \ge 0$}, where $M(t) = N_1(t) + N_2(t) + \dots + N_m(t)$. Also, define the marginal point processes {$N_i(t), t \ge 0$}, $i = 1, 2, \dots, m$, and the corresponding marginal histories of the marginal processes: $H_{it-}, i = 1, 2, \dots, m$. The MPGGP can be defined by generalizing Definition 2.1.

**Definition 4.1. (Multivariate Poisson Generalized Gamma Process)** A multivariate counting process {$N(t), t \ge 0$} is called the multivariate Poisson generalized gamma process (MPGGP) with the set of parameters ($\lambda_i(t), i = 1, 2, \dots, m, v, k, \alpha, l$), where $\lambda_i(t) > 0, \forall t \ge 0, i = 1, 2, \dots, m, v \ge 0, k, \alpha, l > 0$, if

(i) {$N_i(t), t \ge 0$} | ($\Phi = \phi$) ~ $NHP\mathcal{P}(\phi\lambda_i(t))$, $i = 1, 2, \dots, m$, independent;

(ii) $\Phi \sim \mathcal{GG}(v, k, \alpha, l)$.
To state the properties of the MPGPP, we define $\lambda(t) = \sum_{i=1}^{m} \lambda_i(t)$, $\Lambda(t) = \int_{0}^{t} \lambda(v)dv$, $p_i(t) = \lambda_i(t)/\lambda(t)$, and $N_{ui}(t) \equiv N_i(u+t) - N_i(u)$, $i = 1, 2, \dots, m$.

**Proposition 4.2.** Let {$\mathbf{N}(t)$, $t \ge 0$} be the MPGGP with the set of parameters ($\lambda_i(t)$, $i = 1, 2, \dots, m, v, k, \alpha, l$). Then

(i) The pooled process {$M(t)$, $t \ge 0$} is MPGGP($\lambda(t)$, $v, k, \alpha, l$).

(ii) The process {$\mathbf{N}(t)$, $t \ge 0$} is constructed by thinning of {$M(t)$, $t \ge 0$} with thinning probabilities $p_i(t)$, $i = 1, 2, \dots, m$: {$(M_{p_1}(t), M_{p_2}(t), \dots, M_{p_m}(t))$, $t \ge 0$}.

(iii) Given ($\mathcal{H}_{iu-}$, $i = 1, 2, \dots, m$), {$\mathbf{N}_u(t)$, $t \ge 0$}, where $\mathbf{N}_u(t) = (N_{u1}(t), N_{u2}(t), \dots, N_{um}(t))$, is MPGGP($\lambda_i(t+u)$, $i = 1, 2, \dots, m, v, k, \sum_{i=1}^{m} n_i, \alpha + \Lambda(u), l$), where $n_i$ is the realization of $N_i(t-)$, $i = 1, 2, \dots, m$, respectively.

(iv) For any fixed $u \ge 0$, {$\mathbf{N}_u(t)$, $t \ge 0$} is “unconditionally” MPGGP($\lambda_i(t+u)$, $i = 1, 2, \dots, m, v, k, \alpha, l$).

(v) The marginal process {$N_i(t)$, $t \ge 0$} is MPGGP($\lambda_i(t)$, $v, k, \alpha, l$), $i = 1, 2, \dots, m$.

In addition to the basic properties stated in Proposition 4.2, it can be shown that the complete intensity functions are given by:

$$ \lambda_{it} \equiv \lim_{\Delta t \to 0} \frac{P(N_i(t, t + \Delta t) = 1 | \mathcal{H}_{1t-}; \mathcal{H}_{2t-}; \dots; \mathcal{H}_{mt-})}{\Delta t} \\ = \frac{1}{(\alpha + \sum_{i=1}^{m} \Lambda_i(t))} \frac{\Gamma_{\nu}(k + \sum_{i=1}^{m} N_i(t-) + 1, (\alpha + \sum_{i=1}^{m} \Lambda_i(t))l)}{\Gamma_{\nu}(k + \sum_{i=1}^{m} N_i(t-) - 1, (\alpha + \sum_{i=1}^{m} \Lambda_i(t))l)} \lambda_i(t), i = 1, 2, \dots, m. $$

Furthermore, it can also be shown that

$$ P(N_i(u_{i2}) - N_i(u_{i1}) = n_i, i = 1, 2, \dots, m) \\ = \left[ \prod_{i=1}^{m} \frac{(\Lambda_i(u_{i2}) - \Lambda_i(u_{i1}))^{n_i}}{n_i!} \right] \frac{\alpha^{k-v}}{(\alpha + \sum_{i=1}^{m} (\Lambda_i(u_{i2}) - \Lambda_i(u_{i1})))^{k+\sum_{i=1}^{m} n_i - v}} \\ \times \frac{\Gamma_{\nu}(k + \sum_{i=1}^{m} n_i, (\alpha + \sum_{i=1}^{m} (\Lambda_i(u_{i2}) - \Lambda_i(u_{i1})))l)}{\Gamma_{\nu}(k, \alpha l)}. $$

In addition to the above results, other properties obtained in the previous sections could be extended to the multivariate case in similar ways.

**5. Comparison and monotony results with respect to the multivariate likelihood ratio ordering**

The results of the previous sections were dependent on the specific properties of the generalized gamma distribution. We here provide some other results of MPGGPs based on the notion of multivariate likelihood ratio ordering, which hold in the more general setting of a multivariate mixed Poisson process which we now define.

**Definition 5.1.** Let {$\mathbf{N}(t)$, $t \ge 0$} be a multivariate counting process, where $\mathbf{N}(t) = (N_1(t), N_2(t), \dots, N_m(t))$, for all $t \ge 0$ and let $\Phi$ be a non negative random variable, which is assumed to be absolutely continuous with respect to Lebesgue measure. Then {$\mathbf{N}(t)$, $t \ge 0$} is called the multivariate mixed Poisson process (MMPP) with the set of parameters ($\lambda_i(t)$, $i = 1, 2, \dots, m, \Phi$), where $\lambda_i(t) > 0$, $\forall t \ge 0$, $i = 1, 2, \dots, m$, if {$N_i(t)$, $t \ge 0$} |($\Phi = \phi$) ~ $\mathcal{N}\mathcal{H}\mathcal{P}(\phi\lambda_i(t))$, $i = 1, 2, \dots, m$ and are independent.
**Remark 5.2.** Based on its definition, a MPGGP is a specific MMPP, where $\Phi$ has a generalized gamma distribution.

In order to write down the different properties below, we now recall the notion of multivariate likelihood ratio ordering. We refer to Karlin and Rinott [18] for more details.

**Definition 5.3.** Let $X$ and $Y$ be two random vectors on $\mathbb{R}^n$ with respective density $f_X$ and $f_Y$ with respect to a common product measure $\sigma$ ($dX$). (For instance $X$ and $Y$ may be both absolutely continuous or both discrete). Then, $X$ is said to be smaller than $Y$ in the multivariate likelihood ratio ordering (written $X <_{\text{lr}} Y$) as soon as

$$f_X(\mathbf{x})f_Y(\mathbf{y}) \le f_X(\mathbf{x} \wedge \mathbf{y})f_Y(\mathbf{x} \vee \mathbf{y}),$$

for all $\mathbf{x}, \mathbf{y} \in \mathbb{R}^n$, where the minimum $\wedge$ and the maximum $\vee$ are taken componentwise.

In the univariate case, the likelihood ratio order is denoted by $<_{lr}$.

**Definition 5.4.** Let $X$ be a random vector on $\mathbb{R}^n$ with density $f_X$ with respect to a product measure $\sigma$ ($dX$). Then, $X$ is said to be Multivariate Totally Positive property of order 2 (MTP2) as soon as

$$f_X(\mathbf{x})f_Y(\mathbf{y}) \le f_X(\mathbf{x} \wedge \mathbf{y})f_Y(\mathbf{x} \vee \mathbf{y}),$$

for all $\mathbf{x}, \mathbf{y} \in \mathbb{R}^n$, which is just equivalent to $X <_{\text{lr}} X$.

We first provide a result that extends Proposition 4 in Cha and Mercier [9] to the multivariate setting.

**Proposition 5.5.** Let $\{N(t), t \ge 0\}$ be a MMPP with set of parameters $(\lambda_i(t), i = 1, 2, \dots, m, \Phi)$. Then $N(t)$ increases with respect to $t$ in the multivariate likelihood ratio ordering.

*Proof.* Let $0 \le t_1 < t_2$. Let us show that $N(t_1) <_{\text{lr}} N(t_2)$.

Let

$$g_j(n_1, \dots, n_m | \phi) \equiv P(N_i(t_j) = n_i, i = 1, 2, \dots, m | \Phi = \phi) \\ = \prod_{i=1}^{m} P(N_i(t_j) = n_i | \Phi = \phi),$$

for $j = 1, 2$ and $n_i \in \mathbb{N}, i = 1, 2, \dots, m$, where $P(N_i(t_j) = |\Phi = \phi)$ is the Poisson distribution with parameter $\Lambda_i(t_j)$. As this distribution increases with respect to $\Lambda_i(t_j)$ in the likelihood ratio ordering and as $\Lambda_i(t_1) \le \Lambda_i(t_2)$, we derive that $P(N_i(t_1) = |\Phi = \phi) <_{lr} P(N_i(t_2) = |\Phi = \phi)$, and next that $g_1(\cdots|\phi) <_{\text{lr}} g_2(\cdots|\phi)$, as the multivariate likelihood ratio ordering is stable through conjunction (see Shaked and Shanthikumar [24], Theorem 6.E.4(a) page 299).

Using that $\Phi <_{lr} \Phi$, we derive from Theorem 2.4 in Karlin and Rinott [18] that

$$\int_{\mathbb{R}_+} g_1(\cdots|\phi) f_{\Phi}(\phi) d\phi <_{\text{lr}} \int_{\mathbb{R}_+} g_2(\cdots|\phi) f_{\Phi}(\phi) d\phi,$$

which is just equivalent to $N(t_1) <_{\text{lr}} N(t_2)$, and allows to conclude. $\square$

We next show that the marginal increments in a MMPP (taken at possibly different times for each margin) exhibit the MTP2 property.
**Proposition 5.6.** Let {$\mathbf{N}(t)$, $t \ge 0$} be a MMPP with set of parameters $(\lambda_i(t), i = 1, 2, \dots, m, \Phi)$ and let $0 \le u_{i1} \le u_{i2}, i = 1, 2, \dots, m$. Then the random vector $(N_i(u_{i2}) - N_i(u_{i1}), i = 1, 2, \dots, m)$ is MTP2.

*Proof.* The distribution of $(N_i(u_{i2}) - N_i(u_{i1}), i = 1, 2, \dots, m)$ can be written as

$$P(N_i(u_{i2}) - N_i(u_{i1}) = n_i, i = 1, 2, \dots, m) = \int_{\mathbb{R}_+} \left( \prod_{i=1}^{m} g_i(n_i|\phi) \right) f_{\Phi}(\phi) d\phi,$$

for all $n_i \in \mathbb{N}$, $i = 1, 2, \dots, m$, where

$$g_i(\cdot|\phi) = P(N_i(u_{i2}) - N_i(u_{i1}) = \cdot|\Phi = \phi),$$

stands for the Poisson distribution with parameter $(\Lambda_i(u_{i2}) - \Lambda_i(u_{i1})) \phi$. As this distribution increases in the likelihood ordering with $\phi$, the result follows from Property 7.2.18 in Denuit et al. [12]. □

We recall that the MTP2 property is a strong positive dependence property, which entails for instance conditional increasingness in sequence and positive association, see e.g., Belzunce et al. [4]. Hence, such properties are fulfilled by the marginal increments in a MPGGP.

In a similar way, the concept of positive upper orthant dependent multivariate process (PUODMP) was defined in Cha and Giorgio [8]. We refer to this paper for more details. The following positive dependence result now is a direct consequence of Proposition 5.6.

**Corollary 5.7.** A MMPP is a positive upper orthant dependent multivariate process:

$$P(N_i(t_{i2}) - N_i(t_{i1}) > n_i, i = 1, 2, \dots, m) \ge \prod_{i=1}^{m} P(N_i(t_{i2}) - N_i(t_{i1}) > n_i),$$

for all $t_{i2} > t_{i1}$ and $n_i$, $i = 1, 2, \dots, m$

We now come to the extension of Proposition 5 in Cha and Mercier [9] to the multivariate setting.

**Proposition 5.8.** Let {$\mathbf{N}(t)$, $t \ge 0$} and {$\bar{\mathbf{N}}(t)$, $t \ge 0$} be two MMPPs with respective sets of parameters $(\lambda_i(t), i = 1, 2, \dots, m, \Phi)$ and $(\bar{\lambda}_i(t), i = 1, 2, \dots, m, \bar{\Phi})$. Assume that $\Phi <_{\mathbf{lr}} \bar{\Phi}$. Also, let $0 \le u_{i1} \le u_{i2}$, $i = 1, 2, \dots, m$, such that

$$\Lambda_i(u_{i2}) - \Lambda_i(u_{i1}) \le \bar{\Lambda}_i(u_{i2}) - \bar{\Lambda}_i(u_{i1}), \qquad (5.1)$$

for all $i = 1, 2, \dots, m$.

Then, we have the following result:

$$(N_i(u_{i2}) - N_i(u_{i1}), i = 1, 2, \dots, m) <_{\mathbf{lr}} (\bar{N}_i(u_{i2}) - \bar{N}_i(u_{i1}), i = 1, 2, \dots, m), \qquad (5.2)$$

for all $0 \le u_{i1} \le u_{i2}, n_i \in \mathbb{N}$, $i = 1, 2, \dots, m$, where $\mathbf{lr}$ refers to the multivariate likelihood ratio ordering.

*Proof.* We use some arguments from the proof of Theorem 3.8 in Belzunce et al. [5]. See also Theorem 2.7 in Khaledi and Shaked [19] (which is written only for absolutely continuous random variables and hence cannot be applied here).
Let us first write

$$
\begin{align}
g(n_1, \dots, n_m | \phi) &\equiv P(N_i(u_{i2}) - N_i(u_{i1}) = n_i, i = 1, 2, \dots, m | \Phi = \phi) \nonumber \\
&= \prod_{i=1}^{m} P(N_i(u_{i2}) - N_i(u_{i1}) = n_i | \Phi = \phi), \tag{5.3}
\end{align}
$$

where $P(N_i(u_{i2}) - N_i(u_{i1}) = \cdot | \Phi = \phi)$ is the Poisson distribution with parameter $(\Lambda_i(u_{i2}) - \Lambda_i(u_{i1})) \phi$.
Using similar arguments as for the proof of Proposition 5.5, it is clear that

$g(\cdots | \phi) <_{\mathbf{lr}} \bar{g}(\cdots | \phi),$

where $\bar{g}$ is defined in a similar way as (5.3) for the process $\{\bar{N}(t), t \ge 0\}$, because $\Lambda_i(u_{i2}) - \Lambda_i(u_{i1}) \le \bar{\Lambda}_i(u_{i2}) - \bar{\Lambda}_i(u_{i1})$.

As the Poisson distribution with parameter $(\Lambda_i(u_{i2}) - \Lambda_i(u_{i1})) \phi$ also increases with respect to $\phi$ in the likelihood ratio ordering, we derive that $g(n_1, \dots, n_m | \phi)$ is MTP2 in $(n_1, \dots, n_m, \phi)$.

Then, based on Theorem 2.4 in Karlin and Rinott [18], we get that

$$
\int_{\mathbb{R}_+} g(\cdots | \phi) f_{\Phi}(\phi) d\phi <_{\mathbf{lr}} \int_{\mathbb{R}_+} \bar{g}(\cdots | \bar{\phi}) f_{\bar{\Phi}}(\bar{\phi}) d\bar{\phi},
$$

which is just equivalent to (5.2) and achieves the proof. □

As a by-product of the previous proposition, considering $u_{i1} = 0$ and $u_{i2} = t$ for all $i = 1, 2, \dots, m$,
one can see that, if all parameters are fixed except from one, then $\mathbb{N}(t)$ increases in the likelihood
ordering when $k$ increases or $\Lambda(t)$ increases, and when $\alpha$ or $\nu$ decreases.

We next explore the conditions given in Proposition 5.8 to derive the comparison result on a simple
example (BPGGP).

**Example 5.9.** Let $\{\mathbb{N}(t), t \ge 0\}$ and $\{\bar{\mathbb{N}}(t), t \ge 0\}$ be two BPGGPs with sets of parameters $(\lambda_i(t) = \lambda_i, i = 1, 2, \nu, k, \alpha, l = 1)$ and $(\bar{\lambda}_i(t) = \bar{\lambda}_i, i = 1, 2, \bar{\nu}, \bar{k}, \bar{\alpha}, \bar{l} = 1)$, respectively.

For a given $t$, we know from Theorem 2.3 that the joint pdf of $(N_1(t), N_2(t))$ is given by

$$
g_{12}(x_1, x_2) = \mathbb{P}[(N_1(t) = x_1, N_2(t) = x_2)] \\
= \frac{\lambda_1^{x_1} \lambda_2^{x_2} t^{x_1+x_2}}{x_1! x_2!} \frac{\alpha^{k-\nu}}{[\alpha+t(\lambda_1+\lambda_2)]^{k+x_1+x_2-\nu}} \frac{\Gamma_\nu[k+x_1+x_2, \alpha+t(\lambda_1+\lambda_2)]}{\Gamma_\nu(k, \alpha)},
$$

for all $x_1, x_2 \in \mathbb{N}$, with a similar expression for the joint pdf of $(\bar{N}_1(t), \bar{N}_2(t))$ (denoted by $\bar{g}_{12}$).
We set

$$
\begin{align*}
G(x_1, x_2, y_1, y_2) &= g_{12} (x_1 \wedge y_1, x_2 \wedge y_2) \bar{g}_{12} (x_1 \vee y_1, x_2 \vee y_2) - g_{12} (x_1, x_2) \bar{g}_{12} (y_1, y_2), \\
\bar{G}(x_1, x_2, y_1, y_2) &= \bar{g}_{12} (x_1 \wedge y_1, x_2 \wedge y_2) g_{12} (x_1 \vee y_1, x_2 \vee y_2) - g_{12} (x_1, x_2) \bar{g}_{12} (y_1, y_2),
\end{align*}
$$

for all $x_1, x_2, y_1, y_2 \in \mathbb{N}$. Then $(N_1(t), N_2(t)) <_{\mathbf{lr}} [>_{\mathbf{lr}} (\bar{N}_1(t), \bar{N}_2(t))$ if and only if $G[\bar{G}]$ remains non-negative on $\mathbb{R}_+^4$.
<figure><img src="image_4.png" /><figcaption>Figure 1. The quotient $f/\bar{f}$ of the respective pdfs of $\Phi$ and $\bar{\Phi}$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 2. The functions $G(x_1, x_2, y_1, y_2)$ and $\bar{G}(x_1, x_2, y_1, y_2)$ with respect to $(x_2, y_2)$ for $(x_1, y_1) = (3, 10)$ and $(x_1, y_1) = (8, 2)$, respectively.</figcaption></figure>

We take $t=5$, $\lambda_1 = 1 < \bar{\lambda}_1 = 2$, $\lambda_2 = 1.5 < \bar{\lambda}_2 = 3$, $\nu=1$, $k=1$, $\alpha=1$, $\bar{\nu} = 1.25$, $\bar{k} = 0.5$, $\bar{\alpha} = 0.75$.

Then Condition (5.1) on the $\Lambda_i$'s and $\bar{\Lambda}_i$'s given in Proposition 5.8 is true. However $\Phi$ and $\bar{\Phi}$ are not comparable with respect to $\ell r$ ordering. Indeed, the quotient of their respective pdfs ($f$ and $\bar{f}$) is not monotonic, as can be seen in Figure 1.

The functions $G(x_1, x_2, y_1, y_2)$ and $\bar{G}(x_1, x_2, y_1, y_2)$ are next plotted in Figure 2 with respect to $(x_2, y_2)$ for $(x_1, y_1) = (3, 10)$ and $(x_1, y_1) = (8, 2)$, respectively. As can be seen, $G$ and $\bar{G}$ both change sign on $\mathbb{R}_+^4$ and consequently, $(N_1(t), N_2(t))$ and $(\bar{N}_1(t), \bar{N}_2(t))$ are not comparable with respect to the bivariate likelihood ratio ordering (for $t=5$).

Based on the previous example, Condition (5.1) on the $\Lambda_i$'s and $\bar{\Lambda}_i$'s is not sufficient to derive the comparison result in Proposition 5.8, and some additional comparison assumption between $\Phi$ and $\bar{\Phi}$ is required.
We finally come to the comparison between the points of two MMPPs with different parameters. To begin with, we consider the case where the two processes share the same $\lambda_i$'s, $i = 1, \dots, m$.

**Proposition 5.10.** Let $\{\mathbf{N}(t), t \ge 0\}$ and $\{\bar{\mathbf{N}}(t), t \ge 0\}$ be two MMPPs which share the same $(\lambda_i(t), i = 1, 2, \dots, m)$ with different mixture distributions $\Phi$ and $\bar{\Phi}$, respectively. Assume that $\Phi <_{\mathbf{lr}} \bar{\Phi}$. For $i = 1, 2, \dots, m$ and $n \in \mathbb{N}^*$, let $T_{in}$ (resp. $\bar{T}_{in}$) be the n-th point of $\{N_i(t), t \ge 0\}$ (resp. $\{\bar{N}_i(t), t \ge 0\}$).

Then

$$ (\bar{T}_{in_i}, i = 1, 2, \dots, m) <_{\mathbf{lr}} (T_{in_i}, i = 1, 2, \dots, m), $$

for all $n_i \in \mathbb{N}^*$ and all $i = 1, 2, \dots, m$.

*Proof.* Our aim is to use Theorem 3.8 in Belzunce et al. [5]. Let $i \in \{1, 2, \dots, m\}$ be fixed. Let us first show that $[-T_{in_i}|\Phi = \phi]$ increases with respect to $\phi$ in the likelihood ratio ordering.

For that, we will use Theorem 3.7 in Belzunce et al. [3]. Let $\phi_1 \le \phi_2$. Then, the ratio of the respective density functions of $[T_{i1}|\Phi = \phi_1]$ and $[T_{i1}|\Phi = \phi_2]$ (first points in the i-th marginal processes) is

$$ \frac{\phi_1 \lambda_i(t) e^{-\phi_1 \Lambda_i(t)}}{\phi_2 \lambda_i(t) e^{-\phi_2 \Lambda_i(t)}} = \frac{\phi_1}{\phi_2} e^{(\phi_2 - \phi_1) \Lambda_i(t)}, $$

and it is increasing. This implies that $[T_{i1}|\Phi = \phi_2] <_{\mathbf{lr}} [T_{i1}|\Phi = \phi_1]$.

Also, the ratio of the corresponding cumulative hazard rate functions is

$$ \frac{\phi_1 \Lambda_i(t)}{\phi_2 \Lambda_i(t)} = \frac{\phi_1}{\phi_2}, $$

which is constant and hence non decreasing. Based on Theorem 3.7 in Belzunce et al. [3], we derive that $[T_{in_i}|\Phi = \phi_2] <_{\mathbf{lr}} [T_{in_i}|\Phi = \phi_1]$, or equivalently that $[-T_{in_i}|\Phi = \phi_1] <_{\mathbf{lr}} [-T_{in_i}|\Phi = \phi_2]$.

Hence, $[-T_{in_i}|\Phi = \phi]$ increases with respect to $\phi$ in the likelihood ratio ordering.

Also, based on our assumptions, we know that $\Phi <_{\mathbf{lr}} \bar{\Phi}$.

The result now is a direct consequence of Theorem 3.8 in Belzunce et al. [5]. □

A natural question now is: Is it possible to compare the points in two MMPPs which share the same $\Phi$ with different $\lambda_i(t)$'s for $i = 1, 2, \dots, m$? The answer to this question is explored in next example.

**Example 5.11.** Let $\{\mathbf{N}(t), t \ge 0\}$ and $\{\bar{\mathbf{N}}(t), t \ge 0\}$ be two MPGGPs (which are specific MMPPs), which share $(v, k, \alpha, 1) = (0, 1, 1, 1)$, so that $\Phi$ is exponentially distributed with mean 1. Let $\lambda_i(t) = \lambda_i$ and $\bar{\lambda}_i(t) = \bar{\lambda}_i$, $i = 1, 2$, be the corresponding constant baseline intensity functions of the NHPPs, respectively. Then, it is easy to check that the joint density function of $(T_{11}, T_{21})$ is

$$ \begin{aligned} f_{12}(x_1, x_2) &= \int_{\mathbb{R}_+} \lambda_1 \phi e^{-\lambda_1 \phi x_1} \lambda_2 \phi e^{-\lambda_2 \phi x_2} e^{-\phi} d\phi \\ &= \frac{2\lambda_1 \lambda_2}{(\lambda_1 x_1 + \lambda_2 x_2 + 1)^3}, \end{aligned} $$

with a similar expression for the joint density function $\bar{f}_{12}$ of $(\bar{T}_{11}, \bar{T}_{21})$. The point is to see whether $(T_{11}, T_{21})$ and $(\bar{T}_{11}, \bar{T}_{21})$ are comparable with respect to the bivariate likelihood ordering. Let

$$ H(x_1, x_2, y_1, y_2) = f_{12}(x_1 \wedge y_1, x_2 \wedge y_2) \bar{f}_{12}(x_1 \vee y_1, x_2 \vee y_2) - f_{12}(x_1, x_2) \bar{f}_{12}(y_1, y_2), $$

$$ \bar{H}(x_1, x_2, y_1, y_2) = \bar{f}_{12}(x_1 \wedge y_1, x_2 \wedge y_2) f_{12}(x_1 \vee y_1, x_2 \vee y_2) - f_{12}(x_1, x_2) \bar{f}_{12}(y_1, y_2), $$
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The functions $H(x_1, x_2, y_1, y_2)$ and $\bar{H}(x_1, x_2, y_1, y_2)$ with respect to $(x_2, y_2)$ for $(x_1, y_1) = (2, 0.01)$ and $(x_1, y_1) = (0.01, 1)$, respectively.</figcaption></figure>

for all $(x_1, x_2, y_1, y_2) \in \mathbb{R}_+^4$.

Then $(T_{11}, T_{21}) <_{\text{lr}} [>_{\text{lr}}] (\bar{T}_{11}, \bar{T}_{21})$ if and only if $H[\bar{H}]$ remains non-negative on $\mathbb{R}_+^4$. The functions $H(x_1, x_2, y_1, y_2)$ and $\bar{H}(x_1, x_2, y_1, y_2)$ are plotted in Figure 3 with respect to $(x_2, y_2)$ for $(x_1, y_1) = (2, 0.01)$ and $(x_1, y_1) = (0.01, 1)$, respectively, with $\lambda_1 = \bar{\lambda}_1 = \lambda_2 = 1 < \bar{\lambda}_2 = 6$. As can be seen, $H$ and $\bar{H}$ both change sign on $\mathbb{R}_+^4$ and consequently, $(T_{11}, T_{21})$ and $(\bar{T}_{11}, \bar{T}_{21})$ are not comparable with respect to the bivariate likelihood ratio ordering.

Based on this simple example (with constant $\lambda_i$'s and $\bar{\lambda}_i$'s such that $\lambda_1 = \bar{\lambda}_1 = \lambda_2 < \bar{\lambda}_2$), it seems that there is no hope to find conditions under which the points in two MMPPs with different $\lambda_i(t)$'s could be comparable with respect to the multivariate likelihood ratio ordering.

However, it is possible to get comparison results with respect to the weaker usual stochastic ordering. We recall that given two random vectors $\mathbf{X}$ and $\mathbf{Y}$ on $\mathbb{R}^n$, then $\mathbf{X}$ is said to be smaller than $\mathbf{Y}$ in the usual stochastic ordering (written $\mathbf{X} <_{\text{sto}} \mathbf{Y}$) as soon as

$$E[\varphi(\mathbf{X})] \le E[\varphi(\mathbf{Y})],$$

for all non-decreasing function $\varphi: \mathbb{R}^n \to \mathbb{R}$ such that the expectations exist. In the univariate setting (written $X <_{\text{sto}} Y$), it is equivalent to $\bar{F}_X(t) \le \bar{F}_Y(t)$, for all $t \ge 0$.

The multivariate likelihood ratio ordering is known to imply the usual stochastic ordering. See Shaked and Shanthikumar [24] for more details.

We now come to the comparison result.

**Proposition 5.12.** Let $\{\mathbb{N}(t), t \ge 0\}$ and $\{\bar{\mathbb{N}}(t), t \ge 0\}$ be two MMPPs with sets of parameters $(\lambda_i(t), i = 1, 2, \dots, m, \Phi)$ and $(\bar{\lambda}_i(t), i = 1, 2, \dots, m, \bar{\Phi})$, respectively. Assume that $\Phi <_{\text{lr}} \bar{\Phi}$ and $\Lambda_i(t) \le \bar{\Lambda}_i(t)$, for all $t \ge 0$ and all $i = 1, 2, \dots, m$. Using the notations of Proposition 5.10, we have:

$$(\bar{T}_{in_i}, i = 1, 2, \dots, m) <_{\text{sto}} (T_{in_i}, i = 1, 2, \dots, m),$$

for all $n_i \in \mathbb{N}^*$ and all $i = 1, 2, \dots, m$.
Proof. Our aim is to use Theorem 3.1 in Belzunce et al. [5]. We already know from the proof of Proposition 5.10 that $[-T_{in_i}|\Phi = \phi]$ increases with respect to $\phi$ in the likelihood ratio ordering, and hence also in the usual stochastic ordering.

Also, based on $\Lambda_i(t) \le \bar{\Lambda}_i(t)$, for all $t \ge 0$ and all $i = 1, 2, \dots, m$, it is easy to check that the conditional survival functions of $\bar{T}_{i1}$ given $\bar{\Phi} = \phi$ and $T_{i1}$ given $\Phi = \phi$ fulfill

$$ \bar{F}_{\bar{T}_{i1}|\bar{\Phi}=\phi}(t) \equiv \mathbb{P}(\bar{T}_{i1} > t | \bar{\Phi} = \phi) = e^{-\phi \bar{\Lambda}_i(t)} \le e^{-\phi \Lambda_i(t)} = \bar{F}_{T_{i1}|\bar{\Phi}=\phi}(t), $$

for all $t \ge 0$, which means that $[\bar{T}_{i1}|\bar{\Phi} = \phi] <_{sto} [T_{i1}|\Phi = \phi]$. We derive from Theorem 3.1 in Belzunce et al. [3] that $[(\bar{T}_{i1}, \bar{T}_{i2}, \dots, \bar{T}_{in_i}) | \bar{\Phi} = \phi] <_{sto} [(T_{i1}, T_{i2}, \dots, T_{in_i}) | \Phi = \phi]$ and next that $[\bar{T}_{in_i}|\bar{\Phi} = \phi] <_{sto} [T_{in_i}|\Phi = \phi]$ (as the usual stochastic ordering is stable through marginalization), or equivalently that $[-T_{in_i}|\Phi = \phi] <_{sto} [\bar{T}_{in_i}|\bar{\Phi} = \phi]$.

Finally, based on $\Phi <_{lr} \bar{\Phi}$, we can derive that

$$ (-T_{in_i}, i = 1, 2, \dots, m) <_{sto} (-\bar{T}_{in_i}, i = 1, 2, \dots, m), $$

from Theorem 3.1 in Belzunce et al. [5], which allows to conclude. □

**Remark 5.13.** Note that all the results from Propositions 5.5, 5.6 and Corollary 5.7 hold for MPGGPs, as they are specific MMPPs. In order to apply Propositions 5.8, 5.10 and 5.12 for MPGGPs, one can use the following result which provides conditions under which $\Phi <_{lr} \bar{\Phi}$. The arguments are given in the proof of Proposition in Cha and Mercier [9].

**Lemma 5.14.** Let $\Phi \sim \mathcal{GG}(v, k, \alpha, l = 1)$ and $\bar{\Phi} \sim \mathcal{GG}(\bar{v}, \bar{k}, \bar{\alpha}, \bar{l} = 1)$. Then $\Phi <_{lr} \bar{\Phi}$ as soon as one of the following conditions holds:

* $\bar{\alpha} = \alpha, \bar{k} \ge k$ and $\bar{k} - k \ge \bar{v} - v$;

* $\bar{\alpha} < \alpha$ and $(\alpha - \bar{\alpha} + \bar{k} - k + v - \bar{v})^2 - 4(\alpha - \bar{\alpha})(\bar{k} - k) \le 0$;

* $\bar{\alpha} < \alpha$ and $\alpha - \bar{\alpha} + \bar{k} - k + v - \bar{v} \ge 0$.

As a specific case, we can see that, if all parameters are fixed except from one, $\Phi$ increases in the likelihood ordering when $k$ increases, and when $\alpha$ or $v$ decreases.

**Acknowledgment.** The authors thank the reviewers for helpful comments and suggestions.

**Funding statement.** This work was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (Grant Number: 2019R1A6A1A11051177).

**Competing interest.** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## References

[1] Agarwal S.K. & Kalla S.L. (1996). A generalized gamma distribution and its application in reliability. Communications in Statistics – Theory and Methods 25(1): 201–210.

[2] Allen F. & Gale D. (2000). Financial contagion. Journal of Political Economy 108(1): 1–34.

[3] Belzunce F., Lillo R.E., Ruiz J.M. & Shaked M. (2001). Stochastic comparisons of nonhomogeneous processes. Probability in the Engineering and Informational Sciences 15(2): 199–224.

[4] Belzunce F., Martínez-Riquelme C. & Mulero J. (2016). An Introduction to Stochastic Orders, Amsterdam: Elsevier/Academic Press.

[5] Belzunce F., Mercader J.A., Ruiz J.M. & Spizzichino F. (2009). Stochastic comparisons of multivariate mixture models. Journal of Multivariate Analysis 100(8): 1657–1669.

[6] Bowsher C.G. (2006). Modelling security market events in continuous time: intensity based, multivariate point process models. Journal of Econometrics 141(2): 876–912.
[7] Cha J.H. (2014). Characterization of the generalized Pólya process and its applications. *Advances in Applied Probability* 46(4): 1148–1171.

[8] Cha J.H. & Giorgio M. (2016). On a class of multivariate counting processes. *Advances in Applied Probability* 48(2): 443–462.

[9] Cha J.H. & Mercier S. (2021). Poisson generalized Gamma process and its properties. *Stochastics* 93(8): 1123–1140.

[10] Cox D.R. & Lewis P.A.W. (1972). Multivariate point processes. In *Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability*, pp. 401–448.

[11] Daley D.J. (1968). The correlation structure of the output process of some single server queueing systems. *Annals of Mathematical Statistics* 39(3): 1007–1019.

[12] Denuit M., Dhaene J., Goovaerts M. & Kaas R. (2006). *Actuarial Theory for Dependent Risks: Measures, Orders and Models*. Chichester: John Wiley & Sons.

[13] Finkelstein M. (2004). Minimal repair in heterogeneous populations. *Journal of Applied Probability* 41(1): 281–286.

[14] Finkelstein M. (2008). *Failure Rate Modeling for Reliability and Risk*. London: Springer.

[15] Ghitany M.E. (1998). On a recent generalization of gamma distribution. *Communications in Statistics – Theory and Methods* 27(1): 223–233.

[16] Grandell J. (1997). *Mixed Poisson Processes*. Monographs on Statistics and Applied Probability, Vol. 77. London: Chapman & Hall.

[17] Gupta R.C. & Ong S.H. (2004). A new generalization of the negative binomial distribution. *Computational Statistics & Data analysis* 45(4): 287–300.

[18] Karlin S. & Rinott Y. (1980). Classes of orderings of measures and related correlation inequalities. I. Multivariate totally positive distributions. *Journal of Multivariate Analysis* 10(4): 467–498.

[19] Khaledi B.E. & Shaked M. (2010). Stochastic comparisons of multivariate mixtures. *Journal of Multivariate Analysis* 101(10): 2486–2498.

[20] Kobayashi K. (1991). On generalized gamma functions occurring in diffraction theory. *Journal of the Physical Society of Japan* 60(5): 1501–1512.

[21] Partrat C. (1994). Compound model for two dependent kinds of claims. *Insurance: Mathematics and Economics* 15(2): 219–231.

[22] Riordan J. (1937). Moment recurrence relations for binomial, Poisson and hypergeometric frequency distributions. *Annals of Mathematical Statistics* 8(2): 103–111.

[23] Ross S.M. (2003). *Introduction to Probability Models, 8th ed*. San Diego: Academic Press.

[24] Shaked M. & Shanthikumar J.G. (2007). *Stochastic Orders*. New York: Springer.

<footnote>Cite this article: Cha JH. and Mercier S. (2024). An extended class of multivariate counting processes and its main properties. *Probability in the Engineering and Informational Sciences* 1–18. https://doi.org/10.1017/S0269964824000111</footnote>