Arthur Charpentier

Aller au contenu | Aller au menu | Aller à la recherche

Keyword - actuariat

Fil des billets

mercredi 10 mars 2010

Inventaire et Zillmérisation

Au lieu de faire une réponse au commentaire sur les PM (ici), je vais faire un billet. La question était "cettte  provision correspond elle à la PM zillmérisée. Si non, comment la calcule-t-on ? ". Finalement, ça correspond aussi à la suite du billet d’hier sur les termes actuariels (ici). Avant de parler de la PM, il faut revenir un peu davantage sur les primes, et distinguer la prime pure (que l’on a évoquée jusqu’à présent) et la prime commerciale. Ensuite seulement, on pourra parler des calculs de PM.

  • Primes pure et commerciale
Jusqu’à présent, je n’ai pas parlé des frais de gestion ou des frais d’acquisition. Dans le second cas, il s’agit des dépenses engagées pour établir la police, et rémunérer l’agent qui rapporte le contrat. Dans le premier cas, il faut distinguer les frais durant la période où il faut encaisser les primes (période 1), et les frais durant la période où il faut garantir un paiement à l’assuré (période 2). Au cours de la première période, on trouvera les frais d’établissement des quittances, par exemple. Les deux peuvent bien entendu coïncider, auquel cas les frais s’additionne. On supposera que les frais sont proportionnels aux capitaux garantis. On peut aussi supposer qu’il y a des frais de prime commerciale (proportionnels à la prime pure)
Sous ces hypothèses, pour calculer la prime chargée (ou prime commerciale) de l’assuré, on suppose que la valeur actuelle probable des engagements de l’assuré est égale à la valeur actuelle probable des engagements de l’assureur, à laquelle on ajoute la valeur actuelle probable des frais de gestion.
Si on considère un prime unique, on a quelque chose comme
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-01.png
où http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-02.png désigne la prime commerciale unique, http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-03.png désignent les frais d’acquisition, http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-04.png correspond aux frais de gestion des engagements de l’assureur (il n’y a pas de http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-05.png car la prime est versée au début, en une seule fois),http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-06.png l’annuité associée à la garantie proposée par l’assureur (i.e. avec un capital garanti de 1) et r les chargements sur la prime. On suppose ici que tout est payé en début d’année, histoire de simplifier...
Pour une prime annuelle, en notant en minuscule les primes annuelles, et http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-07.png l’annuité associée aux engagements de l’assuré (i.e. payer sa prime annuellement tant qu’il est en vie pendant une période prédéterminée),
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-08.png
soit encore
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-50.png
Pour compliquer un peu les choses, notons qu’il existe une contre-assurance de cette prime commerciale.En effet, certains contrats d’assurance en cas de vie (disons de la retraite pour fixer les choses) où une garantie prévoie un remboursement des primes commerciales versées en cas de décès de l’assuré avant le début de la retraite. Toujours pour fixer les choses, considérons une assurance de capital différé, de durée http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-10.png, souscrite à l’âge http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-11.png. Comme auparavant, on note http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-02.png la prime unique commerciale contenant cette contre assurance, de telle sorte que
La valeur actuelle probable de la contre-assurance correspond à l’annuité d’une assurance temporaire décès de http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-10.png années, i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-12.png
qui peut s’écrire
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-13.png
soit
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-14.png
Si l’on considère une prime annuelle, payée tout au long des n années, on a quelque chose de la forme
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-15.png
où pour rappel (je renvoie au cours à ce sujet) http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-16.png désigne une assurance temporaire décès dont le capital croît d’une unité chaque année, ce qui donne
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-17.png
Ce n’est pas forcément joli, mais comme toujours en actuariat de l’assurance-vie, il suffit d’écrire les choses calmement pour retomber sur des notions qui ont été prédéfinies.
  • PM pure, d’inventaire ou zillmérisée
