Arthur Charpentier

Blog transfert

Arthur Charpentier — Fri, 01 Oct 2010 00:01:00 +0200

As mentioned during the past few weeks, the blog has been transfered (please update links and bookmarks). You will be redirected (shortly) to

http://freakonometrics.blog.free.fr

Statistique de l'assurance STT6705V, partie 5

Arthur Charpentier — Wed, 29 Sep 2010 02:54:00 +0200

Pour la cinquième séance de cours, nous allons continuer sur la modélisation des coûts, à partir des transparents en ligne ici. La première partie portera sur la comparaison du modèle Gamma et du modèle lognormal, et la seconde sur l’écrêtement des gros sinistres.
Il convient de faire attention à la transformation logarithmique. Un modèle de la forme

paraîtra toujours "meilleur" qu’un modèle

(pour peu que les soient bien plus grands que 1). Regardons un jeu de données classique, liant vitesse du véhicule et distance de freinage.
> summary(lm(dist~speed,data=cars))
Call:
lm(formula = dist ~ speed, data = cars)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -17.5791     6.7584 -2.601   0.0123 *
speed         3.9324     0.4155   9.464 1.49e-12 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 15.38 on 48 degrees of freedom
Multiple R-squared: 0.6511,     Adjusted R-squared: 0.6438
F-statistic: 89.57 on 1 and 48 DF, p-value: 1.490e-12

> summary(lm(log(dist)~speed,data=cars))
Call:
lm(formula = log(dist) ~ speed, data = cars)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.67612    0.19614   8.546 3.34e-11 ***
speed        0.12077    0.01206 10.015 2.41e-13 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4463 on 48 degrees of freedom
Multiple R-squared: 0.6763,     Adjusted R-squared: 0.6696
F-statistic: 100.3 on 1 and 48 DF, p-value: 2.413e-13

Le R2 est plus faible, l’écart type des résidus aussi, la log-vraisemblance si on l’avait, etc.
Il ne faut jamais oublier qu’on doit ensuite prendre l’exponentiel sur modèle en log. C’est d’ailleurs ce qui est fait par exemple dans la transformation de Box Cox, où l’on compare les sommes des carrés des résidus sur les , i.e.

où .

> s=summary(lm(log(dist)~speed,data=cars))$sigma
> mean((cars$dist-predict(lm(dist~speed,data=cars)))^2)
[1] 227.0704
> mean((cars$dist-exp(predict(lm(log(dist)~speed,data=cars))+.5*s^2))^2)
[1] 296.2027
Tout ça pour conclure que le modèle linéaire n’est peut être pas si mauvais que ça....

EM and mixture estimation, with R (part 2)

Arthur Charpentier — Sat, 25 Sep 2010 02:56:00 +0200

Following my previous post on optimization and mixtures (here), Nicolas told me that my idea was probably not the most clever one (there).
So, we get back to our simple mixture model,

In order to describe how EM algorithm works, assume first that both and are perfectly known, and the mixture parameter is the only one we care about.

The simple model, with only one parameter that is unknown

Here, the likelihood is

so that we write the log likelihood as

which might not be simple to maximize. Recall that the mixture model can interpreted through a latent variate (that cannot be observed), taking value when is drawn from , and 0 if it is drawn from . More generally (especially in the case we want to extend our model to 3, 4, ... mixtures), and .
With that notation, the likelihood becomes

and the log likelihood

the term on the right is useless since we only care about p, here. From here, consider the following iterative procedure,
Assume that the mixture probability is known, denoted . Then I can predict the value of (i.e. and ) for all observations,

So I can inject those values into my log likelihood, i.e. in

having maximum (no need to run numerical tools here)

that will be denoted . And I can iterate from here.
Formally, the first step is where we calculate an expected (E) value, where is the best predictor of given my observations (as well as my belief in ). Then comes a maximization (M) step, where using , I can estimate probability .

A more general framework, all parameters are now unkown

So far, it was simple, since we assumed that and were perfectly known. Which is not reallistic. An there is not much to change to get a complete algorithm, to estimate . Recall that we had which was the expected value of Z_{1,i}, i.e. it is a probability that observation i has been drawn from .
If , instead of being in the segment was in , then we could have considered mean and standard deviations of observations such that =0, and similarly on the subset of observations such that =1.
But we can’t. So what can be done is to consider as the weight we should give to observation i when estimating parameters of , and similarly, 1- would be weights given to observation i when estimating parameters of .
So we set, as before

and then

and for the variance, well, it is a weighted mean again,

and this is it.

Let us run the code on the same data as before

Here, the code is rather simple: let us start generating a sample
> X1 = rnorm(n,0,1)
> X20 = rnorm(n,0,1)
> Z = sample(c(1,2,2),size=n,replace=TRUE)
> X2=4+X20
> X = c(X1[Z==1],X2[Z==2])
then, given a vector of initial values (that I called and then before),
> s = c(0.5, mean(X)-1, var(X), mean(X)+1, var(X))
I define my function as,
> em = function(X0,s) {
+ Ep = s[1]*dnorm(X0, s[2], sqrt(s[4]))/(s[1]*dnorm(X0, s[2], sqrt(s[4])) +
+ (1-s[1])*dnorm(X0, s[3], sqrt(s[5])))
+ s[1] = mean(Ep)
+ s[2] = sum(Ep*X0) / sum(Ep)
+ s[3] = sum((1-Ep)*X0) / sum(1-Ep)
+ s[4] = sum(Ep*(X0-s[2])^2) / sum(Ep)
+ s[5] = sum((1-Ep)*(X0-s[3])^2) / sum(1-Ep)
+ return(s)
+ }
Then I get , or . So this is it ! We just need to iterate (here I stop after 200 iterations) since we can see that, actually, our algorithm converges quite fast,
> for(i in 2:200){
+ s=em(X,s)
+ }

Let us run the same procedure as before, i.e. I generate samples of size 200, where difference between means can be small (0) or large (4),

Ok, Nicolas, you were right, we’re doing much better ! Maybe we should also go for a Gibbs sampling procedure ?... next time, maybe....

Optimization and mixture estimation, with R

Arthur Charpentier — Thu, 23 Sep 2010 23:09:00 +0200

Recently, one of my students asked me about optimization routines in R. He told me he that R performed well on the estimation of a time series model with different regimes, while he had trouble with a (simple) GARCH process, and he was wondering if R was good in optimization routines. Actually, I always thought that mixtures (and regimes) was something difficult to estimate, so I was a bit surprised...

Indeed, it reminded me some trouble I experienced once, while I was talking about maximum likelihooh estimation, for non standard distribution, i.e. when optimization had to be done on the log likelihood function. And even when generating nice samples, giving appropriate initial values (actually the true value used in random generation), each time I tried to optimize my log likelihood, it failed. So I decided to play a little bit with standard optimization functions, to see which one performed better when trying to estimate mixture parameter (from a mixture based sample). Here, I generate a mixture of two gaussian distributions, and I would like to see how different the mean should be to have a high probability to estimate properly the parameters of the mixture.

The density is here proportional to

The true model is , and being a parameter that will change, from 0 to 4.
The log likelihood (actually, I add a minus since most of the optimization functions actually minimize functions) is
> logvraineg <- function(param, obs) {
+ p <- param[1]
+ m1 <- param[2]
+ sd1 <- param[3]
+ m2 <- param[4]
+ sd2 <- param[5]
+ -sum(log(p * dnorm(x = obs, mean = m1, sd = sd1) + (1 - p) *
+ dnorm(x = obs, mean = m2, sd = sd2)))
+ }
The code to generate my samples is the following,
>X1 = rnorm(n,0,1)
> X20 = rnorm(n,0,1)
> Z = sample(c(1,2,2),size=n,replace=TRUE)
> X2=m+X20
> X = c(X1[Z==1],X2[Z==2])
Then I use two functions to optimize my log likelihood, with identical intial values,
> O1=nlm(f = logvraineg, p = c(.5, mean(X)-sd(X)/5, sd(X), mean(X)+sd(X)/5, sd(X)), obs = X)
> logvrainegX <- function(param) {logvraineg(param,X)}
> O2=optim( par = c(.5, mean(X)-sd(X)/5, sd(X), mean(X)+sd(X)/5, sd(X)),
+ fn = logvrainegX)
Actually, since I might have identification problems, I take either or , depending whether or is the smallest parameter.

On the graph above, the x-axis is the difference between means of the mixture (as on the animated grap above). Then, the red point is the median of estimated parameter I have (here ), and I have included something that can be interpreted as a confidence interval, i.e. where I have been in 90% of my scenarios: the black vertical segments. Obviously, when the sample is not enough heterogeneous (i.e. and rather different), I cannot estimate properly my parameters, I might even have a probability that exceed 1 (I did not add any constraint). The blue plain horizontal line is the true value of the parameter, while the blue dotted horizontal line is the initial value of the parameter in the optimization algorithm (I started assuming that the mixture probability was around 0.2).
The graph below is based on the second optimization routine (with identical starting values, and of course on the same generated samples),

(just to be honest, in many cases, it did not converge, so the loop stopped, and I had to run it again... so finally, my study is based on a bit less than 500 samples (times 15 since I considered several values for the mean of my second underlying distribution), with 200 generated observations from a mixture).
The graph below compares the two (empty circles are the first algorithm, while plain circles the second one),

On average, it is not so bad.... but the probability to be far away from the tru value is not small at all... except when the difference between the two means exceeds 3...
If I change starting values for the optimization algorithm (previously, I assumed that the mixture probability was 1/5, here I start from 1/2), we have the following graph

which look like the previous one, except for small differences between the two underlying distributions (just as if initial values had not impact on the optimization, but it might come from the fact that the surface is nice, and we are not trapped in regions of local minimum).
Thus, I am far from being an expert in optimization routines in R (see here for further information), but so far, it looks like R is not doing so bad... and the two algorithm perform similarly (maybe the first one being a bit closer to the true parameter).

Statistique de l'assurance STT6705V, partie 4b

Arthur Charpentier — Thu, 23 Sep 2010 21:17:00 +0200

Avant toute autre chose, pour revoir le dernier cours, il suffit d’aller ici et là.
Bon, sinon, comme promis, les bases de données pour les projets sont en ligne ici. Le principe est simple. Il y a 28 bases de données, toutes semblables (mais bien sûr différentes), par les numéros ci-dessous. Comme toujours, premier arrivé, premier servi, donc les bases vont être attribuées au fur et à mesure

étape 1: choisir une base

Vous choisissez un numéro et cliquez dessus, afin de me prévenir (par courriel, à l’aide du tableau ci-dessous) de votre choix. Les projets se font par deux, merci de mettre dans le corps du message les noms des deux personnes qui travailleront sur le projet.

étape 2: récupérer les données

une fois que j’aurais confirmé que personne n’a choisi la même base, vous pouvez les récupérer sous R à l’aide du code suivant,
> k=1
> nom=paste("http://perso.univ-rennes1.fr/arthur.charpentier/6705V/UdM-baseC-",
+ k,".txt",sep="")
> baseC=read.table(nom,header=TRUE)
pour la base de contrats (et pour le groupe qui aurait choisi la base 1) et pour la base des sinistres,
> nom=paste("http://perso.univ-rennes1.fr/arthur.charpentier/6705V/UdM-baseS-",
+ k,".txt",sep="")
> baseS=read.table(nom,header=TRUE)
Ensuite, c’est parti, il s’agit de me proposer différents modèles de tarification, et de calculer les primes pures avec les différents modèles, pour une personne parmi les listes des personnes ayant les caractéristiques suivantes,
> client=data.frame(
+ exposition=rep(1,9),
+ zone=c("A","A","A","C","D","E","F","F","F"),
+ puissance=c(6,7,11,6,7,11,6,7,11),
+ agevehicule=c(0,1,5,10,5,1,0,6,10),
+ ageconducteur=c(25,18,55,55,55,40,21,20,18),
+ bonus=c(80,100,50,60,55,50,100,125,100),
+ marque=c(1,2,12,12,12,1,1,1,2),
+ carburant=c("D","E","E","D","D","E","E","D","D"),
+ densite=rep(3000,9),
+ client=rep(baseC$region[1],9))
(si certains modalités ne sont pas présentes dans la base, il faut choisir quelqu’un d’autre.... sur les neuf, il doit bien en avoir un(e) qui pourrait être présent(e) dans votre base de données). En cas de problème, vous avez mon adresse électronique.... Et je reviendrais ultérieurement sur la forme de ce que j’attends.

Soccer, probabilité (et assurance) partie 3

Arthur Charpentier — Thu, 23 Sep 2010 17:54:00 +0200

La dernière fois (il y a quelques mois, ici) on avait évoqué le lien entre cote et probabilités risques neutres (induites). J’avais évoqué les cotes à un instant donné donné, et montré comment on pouvait construire la probabilité risque neutre associée. L’idée est ici de suivre l’évolution des probabilités pendant la coupe du monde, au fur et à mesure que les matchs révèlent de l’information sur le vrai niveau des équipes...
A l’aide des données "World Cup Group A Betting" ("Win Market") sur http://www.oddschecker.com/ (Vincent, alias @Vicnent, avait fait des sauvegarde régulières des pages html), on peut suivre les cotes toutes les 30 minutes, entre le 10 juin et le 22 juin, au sein du groupe de l’équipe de France.

Pour le premier site de pari (bet365), on a les probabilités suivantes, avec les probabilités associées à l’équipe de France (en bleu), à l’Afrique du Sud (en jaune), l’Uruguay (en rouge) et le Mexique (en vert),

On note que les probabilités sont stables en dehors des matchs, autrement dit seul le comportement sur le terrain semble intéresser les parieurs (les matchs ont lieu pendant les date où l’on observe des traits bleus clairs verticaux).
Pour le troisième, on observe quelque chose de très proche (les valeurs manquants signifient que le tableau en ligne était vide, ou que j’ai raté ma lecture du fichier html sous R), avec des variations du même ordre à la fin des matchs,

Pour le sixième on a

pour le neuvième,

et pour le douzième (je n’ai pas affiché tous les graphs)

On notera qu’à la fin, après les seconds matchs plus personnes ne semblait croire à l’équipe de France (en France en tous les cas), sa probabilité de gagner pour les parieurs est resté élevée, et ce, chez tous les sites de paris.... il aurait probablement été intéressant de parier à la fin sur le Mexique (même si je sais que rétrospectivement, il est toujours facile de dire ce qu’il fallait alors faire). Attendre le dernier moment pour faire des paris n’est pas forcément stupide...
La prochaine étape c’est de regarder non pas les pronostiques par poule, mais globalement, sur le futur vainqueur. Nous verrons en particulier comment les matchs des autres influences les probabilités d’une équipe....

Qui peut m'aider à comprendre les sorties de SAS ?

Arthur Charpentier — Thu, 23 Sep 2010 01:05:00 +0200

Je m’étais promis que j’évoquerais une bizarrerie rencontrée avec SAS lors d’une formation.... Écrire ce billet permettra à ceux qui auraient des éléments d’explication de poster un commentaire.
Pour cela, comparons une régression logistique faite avec deux outils différents, sous SAS,

avec la procédure logistique

Le code pour faire une régression logistique ressemble à ça
PROC LOGISTIC DATA=base_logistq;
FORMAT age_soc f2_ageso.;
CLASS sexe_soc age_soc fract_paiemt;
MODEL SPOCAM = sexe_soc age_soc fract_paiemt / selection=stepwise;
RUN; QUIT;
ce qui donne la sortie suivante (je passe l’introduction pour insister sur les coefficients)
                                 The LOGISTIC Procedure

                   Analyse des estimations de la vraisemblance maximum
                                                     Erreur         Khi 2
Paramètre                    DF    Estimation         std       de Wald    Pr > Khi 2

Intercept                     1        1.7833      0.0676      696.9022        <.0001
sexe_soc     Femme            1       -0.2429      0.0619       15.4237        <.0001
age_soc      1_AGESOC_-60     1        0.4578      0.0667       47.1020        <.0001
fract_paiemt Annuel           1        0.6021      0.0997       36.4862        <.0001
fract_paiemt Mensuel          1       -0.5410      0.0842       41.2342        <.0001

avec la procédure genmod (car la régression logistique est un glm)

On peut faire exactement la même chose (théoriquement) en ajustement un modèle GLM,
PROC GENMOD DATA=base_logistq;
FORMAT age_soc f2_ageso.;
CLASS sexe_soc age_soc fract_paiemt;
MODEL SPOCAM = sexe_soc age_soc fract_paiemt / dist = binomial;
RUN;
et la sortie ressemble à ça
                                  The GENMOD Procedure
                      Analyse des résultats estimés de paramètres

                                                  Erreur      Wald 95Limites
Paramètre                     DF   Estimation   standard      de confiance %       Khi 2
Intercept                      1       1.5073     0.1501     1.2131     1.8014    100.85
sexe_soc       Femme           1      -0.4859     0.1237    -0.7284    -0.2434     15.42
sexe_soc       Homme           0       0.0000     0.0000     0.0000     0.0000       .
age_soc        1_AGESOC_-60    1       0.9156     0.1334     0.6542     1.1771     47.10
age_soc        Z_AGESOC_+60    0       0.0000     0.0000     0.0000     0.0000       .
fract_paiemt   Annuel          1       0.6634     0.1770     0.3165     1.0104     14.05
fract_paiemt   Mensuel         1      -0.4798     0.1510    -0.7759    -0.1838     10.09
fract_paiemt   Semestriel      0       0.0000     0.0000     0.0000     0.0000       .
Scale                          0       1.0000     0.0000     1.0000     1.0000

comparaison des deux sorties

Si on regarde l’impact du sexe par exemple, dans la première sortie on peut lire
sexe_soc     Femme            1       -0.2429      0.0619       15.4237        <.0001
alors que dans la seconde sortie, on a
sexe_soc       Femme           1      -0.4859     0.1237    -0.7284    -0.2434     15.42
sexe_soc       Homme           0       0.0000     0.0000     0.0000     0.0000
On dira ce qu’on veut, mais moi je trouve cette différence troublante.... Dans la seconde sortie, le coefficient vaut le double de l’autre....
Alors SAS semble s’y retrouver car si on lui demande d’afficher le score prédit pour un individu au hasard (le premier de la base par exemple), les prédictions sont très proches,
                           fract_                                 proba1_       proba1_
Obs sexe_soc   age_soc paiemt      SPOCAM proba1_logit
1    Homme          71 Annuel         0      0.10242637    0.10241302
Si quelqu’un sait interpréter ce qui est fait avec cette procédure logistique (car R donne la même chose que la sortie GLM), je suis preneur.....

Statistique de l'assurance STT6705V, partie 4

Arthur Charpentier — Wed, 22 Sep 2010 03:15:00 +0200

On continue le cours avec la fin sur la modélisation de la fréquence, où on parlera de surdispersion (évoquée ici par exemple), de quasi-lois (ici ou là), de régression binomiale négative (ici), ou de modèles à inflation de zéros (là par exemple). Et si ça ne suffit pas, on commencera à parler de la modélisation des coûts individuels. Les lois de bases seront les lois lognormales, et gamma,

Les transparents seront bientôt en ligne ici.

Dans le prochain billet, je mettrais les liens vers les bases de données pour le projet d’économétrie. Sinon pour les rappels de la semaine dernière, je renvoie Sinon pour les rappels de la semaine dernière, je renvoie ici et l à.

Espérance de vie à la naissance, ou à 65 ans ?

Arthur Charpentier — Tue, 14 Sep 2010 20:21:00 +0200

Bon, je vais faire un rapide billet pour répondre à un commentaire qui avait été posté ici, qui remettait en cause mes conclusions (ce qui ne me gêne pas, loin de là, c’est comme ça qu’on va faire avancer les choses), et surtout demandait des informations complémentaires. J’avais un peu forcé le trait en affirmant que "une personne de 70 ans vivra à peine plus longtemps - en moyenne - qu’une personne de 70 ans en 1950" (c’est peut être le cas pour les centenaires, mais je n’ai pas de données), mais le but de mon précédant billet était de montrer que les gains sont beaucoup plus faibles que pour les plus jeunes.
J’ai décidé de refaire l’étude, de manière plus claire peut être, et plus complète. Et en utilisant des données de l’INED, qui est une source a priori fiable de données.

Les données (détaillées)

Comme je le disais, tout d’abord, les données: je prends ici les données de l’INED (et pas de http://mortalityorg). On peut trouver en ligne (au format xls) des données de population par âge et par sexe, ici (populations par âge, de 0 à 100 ans, au premier janvier, de 1899 à 1998, selon le territoire couvert par la statistique des décès, ensemble des deux sexes, avec depuis 1946, les données INSEE et avant 1946, des reconstitutions intercensitaires (Vallin, 1973)). Et là, on a les données du nombre de décès par âge et par sexe (décès par âge et par génération, de 1899 à 1997, ensemble des deux sexes, avec pour les moins de 100 ans, depuis 1907, décès hors pertes militaires (SGF, puis INSEE), pour les années de guerre : décès (INSEE + estimations des pertes militaires (Vallin, 1973)) et entre 1899 et 1906 : décès par âge de l’INSEE, répartition entre les triangles de Lexis (Vallin, 1973). Pour les centenaires, jusqu’en 1967 : SGF, INSEE + répartitions estimées par les auteurs, et depuis 1968 : INSEE (extrait d’enregistrements individuels fournis par l’INSEE dans le cadre d’un avenant à la convention INED-INSEE)). Voilà pour les données, ou presque....
Par paresse, et pour pouvoir partager les données plus facilement, j’ai mis en ligne des fichiers csv sur ma page.
Les données (brutes) contiennent des colonnes sans intérêt, que je me suis permis de virer. Sinon dans la base des populations, certaines lignes ont été supprimées car en 1914 il y avait deux populations, sur territoire d’avant 1914 (sans l’Alsace-Lorraine: Moselle, Bas-Rhin, Haut-Rhin), et le territoire sans opérations militaires (manquent, outre l’Alsace-Lorraine: Aisne, Ardennes, Marne, Meurthe-et-Moselle, Meuse, Nord, Oise, Pas-de-Calais, Somme et Vosges). C’est cette dernière qui est retenues dans la base des décès c’est donc cette dernière que l’on utilisera pour calculer l’exposition.
> tabB=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/tabB.csv",
+ sep=";",header=FALSE)
> ANNEE=tabB[,1]
> BASEB=tabB[,seq(2,246,by=2)]
> BASEB=BASEB[,1:100]
> tabC=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/tabC.csv",
+ sep=";",header=FALSE)
> BASEC=tabC[,2:101]
> BASEC=BASEC[-c(16,23,43,48,51,53),]
> BASEC=BASEC[1:nrow(BASEB),]
> AGE=0:99

