# A bivariate failure time model with random shocks and mixed effects

Sophie Mercier, Ha Pham

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

## To cite this version:

Sophie Mercier, Ha Pham. A bivariate failure time model with random shocks and mixed effects. Journal of Multivariate Analysis, 2017, 153, pp.33-51. 10.1016/j.jmva.2016.09.008. hal-01361399

HAL Id: hal-01361399

https://hal.science/hal-01361399

Submitted on 23 Sep 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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>
A bivariate failure time model with random
shocks and mixed effects

Sophie MERCIER*

Laboratoire de Mathématiques et de leurs Applications – Pau (UMR CNRS 5142)

Bâtiment IPRA, Université de Pau et des Pays de l'Adour,
Avenue de l'Université, F-64013 PAU cedex, FRANCE

Tel.: +33 (0)5 59 40 75 37 - Fax: +33 (0)5 59 40 75 55

Email: sophie.mercier@univ-pau.fr

Hai Ha PHAM

Department of Economical Mathematics
Banking University of Hochiminh city

39 Ham Nghi, District 1, Hochiminh city, Vietnam

Email: phamhaiha09@gmail.com

September 5, 2016

Abstract

Two components are considered, which are subject to common ex-
ternal and possibly fatal shocks. The lifetimes of both components are
characterized by their hazard rates. Each shock can cause the immediate
failure of either one or both components. Otherwise, the hazard rate of
each component is increased by a non fatal shock of a random amount,
with possible dependence between the simultaneous increments of the two
failure rates. An explicit formula is provided for the joint distribution of
the bivariate lifetime. Aging and positive dependence properties are de-
scribed, thereby showing the adequacy of the model as a bivariate failure
time model. The influence of the shock model parameters on the bivari-
ate lifetime is also studied. Numerical experiments illustrate and complete
the study. Moreover, an estimation procedure is suggested in a parametric
framework, under a specific observation scheme.

***Keywords:*** Aging properties; Bivariate new better than used; Bivariate non-homogeneous compound Poisson process; Hazard rate order; Hazard rate process; Maximum likelihood estimator; Multivariate total positivity; Positive dependence properties; Reliability; Stochastic order.

<footnote>*Corresponding author</footnote>
# 1 Introduction

Components are considered, which are made dependent through shocks arising from a common external environment. Except for these shocks, the components are assumed to evolve independently. The components may fail either when the shocks occur or between them. Each shock may hence lead to simultaneous failures, also called common cause failures. In practice, there are many possible causes for such shocks. For example, according to [15], they can be due to

> "interfaces, the environment, and major adverse events. The interfaces include power, cooling, material inputs, and external controls. The environment can produce excess temperature, pressure, vibration, impact, noise, and contamination. Events include earthquake, tsunami, hurricane, tornado, flood, and blizzard."

According to [20], such shocks can also be due to "the failure of an external piece of hardware (... such as a pipe leak), or a human error (... such as the miscalibration of an instrument)." Note also that according to [15], "about ten per cent of all failures are usually identified as common cause."

A classical way to model common cause failures is to use the binomial failure rate (BFR) model introduced in [35] (see also [3]), where each unit can fail independently or due to a common shock, with fatality for the different units determined by independent Bernoulli trials, and where the shocks arrive according to an independent homogeneous Poisson process. This model has proved its usefulness in applications; see, e.g., [3, 35] and references therein. However, it has also been seen to be too restrictive in several ways: First, "the assumption that the components will fail independently of each other, given that a shock has occurred, represents a rather serious limitation, and this assumption is often not satisfied in practice" [30, p. 223]. Secondly, the only dependence between components included in the BFR model is due to the possibly synchronized failures whereas, according to [30, p. 213], "even if common cause failures are caused by a common cause, they do not need to occur at the same time. A rather long time between failures does not necessarily mean that there is no dependency between the failure events."

Indeed, one can think that a shock can simultaneously weaken components without necessarily entailing their immediate failures, e.g., think of a temporary high temperature/voltage/pressure... This situation can also reflect the case of some components of an aircraft engine, which are subject to external shocks during takeoff and landing. The shocks can simultaneously increase the deterioration of the components without necessarily leading to their immediate failure. Examples are not restricted to the reliability domain either. In insurance, [11] highlights the necessity of taking into account some dependence between the lifetimes of husband and wife (or twins), because of the common shocks (for
example accidents or contagious diseases) that they may suffer in their common
life. There is thus a need for describing such situations. As will be seen in the
following, the model proposed in this paper does not share with the BFR model
the two above mentioned restrictions, in addition to improving over it on several
other points.

Many other shock models than the BFR have been developed in the reliability literature; see, e.g., [9, 23, 29, 34]. In these writings, common cause failures are typically described through the Marshall-Olkin multivariate exponential distribution family [24] or extensions thereof. In the bivariate case, three independent homogeneous Poisson processes govern the occurrence of shocks, which all are fatal. Two processes impact a single component each while the third one impacts both components simultaneously. This model is fundamental in reliability theory and it remains an important source of inspiration for much research; see, e.g., [4, 5, 13, 18, 21, 22, 36]. In cumulative shock models [25], a shock simultaneously increases some intrinsic characteristics (such as hazard rate, deterioration level, age, etc.) of the components, with possible dependence between the simultaneous random increments.

Recently, an approach for a realistic shock model has been studied in [6]. This article considers a univariate model that takes into account shocks with a mixed effect: a shock can be fatal to the system with a probability depending on the shock's arrival time; a non fatal shock increases the system failure rate of a random increment. An extension of this model to multi-component systems has yet to be considered in full generality. Exceptions are [5] and [27] which both provide an extension to the case of competing soft and sudden failures, where soft failures refer to the reaching of some critical threshold for some degradation level, and sudden failures to accidental failures, characterized by a failure rate. In [27], the authors are mainly interested in the univariate survival function of the competing risks and in the influence of the stressing environment on the system lifetime. As for [5], it focuses on the bivariate survival function of the competing soft and sudden failures, and on the induced conditional properties. The authors also consider a condition-based preventive maintenance policy.

The present paper deals with a complex bivariate failure time model where both lifetimes are characterized by failure rates (only sudden failures). Contrary to [27], the fatality of a shock may differ between the components. Also, unlike the BFR model, the Bernoulli trials governing the fatality for the two components are not independent and even non fatal shocks can entail some dependence between components. Here, a non fatal shock increases the failure rates of the surviving components of a random increment (simultaneous weakening of the components), with possible dependence between simultaneous increments. No assumption is made on the joint distribution of simultaneous increments so that the model includes the possibility of getting increments of the
form $(0, x_2)$, $(x_1, 0)$ or $(x_1, x_2)$ (with $x_1, x_2 > 0$) depending on the realization. Hence, a non-fatal shock can act on a single component or on both of them, simultaneously. As in [5, 27], the probability for a shock to be fatal depends on the shock's arrival time and the occurrence of shocks is modeled through a non-homogeneous Poisson process.

The paper is organized as follows. The model is specified in Section 2. An explicit formula for the joint distribution of the bivariate lifetime is provided in Section 3, as well as aging and positive dependence properties (Bivariate New Better than Used property and Multivariate Total Positivity of order two). The influence of the shock model parameters on the bivariate lifetime is also studied. Numerical experiments are presented in Section 4 which not only illustrate the theoretical results and but also show that stronger results do not hold under the same assumptions. An attempt is made in Section 5 to provide a few elements for parametric estimation of the model under a specific observation scheme. Concluding remarks are given in Section 6.

## 2 The model

Two components are considered. The lifetime of component $i = 1, 2$ is characterized by its intrinsic hazard rate function $h_i(t)$ or by the corresponding cumulative hazard rate function $H_i(t) = \int_0^t h_i(u)du$. Stresses due to the external environment come in the form of shocks, independently of the components' intrinsic deterioration. Except for shocks, the components are assumed to behave independently. The shocks occur at time $T_1, T_2, \dots$ according to a non-homogeneous Poisson process $(N_t)_{t \ge 0}$ with intensity $\lambda(t)$ and cumulative intensity $\Lambda(t) = \int_0^t \lambda(x)dx$. The $n$th shock, which occurs at time $T_n$, increases the hazard rate of the $i$th component by a random amount $V_n^{(i)}$. The bivariate increments $V_1 = (V_1^{(1)}, V_1^{(2)})$, $V_2 = (V_2^{(1)}, V_2^{(2)})$, \dots are assumed to be independent and identically distributed (i.i.d.), as well as independent of the shocks arrival times $(T_n)_{n \ge 1}$. Simultaneous increments $V_n^{(1)}$ and $V_n^{(2)}$ may however be dependent. Furthermore, a shock can be fatal, and can possibly induce the immediate failure of either one or both components. The fatality of a shock does not depend on the system intrinsic deterioration but depends on the shock's arrival time. The following notations are used:

$p_{00}(T_n)$: probability that the shock at time $T_n$ induces the simultaneous failure of both components,

$p_{11}(T_n)$: probability that the shock at time $T_n$ induces no failure at all among the two components,

$p_{01}(T_n)$: probability that the shock at time $T_n$ is fatal only for the first component,
$p_{10}(T_n)$: probability that the shock at time $T_n$ is fatal only for the second component.

Therefore, $p_{00}(\cdot) + p_{01}(\cdot) + p_{10}(\cdot) + p_{11}(\cdot) = 1$ by definition.

The common distribution of the i.i.d. random vectors $V_1 = (V_1^{(1)}, V_1^{(2)}), \dots, V_2 = (V_2^{(1)}, V_2^{(2)}), \dots$ is denoted by $\mu(dv_1, dv_2)$. When subscript $n$ is unnecessary, we drop it and set $V = (V^{(1)}, V^{(2)})$ to be a generic copy of $V_n = (V_n^{(1)}, V_n^{(2)})$. For $j = 1, 2$, the distribution of $V^{(j)}$ is denoted by $\mu_j(dv_j)$.

For $i = 1, 2$, we set $(A_t^{(i)})_{t \ge 0}$ to be the univariate compound Poisson process defined, for all $t \ge 0$, by

$$A_t^{(i)} = \sum_{n=1}^{N_t} V_n^{(i)},$$

where an empty sum equals 0 by convention. We also introduce the bivariate compound Poisson process $(A_t)_{t \ge 0}$, where, for all $t \ge 0$,

$$A_t = (A_t^{(1)}, A_t^{(2)}) = \left( \sum_{n=1}^{N_t} V_n^{(1)}, \sum_{n=1}^{N_t} V_n^{(2)} \right) = \sum_{n=1}^{N_t} V_n.$$

Now, let $\mathcal{F} = \sigma(A_s, s \ge 0)$ be the $\sigma$-field generated by $(A_t)_{t \ge 0}$. Provided that the $i$th ($i = 1, 2$) component is working up to time $t$, the conditional hazard rate of this component at time $t$ given $\mathcal{F}$ is

$$r_i(t) = h_i(t) + A_t^{(i)}.$$

For $i = 1, 2$, let $\tau_i$ be the lifetime of the $i$th component without taking into account the possibility of fatal shocks and let $\xi_i$ be the time of the first fatal shock for the $i$th component. We have

$$ \begin{aligned} E(\mathbf{1}_{\{\tau_i > t\}} | \mathcal{F}) &= \exp\{-H_i(t)\} \exp\left\{-\int_0^t A_s^{(i)} ds\right\} \\ &= \exp\{-H_i(t)\} \exp\left\{-\sum_{k=1}^{N_t} V_k^{(i)} (t - T_k)\right\} \\ &= \exp\{-H_i(t)\} \exp\left\{-\sum_{k=1}^{\infty} V_k^{(i)} (t - T_k)^+\right\} \end{aligned} \quad (1) $$

and

$$E(\mathbf{1}_{\{\xi_i > t\}} | \mathcal{F}) = \prod_{k=1}^{N_t} q_i(T_k)$$

