Ma vie et les bases de données

Ayant remarqué que les bloggeurs avaient souvent, un jour, envie de parler d’eux dans un billet, je vais faire pareil (oui, je sais Igor prétend que je ne fais déjà que ca ici). Je vais donc parler un peu de ma vie en tant que bloggeur.

Ce blog a été créé après un transfert du blog que je tenais à Rennes 1, qui avait été lancé voilà bientôt un peu plus de 2 ans, et je dois avouer (maintenant il y prescription) que j’ai honteusement volé le nom à un extrait de conversation que j’avais avec Francis Kramarz à l’X (reprenant le nom inventé par Stephen Dubner et Steven Levitt, i.e.freakonomics). L’idée du blog était de parler de choses qui m’amusent dans ma vie d’enseignant-chercheur.

L’économie mathématique n’est pas très sexy, et n’intéresse que peu de monde, donc j’évite d’en mettre trop… par contre, j’aime aussi beaucoup les maths appliqués, les statistiques, et l’économétrie. Et l’avantage quand on travaille sur la modélisation, y compris en tant qu’enseignant, c’est qu’on peut illustrer un peu comme on veut. Personnellement je comprends mieux la notion d’optimum de Pareto quand on me parle de se servir de l’eau à la cantine (comme ici), plutôt que de me donner une définition abstraite.

Et il s’avère que ça a plutôt tendance à réveiller les élèves quand on dit que les processus VAR permettent de suivre les mensurations des demoiselles de playboy (comme ici), en tous les cas plus que la modélisation des taux d’intérêt dans la zone euro. Les processus aléatoires, comme  je m’efforce à le dire (ici ou), c’est beaucoup plus marrant appliqué à des problèmes d’alcooliques qu’appliqué à la valorisation d’options financières. Et les chaînes de Markov appliquées à l’analyse de la pluie en Bretagne (comme ici), ça intéresse tous les élèves qui viennent à la fac en vélo, ou sinon je peux aussi parler de consommation de papier toilette dans les toilettes publique (comme ).

Bref, peu importe le sujet, autant que ce soit amusant, l’idée étant d’illustrer une méthode, une technique, un concept abstrait…

Donc pour amuser mes étudiants (et mes lecteurs) je traîne beaucoup sur internet à la recherche de bases de données, afin d’illustrer des méthodes d’estimation. Alors pour ceux qui n’en ont jamais manipulé, les données, c’est un truc assez incroyable…. Vu de loin, ça ne ressemble pas à grand chose, genre ça peut ressembler à ça,

Mais il y a quelques années, on a vu arriver le datamining, et “les gens” ont commencé à retenir une idée un peu générale du genre “vos données valent de l’or”, et depuis, certains refusent de communiquer leurs données. Je ne parle pas de chercheurs qui ont collecté des données, et qui n’ont pas fini de publier dessus (je comprends qu’il puisse exister une relative compétition entre chercheurs), mais je parle de base qui existent…. mais que personne n’exploite. En fait “les gens” ont oublié que s’ils avaient de l’or, il faut peut être des mineurs pour aller le chercher, puis des joaillers pour le retravailler, et enfin des bijoutiers pour rendre ça présentable. Et c’est souvent agaçant de faire face à un refus, de voir ces belles données qui existent, mais qu’on ne peut avoir….

Bref, aujourd’hui je vais faire un billet sur les billets que j’aurais rêvé de faire, mais que je ne fais pas, faute de données….

  • L’argument (magique) de la CNIL pour les bases informatiques

Parmi les arguments qui ressortent sans cesse quand on demande des données, c’est “il faut avoir l’accord de la CNIL”. Bon, c’est l’argument de mauvaise foi par excellence…. On m’a sorti cet argument dans un hôpital, ou je voulais avoir des données sur les accouchements, connaitre la durée de la grossesse en fonction de caractéristiques de la mère (et du nombre d’enfants déjà eu, par exemple). Un hôpital m’a aussi dit “on peut vous ouvrir nos armoires, et vous pourrez noter vous même car ca n’est pas informatisé” (oui oui, j’ai eu ca… et j’ai pu vérifier avec la naissance de ma troisième: tout était sur papier… oups, j’ai dit de quel hôpital il s’agissait).

  • On veut la conclusion de l’étude avant de fournir les études

Toujours sur mes études sur la durée des grossesses, j’ai appris que l’Inserm constituait une grosse base de données. J’ai demande à y avoir accès, et fort gentiment, on m’a demande de déposer un dossier de recherche. J’ai donc mis l’objet de l’étude, le nom des étudiants qui auraient travaille sur la base, une biblio pour expliquer que je ne partais pas à l’aveugle, et la liste des variables qui m’auraient intéressées. Et la, par téléphone, on m’a dit que c’était très bien, mais “on a besoin d’avoir vos conclusions avant de vous donner accès aux données”. La trouve la démarche (supposée scientifique) surprenante. L’argument m’a un peu dépassé, et je n’ai pas donne suite.

  • L’invention d’arguments (farfelus, parfois)