Les taux de décès (population entière)

Commençons par calculer, et tracer les taux de décès, qui sont les ratios entre les nombres de décès et la population totale.
Mon premier point (pour répondre au commentaire) est que si l’on regarde un peu la surface, les gains en terme de taux de mortalité instantanés sont plus importants à la naissance (courbe rouge) qu’à des âges beaucoup plus élevés, comme 90 ans (en bleu) voire à 65 ans (en vert),

Un peu de prospective

Comme on le voit, sur nos données, on est limité car on n’a qu’un siècle de données, et c’est donc délicat de suivre les espérances de vie.
Le plus simple est alors de reprendre ce que j’avais fait (en ligne ici), utilisons les fonctions de Rob Hyndman, en ligne sur son blog (là), mais sur nos données (la dernière fois j’avais utilisé les siennes).
>   library(demography)
Le chargement a nécessité le package : forecast
Le chargement a nécessité le package : fracdiff
This is forecast 2.05
Le chargement a nécessité le package : rainbow
Le chargement a nécessité le package : hdrcde
Le chargement a nécessité le package : locfit
Le chargement a nécessité le package : akima
Le chargement a nécessité le package : lattice
locfit 1.5-6     2010-01-20
Le chargement a nécessité le package : ash
Le chargement a nécessité le package : ks
Le chargement a nécessité le package : KernSmooth
KernSmooth 2.23 loaded
Copyright M. P. Wand 1997-2009
Le chargement a nécessité le package : mvtnorm
hdrcde 2.13 loaded
Attachement du package : ’hdrcde’
Le chargement a nécessité le package : MASS
Le chargement a nécessité le package : pcaPP
Le chargement a nécessité le package : cluster
Le chargement a nécessité le package : ftsa
This is demography 1.03
Il y a eu 13 avis (utilisez warnings() pour les visionner)
Commençons par mettre les données au bon format. Pour des raisons numériques, je suis obligé de me débarrasser des âges élevés pour lesquels j’ai pu avoir une population nulle, ou pour lesquels j’ai des valeurs manquantes. J’ai donc exclu les âges supérieurs à 90 ans,
> donnees = demogdata(data=t(as.matrix(BASEB[,1:90]))/t(as.matrix(BASEC[,1:90])),
+                     pop= t(as.matrix(BASEC[,1:90])),
+                     ages=AGE[1:90],
+                     years=ANNEE,
+                     type="mortality", label="France", name="total", lambda=0)
> donnees
Mortality data for France
    Series: total
    Years: 1899 - 1997
    Ages: 0 - 89
Ensuite, on peut reprendre le code que j’avais fait. Notons qu’on utiliser ici une modélisation à la Lee et Carter, ce qui a été utilisée par l’INED pour extrapoler certaines données dans la table.
>   france.LC2 <- lca(donnees,adjust="dt",series="total")
>   france.fcast <- forecast(france.LC2)
>   L2 <- lifetable(france.fcast)
>   ex2=L2$ex
>   L1=lifetable(donnees,series="total")
>   ex1=L1$ex
Pour l’espérance de vie à la naissance, et son évolution dans le temps, si on regarde les gains moyens année après année, on voit que l’on gagne, en moyenne, 1/4 d’espérance de vie, soit 3 mois par an (comme dans mon précédant billet),
>   age=0
>   ex=c(ex1[age+1,],ex2[age+1,])
>   plot(1899:2047,ex,col="blue")
>   I=(1950:2000)-1898
>   y=ex[I]
>   x=1950:2000
>   lm(y~x)
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept)            x
-421.0384       0.2549
>   abline(lm(y~x),col="red")
>   points(x,y,pch=19,col="red")

On notera qu’en utilisant les autres méthodes d’estimation, j’obtiens la même pente. En revanche, si je regarde l’espérance de vie résiduelle à 65 ans, elle n’est plus que de 0,16 années, soit moins de deux mois chaque années,
>   age=65
>   ex=c(ex1[age+1,],ex2[age+1,])
>   plot(1899:2047,ex,col="blue")
>   I=(1950:2000)-1898
>   y=ex[I]
>   x=1950:2000
>   lm(y~x)
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept)            x
-295.6812       0.1612
>   abline(lm(y~x),col="red")
>   points(x,y,pch=19,col="red")

Dit autrement, entre un bébé qui nait en 1975 et un bébé qui nait en 1995, il y a un gain d’environ 6 ans d’espérance de vie, passant de 81,8 ans à 87,7, alors que pour une personne qui atteignait les 65 ans en 1975, il avait 22 ans à vivre encore, en moyenne, contre seulement 26,5 pour une personne de 65 ans en 1995. Donc je peux maintenir ce que j’affirmais la dernière fois, sur d’autres données.... En fait, si on veut être plus précis, on peut précisément regarder l’évolution de la pente, correspondant au gain moyen d’espérance de vie chaque année,

Bref, de considérables progrès ont été effectués en terme de gains d’espérance de vie, mais beaucoup plus sur les jeunes que sur les personnes plus âgées...

Modèles GLM et variables qualitatives

Arthur Charpentier — Sat, 11 Sep 2010 23:49:00 +0200

Suite à une question lors du cours de mercredi dernier (et comme Brice m’a posé exactement la même par mail tout à l’heure), je vais faire un rapide billet sur les liens entre "modèle(s) GLM sur variable(s) qualitative(s)" (exclusivement) et "moyenne par classe".
Ce que je disais en classe, c’est qu’en tarification, faire des GLM sur des données exclusivement qualitatives, c’était équivalent à calculer la moyenne par modalités (par classe tarifaire finalement). J’en avais déjà parlé dans un billet sur la méthode des marges (car c’est de ça qu’on parle en fait), ici. Mais ici, je voulais illustrer ce point sur un exemple (ou plutôt plusieurs exemples), afin d’expliquer ce que j’entendais pas "c’est la même chose"...

en présence d’une unique variable qualitative

Bon, commençons par simuler un jeu de données pour faire de l’économétrie...
> set.seed(1)
> X = as.factor(sample(c("A","B","C","D"),size=1000,
+      replace=TRUE,prob=c(.4,.3,.2,.1)))
> F=data.frame(X=c("A","B","C","D"))
> M = rep(NA,1000)
> M[X=="A"]=1
> M[X=="B"]=1.5
> M[X=="C"]=2
> M[X=="D"]=2.5
> Z = rexp(1000,rate=M)
> Y = rpois(1000,Z)
> base=data.frame(X,Y)
>   table(base)
   Y
X     0   1   2   3   4   5   6   7   8 11
A 188 102 62 28 12   3   6   2   1   1
B 169 69 22 20   8   2   0   1   0   0
C 140 39   8   6   1   0   0   0   0   0
D 85 17   5   3   0   0   0   0   0   0
Le bon modèle est un modèle de Poisson avec un lien logarithmique.... Regardons déjà ce que sont les moyennes par modalités,
>   mean(Y[X=="A"])
[1] 1.091358
>   mean(Y[X=="B"])
[1] 0.7628866
>   mean(Y[X=="C"])
[1] 0.3969072
>   mean(Y[X=="D"])
[1] 0.3272727
Pour comparer avec des modèles GLM, commençons par le bon modèle,
> predict(reg,newdata=F,type="response")
        1         2         3         4
1.0913580 0.7628866 0.3969072 0.3272727
ce qui est exactement ce que nous avions, en moyenne, par classes. Que se passe-t-il si on se trompe de modèle ? par exemple avec une autre fonction lien,
> reg=glm(Y~X,data=base,family=poisson(link="identity"))
> predict(reg,newdata=F,type="response")
        1         2         3         4
1.0913580 0.7628866 0.3969072 0.3272727
et si on change aussi la loi, là aussi on retombe sur la moyenne par modalité...
> reg=glm(Y~X,data=base,family=gaussian(link="identity"))
> predict(reg,newdata=F,type="response")
        1         2         3         4
1.0913580 0.7628866 0.3969072 0.3272727
On peut aussi tenter une régression binomiale négative par exemple,
> library(MASS)
> reg=glm.nb(Y~X,data=base)
> predict(reg,newdata=F,type="response")
        1         2         3         4
1.0913580 0.7628866 0.3969072 0.3272727
Moralité, on retrouve la même chose avec les deux approches: moyenne par modalité, ou GLM.

en présence d’une unique variable qualitative et d’une exposition différente

En pratique, si on cherche à modéliser la fréquence individuelle de sinistres, il convient de prendre en compte l’exposition (c’est à dire le fait que le nombre de sinistres observés dans la base n’a pas été obtenu sur une année complète, mais parfois quelques mois). Classiquement, le nombre de sinistres sera supposé proportionnel à cette exposition. Considérons la base suivante,
> set.seed(1)
> X = as.factor(sample(c("A","B","C","D"),size=1000,
+      replace=TRUE,prob=c(.4,.3,.2,.1)))
> E = runif(1000)
> F=data.frame(X=c("A","B","C","D"))
> M = rep(NA,1000)
> M[X=="A"]=1
> M[X=="B"]=1.5
> M[X=="C"]=2
> M[X=="D"]=2.5
> Z = rexp(1000,rate=M/E)
> Y = rpois(1000,Z)
> base=data.frame(X,Y)
>   table(base)
   Y
X     0   1   2   3   4   5   6   7
A 287 76 24   7   6   2   2   1
B 222 44 17   8   0   0   0   0
C 156 29   7   2   0   0   0   0
D 90 15   2   2   1   0   0   0
Les moyennes empiriques par modalité doivent être pondérées par ces expositions. Le plus simple est de "normaliser", en modifiant les données pour les passer sur une base annuelle (en divisant par l’exposition), puis en attribuant des poids proportionnels à cette exposition,
>   weighted.mean(Y[X=="A"]/E[X=="A"],E[X=="A"])
[1] 1.001059
>   weighted.mean(Y[X=="B"]/E[X=="B"],E[X=="B"])
[1] 0.6973345
>   weighted.mean(Y[X=="C"]/E[X=="C"],E[X=="C"])
[1] 0.5110974
>   weighted.mean(Y[X=="D"]/E[X=="D"],E[X=="D"])
[1] 0.5758276
Comme l’exposition est supposée intervenir de manière proportionnelle (i.e. multiplicative), il est logique de supposer un lien logarithmique, et si l’on considère une loi de Poisson, on obtient (en faisant une prédicition pour une exposition d’un an)
> base=data.frame(X,Y,E)
> F=data.frame(X=c("A","B","C","D"),E=1)
> reg=glm(Y~X+offset(log(E)),data=base,family=poisson(link="log"))
> predict(reg,newdata=F,type="response")
        1         2         3         4
1.0010593 0.6973345 0.5110974 0.5758276
Bref, encore une loi les résultats sont identiques. Mais, si on fait une régression binomiale négative, on obtient des résultats légèrement différents (0,17%, 0,03%, 1,02% et 3,40% respectivement), mais il ne sont pas égaux,
> reg=glm.nb(Y~X+offset(log(E)),data=base)
> predict(reg,newdata=F,type="response")
        1         2         3         4
0.9993223 0.6975693 0.5163053 0.5562640
Autrement dit, si on doit prendre en compte une exposition, les deux méthodes ne coïncident plus.

en présence de plusieurs variables qualitatives

Pour conclure, regardons ce qui se passe avec deux variables qualitatives,
> set.seed(1)
> X1 = as.factor(sample(c("A","B","C","D"),size=1000,
+       replace=TRUE,prob=c(.4,.3,.2,.1)))
> X2 = as.factor(sample(c("H","F"),size=1000,
+       replace=TRUE,prob=c(.6,.4)))
> M = rep(NA,1000)
> M[(X1=="A")&(X2=="H")==TRUE]=1
> M[(X1=="B")&(X2=="H")==TRUE]=1.5
> M[(X1=="C")&(X2=="H")==TRUE]=2
> M[(X1=="D")&(X2=="H")==TRUE]=2.5
> M[(X1=="A")&(X2=="F")==TRUE]=1.1
> M[(X1=="B")&(X2=="F")==TRUE]=1.4
> M[(X1=="C")&(X2=="F")==TRUE]=2.1
> M[(X1=="D")&(X2=="F")==TRUE]=2.7
> Z = rexp(1000,rate=M)
> Y = rpois(1000,Z)
> base=data.frame(Y,X1,X2)
Bon, on a compris l’idée, je vais me contenter de faire les calculs pour deux types d’individus,
>   mean(Y[(X1=="A")&(X2=="H")==TRUE])
[1] 0.9227642
>   mean(Y[(X1=="C")&(X2=="F")==TRUE])
[1] 0.4473684
Si on fait une régression sur les deux facteurs, on obtient
> F=data.frame(X1=rep(c("A","B","C","D"),2),X2=rep(c("H","F"),each=4))
> reg=glm(Y~X1+X2,data=base,family=poisson(link="log"))
> predict(reg,newdata=F,type="response")
        1         2         3         4         5         6         7         8
0.9081317 0.6647936 0.4915430 0.3998411 0.9723245 0.7117857 0.5262885 0.4281045
Les résultats sont différents des moyennes empiriques, tout simplement car on suppose ici les facteurs multiplicatifs,
> summary(reg)
Call:
glm(formula = Y ~ X1 + X2, family = poisson(link = "log"), data = base)
Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.02807    0.06824 -0.411 0.680851
X1B         -0.31191    0.08759 -3.561 0.000369 ***
X1C         -0.61384    0.11336 -5.415 6.12e-08 ***
X1D         -0.82032    0.15778 -5.199 2.00e-07 ***
X2H         -0.06830    0.07587 -0.900 0.368009
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
    Null deviance: 1581.6 on 999 degrees of freedom
Residual deviance: 1525.5 on 995 degrees of freedom
AIC: 2489.6
Number of Fisher Scoring iterations: 6
> exp(coefficients(reg)[5])
    X2H
0.93398
Autrement dit, pour un homme, je prédis toujours être 6,4% inférieur à une femme. Or il n’y a pas de raison que les deux facteurs interviennent indépendamment... il faut donc plutôt regarder le produit cartésien des modalités,
> predict(reg,newdata=F,type="response")
        1         2         3         4         5         6         7         8
0.9227642 0.6416185 0.5423729 0.3243243 0.9496855 0.7457627 0.4473684 0.5833333
On a alors 8 modalités possibles (le croisement des modalités de chacune des variables), et cette fois, on retombe sur les moyennes par classes.
Moralité, les GLM sur des facteurs peuvent coïncider avec les moyennes empiriques par classes, mais à condition de le faire en suivant quelques règles (en particulier croiser toutes les modalités mais aussi faire attention à l’exposition).

Ramadan, entre le sexe et la faim...

Arthur Charpentier — Fri, 10 Sep 2010 14:35:00 +0200

Dans un précédant billet (ici) je volais un dessin de Martin Vidberg qui parlait de google et de statistiques... Google trend est un outil magique pour noter des tendances et faire des statistiques simples et marrantes. Par exemple, puisque le ramadan touche à sa fin, on peut regarder ce qui motivent les algériens pendant le ramadan. Manifestement, le sexe n’est pas une priorité, que ce soit en anglais,

ou en français,

De manière générale, le mot clé le plus tapé sous google a tendance à se faire tout petit pendant la période du ramadan, alias porno,

En fait, on cherche un peu moins les femmes sur le net,

Par contre, visiblement, on a faim, et on veut bien manger: cuisine explose pendant le ramadan,

On peut aussi regarder les autres pays à forte population musulmane, et observer des tendances similaires, comme le Maroc,

les Émirats Arabes Unis,

ou encore le Pakistan,

(dans ce dernier cas, on observe un fort retour du mot clé sexune fois le ramadan passé... une espèce de manque virtuel). Bref, avec Google, on peut vraiment observer des choses intéressantes... non ?

Le million ! le milllion !

Arthur Charpentier — Thu, 09 Sep 2010 16:08:00 +0200

Hier soir (ou ce matin, je suis perdu avec ce décalage horaire) Christelle me demandait de parler un peu de prévision avec R (ici). Au lieu de renvoyer vers l’aide en ligne, penons un exemple pratique (et simple, voire si possible intéressant): la fréquentation d’un blog (en l’occurrence http://blogperso.univ-rennes1.fr/arthur.charpentier/). Considérons le nombre de pages vues, par jour, selon Google Analytics.
> base=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/million.csv",
+ sep=";",header=TRUE
)> X=base$nombre
> D0=as.Date("08/11/2008","%d/%m/%Y")
> D=D0+1:length(X)
> plot(D,X)
> plot(acf(X,lag=90))
La série a l’allure suivante (oui, le compteur n’a été installé qu’il y a deux ans),

Les autocorrélations montre une forte saisonnalité hebdomadaire, avec ~~moins~~ plus (cf ici) de consultations en semaine que le week-end,

La question que l’on cherche à résoudre est "ai-je une chance d’atteindre le million de pages vues d’ici la fin de l’année ?".

On peut traduire cette question de deux manières

quelle est la probabilité que le 1er janvier, j’ai atteint le million de pages vues,
quelle est la probabilité que la date où le million de pages vues sera atteint soit avant le 1er janvier

Une fois formalisée la question, reste à faire un peu d’économétrie.

modélisation économétrique

En faisant simple et rapide, afin de prendre en compte la corrélation forte avec la semaine précédente, et le fait que l’on s’intéresse à la somme cumulée, on peut considérer un modèle ARIMA

avec un retard d’ordre 7 pour les composantes moyennes mobiles et autorégressives,
avec une série intégrée à l’ordre 1,

L’ajustement se fait de la manière suivante,
> X=cumsum(base$nombre)
> model <- arima(X,c(7 ,   # partie AR
+                     1,    # partie I
+    7))   # partie MA
et la prévision se fait via
> forecast <- predict(model,200)
Ensuite, ce n’est qu’une représentation graphique,
> u=max(D)+1:200
> polygon( c(u,rev(u)), c(forecast$pred - 1.96*forecast$se,
+ rev(forecast$pred + 1.96*forecast$se)), col = "yellow", border=NA)
> lines(u,forecast$pred,col="blue",lwd=2)
> lines(u,forecast$pred- 1.96*forecast$se,col="blue",lty=2)
> lines(u,forecast$pred+ 1.96*forecast$se,col="blue",lty=2)
> abline(v=as.Date("01/01/2011","%d/%m/%Y"),col="red")
> abline(h=1000000,col="red")

Afin de répondre à la question posée, on peut étudier les différentes probabilités envisagées,

probabilité que le 1er janvier, j’ai atteint le million de pages vues

Dans un premier temps, on utilise la normalité des prédictions (en supposant une normalité du bruit) pour obtenir la loi de la prédiction à une date quelconque,
> k=which(u==as.Date("01/01/2011","%d/%m/%Y"))
> x=seq(800000,1100000,by=100)
> y=dnorm(x,mean=forecast$pred[k],sd=forecast$se[k])
> plot(x,y,type="l",lwd=2)
> x0=x[x>=1000000]
> polygon( c(x0,rev(x0)), c(y[x>=1000000],rep(0,length(x0))), col = "yellow",border=NA)
> lines(x,y,type="l",lwd=2)
> abline(v=1000000)

Le cas de la probabilité d’atteindre plus d’un million de visiteur le 1er janvier est alors
> k=which(u==as.Date("01/01/2011","%d/%m/%Y"))
> 1-pnorm(1000000,mean=forecast$pred[k],sd=forecast$se[k])
[1] 0.2604821

problème dual,

Dans un second temps, on peut envisager une autre approche, consistant à se demander quelle pourrait être la loi de la date du jour où le million de pages vues sera atteint,
> P=rep(NA,300)
> for(k in 1:300){
+ P[k]=1-pnorm(1000000,mean=forecast$pred[k],sd=forecast$se[k])}
> plot(max(D)+1:300,P,type="l",lwd=2)
> x=max(D)+1:300
> x0=x[x<=as.Date("01/01/2011","%d/%m/%Y")]
> polygon( c(x0,rev(x0)), c(P[x<=as.Date("01/01/2011","%d/%m/%Y")],rep(0,length(x0))), col = "yellow", border=NA)
> lines(max(D)+1:300,P,type="l",lwd=2)
> abline(v=as.Date("01/01/2011","%d/%m/%Y"))

La probabilité que cette date soit antérieure au 1er janvier est alors
> P[u==as.Date("01/01/2011","%d/%m/%Y")]
[1] 0.2604821
(ce qui, fort heureusement, correspond à la probabilité calculée par le problème primal). Bref, j’ai 1 chance sur 4 d’atteindre le million de pages vues avant la nouvelle année....

sensibilité à un changement de modèle

C’est bien joli tout ça, mais ces calculs sont largement soumis à la contrainte du choix de modèle que j’ai fait arbitrairement au début... on peut se demander si en changeant de modèle, les résultats changent sensiblement, ou pas. Au lieu de tenter un autre modèle ARIMA (voire SARIMA), j’ai préféré changer la série de référence... et me focaliser sur 2010 uniquement. D’un côté j’enlève les premières semaines où le niveau de fréquentation était très faible, de l’autre, je donne un poids très important aux vacances d’été, i.e. la période juin-août, pendant laquelle les internautes semblent moins sensibles à la modélisation économétrique.
Si l’on modélise la fréquentation pour 2010, seulement, on obtient

avec les distributions suivantes, tout d’abord pour la densité du nombre de visiteur atteint au 1er janvier 2011,

et pour la fonction de répartition de la date où sera atteint le million de pages vues,

soit une probabilité de l’ordre de 35%, dans les deux cas,
> P[u==as.Date("01/01/2011","%d/%m/%Y")]
[1] 0.3531648
> k=which(u==as.Date("01/01/2011","%d/%m/%Y"))
> 1-pnorm(1000000,mean=forecast2$pred[k],sd=forecast2$se[k])
[1] 0.3531648

On obtient des probabilités relativement proches avec les deux modèles, et j’aurais envie de croire que l’objectif est envisageable.... "challenge accepted" comme dirait mon ami Barney Stinson (oui, c’est mon ami). Reste juste à trouver des sujets qui attireront du monde en cette fin d’année...

