About the Gumbel-Barnett Archimedean copula

Part 1: About the Gumbel-Barnett generator and its derivatives.

In this article, I discuss an interesting problem encountered while working on Copulas.jl, a Julia package for copula routines that I have developed and continue to maintain. The focus is on the Gumbel-Barnett Archimedean generator, listed as number 4.1.9 in Nelsen’s catalog. The Archimedean copulas constructed from this generator have been known for some time, but are usually restricted to the bivariate case. Its generator is given as follows, with a parameter $\theta$:

$$ \phi(t) = \exp\left(\frac{1 - e^t}{\theta}\right). $$

For practical implementation, it is useful to have explicit formulas for the generator, its inverse, and their derivatives. After some calculations, I was able to derive these expressions.

The inverse generator is given by:

$$ \phi^{-1}(t) = \log\left(1 - \theta \log(t)\right). $$

The first derivative of the inverse is:

$$ \left(\phi^{-1}\right)'(t) = -\frac{\theta}{t \left(1 - \theta \log(t)\right)}. $$

For higher-order derivatives of the generator $\phi^{(k)}(t)$ for $k \in \mathbb{N}$, the computation becomes more involved. Let us define $C(t) = -\frac{e^t}{\theta}$. Note that $C’(t) = C(t)$, which leads to a recursive structure in the derivatives.

The first derivative is:

$$ \phi’(t) = \phi(t) C(t). $$

Applying the chain rule, the second derivative is:

$$ \phi^{(2)}(t) = \left(\phi’(t)\right)’ = \phi’(t) C(t) + \phi(t) C’(t) = \phi(t) \left(C(t)^2 + C(t)\right). $$

Continuing, the third derivative is:

$$ \phi^{(3)}(t) = \left(\phi^{(2)}(t)\right)’ = \phi(t) \left(C(t)^3 + 3C(t)^2 + C(t)\right). $$

A pattern emerges: for each $k$, the $k$-th derivative can be written as

$$ \phi^{(k)}(t) = \phi(t) \sum_{\ell=1}^k a_{k,\ell} C(t)^\ell, $$

where the coefficients $a_{k,\ell}$ satisfy the recursion

$$ a_{k,\ell} = a_{k-1,\ell-1} + \ell, a_{k-1,\ell} $$

with $a_{1,1} = 1$ and $a_{k,0} = 0$. This is easily proved with the chain rule and the product rule.

This result not only provides explicit formulas for implementation but also highlights an elegant connection between copula theory and combinatorics, as this recursion corresponds to Stirling numbers of the second kind! We therefore obtain the following closed-form formula for them:

$$a_{k,\ell} = \frac{1}{\ell!} \sum_{i=0}^\ell (-1)^{\ell -i} \binom{\ell}{i} i^k.$$

Digging a bit more, in fact the polynomial $B_k$ that we want, such that $$B_k(C(t)) = \sum_{\ell=1}^k a_{k,\ell} C(t)^{\ell},$$

is simply called Bell’s polynomial. So our derivative simply is

$$\phi^{(k)}(t) = \exp\left(\frac{1}{\theta} + C(t)\right) * B_k(C(t),…C(t))$$

But now the application of $B_k$ at a vector $(x,…,x)$ is called Touchard’s polynomials $T_k$:

$$\phi^{(k)}(t) = \exp\left(\frac{1}{\theta} + C(t)\right) * T_k(C(t)).$$

But the recursive version is much more interesting for computations, don’t you think? We ended up with the following Julia code to define the Gumbel-Barnett copula:

ϕ(G::GumbelBarnettGenerator, t) = exp((1 - exp(t)) / G.θ)
ϕ⁽¹⁾(G::GumbelBarnettGenerator, t) = -exp(1 - exp(t) / (G.θ + t)) / G.θ
ϕ⁻¹(G::GumbelBarnettGenerator, t) = log1p(-G.θ * log(t))
ϕ⁻¹⁽¹⁾(G::GumbelBarnettGenerator, t) = -G.θ / (t - G.θ * t * log(t))
function ϕ⁽ᵏ⁾(G::GumbelBarnettGenerator, ::Val{k}, t) where k
    A = zeros(k)
    A[1] = 1
    for j in 2:k
        A[2:j] = A[1:j-1] .+ collect(2:j) .* A[2:j]
    end 
    return ϕ(G,t) * dot(A, (-exp(t)/G.θ).^(1:k))
end