Toujours sur mon étude sur les naissances, j’ai voulu avoir l’information à la source, c’est à dire dans les CAF. Elles connaissent le terme prévu, le terme réel entre autres, elles disposent de tout plein de variables explicatives. Bref, lasource de données idéale. J’ai donc contacté une vingtaine de CAF (car les CAF sont départementales), et j’ai eu une vingtaine de réponses différentes allant du “on n’a pas ce genre de données”, “on n’as pas le temps”,  “c’est la Caisse Nationale qui a ces données”, “il faut l’accord de la Caisse Nationale “… Bref, assez souvent, on botte en touche… Un peu comme “mais vous savez, on ne vend pas de fleurs” que m’a sorti Interflora (je voulais savoir si on commandait plus de fleurs a la Saint Valentin ou a la fête des mères, en vain).

  • Nous on fait dans le qualitatif, pas le quantitatif

Il y a quelques années, j’avais vu passer une étude sur l’évolution des mensurations de miss America (et les données sont en ligne ici). Trouvant l’analyse intéressante, j’ai contacte le comité Miss France pour avoir les mensurations des miss (et si possible des prétendantes), et Madame de Fontenay (oui oui, elle m’a écrit) m’a répondu “nous n’avons pas ces informations, les miss ne sont pas choisies pour leurs mensurations mais pour leur personnalité” (je n’invente presque pas – la seule invention vient de la retranscription car j’ai perdu le mail avec tous mes transferts de boites, mais c’est ce qui était dit en substance). Finalement j’ai ressorti mes vieux playboy, ce qui a permis d’avoir encore plus d’observations d’ailleurs….

  • Ah oui, je vois très bien, mais demandez au service ETDDHRC qui s’en occupe

L’argument magique du “c’est pas moi c’est mon frère” pour reprendre La Fontaine (ou ma fille, je ne sais plus trop). Je l’ai eu sur des demandes anciennes, faites à plusieurs reprises (et liées à des sujets qui ont beaucoup fait parler dans la blogosphere)

  • il y a 4 ou 5 ans, j’ai eu accès à un très gros fichier d’expertises de véhicules, consécutives à des accidents de la route, plusieurs centaines de milliers d’accidents analysés… on avait voulu travailler avec des étudiants sur l’impact des radars automatiques, mais surtout sur la fraude à l’assurance… Sauf qu’avoir les accidents c’est bien, mais il fallait normer par le nombre de véhicules en circulation… j’ai donc demande à plusieurs reprises un accès au fichier des cartes grises pour savoir combien de véhicules étant en circulation (marque, modèle)… tous les services se sont renvoyés la balle… Et je n’ai jamais eu la base…
  • j’ai eu la même discussion il y a quelques mois sur les nouvelles plaques… le Ministère de l’Intérieur affiche (ici) les départements les plus demande (car on a maintenant le choix du numéro que l’on pose sur la plaque), mais on ne peut pas croiser le département de résidence et le département de la plaque…. Étude qui m’aurait amusé car j’ai un joli modèle à tester derrière…
  • sinon à la même période, je souhaite analyser l’impact de la taille des classes sur la réussite scolaire (j’ai fait 4 billets sur le sujets, ici ou , mais sur des données israéliennes): j’ai demande au Ministère de l’Éducation Nationale à avoir accès aux résultats individuels (mais anonymes) d’échantillons du test national de CE1. En particulier, j’étais aussi intéresse par la réussite des enfants nés en janvier par rapport à ceux nés en décembre… Tous les services se renvoyaient la balle, et je n’ai jamais eu les données.

J’oublie de préciser, mais a chaque fois, je précise qu’il s’agit d’une étude académique, et que peux m’engager a ne pas diffuser les données.

  • Le service juridique veut que vous signez 10 pages de “disclosure”

Alors ça c’est le pire… c’est arrivé l’autre jour avec un copain actuaire. Je l’avais sollicité pour quelques données, des triangles de liquidation comme on dit. Bon, ce sont des données un peu sensibles, mais je ne demandais pas les montants de provisions constitues, juste des paiements. J’ai un “c’est bon de notre cote, par contre il faudra que tu signes le document joint”, à savoir un document juridique de 10 pages, qui devaient être signées par les présidents de l’Université de Rennes 1 et de l’Université de Montréal (car c’était pour donner un projet à mes étudiants). Ne voulant pas prendre de risques, j’ai laissé tombé….

  • Allez vous servir sur le site (en faisant 5346  requêtes)

