Bug in AlphaStableDistributions?

I could not reconcile the output of a classical sampling algorithm for the Gumbel Copula when using AlphaStableDistributions.jl. I wonder if there is a bug in its simulation algorithm or if I have overseen something?

A classical sampling algorithm for the (nested) Gumbel Copula, see e.g Hofert et al page 4, is based on sampling from a stable distribution. It produces in particular uniformly distributed random numbers u in the unit interval:

using AlphaStableDistributions, Statistics

# Select parameters for calculation
theta = 1.2 # needs to be > 1
nsamples = 10^6

# Algorithm 
# Step (1) 
scal(th) = (cos(π / (2.0 * th)))^th
beta = 1.0
s = AlphaStable(α=1.0 / theta, β=beta, scale=scal(theta), location=0.0)

v = rand(s, nsamples)

# Step (2)
r = -log.(rand(nsamples)) 

# Step (3)
psi(x, th) = exp(-x^(1.0 / th))
u = psi.(r ./ v, theta)

Assumption: The algorithm is specified for the so-called 1-paramaterization in Nolan which is probably the choice of AlphaStableDistributions as it implements the mean, where defined, as the location parameter.

In view of the sample size, mean and var of u are too far away from the theoretical values 0.5 and 0.08333...

julia> (mean(u), var(u))
(0.4578483790138697, 0.10688388000233129)

Replacing the samples v in Step (2) above using the algorithm [Nolan, Theorem 1.19, page 21] the result of the test is as expected:

# sketch for a fix
function stablernd(alpha, beta)  # is S(alpha, beta,1,0;1) distributed
    alpha != 1.0 || error("special case missing")
    t0 = atan(beta * tan((pi * alpha) / 2)) / alpha  # numerical problems if alpha near 1
    Theta = pi * (rand() - 0.5)
    W = -log(rand())
    term1 = sin(alpha * (t0 + Theta)) / (cos(alpha * t0) * cos(Theta))^(1 / alpha)
    term2 = ((cos(alpha * t0 + (alpha - 1) * Theta)) / W)^((1 - alpha) / alpha)
    return term1 * term2

# replace v and recalculate u from above
for k in 1:nsamples
    v[k] = stablernd(1 / theta, beta) * scal(theta)
u = psi.(r ./ v, theta)

Indeed, …

julia> (mean(u), var(u))
(0.5003790541200495, 0.08318313400529477)

Have you opened an issue on the repo? This seems like more of a bug report than a Discourse thread.

This indeed looks like a bug report for AlphaStableDistributions.jl. Note that you could reproduce your first piece of code, including in more dimensions, through :

using Copulas, Statistics
θ = 1.2
d = 3
G = GumbelCopula(d,θ)
nsamples = 10^6
u = rand(G,nsamples)

And indeed I get the same weird (0.4581850334671335, 0.10689564822591738), even for dimensions d > 2, while I totally agree with you that the mean should be a lot closer to 1/2.

This is using the same AlphaStable parametrisation as you are. Edit: The Copulas.jls one is implemented in this file so maybe it requires a separate fix from AlphaStableDistributions.jl.

Thank you for your feedback. I opened an issue on the repo.

1 Like

So both Copulas.jl and AlphaStableDistributions.jl were fixed and registered, so that now after updating your packages you have :

using Copulas, Statistics, Random
θ = 1.2
d = 2
G = GumbelCopula(d,θ)
nsamples = 10^4
u = rand(G,nsamples)

returning (0.5038099870266098, 0.08355682703493884) as it should be.

Thanks !

PS: May I ask what it is you were trying to do with this ?

Many thanks for this fix and the swift response.

P.S.: I detected the bug when I translated some Matlab code I wrote long ago to Julia in order to sample from nested Gumbel copulas - recursively through the tree with theta parameters in the vertices. I found them quite useful in building actuarial toy models.

The nested Archimedean trees copulas are on my todo-list for Copulas.jl. If you have something that works correctly, maybe you could propose it for a PR ? If you dont know how to do such a contribution to an open source repository, you may learn (its a great tool to have) or i could help you.