Freakonometrics

Je poursuis le billet que j'avais commencé il y a quelque temps (ici) sur les copules. En effet, un étudiant me dit avoir utilisé "l'algorithme de Devroye" qui est donné comme exercice dans le livre de Roger Nelsen,

• Un algorithme de génération de la copule de Clayton

Je n'aime pas faire les exos à la place des élèves (surtout quand ce ne sont pas les miens), il faut faire un peu de calculs. En fait, c'est comme ça que Cook et Johson ont obtenu, en 1981, une distribution qui correspondait avec celle obtenue par Clayton.

On notera que ce n'est pas l'utilisation de la représentation frailty de la copule de Clayton...

• Travail sur un algorithme faux, et tests d'ajustement

Malheureusement, dans l'algorithme envoyé, l'étudiant avait permuté les facteurs de scale et de shape de la loi Gamma, ce qui donne des choses assez différentes. Mais sur lesquelles je vais rebondir. On va faire comme si on n'avait pas vu la faute dans l'algorithme, et que l'on se demande si on a généré la copule de Clayton que l'on pense avoir simulée....
Considérons l'algorithme suivant, qui génère effectivement une distribution à marges dans [0,1], on se demande si la copule associé est, ou pas, la copule de Clayton (donc on prend les rangs)
> n=10000
> x1=rexp(n,1)
> x2=rexp(n,1)
> theta=1/2
> x=rgamma(n, shape=1, scale=theta)
> u1=(1+x1/x)^(-theta)
> u2=(1+x2/x)^(-theta)
> x1=rank(u1)/(n+1)
> x2=rank(u2)/(n+1)
> plot(x1,x2)
Vu de loin, cela ressemble à une copule de Clayton.

On peut d'ailleurs essayer de récupérer la densité de la copule,
>  beta.kernel.copula.surface <- function (u,v,bx,by,p) {
+  s <- seq(1/p, len=(p-1), by=1/p)
+  mat <- matrix(0,nrow = p-1, ncol = p-1)
+ for (i in 1:(p-1)) {
+    a <- s[i]
+    for (j in 1:(p-1)) {
+        b <- s[j]
+    mat[i,j] <- sum(dbeta(a,u/bx,(1-u)/bx) * dbeta(b,v/by,(1-v)/by)) / length(u)
+    } }
+ return(data.matrix(mat)) }
> Z= beta.kernel.copula.surface(x1,x2,bx=.1,by=.1,p=26)
>  p=26
> u=seq(1/p, len=(p-1), by=1/p)
> persp(u,u,Z,theta=30,col="green",shade=TRUE)

ou en changeant la fenêtre de l'estimation à noyau
> Z= beta.kernel.copula.surface(x1,x2,bx=.01,by=.01,p=26)

Bref, ça ressemble à du Clayton.... Si on regarde la section diagonale de la copule empirique, on note que l'on n'obtiens ni la copule de Clayton de paramètre 2, ni celle de paramètre 1/2,

> U=seq(.01,1,by=.01)
> D=rep(NA,length(U))
> for(i in 1:length(U)){
+ D[i]=mean((x1<U[i])&(x2<U[i]))
+ }
> U=c(0,U)
> plot(U,c(0,D),cex=.5)
> lines(U,(2*U^(-1/theta)-1)^(-theta),col="red")
> lines(U,(2*U^(-theta)-1)^(-1/theta),col="blue")
Bref, ça ne semble pas correspondre à une des deux copules de Clayton recherchées,

On peut aussi des grandeurs comme la fonction
qui est l'extension en dimension 2 de l'integral probability transform (cf ici par exemple). Ou plutôt la fonction  définie comme
> V=rep(NA,n)
> for(i in 1:n){
+ V[i]=sum((x1<x1[i])&(x2<x2[i]))/(n-1)
+ }
> lambda=sort(V)-(1:n)/(n+1)
> plot(sort(V),lambda)
> u=seq(0,1,by=0.01)
> l=(u^(-theta)-1)/(-theta*u^(-theta-1))
> lines(u,l,col="blue")
> l=(u^(-1/theta)-1)/(-1/theta*u^(-1/theta-1))
> lines(u,l,col="red")

Pour conclure définitivement, je reprendrais des choses que j'avais faites dans ma thèse. En fait, la copule de Clayton est la seule copule continue qui soit invariante par troncature. On peut alors regarder les tau de Kendall ou Rho de Spearman conditionnels. Normalement, ils devraient être constants,
> U=seq(.01,1,by=.01)
> R=K=rep(NA,length(U))
> for(i in 1:length(U)){
+ I=(x1<U[i])&(x2<U[i])
+ R[i]=cor(x1[I],x2[I],method = "spearman")
+ K[i]=cor(x1[I],x2[I],method = "kendall")
+ }
> plot(U,R,cex=.75,type="b")
> abline(h=R[length(R)],col="green")

> plot(U,K,cex=.75,type="b")
> abline(h=K[length(R)],col="purple")

Moralité (même si je n'ai pas tracé les régions de confiance des tests) je pense pouvoir conclure que l'algorithme tel qu'il est là ne génère pas une copule de Clayton de paramètre 2 ou 1/2.