Cet argument m’a été sorti par le Ministère de l’Écologie:. Je voulais connaître la liste de toutes les communes touchées par une catastrophe naturelle (et qui ont bénéficié d’un arrêté), “c’est sur le site internet, vous savez, on est moderne nous”. Cette information est publique (parue au Journal Officiel) mais je voulais une base. J’ai fait un code l’autre jour pour aller faire plusieurs milliers de requetes (sur les matchs de tennis, mais ca prend du temps… et ca n’est pas toujours simple car il faut aller decoder des fichiers html) heureusement, il y a toujours des copains qui ont ces donnees et qui sont prêt à envoyer un fichier xls dans l’heure).

  • On veut bien donner des bases, mais à des journalistes, pas à des universitaires

Pour la première fois, on évoquait la discrétion, mais dans l’autre sens. Classiquement c’est “on vous donne des chiffres, mais faut pas le dire”, mais là c’était “on veut bien les donner, mais faut que ça se sache…”. Heureusement, Twitter est là, et on finira par trouver un point d’entrée…. à suivre donc…
Fort heureusement, certains ont compris leur intérêt, et n’hésite pas à m’envoyer des bases si je leur demande. Ca fait une étude gratuite pour eux, et moi je me serais amusé un peu (ou fait travaillé mes étudiants). Mais ce n’est pas la majorité. En attendant, je suis devenu voleur de données. Beaucoup de données traînent sur internet, et en bricolant un peu on peut les récupérer… L’autre jour un copain (que j’avais sollicité pour aller chercher des bases en faisant automatiquement des requêtes sur un site) me demandait “elles sont à qui ces données ?”. Et je dois avouer que depuis je doute…. Moralité, je suis un peu surpris de l’étymologie du mot données, car rares sont ceux qui en donnent vraiment… et c’est dommage.

Non, les probabilités ne sont pas nos amies…

… et non, nous ne vivons pas le monde merveilleux des Bisounours: la vie n’est pas toujours juste. Les probabilités sont un objet assez curieux à manipuler, pas forcément très intuitif. J’avais évoqué ici et enfin  la notion d’espace probabilisé http://freakonometrics.blog.free.fr/public/perso/espace-probabilise.png… En fait, généralement, les calculs de probabilités sont assez troublants et paradoxaux. On va considérer ici des gens devant prendre des situations en présence d’incertitude, et plus particulièrement deux situations: la difficulté de prendre en compte de l’information (et donc de calculer des probabilités conditionnelles) et le cas où les gens doivent agir en prenant des décisions de manière aléatoire (et donc fixent en même les probabilités).

  • Les probabilités conditionnelles

Avant hier, j’avais mis en ligne (ici) un billet sur le paradoxe dit de Monty Hall, où le candidat à un jeu télévisé doit choisir une porte parmi trois (afin de gagner une voiture), une fois son choix effectué, le présentateur ouvre une des deux autres portes (une qui n’a pas de lot derrière) et demande au candidat s’il souhaite revenir sur son choix. Et la réponse est qu’effectivement, le candidat a intérêt à ouvrir la porte qui est restée ouverte (et que l’animateur n’a pas choisie). L’idée du raisonnement est de formaliser proprement cette information à l’aide de lois conditionnelles.

Cette importance de l’information se retrouve par exemple dans l’histoire des cocus de Bagdad ou des yeux bleus (ici). Sinon, sur les paradoxes de probabilité, je peux aussi renvoyer vers des histoires de des lancers de pièces (ici), mais aussi sur le paradoxe des enveloppes () qui pourrait faire penser à l’histoire de Monty Hall. Mais mon préféré sur les lois conditionnelles reste le paradoxe de Simpson (ici, repris sous une forme plus géométrique )…

  • Les probabilités comme paramètres de décision

Dans tous ces exemples, les probabilités sont connues, et il s’agit simplement d’une mise à jour afin de prendre en compte une nouvelle information. Mais il existe des cas encore plus vicieux, où les probabilités ne peuvent pas être considérées comme exogènes, ou objectives, car elles sont déterminées par les joueurs. Considérons une jeu simple, un peu dans l’esprit du pierre-feuille-ciseau (que j’avais déjà évoqué ici, en rappelant qu’il était fondamental de prendre des décisions de manière aléatoire). Dans le jeu d’aujourd’hui, on a toujours deux joueurs, disons moi (P) et ma fille (F). On cache nos mains dans le dos, et on les sort au même moment, en montant un ou deux doigts,

  • si le nombre de doigts sortis est identique (paires (1,1) ou (2,2)), alors je gagne le nombre de doigts sortis, i.e. 2 ou 4
  • si le nombre de doigts sortis est différent (paires (1,2) ou (2,1)), alors ma fille gagne le nombre de doigts sortis, i.e. 3

La matrice des paiements est alors la suivante

P\F 1 2
1 2 -3
2 -3 4