De la même façon qu’il existe plusieurs formes de chargements et de frais pour les prime, il existera plusieurs sortes de provisions mathématiques. Telle que je l’ai présentée, i.e. différence entre la valeur actuelle probable des engagements de l’assureur, et la valeur actuelle probable des engagements de l’assuré (des primes), il s’agit de la provision mathématique pure. Pour un contrat souscrit à l’age x, et vu après k années, je l’avais notée http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-18.png (c’est la notation officielle). Mais ce n’est pas ce que la réglementation demande de calculer.
"Les provisions mathématiques constituées par les entreprises d’assurance-vie et de capitalisation sont calculées en tenant compte, dans la détermination de l’engagement de l’assuré ou du souscripteur, de la partie des primes devant être versée par l’intéressé et représentative des frais d’acquisition du contrat, lorsque ces frais ont été portés en charge déductible par l’entreprise avant la fin de l’exercice à la clôture duquel la provision est constituée." (Article L331-1 du Code des Assurance, ici)
L’idée est qu’une fois la période de paiement de la prime (voire dès le début en cas de prime unique), la PM pure ne sera pas suffisante, l’assureur ayant davantage de frais.
Il va falloir prendre en compte différentes sortes de frais ou de chargements (évoqués auparavant). Le cas le plus simple est celui des chargements dits de règlement: si pour payer 1€ de prestation, cela coûte http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-20.png€ à l’assureur, la prime annuelle payée par l’assureur doit être multipliée parhttp://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-20.png (de telle sorte que les VAP initiale soient égales). Dans ce cas, un rapide calcul montre que la PM doit être multipliée par http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-20.png également (les engagements de part et d’autre sont multipliés par ce même facteur).
Bon, pour les frais de gestion, c’est un peu plus compliqué. Si on reprend ce qui a été fait ici, on peut montrer que la différence entre  les VAP donne une provision de la forme suivante
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-21.png
si la durée de paiement des primes était initialement http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-24.png, et que http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-10.png est la durée de la garantie. Si http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-24.png=http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-10.png, il n’y a alors pas de provision de gestion, le second terme étant nul. Cette provision existe dès lors que http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-23.png (ce qui arrive au final assez souvent me semble-t-il). On arrive ainsi à la provision d’inventaire, qui est souvent notée http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-22.png. Mais ce n’est pas fini, il y a aussi les frais d’acquisition. Généralement, ces frais sont occasionnés à la signature du contrat, puis amori tout au long de la durée du contrat, via un chargement (constant) sur la prime. Aussi,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-26.png
On pose alors
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-27.png
On baisse ainsi la PM en prenant en compte les chargements pour frais d’acquisition. Cette méthode où l’on prend en compte la valeur actuelle des frais d’acquisition dans que la PM est positive s’appelle la zillmérisation.
Voilà en gros ce que j’ai compris (je suis loin d’être un expert sur le sujet).... Car pour ceux qui en douteraient encore, l’assurance-vie n’est pas du tout mon domaine de prédilection. Promis, c’est mon avant dernier billet sur les PM (le dernier étant une ébauche de correction pour le devoir maison... mais on attendra que tout le monde l’ait rendu). Maintenant s’il y en a qui ont des compléments ou des corrections à apporter, je suis preneur ! Je peux toutefois renvoyer à un billet de Cimon sur le sujet il y a quelque temps maintenant, ici. Pour les historiens, je renvoie en revanche à une traduction (en anglais) du papier d’August Zillmer présentant cette méthode, datant de 1863 () !

mercredi 3 mars 2010

Actuariat, modélisation de la loi composée

Nous avions obtenu, ici et , des modèles pour les lois des fréquence et pour les coûts individuels. Nous avions supposé

  • http://perso.univ-rennes1.fr/arthur.charpentier/latex/tp-sum-13.png suit une loi binomiale négative (ici)
  • http://perso.univ-rennes1.fr/arthur.charpentier/latex/tp-sum-11.png suit une loi mélange de 3 lois (exponentielle tronquée, Dirac et lognormale translatée) ()
  • Utilisation des simulations paramétriques
