2. Background

The competing risks model is often employed in reliability theory and survival analysis, being appropriate to describe the failures of devices or organisms in the presence of different types of risks. Hereafter we describe the basic notions that will be used in the rest of the paper.

Let T be an absolutely continuous non-negative random variable which describes the random failure time of a system subject to two possible causes of failure. We set $\delta=i$ if the failure occurs due to a shock of type $i,\,i=1,2$ . Let N(t) denote the total number of shocks occurring in [0, t], with $t\geq 0$ . Moreover, let $N_i(t)$ denote the counting process representing the number of shocks of type i occurring in [0, t], $i=1,2$ , $N_{1}(t)$ and $N_{2}(t)$ being independent. Clearly we have

\[N(t) = N_1(t) + N_2(t), \quad t\geq 0.\]

Moreover, we assume that the system fails when the sum of shocks of type 1 and type 2 reaches a random threshold M that takes values in $\mathbb{N}=\{1,2,\ldots\}$ . In other words, failures are assumed to occur, due to a single cause, at the first instant in which $N_1(t)+N_2(t)=M$ . The probability distribution and the survival probability of M will be respectively denoted by

(2.1) \begin{equation}p_k=\mathbb{P}(M=k), \quad k\in\mathbb{N},\end{equation}

and

(2.2) \begin{equation}\overline P_k=\mathbb{P}(M>k), \quad k\in\mathbb{N}_{0}=\mathbb{N}\cup\{0\}.\end{equation}

Letting $f_T(t)$ , $t\geq 0$ , denote the probability density function of the failure time (first-hitting time)

\[T=\inf\{t\geq 0\,:\, N(t)=M\},\]

we have

(2.3) \begin{equation}f_T(t)=f_1(t)+f_2(t), \quad t\geq 0,\end{equation}

where $f_i(t)$ is the sub-density defined by

\[f_i(t)=\dfrac{\rm d}{{\rm d}t}\mathbb{P}\{T\leq t,\,\delta=i\},\quad t\geq 0, \ i=1,2.\]

Furthermore, the probability that the failure occurs due to a shock of type i reads

(2.4) \begin{equation}\mathbb{P}(\delta=i)=\int_0^{\infty} f_i(t)\,{\rm d}t,\quad i=1,2.\end{equation}

In order to express the sub-densities $f_i(t)$ , $i=1,2$ , in terms of the joint probability distribution of $(N_1(t),N_2(t))$ , we now introduce the hazard rates

(2.5a) \begin{align}r_1(x_1,x_2;\, t)=\displaystyle\lim_{\tau\to 0^+} \dfrac{1}{\tau}\mathbb{P}\big\{N_1(t+\tau)= x_1+1,\, N_2(t+\tau)= x_2| N_1(t)=x_1,\, N_2(t)=x_2\big\},\end{align}

(2.5b) \begin{align}r_2(x_1,x_2;\, t)=\displaystyle\lim_{\tau\to 0^+} \dfrac{1}{\tau}\mathbb{P}\big\{N_1(t+\tau)= x_1,\, N_2(t+\tau)= x_2+1| N_1(t)=x_1,\, N_2(t)=x_2\big\},\end{align}

with $(x_1,x_2)\in\mathbb{N}^2_0$ and $t\geq 0$ . Given that $x_1$ shocks of type 1 and $x_2$ shocks of type 2 occurred in [0, t], the hazard rate $r_i(x_1,x_2;\, t)$ gives the intensity of the occurrence of a shock of type i immediately after t, with $i=1,2$ . Hence, conditioning on M and recalling (2.1) and (2.5), for $t\geq 0$ and $i=1,2$ , the failure densities can be expressed as

(2.6) \begin{equation}f_i(t)=\sum_{k=1}^{+\infty} p_k \sum_{x_1+x_2=k-1}\mathbb{P}\{N_1(t)=x_1,N_2(t)=x_2\}\,r_i(x_1,x_2;\,t).\end{equation}

A relation similar to (2.6) holds for the survival function of T, denoted by

\[\overline{F}_T(t)=\mathbb{P}\{T>t\}, \quad t\geq 0.\]

Indeed, conditioning on $(N_1(t),N_2(t))$ and recalling (2.2), we obtain

(2.7) \begin{equation}\overline{F}_T(t)=\sum_{k=0}^{+\infty} \overline P_k \sum_{x_1+x_2=k}\mathbb{P}\{N_1(t)=x_1,N_2(t)=x_2\},\quad t\geq 0,\end{equation}

where $\overline P_0=1$ .

3. The model and main results

In this section we introduce the specific stochastic model that will be adopted for the bivariate counting process of interest. Then we extensively study the quantities of interest for the competing risks model under the introduced scheme. Let $\{\mathcal{A}^{\nu}(t)\,:\, t\geq 0\}$ be the stable subordinator, i.e. the non-decreasing Lévy process with Bernštein function $f(u)=u^{\nu},\,\nu\in(0,1)$ . Its Laplace transform is therefore

(3.1) \begin{equation}\mathbb{E}\bigl({\mathrm{e}}^{-u\mathcal{A}^{\nu}(t)}\bigr)={\mathrm{e}}^{-u^{\nu}t},\quad u\geq 0\end{equation}

(see [Reference Applebaum2, Example 1.3.18]). We also recall that the Lévy measure associated with f is

\begin{equation*}\nu({\mathrm{d}} s)=\dfrac{\nu s^{-\nu-1}}{\Gamma(1-\nu)}\,{\mathrm{d}} s.\end{equation*}