Si on suppose que l’on a une chance sur 4 pour chacune des paires, le jeu est alors a somme nulle, puisque
http://freakonometrics.blog.free.fr/public/maths/jeuPF-80.png
Sauf que les probabilités sont fixées par les joueurs, en l’occurrence moi et ma fille. Supposons que je sorte 1 doigt avec probabilité http://freakonometrics.blog.free.fr/public/maths/jeuPF1.png, et 2 avec probabilité http://freakonometrics.blog.free.fr/public/maths/jeuPF2.png(avec http://freakonometrics.blog.free.fr/public/maths/jeuPF-81.png). Ma fille sort 1 doigt avec probabilité http://freakonometrics.blog.free.fr/public/maths/jeuPF3.png, et 2 avec probabilitéhttp://freakonometrics.blog.free.fr/public/maths/jeuPF4.png (avec http://freakonometrics.blog.free.fr/public/maths/jeuPF5.png). La valeur espérée du jeu (de mon point de vue) est
http://freakonometrics.blog.free.fr/public/maths/jeuPF.png
Ma fille cherche http://freakonometrics.blog.free.fr/public/maths/jeuPF15.png qui va minimiser cette valeur, alors que moi, je cherche http://freakonometrics.blog.free.fr/public/maths/jeuPF16.png qui la maximise. Bref, un problème classique de problème minmax de théorie des jeu. L’équilibre (s’il existe) sera http://freakonometrics.blog.free.fr/public/maths/jeuPF17.png où http://freakonometrics.blog.free.fr/public/maths/jeuPF-18.pnget http://freakonometrics.blog.free.fr/public/maths/jeuPF-19.png solution de
http://freakonometrics.blog.free.fr/public/maths/jeuPF-20.png
A p donné, ma fille cherche à trouver
http://freakonometrics.blog.free.fr/public/maths/jeuPF-21.png
soit
http://freakonometrics.blog.free.fr/public/maths/jeuPF-23.png
soit encore
http://freakonometrics.blog.free.fr/public/maths/jeuPF-25.png
en notant que http://freakonometrics.blog.free.fr/public/maths/jeuPF-81.png, on peut réécrire http://freakonometrics.blog.free.fr/public/maths/jeuPF-26.png. Autrement dit, si http://freakonometrics.blog.free.fr/public/maths/jeuPF-27.png, ma fille a toujours intérêt à sortir un doigt (car http://freakonometrics.blog.free.fr/public/maths/jeuPF-82.png) et ce n’est pas un équilibre… pareillement, si http://freakonometrics.blog.free.fr/public/maths/jeuPF-28.png, ma fille a toujours intérêt à sortir deux doigts, et là encore ce n’est pas un équilibre. Autrement dit, le seul équilibre possible est obtenu lorsque http://freakonometrics.blog.free.fr/public/maths/jeuPF-29.png, i.e l’équilibre correspond à http://freakonometrics.blog.free.fr/public/maths/jeuPF-30.png. De manière dual, pour ma part, je cherchais à trouver
http://freakonometrics.blog.free.fr/public/maths/jeuPF-31.png
i.e.
http://freakonometrics.blog.free.fr/public/maths/jeuPF-32.png
soit
http://freakonometrics.blog.free.fr/public/maths/jeuPF-25.png
qui conduit à un équilibre seulement si http://freakonometrics.blog.free.fr/public/maths/jeuPF-33.png. Autrement dit, afin que le jeu soit à l’équilibre, il faut que nous jouions la stratégie suivante: sortir un doigt avec probabilité http://freakonometrics.blog.free.fr/public/maths/jeuPF-40.png (et pas http://freakonometrics.blog.free.fr/public/maths/jeuPF60.png) et deux avec probabilité http://freakonometrics.blog.free.fr/public/maths/jeuPF-41.png. Et la valeur espérée du jeu est alors
http://freakonometrics.blog.free.fr/public/maths/jeuPF-50.png
Moralité, ma fille tu as intérêt à travailler tes maths si tu veux battre un jour ton père ce n’est pas vraiment un jeu “juste“… à condition bien sûr que moi et ma fille soyons ensuite capable de mettre en œuvre cette stratégie… Car si je regarde avec mes yeux de statisticiens, je pense qu’il est assez dur de faire la différence entre http://freakonometrics.blog.free.fr/public/maths/jeuPF-40.png et http://freakonometrics.blog.free.fr/public/maths/jeuPF60.png … La preuve sur cette petite simulation,

> N=100; pv=rep(NA,100000)
> for(s in 1:100000){
+ X=sample(1:2,size=N,prob=c(5/12,7/12),replace=TRUE)
+ P=prop.test(table(X)[1],n=N,p=1/2)
+ pv[s]=P$p.value}
> mean(pv>0.05)
[1] 0.66776

