Category Archives: Actuarial science

Modélisation des tables de mortalité

En cours d’actuariat, nous commencerons à parler des tables de mortalité prospectives.

Côté théorie, je peux renvoyer vers quelques slides de Pierre Devolder, ici. Je renvoie ici pour un court papier de vulgarisation sur le sujet. Sur l’ajustement de lois de mortalité prospectives, je renvoie au très beau mémoire en ligne ici (ou  pour des choses très proches en anglais). Pour l’approche statique, je renvoie également ici pour une méthode simple d’ajustement de loi de Makeham. Pour les aspects dynamique, il y a des documents techniques ici ou .

Sinon pour les données que nous utiliserons en TD, elles sont ici (pour l’exposition, i.e. le nombre de vivants) et  (pour les décès). Il faut retravailler un peu les données, à cause d’une valeur “110+”, mais rien de plus simple,

> DEATH=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/R/FR-Deaths_1x1.txt",header=TRUE)
> EXPOSURE=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/R/FR-Exposures_1x1.txt",header=TRUE)
> DEATH$Age=as.numeric(as.character(DEATH$Age))
Message d'avis :
NAs introduits lors de la conversion automatique 
> DEATH$Age[is.na(DEATH$Age)]=110
> EXPOSURE$Age=as.numeric(as.character(EXPOSURE$Age))
Message d'avis :
NAs introduits lors de la conversion automatique 
> EXPOSURE$Age[is.na(EXPOSURE$Age)]=110

Coté code, je renvoie à un précédant billet sur le sujet, ici, et  pour le package (non officiel) de démographie développé par Rob Hyndman.

Enfin, à la fin du cours, on parlera un peu de modèles de durée. J’essayerai de taper quelques choses dans les jours à venir, mais côté TD, on utilisera ce que fait Frédéric Planchet (qui connaît fort bien le sujet), et qui se trouve en ligne ici. Les données sont en ligne ici, avec un peu de code R pour retraiter les donnéesici et , et un peu de code pour des modèles classiques, comme du modèle de Cox (ici) ou du Kaplan Meier ().

Les codes R utilisés en TD sont les suivants (avec en plus les graphiques obtenus),

> D=DEATH[DEATH$Year==an,]
> E=EXPOSURE[EXPOSURE$Year==an,]
> MU = D[,3:5]/E[,3:5]
> plot(0:110,log(MU[,1]),type="l",col="red")
> lines(0:110,log(MU[,2]),col="blue")

Avec ci-dessus les log des taux de décès instantannés pour l’année 1986.

> PH=matrix(NA,111,111) # ligne X, colonne H
> PF=PH
> for(x in 0:110){
+ PH[x+1,1:(111-x)]=exp(-cumsum(MU[(x+1):111,2]))
+ PF[x+1,1:(111-x)]=exp(-cumsum(MU[(x+1):111,1]))
+ }
> x=0
> plot(1:111,PH[x+1,],ylim=c(0,1),type="l",col="blue")
> lines(1:111,PF[x+1,],col="red")

pour la fonction de survie à la naissance (ci-dessus) ou à 50 ans (ci-dessous),

> somme=function(x){sum(x,na.rm=TRUE)}
> EH=EF=rep(NA,111)
> EH=apply(PH,1,somme)
> EF=apply(PF,1,somme)
> plot(0:110,EH,type="l",col="blue")
> lines(0:110,EF,col="red")

Ce qui correspond aux espérances de vie résiduelles.

> for(k in 1:3){
+ Z=MU[,k]
+ Z[is.nan(Z)==TRUE]=NA
+ MU[,k]=Z
+ }
> AGE=unique(DEATH$Age)
> ANNEE=unique(DEATH$Year)
> MUF=matrix(MU[,1],length(AGE),length(ANNEE))
> MUH=matrix(MU[,2],length(AGE),length(ANNEE))
> persp(AGE[1:100],ANNEE,log(MUH[1:100,]),
+ theta=-30,col="light green",shade=TRUE)
> persp(AGE[1:100],ANNEE,log(MUH[1:100,107:1]),
+ theta=-30,col="light green",shade=TRUE)
> image(AGE,ANNEE,log(MUH))

On a retourné le graph en mettant les années à l’envers, mais on peut aussi regarder la surface pour les femmes, ainsi que l’allure des courbes de niveau

> persp(AGE[1:100],ANNEE,log(MUF[1:100,107:1]),
+ theta=-30,col="light green",shade=TRUE)
> image(AGE,ANNEE,log(MUF))

Mais nous n’avons ici qu’une lecture statique des tables, et nous ne suivons pas vraiment un individu: en même temps que la personne vieilli, le temps passe…. Considérons une personne née en 1900, et regardons (ex post) l’évolution de sa probabilité de décès (en regardant non pas la probabilité de décéder à 18 d’une personne de 18 ans en 1900, mais en tenant compte du fait que la personne a fêté ses 18 ans en 1918).

> MUFP=MUHP=rep(NA,111)    # MU PROSPECTIF
> for(i in 1:111){
+     MUHP[i]=MUH[i,colonne+i]
+     MUFP[i]=MUF[i,colonne+i]
+ }
> plot(0:110,log(MUHP),col="purple",type="l")
> D=DEATH[DEATH$Year==annee.naissance,]
> E=EXPOSURE[EXPOSURE$Year==annee.naissance,]
> MU = D[,3:5]/E[,3:5]
> lines(0:110,log(MU[,2]),col="blue")

On peut aussi l’ajustement par un modèle à la Lee-Carter (via la fonction fournie par LifeMetrics)

> source("https://perso.univ-rennes1.fr/arthur.charpentier/fitModels.r")
> DEATH=DEATH[DEATH$Age<90,] # virer age>90
> EXPOSURE=EXPOSURE[EXPOSURE$Age<90,] # virer age>90
> XV=unique(DEATH$Age)
> YV=unique(DEATH$Year)
> ETF=t(matrix(EXPOSURE[,3],length(XV),length(YV)))
> DTF=t(matrix(DEATH[,3],length(XV),length(YV)))
> ETH=t(matrix(EXPOSURE[,4],length(XV),length(YV)))
> DTH=t(matrix(DEATH[,4],length(XV),length(YV)))
> WA=matrix(1,length(YV),length(XV))
> LCF = fit701(xv=XV,yv=YV,etx=ETF,dtx=DTF,wa=WA)
-1166100 -> -263940.2  ->-204175.2 
-188375.9 -> -185822.4  ->-184403.5 
-183223.3 -> -183069.8  ->-182350.3 
-181701.1 -> -181600.6  ->-181145.5 
-180729.6 -> -180671.5  ->-180371.8 
-180101.9 -> -180068.6  ->-179871.2 
> LCH = fit701(xv=XV,yv=YV,etx=ETH,dtx=DTH,wa=WA)
-2564403 -> -1028776  ->-953629.9 
-881334.2 -> -868279.3  ->-834413.6 
-804946.8 -> -802817.2  ->-786343 
-772775 -> -772007.5  ->-764501.7 
-758653.8 -> -758399  ->-755212 
-752931 -> -752851  ->-751620.2 
-750811.1 -> -750786.7  ->-750351.4 
> plot(LCF$x,LCF$beta1,type="l",col="red")
> lines(LCH$x,LCH$beta1,col="blue")

> plot(LCF$x,LCF$beta2,type="l",col="red")
> lines(LCH$x,LCH$beta2,col="blue")

> plot(LCF$y,LCF$kappa2,type="l",col="red")
> lines(LCH$y,LCH$kappa2,col="blue")

> MUF      = DTF/ETF
> MUH      = DTH/ETH
> MUFMODEL = MUF; MUHMODEL = MUH
> for(i in 1:nrow(MUF)){
+ for(j in 1:ncol(MUF)){
+ MUFMODEL[i,j] = exp(LCF$beta1[j]
+ +LCF$beta2[j]*LCF$kappa2[i])
+ MUHMODEL[i,j] = exp(LCH$beta1[j]
+ +LCH$beta2[j]*LCH$kappa2[i])
+ }}

soit, pour avoir un dessin,

> persp(LCH$y,LCH$x,log(MUH),
+ theta=45,col="light green",shade=TRUE)
> persp(LCH$y,LCH$x,log(MUHMODEL),
+ theta=45,col="light blue",shade=TRUE)

On peut ainsi comparer les taux bruts (en vert) avec les taux estimés (en bleu)

mais aussi regarder  se trouvent les erreurs le plus grandes,

Peut on faire l’économie du formalisme quand on parle d’extrêmes ?

Tous les blogs économiques saluent la parution en poche du joli petit livre de Daniel Zajdenweber, Economie des Extrêmes. En particulier, beaucoup de monde salue ce livre qui explique simplement des choses complexes…. Par exemple Alexandre dès 2001 “Passé le premier chapitre, un peu ardu, et qui nécessite du lecteur des connaissances de base en statistique et probabilités (notion de lois de probabilité, d’espérance, de variance…) qui décrit en termes littéraires les caractéristiques de ces lois, l’auteur applique ces résultats à un grand nombre de phénomènes concrets, et en tire les conséquences“. Mais peut-on parler d’économie des extrêmes sans être technique ?

