statsmodels: Copula density not correct for d>2

It seems that the pdf of some Archimedian Copulas are incorrect for d>2. They do not throw an error or warning in this case. The formulas look generic for any dimension at first but are restricted to d=2. At least for Gumbel and Clayton. I compared with this: https://arxiv.org/pdf/1207.1708.pdf For Frank there is a dim check (“pdf is currently only available for bivariate copula”).

from statsmodels.distributions.copula.archimedean import GumbelCopula
c=GumbelCopula()
c.pdf([.7,.7,.7],20)
> array([3.97063356e-08])

For example in this case, a high value for the density would be expected.

Also it would be good to know the source where the formulas have been taken from.

About this issue

Original URL
State: open
Created a year ago
Comments: 21 (21 by maintainers)

Most upvoted comments

parking this here for now

Sterling code is from Rosetta Code converted to class for cache polylog Li code formulas from Wikipedia (I also have a brute force polynomial series function for Li, but finding truncation point is non-trivial, and not needed for negative integer values)

n = k_dim, but Li argument is -n (defined for first argument is negative integer)

class Sterling1():
    
    def __init__(self):
        self._cache = {}
    
    def __call__(self, n, k):
        key = str(n) + "," + str(k)

        if key in self._cache.keys():
            return self._cache[key]
        if n == k == 0:
            return 1
        if n > 0 and k == 0:
            return 0
        if k > n:
            return 0
        result = sterling1(n - 1, k - 1) + (n - 1) * sterling1(n - 1, k)
        self._cache[key] = result
        return result
    
    def clear_cache(self):
        self._cache = {}

sterling1 = Sterling1()


class Sterling2():
    
    def __init__(self):
        self._cache = {}
    
    def __call__(self, n, k):
        key = str(n) + "," + str(k)

        if key in self._cache.keys():
            return self._cache[key]
        if n == k == 0:
            return 1
        if (n > 0 and k == 0) or (n == 0 and k > 0):
            return 0
        if n == k:
            return 1
        if k > n:
            return 0
        result = k * sterling2(n - 1, k) + sterling2(n - 1, k - 1)
        computed[key] = result
        return result

    def clear_cache(self):
        self._cache = {}
        
sterling2 = Sterling2()


def li3(z):
    return z * (1 + 4 * z + z**2) / (1 - z)**4


def li4(z):
    return z * (1 + z) * (1 + 10 * z + z**2) / (1 - z)**5


def lin(n, z):
    """Polylogarithm for negative integer
    """
    if np.size(z) > 1:
        z = np.array(z)[..., None]
        # raise ValueError("currently only for scalar z")
        
    k = np.arange(n+1)
    st2 = np.array([sterling2(n + 1, ki + 1) for ki in k]) 
    res = (-1)**(n+1) * np.sum(factorial(k) * st2 * (-1 / (1 - z))**(k+1), axis=-1)
    return res

josef-pkt on Jan 24, 2023