We denote the density of $\mathcal{A}^{\nu}(t)$ by $f_{\mathcal{A}^{\nu}(t)}$ , that is (see [Reference Uchaikin and Zolotarev23]),

\begin{equation*}f_{\mathcal{A}^{\nu}(t)}(y,t)\,{\mathrm{d}} y=\mathbb{P}(\mathcal{A}^{\nu}(t)\in{\mathrm{d}} y)=\dfrac{1}{\pi}\sum_{k=1}^{\infty} ({-}1)^{k+1} \dfrac{\Gamma(\nu k +1)}{k!}\, \dfrac{t^k}{y^{\nu k+1} } \sin (\pi \nu k)\,{\mathrm{d}} y,\quad y>0.\end{equation*}

Furthermore, the joint density function, defined as

\begin{equation*}f_{\mathcal{A}^{\nu}(t_{1}),\,\mathcal{A}^{\nu}(t_{2})}(y_{1},t_{1};\,y_{2},t_{2})\,{\mathrm{d}} y_{1}\,{\mathrm{d}} y_{2}=\mathbb{P}(\mathcal{A}^{\nu}(t_{1})\in{\mathrm{d}} y_{1},\,\mathcal{A}^{\nu}(t_{2})\in{\mathrm{d}} y_{2}),\end{equation*}

is given by

(3.2) \begin{equation}f_{\mathcal{A}^{\nu}(t_{1}),\,\mathcal{A}^{\nu}(t_{2})}(y_{1},t_{1};\,y_{2},t_{2})\,{\mathrm{d}} y_{1}\,{\mathrm{d}} y_{2}=f_{\mathcal{A}^{\nu}(t_{2}-t_{1})}(y_{2}-y_{1},t_{2}-t_{1})f_{\mathcal{A}^{\nu}(t_{1})}(y_{1},t_{1})\,{\mathrm{d}} y_{1}\,{\mathrm{d}} y_{2},\end{equation}

since the process $\mathcal{A}^{\nu}(t)$ has independent and stationary increments. We assume that the two kinds of shocks affect the system according to a bivariate space-fractional Poisson process, which is a particular case of the more general multivariate space–time fractional Poisson process considered in [Reference Beghin and Macci4]. It is the time-change of a bivariate vector of independent classical homogeneous Poisson processes where the time change is an independent stable subordinator, defined hereafter.

Definition 3.1. Let $\{\{{\mathcal P}_{i}(t)\,:\, t\geq 0\}\,:\, i\in\{1,2\}\}$ be two independent homogeneous Poisson processes with intensities $\lambda_{1}$ and $\lambda_{2}$ respectively. Then, for $\nu\in(0,1]$ , we consider the bivariate process

\begin{equation*}\bigl(N_{1}^{\nu}(t),N_{2}^{\nu}(t)\bigr)\,:\!=\, \big({\mathcal P}_{1}(\mathcal{A}^{\nu}(t)),{\mathcal P}_{2}(\mathcal{A}^{\nu}(t))\big),\quad t\geq 0,\end{equation*}

where ${\mathcal P}_{i}(t),\,i=1,2$ , and $\mathcal{A}^{\nu}(t)$ are independent.

We remark that the processes $\bigl\{\bigl\{N_{i}^{\nu}(t)\,:\, t\geq 0\bigr\}\,:\, i\in\{1,2\}\bigr\}$ are conditionally independent given $\mathcal{A}^{\nu}(t)$ for $0<\nu<1$ . Unlike the classical Poisson process, they are non-renewal processes, as noted in [Reference Garra, Orsingher and Scavino12] and [Reference Polito and Scalas22]. In addition, they are point processes with stationary independent increments generalizing the Poisson process, in that they perform non-unitary jumps with non-negligible probability. We refer the reader to [Reference Orsingher and Polito20] and [Reference Orsingher and Toaldo21] for more details. Note that for $\nu=1$ we have $\mathcal{A}^{1}(t)=t$ , so they are independent in this case. The state probabilities of the process in Definition 3.1 read, for $x_{1},\,x_{2}\in\mathbb{N}_{0}$ (see [Reference Beghin and Macci4, equation (6)]),

(3.3) \begin{align}\mathbb{P}\bigl\{N_{1}^{\nu}(t)=x_{1},\,N_{2}^{\nu}(t)=x_{2}\bigr\}&=\dfrac{\lambda_{1}^{x_{1}}\lambda_{2}^{x_{2}}}{(\lambda_{1}+\lambda_{2})^{x_{1}+x_{2}}}\dfrac{({-}1)^{x_{1}+x_{2}}}{x_{1}!x_{2}!}\nonumber\\&\quad \times{}_1 \Psi_1\biggl[\begin{matrix}(1,\nu)\\(1-(x_{1}+x_{2}),\nu)\end{matrix} \biggm| ({-}(\lambda_{1}+\lambda_{2})^{\nu}t)\biggr],\end{align}

where

\begin{equation*}{}_p \Psi_q\biggl[\begin{matrix}(a_{l},\alpha_{l})_{1,p}\\[4pt](b_{l},\beta_{l})_{1,q}\end{matrix} \biggm| z\biggr]\,:\!=\, \sum_{k=0}^{+\infty}\dfrac{\prod_{i=1}^{p}\Gamma(a_{l}+\alpha_{l}k)}{\prod_{j=1}^{q}\Gamma\big(b_{j}+\beta_{j}k\big)}\dfrac{z^{k}}{k!}\end{equation*}