Statistique de l'assurance STT6705V, partie 2

Arthur Charpentier — Mon, 06 Sep 2010 23:24:00 +0200

Après un long week end - fête du travail oblige - alors que les étudiants Rennais reprenaient le chemin de la fac, on va poursuivre le cours de statistique de l’assurance. Ce premier cours permettra de présenter un peu la théorie des modèles GLM, les modèles linéaires généralisés,

le but sera de survoler la théorie, en insistant sur leur utilisation en assurance. Les transparents seront mis en ligne ici mercredi, un peu avant le cours.
Pour les étudiants Rennais, j’avais envoyé un lien (ou plutôt fait envoyé, donc consulté vos mails) pour voir les 3 heures d’introduction qui avait été filmées la semaine dernière (cf ici). Pour des raisons de droits d’auteur, je crois que je ne suis pas autorisé à donner le lien (mais je vais me renseigner).

Je vous encourage à aller passer 3 heures devant votre ordinateur avant le cours de mercredi¹....

côté programmation, j’ai mis un peu de code R en ligne afin d’importer les bases que l’on utilisera au cours des premières session, et les mettre un peu en fome.....
> sinistreUdM <- read.table("http://perso.univ-rennes1.fr/arthur.charpentier/sinistreUdM.txt",
+ header=TRUE,sep=";")
> sinistres=sinistreUdM[sinistreUdM$garantie=="1RC",]
> nrow(sinistres)
> contratUdM <- read.table("http://perso.univ-rennes1.fr/arthur.charpentier/contratUdM.txt",
+ header=TRUE,sep=";")
A partir de ces deux bases, on peut récupérer le nombre de sinistres, par police,
> T=table(sinistres$nocontrat)
> T1=as.numeric(names(T))
> T2=as.numeric(T)
> nombre1 = data.frame(nocontrat=T1,nbre=T2)
> I = contratUdM$nocontrat%in%T1
> T1=contratUdM$nocontrat[I==FALSE]
> nombre2 = data.frame(nocontrat=T1,nbre=0)
> nombre=rbind(nombre1,nombre2)
.... et voilà, on peut commencer,
> base = merge(contratUdM,nombre)
> head(base)
nocontrat exposition zone puissance agevehicule ageconducteur bonus
1        27       0.87    C         7           0            56    50
2       115       0.72    D         5           0            45    50
3       121       0.05    C         6           0            37    55
4       142       0.90    C        10          10            42    50
5       155       0.12    C         7           0            59    50
6       186       0.83    C         5           0            75    50
marque carburant densite region nbre
1     12         D      93     13    0
2     12         E      54     13    0
3     12         D      11     13    0
4     12         D      93     13    0
5     12         E      73     13    0
6     12         E      42     13    0

¹ de toutes façons, la rentrée des séries télés n’a pas commencé.... il faudra attendre 2 semaines avant que How I Met Your Mother ne reprenne, ainsi que House, Dexter ou encore The Big Bang Theory . Bref, pas d’excuse (pour l’instant), il n’y a rien d’autre à la télé....

Statistique de l'assurance STT6705V, partie 1

Arthur Charpentier — Tue, 31 Aug 2010 04:05:00 +0200

Le premier cours commencera mercredi, pavillon Roger Gaudry, salle S116. Le premier cours ne sera pas transmis en visio-conférence (la rentrée des masters à Rennes n’ayant lieu que la semaine prochaine). Ce cours proposera une introduction au cours qui aura lieu tous les mercredis matins pendant la session d’automne. Les transparents du premier cours sont en ligne ici. Pour les prochains cours, je mettrais en avance les lignes le code que j’utiliserais.

Sinon pour le barème, il y aura deux mini-projets (25-25) de provisionnement et de mortalité, et un plus gros (50) de tarification. Le cours se basera sur R, mais pour les amateurs de SAS (ou autre), je suis ouvert....

Passer des cassettes aux mp3

Arthur Charpentier — Sun, 22 Aug 2010 22:25:00 +0200

Plus le temps passe, plus j’ai l’impression d’être un dinosaure.... Lorsque j’étais jeune, j’écoutais de la musique sur des cassettes¹. Je me souviens même avoir été enchanté par l’arrivée du premier walkman réversible: on n’était alors plus obligé de tourner la cassette au bout de 20 minutes. Quand soudain, à la fin des années 80 est arrivé le compact disc, et les lecteurs avec leur fonction Random, ou Shuffle.

J’aurais du me douter que je tournerais mal parce que je me souviens avoir joué avec la platine de mes parents, en essayant de voir si la première chanson pouvait passer en premier avec une lecture aléatoire, mais aussi la seconde pouvait passer juste après. Bon, avec le recul, je me rends compte que j’ai du beaucoup user la platine parce qu’avec un disque de 10 titres, il faut 90 (ou 100) essais (en moyenne) pour qu’une lecture aléatoire commence effectivement par 1 et 2. C’est d’ailleurs à ce moment que j’ai aussi découvert les tirages avec et sans remises: avec certaines platines, on ne pouvait pas obtenir deux fois le même titre, alors que d’autres oui.
Puis sont arrivés les lecteurs mp3, où on n’avait plus 10 titres tirés au hasard, mais 200 voire 1000. Et l’autre jour en courant, je me suis posé les mêmes questions que lorsque j’étais tout petit (oui, je suis un dinosaure, mais j’étais jeune à cette époque): si le lecteur propose une fonction random avec remise, combien de chansons dois-je attendre, en moyenne, avant de réentendre une chanson que j’avais déjà entendu ?

De la probabilité d’entendre deux fois une chanson ?

Pour faire un peu formel, supposons que mon lecteur contienne chansons.
La probabilité de n’avoir aucune répétition sur lectures (on suppose ) est

qui correspond aussi au nombre d’arrangements de éléments parmi , divisé par , soit

Si le ratio est faible (ou n est grand), on peut utiliser l’approximation et donc la probabilité s’approche par

soit tout simplement - toujours en cherchant une approximation - une probabilité de n’avoir aucune répétition de

Si on se donne une probabilité de n’avoir aucune répétition de , alors

Avec =200 chansons sur mon lecteur mp3, en jouant 17 titres, j’ai environ 1 chance sur 2 d’avoir au moins un titre joué deux fois. Et 17 chansons, ça correspond à ce que j’écoute en faisant mon jogging.... Bref, le temps de faire mon jogging, j’ai 1 chance sur 2 d’entendre un titre deux fois.
Mais que se passerait-il si je tenais le raisonnement inverse ?

Problème dual et jogging interminable....

Autrement dit, combien de temps vais-je courir si je décide de m’arrêter lorsque j’entends (pour la première fois) un titre pour la seconde fois ? En fait, Donald Knuth en parle un peu dans son livre the art of computer programming: on peut montrer que l’espérance du nombre de titres qu’il faudra jouer pour avoir une première répétition parmi n titres est de la forme

qui admet un développement asymptotique

soit, toujours avec 200 titres dans mon lecteur, une vingtaine de titres à jouer, en moyenne. Ce qui rallongera un peu mon jogging, finalement...
> n=200
> sqrt(2*n*log(2))
[1] 16.65109
> 1+sqrt(pi*n/2)-1/3+1/12*sqrt(pi/2*n)-4/135/n
[1] 19.8681
Formellement, le premier correspond à une médiane, alors que le second à une espérance. Ce qui explique la différence.
Notons que si mon jogging classique dure la durée de 17 titres, la probabilité d’avoir un jogging deux fois plus long est quand même d’une chance sur 10,
> k=31
> exp(-k*(k-1)/(2*n))
[1] 0.09778344
> sqrt(2*n*log(1/.01))
[1] 42.91932
autrement dit, j’ai une chance sur 100 (ce qui devrait donc arriver une fois dans l’année, en moyenne) de courir pendant 43 chansons... qui vont me paraître très très longues....
¹ je passais mon temps l’oreille sur la radio (car j’ai connu les premières radios libres, comme on disait alors), prêt à enregistrer dès qu’une super chanson passait.. je me souviens avoir guetter des jours durant le passage de you spin me around (ici) que mes élèves (qui n’ont pas connu cette époque) ou les cinéphiles connaissent probablement seulement sous la version right round (là) qui a abondamment samplé le thème. Oui, je suis un dinosaure..

Des étés pluvieux en Bretagne ? une réalité statistique...

Arthur Charpentier — Sat, 21 Aug 2010 15:32:00 +0200

Pour compléter le précédant billet (ici) on peut se demander en quoi la Bretagne est différente des autres régions françaises...

Nous avions vu ici le niveau de précipitation moyen, jour après jours pendant les mois d’été, en Bretagne. A Rennes. En revanche, à Paris on obtient la moyenne suivante,

que l’on peut comparer à Marseille,

ou encore à Strasbourg,

Sur la figure ci-dessous, on voit que la probabilité d’avoir de la pluie à Paris (au sens au moins 0.1 mm d’eau dans la journée, en trait gras bleu, au moins 2 mm d’eau dans la journée, en trait bleu) est supérieure à la probabilité d’avoir de la pluie à Rennes (respectivement en bleu clair gras, et en bleu clair fin)

On est certes très au dessus de Marseille,

mais très en dessous de Strasbourg,

Mais au delà des lois marginales, ces villes sont différentes de la Bretagne si l’on regarde les matrices de transition.

Transition d’un jour sur l’autre

Pour Rennes, si on regarde jour après jour, on obtient

		jour
		beau temps	pluie
jour	beau temps	1955	612	2567
jour	pluie	606	723	1329

ce qui donne les probabilités de transition suivantes,

		jour
		beau temps	pluie
jour	beau temps	76,15 %	23,85 %
jour	pluie	45,60 %	54,40 %

Pour Paris, la probabilité de transition jour après jour a la forme suivante

		jour
		beau temps	pluie
jour	beau temps	2689	959	3648
jour	pluie	946	1466	2412

ce qui donne les probabilités de transition suivantes,

		jour
		beau temps	pluie
jour	beau temps	73,71 %	26,29 %
jour	pluie	39,22 %	60,78 %

Pour Marseille, la probabilité de transition jour après jour a la forme suivante

		jour
		beau temps	pluie
jour	beau temps	2527	375	2902
jour	pluie	362	216	578

ce qui donne les probabilités de transition suivantes,

		jour
		beau temps	pluie
jour	beau temps	87,08 %	12,92 %
jour	pluie	62,63 %	37,37 %

Pour Strasbourg, la probabilité de transition jour après jour a la forme suivante

		jour
		beau temps	pluie
jour	beau temps	31	128	159
jour	pluie	132	1464	1596

ce qui donne les probabilités de transition suivantes,

		jour
		beau temps	pluie
jour	beau temps	19,50 %	80,50 %
jour	pluie	8,27 %	91,73 %

Transition d’une semaine sur l’autre

Si en revanche on regarde les matrices de transition semaine par semaine, on a des résultats assez différents. Une bonne semaine signifie aucun jour avec plus de 2 dm de pluie.
Pour Rennes, si on regarde semaine après semaine

		jour
		beau temps	pluie
jour	beau temps	379	25	404
jour	pluie	26	7	33

ce qui donne les probabilités de transition suivantes,

		jour
		beau temps	pluie
jour	beau temps	93,81 %	6,19 %
jour	pluie	78,79 %	21,21 %

Pour Paris, la probabilité de transition semaine après semaine a la forme suivante

		jour
		beau temps	pluie
jour	beau temps	576	46	622
jour	pluie	53	4	57

ce qui donne les probabilités de transition suivantes,

		jour
		beau temps	pluie
jour	beau temps	92,60 %	7,40 %
jour	pluie	92,98 %	7,02 %

Pour Marseille, la probabilité de transition semaine après semaine a la forme suivante

		jour
		beau temps	pluie
jour	beau temps	274	59	333
jour	pluie	47	9	56

ce qui donne les probabilités de transition suivantes,

		jour
		beau temps	pluie
jour	beau temps	82,28 %	17,72 %
jour	pluie	83,93 %	16,07 %

Pour Strasbourg, la probabilité de transition semaine après semaine a la forme suivante

		jour
		beau temps	pluie
jour	beau temps	1494	614	2018
jour	pluie	613	939	1552

ce qui donne les probabilités de transition suivantes,

		jour
		beau temps	pluie
jour	beau temps	70,87 %	29,13 %
jour	pluie	39,50 %	60,50 %

Les tests du chi deux, d’indépendance d’une semaine sur l’autre donnent

à Rennes, une statistique du chi-deux de 8,054, soit une p-value de 0,45%
à Paris, une statistique du chi-deux de 0,025, soit une p-value de 87,26%
à Marseille, une statistique du chi-deux de 0,012, soit une p-value de 91,24%
à Strasbourg, une statistique du chi-deux de 0,7649, soit une p-value de 38,18%

autrement dit l’hypothèse d’indépendance est acceptée partout, sauf à Rennes....

Moralité ?

De manière assez paradoxale, on prétend que la Bretagne a un temps changeant, et pour reprendre le titre du précédant billet, effectivement, en Bretagne, il peut faire beau plusieurs fois par jour. Mais sur le long terme, d’une semaine sur l’autre, le temps est au contraire très corrélé, contrairement aux autres régions. A Paris, Marseille ou Strasbourg, qu’il ait fait beau, ou qu’il ait plu la semaine précédente, cela n’apporte aucune information sur la probabilité d’avoir de la pluie la semaine où l’on vient en vacances.... Mais pas en Bretagne: manifestement, il existe donc des étés pourris, où il pourra pleuvoir toutes les semaines, et des étés superbes où il ne pleut jamais....

En Bretagne, il fait beau plusieurs fois par jour

Arthur Charpentier — Sat, 21 Aug 2010 09:27:00 +0200

bon, et à l’occasion il peut pleuvoir un peu.... Il peut donc être intéressant pour planifier un peu ses vacances de calculer la probabilité d’avoir de la pluie.

Les calculs ont été fait sur les données de pluviométrie à Rennes, en ligne sur le site eca&d, ici (données de qualité gratuites).

Probabilité d’avoir de la pluie, pendant les vacances d’été

Considérons ici la série du niveau de précipitation par jour, en par 0,1 mm,

et la probabilité d’avoir de la pluie dans la journée (régression logistique),

Bref, il peut pleuvoir à Rennes l’été. Mais ça n’aide pas vraiment pour planifier ses vacances. Car si ça se trouve, il y a des étés sans pluie, et des étés où il pleut sans cesse.

Dynamique et matrice de transition (par jour et par semaine)

Autrement dit, au lieu de regarder les lois marginales (comme la probabilité d’avoir de la pluie dans la journée), on peut s’intéresser à la dynamique de la série, modélisée sous la forme d’une chaîne de Markov. Si on regarde jour après jour, avec les 30 mois de juillet et août entre 1980 et 2009,

		jour
		beau temps	pluie
jour	beau temps	871	298	1169
jour	pluie	292	335	627

ce qui donne les probabilités de transition suivantes,

		jour
		beau temps	pluie
jour	beau temps	75,51 %	25,49 %
jour	pluie	46,57 %	53, 43 %

autrement, s’il fait beau aujourd’hui, on a 3 chances sur 4 d’avoir du beau temps demain.
Si on regarde semaine après semaine, où l’intérêt sont les semaines sans pluie,

		semaine
		beau temps	pluie
semaine	beau temps	13	26	39
semaine	pluie	23	140	163

avec là aussi les probabilités de transition suivantes,

		semaine
		beau temps	pluie
semaine	beau temps	33,33 %	66,67 %
semaine	pluie	14,11 %	85,89 %

Si on regarde semaine après semaine, où l’intérêt sont les semaines avec six jours sans pluie (on s’autorise une journée de pluie),

		semaine
		beau temps	pluie
semaine	beau temps	36	36	72
semaine	pluie	30	100	130

ce qui donne les probabilités de transition suivantes,

		semaine
		beau temps	pluie
semaine	beau temps	50,00 %	50,00 %
semaine	pluie	23,08 %	76,92 %

Bon, maintenant, si on pense qu’avoir 2 mm de pluie dans la journée, c’est juste un peu d’humidité dans l’air, les matrices de transitions sont sensiblement différentes, dans le cas où on s’autorise une journée dans la semaine avec un peu d’humidité pour parler de beau temps,

		semaine
		beau temps	pluie
semaine	beau temps	106	36	142
semaine	pluie	31	29	60

ce qui donne les probabilités de transition suivantes,

		semaine
		beau temps	pluie
semaine	beau temps	75,65 %	25,35 %
semaine	pluie	51,67 %	48,33 %

Autrement dit, on retrouve une matrice proche de celle que nous avions sur les données journalière: s’il a fait beau la semaine avant de venir en Bretagne, on a 3 chances sur 4 d’avoir du beau temps. En revanche, s’il a fait mauvais, on a une chance sur deux d’avoir beau temps la semaine suivante.
Si notre définition de pluie est encore plus laxiste (il faut qu’il y ait eu un déluge une journée dans la semaine, à savoir plus de 2 cm d’eau dans la journée), alors cette fois, on obtient,

		semaine
		beau temps	pluie
semaine	beau temps	171	12	183
semaine	pluie	14	5	130

ce qui donne les probabilités de transition suivantes,

		semaine
		beau temps	pluie
semaine	beau temps	93,44 %	6,56 %
semaine	pluie	73,68 %	26,32 %

Autrement dit, s’il y a eu au moins une très mauvaise journée la semaine passée, on a 1 chance sur 4 d’en avoir également une la semaine suivante. En revanche, s’il a fait beau tout le temps, on a presque 95% de chances d’avoir du beau temps la semaine suivante.
Moralité, quelle que soit la définition retenue pour définir le beau temps, le temps à Rennes n’est pas indépendant d’une semaine sur l’autre: il y a manifestement des étés pluvieux, et d’autre non.

Allez, on va fermer (temporairement...)

Arthur Charpentier — Sun, 08 Aug 2010 01:26:00 +0200

Vitesse et comportement dangereux au volant (partie 3)

Arthur Charpentier — Sat, 07 Aug 2010 02:42:00 +0200

Après plusieurs billets sur les accidents de la route (ici et là), parlons un peu de vitesse des véhicules. J’ai pu récupérer des données via l’inrets, d’observations près de Saclay en région parisienne, sur une 2x2 limitée à 110 km/h.

J’ai ainsi le passage de plusieurs millions de véhicules, sur quelques jours, incluant tous les véhicules. Mais si j’ai signé des papiers quant à la confidentialité des données, je dois pouvoir mettre des dessins construits à partir de ces données...

A quelle vitesse roulent les conducteurs parisiens ?

Comme pour les accidents, distinguons semaine et week end. En semaine, on retrouve les pics du matin et du soir, qui poussent les gens à ralentir. La courbe en noir est la vitesse moyenne des véhicules qui passent sur ce tronçon, et la courbe en bleu, la vitesse des 5% des véhicules les plus rapides,

Ces pics sont légèrement atténués le week end,

avec ici un trafic très dense tout l’après midi, entre 15 et 20 heures.

Quelle proportion de véhicules dépasse la limite autorisée ?

Comme pour les accidents, distinguons semaine et week end. En semaine, on retrouve les pics du matin et du soir, avec en gras la proportion des véhicules qui dépasse 110 km/h, et en trait fin, la proportion qui dépasse 115 km/h (soit 5% de plus),

avec des choses finalement assez proches le week end,

Autrement dit, aux heures de pointe, peu de monde dépasse la vitesse autorisée, mais vers midi, ou surtout la nuit, environ un véhicule sur quatre dépasse la vitesse autorisée (mais seulement un sur six la dépasse de plus de 5%).

Vitesse et respect des distances de sécurité

La base est très complète, puisque j’ai tous les véhicules qui sont passés près de Saclay, avec l’heure d’observation, la vitesse, ce qui permet d’obtenir la distance entre deux véhicules qui se suivent (en supposant que le premier véhicule maintienne sa vitesse constante).
Rappelons qu’il existe un lien théorique entre la vitesse d’un véhicule et la distance de freinage (expliqué ici ou là). Pour faire simple, la distance de freinage est la somme d’une distance parcourue avant de réagir , et d’une distance de freinage proprement dite ,

Si on suppose que la personne met une seconde à réagir, la distance de réaction est simplement la vitesse (exprimée en m/h), i.e.

Quant au freinage à prorement parler, rappelons que les équations qui charactérisent le mouvement d’un objet sont de la forme

Aussi, en mettant la deuxième au carré, et en combinant avec la troisième, on a

avec (on cherche à arrêter le véhicule) et , aussi, avec , on obtient quelque chose qui peut s’écrire

à condition d’avoir des vitesses en km/h. Aussi en agglomération, il faut 23 m pour s’arrêter. Et sur notre route à Saclay, pour une personne qui roule à 110 km/h, il faut compter 110 m (ce qui correspond au message véhiculé par la sécurité routière).
Sur nos données, on peut représenter en rouge les individus qui ne respectent pas la distance de sécurité.

La courbe bleue est la distance moyenne des véhicule en fonction de la distance qui les sépare du véhicule qui les précède,
Aussi, en faisant une régression logistique (lissée), on obtient que 60% des conducteurs qui roulent entre 50 et 110 km/h ne respectent pas les distances de sécurité.

Quand surviennent les accidents de la route ? (partie 2)

Arthur Charpentier — Fri, 06 Aug 2010 23:33:00 +0200

Je poursuis un peu mon étude des accidents de la route (ou pour être plus précis, des accidents corporels, ayant causé des blessés, et ayant fait l’objet d’un rapport de police).

J’avais parlé ici de la tendance de long terme, sur 6 ou 7 ans, montrant qu’il y avait une tendance à la baisse aussi bien du nombre de tués que du nombre de blessés grave. On avait observé que les cycles annuels avaient tendance à diminuer également. Mais qu’en est-il des cycles à court terme ?
En particulier, on peut se demander s’il n’existe pas des cycles dans la journée. En semaine, l’hiver (entre novembre et début février), on observe les tendances suivantes,

(les points sont les moyennes brutes, par tranche d’une demi heure, et la courbe est un lissage de ces points). On distinguera les années 2002-2004 (en rouge) et 2004-2007 (en bleu) Si l’on regarde les semaines l’été (entre mai et début août), on a

Notons que si l’on compte non plus les accidents, mais les blessés graves ou les tués (il y a généralement plusieurs victimes dans ces accidents), on observe les cycles suivants,