Simuler une loi binomiale négative peut se faire en se souvenant que la loi binomiale négative est une loi Poisson-Gamma (mélange Poisson avec un paramètre Gamma). On peut alors simplement simuler une loi Gamma, puis simuler une loi de Poisson.  Mais faisons plus simple, via la fonction
> Ns=rnbinom(nsimul,.3,mu=.12)

Avant de simuler la charge totale, il nous faut quelques lignes de codes pour générer une loi exponentielle tronquée,
> rexpt=function(n,t,m){
+ ze=rexp(2*n,m)
+ z=ze[ze<t]
+ while(length(z)<n){
+ ze=rexp(2*n,m)
+ z=c(z,ze[ze<t])
+ }
+ return(z[1:n])
+ }

On peut alors enfin simuler la charge totale simplement via
> Ns=rnbinom(nsimul,.3,mu=.12)
> Z =sample(c("bleu","vert","rouge"),sum(Ns)+1,replace=TRUE,prob=c(.25,.4,.35))
> X = (Z=="vert")*1170 +
+     (Z=="bleu")*(1170+rlnorm(sum(Ns)+1,7.07,1.55)) +
+     (Z=="rouge")*rexpt(sum(Ns)+1,1170,1/400)

On commence par regarder, par police, combien il y a de sinistres, puis on simule autant de sinistres qu’il en faut. Comme on a un mélange, on simule le paramètre du mélange, i.e. le type de sinistre.
Enfin, le vecteur de la charge totale est ici obtenu via le code suivant,
> S=rep(0,nsimul)
> s=0
> for(i in (1:nsimul)[Ns>0]){
+ S[i]=sum(X[(s+1):(s+Ns[i])])
+ s=s+Ns[i]
+ }

On commence par vérifier que les ordres de grandeurs coïncident avec ce que nous dit la théorie
> Nemp=read.table("D:\\r-data\\base-nb.txt",header=TRUE)$x
> Xemp=read.table("D:\\r-data\\base-cout.txt",header=TRUE)$x
> mean(Ns)
[1] 0.120989
> mean(Nemp)
[1] 0.1216735
> mean(X)
[1] 1831.396
> mean(Xemp)
[1] 2881.173
> mean(S)
[1] 221.5804
> mean(Ns)* mean(X)
[1] 221.5788
On notera une différence significative sur les coûts moyens: on rate ici les très gros couts de sinistres.
Si l’on suppose l’ajustement valide, on en déduit alors une estimation d’un quantile élevé simplement
> quantile(S,.99)
     99%
3984.903
ou encore la probabilité qu’une police coute plus de 15000 euros sur un an,
> sum(S>15000)/length(S)
[1]  0.00179
  • Utilisation des simulations paramétriques et nonparamétriques
Comme tenu de la différence observée sur le coût moyen des sinistres individuels, on peut aussi tenter de bootrapper dans les charges observées,
> Sb=rep(0,nsimul)
> for(i in (1:nsimul)[Ns>0]){
+ Sb[i]=sum(sample(Xemp,size=Ns[i],replace=TRUE))
+ }
>  mean(Sb)
[1] 336.3616
>  mean(Nemp)*mean(Xemp)
[1] 350.5623


Cette fois, l’ordre de grandeur paraît plus acceptable (on retrouve la prime pure mentionnée ici). Et si on regarde les quantiles et les probabilités d’évnènement rare,
> quantile(Sb,.99)
     99%
3602.766
> sum(Sb>15000)/length(Sb)
[1] 0.001403
On peut visualiser ces deux distributions ci-dessous, avec en rouge la fonction de survie de variable simulée, suivant des lois paramétriques, et en bleu la fonction de survie lorsqu’on boostrappe les coûts des sinistres,
  • Approximation Normal Power
A l’aide de l’algorithme précédant, on peut récuper les premiers moments de la charge totale,
> library(fBasics)
Le chargement a nécessité le package : timeDate
Le chargement a nécessité le package : timeSeries
> mean(S)
[1] 221.5804
> var(S)
[1] 5781070
> skewness(S)
[1] 125.4405
attr(,"method")
[1] "moment"
> kurtosis(S)
[1] 32134.53
attr(,"method")
[1] "excess"