is the generalized Wright function

\[\big(z\in\mathbb{C}, a_{i},b_{j}\in\mathbb{C}, \alpha_{i},\beta_{j}\in\mathbb{R}, \alpha_{i},\beta_{j}\neq 0;\, i=1,2,\ldots,p;\,j=1,2,\ldots,q\big),\]

under the convergence condition

\begin{equation*}\sum_{k=1}^{q}\beta_{k}-\sum_{h=1}^{p}\alpha_{h}>-1.\end{equation*}

For details see e.g. [Reference Kilbas, Srivastava and Trujillo16, equations (1.11.14), (1.11.15)]. We observe that the joint distribution of the process $\bigl(N_{1}^{\nu}(t),N_{2}^{\nu}(t)\bigr)$ is exchangeable. Furthermore, by resorting to the expression of the probability generating function given in equation (5) of [Reference Beghin and Macci4], one can easily check that the expected value of each of the two components of the process is infinite, thus making our model suitable to describe bursty dynamics. In order to evaluate the hazard rates, we state the following lemma.

Lemma 3.1. Fix $m\in\mathbb{N}_{0}$ . Then

(3.4) \begin{equation}\dfrac{{\mathrm{d}} ^{m} }{{\mathrm{d}} u^{m}}\bigl( {\mathrm{e}}^{-u^{\nu}t} \bigr)=u^{-m}\,{\mathrm{e}}^{-u^{\nu}t}\sum_{k=0}^{m}\dfrac{(u^{\nu}t)^{k}}{k!}\sum_{j=0}^{k}({-}1)^{j}\binom{k}{j}(\nu j)_{m},\end{equation}

where $(x)_{m}=x(x-1)\cdots(x-m+1)$ denotes the falling factorial.

Proof. We recall that if a and b are functions with a sufficient number of derivatives, then Hoppe’s formula for the m-fold derivative of a composition of functions [Reference Johnson15] states that

\begin{equation*}\dfrac{{\mathrm{d}} ^{m}}{{\mathrm{d}} t^{m}}a(b(t))=\sum_{k=0}^{m}\dfrac{1}{k!}\dfrac{{\mathrm{d}} ^{k}}{{\mathrm{d}} b^{k}}a(b)A_{m,k}(b(t)),\end{equation*}

where

\begin{equation*}A_{m,k}(b(t))=\sum_{j=0}^{k}\binom{k}{j}({-}b(t))^{k-j}\dfrac{{\mathrm{d}} ^{m} }{{\mathrm{d}} t^{m}}(b(t))^{j}.\end{equation*}

For $b(t)\,:\!=\, -u^{\nu}t$ and $a(b)\,:\!=\, {\mathrm{e}}^{b}$ , we have

\begin{equation*}\dfrac{{\mathrm{d}} ^{m}}{{\mathrm{d}} t^{m}}(b(t))^{j}=({-}t)^{j}(\nu j)_{m}u^{\nu j-m},\quad 0\leq j\leq m.\end{equation*}

Therefore

\begin{equation*}\dfrac{{\mathrm{d}} ^{m} }{{\mathrm{d}} u^{m}}\bigl( {\mathrm{e}}^{-u^{\nu}t} \bigr)=\sum_{k=0}^{m}\dfrac{{\mathrm{e}}^{-u^{\nu}t}}{k!}\sum_{j=0}^{k}\binom{k}{j}(u^{\nu}t)^{k-j}({-}t)^{j}(\nu j)_{m}u^{\nu j-m},\end{equation*}

and, after simplifying, we obtain (3.4).

We now derive the expression of the hazard rates (2.5) of the process introduced in Definition 3.1.

Proposition 3.1. Let $i\in\{1,2\}$ . Under the assumptions of the model in Definition 3.1, for all $x_{1},\,x_{2}\in\mathbb{N}_{0}$ and $t\geq 0$ , we have

(3.5) \begin{align}r_i(x_1,x_2;\, t)&=\biggl({}_1 \Psi_1\biggl[\begin{matrix}(1,\nu)\\(1-(x_{1}+x_{2}),\nu)\end{matrix} \biggm| ({-}(\lambda_{1}+\lambda_{2})^{\nu}t)\biggr]\biggr)^{-1}\notag \\&\quad \times \lambda_{i}\nu(\lambda_{1}+\lambda_{2})^{\nu-1}\,{\mathrm{e}}^{-(\lambda_{1}+\lambda_{2})^{\nu}t}\sum_{k=0}^{x_{1}+x_{2}}\dfrac{[(\lambda_{1}+\lambda_{2})^{\nu}t]^{k}}{k!}\sum_{j=0}^{k}({-}1)^{j}\binom{k}{j}(\nu j)_{x_{1}+x_{2}}.\end{align}

Proof. Fix $i=1$ and let $t\geq 0$ . To begin with, we observe that (2.5a) can be equivalently rewritten as

(3.6) \begin{equation}\displaystyle\lim_{\tau\to t} \dfrac{\mathbb{P}\bigl\{N_1^{\nu}(\tau)= x_1+1,\, N_2^{\nu}(\tau)= x_2,\, N_1^{\nu}(t)=x_1,\, N_2^{\nu}(t)=x_2\bigr\}}{(\tau-t)\mathbb{P}\bigl\{N_1^{\nu}(t)=x_1,\, N_2^{\nu}(t)=x_2\bigr\}}.\end{equation}

