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

Pour avoir un cycle de semaines entières, il faut compter 400 ans, soit 20871 semaines (ce que l’on appelle un cycle grégorien). Mais en 400 ans, on a 4800 mois, soit 4800 “13”… qui n’est pas divisible par 7, ce qui tend à laisser penser que certains jours de la semaine tomberont (un peu) plus souvent que les autres… Et Brown a montré que c’était le cas (dès 1933, dans le American Mathemical Monthly1). Pour les amateurs, Jean-Luc Nothias avait fait un article en mars dernier sur le sujet (







1 B. H. Brown, “Solution to Problem E36”, American Mathematical Monthly, vol. 40, issue 10 (1933), p. 607 Pour faire simple, un cycle grégorien correspond à 400 années, soit 136097 jours, ou encore 20 871 semaines. Si on regarde la distribution des jours on voit que samedi et jeudi tombent 684 fois, lund et mardi 685, mercredi et dimanche 687, et vendredi 688. Bon, d’accord, je pinaille…




est un équilibre de Nash si aucun des jours n’a intérêt de dévier dans sa stratégie, i.e.


, et on note que
ainsi que
. Aussi, pour toute stratégie mixte,


,
. Cette suite est sera une
, et pour tout



et de matrice de transition P, pour tout
,

,







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 (
Lors du dîner de gala des
le montant de la prime associée à la kième affaire (si on utilisait le joker). Il faut alors arbitrer, chaque jour k, entre
la valeur du joker à la date k. Alors


(le dernier jour, si on a le joker, on l’utilise). Autrement dit, on devrait y arriver par induction backward… Et la résolution dépend de la loi des montants des affaires.
,











en 
le facteur d’actualisation entre les dates
et 
est la filtration naturelle, et
est une probabilité risque neutre, sous laquelle la valeur actualisée de l’actif est une martingale, i.e.




































et 
où
.






et 

).






et vérifiant le système
et


















et
désignent les quantités des deux biens qui sont disponibles, respectivement. Classiquement, on résout ce programme d’optimisation à l’aide du lagrangien

















