Class for Pareto Bernstein functions
Source:R/s4-ParetoBernsteinFunction.R
ParetoBernsteinFunction-class.RdFor the Pareto-jump compound Poisson process with index \(0 < \alpha < 1\) and cutoff point \(x0\), the corresponding Bernstein function is $$ \psi(x) = 1 - e^{-x x_0} + (x_0 x)^\alpha \Gamma(1-\alpha, x_0 x) , 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_{x_0}^\infty e^{-ux} (1-e^{-u})^k
\alpha \frac{{x_0}^\alpha}{t^{1+\alpha}} du,
x>0, k>0 .
$$
The Pareto Bernstein function has the Lévy density \(\nu\): $$ \nu(du) = \alpha \frac{x_0^\alpha}{u^{\alpha + 1}}, \quad u > x_0 . $$
The Pareto Bernstein function, in combination with a linear Bernstein function can be used to approximate the Bernstein function of an \(\alpha\)-stable subordinator, see Sec. 5.3 of (Loroño 2015) .
References
Loroño LF (2015). Selected topics in financial engineering: first-exit times and dependence structures of Marshall-Olkin Kind. Ph.D. thesis, University of the Basque Country.
See also
getLevyDensity(), 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,
InverseGaussianBernsteinFunction-class,
LevyBernsteinFunction-class,
LinearBernsteinFunction-class,
PoissonBernsteinFunction-class,
ScaledBernsteinFunction-class,
SumOfBernsteinFunctions-class
Other Levy Bernstein function classes:
AlphaStableBernsteinFunction-class,
CompleteBernsteinFunction-class,
ExponentialBernsteinFunction-class,
GammaBernsteinFunction-class,
InverseGaussianBernsteinFunction-class,
LevyBernsteinFunction-class,
PoissonBernsteinFunction-class
Other Algebraic Bernstein function classes:
AlphaStableBernsteinFunction-class,
ExponentialBernsteinFunction-class,
InverseGaussianBernsteinFunction-class
Examples
# Create an object of class ParetoBernsteinFunction
ParetoBernsteinFunction()
#> An object of class "ParetoBernsteinFunction"
#> (invalid or not initialized)
ParetoBernsteinFunction(alpha = 0.2, x0 = 1e-2)
#> An object of class "ParetoBernsteinFunction"
#> - alpha: 0.2
#> - x0: 0.01
# Create a Lévy density
bf <- ParetoBernsteinFunction(alpha = 0.7, x0 = 1e-2)
levy_density <- getLevyDensity(bf)
integrate(
function(x) pmin(1, x) * levy_density(x),
lower = attr(levy_density, "lower"),
upper = attr(levy_density, "upper")
)
#> 0.1093688 with absolute error < 7.4e-05
# Evaluate the Bernstein function
bf <- ParetoBernsteinFunction(alpha = 0.3, x0 = 1)
calcIterativeDifference(bf, 1:5)
#> [1] 0.9421122 0.9865239 0.9963303 0.9989241 0.9996704
# Calculate shock-arrival intensities
bf <- ParetoBernsteinFunction(alpha = 0.8, x0 = 1e-2)
calcShockArrivalIntensities(bf, 3)
#> [1] 0.03709606 0.03709606 0.00846491 0.03709606 0.00846491 0.00846491 0.02132417
calcShockArrivalIntensities(bf, 3, tolerance = 1e-4)
#> [1] 0.037096063 0.037096063 0.008464909 0.037096063 0.008464909 0.008464909
#> [7] 0.021324170
# Calculate exchangeable shock-arrival intensities
bf <- ParetoBernsteinFunction(alpha = 0.4, x0 = 1e-2)
calcExShockArrivalIntensities(bf, 3)
#> [1] 0.04823193 0.02054932 0.14003630
calcExShockArrivalIntensities(bf, 3, tolerance = 1e-4)
#> [1] 0.04823193 0.02054932 0.14003630
# Calculate exchangeable shock-size arrival intensities
bf <- ParetoBernsteinFunction(alpha = 0.2, x0 = 1e-2)
calcExShockSizeArrivalIntensities(bf, 3)
#> [1] 0.12750303 0.07180652 0.37062172
calcExShockSizeArrivalIntensities(bf, 3, tolerance = 1e-4)
#> [1] 0.12750303 0.07180651 0.37062172
# Calculate the Markov generator
bf <- ParetoBernsteinFunction(alpha = 0.6, x0 = 1e-2)
calcMDCMGeneratorMatrix(bf, 3)
#> [,1] [,2] [,3] [,4]
#> [1,] -0.2257529 0.1306000 0.04112552 0.05402745
#> [2,] 0.0000000 -0.1822196 0.11448366 0.06773596
#> [3,] 0.0000000 0.0000000 -0.12497779 0.12497779
#> [4,] 0.0000000 0.0000000 0.00000000 0.00000000
calcMDCMGeneratorMatrix(bf, 3, tolerance = 1e-4)
#> [,1] [,2] [,3] [,4]
#> [1,] -0.2257529 0.1306000 0.04112552 0.05402745
#> [2,] 0.0000000 -0.1822196 0.11448366 0.06773596
#> [3,] 0.0000000 0.0000000 -0.12497779 0.12497779
#> [4,] 0.0000000 0.0000000 0.00000000 0.00000000