We concentrate our attention on the numerator. From Definition 3.1 we know that the governing fractional bivariate process can be considered as a two-dimensional vector of homogeneous Poisson processes stopped at a random time $\mathcal{A}^{\nu}(t)$ . Therefore we have

(3.7) \begin{align}&\mathbb{P}\bigl\{N_{1}^{\nu}(\tau)=x_{1}+1,\,N_{2}^{\nu}(\tau)=x_{2},\,N_{1}^{\nu}(t)=x_{1},\,N_{2}^{\nu}(t)=x_{2}\bigr\}\nonumber\\&\quad =\int_{0}^{+\infty}\int_{0}^{v}\mathbb{P}\big\{N_{1}(v)=x_{1}+1,\,N_{2}(v)=x_{2},\,N_{1}(u)=x_{1},\,N_{2}(u)=x_{2}\big\}\nonumber\\&\quad \quad{}\times \mathbb{P}(\mathcal{A}^{\nu}(\tau)\in{\mathrm{d}} v,\mathcal{A}^{\nu}(t)\in{\mathrm{d}} u)\nonumber\\&\quad \overset{\text{by (3.2)}}{=}\int_{0}^{+\infty}\int_{0}^{v}\mathbb{P}\big\{\mathcal{P}_{1}(v)=x_{1}+1,\,\mathcal{P}_{2}(v)=x_{2},\,\mathcal{P}_{1}(u)=x_{1},\,\mathcal{P}_{2}(u)=x_{2}\big\}\nonumber\\&\quad \quad{}\times f_{A^{\nu}(\tau-t)}(v-u,\tau-t)f_{A^{\nu}(t)}(u,t)\,{\mathrm{d}} v\,{\mathrm{d}} u.\end{align}

The two homogeneous Poisson processes $\mathcal{P}_{1}(t)$ and $\mathcal{P}_{2}(t)$ are independent, and hence

(3.8) \begin{align}&\mathbb{P}\big\{\mathcal{P}_{1}(v)=x_{1}+1,\,\mathcal{P}_{2}(v)=x_{2},\,\mathcal{P}_{1}(u)=x_{1},\,\mathcal{P}_{2}(u)=x_{2}\big\}\nonumber\\&\quad =\mathbb{P}\big\{\mathcal{P}_{1}(v-u)=1,\,\mathcal{P}_{2}(v-u)=0\big\}\mathbb{P}\big\{\mathcal{P}_{1}(u)=x_{1},\,\mathcal{P}_{2}(u)=x_{2}\big\}\notag \\&\quad =\dfrac{{\mathrm{e}}^{-(\lambda_{1}+\lambda_{2})v}\lambda_{1}^{x_{1}+1}\lambda_{2}^{x_{2}}u^{x_{1}+x_{2}}(v-u)}{x_{1}!x_{2}!}.\end{align}

We put (3.8) into (3.7) and get

\begin{align*}&\mathbb{P}\bigl\{N_{1}^{\nu}(\tau)=x_{1}+1,\,N_{2}^{\nu}(\tau)=x_{2},\,N_{1}^{\nu}(t)=x_{1},\,N_{2}^{\nu}(t)=x_{2}\bigr\}\\&\quad =\dfrac{\lambda_{1}^{x_{1}+1}\lambda_{2}^{x_{2}}}{x_{1}!x_{2}!}\int_{0}^{+\infty}\,{\mathrm{e}}^{-(\lambda_{1}+\lambda_{2})v}\,{\mathrm{d}} v\int_{0}^{v}u^{x_{1}+x_{2}}(v-u)f_{A^{\nu}(\tau-t)}(v-u,\tau-t)f_{A^{\nu}(t)}(u,t)\,{\mathrm{d}} u\\&\quad =\dfrac{\lambda_{1}^{x_{1}+1}\lambda_{2}^{x_{2}}}{x_{1}!x_{2}!}\int_{0}^{+\infty}u^{x_{1}+x_{2}}f_{A^{\nu}(t)}(u,t)\,{\mathrm{d}} u\int_{u}^{+\infty}\,{\mathrm{e}}^{-(\lambda_{1}+\lambda_{2})v}(v-u)f_{A^{\nu}(\tau-t)}(v-u,\tau-t)\,{\mathrm{d}} v,\end{align*}

by changing the order of integration. By the change of variable $v-u=y$ , we obtain

(3.9) \begin{align}&\mathbb{P}\bigl\{N_{1}^{\nu}(\tau)=x_{1}+1,\,N_{2}^{\nu}(\tau)=x_{2},\,N_{1}^{\nu}(t)=x_{1},\,N_{2}^{\nu}(t)=x_{2}\bigr\}\nonumber\\&\quad =\dfrac{\lambda_{1}^{x_{1}+1}\lambda_{2}^{x_{2}}}{x_{1}!x_{2}!}\,\mathcal{I}(x_{1}+x_{2},t)\cdot \mathcal{I}(1,\tau-t),\end{align}

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

\begin{equation*}\mathcal{I}(h,t)\,:\!=\, \int_{0}^{+\infty}\,{\mathrm{e}}^{-(\lambda_{1}+\lambda_{2})u}u^{h}f_{A^{\nu}(t)}(u,t)\,{\mathrm{d}} u.\end{equation*}

It turns out that, due to (3.1) and (3.4),

