Digital resources in the Social Sciences and Humanities OpenEdition Our platforms OpenEdition Books OpenEdition Journals Hypotheses Calenda Libraries OpenEdition Freemium Follow us

Les députés sont-ils à l’image de la population

Beaucoup de choses ont été écrites sur le fait que les députés ne sont pas vraiment le reflet de la population, que ce soit en terme de profession, de sexe, d’origine, d’age, etc. La liste pourrait être longue. Il y a plusieurs mois, j’avais commencé à regarder le profil des députés, par age. En effet, le site http://assemblee-nationale.fr/ permet d’accéder à des données sur tous les députés, depuis la Révolution. Y compris leur date de naissance. En croisant ces données avec des données de population, par exemple via http://www.mortality.org/, on peut comparer la répartition des ages des députés, avec la répartition des ages de la population.

Pour les amateurs, le code pour récupérer les données (ou au moins les dates de naissance des députés) ressemble à

N=2002
URL=paste("http://www.assemblee-nationale.fr/
sycomore/result.asp?radio_dept=tous_departements&
regle_nom=est&Nom=&departement=&
choixdate=intervalle&D%C3%A9butMin=01%2F01%2F",
N,"&FinMin=31%2F12%2F",N,"&Dateau=&legislature=",
s,"&choixordre=chrono&Rechercher=
Lancer+la+recherche",sep="")
HTML=scan(URL,what="character")

k=which(HTML=="class=\"titre\">Né")
vHTML=HTML[k:length(HTML)]
vk=which(substr(vHTML,1,7)=="> ")
liste=vHTML[vk]
naissance=liste[seq(1,length(liste),by=2)]
NAISSANCE=as.Date(substr(naissance,8,17),
"%d/%m/%Y")

Maintenant, pour être tout à fait honnête, je ne suis pas certain de ce qui est vraiment renvoyé, et j’ai des doutes que cela correspondent réellement à la requête faite. En effet, même si je demande à avoir la liste des députés après l’élection, j’ai trop de monde… mais peut-être est-ce du aux décès éventuels, et il est possible que l’ensemble des députés qui ont siégé pendant la mandature apparaissent dans le résultat de la requête.

Sur la figure suivante on voit, sur plusieurs élections depuis plus de 100 ans, comment les deux distributions se déforment, avec en rouge la distribution de l’age des députés, et en bleu, la distribution de la population française, dans son ensemble (population de plus de 18 ans)

Si on veut tout suivre sur un graphique, au lieu de se regarder une animation, on peut représenter les différents quantiles (10%, 25%, 75% et 90%, retenus sur la population de plus de 18 ans, et l’age médian, au centre), avec la population française l’année de l’élection,

et l’ensemble des élus au parlement,

Si on veut faciliter la comparaison, on peut se contenter de visualiser l’évolution des ages moyens,

ou encore, du ratio (en % de différence) entre l’age moyen des députés, et celui de l’ensemble de la population.

Sur ce graphique, on voit que depuis 30 ans, l’age moyen des députés croit plus vite que celui de la population: la population français vieilli, mais moins que ses députés… La gérontocratie perdure donc en France. En espérant que cela ne débouche pas sur le clash générationnel que l’on semble observer ces temps-ci au Québec…

A la recherche des groupes de trolls sur Twitter

Il y a quelques jours @olihb avait mis en ligne sur son blog http://olihb.com/ une très jolie carte permettant de visualiser les personnes (ou les comptes) impliquées sur Twitter sur les manifestations récentes au Québec (la carte était basée sur gephi). Beaucoup de monde twittait avec des hashtags, et l’idée était d’utiliser les plus classiques et les plus populaires pour identifier les comptes actifs (#ggi#manifencours, etc). Mais la semaine passée j’avais l’impression que le débat s’était durci, avec beaucoup plus de trolls dans les discussions.

Bref, le principe est que les trolls évoluent souvent en petits groupes, donc avec @3wen on a voulu voir si on pouvait identifier facilement les groupes de trolls (sans pour autant faire du sentiment analysis, juste visualiser des personnes qui se suivent beaucoup entre elles, les connivences entre twittos, i.e. du suivi mutuel).

La méthodologie est simple. On a deux bases à notre disposition. Un première contient tous les tweets contenant #ggi entre le 7 juin à 1 heure du matin et le 8 à minuit 45 (soit 14,700 tweets). La seconde contient tous les tweets contenant #ggi entre le 12 juin à 23 heures 30 et le 14 juin à 17 heures (soit 14,784 tweets). Ah oui, les heures sont GMT. Plus précisément, on a la fréquence d’arrivée suivante

Les durées d’observations sont courtes, mais l’idée est que si une personne a tweetté au moins une fois pendant ces périodes, on les a dans notre base, et ça suffit pour les suivre. On ne va pas étudier l’intensité de la mise en ligne de tweets. Je ne rentre pas dans le détail du code, mais c’est basé sur les fonctions présentées dans un précédant billet, et toujours avec la même librairie,

require("RJSONIO")

Sur des deux bases, i.e. 29,484 tweets, on récupère un ensemble de 6,556 twittos uniques (mais répartis dans deux listes: ceux qui ont tweeté le 7 et ceux qui ont tweeté le 14). Pour toutes ces personnes, on est allé voir qui les suivaient. Sur le principe, c’est simple, mais on fonctionne avec des API qui imposent des limites horaires. Bref, il faut bricoler. Pour une personne, on récupère avec la fonction suivante la liste de ses followers,

recuperefollowers=function(id,nbrequetes){
followers=try(scan(paste(
"http://api.twitter.com/1/followers/ids.json?cursor=-1&user_id=",
id,sep=""),what = "character", encoding="latin1"))
if(is.null(attr(followers,"condition"))){
followers=paste(followers[1:length(followers)],collapse=" ")
followers=fromJSON(followers, method = "C")
id_followers=followers$id
next_cursor=followers$next_cursor_str
nbrequetes=nbrequetes+1
plusDe5000=FALSE
while(nbrequetes<150 & next_cursor!="0" &
is.null(attr(followers,"condition"))){
plusDe5000=TRUE
followers=scan(paste(
"http://api.twitter.com/1/followers/ids.json?cursor=",
next_cursor,"&user_id=",id,sep="")
,what = "character", encoding="latin1")
followers=paste(followers[1:length(followers)],collapse=" ")
followers=fromJSON(followers, method = "C")
id_followers=c(id_followers,followers$id)
next_cursor=followers$next_cursor_str
nbrequetes=nbrequetes+1
}
if(plusDe5000 & nbrequetes>=149){
return(c(list(),nbrequetes,TRUE))
}else{
return(c(list(id_followers),nbrequetes,FALSE))
}}else{
nbrequetes=nbrequetes+1
print(paste("Probleme rencontre avec ",id,sep=""))
return(c(list("PROBLEME"),nbrequetes,FALSE))}}

On peut ensuite lancer cette fonction sur tous les identifiants qu’on a récupéré, et on stocke tout le monde dans une (petite) base

temp=try(recuperefollowers(
lesIDAParcourir[unID_index],nbrequetes))
lesFollowers=temp[[1]]
resul=cbind(rep(lesIDAParcourir[unID_index]
,length(lesFollowers)),lesFollowers)
nbrequetes=temp[[2]]
refaire=temp[[3]]
Sys.sleep(runif(1,1,1.5))
compte=compte+1
write.table(resul,paste("auto_resul_",
compte,".txt",sep=""),sep=";")

Si on veut que ça tourne, il faut juste faire patienter un heure après avoir fait 150 requêtes. Bref, dans la boucle, on met une petite fonction qui temporise,

if(nbrequetes>148){
nbrequetes=1
Sys.sleep(60*60+60*trunc(runif(1,2,3))) }

Une fois créées les 6,500 bases de nos 6,500 comptes, contenant la liste de tous les followers, on va les agréger (fort heureusement, on n’avait pas les  gros comptes avec plus d’un million de followers). La liste se récupère avec

N=list.files("id/petits")

D’ailleurs, pour recoller les 6,500 bases ensemble, Ewen suggérait le code suivant, beaucoup beaucoup plus rapide qu’une gros boucle,

recupere_petits_fichiers=function(x){
temp=read.table(paste("id/petits/",x,sep=""),
header=TRUE,sep=";",comment.char="",
check.names=FALSE,colClasses=
c("character","character"))
colnames(temp)=list("id_twittos","id_follower")
head(temp)

On fait ensuite une récupération des followers pour la première base de tweets, et une autre pour ceux de la seconde. Les deux listes sont obtenues avec le code suivant

present_avant_temp=temp[which(
temp$id_twittos%in%lesID_avant),]
if(nrow(present_avant_temp)>0){
present_avant_temp=
present_avant_temp[which(
present_avant_temp$id_follower%in%lesID_avant),]}
else{
present_avant_temp=NULL}

present_apres_temp=
temp[which(temp$id_twittos%in%lesID_apres),]
if(nrow(present_apres_temp)>0){present_apres_temp=
present_apres_temp[
which(present_apres_temp$id_follower%in%lesID_apres),]}
else{present_apres_temp=NULL}
return(c(list(present_avant_temp),
list(present_apres_temp),
list(nbfollowers=nrow(temp))))

Pour constituer notre grosse base, la ruse est d’utiliser non pas une boucle sur tous les petites bases, mais

resul=lapply(N,recupere_petits_fichiers)
liens_avant=data.frame(do.call(
"rbind",lapply(resul,function(x) x[[1]])),
stringsAsFactors=FALSE)
liens_apres=data.frame(do.call(
"rbind",lapply(resul,function(x) x[[2]])),
stringsAsFactors=FALSE)

On a ainsi récupéré une liste avec plusieurs centaines de milliers de comptes Twitter. On a retenu seulement ceux qui avaient twitté avec le hashtag #ggi. On a ainsi pu constituer, pour nos 6,556 twittos une base de personnes qui tweettent, et de personnes qui les suivent (et qui ont twitté pendant cette – courte – période avec le hastag #ggi). Bref, une grosse matrice de correspondance de qui suit qui. On a pu le faire sur les deux bases. C’est précisément la différence entre les deux bases que l’on a cherché à visualiser, et gephi semble permettre de comparer les bases dans le temps. A gauche, on a le nuage de connexions entre les comptes pour la première période, et à droite pour la seconde.

Bon, le premier est plus compact, mais on voit qu’il y a des différences, en particulier sur les bords, que ce soit en haut

ou à droite du nuage,

Dans les deux cas, on voit que plein de petits comptes sont venus se greffer au nuage, souvent en suivant quelqu’un qui le suit en retour. Et parfois, ce sont les seuls connexions qu’ils ont (qui utilisent le hashtag #ggi). Ce sont ces comptes qui pourraient constituer les nouveaux trolls qu’on espérait identifier. L’avantage est qu’on a les noms de tous ces comptes (on ne les a pas mis sur le graphique), mais dans un second temps, on va aller regarder ce qu’on fait précisément ces comptes. Peut-être interprétons nous mal ces graphiques, qui représentent une arrivée de nouveaux comptes, et que cette dynamique est la même pour tous les nouveaux comptes, qui commencent par de faibles connexions, et ensuite vont s’étoffer au fur et à mesure. Ce qu’il faut que l’on comprenne aussi c’est pourquoi ils sont tous aussi rapprochés. Ils semblent être suivis par des personnes qui sont suivis pas les mêmes personnes. Ce qui laisse à penser que l’on a effectivement identifié une certaine catégorie de comptes.
On peut reprendre les graphs avec une résolution plus fine, avant

et après

ou alors en changeant les options, avant

à comparer avec après

On voit peut être encore plus clairement sur ces deux dernières images les nouveaux comptes se greffer sur le graph existant. Si la qualité ne suffit pas, des graphiques (très) haute résolution sont téléchargeables (15Mo chacun), avant etaprès (avec les identifiants cette fois). Sinon, en promenant la souris sur le dessin dessous, on voit encore mieux la différence,
 

Maintenant – pour être complètement honnête – sur la lecture des graphs, j’avoue ne pas avoir tout bien compris. Les représentations semblent claires, mais j’aurais bien aimé être certain d’avoir compris la construction proposée par gephi avant de les mettre en ligne. Car si on regarde sur moins ce compte, on comprend que l’algorithme permet effectivement de bien mettre en avant les réseaux. Parmi les illustrations les plus frappantes, la figure suivante montre des interconnections avec Facebook à gauche, et purement aléatoire à droite (illustration trouvé sur le site de gephi),

La grande difficulté est alors de représenter spatialement les groupes, comme le notent Gastner et Newman dans The Spatial Structure of Networks, Watts et Strogatz dansCollective dynamics of ‘small-world’ networksou encore Nisha et Venkatesh dans Small worlds: How and why. Mais pour etre tout à fait complet, je renverrais aussi vers la thèse de doctorat de Jure Leskovec, intituléeDynamics of large networks ainsi qu’au livre (intégralement téléchargeable) Networks, Crowds, and Markets: Reasoning About a Highly Connected World de David Easley et Jon Kleinberg. Bref, on a découvert un outils probablement très riche, mais il va nous falloir du temps pour comprendre ce qu’il fait vraiment.

Pricing options on multiple assets

I am a big fan of trees. It is a very nice way to see how financial pricing works, for derivatives. An with a matrix-based language (R for instance), it is extremely simple to compute almost everything. Even options multiple assets. Let us see how it works. But first, I have to assume that everyone knows about trees, and risk neutral probabilities, and is familiar with standard financial derivatives. Just in case, I can upload some old slides of the first course on asset pricing we gave a few years ago at École Polytechnique.

Let us get back on the pricing of (European) call options, with trees.The idea is simple. We have to fix the number of periods. Let us start with only one (as described in the slides above). The stock has price and can go either up, and then have price or go down, and have price . And the fundamental theorem of asset pricing says that we do not really care about probabilities of going up, or down. Assuming that we can buy or sell that stock, and that a risk free asset is available on the market, it is possible to price any contingent financial product, like a financial option. Since we know the final value of the option when the stock goes either up, or down, it is possible to replicate the payoff of that option using the stock and the risk free asset. And we can prove that the price of the option is simply

where the probability is the so-called risk neutral probability

So, we’ve done it here with only one single period, but it is possible to extend it to multiperiods. The idea is to keep that multiplicative representation of possible values of the stock, and to get a recombinant tree. At step 2, the stock can take only three different values: went up twice, went down twice, or went up and down (or the reverse, but we don’t care: this is the point of recombining). If we write things down, then we can prove that

for some probability parameter (the so-call risk neutral probability, if it is unique). But we do not really care about those closed formula, the goal is to write an algorithm which computes the tree, and return the price of a call option (say). But before starting, we have to make a connection between that model with up and down prices, and the parameters of the Black-Scholes diffusion, for the stock price. The idea is to identify the first and the second moment, i.e.

(where, under the risk neutral probability, the trend is the risk free rate) and

The code might look like that

n=5; T=1; r=0.05; sigma=.4;S=50;K=50
price=function(n){
u.n=exp(sigma*sqrt(T/n));
d.n=1/u.n
p.n=(exp(r*T/n)-d.n)/(u.n-d.n)
SJ=matrix(0,n+1,n+1)
SJ[1,1]=S
for(i in(2:(n+1)))
{for(j in(1:i)){SJ[i,j]=S*u.n^(i-j)*d.n^(j-1)}}
OPT=matrix(0,n+1,n+1)
OPT[n+1,]=(SJ[n+1,]-K)*(SJ[n+1,]>K)
for(i in(n:1))
{for(j in(1:i)){OPT[i,j]=exp(-r*T/n)*(OPT[i+1,j]*p.n+
(1-p.n)*OPT[i+1,j+1])}}
return(OPT[1,1])
}

We can plot the evolution of the price, as a function of the number of time periods (or subdivision of the time interval, from now till maturity of the European option),

N=10:400
V=Vectorize(price)(N)
plot(N,V,type="l")

Note that we can compare with the Black-Scholes price of this call option, given by

where

and

d1=1/(sigma*sqrt(T))*(log(S/K)+(r+sigma^2/2)*T)
d2=d1-sigma*sqrt(T)
BS=S*pnorm(d1)-K*exp(-r*T)*pnorm(d2)
abline(h=BS,lty=2,col="red")

The code is clearly not optimal, but at least, we see what’s going on. For instance, we do not need a matrix when we calculate using backward recursions the price of the option. We can just keep a single vector. But this matrix is nice, because we can use it to price American options. For instance, with the code below, we compare the price of an American put option, and the price of European put option.

price.american=function(n,opt="put"){
u.n=exp(sigma*sqrt(T/n)); d.n=1/u.n
p.n=(exp(r*T/n)-d.n)/(u.n-d.n)
SJ=matrix(0,n+1,n+1)
SJ[1,1]=S
for(i in(2:(n+1)))
{for(j in(1:i)) {SJ[i,j]=S*u.n^(i-j)*d.n^(j-1)}}
OPTe=matrix(0,n+1,n+1)
OPTa=matrix(0,n+1,n+1)
if(opt=="call"){
OPTa[n+1,]=(SJ[n+1,]-K)*(SJ[n+1,]>K)
OPTe[n+1,]=(SJ[n+1,]-K)*(SJ[n+1,]>K)
}
if(opt=="put"){
OPTa[n+1,]=(K-SJ[n+1,])*(SJ[n+1,]<K)
OPTe[n+1,]=(K-SJ[n+1,])*(SJ[n+1,]<K)
}
for(i in(n:1))
{
for(j in(1:i))
{if(opt=="call"){
OPTa[i,j]=max((SJ[i,j]-K)*(SJ[i,j]>K),
exp(-r*T/n)*(OPTa[i+1,j]*p.n+
(1-p.n)*OPTa[i+1,j+1]))}
if(opt=="put"){
OPTa[i,j]=max((K-SJ[i,j])*(K>SJ[i,j]),
exp(-r*T/n)*(OPTa[i+1,j]*p.n+
(1-p.n)*OPTa[i+1,j+1]))}

OPTe[i,j]=exp(-r*T/n)*(OPTe[i+1,j]*p.n+
(1-p.n)*OPTe[i+1,j+1])}}
priceop=c(OPTe[1,1],OPTa[1,1])
names(priceop)=c("E","A")
return(priceop)}

It is possible to compare those price, obtained on trees, with prices given by closed (approximated) formulas.

> d1=1/(sigma*sqrt(T))*(log(S/K)+(r+sigma^2/2)*T)
> d2=d1-sigma*sqrt(T)
> (BS=-S*pnorm(-d1)+K*exp(-r*T)*pnorm(-d2)  )
[1] 6.572947
> N=10:200
> M=Vectorize(price.american)(N)
> plot(N,M[1,],type='l',col='blue',ylim=range(M))
> lines(N,M[2,],type='l',col='red')
> abline(h=BS,lty=2,col='blue')
> library(fOptions)
> (am=BAWAmericanApproxOption(TypeFlag =
+ "p", S = S,X = K, Time = T, r = r,
+ b = r, sigma =sigma)@price)
[1] 6.840335
> abline(h=am,lty=2,col='red')

Another great thing with trees, is that it becomes possible to plot to region where it is optimal to exercise our right to sell the stock.

Let us move now to a model with two assets, as suggested by Rubinstein (1994). First, observe that a discretization of two independent Brownian motions will be based on two independent random walk, taking values

i.e. both went up (NW), both went down (SE), and one went up while the other went down (either NE or SW). With independent and symmetric random walks, the probabilities will be respectively 1/4. An if we move one step foreward, we have the following tree.

Here it is still recombining. But the size will increase much faster than in the univariate case. Now, assume that there might be some correlation. Then one can consider the following values, to have a specific correlation,

And again, the idea is then to identify the first two moments. This gives us the following system of equations for the four respective (risk neutral) probabilities

For those willing to do the maths, please do. The answer should be

and for the last one

The code here looks like that

price.spead=function(n){
T=1; r=0.05; K=0
S1=105
S2=100
sigma1=0.4
sigma2=0.3
rho=0.5
u1.n=exp(sigma1*sqrt(T/n)); d1.n=1/u1.n
u2.n=exp(sigma2*sqrt(T/n)); d2.n=1/u2.n

v1=r-sigma1^2/2; v2=r-sigma2^2/2
puu.n=(1+rho+sqrt(T/n)*(v1/sigma1+v2/sigma2))/4
pud.n=(1-rho+sqrt(T/n)*(v1/sigma1-v2/sigma2))/4
pdu.n=(1-rho+sqrt(T/n)*(-v1/sigma1+v2/sigma2))/4
pdd.n=(1+rho+sqrt(T/n)*(-v1/sigma1-v2/sigma2))/4
k=0:n
un=matrix(1,n+1,1)
SJ= (S1 * d1.n^k * u1.n^(n-k-1)) %*% t(un) -
un %*%t(S2 * d2.n^k * u2.n^(n-k-1))
OPT=(SJ)*(SJ>K)
for(k in(n:1))
{
OPT0=matrix(0,k,k)
for(i in(1:k))
{
for(j in(1:k))
{OPT0[i,j]=(OPT[i,j]*puu.n+OPT[i+1,j]*pdu.n+
OPT[i,j+1]*pud.n+OPT[i+1,j+1]*pdd.n)*exp(-r*T/n)}}
OPT=OPT0}
return(OPT[1,1])}

If we look at the details, consider two periods, like on the figure above, the are nine values for the spread,

> n=2
> SJ
[,1]      [,2]       [,3]
[1,]  32.02217  84.86869 119.443578
[2,] -47.84652   5.00000  39.574891
[3,] -93.20959 -40.36308  -5.788184

and the payoff of the option is here

> OPT
[,1]     [,2]      [,3]
[1,] 32.02217 84.86869 119.44358
[2,]  0.00000  5.00000  39.57489
[3,]  0.00000  0.00000   0.00000

So if we go backward of one step, we have the following square of values

> k=n
> OPT0<-matrix(0,k,k)
> for(i in(1:k))
+ {
+   for(j in(1:k))
+   {
+     OPT0[i,j]=(OPT[i,j]*puu.n+OPT[i+1,j]*pdu.n+
+ OPT[i,j+1]*pud.n+OPT[i+1,j+1]*pdd.n)*exp(-r*T/n)
+ }
+ }
> OPT0
[,1]      [,2]
[1,] 22.2741190 58.421275
[2,]  0.5305465  5.977683

The idea is then to move backward once more,

> OPT=OPT0
> OPT0<-matrix(0,k,k)
> for(i in(1:k))
+ {
+   for(j in(1:k))
+   {
+     OPT0[i,j]=(OPT[i,j]*puu.n+OPT[i+1,j]*pdu.n+
+ OPT[i,j+1]*pud.n+OPT[i+1,j+1]*pdd.n)*exp(-r*T/n)
+ }
+ }
> OPT0
[,1]
[1,] 16.44106

Here calculations are much (much) longer,

> price.spead(250)
[1]  15.66496

and again, it is possible to use standard approximations to compare that price with a more standard one,

> (sp=SpreadApproxOption(TypeFlag =
+ "c", S1 = 105, S2 = 100, X = 0,
+ Time = 1, r = .05, sigma1 = .4,
+ sigma2 = .3, rho = .5)@price)
[1]  15.65077

Well, playing with trees is nice, but it might not be optimal for complex products. Next time, we’ll discuss other techniques…

Date of death, birthday and Elvis Presley

10 days ago, a study published on http://www.annalsofepidemiology.org/ mentioned that “Death has a preference for birthdays” (as claimed in the title). The conclusion of the paper is that, in general, birthdays do not evoke a postponement mechanism but appear to end up in a lethal way more frequently than expected (“anniversary reaction”). Well, this is not new, and several previous articles have mentioned that point, e.g. Angermeyer et al. (1987).

I found the idea interesting since in demography, there is a large literature trying to extrapolate death rates from discrete to continuous time. Extrapolation are usually extremely smooth. But none of them integrate that aspect of mortality precisely on the birthday. The problem is that it is rather difficult to say something since datasets with individual observations are rare, online.

But yesterday, @coulmont sent me a tweet mentioning a website. I do not know if this is legal (even if some explanations are given), but I will mention courtesy of http://ssdmf.info/. It is a so-called Social Security Death Master File, containing individual informations about deaths in the US, as well as geographic information (as described on http://www.ssa.gov/), for people having a social security number.

With R, it is possible to work on those files (even they are huge, with tens of millions observations). For instance, we can check who is inside.

> elvis=scan("ssdm2",skip=22371720,n=1,what="character",sep=",")
> elvis
[1] " 409522002PRESLEY         ELVIS     0800197701081935  "

If you believe that Elvis is dead, you might agree that this database can be accurate (or at least, not too bad). And further, we can see here how to read the result: Elvis was born on January 8, 1935 (8 last digits), and died on August 16, 1977 (8 digits before). Obviously here, there are some problems with the dataset (we do not have the day of the death of Elvis). So here, we remove all the observations that do not give us proper dates. Then, the idea is to assume that the person died in 2000 (or any year since the point is to focus on days and months). Then, we count the number of days between the day of death and the birthday in 2001 (that would have been after) and the one in 2000 (that was either before or after the death), so that we can derive the number of days after the birthday,

dates=substr(base,66,81)
death=as.Date(substr(dates,1,8),"%m%d%Y")
birth=as.Date(substr(dates,9,16),"%m%d%Y")
indice=is.na(death)|is.na(birth)
mean(indice)
mdeath=substr(dates,1,2)
ddeath=substr(dates,3,4)
mbirth=substr(dates,9,10)
dbirth=substr(dates,11,12)
indice=which(ddeath!="00")
birth1=as.Date(paste(mbirth[indice],
dbirth[indice],"2000",sep=""),"%m%d%Y")
birth2=as.Date(paste(mbirth[indice],
dbirth[indice],"2001",sep=""),"%m%d%Y")
death=as.Date(paste(mdeath[indice],ddeath[indice],
"2000",sep=""),"%m%d%Y")
k=length(indice)
diffday=cbind((as.numeric(death-birth1))[1:k],
(as.numeric(death-birth2))[1:k])
DIFF=apply(diffday,1,function(x) {min(x[x>=0])})

What we have here is the number of days following the previous birthday. If we look at the distribution of that number of days, we obtain

counts=table(DIFF)
plot(as.numeric(names(counts)),
as.numeric(counts))
counts["0"]/(mean(counts[100:200]))
> counts["0"]/(mean(counts[100:200]))
0
1.121261

Thus, the death excess on the day of birth was around 12%, which is rather close to the one obtained from the Swiss mortality statistics 1969–2008 (in Ajdacic-Gross et al. (2012)). Note that here, we just play with a small subset of the entire dataset,

That database is probably extremely interesting, except that it suffers a huge selection bias, since only dead people are in that database. So it might be useless if we wish to study life expectancy of people named Bill versus people named Georges (that was something I wanted to investigate initially). But we’ll see what else we can do with it (since Ewen have been able to write some code to go through that huge dataset).

Do you still have time to sleep ?

Last week, @3wen (Ewen) helped me to write nice R functions to extract tweets in R and build datasets containing a lot of information. I’ve tried a couple of time on my own. Once on tweet contents, but it was not convincing and once on the activity on Twitter following an event (e.g. the death of someone famous). I have to admit that I am not a big fan of databases that can be generated using standard function to study tweets. For instance, we can only extract tweets, notre-tweets (which is also an important indicator of tweet-activity). @3wen suggested to use

require("RJSONIO")

The first step is to extract some information from a tweet, and store it in a dataset (details can be found on https://dev.twitter.com/)

obtenir_ligne <- function(unTweet){
date_courante=unTweet$created_at
id_courant=unTweet$id_str
text=unTweet$text
nb_followers=unTweet$user$followers_count
nb_amis=unTweet$user$friends_count
utc_offset=unTweet$user$utc_offset
listeMentions=unTweet$entities$user_mentions
return(c(list(c(id_courant,date_courante,text,
nb_followers,nb_amis,utc_offset)),
list(do.call("rbind",lapply(listeMentions,
function(x,id_courant) c(id_courant,
x$screen_name),unTweet$id_str)))))
}

Now that we  have the code to extract information from one tweet, let us find several tweets, from one user, say my account,

nom="Freakonometrics"

The (small) problem here, is that we have a limitation: we can only get 100 tweets per call of the function

n=100
tweets_courants=scan(paste(
"http://api.twitter.com/1/statuses/user_timeline.json?
include_entities=true&include_rts=true&screen_name=
",nom,"&count=",n,sep=""),what = "character",
encoding="latin1")
tweets_courants=paste(tweets_courants[
1:length(tweets_courants)],collapse=" ")
tweets_courants=fromJSON(tweets_courants,
method = "C")

Then, we use our function to build a database with 100 lines,

extracTweets <- lapply(tweets_courants,
obtenir_ligne)
mentions=do.call("rbind",lapply(extracTweets,
function(x) x[[2]]))
colnames(mentions)=list("id","screen_name")
res=t(sapply(extracTweets,function(x) x[[1]]))
colnames(res) <- list("id","date","text",
"nb_followers","nb_amis","utc_offset")

The idea then is simply to use a loop, based on the latest id observed

dernier_id=tweets_courants[[length(
tweets_courants)]]$id_str

So, here we go,

compteurLimite=100

while(compteurLimite<4100){
tweets_courants=scan(paste(
"http://api.twitter.com/1/statuses/user_timeline.json?
include_entities=true&include_rts=true&screen_name=
",nom,"&count=",n,"&max_id=",dernier_id,sep=""),
what = "character", encoding="latin1")
tweets_courants=paste(tweets_courants[
1:length(tweets_courants)],collapse=" ")
tweets_courants=fromJSON(tweets_courants,
method = "C")

extracTweets <- lapply(tweets_courants[
2:length(tweets_courants)],obtenir_ligne)
mentions=rbind(mentions,do.call("rbind",
lapply(extracTweets,function(x) x[[2]])))
res=rbind(res,t(sapply(extracTweets,function(x) x[[1]])))
t(sapply(extracTweets,function(x) x[[1]]))
dernier_id=tweets_courants[[length(
tweets_courants)]]$id_str
compteurLimite=compteurLimite+100
}

resFreakonometrics=res=
data.frame(res,stringsAsFactors=FALSE)

All the information about my own tweets (and re-tweets) are stored in a nice dataset. Actually, we have even more, since we have extracted also names of people mentioned in tweets,

mentionsFreakonometrics=
data.frame(mentions)

We can look at people I mention in my tweets

gazouillis=sapply(split(mentionsFreakonometrics,
mentions$screen_name),nrow)
gazouillis=gazouillis[order(gazouillis,
decreasing=TRUE)]

plot(gazouillis)
plot(gazouillis,log="xy")
> gazouillis[1:20]
tomroud freakonometrics       adelaigue       dmonniaux
155              84              77              56
J_P_Boucher         embruns      SkyZeLimit        coulmont
42              39              35              31
Fabrice_BM            3wen          obouba          msotod
31              30              29              27
StatFr     nholzschuch        renaudjf        squintar
26              25              23              23
Vicnent        pareto35        romainqc        valatini
23              22              22              22

If we plot those frequencies, we can clearly observe a standard Pareto distribution,

Now, let us spend some time with dates and time of tweets (it was the initial goal of this post)… One more time, there is a (small) technical problem that we have to deal with: language. We need a function to convert date in English (on Twitter) to dates in French (since I have a French version of R),

changer_date_anglais <- function(date_courante){
mois <- c("Jan","Fév", "Mar", "Avr", "Mai",
"Jui", "Jul", "Aoû", "Sep", "Oct", "Nov", "Déc")
months <- c("Jan", "Feb", "Mar", "Apr", "May",
"Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
jours <- c("Lun","Mar","Mer","Jeu",
"Ven","Sam","Dim")
days <- c("Mon","Tue","Wed","Thu",
"Fri","Sat","Sun")
leJour <- substr(date_courante,1,3)
leMois <- substr(date_courante,5,7)
return(paste(jours[match(leJour,days)]," ",
mois[match(leMois,months)],substr(
date_courante,8,nchar(date_courante)),sep=""))
}

So now, it is possible to plot the times where I am online, tweeting,

DATE=Vectorize(changer_date_anglais)(res$date)
DATE=sapply(resSkyZeLimit$date,
changer_date_anglais,simplify=TRUE)

DATE2=strptime(as.character(DATE),
"%a %b %d %H:%M:%S %z %Y")
lt= as.POSIXlt(DATE2, origin="1970-01-01")
heure=lt$hour+lt$min/60
plot(DATE2,heure)

On this graph, we can see that I am clearly not online almost 6 hours a day (or at least not on Twitter). It is possible to visualize more precisely the period of the day where I might be on Twitter,

hist(heure,breaks=0:24,col="light green",proba=TRUE)
X=c(heure-24,heure,heure+24)
d=density(X,n = 512, from=0, to=24,bw=1)
lines(d$x,d$y*3,lwd=3,col="red")

or, if we want to illustrate with some kind of heat plot,

Note that we did it for my Twitter account, but we can also run the code on (almost) anyone on Twitter. Consider e.g. @adelaigue. Since Alexandre is tweeting in France, we have to play with time-zones,

res=extractR("adelaigue")
DATE=Vectorize(changer_date_anglais)(res$date)
DATE2=strptime(as.character(DATE),
"%a %b %d %H:%M:%S %z %Y",tz = "GMT")+2*60*60

or I can also look at @skythelimit who’s usually twitting from Singapore (I am in Montréal). I can seen clearly when we might have overlaps,

res=extractR("skythelimit")

Nice isn’t it. But it is possible to do much better… for instance, for those who do not ask specifically not to be Geo-located, we can see where they do tweet during the day, and during the night… I am quite sure a dozen posts with those functions can be written…

La pyramide des âges, au Canada

L’autre jour, Statistique Canada publiait une jolie animation interactive permettant de visualiser la déformation de la “pyramide” des âges (évoquée par exemple surhttp://www.lapresse.ca/). Bon, j’ai des données assez proches, tirées dehttp://www.mortality.org/).

> pop=read.table(
+ "http://freakonometrics.blog.free.fr/PopulationCanada.txt",
+ header=TRUE,skip=2)
> pop$Age=as.numeric(as.character(pop$Age))
> pop$Year=as.numeric(as.character(pop$Year))
> pop=pop[is.na(pop$Age)==FALSE,]

On peut faire assez facilement des pyramides des ages, par année, avec le code suivant,

> library(plotrix)
> seuils=seq(0,110,by=10)
> pop$tranche=cut(pop$Age,seuils, right = FALSE)
> an=1955
> base=pop[pop$Year==an,]
> women=aggregate(base$Female,
+ list(base$tranche),sum)[,2]
> men=aggregate(base$Female,
+ list(base$tranche),sum)[,2]
> nom=as.character(unique(pop$tranche))
> pyramid.plot(men/sum(men)*100,
+ women/sum(women)*100,labels=nom,gap=2,
+  lxcol=c("blue","blue","purple","purple","purple",
+ "purple","red","red","red","red","red"),
+  rxcol=c("blue","blue","purple","purple","purple",
+ "purple","red","red","red","red","red"))

http://freakonometrics.blog.free.fr/public/perso6/pyramide-ages-canada2.gif

Au lieu de faire des barres horizontales, on peut cumuler, par tranche d’age, et regarder l’évolution dans le temps des proportions “moins de 20 ans” (bleu), “entre 20 et 60 ans” (mauve) et “plus de 60 ans” (rouge).

> seuils=c(0,20,60,110)
> pop$tranche2=cut(pop$Age,seuils, right = FALSE)
> YEAR=1921:2008
> totaux=matrix(NA,length(YEAR),3)
> for(i in 1:length(YEAR)){
+ base=pop[pop$Year==YEAR[i],]
+ totaux[i,]=aggregate(base$Total,
+ list(base$tranche2),sum)[,2]
+ }
> stotaux=apply(totaux,1,sum)
> X=totaux[,1]/stotaux
> Y=X+totaux[,2]/stotaux

Si on regarde bien la pyramide, on s’aperçoit qu’après guerre, on a une génération importante, qui se déplace ensuite progressivement vers le haut, les baby-boomers. Pour visualiser encore davantage cette génération, on peut l’isoler sur le graphique ci-dessous, en suivant la cohorte née entre 1945 et 1950.

Mais on reviendra un peu plus en détails très bientôt sur cette génération (peut-être plutôt sur des données françaises cette fois afin de limiter les erreurs d’interprétation). A suivre…

Claims reserving and IBNR with R

Following previous posts on life contingencies and longevity and mortality models, I upload additional material for the short course at the 6th R/Rmetrics Meielisalp Workshop & Summer School on Computational Finance and Financial Engineering organized by ETH Zürich, https://www.rmetrics.org/. The third part of the talk (on Actuarial models with R) will be dedicated to IBNR and claims reserving. A complete set of slides can be downloaded from the blog, but again, only some part will be presented. Note that the slides start with a parallel between mortality tables (in life insurance) and payment triangles (in non-life insurance).

Once again, the codes are from a book on actuarial science in R, written with Christophe Dutang (so far in French) that should appear, some day… The code used in the slides above are based on the following datasets,

> source("https://perso.univ-rennes1.fr/arthur.charpentier/ + bases.R")

We will built our own functions to derive all quantities. One function used can be found here

> source("https://perso.univ-rennes1.fr/arthur.charpentier/ + merz-wuthrich-triangle.R")

Finally, note that most of the code can be found in the following library

> library(ChainLadder)

Longevity and mortality dynamics with R

Following the previous post on life contingencies and actuarial models in life insurance, I upload additional material for the short course at the 6th R/Rmetrics Meielisalp Workshop & Summer School on Computational Finance and Financial Engineering organized by ETH Zürich, https://www.rmetrics.org/. The second part of the talk (on Actuarial models with R) will be dedicated to longevity and mortality. A complete set of slides can be downloaded from the blog, but again, only some part will be presented.

As mentioned earlier, the codes are from a book on actuarial science in R, written with Christophe Dutang (so far in French) that should appear, some day… The code used in the slides above can be downloaded from here, and datasets are the following,

> DEATH <- read.table(
+ "http://freakonometrics.free.fr/Deces-France.txt",
+ header=TRUE)
> EXPO  <- read.table(
+ "http://freakonometrics.free.fr/Exposures-France.txt",
+ header=TRUE,skip=2)

For additional resources, I will use Rob Hyndman‘s package on demography, Heather Turner and David Firth’s package on generalized nonlinear models (e.g. the slides of the short course Heather gave in Rennes at the UseR! conference in 2009), as well as functions developed by JPMorgan’s LifeMetrics (functions are  fully documented in the LifeMetrics Technical Document). All those functions can be obtained using

> library(demography)
> library(gnm)
> source("http://freakonometrics.free.fr/fitModels.R")

Life contingencies with R

I will be giving in less than four weeks a short course at the 6th R/Rmetrics Meielisalp Workshop & Summer School on Computational Finance and Financial Engineering organized by ETH Zürich, https://www.rmetrics.org/. The talk will be on Actuarial models with R, and first part will be dedicated to life insurance. A complete set of slides can be downloaded from the blog, but in the talk, only some part will be presented.

The codes are from a book on actuarial science in R, written with Christophe Dutang (so far in French) that should appear, some day… The code used in the slides can be downloaded from here, and datasets are the following,

> TD <- read.table(
+ "https://perso.univ-rennes1.fr/arthur.charpentier/TD8890.csv",sep=";",header=TRUE)
> TV <- read.table(
+ "https://perso.univ-rennes1.fr/arthur.charpentier/TV8890.csv",sep=";",header=TRUE)

For additional resources, I recommend Emiliano’s website, http://www.math.uconn.edu/, with great lectures on life insurance mathematics, and the (new) lifecontinfencies vignette on http://cran.r-project.org/,

> library(lifecontingencies)

French dataset: population and GPS coordinates

A short post today based on recent work by @3wen (Ewen Gallic, graduate Student in Rennes, spending a year in Montreal). Since we were working on a detailed French dataset (per commune), we needed a dataset containing a list allcommunes, with population and location. GPS coordinates were extracted from Google, using the following php file, inspired by http://www.andrew-kirkpatrick.com/ on Google geocoding api with php webpage. Population was interpolated from INSEE’s datasets, i.e. http://www.insee.fr/ (since data are over a 35 year period, from 1975 to 2010, changes have been taken into account as carefully are possible – e.g. merges and splits of cities – based on thatdescription). A spline model has been used for all cities (with three degrees of freedom, and null and negative interpolation became one, since we’ll be using loglinear models afterwards). Names are from that dataset, still on INSEE’s website, http://www.insee.fr/.

A zipped file can be downloaded here popfr19752010.zip, but it is also possible to use the code below (it is a 24Mo dataset). Since it was hard to find such a dataset online (different files can be found, but we found none with population and location), we have decided to upload that dataset. Please let us know if there are problems with those data…

> base=read.csv(
+ "http://freakonometrics.free.fr/popfr19752010.csv",
+ header=TRUE)

Using that code, it is possible to locate all the communes in France (metropolitan), for instance

> library(maps)
> map("france")
> points(base$long,base$lat,cex=.1,col="red",pch=19)
> points(base$long,base$lat,cex=2*base$pop_2010/
+ max(base$pop_2010),col="blue",pch=19)

Several additional lines of code on that dataset (and also others) will be uploaded, soon.

Cette oeuvre est mise à disposition sous licence Paternité – Partage à l’Identique 3.0 non transposé. Pour voir une copie de cette licence, visitez http://creativecommons.org/. Date : 24 mai 2012, par Ewen GALLIC. Sources : INSEE, API Google Maps v3 et GeoHack (coordonnées GPS), propres calculs (estimation de population à partir des données INSEE).

  • reg : code region INSEE (character)
  • dep : code departement INSEE (character, corse 201 et 202 au lieu de 2A et 2B)
  • com : code commune INSEE (character)
  • article : article du nom de la commune (character)
  • com_nom : nom de la commune (character)
  • long : longitude (numeric)
  • lat : latitude (numeric)
  • pop_i : estimation de la population à la date i (ramenée à 1 si <=0), i=1975,…,2010 (numeric)

Qui ne voit pas la présence du diable dans la loi ?

Je suis de plus en plus dérouté par le mélange des genres dans l’espace médiatique, de ces prétendus spécialistes qui nous expliquent des choses qu’ils ont du découvrir 10 minutes avant de commencer à taper un article. Ce matin, des éditorialistes se sentaient autoriser à fabuler sur un sondage conduit sur internet auprès de 800 personnes, et j’avoue que cela m’a agacé… J’ignore si j’ai davantage de légitimité, mais je me suis dit que la loi 78 adoptée vendredi (au Québec) pourrait être un prétexte intéressant pour écrire sur un sujet d’actualité (chose que j’hésite d’ordinaire à faire, car ce n’est pas mon métier, loin de là). Pourquoi ne pourrais-je pas, moi aussi, faire mon malin et analyser un texte de loi ? On va donc regarder cette loi dans le détail, à partir du texte en ligne sur http://profscontrelahausse.org/ (ou pour être complètement honnête, le début seulement de la loi)

L=paste("Dans la presente loi, a moins",
"que le contexte nindique",
"un sens different, on entend par:",
"association d’etudiants}}:une",
"association ou un regroupement ",
"d’associations de niveau postsecondaire",

etc

"la qualite de l’enseignement, les ",
"services requis de façon a ",
"tenir compte des circonstances particulieres",
"resultant de l’interruption ",
"de la session d’hiver de l’annee 2012 ou ",
"de la session d’ete de l’annee 2012.")

L2=substring(toupper(L),1:nchar(L),1:nchar(L))
alphabet=c("A","B","C","D","E","F","G"
,"H","I","J","K","L","M","N","O",
"P","Q","R","S","T","U","V","W","X","Y","Z")
L2=L2[L2%in%alphabet]
ML=t(matrix(L2,nc,nl))

On peut visualiser le début du texte ci-dessous (l’idée est de ne garder que les lettres, pour faire simple, mais on peut tenter plus largement en gardant la ponctuation, et les chiffres, ça ne change rien aux algorithmes que nous allons développer)

Commencons pas des trivialités: le texte de cette loi est fondamentalement inégalitaire: ce sont toujours les mêmes qui sortent gagnants, et les mêmes qui sortent perdants,

> sort(table(L))[c(1,2,3,4,5,20,21,22,23)]
L
Y   J   X   H   F   N   A   S   E
7   8  10  24  26 355 386 402 820
On retrouve ici reproduites les inégalités observées dans presque tous les textes (sauf peut être La Disparition, mais c’est une autre histoire): N, A, S et E en tête, et toujours les mêmes petites lettres qui se trouvent honteusement oubliées…

Maintenant, si on prend le temps de regarder en détail le contenu de la loi, force est de constater qu’on y voit la présence du Malin constamment ! Pour reprendre une idée de Bahye ben Asher ibn Halawa (רבינו בחיי) on peut tenter d’utiliser le code de la ThorahJulien Prévieux avait utilisé cette technique il y a quelques temps pour retrouver toute une terminologie associée aux crashs boursiers dans une page du Capital de Karl Marx, prise au hasard, ou aux paradis fiscaux dans la Richesse des Nations d’Adam Smith. Les amateurs de la théorie de Ramsey prétendront que c’est une supercherie, et ils auront probablement raison. Mais après avoir joué l’autre jour sur la recherche de schémas graphiques dans des images (en l’occurrence des rayures rouges et blanches, afin de repérer Charlie), j’ai eu envie de jouer à chercher des mots dans une matrice de lettres.

Commençons par nous donner un mot, au hasard, e.g.`
MOT="DEMON"
Le but est de voir si ce mot figure, ou pas, dans la matrice de lettre, suivant un schéma simple (c’est l’idée du code de la Thorah). Ma stratégie est ici assez simple: il faut commencer par chercher la lettre la moins présente dans la matrice, qui figure dans le mot. En fait, dans DEMON, je vais me focaliser sur les lettres du centre (en excluant la première et la dernière), i.e. DEMON
lettrelaplusrare=function(mot=MOT){
lettres=substring(mot,2:(nchar(mot)-1),
2:(nchar(mot)-1));
fL=TL[lettres];
return(names(fL)[which.min(fL)])}
rare=lettrelaplusrare()
On trouve ainsi la lettre qui va servir de base dans notre algorithme. A la rigueur, je peux faire une petite fonction qui va colorer une lettre dans la matrice,
dessinlettre=function(ltr="A",clr="yellow"){
which(L==ltr);
CLR=rep(NA,length(ML));
CLR[which(as.vector(ML)==ltr)]=clr;
dessin(ML,CLR,plt=FALSE)}
mais pour le moment, ça ne sert à rien. Dans un second temps, on regarde les lettres qui entourent la lettre la plus rare (c’est pour ça que j’ai exclus les lettres au bord, afin d’être certain d’avoir toujours une lettre avant, et une après),
lettres=substring(MOT,1:(nchar(mot)),
1:(nchar(mot)))
avant=lettres[which(lettres==rare)-1]
apres=lettres[which(lettres==rare)+1]
avant=avant[1]
apres=apres[1]
dessinlettre(rare,"yellow")
dessinlettre(avant,"green")
dessinlettre(apres,"pink")
On a ainsi trois lettres importantes: une relativement rare au centre, en jaune, une avant, en vert, et une après en rose. C’est ce qu’on peut visualiser sur le dessin ci-dessous,
Bon, ensuite on rentre dans l’algorithme, un peu lourd à mon gout: on va étudier tous les triplets possibles de ces trois lettres, en identifiants ceux pour lesquels la lettre au centre est précisément au milieu, coincée entre les deux autres lettres. C’est le principe de ce code de la Thorah. La première étape est de récupérer les coordonnées (dans la matrice) de ces trois lettres,
lrr=which(as.vector(ML)==rare)
lav=which(as.vector(ML)==avant)
lap=which(as.vector(ML)==apres) 
LIGrr=(lrr-1)%%nrow(ML)+1
COLrr=(lrr-1)%/%nrow(ML)+1
LIGav=(lav-1)%%nrow(ML)+1
COLav=(lav-1)%/%nrow(ML)+1
LIGap=(lap-1)%%nrow(ML)+1
COLap=(lap-1)%/%nrow(ML)+1
ptrr=cbind(LIGrr,COLrr)
ptav=cbind(LIGav,COLav)
ptap=cbind(LIGap,COLap)
Ensuite, on balaye l’ensemble des triplets (et ça peut être long),
pointsalignes=function(x,y,z){
((y[1]+2*(x[1]-y[1]))==z[1])&((y[2]+
2*(x[2]-y[2]))==z[2])}
LEQUEL=rep(NA,3)
for(a in 1:length(LIGrr)){
for(b in 1:length(LIGav)){
for(c in 1:length(LIGap)){
if(pointsalignes(ptrr[a,],ptav[b,],ptap[c,])==
TRUE){LEQUEL=cbind(LEQUEL,c(a,b,c))}
}}}
Parmi les triplets, on retient seulement ceux pour lesquels la lettre du milieu est précisément au centre, entre les deux autres. Enfin arrive la dernière étape pénible: regarder suivant la suite logique ce que sont les mots que l’on obtient (une fois placées ces trois lettres, la récurrence est complètement décrite). On retient alors les quintuplets correspondant au mot DEMON,
V=LEQUEL
for(k in 2:ncol(V)){
CLR=rep(NA,length(ML));
i1=ptrr[V[1,k],];
CLR[(i1[2]-1)*ncol(ML)+i1[1]]=clr;
lettres=substring(MOT,1:(nchar(mot)),
1:(nchar(mot)))
milieu=which(lettres==rare)
obtenu=ML[i1[1],i1[2]]
for(u in 1:(milieu-1)){
i2=ptav[V[2,k],]+(u-1)*(ptav[V[2,k],]-ptrr[V[1,k],]);
indice=(i2[2]-1)*ncol(ML)+i2[1]
if((i2[2]<=ncol(ML))&(i2[1]>=1)&(i2[1]<=
nrow(ML))&(i2[1]>=1)){CLR[indice]=clr;
obtenu=paste(ML[i2[1],i2[2]],obtenu,sep="")}}
for(u in (milieu+1):nchar(MOT)){
i3=ptap[V[3,k],]+(u-milieu-1)*(ptap[V[3,k],]-
ptrr[V[1,k],]);
indice=(i3[2]-1)*ncol(ML)+i3[1]
if((i3[2]<=ncol(ML))&(i3[1]>=1)&(i3[1]<=
nrow(ML))&(i3[1]>=1)){CLR[indice]=clr;
obtenu=paste(obtenu,ML[i3[1],i3[2]],sep="")}}
if(obtenu==MOT)
{print(paste(k,"  ",obtenu))}}
Reste à faire tourner cette boucle. Et ca ne manque pas, DEMON apparaît effectivement dans le texte,
On peut aller vérifier, pour ceux qui douteraient,

Coïncidence ? Peut-être…. alors tentons un autre mot, comme MALIN, qui apparait aussi,

ou alors SATAN… ce dernier apparaît huit fois dans le début du texte !

Étonnant non…. en tous les cas, je me suis bien amusé à coder ce petit algorithme, et de voir qu’il fonctionnait aussi bien ! Les amateurs de probabilités pourront se lancer dans des calculs, ou consulter quelques articles qui évoquent cette technique, comme McKay,Bar-Natan, Bar-Hillel & Kalai (2001), et tous les articles évoqués dans les références… Quant au juge qui devra se prononcer sur ma condamnation pour blasphème et pour hérésie, il va de soi que ce billet est une blague. Le diable n’existe pas ! On n’est plus au moyen-age quand même….

Births and week-ends, in France

This week, I have seen on the internet (sorry, I cannot find proper references) the graph produced here on the right: which birthday is most likely ? The fact that I have no further information is important, since I do not know in which country such a graph was obtained. At least, I know it should not be France…

In France, I have already mentioned that there is a strong week-end effect: nowadays, there is 25% less deliveries during week-ends than during the week. Calot (1981) observed already that there were less deliveries on Sundays. This has been confirmed more recently, e.g. in http://www.lepoint.fr/ or http://www.prepabl.fr/, with a significant difference between week days, and week-ends. Here  is the number of birth per day, over 40 years, with in blue the average trend during the week, and in red, during week-ends,

naissance=read.table(
"http://freakonometrics.free.fr/naissanceFR2.txt")
attach(naissance)
date=as.Date(date)
plot(date, nbre,cex=.5)
t2=as.POSIXlt(date)
jour=t2$wday
X=naissance$date
Y=naissance$nbre
J=jour
df=data.frame(X,Y,J)
library(splines)
regs=lm(Y~bs(X,df=20),data=df[jour%in%c(0,6),])
Yp=predict(regs,newdata=df)
lines(X,Yp,col="red",lwd=3)
regs=lm(Y~bs(X,df=20),data=df[jour%in%1:5,])
Yp=predict(regs,newdata=df)
lines(X,Yp,col="blue",lwd=3)

If we look at the evolution of the ratio week-ends over weeks days, we have the following graph

t2=as.POSIXlt(date)
jour=t2$wday
jour=jour[1:(1982*7)]
nbre2=jour
for(i in 1:1982){
taux=sum(nbre[6:7+7*(i-1)])/
sum(nbre[1:5+7*(i-1)])/2*5
nbre2[1:5+7*(i-1)]=nbre[1:5+7*(i-1)]*taux
nbre2[6:7+7*(i-1)]=nbre[6:7+7*(i-1)]
nbre2[1:7+7*(i-1)]=
mean(nbre[1:7+7*(i-1)])/mean(nbre2[1:7+7*(i-1)])*
nbre2[1:7+7*(i-1)]
}
nbretaux=jour
for(i in 1:1982){
taux=sum(nbre[6:7+7*(i-1)])/
sum(nbre[1:5+7*(i-1)])/2*5
nbretaux[1:7+7*(i-1)]=taux
}
plot(date[1:length(nbre2)],nbretaux)
X= date[1:length(nbre2)]
Y=nbretaux
library(splines)
reg=lm(Y~bs(X,df=20))
Yp=predict(reg)
lines(X,Yp,col="red",lwd=3)

In the beginning of the 70’s, during week-ends, there were 5% less deliveries, but 25% less around 2000. It is then possible to produce the same kind of graphs as the one above, per year of birth. And here, we clearly observe the importance of the week end effect (maybe also because of color choice)

naissance=read.csv(
"http://freakonometrics.free.fr/naissanceFR.csv",
sep=";")
M=as.matrix(naissance[,3:ncol(naissance)])
BIRTH=as.vector(t(M))
YEAR=rep(1968:2005,each=12*31)
MONTH=rep(rep(1:12,each=31),38)
DAY=rep(1:31,12*38)
X=NA
for(y in 1968:2005){
sbase=base[YEAR==y,]
X=c(X,sbase$BIRTH/sum(sbase$BIRTH,
na.rm=TRUE))
}
base=data.frame(YEAR,MONTH,DAY,
BIRTH,BIRTHDAYPROB=X[-1])

m1=min(base$BIRTHDAYPROB,na.rm=TRUE)
m2=max(base$BIRTHDAYPROB,na.rm=TRUE)
y=1980
colr=rev(heat.colors(100))
sbase=base[YEAR==y,]
plot(0:1,0:1,col="white",xlim=c(-1,12),
ylim=c(-31,1),axes=FALSE,xlab=
paste("Naissance en",y,sep=" "),ylab="")
for(x in 1:nrow(sbase)){
a=sbase$MONTH[x];b=sbase$DAY[x]
polygon(c(a-.9,a-.9,a-.1,a-.1),-c(b-.9,b-.1,
b-.1,b-.9),col=colr[(sbase$BIRTHDAYPROB[x]-m1)/
(m2-m1)*100],border=NA)
}
text((1:12)-.5,.5,c("J","F","M","A","M","J","J",
"A","S","O","N","D"),cex=.7)
text(-.5,-(1:31)+.5,1:31,cex=.7)

"sendo l'intento mio scrivere cosa utile a chi la intende…"