La distribution de S est alors
>  m=mean(S)
>  s=sd(S)
>  (g=as.numeric(skewness(S)))
[1] 125.4405
>  mean((S-mean(S))^3)/sd(S)^3
[1] 125.4405
> x <- seq(0,30000,length=500)
> z=(x-m)/s
> y=1-pnorm(sqrt(9/g^2 + 6*z/g + 1) - 3/g) 

qui n’est pas vraiment réaliste...
On voit que pour une police sur 20, la charge sinistre dépasse les 100000 euros. Cela n’est pas réaliste ! La moyenne devrait alors forcément excéder 5000.
  • Algorithme de Panjer
Histoire de conclure, on peut aussi regarder ce qu’aurait donné l’algorithme de Panjer. Commençons par discrétiser la distribution de la charge individuelle.
> F=function(x){
+ y=NA
+ if(x<1170){y=pexp(x,1/400)/pexp(1170,1/400)*.35}
+ if(x==1170){y=.75}
+ if(x>1170){y=.75+plnorm(x-1170,7.07,1.55)*.25}
+ return(y)
+ }
> x  = seq(0,100000,length=1000)
> Vf = rep(NA,length(x)-1)
> for(i in 1:(length(x)-1)){
+   Vf[i] = F(x[i+1])-F(x[i])
+ }
> sum(Vf)
[1] 0.9994685

Il faut faire attention lors de la censure, car la discrétisation sera la même pour la charge totale. Une fois discrétisée la loi de la charge individuelle, on travaille sur la loi de la fréquence. Avec les notations de la classe de Sundt, on a
> size=.3
> mu=.12
> a=1-size/(size+mu)
> b=(size-1)*a
> a
[1] 0.2857143
> b
[1] -0.2
L’algorithme de Panjer s’écrit alors
>  G=rep(NA,length(Vf))
>  G[1]=dnbinom(0,size=.3,mu=.12)
>  for(i in 1:(length(G)-1)){
+  G[i+1]=sum((a+b*(1:i)/i)*Vf[1:i]*G[i-(1:i)+1])
+  }
>  sum(G)
[1] 0.9999346
et l’on peut visualiser la fonction de survie avec le code suivant
> lines(x[2:1000],1-cumsum(G),lwd=4)

Fort logiquement, ces méthodes (simulations et Panjer) donnent des grandeurs très proches en terme de quantiles.

mardi 2 mars 2010

Actuariat, calcul(s) de PM (partie 1)

Suite aux deux billets sur la modélisation des sinistres en assurance auto (ici et ), un billet pour décrire le contrat d’assurance-vie sur lequel on travaillera en salle machine mercredi.
Pour s’entraîner, on va considérer un contrat d’assurance sur deux têtes, signé par un couple d’âge 50 ans (pour l’homme) et 52 ans (pour la femme). On supposera leur durées de vies résiduelles modélisées par les tables TV88-90 et TD88-90 respectivement. On supposera dans un premier temps que les durées de vie sont indépendantes (hypothèse que l’on relâchera éventuellement par la suite). Il s’agit d’un contrat d’assurance-vie proposant simplement en garantie en cas de vie (ça n’a rien de plus compliqué si l’on rajoute une garantie décès, on utilise l’additivité des VAP). Le couple cotise (annuellement) pendant 20 ans, soit en payant une prime http://perso.univ-rennes1.fr/arthur.charpentier/latex/tp-sum-21.png si les deux époux sont en vie, soit en payant une prime http://perso.univ-rennes1.fr/arthur.charpentier/latex/tp-sum-22.pngsi l’un des époux décède. En contrepartie, l’assureur verse 10 000 euros si les deux époux sont en vie, et 60% de ce montant au conjoint survivant.
Comme dans le cours, le plus simple pour comprendre ce contrat est sûrement de faire un petit dessin. On met en rouge les flux positifs pour l’assureur (il encaisse de l’argent) et en vert les flux négatifs (il dépense de l’argent). Dans un premier temps, les deux conjoints décèdent pendant leur retraite,