where $q_1(\cdot) = p_{11}(\cdot) + p_{10}(\cdot)$ and $q_2(\cdot) = p_{11}(\cdot) + p_{01}(\cdot)$.

Given $\mathcal{F}$, the random variables $\tau_1$ and $\tau_2$ are assumed to be conditionally independent, and conditionally independent of $\xi_1$ and $\xi_2$. However, based on the fact that the two components may fail simultaneously at each shock, $\xi_1$ and
$\xi_2$ are not conditionally independent given $\mathcal{F}$. Given $\mathcal{F}$, their conditional joint survival function is

$$ E(\mathbf{1}_{\{\xi_1 > s\}} \mathbf{1}_{\{\xi_2 > t\}} | \mathcal{F}) = \begin{cases} \prod_{i=1}^{N_t} p_{11}(T_i) \prod_{i=N_t+1}^{N_s} q_1(T_i) & \text{if } s \ge t, \\ \prod_{i=1}^{N_s} p_{11}(T_i) \prod_{i=N_s+1}^{N_t} q_2(T_i) & \text{if } s < t \end{cases} \qquad (2) $$

where an empty product is equal to 1 by convention.

Letting $Y = (Y_1, Y_2)$ denote the bivariate lifetime of the two components, it is easy to see that $Y_i = \min(\tau_i, \xi_i)$ for $i = 1, 2$.

# 3 Theoretical results

## 3.1 Joint survival function

**Proposition 1** The joint survival function of $Y = (Y_1, Y_2)$ is given by

$$ \bar{F}_Y(s,t) = e^{-H_1(s)-H_2(t)-\Lambda\{\max(s,t)\}} \exp\left\{\int_0^{\min(s,t)} \tilde{\mu}(s-w, t-w)p_{11}(w)\lambda(w)dw + \mathbf{1}_{\{t\le s\}}\int_t^s \tilde{\mu}_1(s-w)q_1(w)\lambda(w)dw + \mathbf{1}_{\{t>s\}}\int_s^t \tilde{\mu}_2(t-w)q_2(w)\lambda(w)dw\right\} \quad (3) $$

for all $s, t \ge 0$, where $\tilde{\mu}$ stands for the Laplace transform of the bivariate distribution $\mu$ of $V = (V^{(1)}, V^{(2)})$, defined, for all $x_1, x_2 \ge 0$, by

$$ \tilde{\mu}(x_1, x_2) = \iint_{\mathbb{R}_+^2} e^{-x_1 v_1 - x_2 v_2} \mu(dv_1, dv_2) = E(e^{-x_1 V^{(1)} - x_2 V^{(2)}}), $$

and for $i = 1, 2$, the symbol $\tilde{\mu}^{(i)}$ stands for the Laplace transform of the univariate distribution $\mu^{(i)}$ of $V^{(i)}$ defined, for all $x_i \ge 0$, by

$$ \tilde{\mu}_i(x_i) = \int_0^\infty e^{-x_i v_i} \mu_i(dv_i) = E(e^{-x_i V^{(i)}}). $$

**Proof.** We have

$$ \begin{aligned} \bar{F}_Y(s,t) &= \Pr(Y_1 > s, Y_2 > t) \\ &= \Pr(\tau_1 > s, \tau_2 > t, \xi_1 > s, \xi_2 > t) \\ &= E\{E(\mathbf{1}_{\{\tau_1 > s\}} | \mathcal{F}) E(\mathbf{1}_{\{\tau_2 > t\}} | \mathcal{F}) E(\mathbf{1}_{\{\xi_1 > s\}} \mathbf{1}_{\{\xi_2 > t\}} | \mathcal{F})\} \end{aligned} $$

based on the conditional independence of $\tau_1, \tau_2$ and $(\xi_1, \xi_2)$ given $\mathcal{F}$.
We first consider the case $s \ge t$. Substituting (1) and (2) into the previous equation, we get

$$
\begin{align}
\bar{F}_Y(s,t) &= e^{-H_1(s)-H_2(t)} \mathrm{E} \left\{ e^{-\sum_{i=1}^{\infty} V_i^{(1)}(s-T_i)^+} e^{-\sum_{i=1}^{\infty} V_i^{(2)}(t-T_i)^+} \prod_{i=1}^{N_t} p_{11}(T_i) \prod_{i=N_{t+1}}^{N_s} q_1(T_i) \right\} \tag{4} \\
&= e^{-H_1(s)-H_2(t)} \mathrm{E} \left\{ e^{-\sum_{i=1}^{\infty} \psi_{s,t}(V_i^{(1)}, V_i^{(2)}, T_i)} \right\}
\end{align}
$$

where

$$
\psi_{s,t}(v_1, v_2, w) = v_1(s-w)^+ + v_2(t-w)^+ - \ln p_{11}(w) \mathbf{1}_{[0,t]}(w) - \ln q_1(w) \mathbf{1}_{(t,s]}(w) \tag{5}
$$

for all $v_1, v_2, w \in \mathbb{R}_+$.

The sequence $(V_i^{(1)}, V_i^{(2)}, T_i)_{i \ge 1}$ can be seen as the points of a Poisson ran-
dom measure with intensity $\nu$($dv_1, dv_2, dw$) = $\mu$($dv_1, dv_2$)$\lambda$($w$)$dw$; see, e.g.,
Corollary 3.5 in Chapter 6 of [7]. The formula for Laplace functionals of a
Poisson random measure (Theorem 2.9 in Chapter 6 of [7]) yields

$$
\mathrm{E}\{e^{-\sum_{i=1}^{\infty} \psi_{s,t}(V_i^{(1)}, V_i^{(2)}, T_i)}\} = \exp \left[ - \iiint_{\mathbb{R}_+^3} \{1 - e^{-\psi_{s,t}(v_1, v_2, w)}\} \nu dv_1 dv_2 dw \right].
$$

Substituting $\psi_{s,t}$ by its expression given in (5), we get

$$
\begin{align*}
& \iiint_{\mathbb{R}_+^3} \{1 - e^{-\psi_{s,t}(v_1, v_2, w)}\} \nu(dv_1, dv_2, dw) \\
&= \Lambda(s) + \int_0^t \left\{ \int_0^\infty \int_0^\infty e^{-v_1(s-w)-v_2(t-w)} \mu(dv_1, dv_2) \right\} p_{11}(w) \lambda(w) dw \\
&\quad + \int_t^s \left\{ \int_0^\infty \int_0^\infty e^{-v_1(s-w)} \mu(dv_1, dv_2) \right\} q_1(w) \lambda(w) dw \\
&= \Lambda(s) + \int_0^t \tilde{\mu}(s-w, t-w) p_{11}(w) \lambda(w) dw + \int_t^s \tilde{\mu}_1(s-w) q_1(w) \lambda(w) dw.
\end{align*}
$$

The case $s < t$ is obtained by symmetry and this provides Formula (3). $\blacksquare$

**Remark 1** One can easily derive the univariate lifetimes of both components.
For example, the survival function for the first component is

$$
\bar{F}_{Y_1}(s) = \bar{F}_Y(s, 0) = e^{-\Lambda(s)-H_1(s)+\int_0^s \tilde{\mu}_1(s-w)q_1(w)\lambda(w)dw}.
$$

This formula is stated in Theorem 1 from [6] and thus Proposition 1 appears as
a generalization of this theorem. Note, however, that the present use of Poisson
random measures allows for a shorter proof.
**Remark 2** Based on Formula (3), it is clear that the bivariate survival function $\bar{F}_Y$ will generally not be differentiable at a point of the main diagonal $\Delta = \{(x, x) : x \in \mathbb{R}_+\}$. To be more specific, the distribution of $Y$ admits a non absolutely continuous part with support $\Delta$ and an absolutely continuous part for the remaining of the distribution, which writes:

$$ \mathrm{Pr}_Y(ds, dt) = f_Y(s, s) ds \delta_s(dt) + f_Y(s, t) ds dt $$

where $\delta_s(dt)$ stands for the Dirac mass at point $s$ and where

$$ f_Y(s,t) = \begin{cases} \frac{\partial^2}{\partial s \partial t} \bar{F}_Y(s,t) & \text{if } s \neq t \\ \lim_{h \to 0^+} \frac{1}{2h} \mathrm{Pr} \left( s-h < Y_1 \leq s+h, s-h < Y_2 \leq s+h \right) & \text{if } s=t. \end{cases} $$

In the most general case, $f_Y(s, t)$ can be computed through the following:

$$ f_Y(s, s) = \lim_{h \to 0^+} \frac{1}{2h} \left\{ \bar{F}_Y(s+h, s+h) - \bar{F}_Y(s-h, s+h) - \bar{F}_Y(s+h, s-h) + \bar{F}_Y(s-h, s-h) \right\} $$

when $s = t$, and through numerical differentiation when $s \neq t$.

We now provide two examples where an explicit expression is available for the joint survival function.

**Example 1** Let $h_1 = h_2 = 0$ and assume that for all $i, j = 0, 1, p_{i,j}$ is constant, as is the size increment, viz. $V = (v_1, v_2) \in \mathbb{R}_+^2$. This leads to

$$ \bar{F}_Y(s,t) = e^{-\Lambda\{\max(s,t)\}} \exp \left\{ p_{11} \int_0^{\min(s,t)} e^{-v_1(s-w)-v_2(t-w)} \lambda(w) dw + \mathbf{1}_{\{t \le s\}} q_1 \int_t^s e^{-v_1(s-w)} \lambda(w) dw + \mathbf{1}_{\{t > s\}} q_2 \int_s^t e^{-v_2(t-w)} \lambda(w) dw \right\}. $$

Specializing to $\lambda(w) = e^{\alpha w}$ with $\alpha \in \mathbb{R} \setminus \{0, -v_1, -v_2, -v_1 - v_2\}$, we easily get

$$ \bar{F}_Y(s,t) = \exp \left\{ -\frac{e^{\alpha \max(s,t)} - 1}{\alpha} + p_{11} e^{-v_1 s - v_2 t} \frac{e^{(v_1+v_2+\alpha) \min(s,t)} - 1}{v_1 + v_2 + \alpha} \right\} \\ \times \exp \left\{ \mathbf{1}_{\{t \le s\}} q_1 e^{-v_1 s} \frac{e^{(v_1+\alpha)s} - e^{(v_1+\alpha)t}}{v_1 + \alpha} + \mathbf{1}_{\{t > s\}} q_2 e^{-v_2 t} \frac{e^{(v_2+\alpha)t} - e^{(v_2+\alpha)s}}{v_2 + \alpha} \right\}. $$

**Example 2** Let $h_1 = h_2 = 0$ and assume that for all $i, j = 0, 1, p_{i,j}$ is constant, as is $\lambda$. The size $V = (V^{(1)}, V^{(2)})$ of an increment has the Marshall-Olkin distribution with parameters $(\lambda_1, \lambda_2, \lambda_{12})$, so that its bivariate survival function
is given, for all $x_1, x_2 \ge 0$, by

$$ \bar{F}_{(V^{(1)}, V^{(2)})}(x_1, x_2) = e^{-\lambda_1 x_1 - \lambda_2 x_2 - \lambda_{12} \max(x_1, x_2)}; $$

see [24]. Using [13, Theorem 3.4] for the moment-generating function of the Marshall-Olkin bivariate exponential distribution with parameter $(\alpha, \lambda_1, \lambda_2, \lambda_{12})$ and letting $\alpha \to \infty$, we get

$$ \tilde{\mu}(s, t) = \frac{1}{\gamma + s + t} \left( \frac{\lambda_1 \gamma_2}{\gamma_2 + t} + \frac{\lambda_2 \gamma_1}{\gamma_1 + s} + \lambda_{12} \right) $$