\begin{align*}\mathcal{I}(h,t)&=({-}1)^{h}\dfrac{{\mathrm{d}} ^{h}}{{\mathrm{d}} u^{h}}\mathbb{E}\big({\mathrm{e}}^{-u\mathcal{A}^{\nu}(t)}\big)\biggr\rvert_{u=\lambda_{1}+\lambda_{2}}\nonumber\\&=({-}1)^{h}\dfrac{{\mathrm{d}} ^{h} }{{\mathrm{d}} u^{h}} \,{\mathrm{e}}^{-u^{\nu}t}\biggr\rvert_{u=\lambda_{1}+\lambda_{2}}\nonumber\\&=\dfrac{({-}1)^{h}}{(\lambda_{1}+\lambda_{2})^{h}}\,{\mathrm{e}}^{-(\lambda_{1}+\lambda_{2})^{\nu}t}\sum_{k=0}^{h}\dfrac{[(\lambda_{1}+\lambda_{2})^{\nu}t]^{h}}{h!}\sum_{j=0}^{k}({-}1)^{j}\binom{k}{j}(\nu j)_{h}.\end{align*}

By exploiting the previous expression, we evaluate (3.9) as follows:

(3.10) \begin{align}&\mathbb{P}\bigl\{N_{1}^{\nu}(z)=x_{1}+1,\,N_{2}^{\nu}(z)=x_{2},\,N_{1}^{\nu}(t)=x_{1},\,N_{2}^{\nu}(t)=x_{2}\bigr\}\nonumber\\&\quad =\dfrac{\lambda_{1}^{x_{1}+1}\lambda_{2}^{x_{2}}}{x_{1}!x_{2}!}\nu(z-t)(\lambda_{1}+\lambda_{2})^{\nu-1-x_{1}-x_{2}}\,{\mathrm{e}}^{-(\lambda_{1}+\lambda_{2})^{\nu}z}({-}1)^{x_{1}+x_{2}} \notag \\&\quad\quad \times\sum_{k=0}^{x_{1}+x_{2}}\dfrac{[(\lambda_{1}+\lambda_{2})^{\nu}t]^{k}}{k!}\sum_{j=0}^{k}({-}1)^{j}\binom{k}{j}(\nu j)_{x_{1}+x_{2}}.\end{align}

The hazard rate (3.6) can now be easily computed by taking into account (3.3) and (3.10), thus getting the desired result. The case $i=2$ can be treated analogously, and this concludes the proof.

Remark 3.1. We observe that (3.5) with $\nu=1$ meets a known formula for the non-fractional case (see [Reference Di Crescenzo and Longobardi9, Section 3]), i.e. $r_{i}(x_{1},x_{2};\,t)=\lambda_{i},\,i\in\{1,2\}$ . Indeed, since $(j)_{x_{1}+x_{2}}=0$ when $j<x_{1}+x_{2}$ , the sum

\begin{equation*}\sum_{k=0}^{x_{1}+x_{2}}\dfrac{[(\lambda_{1}+\lambda_{2})t]^{k}}{k!}\sum_{j=0}^{k}({-}1)^{j}\binom{k}{j}( j)_{x_{1}+x_{2}}\end{equation*}

reduces to

\begin{equation*}\dfrac{[(\lambda_{1}+\lambda_{2})t]^{x_{1}+x_{2}}}{(x_{1}+x_{2})!}({-}1)^{x_{1}+x_{2}}(x_{1}+x_{2})_{x_{1}+x_{2}}=[{-}(\lambda_{1}+\lambda_{2})t]^{x_{1}+x_{2}}.\end{equation*}

Moreover, the Fox–Wright function ${}_1 \Psi_1$ simplifies to

\begin{align*}{}_1 \Psi_1\biggl[\begin{matrix}(1,1)\\[3pt](1-(x_{1}+x_{2}),1)\end{matrix} \biggm| -(\lambda_{1}+\lambda_{2})t\biggr]&=\sum_{h=0}^{+\infty}\dfrac{\Gamma(h+1)}{\Gamma(1-x_{1}-x_{2}+h)}\dfrac{[{-}(\lambda_{1}+\lambda_{2})t]^{h}}{h!}\\&=\sum_{h=x_{1}+x_{2}}^{+\infty}\dfrac{[{-}(\lambda_{1}+\lambda_{2})t]^{h}}{(h-x_{1}-x_{2})!}\\&=[{-}(\lambda_{1}+\lambda_{2})t]^{x_{1}+x_{2}}\,{\mathrm{e}}^{-(\lambda_{1}+\lambda_{2})t}.\end{align*}

The second equality holds because the summands with $h<x_{1}+x_{2}$ are equal to 0. In the light of this, after some manipulations, the hazard rate (3.5) reduces to $\lambda_{i}$ , i.e. the hazard rate in the classical case.

Hereafter we present the formulas for the failure densities (2.6) and for the survival function of T given in (2.7).

Proof. (i) We substitute in (2.6) the expression of the state probability (3.3) and of the hazard rate (3.5), thus getting, after some computations,