pour l’hiver, alors que l’été,

Autrement dit, les tendances sont les mêmes: les accidents du matin ne sont pas moins meutriers que ceux du soir... Les week-ends, les tendances sont assez différentes entre l’hiver et l’été, avec respectivement

Si l’on regarde les blessés graves et les tués, on a les cycles suivants,

Nombre de morts sur la route en baisse... et ? (partie 1)

Arthur Charpentier — Thu, 05 Aug 2010 00:22:00 +0200

Tout le monde semble se réjouir de la baisse du nombre de morts sur la route depuis quelques années (longtemps attribué à la présence de radars automatiques, mais qu’on pourrait aussi associer à l’instauration du permis à point). Et c’est tant mieux. Mais sans vouloir jouer les cyniques, et en tous les cas c’est une statistique qui intéresse probablement davantage les assureurs, il peut sembler intéressant de s’intéresser aussi aux accidents corporels graves....

Le principal soucis méthodologique est qu’au 1er janvier 2005, la définition d’accidents graves ou encore de décès a changé. Avant 2005, un tué, était une victime décédée sur le coup ou dans les six jours qui suivent l’accident. A partir de 2005, les 6 jours passent à 30 jours. De même, un blessé grave était un blessé dont l’état nécessitait plus de six jours d’hospitalisation. Après 2005, cette classe disparait, et seuls restent les blessés qui ont passé plus de 24 heures à l’hôpital.... Bref, on est un peu bloqué pour vraiment analyser les tendances sur une longue période...

Nombre de morts sur la route

Le graphique suivant montre le nombre quotidien de morts sur les routes, passant d’une vingtaine en 2002 à une douzaine (en 2008).

On notera qu’il n’y a pas eu de réelle rupture en 2005, malgré la nouvelle définition de décès.

Nombre de blessés graves

En revanche, les statistiques sur les blessés (graves) ont connu une réelle rupture en 2005.

On peut essayer de corriger afin de rendre les séries comparables (au moins en terme d’ordre de grandeur),

Au delà de la tendance (et de la rupture de la forte baisse observée en 2002-2004), on notera que le cycle annuel est plus faible. De la même manière que les cycles dans les naissances semblent avoir diminué (comme je l’évoquais ici), on notera que le cycle dans les accidents graves (mais aussi dans une moindre mesure les décès) a fortement diminué...

Propriétaire ou locataire ?

Arthur Charpentier — Wed, 04 Aug 2010 22:29:00 +0200

Il existe une mythologie aujourd’hui bien ancrée, consistant à penser qu’il faut être propriétaire de son logement. C’est d’ailleurs le message que l’État (c’est à dire nous tous) véhicule sur différentes affiches, comme celle à droite. Tout le monde a déjà entendu l’adage prétendant que "payer un loyer, c’est jeter de l’argent par les fenêtres"....mais personnellement, j’ai l’impression que rembourser un crédit, c’est aussi ’"jeter de l’argent par les fenêtres". Dans le premier cas, c’est le propriétaire qui est derrière la fenêtre, alors que dans le second cas c’est une banque (mais aussi un notaire, et l’état qui récupère beaucoup d’impôt).
Visiblement, l’État semble avoir pris le parti qu’enrichir les propriétaires, c’est mal, alors qu’enrichir son banquier c’est bien... Peut être est-ce mieux dans l’intérêt collectif (car c’est le but ultime de l’État me semble-t-il) d’être tous propriétaires ? Mais qu’en est-il de l’intérêt individuel ? Ai-je vraiment intérêt à préférer être propriétaire plutôt que locataire ?
On peut essayer de faire un petit modèle simple, voire simpliste... Considérons quelqu’un possédant un apport personnel de 100 000 € (je mets des sommes rondes pour simplifier le modèle), et souhaitant loger dans une maison qui vaut 300 000 €. Il a deux possibilités,

devenir propriétaire du bien
être locataire du bien (on suppose que le choix est possible pour le même bien)

Dans le premier cas, on va supposer que qu’il est possible d’avoir un crédit sur 20 ans pour un taux de 5,5%. On va supposer qu’il y a 10% de frais d’acquisition, entre l’agent immobilier et le notaire. On oubliera les charges annuelles en tant que propriétaire, et les éventuels crédit d’impôt sur les intérêts.
Dans le second cas, il paye un loyer de 1 000 € par mois, éventuellement revalorisé (+2% par an par exemple), mais il a la possibilité d’épargner (en plus de capitaliser ce qui constituait l’apport initial s’il achetait).
En fait, dans les deux cas, il épargne. On peut partir du fait qu’il a 2 000 € allouer au logement (ce montant peut également être revalorisé). Dans le second cas, il peut épargner 1 000 € par mois, et dans le premier, la différence entre les mensualités qu’il doit à la banque et les 2 000 €. On suppose que les taux de placement sont à 4,5%.
Au bout de 15 ans, la maison vaut 500 000 €. Notons que la revalorisation à 2% est plus faible que la croissance de la valeur du bien (passer de 300 à 500 en 20 ans correspond à une valorisation annuelle de 2,5%). Le propriétaire a-t-il vraiment fait une bonne affaire ?
> capital=100000
> revenu =2000
> loyer =1000
> i1     =.055
> i2     =.045
> i3     =.02
> T      =20
> maison1=300000
> maison2=500000
> mensualite=(maison1*1.1-capital)/sum((1+i1)^(-(1:(T*12))/12))
> mensualite
[1] 1564.789
> (1+i3)^T
[1] 1.485947
> locataire    = capital*(1+i2)^T+sum((1+i3)^((1:(T*12))/12)*(revenu-loyer)*(1+i2)^((1:(T*12))/12))
> proprietaire = maison2+sum((revenu-mensualite)*(1+i3)^((1:(T*12))/12)*(1+i1)^((1:(T*12))/12))
> locataire
[1] 728274.5
> proprietaire
[1] 738245.5
> (proprietaire-locataire)/proprietaire
[1] 0.01350635
Autrement dit pour être propriétaire, les mensualités sont de 1 564 € par mois (constantes), alors que le locataire aurait eu un loyer de 1 000 € initialement, revalorisé de 2% par an, soit 1 485 € par mois sur la fin. Bref, au final, le locataire, qui a mis 1 000 € par mois (revalorisés là aussi à 2%) se retrouve à la tête de 728 milliers € alors que le propriétaire a vu son bien se valoriser (et atteindre 500 milliers €), mais il a également pu épargner un peu. Bref, il se retrouve à la tête d’un patrimoine de 738 milliers €. Ce qui est comparable.... Moralité, on jette autant d’argent par les fenêtres dans les deux cas, mais du point de vue de la personne qui loge dans l’appartement, ça se vaut...
Mais au delà des montants, on peut surtout étudier l’impact des différents paramètres....
Le graphique ci-dessous montre l’impact du capital initial: s’il est trop faible, être locataire est beaucoup plus intéressant,

(les ordonnées positives signifie que la richesse en tant que propriétaire excède celle que l’on aurait en tant que locataire). Autrement dit, avec un faible capital initial, il sera coûteux de devenir propriétaire. Ici, il faut disposer d’au moins 1/3 du prix de la maison en apport. Mais grosso modo, plus on est riche, plus on peut avoir intérêt à devenir propriétaire. On peut aussi regarder l’impact du taux d’emprunt,

(où l’on retrouve que si le taux d’emprunt est trop élevé, je n’ai aucun intérêt à emprunter) ou encore du revenu que l’on souhaite allouer au logement,

Pour comparer l’impact de la durée, on va supposer que la valorisation du loyer et du salaire reste à 2%, mais que la valorisation du bien immobilier est de 2,5% par an.

Certes, le locataire épargne moins sur une durée courte, mais le locataire a de telles mensualités à payer qu’il vaut mieux épargner. Bon, bien sûr la difficulté est de supposer qu’on peut faire varier ces paramètres indépendamment les uns des autres, mais j’ai l’impression que cela permet de mieux comprendre qui peut être intéressé pour devenir propriétaire.... en l’occurrence ceux qui ont le temps, et ceux qui ont de l’argent....
Bref, je laisse ceux qui le souhaite améliorer le modèle car il est ici très (trop ?) simpliste... même s’il laisse à penser que le choix d’acheter ou de louer n’est pas aussi simple qu’il y paraît, loin de là.....

Importer une base SAS sous R

Arthur Charpentier — Tue, 03 Aug 2010 15:30:00 +0200

Depuis pas mal de temps, je reçois des bases envoyées par différentes personnes, pour faire des stats, et si la plupart ont la bonne idée de me les envoyer directement en format csv, malheureusement, je reçois parfois des bases sas dont je ne sais trop que faire (car je n’ai pas sas).
La library(foreign) de R propose d’importer des bases SAS qui sont au format xport (ici). Mais généralement, les données SAS sont au format sas7bdat. Pendant l’année, j’avais l’habitude d’aller squatter la salle informatique qui dispose de SAS, pour ouvrir SAS et exporter les données en csv. Mais avec la fermeture estivale de la faculté, j’ai été un peu bloqué. J’ai alors découvert la version online de SAS sur l’ent de l’université. Malheureusement, il faut que les données soient dans un répertoire sur le réseau, ce qui nécessite au préalable d’envoyer des données sur le réseau au lieu de pointer sur un répertoire local (ce qui peut être loin, voir impossible si la base est trop grosse). Fort heureusement, l’autre jour, j’ai découvert,

On peut en effet télécharger gratuitement sur le site de www.sas.com un lecteur de base, qui permet ensuite d’exporter la base en csv... et ensuite de l’importer facilement sous R... Trop facile !

Tempête et inondation: 1,5 milliards d'euros pour Xynthia

Arthur Charpentier — Fri, 30 Jul 2010 21:43:00 +0200

Les premiers rapports sur Xynthia, ses conséquences mais aussi ses causes commencent à être diffusés, comme ici par exemple. Pour rappels, dans la nuit du 27 au 28 février, la tempête Xynthia a balayé plusieurs pays européens, causant de très nombreux dégâts en France, principalement en Charente Maritime et en Vendée. A l’ile de Ré, des rafales de 160 km/h ont été mesurées. Mais les effets de la tempête ont été décuplés par la présence de polders, l’océan reprenant alors en quelques heures des territoires que l’homme croyait s’être approprié depuis quelques siècles. Comble de malchance, cette tempête est survenu lors des grandes marées, ce qui a de nouveau accentué les effets : afin de compenser la baisse de pression atmosphérique, le niveau des océans monte davantage, créant ainsi une surcote, et causant des inondations dès que le niveau des digues est atteint.

Ces inondations ont causé environ 420 000 sinistres, pour un coût de 700 millions € (soit un coût moyen par sinistre de l‘ordre de 1700 €). 22 000 sinistres ont été des dommages matériels de véhicules (pour un coût total de l’ordre de 35 millions €), 69 000 sinistres ont été des dommages de professionnels (pour 255 millions €) et près de 330 000 sinistres ont été des dommages à des particuliers, pour mon montant total de 425 millions d’euros. Le risque inondation est couvert par le régime dit cat nat, instauré par la loi du 13 juillet 1982, qui est obligatoire dans tout contrat de dommages aux biens, et qui couvre tous les dommages causées par une « intensité anormale d’un agent naturel », excepté les tempêtes. Si une commune touchée est classée comme été victime d’une catastrophe naturelle, les habitants peuvent alors faire jouer la garantie. Toutefois, si la commune n’a pas fait de Plan de Prévention des Risques,. Le montant des primes d’assurance affectée par ces catastrophes est fixé par l’Etat, ou plutôt la part de la prime habitation (pour les particulier) qui sera affectée à la couverture du risque, en l‘occurrence 12%: peu importe le risque réel, l’Etat impose un mécanisme de solidarité entre les assurés. Les assureurs touchent ainsi environ 1,3 milliards € au titre de l’ensemble des risques dit cat nat (ce qui inclus la sécheresse ou les mouvements de terrain, en plus des inondations). Les assureurs peuvent également souscrire un contrat de réassurance auprès de la Caisse Centrale de Réassurance (CCR, possédant une garantie illimitée de l‘Etat) en quote-part, conduisant à un partage pour moitié entre l’assureur et le réassurance. La CCR offre de plus une clause de limitation des pertes, dans le cas où un sinistre couterait plus du double des primes encaissées. Compte tenu du rôle central joué par la CCR, c’est cette dernière qui fournit généralement les statistiques sur les évènements naturels.
Les tempêtes ont fait l’objet d’une extension obligatoire, en juin 1990, dans tous les contrats dommage, avec une couverture qui incluse généralement la grêle et le poids de la neige. La tempête en elle même a causé environ 35 000 sinistres, pour un coût de 700 millions € (soit un coût moyen par sinistre de l‘ordre de 20 000 €). Environ 10 500 de ces sinistres sont liés à des couvertures de véhicules, avec un coût moyen de l’ordre de 5 000 €. 5 500 ont causé des dommages à des professionnels, pour un coût total de l’ordre de 250 millions €, alors que 400 millions ont été affecté à 19 000 sinistres de particuliers. Cette garantie est certes elle aussi obligatoire, mais l’Etat ne propose plus d’intervenir, et ce sont les réassureurs privés qui permettent aux assureurs de se couvrir contre les très grosses tempêtes (on aura en mémoire celles de décembre 1999, dont le montant réévalué atteint les 10 milliards €). Notons que pour couvrir le risque tempêtes, les assureurs collectent environ 1 milliard € par an. Dans le cas de la garantie tempête, la principale source de données est la Fédération Française des Sociétés d’Assurance, qui collecte les données des différentes mutuelles et compagnies d’assurance.
Cet évènement a marqué les esprits non seulement à cause du nombre élevé de victimes, mais aussi car il rappelle deux évènements très similaires. Le premier est survenu dans la nuit du 31 janvier au 1er février 1953, au Pays Bas, où une tempête a traversé le pays à quelques jours des grandes marées (ici). Les digues subirent plus de 400 brèches, et plus de 500 km de digues furent à reconstruire. Mais surtout, 1 836 personnes trouvèrent la mort, ce qui, ramené à la population des Pays-Bas de l’époque (de l’ordre de 10 millions d’habitants), correspondrait – aujourd’hui – à une catastrophe causant la mort de plus de 10 000 personnes en France (ou plus de 60 000 aux Etats Unis). Le second est l’ouragan Katrina qui a détruit la Nouvelle Orléans le 28 août 2005, où le raz de marée qui est survenu au même moment a causé davantage de dégâts que les rafales de vent (pourtant très fort). Ce risque avait été évoqué par des ingénieurs quelques mois plus tôt.
Le parallèle avec Katrina est d’autant plus saisissant que, comme à la Nouvelle Orléans, on retrouve des rapports d’experts qui évoquent la catastrophe quelques mois plus tôt. En particulier, l’idée d’avoir un « Plan de Prévention des Risques de submersion marine dans l’estuaire du Lay sur les communes de La Faute-sur-Mer et de l’Aiguillon-sur-Mer, où la conjonction de deux phénomènes, de crue dans l’estuaire du Lay et de submersion marine, pourrait avoir un impact très important sur les zones densifiées à l’arrière d’un réseau de digues vieillissant. » (ici). A La Faute-sur-Mer, « la rupture des digues sur ce secteur engendrerait des dégâts majeurs aux biens et aux personnes », ajoute le rapport en citant l’exemple des tempêtes Lothar et Martin qui « ont démontré que les zones côtières pouvaient être submergées par la mer, notamment en zone estuarienne ».

(le graphique a été volé sur le site de 20min, ici). En effet, depuis les tempêtes de décembre 99, les tempêtes majeures en Europe sont jugées « rares », mais plus « exceptionnelles ». On a même vu des cyclones tropicaux se diriger vers l’Europe (au lieu de se diriger vers les Amériques) comme la tempête Grace en octobre 2009 (ici). Autrement dit le risque existe. Et comme toujours quand on fait face à un risque, deux solutions se posent à nous : accepter de prendre le risque, quitte éventuellement à transférer la conséquence financière à un tiers (un assureur le plus souvent), ou tout mettre en œuvre pour le réduire (on parlera alors de prévention). Mais encore faut-il toutefois avoir conscience de l’existence de ces risques, en l’occurrence de risques à long terme. Et ce n’est pas essentiel si l’on mène des politiques à court terme.

L'équipe de France de foot devrait jouer seulement en semaine

Arthur Charpentier — Sat, 24 Jul 2010 22:48:00 +0200

Il y a quelques semaines, le service de centralisation des statistiques de l’assurance accidents LAA, par l’intermédiaire de Stefan Scholz Odermatt, a publié une étude sur la survenance d’accident automobile les jours de match de foot.

Selon l’étude (mentionnée ici ou là) le nombre d’accidents les jours où il y a match est supérieur aux jours sans match,

Ayant réussi à récupérer une base presque exhaustive d’accidents corporels en France, je peux faire la même étude, entre janvier 2002 et décembre 2007. Les conclusions sont moins flagrantes que sur le cas suisse, en particulier les jours de semaine, où manifestement, personne ne s’intéresse aux matchs (où en tous les cas pas assez pour être moins vigilent au volant).

En revanche, le week end, l’effet devient plus prononcé

avec davantage d’accidents dans la journée où il y a un match de l’équipe de France.
Moralité, l’équipe de France de foot devrait jouer seulement en semaine, à chaque match, cela éviterait une vingtaine d’accidents de la route à chaque fois (et je ne parle que d’accidents corporels, ayant causé des blessures et ayant donné lieu à un constat de police) ! Comme le disait un copain à qui je racontais ça, "l’équipe de France devrait jouer, ça serait un début".... mais c’est un autre sujet.

Estimateur(s) sans biais ?

Arthur Charpentier — Fri, 23 Jul 2010 23:17:00 +0200

En fin de semaine j’étais à Lyon à l’école d’été organisée par Stéphane, pour parler bootstrap et estimation de l’incertitude associée à l’estimation du montant des réserves (pour faire simple, les slides sont en ligne ici). Je ferais un billet plus détaillé (avec les codes pour répondre à la demande que m’ont faites plusieurs personnes). A la fin de l’exposé, on discutait des estimateurs sans biais avec quelqu’un qui se plaignait du biais à distance fini de certains estimateurs, et qui me demandait "vous ne pouvez pas trouver un estimateur sans biais du msepc ? ". Si à la rigueur la recherche des propriétés d’estimateur peut être un passe temps occasionnel, je n’ai pas eu la présence d’esprit de lui faire remarquer qu’il existe des paramètres dont il est impossible de construire un estimateur sans biais... étonnant, non ?

Prenons la loi de Poisson de paramètre . Alors

Supposons que l’on cherche à estimer . S’il existe une statistique qui soit un estimateur sans biais de , alors

Or est une fonction de , donc

Cela signifie tout simplement que

Ça c’est bien gênant parce que la fonction n’admet pas de développement en série entière (je renvoie à des cours de prépa, ou là par exemple). Moralité, il existe des grandeurs dont on ne peut pas construire un estimateur sans biais.... Qui l’eut cru ?

Bonus-malus et non déclaration de sinistres

Arthur Charpentier — Fri, 23 Jul 2010 00:34:00 +0200

Le mécanisme bonus-malus a de très nombreuses vertus, dont celui de renforcer la solidarité. Mais il semble qu’il incite aussi à ne pas déclarer certains petits sinistres à son assureur. C’est ce que nous allons montrer ici.
La distribution des personnes ayant un bonus de 50% (ce qui est le niveau le plus bas que l’on puisse atteindre, en théorie) a la forme suivante, en fonction de l’âge du conducteur,
> library(splines) ; library(pscl)
> reg <- glm((bonus==50)~bs(ageconducteur),data=sinistres,family=binomial)
> base <- data.frame(ageconducteur=seq(18,80))
> y=predict(reg,newdata=base)
> plot(seq(18,80),y)

Comme je le notais déjà ici, une personne qui a un niveau de bonus bas peut être incitée à ne pas déclarer un petit sinistre à son assureur (et de proposer un arrangement à l’amiable avec la contrepartie).
Cette sur-représentation des 0 dans la base pour les très bas niveaux de bonus peut être prise en compte à l’aide des modèles dits zero-inflated.
Classiquement, on supposait que

dans le cas d’un modèle de Poisson. On va supposer ici que l’assuré peut décider de ne pas déclarer certains sinistres. Autrement dit

et pour

On peut considérer un modèle logistique pour modéliser cette probabilité de non-déclaration,

alors que pour le modèle de Poisson

En fait, si l’on suppose que l’impact d’une variable n’est pas linéaire, on peut introduire des splines pour estimer la transformation optimale,
> library(splines) ; library(pscl)
> reg=zeroinfl(nombre~bs(ageconducteur,df=4) | bs(ageconducteur), data = nombre,
+ dist = "poisson",link="logit",offset=log(exposition))
> age=seq(18,80)
> DT=data.frame(ageconducteur=age,exposition=1)
> Y=predict(reg,newdata=DT,type="zero")
> plot(age,Y)
Ce qui permet de ne faire une prédiction que sur la composante d’inflation zero. Sur une base de données sur laquelle je devrais revenir à la rentrée, on obtient la tendance suivante,

Malheureusement, l’interprétation est plus délicate, car avec une régression binomiale négative, qui autorise plus de variance, on obtient des ordres de grandeur très différents
> reg=zeroinfl(nombre~bs(ageconducteur,df=4) | bs(ageconducteur), data = nombre,
+ dist = "negbin",link="logit",offset=log(exposition))

(on oubliera le comportement de bord pour les âges élevés, peu de personnes âgées appartenant à ce portefeuille, les résultats sont peu robustes à droite).
On retrouve toutefois une fonction croissante (au moins entre 20 et 60 ans), ce qui peut être relié avec la distribution du bonus en fonction de l’âge: plus on est âgé, plus on a de chance d’avoir un très bon bonus, et plus on a de chances de ne pas déclarer un sinistre à son assureur...
D’ailleurs, si on fait la régression directement sur le niveau de bonus, et plus sur l’âge, on a l’impact suivant
> reg=zeroinfl(nombre~bs(bonus,df=4) | bs(bonus), data = nombre,
+ dist = "poisson",link="logit",offset=log(exposition))

Moralité, c’est bien le niveau élevé de bonus qui incite les assurés à ne pas déclarer de sinistres à leurs assureurs (et pas forcément un effet d’âge que l’on pourrait associer à de l’Alzheimer).