where $\gamma = \lambda_1 + \lambda_2 + \lambda_{12}$, $\gamma_1 = \lambda_1 + \lambda_{12}$, $\gamma_2 = \lambda_2 + \lambda_{12}$ and, for $i = 1, 2$,

$$ \tilde{\mu}_i(s) = \frac{\gamma_i}{\gamma_i + s}. $$

For $s \le t$, easy (but tedious) computations yield

$$ \begin{aligned} \bar{F}_Y(s,t) &= e^{-\lambda t + p_{11}\lambda \int_0^s \tilde{\mu}(s-w,t-w) dw + q_2\lambda \int_s^t \tilde{\mu}_2(t-w) dw} \\ &= e^{-\lambda t} \left( \frac{\gamma_2 + t}{\gamma_2 - s + t} \frac{\gamma - s + t}{\gamma + s + t} \right)^{p_{11}\lambda \frac{\lambda_1 \gamma_2}{\gamma - 2\gamma_2 + s - t}} \left( \frac{\gamma_1 + s}{\gamma_1} \frac{\gamma - s + t}{\gamma + s + t} \right)^{p_{11}\lambda \frac{\lambda_2 \gamma_1}{\gamma - 2\gamma_1 - s + t}} \\ & \quad \times \left( \frac{\gamma + s + t}{\gamma - s + t} \right)^{p_{11}\lambda \frac{\lambda_2}{2}} \left( \frac{\gamma_2 - s + t}{\gamma_2} \right)^{q_2 \lambda \gamma_2} \end{aligned} $$

(with a similar expression for $s > t$). Furthermore, one has, for $i = 1, 2$,

$$ \bar{F}_{Y_i}(s) = e^{-\lambda s} \left( \frac{\gamma_i + s}{\gamma_i} \right)^{q_i \lambda \gamma_i}. $$

## 3.2 A bivariate aging property

Let us first recall that the univariate lifetime $Y_i$ of a component with cumulative hazard rate $H_i$ has the New Better than Used (NBU) property if and only if $H_i$ is super-additive, i.e., $H_i(x_i + y_i) \ge H_i(x_i) + H(y_i)$ for all $x_i, y_i \ge 0$, which is equivalent to $\bar{F}_{Y_i}(x_i + y_i) \le \bar{F}_{Y_i}(x_i)\bar{F}_{Y_i}(y_i)$ for all $x_i, y_i \ge 0$. This property means that the lifetime of a new component is stochastically larger than the lifetime of an older component (see, e.g., [19]) and it is a so-called aging property. In the bivariate setting, there are different ways to define a Bivariate New Better than Used (BNBU) property (see Section 8.5 in [19] for example). We provide here sufficient conditions under which $Y$ has the BNBU property in the following sense:

$$ \bar{F}_Y(s_1 + s_2, t_1 + t_2) \le \bar{F}_Y(s_1, t_1)\bar{F}_Y(s_2, t_2) \quad (6) $$

for all $s_1, s_2, t_1, t_2 \ge 0$ such that $(s_1 - t_1)(s_2 - t_2) \ge 0$.
**Proposition 2** Assume that $H_1, H_2, \Lambda$ are super-additive and that the functions $p_{11}(w)\lambda(w)$ and $q_i(w)\lambda(w)$ ($i = 1,2$) are non-increasing. Then $Y$ is BNBU.

**Proof.** For any $s \le t$, let us write

$$ \bar{F}_Y(s,t) = e^{-H_1(s)-H_2(t)-\Lambda(t)+A(s,t)+B(s,t)}, \quad (7) $$

where

$$ A(s,t) = \int_{0}^{s} \tilde{\mu}(s-w, t-w)p_{11}(w)\lambda(w)dw, $$

$$ B(s,t) = \int_{s}^{t} \tilde{\mu}_{2}(t-w)q_{2}(w)\lambda(w)dw. $$

The objective is to prove (6). Let us first suppose that $s_i \le t_i$ for $i = 1,2$. Because $H_1, H_2$ and $\Lambda$ are super-additive, we already know that

$$ e^{-H_1(s_1)-H_2(t_1)-\Lambda(t_1)}e^{-H_1(s_2)-H_2(t_2)-\Lambda(t_2)} \ge e^{-H_1(s_1+s_2)-H_2(t_1+t_2)-\Lambda(t_1+t_2)}. $$

Thus, using (7) for $(s_1, t_1)$, $(s_2, t_2)$ and $(s_1 + s_2, t_1 + t_2)$, we can see that it is sufficient to prove that

$$ A(s_1, t_1) + A(s_2, t_2) \ge A(s_1 + s_2, t_1 + t_2), \quad (8) $$

$$ B(s_1, t_1) + B(s_2, t_2) \ge B(s_1 + s_2, t_1 + t_2). \quad (9) $$

Writing $\int_{0}^{s_1+s_2} = \int_{0}^{s_1} + \int_{s_1}^{s_2}$, we have

$$ A(s_1 + s_2, t_1 + t_2) = a(s_1, s_2, t_1, t_2) + b(s_1, s_2, t_1, t_2) \quad (10) $$

with

$$ \begin{aligned} a(s_1, s_2, t_1, t_2) &= \int_{0}^{s_1} \tilde{\mu}(s_1 + s_2 - w, t_1 + t_2 - w)p_{11}(w)\lambda(w)dw \\ &\le \int_{0}^{s_1} \tilde{\mu}(s_1 - w, t_1 - w)p_{11}(w)\lambda(w)dw \\ &= A(s_1, t_1) \end{aligned} \quad (11) $$
(by non-increasingness of $\tilde{\mu}$) and

$$ 
\begin{aligned}
b(s_1, s_2, t_1, t_2) &= \int_{s_1}^{s_1+s_2} \tilde{\mu}(s_1 + s_2 - w, t_1 + t_2 - w) p_{11}(w) \lambda(w) dw \\
&= \int_{0}^{s_2} \tilde{\mu}(s_2 - u, t_1 + t_2 - s_1 - u) p_{11}(u + s_1) \lambda(u + s_1) du \\
&\le \int_{0}^{s_2} \tilde{\mu}(s_2 - u, t_2 - u) p_{11}(u) \lambda(u) du \\
&= A(s_2, t_2)
\end{aligned} 
\quad (12) $$

(setting $u = w - s_1$ in the second line, and using the non-increasingness of both $p_{11}(w)\lambda(w)$ and $\tilde{\mu}$ for the inequality). Gathering (10), (11) and (12) provides (8).

Inequality (9) can be proved in a similar way by writing

$$ \int_{s_1+s_2}^{t_1+t_2} = \int_{s_1+s_2}^{t_1+s_2} + \int_{t_1+s_2}^{t_1+t_2} $$

in $B(s_1 + s_2, t_1 + t_2)$ and showing that the first (resp. second) term is smaller than $B(s_1, t_1)$ (resp. $B(s_2, t_2)$) under non increasingness of $q_2(w)\lambda(w)$.

The case $s_i \ge t_i$, $i = 1, 2$ is similar, given the non-increasingness of $q_1(w)\lambda(w)$.

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

These assumptions are quite natural to guarantee that $Y$ is BNBU. Indeed, the super-additivity of $H_1$ and $H_2$ states that without taking into account the effect of the external environment, both components are NBU. The super-additivity of $\Lambda$ means that the shocks are more and more frequent. The decreasing property of the functions $p_{11}(w)\lambda(w)$, $q_i(w)\lambda(w)$, $i = 1, 2$ implies that non-fatal shocks appear more and more rarely, and consequently, as shocks are more and more frequent, that the probability of a fatal shock is increasing.

**Remark 3** Specializing to $s = t$, it is easy to see that $\bar{F}_Y(s, s)$ is the probability that both components are still working at time $s$, so that $\bar{F}_Y(s, s)$ appears as the reliability of a two-unit series system at time $s$. Then, writing the BNBU property for $s_i = t_i$ ($i = 1, 2$), it can be seen that, under the conditions of Proposition 2, the lifetime of the series system has the (univariate) NBU property.

Now, recall that a stronger BNBU property is provided by

$$ \bar{F}_Y(s_1 + s_2, t_1 + t_2) \le \bar{F}_Y(s_1, t_1) \bar{F}_Y(s_2, t_2) \quad (13) $$

for all $s_1, s_2, t_1, t_2 \ge 0$. One may consequently wonder whether this stronger BNBU property (or the even stronger property from [26]) could be valid under the assumption of Proposition 2. The answer is negative as will be shown in Example 4; see Section 4.
### 3.3 A positive dependence property

Given a bivariate random vector $Y = (Y_1, Y_2)$, we recall that $\bar{F}_Y$ is multivariate totally positive of order 2 (MTP2) if

$$ \forall_{x,y \in \mathbb{R}_+^2} \bar{F}_Y(y) \bar{F}_Y(x) \le \bar{F}_Y(x \vee y) \bar{F}_Y(x \wedge y) $$

or equivalently if

$$ \bar{F}_Y(x_1, x_2) \bar{F}_Y(y_1, y_2) \le \bar{F}_Y(x_1, y_2) \bar{F}_Y(y_1, x_2) $$

for all $x_1 \le y_1$ and $x_2 \ge y_2$.

We also recall that $Y_2$ is said to be Right-Tail-Increasing in $Y_1$, written RTI($Y_2|Y_1$), as soon as $\mathrm{Pr}(Y_2 > x_2|Y_1 > x_1)$ is non-decreasing in $x_1$ for all $x_2 \ge 0$. Finally, $Y_2$ is said to be Left-Tail-Decreasing in $Y_1$, written LTD($Y_2|Y_1$), as soon as $\mathrm{Pr}(Y_2 \le x_2|Y_1 \le x_1)$ is non-increasing in $x_1$ for all $x_2 > 0$. Both RTI($Y_2|Y_1$) and LTD($Y_2|Y_1$) properties are positive dependence properties, which imply association and positive quadrant dependence of $Y$; see, e.g., [19, Chapter 9] for more details on these different notions.

**Theorem 3** $\bar{F}_Y$ is MTP2.

The proof of Theorem 3 is long and technical, and hence postponed to Appendix A.

As is well known, the MTP2 property of $\bar{F}_Y$ entails that both RTI($Y_2|Y_1$) and RTI($Y_1|Y_2$) properties are true; see, e.g., Theorem 8.5 in [14]. However, as will be shown in Example 5 from Section 4, the LTD($Y_2|Y_1$) property is not always true. As the MTP2 property of $F_Y$ implies that both LTD($Y_2|Y_1$) and LTD($Y_1|Y_2$) properties are true, $F_Y$ is not always MTP2, again by Theorem 8.5 in [14].

### 3.4 Influence of the shock model parameters on the bivariate lifetime

We now study the influence of different parameters on the bivariate lifetime. To this end, two similar systems are considered ($S$ and $\bar{S}$, say), with identical parameters except from one. An upper bar is added to all quantities referring to the second system. (For instance, we use $\bar{\lambda}(w)$ for the second system). The bivariate lifetimes $Y$ and $\bar{Y}$ are next compared, using different stochastic orders whose definitions are now recalled.

Let $X = (X_1, X_2)$ and $Y = (Y_1, Y_2)$ be two bivariate random vectors. Then $X$ is said to be smaller than $Y$ with respect to the Upper Orthant (UO) order
(denoted by $X \le_{uo} Y$) if

$$ \forall (x_1, x_2) \in \mathbb{R}_+^2 \quad \bar{F}_X(\mathbf{x}) \le \bar{F}_Y(\mathbf{x}), $$

or equivalently (in the two-dimensional case) if

$$ \forall (x_1, x_2) \in \mathbb{R}_+^2 \quad F_X(\mathbf{x}) \le F_Y(\mathbf{x}). $$

