Category Archives: reserving

Provisionnement et règle de trois

Longtemps, dans les formations d’actuariat, le provisionnement était un sujet traité en 2 heures, car malgré son importance stratégique dans les résultats d’une entreprise d’assurance, il ne s’agissait que d’une “règle de trois“.

  • De l’importance des provisions en assurance

D’un point de vue comptable, les provisions pour sinistres à payer sont une quantification des engagements pris par l’assureur envers ses assurés. En reprenant des vieux tableaux dans l’Argus, avec une compagnie d’assurance spécialisée en grands risques (AXA Corporate Solutions), une mutuelle (MAIF), et une autre société anonyme (Groupama), on obtient

Pour les grands risques, on notera que le montant de provision (plus de 5 milliards) est plus de 3 fois supérieur au montant des primes encaissées (1,5 milliards). Aussi, si l’on pense qu’il faut augmenter de 5% le montant des provisions (ce qui pourrait paraître acceptable à un statisticien) il faudrait augmenter le chiffre d’affaire de 15%.

  • Histoire de la “règle de trois

Le Petit Robert prétendait que l’expression « règle de trois » remontrait à 1538, même si le résultat sous-jacent est tout simplement le théorème de Thalès… Les manuels d’arithmétiques expliquaient ainsi que “la règle de trois est la méthode servant à résoudre les problèmes dans lesquels, étant données des grandeurs pouvant former les trois termes d’une proportion, on demande de trouver le quatrième terme“. On peut retrouver beaucoup plus d’infos sur le site images des mathématiques (ici) ou (), où on trouve d’ailleurs cette illustration,

La “règle de trois” est aussi le premier souvenir de raisonnement mathématique exposée par Paul Lévy dans sa biographie (“quelques aspects de la pensée d’un mathématicien“), lorsqu’il raconte qu’à 7 ans, il tenait la main d’une fille de 12 ans, il était surpris de lever la main. En effet, certes la jeune fille était plus grande que lui, mais son bras devait également être plus grand. Il compris alors le principe des règles de proportionnalité.
Bref, comme toujours, les actuaires ont inventé un nouveau mot pour la “règle de trois” (je renvoie à un vieux papier de Claude Bébéar sur la novlangue utilisée par les actuaires, ici, mais je reviendrais longuement sur le jargon une autre fois).

  • Les méthodes de provisionnement: Chain Ladder

J’avais évoqué il y a quelques temps les travaux d’Eugène Astesan (ici), précurseur de la méthode Chain Ladder, même si l’estimateur qu’il avait alors proposé n’était pas celui de variance minimale dans la classe des estimateurs. Les premiers travaux sur les modèles stochastiques de provisionnement ont chercher à justifier la règle de trois, par exemple via le papier de Thomas Mack “which stochastic model is underlying the chain ladder method” (ici).
Pour venir deux minutes sur la méthode “chain ladder“, par exemple à l’aide du triangle suivant, on peut estimer les montants à payer (dans la partie inférieure) via une simple règle de trois. Considérons le triangle suivant, tiré du livre de Partar, Lecoeur, Nessi,  Nisipasu et Reiz

Pour la seconde colonne, on fait une simple règle de trois, en calculant les sommes sur chacune des colonnes,

Aussi, ici,

i.e.

117,945 = 56,762 x 440,438 / 211,961

On peut ensuite passer à la troisième colonne, puis la quatrième, etc. La règle de trois donne ici

Voici pour la méthode dite “chain ladder“. Et Thomas Mack a été le premier à essayer de quantifier l’erreur de prédiction commise en utilisant cet estimateur. En particulier, les travaux théoriques des années 90 ont essayer d’étudier l’incertitude associée à l’estimation de la charge ultime, calculée aujourd’hui. Comme le notait un article paru dans l’Argus, les compagnies d’assurance qui utilisaient ces techniques étaient encore rares il y a 5 ans (et celles qui le faisaient cherchaient également à le faire savoir, ici). Mais les derniers documents (ici ou ) partent sur un objectif radicalement différent, que j’exposerais d’ici quelques jours, en revenant sur la novlangue des récents papiers sur le provisionnement..

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