autrement dit, en supposant que ma fille et moi arrivions a jouer 100 fois à ce jeu, dans 2/3 des cas je n’arrive pas à faire la différence entre une probabilité valant http://freakonometrics.blog.free.fr/public/maths/jeuPF60.png et une probabilité valant http://freakonometrics.blog.free.fr/public/maths/jeuPF-40.png.

Tennis and risk management

As mentioned already here, while we were going to Québec City for the workshop, we had interesting discussions in the car, and Maciej mentioned an article recently published in The Actuary,

Hence, I wanted to discuss (extremely) rare event probabilities in tennis. The story is simple: in June 2010, at Wimbledon, Nicolas Mahut and John Isner have played the longest match ever. 980 points, 11 But first of all, we need a dataset. Thanks to Duncan Murdoch, I have been able to run a short code to build up a dataset:

CITIES=c("berlin","madrid","paris","rolandgarros","wimbledon","sydney",
"beijing","shanghai","singapore","tokyo","melbourne","melbourne-indoor")
YEARS=1970:2009
BASE0=data.frame(YEAR=NA,TRNMT=NA,LENGTH=NA,SETS=NA)
for(i in 1:length(CITIES)){
for(j in 1:length(YEARS)){
city=CITIES[i]
year=YEARS[j]
localization = paste("http://www.resultsfromtennis.com/",
year,"/atp/",city,".html",sep="")
essai = try(readLines(localization), silent=TRUE)
ERROR404=FALSE
if(inherits(essai, "try-error")){ERROR404=TRUE}
if(ERROR404==FALSE){
B=scan(localization,"character")
SETS=NA
LENGTH=NA
if(length(B)>270){
I=(substr(B,1,10)=="class=rez>")
sum(I)
X0=B[I]
X3=as.numeric(substr(X0,11,13))
X2=as.numeric(substr(X0,11,12))
X1=as.numeric(substr(X0,11,11))
X0=X3
X0[is.na(X3)==TRUE]=X2[is.na(X3)==TRUE]
X0[is.na(X2)==TRUE]=X1[is.na(X2)==TRUE]
JL=c(which(substr(B,1,9)=="class=nl>"),length(B))
IL=which(substr(B,1,10)=="class=rez>")
IC=cut(IL,JL)
base=data.frame(IC,X0)
LENGTH=as.numeric(tapply(X0,IC,sum))
SETS=as.numeric(tapply(X0,IC,length))/2}
BASE=data.frame(YEAR=year,TRNMT=city,LENGTH,SETS)
BASE0=rbind(BASE0,BASE)}}}
write.table(BASE0,"BASE-TENNIS-TOTAL.txt")

Here I consider only tournaments where players have to win 3 sets (and actually more tournaments than those in the code above), and I have something like a bit more than 72,000 matches,

> I=is.na(TENNIS$LENGTH)==FALSE
> BT=TENNIS[I,]
> nrow(BT)
[1] 72754
> maxr=function(x){max(x,na.rm=TRUE)}
> T=paste(BT$TRNMT,BT$YEAR)
> DUREE=tapply(BT$SETS,T,maxr)
> LISTE=names(DUREE[DUREE>3])
> BT=BT[T%in%LISTE,]

so, if we look briefly at matches over 35 years, we have the following boxplot (one boxplot per year),

The red line being the epic Isner-Mahut match in June 2010 (4-6, 6-3, 7-6, 6-7, 70-68, i.e. 183 games, here for the score card).

If we study theory (e.g. from Paul Newton and Kamran Aslam), a lot of results can be obtained for the expected value of the number of games, but if we want to study extremely rare events, we should generate Markov chains (with a lot of generation since the probability should be extremely small). But how many ? Consider below matches with more than 50 games,

The tail plot (over 50), i.e. the log-log Pareto plot indicates that it will be difficult to study tails,

and similarly with the Hill plot (assuming that tails are Pareto type….)

Anyway, if we want to study tails, we should consider a threshold high enough. For instance, with a threshold at 68 (we keep only 24 match), we have

> seuil=68+0.25
> GPD1=gpd(X,seuil,method = "ml")
> GPD2=gpd(X,seuil,method = "pwm")
>
> xi=GPD1$par.ests[1]
> mu=seuil
> beta=GPD1$par.ests[2]
> x=180
> P=exp((-1/xi)*log(1 + (xi * (x - mu))/beta))
> as.numeric((1-GPD1$p.less.thresh)*P)
[1] 5.621281e-09
>
> xi=GPD2$par.ests[1]
> mu=seuil
> beta=GPD2$par.ests[2]
> x=180
> P=exp((-1/xi)*log(1 + (xi * (x - mu))/beta))
> as.numeric((1-GPD2$p.less.thresh)*P)
[1] 3.027095e-09

