Freakonometrics

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
où  désigne la prime commerciale unique,  désignent les frais d'acquisition,  correspond aux frais de gestion des engagements de l'assureur (il n'y a pas de  car la prime est versée au début, en une seule fois), 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  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),
soit encore
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 , souscrite à l'âge . Comme auparavant, on note  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  années, i.e.
qui peut s'écrire
soit
Si l'on considère une prime annuelle, payée tout au long des n années, on a quelque chose de la forme
où pour rappel (je renvoie au cours à ce sujet)  désigne une assurance temporaire décès dont le capital croît d'une unité chaque année, ce qui donne
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  (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 € à l'assureur, la prime annuelle payée par l'assureur doit être multipliée par (de telle sorte que les VAP initiale soient égales). Dans ce cas, un rapide calcul montre que la PM doit être multipliée par  é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
si la durée de paiement des primes était initialement , et que  est la durée de la garantie. Si =, il n'y a alors pas de provision de gestion, le second terme étant nul. Cette provision existe dès lors que  (ce qui arrive au final assez souvent me semble-t-il). On arrive ainsi à la provision d'inventaire, qui est souvent notée . 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,
On pose alors
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 (là) !

Suite aux deux billets sur la modélisation des sinistres en assurance auto (ici et là), 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  si les deux époux sont en vie, soit en payant une prime si 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,
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
En notant que , et  (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,
Pour aller plus loin, il faut expliciter les probabilités jointes, en particulier  (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)

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. 
• la loi de Poisson, i.e. 
• la loi binomiale négative, i.e. 

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
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

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 à  (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....

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,

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
où
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
L'optimum  est alors la solution du système 3 équations et 3 inconnues , correspondant aux trois conditions du 1er ordre,
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)