\begin{align*}f_{i}(t)&=\sum_{k=1}^{+\infty}p_{k}\sum_{x_{1}+x_{2}=k-1}\dfrac{\lambda_{1}^{x_{1}}\lambda_{2}^{x_{2}}({-}1)^{x_{1}+x_{2}}}{(\lambda_{1}+\lambda_{2})^{x_{1}+x_{2}}x_{1}!x_{2}!}\lambda_{i}\nu(\lambda_{1}+\lambda_{2})^{\nu-1}\,{\mathrm{e}}^{-(\lambda_{1}+\lambda_{2})^{\nu}t}\\&\quad{}\times\sum_{h=0}^{x_{1}+x_{2}}\dfrac{[(\lambda_{1}+\lambda_{2})^{\nu}t]^{h}}{h!}\sum_{j=0}^{h}({-}1)^{j}\binom{h}{j}(\nu j)_{x_{1}+x_{2}}\\&=\lambda_{i}\nu(\lambda_{1}+\lambda_{2})^{\nu-1}\,{\mathrm{e}}^{-(\lambda_{1}+\lambda_{2})^{\nu}t}\sum_{k=1}^{+\infty}\dfrac{p_{k}}{(k-1)!}\dfrac{({-}1)^{k-1}}{(\lambda_{1}+\lambda_{2})^{k-1}}\\&\quad{}\times\sum_{x_{1}=0}^{k-1}\dfrac{(k-1)!}{x_{1}!(k-1-x_{1})!}\lambda_{1}^{x_{1}}\lambda_{2}^{k-1-x_{1}}\sum_{h=0}^{k-1}\dfrac{[(\lambda_{1}+\lambda_{2})^{\nu}t]^{h}}{h!}\sum_{j=0}^{h}({-}1)^{j}\binom{h}{j}(\nu j)_{k-1}.\end{align*}

The thesis follows from a straightforward application of the binomial theorem.

(ii) The second desired result can be obtained similarly, by considering (2.7) and (3.3) and by performing the same kind of computation.

We conclude this section by deriving the distribution of the ith cause of failure (2.4).

Proposition 3.3. Let $i=1,2$ . Under the assumptions of the model specified in Definition 3.1, we have

(3.13) \begin{equation}\mathbb{P}(\delta=i)=\dfrac{\lambda_{i}}{\lambda_{1}+\lambda_{2}}\nu\sum_{k=1}^{+\infty}({-}1)^{k-1}\dfrac{p_{k}}{(k-1)!}\sum_{j=0}^{k-1}({-}1)^{j}\binom{k}{j+1}(\nu j)_{k-1}.\end{equation}

Proof. With reference to (2.4) and to (3.11), we have

\begin{align*}\mathbb{P}(\delta=i)&=\lambda_{i}\nu(\lambda_{1}+\lambda_{2})^{\nu-1}\sum_{k=1}^{+\infty}({-}1)^{k-1}\dfrac{p_{k}}{(k-1)!}\sum_{h=0}^{k-1}\dfrac{(\lambda_{1}+\lambda_{2})^{\nu h}}{h!}\\&\quad{}\times\int_{0}^{+\infty}\,{\mathrm{e}}^{-(\lambda_{1}+\lambda_{2})^{\nu}t}t^{h}\,{\mathrm{d}} t\sum_{j=0}^{h}({-}1)^{j}\binom{h}{j}(\nu j)_{k-1}.\end{align*}

We apply formula 3.351.3 of [Reference Gradshteyn and Ryzhik13], and get

\begin{equation*}\mathbb{P}(\delta=i)=\dfrac{\lambda_{i}}{\lambda_{1}+\lambda_{2}}\nu\sum_{k=1}^{+\infty}({-}1)^{k-1}\dfrac{p_{k}}{(k-1)!}\sum_{h=0}^{k-1}\sum_{j=0}^{h}({-}1)^{j}\binom{h}{j}(\nu j)_{k-1}.\end{equation*}

By changing the order of the finite summations and using the hockey stick identity, the thesis immediately follows.

A typical problem in this area is to investigate the dependence between the failure time T and the cause of failure $\delta$ . For instance, Di Crescenzo and Meoli [Reference Di Crescenzo and Meoli10] presented a stochastic model for competing risks involving the Mittag–Leffler distribution, and inspired by fractional random growth phenomena, where T and $\delta$ prove to be independent.

However, in the present context, formula (3.13) makes it clear that the failure time T and the type of failure $\delta$ are in general dependent. In fact, due to (2.3) and (3.11), the probability density function of the failure time T is, for $t\geq 0,$

(3.14) \begin{align}f_{T}(t)&=\nu(\lambda_{1}+\lambda_{2})^{\nu}\,{\mathrm{e}}^{-(\lambda_{1}+\lambda_{2})^{\nu}t}\sum_{k=1}^{+\infty}({-}1)^{k-1}\dfrac{p_{k}}{(k-1)!} \notag \\&\quad\times\sum_{h=0}^{k-1}\dfrac{[(\lambda_{1}+\lambda_{2})^{\nu}t]^{h}}{h!}\sum_{j=0}^{h}({-}1)^{j}\binom{h}{j}(\nu j)_{k-1}.\end{align}

For there to be independence, it must be $f_{i}(t)=f_{T}(t)\mathbb{P}(\delta=i)$ , or, equivalently, $\mathbb{P}(\delta=i)={{{\lambda_{i}}/{(\lambda_{1}+\lambda_{2})}}}$ (see (3.11) and (3.14)). In general this is not true, as can easily be checked in some special cases.

4. Special cases

With reference to the stochastic model treated so far, in this section we present some examples for the probability distribution of the random lifetime T by specializing the distribution of the threshold M by means of (2.2). The starting point for our investigation is the following expansion, proved in Theorem 2.2 of [Reference Orsingher and Polito20]:

(4.1) \begin{equation}{\mathrm{e}}^{-\lambda^{\alpha}(1-u)^{\alpha}t}=\sum_{k=0}^{+\infty}u^{k}\dfrac{({-}1)^{k}}{k!}{}_1 \Psi_1\biggl[\begin{matrix}(1,\nu)\\(1-k,\nu)\end{matrix} \biggm| -\lambda^{\nu}t\biggr],\quad t\geq 0.\end{equation}