Histoire que mon message ne soit pas déformé, je trouve passionnant ce petit livre introductif à la problématique des risques extrêmes (qui est un de mes dadas depuis quelques années) mais j’espère qu’il servira d’encouragement à une lecture d’ouvrages plus détaillés sur le sujet. Car la vulgarisation a des limites que l’on atteint vite quand on parle de sujets aussi complexes.

L’exemple que j’ai le plus étudié est celui des sinistres de perte d’exploitation (longuement évoqué par Daniel Zajdenweber dans son livre). Il y a quelques années j’avais utilisé cette partie du livre comme base pour faire un sujet d’examen pour le cours de “réassurance et grands risques” que je donnais alors à l’ENSAE1. Et malheureusement, mes compétences littéraires sont très limitées, donc je vais faire des maths. Dans le livre, la figure suivante est présentée,

qui correspond effectivement à la fonction tracée dès 1925 par Karl Gustav Hagstroem (j’avais souligné (ici) ses travaux précurseurs où l’intérêt de la loi de Pareto pour modéliser les très grands riques apparaissait pour la première fois). C’est en effet assez naturel: si on a une loi de Pareto, i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z01.png

alors on pourrait écrire, en passant au logarithme

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z02.png

Si on représente la version empirique, c’est à dire le nuage de points

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z03.png

alors pour une loi de Pareto, les points devraient être alignés suivant une droite, et la pente doit correspondre au paramètre de la fonction puissance. C’est visiblement l’idée exploitée ici.

Autrement dit, les pointillés ne sont un intervalle de confiance, mais juste un outils graphique pour se demander si la pente vaut 1, ou pas. Daniel Zajdenweber affirme que la pente doit ici être -1.

Le fait que la valeur soit unitaire ou pas a en effet un impact très important en terme d’assurabilité du risque de perte d’exploitation. Rappelons que pour une variable positive (et c’est le cas ici). Et si on a une telle loi de Pareto (de puissance unitaire), alors la prime pure d’un traité de réassurance, couvrant entre m et M s’écrit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z06.png

soit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z07.png

ce qui correspond aux calculs de Daniel Zajdenweber… Mais encore une fois “l’absence d’espérance mathématique de la distribution des sinistres” est une conclusion très forte sur laquelle on peut essayer de revenir.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z04.png

aussi ici

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z05.png

autrement dit, l’espérance est finie si la pente est strictement plus grande (en valeur absolue) que 1. Si la pente est inférieure (ou égale) à 1, le risque n’est pas assurable ! Ce qui est une conclusion très très forte pour les assureurs.
J’ai donc demandé à la FFSA la base de données utilisée ici, et pour éviter des problèmes d’inflation des coûts de sinistres entre 1992 et 2000. Si je prends tous les sinsitres, on obtient l’ajustement de Pareto suivant

soit une pente (en valeur absolue de 1.47). Mais encore une fois, l’ajustement de Pareto se fait sur les grands sinistres. Hill a proposé un estimateur très populaire pour estimer ce coefficient, où on ne prend en compte que les k observations les plus grandes, et on regarde l’estimation de la pente du graphique de Pareto pour ces quelques valeurs. On représente alors l’estimation en fonction du nombre de grands sinistres, ou du seuil définissant les graphs sinistres. Numériquement, en posant

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z08.png

on peut écrire comme estimateur de la pente

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z09.png

soit, en simplifiant le numérateur,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z10.png

tel que l’a construit Hill en 1975. Graphiquement, on a ici

Bref, la question est ce savoir si on atteint la valeur 1 pour les grands sinistres. Graphiquement, on a malgré tout envie de rejeter cette hypothèse.
Une solution peut être de faire un test statistique, basé sur de ratio de vraisemblance, comme le suggèrent Reiss & Thomas (2001) or Coles (2001). En fait, on peut même utiliser d’autres estimateurs que celui de Hill, comme celui obtenu en faisant un ajustement de loi GPD (Pareto généralisée
) sur la loi des Excès, ou une loi GEV sur des maximas par blocs (Generalized Extreme Value). On introduit alors la statistique de test suivante

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z12.png

et on regarde les p-value (ainsi que la correction de Bartlett à droite),

On peut aussi, plus simplement, estimer plusieurs coefficients de pentes pour des seuils différents, et regarder la borne supérieure de l’intervalle de confiance,

Bref, même si avec un des ajustements de loi GPD on hésite à retenir une pente unitaire, la plupart des tests rejettent cette hypothèse, et donc le risque de perte d’exploitation semble être assurable, d’espérance mathématique finie. Bref, les dessins c’est très bien pour faire passer une idée, mais ne retenir que ça pour en tirer des conclusions aussi fortes me laisser sceptique….

Incertitude(s) dans les triangles de paiements

Un petit billet aujourd’hui sur les méthodes de provisionnements, et sur les calculs d’incertitude à 1 an, ou à l’ultime, de montants de provisions pour sinistres restant à payer…. Ce sujet avait été évoqué ici (surtout dans les commentaires à l’époque).

  • Mack et incertitude à long terme 

Le mse dans le cadre du modèle de Mack est proposé dans la library(ChainLadder), via la fonction MackChainLadder(). Sur le triangle utilisé dans les exemples et formations en provisionnement, on obtient les mse estimés suivants

> source("https://perso.univ-rennes1.fr/arthur.charpentier/bases.R")
> library(ChainLadder)
> MackChainLadder(PAID)
MackChainLadder(Triangle = PAID)

  Latest Dev.To.Date Ultimate    IBNR Mack.S.E CV(IBNR)
1  4,456       1.000    4,456     0.0    0.000      NaN
2  4,730       0.995    4,752    22.4    0.639   0.0285
3  5,420       0.993    5,456    35.8    2.503   0.0699
4  6,020       0.989    6,086    66.1    5.046   0.0764
5  6,794       0.978    6,947   153.1   31.332   0.2047
6  5,217       0.708    7,367 2,149.7   68.449   0.0318

              Totals
Latest:    32,637.00
Ultimate:  35,063.99
IBNR:       2,426.99
Mack S.E.:     79.30
CV(IBNR):       0.03

Mais il s’agit fondamentalement d’une formule sur l’estimation ultime de la charge. Pour une année de survenance donnée, la formule dérivée par Thomas Mack ressemble à ça

(honteusement emprunté dans des slides de Mario Wüthrich, par paresse, ici), où pour rappel,

  • Muerz & Wüthrich et incertitude à court terme

Pour faire simple, il existe une formule fermée, un peu complexe, et une simplification, reposant sur l’idée que, si on note

alors, Michael Merz et Mario Wüthrich ont montré que

Le code sur le triangle mentionné ici donne le résultat suivant (ce code est parfois appelé dans divers documents MackMerzWuthrich.R)

> source("https://perso.univ-rennes1.fr/arthur.charpentier/merz-wuthrich-triangle.R")
> MSEP_Mack_MW(PAID,0)
   MSEP Mack MSEP observable approchée MSEP observable exacte
1  0.0000000                  0.000000               0.000000
2  0.6393379                  1.424131               1.315292
3  2.5025153                  2.543508               2.543508
4  5.0459004                  4.476698               4.476698
5 31.3319292                 30.915407              30.915407
6 68.4489667                 60.832875              60.832898
7 79.2954414                 72.574735              72.572700

La dernière ligne (notée ici 7 correspond au total – toutes années de survenance confondues). Ici, on regarde l’erreur à un an. Et Mario Wüthrich concluait – lors d’une jolie école d’été à Lyon il y a quelques années – en notant le point suivant

Sur le triangle étudié ici, on était à 91%. Mais supposer que l’on soit inférieur à 1 est un point qu’il me semble intéressant de discuter…

  • incertitude plus grande à court terme ?

Dans un billet datant un peu (ici), on laissait entendre qu’il était sûrement possible que l’incertitude à court terme excède l’incertitude à long terme.
Arnaud Lacoume a trouvé un vrai triangle dans lequel l’incertitude de court terme excède l’incertitude à long terme (détaillé ici dans une feuille xls). Pour information, il s’agit du risque GAV, garanties accidents de la vie, bref, rien de très exotique…

> Triangle = read.table("https://perso.univ-rennes1.fr/arthur.charpentier/GAV-triangle.csv",sep=";")

Si on utilise le code gentiment offert par Arnaud, on note qu’effectivement, l’incertitude à 1 an est grande, très grande…

> MSEP_Mack_MW(Triangle,0)
    MSEP Mack MSEP observable approchée MSEP observable exacte
1        0.00                       0.0                    0.0
2    12459.74                  129692.2               152605.9
3   209431.14                  214136.5               214419.6
4   258003.38                  198072.3               198773.0
5  3055297.40                 3048489.5              3065525.1
6 58429393.29                57056117.3             67375794.0
7 58669646.13                57301552.4             67586106.6

