Class for Inverse Gaussian Bernstein function
Source:R/s4-InverseGaussianBernsteinFunction.R
InverseGaussianBernsteinFunction-class.RdFor the inverse Gaussian Lévy subordinator with \(\eta > 0\), the corresponding Bernstein function is the function $$ \psi(x) = \sqrt{2x + \eta^2} - \eta, x>0. $$
Details
For the inverse Gaussian Bernstein function, the higher-order alternating
iterated forward differences are not known in closed-form, 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{1}{\sqrt{2\pi}
u^{3/2}} e^{-\frac{1}{2}\eta^2 u} du, x>0, k>0.
$$
This Bernstein function can be found on p. 309 in (Mai and Scherer 2017) . Furthermore it is a transformation of no. 2 in the list of complete Bernstein functions in Chp. 16 of (Schilling et al. 2012) .
The inverse Gaussian Bernstein function has the Lévy density \(\nu\): $$ \nu(du) = \frac{1}{\sqrt{2 \pi u^3}} \operatorname{e}^{-\frac{1}{2} \eta^2 u} , \quad u > 0 , $$
and it has the Stieltjes density \(\sigma\): $$ \sigma(du) = \frac{ \sin(\pi / 2) }{ \pi } \cdot \frac{ \sqrt{2 x - \eta^2} }{ x } , \quad u > \eta^2 / 2 . $$
References
Mai J, Scherer M (2017).
Simulating copulas: stochastic models, sampling algorithms and applications, Series in Quantitative Finance, 2 edition.
World Scientific.
doi:10.1142/10265
.
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,
GammaBernsteinFunction-class,
LevyBernsteinFunction-class,
LinearBernsteinFunction-class,
ParetoBernsteinFunction-class,
PoissonBernsteinFunction-class,
ScaledBernsteinFunction-class,
SumOfBernsteinFunctions-class
Other Levy Bernstein function classes:
AlphaStableBernsteinFunction-class,
CompleteBernsteinFunction-class,
ExponentialBernsteinFunction-class,
GammaBernsteinFunction-class,
LevyBernsteinFunction-class,
ParetoBernsteinFunction-class,
PoissonBernsteinFunction-class
Other Complete Bernstein function classes:
AlphaStableBernsteinFunction-class,
ExponentialBernsteinFunction-class,
GammaBernsteinFunction-class
Other Algebraic Bernstein function classes:
AlphaStableBernsteinFunction-class,
ExponentialBernsteinFunction-class,
ParetoBernsteinFunction-class
Examples
# Create an object of class InverseGaussianBernsteinFunction
InverseGaussianBernsteinFunction()
#> An object of class "InverseGaussianBernsteinFunction"
#> (invalid or not initialized)
InverseGaussianBernsteinFunction(eta = 0.3)
#> An object of class "InverseGaussianBernsteinFunction"
#> - eta: 0.3
# Create a Lévy density
bf <- InverseGaussianBernsteinFunction(eta = 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.023005 with absolute error < 9.2e-06
# Create a Stieltjes density
bf <- InverseGaussianBernsteinFunction(eta = 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 with absolute error < 9e-05
# Evaluate the Bernstein function
bf <- InverseGaussianBernsteinFunction(eta = 0.3)
calcIterativeDifference(bf, 1:5)
#> [1] 1.145683 1.722375 2.167793 2.544293 2.876476
# Calculate shock-arrival intensities
bf <- InverseGaussianBernsteinFunction(eta = 0.8)
calcShockArrivalIntensities(bf, 3)
#> [1] 0.4227538 0.4227538 0.1065044 0.4227538 0.1065044 0.1065044 0.1890450
calcShockArrivalIntensities(bf, 3, method = "stieltjes")
#> [1] 0.4227538 0.4227538 0.1065044 0.4227538 0.1065044 0.1065044 0.1890450
calcShockArrivalIntensities(bf, 3, tolerance = 1e-4)
#> [1] 0.4227538 0.4227538 0.1065044 0.4227538 0.1065044 0.1065044 0.1890450
# Calculate exchangeable shock-arrival intensities
bf <- InverseGaussianBernsteinFunction(eta = 0.4)
calcExShockArrivalIntensities(bf, 3)
#> [1] 0.4423269 0.1275870 0.3721928
calcExShockArrivalIntensities(bf, 3, method = "stieltjes")
#> [1] 0.4423269 0.1275870 0.3721928
calcExShockArrivalIntensities(bf, 3, tolerance = 1e-4)
#> [1] 0.4423269 0.1275871 0.3721929
# Calculate exchangeable shock-size arrival intensities
bf <- InverseGaussianBernsteinFunction(eta = 0.2)
calcExShockSizeArrivalIntensities(bf, 3)
#> [1] 1.3429981 0.4020703 0.5125728
calcExShockSizeArrivalIntensities(bf, 3, method = "stieltjes")
#> [1] 1.3429981 0.4020703 0.5125728
calcExShockSizeArrivalIntensities(bf, 3, tolerance = 1e-4)
#> [1] 1.3429981 0.4020703 0.5125728
# Calculate the Markov generator
bf <- InverseGaussianBernsteinFunction(eta = 0.6)
calcMDCMGeneratorMatrix(bf, 3)
#> [,1] [,2] [,3] [,4]
#> [1,] -1.921904 1.301528 0.3539682 0.2664076
#> [2,] 0.000000 -1.488061 1.1036643 0.3843970
#> [3,] 0.000000 0.000000 -0.9362291 0.9362291
#> [4,] 0.000000 0.000000 0.0000000 0.0000000
calcMDCMGeneratorMatrix(bf, 3, method = "stieltjes")
#> [,1] [,2] [,3] [,4]
#> [1,] -1.921904 1.301528 0.3539682 0.2664076
#> [2,] 0.000000 -1.488061 1.1036643 0.3843970
#> [3,] 0.000000 0.000000 -0.9362291 0.9362291
#> [4,] 0.000000 0.000000 0.0000000 0.0000000
calcMDCMGeneratorMatrix(bf, 3, tolerance = 1e-4)
#> [,1] [,2] [,3] [,4]
#> [1,] -1.921904 1.301528 0.3539682 0.2664076
#> [2,] 0.000000 -1.488061 1.1036643 0.3843970
#> [3,] 0.000000 0.000000 -0.9362291 0.9362291
#> [4,] 0.000000 0.000000 0.0000000 0.0000000