Davantage de régularité dans les naissances ?

Arthur Charpentier — Tue, 20 Jul 2010 00:50:00 +0200

L’autre jour, l’INED avait fait beaucoup de bruit en publiant une étude expliquant que le pic des naissance classiquement observé en mai semblait avoir disparu (ici ou là). On va essayer de vérifier ce point, à partir de toutes les naissances observées en France entre 1968 et 2005. Ca sera le billet sexe de l’été...

Le nombre quotidien de naissances en France

La série brute ressemble à ça,

On pourrait à la rigueur extraire les cycles directement sur cette série, mais ça ne serait pas forcément convainquant (même si visuellement on retrouve ce qui était évoqué),

mais la première chose que l’on observe est que la différence entre le nombre de naissances le week end (en rouge) et en semaine (en bleu) s’est fortement creusée dans les années 80-90.

On peut regarder le ratio du nombre de naissance le week end, par rapport aux jours de la semaine précédant le week end,

et on peut alors "corriger" ou "lisser" la série, de manière à maintenir la moyenne hebdomadaire, mais en transférant quelques naissances de la semaine au week end.

Sur cette série lissée, on peut extraire une tendance de long terme.

C’est sur le "bruit" autour de cette tendance que l’on peut essayer d’étudier un éventuel cycle annuel. Notons que l’on retrouve les tendances sur la série brute

Série corrigée de la tendance de long terme

La série suivante permet de voir l’évolution du "cycle" (série brute à laquelle on retranche la tendance observée ci-dessus) à la fin des années 60 (en rouge), des années 70 (en orange), des années 80 (en jaune), voire des années 90 (en presque blanc),

Sur les années 68-78, on a la tendance suivante

alors que pour les années 95-2005, cela donne

Au final, les deux tendances sont les suivantes, avec en rouge les années 70 et en bleu les années 2000.

Effectivement, le pic classiquement observé au mois de mai semble avoir disparu, et la série aujourd’hui semble beaucoup plus lisse. Et comme l’INED, on note qu’il y a effectivement un mini-pic au mois de septembre.

Des naissances au comportement sexuel

Tout comme l’INED, on peut relier ça à l’activité sexuelle des français, que l’on peut visualiser sur le graphique ci-dessous (en décalant tout simplement de 9 mois),

Bref, l’activité sexuelle semblait s’intensifier pendant les vacances scolaires, en particulier début août (mais aussi plus généralement pendant les vacances scolaires, que ce soit Pâques, les vacances d’été, ou les vacances de Noël d’ailleurs). Ou peut être pourrait-on conclure à la rigueur que l’été, on sort moins couvert.... car on observe les comportements sexuels associés à une naissance 9 mois plus tard (ce qui ne doit pas correspondre à l’unique activité sexuelle).
Les français semblent avoir intensifié leur activité sexuelle pendant une période plus fraîche (novembre-mars), au détriment de la période estivale. Un effet du réchauffement climatique (il fait trop chaud pour batifoler l’été) ? Un effet des 35 heures et des RTT (et des vacances qui sont plus réparties dans l’année qu’avant) ?

L’été s’ra chaud ?

En fait, au lieu de penser en terme de saison, on peut penser en terme de température (j’ai pris la température quotidienne minimale à Paris, c’est à dire la nuit... ).

c’est à dire en zoomant

On peut aussi chercher à étudier davantage ce qui se passait il y a 30 ans et maintenant (comme dans la partie précédente). Il y a 30 ans (courbe rouge) on observait un (très) fort pic estival. Ce pic a disparu aujourd’hui (courbe bleue), avec un comportement comparable, qu’il fasse 0 ou 20 degrés dehors.

Quid de l’impact des bourses ?

Mais la recherche de la variable qui pourrait expliquer cette saisonnalité ou ces pics est sans fin. Beaucoup de monde évoque la bonne santé économique pour expliquer les booms de natalité.... ça marche peut être à long terme, mais à court terme ?
En fait, si on croisse l’activité sexuelle avec le niveau de la bourse, on peut presque retrouver ce genre de comportement,

une bourse en très forte hausse semble provoquer un (léger) pic de natalité 9 mois plus tard, alors qu’un effondrement de la bourse se traduit plutôt par une chute (légère certes) 9 mois plus tard...

Qu'est ce que la modélisation économétrique (4) ?

Arthur Charpentier — Fri, 16 Jul 2010 17:31:00 +0200

Suite à une discussion au téléphone l’autre jour avec Elsa, j’ai voulu rajouter un petit billet - un peu technique - sur la multicolinéarité. Ou comment interpréter le signe (et la valeur) d’un coefficient dans une régression multiple.

Interpréter une régression avec des variables explicatives corrélées

Assez souvent en économétrie, les variables explicatives peuvent être corrélées. Ce qui peut introduire des erreurs d’interprétation. Considérons le petit modèle suivant, où l’on dispose de deux variables explicatives.
> library(mnormt)
> set.seed(1)
> Z=rmnorm(200,c(0,0),matrix(c(1,-.8,-.8,1),2,2))
> Y=1+2*Z[,1]+.5*Z[,2]+rnorm(200)
Si l’on se contente d’analyses univariées, on notera que semble décroitre avec .
> summary(lm(Y~Z[,2]))
Call:
lm(formula = Y ~ Z[, 2])
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   1.0229     0.1132   9.038   <2e-16 ***
Z[, 2]       -1.0717     0.1171 -9.155   <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.6 on 198 degrees of freedom
Multiple R-squared: 0.2974,     Adjusted R-squared: 0.2938
F-statistic: 83.81 on 1 and 198 DF, p-value: < 2.2e-16

Pourtant, si l’on contrôle à l’aide la variable , on s’aperçoit que l’impact de sur va plutôt dans l’autre sens, avec un signe positif, et significativement positif,
> summary(lm(Y~Z[,1]+Z[,2]))
Call:
lm(formula = Y ~ Z[, 1] + Z[, 2])
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.95647    0.07597 12.591 < 2e-16 ***
Z[, 1]       2.04287    0.13082 15.615 < 2e-16 ***
Z[, 2]       0.45621    0.12541   3.638 0.000351 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.073 on 197 degrees of freedom
Multiple R-squared: 0.686,      Adjusted R-squared: 0.6828
F-statistic: 215.2 on 2 and 197 DF, p-value: < 2.2e-16
Graphiquement, on a le modele suivant

autrement dit, croit fortement avec , mais à fixé, on note que est croissant avec . L’intuition derrière est qu’il faut regarder l’impact de sur pour des individus qui se ressemblent - ou identique - c’est à dire avec des valeurs de proches. Si l’on considère le sous-échantillon des petites valeurs de , on obtient
> I=Z[,1]<(-1)
> plot(Z[,2],Y,xlab="",ylab="")
> abline(lm(Y[I]~Z[I,2]),col="red",lwd=2)

alors que pour les grandes valeurs de ,

Bref, l’impact négatif semble ici positif. En considérant un échantillon plus grand (en simulant 2000 valeurs au lieu de 200), on obtient des choses ressemblant à ça,

Interprétation sur le problème des tailles des classes

Si l’on revient au sujet de l’impact de la taille de la classe (dans les écoles) sur les résultats scolaires, on pourrait être dans la même configuration. Sans variable de contrôle, la taille de la classe impact de manière positive les résultats scolaires: plus la classe est nombreuse, meilleures sont les résultats. Or comme le notaient Angrist et Lavy (ici), la taille de la classe est très corrélée à des variables socio-économiques. Si l’on arrive à trouver de bonnes variables de contrôle, on peut parfaitement obtenir un impact négatif, comme dans l’exemple simulé ci-dessus...

Tails of copulas, une lecture graphique

Arthur Charpentier — Thu, 08 Jul 2010 04:26:00 +0200

Suite à une formation que je faisais en fin de semaine à Brest (les slides sont ici et là), je voulais revenir sur les histoires de tails of copulas, pour reprendre le titre de l’article (ici) de Gary Venter (et qui correspond à des choses que j’avais pu présenter il y a quelques années à Berlin, les slides étant en ligne ici).

Quantifier la dépendance de queue

L’idée est de noter qu’il est noter qu’il existe deux manières de quantifier la dépendance de queue. La première est liée à l’approche de Joe (1990, ici, ou 1997 pour le livre), qui a introduit un (strong) tail dependence index. Par exemple pour la queue inférieure,

soit

La seconde est liée à une idée que l’on retrouve dans les travaux de Janet Heffernan, Stuart Coles ou Jonathan Tawn. L’intuition est la suivante (on peut la retrouver en ligne ici). Si et ont la même loi et que l’on suppose les variables indépendantes, alors

En revanche, si les variables sont comonotones (c’est à dire égales comme on suppose les lois identiques),

Aussi, on peut supposer qu’il existe un indice tel que

Le soucis est que le cas d’indépendance correspond à =2, alors que le cas de dépendance forte correspond au cas =1. Il est alors usuel de faire une transformation affine pour se ramener sur [0,1], et que la force de la dépendance soit croissante avec , e.g.

Posons alors

qui pourra être interprété comme un (weak) tail dependence index.
Bref, ces deux mesures donnent de l’information sur le comportement dans les queues de distribution.

Les fonctions de concentration dans les queues

L’idée est de noter qu’il est possible d’étudier ces fonctions afin de mieux comprendre le comportement dans les queues. En s’inspirant de Gary Venter, on peut définir

pour étudier le comportement dans la queue inférieure, et

pour la queue supérieure,où est la copule de survie associée à , au sens où

Cet outil permettra de modéliser la dépendance forte. On peut également poser, afin d’étudier la dépendance faible,

Application statistique

L’idée est de noter qu’il est facile d’estimer ces fonctions. Ces outils peuvent être utiles pour mieux comprendre le comportement dans les queues.
Par exemple pour une copule Gaussienne de corrélation 0,5, on a la forme théorique suivante pour les fonctions de concentration (au sens fort)

Statistiquement, il est possible d’estimer ces quantités en comptant simplement le nombre d’observations dans le coin inférieur gauche, ou le coin supérieur droit. Si on dispose d’un échantillon, on peut alors regarder ce que donnent les versions

Pour un échantillon de taille n=500, on obtient les intervalles de confiance à 90% de la forme suivante,

Le code R ressemble à ça
> library(evd); data(lossalae)
> cor(lossalae,method="spearman")
Loss ALAE
Loss 1.000000 0.451872
ALAE 0.451872 1.000000
avec le code suivant pour la version empirique,
> library(evd); data(lossalae)
> z=seq(0,.5,by=.001)
> U=rank(v[,1])/(nrow(v)+1)
> V=rank(v[,2])/(nrow(v)+1)
> Lemp=rep(NA,length(z))
> Remp=rep(NA,length(z))
> for(i in 1:length(z)){
+ Lemp[i]=sum((U<=z[i])&(V<=z[i]))/sum(U<=z[i])
+ Remp[i]=sum((U>=1-z[i])&(V>=1-z[i]))/sum(U<=z[i])
+ }
et pour la version théorique,
> Lg=(pcopula(copclayton,cbind(z,z)))/(z)
> Rg=((1-2*(1-z)+pcopula(copclayton,cbind(1-z,1-z))))/(z)
> plot(c(1-z,z),c(Lg,Rg))
De plus, on a des fonctions similaires pour la dépendance au sens faible, avec le code suivant pour la version théorique,
> Lg=log(pcopula(cop,cbind(z,z)))/log(z)
> Rg=log((1-2*(1-z)+pcopula(cop,cbind(1-z,1-z))))/log(z)
> Lg=1/Lg*2-1
> Rg=1/Rg*2-1
et celui là pour la version empirique
> z=seq(0,.5,by=.001)
> v <- lossalae
> U=rank(v[,1])/(nrow(v)+1)
> V=rank(v[,2])/(nrow(v)+1)
> Lemp=rep(NA,length(z))
> Remp=rep(NA,length(z))
> for(i in 1:length(z)){
+ Lemp[i]=log(mean((U<=z[i])&(V<=z[i])))/log(mean(U<=z[i]))
+ Remp[i]=log(mean((U>=1-z[i])&(V>=1-z[i])))/log(mean(U<=z[i]))
+ }
> Lemp=1/Lemp*2-1
> Remp=1/Remp*2-1
Bref, on peut utiliser ces fonctions sur des vrais échantillons. Considérons l’exemple classique loss-alae (où l’on couple les frais dans des sinistres assurés, et les frais payés par l’assureur). On souhaite ajuster une copule, sans trop savoir laquelle. On peut commencer par étudier la dépendance forte, et comparer avec une copule Gaussienne. La copule Gaussienne de référence possède ici le même rho de Spearman que l’échantillon dont on dispose,
> library(evd); data(lossalae)
> cor(lossalae,method="spearman")
Loss ALAE
Loss 1.000000 0.451872
ALAE 0.451872 1.000000
> library(copula)
> paramgauss=.47
> paramclayton=.9
> paramgumbel=1.45
> copgauss=normalCopula(paramgauss)
> copclayton=claytonCopula(paramclayton, dim = 2)
> copgumbel=gumbelCopula(paramgumbel, dim = 2)
On obtient ici

La courbe verte est l’intervalle de confiance (ponctuel) à 95% pour une copule Gaussienne et un échantillon de même taille. On voit qu’on modélise mal la structure de dépendance. Avec une copule duale de Clayton, on obtient

et enfin pour une copule de Gumbel,

Bref, la copule de Gumbel semble réellement bien adaptée... Si on creuse en étudiant la dépendance au sens faible, on peut valider là aussi ce modèle. En effet, si la référence est la copule Gaussienne,

ou pour une copule de Clayton,

alors qu’une copule de Gumbel donnerait

Transfert du blog, suite

Arthur Charpentier — Thu, 01 Jul 2010 19:18:00 +0200

Autant l’admettre rapidement, transférer un blog n’est pas aussi simple que j’aurais pu le penser... Certes, l’essentiel des billets sont désormais lisibles sur freakonometrics.free.fr, mais les catégories n’ont pas été transférées, ainsi que les commentaires. Et comme certains commentaires qui m’ont été fait sont plus intéressants que les billets eux-mêmes, je trouve que ça manque. Mais surtout le référencement a du mal à suivre...

Bref, autant faire simple, et je vais continuer à nourrir les deux blogs en parallèles (quite probablement à ne pas mettre forcément tous les billets en ligne ici).

Econometrica et Ouest France: même combat ?

Arthur Charpentier — Fri, 18 Jun 2010 09:45:00 +0200

Un court billet aujourd’hui pour présenter deux ou trois petits graphiques (je n’ai malheureusement pas encore trouvé comment aller plus loin sur l’analogie...).

La recherche et l’impact factor des revues

Dans le système d’évaluation des chercheurs (je pourrais même dire des enseignants chercheurs) on est évalué tout simplement en fonction des revues dans lesquelles on arrive à placer nos papiers. Les revues sont en effet notées, classées, le classement le plus connu étant l’impact factor (j’en avais parlé l’an dernier lors du versement de primes à des chercheurs à Lyon, ici). Quelques revues ont un facteur d’impact très élevé, et ensuite se trouve des centaines voire des milliers de revues spécialisées. La distribution de ce facteur d’impact a l’allure suivante,

Comme l’ont noté plusieurs auteurs (ici ou là), la loi de dite de Zipf pourrait bien reproduire cette distribution (j’en avais parlé il y a plus d’un an dans un billet sur la loi de George Kingsley Zipf (ici), que j’avais fait un peu après avoir joué avec la loi de Benford (ici et là)).

La presse quotidienne, hebdomadaire ou mensuelle, et les nombres de tirages

Le monde de la recherche essaye régulièrement de donner l’illusion d’une brillante autorégulation, avec une reconnaissance par ses pairs, des relecteurs anonymes des articles soumis dans les revues, etc. Mais le monde de la recherche n’est pas si éloigné que ça des querelles d’épiciers sur le tirages de la presse.... En fait, si on regarde ce qu’on appelle la diffusion de la presse payante technique et professionnelle (ici) et la presse payante grand public (là), on retrouve à très peu de choses prêt la même loi...

Bref, plus publier dans une revue à gros tirage veut dire qu’on est un meilleur chercheur... ou quelque chose du genre. Le code est ici pour ceux qui doutent,
> presse=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/presses.csv",sep=";",header=TRUE)
# presse=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/Science_Journal_Ranking_Version_2003.csv",sep=";",header=TRUE)
> X=as.numeric(as.character(presse$Diffusion))
> X=X[is.na(X)==FALSE]
> Xs=sort(X)
> plot(length(X):1,Xs,col="blue",xlab="Ranking",ylab="Tirage (diffusion France payéé)")
> plot(log(length(X):1),Xs,col="blue",xlab="Ranking (log)",ylab="Tirage (diffusion France payéé)")

Liens avec la théorie des valeurs extrêmes

En fait, la loi de Zipf est tout simplement une loi de Pareto dans le cas discret (oui, les chercheurs ont un don pour donner des noms différents à des objets identiques, ça permet de faire plusieurs papiers sans que personne ne se doute de rien... pour ceux qui en doute, il suffit de croiser sous Google finance avec Zipf, Pareto, Power law, on va retrouver plein de monde qui raconte exactement la même chose). Bref, la loi de Pareto donne des choses très proches,

y compris la loi exponentielle qui est un cas limite de la loi de Pareto

On peut d’ailleurs aller un peu plus loin, en estimant le paramètre de la puissance, avec l’estimateur de Hill (discuté là),

ce qui confirme que les lois sont proches (des lois puissance), même si dans le cas de la presse française, on est sur un indice proche de 1.5 (à droite), alors que la presse académique est de l’ordre de 2 (à gauche), avec toutefois une explosion du côté des extrêmes (les 30 meilleures revues ayant visiblement un impact beaucoup trop important).

Le paradoxe de Saint Pétersbourg, partie -2

Arthur Charpentier — Thu, 17 Jun 2010 22:17:00 +0200

(oui, je commence à numéroter comme les Donjons, ou comme lorsque je commence des sections ou des sous-sections quand je fais cours au tableau). J’avais fait l’autre jour un billet sur le paradoxe de Saint-Pétersbourg, expliquant que, bien que l’espérance de ce jeu soit infinie, ça reste un jeu qui s’achève très rapidement, car en moyenne, on n’a le droit de ne faire que deux lancers...
Mais avant d’aller plus loin, un petit retour en arrière s’impose. En fait, j’ai toujours présenté en cours la théorie de l’espérance de l’utilité comme une réponse à ce paradoxe. Je me rend compte, après avoir cherché un peu, que c’est un peu plus compliqué ça. Installons nous confortablement, je vais raconter la petite histoire....

En 1738, Daniel Bernoulli a publié un ouvrage intitulé Specimen theoriae novae de mensura sortis (que nous pourrons traduire théorie sur la mesure du risque, dont la traduction est longuement commentée dans un document de Pierre Charles Pradier, ici, qui note que sortis est lié au sort, et donc au hasard, au risque). C’est dans cet ouvrage qu’il introduit le concept d’espérance morale.
Toute cette théorie a été reprise dans la théorie analytique des probabilités publiée par Laplace en 1812,

La construction de la fortune morale

Pour comprendre le cheminement de Daniel Bernoulli, supposons qu’une personne passe d’une somme de monnaie à une somme . Il pense que la valeur relative de cet incrément doit être proportionnel à l’incrément , mais inversement proportionnel à la richesse . Autrement dit, la variation de la valeur vérifie

ou encore

Pour reprendre la terminologie de Laplace, est la fortune physique et y la fortune morale. Dans la langue de Daniel Bernoulli, est appelé emolunmentum, et summa bonorum. Introduisons alors un peu d’aléa. La personne commence avec une dotation physique . Il peut alors gagner avec une probabilité . Daniel Bernoulli propose alors de poser

Daniel appelle emolunmentum medium. En conservant une expression de la forme

il convient de définir la fortune physique comme

Pour reprendre la terminologie de Laplace, est "l’accroissement de la fortune physique qui procurerait à l’individu le même avantage moral qui résulte pour lui, de son expectative", ou encore, selon Daniel Bernoulli est le "lucurm legitime expectandum seu sors quoesita".
Il considère alors deux états possibles, de probabilités respectives . La fortune physique est alors

On supposera que le jeu est juste, autrement dit , ou avec les notations de Laplace,

}

Daniel Bernoulli (puis plus tard Laplace) montrèrent que même si le jeu était juste, la fortune physique était toujours plus petite que . En fait, la démonstration est simple, car comme on a une mesure de probabilité, , et donc

Autrement dit, on doit simplement montrer que

L’astuce pour montrer ce dernier résultat consiste à noter que le terme de droite peut s’écrire

ce que revient à comparer une moyenne géométrique et une moyenne arithmétique (ce que j’avais évoqué ici avec des dessins).
A partir de cette théorie, Daniel Bernoulli propose ensuite une application en assurance. Un armateur possède une fortune physique , et il espère toucher x si le bateau arrive à bon port (ce qui surviendra avec probabilité ). S’il ne s’assure pas, sa fortune physique sera

En invoquant encore une fois cette comparaison entre moyenne géométrique et moyenne arithmétique, on en déduit que

Daniel Bernoulli en déduit que l’armateur a intérêt à s’assurer dès lors que la prime d’assurance ne dépasse pas (qui est la prime pure en terminologie actuarielle).
En notant e la prime d’assurance, Bernoulli montre qu’il faut chercher ce que nous appellerions un équivalent certain, solution de

Il essaye alors de résoudre ce problème comme une équation en (la prime d’asssurance e étant supposée donnée). Mais je m’égare un peu....

Le paradoxe de Saint Petersbourg

Une fois présenté son problème d’assurance, Daniel Bernoulli reprend un problème qu’il avait soumis à Pierre Rémond de Montmort en septembre 1713 (ici).

Dans ce jeu, un joueur lance un pièce en l’air. Si "face" apparaît au premier lancer, il touche 1 shilling. S’il tombe sur "pile", il peut relancer la pièce. Si "face" apparaît au second lancer, il touche 2 shilling. S’il tombe sur "pile", il peut relancer la pièce. Si "face" apparaît au troisième lancer, il touche 4 shilling. S’il tombe sur "pile", il peut relancer la pièce, etc. A chaque fois qu’il fait "face", il double ses gains.
L’espérance de gain est alors