On note ici que le MSE exact à un an (car Merz et Wüthrich proposent aussi une formule approchée, basée sur une approximation qui contraint l’estimation à un an à être inférieure à l’estimation à horizon ultime) est – si on regarde toutes années de survenance confondues – 15% plus élevé que le MSE à ultime. Avec un résultat du même ordre pour la dernière année (ce n’est donc pas un problème d’agrégation). On peut jouer à simuler des triangles, via un processus de déroulé de paiements.

Regardons par exemple le code suivant pour simuler des triangles. L’idée est ici que l’on n’a pas d’incertude sur la valeur à l’ultime,

> M=matrix(NA,10,10)
> for(i in 1:10){
+ X=rexp(10)
+ Y=cumsum(X)/sum(X)
+ M[i,1:(11-i)]=Y[1:(11-i)]*100}
> Triangle = M

On tire des variables positives, car on veut un modèle à incréments positifs. Par exemple, une des simulations a donné le triangle suivant

> Triangle
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,] 35.4 36.3 38.2 47.6 50.7 66.1 71.3 72.1 98.7   100
 [2,] 15.5 20.2 22.3 23.7 63.9 68.2 75.9 79.3 95.1    NA
 [3,]  4.1 30.3 32.3 33.3 49.4 53.1 62.8 77.1   NA    NA
 [4,]  1.1 19.4 38.9 41.9 48.0 66.5 68.3   NA   NA    NA
 [5,]  0.3  7.2 10.8 17.8 28.1 30.1   NA   NA   NA    NA
 [6,] 15.6 32.2 46.0 50.6 68.1   NA   NA   NA   NA    NA
 [7,]  4.0  9.9 10.6 15.8   NA   NA   NA   NA   NA    NA
 [8,]  6.1 24.0 29.4   NA   NA   NA   NA   NA   NA    NA
 [9,]  0.0  0.5   NA   NA   NA   NA   NA   NA   NA    NA
[10,]  2.3   NA   NA   NA   NA   NA   NA   NA   NA    NA

En faisant le calcul sur les mesures d’incertitude à la Merz & Würthrich, on obtient une très grande incertitude à court terme,

> MSEP_Mack_MW(Triangle,0)
   MSEP Mack MSEP observable approchée MSEP observable exacte
1   0.000000                  0.000000               0.000000
2   6.820602                 12.776729             118.050196
3  13.377536                 13.119557              13.126625
4  17.415678                 13.638880              13.651131
5  11.797752                  5.574916               5.580840
6  26.767233                 18.407916              18.429069
7  27.184076                 23.725348              23.773581
8  44.100071                 16.410374              16.444281
9   7.227724                  3.450343               3.460982
10 54.012270                 49.443633              49.593284
11 99.184167                 85.042393             145.008282

Ici je suis 50% plus élevé ! Bon, mais c’est beau de montrer pourquoi une idée partagée par beaucoup de monde n’est pas toujours vraie…mais il reste à comprendre plus en détails le pourquoi … à suivre donc !

En attendant, je mets ici le code utilisé pour quantifier cette incertitude à un

### MERZ & WUTHRICH

MWChainLadder <- function(Trian, bool = 0, warning = FALSE){    

##  variables
I <- nrow(Trian)
J <- ncol(Trian)
diag <- matrix(NA,I,1)
diag_inv <- matrix(NA,I,1)
S_I <- matrix(NA,1,J)
S_II <- matrix(NA,1,J)
Phi <- matrix(NA,I,1)
Psi <- matrix(NA,I,1) 
Delta <- matrix(NA,I,1)
Gamma <- matrix(NA,I,1)
Lambda <- matrix(NA,I,1)
Upsilon <- matrix(NA,I,1)
cov_obs <- matrix(0,I,J)
cov_reel <- matrix(0,I,J)
msep_obs <- matrix(NA,1,J+1)
msep_reel <- matrix(NA,1,J+1)
Delta[1] <- Phi[1] <- Psi[1] <- Upsilon[1] <- Lambda[1] <- 0

Mack<-MackChainLadder(Trian)
Mack$sigma[J-1] <- sqrt(min(Mack$sigma[J-2]^4/
Mack$sigma[J-3]^2,min(Mack$sigma[J-2]^2,Mack$sigma[J-3]^2)))
#plot(Mack)

for (i in 1:I){
    diag[i] = Trian[i,I-i+1]
    diag_inv[i] = Trian[I-i+1,i]
}

for (j in 1:J){
    S_I[j] <- sum(Trian[1:(I-j),j])
    S_II[j] <- sum(Trian[1:(I-j+1),j])
}
S_I[I] <- 0

Delta[2] <- Mack$sigma[I-1]^2/(S_I[I-1]*(Mack$f[I-1])^2)
Phi[2] <- 0
Psi[2] <- Mack$sigma[I-1]^2/(diag[2]*(Mack$f[I-1])^2)
Upsilon[2] <- Mack$sigma[I-1]^2/(S_II[I-1]*(Mack$f[I-1])^2)
Lambda[2] <- diag[2]*Mack$sigma[I-1]^2/((Mack$f[I-1]^2)*S_II[I-1]*S_I[I-1])
for (i in 3:I){
    Delta[i] <- Mack$sigma[I-i+1]^2/(S_I[I-i+1]*(Mack$f[I-i+1])^2) + sum(
            (diag_inv[(I-i+2):(J-1)]/S_II[(I-i+2):(J-1)])^2*
         Mack$sigma[(I-i+2):(J-1)]^2/(S_I[(I-i+2):(J-1)]*(Mack$f[(I-i+2):(J-1)])^2))        
    Phi[i] <- sum((diag_inv[(I-i+2):(J-1)]/S_II[(I-i+2):(J-1)])^2*
           Mack$sigma[(I-i+2):(J-1)]^2/(diag_inv[(I-i+2):(J-1)]*(Mack$f[(I-i+2):(J-1)])^2))    
    Psi[i] <- Mack$sigma[I-i+1]^2/(diag[i]*(Mack$f[I-i+1])^2)
    Upsilon[i] <- Phi[i] + Mack$sigma[I-i+1]^2/(S_II[I-i+1]*(Mack$f[I-i+1])^2)

Lambda[i] <- diag[i]*Mack$sigma[I-i+1]^2/((Mack$f[I-i+1]^2)*S_II[I-i+1]*S_I[I-i+1]) + sum(
            (diag_inv[(I-i+2):(J-1)]/S_II[(I-i+2):(J-1)])^2*
           Mack$sigma[(I-i+2):(J-1)]^2/(S_I[(I-i+2):(J-1)]*(Mack$f[(I-i+2):(J-1)])^2))
}
Gamma = Phi + Psi

#MSEP par acc. year
for (i in 1:I){
    msep_obs[i] = Mack$FullTriangle[i,J]^2 * (Gamma[i] + Delta[i])
    msep_reel[i] = Mack$FullTriangle[i,J]^2 * (Phi[i] + Delta[i])
}

#Covariance
for (i in 2:(I-1)){ 
        for (k in (i+1):I){
        cov_obs[i,k] <- Mack$FullTriangle[i,J]*Mack$FullTriangle[k,J]*(Upsilon[i] + Lambda[i])
        cov_reel[i,k] <- Mack$FullTriangle[i,J]*Mack$FullTriangle[k,J]*(Phi[i] + Lambda[i])    
    }
}

#MSEP aggregated
msep_obs[I+1] = sum(msep_obs[2:I]) + 2*sum(cov_obs)
msep_reel[I+1] = sum(msep_reel[2:I]) + 2*sum(cov_reel)

##################################################################################################################
## EXACTES
##################################################################################################################
#DÈclaration des variables
facteur <- matrix(1,I,1)
Phi_exact <- matrix(NA,I,1)
Psi_exact <- matrix(NA,I,1)
Gamma_exact <- matrix(NA,I,1)
Upsilon_exact <- matrix(NA,I,1)
cov_obs_exact <- matrix(0,I,J)
cov_reel_exact <- matrix(0,I,J)
msep_obs_exact <- matrix(NA,1,J+1)
msep_reel_exact <- matrix(NA,1,J+1)
Phi_exact[1] <- Upsilon_exact[1]<- Gamma_exact[1] <- 0

# i=2
Phi_exact[2] <- Psi_exact[2] <- 0
Gamma_exact[2] <- Mack$sigma[I-1]^2/(diag[I-1]*(Mack$f[I-1])^2)
Upsilon_exact[2] <- Mack$sigma[I-1]^2/(S_II[I-1]*(Mack$f[I-1])^2)

# loop for prior
for (i in 3:I){
    facteur[i] <-prod(1 + diag_inv[(I-i+2):(J-1)]*(Mack$sigma[(I-i+2):(J-1)]^2)/((S_II[(I-i+2):(J-1)]*
            Mack$f[(I-i+2):(J-1)])^2))
    Phi_exact[i] <- (1 + Mack$sigma[I-i+1]^2/(diag_inv[I-i+1]*(Mack$f[I-i+1])^2))* (facteur[i]-1)    
    Psi_exact[i] <- (1 + Mack$sigma[I-i+1]^2/(S_II[i]*(Mack$f[I-i+1])^2)) * Phi[i] /
                (1 + Mack$sigma[I-i+1]^2/(diag_inv[I-i+1]*(Mack$f[I-i+1])^2))
    Upsilon_exact[i] <- (1 + Mack$sigma[I-i+1]^2/(S_II[I-i+1]*(Mack$f[I-i+1])^2))* facteur[i]-1    
    Gamma_exact[i] <- (1 + Mack$sigma[I-i+1]^2/(diag_inv[I-i+1]*(Mack$f[I-i+1])^2))* facteur[i]-1    
}

#MSEP per accident year
for (i in 1:I){
    msep_obs_exact[i] = Mack$FullTriangle[i,J]^2 * (Gamma_exact[i] + Delta[i])
    msep_reel_exact[i] = Mack$FullTriangle[i,J]^2 * (Phi_exact[i] + Delta[i])
}

#Covariance
for (i in 2:(I-1)){ 
    for (k in (i+1):I){
        cov_obs_exact[i,k] <- Mack$FullTriangle[i,J]*Mack$FullTriangle[k,J]*(Upsilon_exact[i] + Lambda[i])
        cov_reel_exact[i,k] <- Mack$FullTriangle[i,J]*Mack$FullTriangle[k,J]*(Psi_exact[i] + Lambda[i])    
    }
}

#MSEP ag.
msep_obs_exact[I+1] <- sum(msep_obs_exact[2:I]) + 2*sum(cov_obs_exact)
msep_reel_exact[I+1] <- sum(msep_reel_exact[2:I]) + 2*sum(cov_reel_exact)

##################################################################################################################
##Outputs Fonction 
##################################################################################################################
msep_Mack <- array(0,c(1,I+1))
msep_Mack[1:I] <- Mack$Mack.S.E[,I]
msep_Mack[I+1] <- Mack$Total.Mack.S.E
Vari <- array(0,c(1,I+1))
for (i in 1:I){
    Vari[i] <- Mack$FullTriangle[i,J]^2 * Psi[i]
}
Vari[I+1] <- sum(Vari[1:I])

# complete vision
if(bool==1){    
    result <- cbind(t(msep_Mack), t(sqrt(Vari)), t(sqrt(msep_obs)), t(sqrt(msep_obs_exact)), t(sqrt(msep_reel)), t(sqrt(msep_reel_exact)))
    result <- as.data.frame(result)
    l <- list("Variance CDR", "MSEP Mack","MSEP obs. approx.","MSEP obs. exact", "MSEP real approx.","MSEP real exact")
    names(result) <- l
    return(result)
    }

#si l'on souhaite uniquement l'aspect solvabilité
if(bool==0){
    result <- cbind(t(msep_Mack), t(sqrt(msep_obs)), t(sqrt(msep_obs_exact)))
    result <- as.data.frame(result)
    l <- list("MSEP Mack","MSEP obs. approx","MSEP obs. exact")
    names(result) <- l
    return(result)
    }
}