Dans un second temps, un des conjoint décède avant le début de la période de retraite,

  • Calculs des valeurs actuelles probables
La valeur actuelle probable de l’assureur s’écrit comme l’espérance de la somme actualisée des engagements de l’assureur,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/tp-sum-23.png
que l’en sépare entre les engagements versés si une ou deux personnes sont en vie. Si l’on regarde les engagements de l’assuré, alors
http://perso.univ-rennes1.fr/arthur.charpentier/latex/tp-sum-24.png
En notant que http://perso.univ-rennes1.fr/arthur.charpentier/latex/tp-sum-25.png, et http://perso.univ-rennes1.fr/arthur.charpentier/latex/tp-sum-26.png (i.e. les deux membres du couple sont en vie à la signature du contrat).
On calcule la prime annuelle \pi en notant que les deux VAP doivent être égaux,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/tp-sum-27.png
Pour aller plus loin, il faut expliciter les probabilités jointes, en particulier http://perso.univ-rennes1.fr/arthur.charpentier/latex/tp-sum-28.png (les autres découlant de cette dernière). Le début du code est
> Lxv=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/TV8890.csv",header=TRUE,sep=";")
> Lxd=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/TD8890.csv",header=TRUE,sep=";")
On poursuivra en salle machine pour calculer tous ces coefficients, mais aussi calculer les provisions mathématiques..... à suivre donc.
PS et je ferais un jour un billet sur les autres tables utilisées, car je continue à utiliser des tables qui sont dépassées (mais qui, à mon avis, m’enlève rien à la modélisation)

Actuariat, ajustement de lois pour les coûts individuels

Suite de notre billet sur la modélisation des sinistres en assurance auto (ici). Pour commencer, on peut essayer des visualiser les coûts des sinistres, voire le logarithme des coûts de sinistres,
> X=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/base-cout.txt",header=TRUE)$x
> plot(ecdf(X))
> plot(ecdf(log(X)))

On note que l’ajustement sera plus compliqué que l’ajustement de la fréquence. Si l’on regarde les premiers moments, on obtient
> mean(X)
[1] 2881.173
> sd(X)
[1] 39227.22

i.e. des coûs avec une très grande variance.
Classiquement, on peut commencer par faire des ajustements (et des tests d’ajustements) des lois les plus classiques, loi lognormale, loi Gamma, loi exponentielle, loi de Pareto, loi de Weibull, loi log-Gamma, etc....
Commençons par exemple par une loi exponentielle, dont l’estimation par la méthode des moments et par la méthode du maximum de vraisemblance coïncident là encore.
> x<-seq(0,20000,100)
> plot(x,pexp(x,1/mean(X)),type="l",col="red",lwd=3,
+ main="Ajustement d’une loi exponentielle")
> plot(ecdf(X),add=TRUE)
> ks.test(X,"pexp", 1/mean(X))
        One-sample Kolmogorov-Smirnov test
data:  X
D = 0.4274, p-value < 2.2e-16
alternative hypothesis: two-sided 

Bref, l’ajustement n’est pas terrible, mais visuellement, on pouvait s’y attendre. Afin d’améliorer l’ajustement, on peut envisager une loi à deux paramètres. Par exemple, la loi de Weibull, même si n’importe quelle loi ferait l’affaire, a priori.
> fitdistr(X,"weibull")
      shape          scale   
  6.742917e-01   1.513545e+03
 (1.048903e-02) (6.566297e+01)
> ks.test(X,"pweibull", shape=.675,scale=1500)
        One-sample Kolmogorov-Smirnov test