I.e. the probability that one match last more than 183 games is 1 chance over a billion… With, say, 2500 match per year, that gives us a return period of 400 years. So yes, we might say that this way a rare event… So perhaps, generating several billions of chains, it should be possible to get a more precise estimation of the probability to play 183 games in a single match…

I really need to find hot (and sexy) topics

50 days ago (here), I was supposed to be very optimistic about the probability that I could reach a million viewed pages on that blog (over a bit more than two years). Unfortunately, the wind has changed and today, the probability is quite low…

 base=read.table("millionb.csv",sep=";",header=TRUE)
X1=cumsum(base$nombre)
base=read.table("million2b.csv",sep=";",header=TRUE)
X2=cumsum(base$nombre)
X=X1+X2
 D=as.Date(as.character(base$date),"%m/%d/%Y")
kt=which(D==as.Date("01/06/2010","%d/%m/%Y"))
D0=as.Date("08/11/2008","%d/%m/%Y")
D=D0+1:length(X1)
P=rep(NA,(length(X)-kt)+1)
for(h in 0:(length(X)-kt)){
model  <- arima(X[1:(kt+h)],c(7 1,7),method="CSS") 
 forecast <- predict(model,200)
u=max(D[1:kt+h])+1:300
k=which(u==as.Date("01/01/2011","%d/%m/%Y"))
(P[h+1]=1-pnorm(1000000,forecast$pred[k],forecast$se[k]))
}
plot( D[length(D)-length(P)]+1:220,c(P,rep(NA,220-length(P))),
ylab="Probability to reach 1,000,000",xlab="",
type="l",col="red",ylim=c(0,1))
So, I guess my posts on multiple internal rates of return, or Young’s inequality will have to wait next year… I really need to find some more sexy post to attract readers.. Challenge accepted !

Probabilities, and opening doors (or boxes)

Recently, while we were in the car to Québec city with some PhD students, someone mentioned  the Monty Hall paradox, and we discussed possible extensions… The Monty Hall paradox is usually presented from tv show,

Craig F. Whitaker wrote the problem as follows, in Parade Magazine, September 1990, « Suppose you’re on a game show, and you’re given the choice of three doors: Behind one door is a car; behind the others, goats. You pick a door, say No. 1, and the host, who knows what’s behind the doors, opens another door, say No. 3, which has a goat. He then says to you, “Do you want to pick door No. 2?” Is it to your advantage to switch your choice? ». Actually, Bertrand proposed the same problem, but with boxes instead of doors…. but the problem was the same.

Assume that the candidate chooses door 1 (without loss of generality since the problem is clearly symmetric).  The probability that the car is behind door 2 is http://freakonometrics.blog.free.fr/public/maths/montyh52.png. The animator can either open door 2 or door 3:

  • if the car is behind door 3, he has to open door 2,
  • if the car is behind door 2, he has to open door 3,
  • if the car is behind door 1, he can open door 2 or 3, and we assume that the opening is equiprobable,

 

Assume that the animator says “the second box is empty”, what should the candidate do ?
To formalize the problem let http://freakonometrics.blog.free.fr/public/maths/montyh1.png denote the event that the car is behind door http://freakonometrics.blog.free.fr/public/maths/montyh2.png, and http://freakonometrics.blog.free.fr/public/maths/montyh3.png the event that the animator opens door http://freakonometrics.blog.free.fr/public/maths/montyh4.png. So, if he opens door 2, the probability that the car is behind door 3 is
http://freakonometrics.blog.free.fr/public/maths/montyh5.png
where
http://freakonometrics.blog.free.fr/public/maths/montyh7.png
from the previous discussion, since he cannot open door 1 (the candidate chose it) and the cannot open door 3 (since the car is behind). Further
http://freakonometrics.blog.free.fr/public/maths/montyh8.png
from equiprobability. And for http://freakonometrics.blog.free.fr/public/maths/montyh9.png we get, similarly
http://freakonometrics.blog.free.fr/public/maths/montyh10.png
Thus,
http://freakonometrics.blog.free.fr/public/maths/montyh11.png
while
http://freakonometrics.blog.free.fr/public/maths/montyh12.png
So the optimal strategy is to open the third door (even if I chose the first one)… It is usually seen as a paradox, but if you consider a much larger number of doors (say 4),

and that the animator opens 2 doors, then should we still change, and open the door that is still closed ? The higher the number of doors, the higher the probability to have something behind the other door…

