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,
u are too far away from the theoretical values
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)
julia> (mean(u), var(u)) (0.5003790541200495, 0.08318313400529477)