data:  X
D = 0.1912, p-value < 2.2e-16
alternative hypothesis: two-sided
qui est certes, moins pire qu’avant, mais toujours pas génial.... On peut visualiser cet ajustement
> x<-seq(0,20000,100)
> plot(x,pweibull(x,shape=.675,scale=1500),type="l",col="red",lwd=3,
+ main="Ajustement d’une loi de Weibull")
> plot(ecdf(X),add=TRUE)

Comme visiblement l’ajustement par une loi simple ne donne rien, on peut tenter un mélange, permettant en particulier de prendre en compte le caractère multimodal de la distribution, et des nombreux doublons.
Sur 1307 sinistres1, certains montants apparaissent beaucoup
> sum(X==1204)
[1] 216
> sum(X==1128.12)
[1] 151
> sum(X==1172)
[1] 88
> sum(X==1128)

[1] 34
> sum((X>=1128)&(X<=1204))
[1] 507
soit presque 40% des sinistres. On peut donc, un peu abusivement, supposer que l’on a une masse de Dirac à 1170, pour un poids de 40%
> mean(X[(X>=1128)&(X<=1204)])
[1] 1169.583

Pour les coûts de sinistres inférieurs, on peut tenter un ajustement par une loi expontielle tronquée, et au délà une loi de Weibull translatée, ou lognormale, par exemple.
Si l’on cherche le maximum de vraisemblance pour une loi expontielle tronquée, on obtient
> library(stats4)
> ll<-function(theta){
+ mu=1128
+ x=X[X<mu]
+ sum(log(dexp(x,1/theta)/pexp(mu,1/theta)))
+ }
> nlm(ll,p=mean(X[X<1128]))
$minimum
[1] -3476.982
$estimate
[1] 392.5433
$gradient
[1] 0.3038019


Pour la partie supérieure, tentons un ajustement par les deux lois mentionnées
> fitdistr(X[X>1204]-1204,"lognormal")
    meanlog       sdlog  
  7.06796351   1.55143837
 (0.08927525) (0.06312713)
> fitdistr(X[X>1204]-1204,"weibull")
      shape          scale   
     0.5352695   2461.9708666
 (   0.0192046) ( 271.8158455)
Graphiquement, pour les sinistres dépassant 1204 euros, on obtient,

où la courbe bleue est la loi lognormale, et la courbe rouge la loi de Weibull. Bref, on retiendra la loi blue pour 23% des sinistres (les plus gros). Au final, on a ajusté un triple mélange,
  • 40% des sinistres coûtent 1170 euros (type vert)
  • 25% des sinistres suivent une loi lognormale translatée http://perso.univ-rennes1.fr/arthur.charpentier/latex/tp-sum-05.png+1170 (type bleu)
  • 35% des sinistres suivent un loi expontielle tronquée http://perso.univ-rennes1.fr/arthur.charpentier/latex/tp-sum-06.png (type rouge)
Graphiquement l’ajustement est le suivant

avec le code
> plot(x[x<1170],pexp(x[x<1170],1/400)/pexp(1170,1/400)*.35,type="l",col="red",lwd=3,
+ main="Ajustement d’une loi mélange",ylim=c(0,1),xlim=range(x))
> plot(ecdf(X),add=TRUE,col="grey")
> lines(x[x>1170],plnorm(x[x>1170]-1170,7.0679,1.55)*0.25+.75,lwd=3,col="blue")
> lines(x[x<1170],pexp(x[x<1170],1/400)/pexp(1170,1/400)*.35,type="l",col="red",lwd=3)
> segments(1170,.35,1170,.75,lwd=3,col="green")
En creusant un peu, on devrait pouvoir faire mieux, mais en 30 minutes, ce n’est pas trop mal... on peut écrire facilement la "densité" de cette loi, mais aussi un algorithme de simulation.
Notons qu’il aurait également été possible de supprimer la masse de Dirac, et d’ajuster une loi sur les sinistres plus ou moins grands que 1170 euros,
> X=X[(X<1128)|(X>1204)]
> fitdistr(X,"weibull")
      shape          scale   
  5.770185e-01   1.427570e+03
 (1.259671e-02) (9.126271e+01)
