Practicum 13

Ce TP servira à illustrer les différentes fonctions du progiciel R copula, par le biais d’une analyse d’un modèle de copule sur le jeu de données LOSS-ALAE. Pour des illustrations de plusieurs copules, consultez la copulathèque.

Les données sont tirées de Frees, E & Valdez A.(1998). Understanding relationships using copulas, North American Actuarial Journal, 2, 1–26.

1500 réclamations générales choisies aléatoirement à partir du bureau des services d’assurance

Notez que 34 observations sont censurées à droite, c’est-à-dire que la réclamation égale le montant limite de la police d’assurance. Shih & Louis (1995) étendent la méthode de pseudo-vraisemblance à des données bivariées censurées, en utilisant l’estimateur de Kaplan-Meier et nous procédons similairement. Le diagramme de Kendall illustre sans surprise la dépendence positive entre réclamations et provisions.

library(copula)
data(loss)
n <- nrow(loss)
#Nuage de point des données, valeurs censurées en rouge
plot(log(loss$loss), log(loss$alae),pch=20,xlab="log des paiements d'indemnités", 
     ylab="log des provisions ajustées pour sinistres survenus",main="Données LOSS-ALAE")
points(log(loss$loss[which(loss$censored==1)]),log(loss$alae[which(loss$censored==1)]),col=2,pch=20)

#Diagramme de Kendall
par(pty="s")
lcopula::K.plot(loss[,1:2],pty="s")

#Calcul du tau de Kendall et du rho de Spearman
#tau <- cor(loss[,1:2], method="kendall")[1,2]
tau <- pcaPP::cor.fk(loss[,1:2])[1,2] 
#Cet algorithme est de complexité O(nlog(n)), comparativement à l'implémentation classique qui est O(n^2)
rho <- cor(loss[,1:2], method="spearman")[1,2]
tau; rho
## [1] 0.3154175
## [1] 0.451872
#Présence de doublons en plus des données censurées: 65 pour ALAE, 925 pour LOSS
U.sp <- loss[,1:2] #copie de la matrice de données
U.sp[,2] <- rank(loss[,2],ties.method="random")/(n+1) #pseudo-observations pour ALAE
library(survival)
#Estimateur de Kaplan-Meier pour les données censurées à droite
KM <- survfit(Surv(time=loss$loss, time2=loss$limit, type="interval",event=-loss$censored+1)~1)
censoredECDF <- stepfun(x=KM$time,y=c(0,1-KM$surv))
U.sp[,1] <- censoredECDF(loss$loss)/(n+1)*n
U.sp <- as.matrix(U.sp)

La fonction pobs de copula automatise le processus, et prend une matrice comme argument. Notez qu’en raison des doublons, il peut être nécessaire de considérer les égalités en fournissant une méthode à l’argument ties.method de rank.

U.sp <- pobs(data,ties.method="random")
#Pas approprié ici en raison des données censurées

Les données LOSS-ALAE sont fortement asymétriques et ont une queue lourde. Deux distributions typiques utilisées en actuariat pour ce type de données sont les distributions de Burr et de Pareto. L’estimation multivariée est possible en utilisant la fonction mvdc, qui sert à la création d’une fonction multivarié construite à partir d’une copule. Comme pour la plupart des objets de classe copula, des arguments doivent êtres fournis (qu’ils soient optionnels ou pas). On peut aussi estimer par la méthode des moments le paramètre; dans la fonction iTau, on calibre le paramètre de gumbelCopula au \(\tau\) de Kendall observé dans l’échantillon. Les paramètres sont tirés du livre de Joe (2014).

#Modèles paramétriques tirées de actuar
library(actuar,warn.conflicts=F)
#Paramètre correspond à la méthode de moments pour la copule de Gumbel
theta.hat <- iTau(gumbelCopula(), tau)
LOSS.Mvd <- mvdc(copula = gumbelCopula(2, dim=2),
                 margins = c("pareto","burr"), paramMargins = list(list(shape = 1.33, scale = 1.66),
                                                                   list(shape1 = 2, shape2=1, scale = 1.5)))