We examine four special cases for the distribution of the threshold M.

1. (1) M is geometrically distributed with parameter p, that is to say,

   \begin{equation*}\overline P_k=(1-p)^{k},\quad 0<p\leq 1.\end{equation*}
   To derive the distribution of T, put the previous survival probability in (3.12); then set $\lambda\,:\!=\, \lambda_{1}+\lambda_{2}$ and $u\,:\!=\, 1-p$ . Due to (4.1), we thus find that the random time T is exponentially distributed with parameter $[p(\lambda_{1}+\lambda_{2})]^{\nu}$ , that is,
   \begin{equation*}\overline{F}_T(t)={\mathrm{e}}^{-[p(\lambda_{1}+\lambda_{2})]^{\nu}t},\quad t\geq 0.\end{equation*}

2. (2) M is logarithmically distributed with parameter p. The survival probability of M reads

   \begin{equation*}\overline P_k=-\dfrac{B(p;\,k+1,0)}{\ln(1-p)},\quad 0<p<1,\end{equation*}
   where
   \[B(x;\,a,b)=\int_{0}^{x}u^{a-1}(1-u)^{b-1}\,{\mathrm{d}} u\]
   is the incomplete beta function. From (3.12) and (4.1), we get, for $t\geq 0$ ,
   \begin{align*}\overline{F}_T(t)&=\sum_{k=0}^{+\infty}\dfrac{({-}1)^{k}}{k!}\biggl({-}\dfrac{B(p;\,k+1,0)}{\ln(1-p)}\biggr){}_1 \Psi_1\biggl[\begin{matrix}(1,\nu)\\(1-k,\nu)\end{matrix} \biggm| -(\lambda_{1}+\lambda_{2})^{\nu}t\biggr]\\&=-\dfrac{1}{\ln(1-p)}\sum_{k=0}^{+\infty}\dfrac{({-}1)^{k}}{k!}\int_{0}^{p}z^{k}(1-z)^{-1}{}_1 \Psi_1\biggl[\begin{matrix}(1,\nu)\\(1-k,\nu)\end{matrix} \biggm| -(\lambda_{1}+\lambda_{2})^{\nu}t\biggr]\,{\mathrm{d}} z\\&=-\dfrac{1}{\ln(1-p)}\int_{0}^{p}(1-z)^{-1}\,{\mathrm{e}}^{-(\lambda_{1}+\lambda_{2})^{\nu}t(1-z)^{\nu}}\,{\mathrm{d}} z.\end{align*}
   Due to the Taylor expansion of the exponential function, we have
   \begin{align*}\overline{F}_T(t)&=1-\dfrac{1}{\ln(1-p)}\sum_{k=1}^{+\infty}\dfrac{[{-}(\lambda_{1}+\lambda_{2})^{\nu}t]^{k}}{k!}\int_{0}^{p}(1-z)^{-1+\nu k}\,{\mathrm{d}} z\\&=1-\dfrac{1}{\nu\ln(1-p)}\Biggl[\sum_{k=1}^{+\infty}\dfrac{[{-}(\lambda_{1}+\lambda_{2})^{\nu}t]^{k}}{k\cdot k!}-\sum_{k=1}^{+\infty}\dfrac{[{-}(\lambda_{1}+\lambda_{2})^{\nu}t(1-p)^{\nu}]^{k}}{k\cdot k!}\Biggr].\end{align*}
   In order to simplify the previous expression, we make use of the series expansion of the exponential integral
   \[E_{1}(z)=\int_{1}^{+\infty}\frac{{\mathrm{e}}^{-uz}}{u}\,{\mathrm{d}} u\]
   (see [Reference Abramowitz and Stegun1, 5.1.11]). The survival function of T finally reads
   \begin{equation*}\overline{F}_T(t)=\dfrac{E_{1}[(\lambda_{1}+\lambda_{2})^{\nu}t]-E_{1}[(\lambda_{1}+\lambda_{2})^{\nu}(1-p)^{\nu}t]}{\nu\ln(1-p)},\quad t\geq 0.\end{equation*}

The next two cases can be evaluated similarly, and thus we omit the calculations.

1. (3) If M is such that

   \begin{equation*}\overline P_k=\dfrac{\gamma(k+1;\,a,p)}{{\mathrm{e}}^{-a}-{\mathrm{e}}^{-p}},\quad 0\leq a<p\leq 1,\end{equation*}
   where
   \[\gamma(a,x_{1},x_{2})=\int_{x_{1}}^{x_{2}}\,{\mathrm{e}}^{-u}u^{a-1}\,{\mathrm{d}} u\]
   is the generalized incomplete gamma function, then
   \begin{equation*}\overline{F}_T(t)=\dfrac{1}{{\mathrm{e}}^{-a}-{\mathrm{e}}^{-p}}\int_{a}^{p}\,{\mathrm{e}}^{-z-(\lambda_{1}+\lambda_{2})^{\nu}t(1-z)^{\nu}}\,{\mathrm{d}} z,\quad t\geq 0.\end{equation*}

