Skip to contents

For the \(\alpha\)-stable Lévy subordinator with \(0 < \alpha < 1\), the corresponding Bernstein function is the power function with exponent \(\alpha\), i.e. $$ \psi(x) = x^\alpha, \quad x>0. $$

Details

For the \(\alpha\)-stable 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 \alpha \frac{1}{\Gamma(1-\alpha) u^{1+\alpha}} du, x>0, k>0 . $$

This Bernstein function is no. 1 in the list of complete Bernstein functions in Chp. 16 of (Schilling et al. 2012) .

The \(\alpha\)-stable Bernstein function has the Lévy density \(\nu\): $$ \nu(du) = \frac{\alpha}{\Gamma(1-\alpha)} u^{-1 - \alpha} , \quad u > 0 , $$ and it has the Stieltjes density \(\sigma\): $$ \sigma(du) = \frac{\sin(\alpha \pi)}{\pi} u^{\alpha - 1}, \quad u > 0 . $$

Slots

alpha

The index \(\alpha\).

References

Schilling RL, Song R, Vondracek Z (2012). Bernstein functions, 2 edition. De Gruyter. doi:10.1515/9783110269338 .

Examples

# Create an object of class AlphaStableBernsteinFunction
AlphaStableBernsteinFunction()
#> An object of class "AlphaStableBernsteinFunction"
#> 	 (invalid or not initialized)
AlphaStableBernsteinFunction(alpha = 0.5)
#> An object of class "AlphaStableBernsteinFunction"
#> - alpha: 0.5

# Create a Lévy density
bf <- AlphaStableBernsteinFunction(alpha = 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.114243 with absolute error < 1.4e-05

# Create a Stieltjes density
bf <- AlphaStableBernsteinFunction(alpha = 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 < 8.4e-06

# Evaluate the Bernstein function
bf <- AlphaStableBernsteinFunction(alpha = 0.3)
calcIterativeDifference(bf, 1:5)
#> [1] 1.000000 1.231144 1.390389 1.515717 1.620657

# Calculate shock-arrival intensities
bf <- AlphaStableBernsteinFunction(alpha = 0.8)
calcShockArrivalIntensities(bf, 3)
#> [1] 0.66712356 0.66712356 0.07397757 0.66712356 0.07397757 0.07397757 0.18492131
calcShockArrivalIntensities(bf, 3, method = "stieltjes")
#> [1] 0.66712356 0.66712356 0.07397757 0.66712356 0.07397757 0.07397757 0.18492131
calcShockArrivalIntensities(bf, 3, tolerance = 1e-4)
#> [1] 0.66712356 0.66712356 0.07397757 0.66712356 0.07397757 0.07397757 0.18492130

# Calculate exchangeable shock-arrival intensities
bf <- AlphaStableBernsteinFunction(alpha = 0.4)
calcExShockArrivalIntensities(bf, 3)
#> [1] 0.23233766 0.08717025 0.59332184
calcExShockArrivalIntensities(bf, 3, method = "stieltjes")
#> [1] 0.23233766 0.08717025 0.59332184
calcExShockArrivalIntensities(bf, 3, tolerance = 1e-4)
#> [1] 0.23233766 0.08717027 0.59332185

# Calculate exchangeable shock-size arrival intensities
bf <- AlphaStableBernsteinFunction(alpha = 0.2)
calcExShockSizeArrivalIntensities(bf, 3)
#> [1] 0.2910978 0.1549973 0.7996359
calcExShockSizeArrivalIntensities(bf, 3, method = "stieltjes")
#> [1] 0.2910978 0.1549973 0.7996359
calcExShockSizeArrivalIntensities(bf, 3, tolerance = 1e-4)
#> [1] 0.2910978 0.1549974 0.7996358

# Calculate the Markov generator
bf <- AlphaStableBernsteinFunction(alpha = 0.6)
calcMDCMGeneratorMatrix(bf, 3)
#>           [,1]      [,2]       [,3]      [,4]
#> [1,] -1.933182  1.252396  0.2947533 0.3860323
#> [2,]  0.000000 -1.515717  1.0314331 0.4842834
#> [3,]  0.000000  0.000000 -1.0000000 1.0000000
#> [4,]  0.000000  0.000000  0.0000000 0.0000000
calcMDCMGeneratorMatrix(bf, 3, method = "stieltjes")
#>           [,1]      [,2]       [,3]      [,4]
#> [1,] -1.933182  1.252396  0.2947533 0.3860323
#> [2,]  0.000000 -1.515717  1.0314331 0.4842834
#> [3,]  0.000000  0.000000 -1.0000000 1.0000000
#> [4,]  0.000000  0.000000  0.0000000 0.0000000
calcMDCMGeneratorMatrix(bf, 3, tolerance = 1e-4)
#>           [,1]      [,2]       [,3]      [,4]
#> [1,] -1.933182  1.252396  0.2947533 0.3860323
#> [2,]  0.000000 -1.515717  1.0314331 0.4842834
#> [3,]  0.000000  0.000000 -1.0000000 1.0000000
#> [4,]  0.000000  0.000000  0.0000000 0.0000000