start <- c(1.33,1.66,2,1,1.5,theta.hat)
bivfit <- fitMvdc(data=as.matrix(loss[,1:2]/10000), mvdc=LOSS.Mvd, start = start)
print(bivfit)
## The Maximum Likelihood estimation is based on 1500 observations.
## Margin 1 :
##          Estimate Std. Error
## m1.shape    1.226      0.068
## m1.scale    1.593      0.147
## Margin 2 :
##           Estimate Std. Error
## m2.shape1   1.5980      0.163
## m2.shape2   1.1165      0.038
## m2.scale    0.9488      0.131
## Copula:
##       Estimate Std. Error
## param    1.444      0.033
## The maximized loglikelihood is -4505.904 
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient 
##       51       14

Lorsque l’on utilise la fonction optim, il est de convenance de passer les arguments transformés pour que le domaine d’optimisation soit non-contraint. Cela aide l’algorithme si certains paramètres sont situés sur la bordure de l’espace admissible (en insérant exp si les paramètres doivent être positifs, par exemple. Dans ce cas, les erreurs standards doivent être ajustées par la méthode \(\Delta\). Notez que par défaut, l’optimisation correspond à minimiser la fonction objective. Vous pouvez maximiser en fournissant l’argument control=list(fnscale = -1). Le vecteur de départ de l’algorithme est souvent spécifié à l’aide d’estimateurs de moments, qui sont rapides à calculer. Dans le cas présent, on utilise les EMV issus d’un calcul prélimaire de l’ouvrage de référence. Pour les fonctions standards, une fonction enveloppe pour optim est fitdistr.

La distribution logistique est un cas spécial de la loi de Burr quand le paramètre de forme shape2 est 1. Puisque les lois sont imbriquées, un test de ratio de vraisemblances est possible. La log-vraisemblance est presque identique pour les deux modèles.

#Définitions de fonctions de log-vraisemblance (-ll) pour les lois Burr et Pareto, suggérées par Joe (2015)
loglik.marg.bur <- function(b, x){ -sum(dburr(x, shape1= b[1], shape2=b[2], scale= b[3], log=T))}
#Estimateurs de maximum de vraisemblances, variables censurées
#Notez que l'estimateur paramétrique est légèrement différent du fait que les données soient censurées
loglik.marg.par.cens  <- function(b, x, y){
  -sum(dpareto(x[which(y==0)], shape=b[1], scale=b[2], log=T))-
    sum(log(1-ppareto(x[which(y==1)], shape=b[1], scale=b[2])))
  }
loglik.marg.bur.cens <- function(b, x, y){
  -sum(dburr(x[which(y==0)], shape1= b[1], shape2=b[2], scale= b[3], log=T))-
    sum(log(1-pburr(x[which(y==1)], shape1= b[1], shape2=b[2], scale= b[3])))
  }

#Division par 10000 pour préserver le paramètre d'échelle dans un intervalle acceptable
opt1burr <- optim(par=c(1.32, 1.67, 1),fn=loglik.marg.bur.cens, x=loss$loss/10000, y=loss$censored, hessian=T)
opt1par  <- optim(par=c(1.32, 1.67),fn=loglik.marg.par.cens, x=loss$loss/10000,y=loss$censored,hessian=T)
#LRT: calcul de la valeur p
pchisq(2*(opt1par$value-opt1burr$value),df=1)
## [1] 0.7514158
#On ne rejette pas l'hypothèse nul que le modèle simplifié est adéquat
par1 <- opt1par$par
opt2burr <- optim(par=c(2.12, 1.10, 1.42),fn=loglik.marg.bur, x=loss$alae/10000, hessian=T)
par2 <- opt2burr$par
#Erreurs standard
se.par1 <- sqrt(diag(solve(opt1par$hessian)))
se.par2 <- sqrt(diag(solve(opt2burr$hessian)))
param <- rbind(c(par1,par2),c(se.par1,se.par2))
colnames(param) <- c("shape","scale","shape1","shape2","scale")
rownames(param) <- c("Parameter","St.error")
#Paramètres et erreurs standards pour les lois marginales
param
##                shape     scale    shape1    shape2     scale
## Parameter 1.13459895 1.4440678 1.7047344 1.0990071 1.0334128
## St.error  0.06602523 0.1386276 0.2008986 0.0396571 0.1657517
#Transformation des marges à l'échelle uniforme
U.p <- cbind(ppareto(q=loss$loss/10000, shape=par1[1], scale=par1[2]),
             pburr(q=loss$alae/10000, shape1=par2[1],shape2=par2[2],scale=par2[3]))
par(mfrow=c(1,2),pch=20,pty="s",mar=c(4.5,5,1,1))
plot(U.p, xlab=expression(widehat(U)[loss]),ylab=expression(widehat(U)[alae]),main="Paramétrique")
plot(U.sp, xlab=expression(widehat(U)[loss]),ylab=expression(widehat(U)[alae]),main="Semiparamétrique")

On procède maintenant à l’estimation du paramètre de la copule. Ici, on a deux choix: fournir les pseudo-observations (nécessaire en raison de la censure et des duplicatas) ou encore transformer non-paramétriquement les données en pseudo-observations en utilisant l’inférence pour les marges (mpl), via la fonction pobs. L’utilisation de cette méthode dans le contexte est uniquement illustrative. Une façon heuristique de juger de notre modèle est de comparer la dépendence mesurée et celle du modèle postulé. Si les estimés sont radicalement différents, c’est mauvais signe.

#Inférence pour les marges (IFM)
fit.ipm <- fitCopula(copula=LOSS.Mvd@copula, data=U.p, method="ml", start = theta.hat)
#Estimation de la pseudo-vraisemblance
fit.pv <- fitCopula(copula=gumbelCopula(), data=U.sp, method="mpl")
#Essai avec un autre modèle de valeurs extrêmes
fit.pv2 <- fitCopula(copula=galambosCopula(),data=U.sp, method="mpl")
#Imprimer les résultats de l'optimisation
fit.ipm
## fitCopula() estimation based on 'maximum likelihood'
## and a sample of size 1500.
##       Estimate Std. Error z value Pr(>|z|)    
## param  1.45417    0.02918   49.84   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## The maximized loglikelihood is  206.2 
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient 
##       16        3
fit.pv
## fitCopula() estimation based on 'maximum pseudo-likelihood'
## and a sample of size 1500.
##       Estimate Std. Error z value Pr(>|z|)    
## param  1.45821    0.03284   44.41   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## The maximized loglikelihood is  208 
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient 
##       15        3
fitted <- mvdc(copula=fit.pv@copula,margins = c("pareto","burr"), 
               paramMargins = list(list(shape=par1[1], scale=par1[2]), 
                                   list(shape1=par2[1],shape2=par2[2],scale=par2[3])))
#Valeur de tau et de rho pour le modèle estimé
tau(fit.pv@copula); rho(fit.pv@copula)
## [1] 0.3142272
## [1] 0.4498425
par(mfrow=c(1,1),pch=20)
plot(log(rMvdc(fitted, n=n)),xlab="LOSS simulées", ylab="ALAE simulées")

#Densité de la copule de Gumbel - 
sum(dCopula(u=U.p, copula = fit.pv@copula, log=T))
## [1] 206.2315

Il reste maintenant à tester l’adéquation de noter modèle, sélectionné il faut dire sur la base de quelques faits empiriques et au vu des modèles sélectionnés dans les articles qui ont analysé le jeu de données LOSS-ALAE. Malheureusement, la plupart des tests sont basés sur les pseudo-observations et une des hypothèse pour leur validité ne tient pas (dans le cas présent la continuité absolue de la loi des données). Nous illustrons les tests en reprenant les observations (moins les données censurées) et en ordonnant de façon aléatoire les égalités dans les données.

Le code ci-dessous teste l’adéquation du modèle au données, l’indépendence des données et l’hypothèse qui veut que la copule sous-jacente soit une de valeur extrême.

#Test d'adéquation basé sur l'autoamorçage par la méthode des multiplicateurs
#Tiré de Yan et Kojadinovic (2010), JSS
Loss.ncens <- subset(loss, censored == 0, select = c("loss", "alae"))
pseudoLoss <- sapply(Loss.ncens, rank, ties.method = "random") /(nrow(Loss.ncens) + 1)

#Test d'indépendence
empsamp <- indepTestSim(nrow(pseudoLoss), p = 2, N = 1000, verbose=F)
indepTest(pseudoLoss, empsamp)
## 
## Global Cramer-von Mises statistic: 2.80509 with p-value 0.0004995005 
## Combined p-values from the Mobius decomposition:
##    0.0004995005  from Fisher's rule,
##    0.0004995005  from Tippett's rule.
#Rejet clair de l'hypothèse nulle, évident visuellement
#Test d'adéquation
lossGof.gumbel.mult <- gofCopula(gumbelCopula(), pseudoLoss, method="Sn",estim.method = "mpl", simulation = "mult")
lossGof.gumbel.mult
## 
##  Multiplier bootstrap goodness-of-fit test with 'method'="Sn",
##  'estim.method'="mpl"
## 
## data:  x
## statistic = 0.026482, parameter = 1.4248, p-value = 0.1334
#gofTest <- gofCopula(copula=fit.pv@copula, x=as.matrix(loss[,1:2]), N = 1000, estim.method="itau")
#Nous avons postulé une copule de valeur extrême. Cette hypothèse est-elle valide?
evTestK(pseudoLoss)
## 
##  Test of bivariate extreme-value dependence based on Kendall's
##  process with argument 'method' set to "fsample"
## 
## data:  pseudoLoss
## statistic = -0.073288, p-value = 0.9416
#Pas d'évidence pour rejeter

Pour représenter graphiquement ou créer des objets de class copule, on peut utiliser par exemple les commandes suivantes

par(mfrow=c(1,2), pty="s",pch=20)
nc <- normalCopula(param=0.5)
plot(rCopula(1000,nc),xlab="U",ylab="V")
dCopula(c(0.3,0.4), nc)
## [1] 1.192296
pCopula(c(0.3,0.4), nc)
## [1] 0.1918907
contour(nc, dCopula, xlim=c(0.05,0.95),
        ylim=c(0.05,0.95))

graphics.off()
persp(nc, dCopula)

Question 4

  1. Charger les données, faites une analyse exploratoire (représentation graphique et statistiques descriptives) et tester pour la dépendance extrémale et l’échangeabilité).