Les méchanismes bonus-malus

L’outil de base pour l’étude des systèmes bonus-malus est l’utilisation des chaînes de Markov, dès lors que l’on travaille sur des classes de bonus. Pour plus de compléments sur les chaînes de Markov, je renvoie à des vieilles notes de cours (inachevée, et sûrement remplis de coquillettes) ici.

  • Les chaînes de Markov

L’outil de base pour l’étude des systèmes bonus-malus est l’utilisation des chaînes de Markov, dès lors que l’on travaille sur des classes de bonus. Pour plus de compléments sur les chaînes de Markov, je renvoie à des vieilles notes de cours (inachevée, et sûrement remplis de coquillettes) ici.
Pour reprendre les grandes idées, on considère une suite de variables aléatoires à valeurs dans https://perso.univ-rennes1.fr/arthur.charpentier/latex/calE.png, https://perso.univ-rennes1.fr/arthur.charpentier/latex/suiteXn.png. Cette suite est sera une chaîne de Markov si pour tout https://perso.univ-rennes1.fr/arthur.charpentier/latex/ngeq1.png, et pour tout

https://perso.univ-rennes1.fr/arthur.charpentier/latex/suitei.png

telle que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/markov-proba1.png

on ait

https://perso.univ-rennes1.fr/arthur.charpentier/latex/proab1.png
= https://perso.univ-rennes1.fr/arthur.charpentier/latex/proab2.png

On parlera d’absence de mémoire (au sens où seule la dernière valeure intervient dans le calcul de la loi conditionnelle, et pas l’intégralité de la trajectoire). Si https://perso.univ-rennes1.fr/arthur.charpentier/latex/suiteXn.png est une chaîne de Markov de distribution initiale https://perso.univ-rennes1.fr/arthur.charpentier/latex/lambdaV.png et de matrice de transition P, pour tout https://perso.univ-rennes1.fr/arthur.charpentier/latex/nkentier.png,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/proba-markov-statioonnaire.png
https://perso.univ-rennes1.fr/arthur.charpentier/latex/pr135.png

Pour étudier l’évolution sur le long terme, on peut utiliser l’équation de Chapman-Kolmogorov: si https://perso.univ-rennes1.fr/arthur.charpentier/latex/suiteXn.png  est une chaîne de Markov de distribution initiale  https://perso.univ-rennes1.fr/arthur.charpentier/latex/lambdaV.png et de matrice de transition P,  pour tout
https://perso.univ-rennes1.fr/arthur.charpentier/latex/nkentier.png, pour tout https://perso.univ-rennes1.fr/arthur.charpentier/latex/zvze54v6.png,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/chapman-kol.png

Un corrolaire important est que si

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ezjpjio.png

alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcnaoli.png

Bref, on peut facilement étudier l’évolution à long terme de la chaîne de Markov en faisant un peu de calcul matriciel, i.e. des calculs de puissances de matrices.
En fait, formellement, on peut soit étudier des matrices, soit des graphes. Par exemple,

Si on veut étudier l’évolution sur deux périodes, il suffit de calculer le carré de la matrice, i.e.

  • Les systèmes bonus-malus par classes

Les tableaux ci-dessous sont tirés des annexes de la bible de Jean Lemaire, Bonus-Malus systems in automobile insurance. Le plus simple est probablement celui de Hong Kong, qui est Markovien (à l’ordre 1),

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.BM-hong-kong-B-9_m.jpg

L’exemple ci-dessous est celui qui a existé en Allemagne dans les années 80,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.BM-germany-B-7_m.jpg

puis celui qui existait en Finlande,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.BM-finnish-B-6_m.jpg

et enfin le système italien du début des années 90,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.BM-italy-B-11_m.jpg

Pour tous ces sytèmes, une fois spécifiée la loi du nombre de sinistres par an, il est possible de construire la matrice de transition des classes de bonus sur un an, ou plus.
En particulier, l’étude de la loi invariante de ce genre de sytème (que je mentionnais en introduction) permet d’étudier la viabilité du système à plus ou moins long terme, et la vitesse de la convergence vers l’état stationnnaire.

  • Le système bonus-malus français

En France, le mécanisme de bonus-malus peut paraître un peu plus compliqué, et il fonctionne de la manière suivante.

Le calcul du bonus-malus s’ef
fectue en fonction des sinistres survenus (ou non) pendant la période de douze mois consécutifs précédant de deux mois l’échéance annuelle du contrat d’assurance.
Seuls les accidents responsables donnent lieu à malus. C’est à dire qu’en cas de foce majeur, ou si l’accident est entièrement imputable à un tiers, alors le conducteur n’aura pas de malus. Notons aussi que Le malus est rattaché au contrat d’assurance. Autrement dit, si quelqu’un cause un accident avec un véhicule assuré par un tiers, il se supportera pas le malus.
La prime de référence est la prime qui sera utilisée pour ensuite actualiser le montant de prime en tenant compte du coefficient de majoration ou minoration. Il s’agit de la prime technique calculée par l’assureur en fonction de charactéristiques a priori, comme la puissance du véhicule véhicule, la zone géographique de circulation, l’usage socioprofessionnel ou encore le kilométrage parcouru….
A partir de cette prime de référence, on calcule une prime qui tient compte de l’expérience passée, via un coefficient multiplicatif.
La valeur de référence est de 1 pour le taux de majoration-minoration. Ce taux subie une réduction – corrrespondant à un bonus – pour une année sans accident, d’un niveau de 5 %. Toutefois, la valeur du coefficient est arrêté à la deuxième décimale et arrondi par défaut. Au final, le taux prend donc un nombre fini de valeurs. Pour résumer schématiquement