For instance, with 4 doors, or boxes, if the candidate still chose the first door, and that the animator opens doors 2 and 3, then, the probability that the car is behind the fourth one is
http://freakonometrics.blog.free.fr/public/maths/montyh15.png
i.e.
http://freakonometrics.blog.free.fr/public/maths/montyh50.png
http://freakonometrics.blog.free.fr/public/maths/montyh52.png that appears at the denominator since we focus on the pair http://freakonometrics.blog.free.fr/public/maths/montyh54.png) while
http://freakonometrics.blog.free.fr/public/maths/montyh55.png
Once again, we have that
http://freakonometrics.blog.free.fr/public/maths/monthyh56.png
i.e. the opening of a doors bring us no information about our choice, but it will after conditional probability for the remaining door.
More generally, with n doors, if the animator opens http://freakonometrics.blog.free.fr/public/maths/monthy20.png doors,

then
http://freakonometrics.blog.free.fr/public/maths/montyh60.png
while
http://freakonometrics.blog.free.fr/public/maths/montyh62.png
And here the result is even more intuitive: we have to open the door that was left closed. Actually, it is possible to see that it can be extend to the case where the are http://freakonometrics.blog.free.fr/public/maths/montyh75.png doors (or boxes), the candidate chooses http://freakonometrics.blog.free.fr/public/maths/montyh73.png doors, and the animator opens http://freakonometrics.blog.free.fr/public/maths/montyh71.png (out of the http://freakonometrics.blog.free.fr/public/maths/montyh95.png remaining doors). Then, behind each door chosen by the candidate, the probability does not change
http://freakonometrics.blog.free.fr/public/maths/montyh70.png
where i goes from 1 to http://freakonometrics.blog.free.fr/public/maths/montyh73.png, while
http://freakonometrics.blog.free.fr/public/maths/montyh80.png
where i goes from http://freakonometrics.blog.free.fr/public/maths/montyh76.png to http://freakonometrics.blog.free.fr/public/maths/montyh75.png.

Is it that stupid to make extremely long term forecast when studying mortality ?

I received recently a comment by FCA (here) who raised an important question, about forecast in dynamic mortality models. (S)he mentioned that from his(her) point of view, the econometric models I considered were “good to predict for the next, say, 3 or 4 years. Not for the next 50 years…”. Which was the message I tried to stress last year in a conference about retirement in France (here). But from a quantitativepoint of view, how inconsistent were forecasts made 35 years ago, or 60 years ago ?

Consider here the Lee Carter model, obtained on the periods 1816-1950 (in black below), 1816-1975 (in red) and 1816-2000 (in blue), unfortunately, it is difficult to compare http://freakonometrics.blog.free.fr/public/maths/viekt.png‘s since we have identifiability problems here. Nevertheless, we if consider affine transformation so that  http://freakonometrics.blog.free.fr/public/maths/viekt.png‘s are equal in 1900 and 1950 (say), we obtain

On that graph, we considered an ETS (AAN) forecast. If we do not consider the entire series for forecasting, but only observations following WWI (1945), we obtain

For sketches of the R code,

T=1980
base0=data.frame(D,E,A,Y,a=as.factor(A),
y=as.factor(Y))
base=base0[base0$Y<=T,]
LC2=gnm(D~a+Mult(a,y),offset=log(E),family=
poisson,data=base)
A=LC2$coefficients[1]+LC2$coefficients[2:110]
B=LC2$coefficients[111:220]
K0=LC2$coefficients[221:length(LC2$coefficients)]
Y=as.numeric(K0)
K1=c(K0,forecast(ets(Y,model="AAN"),h=240)$mean)
K2=c(K0,forecast(auto.arima(Y,allowdrift=TRUE),h=240)$mean)
MU=matrix(NA,length(A),length(K1))
MU1=MU2=MU
for(i in 1:length(A)){
for(j in 1:length(K1)){
MU1[i,j]=exp(A[i]+B[i]*K1[j])
MU2[i,j]=exp(A[i]+B[i]*K2[j])
}}
x=40
s=seq(0,109-x-1)
t=2000
Pxt1=cumprod(exp(-diag(MU1[x+1+s,t+s-base1$Year[1]-1])))
Pxt2=cumprod(exp(-diag(MU2[x+1+s,t+s-base1$Year[1]-1])))
r=.035
m=70
h=seq(0,39)
V1=1/(1+r)^(m-x+h)*Pxt1[m-x+h]
V2=1/(1+r)^(m-x+h)*Pxt2[m-x+h]
M=cbind(V1,V2)
apply(M,2,sum)