> fitdistr(X,"lognormal")
    meanlog       sdlog  
  6.54117332   1.51491685
 (0.05356040) (0.03787292)


Sinon plus généralement, pour estimer par maximum de vraisemblance les paramètres d’un mélange, rien de plus simple, via une routine d’optimisation
Par exemple un mélange de deux lois lognormales, dont une translatée,
> library(stats4)
> ll<-function(theta) {
+ x<-X
+ p =theta[1]
+ m1=theta[2]
+ s1=theta[3]
+ m2=theta[4]
+ s2=theta[5]
+ m =theta[6]
+ -sum(log(
+ p*dlnorm(x,m1,s1)+
+ (1-p)*dlnorm(x-m,m2,s2)*(x>m) )) }
> nlm(ll,p=c(.5,4,.5,6,1.5,2000))
$minimum
[1] 6695.37
$estimate
[1] 9.853219e-01 6.524964e+00 1.519144e+00 6.000751e+00 9.597966e-02
[6] 1.806010e+03

On utilise plutôt http://perso.univ-rennes1.fr/arthur.charpentier/latex/tp-sum-07.png dans la fonction d’optimisation car les algorithmes d’optimisation permettent surtout de minimiser des fonctions. Sinon il existe aussi la library(mixdist), basé sur de l’algorithme EM.
Pour la suite, on verra en salle machine mercredi....
1 Pour reprendre une remarque qui m’a été faite, la base des coûts contenait 1317 sinistres
> length(X)
[1] 1317
ce qui est moins que le nombre effectif de sinistres observés dans la base des fréquences,
> sum(table(N)[2:9])
[1] 3794
i.e. si on exclue les polices n’ayant eu aucun sinistre, il en reste 3794. Ceci ne change rien à la généralité de notre analyse, dès lors que l’on suppose que les 1317 sinistres sont représentatifs.

Actuariat, ajustement de lois pour les fréquences

A la suite des remarques des élèves de M2 lors de la réunion pédagogique, je rajoute un TP dans le cours d’actuariat de M1, qui aura lieu mercredi après midi. Nous ferons un peu ensemble ce que j’ai demandé pour les DM sur d’autres bases, et d’autres contrats. Pour commencer, travaillons sur le DM2, i.e. le modèle collectif. Dans un premier temps, un premier billet pour modéliser la fréquence de sinistres.
Rappelons que nous avons présenté 3 lois de base pour modéliser la fréquence,

  • la loi binomiale, i.e. http://perso.univ-rennes1.fr/arthur.charpentier/latex/tp-sum-01.png
  • la loi de Poisson, i.e. http://perso.univ-rennes1.fr/arthur.charpentier/latex/tp-sum-02.png
  • la loi binomiale négative, i.e. http://perso.univ-rennes1.fr/arthur.charpentier/latex/tp-sum-03.png
Afin d’opérer rapidement une première sélection, on peut calculer espérance et variance empiriques,
> N=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/base-nb.txt",header=TRUE)$x
> mean(N)
[1] 0.1216735
> var(N)
[1] 0.1831262

Avant de se lancer dans toute modélisation, regardons également la distribution du nombre de sinistres par contrat
> (freq.empirique = table(N))
N
    0     1     2     3     4     5     7    16    21
35549  3008   640   119    22     2     1     1     1

En l’occurence, la loi binomiale négative semblerait plus adaptée. Mais auparavant, regardons les autres lois.Dans le cas de la loi de Poisson, l’estimateur du maximum de vraisemblance et l’estimateur par la méthode des moments coïncide. Le paramètre est ici
> lambda=mean(N)
Si l’on compare la distribution empirique avec la loi de Poisson, on obtient
> (freq.theorique = length(N)*dpois(as.numeric(names(freq.empirique)),lambda))
[1] 3.483576e+04 4.238589e+03 2.578619e+02 1.045832e+01 3.181251e-01
[6] 7.741478e-03 2.728767e-06 3.841856e-24 4.195623e-35