1 année sans accident 100 x 0.95 = 0.95 de Bonus
2 années sans accident 0.95 x 0.95 = 0.90 de Bonus
3 années sans accident 0.90 x 0.95 = 0.85 de Bonus
4 années sans accident 0.85 x 0.95 = 0.80 de Bonus
5 années sans accident 0.80 x 0.95 = 0.76 de Bonus
6 années sans accident 0.76 x 0.95 = 0.72 de Bonus
7 années sans accident 0.72 x 0.95 = 0.68 de Bonus
8 années sans accident 0.68 x 0.95 = 0.64 de Bonus
9 années sans accident 0.64 x 0.95 = 0.60 de Bonus
10 années sans accident 0.60 x 0.95 = 0.57 de Bonus
11 années sans accident 0.57 x 0.95 = 0.54 de Bonus
12 années sans accident 0.54 x 0.95 = 0.51 de Bonus
13 années sans accident 0.51 x 0.95 = 0.48 plafonné à 0.50

La réduction maximale est donc de 0.50, soit 50% de réduction sur la prime de référence.
En revanche, chaque accident responsable donne lieu à une majoration – on palera de malus – de 25 % du coefficient par accident. Toutefois, cette majoration est réduite de moitié (12,5 %) si la responsabilité n’est que partiellement engagée.Là encore, le coefficient obtenu est arrêté à la deuxième décimale et arrondi par défaut. Aussi, avec deux accidents responsables dans l’année, le coefficient sera multiplié par un coefficient de 1,56. Afin d’éviter l’absence d’assurance, on notera qu’en aucun cas, le coefficient de réduction-majoration ne peut excéder 3,50.
Pour les jeunes conducteurs, l’article A. 335-9-1 du code des assurances prévoit que l’assureur a la possibilité de majorer la prime de référence (dans des limites fixées par la loi là encore).
Notons enfin qu’il existe des avantages pour les bons conducteurs: si un conducteur est au bonus maximal depuis au moins trois ans, il conserve son bonus après le premier accident responsable. Notons que c’est le genre de disposition qui figure dans la loi et que certains assureurs utilisent à des fins commerciales (comme le montre l’affiche ci-dessus ou ci-dessous, qui a été analysée sur le blog de Gizmo (ici)).

Il existe aussi une règle dite de descente rapide qui peut atténuer les effets des malus: après deux années consécutives sans accident, votre coefficient ne peut être supérieur à 1.
Ce qui rend le système intéressant, c’est qu’il est obligatoire, en France. Aussi, en cas de changement d’assureur, le coefficient de bonus-malus sera recalculé à partir d’un relevé d’information: l’assureur est tenu de fournir ce document à la résiliation du contrat (contenant des informations sur le conducteur et le véhicule, ainsi que le nombre des sinistres survenus au cours des cinq années précédant l’établissement du relevé d’information, ainsi que leur nature et la part de responsabilité retenue)
Davantage d’information se trouve ici, (ce document propose d’ailleurs une petite simulation à la fin sur le coût d’un accident) ou encore (pour une approche plus actuarielle), voire (pour une approche plus économique).

De l’interprétation d’un effet nonlinéaire en régression

Je vais poursuivre un peu le TD de statistique de l’actuariat 2, où nous avions parlé tarification et régression Poissonnienne.

  • Analyser un effet nonlinéaire

Dans un premier temps, nous avions regardé un modèle simple, où l’on voulait voir si l’âge du conducteur expliquait la fréquence de sinistres. Bref, on peut utiliser une régression log-Poisson, qui donne la prédiction suivante, en fonction de l’âge,

Mais l’âge est une variable un peu particulière car c’est une fausse variable continue: seule les valeurs entières apparaissent. Comme on chercher à calculer

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-glm-nb-1.png

notons qu’il est possible de considérer l’estimateur empirique naturel de cette grandeur, i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-glm-nb-2.png

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-glm-nb-3.png

Graphiquement, on obtient la courbe suivante,

Le côté ératique de la série à droite et à à gauche vient de la faible représentation dans la base, comme en atteste l’histogramme de l’âge du conducteur,

Je peux en profiter pour revenir à un commentaire posté ici sur les modèles GAM. En effet, au lieu de regarder la fréquence en fonction de l’âge, on peut s’intéresser à la la différence par rapport à la fréquence globale.

> m=mean(base$nbsin)
> Y=predict(REG,newdata=data.frame(age=X,expo=E),type = "response")
> lines(X,Y/m-1,col="red")

On obtient alors le graphique suivant (en superposant à l’ajustement obtenu par GAM)

Notons que si l’on ajuste un modèle log-Poisson uniquement sur la tranche d’âge [50-70], on est très proche de ce que propose le modèle GAM, d’où l’interprétation de l’ajustement GAM comme un ajustement local.

 

  • Changer de modèle (et les comparer)

On a choisi un modèle log-Poisson, mais on se doute que d’autres modèles seraient possible. En particulier dans notre base, les assurés ont eu 0 ou 1 sinistre l’année passée. Autrement dit, des lois binomiales peuvent être ajustée,

> REG0 = glm(nbsin~age,base,
+  family=poisson(link="log"),
+  offset=expo)
> REG1a = glm(nbsin~age,base,
+  family=binomial(link = "logit"),
+  offset=expo)
> REG1b = glm(nbsin~age,base,
+  family=binomial(link = "probit"),
+  offset=expo)
> REG1c = glm(nbsin~age,base, 
+  family=gaussian(link = "identity"),
+  offset=expo)

Le critère AIC est souvent invoqué pour juger de l’adéquation d’un modèle. Il vaut ici

> AIC(REG0)
[1] 36470.7
> AIC(REG1a)
[1] 36044.45
> AIC(REG1b)
[1] 37304.02
> AIC(REG1c)
[1] 32549.83

Mais on peut également regarder le critère BIC,

> AIC(REG0, k = log(nrow(base)))
[1] 36488.81
> AIC(REG1a, k = log(nrow(base)))
[1] 36062.57
> AIC(REG1b, k = log(nrow(base)))
[1] 37322.14
> AIC(REG1c, k = log(nrow(base)))
[1] 32577

Autrement dit, en utilisant le règle “the smaller, the better” on retiendrait le dernier modèle qui est un modèle Gaussien, manifestement assez mauvais: sur le graphique ci-dessous, la régression log-Poisson est en mauve, presque confondu avec la régression logistique1, en bleu, la régression probit en vert et le modèle linéaire (Gaussien) en rouge,

Bref, ce dernier modèle n’est pas convainquant… Par contre, notons que si l’on compare le modèle log-Poisson et logit, ce dernier conduit à tarifer moins cher les âges extrêmes,

Autrement dit, au lieu de comparer des modèles actuariels sur la base d’outils statistiques, et uniquement de ces outils, il peut être instructif d’essayer de comprendre qui verra sa prime baisser avec tel ou tel modèle… Mais je reviendrais plus longuement, une autre fois, sur les histoires de choix de modèles…..

Crédibilité et calcul bayésien

On a commencé en cours la tarification a posteriori, en commençant par la tarification bayésienne, et en présentant ensuite l’approche d’Hans Bühlmann comme un cas particulier en projetant dans un sous espace de combinaisons linéaires (un peu comme l’économétrie linéaire était présenté en cours comme un cas particulier de calcul d’espérance conditionnelle).La classification a priori sera poursuivie un TD, sur des bases de données, afin d’orienter davantage le cours sur la pratique des GLM.
L’utilisation des techniques bayésiennes en actuariat ne s’est pas faite sans peine. Pour ceux qui souhaitent s’en convaincre, il suffit de relire les comptes rendus de la conférence à Trieste en 1962. Dans le discours de Bruno de Finetti (publié dans l’ASTIN Bulletin), on retrouve un point historique, présentant les nouveaux travaux d’Hans Bühlmann,

Et Bruno de Finetti de conclure en rappelant les apports d’Arthur Bailey en 1950,

  • La formule de Bayes, ou apprendre par l’expérience

La formule de Bayes ne dit pas grand chose quand on le regarde sous la forme suivante

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-bayes-01.png

A la rigueur on peut aussi la réécrire en utilisant deux fois cette formule

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-bayes-02.png

Bon, là encore, on ne voit peut être toujours pas ce qu’on peut en faire…. L’idée géniale est de faire la lecture suivante

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-bayes-03.gif

Autrement dit la loi a posteriori du paramètre, conditionnellement aux observations (en bleu) est proportionnelle au produit de  la loi a priori du paramètre (en vert) et de la vraisemblance des observations (en rouge, qui sera un simple produit de lois si on suppose l’indépendance conditionnelle au paramètre). Le paramètre en noir est juste un facteur de normalisation, histoire que la loi soit effectivement une mesure de probabilité. D’où parfois l’écriture

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-bayes-4.gif