Also, $X$ is said to be smaller than $Y$ in the sense of the Weak Hazard Rate order ($X \le_{WHR} Y$) if $\bar{F}_Y(\mathbf{x})/\bar{F}_X(\mathbf{x})$ is non-decreasing with respect of $\mathbf{x} \in \{\mathbf{y} \in \mathbb{R}_+^2 : \bar{F}_Y(\mathbf{y}) > 0\}$. The random variable $X$ is said to be weaker than $Y$ in the Hazard Rate sense ($X \le_{HR} Y$) if

$$ \forall_{\mathbf{x}, \mathbf{y} \in \mathbb{R}_+^2} \quad \bar{F}_Y(\mathbf{y}) \bar{F}_X(\mathbf{x}) \le \bar{F}_Y(\mathbf{x} \lor \mathbf{y}) \bar{F}_X(\mathbf{x} \land \mathbf{y}). $$

The HR order implies the WHR order.

Finally, $X$ is said to be smaller than $Y$ in the sense of the bivariate Laplace transform order (denoted by $X \le_L Y$) if

$$ \forall_{\mathbf{s}=(s_1, s_2) \in \mathbb{R}_+^2} \quad \mathcal{L}_X(\mathbf{s}) = E(e^{-s_1 X_1 - s_2 X_2}) \ge \mathcal{L}_Y(\mathbf{s}) = E(e^{-s_1 Y_1 - s_2 Y_2}); $$

see [31]. Note in passing that [8] uses the reverse inequality.

The following result now shows that, as expected, the more frequent the shocks, the shorter the bivariate lifetime $Y$.

**Proposition 4** Let us consider two systems $S$ and $\bar{S}$ with the same parameters except from the intensity of the non-homogeneous Poisson process. Assume that $\lambda(w) \le \bar{\lambda}(w)$ for all $w \ge 0$. Then $\bar{Y}$ is smaller than $Y$ in the sense of the Hazard Rate order ($\bar{Y} \le_{HR} Y$).

**Proof.** We first show that $\bar{Y} \le_{WHR} Y$. Note that $\bar{F}_{\bar{Y}}(\mathbf{x})/\bar{F}_{\bar{Y}}(\mathbf{x}) = \bar{F}_{\bar{Y}}(\mathbf{x})$, where $\tilde{Y}$ stands for the bivariate lifetime of a similar system to $S$ and $\bar{S}$ with identical parameters, except from the intrinsic failure rates with $\tilde{h} = 0$ and the intensity of shocks, with $\tilde{\lambda} = \bar{\lambda} - \lambda$. As $\bar{F}_{\bar{Y}}(\mathbf{x})$ is non-increasing in $\mathbf{x}$, it is clear that $\bar{F}_{\bar{Y}}(\mathbf{x})/\bar{F}_{\bar{Y}}(\mathbf{x})$ is non-decreasing in $\mathbf{x}$ and $\bar{Y} \le_{WHR} Y$. Now, remembering that $\bar{F}_{\bar{Y}}$ is MTP2 from Theorem 3, we conclude that $\bar{Y} \le_{HR} Y$ from Theorem 6.D.1 in [32].

**Remark 4** As the HR order implies the UO order [32, Eq. (6.G.10)] and as the UO order implies the Laplace transform order (consequence of Theorem 6.G.14 in [32] in the bivariate case), we deduce from the previous result that if $\lambda(w) \le \bar{\lambda}(w)$ for all $w \ge 0$, then $\bar{Y} \le_{uo} Y$ and $\bar{Y} \le_L Y$ as well. Based on [12], we also conclude that

$$ \forall_{\mathbf{x} \in \mathbb{R}_+^2} \quad m_{\bar{Y}}(\mathbf{x}) = E(\bar{Y} - \mathbf{x}|\bar{Y} > \mathbf{x}) \le m_Y(\mathbf{x}) = E(Y - \mathbf{x}|Y > \mathbf{x}), $$
where $m_{\bar{Y}}(\mathbf{x})$ and $m_Y(\mathbf{x})$ are the bivariate mean residual lifetimes of $\bar{Y}$ and $Y$, respectively.

The next result shows that the larger the fatality of shocks, the shorter the bivariate lifetime $Y$.

**Proposition 5** For each $t \ge 0$, let $U(t)$ with distribution $\text{Pr}\{U(t) = (i, j)\} = p_{ij}(t)$ for all $i, j \in \{0, 1\}$ and let $\bar{U}(t)$ be defined in the same way with respect to the family $\{\bar{p}_{ij}(t)\}_{i,j \in \{0, 1\}}$. Let us assume that $U(t) \le_{UO} \bar{U}(t)$ for all $t \ge 0$. Then $Y \le_{UO} \bar{Y}$.

**Proof.** As $U(t) \le_{UO} \bar{U}(t)$, we know that for all $i, j \in \{0, 1\}$,

$$ \text{Pr}\{U(t) \ge (i, j)\} \le \text{Pr}\{\bar{U}(t) \ge (i, j)\}. $$

Taking $(i, j) = (1, 1)$ (resp. $(1, 0)$, $(0, 1)$) this yields $p_{11}(t) \le \bar{p}_{11}(t)$ (resp. $p_{11}(t) + p_{10}(t) = q_1(t) \le \bar{q}_1(t)$, $p_{11}(t) + p_{01}(t) = q_2(t) \le \bar{q}_2(t)$). Therefore, one has $\bar{F}_Y(\mathbf{x}) \le \bar{F}_{\bar{Y}}(\mathbf{x})$ for all $\mathbf{x} \in \mathbb{R}_+^2$. Hence $Y \le_{UO} \bar{Y}$. ■

Example 7 from Section 4 shows that even if $U(t) \le_{WHR} \bar{U}(t)$ for all $t \ge 0$, then $Y$ is not always smaller than $\bar{Y}$ in the WHR sense.

The next and last result shows that the larger the increments of hazard rates are at jump times, the shorter the bivariate lifetime $Y$ is.

**Proposition 6** Assume $V$ to be smaller than $\bar{V}$ in the sense of the bivariate Laplace transform order ($V \le_{\mathcal{L}} \bar{V}$). Then $\bar{Y} \le_{UO} Y$.

**Proof.** Based on $V \le_{\mathcal{L}} \bar{V}$, we have $\tilde{\mu}(x_1, x_2) \ge \tilde{\mu}(x_1, x_2)$ and $\tilde{\mu}_i(x_i) \ge \tilde{\mu}_i(x_i)$ for all $x_i \ge 0$ and $i = 1, 2$. Inserting these inequalities into Eq. (3) easily provides the result. ■

Let us assume that $V$ and $\bar{V}$ have identical marginal distributions and remark that $Y$ and $\bar{Y}$ then share the same property. Under this assumption and in the present bivariate case, the concordance order is known to boil down to the UO order; see [28] for more details on this notion. Also, the concordance order measures the dependence between the margins of a random vector (with given marginal distributions). In that case, the previous result means that the larger the increments of hazard rates are at jump times (with respect to $\le_{\mathcal{L}}$), the less dependent the marginal lifetimes are. As the Laplace transform order is implied by the concordance (or UO) order, we also derive that the more dependent the increments of hazard rates are ($V \le_{\mathcal{C}} \bar{V}$), the less dependent the marginal lifetimes are ($\bar{Y} \le_{\mathcal{C}} Y$).

**Remark 5** As already noticed in Remark 3, $\bar{F}(t, t)$ corresponds to the reliability of a series system at time $t$. Then, all of the above comparison results entail similar ones for the series system: if $\lambda(w)$ increases, then the lifetime of the
<figure><img src="image_4.png" /><figcaption>Figure 1: Example 3, Joint survival function of Y</figcaption></figure>

series system decreases with respect to the univariate HR order; if $U(t)$ increases with respect to the UO order, or if $V$ decreases with respect to the Laplace order, then the lifetime of the series system increases with respect to the usual stochastic (or UO) order.

# 4 Numerical experiments

In all the following experiments, we take $h_i = 0$ for $i = 1, 2$. In Examples 4, 6, 7 and 8, $V = (V^{(1)}, V^{(2)})$ has a Marshall–Olkin distribution with parameters $(\lambda_1, \lambda_2, \lambda_{12})$ as in Example 2. In the other examples, we take $V^{(i)} = U^{(i)} + U^{(3)}$, $i = 1, 2$, where $U^{(1)}, U^{(2)}$ and $U^{(3)}$ are independent and $U^{(i)}$ is gamma or exponentially distributed for $i = 1, 2, 3$, where the parametrization of the gamma distribution is such that

$$f(x) = \frac{a^b}{\Gamma(a)} x^{a-1} e^{-bx} \mathbf{1}_{\{x>0\}}.$$

**Example 3** The parameters are: $\lambda(x) = 2x$, $p_{11}(x) = p_{01}(x) = p_{10}(x) = e^{-x}/4$ and the $U^{(i)}$'s ($i = 1, 2, 3$) are gamma distributed with respective parameters (1, 1), (2, 1) and (3, 1). The joint survival function of Y is displayed in Figure 1. As expected (see Remark 2), the joint survival function $\bar{F}_Y(s, t)$ is continuous on $\mathbb{R}_+^2$ but it is not differentiable at points $(s, s)$ of the main diagonal.

**Example 4** The parameters are $p_{01}(x) = p_{10}(x) = p_{11}(x) = 1/4$, $\lambda(x) = 2$ and $V$ is Marshall–Olkin distributed with parameters (1, 1, 1). The assumptions of Proposition 2 are verified. We take $s_2 = 0.10 < t_2 = 0.14$ and we set

$$D(s_1, t_1) = \bar{F}_Y(s_1, t_1)\bar{F}_Y(s_2, t_2) - \bar{F}_Y(s_1 + s_2, t_1 + t_2).$$
<figure><img src="image_5.png" /><figcaption>Figure 2: Example 4, BNBU property, case $s_1 < t_1$ (left) and without condition (right)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 3: Example 5, MTP2 property</figcaption></figure>

The function $D(s_1, t_1)$ is plotted in Figure 2 in the case where $s_1 < t_1$ (left) and also removing this condition (right). We observe that, as expected, $D(s_1, t_1)$ remains non-negative on the left plot, which shows that Y is BNBU in the sense of Proposition 2. However, this is no longer true when condition $s_1 < t_1$ is removed. Hence, Y is not BNBU in the stronger sense provided by (13).

**Example 5** The parameters are: $\lambda(x) = e^x$, $p_{11}(x) = p_{01}(x) = p_{10}(x) = e^{-x}/4$. Also, the $U^{(i)}$'s ($i = 1, 2, 3$) are exponentially distributed with respective means 1, 5 and 6. Taking $\mathbf{y} = (0.1, 0.675)$, the difference

$$d(\mathbf{x}) = \bar{F}_Y(\mathbf{y})\bar{F}_Y(\mathbf{x}) - \bar{F}_Y(\mathbf{x} \vee \mathbf{y})\bar{F}_Y(\mathbf{x} \wedge \mathbf{y})$$

is plotted in Figure 3. We observe that it remains non-positive, which is coherent with the MTP2 property from Theorem 3. Figure 4 (left) shows the right tail $RT(x_1) = \bar{F}_Y(x_1, x_2)/\bar{F}_{Y_1}(x_1)$ with respect to $x_1$ for various values of $x_2$. Whatever $x_2$ is, we observe that the right tail is always increasing ($RTI(Y_2|Y_1)$ property), which is coherent with the fact that $\bar{F}_Y$ is MTP2 (see the lines follow-
<figure><img src="image_7.png" /><figcaption>Figure 4: Example 5, RTI and not LTD property</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5: Example 6, $d(\mathbf{x})$ is non positive - influence of intensity</figcaption></figure>

ing Theorem 3). However, the left tail $LT(x_2) = F_Y(x_1, x_2)/F_{Y_2}(x_2)$ is plotted in Figure 4 (right) for $x_1 = 1.2$ and we observe that it is not monotonic. The LTD $(Y_2|Y_1)$ property is consequently not true and the MTP2 property cannot hold for $F_Y$ in a general setting.