Actually, it is not that bad…. even if it is only a qualitative intuition. Again, I am not a demographer, and my interest is more on actuarial science… so if we look at the estimation of annuities (still the same insurance contract, as here) for some insured of age 40 in 2000, we get the following graph (where forecasts http://freakonometrics.blog.free.fr/public/maths/viekt.png‘s were obtained on the complete series, i.e. from 1816 until the year we consider),

(here it means that in 1900, I had to forecast mortality for someone of age 40 in 2000… so we had to forecast mortality with a 150 year horizon). Obviously, even if we are able to forecast improvement of mortality rates, it is not enough since it looks like, each year, improvement are alway higher than what what expected. Note that if we run it twice (since there might be problem with initial values in the econometric procedure) we obtain something similar,

So, the output is consistent. And if we change the way we predict future values, e.g. on focusing only on the past 50 years, i.e.

K1=c(K0,forecast(ets(Y[(length(Y)-50):length(Y)],
model="AAN"),h=240)$mean)
K2=c(K0,forecast(auto.arima(Y[(length(Y)-50):length(Y)],
allowdrift=TRUE),h=240)$mean)

we obtain the following graph for the annuity associated to an insurance contract sold in 2000,

so that relative changes compared with 1980 are (in %)

Hence, over a bit more than 25 years, we underestimated annuities of 25%. We if start to take into account possible investments, it is not so bad, I think….  don’t you think ?

 

Finding roots of functions in actuarial science

The following simple code can be used to find roots of functions (based on the secant algorithm),

secant=function(fun, x0, x1, tolerence=1e-07, niter=500){
for ( i in 1:niter ) {
	x2 <- x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0))
	if (abs(fun(x2)) < tolerence)
		return(x2)
	x0 <- x1
	x1 <- x2
}}

It can be interesting in actuarial science, e.g. to find the actuarial rate so that to present values are equal. For instance, consider the following capital, given only if the insured is still alive (this example was initially considered here). We would like to find the rate so that the probable discounted value is 600,

> Lx=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/TV8890.csv",
+ header=TRUE,sep=";")
> capital=c(100,100,125,125,150,150)
> n=length(capital)
> x=0.035
> X=45
> f=function(x){
+ capital.act=capital*(1/(1+x))^(1:n)
+ PROBA=Lx[((Lx[,1]>X)*(Lx[,1]<=(X+n)))==1,2]/Lx[(Lx[,1]==X)==1,2]
+ return(sum(capital.act*PROBA))}
>
> f1=function(x){f(x)-600}
> secant(f1,0,0.1)
[1] 0.06022313
> f(0.06022313)
[1] 600

*

Millenium bridge, endogeneity and risk management

In less than 48 hours, last week two friends mentioned the Millennium Bridge as an illustration of a risk management concept. There are several documents with that example, here (for the initial idea of using the Millennium Bridge to illustrate issues in risk management) here or there, e.g.

When we mention resonance effects on bridges, we usually thing of the Tacoma Narrows Bridge (where strong winds set the bridge oscillating) or the Basse-Chaîne Bridge (in France, which collapsed on April 16, 1850, when 478 French soldiers marched across it in lockstep). In the first case, there is nothing we can do about it, but for the second one, this is why soldiers are required to break step on bridges.

But for the Millennium bridge, a ‘positive feedback‘ phenomenon (known as Synchronous Lateral Excitation in physics) has been observed: the natural sway motion of people walking caused small sideways oscillations in the bridge, which in turn caused people on the bridge to sway in step, increasing the amplitude of the oscillations and continually reinforcing the effect. That has been described in a nice paper in 2005 (here). In the initial paper by Jon Danielsson and Hyun Song Shin, they note that “what is the probability that a thousand people walking at random will end up walking exactly in step? It is tempting to say “close to zero”, or “negligible”. After all, if each person’s step is an independent event, then the probability of everyone walking in step would be the product of many small numbers – giving us a probability close to zero. Presumably, this is the reason why Arup – the bridge engineers – did not take this into account. However, this is exactly where endogenous risk comes in. What we must take into account is the way that people react to their environment. Pedestrians on the bridge react to how the bridge is moving. When the bridge moves under your feet, it is a natural reaction for people to adjust their stance to regain balance. But here is the catch. When the bridge moves, everyone adjusts his or her stance at the same time. This synchronized movement pushes the bridge that the people are standing on, and makes the bridge move even more. This, in turn, makes the people adjust their stance more drastically, and so on. In other words, the wobble of the bridge feeds on itself. When the bridge wobbles, everyone adjusts their stance, which sets off an even worse wobble, which makes the people adjust even more, and so on. So, the wobble will continue and get stronger even though the initial shock (say, a gust of wind) has long passed. It is an example of a force that is generated and amplified within the system. It is an endogenous response. It is very different from a shock that comes from a storm or an earthquake which are exogenous to the system.

And to go further, they point out that this event is rather similar to what is observed in financial markets (here) by quoting The Economist from October 12th 2000 “So-called value-at-risk models (VaR) blend science and art. They estimate how much a portfolio could lose in a single bad day. If that amount gets too large, the VAR model signals that the bank should sell. The trouble is that lots of banks have similar investments and similar VAR models. In periods when markets everywhere decline, the models can tell everybody to sell the same things at the same time, making market conditions much worse. In effect, they can, and often do, create a vicious feedback loop.