I don’t know if it’s because of Halloween (or the new season of walking dead), but I had weird thoughts recently. For instance, I wondered if zombies walk was a random walk… Note that I am not the only one: there are papers on that topic on Wired or even Nature (at least on one of their blogs). Before modeling that walk, it might be better to describe the process, and to generate zombies walk through simulations.
Monthly Archives: October 2011
Régression sur des variables catégorielles
Petit complément par rapport au cours de mardi. On avait évoqué tout d’abord la lecture des sorties lorsque l’on régresse sur des variables catégorielles (des facteurs). Commençons par supprimer la constante de la régression
> reg0=glm(nbre~0+zone,offset=log(exposition),data=base, + family=poisson(link="log")) > summary(reg0) Call: glm(formula = nbre ~ 0 + zone, family = poisson(link = "log"), data = base, offset = log(exposition)) Deviance Residuals: Min 1Q Median 3Q Max -0.5717 -0.3968 -0.2996 -0.1547 12.6722 Coefficients: Estimate Std. Error z value Pr(>|z|) zoneB -2.54187 0.06287 -40.43 <2e-16 *** zoneA -2.54912 0.05285 -48.23 <2e-16 *** zoneC -2.38525 0.03753 -63.56 <2e-16 *** zoneD -2.13454 0.03878 -55.05 <2e-16 *** zoneE -2.00204 0.03965 -50.49 <2e-16 *** zoneF -2.06932 0.11547 -17.92 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 50966 on 50000 degrees of freedom Residual deviance: 15692 on 49994 degrees of freedom AIC: 20800 Number of Fisher Scoring iterations: 6 > predict(reg0,newdata=data.frame( + zone=c("A","B","C","D","E"),exposition=rep(1,5))) 1 2 3 4 5 -2.549120 -2.541870 -2.385253 -2.134543 -2.002044
On voit que toutes les modalités sont présentes, et toutes sont significatives. Si on régresse sur la constante, il faudra supprimer une modalité pour rendre le modèle identifiable. On peut forcer pour que la modalité de référence soit la seconde,
> base$zone=relevel(base$zone,"B") > regB=glm(nbre~zone,offset=log(exposition),data=base, + family=poisson(link="log")) > summary(regB) Call: glm(formula = nbre ~ zone, family = poisson(link = "log"), data = base, offset = log(exposition)) Deviance Residuals: Min 1Q Median 3Q Max -0.5717 -0.3968 -0.2996 -0.1547 12.6722 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.54187 0.06287 -40.431 < 2e-16 *** zoneA -0.00725 0.08213 -0.088 0.929661 zoneC 0.15662 0.07322 2.139 0.032432 * zoneD 0.40733 0.07387 5.514 3.50e-08 *** zoneE 0.53983 0.07433 7.263 3.80e-13 *** zoneF 0.47255 0.13148 3.594 0.000325 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 15809 on 49999 degrees of freedom Residual deviance: 15692 on 49994 degrees of freedom AIC: 20800 Number of Fisher Scoring iterations: 6 > predict(regB,newdata=data.frame( + zone=c("A","B","C","D","E"),exposition=rep(1,5))) 1 2 3 4 5 -2.549120 -2.541870 -2.385253 -2.134543 -2.002044
On notera que les prédictions ne changent pas. On peut aussi choisir la première comme modalité de référence,
> base$zone=relevel(base$zone,"A") > reg=glm(nbre~zone,offset=log(exposition), > data=base,family=poisson(link="log")) > summary(reg) Call: glm(formula = nbre ~ zone, family = poisson(link = "log"), data = base, offset = log(exposition)) Deviance Residuals: Min 1Q Median 3Q Max -0.5717 -0.3968 -0.2996 -0.1547 12.6722 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.54912 0.05285 -48.232 < 2e-16 *** zoneB 0.00725 0.08213 0.088 0.929661 zoneC 0.16387 0.06482 2.528 0.011471 * zoneD 0.41458 0.06555 6.324 2.54e-10 *** zoneE 0.54708 0.06607 8.280 < 2e-16 *** zoneF 0.47980 0.12699 3.778 0.000158 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 15809 on 49999 degrees of freedom Residual deviance: 15692 on 49994 degrees of freedom AIC: 20800 Number of Fisher Scoring iterations: 6
Le fait que la seconde modalité ne soit pas significative se lit par rapport à la modalité de référence (en l’occurrence la première): non significatif signifie alors non significativement différente. Autrement dit, on peut regrouper les modalités en une seule.
> base$zonesimple=base$zone > base$zonesimple[base$zone%in%c("A","B")]="A" > reg=glm(nbre~zonesimple,offset=log(exposition), + data=base,family=poisson(link="log")) > summary(reg) Call: glm(formula = nbre ~ zonesimple, family = poisson(link = "log"), data = base, offset = log(exposition)) Deviance Residuals: Min 1Q Median 3Q Max -0.5717 -0.3959 -0.2989 -0.1547 12.6722 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.54612 0.04046 -62.937 < 2e-16 *** zonesimpleC 0.16087 0.05518 2.915 0.00355 ** zonesimpleD 0.41158 0.05604 7.345 2.06e-13 *** zonesimpleE 0.54408 0.05665 9.605 < 2e-16 *** zonesimpleF 0.47681 0.12235 3.897 9.74e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 15809 on 49999 degrees of freedom Residual deviance: 15692 on 49995 degrees of freedom AIC: 20798 Number of Fisher Scoring iterations: 6
On note qu’avec ce regroupement, les autres modalités sont sensiblement différentes. On peut aussi retenir la troisième comme modalité de référence
> base$zonesimple=relevel(base$zonesimple,"C") > reg=glm(nbre~zonesimple,offset=log(exposition), + data=base,family=poisson(link="log")) > summary(reg) Call: glm(formula = nbre ~ zonesimple, family = poisson(link = "log"), data = base, offset = log(exposition)) Deviance Residuals: Min 1Q Median 3Q Max -0.5717 -0.3959 -0.2989 -0.1547 12.6722 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.38525 0.03753 -63.557 < 2e-16 *** zonesimpleA -0.16087 0.05518 -2.915 0.00355 ** zonesimpleD 0.25071 0.05396 4.646 3.39e-06 *** zonesimpleE 0.38321 0.05460 7.019 2.24e-12 *** zonesimpleF 0.31593 0.12142 2.602 0.00927 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 15809 on 49999 degrees of freedom Residual deviance: 15692 on 49995 degrees of freedom AIC: 20798 Number of Fisher Scoring iterations: 6
Comme toutes les modalités semblent significatives, on peut tenter de prendre comme modalité de référence une des dernières (dont les estimations des coefficients donnent des résultats très proches)
> base$zonesimple=relevel(base$zonesimple,"F") > reg=glm(nbre~zonesimple,offset=log(exposition), + data=base,family=poisson(link="log")) > summary(reg) Call: glm(formula = nbre ~ zonesimple, family = poisson(link = "log"), data = base, offset = log(exposition)) Deviance Residuals: Min 1Q Median 3Q Max -0.5717 -0.3959 -0.2989 -0.1547 12.6722 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.06932 0.11547 -17.921 < 2e-16 *** zonesimpleC -0.31593 0.12142 -2.602 0.00927 ** zonesimpleA -0.47681 0.12235 -3.897 9.74e-05 *** zonesimpleD -0.06522 0.12181 -0.535 0.59232 zonesimpleE 0.06727 0.12209 0.551 0.58161 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 15809 on 49999 degrees of freedom Residual deviance: 15692 on 49995 degrees of freedom AIC: 20798 Number of Fisher Scoring iterations: 6
Au vue de cette dernière sortie, on peut tenter de fusionner toutes les dernières classes ensembles
> base$zonesimple[base$zone%in%c("D","E","F")]="F" > reg=glm(nbre~zonesimple,offset=log(exposition), + data=base,family=poisson(link="log")) > summary(reg) Call: glm(formula = nbre ~ zonesimple, family = poisson(link = "log"), data = base, offset = log(exposition)) Deviance Residuals: Min 1Q Median 3Q Max -0.5660 -0.3959 -0.3004 -0.1547 12.5929 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.07182 0.02696 -76.853 < 2e-16 *** zonesimpleC -0.31344 0.04621 -6.783 1.18e-11 *** zonesimpleA -0.47431 0.04861 -9.757 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 15809 on 49999 degrees of freedom Residual deviance: 15698 on 49997 degrees of freedom AIC: 20800 Number of Fisher Scoring iterations: 6
Bon, formellement, regrouper deux modalités (i.e. décréter que deux variables sont non significatives simultanément) demande un peu plus qu’un test de Student, ou que deux tests de Student…. Si on remonte un peu en arrière, on aurait pu faire un test multiple avant de fusionner les trois modalités (un test de type Fisher, ou une ANOVA)
> base$zonesimple=relevel(base$zonesimple,"F") > reg=glm(nbre~zonesimple,offset=log(exposition), + data=base,family=poisson(link="log")) > library(car) > linearHypothesis(reg,c("zonesimpleD=0","zonesimpleE=0")) Linear hypothesis test Hypothesis: zonesimpleD = 0 zonesimpleE = 0 Model 1: restricted model Model 2: nbre ~ zonesimple Res.Df Df Chisq Pr(>Chisq) 1 49997 2 49995 2 5.7073 0.05763 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Manifestement, on peut accepter l’hypothèse que ces trois catégories n’en font qu’une. La zone géographique peut alors se découper en trois grandes zones, et pas en six. On notera que cela correspond à ce que propose un arbre de régression
Consecutive number and lottery
Recently, I have been reading odd things about strategies to win at the lottery. E.g.
or
I wrote something a long time ago, but maybe it would be better to write another post. First, it is easy to get data on the French lotteries, including draws, number of winners and gains,
loto=read.table("http://freakonometrics.blog.free.fr/public/ data/loto.csv",sep=";",header=TRUE) balls=loto[,c("boule_1","boule_2","boule_3", "boule_4","boule_5","boule_6")] q=function(x){quantile(x,(0:5)/5)} sortballs=balls consec=balls[,-1] sortconsec=consec for(i in 1:nrow(balls)){sortballs[i,]=q(balls[i,]) consec[i,]=sortballs[i,2:6]-sortballs[i,1:5] sortconsec[i,]=sort(consec[i,])} winner1=loto[,"nombre_de_gagnant_au_rang1"] gain1=as.numeric(as.character(loto[,"rapport_du_rang1"])) winner2=loto[,"nombre_de_gagnant_au_rang2"] gain2=as.numeric(as.character(loto[,"rapport_du_rang2"])) winner3=loto[,"nombre_de_gagnant_au_rang3"] gain3=as.numeric(as.character(loto[,"rapport_du_rang3"])) winner4=loto[,"nombre_de_gagnant_au_rang4"] gain4=as.numeric(as.character(loto[,"rapport_du_rang4"])) winner5=loto[,"nombre_de_gagnant_au_rang5"] gain5=as.numeric(as.character(loto[,"rapport_du_rang5"])) which1=(sortconsec[,1]==1) which2=(sortconsec[,2]==1) which3=(sortconsec[,3]==1) which4=(sortconsec[,4]==1) which5=(sortconsec[,5]==1)
There several ways to defining “winning at the lottery” (2 out of 6, 3 out of 6, 4 out of 6, etc) and to define “having consecutive numbers” (it can be 2 out of 6, 3 out of 6, etc). For instance,
It is also possible to compare the number of winners obtained for medium winners (3 out of 6, so called vainqueur de rang 4) when there were 2 consecutive numbers
> t.test(winner4[which1==TRUE],winner4[which1==FALSE]) Welch Two Sample t-test data: winner4[which1 == TRUE] and winner4[which1 == FALSE] t = -3.2132, df = 4792.491, p-value = 0.001321 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -6864.430 -1662.123 sample estimates: mean of x mean of y 33887.82 38151.10
With a simple mean comparison test, we have that there is a significant difference between average number of winners when there were at least 2 consecutive numbers out of 6 balls drawn. And actually, the average number of winners was lower when there are consecutive numbers. And if we look at the average gain, we have also a significant difference
> t.test(gain4[which1==TRUE],gain4[which1==FALSE]) Welch Two Sample t-test data: gain4[which1 == TRUE] and gain4[which1 == FALSE] t = 5.8926, df = 3675.361, p-value = 4.143e-09 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 11.06189 22.09337 sample estimates: mean of x mean of y 173.9788 157.4012
Here we see that if we play consecutive numbers, on average, the gain is larger. Perhaps it would be better to look at that on graphs.
which0=which1 WIN=c(mean(winner5[which0==TRUE]),mean(winner5[which0==FALSE]), mean(winner4[which0==TRUE]),mean(winner4[which0==FALSE]), mean(winner3[which0==TRUE]),mean(winner3[which0==FALSE]), mean(winner2[which0==TRUE]),mean(winner2[which0==FALSE]), mean(winner1[which0==TRUE]),mean(winner1[which0==FALSE])) MWIN=matrix(WIN,2,5) plot(1:5,MWIN[1,],type="b",col="red",log="y", ylim=c(1,1000000),xlab="two consecutive numbers", ylab="number of winners (log scale)") lines(1:5,MWIN[2,],type="b",col="blue",pch=4)
If we focus on the case where “having consecutive numbers” means two consecutive numbers, we have below the number of winners, with first rank (6 out of 6), then second rank (5 out of 6), etc,
Note that the y-axis is on a log scale, and that draws with consecutive balls are in red, and no consecutive balls are in blue. If we focus on average gains, curves are in opposite order,
But if we consider the case of three consecutive balls, we have, for the number of winners,
or for average gains
Here it starts to get slightly different: there are more “big winners” when there are at least three consecutive numbers. And with four consecutive numbers, it is clearly the opposite
Here we see that there are much more winners with four consecutive numbers (actually, it might be a triplet and a pair). So I have to confess that I am not convinced by the conclusion: actually a lot a people pick consecutive numbers… Actually, if we look at draws where there were the more winners, we can clearly see that a lot of players like consecutive numbers (perhaps not has much as playing birthdays, since most numbers are lower than 31),
loto[loto$"nombre_de_gagnant_au_rang1">50, c("combinaison_gagnante_en_ordre_croissant", "nombre_de_gagnant_au_rang1","date_de_tirage")] combinaison_gagnante gagnant_au_rang1 3189 2-4-13-16-28-31 103 3475 1-5-9-10-12-25 64 4018 4-5-7-14-15-17 63 4396 26-27-28-35-36-37 64 4477 7-11-15-27-33-44 53 4546 2-9-12-14-19-24 60 4685 2-8-10-12-14-16 96 date_de_tirage 3189 19930626 3475 19920212 4018 19880504 4396 19840919 4477 19830914 4546 19820519 4685 19790919
On September 1979, there were 5 even consecutive numbers (ok, it can not be considered strictly as consecutive numbers) and 96 winners with 6 numbers out of 6 ! And if we look at the others, even if they are not strictly consecutive, there is a lot of regularity. So I believe that picking consecutive numbers might not be a great strategy if you want to win a lot of money !
On tourne en rond
Hier soir, nous sommes allé mangé au restaurant, en famille, tous les 5. Le serveur nous a placé à une table ronde. Tout le monde s’est installé. Alors que nous commencions l’apéro, ma fille s’est exclamé “ça alors, c’est dingue, on est rangé par age…“. Vraiment ? c’est une coïncidence exceptionnelle ?… allons, quelle est la probabilité que si 5 personnes se mettre autour d’une table ronde, les personnes soient ordonnées par age ? Je vais laisser un peu de temps pour trouver la solution…
Personnellement je dirais 1 chance sur 12….
oui, car 12, c’est 2 fois 3 fois 2..
Non ? bon, la réponse en images alors…
Bon, et je remets en jeu le titre de Sheldon d’or au premier qui reconnaît le morceau…
Dictées, et mathématiques
Tous les soirs, c’est la guerre avec mon fils pour faire les devoirs… et surtout pour préparer la dictée du vendredi. Car pour se préparer, on fait trois dictées à la maison, avant la dictée en classe (le mardi soir, le mercredi soir et le jeudi soir). Le truc, c’est qu’avec sa mère, on n’adore pas faire la dictée, donc la règle tacite est “celui qui a fait la dictée un soir ne la fera pas le lendemain…“.
Mon fils, qui cherchait à négocier un cinéma en fin de semaine, proposa la règle suivante “si je fais moins de trois fautes au moins une fois je peux aller au cinéma samedi !…“. Afin de ne pas sembler trop laxiste, il a fallu que je demande un peu plus, “euh… si deux soirs de suite, tu fais moins de trois fautes, tu as gagné !…“. Bon, il faut reconnaître que j’ai toujours tendance à faire des dictées avec des phrases plus alambiquées, et la probabilité de faire plus de trois fautes est plus grande s’il fait une dictée avec son père que s’il fait une dictée avec sa mère. Avant-hier, comme nous étions mardi soir, je lui demande “maintenant que tu connais la règle, avec qui tu commences ?“. Alors, quelle aurait du être la stratégie optimale de mon fils…?
Oui, il faut que tout le monde cherche, alors on va poursuivre encore un peu l’interlude musical (d’ailleurs, mon admiration sans bornes au premier qui reconnaît les deux morceaux: je lui remettrais personnellement un Sheldon d’or)…
Faisons le calcul: soit 
la probabilité de faire moins de trois fautes avec son père, et 
la probabilité de faire moins de trois fautes avec sa mère. S’il commence avec son père (i.e. puis sa mère puis enfin son père, sur les trois soirs), la probabilité de gagner son cinéma est
En revanche, en commençant avec sa mère, elle s’écrit
Or comme 
(la probabilité de faire moins de trois fautes est plus faible avec son père), alors 
, et donc, la seconde stratégie est moins bonne: il a donc intérêt à faire deux fois la dictée avec moi, même si c’est avec moi qu’il fait le plus de fautes… Étonnant non ? (en fait pas tant que ça si on note que pour gagner, il est indispensable qu’il ait moins de trois fautes à la dictée du mercredi… donc il est naturel de penser qu’il faut que cette dictée se fasse avec sa mère)
Short selling, volatility and bubbles
Yesterday, I wrote a post (in French) about short-selling in financial market since some journalists claimed that it was well-known that short -selling does increase volatility on financial market. Not only in French speaking journals actually, since we can read on http://www.forbes.com that « in a market with restrictions on short-selling, volatility is reduced ». But things are not that simple. For instancehttp://www.optionsatoz.com/ explains it from a theoretical point of view. But we can also look at the data. For instance, we can compare the stock price of Air China, exchanged in Shanghai in blue (where short-selling is forbidden) and in Hong Kong in rouge (where short-selling is allowed), since @Igor gave me the tickers of those stocks
library(tseries) X<-get.hist.quote("0753.HK") Y<-get.hist.quote("601111.SS") plot(Y[,4],col="blue",ylim=c(0,30)) lines(X[,4],col="red")
But as @alea_ pointed out, one asset is expressed here in Yuans renminbi, and the other one in HK dollars. So I downloaded the exchange rate fromhttp://www.oanda.com/
Z=read.table("http://freakonometrics.blog.free.fr/public/ data/change-cny-hkd.csv",header=TRUE,sep=";",dec=",") D=as.Date(as.character(Z$date),"%d/%m/%y") z=as.numeric(Z$CNY.HKD) plot(D,z,type="l") X2=X[,4] for(t in 1:length(X2)){ X2[t]=X2[t]*z[D==time(X2[t])]} X2=X[,4] plot(Y[,4],col="blue",ylim=c(0,30)) lines(X2,col="red")
Now both stocks are expressed in the same currency. To compare returns volatility, a first idea can be to use GARCH models,
RX=diff(log(X2)) RY=diff(log(Y[,4])) Xgarch = garch(as.numeric(RX)) SIGMAX=predict(Xgarch) Ygarch = garch(as.numeric(RY)) SIGMAY=predict(Ygarch) plot(time(Y)[-1],SIGMAY[,1],col="blue",type="l") lines(time(X2)[-1],SIGMAX[,1],col="red")
But volatility is here too eratic. So an alternative can be to use exponentially-weighted moving averages, where simple recursive relationships are considered
or equivalently
The code is not great, but it is easy to understand,
moy.ew=function(x,r){ m=rep(NA,length(x)) for(i in 1:length(x)){ m[i]=weighted.mean(x[1:i], rev(r^(0:(i-1))))} return(m)} sd.ew=function(x,r,m){ sd=rep(NA,length(x)) for(i in 1:length(x)){ sd[i]=weighted.mean((x[1:i]-m[i])^2, rev(r^(0:(i-1))))} return(sd)} q=.97 MX=moy.ew(RX,q) SX=sd.ew(RX,q,MX) MY=moy.ew(RY,q) SY=sd.ew(RY,q,MY) plot(time(Y)[-1],SY,col="blue",type="l") lines(time(X2)[-1],SX,col="red")
And now we have something less erratic, so we can focus now on the interpretation.
It is also possible to look on the difference between those two series of volatility, areas in blue means that in Shanghai (again, where short-selling is forbidden) returns are more volatile than in Hong Kong, and areas in red are periods where returns are more volatile in Hong Kong,
a=time(X2)[which(time(X2)%in%time(Y))] b=SY[which(time(Y)%in%time(X2))]- SX[which(time(X2)%in%time(Y))] n=length(a) a=a[-n];b=b[-n] plot(a,b,col="black",type="l") polygon(c(a,rev(a)),c(pmax(b,0),rep(0,length(a))), col="blue",border=NA) polygon(c(a,rev(a)),c(pmin(b,0),rep(0,length(a))), col="red",border=NA)
So clearly, there is nothing clear that can be said… Sometimes, volatility is higher in Hong Kong, and sometimes, it is higher in Shanghai. But if we look at the price, instead of looking at volatility,
a=time(X2)[which(time(X2)%in%time(Y))] b=as.numeric(Y[which(time(Y)%in%time(X2)),4])- as.numeric(X2[which(time(X2)%in%time(Y))]) n=length(a) a=a[-n];b=b[-n] plot(a,b,col="black",type="l") polygon(c(a,rev(a)),c(pmax(b,0),rep(0,length(a))), col="blue",border=NA) polygon(c(a,rev(a)),c(pmin(b,0),rep(0,length(a))), col="red",border=NA)
Here, it looks like bans on short-selling creates bubbles. Might not not be a goodthing.
La vente à découvert, un moindre mal ?
Doit-on penser, comme l’avait dit Amartya Sen, que l’économie est une science morale ? La science économique (comme science ‘humaine’) a-t-elle un contenu moral per se, contrairement aux sciences ‘dures’ ? Si la morale est un choix de ce qui est bien et de ce qui est mal, on admettra que la loi de la gravitation universelle n’est bien ni mal. Pareil pour la division des cellules, ou le théorème de Pythagore. Pourtant en économie, il est fréquent d’évoquer une raison morale comme justification de tout et n’importe quoi… peut-être encore plus dès que l’on parle de finance…
Par exemple dans La Presse datant de fin septembre, on apprend que « les lois de l’offre et de la demande sont malheureusement faussées par le short selling […] Pourquoi permet-on de telles pratiques? Comme vous le savez la très grande proportion des short sellers ce sont des spéculateurs […] Malheureusement les autorités des marchés font la sourde oreille devant le gros bon sens d’une telle demande. Leur argument ? Oui! le «short selling» crée de la volatilité artificielle, mais l’interdire brimerait le libre jeu de l’offre et de la demande ».
Bref, la vente à découvert c’est mal, il faudrait l’interdire ! En Europe, La Tribune ou Le Monde en ont aussi parler, avec des arguments similaires.
Bruno Rémillard, que j’ai croisé en fin de semaine à un séminaire, me disait avoir envoyé à La Presse une lettre afin d’apporter des précisions (réponse écrite avec Nicolas Papageorgiou ), lettre restée pour l’instant sans réponse (et non publiée). Il reprenait en partie des résultats d’Alexandre Roch qui montrait que l’interdiction de vente à découvert pouvait générer des bulles spéculatives. Par exemple sur le graphique ci-dessus, emprunté à Alexandre, on a comparé les valeurs du titre d’Air China, échangé à la fois à Shanghai en bleu (où la vente à découvert est interdite) mais également à Hong Kong en rouge (où la vente à découvert est permise).
Je pense qu’il faudrait regarder plus attentivement pour montrer que le titre est plus volatile sur le marché où la vente à découvert n’est pas autorisée. Mais surtout, interdire la vente à découvert signifie encourager la création de bulles. Je ne suis pas sur que le jeu en vaille la chandelle….
De la difficulté de générer du hasard
J’avais parlé cet été de l’utilisation des chiffres du us treasury balance comme générateur (involontaire) de nombres aléatoires. Cet été, un débat intéressant a secoué le Québec, suite à la plainte déposée par deux joueurs (compulsifs) de l’Extra. L’histoire est raconté en détails sur le site de canoe.ca, ou de lapresse.ca. Le jeu est assez simple. L’histoire un peu moins (c’est pour cela que j’ai préféré renvoyer vers deux sources racontant la même histoire). Pour essayer d’expliquer simplement ce qui s’est passé revenons un peu sur le jeu. Tous les jours, Loto-Québec tire au hasard un nombre de 7 chiffres.
Comme on le voit ci-dessus, le but est d’avoir les premiers (ou les derniers) chiffres bons.
Quand on achète un ticket, un numéro nous est attribué… et on espère qu’il sera tiré le soir. Mais le soucis est que les joueurs qui intentent un procès n’achetaient pas un mais dix tickets. Et étrangement, à chaque fois, ils avaient un ticket gagnant, et un seul.
En fait, quand on achetait dix combinaisons d’extra, il semble que le premier chiffre sans remise parmi 0 à 9, les cinq suivants avec remise parmi 0 à 9 et le dernier sans remise parmi 0 à 9 (pour reprendre la formulation qu’a utilisé Christian Léger de l’UdM pour me l’expliquer). Sur un achat, on ne voit rien, on a des nombres uniforme… mais avec répétition, ça se complique.
En fait, c’est toujours quand on étudie la dynamique des générateurs de nombres aléatoires que l’on a des surprises.
Par exemple sur les tirages de l’extra (et non pas les numéros des combinaisons des billets, car on ne les trouve pas sur internet, il faudrait acheter les billets pour refaire l’expérience). On peut récupérer l’historique des tirages en ligne, sur le site de loto-quebec.com.
Quand on achète un ticket, un numéro nous est attribué… et on espère qu’il sera tiré le soir. Mais le soucis est que les joueurs qui intentent un procès n’achetaient pas un mais dix tickets. Et étrangement, à chaque fois, ils avaient un ticket gagnant, et un seul.
En fait, quand on achetait dix combinaisons d’extra, il semble que le premier chiffre sans remise parmi 0 à 9, les cinq suivants avec remise parmi 0 à 9 et le dernier sans remise parmi 0 à 9 (pour reprendre la formulation qu’a utilisé Christian Léger de l’UdM pour me l’expliquer). Sur un achat, on ne voit rien, on a des nombres uniforme… mais avec répétition, ça se complique.
En fait, c’est toujours quand on étudie la dynamique des générateurs de nombres aléatoires que l’on a des surprises.
Par exemple sur les tirages de l’extra (et non pas les numéros des combinaisons des billets, car on ne les trouve pas sur internet, il faudrait acheter les billets pour refaire l’expérience). On peut récupérer l’historique des tirages en ligne, sur le site de loto-quebec.com.
> extra=read.table("http://freakonometrics.blog. free.fr/public/data/extra-loto-quebec.csv",sep=";",header=TRUE) > n=1940 > extra=extra[1:(n-1),1:2] > X=extra[,2]/10000000
Si on regarde les 1940 derniers tirages (car avant, il semble qu’il n’y avait pas 7 mais 6 chiffres), on a la série suivante,
Et un test d’uniformité semble confirmer l’intuition que les nombres sont uniformément répartis,
Mais générer du hasard, ça n’est pas que ça… Et sur cette série, si on fait des tests sur l’indépendance des observations, on est un peu surpris,
> library(lawstat) > runs.test(X, plot.it = FALSE, alternative = "two.sided") Runs Test - Two sided data: X Standardized Runs Statistic = -1.8, p-value = 0.07273 > runs.test(X, plot.it = FALSE, alternative = "positive.correlated") Runs Test - Positive Correlated data: X Standardized Runs Statistic = -1.8, p-value = 0.03637 > library(FitAR) > LjungBoxTest(qnorm(X), k=0, lag.max=12, + StartLag=1, SquaredQ=FALSE) m Qm pvalue 1 7.2 0.0075 2 9.2 0.0102 3 10.9 0.0125 4 10.9 0.0276 5 15.3 0.0090 6 15.4 0.0171 7 16.1 0.0244 8 16.3 0.0380 9 16.7 0.0538 10 20.6 0.0240 11 20.8 0.0357 12 20.8 0.0533
ACT2040, introduction aux modèles linéaires généralisés
On commencera ce mardi les GLM, après avoir introduit les lois exponentielles (qui ont du être revues en démonstration vendredi dernier). La notation utilisée sera que la loi (densité ou fonction de probabilité) de 
sera de la forme
Pour un complément plus exhaustif, je renvoie au chapitre en ligne.
- Le modèle linéaire (Gaussien)
Le modèle de base est le modèle Gaussien que l’on avait revu au dernier cours,
> X=c(1,2,3,4) > Y=c(1,2,5,6) > base=data.frame(X,Y) > reg1=lm(Y~1+X,data=base) > nbase=data.frame(X=seq(0,5,by=.1)) > Y1=predict(reg1,newdata=nbase)
Pour une prédiction (unique), on obtient la prédiction suivante
Le code pour une telle représentation est le suivant
> plot(X,Y,pch=3,cex=1.5,lwd=2,xlab="",ylab="") > lines(nbase$X,Y1,col="red",lwd=2) > u=2 > mu=predict(reg1)[2] > sigma=summary(reg1)$sigma > y=seq(0,7,.05) > loi=dnorm(y,mu,sigma) > segments(u,y,loi+u,y,col="light green") > lines(loi+u,y) > abline(v=u,lty=2) > points(X[2],Y[2],pch=3,cex=1.5,lwd=2) > points(X[2],predict(reg1)[2],pch=19,col="red") > arrows(u-.2,qnorm(.05,mu,sigma), + u-.2,qnorm(.95,mu,sigma),length=0.1,code=3,col="blue")
On peut multiplier les prédictions, en se basant sur l’hypothèse d’homoscédasticité (la variance sera alors constante)
Mais on peut aller plus loin
- Le modèle linéaire généralisé
Plusieurs modèles peuvent etre estimés, en changeant la loi de la variable à expliquer, et la fonction lien,
> reg2=glm(Y~1+X,data=base,family=poisson(link="identity")) > Y2=predict(reg2,newdata=nbase,type="response") > reg3=glm(Y~1+X,data=base,family=poisson(link="log")) > Y3=predict(reg3,newdata=nbase,type="response") > reg4=glm(Y~1+X,data=base,family=gaussian(link="log")) > Y4=predict(reg4,newdata=nbase,type="response") > sigma=sqrt(summary(reg4)$dispersion)
Pour le modèle Poissonnien avec un lien identité, on obtient
On obtient ainsi une variance qui augmente avec la prédiction,
Pour une régression de Poisson avec un lien logarithmique,
i.e. pour nos quatre prédictions
On peut comparer avec une prédiction d’un modèle Gaussien avec un lien logarithmique,
i.e. pour les quatre prédictions