#Charger les données
library(lcopula)
data(nutrient)
dat <- as.matrix(nutrient[-690,3:4])

Vous pouvez utiliser les fonctions evTestK et exchTest pour tester les hypothèses.

  1. Choisir un modèle adéquat pour les données.

Les classes evCopula, archmCopula et ellipCopula contiennent des copules de valeurs extrêmes, archimédiennes et elliptiques. D’autres copules sont disponibles (par exemple fgmCopula ou plackettCopula). Une plus grande variété est disponible en chargeant le paquet VineCopula. Utiliser la copulathèque pour vous guider dans votre choix. Les copules de valeurs extrêmes ont de la dépendence extrémale, les copules archimédiennes sont échangeables, les copules elliptiques ont la même dépendence de queue inférieure et supérieure, etc.

La fonction BiCopSelect du paquet VineCopula estime une série de modèles et retourne le meilleur choix selon un critère d’information. Cette façon de faire n’est pas justifiée, mais peut vous aider à choisir un modèle.

  1. Estimer les paramètres de la copule en assumant un modèle paramétrique pour les marges et en utilisant la technique d’inférence pour les marges ou de maximum de vraisemblance.

  2. Estimer les paramètres à l’aide du pseudo-maximum de vraisemblance ou de la méthode des moments.

