rpextmo: Simulate from parametrized families of extendible MO...
In hsloot/rmo: A package for the Marshall-Olkin distribution

rpextmo

R Documentation

Simulate from parametrized families of extendible MO distributions

Description

Draws n iid d-variate samples from a parametrized family of extendible MO distributions.

Usage

rpextmo(
  n,
  d,
  a = 0,
  b = 0,
  gamma = 1,
  eta = NULL,
  family = c("Armageddon", "Poisson", "Pareto", "Exponential", "AlphaStable",
    "InverseGaussian", "Gamma"),
  method = c("MDCM", "LFM", "AM", "ESM")
)

Arguments

`n`	An integer for the number of samples.
`d`	An integer for the dimension.
`a`	A non-negative double representing the killing-rate `a` of the Bernstein function.
`b`	A non-negative double representing the drift `b` of the Bernstein function.
`gamma`	a positive double representing the scaling factor of the integral part of the Bernstein function.
`eta`	A numeric vector representing the distribution family's parameters, see the Details section.
`family`	A string representing the parametrized family. Use `"Armageddon"` for the Armageddon family, `"Poisson"` for the Poisson family, `"Pareto"` for the Pareto family, `"Exponential"` for the Exponential family, `"AlphaStable"` for the `\alpha`-stable family, `"InverseGaussian"` for the Inverse-Gaussian family, `"Gamma"` for the Gamma family, see the Details section.
`method`	A string representing which sampling algorithm should be used. Use `"MDCM"` for the Markovian death-set model, `"LFM"` for the Lévy–frailty model, `"AM"` for the Arnold model, and `"ESM"` for the exogenous shock model (in case of the Armageddon family, the algorithm is optimized to consider only finite shocks). We recommend using the ESM only for small dimensions; the AM can be used up to dimension `30`.

Details

A parametrized ext. MO distribution is a family of ext. MO distributions, see rextmo(), corresponding to Bernstein functions of the form

\psi{(x)} = 1_{\{ x > 0 \}} a + b x + \gamma \cdot \int_{0}^{\infty}{ {[1 - e^{-ux}]} \nu{(\mathrm{d}u)} }, \quad x \geq 0 ,

\psi{(x)} = 1_{\{ x > 0 \}} a + b x + \gamma \cdot \int_{0}^{\infty}{ \frac{x}{x + u} \sigma{(\mathrm{d}u)} }, \quad x \geq 0 ,

where a, b \geq 0. The \nu and \sigma represent the Lévy measure and Stieltjes measure, respectively. At least one of the following conditions must hold: a > 0, b > 0, or \nu \not\equiv 0 (resp. \sigma \not\equiv 0).

Families

All implemented families are listed in the following; some re-combinations are possible, see ScaledBernsteinFunction, SumOfBernsteinFunctions, and CompositeScaledBernsteinFunction.
Armageddon family: We have for \nu = \sigma \equiv 0 the Bernstein function \psi:

\psi{(x)} = 1_{\{ x > 0\}} a + b x , \quad x \geq 0 ,

see ConstantBernsteinFunction and LinearBernsteinFunction.
Poisson family: We have for \eta > 0 the Bernstein function \psi::

\psi{(x)} = 1_{\{ x > 0\}} a + b x + \gamma \cdot {[1 - e^{-\eta x}]}, \quad x \geq 0 ,

with (discrete) Lévy measure \nu:

\nu{(\mathrm{d}u)} = \delta_{\{ \eta \}}{(\mathrm{d}u)} ,

see PoissonBernsteinFunction.
Pareto family: We have \eta \in \mathbb{R}^2 with \eta_1 \in {(0, 1)}, \eta_2 > 0 a Bernstein function with Lévy measure \nu:

\nu{(\mathrm{d}u)} = \eta_{1} \eta_{2}^{\eta_{1}} \cdot u^{-\eta_{1}-1} 1_{\{ u > \eta_{2}\}} \mathrm{d}u ,

see ParetoBernsteinFunction.
Exponential family: We have for \eta > 0 the Bernstein function \psi:

\psi{(x)} = 1_{\{ x > 0\}} a + b x + \gamma \cdot \frac{x}{x + \eta} , \quad x \geq 0 ,