Bref, on a réussit à inverser les deux termes dans la probability (on retrouve le terme “inverse probability” dans la littérature anglo-saxonne).
En cours, on avait vu que si on supposait le nombre annuel de sinistres Poisson, conditionnellement au paramètre, avec une indépendance (toujours conditionnellement à ce paramètre), et que l’on suppose que la loi a priori du paramètre est une loi Gamma, alors la loi a posteriori se calcule facilement et c’est une loi Gamma.
Considérons la petite simulation suivante: on suppose que le nombre de sinistres suit une loi de Poisson, et on cherche à estimer le paramètre de cette loi. On commence par une loi a priori Gamma, centrée sur 5, et on regarde évoluer la loi a posteriori au fur et à mesure que l’on observe les sinistres arriver (ou les années passer en l’occurence).

On peut également considerer l’évolution de l’espérance a posteriori du paramètre, ainsi qu’un intervalle de confiance (puisqu’on connaît la loi)

La “magie” des méthodes bayésiennes est que peu importe la loi a priori, on va converger vers la vraie valeur. Par exemple sur l’animation ci-dessous, on part d’une loi a priori de moyenne 2,

avec là aussi l’évolution de la moyenne a posteriori avec un intervalle de confiance,

La morale de l’histoire est que comme la loi a priori n’influence pas (j’exagère un peu car elle va influencer la vitesse de la convergence, ce qui n’est pas rien en pratique), autant prendre quelque chose de simple pour les calculs: la loi conjuguée de la loi exponentielle. La loi conjuguée d’une loi de la famille exponentielle de la forme

s’écrit

où et sont des constantes, et où  est une constante de normalisation (afin que cette fonction soit effectivement une mesure de probabilité).

  • L’hérésie de la statistique bayésienne

La principale hérésie des statisticiens bayésien a été de proposer une lecture des probabilités d’un point de vue subjectif. En particulier Jimmie Savage et Bruno de Finetti partait de l’hypothèse que chacun avait à sa disposition son propre espace de probabilité.C’est ce qui lui permet de se construire un a priori. Et au fur et à mesure que le temps passe, on apprend, et on peut dériver une loi a posteriori. Par exemple  je pars d’un a priori comme quoi un enfant sur deux qui nait est un garçon, en France. Je peux ensuite récupérer des données, et affiner ma connaissance. C’est ce qu’avait fait Laplace il y a fort fort longtemps….
L’hérésie est ici un
problème philosophique lié au concept même de probabilité. Mais on peut noter qu’à l’époque dans les années 60, les débats étaient assez virulants. Et pour Jimmie Savage et Bruno de Finetti, il n’est pas nécessaire d’avoir le formalisme mathématique d’Andreï Kolmogorov pour appréhender la notion de probabilité (c’est à dire la mesure d’ensembles dans un espace abstrait d’évènements, voir ici ou ). On notera que pour John Keynes (qui avait écrit un “traité de probabilités“) il existait une probabilité subjective, mais davantage liée à une information donnée par l’éducation. C’est un problème culturel plus que personnel. Mais finalement, il a surtout fallu attendre les travaux de Daniel Kahneman et Amos Tversky dans les années 70 pour que des psychologues se penchent sur la construction individuelle des probabilités.

  • La crédibilité au sens de Bühlmann

L’idée d’Hans Bühlmann est de proposer une formule fermée qui approche l’espérance conditionnelle, tout simplement en utilisant celle qui avait été obtenue dans le cas bayésien. En tarification, on s’intéresse à des primes pures, et il est finalement assez naturel de chercher à calculer

Hans Bühlmann a suggéré d’approcher cette grandeur par

avec classiquement

est appelé facteur de crédibilité, avec , et où respectivement,


Bref, on a un facteur de crédibilité qui dépend de la variance intra et de la variance inter, pour reprendre des termes connus en statistique. L’énorme avantage est que cette formule peut très facilement s’utiliser sur des données. En particulier, considérons un tableau de la forme suivante

où on regarde si un contrat a – ou n’a pas – eu de sinistre chacune des années. Formellement, on dispose d’un tableau d’observations,

(on notera que le nombre d’années d’observations pour chacun des contrats peut être différent). On peut alors obtenir les estimateurs empiriques des moments nécéssaires, i.e. on pose

correspondant à la moyenne empirique pour un contrat,

et

Le facteur de crédibilité pour le contrat i est alors

où tout naturellement

Formation ERM de l’Institut des Actuaires

Intervention jeudi après midi dans le cadre de la formation ERM (Enterprise Risk Management) de l’Institut des Actuaires, sur le thème “everything you wanted to ask about tail dependence in risk management and you’re afraid to ask“.
Comme Anne-Laure Fougères est déjà intervenue sur les extrêmes, et que Stéphane Loisel ont présenté les copules, le cours sera une discussion autour des diverses recommandations du Groupe Consultatif, des documents du CEIOPS, et de divers documents qui circulent sur le calcul du SCR (Solvency Capital Requierment) de l’agrégation des risques. Ou “everything you wanted to ask about tail dependence in risk management and you’re afraid to ask“… Les slides sont en ligne ici, et ne contienent que les éléments illustratifs du cours.

Les généralités sont généralement fausses

Suite à un commentaire sur un vieux billet (ici), concomitante avec une demande d’intervention pour la formation ERM proposée par l’Institut des Actuaires, je me suis replongé dans les calculs de SCR (solvency capital requierement). Plus je lis la documentation technique, plus je suis surpris par l’importance de l’hypothèse de normalité ! J’ai repris mon Loss Model, qui est le livre de base pour l’examen d’actuariat non-vie dans les pays anglo-saxons, et j’ai vu des dizaines de lois bizarres… mais rien sur la loi normale…

  • de l’hypothèse de normalité dans la formule d’agrégation

Comme je le disais dans un vieux billet (ici), on ne peut montrer la validité de la formule de calcul du SCR que dans le cas Gaussien (ou plus précisément dans le cas elliptique, voire ici). En effet, la preuve est assez simple. Je renvois vers un papier tout chaud de Laurent Devineau et Stéphane Loisel sur le sujet, ici.

  • où l’on continue à entrevoir l’hypothèse de normalité

En fait l’hypothèse de normalité est vraiment partout dans les calculs de SCR. Mais bien sûr, sans le dire. Par exemple dans un papier de la SCOR,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.doc14-1_m.jpg

Comme je le disais l’autre jour (ici), le problème des documents “professionels” (par opposition à “académiques“), c’est que comme on veut que le CEO le lise, on évite de parler d’hypothèses, ou pire, de mettre un peu de maths…. Bref, on apprend que “generally” (terme d’une grande rigueur scientifique, que je recaserai un jour dans mes démonstrations), “the expected shortfall[…] at the 99 % level corresponds quite closely to the[…] value-at-risk at a 99.6% level “. Mais on retrouve la même phrase dans un rapport de l’edhec

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.SCR-edhec_m.jpg

Là aussi, aucune source… Damned ! En creusant un peu, Actuaris m’a mis sur la piste (ici), car on apprend que c’est l'”autorité suisse” qui a énoncé ce résultat que je qualifierais de surprenant1.

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.actuaris-var_m.jpg

Effectivement, on peut retrouver le document originel,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.doc13-1_m.jpg

J’avoue avoir essayé – sans grand succès – de lire l’annexe D. Je n’ai pas tout compris, mais j’y ai lu un “assuming that the distribution functions follows a normal law”,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.doc13-2_m.jpg

Donc si j’ai bien tout compris, “generally”, “the distribution functions follows a normal law” ! Bon, je crois que je vais pouvoir envoyer un mail à Stuart, Harry et Gordon pour leur expliquer qu’ils peuvent virer des pages à leur livre ! Ou peut-être Jean-Luc Besson (pour lequel j’ai une immense estime !) qui présente beaucoup de lois utiles en assurance dans l’ouvrage qu’il a coécrit avec Christian Partrat…
Mais peut-être a-t-il raison…? On peut se dire qu’en actuariat, on avait l’habitude de présenter l’approximation “normal power“, et peut être que ce résultat est vrai pour d’autres lois….
Regardons vite-fait sur quelques lois….. Par exemple, pour la loi lognormale, si on regarde en fonction de l’écart-type de la loi normale sous-jacente,

Effectivement, on colle au quantile à 99.6%. L’ordonnée est en échelle logarithmique pour que le dessin soit présentable, avec en pointillé la VaR (formule exacte, à 99.6% en rouge, et à 99.8% en bleu), et les ronds sont des estimations de la TVaR à 99%. Par contre, par exemple pour la loi lognormale, si on regarde en fonction de l’écart-type de la loi normale sous-jacente, plus la volatilité est grande, plus l’approximation est mauvaise.

Bon, essayons une loi Gamma (classique en assurance, car c’est une loi de la famille exponentielle, et souvent comparée à la loi lognormale).

Visiblement, là aussi l’approximation semble marcher….  mais rappelons que la loi Gamma est une loi à queue fine… Allons-y franchement, comparons avec le cas d’une loi de Pareto…