**Example 6** The parameters are $p_{01}(x) = p_{10}(x) = p_{11}(x) = e^{-x}/4$ and $V$ is Marshall-Olkin distributed with parameters $(1, 2, 1)$. We consider $\lambda(x) = x \le \bar{\lambda}(x) = 2x$, so that the assumptions of Proposition 4 hold true. Taking $\mathbf{y} = (0.36, 0.2)$, the difference

$$d(\mathbf{x}) = \bar{F}_Y(\mathbf{y})\bar{F}_{\bar{Y}}(\mathbf{x}) - \bar{F}_Y(\mathbf{x} \vee \mathbf{y})\bar{F}_{\bar{Y}}(\mathbf{x} \wedge \mathbf{y})$$

is always non-positive (see Figure 5). This means that $\bar{Y}$ is weaker than $Y$ in the HR sense, in accordance with Proposition 4.
<figure><img src="image_9.png" /><figcaption>Figure 6: Example 7 - Influence of fatality of shocks</figcaption></figure>

**Example 7** This example illustrates the influence of the fatality of shocks on the bivariate lifetime as in Proposition 5. Let $Y$ and $\bar{Y}$ be the lifetimes corresponding to $p_{11}(x) = p_{10}(x) = p_{01}(x) = 1/6$ and $\bar{p}_{11}(x) = \bar{p}_{10}(x) = \bar{p}_{01}(x) = 1/6$. The $V^{(i)}$'s have a Marshall-Olkin distribution with parameters (4, 4, 1) and the intensity is $\lambda(x) = 1$. We have $U \le_{UO} \bar{U}$ and, as expected from Proposition 5, the difference $D = \bar{F}_{\bar{Y}} - \bar{F}_{Y}$ is always non-negative; see Figure 6 (left). Accordingly, the bivariate lifetime $Y$ is larger than $\bar{Y}$ in the sense of the UO order. However, taking $x_1 = 0.7$, Figure 6 (right) shows that $r(x_1, x_2) = \bar{F}_{\bar{Y}}(x_1, x_2)/\bar{F}_{Y}(x_1, x_2)$ is not monotone with respect to $x_2$. It implies that $Y$ and $\bar{Y}$ are not comparable with respect to the WHR order (even though $U \le_{WHR} \bar{U}$ is valid).

**Example 8** This example illustrates the influence of the size of $V$, as measured by the Laplace transform order; see Proposition 6. All parameters are the same as in Example 4 except from a Marshall-Olkin distribution with parameters (1, 1, 1) for $V$ and with parameters (1, 1, 2) for $\bar{V}$. It is then easy to check that $\bar{V} \le_{UO} V$, and hence that $\bar{V} \le_{L} V$. We observe in Figure 7 (left) that, as expected, the difference $D = \bar{F}_{\bar{Y}} - \bar{F}_{Y}$ is always non-negative, so that $Y$ is smaller than $\bar{Y}$, in the sense of the upper orthant order. Also, Figure 7 (right) plots $r(x_1, x_2) = \bar{F}_{\bar{Y}}(x_1, x_2)/\bar{F}_{Y}(x_1, x_2)$ with respect to $x_2$ for $x_1 = 0.725$. As it is not monotonic, $Y$ and $\bar{Y}$ are not comparable with respect to the WHR order.

**Example 9** As a last example, we look at the influence of the dependence between $V^{(1)}$ and $V^{(2)}$ on the bivariate lifetime, when the marginal distributions of $V$ are fixed. We take $p_{1,1}(x) = 1$, $\lambda(x) = 2x$. The $U^{(i)}$'s ($i = 1, 2, 3$) are gamma distributed with respective parameters (3, 1), (3, 1) and (0, 1), and the parameters for the $\bar{U}^{(i)}$'s are (0, 1), (0, 1) and (3, 1), which leads to identical marginal
<figure><img src="image_10.png" /><figcaption>Figure 7: Example 8 - Influence of the size of V</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 8: Example 8 - Influence of the dependence between $V^{(1)}$ and $V^{(2)}$</figcaption></figure>

distributions for $V$ and $\bar{V}$. It is easy to check that $V \le_{\mathcal{L}} \bar{V}$. As expected, the difference $D = \bar{F}_{\bar{Y}} - \bar{F}_{Y}$ is always non-negative; see the left panel of Figure 8. However, for $\mathbf{y} = (0.7, 0.7)$, the difference

$$d(\mathbf{x}) = \bar{F}_{Y}(\mathbf{y})\bar{F}_{\bar{Y}}(\mathbf{x}) - \bar{F}_{Y}(\mathbf{x} \vee \mathbf{y})\bar{F}_{\bar{Y}}(\mathbf{x} \wedge \mathbf{y})$$

is not always of the same sign, which shows that $Y$ and $\bar{Y}$ are not comparable with respect to the HR order (nor to the WHR order, due to the MTP2 property of $\bar{F}_{Y}$).
# 5 Estimation procedure

Assuming a parametric framework, we briefly suggest here some possible approaches for estimating the model parameters, according to the observation scheme. The estimation procedures are illustrated in Example 2, where we recall that the bivariate failure rate increments have a Marshall-Olkin distribution. For each estimation method, 500 independent sets of $n_{MC}$ independent bivariate lifetimes are generated. The parameters are estimated for each set of $n_{MC}$ bivariate data. This provides 500 estimates for each parameter, from which we report the mean, the standard deviation and the 2.5% and 97.25% quantiles ($q_{0.025}$ and $q_{0.975}$).

## 5.1 Observation of failure times

Based on Remark 2, it is possible to compute the density function $f_Y(s, t)$ of the distribution of $Y$ with respect to $ds \, dt + ds \, \delta_s(dt)$ (at least numerically). It is hence possible to use the standard maximum likelihood (ML) method. For $i = 1, 2$, the marginal density $f_{Y_i}(s) = -\partial \bar{F}_{Y_i}(s)/\partial s$ of the distribution of $Y_i$ with respect to Lebesgue measure is also available and can be used to compute ML estimators for the marginal parameters. We consequently suggest the following two-step procedure:

1. Estimation of $\lambda$ and of the marginal parameters from the marginal data. In Example 2, this provides $(\hat{\lambda}^{(i)}, \hat{\gamma}_i, \hat{q}_i)$ for $i = 1, 2$ where $\hat{\lambda}^{(i)}$ for $i = 1, 2$ refer to the two different estimation results based on the two marginal data sets.

2. New estimation of $\lambda$ (starting from $(\hat{\lambda}^{(1)} + \hat{\lambda}^{(2)})/2$ as an initial guess) and of the remaining parameters from the bivariate data. In Example 2, this provides $(\hat{\lambda}, \hat{\gamma}, \hat{p}_{11})$.

Results are provided in Table 1 for $n_{MC} = 500$ and one parameter set. Note that this parameter set provides a probability $Q_{11}$ that the failures of the two components are simultaneous (and consequently due to a shock) of about 10.5%. The probability that the failure of one component ($Q_1$ and $Q_2$) is due to a shock is about 22% for each of them. As can be seen in Table 1, the ML method does not provide reliable results, even for the marginal parameters. This is not very surprising based on the large number of parameters to estimate. Even if enlarging the number $n_{MC}$ of observations provides better results (see Table 2 where $n_{MC} = 1000$), it seems that this method requires too large a sample size for the results to be reliable in practice.

However, in a real life application, expert advice (or the data themselves) can provide additional information that can be of great help for estimation
<table><caption>Table 1: Estimation results for ML estimates based on failure data with $n_{MC} = 500$; computing time ≈ 399 c.p.u. time; $Q_{11} \approx 10.5\%$, $Q_1 \approx 22.5\%$, $Q_2 \approx 22.5\%$.</caption><thead><tr><th>Parameter</th><th>True value</th><th>Mean (std)</th><th>[q<sub>0.025</sub>, q<sub>0.075</sub>]</th></tr></thead><tbody><tr><td>q<sub>1</sub></td><td>0.8</td><td>0.8011(0.1889)</td><td>[0.3910, 1]</td></tr><tr><td>γ<sub>1</sub></td><td>0.15</td><td>0.3144(0.4314)</td><td>[0.0500, 1.9999]</td></tr><tr><td>q<sub>2</sub></td><td>0.8</td><td>0.8109(0.1803)</td><td>[0.4013, 1]</td></tr><tr><td>γ<sub>2</sub></td><td>0.15</td><td>0.2939(0.4037)</td><td>[0.0500, 1.9999]</td></tr><tr><td>λ</td><td>0.4</td><td>0.3964 (0.0566)</td><td>[0.2808, 0.5267]</td></tr><tr><td>p<sub>11</sub></td><td>0.7</td><td>0.6452(0.2732)</td><td>[0.0773, 1]</td></tr><tr><td>γ</td><td>0.25</td><td>0.4908(0.6086)</td><td>[0.0825, 2.3155]</td></tr></tbody></table>

<table><caption>Table 2: Estimation results for ML estimates based on failure data with $n_{MC} = 1000$, computing time ≈ 708 c.p.u. time; $Q_{11} \approx 10.5\%$, $Q_1 \approx 22.5\%$, $Q_2 \approx 22.5\%$.</caption><thead><tr><th>Parameter</th><th>True value</th><th>Mean (std)</th><th>[q<sub>0.025</sub>, q<sub>0.075</sub>]</th></tr></thead><tbody><tr><td>q<sub>1</sub></td><td>0.8</td><td>0.8115(0.1457)</td><td>[0.4954, 1]</td></tr><tr><td>γ<sub>1</sub></td><td>0.15</td><td>0.2149(0.2335)</td><td>[0.0500, 0.7759]</td></tr><tr><td>q<sub>2</sub></td><td>0.8</td><td>0.8108(0.1460)</td><td>[0.5161, 1]</td></tr><tr><td>γ<sub>2</sub></td><td>0.15</td><td>0.2111(0.2272)</td><td>[0.0500, 0.7528]</td></tr><tr><td>λ</td><td>0.4</td><td>0.3983 (0.0324)</td><td>[0.3316, 0.4719]</td></tr><tr><td>p<sub>11</sub></td><td>0.7</td><td>0.6880(0.2522)</td><td>[0.0773, 1]</td></tr><tr><td>γ</td><td>0.25</td><td>0.3371(0.3295)</td><td>[0.0931, 1.3477]</td></tr></tbody></table>

purposes. As an example, it can be considered that simultaneous failures are impossible or on the contrary that failures induced by shocks always induce simultaneous failures (when both components are still alive). Also, it may be possible that the two components could be considered as identical. All these situations allow to reduce the number of parameters to be estimated and will consequently improve the estimation results. Another situation that typically fits most examples we have in mind corresponds to the case where the times of shocks are observed. Contrary to the previously described situations, this one requires a specific estimation procedure which we now describe.

## 5.2 Observation of both shocks and failure times

We only consider the case where $h_1 = h_2 = 0$, which means that only the influence of shocks is estimated. Here, the shock times are observed, which allows to use the standard ML method for estimating the Poisson process parameters. This first step is not developed in the following, as it is standard. In the second step, the other parameters are estimated from the failure data, based on a conditional likelihood function given the shock times.

Starting again from (4), one easily obtains that the conditional survival
function $\bar{F}_{Y|\{T_n\}}(\cdot, \cdot|\{t_n\})$ of Y given $\{T_n = t_n, \forall n \in \mathbb{N}^*\}$ is given by