with Lévy measure \nu:

\nu{(\mathrm{d}u)} = \eta e^{-\eta u} \mathrm{d}u ,

and with (discrete) Stieltjes measure \sigma:

\sigma{(\mathrm{d}u)} = \delta_{\{ \eta \}}{(\mathrm{d}u)} ,

see ExponentialBernsteinFunction.
\alpha-stable family: We have for \eta \in {(0, 1)} the Bernstein function \psi:

\psi{(x)} = 1_{\{ x > 0\}} a + b x + \gamma \cdot x^{\eta} , \quad x \geq 0 ,

with Lévy measure \nu:

\nu{(\mathrm{d}u)} = \frac{\eta}{\Gamma{(1 - \eta)}} \cdot u^{-\eta-1} \mathrm{d}u ,

and with Stieltjes measure \sigma:

\sigma{(\mathrm{d}u)} = \frac{\sin{(\eta \pi)}}{\pi} \cdot u^{\eta - 1} \mathrm{d}u ,

see AlphaStableBernsteinFunction.
Inverse-Gaussian family: We have for \eta > 0 the Bernstein function \psi:

\psi{(x)} = 1_{\{ x > 0\}} a + b x + \gamma \cdot {\left[ \sqrt{2 x + \eta^2} - \eta \right]}, \quad x \geq 0 ,

with Lévy measure \nu:

\nu{(\mathrm{d}u)} = \frac{1}{ \sqrt{2 \pi} } \cdot \frac{ e^{-\frac{1}{2} \eta^2 u} }{ \sqrt{u^3} } \mathrm{d}u ,

and with Stieltjes measure \sigma:

\sigma{(\mathrm{d}u)} = \frac{\sin{(\pi / 2)}}{\pi} \cdot \frac{\sqrt{2 u - \eta^2}} {u} 1_{\{ u > \eta^2 / 2 \}} \mathrm{d}u ,

see InverseGaussianBernsteinFunction.
Gamma family: We have for \eta > 0 the Bernstein function \psi:

\psi{(x)} = 1_{\{ x > 0\}} a + b x + \gamma \cdot \log{\left( 1 + \frac{x}{\eta} \right)} , \quad x \geq 0 ,

with Lévy measure \nu:

\nu{(\mathrm{d}u)} = e^{-\eta u} u^{-1} \mathrm{d}u ,

and with Stieltjes measure \sigma:

\sigma{(\mathrm{d}u)} = u^{-1} 1_{\{ u > \eta \}} \mathrm{d}u ,

see GammaBernsteinFunction.

Simulation algorithms

The MDCM, AM, and ESM simulation algorithms for the exchangeable Marshall–Olkin distribution can be used. For this, the corresponding Bernstein function is passed to rextmo(). An exception is the ESM for the Armageddon family which uses an optimized version considering only finite shock-times.
The Lévy-frailty model (LFM) simulates the elements of the random vector as first-hitting times of a compound Poisson subordinator \Lambda into sets (E_i, \infty) for iid unit exponential random variables. Here, the subordinator is a linear combination of a pure-drift subordinator, a pure-killing subordinator, and a pure-jump compound Poisson subordinator, i.e.

\Lambda_{t} = \infty \cdot 1_{\{ \epsilon > a t \}} + b t + \sum_{j=1}^{N_{\gamma t}} X_{j} , \quad t \geq 0,

where \epsilon is a unit exponential rv, n is a Poisson process, and X_{1}, X_{2}, \ldots are iid jumps from the corresponding jump distribution, see \insertCite@see pp. 140 sqq. @Mai2017armo.

Value

rpextmo returns a numeric matrix of size n x d. Each row corresponds to an independently and identically (iid) distributed sample from a d-variate parametrized extendible Marshall–Olkin distribution with the specified parameters.

References

\insertAllCited

Examples

## Armageddon
rpextmo(10, 3, a = 0.2, b = 0.5)
## comonotone
rpextmo(10, 3, a = 1)
## independence
rpextmo(10, 3, b = 1)