2. (4) If M is such that

   \begin{equation*}\overline P_k=\dfrac{\operatorname{Si}(k+1;\,a,p)}{\cos a-\cos p}, \quad 0\leq a<p\leq 1,\end{equation*}
   where
   \[\operatorname{Si}(a,x_{1},x_{2})=\int_{x_{1}}^{x_{2}}u^{a-1}\sin u\,{\mathrm{d}} u\]
   is the generalized sine integral, then
   \begin{equation*}\overline{F}_T(t)=\dfrac{1}{\cos a-\cos p}\int_{a}^{p}\sin z\,{\mathrm{e}}^{-(\lambda_{1}+\lambda_{2})^{\nu}t(1-z)^{\nu}}\,{\mathrm{d}} z,\quad t\geq 0.\end{equation*}

In Figure 1 we provide some plots of the survival function of T with reference to the four cases analyzed, for various choices of the parameters involved. In cases (3) and (4) the relevant integrals have been evaluated by means of routines available in MATLAB $^{\circledR}$ . We observe that in each of the four cases the survival function of T is increasing in $\nu$ . Moreover, the distributions that refer to cases (2), (3), and (4) have thinner tails than the exponential distribution, which refers to case (1).

Figure 1. Plots of the survival function of T with $\lambda_{1}+\lambda_{2}=1,$ $a=0,\,p=0.5$ , when the distribution of M is from left to right, from top to bottom, geometric, logarithmic, $p_k=\frac{\gamma(k;\,a,p)}{{\mathrm{e}}^{-a}-{\mathrm{e}}^{-p}}-\frac{\gamma(k+1;\,a,p)}{{\mathrm{e}}^{-a}-{\mathrm{e}}^{-p}}$ , $p_k=\frac{\operatorname{Si}(k;\,a,p)}{\cos a-\cos p}-\frac{\operatorname{Si}(k+1;\,a,p)}{\cos a-\cos p}$ .

5. Ageing notions

In this section we discuss some ageing characteristics of the random life-length T with respect to some of the special cases analyzed in Section 4. We briefly recall the main issues which we will refer to. See [Reference Marshall and Olkin17] or [Reference Navarro19] for a comprehensive coverage of the subject matter. A random lifetime T of an item is said to be NBU [NWU] (new better [worse] than used) if $\overline{F}_{T}(t+x)\leq [{\geq}]\, \overline{F}_{T}(t)\overline{F}_{T}(x)$ for all $x,\,t\geq 0$ . This means that the lifetime of a used item of age t is stochastically less [greater] than the lifetime of a brand new item. Moreover, if T has a density $f_{T}(t)$ for which the hazard rate $f(t)/\overline{F}_T(t)$ is increasing [decreasing] in t, the failure time T is said to have an increasing [decreasing] hazard rate (IHR [DHR]). Proposition B.8.a. of [Reference Marshall and Olkin17] states that if the density $f_{T}(t)$ is log-concave, then T is IHR, while if $f_{T}(t)$ is log-convex on $[0,+\infty)$ , then T is DHR. Clearly, in case (1) of Section 4, the density of T, being exponential, is both log-concave and log-convex.

The following result refers to case (2) of Section 4.

Proof. We consider case (2) of Section 4. By differentiation, the density of T reads

\begin{equation*}f_{T}(t)=\dfrac{{\mathrm{e}}^{-(\lambda_{1}+\lambda_{2})^{\nu}(1-p)^{\nu}t}-{\mathrm{e}}^{-(\lambda_{1}+\lambda_{2})^{\nu}t}}{t\nu\log(1-p)^{-1}}, \quad t>0.\end{equation*}

With reference to [Reference Guo and Qi14], set $a\,:\!=\, {\mathrm{e}}^{-(\lambda_{1}+\lambda_{2})^{\nu}t}$ and $b\,:\!=\, {\mathrm{e}}^{-(\lambda_{1}+\lambda_{2})^{\nu}(1-p)^{\nu}t}$ . It turns out that $b>a>0$ , so

\begin{equation*}f_{T}(t)=\dfrac{g_{(a,b)}(t)}{\nu\log(1-p)^{-1}},\quad t>0,\end{equation*}

where

\begin{equation*}g_{a,b}(t) =\begin{cases}\dfrac{b^{t}-a^{t}}{t} & \text{if $ t \neq 0$,} \\[5pt]b-a & \text{if $ t=0$.}\end{cases}\end{equation*}

Since ${\nu\log(1-p)^{-1}}>0$ and $g_{(a,b)}(t)$ is logarithmically convex on $(0,+\infty)$ (see [Reference Guo and Qi14, Section 2]), the density $f_{T}(t)$ is logarithmically convex as well. Therefore T is DHR.

It is well known that the NBU [NWU] property is implied by the IHR [DHR] property. The NBU ageing notion and its dual, the NWU one, have been expressed in the context of competing risks models in [Reference Di Crescenzo and Longobardi8]. In this case they are denoted as NBU $^{*}$ and NWU $^{*}$ and are defined as follows. Let $i\in\{1,2\}$ and let $X_i$ be the random variable that describes the lifetime of the item when $\delta=i$ . We have that $X_{i}$ is NBU $^{*}$ [NWU $^{*}$ ] if and only if $\mathbb{P}(T>t+x,\delta=i)\leq[{\geq}]\,\overline{F}_{T}(t)\mathbb{P}(T>x,\delta=i),\,x,t\geq 0.$ Namely, NBU $^{*}$ [NWU $^{*}$ ] expresses the positive [negative] ageing notion for a specific cause of failure. Finally, in the following proposition we consider the NBU $^{*}$ ageing notion in a special case. It is an immediate consequence of Remark 3.2 of [Reference Di Crescenzo and Longobardi8].