$$
\begin{align}
\bar{F}_{Y|\{T_n\}}(s, t|\{t_n\}) & = \prod_{i=1}^{n_t} p_{11}(t_i) \prod_{i=n_t+1}^{n_s} q_1(t_i) \mathrm{E}\{e^{-\sum_{i=n_t+1}^{n_s} V_i^{(1)}(s-t_i)}\} \mathrm{E}\left[e^{-\sum_{i=1}^{n_t} \{V_i^{(1)}(s-t_i)+V_i^{(2)}(t-t_i)\}}\right] \nonumber \\
& = \prod_{i=1}^{n_t} \{p_{11}(t_i)\tilde{\mu}(s-t_i, t-t_i)\} \prod_{i=n_t+1}^{n_s} \{q_1(t_i)\tilde{\mu}_1(s-t_i)\} \tag{14}
\end{align}
$$

for $s \ge t$ (with a similar expression for $s < t$), where $n_s$ (resp. $n_t$) stands for the observation of $N_s$ (resp. $N_t$). Also:

$$
\bar{F}_{Y_k|\{T_n\}}(s|\{t_n\}) = \prod_{i=1}^{n_s} \{q_k(t_i)\tilde{\mu}_k(s-t_i)\}
$$

for $k = 1, 2$. One easily deduces that the conditional distribution of $Y_k$ given $\{T_n = t_n : n = 1, 2, \dots\}$ has the following density with respect to $ds + \sum_{j=1}^{+\infty} \delta_{t_j}(ds)$:

$$
f_{Y_k|\{T_n\}}(s|\{t_n\}) = \begin{cases} -\prod_{i=1}^{n_s} \{q_k(t_i) \tilde{\mu}_k'(s-t_i)\} & \text{if } s \notin \{t_n\}, \\ (1-q_k(t_n)) \prod_{i=1}^{n-1} \{q_k(t_i) \tilde{\mu}_1(t_n-t_i)\} & \text{if } s = t_n \text{ with } n \ge 1, \end{cases}
$$

where $\tilde{\mu}_k'$ stands for the derivative of $\tilde{\mu}_k$. This enables one to write down the conditional (log)-likelihood function for the marginal data.

Now, setting $f_Y(s, t|\{t_n\})$ to be the bivariate density of the conditional distribution of Y given $\{T_n = t_n : n \ge 1\}$ with respect to

$$
ds \, dt + \sum_{j=1}^{+\infty} ds \, \delta_{t_j}(dt) + \sum_{j=1}^{+\infty} \delta_{t_j}(ds) \, dt + \sum_{j=1}^{+\infty} \sum_{k=1}^{+\infty} \delta_{(t_j, t_k)}(ds, dt),
$$

we first have

$$
f(s, t | \{t_n\}) = \frac{\partial^2}{\partial s \partial t} \bar{F}_{Y|(T_n)}(s, t | \{t_n\})
$$

for all $s, t \in \mathbb{R}_+ \setminus \{t_n\}$.

Starting again from (14), we observe that

$$
\begin{align*}
& \Pr(Y_1 > s, Y_2 = t_k | \{T_n = t_n\}) \\
&= \bar{F}_{Y|(T_n)}(s, t_k^- | \{t_n\}) - \bar{F}_{Y|(T_n)}(s, t_k | \{t_n\}) \\
&= \begin{cases}
p_{10}(t_k) \tilde{\mu}_1 (s-t_k) \prod_{i=1}^{k-1} \{p_{11}(t_i) \tilde{\mu} (s-t_i, t_k-t_i)\} & \text{if } s \ge t_k, \\
\times \prod_{i=k+1}^{n_s} \{q_1(t_i) \tilde{\mu}_1 (s-t_i)\} & \\
\{1-q_2(t_k)\} \prod_{i=1}^{n_s} \{p_{11}(t_i) \tilde{\mu} (s-t_i, t_k-t_i)\} & \\
\times \prod_{i=n_s+1}^{k-1} \{q_2(t_i) \tilde{\mu}_2 (t_k-t_i)\} & \text{if } s < t_k.
\end{cases}
\end{align*}
$$
<table><caption>Table 3: Estimation results for ML estimates based on failure data and shock times for $n_{MC} = 100$; $\lambda = 0.4$; 399 c.p.u. time; $Q_{11} \approx 10.5\%$, $Q_1 \approx 22.5\%$, $Q_2 \approx 22.5\%$</caption><thead><tr><th>Parameter</th><th>True value</th><th>Mean (std)</th><th>[q<sub>0.025</sub>, q<sub>0.075</sub>]</th></tr></thead><tbody><tr><td>q<sub>1</sub></td><td>0.8</td><td>0.8003(0.0363)</td><td>[0.7248, 0.8696]</td></tr><tr><td>γ<sub>1</sub></td><td>0.15</td><td>0.1534(0.0267)</td><td>[0.1056, 0.2108]</td></tr><tr><td>q<sub>2</sub></td><td>0.8</td><td>0.7992(0.0354)</td><td>[0.7350, 0.8704]</td></tr><tr><td>γ<sub>2</sub></td><td>0.15</td><td>0.1533(0.0276)</td><td>[0.1065, 0.2091]</td></tr><tr><td>p<sub>11</sub></td><td>0.7</td><td>0.7480(0.0503)</td><td>[0.6399, 0.8296]</td></tr><tr><td>γ</td><td>0.25</td><td>0.2523(0.0578)</td><td>[0.1529, 0.3624]</td></tr></tbody></table>

For $s \in \mathbb{R}_+ \setminus \{t_n\}$, we can get

$$f(s, t_k | \{t_n\}) = - \frac{\partial}{\partial s} \Pr(Y_1 > s, Y_2 = t_k | \{T_n = t_n\})$$

and a similar expression for $f(t_k, t)$ and $t \in \mathbb{R}_+ \setminus \{t_n\}$. Finally,

$$f(t_j, t_k | \{t_n\}) = \Pr(Y_1 > t_j^-, Y_2 = t_k | \{T_n = t_n\}) - \Pr(Y_1 > t_j, Y_2 = t_k | \{T_n = t_n\}) \\ = \begin{cases} p_{00}(t_j) \prod_{i=1}^{j-1} \{p_{11}(t_i)\tilde{\mu}(t_j - t_i, t_j - t_i)\} & \text{if } j = k, \\ p_{01}(t_j) (1 - q_2(t_k)) \tilde{\mu}_2(t_k - t_j) \prod_{i=1}^{j-1} \{p_{11}(t_i)\tilde{\mu}(t_j - t_i, t_k - t_i)\} \\ \times \prod_{i=j+1}^{k-1} \{q_2(t_i)\tilde{\mu}_2(t_k - t_i)\} & \text{if } j < k, \end{cases}$$

with a similar expression for $j > k$. It is then possible to write the conditional (log)-likelihood function for the bivariate data.

We now provide results for Example 2 with $n_{MC} = 100$ (and 500 replications). The estimation results are displayed in Table 3. Comparing to the estimation results from Table 1, one can see that these estimations are much better, even though the data size is 5 times smaller in the present case than in Table 1).

Though the procedure clearly deserves a more thorough study (with more numerical experiments), it seems possible to estimate the model parameters when shock times are observed. In case of unobserved shock times, a possibility might be to consider them as masked data and to use an EM algorithm as in, e.g., [10] or [17], where the authors use such a method for estimating the parameters of multivariate extensions either of the modified Sarhan-Balakrishnan or the Marshall-Olkin class of distributions. However, the development of an EM (or SEM) algorithm requires much more work in our present setting.

# 6 Concluding remarks

We proposed here a bivariate random shock model with competing failure modes. The model takes into account different kinds of dependence between
components arising from an external environment. The joint survival function of the bivariate lifetime is obtained explicitly. Conditions under which the bivariate life time has a BNBU property are provided. It is also numerically observed that under such conditions, the bivariate system lifetime is not BNBU in any stronger sense. Note that there exist many other multivariate aging properties in the literature, that have not been considered here and hence require further study. As an example, it would be of interest to investigate aging properties based on conditional distributions such as those developed in [1, 2, 33].

A strong positive dependence property (the MTP2 property) is also proved in the paper for the bivariate survival function, without any additional assumption. This entails that the RTI property is valid, too. The dual LTD property is however observed not to hold in a general setting. This shows that neither the MTP2 property for the bivariate cumulative distribution function nor the stochastic increasingness property of one lifetime with respect to the other [19] can hold in a similar general setting.

The influence of the shocks parameters is also studied. It is proved that the more frequent the shocks, the smaller the bivariate lifetime, in the sense of the (strong) hazard rate order. It is also showed that the smaller the fatality of shocks, the larger the bivariate lifetime. Finally, the larger the increments of hazard rates, the larger the bivariate lifetime, in the sense of the upper orthant order (but not in the sense of the weak hazard rate order). As a by-product, the more dependent the increments are, the less dependent the marginal lifetimes are.

The proposed model has thus many desirable aging and positive dependence properties to be used as a bivariate lifetime in reliability. A next step might be to propose and study preventive maintenance policies to enlarge the components lifetimes.

Finally, a parametric estimation procedure has been proposed, which seems to well behaved in case of observed shock times. As already mentioned, however, the procedure needs to be studied more thoroughly. Also, the case of non observed shock times requires further investigation such as (maybe?) the development of an EM (or SEM) algorithm.

# Appendix A

In this appendix, we show Theorem 3, namely we prove that $\bar{F}_Y$ is MTP2. Let us first write

$$ \bar{F}_Y(x_1, x_2) = e^{-H_1(x_1)-H_2(x_2)} \exp \{g(x_1, x_2)\} $$
with

$$g(x_1, x_2) = A(x_1, x_2) + B_1(x_1, x_2) + B_2(x_1, x_2)$$

where

$$A(x_1, x_2) = \int_{0}^{\min(x_1, x_2)} \{-1 + \tilde{\mu}(x_1 - w, x_2 - w)p_{11}(w)\} \lambda(w) dw,$$

$$B_1(x_1, x_2) = \mathbf{1}_{\{x_2 \le x_1\}} \int_{x_2}^{x_1} \{-1 + \tilde{\mu}_1(x_1 - w)q_1(w)\} \lambda(w) dw,$$

$$B_2(x_1, x_2) = \mathbf{1}_{\{x_2 > x_1\}} \int_{x_1}^{x_2} \{-1 + \tilde{\mu}_2(x_2 - w)q_2(w)\} \lambda(w) dw.$$

The function $(x_1, x_2) \mapsto e^{-H_1(x_1)-H_2(x_2)}$ is clearly MTP2, because both $H_1$ and $H_2$ are non-decreasing. As the product of MTP2 functions is MTP2 [16, Prop. 3.3], it is sufficient to show that

$$(x_1, x_2) \mapsto G(x_1, x_2) = \exp\{g(x_1, x_2)\}$$

is MTP2. Note that we have included the term $-\Lambda\{\max(x_1, x_2)\}$ in the function $g$ because the function $(x_1, x_2) \mapsto \exp\{-\max(x_1, x_2)\}$ is not MTP2 so that $(x_1, x_2) \mapsto \exp[-\Lambda\{\max(x_1, x_2)\}]$ will generally not be MTP2 either.

Let $x_1 \le y_1$ and $x_2 \ge y_2$ be fixed. One needs to prove that

$$g(x_1, x_2) + g(y_1, y_2) \le g(x_1, y_2) + g(y_1, x_2),$$

namely

$$A(x_1, x_2) + B_1(x_1, x_2) + B_2(x_1, x_2) + A(y_1, y_2) + B_1(y_1, y_2) + B_2(y_1, y_2) \\ \le A(x_1, y_2) + B_1(x_1, y_2) + B_2(x_1, y_2) + A(y_1, x_2) + B_1(y_1, x_2) + B_2(y_1, x_2). \tag{15}$$

Let us write all $A(z_1, z_2)$ terms in (15) as

$$A(z_1, z_2) = A_1(z_1, z_2) + A_2(z_1, z_2)$$

with