On peut aussi utiliser différents tests pour vérifier la qualité de l’ajustement
> library(vcd)  
Le chargement a nécessité le package : MASS
Le chargement a nécessité le package : grid
Le chargement a nécessité le package : colorspace
> gof = goodfit(N,type= "poisson",method= "ML")
> plot(gof,main="Ajustement d’une loi de Poisson")

Le test évoqué ici est celui du chi-deux, qui formellement donne les résultats suivants (en bricolant à la main pour obtenir une distribution empirique),
> x=seq(0,25)
> y=table(N)/length(N)
> freq.empirique=rep(0,26)
> freq.empirique[1:6]=y[1:6]
> freq.empirique[8]=y[7]
> freq.empirique[16+1]=y[8]
> freq.empirique[22]=y[9]
>  (freq.theorique1 = dpois(x,lambda))
 [1] 8.854374e-01 1.077343e-01 6.554202e-03 2.658242e-04 8.085939e-06
 [6] 1.967689e-07 3.990259e-09 6.935839e-11 1.054885e-12 1.426128e-14
[11] 1.735219e-16 1.919365e-18 1.946132e-20 1.821482e-22 1.583044e-24
[16] 1.284096e-26 9.765029e-29 6.989089e-31 4.724371e-33 3.025425e-35
[21] 1.840570e-37 1.066422e-39 5.897966e-42 3.120114e-44 1.581813e-46
[26] 7.698588e-49 
> chisq.test(x=freq.empirique,p=freq.theorique1)
        Chi-squared test for given probabilities
data:  freq.empirique
X-squared = 6.058094e+29, df = 25, p-value < 2.2e-16

Bref, on a du mal à accepter l’ajustement de la loi de Poisson (mais on s’y attendait un peu).
Pour la loi binomiale négative, les méthodes classiques (méthode des moments et méthode du maximum de vraisemblance) n’ont pas de solution analytique. Il faut donc faire des calculs. Pour le maximum dde vraisemblance, on peut utiliser le code suivant,
> library(MASS)
> fitdistr(N,"negative binomial")
      size           mu    
  0.287203791   0.121677345
 (0.013202351) (0.002098181)

Mais on peut aussi tenter un test de goodness of fiit,
> gof = goodfit(N,type= "nbinomial",method= "ML")
> plot(gof,main="Ajustement d’une loi Binomiale Négative")

On peut là aussi faire un test du chi-deux,
> x=seq(0,100)
> y=table(N)/length(N)
> freq.empirique=rep(0,101)
> freq.empirique[1:6]=y[1:6]
> freq.empirique[8]=y[7]
> freq.empirique[16+1]=y[8]
> freq.empirique[22]=y[9]
> parametres= fitdistr(N,"negative binomial")$estimate
>  (freq.theorique2 = dnbinom(x,parametres[1],mu=parametres[2]))
  [1] 9.035266e-01 7.722249e-02 1.479019e-02 3.355599e-03 8.206336e-04
  [6] 2.093949e-04 5.491026e-05 1.467661e-05 3.978407e-06 1.090153e-06
[...]
 [96] 1.104459e-52 3.262314e-53 9.636853e-54 2.846935e-54 8.411084e-55
[101] 2.485180e-55
> chisq.test(x=freq.empirique,p=freq.theorique2)
        Chi-squared test for given probabilities
data:  freq.empirique
X-squared = 2237.183, df = 100, p-value < 2.2e-16

autrement dit, on accepte là aussi difficilement l’hypothèse de loi binomiale négative, même si la distance du chi-deux est beaucoup plus faible qu’avec la loi de Poisson. Bref, comme il faut retenir un modèle, on acceptera le modèle suivant
http://perso.univ-rennes1.fr/arthur.charpentier/latex/tp-sum-04.png
dont la  moyenne théorique vaut 0,121, et la variance théorique 0,17, ce qui n’est pas trop éloigné de
> mean(N)
[1] 0.1216735
> var(N)
[1] 0.1831262

- page 2 de 9 -