The course on risk and insurance in Luminy, starts at 10.30 on Friday (here) instead of Thursday (I switched with Patrice Bertail). The slides can be found here,
Then, it will be time to leave Marseille
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).
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,
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…
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) } }