stable
R packageThis package is intended to be the developmental version to the CRAN version of Jim Lindsey’s stable. The .zip files listed on his homepage have been listed as version 1.0 since 2005. For the subsequent maintenance on this github and CRAN, we will start at 1.1.
To compare this version with the static v1.0 files on Jim Lindsey’s Homepage, it may be useful to use the compare page for this repo’s two branches.
stabledist
R packageIn brief, the parameters have different names and are transformations for each other. First, the names:
stabledist | stable |
---|---|
alpha | tail |
beta | skew |
gamma | disp |
delta | loc |
If you read the Lambert and Lindsey (1999 JRSS-C) PDF in this repo,
be aware that location is given the greek letter gamma and scale is
given the greek letter delta. The Nolan PDF does the opposite and is
used for stabledist
.
[Swihart 2022 update, see references in ?dstable
:] In
this README we detail how to make equivalent calls to those of
‘stabledist’ (i.e., Nolan’s 0-parameterization and 1-parameterization as
detailed in Nolan (2020)). See github for Lambert and Lindsey 1999
JRSS-C journal article, which details the parameterization of the Buck
(1995) stable distribution which allowed a Fourier inversion to arrive
at a form similar to but not exactly the
The Nolan (2020) reference is a textbook that provides an accessible and
comprehensive summary of stable distributions in the 25 years or so
since the core of this R package was made and put on CRAN.
The Buck (1995) parameterization most closely resembles the Zolotarev B
parameterization outlined in Definition 3.6 on page 93 of Nolan (2020) –
except that Buck (1995) did not allow the scale parameter to multiply
with the location parameter.
This explains why the Zolotarev B
entry in Table 3.1 on
page 97 of Nolan (2020) has the location parameter being multiplied by
the scale parameter whereas in converting the Lindsey and Lambert (1999)
to Nolan 1-parameterization the location parameter stays the same.
To be clear, and are evaluated by numerically integrating the inverse
Fourier transform. The code works reasonably for small and moderate
values of x, but will have numerical issues in some cases large x(such
as values from being greater than 1 or or not being monotonic). The
arguments , , , and can be adjusted to improve accuracy at the cost of
speed, but will still have limitations. Functions that avoid these
problems are available in other packages (such as and ) that use an
alternative method (as detailed in Nolan 1997)
distinct from directly numerically integrating the Fourier inverse
transform. See last example in this README.
For some values for some distributions things match up nicely, as we see with Normal and Cauchy:
<- 3
q ::pstable(q, tail =2, skew=0, disp =1, loc =0)
stable#> [1] 0.9830526
::pstable(q, alpha=2, beta=0, gamma=1, delta=0)
stabledist#> [1] 0.9830526
<- 3
q ::pstable(q, tail =1, skew=0, disp =1, loc =0)
stable#> [1] 0.8975836
::pstable(q, alpha=1, beta=0, gamma=1, delta=0)
stabledist#> [1] 0.8975836
However, to make stable
equivalent to
stabledist
in general, some transformations are needed.
Please see the following examples. Between stabledist
and
stable
, the alpha
is equivalent to
tail
and the delta
is equivalent to
loc
with no transformation. For the beta
(skew
) and gamma
(disp
)
parameters, a transformation is needed to get equivalent calls. Note
differences still may exist to numerical accuracy.
<- 0.9
q
# nolan pm=1 parameters:
<- 0.5
a <- 1
b <- .25
c <- 0.8
d
# lindsey-(3) page 415 conversion:
# tail/alpha and location stay the same
<- a
a3 <- d
d3 # the others require calcs:
<- cos(pi/2 * a)^2 + (-b)^2*sin(pi/2 * a)^2
DEL2 <- sqrt(DEL2) * sign(1-a)
DEL <- min(a, 2-a)
eta_a # the lindsey-(3) beta:
<- 2/(pi*eta_a)*acos( cos(pi/2 * a) / DEL )
b3 # the lindsey-(3) scale:
<- ( (DEL*c^a) / cos(pi/2 * a) )^(1/a)
c3
::pstable(q, tail =a, skew=b3, disp =c3, loc =d)
stable#> [1] 0.1154242
::pstable(q, alpha=a, beta=b , gamma=c , delta=d, pm=1)
stabledist#> [1] 0.1138462
::plevy(q, m=d, s=c)
rmutil#> [1] 0.1138463
# more accuracy!!!!?!
::pstable(q, tail =a, skew=b3, disp =c3, loc =d, eps = 0.13*1e-7)
stable#> [1] 0.1138786
<- 0.9
q
# nolan pm=1 parameters:
<- 0.5
a <- 1
b <- .25
c <- 0.8
d
# lindsey-(3) page 415 conversion:
# tail/alpha and location stay the same
<- a
a3 <- d
d3 # the others require calcs:
<- cos(pi/2 * a)^2 + (-b)^2*sin(pi/2 * a)^2
DEL2 <- sqrt(DEL2) * sign(1-a)
DEL <- min(a, 2-a)
eta_a # the lindsey-(3) beta:
<- 2/(pi*eta_a)*acos( cos(pi/2 * a) / DEL )
b3 # the lindsey-(3) scale:
<- ( (DEL*c^a) / cos(pi/2 * a) )^(1/a)
c3
::dstable(q, tail =a, skew=b3, disp =c3, loc =d)
stable#> [1] 1.806389
::dstable(q, alpha=a, beta=b , gamma=c , delta=d, pm=1)
stabledist#> Warning in uniroot(function(th) log(g(th)), lower = l.th, upper = u.th, : -
#> Inf replaced by maximally negative value
#> Warning in uniroot(function(th) log(g(th)), lower = l.th, upper = u.th, : -
#> Inf replaced by maximally negative value
#> Warning in .integrate2(g1, lower = a, upper = b, subdivisions =
#> subdivisions, : roundoff error is detected in the extrapolation table
#> [1] 1.807224
::dlevy(q, m=d, s=c)
rmutil#> [1] 1.807224
<- .3
p
# nolan pm=1 parameters:
<- 0.5
a <- 1
b <- .25
c <- 0.8
d
# lindsey-(3) page 415 conversion:
# tail/alpha and location stay the same
<- a
a3 <- d
d3 # the others require calcs:
<- cos(pi/2 * a)^2 + (-b)^2*sin(pi/2 * a)^2
DEL2 <- sqrt(DEL2) * sign(1-a)
DEL <- min(a, 2-a)
eta_a # the lindsey-(3) beta:
<- 2/(pi*eta_a)*acos( cos(pi/2 * a) / DEL )
b3 # the lindsey-(3) scale:
<- ( (DEL*c^a) / cos(pi/2 * a) )^(1/a)
c3
::qstable(p, tail =a, skew=b3, disp =c3, loc =d)
stable#> [1] 1.031301
::qstable(p, alpha=a, beta=b , gamma=c , delta=d, pm=1)
stabledist#> [1] 1.032735
::qlevy(p, m=d, s=c)
rmutil#> [1] 1.032733
<- -1.97
q
# nolan pm=1 parameters:
<- 0.8
a <- 0
b <- 1
c <- 0
d
# lindsey-(3) page 415 conversion:
# tail/alpha and location stay the same
<- a
a3 <- d
d3 # the others require calcs:
<- cos(pi/2 * a)^2 + (-b)^2*sin(pi/2 * a)^2
DEL2 <- sqrt(DEL2) * sign(1-a)
DEL <- min(a, 2-a)
eta_a # the lindsey-(3) beta:
<- 2/(pi*eta_a)*acos( cos(pi/2 * a) / DEL )
b3 # the lindsey-(3) scale:
<- ( (DEL*c^a) / cos(pi/2 * a) )^(1/a)
c3
::pstable(q, tail =a, skew=b3, disp =c3, loc =d)
stable#> [1] 0.1722953
::pstable(q, alpha=a, beta=b , gamma=c , delta=d)
stabledist#> [1] 0.1722945
<- -1
q
# nolan pm=1 parameters:
<- 1.3
a <- 0.4
b <- 2
c <- 0.75
d
# lindsey-(3) page 415 conversion:
# tail/alpha and location stay the same
<- a
a3 <- d
d3 # the others require calcs:
<- cos(pi/2 * a)^2 + (-b)^2*sin(pi/2 * a)^2
DEL2 <- sqrt(DEL2) * sign(1-a)
DEL <- min(a, 2-a)
eta_a # the lindsey-(3) beta:
<- -sign(b)*2/(pi*eta_a)*acos( cos(pi/2 * a) / DEL )
b3 # the lindsey-(3) scale:
<- ( (DEL*c^a) / cos(pi/2 * a) )^(1/a)
c3
::pstable(q, tail =a, skew=b3, disp =c3, loc =d)
stable#> [1] 0.4349168
::pstable(q, alpha=a, beta=b , gamma=c , delta=d, pm=1)
stabledist#> [1] 0.4348957
::dstable(q, tail =a, skew=b3, disp =c3, loc =d)
stable#> [1] 0.1454112
::dstable(q, alpha=a, beta=b , gamma=c , delta=d, pm=1)
stabledist#> [1] 0.1454111
sd2s
and s2sd
<- -1
q # nolan pm=1 parameters:
<- 1.3
a <- -0.4
b <- 2
c <- 0.75
d # sd2s takes nolan (stabledist) parameters and returns lindsey (stable)
<- stable::sd2s(alpha=a, beta=b, gamma=c, delta=d)
s ::pstable(q, tail = s$tail, skew=s$skew, disp = s$disp, loc = s$loc)
stable#> [1] 0.196531
::pstable(q, alpha=a, beta=b , gamma=c , delta=d, pm=1)
stabledist#> [1] 0.1965513
# s2sd takes lindsey (stable) parameters and returns nolan (stabledist)
<- stable::s2sd(tail = s$tail, skew=s$skew, disp = s$disp, loc = s$loc)
sd ::pstable(q, alpha=sd$alpha, beta=sd$beta , gamma=sd$gamma , delta=sd$delta, pm=1)
stabledist#> [1] 0.1965513
<- -1
q # nolan pm=1 parameters:
<- 1.3
a1 <- -0.4
b1 <- 2
c1 <- 0.75
d1 # for a1 != 1
<- d1 + b1*c1*tan(pi*a1/2)
d0
# Calculate d0 by hand or use pm1_to_pm0():
# Convert to nolan pm=0 parameters:
<- stable::pm1_to_pm0(a1,b1,c1,d1)
pm0 <- pm0$a0
a0 <- pm0$b0
b0 <- pm0$c0
c0 <- pm0$d0
d0 # check:
::pstable(q, alpha=a1, beta=b1 , gamma=c1 , delta=d1, pm=1)
stabledist#> [1] 0.1965513
# only change delta=d0 for pm=0
::pstable(q, alpha=a1, beta=b1 , gamma=c1 , delta=d0, pm=0)
stabledist#> [1] 0.1965513
::pstable(q, alpha=a0, beta=b0 , gamma=c0 , delta=d0, pm=0)
stabledist#> [1] 0.1965513
::dstable(q, alpha=a1, beta=b1 , gamma=c1 , delta=d1, pm=1)
stabledist#> [1] 0.0572133
# only change delta=d0 for pm=0
::dstable(q, alpha=a1, beta=b1 , gamma=c1 , delta=d0, pm=0)
stabledist#> [1] 0.0572133
::dstable(q, alpha=a0, beta=b0 , gamma=c0 , delta=d0, pm=0)
stabledist#> [1] 0.0572133
<- -1
q # nolan pm=1 parameters:
# a1 <- 1.3
# b1 <- 0.4
# c1 <- 2
# d1 <- 0.75
<- 1.3
a1 <- .5
b1 <- 1
c1 <- 0
d1 # for a1 != 1
<- d1 + b1*c1*tan(pi*a1/2)
d0
<- stable::sd2s(alpha=a1, beta=b1, gamma=c1, delta=d1)
s ::stable.mode(tail = s$tail, skew=s$skew, disp = s$disp, loc = s$loc)$ytilde
stable#> [1] -1.13224
*stabledist::stableMode(alpha=a1, beta=b1)+d0
c1#> [1] -1.133257
<- seq(-2.5,2.6,0.001)
xran <- stabledist::dstable(xran, alpha=a1, beta=b1, gamma=c1, delta=d1, pm=1)
ysd #plot(xran, ysd)
== max(ysd)]
xran[ysd #> [1] -1.133
<- stable::dstable(xran, tail = s$tail, skew=s$skew, disp = s$disp, loc = s$loc)
ys #points(xran, ys, col="blue")
== max(ys)]
xran[ys #> [1] -1.133
<- stable::s2sd(1.5, 0.5, 1/sqrt(2), 0)
param_conv
param_convhead(stable::dstable(q, tail =1.5, skew=0.5, disp =1/sqrt(2), loc = 0))
head(stabledist::dstable(q, alpha=param_conv$alpha, beta=param_conv$beta , gamma=param_conv$gamma , delta=param_conv$delta, pm=1))
plot(q,stable::dstable(q, tail =1.5, skew=0.5, disp =1/sqrt(2), loc = 0), type="s")
plot(q,stabledist::dstable(q, alpha=param_conv$alpha, beta=param_conv$beta , gamma=param_conv$gamma , delta=param_conv$delta, pm=1), type="s")
<- stable::s2sd(1.5, 0.5, 1/sqrt(2), 0)
param_conv
param_convplot(q,stable::pstable(q, tail =1.5, skew=0.5, disp =1/sqrt(2), loc = 0), type="s")
plot(q,stabledist::pstable(q, alpha=param_conv$alpha, beta=param_conv$beta , gamma=param_conv$gamma , delta=param_conv$delta, pm=1), type="s")