Choisissez une famille pour les distributions marginales (par exemple lognormale ou gamma). Vérifier l’adéquation à l’aide d’un diagramme QQ. La fonction MASS::fitdistr permet d’estimer quelques modèles courants. Autrement, utiliser optim, avec par exemple pour la distribution mydist la fonction -sum(dmydist(..., log=TRUE)).

Les fonctions fitCopula avec method respectivement ml pour l’inférence pour les marges (vous devez d’abord transformer les données à l’aide de l’estimé de votre modèle), mpl pour les données pobs(), itau et irho pour la méthode des moments (copules archimédiennes et elliptiques). Vous pouvez aussi utiliser fitMvdc pour l’estimation paramétrique en une étape.

  1. Créer un échantillon du modèle estimé et comparer graphiquement sur \([0,1]^2\) ce dernier aux pseudo-observations. Répéter le diagnostic graphique, cette fois à l’échelle originale des données.

La fonction rCopula() vous permettra de générer des données du modèle. Utiliser la transformation des quantiles et la transformées intégrale de probabilité pour passer d’une échelle à l’autre.

  1. Superposer aux données un contour correspondant à la densité du modèle estimé.

Voir contour et l’exemple des données LOSS-ALAE.

  1. Faire un test d’adéquation.

La fonction gofCopula suffira (il y a des variantes selon votre choix de copules). L’option simulation="mult" peut être utilisée pour accélérer la procédure d’autoamorçage, attendu que la taille de l’échantillon est suffisante ici.