cette fois, l’abscisse se lit dans l’autre sens: à gauche, c’est le cas le plus risqué (avec une variance infinie entre 1 et 2). On notera que pour les risques les plus importants, là aussi on sous-estime la TVaR en considérant la VaR (à un seuil plus important). On retrouve d’ailleurs exactement le même type de graphique avec une loi de Student,

Bref, on me dira que je pinaille, n’empêche qu’à mon avis, les compangies sous-estiment “généralement” la TVaR (à supposer que ça soit la mesure de risque qui les intéresse) en prenant une hypothèse de normalité…

1 en tant qu’enseignant, ce genre de conclusion me gêne toujours, car je me fatigue à faire des cours sur la mesure de risque, j’essaye de présenter des notions aussi proprement que possible, en expliquant que les hypothèses c’est important, qu’il faut définir différentes notions de mesures de risques (entre les convexes, les cohérentes, les dynamiques, etc), et là tout d’un coup, on nous apprend que “generally” c’est kif kif….

Grippe, cravate et (ré)assureurs

Entendu sur France Info, dans la revue de presse de Jean Christophe Martin, ce matin, ici, “parmi ces conseils, une bonne nouvelle pour les allergiques à la cravate : en plus des recommandations classiques, se laver les mains et se dire bonjour sans bisous, le Parisien signale que le groupe Scor, 5ème société mondiale de réassurance, demande aussi à ses salariés d’éviter de porter une cravate. Une consigne qui a eu un gros succès dans les bureaux de la Défense, au siège du groupe. Et oui, contrairement aux chemises, on ne fait pas souvent nettoyer sa cravate, qui peut se transformer en boîte à postillons et en bouillon de culture… Le nouveau mot d’ordre dans les bureaux : à bas la cravate, vive la grippe A...”. Chez AXA, on apprend que c’est la bise, qui est interdite (ici).

Quand les assureurs souhaiteraient faire la pluie et le beau temps

Damned, je me suis fais devancé par un très très beau billet d’Anne (ici) sur l’assurance proposée par AXA et AON. Ca m’apprendra à ne plus lire régulièrement Le Monde

Tout comme Anne, je suis un peu agacé par le nom donné sur le site d’AON: un assureur ne peut pas garantir la pluie et le beau temps, de même qu’il ne peut pas garantir que vous ne tomberez pas malade ou que vous n’aurez pas d’accident de voiture dans l’année ! Ou alors c’est autre chose que de l’assurance…

Je suis également un peu sceptique car ces risques ne sont pas vraiment indépendants, surtout au niveau local (dans un même département, il y a des chances pour que tout le monde soit sinistré la même semaine). Peut être qu’au niveau national, il existe des dipôles climatiques qui font que s’il pleut dans un département, il y a de très fortes chances pour qu’il ne pleuve pas dans un autre. C’était une étude qu’on souhaitait mener avec des collègues de Paris 6, et nous n’avions riens trouvé sur le sujet. Mais les gens d’AXA et de Metnext savent ce qu’ils font, je pense… Car il faut reconnaître que les garanties (ici) excluent beaucoup de choses si on regarde en détails le contrat, en particulier quand on lit qu’il faut “3 heures consécutives” de pluie, “entre 8 heures et 20 heures“…etc.
En tous les cas AON a l’honnêteté d’annoncer clairement la couleur sur son site (ici). Tout d’abord avec de joli petit logo que j’ai reproduit ci dessous, avec cette superbe accroche (évoquée en détails ici) cette garantie est, pour les voyagistes, “une source de revenus complémentaires“…
Sinon il faut reconnaître également que les assureurs offraient déjà ce genre de garanties, mais seulement dans les grands évènements, comme des championnats du monde de football, des tournois de tennis, ou des concerts de grosses stars….pour la première fois, la garantie n’est plus pour Michael Jackson qui souhaiterait remplir des stades, mais pour Robert et Raymonde qui voudraient faire du camping au sec….

Un petit billet sur la racine carrée

En lisant le Canard l’autre jour, je suis tombé sur ce “joli” papier sur la règle de la racine carrée.

Bon, en l’occurence, c’était juste pour évoquer la jolie forme (i.e. la forme du symbole mathématique, aucunement la fonction qui se trouve derrière).
En fait, je voulais en profiter pour revenir sur la square too rule for reinsurance pour reprendre l’expression de Michael Powers. Je connaissais (plus ou moins) une loi de la racine carré en assurance, correspondant au théorème central limite (que j’ai évoqué ici via un papier d’Alain Tosetti), que l’on connait en statistique sous la forme(sous quelques hypothèses)

Mais ici le résultat est de toute autre nature. En faisant un petit modèle de théorie des jeux1, Michael Powers et Martin Shubik avaient étudié, en 2001, la concentration dans le marché de la réassurance (ici). Ils avaient alors élaboré cette square too rule for reinsurance qu’ils ont par la suite validé sur des données (ici).
En supposant que les assureurs ont des fonctions d’utilité CARA (homogènes) alors que les réassureurs sont neutres au risque, que les risques soient indépendants et que les couvertures soient proportionnelles, ils montrent que, pour éviter la saturation du marché, le nombre optimal de réassureur doit être en racine carrée du nombre d’assureurs. En fait, formellement, ce n’est pas trop une “racine carrée“, mais le résultat d’un programme d’optimisation de la forme suivante,

mais si on résout numériquement, on tombe sur une fonction racine carrée, ou presque,

(en tous les cas, la solution au programme semble unique, et asymptotiquement, on a effectivement une racine carrée du nombre d’assureurs). Venezian, Viswanathan et Juca ont testé (ici) cette relation en 2005 sur les données OCDE (que l’on retrouve dans les statistical yearbooks), c’est à dire 18 à 20 pays sur une dizaine d’années. Ils ont alors testé un modèle économétrique simple en passant au log, i.e.

et en étudiant l’évolution de la pente (correspondant à la puissance sur les nombres), les auteurs obtiennent le graphique suivant,

Sauf que derrière se cache peut être aussi des effets prix, ce qui est corrigé par la suite en prenant un modèle de la forme suivante

où le prix est obtenu en utilisant un modèle à la Madsen et Pedersen (2003, ici puis , voire ici ou pour des slides), qui a l’allure suivante

Dans ce cas en utilisant des GMM2 pour estimer les coeffiicients, Michael Powers et Martin Shubik ont obtenu un coefficient de l’ordre de 0,477, ce qui validerait l’hypothèse de loi en “racine carrée“.
Remarque: en notant que la rétrocession était un cas particulier de réassurance, ils évoquent une loi en puissance 1/4 pour le nombre de rétrocessionnaires.
1 je ferais sûrement un billet sur l’utilisation de la théorie des jeux en assurance, et en réassurance puisque je dois faire un exposé fin août à Grenoble sur ce sujet aux journées AST&Risk. A suivre donc….
2 promis, je ferais un cours billet sur les GMM…. ce n’est pas la première fois que j’en parle, mais pour plus d’information, je peux renvoyer ici, pour wikipedia, ou pour des slides de Bernard Salanié, voire enfin ici pour le chapitre dans l’incroyable livre de Russell Davidson et James G. MacKinnon (Estimation et Inférence en Économétrie).

Quantile de sommes, ou somme de quantiles ?

Je reprends ici un commentaire que j’ai longtemps entendu dans mon expérience d’actuaire dans une vie antérieure (et qu’on peut lire – entre les lignes le plus souvent – dans certaines réactions de risk managers). Les quantiles sont la base des mesures de risques, en particulier via la VaR en finance, mais quand j’étais actuaire, on calculait des 90/10 pour désigner des quantiles à 90% des pertes. Donc tout ça n’est pas spécialement nouveau. Je rappelle que la fonction quantile est simplement  où

Les quantiles sont calculés par type de risque (en gestion des risques), et forcément à la fin, on souhaite un quantile de la somme des risques, et il a longtemps été coutume de prendre la somme des quantiles. Je vais donc discuter un peu le résultat suivant,

  • L’égalité n’est pas vraie dans le cas indépendant

Suite à une espère de déformation professionnelle (à cause de la variance), beaucoup de monde pense que

mais ce résultat est faux…. Un exemple est d’ailleurs donné ci-dessous dans le cas de variables exponentielles.

  • L’égalité est vraie dans le cas comonotone…

En fait, le résultat précédant est vrai dans le cas comonotone. La comonotonie correspond au cas de dépendance parfaite positive, de corrélation maximale (si la corrélation existe), i.e.  avec  une fonction strictement croissante. Alors

  • … mais ce n’est pas un “worst case

Seconde bizarrerie, la comonotonie n’est pas le pire des cas. La comonotonie est souvent envisagée comme un worst case scenario mais il n’en est rien (cette idée se retrouve dans les QIS 3 par exemple, ici).
Formellement, si  (la classe de Fréchet), alors

et

En passant en dimension 2 (cas où la borne inférieure est effectivement une vraie fonction de répartition), si et  désignent respectivement des version anticomonotone et comonotone de , alors on peut se demander, en considérant une mesure de risque quelconque  si