qui est une somme infinie.
Pour résoudre le problème, Daniel Bernoulli calcule la fortune physique i.e.

qui est finie dès lors que est finie. Par exemple avec , on obtient 2. Pour ceux qui veulent s’en convaincre (numériquement)
> x0=0
> prod((x0+2^(0:100))^(1/2^(1:101)))-x0
[1] 2
Pour une richesse initiale de 1000, Bernoulli montre que la fortune physique est alors de 6,
> x0=1000
> prod((x0+2^(0:100))^(1/2^(1:101)))-x0
[1] 5.972253
Bref, Daniel Bernoulli avait ainsi montré que ses calculs de fortune physique permettait d’expliquer pourquoi les personnes jouant à son jeu étaient près à jouer une somme assez faible (en tous les cas comparé à la valeur espérée de gain).

Montmort, Bernoulli(s) et Cramer

En fait, si Montmort avait évoqué ce paradoxe dès 1713, on peut retrouver également ce problème soulevé en 1728 par Gabriel Cramer, toujours dans des lettres échangées avec Bernoulli, mais cette fois avec Nicolas, et pas Daniel, ici,

A la lecture de tout cela, on se rend compte que le paradoxe de Saint-Pétersbourg a été proposé à Bernoulli par Montmort et Cramer, que Cramer introduit également la notion d’espérance morale, bien avant Bernoulli.... L’histoire est souvent bien injuste dès qu’on se penche sur les héritages...

Transfert du blog....

Arthur Charpentier — Fri, 11 Jun 2010 16:14:00 +0200

Le blog devrait progressivement migrer depuis

http://blogperso.univ-rennes1.fr/arthur.charpentier/index.php/

vers cette nouvelle adresse,

http://freakonometrics.blog.free.fr/index.php/

Encore quelques réglage à faire, et comprendre comment faire proprement du dotclear avec free... En attendant, je vais probablement mettre les deux blogs à jour de manière sumultanée. Encore quelques semaines avant l’inauguration.... et un basculement définitif probablement pour la rentrée...

Paradoxe de Simpson expliqué par un dessin

Arthur Charpentier — Thu, 10 Jun 2010 22:36:00 +0200

Il y a plusieurs mois, j’avais fait un billet sur le paradoxe de Simpson (ici). Pour reprendre un exemple simple, supposons que l’on a deux hôpitaux, affichant les statistiques suivantes, pour les personnes "saines",

hôpital	total	survivants	décès	taux de survie
hopital A	600	590	10	98%	x
hopital B	900	870	30	97%

et pour les personnes "malades" (on suppose vraiment que le critère est identique pour les deux hôpitaux)

hôpital	total	survivants	décès	taux de survie
hopital A	400	210	190	53%	x
hopital B	100	30	70	30%

Bref, peu importe son état de santé, on a toujours intérêt à choisir l’hôpital A. D’où le paradoxe. En revanche, si on combine, on obtient "globalement",

hôpital	total	survivants	décès	taux de survie
hopital A	1000	800	200	80%
hopital B	1000	900	100	90%	x

on note que l’on a intérêt à choisir l’hôpital B. Étonnant non ? (l’exemple vient de Marco Scarsini, et je n’ai jamais trouvé plus élégant).
J’évoque souvent ce paradoxe en économétrie, et c’est ce que j’avais fait dans mon précédant billet (ici). En fait, Roger Nelsen propose un petit exemple simple permettant de comprendre comment on peut obtenir un tel paradoxe. L’idée est de noter que

et pourtant

Cela n’a a priori rien d’exceptionnel... On peut parfaitement avoir

tout en ayant,

C’est ce que l’on peut visualiser ci-dessous,

Merci Roger....

Dessins et preuves ?

Arthur Charpentier — Thu, 10 Jun 2010 09:10:00 +0200

Plusieurs personnes m’ont laissé des messages et des commentaires suites aux derniers dessins, disant en substance "c’est bien joli, mais est-ce que c’est vraiment une preuve ? ". Etant incompétent en épistémologie, je n’aurais pas la prétention d’apporter des lumières... En revanche, Olivier m’a forwardé (ici) le document en ligne là (on l’on retrouve page 4 un dessin plus beau que celui que j’avais pu faire pour montrer que la loi géométrique est une vraie loi de proba), et surtout, il m’a fait découvrir ici la notion de "preuve sans mots", où l’on retrouve une très jolie preuve animée (je me rends compte que l’animation aide vraiment à la compréhension, c’est ça que j’aurais du faire depuis le début) du théorème de Pythagore (que je continue de mentionner en cours, comme là, lorsque je dois parler de formule de décomposition de la variance). Bref, sur la page wiki est mentionné des ouvrages de Roger Nelsen, que je dévore depuis ce matin.

(ainsi que quelques livres avec Claudi Alsina). Pour la petite histoire, c’est justement Roger qui m’avait expliqué comme obtenir ces propriétés de la loi géométrique, à l’époque où je faisais un short course à Samos, en , sur advanced statistical methods in insurance. J’avais évoqué l’importance de la loi géométrique dans le concept de période de retour, et Roger m’avait expliquer comment retrouver l’histoire de l’espérance de la loi géométrique sans équations....

J’ai retrouvé plusieurs dessins que j’avais fait dans mes cahiers de prépas (que je commence à ressortir ces derniers temps, depuis que je dois revenir sur quelques concepts fondamentaux de convexité et de topologie... et qui m’avait servi pour les dessins que j’avais déjà pu faire, ici en particulier), et surtout j’en ai découvert un paquet d’autres, que je ne pourrais pas m’empêcher d’évoquer dans de futurs billets... D’ailleurs, je me suis rendu compte qu’on pouvait retrouver un bon nombre de ces exemples directement sur jstor, là. Par exemple, il propose une interprétation de la comparaison de plusieurs moyennes (arithmétique, géométrique, etc), plus complète que ce que j’avais pu faire,

Sur la preuve par le dessin, Jean-Paul Delahaye apporte les éclaircissements suivant, ici, en citant David Hilbert qui disait : « les règles pour écrire une démonstration doivent être si claires que si quelqu’un vous propose une démonstration, alors une procédure mécanique doit pouvoir assurer que la démonstration est correcte, c’est-à-dire qu’elle obéit aux règles. ». Il détaille par ailleurs cet argument, "ce que les Grecs – par exemple Euclide – considéraient comme une démonstration est aujourd’hui encore, moyennant quelques adaptations dans le style et les notations, toujours considéré comme une démonstration. En revanche, nous avons beaucoup avancé sur les notations, la structure logique des arguments, et la compréhension de ce qu’il est possible d’attendre de la méthode axiomatique. Nos notations et la logique mathématique ont tellement progressé que, contrairement aux Grecs, nous savons formaliser les démonstrations que nous écrivons."
Bref, pour ma part, je trouve cela tellement élégant (et surtout je sais que Roger est beaucoup plus compétent en maths que moi) que je ne m’empêcherais pas de mettre ici le compte rendu de mes lectures ! Merci Roger, et Olivier....

Discussions géométriques du soir (partie 2)

Arthur Charpentier — Mon, 07 Jun 2010 23:26:00 +0200

J’avais raconté l’autre jour l’aventure de mon fils qui souhaitait faire des triangles avec des bouts de bois (ici). Il a récidivé ce soir. Sauf qu’il n’avait qu’un grand baton, et que j’avais fermé le portillon qui l’empêchait d’aller dans la forêt en chercher (la nuance entre romancer et devenir mythomane doit être mince). Je lui ai suggéré de le casser en trois portions, puis de constituer un triangle avec les trois morceaux restant. Manque de bol, l’expérience rate (comme l’autre fois). Et me souvenant de mes pénibles calculs de la dernière fois, je me suis rappelé qu’il avait une chance sur deux de pouvoir effectivement construire un triangle... sauf que la manière d’obtenir les trois bouts de bois me semble différente... donc peut être que ce 1/2 n’est plus valide ?
En fait, je crois qu’il faut définir clairement comment on va couper les bouts de bois.

méthode 1: le coupage simultané

Je suppose que je peux couper d’un coup le bout de bois en deux. Formellement, je tire deux points sur l’intevalle correspondant à mon bout de bois, les tirages étant uniformes, et indépendants. On arrive alors à trois morceaux, , et , en sachant que

la longueur du bout de bois dont je diposais initialement. On supposera que =1 pour simplifier (ce qui ne changera rien au résultat final, on s’en doute). L’ensemble des points que je peux obtenir peuvent alors se visualiser comme des points du simplexe,

autrement dit, on peut les représenter dans le triangle équilatéral ci-dessous.
Or on ne peut constituer un triangle que si la plus grande des longeurs n’excède pas la somme des deux autres. Supposons que soit la plus grande longueur (et qu’elle soit visualisée graphiquement par rapport à la distance au côté en bas), le point rouge ci-dessous (correpondant à un triplet de longueurs possibles) ne convient pas. En fait, dès lors que l’on entre dans l’aire rosée, on ne peut plus constituer un triangle...
Moralité, à moins de tomber dans le triangle jaune, on ne peut pas constituer de triangle.... Ce qui veut dire qu’on n’a plus qu’une chance sur quatre !

méthode 2: le coupage séquentiel

Une autre méthodde consister à couper d’abord en deux (en tirant un point au hasard sur mon intervalle), puis de recouper en deux le plus grand morceau... Damned, mon fils a de bonne idée, mais ça sent les calculs à n’en plus finir ça...Ou alors va falloir continuer à essayer de faire des dessins, car mine de rien, c’est plus simple que les calculs de la dernière fois. Bon, la longueur de mon bout de bois est la hauteur de mon triangle. Je commence par le couper en deux, la distance au côté en bas devant être la plus petite des deux.

Une fois fixé la taille de ce petit bout, il me reste à couper le plus grand morceaux en deux. Pour cela, je vais me déplacer le long du segment rouge ci-dessous,

La position sur ce petit segment me donnant ensuite les deux autres longueurs des morceaux. Or rappelons que pour que nos bouts de bois forments un triangle, il faut forcément appartenir au triangle jaune,
Pour résumé, nous avions exclus la partie supérieure du triangle, et conditionnellement au fait que nous sommes dans le lozange inférieur, nous souhaitons être dans le triangle jaune. Arg, cette fois-ci la probabilité passe à une chance sur trois !

couper d’un coup comme découpage séquentiel aléatoire ?

Mais attends deux minutes mon bonhomme (oui, je parle à mon fils). Si j’avais choisi de découper le plus petit morceau, on sait qu’il est impossible de constituer un triangle car par définition Supposons que l’on casse le baton en deux, et que l’on tire ensuite au hasard (avec probabilité 1/2 et 1/2) quelle partie sera alors coupée en deux, on doit se ramener au cas précédant non ? La probabilité de constituer un triangle doit alors être

Damned, mais ça ne fait pas un quart ça ! Bon, je crois que je vais prendre une aspirine et aller demander à mon fils d’aller ranger ses bâtons !

Comprendre les tableaux croisés, 1

Arthur Charpentier — Mon, 07 Jun 2010 11:12:00 +0200

Je voulais faire quelques billets rapides sur les tableaux croisés, suite à des demandes consécutives à mon billet sur le sexe des bébés (ici). Pour commencer, je vais m’intéresser au cas où les variables supposées indépendantes (d’autres billets seront l’occasion de tester d’autres types d’hypothèses). Histoire de trouver des données un peu originales, on peut traîner sur internet. Etant tombé (presque) par hasard sur le fichier pdf ici, j’ai voulu m’intéresser à l’origine des enseignants-chercheurs recrutés lors de la campagne 2008 (comme le dit le titre du document en ligne sur le site du ministère). Prenons deux variables catégorielles géographiques, par exemple

, ville (au sens large) où l’on a soutenu sa thèse
, ville (au sens large) où l’on a obtenu un poste de maître de conférence

Compte tenu des règles de non-recrutement local, on devrait pouvoir supposer que ces variables sont indépendantes. En fait, je ne raisonne pas ici par université, mais par ville, ce qui correspond à un recrutement local au sens large (passer de Paris 6 à Paris 7 sera ici considéré comme "local").
L’étude a été faite à partir des données disponibles ici et là. Petite précaution d’usage: ce sont les données de la campagne 2008, avant que les comités de sélection ne soient mis en place. Il conviendrait de refaire l’étude une fois les nouvelles données publiques. On a deux bases, une en droit-économie-gestion, et une en sciences. J’insiste sur le fait que l’on a des données sur une année seulement, donc avec peu d’observations... Mais comme l’indique le titre, je vais parler des outils visualisation des tableaux croisés, pas faire un billet sur le recrutement local. Voilà pour les précautions d’usage...
Pour visualiser ici, j’ai fait des regroupements, car avec 25 villes, ça ne tient pas dans mon billet (mais par la suite, l’analyse sera plus fine). Les données de base sont les suivantes, les lignes étant l’endroit où un doctorant a soutenu sa thèse (en économie ou droit), et en colonnes, l’endroit où a été recruté un maître de conférence,

	Rennes	Province (autre)	Région parisienne	total
Rennes	5	8	0	13
Province (autre)	3	124	23	150
Région parisienne	1	13	49	63
total	9	145	72	226

Classiquement, on peut retraduire ces chiffres (qui sont des comptages) sous la forme de probabilités,

	Rennes	Province (autre)	Région parisienne	total
Rennes	2,21%	3,54%	0,00%	5,75%
Province (autre)	1,32%	54,86%	10,18%	66,37%
Région parisienne	0.44%	5,75%	21,68%	27,88%
total	3,98%	64,15%	31,86%	100,00%

Les probabilités qui apparaissent sur les totaux sont les probabilités marginales. Aussi, au bout de la première ligne on trouve que 5,75% des doctorants qui ont eu un poste en 2008 avaient soutenu leur thèse à Rennes (ou dans la région, i.e. les universités de Rennes 1, de Rennes 2 et de Brest). De manière similaire, en bas de la première colonne, 3,98% des maîtres de conférences recrutés en 2008 ont été recrutés à Rennes.
Deux autres tableaux classiques sont alors les profils ligne et profils colonnes, c’est à dire les probabilités conditionnelles,

	Rennes	Province (autre)	Région parisienne	total
Rennes	38,46%	61,54%	0,00%	100,00%
Province (autre)	2,00%	82,67%	15,33%	100,00%
Région parisienne	1,58%	20,64%	77,78%	100,00%
total	3,98%	64,15%	31,86%	100,00%

par ligne, autrement dit, pour un étudiant qui a soutenu sa thèse à Rennes (et qui a eu un poste de Maître de Conférences), il avait 38,46% de chance de rester à Rennes (ou plutôt en Bretagne), et 61,54% de chance d’aller dans une autre ville de province. Côté colonnes,

	Rennes	Province (autre)	Région parisienne	total
Rennes	55,56%	5,52%	0,00%	5,75%
Province (autre)	33,33%	85,52%	31,94%	66,37%
Région parisienne	11,11%	8,96%	68,06%	27,88%
total	100,00%	100,00%	100,00%	100,00%

Ici, on lit que pour un Maître de Conférences recruté sur la région parisienne, il y avait 68,06% de chances pour qu’il ait fait en région parisienne.
Bref, à partir de ces tableaux, on peut essayer de tester l’hypothèse d’indépendance, et si cette dernière n’est pas vérifiée, on va essayer de comprendre pourquoi.
Mais avant, rappelons que l’hypothèse d’indépendance se traduit par le fait que

La valeur de gauche est la probabilité jointe (que l’on retrouve au coeur du second tableau), alors que les deux autres sont les probabilités marginales, qui apparaissent en bas, et à droite (dans les totaux). Bref, on pourrait construire le tableau des probabilités que l’on devrait avoir sous hypothèse d’indépendance,

	Rennes	Province (autre)	Région parisienne	total
Rennes	0,22%	3,69%	1,83%	5,75%
Province (autre)	2,64%	42,58%	21,14%	66,37%
Région parisienne	1,11%	17,88%	8,88%	27,88%
total	3,98%	64,15%	31,86%	100,00%

que l’on peut retraduire en terme d’effectif espéré, toujours sous l’hypothèse d’indépendance. Pour ceux qui connaissent les copules, on change juste la structure de dépendance ici, pas les lois marginales (ou les effectifs marginaux) qui restent inchangés.

	Rennes	Province (autre)	Région parisienne	total
Rennes	0,5	8,3	4,2	13
Province (autre)	6,0	96,2	47,8	150
Région parisienne	2,5	40,5	20	63
total	9	145	72	226

L’outils de base pour comparer ces données est la distance du chi-deux.
La statistique du chi-deux est

ou encore

Karl Pearson avait été un des premiers à étudier ces grandeurs, sous des hypothèses de normalité (des sommes de lois normales indépendantes au carré se ramenant à des lois du chi-deux).
On notera qu’une approche alternative est de supposer que les nombres (ou les effectifs) suivent une loi de Poisson. On peut alors se demander si

Si on suppose avoir une loi de Poisson, on peut s’intéresser aux résidus de Pearson,

autrement dit, on retrouve des grandeurs très proches de la contribution au chi-deux (en prenant le carré). On aussi le signe qui nous indique s’il y a trop peu de monde, ou au contraire par assez....
Si je travaille sur la base totale, soit environ 25 "villes", on obtient
> cont1=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/contingence-universite_1.csv",
+ header=TRUE,sep=";")
> cont2=cont1
> for(i in 1:ncol(cont2)){
+ for(j in 1:nrow(cont2)){
+ if(is.na(cont2[j,i]==TRUE)){cont2[j,i]=0} } }
> X=as.matrix(cont2[1:nrow(cont2),2:ncol(cont2)])
> rownames(X)=as.character(cont2[,1])
> colnames(X)=names(cont2[2:ncol(cont2)])
> X=X[,-2]
> X=X[-2,]
> X=X[,-9]
> X=X[-9,]
> Z=X[-nrow(X),]
> Z=Z[,-ncol(Z)]
> Z=Z[-nrow(Z),]
> chisq.test(Z)
        Pearson’s Chi-squared test
data: Z
X-squared = 1342.546, df = 576, p-value < 2.2e-16
La base est assez sale, j’ai fait un copier/coller rapide du fichier pdf (et je vire les Antilles et la Réunion car après je voudrais faire un peu de visualisation sur une carte). Autrement dit, l’hypothèse d’indépendance ne semble pas validée par le test du chi-deux. Essayons maintenant de comprendre un peu mieux les interactions dans notre carré...
Si on regarde les problèmes de recrutement local, il faut regarder ce qui se passe sur la diagonale de la matrice des résidus,
> CT=chisq.test(Z)
Message d’avis :
In chisq.test(Z) : l’approximation du Chi-2 est peut-être incorrecte
> Dia=diag(CT$residuals)
> names(Dia)=rownames(Z)
> Dia
AIX-MARSEILLE        BESANCON        BORDEAUX            CAEN        CLERMONT
     6.43901593      3.26116568      8.37261202     10.53607373      5.12692873
         CORSE           DIJON        GRENOBLE           LILLE         LIMOGES
     8.56426340      4.10931024      9.18869705      7.59456681      8.56426340
          LYON     MONTPELLIER           NANCY          NANTES            NICE
     3.16119625      4.18366071      8.84765951      2.67087742      7.86973591
       ORLEANS        POITIERS           REIMS          RENNES           ROUEN
     1.79480842      5.81144228     -0.06651901      6.22963060      6.57435350
    STRASBOURG        TOULOUSE           PARIS         CRETEIL      VERSAILLES
     5.64850464      4.12289099      5.46092001     -0.22061860      4.01351350
Visuellement, sur une carte de France, on obtient la distribution ci-dessous, pour les postes en économie et droit,

et ci-dessous pour les postes en sciences et pharmacie,

Si on regarde le tableau des résidus on obtient les valeurs suivantes sur la ligne correspondant à Rennes, i.e. les endroits où vont les étudiants ayant soutenu leur thèse à Rennes,

avec en bleu les valeurs positives, et en rouge les valeurs négatives. A droite sont représentées les lignes, i.e. les endroits d’où viennent les Maîtres de Conférences recrutés à Rennes.
Au lieu de regarder ville par ville, on peut aussi faire une petite animation, soit sur les endroits où sont partis les doctorants,

ou les endroits d’où viennent les Maîtres de Conférence recrutés,

On peut également s’il existe des interactions, ou des échanges d’une ville à l’autre. Dans ce cas, les résidus croisées et doivent être positifs.

Bref, au delà du "recrutement local" on voit également apparaître de fortes interactions régionales (Rennes-Nantes, Toulouse-Bordeaux, Orléans-Paris).

Avant de conclure, il existe une autre technique permettant de mieux visualiser les interactions, c’est l’analyse factorielle des correspondances (je renvois ici pour des slides). Si on regarde les 2 premiers axes principaux, en sciences (à gauche) et en droit et économie (à droite), on obtient les projections suivantes des modalités, avec en bleu les villes où les étudiants ont soutenu leur thèse, en rouge les villes où les postes de Maître de Conférences ont été obtenus,

et sur les les 2 axes suivants

Bref, comme souvent en statistique, afin de mieux comprendre un ajustement de modèle (ou comme ici lors la réalisation d’un test d’hypothèse) la visualisation des résidus est la technique de base (c’était l’idée d’ailleurs du commentaire que j’avais fait là).
Voilà pour la première partie, sur les modèles où l’on suppose l’indépendance. La prochaine fois, on verra d’autres formes d’hypothèses.

Le paradoxe de Saint Péterbourg, partie 1

Arthur Charpentier — Fri, 04 Jun 2010 21:50:00 +0200

Pour répondre à une question sur mon précédant billet (ici), je vais revenir sur un paradoxe assez classique, le paradoxe de Saint Petersbourg. Mais avant de parler du paradoxe, et de ses implications en théorie de la décision dans l’incertain, je voulais présenter le jeu, et en profiter pour jouer à faire des dessins puisque j’ai cru remarquer que j’avais des amateurs d’explications géométriques.
Le jeu est simple, c’est un jeu de pile ou face répété, le jeu s’arrêtant à la sortie du premier "face" (et on le verra par la suite, l’idée est de doubler ses gains chaque fois que "pile" sort).
Aussi, le temps d’arrêt du jeu est dont la loi est simplement

