Quand il est optimal de ne rien faire…

Tous les bloggers d’économie se sentent régulièrement obligés de parler de ce film culte… et je n’y couperais par aujourd’hui.
Ce week end, en sortant de l’exposition sur les mammouths à Rennes avec mon fils, on est passé à la fête de la science, voir si on pouvait faire quelques casses têtes à un stand tenu par des collègues de Beaulieu. Et là, un petit problème de théorie des jeux m’a été soumis, et je ne pouvais pas passer à côté de l’occasion (surtout quand on connaît la morale) d’en parler ici.
Tout le monde connaît la scène du “duel à trois” où le bon, la brute et le truand se font face… Pour les besoin du jeu, on va faire quelques hypothèses (c’est le principe de base d’un jeu: il faut une règle),

  • tout d’abord, on suppose que l’ordre dans lequel ils tirent a été fixé (et qu’ils vont s’y tenir, sans tricher (même Tuco): le bon, puis le truand, et enfin la brute,
  • ils peuvent tirer (et doivent tirer) tant qu’ils ont des balles, et tant qu’ils sont en vie, le jeu s’arrêtant quand plus personne n’a de balle, ou qu’il n’y a qu’un survivant1
  • on connaît les probabilités que chacun rate sa cible:  le bon (alias Blondin) est un peu nul (oui, il faut un peu d’imagination) car il a 1 chance sur 3 d’atteindre sa cible; le truand (Tuco) est un peu moins mauvais car il a 1 chance sur 2; enfin la brute (Sentenza) est super fort car il ne rate jamais sa cible.
  • ah oui, et les personnages ne sont pas rancuniers… ce n’est pas parce que quelqu’un nous tire dessus qu’on voudra se venger quand ça sera notre tour (et si on a la chance d’être en vie)

..La question est simple “que doit faire Clint, alias Blondin ?“. Raisonnons simplement en étudiant les scénarios,

  • s’il tire sur Tuco, il a 1 chance sur 3 de le tuer. S’il est chanceux, il le tue. Alors on retrouve avec un duel entre la brute et le bon. Et comme c’est à la brute de tirer et qu’il ne rate jamais son coup… ce n’est pas franchement malin,
  • s’il tire sur la brute, supposons qu’il le tue. Alors c’est au truand de tirer, et il a 1 chance sur 2 de gagner le dual, alors que le bon a 1 chance sur 3… Bref, les probas ne sont pas franchement en sa faveur…

mais on a regarder ce qui se passait si, par malheur, Blondin tuait quelqu’un… peut être devrait on se demander ce qu’il se passait s’il rate son coup

  • dans ce cas, peu importe qui a été visé car personne n’est rancunier. Dans ce cas, c’est à Tuco de tirer.S’il tire sur Clint, et qu’il le tue, ça sera à la brute de tirer, et alors il n’a aucune chance de s’en sortir vivant. Il est donc impératif que Tuco tire sur Sentenza.
  • si Tuco rate son coup, alors la brute a tout intérêt à tuer d’abord le truand. En effet, il risque plus de mourir s’il laisse le truand en vie que la brute. Donc, le truand meurt car Sentenza est super fort….
  • Moralité, après un tour, c’est au bon de tirer, il ne peut que tirer sur la brute et espérer réussir son coup. Sinon, c’en est fait de lui…

Bref, en ratant son premier tir, il augmente ses chances de gagner…. je laisse les personnes intéressées faire les maths, mais ce petite exemple montre clairement que, dans la vie, il est parfois optimal de ne rien faire !
1 a priori on peut supposer aussi un jeu à la roulette russe, avec une balle chacun, qui doit arriver à la même conclusion…

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

Vendredi 13 et probabilités….

Je ne sais pas si c’est la proximité avec Halloween, mais je trouvais ce vendredi 13 qui s’annonce avait quelque chose de spécial.. Et tant pis si après mes études sur les jours de naissance les jours de pleine lune, je passe définitivement pour un mystique. Disons plutôt que j’aime bien interroger les croyances de grands mères…. de la même manière que j’aime questionner les pratiques des actuaires (souvent pleines de bon sens, mais parfois un peu bancales).
Le problème avec les vendredi 13 est que l’on ne sait jamais si c’est un jour qui porte chance, ou au contraire un jour qui doit porter malheur… Dans le doute, je n’ai pas pu m’empêcher de mettre quelques trèfles à quatre feuilles, des brin de muguet et un pompon de marin.. Comme disait un très grand scientifique, “je n’y crois pas, mais il parait que ça marche même pour ceux qui n’y croient pas“. Je ne vais donc pas me priver !

Bon notons tout d’abord qu’avoir un vendredi 13 n’a rien d’exceptionnel… il y en a déjà au moins un par an. Mais ce qui est troublant c’est que le 13 d’un mois a plus de chances de tomber un vendredi qu’un autre jour. 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 (ici).
De plus, certains auteurs ont tenté des études statistiques sur le sujet. Par exemple, en 1993, le British Medical Journal avait publié un article sur le sujet (“Is Friday the 13th Bad for Your Health?” , ici) “Friday 13th is unlucky for some. The risk of hospital admission as a result of a transport accident may be increased by as much as 52 percent. Staying at home is recommended.” Bon, on va essayer de jeter un oeil….

  • Visites chez le médecin (généraliste) les vendredi 13

Les données a notre disposition (gentillement données par la Caisse d’Assurance Maladie) correspondent au nombre de consultations d’un médecin généraliste, par jour… Bon, ces données présentent un soucis technique: les données arrivent avec un (léger) délai, et les données du vendredi arrivent regroupées avec celles du samedi et du dimanche le mardi suivant. On ne peut donc pas comparer un vendredi au jeudi précédant par exemple… Pour éviter ce genre de soucis, la stratégie la plus simple consiste à comparer le nombre de consultations chez le médecin….

avec en rouge les vendredi 13, et en jaune la moyenne des vendredi 6 et 20…. Vu de loin, on consulte un peu plus les vendredi 13 que les vendredi avant et après. Si on fait un test de Student (quelle que soit l’hypothèse alternative) on rejette l’hypothèse comme quoi on conculte plus son médecin un vendredi 13…

> t.test(V13,V)
        Welch Two Sample t-test
data:  V13 and V
t = 0.6126, df = 9.052, p-value = 0.5552
alternative hypothesis: true difference in means is not equal to 0

Mais peut être est-ce par qu’une consultation chez le médecin généraliste n’est pas associé à de la “malchance” ?

  • Quid des catastrophes naturelles ?

Cherchons autre chose que l’on pourrait associer à cette prétendue “malchance“. Si on regarde les communes classées catastrophes naturelles en France, entre août 1982 et mars 2008, il y a eu 131 285 classements de communes comme touchées par une catastrophe naturelle.(et paru au Journal Officiel, ces données sont publiques). Seuls 393 classements ont été  fait un vendredi 13, ce qui pourrait sembler très faible… En fait, cette série est très ératique.

Le graphique montre l’évolution du nombre de classements les mois où il y a eu un vendredi 13. Statistiquement, on pourrait noter qu’il y a beaucoup plus de communes touchées par une catastrophe naturelle un mercredi 11 par exemple… Bref, la base semble inutilisable pour cette étude… simplement car un cinquième des classements concernent les tempêtes de décembre 1999 (ou plutôt les inondations qui ont suivies, car les tempêtes sont exclues du régime cat nat, en métropole).

  • Quid des gros sinistres de perte d’exploitation et des incendies entreprises…

J’aurais préféré avoir d’autres données d’assureurs, mais comme toujours, lorsque j’ai demandé les données, on m’a expliqué que c’étaient des données sensibles, et que je ne pourrais rien avoir ! Fort heureusement, il me reste quelques jeux de données…. en particulier une base de gros sinistres de sinistres en perte d’exploitation en France.
Cette fois, ma base va de début 1985 à fin 2000, avec en tout 2387 sinistres (cette fois l’interpréation devrait être plus simple car les données sont identifiées comme des sinistres). Au final 17 sinistres majeurs2 se sont produits un vendredi 13. Cette fois, le chiffre est significativement grand, car par exemple, seulement 4 se sont produit un jeudi 12, et 9 un samedi 14…

Là encore, si on regarde ce qui s’est passé les mois où il y avait un vendredi 13, on note qu’il y a eu beaucoup de sinistres ces jours là… Envie de faire un test ? On a en tout 2387 sinistres à faire tomber sur 5824 jours, parmi lesquels, on compte 27 vendredi 13. Bref, en supposant les sinistres indépendants, chaque sinistre à une probabilité de 0,463% de tomber un vendredi 13. Le nombre de sinistres tombant un vendredi 13 est donc supposé suivre une loi binomiale

On pourrait espérer avoir 11 sinistres qui tombent un vendredi 13 (et on en observe 17). Si on calcule le quantile à 95%, on tombe sur 17. En fait, on avait environ 1 chance sur 20 d’avoir autant de sinistres ces jours là…. ce qui correspond exactement au seuil que l’on retient pour juger de la significativité d’un test (comme je l’évoquais un peu ici). Bref, on ne peut pas vraiment conclure….
En fait, comme il s’agissait de sinistres d’entreprises, je pensais qu’il y aurais un biais, car (par définition) les vendredi 13 tombent un vendredi, et que les jours de semaines devaient avoir plus de sinistres que les jours de week end. Mais ce n’est pas le cas

Bref, on est dans la zone d’incertitude des tests statistiques….
Mais on peut aussi utiliser une base d’incendies intreprises que j’ai sous le coude depuis quelques années, avec des sinistres survenus entre 1982 et 1996. Cette base contient 9613 gros sinistres: En utilisant la même méthode que pour les sinistres de perte d’exploitation, on note que beaucoup d’incendies sont survenus un vendredi 13 dans cette base: 52. C’est un peu plus que ce que nous aurions pu espérer, à savoir 45,62. Là encore, en notant que le nombre d’incendies devrait suivre une loi binomiale, notons qu’il y avait 15,4% de chances d’atteindre 52 sinistres, ce qui n’est – cette fois-ci – pas particulièrement rare.

Par contre, sur cette base, l’effet week end semble plus prononcé, ce qui semble aller dans le sens que ces 52 sinistres, c’est loin d’être exceptionnellement elévé pour un jour de semaine.

On pourrait aussi parler de karma, ou de chose dans le genre….

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…

2 pour les plus curieux, ou les plus sceptiques, je mets ici les codes pour ceux qui veulent jouer avec les données… c’est un exercice amusant pour jouer avec la manipulation de dates sous R…

> B1=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/data-cat-nat.txt",header=TRUE)
> B2=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/data-pe.txt",header=TRUE)
> B4=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/data-incendies.txt",header=TRUE)
> B3=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/date-v13.txt",header=TRUE) 
> B2$DATEPE=as.character(B2$DATEPE)
> x2=as.Date(B2$DATEPE,"%d/%m/%Y")   
> B3$date=as.character(B2$date)
> x3=as.Date(B3$date,"%d/%m/%Y")
> p2=sum((x3>=min(x2))&(x3<=max(x2)))/as.numeric(max(x2)-min(x2))
> n2=length(x2);  n2*p2
> qbinom(.05,n2,p2); qbinom(.95,n2,p2)
> pbinom(sum(x2%in%x3)  ,n2,p2)
> NOM=weekdays(x3[1]+1:7-5)
> S=rep(NA,7)
> for(i in 1:7){S[i]=sum(weekdays(x2)==NOM[i])}
> barplot(S,names.arg=NOM,col="light blue")


voilà pour les sinistres en perte d’exploitation… le reste n’est qu’une variation sur le même thème…

Jeu de mains, jeu de malins ?

Ce week-end – merci à Serge qui m’a raconté cette histoire – beaucoup de mes illusions de jeunesse se sont envolées…. Mais gardons un peu de suspens pour que mes lecteurs lisent le billet en entier….

  • Pierre, ciseaux ou feuille

Le jeu “pierre, ciseaux, feuille” semble être né en Chine (soit sous le joli nom 五雜俎, alias wuzazu, datant de la fin de période Ming, soit 手勢令 i.e.shoushiling) ou au Japon (soit sous la forme 一二三 , i.e. Hi-fu-mu, c’est à dire un, deux, trois,soit じゃんけんぽん ou plus simplement じゃんけん, alias janken  , lui même basé sur deux anciens jeux, 数拳 – su ken – et 三すくみ拳 -saan sukumi ken). Pour rappel, la pierre (ぐう) bat les ciseaux (ちょき) qui bat à la feuille (ぱあ) qui à son tour bat la pierre. Bref, ce jeux vient d’orient !
Pour preuve, en avril 2005, Takashi Hashiyama, alors président de Maspro Denkoh Corporation, une entreprise d’électronique située à Nagoya, au Japon, n’arrivait pas à savoir chez qui (entre Christie’s et Sotheby’s) il allait vendre sa collection d’art (estimée à plus de 20 millions de dollars). N’arrivant pas à choisir, le choix s’est fait par “pierre, ciseaux, feuille” (ici).

  • Le théorème min-max

Pour comprendre ce jeu, il faut revenir un peu sur des résultats classiques de théorie des jeux. La matrice de gain pour le joueur X est ici

GAINS
X joue « pierre »
X joue « feuille »
X joue « ciseaux »
Y joue « pierre »
0
1
-1
Y joue « feuille »
-1
0
1
Y joue « ciseaux »
1
-1
0

et de manière symétrique, on peut définir la matrice de gain pour le joueur Y. Notons que l’on suppose ici un “jeu à somme nulle” (i.e. ce que perd X,  Y le gagne, et réciproquement). Cette seconde matrice de gains est alors l’opposé de celle-ci.
La stratégie de base du jeu a été énoncée par Italo Calvino sous la forme suivante, “tu sais que ce que tu peux espérer de mieux est d’éviter le pire”.
Intuitivement, supposons que X souhaite jouer “pierre” on peut penser que Y reconstitue les anticipations de X, et donc joue “feuille”, mais à son tour X doit modifier ses anticipations et donc jouer “ciseaux”, plutôt que “feuille”, etc. Bref, rien ne se stabilise et aucune stratégie pure ne s’impose.
Mais s’il n’y a pas de stratégie pure, tentons des stratégies mixtes, i.e. choisissons intelligemment une mesure de probabilités sur l’ensemble des stratégies pures. En pratique, on jouera alors aléatoirement, en pondérant chaque tirage suivant la mesure de probas choisie. Une stratégie mixte (pour X) est alors caractérisée par un vecteur du simplexe,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/minimax-1.png

La  stratégie mixte

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/minimax-2.png

consiste à décider au  hasard quel coup on jouera avec équiprobabilité.
Formellement, notons que le tirage équiprobable – par les deux joueurs – est un équilibre de Nash. En effet, si A désigne la matrice de gain pour X (on est sur des jeux à somme nulle), John von Neumann a montré ce que l’on appelle aujourd’hui le théorème du minimax (en fait, Emile Borel avait également obtenu ce résultat) en 1926, sous la forme suivante: il y a un équilibre de Nash associé à A si et seulement si

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/minimax-3.png

L’intuition est assez évidente (mais la preuve de ce théorème du point fixe est un peu technique): https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/minimax-4.png est un équilibre de Nash si aucun des jours n’a intérêt de dévier dans sa stratégie, i.e.

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/minimax-5.png

et

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/minimax-6.png

On recherche alors un point selle dans la matrice des gains.
Dans le cas de notre jeu de “pierre, ciseaux, feuille“, notons que si

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/minimax-7.png

alors  https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/minimax-8.png, et on note que  https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/minimax-9.png ainsi que https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/minimax-10.png. Aussi, pour toute  stratégie mixte,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/minimax-11.png

et de même

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/minimax-12.png

Bref, cette stratégie est une équilibre de Nash.

  • Pierre, ciseaux, feuille, ou comment trouver une stratégie

Bref, comme on vient de le noter, la meilleure stratégie à adopter quand on jour à “pierre, caillou, feuille” c’est de ne pas avoir de stratégie… ce qui revient à dire que l’issue du jeu est purement aléatoire. Pour Kenneth Binmore (ici) la méthode du minimax consiste à dire que face à un “bon” joueur, il n’y a rien de mieux à faire de que choisir la stratégie mixte du théorème de von Neumann.
Bon, si on réfléchit deux minutes, on sait qu’on ne joue pas une fois à ce jeux, mais que l’on répète le jeux quelques fois (avant de se lasser). C’est ce qu’a décrit Julia Robinson en 1950. L’idée était d’avoir une approche que l’on peut qualifier de bayésienne, consistant à utiliser pour la stratégie de l’adversaire les probabilités observées.

  • Un champion du monde de “pierre, ciseaux, feuille” ?

Maintenant que j’ai fini avec les rappels de théorie, revenons à la discussion de samedi soir…. J’ai appris qu’il existait un champion du monde de “pierre, ciseaux, feuille” ! Serge prétendait – samedi soir – qu’un joueur avait été vainqueur deux années de suite (ce qui m’a vraiment fait douté sur le côté aléatoire du jeu… malheureusement pour mon histoire, ce n’est qu’une légende, car les champions du monde entre 2002 et 2008 ont été Peter Lovering, Rob Krueger, Lee Rammage, Andrew Bergel, Bob Cooper, Andrea Farina et Monica Martinez).
Comment peut-il y avoir un champion du monde à un jeu qui repose exclusivement sur du hasard ? C’est comme si on inventait un champion du monde de lancer de pièces (ou de dés), qui serait récompensé pour sa régularité1 !
En creusant un peu, j’ai aussi découvert qu’il existait un livre qui apprend à gagner au jeu (je ne l’ai encore lu, mais promis, si la BU le commande, je le dévorerais et ferais un compte rendu détaillé ici même).