En fait, Andre Tchen a montré en 1980 un résultat qui montre que la comonotonie est effectivement un worst case, mais seulement dans certains cas très limités.
Si est une fonction supermodulaire, i.e.

pour tout  et . Alors dans ce cas, pour tout ,

(la preuve se trouve ici). Je renvoie aux papiers de Michel Denuit, de Jan Dhaene ou de Marc Goovarts (ici ou ) pour des applications en tarification sur deux têtes, par exemple, ou sur les primes stop-loss en réassurance

  • Borner (numériquement) le quantile d’une somme

Commençons par un exemple, avec deux lois exponentielles, car c’est le plus simple, et surtout on peut calculer explicitement les bornes (dans le cas général, on se contentera de méthode numérique comme l’avait fait Khalil ici).
Si  et , où, pour être plus précis,

l’encadrement

est valide pour tout  quelle que soit la structure de dépendance entre et , où

Aussi, en prenant les inverses (car ces fonctions sont strictement croissante), on obtient une inégalité en terme de quantile, ou de Value-at-Risk,

pour tout .
Rappelons que dans le cas de variables indépendantes  alors que dans le cas comonotone . La figure ci-dessous montre les valeurs possible pour  la fonction de répartition de la somme, avec le cas indépendant en bleu, et le cas comotonone en rouge.

On peut aussi visualiser les valeurs possibles pour les quantiles de la somme. Dans un cas général, Ce problème a été traité par Makarov en 1981, et par Frank, Nelsen et Schweizer en 1997 (ici), tout en restant en dimension 2. En fait, ils ont travaillé sur les bornes de la fonction de répartition d’une somme mais c’est pareil, comme on l’a vu juste auparavant. Ils ont montré le résultat suivant: si  alors pour tout

désigne la copule anticomonotone (qui est une copule en dimension 2),

en posant , et

On peut retrouver des résultats plus fin en creusant dans la direction de l’arithmetic probability, par exemple avec la thèse de Williamson (ici, cette thèse est remarquable, et trop peu citée, malheureusement). Je ne parle que de la dimension 2, l’extension en dimension supérieure est délicate (je renvoie aux travaux de Paul Embrechts sur le sujet, ici ou , par exemple).

Sommer des variables binomiales négatives

Il existe trois grandes familles de lois utilisées pour modéliser les variables de comptage,

  • la loi binomiale, 
  • la loi de Poisson, 
  • la loi binomiale négative, 

Ces trois lois sont les trois dans la famille de Sundt (ici), au sens où

En plus ces trois lois sont dans la famille exponentielle (et peuvent servir en économétrie via les GLM). Pour la loi binomiale négative, il existe des dizaines de manières de l’écrire. On peut remonter à Blaise Pascal en 1679 pour sa première utilisation. Pour reprendre la présentation de Georges Polya dans le contexte des urnes, on a des boules blanches et noires dans une urne, et on fait des tirages avec remise. Il y a une proportion p de boules noires. Le nombre de boules blanches que l’on a tiré avant de tirer la ième boule noire suit une loi binomiale négative. Mais ça ne marche que si  est un entier. On pourra envisager une écriture plus générale par la suite.
On retiendra aussi l’approche de Greenwood et Yule (1920), présentant la loi binomiale négative comme un mélange de lois de Poisson. Supposons que   suive une loi de Poisson  et que  suit une loi Gamma, de densité

Alors la loin (non conditionnelle) de X est une loi binomiale négative,

L’interprétation via les lois mélanges est assez élégante (j’avais parlé de la surdispersion de la loi de Poisson en pratique ici), et notons d’ailleurs que Greenwood et Yule avaient parlé de “accident proneness” pour décrire la loi sous cette forme
Bref, la loi binomiale négative de paramètres  et  a pour loi de probabilité

.

Remarque: un cas particulier de la loi binomiale négative est la loi ditegéométrique, obtenue quand . J’ai évoqué cette loi ici, afin d’expliquer la notion de période de retour (et l’interprétation d’un probabilité en terme de durée). Feller avait d’ailleur noté que la loi binomiale négative pouvait être vue comme une convolution de lois géométriques.
Mais je m’écarte du sujet…. J’ai reçu un mail cet après midi me demandant (je vire certaines parties pour garder la question) “la loi binomiale négative est elle additive ?”. Je comprends cette question comme la stabilité de la famille de loi par sommation: la somme de variables binomiales négatives (indépendantes) est-elle encore une loi binomiale négative ?
Notons que la loi binomiale négative est indéfiniment divisible. Mais il existe un autre résultat sur l’additivité qui nous sera encore plus utile. Liiders a montré en 1934, puis Quenouille en 1949, que la loi binomiale négative est un loi Poisson composée, lorsque les sauts sont des variables logarithmiques indépendantes (ce qui a aussi été appelé processus de Arfwedson).

où  les  sont i.i.d. de loi

La fonction génératrice est alors

où  correspond à la probabilité, et la puissance est le second paramètre.
En utilisant ce résultat, on peut en déduire que si deux variables ont le même paramètre  c’est à dire la même probabilité sous jacente, alors la loi de la somme (de deux variables indépendantes) est toujours une loi binomiale négative, dont  le paramètre sera lié à la somme de deux lois de Poisson indépendantes: ça sera donc la somme des deux autres paramètres. Aussi, si

 et 

sont deux variables indépendantes, alors

Une autre solution est de faire les calculs…(si quelqu’un veut s’y coller, je suis preneur). Enfin, si on se place dans le cas où les paramètres  sont entiers, l’interprétation en terme d’urne rend ce résultat trivial.
Remarque: on doit également obtenir la stabilité par sommation en relachant l’hypothèse d’indépendance, avec des l’indépendance conditionnellement, à condition que le paramètre de mélange (permettant d’obtenir l’indépendance conditionnelle) coïncide avec le paramètre de mélange permettant d’écrire la loi binomiale négative comme un mélange Poisonnien.
Greenwood, M., and Yule, G. U. (1920). An inquiry into the nature of frequency distributions representative of multiple happenings with particular reference to the occurrence of multiple attacks of disease or of repeated accidents, Journal of the Royal Statistical Society, Series A, 83, 255-279.
Quenouille, M. H. (1949). A relation between the logarithmic, Poisson, and negative binomial series, Biometrics, 5, 162- 164.

“too big to fail” et mesures de risque

Je suis tombé ce matin sur un billet sur “too big to fail” sur le blog de rationalité limitée (ici, et plus récemment ). On retrouve cette idée sur le blog de Paul Krugman (ici) ou sur ce qui s’est passé pour AIG, l’assureur américain, ici ou dans des journaux français, ici avec une interview de Jean-Hervé Lorenzi. Bref, c’est un sujet incroyablement à la mode… et il n’y a pas besoin d’une grande culture historique pour savoir qu’il n’existe aucune entreprise ou civilisation qui soit “too big to fail”. Sur ce dernier point, je renvoie vers le très bon livre de Jared Diamond, sur la disparition des grandes civilisations du passé.Mais on s’écarte un peu du sujet…
Le billet note fort justement que “l’idée est intéressante mais comme toujours se pose la question des critères et des méthodes de mesure”. Bref, tout ça m’a fait penser à un commentaire que m’avait fait un jour la personne qui m’a tout appris sur l’assurance, à savoir Gilles Bénéplanc1. Un grand débat sur les mesures de risque, et leur axiomatique (c’est à dire leurs propriétés souhaitables), porte sur l’axiome dit desous-additivité. Formellement, une mesure de risque est sous-additive si

autrement dit, la fusion diminue le risque. On peut aisément montrer que cet axiome peut d’ailleurs s’écrire

Il est connu que la Value-at-Risk (correspondant à un quantile de la distribution des pertes) n’est pas sous-additivite, hormis certains cas particuliers.

Et la question qui a déchiré nombre d’économistes et de mathématiciens a été “est-ce vraiment génant que la VaR ne soit pas sous-additive ?”. Gilles Bénéplanc avait fait une remarque plein de bon sens à savoir que ces questions visant à comparer

et

étaient en fait stériles car on ne devraient pas fixer les mêmes seuils (les probabilités de faire faillite si on veut interpréter simplement) pour les petites entités et l’entité agrégée. On peut accepter des ruines individuels dans un marché très concurrentiel, mais pas si une entreprise est “too big to fail”…

Il serait donc temps que la réglementation (Solvabilité II pour les compagnies d’assurance) prennent en compte cette idée de “too big to fail”, ça induirait des horizon de ruine différents pour les petites entreprises et les très gros groupes.

1 Gilles Bénéplanc a été co-encadrant de mon mémoire d’actuariat, avec François Bucchini, et comme je l’avais dit lors de la remise du prix SCOR voilà un peu plus de 2 ans, c’est à eux que je dois tout ce que je sais en assurance. Gilles avait fait une jolie thèse (encadrée par Jean Charles Rochet) sur Solvabilité des Entreprises d’Assurances et Financement des Grands Risques, et est aussi éditeur de la revue Risques.