The Marshall-Olkin distribution
Source:vignettes/The-Marshall-Olkin-distribution.Rmd
The-Marshall-Olkin-distribution.RmdDISCLAIMER: This article is a brief work of reference about Marshall–Olkin distributions. For a detailed mathematical treatment of this class, we refer to the references at the end of this document.
Marshall–Olkin distributions
Definition
A d-variate random vector \(\boldsymbol{\tau} = {( \tau_{1}, \ldots, \tau_{d} )}\) has a Marshall–Olkin (MO) distribution with parameters \(\lambda_{I} \geq {0}\), \(\emptyset \neq I \subseteq {\{ 1, \ldots, d \}}\), called shock arrival intensities, if and only if it has the multivariate survival function
\[ \bar{F} (\boldsymbol{t}) = \exp{\left \{ -\sum_{\emptyset \neq I \subseteq {\{ 1, \ldots, d \}}}{ \lambda_I \max_{i \in I}{ t_{i} } } \right \}}, \quad \boldsymbol{t} = {( t_{1}, \ldots, t_{d} )} \geq 0 , \]
and the rates for the marginal exponential distributions fulfill
\[ c_{i} := \sum_{I \ni i} \lambda_I > 0 , \quad \forall i \in {\{ 1, \ldots, d \}} . \]
The exogenous shock model
A d-variate Marshall–Olkin distribution has a natural stochastic model: consider independent exponentially distributed random variables \(E_{I} \sim \mathrm{Exp}{(\lambda_{I})}\), \(\emptyset \neq I \subseteq {\{ 1, \ldots, d \}}\), whose rates are the shock arrival intensities, and define
\[ \tau_{i} := \min{\Big\{ E_{I} \ : \ I \ni i \Big\}} , \quad i \in {\{ 1 , \ldots , d \}} . \]
The Arnold model
A d-variate Marshall–Olkin distribution has another stochastic model, which is strongly linked to continuous time, homogeneous Markovian processes. For this, consider the transformed parameters: the total shock arrival intensity
\[ \lambda := \sum_{\emptyset \neq I \subseteq {\{ 1, \ldots, d \}}}{ \lambda_{I} } \]
and the shock arrival probabilies
\[ p_{I} := \frac{ \lambda_{I} }{ \lambda }, \quad \emptyset \neq I \subseteq {\{ 1, \ldots, d \}} . \]
Let the following two sequences be independent:
- \(\epsilon_{1}, \epsilon_{2}, \ldots \sim \mathrm{Exp}{(\lambda)}\) iid.
- \(Y_{1}, Y_{2}, \ldots\) iid in \(\mathcal{P}{(\{1, \ldots, d\})} \setminus \emptyset\) with \(\mathbb{P}{( Y_{1} = I )} = p_{I}\), \(\emptyset \neq I \subseteq {\{ 1, \ldots , d\}}\).
Now, define \(\boldsymbol{\tau}\) by
\[ \tau_{i} := \sum_{j=1}^{\min{\{ k \in \mathbb{N} : i \in Y_{k} \}}}{ \epsilon_{j} } , \quad i \in {\{ 1, \ldots, d \}} . \]
The Markovian death-set model
A d-variate Marshall–Olkin distribution has yet another stochastic model through the associated death-set process, indicating which entities are dead, which is a Markov process with infinitesimal generator matrix \(Q = {( q_{I, J} : I, J \subseteq {\{ 1 , \ldots , d \}} )}\) defined by
\[ q_{I, J} = \begin{cases} -\sum_{I \subsetneq K \subseteq {\{ 1 , \ldots , d \}}} { q_{I, K} } , & I = J , \\ \sum_{L \subseteq I}{ \lambda_{L \cup {(J \setminus I )}} } , & I \subsetneq J , \\ 0 , & \text{else} . \end{cases} \]
Given a sample of \(Z\), simulated until reaching its absorbing state and associated waiting times \(W_{1}, W_{2}, \ldots\) and corresponding states \(Z_{W_{1}}, Z_{W_{1} + W_{2}}, \ldots\), define
\[ \tau_{i} := \sum_{j=1}^{ \min{\{ k : Z_{W_{1} + \cdots + W_{k}} \ni i \}} }{ W_{j} } , \quad i \in {\{ 1 , \ldots , d \}} . \]
Exchangeable Marshall–Olkin distributions
Exchangeability
A random vector is called exchangeable if a reordering of the components does not change its law. An equivalent definition is that the survival function, distribution function, or density (if it exists) does not change if the arguments are permutated.
Exchangeable shock arrival intensities
A Marshall–Olkin distribution is exchangeable if and only if
\[ {\lvert I \rvert} = {\lvert J \rvert} \quad \Rightarrow \quad \lambda_{I} = \lambda_{J} , \quad \forall \emptyset \neq I, J \subseteq {\{ 1 , \ldots , d \}} . \]
For this reason, \(\lambda_i := \lambda_{I}\), for \(\emptyset \neq I \subseteq {\{ 1, \ldots, d \}}\) with \(i = {\lvert I \rvert}\), \(i \in {\{ 1 , \ldots , d \}}\), are called (unscaled) exchangeable shock arrival intensities.
The Markovian death-counting model
A d-variate exchangeable Marshall–Olkin distribution has a stochastic model through the associated death-counting process \(Z^\ast\), indicating how many entities are dead, which is a Markov process with infinitesimal generator matrix \(Q^\ast = {( q^\ast_{i, j} : 0 \leq i, j \leq d )}\) defined by
\[ q^\ast_{i, j} = \begin{cases} -\sum_{l=1}^{d-i}{ q^\ast_{i, l} } , & i = j , \\ \binom{d-i}{j-i} \sum_{k=0}^{i}{ \binom{i}{k} \lambda_{k + {(j-i)}} } , & i < j , \\ 0 , & \text{else} . \end{cases} \]
The first row of this infinitesimal generator matrix consists of the (scaled) exchangeable shock-size arrival intensities \(\eta_{i} = \binom{d}{i} \lambda_{i}\), \(0 < i \leq d\), and the remaining rows can be calculated with the recursion
\[ q^\ast_{i+1, j+1} = \frac{d-j}{d-i} q^\ast_{i, j} + \frac{j+1-i}{d-i} q^\ast_{i, j+1} , \quad 0 \leq i < j < d . \]
Given a sample of \(Z^\ast\), simulated until reaching its absorbing state and associated waiting times \(W_{1}, W_{2}, \ldots\) and states \(Z^\ast_{W_{1}}, Z^\ast_{W_{1} + W_{2}}, \ldots\), and an independent uniform permutation \(\Pi\) on \({\{ 1 , \ldots, d \}}\), define
\[ \tau_{i} := \sum_{j = 1}^{\min{\{ k : Z^\ast_{W_{1} + \cdots + W_{k}} \geq \Pi{(i)} \}}}{ W_{j} } , \quad i \in {\{ 1 , \ldots , d \}} . \]
Extendible Marshall–Olkin distributions
Exchangeable sequence \(\tau_{1}, \tau_{2}, \ldots\) with finite margins of Marshall–Olkin type are identified with a so-called Bernstein function \(\psi\) and associated Lévy subordinator \(\Lambda\) such that
\[ \tau_{i} := \inf{\{ t > 0 : \Lambda_{t} \geq E_{i} \}} , \quad i \in {\{ 1 , \ldots , d \}} , \] for iid unit exponentially distributed random barrier variables \(E_{1}, \ldots, E_{d}\); the Bernstein function is the Lévy subordinator’s Laplace exponent, and it holds that
\[ \eta_{i} = \binom{d}{i} {(-1)}^{i-1} \Delta^{i} \psi{(d-i)}, \quad i \in {\{ 1 , \ldots , d \}}. \]
Bernstein functions
Bernstein functions are the Laplace exponents of (killed) Lévy subordinators. They are associated with a Lévy triplet \({(a, b, \nu)}\), with a killing rate \(a \geq 0\), a drift \(b \geq 0\), and a pure-jump Lévy measure \(\nu\) on \({(0, \infty)}\) with \(\int_{0}^{\infty} {(1 \wedge u)} \nu{(du)} < \infty\):
\[ \psi{(x)} = a 1_{\{ x > 0 \}} + b x + \int_{0}^{\infty}{ {(1 - e^{u x})} \nu{(du)}} , \quad x \geq 0 . \]
Complete Bernstein functions are also associate with a Stieltjes triples \({( a , b , \sigma )}\), with \(a, b\) from the Lévy triplet, and a Stieljtes measure \(\sigma\) on \({(0, \infty)}\) with \(\int_{0}^{\infty} {( 1 + u)}^{-1} \sigma{(du)} < \infty\):
\[ \psi{(x)} = a 1_{\{ x > 0 \}} + b x + \int_{0}^{\infty}{ \frac{x}{x + u} \sigma{(du)} } , \quad x \geq 0 . \]
Calculating the exchangeable shock-size arrival intensities
Numerically stable approximations of the exchangeable shock-size arrival intensities \(\eta_{i}\), \(i \in {\{ 1 , \ldots , d\}}\), can be derived from the Bernstein function’s integral representations.
-
For a Bernstein function with Lévy triplet \({(a, 0, 0 )}\), they can be calculated as
\[ \eta_{i} = \binom{d}{i} a 1_{\{ i = d \}} , \quad i \in {\{ 1 , \ldots , d \}} . \]
-
For a Bernstein function with Lévy triplet \({(0, b, 0 )}\), they can be calculated as
\[ \eta_{i} = \binom{d}{i} b 1_{\{ i = 1 \}} , \quad i \in {\{ 1 , \ldots , d \}} . \]
-
For a Bernstein function Lévy triplet \({(0, 0, \nu)}\), they can be calculated as
\[ \eta_{i} = \int_{0}^{\infty}{ \binom{d}{i} {(1 - e^{-u})}^{i} e^{-u {(d - i)}} \nu{(du)} } , \quad i \in {\{ 1 , \ldots , d \}} . \]
-
For a complete Bernstein function with Stieltjes triples \({(0, 0, \sigma )}\), they can be calculated as
\[ \eta_{i} = \int_{0}^{\infty}{ \binom{d}{i} u B{(1 + i, d-i + u)} \sigma{(du)} } , \quad i \in {\{ 1 , \ldots , d \}} . \]
Finally, linear combinations are associated with linear combinations of their Lévy (or Stieltjes) triples with the same coefficients. Similarly, their exchangeable shock-size arrival intensities can be calculated as a linear combination of the individual exchangeable shock-size arrival intensities with these coefficients.
The Lévy frailty model
A d-variate extendible Marshall–Olkin distribution with associated (killed) compound Poisson subordinator with drift \[ \Lambda_{t} = \infty \cdot 1_{\{ \epsilon \leq a t \}} + b t + \sum_{j = 1}^{N_{t}}{ X_{j} } , \quad t \geq 0 , \]
with a unit exponenential random variable \(\epsilon\) and compound Poisson process \(N_{t} = \max{\{ k : W_{1} + \cdots + W_{k} \leq t \}}\), define
\[ \tau_{i} := W_{1} + \cdots + W_{K_{i}} + \min{\left \{ W_{K_{i} + 1} , \frac{E_{i} - \Lambda_{W_{1} + \cdots + W_{K_{i}}}}{b} , \frac{\epsilon}{a} - {( W_{1} + \cdots + W_{K_{i}} )} \right \}} , \quad i \in {\{ 1 , \ldots , d \}} , \]
where
\[ K_{i} := \max{\left \{ k : \Lambda_{W_{1} + \cdots + W_{k}} < E_{i} \right \}} , \quad i \in {\{ 1 , \ldots , d \}} . \]
Parametric extendible families
Constant Bernstein function
Suppose \(a \geq 0\) and consider the constant Bernstein function
\[ \psi{(x)} = a 1_{\{ x > 0 \}} , \quad x \geq 0 . \]
Linear Bernstein function
Suppose \(b \geq 0\) and consider the linear Bernstein function
\[ \psi{(x)} = b x , \quad x \geq 0 . \]
Armageddon family
Suppose \(a, b \geq 0\) and consider the Armageddon Bernstein function
\[ \psi{(x)} = a 1_{\{ x > 0 \}} + b x , \quad x \geq 0 , \]
which is the sum of a constant Bernstein function and a linear Bernstein function.
-
The corresponding Lévy subordinator has an exponential killing rate \(a\), a linear drift \(b\), and no compound Poisson component:
\[ \Lambda_{t} = \infty \cdot 1_{\{ \epsilon \leq a t \}} + b t , \quad t \geq 0 . \]
-
The corresponding shock arrival intensities are
\[ \lambda_{I} = \begin{cases} b , & {\lvert I \rvert} = 1 , \\ a , & I = {\{ 1 , \ldots, d \}} , \\ 0 , & \text{else} , \end{cases} \]
such that a simplified exogenous shock model can be used.