$$A_1(z_1, z_2) = \int_{0}^{\min(z_1, z_2)} \{-1 + \tilde{\mu}(z_1 - w, z_2 - w)p_{11}(w)\} \lambda(w) dw,$$

$$A_2(z_1, z_2) = \int_{\min(z_1, z_2)}^{\min(z_1, z_2)} \{-1 + \tilde{\mu}(z_1 - w, z_2 - w)p_{11}(w)\} \lambda(w) dw.$$

We will prove that

$$A_1(x_1, x_2) + A_1(y_1, y_2) \le A_1(x_1, y_2) + A_1(y_1, x_2) \tag{16}$$
and

$$ 
\begin{aligned}
& A_2(x_1, x_2) + B_1(x_1, x_2) + B_2(x_1, x_2) + A_2(y_1, y_2) + B_1(y_1, y_2) + B_2(y_1, y_2) \\
& \leq B_1(x_1, y_2) + B_2(x_1, y_2) + A_2(y_1, x_2) + B_1(y_1, x_2) + B_2(y_1, x_2)
\end{aligned} 
\quad (17) $$

(note that $A_2(x_1, y_2) = 0$), which will provide (15) by summation. To prove (16), we first note that it is equivalent to

$$ 
\begin{aligned}
& \int_{0}^{\min(x_1, y_2)} \{ \tilde{\mu}(x_1 - w, y_2 - w) + \tilde{\mu}(y_1 - w, x_2 - w) \\
& \qquad - \tilde{\mu}(y_1 - w, y_2 - w) - \tilde{\mu}(x_1 - w, x_2 - w) \} p_{11}(w) \lambda(w) dw \geq 0.
\end{aligned} 
\quad (18) $$

As $x_1 \leq y_1$ and $x_2 \geq y_2$, we have

$$ 
\begin{aligned}
& \tilde{\mu}(x_1 - w, y_2 - w) + \tilde{\mu}(y_1 - w, x_2 - w) - \tilde{\mu}(y_1 - w, y_2 - w) - \tilde{\mu}(x_1 - w, x_2 - w) \\
& = E[\{e^{-(x_1 - w)V^{(1)}} - e^{-(y_1 - w)V^{(1)}}\}\{e^{-(y_2 - w)V^{(2)}} - e^{-(x_2 - w)V^{(2)}}\}] \\
& \geq 0
\end{aligned} 
$$

for all $w \in [0, \min(x_1, y_2)]$. We deduce that inequality (18) is true, as well as inequality (16). We now come to inequality (17) and we distinguish between different cases.

**Case 1.** Assume that $y_1 \leq y_2$. Hence: $x_1 \leq y_1 \leq y_2 \leq x_2$. Noting that $A_2(x_1, x_2) = B_1(x_1, x_2) = B_1(y_1, y_2) = B_1(x_1, y_2) = B_1(y_1, x_2) = 0$, we must prove that

$$ B_2(x_1, x_2) + A_2(y_1, y_2) + B_2(y_1, y_2) \leq B_2(x_1, y_2) + A_2(y_1, x_2) + B_2(y_1, x_2). $$

This can be rewritten as

$$ 
\begin{aligned}
& \int_{x_1}^{x_2} \{-1 + \tilde{\mu}_2(x_2 - w)q_2(w)\} \lambda(w) dw + \int_{x_1}^{y_1} \{-1 + \tilde{\mu}(y_1 - w, y_2 - w) p_{11}(w)\} \lambda(w) dw \\
& + \int_{y_1}^{y_2} \{-1 + \tilde{\mu}_2(y_2 - w)q_2(w)\} \lambda(w) dw \\
& \leq \int_{x_1}^{y_2} \{-1 + \tilde{\mu}_2(y_2 - w)q_2(w)\} \lambda(w) dw + \int_{x_1}^{y_1} \{-1 + \tilde{\mu}(y_1 - w, x_2 - w) p_{11}(w)\} \lambda(w) dw \\
& + \int_{y_1}^{x_2} \{-1 + \tilde{\mu}_2(x_2 - w)q_2(w)\} \lambda(w) dw,
\end{aligned} 
$$
which may be simplified to

$$ 
\begin{aligned}
& \int_{x_1}^{y_1} \tilde{\mu}_2(x_2 - w)q_2(w) \lambda(w) dw + \int_{x_1}^{y_1} \tilde{\mu}(y_1 - w, y_2 - w) p_{11}(w) \lambda(w) dw \\
& \le \int_{x_1}^{y_1} \tilde{\mu}_2(y_2 - w)q_2(w) \lambda(w) dw + \int_{x_1}^{y_1} \tilde{\mu}(y_1 - w, x_2 - w) p_{11}(w) \lambda(w) dw.
\end{aligned} 
\quad (19) $$

Based on $y_2 \le x_2$, we have $\tilde{\mu}_2(y_2 - w) - \tilde{\mu}_2(x_2 - w) \ge 0$ for all $w \in [x_1, y_1]$. Using also that $q_2 \ge p_{11}$ and $x_1 \le y_1$, we get

$$ 
\begin{aligned}
& \int_{x_1}^{y_1} [\{\tilde{\mu}(y_1 - w, x_2 - w) - \tilde{\mu}(y_1 - w, y_2 - w)\} p_{11}(w) + \{\tilde{\mu}_2(y_2 - w) - \tilde{\mu}_2(x_2 - w)\} q_2(w)] \lambda(w) dw \\
& \ge \int_{x_1}^{y_1} \{\tilde{\mu}(y_1 - w, x_2 - w) - \tilde{\mu}(y_1 - w, y_2 - w) + \tilde{\mu}_2(y_2 - w) - \tilde{\mu}_2(x_2 - w)\} p_{11}(w) \lambda(w) dw \\
& = \int_{x_1}^{y_1} E \left[ \left\{ 1 - e^{-(y_1 - w)V^{(1)}} \right\} \left\{ e^{-(y_2 - w)V^{(2)}} - e^{-(x_2 - w)V^{(2)}} \right\} \right] p_{11}(w) \lambda(w) dw \\
& \ge 0.
\end{aligned} 
$$

Therefore, inequality (19) is true, which concludes this case.

**Case 2.** Assume that $x_1 \le y_2 \le y_1 \le x_2$. We must prove that

$$ B_2(x_1, x_2) + A_2(y_1, y_2) + B_1(y_1, y_2) \le B_2(x_1, y_2) + A_2(y_1, x_2) + B_2(y_1, x_2). $$

This may be written as

$$ 
\begin{aligned}
& \int_{x_1}^{x_2} \{-1 + \tilde{\mu}_2(x_2 - w)q_2(w)\} \lambda(w) dw + \int_{x_1}^{y_2} \{-1 + \tilde{\mu}(y_1 - w, y_2 - w) p_{11}(w)\} \lambda(w) dw \\
& + \int_{y_2}^{y_1} \{-1 + \tilde{\mu}_1(y_1 - w)q_1(w)\} \lambda(w) dw \\
& \le \int_{x_1}^{y_2} \{-1 + \tilde{\mu}_2(y_2 - w)q_2(w)\} \lambda(w) dw + \int_{x_1}^{y_1} \{-1 + \tilde{\mu}(y_1 - w, x_2 - w) p_{11}(w)\} \lambda(w) dw \\
& + \int_{y_1}^{x_2} \{-1 + \tilde{\mu}_2(x_2 - w)q_2(w)\} \lambda(w) dw
\end{aligned} 
$$

which may then be simplified to

$$ 
\begin{aligned}
& \int_{x_1}^{y_1} \tilde{\mu}_2(x_2 - w)q_2(w) \lambda(w) dw + \int_{x_1}^{y_2} \tilde{\mu}(y_1 - w, y_2 - w) p_{11}(w) \lambda(w) dw \\
& + \int_{y_2}^{y_1} \{-1 + \tilde{\mu}_1(y_1 - w)q_1(w)\} \lambda(w) dw \\
& \le \int_{x_1}^{y_2} \tilde{\mu}_2(y_2 - w)q_2(w) \lambda(w) dw + \int_{x_1}^{y_1} \tilde{\mu}(y_1 - w, x_2 - w) p_{11}(w) \lambda(w) dw.
\end{aligned} 
\quad (20) $$
Splitting the domain of integration into $(x_1, y_2)$ and $(y_2, y_1)$ and considering the two terms separately, we first have to prove that

$$ 
\begin{aligned}
& \int_{x_1}^{y_2} \tilde{\mu}_2(x_2 - w)q_2(w) \lambda(w) dw + \int_{x_1}^{y_2} \tilde{\mu}(y_1 - w, y_2 - w) p_{11}(w) \lambda(w) dw \\
& \leq \int_{x_1}^{y_2} \tilde{\mu}_2(y_2 - w)q_2(w) \lambda(w) dw + \int_{x_1}^{y_2} \tilde{\mu}(y_1 - w, x_2 - w) p_{11}(w) \lambda(w) dw,
\end{aligned} 
\quad (21) $$

which is true, using similar arguments as for (19). Looking at the second integral, we must next show that

$$ 
\begin{aligned}
& \int_{y_2}^{y_1} \tilde{\mu}_2(x_2 - w)q_2(w) \lambda(w) dw + \int_{y_2}^{y_1} \{-1 + \tilde{\mu}_1(y_1 - w)q_1(w)\} \lambda(w) dw \\
& \leq \int_{y_2}^{y_1} \tilde{\mu}(y_1 - w, x_2 - w) p_{11}(w) \lambda(w) dw.
\end{aligned} 
\quad (22) $$

Noting that

$$ 
\begin{aligned}
& \tilde{\mu}(y_1 - w, x_2 - w) p_{11}(w) - \tilde{\mu}_2(x_2 - w)q_2(w) + 1 - \tilde{\mu}_1(y_1 - w)q_1(w) \\
& = \{\tilde{\mu}(y_1 - w, x_2 - w) - \tilde{\mu}_2(x_2 - w) - \tilde{\mu}_1(y_1 - w)\} p_{11}(w) - \tilde{\mu}_2(x_2 - w)p_{01}(w) \\
& \quad + 1 - \tilde{\mu}_1(y_1 - w)p_{10}(w) \\
& = E[\{1 - e^{-(y_1 - w)V^{(1)}}\}\{1 - e^{-(x_2 - w)V^{(2)}}\} - 1]p_{11}(w) - \tilde{\mu}_2(x_2 - w)p_{01}(w) \\
& \quad + 1 - \tilde{\mu}_1(y_1 - w)p_{10}(w) \\
& \geq -p_{11}(w) - p_{01}(w) + 1 - p_{10}(w) \\
& \geq 0,
\end{aligned} 
$$

we can conclude that (22) is true, so that (20) is also true, upon summing (21) and (22).

**Case 3.** Assume that $x_1 \leq y_2 \leq x_2 \leq y_1$. We have to prove that

$$ B_2(x_1, x_2) + A_2(y_1, y_2) + B_1(y_1, y_2) \leq B_2(x_1, y_2) + A_2(y_1, x_2) + B_1(y_1, x_2). $$

After simplification, one gets

$$ 
\begin{aligned}
& \int_{x_1}^{x_2} \tilde{\mu}_2(x_2 - w)q_2(w) \lambda(w) dw + \int_{x_1}^{y_2} \tilde{\mu}(y_1 - w, y_2 - w) p_{11}(w) \lambda(w) dw \\
& + \int_{y_2}^{x_2} \{-1 + \tilde{\mu}_1(y_1 - w)q_1(w)\} \lambda(w) dw \\
& \leq \int_{x_1}^{y_2} \tilde{\mu}_2(y_2 - w)q_2(w) \lambda(w) dw + \int_{x_1}^{x_2} \tilde{\mu}(y_1 - w, x_2 - w) p_{11}(w) \lambda(w) dw.
\end{aligned} 
\quad (23) $$
Here again, we use $\int_{x_1}^{x_2} = \int_{x_1}^{y_2} + \int_{y_2}^{x_2}$ and we first look at the $\int_{x_1}^{y_2}$ terms. We obtain

