Freakonometrics

To content | To menu | To search

Teaching › actuariat - M1-09/10

Entries feed - Comments feed

Thursday, March 18 2010

Toujours des calculs de PM

A la demande générale, un court billet sur les calculs de PM (oui, encore un après celui sur la Zillmérisation - ici - mais surtout le calcul détaillé dans le cas d'une temporaire décès - ici et ). Mais après l'assurance en cas de décès, on va parler aujourd'hui de l'assurance en cas de vie. Une rente en quelque sorte, ou une retraite. Avant tout, suite à un commentaire (posté ici), je vais juste remettre un lien vers les notations actuarielles de base (ici), emprunté à l'ouvrage de Pierre Petauton. J'ai rajouté des tables de la Society of Actuaires, et quelques articles de presse sur l'assurance vie, et la modélisation.

  • Description du contrat, et calcul de la prime annuelle
On considère un assuré d'âge x, cotisant pendant m années pour sa retraite, et touchant au bout de n années de cotisation une rente annuelle d'un montant C, payé tous les ans à terme échu s'il est en vie, jusqu'à son décès (i.e. une annuité viagère). La prime pure unique (correspondant à la valeur actuelle probable des engagements de l'assureur) est de la forme
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM-rente-01.png
qui peut se réécrire, en utilisant les notations usuelles,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM-rente-02.png
(i.e. la valeur actuelle probable d'une annuité viagère différée de n années). Si l'assuré paye une prime annuelle constante pendant ces n années, en début d'année, alors la prime est
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM-rente-03.png
On peut alors passer au calcul de la provision mathématique, en notant qu'il faudra distinguer les n premières années (période où l'assuré paye sa prime) et les dernières (période où l'assureur verse la rente).
  • Calcul de la PM, méthode prospective
Comme toujours pour le calcul des PM, on se place au bout de k années. Si k<n (l'assuré paye encore sa prime), en faisant la différence entre les engagements restants de l'assureur et ceux de l'assuré, on obtient
L'image “http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM-rente-04.png” ne peut être affichée car elle contient des erreurs.
Si en revanche on suppose que k> n (seul l'assureur a encore des engagements) alors
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM-rente-05.png
Tout simplement. En effet, dans le premier cas, l'assuré à vieilli, et il a moins de versements à venir (c'est la partie de droite). Pour l'assureur, il s'agit toujours d'une annuité différée (de moins en moins différée). Dans le second cas, l'assureur doit verser une rente viagère.
  • Calcul de la PM, méthode rétrospective
Là aussi, il faut distinguer suivant la valeur de k. Si k<n, on obtient simplement que
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM-rente-06.png
puisque sur cette période, seule l'assuré a pris des engagements. Pour rappels (je repense au commentaire posté ici), http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM-rente-07.png est la vap du capital diffféré, i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM-rente-08.png
Pour la seconde période, si k> n,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM-rente-09.png
avec à gauche un terme constant, les engagements de l'assuré entant passés, et à droite les engagements qu'avait pris l'assureur, i.e. les k-n années qui ont suivi l'année n.
  • Calcul de la PM, méthode par récurence
Allez, pour prouver que je suis capable de faire des calculs de PM, on va tenter le dernier, à savoir la méthode par récurence... J'ai l'impression que si k<n, on aurait quelque chose de la forme
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM-rente-10.png
alors que si k> n
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM-rente-15.png
  • Calcul numérique de la PM
Bon, c'est bien joli d'écrire des formules incompréhensibles (comme dirait le Canard Enchaîné, ici), maintenant il va falloir taper le code, si possible sans erreurs. Pour reprendre une partie du code que j'avais fait la dernière fois (ici),
>  L=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/TD8890.csv",
+  sep=";",header=TRUE)
>  Lx=L$Lx
>  m=length(Lx)
>  p=matrix(0,m,m); d=p
>  for(i in 1:m){
+  p[1:(m-i+1),i]=Lx[1+i:m]/Lx[i]
+  d[1:(m-i+1),i]=(Lx[i:(m)]-Lx[i:(m)+1])/Lx[i]}
>  diag(d[m:1,])=0
>  diag(p[m:1,])=0
>  q=1-p
>  a=E=A=adiff=matrix(0,m,m)
>  r=3/100
>  for(i in 1:(m-1)){
+  a[,i]=cumsum(1/(1+r)^(0:(m-1))*c(1,p[1:(m-1),i]))
+  }
>  for(i in 1:(m-1)){
+  #for(j in 1:(m-i)){
+  A[,i]=cumsum(1/(1+r)^(1:m)*d[,i])
+  }
>  for(i in 1:m){
+  E[,i]=(1/(1+r)^(1:m))*p[,i]
+  }

Par rapport à cet ancien code, je vais juste définir la fonction suivante pour calculer http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM-rente-16.png, en utilisant le fait que
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM-rente-17.png
>  for(i in 1:(m-1)){
+  adiff[(1+0:(m-i-1)),i]=E[(1+0:(m-i-1)),i]*a[m,1+i+(0:(m-i-1))]
+  }

On en déduit alors la prime annuelle, en utilisant le fait que
>  a[n,x]
[1] 19.16241
>  sum(1/(1+r)^(0:(n-1))*c(1,p[1:(n-1),x]) )
[1] 19.16241
>  adiff[n,x]
[1] 4.141035
>  sum(1/(1+r)^((n):m)*p[(n):m,x] )
[1] 4.141035
>  (prime=adiff[n,x] / (a[n,x]))
[1] 0.2161019
>  sum(1/(1+r)^((n):m)*p[(n):m,x] )/sum(1/(1+r)^(0:(n-1))*c(1,p[1:(n-1),x]) )
[1] 0.2161019
(les calculs de part et d'autres sont là simplement pour vérifier). Pour le calcul de la PM, je n'ai programmé que les deux premières méthodes,
>  VR=VP=rep(NA,m-x)
>  ### méthode rétrospective
>  adiff[n-1,x+1]-prime*a[n-1,x+1]
[1] 0.2230330
>  adiff[n-2,x+2]-prime*a[n-2,x+2]
[1] 0.453264
>  (VR[1:(n-1)]=diag(adiff[n-(1:(n-1)),x+(1:(n-1))] - a[n-(1:(n-1)),x+(1:(n-1))]*prime))
 [1]  0.2230330  0.4532640  0.6909852  0.9365169  1.1901967  1.4523249  1.7233908  2.0039135  2.2942649
[10]  2.5954383  2.9077644  3.2316703  3.5678273  3.9170547  4.2805038  4.6597691  5.0559545  5.4702573
[19]  5.9044436  6.3600327  6.8391596  7.3447784  7.8776815  8.4400483  9.0368831  9.6701276 10.3447773
[28] 11.0642559 11.8332036
>  (VR[n:(m-x)]=a[m,x+n:(m-x)])
 [1] 12.656739 12.259522 11.858376 11.451953 11.045177 10.637713 10.230173  9.822171  9.421634  9.023840
[11]  8.630055  8.238339  7.854197  7.474624  7.105806  6.748102  6.400061  6.061537  5.743411  5.437153
[21]  5.143796  4.864296  4.596681  4.338471  4.097980  3.879335  3.672003  3.474060  3.288933  3.123415
[31]  2.962928  2.804341  2.631413  2.464020  2.272099  2.140384  2.023163  1.911479  1.791177  1.673873
[41]  1.510665  1.277393  1.000000  0.000000
>  plot(x+0:(m-x),c(0,VR),col="red",xlab="âge de l'assuré",ylab="montant de la PM")
Pour la méthode prospective, on utilisera le fait que
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM-rente-18.png
>  ### méthode prospective
>  points(x+0:(m-x),c(0,VP),pch=4,col="blue")
>  VP[1:(n)]=a[1:(n),x]*prime/E[1:(n),x]
>  VP[(n+1):(m-x)]=(a[n,x]*prime- (adiff[(n+1),x]-adiff[(n+1):(m-x)+1,x]))/E[(n+1):(m-x),x]
>  points(x+0:(m-x),c(0,VP),pch=4,col="blue")
Avec le graphique suivant pour les deux méthodes, qui donnent ici la même chose,


On notera que si on change l'âge de l'assuré, qui aura maintenant 45 ans, mais qui part toujours à la retraite à 65 ans, la PM change entre 45 et 65 (il n'y a rien entre 35 et 45 puisque personnne n'a d'engagement, et surtout ne change pas une fois atteint l'âge de la retaite, les engagements de l'assureur étant identiques pour ces deux contrats.

Pour ceux qui n'auraient pas noté que la plupart des contrats sont des combinaisons linéaires de contrats plus simples (que l'on peut toujours décomposer sous la forme de contrats simples d'assurance décès, ou de rentes plus ou moins différées), et que l'espérance mathématique (et donc les valeurs actuelles probables) est un opérateur linéaire, tous les calculs de PM peuvent se faire avec les morceaux de code que j'ai mis en ligne. Mais je n'en dis pas davantage, sinon je risque de me faire incendier par les consultants qui verront d'un mauvais oeil que je mette en ligne du code gratuitement (et tapé en 2 heures) alors qu'ils facturent ça généralement une fortune...

Wednesday, March 10 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 () !

Wednesday, March 3 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:\ -data\\base-nb.txt",header=TRUE)$x
> Xemp=read.table("D:\ -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.

Tuesday, March 2 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

Thursday, February 18 2010

Sommes et assurance décès

Un rapide billet suite à une question qui devrait intéresser les élèves qui travaillent sur leur projet d'assurance-vie: "pourquoi a-ton une somme dans la VAP d'un contrat d'assurance-décès ?" (ou pour compléter "parce dans des rentes, je comprends bien, ce sont les somme des différents flux").
Effectivement, je pense que les deux sommes sont assez différentes. Dans le cas de l'assurance décès (en fait j'avais évoqué tout ça dans un précédant billet, ici pour la théorie et pour la mise en œuvre pratique), supposons que l'on verse 1€ lors du décès d'une personne d'âge http://perso.univ-rennes1.fr/arthur.charpentier/latex/somme-dc-0.png à la souscription d'un contrat, le flux actualisé est la variable aléatoire

http://perso.univ-rennes1.fr/arthur.charpentier/latex/somme-dc-1.png
où http://perso.univ-rennes1.fr/arthur.charpentier/latex/somme-dc-2.png est la durée de vie résiduelle de l'individu. La valeur actuelle probable est l'espérance mathématique de ce flux, i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/somme-dc-3.png
Pour effectuer ce calcul, on utilise la formule des probabilités totales, en prenant comme événement conditionnant la date du décès, i.e.
L'image “http://perso.univ-rennes1.fr/arthur.charpentier/latex/somme-dc-4.png” ne peut être affichée car elle contient des erreurs.
 Bref, on peut écrire
http://perso.univ-rennes1.fr/arthur.charpentier/latex/somme-dc-5b.png
ou quelque chose du genre, d'où la somme. Notons que si les décès surviennent en milieu d'année, comme l'évoque Cimon (ici), on aurait
http://perso.univ-rennes1.fr/arthur.charpentier/latex/somme-dc-6.png
si l'on suppose que les décès surviennent

Monday, February 1 2010

Calculer une provision mathématique (2)

Après un petit clin d'œil à Michel Rabagliati qui a gagné le prix du public hier à Angoulême.... Je vais poursuivre le billet commencé la semaine dernière (ici), histoire de montrer que ce n'est pas si compliqué que ça de calculer des PM avec R. Compte tenu des doubles indices dans toutes les fonctions, nous allons créer des matrices afin de stocker les différentes grandeurs.
> L=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/TD8890.csv",
+ sep=";",header=TRUE)
> Lx=L$Lx
> m=length(Lx)
> p=matrix(0,m,m); d=p
> for(i in 1:m){
+ p[1:(m-i+1),i]=Lx[1+i:m]/Lx[i]
+ d[1:(m-i+1),i]=(Lx[i:(m)]-Lx[i:(m)+1])/Lx[i]}
> diag(d[m:1,])=0
> diag(p[m:1,])=0
> q=1-p
aussi, p[10,40] correspondra à http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM80.png (attention, cf commentaire plus bas). Je stocke ainsi la probabilité de décéder au bout de k années, ou bien exactement dansk années. On peut alors calculer les différentes grandeurs,
> a=E=A=matrix(0,m,m)
> r=3/100
> for(i in 1:(m-1)){
+ a[,i]=cumsum(1/(1+r)^(0:(m-1))*c(1,p[1:(m-1),i]))
+ }
> for(i in 1:(m-1)){
+ #for(j in 1:(m-i)){
+ #A[j,i]=sum(1/(1+r)^(1:j)*p[1:j,i]*q[1,i-1+1:j])
+ A[,i]=cumsum(1/(1+r)^(1:m)*d[,i])
+ }
> for(i in 1:m){
+ E[,i]=(1/(1+r)^(1:m))*p[,i]
+ }

A partir de ces constantes (ou de ces matrices de constantes), on peut déjà calculer la prime annuelle du contrat décès que nous avions considéré,
> x=50; n=30
> prime=A[n,x]/a[n,x]
> sum(prime/(1+r)^(0:(n-1))*c(1,p[1:(n-1),x]))
[1] 0.3116454
> sum(1/(1+r)^(1:n)*d[1:n,x])
[1] 0.3116454

De plus, on peut calculer les provisions mathématiques à l'aide des différentes méthodes (et vérifier qu'elles donnent bien le même résultat).
> ### méthode rétrospective
> VR=(prime*a[1:n,x]-A[1:n,x])/E[1:n,x]
> plot(0:n,c(0,VR),col="red")
> ### méthode prospective
> VP=diag(A[n-(0:(n-1)),x+(0:(n-1))])-prime*diag(a[n-(0:(n-1)),x+(0:(n-1))])
> points(0:n,c(VP,0),pch=4,col="blue")
> # méthode itérative
> VI=0
> for(k in 1:n){
+ VI=c(VI,(VI[k]+prime-A[1,x+k-1])/E[1,x+k-1])
+ }
> points(0:n,VI,pch=5,col="green")
Pour un contrat signé par une personne de 50 ans, pour une couverture de 30 ans, au taux d'actualisation de 3%, on obtient les provisions mathématiques suivantes,

(les trois vecteurs coïncident, presque par miracle). Si la personne avait 45 ans, les provisions auraient l'allure suivante (en bleu),

Enfin, si on regarde la sensibilité au taux d'actualisation, un taux de 5% donnerait les provisions suivantes (en vert),

Ben voilà. La prochaine étape ça sera de distribuer les devoirs maison....

Friday, January 29 2010

Calculer une provision mathématique (1)

Bon, je vais essayer de clarifier ici une partie du cours que, tous les ans, je rate: les méthodes de calculs de la provision mathématique. Après quelques éléments de droit, je pense que le plus simple est d'illustrer le calcul de la PM sur un cas simple, comme en cours, i.e. une temporaire décès avec prime annuelle. Cet exemple avait été repris ici, très (trop) brièvement. Sinon je peux renvoyer vers les slides de Pierre Devolder, ici,

  • Valorisation d'une assurance temporaire décès
Le "principe fondamental de valorisation", on doit avoir
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM01.png
en faisant une valorisation à la date 0, i.e. la date de souscription du contrat.
Pour l'assuré, il souhaite payer une prime annuelle constante http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM02.png, noté plus simplement http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM03.png, tant qu'il est en vie i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM04.png

http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM05.png
(le http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM06.png car le paiement se faisant ici en début de période).  De même,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM07.png

http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM08.png
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.losser-france_s.jpg (l'indemnité étant versé par l'assureur à terme échu). J'utilie ici les notations officielles selons les normes françaises que l'on retrouve dans les ouvrages de Pierre Petauton et Christian Hess (ormis peut être pour le taux d'actualisation qui est généralement noté i). Il existe des nuances entre les notations françaises et les notations internationales. Je renvoie ici ou pour les standards internationaux. Parmi les notations internationales, on notera que ce que j'ai noté A ici est souvent noté \bar{A} dans la littérateure anglosaxonne (le A signifiant un paiement en milieu d'année).
Berf, pour revenir à notre prime annuelle, on en déduit que
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM09.png
  • Définition de la provision mathématique
Notons http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM10.png la valeur actuelle probable, en http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM11.png, des engagements de l'assurés pour la période http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM12.png. Aussi, http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM13.png sera la valeur actuelle probable, en 0, des k premières primes annuelles. Et on notera http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM14.png la valeur actuelle probable, en 0, des engagements de l'assurés pour la période http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM15.png, i.e. la valeur actuelle probable des n-k dernières primes annuelles.
De manière analogue, notons http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM16.png la valeur actuelle probable, en http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM11.png, des engagements de l'assureur pour la période http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM12.png.
Comme tenu du principe fondamental de valorisation, pour un contrat arrivant à échéance au bout de n années, on doit avoir
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM20.png
pour un contrat soucrit à la date 0 et tel qu'il n'y a plus d'engagement de part et d'autre passé n années.
Aussi
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM21.png
avec, de manière générale
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM22.png
et
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM23.png
(d'où le principe d'inversion du cycle de production de l'assurance).
La provision mathématique (pures) de l'année k sera noté http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM30.png si elles sont actualisée à la date t. La référence étant L'image “http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM31.png” ne peut être affichée car elle contient des erreurs. (i.e. on actualise en k). C'est ce qui a été noté http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM32.png dans le cours. On  définit http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM33.png par
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM34.png
Cette définition sera dite rétrospective (car on se place sur la période antérieure à k).On peut aussi écrire, de manière équivalente (compte tenu du principe de valorisation)
http://perso.univ-rennes1.fr/arthur.charpentier/latex/MP000.png
Cette définition sera dite prospective (car on se place sur la période postérieure à k).
Enfin, il existe une dernière méthode, correspondant à une simple mise à jour, i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM36.png
Cette méthode sera dite itérative,voire en l'occurence itérative ascendante, car on initialise avec http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM37.png. Mais il sera aussi possible de construire une méthode itérative descendante, récursive, commençant à la fin du contrat. Mais n'essayons pas de compliquer les choses (j'ai déjà suffisement de mal à m'y retrouver moi même).
  • Un peu de législation
Comme l'explique l'article R331-3 du Code des Assurance (ici), "Provision mathématique : différence entre les valeurs actuelles des engagements respectivement pris par l'assureur et par les assurés, à l'exception, pour les contrats mentionnés à l'article L. 142-1, des engagements relatifs à la provision de diversification" ou ici pour le lexique de l'ACAM. Le document ici reprend aussi une définition de ces provisions,
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/2008_11_30_Le_Monde_Gerner-1.jpg
  • La vision du Bowers et al.
Dans l'actuariat américain, le calcul des benefit reserves se fait en introduisant http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM40.png, correspondant à la valeur actualisée à la date t des gains (ou des pertes) futurs de l'assureur (obtenus comme différence entre les engagements de part et d'autre). La provision mathématique est alors "l'espérance de la valeur actualisée à la date k des gains futurs de l'assureur conditionnelemnt au fait que (x) soit en vie". La forme concise est alors
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM41.png
C'est cette définition que nous avions utilisé en cours.
  • La méthode prospective
Nous avions commencé, en cours, par la méthode prospective, car c'est la plus naturelle (de mon point de vue) compte tenu de la législation. Nous avions posé
http://perso.univ-rennes1.fr/arthur.charpentier/latex/MP000.png
Notons que
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM51.png
où le terme de droite désigne la valeur actuelle probable d'un capital différé, relatif au versement d'un euro dans n année, conditionnée par la survie de (x), i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM52.png
Si l'on se place à la date k (car c'est le plus simple, mais l'assuré a alors l'âge x+k), notons que la différence entre les valeurs actuelles probables des engagements des deux parties donne, simplement
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM53.png
car d'un côté, on a une temporaire décès sur les n-k années restantes pour un assuré d'âge x+k, et de l'autre, l'assuré a pris l'engagement de verser sa prime (qui reste inchangée) pendant n-k années s'il vit. Aussi,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM54.png
où l'on considère des assurances décès différées. On peut aussi écrire
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM55.png
  • La méthode rétrospective
La méthode rétrospective n'est pas beaucoup plus compliquée. On écrit simplement
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM60.png
i.e. http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM61.png. Or http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM62.png, et donc
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM63.png
  • La méthode itérative
L'idée est ici de décrire la variation de la provision mathématique entre deux dates en fonction des variation des engagements de part et d'autre. D'un côté il y a le paiement de la prime (en début de période, donc pas de problème d'actualisation et de non paiement), et de l'autre, une assurance décès sur un an. Aussi
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM70.png
Or http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM72.png, ce qui donne, finalement
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM73.png
avec la convention que la première PM est nulle (de part notre principe fondamental de valorisation).
  • Mise en oeuvre pratique
Bon, c'est bien joli ces formules, mais il serait bon de voir si on peut faire des calculs. J'avais fait en cours les calculs dans une feuille excel (comme ici), mais avec R on devrait pouvoir y arriver facilement.... ça sera l'objet d'un billet qui arrivera rapidement (normalement).

Thursday, January 21 2010

Modéliser des durées de vies humaines (1)

Dans le TD de ce matin, nous avons commencé par le modèle d'Abraham De Moivre, avant de présenter celui de Benjamin Gompertz. Je vais continuer à illustrer historiquement le cours d'assurance-vie, pour revenir sur les grands modèles paramétriques classiques...

  • Le modèle de De Moivre
Le modèle d'Abraham De Moivre peut paraître trop simple, mais il n'en est pas moins relativement riche, et intéressant. L'idée de De Moivre était que la durée de vie résiduelle pour un individu d'âge x était uniformément distribuée entre 0 et w-x, w étant la durée de vie limite d'un individu. Le taux de hasard est alors une fonction croissante de l'âge, ce qui est une hypothèse relativement générale (en particulier passé 10 ans). Rappelons qu'en 1724, c'était un blasphème de penser qu'il existait des lois (autres que divines) pour modéliser la vie humaine.
De Moivre a pu acquérir une légitimité sur la modélisation des durées de vies humaines après avoir quelque temps en Angleterre, à côtoyer en particulier Edmond Halley. L'autre légitimité a été acquise sur son lit de mort, en 1754, à Londres. La légende raconte que De Moivre avait remarqué qu'il dormait, chaque nuit, 15 minutes de plus. En faisant un petit calcul sur les suites arithmétique, il avait prédit qu'il mourrait le jour où il dormirait 24 heures. Et la légende prétend qu'il ne s'est pas trompé (ici).
  • Le modèle de Gompertz
Un siècle plus tard, Benjamin Gompertz proposait un autre modèle pour modéliser la durée de vie humaine. Si on retient souvent la forme simple sur le taux de hasard, notons qu'il est possible de caractériser cette loi en notant que
Il est alors possible d'ajuster cette loi de manière assez simple, en demandant à ce qu'elle passe par 3 points particuliers (comme il y a 3 paramètres). En cours, nous avions mis des contraintes sur les âges 50, 60 et 70. Le code est le suivant (pour une ajustement de la table TD88-90),
> L=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/TD8890.csv",
+ sep=";",header=TRUE)
> h=10;  x1=50
> x2=x1+h;  x3=x2+h
> c10=log(L$Lx[L$Age==x2]/L$Lx[L$Age==x3])/log(L$Lx[L$Age==x1]/L$Lx[L$Age==x2])
> c=c10^(1/(x2-x1))
> g=exp(log(L$Lx[L$Age==x1]/L$Lx[L$Age==x2])/c^x1/(1-c^(x2-x1)))
> k=L$Lx[L$Age==x1]/(g^(c^x1))
> LG=k*g^(c^(L$Age))
> plot(L$Age,LG,type='l',col="blue",xlab="Age",ylab="Nombre de survivants Lx")
> lines(L$Age,L$Lx,col="red")
L'ajustement se visualise sur la figure ci-dessous,

avec l'erreur observée ci-dessous.

Mais on peut tenter un ajustement sur les âges 30, 50 et 70,

Saturday, January 16 2010

Modéliser les nombres de sinistres

La surdispersion est souvent évoquée pour les lois binomiales ou de Poisson. J'avais évoqué le fait qu'empiriquement, la variance d'une variable de nombre est souvent plus élevée que l'espérance (alors que pour une loi de Poisson, les deux sont supposés égaux, j'en parlais en particulier ici). La sur-dispersion signifie que la variance de Y dépasse la variance attendue. Considérons un modèle binomial par exemple, i.e. Y prend deux valeurs, 0 ou  1.

On sait que
et que
Mais il est possible qu'en pratique, la variance soit différente (souvent supérieure) de cette valeur. Assez souvent, il existe des classes (clusters) non identifiées par les variables à notre disposition, . Comme nous l'avions mentionné ici (par exemple), l'hétérogénéité non-observée entraine de la surdispersion.
Il existe des tests de surdispersion, détaillé dans le cas Poissonnien dans la section 3.4 du livre de Colin Cameron et Pravin Trivedi, Regression Analysis of Count Data. L'idée la plus simple est de noter que la loi de Poisson est un cas particulier de la loi binomiale négative, dans le cas où http://perso.univ-rennes1.fr/arthur.charpentier/latex/offset0.png. Il suffit de faire un simple test (ratio de vraisemblance par exemple) de nullité d'un coefficient. Le hic est qu'on teste si un coefficient est sur le bord des valeurs possibles, mais Moran (1971) par exemple propose des solutions. Sur l'estimation du paramètre de surdispersion, j'en ai parlé ici (en présentant la loi quasiPoisson).
Sinon parmi les notions classiques sur la modélisation des nombres de sinistres, il y a la variable offset. Pour expliquer cette notion, on suppose que l'on rajoute un terme au prédicteur linéaire, i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/offset01.png
de telle sorte que
http://perso.univ-rennes1.fr/arthur.charpentier/latex/offset02.png
L'exemple classique est obtenu dans le cas où le nombre suit une loi de Poisson et que l'on prend la fonction de lien canonique, i.e. la fonction log,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/offset03.png
Aussi, si l'on considère http://perso.univ-rennes1.fr/arthur.charpentier/latex/offset04.png on obtient un effet multiplicatif. Et http://perso.univ-rennes1.fr/arthur.charpentier/latex/offset05.png est alors l'exposition, i.e. la durée de la période de couverture (une personne restant 6 mois au lieu d'un an devrait avoir deux fois moins de sinistres, en moyenne, qu'une personne ayant des caractéristiques identiques. Au niveau des unités, notre référence sera une exposition de 1. Il convient alors de normaliser l'exposition afin qu'une année corresponde à 1 (pour des contrats d'une période d'un an). On parlera parfois du "nombre d'années police" Sous R, on peut utiliser le code suivant
>  glm(Y~X1+X2,family=poisson(link=log),offset=log(E))
Dans le cas d'une variable binomiale (y-a-t-il eu des sinistres, ou pas), par exemple avec un modèle logistique, on peut introduire la variable offset via la commande suivante sous R,
>  glm(Y~X1+X2+offset(E),family=binomial)

Thursday, January 7 2010

Réassurance optimale, une introduction

Je vais revenir un peu sur l'exercice 3 du TD de ce matin (ici). Le but était de voir comment lier la probabilité de ruine, le montant de fonds propres, et le niveau des primes. Nous relierons tout ça de manière formelle par la suite (dans la partie du cours abordant plus spécifiquement les processus de comptage et la probabilité de ruine dans un modèle à la Cramér-Lundberg), mais le modèle simple présenté ici donne une première idée. En particulier, nous allons évoquer l'utilisation de la réassurance non-proportionnelle pour réduire la probabilité de ruine.
Dans le cas où un traité de réassurance proportionelle (avec rétention ) pour la classe 2 la  probabilité de ruine s'écrit

avec et donc, toujours en supposant  et  indépendants,
Le bénéfice moyen global est alors
En élevant au carré, trouver  en fonction de T revient à résoudre une équation du second degré
On considère les paramètres suivants,
> m1=800
> m2=900
> s1=200
> k1=.04
> s2=300
> k2=.06
> R=800
> T=4

T=4 si et seulement si  et  (une seule des racines de l'équation està valeur dans [0,1]).
> f2=function(x2){
+ T^2*(s1^2+x2^2*s2^2)-(R+k1*m1+x2*k2*m2)^2  }
> uniroot(f2 , c(0,1))
$root
[1] 0.2244363

On peut visualiser ça sur la figure ci-dessous,
> X=seq(-.2,1.2,by=0.0025)
> Y2=f2(X)
> plot(X,Y2,type="l",ylim=range(c(0,Y1)),col="red",lwd=2)

On peut faire pareil pour un programme de réassurance uniquement sur la première branche, mais numériquement il n'est pas possible même en cédant tout le risque) d'atteindre la valeur souhaitée,
> f1=function(x1){
+ T^2*(x1^2*s1^2+s2^2)-(R+x1*k1*m1+k2*m2)^2  }
> uniroot(f1 , c(0,1))
Erreur dans uniroot(f1, c(0, 1)) :
  les valeurs de f() aux points extrêmes ne sont pas de signe opposé

Dans le cas où un traité de réassurance proportionnelle (avec rétention ) existe pour les deux classes, la probabilité de ruine s'écrit
comme précédemment, saut que cette fois  et donc,
qui peut s'inverser afin d'exprimer en fonction de T,
Le bénéfice moyen global est alors

T=4 si et seulement si  et ,
> f=function(x){
+ T^2*(x^2*s1^2+x^2*s2^2)-(R+x*k1*m1+x*k2*m2)^2  }
> uniroot(f , c(0,1))
$root
[1] 0.5898778
Dans le cas où deux traités de réassurance proportionelle (avec rétentions  et ) existent pour les deux classes de risques, et
Le bénéfice moyen global est alors . On a un degré de liberté supplémentaire, que l'on peut utiliser pour maximiser le profit espéré de l'assureur. Les coefficients   et  optimaux sont les solutions du programme d'optimisation
qui peut se réécrire%
Le Lagrangien du programme précédant s'écrit

i.e.
L'optimum  est alors la solution du système 3 équations et 3 inconnues , correspondant aux trois conditions du 1er ordre,

(je ferai un billet un jour sur l'interprétation du multiplicateur de Lagrange...).Des deux premi\`{e}res équations, puisque
on peut écrire,
soit
et donc
ou encore
En substituant dans la dernière équation, on en déduit alors
soit
ou encore

De façon similaire, s'obtient à l'aide d'une équation de degré deux analogue à celle ci-dessus.
La résolution de l'équation ci-dessus donne  et . Le bénéfice moyen est alors de 51,25.
Numériquement, je j'ai pas calculé la valeur optimale, mais l'ensemble des traités possible de réassurance sont dans la partie inférieure du graphique ci-dessous. Il serait possible de calculer le bénéfice le long de la frontière bleue (correspondant au cas où la contrainte est saturée, i.e. la probabilité de ruine est exactement celle qui était requise)