... etc. On retrouve ainsi une loi géométrique de paramètre 1/2. L’espérance du nombre de lancers que l’on ferra avant que le jeu ne s’arrête est fini (alors qu’on le verra le gain espéré est infini), et vaut l’inverse de cette probabilité, i.e. 2. Pour ceux qui ont oublié la formule (que l’on obtient de manière assez calculatoire, en faisant une dérivation dans une série), on peut refaire le calcul,

Personnellement, je ne sais pas calculer cette somme (infinie), mais comme dans mon précédant billet (ici), il existe une somme que l’on peut calculer simplement en faisant un petit dessin,

J’ai 1 carré (mauve) de surface 1, 2 (bleus ) de surface 1/2, puis 3 (verts) de surface 1/4, puis 4 (jaunes) de surface 1/8, puis 5 (rouges) de surface 1/16, 6 (oranges) de surface 1/32, etc. Autrement dit

Or l’espérance que l’on cherche à calculer, c’est cette somme, à un facteur 2 près.... Aussi, en moyenne, on peut espérer faire 2 lancers (la moitié de la surface) ! Ce qui correspond précisément à l’espérance d’une loi géométrique de paramètre 1/2. Décidément, on peut faire plein de choses avec des petits dessins..

Zénon, sommes infinies et loi géométrique

Arthur Charpentier — Fri, 04 Jun 2010 15:02:00 +0200

Ayant parlé plusieurs fois de Luc Chatel sur mon blog je n’ai pas pu m’enpêcher de mettre l’image qui est à la une du Canard de cette semaine.
Pour ceux qui se souviennent de leurs cours de philo au lycée, c’est le paradoxe de Zénon (ou presque). Pour cela il faut fermer les yeux et imaginer une flèche en vol (en fait non, je viens de faire l’expérience, si on ferme les yeux, on ne peut pas lire la suite, ça ne marche qu’en cours cette histoire). La flèche est lancée à 10 mètre de sa cible. Il faut un certain temps, non nul, à cette pierre pour parcourir la moitié de la distance qui la sépare de la cible. Ensuite, il lui reste encore 5 mètres à parcourir, dont elle accomplit d’abord la moitié, 2,5 mètres, ce qui lui prend un certain temps. Puis elle parcourt la moitié, puis la moitié, et cela - à chaque fois - avec un temps non nul. Zénon (Ζήνων, parfois appelé aussi Zénon d’Elée) en conclut que la pierre ne pourra frapper l’arbre qu’au bout d’un temps infini, c’est-à-dire jamais (ou pour la version illustrée, je peux renvoyer à un extrait de Kid Paddle, ici). En fait, Zénon introduit la notion d’infini divisibilité (on dira qu’un élément est indéfiniment divisible, c’est d’ailleurs surprenant ce lien en français entre l’infini et l’indéfini), qui philosophiquement, est un complexe qui a mis du temps à s’imposer: une somme infinie d’éléments non nuls peut être finie (je renvois ici par exemple).

Formellement, l’histoire de la flèche revient à se demander si on peut calculer quelquechose qui ressemble à

soit, de manière plus formelle,

on retrouve la somme d’une série géométrique (au premier terme près), dont on sait que

Un application intéressante est que la loi géométrique, vérifiant

est effectivement une loi de probabilité,

(les termes sont positifs, et la somme vaut 1). Il existe une démonstration élégante de ce dernière résultat, qui tient en un dessin (à condition d’accepter ce principe un peu fractal d’infini divisibilité):

Je reviendrais un jour sur les lois de probas infiniment divisibles, car elles sont très très intéressantes de par certaines propriétés....

Les abaques comme outils de pédagogie statistique

Arthur Charpentier — Fri, 04 Jun 2010 10:18:00 +0200

Comme je l’avais mentionné à propos d’un vieux papier paru dans Le Canard, c’est souvent compliqué de communiquer sur les retraites. C’est un sujet complexe et assez rapidement technique. Historiquement, les abaques sont les tables de calculs,

Mais il faut admettre que retrouve un peu n’importe quoi dans ce terme d’abaque.... Certaines sont évoquées ici ou là, voire là pour quelque chose de très joli. On notera qu’il existe une idée sous jacente commune, à savoir construire un outil graphique permettant de visualiser comment 3 paramètres peuvent être liés (générallement, ils sont liés par une fonction complexe).

Abaques et retraites

Parmi les outils que l’on retrouve dans les documents du COR, j’ai été surpris de voir ressortir l’utilisation des abaques (par exemple ici). On nous en fait même la promotion en vidéo,

Pour la petite histoire, j’avais évoqué son utilisation dans un graphique ici sur les problèmes d’échantillonnage. C’est effectivement un outil graphique intéressant, à condition de savoir le décoder...
L’idée est simple, à savoir lier les trois paramètres essentiels (ou supposés comme tels) du bilan d’un système de retraite, à savoir

le niveau des prélèvements,
le niveau de l’âge effectif de départ en retraite
le niveau des pensions.

La difficulté est de trouver une forme graphique permettant de relier ces trois niveaux (la dimension 3 étant souvent trop abstraite). Par exemple dans les sondages, pour construire un intervalle de confiance, les abaques sont utilisées (comme ici) pour relier

le niveau des intervalles de confiance (5%-95% ou 10%-90%)
le niveau de probabilité attendu dans le sondage
la taille de l’échantillon interrogé

On peut d’ailleurs trouver d’autres exemple en statistiques ici ou là. Mais pour revenir au problème des retraite, sur le graphique ci-dessous (appendix 1 ici), en abscisse on met le taux de remplacement (c’est à dire le ratio de la pension moyenne par rapport au revenu moyen) et en ordonnée on met le niveau d’un hausse du taux de cotisation (en points de cotisation). chaque ligne représente un âge moyen de départ en retraite effectif.

On se fixe alors des scénarios d’équilibre et on regarde la manière dont se déplacent les équilibres. J’avais prévu toute une analyse technique, mais je viens de me rendre compte qu’Antoine a mis en ligne un billet remarquable sur le sujet sur http://www.ecopublix.eu/ (ici), que je ne pourrais jamais égaler en clareté. Donc je vais proposer une autre construction d’abaque afin de montrer son intérêt pédagogique.

Comprendre l’actualisation

Sur le site http://images.math.cnrs.fr/ (ici), Xavier Caruso explique son émerveillement devant la difficulté de faire des calculs d’actualisation. En un sens, ça me rassure de voir que les calculs actuariels ne paraissent pas triviaux à tout le monde (et surtout plein de contre-intuitions). Xavier essaye d’expliquer la différence entre

ouvrir un compte rémunéré à 4,5% sur lequel on dépose 500 euros chaque mois pendant 10 ans
ouvrir un compte rémunéré à 4,5% sur lequel on dépose 60000 euros initialement sans plus jamais l’alimenter ensuite.

Si effectivement la somme des 120 versements de 500 euros correspond aux 60000 euros, la valeur intégrant les intérêts n’est pas du tout la même au bout de 10 ans, car dans le premier car, comme le disent les banquiers, "les intérêts se capitalisent". Bon, si Xavier fait des maths (et donc des développements limités par avoir des formules fermés), on va plutôt faire ici un peu de calculs numériques. Tout d’abord, pour visualiser l’évolution des montants disponibles sur les comptes, rien de plus simple,
> tps=1/12*(1:120)
> plot(tps,60000*(1+0.045)^tps,ylim=c(0,100000))
> lines(tps,cumsum(500*(1+0.045)^(tps)))
les montants au bout de 10 ans étant
> sum(500*(1+0.08525)^(1/12*(1:120)))
[1] 93180.33
> sum(500*(1+0.045)^(1/12*(1:120)))
[1] 75514.32
Ce qui confirme l’intuition que nous avions. Comme le note Xavier, il faudrait un taux d’intérêt deux fois plus grand pour le placement mensuel, la preuve
> sum(500*(1+0.09)^(1/12*(1:120)))
[1] 95543
> sum(500*(1+0.08525)^(1/12*(1:120)))
[1] 93180.33
Xavier obtient cette valeur en effectuant des développements limités. Il résume ça à l’aide du dessin

On peut aller plus loin, en notant qu’à taux sur le placement mensuel donné, ainsi que la maturité, le taux sur le placement bloqué qui rapporte la même valeur à échéance est donné par

> ((sum(500*(1+0.09)^(1/12*(1:120))))/(500*12*10))^(1/10)-1
[1] 0.04762238
(dans le cas d’un placement rapportant 9%). En fait, si on regarde le montant dont on dispose sur les comptes à chaque date, on obtient précisément les valeurs suivantes, en fonction du temps,

on retrouve effectivement qu’à 4,5%, les deux placements ne donnent pas le même montant à échéance (ici 10 ans). Alors que si le placement alimenté tous les mois rapportait 8,525%, on aurait exactement la même somme au bout de 10 ans

ou si on compare avec un placement à 9% (le double, ce qui est obtenu par développement limité), on obtient

Pour trouver la valeur exacte du multiplicateur en fonction du taux offert sur le compte bloqué et la maturité, la fonction est tout simplement
> taux=4.5/100
> T=10
> f=function(k){(1*12*T)*(1+taux)^T-sum(1*(1+k*taux)^(1/12*(1:(12*T)))) }
> uniroot(f,interval=c(0,20))$root
[1] 1.894332
L’idée d’utiliser une abaque peut s’avérer intéressante, car on essaye ici de relier

le taux de rendement de l’argent bloqué
le taux proposé pour le placement alimenté tous les mois
la durée envisagée du placement

Je peux mettre en abscisse le taux du placement bloqué. En ordonné, je met le multiplicateur utilisé pour le placement alimenté tous les mois. A maturité donné (par exemple 10 ans en rouge), on peut regarder la valeur du multiplicateur pour que les valeurs à échéance coïncident,

On retrouve directement que le taux équivalent pour un placement alimenté tous les mois doit être 1,89 fois le taux offert sur le placement bloqué si on se fixe une échéance de 10 ans pour égaliser les placements. Le code est simplement le suivant
> tx=seq(0.005,.05,by=.001)
> tps=seq(5,25,by=.25)
> M=matrix(NA,length(tps),length(tx))
> for(i in 1:length(tps)){
+ for(j in 1:length(tx)){
+ T=tps[i]
+ taux=tx[j]
+ f=function(k){
+ (1*12*T)*(1+taux)^T-sum(1*(1+k*taux)^(1/12*(1:(12*T))))
+ }
+ M[i,j]=uniroot(f,interval=c(0,2))$root
+ }}
> contour(tps,tx,M,lwd=2,
+ xlab="Durée du placement (en années)",ylab="Taux (du placement bloqué)")
> abline(v=seq(5,25,by=1),lty=2,col="grey")
> abline(h=seq(0,.05,by=.005),lty=2,col="grey")
Bref, l’actualisation est un sujet délicat (même ce n’est que des calculs de sommes et de puissances). Et faire des petits dessins permet souvent de mieux comprendre...

Qu'est ce que la modélisation économétrique (3) ?

Arthur Charpentier — Thu, 03 Jun 2010 15:00:00 +0200

Compte tenu du nombre de billets qui agitent la blogosphère en ce moment, je me suis senti obligé de reprendre la plume pour poursuivre le billet précédant (ici). Initialement, le but était de montrer aux élèves comment mener une étude économétrique et de répondre à une question simple (en l’occurence "la taille d’un classe a-t-elle un impact positif ou négatif sur les résultats scolaires ? "). Malheureusement, j’ai été rattrapé par l’actualité. Le premier billet montrait qu’en faisant un modèle de régression simple, l’effet semblait positif: plus la taille est grande, meilleure est la moyenne. Dans le second billet, nous avions noté que cet effet pouvait cacher quelquechose, comme des conditions socioéconomiques. Essayons donc d’aller un peu plus loin....

Régression multiple

En faisant une estimation par moindres carrés sur les deux variables, on obtient
> summary(lm(avgverb~tipuach+classize,data=base0))
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 80.26363    0.74472 107.777   <2e-16 ***
tipuach     -0.34994    0.01077 -32.486   <2e-16 ***
classize    -0.03146    0.02221 -1.416    0.157
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
L’impact de l’indice socio-économique est clairement négatif, et si l’effet taille a un impact négatif, ce dernier n’est pas significatif. Il va falloir essayer d’aller un peu plus loin....

Approche nonlinéaire

On peut ragarder rapidement ce qu’aurait donné une régression nonlinéaire
> library(mgcv)
> reg=gam(avgverb~s(tipuach,classize),data=base0)
> reg=gam(avgmath~s(tipuach,classize),data=base0)
soit visuellement, les prédictions suivantes pour l’épreuve de lecture

et la surface ci-dessous pour l’épreuve de maths,

Bref, l’effet ne parait pas forcément simple, même si l’on intuite une évolution croissante. Mais depuis le début, nous omettons un point important sur la base de Joshua Angrist et Victor Lavy: dans ces écoles israëliennes, la règle de Maimonides est supposer s’appliquer... et jusqu’à présent je n’en ai jamais parlé !

Prise en compte de la règle de Maimonides

Compte tenu de la règle de Maimonides, la taille des classes en fonction de la taille de l’école n’est pas du tout linéaire, comme le montre la figure ci-dessous,

En fait, si on défini le nombre théorique d’élèves dans les classes, avec cette méthode (ce que j’avais détaillé dans la section 3 ici), on obtient des résultats qui pourrait ressemblerà ce que nous avions observé sur le nombre réel d’élèves dans les classes,
> reg=lm(avgverb~func1,data=base0)
> summary(reg)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 70.62076    0.88004   80.25 < 2e-16 ***
func1        0.12159    0.02789    4.36 1.37e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
On retrouve un effet positif (significatif) qui se visualise sur le graphique ci-dessous,

et surtout, le modèle devient particulièrement intéressant si on régresse sur cette effectif théorique et le "percent disadantaged" (PD, décrit dans le précédant billet, là)
> reg=lm(avgverb~func1+tipuach,data=base0)
> summary(reg)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 82.91057    0.78669 105.392 < 2e-16 ***
func1       -0.11160    0.02321 -4.809 1.63e-06 ***
tipuach     -0.35943    0.01050 -34.229 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
où cette fois l’impact de la taille de la classe est clairement significatif, et négatif: à charactéristiques socio-économique identiques, si la taille des classes suivait la règle de Maimonides, la taille de la classe aurait un impact négatif sur les résultats scolaires....
Une solution est d’utiliser une régression par variables instrumentales où l’instrument permet précisément de prendre en compte cette discontinuité. L’idée est ici qu’il existe une variable telle que

On voit alors apparaître naturellement l’idée d’utiliser les variables instrumentales pour corriger d’un éventuel biais lors de l’estimation des coefficients lors de la régression (comme cela est développé dans le chapitre 25.7 du microeconometrics de Colin Cameron et Pravin Trivedi). Et c’est précisément ce qui survient dans nos écoles compte tenu de la règle de Maimonides, On utilise la variable que nous avions définie comme le nombre théorique d’élèves dans la classe !
Damned, l’étude n’est donc toujours pas finie, il va falloir creuser encore davantage (surtout qu’il reste un paquet de variables dans cette base)... mais pour les impatients, la conclusion à laquelle arrivent Joshua Angrist et Victor Lavy est résumé dans le tableau ci-dessous,

Espérance de vie et retraite

Arthur Charpentier — Thu, 03 Jun 2010 11:49:00 +0200

Dans un billet qui a presque un an (ici), j’étais revenu sur l’idée reçue comme quoi "tous les ans, on gagne un trimestre d’espérance de vie" . J’avais du alors expliquer que c’était effectivement le cas, en tous les cas pour l’espérance de vie à la naissance. C’est ce que raconte le graphique ci-dessous, avec une projection par la méthode de Lee & Carter, pour les femmes uniquement

> library(demography)
> france.LC1 <- lca(fr.mort,adjust="e0",series="female",years=c(1900,2040))
> france.fcast <- forecast(france.LC1)
> L2 <- lifetable(france.fcast)
> ex2=L2$ex
> L1=lifetable(fr.mort,series="female")
> ex1=L1$ex
> age=0
> exF=c(ex1[age+1,],ex2[age+1,])
> plot(1816:2056,ex,col="blue")
On peut alors faire une régression pour quantifier un peu mieux ce qui s’est passé au cours des 50 dernières années.
> I=(1950:2000)-1815
> y=exF[I]
> x=1950:2000
> lm(y~x)
Coefficients:
(Intercept)            x
-438.7249       0.2611
Avec un pente de 0.26, effectivement, pour une année de plus (entre 1980 et 1981 par exemple), l’espérance de vie à la naissance gagne 0,26 année, soit un trimestre. Pour les hommes, on a
> I=(1950:2000)-1815
> y=exH[I]
> x=1950:2000
> abline(lm(y~x),col="red")
> points(x,y,pch=19,col="red")
> points(x,y,pch=19,col="red")
> lm(y~x)
Coefficients:
(Intercept)            x
-357.9901       0.2164
que l’on pourrait trouver très proche.

J’ai été un peu surpris de voir cet argument du trimestre gagné avancé lors du débat sur les retraites. Pour commencer, j’ai été surpris de voir Xavier Bertrand prétendre que "Quand vous preniez votre retraite à 60 ans en 1982, vous aviez dix ans d’espérance de vie, aujourd’hui vous avez vingt ans d’espérance de vie. Cette formidable bonne nouvelle, il faut la financer". Vérifions puisqu’on a les codes, il suffit de refaire tourner le programme avec
> age=60
Le graphique est alors

pour les femmes, et pour les hommes, on a

Côté chiffres, pour les femmes,
> exF["1982"]
    1982
22.69534
et pour les hommes,
> exH["1982"]
    1982
17.65394
Bref, j’ai du mal à trouver d’où sortent ces dix ans. Quant à aujourd’hui, pour les femmes,
> exF["2010"]
    2010
26.44607
et pour les hommes,
> exH["2010"]
    2010
21.11194
autrement dit on a dépassé les 20 ans.
Si on calcule la pente de la régression, on obtient pour les femmes,
> lm(y~x)
Coefficients:
(Intercept)            x
-281.0926       0.1533
et pour les hommes,
> lm(y~x)
Coefficients:
(Intercept)            x
-192.7518       0.1063
Autrement dit le gain qui est de l’ordre du trimestre par an à la naissance est plutôt de l’ordre de 2 mois pour les femmes et de 5 semaines pour les hommes. En gros, le gain est divisé par deux entre la naissance et 60 ans. Le graphique ci-dessous montre le gain annuel d’espérance de vie en fonction de l’âge (je suis passé en jours en ordonnées), avec les femmes (en rouge) et les hommes (en bleu)

Bref, si le gain a été très important à la naissance (la surmortalité infantile ayant fortement chuté depuis les années 50), il n’est pas aussi important par la suite, et chute même passé 55 ans ! Autrement dit un bébé vivra - en moyenne - plus longtemps qu’un bébé né en 1950. Mais une personne de 70 ans vivra à peine plus longtemps - en moyenne - qu’une personne de 70 ans en 1950. Bref, parler de ce gain d’un trimestre par an dans un débat sur les retraites (comme cela est fait ici par une députée de Meurthe-et-Moselle) n’a pas de sens.

Foot, probabilité, et assurance, partie 2

Arthur Charpentier — Thu, 03 Jun 2010 10:05:00 +0200

Continuons le billet commencé ici, où j’essayais noté que les enseignes qui se sont engagé à faire des cadeaux en cas de victoire de l’équipe de France s’étaient couvertes via des contrats d’assurance. En fait, elles peuvent également se couvrir sur les marchés (en l’occurrence les sites de paris en ligne). Il existe plusieurs façons de faire des paris. La première est de faire (ou d’acheter) un pari sur un événement (la France ira en finale de la coupe du monde, par exemple). Ce pari est assorti d’une cote que défini le bookmaker. A une cote de 20 contre 1, si je paris 1 euro, je gagne 20 euros en cas de qualification pour la finale.
Dans ce cas particulier, le bookmaker est un agent important dans le marché, car il fixe le prix. Le site oddschecker.com (ici) permet de voir les cotes offertes par différents bookmakers.

Mais il existe aussi des marchés de paris, sans bookmakers. Il s’agit réellement d’un marché, où les parieurs prennent des positions les uns contre les autres, le bookmaker ne servant que d’intermédiaire et d’animateur sur ce marché. C’est le cas sur betfair.com (là) qui est un site d’échange de paris.

Pour commencer simplement, imaginons que nous sommes quelques jours avant une finale opposant la Corée du Sud aux Etats-Unis. Comme il ne peut y avoir qu’un vainqueur il existe des contrats symétrique, où l’on peut gagner 20 euros si la Corée gagner (et rien si elle perd) et un autre où l’on peut gagner 20 euros si la Corée perd. Bref, soit je suis acheter sur la Corée, soit je suis vendeur. Le prix de ces contrats indiquent la probabilité qu’a la Corée de gagner. On parle aussi de back or lay, deux évènements étant possibles.

Lecture d’une cote

Regardons maintenant le match d’ouverture, opposant l’Afrique du Sud au Mexique. Les cotes sont les suivantes, "Mexique 6/4, Afrique du Sud 15/8 et Draw (match nul) 11/5" sur un des sites. Le premier nombre indique ce que l’on gagne (net) en misant le second, si mon pari est gagnant. Autrement dit, si je paris 5 euros sur un match nul, et que match nul il y a, alors mon bénéfice sera de 11 (que je touche en plus de ma mise). La cote est alors de 1+11/5, soit 16/5 (j’avais fait ici un billet sur le lien entre cote et probabilité). Elle indique non plus le bénéfice, mais le gain: en misant 5 euros, j’empoche 16 euros. On parlera de paris 1X2.

Pour aller plus loin, on peut aussi parier sur le nombre de buts marqués. On parlera de paris over/under.

Chez le même bookmaker, si je paris 1 euro que la France gagne la coupe du monde, je fais un bénéfice de 18 euros.

Cote ou probabilité ?

Lorsque j’avais essayé d’expliquer l’origine de la loi logistique (ici), j’avais expliqué qu’il existait un lien simple entre probabilité et cote, la cote étant la probabilité que l’évènement se réalise, divisé par la probabilité qu’il ne se réaliste pas. Malheureusement, les choses sont un peu plus complexes dans les paris...
Zoomons un peu sur le match d’ouverture,

La cote est la valeur offerte la par le bookmarker, à laquelle on ajoute 1. La "probabilité" est alors l’inverse de cette cote.

