Is it possible to compose a copula of multiple variables where a single variable is dependent on multiple other variables?

Python usage
1

Suppose I have 3 variables x1, x2, x3. x1 and x2 are correlated by means of gumbel copula. x1 and x3 are correlated by means of another gumbel copula. x2 is independent of x1 and x3.

(my actual case has 2 more variables, which are both independent of x1, x2, and x3. But I left that out for this example.)

How would I define a composed copula in that case? I already found the ComposedCopula class, but combines independent copula, which is clearly not the case in my example (because of the shared variable x1).

Code I have so far based on ComposedCopula:

x1 = ot.Normal(0, 1)
x2 = ot.Normal(1, 2)
x3 = ot.Normal(2, 4)
cop1 = ot.GumbelCopula(2)
cop2 = ot.GumbelCopula(4)
jointcop = ot.ComposedCopula([cop1, cop2])
compdist = ot.ComposedDistribution([x1, x2, x1, x3], jointcop)

This leads to a composed copula with dimension 4, where the first and third dimension refer to x1. Which, are independent of eachother. So clearly this is not correct…

2

Hi,

You cannot have X1 and X2 linked by a Gumbel copula and X2 independent of X1 (excepted if theta=1 in the Gumbel copula, in which case it is the independent copula). A copula build a symmetrical link between random variables: you cannot have X1 linked to X2 in a given way and X2 linked to X1 in another way.

Could you double check the description of your need? I will be happy to help.

Cheers

Régis

3

Hi Régis,

Thank you for your reply! And you are absolutely right, I made a mistake in my explanation. What I meant to say is:

  • x1 and x2 are correlated by means of gumbel copula
  • x1 and x3 are correlated by means of another gumbel copula
  • x2 is independent of x3

Put in a different way, this is what I’m trying to achieve in a table form (g1 is the dependency by means of the first gumbel copula, g2 the second):

     +----+----+----+
     | x1 | x2 | x3 |
+====+====+====+====+
| x1 | 1  | g1 | g2 |
+----+----+----+----+
| x2 | g1 | 1  | 0  |
+----+----+----+----+
| x3 | g2 | 0  | 1  |
+----+----+----+----+
4

to elaborate on my first reply, with the code with composedcopula:

x1 = ot.Normal(0, 1)
x2 = ot.Normal(1, 2)
x3 = ot.Normal(2, 4)
cop1 = ot.GumbelCopula(2)
cop2 = ot.GumbelCopula(4)
jointcop = ot.ComposedCopula([cop1, cop2])
compdist = ot.ComposedDistribution([x1, x2, x1, x3], jointcop)

I get a dependency structure similar to the following

     +----+----+----+----+
     | x1a| x2 | x1b| x3 |
+====+====+====+====+====+
| x1a| 1  | g1 | 0  | 0  |
+----+----+----+----+----+
| x2 | g1 | 1  | 0  | 0  |
+----+----+----+----+----+
| x1b| 0  | 0  | 1  | g2 |
+----+----+----+----+----+
| x3 | 0  | 0  | g2 | 1  |
+----+----+----+----+----+

where x1a and x1b should be the same variable (or perfectly correlated), but are now independent.

5

Sorry for the delay, I was pretty busy these last weeks!
I better understand your problem, and it appears to be highly nontrivial. It is not a matter of how to do it with OpenTURNS, but rather if it is possible to do it at all (from a mathematical point of view). And the bad news is that you cannot have (x1,x2) linked by a Gumbel copula, (x1,x3) linked by a Gumbel copula and (x2,x3) independent. The demonstration is quite simple:

  1. If x2 and x3 are independent, in particular they are independent given x1
  2. In this case, the conditional copula c(x1,x2,x3|x1)=c(x1,x2,x3)/c(x1)=c(x1,x2,x3) as c(x1)=1 (an univariate copula is always the uniform distribution over [0,1]) is the product of the two conditional copulas c(x1,x2|x1)=c(x1,x2)/c(x1)=c(x1,x2) and c(x1,x3|x1)=c(x1,x3)/c(x1)=c(x1,x3)
  3. So we must have c(x1,x2,x3)=c(x1,x2)c(x1,x3)=Gumbel(theta12).pdf(x1,x2)xGumbel(theta13).pdf(x1,x3)
  4. The copula of (x2,x3) is then \int_0^1c(x1,x2)c(x1,x3)dx1, which is not equal to 1 in general (except if theta12=theta13=1

You can play with the following script to see what’s happen when the parameters change:

import openturns as ot

class ComplexCopula(ot.PythonDistribution):
    def __init__(self, copula01, copula02, N=128):
        super(ComplexCopula, self).__init__(3)
        self.copula01_ = copula01
        self.copula02_ = copula02
        self.range_ = ot.Interval(3)
        self.algo_ = ot.GaussLegendre([N]*3)

    def getRange(self):
        return self.range_

    def getRealization(self):
        X01 = self.copula01_.getRealization()
        X2 = self.copula02_.computeConditionalQuantile(ot.RandomGenerator.Generate(), [X01[0]])
        return [X01[0], X01[1], X2]

    def computeCDF(self, X):
        interval = ot.Interval([0.0]*3, X).intersect(self.range_)
        if interval.isEmpty():
            return 0.0
        if interval == self.range_:
            return 1.0
        
        def kernel(X):
            return [self.computePDF(X)]
        
        integrand = ot.PythonFunction(3, 1, kernel)
        return self.algo_.integrate(integrand, interval)[0]

    def computePDF(self, X):
        X01 = [X[0], X[1]]
        X02 = [X[0], X[2]]
        return self.copula01_.computePDF(X01) * self.copula02_.computePDF(X02)

    def isCopula(self):
        return True

x1 = ot.Normal(0, 1)
x2 = ot.Normal(1, 2)
x3 = ot.Normal(2, 4)
cop1 = ot.GumbelCopula(2)
cop2 = ot.GumbelCopula(4)

c123 = ot.Distribution(ComplexCopula(cop1, cop2))
print("PDF=", c123.computePDF([0.5,0.5,0.5]))
print("CDF=", c123.computeCDF([0.5,0.5,0.5]))
print("sample=", c123.getSample(10))
sample = c123.getSample(10000)
for i in range(3):
    ot.Show(ot.HistogramFactory().build(sample.getMarginal(i)).drawPDF())
ot.Show(ot.VisualTest.DrawPairs(sample))

x123 = ot.ComposedDistribution([x1, x2, x3], c123)
print("PDF=", x123.computePDF([2,3,4]))
print("CDF=", x123.computeCDF([2,3,4]))
print("sample=", x123.getSample(10))
ot.Show(ot.VisualTest.DrawPairs(x123.getSample(1000)))

and the pair plot of the joint copula is

Figure_1
Figure_11200×961 148 KB

which shows clearly a dependence between x2 and x3 (lower right plot).

Cheers

Régis

1 Like