Class for Gamma Bernstein functions
Source:R/s4-GammaBernsteinFunction.R
GammaBernsteinFunction-class.RdThe Gamma Bernstein function, is the Bernstein function of a subordinator with a (scaled) Gamma distribution. The representation is for \(a > 0\) $$ \psi(x) = \log(1 + \frac{x}{a}), x > 0. $$
Details
For this Bernstein function, the higher-order alternating iterated forward
differences are known in closed form but cannot be evaluated numerically
without the danger of loss of significance. But we can use numerical
integration (here: stats::integrate()) to approximate it with the
following representation:
$$
{(-1)}^{k-1} \Delta^{k} \psi(x)
= \int_{0}^{\infty} e^{-ux} {(1 - e^{-u})}^{k}
\frac{e^{-au}}{u} du, x>0, k>0.
$$
This Bernstein function is no. 26 in the list of complete Bernstein functions in Chp. 16 of (Schilling et al. 2012) .
The Gamma Bernstein function has the Lévy density \(\nu\): $$ \nu(du) = \frac{\operatorname{e}^{-a u}}{u}, \quad u > 0 , $$
and it has the Stieltjes density \(\sigma\): $$ \sigma(du) = u^{-1} du, u > a . $$
References
Schilling RL, Song R, Vondracek Z (2012). Bernstein functions, 2 edition. De Gruyter. doi:10.1515/9783110269338 .
See also
getLevyDensity(), getStieltjesDensity(),
calcIterativeDifference(), calcShockArrivalIntensities(),
calcExShockArrivalIntensities(), calcExShockSizeArrivalIntensities(),
calcMDCMGeneratorMatrix(), rextmo(), rpextmo()
Other Bernstein function classes:
AlphaStableBernsteinFunction-class,
BernsteinFunction-class,
CompleteBernsteinFunction-class,
CompositeScaledBernsteinFunction-class,
ConstantBernsteinFunction-class,
ConvexCombinationOfBernsteinFunctions-class,
ExponentialBernsteinFunction-class,
InverseGaussianBernsteinFunction-class,
LevyBernsteinFunction-class,
LinearBernsteinFunction-class,
ParetoBernsteinFunction-class,
PoissonBernsteinFunction-class,
ScaledBernsteinFunction-class,
SumOfBernsteinFunctions-class
Other Levy Bernstein function classes:
AlphaStableBernsteinFunction-class,
CompleteBernsteinFunction-class,
ExponentialBernsteinFunction-class,
InverseGaussianBernsteinFunction-class,
LevyBernsteinFunction-class,
ParetoBernsteinFunction-class,
PoissonBernsteinFunction-class
Other Complete Bernstein function classes:
AlphaStableBernsteinFunction-class,
ExponentialBernsteinFunction-class,
InverseGaussianBernsteinFunction-class
Examples
# Create an object of class GammaBernsteinFunction
GammaBernsteinFunction()
#> An object of class "GammaBernsteinFunction"
#> (invalid or not initialized)
GammaBernsteinFunction(a = 2)
#> An object of class "GammaBernsteinFunction"
#> - a: 2
# Create a Lévy density
bf <- GammaBernsteinFunction(a = 0.7)
levy_density <- getLevyDensity(bf)
integrate(
function(x) pmin(1, x) * levy_density(x),
lower = attr(levy_density, "lower"),
upper = attr(levy_density, "upper")
)
#> 1.092933 with absolute error < 3.1e-05
# Create a Stieltjes density
bf <- GammaBernsteinFunction(a = 0.5)
stieltjes_density <- getStieltjesDensity(bf)
integrate(
function(x) 1/(1 + x) * stieltjes_density(x),
lower = attr(stieltjes_density, "lower"),
upper = attr(stieltjes_density, "upper")
)
#> 1.098612 with absolute error < 8.9e-13
# Evaluate the Bernstein function
bf <- GammaBernsteinFunction(a = 0.3)
calcIterativeDifference(bf, 1:5)
#> [1] 1.466337 2.036882 2.397895 2.662588 2.871680
# Calculate shock-arrival intensities
bf <- GammaBernsteinFunction(a = 0.8)
calcShockArrivalIntensities(bf, 3)
#> [1] 0.3053816 0.3053816 0.1364511 0.3053816 0.1364511 0.1364511 0.2326464
calcShockArrivalIntensities(bf, 3, method = "stieltjes")
#> [1] 0.3053816 0.3053816 0.1364511 0.3053816 0.1364511 0.1364511 0.2326464
calcShockArrivalIntensities(bf, 3, tolerance = 1e-4)
#> [1] 0.3053816 0.3053816 0.1364511 0.3053816 0.1364511 0.1364511 0.2326464
# Calculate exchangeable shock-arrival intensities
bf <- GammaBernsteinFunction(a = 0.4)
calcExShockArrivalIntensities(bf, 3)
#> [1] 0.3483067 0.1906898 0.5230767
calcExShockArrivalIntensities(bf, 3, method = "stieltjes")
#> [1] 0.3483067 0.1906898 0.5230767
calcExShockArrivalIntensities(bf, 3, tolerance = 1e-4)
#> [1] 0.3483067 0.1906898 0.5230767
# Calculate exchangeable shock-size arrival intensities
bf <- GammaBernsteinFunction(a = 0.2)
calcExShockSizeArrivalIntensities(bf, 3)
#> [1] 1.1240803 0.6943271 0.9541813
calcExShockSizeArrivalIntensities(bf, 3, method = "stieltjes")
#> [1] 1.1240803 0.6943271 0.9541813
calcExShockSizeArrivalIntensities(bf, 3, tolerance = 1e-4)
#> [1] 1.1240803 0.6943271 0.9541813
# Calculate the Markov generator
bf <- GammaBernsteinFunction(a = 0.6)
calcMDCMGeneratorMatrix(bf, 3)
#> [,1] [,2] [,3] [,4]
#> [1,] -1.791759 0.9762672 0.4802562 0.3352360
#> [2,] 0.000000 -1.4663371 0.9710156 0.4953214
#> [3,] 0.000000 0.0000000 -0.9808293 0.9808293
#> [4,] 0.000000 0.0000000 0.0000000 0.0000000
calcMDCMGeneratorMatrix(bf, 3, method = "stieltjes")
#> [,1] [,2] [,3] [,4]
#> [1,] -1.791759 0.9762672 0.4802562 0.3352360
#> [2,] 0.000000 -1.4663371 0.9710156 0.4953214
#> [3,] 0.000000 0.0000000 -0.9808293 0.9808293
#> [4,] 0.000000 0.0000000 0.0000000 0.0000000
calcMDCMGeneratorMatrix(bf, 3, tolerance = 1e-4)
#> [,1] [,2] [,3] [,4]
#> [1,] -1.791759 0.9762672 0.4802562 0.3352360
#> [2,] 0.000000 -1.4663371 0.9710156 0.4953214
#> [3,] 0.000000 0.0000000 -0.9808293 0.9808293
#> [4,] 0.000000 0.0000000 0.0000000 0.0000000