rpextmo(10, 3, a = 0.2, b = 0.5, method = "ESM")
## comonotone
rpextmo(10, 3, a = 1, method = "ESM")
## independence
rpextmo(10, 3, b = 1, method = "ESM")

rpextmo(10, 3, a = 0.2, b = 0.5, method = "LFM")
## comonotone
rpextmo(10, 3, a = 1, method = "LFM")
## independence
rpextmo(10, 3, b = 1, method = "LFM")

rpextmo(10, 3, a = 0.2, b = 0.5, method = "MDCM")
## comonotone
rpextmo(10, 3, a = 1, method = "MDCM")
## independence
rpextmo(10, 3, b = 1, method = "MDCM")

rpextmo(10, 3, a = 0.2, b = 0.5, method = "AM")
## comonotone
rpextmo(10, 3, a = 1, method = "AM")
## independence
rpextmo(10, 3, b = 1, method = "AM")

rpextmo(10, 3, a = 0.2, b = 0.5, family = "Armageddon")
## comonotone
rpextmo(10, 3, a = 1, family = "Armageddon")
## independence
rpextmo(10, 3, b = 1, family = "Armageddon")

rpextmo(
  10, 3,
  a = 0.2, b = 0.5,
  family = "Armageddon",
  method = "ESM"
)
## comonotone
rpextmo(
  10, 3,
  a = 1,
  family = "Armageddon",
  method = "ESM"
)
## independence
rpextmo(
  10, 3,
  b = 1,
  family = "Armageddon",
  method = "ESM"
)

rpextmo(
  10, 3,
  a = 0.2, b = 0.5,
  family = "Armageddon",
  method = "LFM"
)
## comonotone
rpextmo(
  10, 3,
  a = 1,
  family = "Armageddon",
  method = "LFM"
)
## independence
rpextmo(
  10, 3,
  b = 1,
  family = "Armageddon",
  method = "LFM"
)

rpextmo(
  10, 3,
  a = 0.2, b = 0.5,
  family = "Armageddon",
  method = "MDCM"
)
## comonotone
rpextmo(
  10, 3,
  a = 1,
  family = "Armageddon",
  method = "MDCM"
)
## independence
rpextmo(
  10, 3,
  b = 1,
  family = "Armageddon",
  method = "MDCM"
)

rpextmo(
  10, 3,
  a = 0.2, b = 0.5,
  family = "Armageddon",
  method = "AM"
)
## comonotone
rpextmo(
  10, 3,
  a = 1,
  family = "Armageddon",
  method = "AM"
)
## independence
rpextmo(
  10, 3,
  b = 1,
  family = "Armageddon",
  method = "AM"
)

## Poisson
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "Poisson"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "Poisson",
  method = "ESM"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "Poisson",
  method = "LFM"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "Poisson",
  method = "MDCM"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "Poisson",
  method = "AM"
)

## Pareto
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = c(0.5, 1e-4), family = "Pareto"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = c(0.5, 1e-4), family = "Pareto",
  method = "ESM"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = c(0.5, 1e-4), family = "Pareto",
  method = "LFM"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = c(0.5, 1e-4), family = "Pareto",
  method = "MDCM"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = c(0.5, 1e-4), family = "Pareto",
  method = "AM"
)

## Exponential
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "Exponential"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "Exponential",
  method = "ESM"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "Exponential",
  method = "LFM"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "Exponential",
  method = "MDCM"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "Exponential",
  method = "AM"
)

## Alpha-Stable
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "AlphaStable"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "AlphaStable",
  method = "ESM"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "AlphaStable",
  method = "MDCM"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "AlphaStable",
  method = "AM"
)

## Inverse Gaussian
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "InverseGaussian"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "InverseGaussian",
  method = "ESM"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "InverseGaussian",
  method = "MDCM"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "InverseGaussian",
  method = "AM"
)

## Gamma
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "Gamma"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "Gamma",
  method = "ESM"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "Gamma",
  method = "MDCM"
)
rpextmo(
  10, 3,
  a = 0.2, b = 0.5, gamma = 2,
  eta = 0.5, family = "Gamma",
  method = "AM"
)

hsloot/rmo documentation built on May 1, 2024, 4:28 a.m.