résultat	Bookie Offers	Cote	"probabilités"
Mexico	11/8	2.375	42.10%
South Africa	13/8	2.623	38.09%
Draw	11/5	3.200	31.25%
		total	111.45%

A-t-on raté quelquechose ? car visiblement la "probabilité" n’est pas une vraie probabilité... En fait, supposer que l’on puisse construire une mesure de probabilité à partir des cotes doit pouvoir être relié au "théorème fondamental d’asset pricing" (et de probabilités risques neutres mentionnés tout à l’heure)

Relecture des cotes à l’aide du théorème fondamental

Autrement dit, on suppose qu’il existe une mesure telle que

où est l’ensemble des gains possibles, et un équivalent certain.
On suppose qu’il existe un ensemble d’évènements (mutuellement exclusif comme on dit, formant une partition de ), un seul pouvant survenir (en l’occurrence l’équipe qui reçoit gagne, ou l’équipe qui reçoit perd, ou il y a match nul). On considère un bookmaker qui accepte de payer si l’évènement i survient à un joueur qui aura parié 1 euro. On pourra légitimement pense que . Notre souhait est de relier ces et , si un lien quelconque pouvait exister. Ce problème avait été soulevé par Ramsey ou de Finetti.
En notant le montant total placé sur l’évènement . On pourra dire que l’on est à l’équilibre si est constant. Autrement dit, le ratio des cotes doit être l’inverse des ratios de montants placés,

On peut ainsi s’intéresser à

qui est la fraction de l’argent qui sera effectivement rétribuée aux parieurs, ou sera la part gardée par le bookmaker. On peut penser qu’un bookmaker fixe , et alors

En fait, k peut être vu comme un taux d’actualisation (l’équivalent d’un en mathématique financière). On peut d’ailleurs noter que si on note le produit

alors la stratégie consistant à parier

sur l’évènement , pour tout les évènements, aura un coût unitaire (j’ai ici constitué mon portefeuille si on revient à la théorie de l’arbitrage). Aussi

En l’occurrence, sur le match d’ouverture de la coupe du monde, on obtient de l’ordre de 90% (l’inverse de la somme de la dernière colonne). Autrement dit, le bookmaker prend de l’ordre de 10%.
A partir de là, on peut en déduire de ces prix d’Arrow Debreu une "probabilité implicite" qui est la probabilité risque neutre, ou la croyance des investisseurs dans la survenance des évènements. Pour cela, on pourrait poser

sauf que dans ce cas, compte tenu de la dernière relation, . On définie alors la probabilité implicite des parieurs sous la forme

On a réussit à construire une mesure de probabilité , qui correspond à l’intuition des parieurs... mais qui n’a rien à voir avec la probabilité réelle. Si l’on compare les bookmakers, on note que les probabilités induites sont très proches,

et si on compare la part prise par le bookmaker sur ce match, on est là aussi sur des ordres de grandeurs comparables,

Justin Wolfers et Eric Zitzewitz, dans Prediction Markets paru dans le Journal of Economic Perspectives de 2004 (ici) expliquent d’ailleurs tout cela très clairement. Je peux aussi renvoyer à un document de travail datant de novembre dernier de Egon Franck, Erwin Verbeek, et Stephan Nüesch (en ligne ici).
J’ai parlé tout à l’heure de probabilité risque neutre, mais cela sous-entend que l’on puisse supposer les marchés efficients, et qu’il n’y ait absence d’opportunité d’arbitrage. C’est précisément ce qu’avait regardé Steven Levitt et Richard Gill ici, dans un vieux papier, évoqué également là. Pour résumer les conclusions, il n’y a pas d’opportunités d’arbitrage sur ces marchés, et il existe des "market makers", c’est à dire des parieurs misant énormément, ce qui accroît considérablement la liquidité de ces marchés. Mais cette conclusion sur l’absence d’opportunité d’arbitrage a été largement remise en cause en pratique. Pour aller plus loin, on peut aussi regarder un article de Moris Eaton et David Freedman, Dutch Book against some ’Objective’ Prior, paru dansBernoulli en 2004, en ligne là. C’est marrant car cette relecture Bayésienne des jeux de pari peut être étudier en pratique, pour la coupe du monde de foot sur worldcup.bayesialab.com (là, merci Avner pour le lien). Je vais continuer de suivre ça lorsque les matchs commençerons, afin de voir les mises à jour des paris au fur et à mesure de l’évolution des matchs... à suivre donc.

Des dangers de la corrélation (de Pearson)

Arthur Charpentier — Wed, 02 Jun 2010 23:52:00 +0200

Suite à mon billet sur les copules empiriques, un ancien collègue m’a fait remarqué que je parlais peu des problèmes de corrélation ou de copules, voire de tail dependence, alors que ça a été mes sujets de prédilection pendant quelques années. Histoire de le faire mentir, je vais faire un billet pour raconter une aventure étrange que j’ai vécue pas plus tard qu’avant hier sur un problème de corrélation.
La corrélation au sens de Pearson est utilisée partout, par exemple sur les séries temporelles quand on regarde les autocorrélations.
Pour montrer un petit exemple et faire un peu réfléchir, considérons un petit exemple simple (qui devrait faire réfléchir les banquiers qui ont à faire face à des problèmes similaires). Considérons deux villes ou passe la Seine, l’une étant - on s’en doute - en amont de l’autre.
Il existe des cartes construites par des hydrologues de vitesse d’avancée moyenne de l’eau.

En tant que statisticien, je me disait qu’il serait possible d’étudier les séries et où serait a priori proche de la valeur indiquée par les hydrologues. Dans un monde idéal, pourrait aussi être un processus stochastique, mais je débute dans les problèmes de changements de temps.

Naturellement, comme tout statisticien qui a pratiqué (et même enseigné) les séries temporelles, je voulais utiliser la corrélation entre les séries, pour trouver la valeur de h optimale. A priori, résoudre

devrait être une bonne idée.
Compte tenu des données dont je dispose (en l’occurrence les points sont séparés par une journée sur la figure ci-dessus), en extrapolant un peu à l’occasion (de manière linéaire), je peux regarder des décalages horaires. Étrangement, j’obtiens la courbe suivante pour l’évolution de la corrélation,

Autrement dit, la corrélation (au sens de Pearson) donne une intuition assez étrange, à savoir qu’il pourrait être intéressant de regarder la série décalée de l’ordre d’une semaine.
En fait, si j’avais regardé la corrélation au sens de Spearman, j’obtiens des résultats beaucoup plus en accord avec l’intuition physique,

Étonnant, non ? Moralité, il faut toujours se méfier de la corrélation !

Qu'est ce que la modélisation économétrique (2) ?

Arthur Charpentier — Wed, 02 Jun 2010 14:10:00 +0200

En début de semaine, Luc Chatel semblait noter "les études les plus récentes indiquent que la diminution des effectifs dans les classes n’a pas d’effet avéré sur les résultats des élèves". Tout a été repris, et détaillé ici, avec d’ailleurs une copie de la fiche donnée aux inspecteur, là.
Je suis un peu surpris car un peu de bon sens (et peut être d’expérience) me laisser penser le contraire. Je me demandais qui avait bien pu conclure à une chose pareil. J’ai traîné un peu sur le net, je suis tombé sur le rapport de Thomas Piketty et Mathieu Valdenaire par exemple (ici pour le rapport et là pour les slides), qui avait fait de l’économétrie sur le panel primaire de 1997, montrant par exmple que la réduction d’un élève par classe en CE1 augmentait de 0,7 points les résultats en maths en début de CE1. Bref, sur le lien effectifs et résultats, j’avais l’impression que l’effet était avéré. Les papiers de Robert Gary Bobo (ici) allaient dans le même sens. Mais qui donc avait pu pondre ces "études les plus récentes" ?

Après avoir cherché un peu, j’ai fini par trouver un élément de réponse sur mon blog (ici): dans un billet sur la modélisation économétrique, j’avais commencé à reprendre le papier de Joshua Angrist et Victor Lavy, "using Maimonides’ rule to estimate the effect the effect of class size on scholastic achivement" (publié que le QJE en 1999, ici). Malheureusement mon billet était long, j’avais fait une régression simple, et conclu à un effet croissant de l’effectif (ici en abscisse) sur les résultats scolaires: plus grande est la classe, meilleurs sont les résulats !

Damned, c’est moi le fautif ? Je devrais me rejouir que Luc lise mon blog (maintenant que je sais qu’il lit mon blog, je peux l’appeler par son prénom), sauf que mon billet se terminait par une mise en garde. Ce que l’on mesure est probablement un effet caché, la mesure d’autre chose.... Afin de justifier ce point, je vais enfin me lancer (avec 10 mois de retard) dans la rédaction de la suite !
En fait, comme le notent Joshua Angrist et Victor Lavy, une variable importante dans l’analyse est le "school level index of student’s socioeconomic status", appelé dans le texte original "percent disadantaged" (PD). Cet indice est suffisement sérieux pour que le Ministère de l’Education l’utilise. Il est fonction du niveau d’éducation du père, et de la taille de la famille de l’enfant scolarisé. On le voit très clairement, plus l’indice est élevé, i.e. plus l’école est situé dans un endroit défavorisé, moins bons sont les résultats scolaires.

Bref, pour ceux qui concluaient après le précédant billet qu’il fallait faire des classes les plus grands possibles pour avoir de meilleurs résultats, ici on pourrait conclure qu’il faut interdire aux enfants de familles nombreuses, ou celles dont le père n’a pas fait de longues études, d’aller à l’école. C’est un peu stupide. Mais on se doute que tout cela est sûrement très lié. Moralité, il faudra faire une étude plus poussée pour voir si, à contexte socio-économique identique, la taille des classes et le résultat aux tests sont positivement ou négativement corrélés... à suivre donc....

Arrêtez donc vos singeries !

Arthur Charpentier — Wed, 02 Jun 2010 00:23:00 +0200

Tout le monde connaît l’expérience où l’on compare les gains financiers obtenus par un traders expérimenté et un singe. Même si beaucoup de sites évoquent "la célèbre expérience", je n ’en ai pas trouvé beaucoup - en France - qui donnent une référénce claire. Cette expérience était née d’une théorie de Burton Malkiel dans son livre a random walk down wall street, qui affirmait "a blindfolded monkey throwing darts at a newspaper’s financial pages could select a portfolio that would do just as well as one carefully selected by experts". En 1988, le Wall Street Journal a pris Burton Malkiel au mot, et a tenté l’expérience.Mais ils n’utilisèrent pas des singes, et ce furent des journalistes qui tiraient au hasard les actions sur lesquelles ils investiraient (alors que parallèlement, les mêmes sommes étaient confiées à des vrais inverstisseurs). Après 6 mois d’expérience, ils firent le bilan. "On October 7, 1998 the Journal presented the results of the 100th dartboard contest. So who won the most contests and by how much? The pros won 61 of the 100 contests versus the darts. That’s better than the 50% that would be expected in an efficient market. On the other hand, the pros losing 39% of the time to a bunch of darts certainly could be viewed as somewhat of an embarrassment for the pros. Additionally, the performance of the pros versus the Dow Jones Industrial Average was less impressive. The pros barely edged the DJIA by a margin of 51 to 49 contests. In other words, simply investing passively in the Dow, an investor would have beaten the picks of the pros in roughly half the contests (that is, without even considering transactions costs or taxes for taxable investors). The pro’s picks look more impressive when the actual returns of their stocks are compared with the dartboard and DJIA returns. The pros average gain was 10.8% versus 4.5% for the darts and 6.8% for the DJIA." Il semble que l’Expansion en France ait aussi voulu tenter l’expérience.

La conclusion de cette expérience est qu’en moyenne, il est dur de battre le hasard ! Une autre conclusion peut aussi être qu’en mettant 500 singes dans une salle de marchés, au bout d’un certain temps on peut en trouver un qui sera aussi riche que Georges Soros, comme le disait Jean-Philippe Bouchaud. On peut parler de grandes déviations éventuellement...
Cette histoire peut prêter à sourire, effectivement, mais un think tank plutôt conservateur outre atlantique (le National Center for Public Policy Research) envisage de faire la même chose pour prédire les ouragans: au lieu d’écouter les climatologues, ils envisagent de faire appel à un singe (comme mentionné ici ou là). Même si les prévisions des chercheurs ne sont pas toujours très fiables (en particulier sur les catastrophes comme les séismes ou les ouragants), elles reposent sur des éléments rationnels, avec souvent des intervalles de confiance.... Oui, je suis le premier à râler quand Météo France prévoit un beau week end et qu’au final il pleut, mais je suis agacé parces démagogues qui critiquent le travail des scientifiques (surtout qu’en l’occurence, je pense que le NCPPR confond météo et climat, mais c’est un autre débat). Il y a un an, la communauté scientifique essayait de faire comprendre qu’elle n’était pas un parasite, voilà que maintenant il faut que l’on se batte pour faire comprendre qu’on fait des choses qu’un chimpanzé ne saurait pas faire ! On est décidément tombé bien bas...

Copules et processus empiriques

Arthur Charpentier — Tue, 01 Jun 2010 13:07:00 +0200

Tarek Zari a soutenu sa thèse au début du mois, présentant une "contribution à l’étude du processus empirique de copule", et sa thèse est en ligne ici. Je mets aussi une copie de ses slides là. Historiquement, il semble que Frits Ruymgaart a été le premier a parler de processus empirique de copules, en 1973 (sa thèse est en ligne ici).

Paul Deheuvels avait également introduit la notion en copule empirique dès 1979 sous le nom de "fonction de dépendance empirique". A la même époque, Ludger Rüschendorf proposait également une étude asymptotique des processus empiriques de copules (ici en 1976), ou encore Gäenssler et Stute dans leur seminar on empirical processes et Winfried Stute dans les années 80 (là). Une revue de la littérature sur les processus empiriques multivariés a été publié à cette époque, en ligne là. Depuis Jean-David Fermanian a publié un papier ici sur la convergence faible, et Paul Deheuvels ou Ludger Rüschendorf ont publié énormément de choses, en particulier sur la vitesse de convergence...

Fibonacci: des lapins à l'actualisation

Arthur Charpentier — Tue, 01 Jun 2010 09:48:00 +0200

Depuis qu’on peut facilement commander des livres d’occasion outre-atlantique, j’avoue en abuser un peu... La dernière commande que j’ai faite est le "origins of value" de Richard Michod, Lynn Nadel et Michael Hechter. Cet ouvrage, qui est une compilation d’articles d’histoire de la finance, contient un chapitre passionnant sur Leonardo Pisano alias Fibonacci. Pour remettre dans le contexte, Fibonacci vivait il y a fort fort longtemps, autour de 1200 (300 ans avant la renaissance et Léonard de Vinci). Je peux renvoyer ici par exemple pour une analyse de quelques uns de ses ouvrages.
Pour moi, Fibonnaci, c’était l’histoire de la reproduction des lapins, c’est à dire, l’étude des suites récurrentes pour ceux qui se souviennent de leurs cours de lycée, avec des suites de la forme

Le problème que l’on cherche à résoudre - posé par Fibonacci - est "partant d’un couple, combien de couples de lapins obtiendrons-nous après un nombre donné de mois sachant que chaque couple produit chaque mois un nouveau couple, lequel ne devient productif qu’après deux mois ? ". Pour ceux qui se souviennent de la résolution, l’équation caractéristique de cette relation de récurrence est un polynôme de degré 2,

dont les solutions sont

Le premier nombre est souvent connu comme étant le nombre d’or (ici). Bref, Fibonacci a été recyclé abondamment par l’utilisation ésotérique de cette racine de polynôme de degré 2. Aussi, lorsque j’ai vu évoqué le nom de Fibonacci dans cet ouvrage de finance, j’ai eu peu peur, d’autant plus que certaines personnes faisant de l’analyse chartiste utilisent ce nombre d’or (comme l’évoque rapidement certains site, comme http://www.forexfibonacci.com/ ou http://www.fibonaccitrader.com/). Aussi, parler de Fibonacci en finance, c’est un peu comme parler de Nostradamus dans les modèles de prévisions.

Que nenni mon ami ! Fibonacci est l’auteur d’un Liber Abaci, correspondant à un livre des calculs. En fait, Leonardo s’est énormément inspiré des mathématiques arabes, alors que la norme était plutôt l’arithmétique romaine. Pour ceux qui ont essayé un jour dans leur vie de lire des nombres latins (même dans Astérix), la logique est assez déroutante ! Par exemple 888 s’écrit DCCCLXXXVIII... Essayez avec cette écriture de poser une addition ! Les nombres arabes - eux - utilisaient le système décimale, ce qui a fait gagné énormément de temps de calcul (on peut aussi remonter à l’arithmétique indienne, mais on sort largement du cadre de mon billet). Fibonacci évoque d’ailleurs constamment l’utilisation de la règle de trois (sans l’appeler ainsi, j’étais revenu ici sur l’origine de cette règle).
Dans son livre d’arithmétique des nombres arabes, Fibonacci montre une application intéressante: le calculs d’actualisation et des marges commerciales. William Goetzmann reprend tout cela dans un article publié en 2003, " Fibonacci and the Financial Revolution" en ligne ici. Il faut aller pour cela dans le chapitre 12 du Liber Abaci, où Fibonacci parle de partage du profit (ce que nous appellerions l’allocation de capital, et que l’on retrouve dans les travaux sur les jeux coopératifs, en particulier les travaux de Shapley), mais aussi d’actualisation de flux futurs, et surtout, de réflexions sur la différence entre un taux annuel et un taux trimestriel. La version en anglaise est la suivante

ou encore

Dans la version originale, cela donne

Quand on pense qu’auparavant l’arithmétique se faisait en chiffre romain, et n’autorisait pas de raisonnements aussi subtils, on peut légitimement penser que Fibonacci est révolutionné le calcul actuariel (certes, en se contentant de traduire des textes arabes). Si on cherche à remonter encore dans le temps, on peut également citer un mathématicien indien, Brahmagupta, ब्रह्मगुप्त, vers 650, qui semble également avoir beaucoup inspiré Fibonacci.

Histoire éthylique, partie 4

Arthur Charpentier — Mon, 31 May 2010 17:05:00 +0200

suite aux pressions générales, je vais reprendre mes discussions d’alcoolique.... ou plutôt reprendre des classiques de finance, en expliquant que ce sont simplement des problèmes que se posent les amateurs de boissons fortes (de là à conseiller plutôt de recruter dans les bars qu’à la sortie des grandes écoles...).
Bref, avant d’avoir entamé sa marche aléatoire dans la rue de la soif (ici, correspondant aux problèmes d’options à barrière traduit en termes financier), puis d’avoir un soucis avec ses clés (là), puis la maréchaussée (ici), notre héros (car on peut maintenant l’appeler un héros après 4 billets qui lui sont consacrés) avait du choisir son bar... Le problème est loin d’être simple. Il y a 20 bars dans la rue (disons pour faire quelque chose de plus formel). Il arrive de la place sainte Anne, et là, il souhaite choisir le bar le moins cher. Le soucis est qu’il n’a pas le droit de faire demi-tour¹ et il ne connaît pas les prix pratiqués dans les différents bars. Il part avec un a priori qui est que le prix d’une pinte est compris entre 3 et 6 euros, que le prix est uniformément réparti entre ces deux prix, et que les prix sont indépendants d’un bar à l’autre. Pour les financiers, il a une option (de commander une bière), et peut l’exercer quand il le souhaite. Une option américaine en quelque sorte. Supposons qu’on soit arrivé au bar . On peut soit payer (qui est supposé aléatoire, uniformément distribué et indépendant des autres bars), soit espérer que l’on puisse payer moins cher plus loin,
Soit la valeur de cette option, alors

i.e.

soit

où désigne la loi du prix de la bière (soit ici une loi uniforme) avec une condition terminale de la forme

car il a soif, et ne quitera pas la rue sans avoir bu un verre !
Classiquement, par backward induction, on peut résoude ce programme, à partir de la loi de . Posons . Alors

soit simplement

soit enfin

Je laisse les plus courageux simplifier les calculs. La "frontière d’exercice" est alors obtenue par récurence. Numériquement, le code est alors simplement
> n=20
> u=rep(NA,n)
> b=6;a=3
> u[n]=(b+a)/2
> for(k in (n-1):1){
+ u[k]=1/(b-a)*(u[k+1]*(b-u[k+1])+(u[k+1]^2-a^2)/2)
+ }

Dès qu’on atteind la barrière, on s’assoit au bar. On note que plus on avance dans la rue, moins on est exigent: au tout début, on ne s’assoit pas à moins de 3 euros 30... mais plus on avance, plus on relève le seuil d’exigence. Le calcul sous forme intégrale donne ici
> u=rep(NA,n)
> b=6;a=3
> u[n]=(b+a)/2
> for(k in (n-1):1){
+ g=function(x){pmin(x,rep(u[k+1],length(x)))/(b-a)}
+ u[k]=integrate(g,lower=a,upper=b)$value
+ }
J’avais déjà abordé ce problème dans un précédant billet, sur les options américaines, mais on peut maintenant aller un peu plus loin... que se passe-t-il si on suppose que les prix sont discret (par exemple par tranches de 50 centimes ou 1 euro) ? L’avantage avec ces méthodes numériques est que l’on peut très facilement enlever des hypothèses, par exemple ici on aurait
> h=2
> K=(b-a)*h+1
> PRIX=seq(a,b,by=1/h)
> u2=rep(NA,n)
> b=6;a=3
> u2[n]=(b+a)/2
> for(k in (n-1):1){
+ g=function(x){pmin(x,rep(u[k+1],length(x)))}
+ u2[k]=sum(g(PRIX)*1/K)}
pour des seuils à 1 euros (les seuls prix possibles étant 3,4,5 ou 6 euros).

Ou la frontière suivante si les prix varient par tranche de 50 centimes.

Compte tenu de la discrétisation, notons que la vraie frontière devient alors ici

Bref, comme toujours, les problèmes d’alcooliques rejoignent les problèmes d’exercice optimal d’options américaines, problème classique en finance de marché...
¹ pour rendre cette histoire crédible, à chaque bar rencontré il demande le prix d’une pinte. S’il estime que c’est trop cher, il s’exclame "mais c’est bien trop cher ici !" et s’en va. Sinon il commande et s’installe. Cette exclamation rend improbable - à ses yeux - l’idée de revenir finallement s’intaller au bar....