Competing risks and shock models governed by a generalized bivariate Poisson process

Abstract

We propose a stochastic model for the failure times of items subject to two external random shocks occurring as events in an underlying bivariate counting process. This is a special formulation of the competing risks model, which is of interest in reliability theory and survival analysis. Specifically, we assume that a system, or an item, fails when the sum of the two types of shock reaches a critical random threshold. In detail, the two kinds of shock occur according to a bivariate space-fractional Poisson process, which is a two-dimensional vector of independent homogeneous Poisson processes time-changed by an independent stable subordinator. Various results are given, such as analytic hazard rates, failure densities, the probability that the failure occurs due to a specific type of shock, and the survival function. Some special cases and ageing notions related to the NBU characterization are also considered. In this way we generalize certain results in the literature, which can be recovered when the underlying process reduces to the homogeneous Poisson process.

1. Introduction

Consider a shock model for a system subject to two types of shocks that occur according to a generalized bivariate Poisson process $\big(N_1^{\nu}(t),N_2^{\nu}(t)\big)$ , $t\geq 0$ . Namely, the underlying process is a bivariate counting process. Its components are independent homogeneous Poisson processes with time changed by means of a suitable non-decreasing Lévy process, i.e. the stable subordinator. For this model we analyze the related competing risks model. Specifically, we deal with the pair $(T,\delta)$ , where T is the failure time of the system and $\delta$ is the failure cause. In other terms, since $N_i^{\nu}(t)$ describes the number of shocks of type i that have arrived up to time t, $i=1,2$ , then T is the first-hitting time of the total number of shocks $N_1^{\nu}(t)+N_2^{\nu}(t)$ through a random integer-valued threshold, say M. The latter represents the number of fatal shocks that cause the failure of the system. Moreover, $\delta=i$ means that the cause of the failure is a shock of type i, $i=1,2$ , i.e. an instantaneous jump of $N_i^{\nu}(t)$ by which $N_1^{\nu}(t)+N_2^{\nu}(t)=M$ for the first time.

In many contexts of applied probability, counting processes are considered to model occurrences of shocks, claims, or other kinds of events. Besides the standard cases based on the classical Poisson process, as an example we recall here some recent developments centered on other kinds of processes. For instance, Cha and Giorgio [6] developed a new class of bivariate counting processes having a ‘marginal regularity’ property and allowing simultaneous occurrences of two types of events, and considered some applications to a shock model. Further kinds of underlying counting processes to shock models were investigated recently by Cha and Finkelstein [5] in the case of a generalized Pólya process, and by Di Crescenzo and Pellerey [11] in the case of a geometric counting process. In the present paper, the underlying model is based on a time-changed bivariate Poisson process, which is also known in the literature as a space-fractional Poisson process. This choice is motivated by the intention of constructing a more general model of the Poisson process, aiming to describe even more complex phenomena. For recent developments in this research area, see e.g. [18] and [22].

The classical competing risks model is largely adopted in reliability theory and survival analysis to describe system failures; see e.g. [3] and [7]. In some instances the risks can be identified with certain first-hitting times of stochastic processes. For instance, the paper [9] is concerned with the case when the risks are due to shocks governed by a bivariate homogeneous Poisson process with independent components.

This paper is oriented to the more general setting, which includes the presence of an underlying stable subordinator, say $\mathcal{A}^{\nu}(t)$ , $t\geq 0$ . The latter describes the time change of the bivariate Poisson process leading to a bivariate space-fractional Poisson process. This entails a more flexible and complex model, due to the presence of the index of stability $\nu\in(0,1]$ . The analysis is centered on the determination of the hazard rates of the two causes of failure. Such functions allow us to find the failure densities of the model, i.e. the sub-densities that describe the first-hitting times through M.

Another function of interest is the survival function of the system. In this case it is found to be expressed in terms of the Fox–Wright function. Moreover, the probabilities of the causes of failures are expressed in series form, depending on the distribution of the threshold M. Hence the specification of the distribution of M is essential to obtain closed-form expressions of the survival function of T in suitable cases.

The analysis of the model also involves the determination of the ageing characteristics of the random life-length T. In particular, we refer to certain typical ageing notions, i.e. the new better than used and new worse than used characteristics, which are defined properly within the competing risks model (see [8]).

The structure of the paper is as follows. In Section 2 we recall some basic facts regarding bivariate shock models. Then, in Section 3, we describe the bivariate space-fractional Poisson process, focusing on its main characteristics, and obtain expressions for the hazard rate, the failure density, the survival function, and the failure probability. Section 4 is devoted to the analysis of some special cases: we give closed-form expressions for the survival function of the random lifetime. In the final section we discuss some ageing notions with regard to the NBU and NBU $^{*}$ characterization.

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 [2, 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 [23]),

\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 [4]. 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 [12] and [22]. 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 [20] and [21] 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 [4, 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. [16, 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 [4], 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 [15] 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 [9, 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).

Proposition 3.2. Under the assumptions of the model in Definition 3.1, for $t\geq 0$ and $i=1,2,$ we have the following results.

1. (i) The failure densities can be expressed as

   (3.11) \begin{align}f_{i}(t)&=\lambda_{i}\nu(\lambda_{1}+\lambda_{2})^{\nu-1}\,{\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}

2. (ii) The survival function of T reads

   (3.12) \begin{equation}\overline{F}_T(t)=\sum_{k=0}^{+\infty}\dfrac{({-}1)^{k}}{k!}\overline P_k\,{}_1 \Psi_1\biggl[\begin{matrix}(1,\nu)\\(1-k,\nu)\end{matrix} \biggm| -(\lambda_{1}+\lambda_{2})^{\nu}t\biggr].\end{equation}

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 [13], 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 [10] 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 [20]:

(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 [1, 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 [17] or [19] 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 [17] 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.

Proposition 5.1. If M has a logarithmic distribution with parameter $p,\,0<p< 1$ , then T is DHR.

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 [14], 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 [14, 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 [8]. 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 [8].

Proposition 5.2. If M has a geometric distribution with parameter $p,\,0<p\leq 1$ , then at most one of the following statements holds:

1. (i) $X_{1}$ and $X_{2}$ are simultaneously NBU $^{*}$ and NWU $^{*}$ ,

2. (ii) one of the variables $X_{1}$ and $X_{2}$ is NBU $^{*}$ and the other is NWU $^{*}$ .

Proof. We consider case (1) of Section 4. The result is an immediate consequence of Remark 3.2 of [8], since T is exponentially distributed.

6. Conclusions

Motivated by applications in reliability theory, medical research, insurance, and economics, to name just a few, in this paper we have studied a counterpart of a known bivariate stochastic shock model based on the space-fractional Poisson process. First we described the general setting of the model, then we obtained the hazard rates and showed a wide range of related results. Next we obtained explicit formulae for the survival function of the random lifetime in four special cases. Finally we discussed ageing notions with respect to the NBU and NBU $^{*}$ characterizations. For the sake of simplicity, we restricted our attention to an underlying bivariate space-fractional Poisson process, but our analysis can be easily generalized by considering a model based on the multivariate space-fractional Poisson process [4]. Follow-up research might relate to other ways of generalizing the underlying bivariate counting process. Indeed, we might consider the composition of a two-dimensional vector of independent classical Poisson processes with an independent tempered stable subordinator, with an independent subordinator associated with a general Bernštein function, or with the components of an independent bivariate stable subordinator.