$$ \begin{aligned} & \int_{x_1}^{y_2} \tilde{\mu}_2(x_2 - w)q_2(w) \lambda(w) dw + \int_{x_1}^{y_2} \tilde{\mu}(y_1 - w, y_2 - w) p_{11}(w) \lambda(w) dw \\ & \le \int_{x_1}^{y_2} \tilde{\mu}_2(y_2 - w)q_2(w) \lambda(w) dw + \int_{x_1}^{y_2} \tilde{\mu}(y_1 - w, x_2 - w) p_{11}(w) \lambda(w) dw, \end{aligned} \quad (24) $$

which is clear, using similar arguments as for (19). We now look at the $\int_{y_2}^{x_2}$ remaining terms in (23), which may be written as

$$ \begin{aligned} & \int_{y_2}^{x_2} \tilde{\mu}_2(x_2 - w)q_2(w) \lambda(w) dw + \int_{y_2}^{x_2} \{-1 + \tilde{\mu}_1(y_1 - w)q_1(w)\} \lambda(w) dw \\ & \le \int_{y_2}^{x_2} \tilde{\mu}(y_1 - w, x_2 - w) p_{11}(w) \lambda(w) dw \end{aligned} \quad (25) $$

and may be treated in a similar way as (22), which allows to conclude this case.

**Case 4.** Assume that $y_2 \le x_1 \le x_2 \le y_1$. We have to prove that

$$ A_2(x_1, x_2) + B_2(x_1, x_2) + B_1(y_1, y_2) \le B_1(x_1, y_2) + A_2(y_1, x_2) + B_1(y_1, x_2). $$

After simplification, this reduces to

$$ \begin{aligned} & \int_{y_2}^{x_1} \tilde{\mu}(x_1 - w, x_2 - w) p_{11}(w) \lambda(w) dw + \int_{x_1}^{x_2} \{-1 + \tilde{\mu}_2(x_2 - w)q_2(w)\} \lambda(w) dw \\ & + \int_{y_2}^{x_2} \tilde{\mu}_1(y_1 - w)q_1(w) \lambda(w) dw \\ & \le \int_{y_2}^{x_1} \tilde{\mu}_1(x_1 - w)q_1(w) \lambda(w) dw + \int_{y_2}^{x_2} \tilde{\mu}(y_1 - w, x_2 - w) p_{11}(w) \lambda(w) dw \end{aligned} \quad (26) $$

Here again, we use $\int_{y_2}^{x_2} = \int_{y_2}^{x_1} + \int_{x_1}^{x_2}$ and we first look at the $\int_{y_2}^{x_1}$ terms. This yields

$$ \begin{aligned} & \int_{y_2}^{x_1} \tilde{\mu}(x_1 - w, x_2 - w) p_{11}(w) \lambda(w) dw + \int_{y_2}^{x_1} \tilde{\mu}_1(y_1 - w)q_1(w) \lambda(w) dw \\ & \le \int_{y_2}^{x_1} \tilde{\mu}_1(x_1 - w)q_1(w) \lambda(w) dw + \int_{y_2}^{x_1} \tilde{\mu}(y_1 - w, x_2 - w) p_{11}(w) \lambda(w) dw \end{aligned} $$

which is clear, using similar arguments as for (19). We now look at the $\int_{x_1}^{x_2}$ terms.
remaining terms in (26), namely

$$ \int_{x_1}^{x_2} \{-1 + \tilde{\mu}_2(x_2 - w)q_2(w)\} \lambda(w) dw + \int_{x_1}^{x_2} \tilde{\mu}_1(y_1 - w)q_1(w) \lambda(w) dw \\ \le \int_{x_1}^{x_2} \tilde{\mu}(y_1 - w, x_2 - w) p_{11}(w) \lambda(w) dw $$

and may be handled in a similar way as (22).

**Case 5.** Assume that $y_2 \le x_2 \le x_1 \le y_1$. We have to prove that

$$ A_2(x_1, x_2) + B_1(x_1, x_2) + B_1(y_1, y_2) \le A_2(y_1, x_2) + B_1(x_1, y_2) + B_1(y_1, x_2) $$

or equivalently that

$$ \int_{y_2}^{x_2} \tilde{\mu}(x_1 - w, x_2 - w) p_{11}(w) \lambda(w) dw + \int_{y_2}^{x_2} \tilde{\mu}_1(y_1 - w) q_1(w) \lambda(w) dw \\ \le \int_{y_2}^{x_2} \tilde{\mu}(y_1 - w, x_2 - w) p_{11}(w) \lambda(w) dw + \int_{y_2}^{x_2} \tilde{\mu}_1(x_1 - w) q_1(w) \lambda(w) dw. $$

This may be proved in a similar way as (19).

**Case 6.** Assume that $y_2 \le x_1 \le y_1 \le x_2$. We have to prove that

$$ A_2(x_1, x_2) + B_2(x_1, x_2) + B_1(y_1, y_2) \le B_1(x_1, y_2) + B_2(y_1, x_2) + A_2(y_1, x_2) $$

which may be simplified into

$$ \int_{y_2}^{x_1} \tilde{\mu}(x_1 - w, x_2 - w) p_{11}(w) \lambda(w) dw + \int_{x_1}^{y_1} \{-1 + \tilde{\mu}_2(x_2 - w)q_2(w)\} \lambda(w) dw \\ + \int_{y_2}^{y_1} \tilde{\mu}_1(y_1 - w) q_1(w) \lambda(w) dw \\ \le \int_{y_2}^{x_1} \tilde{\mu}_1(x_1 - w) q_1(w) \lambda(w) dw + \int_{y_2}^{y_1} \tilde{\mu}(y_1 - w, x_2 - w) p_{11}(w) \lambda(w) dw. $$

This case may be proved in a similar way as Case 3., which achieves this proof.

## Acknowledgments

Hai Ha Pham would like to thank the Vietnam Institute for Advanced Study in Mathematics (VIASM) for its support. Both authors warmly thank the referees, Associate Editor and Editor for their thorough review of the paper, their insightful comments and suggestions, as well as their help in correcting linguistic and typographical mistakes. All this led to a better paper.
References

[1] J.P. Arias-Nicolás, F. Belzunce, O. Núñez-Barrera, and A. Suárez-Llorens, *A multivariate IFR notion based on the multivariate dispersive ordering*, Applied Stochastic Models in Business and Industry **25** (2009), 339-358.

[2] A. Arriaza, F. Belzunce, J. Mulero, and A. Suárez-Llorens, *On a new multivariate IFR aging notion based on the standard construction*, Applied Stochastic Models in Business and Industry **32** (2016), 292-306.

[3] C.L. Atwood, *The binomial failure rate common cause model*, Technometrics **28** (1986), 139-148.

[4] I. Bayramoglu and M. Ozkut, *The reliability of coherent systems subjected to Marshall-Olkin type shocks*, IEEE Transactions on Reliability **64** (2015), 435-443.

[5] J.H. Cha and I.T. Castro, *A stochastic failure model with dependent competing risks and its applications to condition-based maintenance*, Journal of Applied Probability **52** (2015), 558-573.

[6] J.H. Cha and J. Mi, *On a stochastic survival model for a system under randomly variable environment*, Methodology and Computing in Applied Probability **13** (2011), 549-561.

[7] E. Çınlar, *Probability and Stochastics*, Springer, New York, 2011.

[8] M. Denuit, *Laplace transform ordering of actuarial quantities*, Insurance: Mathematics and Economics **34** (2001), 83-102.

[9] M. Finkelstein and J.H. Cha, *Stochastic modelling for reliability*, Springer, 2013.

[10] M. Franco, D. Kundu, and J.-M. Vivo, *Multivariate extension of modified sarhan-balakrishnan bivariate distribution*, Journal of Statistical Planning and Inference **141** (2011), 3400-3412.

[11] E.W. Frees, J. Carrière, and E.A. Valdez, *Annuity valuation with dependent mortality*, The Journal of Risk and Insurance **63** (1996), 229-261.

[12] T. Hu, B. Khaledi, and M. Shaked, *Multivariate hazard rate orders*, Journal of Multivariate Analysis **84** (2003), 173-189.

[13] A. Jamalizadeh and D. Kundu, *Weighted Marshall-Olkin bivariate exponential distribution*, Statistics: A Journal of Theoretical and Applied Statistics **47** (2013), 917-928.

[14] H. Joe, *Dependence modeling with copulas*, CRC Press, 2014.
[15] H.W. Jones, *Common cause failures and ultra reliability*, 42nd International Conference on Environmental Systems, 2012, p. 3602.

[16] S. Karlin and Y. Rinott, *Classes of orderings of measures and related correlation inequalities. i. multivariate totally positive distributions*, Journal of Multivariate Analysis **10** (1980), 467–498.

[17] D. Karlis, *ML estimation for multivariate shock models via an algorithm*, Annals of the Institute of Statistical Mathematics **55** (2003), 817–830.

[18] D. Kundu and R.D. Gupta, *Bivariate generalized exponential distribution*, Journal of Multivariate Analysis **100** (2009), 581–593.

[19] C.D. Lai and M. Xie, *Stochastic aging and dependence for reliability*, Springer, New York, 2006.

[20] M. Lai and J. Yuan, *Periodic replacement model for a parallel system subject to independent and common cause shock failures*, Reliability Engineering & System Safety **31** (1991), 355–367.

[21] H. Li and S.H. Xu, *Stochastic bounds and dependence properties of survival times in a multicomponent shock model*, Journal of Multivariate Analysis **76** (2001), 63–89.

[22] X. Li and F. Pellerey, *Generalized Marshall-Olkin distributions and related bivariate aging properties*, Journal of Multivariate Analysis **102** (2011), 1399–1409.

[23] F. Mallor and J. Santos, *Classification of shock models in system reliability*, Monografías del Semin. Matem. García de Galdeano **27** (2003), 405–412.

[24] A.W. Marshall and I. Olkin, *A multivariate exponential distribution*, Journal of the American Statistical Association **62** (1967), 30–44.

[25] A.W. Marshall and M. Shaked, *Multivariate shock models for distributions with increasing hazard rate average*, The Annals of Probability **7** (1979), 343–358.

[26] —, *A class of multivariate new better than used distributions*, The Annals of Probability **10** (1982), 259–264.

[27] S. Mercier and H.H. Pham, *A random shock model with mixed effect, including competing soft and sudden failures, and dependence*, Methodology and Computing in Applied Probability **18** (2016), 377–400.

[28] A. Müller and D. Stoyan, *Comparison methods for stochastic models and risks*, Wiley, Chichester, 2002.
[29] T. Nakagawa, *Shock and damage models in reliability theory*, Springer, 2007.

[30] M. Rausand and A. Høyland, *System reliability theory. Models, statistical methods and applications*, 2nd ed., Hoboken, NJ: Wiley, 2004.

[31] M. Shaked, *Stochastic comparisons of multivariate random sums in the laplace transform order, with applications*, Statistics & Probability Letters **34** (2007), 1339–1344.

[32] M. Shaked and J.G. Shanthikumar, *Stochastic orders*, Springer, New York, 2007.

[33] __________, *Multivariate conditional hazard rate functions-an overview*, Applied Stochastic Models in Business and Industry **31** (2015), 285–296.

[34] N.D. Singpurwalla, *Survival in dynamic environments*, Statistical Science **10** (1995), 86–103.

[35] W.E. Vesely, *Estimating common cause failure probabilities in reliability and risk analysis: Marshall-Olkin specializations*, Nuclear systems reliability engineering and risk assessment (1977), 314–341.

[36] R. Wang, *A reliability model for multivariate exponential distributions*, Journal of Multivariate Analysis **98** (2007), 1033–1042.