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
end
# replace v and recalculate u from above
for k in 1:nsamples
v[k] = stablernd(1 / theta, beta) * scal(theta)
end
u = psi.(r ./ v, theta)
Indeed, …
julia> (mean(u), var(u))
(0.5003790541200495, 0.08318313400529477)