Test the copula

from __future__ import print_function
import openturns as ot
import openturns.viewer as viewer
from matplotlib import pylab as plt

Copula fitting test using Kendall plot

We first perform a visual goodness-of-fit test for a copula with the Kendall plot test.

We create two samples of size 500.

N = 500
dist1 = ot.ComposedDistribution([ot.Normal()] * 2, ot.GumbelCopula(3.0))
sample1 = dist1.getSample(N)
dist2 = ot.ComposedDistribution([ot.Normal()] * 2, ot.ClaytonCopula(0.2))
sample2 = dist2.getSample(N)

We change the parameter for the evaluation of E(Wi) thanks to the ResourceMap :

ot.ResourceMap.SetAsUnsignedInteger('VisualTest-KendallPlot-MonteCarloSize', 25)

We can test a specific copula model for a given sample,

copula_test = ot.GumbelCopula(3)
graph = ot.VisualTest.DrawKendallPlot(sample1, copula_test)
view =viewer.View(graph)
Kendall Plot

or test whether two samples have the same copula model

graph = ot.VisualTest.DrawKendallPlot(sample1, sample2)
view = viewer.View(graph)
Kendall Plot

The first test gives a positive result as the blue curve is near the identity line which is not the case for the second test.

Graphical copula validation

In this paragraph we visualize an estimated copula versus the data in the rank space.

First we create data

marginals = [ot.Normal()] * 2
dist = ot.ComposedDistribution(marginals, ot.ClaytonCopula(3))
N = 500
sample = dist.getSample(N)

We build a estimate copula from the previous sample using the ClaytonCopulaFactory :

estimated = ot.ClaytonCopulaFactory().build(sample)

We represent data as a cloud in the rank space :

ranksTransf = ot.MarginalTransformationEvaluation(marginals, ot.MarginalTransformationEvaluation.FROM)
rankSample = ranksTransf(sample)
rankCloud = ot.Cloud(rankSample, 'blue', 'plus', 'sample')

We can plot the graph with rank sample and estimated copula :

myGraph = ot.Graph('Parametric estimation of the copula', 'X', 'Y', True, 'topleft')

and draw the iso-curves of the estimated copula

minPoint = [0.0]*2
maxPoint = [1.0]*2
pointNumber = [201]*2
graphCop = estimated.drawPDF(minPoint, maxPoint, pointNumber)
contour_estCop = graphCop.getDrawable(0)
# Erase the labels of the iso-curves
# Change the levels of the iso-curves
nlev = 21
levels = [0.0] * nlev
for i in range(nlev):
    levels[i] = 0.25 * nlev / (nlev - i)
# Change the legend of the curves
contour_estCop.setLegend('Gumbel copula')
# Change the color of the iso-curves
# Add the iso-curves graph into the cloud one
view = viewer.View(myGraph)
Parametric estimation of the copula

Display figures


Total running time of the script: ( 0 minutes 0.370 seconds)

Gallery generated by Sphinx-Gallery