Note the use of Val{k}() in our interface that is supposed to make as many computations as possible happen at compile time… I am sure we could have this recursion there too with a bit more work. Not today.

Part 2: Which dimensions can the Gumbel-Barnett copula really have?

It is well known that the Archimedean construction yields a true copula in dimension $d$ if and only if the generator is $d$-monotone. Is it the case here? The generator has a parameter $\theta$, what are its potential values?

Let us look at the signs of the derivatives.

1. The generator itself is always positive.

2. Its derivative $\phi’(t) = \phi(t) C(t)$ is clearly negative if and only if $\theta > 0$.

3. The second one, however… $\phi’'(t) = \phi(t) \left(C(t) + C^2(t)\right)$… will be positive if $C(t)^2 + C(t)$ is. Since $C(t) <0$, we require $1 + C(t) < 0$ which means $1 - \frac{e^t}{\theta} < 0$ or $t \ge \log(\theta)$ for all $t \in \mathbb R_+$, hence $\theta \le 1$

So the bivariate model only works for $\theta \in [0,1]$, which is an already known result.

What about higher-dimensional models? Let us look at the sign of the next derivatives, recalling that

$$\phi^{(k)}(t) = \exp\left(\frac{1}{\theta} + C(t)\right) T_k(C(t))$$

These functions are positive if and only if $T_k(C(t))$ is. Recall that Touchard’s polynomials have all their roots real and negative, while $C(t)$ is negative, so this does not look very good.

At least $C(t) = -\frac{e^t}{\theta} < \frac{-1}{\theta}$ for all $t \in \mathbb R_+$, that is certain. If we can guarantee that there are no roots of $T_k$ positioned to the left of $\frac{-1}{\theta}$, the sign of $T_k$ would be stabilized after this last root, and the associated $k$-dimensional Archimedean copula will exist.

Let us therefore look at Touchard’s roots. Here are the first few leftmost roots:

So for example, in dimension 3, we would require $\frac{-1}{\theta} < \frac{-3 - \sqrt{5}}{2}$ which means $\theta < \frac{2}{3+\sqrt{5}}$.

In greater dimensions, getting exact roots is complicated. Fortunately, Wikipedia exists and highlights the following lower bound on the roots’ positions: for every root $x^*$ of $T_k$, we have:

$$\lvert x^*\rvert > E(n) = \frac1n\binom{n}{2}+\frac{n-1}{n}\sqrt{\binom{n}{2}^2-\frac{2n}{n-1}\left(\binom{n}{3}+3\binom{n}{4}\right)}.$$

The first values of this bound are:

Hence, we can accept parameters $\theta \in [0, 1/E(n)]$ for each dimension $n$, and finally complete our Gumbel-Barnett implementation by setting the max_monotony of the generator as follows:

function _inv_E(n::Int)
    n==2  && return 1.0
    n==3  && return 0.380
    n==4  && return 0.216
    n==5  && return 0.145
    n==6  && return 0.106
    n==7  && return 0.082
    n==8  && return 0.066
    n==9  && return 0.055
    n==10 && return 0.046
    C2 = binomial(n, 2)
    C3 = binomial(n, 3)
    C4 = binomial(n, 4)
    term1 = C2 / n
    discr = C2^2 - (2n / (n - 1)) * (C3 + 3*C4)
    term2 = (n - 1)/n * sqrt(discr)
    return 1/(term1 + term2)
end

function max_monotony(G::GumbelBarnettGenerator)
    G.θ == 0 && return Inf
    n = 3
    last_valid_n = 2 # minimum 2. 
    while n <= 1e3 # we look until 1000, with E(1000) \approx 20 000 so \theta < 0.00005
        val = _inv_E(n)
        if val > G.θ
            last_valid_n = n
            n += 1
        else
            break
        end
    end
    return last_valid_n
end

end

Recall from Copulas.jl’s docs that this max_monotony function is part of the interface and must be specified. In fact, specifying only ϕ(G, t) and max_monotony(G) is enough for the object G to behave as a valid Archimedean generator; this is how pleasing our interface can be sometimes…

Anyway, we now have a version of the Gumbel-Barnett Archimedean Copula that is not restricted to two dimensions, but which can handle more dimensions—if and only if its parameter is small enough, of course. I believe that this is the only implementation of multivariate Gumbel-Barnett copulas around… It will be released in Copulas.jl soon :)

PS: If you want to hear more about how I’m trying to invert the kth derivative of the generator to make the Rosenblatt transformation work, you can take a look at this discourse thread.