Freakonometrics

To content | To menu | To search

Friday, December 21 2012

The Maya prediction was correct...

...it is the end of a World. This blog is now officially dead. See you from now on on my new blog,

http://freakonometrics.hypotheses.org/

Sunday, January 9 2011

From one extreme (0) to another (1): challenge failed, but who cares...

Just after arriving in Montréal, at the beginning of September, I discussed statistics of my blog, and said that it might be possible - or likely - that by new year's Eve, over a million page would have been viewed on my blog (from Google's counter, here). By the end of October (here) I was very optimistic, but mi-December (here) the challenge was likely to be failed. An indeed, the million page target was hit one week after, on January 8th,

base=read.table("http://freakonometrics.blog.free.fr/public/data/million1.csv",
sep="\t",header=TRUE)
X1=cumsum(base$nombre)
X0=X1
base=read.table("http://freakonometrics.blog.free.fr/public/data/million2.csv",
sep="\t",header=TRUE)
X2=cumsum(base$nombre)
X=X1+X2
 
D0=as.Date("08/11/2008","%d/%m/%Y")
D=D0+1:length(X1)
plot(D,X1,xlim=c(as.Date("08/06/2010","%d/%m/%Y"),as.Date("08/02/2011","%d/%m/%Y")),
ylim=c(800000,1050000))
abline(h=1000000,col="red")
abline(v=as.Date("01/01/2011","%d/%m/%Y"),col="red")
points(D,X,col="blue")

Again, the black points were from the previous blog (http://blogperso.univ-rennes1.fr/arthur.charpentier/) which was transferred to that new one (http://freakonometrics.blog.free.fr) this Autumn. So I just sum up the stats to get the blue points. At each date, I fit an ARIMA, and use it to make forecast the total number of pages viewed on January 1st, and calculate the probability to reach a million page viewed at that date (using a Gaussian ARIMA model). Actually, here, I changed a little bit the challenge, and asked "what would have been the probability to reach a million page viewed on January 1st, and on January 8th" ?

kt=which(D==as.Date("01/06/2010","%d/%m/%Y"))
Xbase=X
X=X1+X2
P1=P2=rep(NA,(length(X)-kt)+7)
for(h in 0:(length(X)-kt+7)){
model <- arima(X[1:(kt+h)],c(7 ,1,7),method="CSS")
forecast <- predict(model,200)
u=max(D[1:kt+h])+1:300
if(min(u)<=as.Date("01/01/2011","%d/%m/%Y")){
k=which(u==as.Date("01/01/2011","%d/%m/%Y"))
(P1[h+1]=1-pnorm(1000000,forecast$pred[k],forecast$se[k]))}
k=which(u==as.Date("08/01/2011","%d/%m/%Y"))
(P2[h+1]=1-pnorm(1000000,forecast$pred[k],forecast$se[k]))
}
The red curve is the probability to reach 1 million viewed on January 1st (as done earlier, using an ARIMA projection). The blue one is the probability to reach 1 million viewed one week after, on January 8th.

 and here is the difference between probabilities,

The flat part at the beginning of November corresponds to the bump that was observed on the initial graph. But then, the slope was too low, and in December, the challenge was failed... Obviously, looking at statistics during a blog migration is not a bright idea...

Friday, January 7 2011

Lecture de blog(s): quelques précautions

En finissant mon billet sur le million de pages vues en deux ans sur ce blog, je me suis rendu compte de l'énormité de la chose.... et je me suis dit qu'il convenait de faire une mise en garde (ça en fait des lecteurs potentiels pour mes élucubrations).
L'écriture d'un blog est un exercice particulier. Contrairement à une page ouebe régulièrement mise à jour, dans un blog, on ne revient pas en arrière. Ou plutôt pas complètement....Et les billets s'accumulent.
Si un billet contient des fautes, je peux le reprendre, et je raye les fautes. Mais généralement, des lecteurs font des remarques qui éclairent le billet, et le font gagner en clarté.
Contrairement aux blogs très fréquentés où les commentateurs sont là pour critiquer, souvent de manière malsaine, ou pour discuter entre eux, il y a peu de commentaires sur mon blog, et souvent, je trouve qu'ils corrigent bien mes imprécisions ou mes coquilles (comme ici). 
Moralité, je recommande aux lecteurs intéressés de ne pas se contenter des billets, mais d'aller faire un tour sur les commentaires.
Voilà, maintenant au boulot, il y a des billets à finir !

Wednesday, October 27 2010

Prétentieux, va !

Depuis quelques semaines, mes billets sont presque exclusivement en anglais... En fait, même si je fais cours en français, et que la plupart de mes collègues ici parlent français, je me rends compte qu'on parle naturellement anglais.... en tous les cas plus naturellement qu'en France. Les séminaires sont toujours en anglais, dès qu'on discute à plus de 3, la probabilité qu'une personne ne parle pas français est tellement grande qu'on parle le plus souvent en anglais....

Finalement, comme toujours quand je ne suis pas en France, je trouve naturel d'écrire en anglais... (ok, d'essayer d'écrire quelque chose qui ressemble à de l'anglais). Mais hier soir, ma femme m'a fait remarqué que c'était prétentieux d'écrire en anglais... Et effectivement, quand je traîne sur les blogs francophones, presque tous sont exclusivement en français.... Bref, je me remets à douter.
J'ai vu l'autre jour que David Monniaux s'en prenait sur son blog à un édito de Jacques Julliard (ici), et d'ailleurs David publie de plus en plus en anglais ces derniers temps (en tous les cas dès que l'on quitte les sujets franco-français). J'ai vu aussi que de plus en plus de doctorants blogueurs se lancent en anglais...
Le débat reste ouvert, et je suis preneur de toutes les remarques....

A million ? what are the odds...

50 days ago, I published a post, here, on forecasting techniques. I was wondering what could be the probability to have, by the end of this year, one million pages viewed (from Google Analytics) on this blog. Well, initially, it was on my blog at the Université de Rennes 1 (http://blogperso.univ-rennes1.fr/arthur.charpentier/), but since I transfered the blog, I had to revise my code. Initially, I had that kind of graphs,

and when I look at the cumulative distribution of the number of pages viewed on January first, I had

while for the distribution of the time I should read this million (the dual problem), I obtained

and I said that I should have around 35% chance to reach the million pages viewed by the end of this year.
Here is the updated graph, with the blog à Université de Rennes 1 (still in black) and the one here (in blue, where I add the two blogs together).

Actually, I decided to look at the evolution of the probability to reach the million by New Year's Eve...

The code looks like that,
> base=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/million2.csv",
+ sep=";",header=TRUE)
> X2=cumsum(base$nombre)
> X=X1+X2
> kt=which(D==as.Date("01/06/2010","%d/%m/%Y"))
> D0=as.Date("08/11/2008","%d/%m/%Y")
> D=D0+1:length(X1)
> P=rep(NA,(length(X)-kt)+1)
> for(h in 0:(length(X)-kt)){
+ model  <- arima(X[1:(kt+h)],c(7 ,   # partie AR
+                     1,    # partie I
+                     7),method="CSS")   # partie MA
+ forecast <- predict(model,200)
+ u=max(D[1:kt+h])+1:300
+ k=which(u==as.Date("01/01/2011","%d/%m/%Y"))
+ (P[h+1]=1-pnorm(1000000,forecast$pred[k],forecast$se[k]))
+ }
It has been a bit tricky, since I wanted an automatic fit of the ARIMA process, meaning that I had to assess a priori the orders of the ARIMA process. And I had numerical problems, since we got non stationary AR part at least at one period of time considered.... So finally I used here the CSS method which uses conditional-sum-of-squares to find starting values in the optimization procedure.

Actually, if we consider a classical descritption of traders, it looks like I act as a trader (dealing with millions and forgetting about real people): it is the same here, I do not know what a million means, I cannot imagine 250,000 visitors looking at that blog... But I can still do the maths. Anyway, a million is huge when I start to think about it... but perhaps I should not... I cannot possibility imagine that so many people might find interesting my mathematical lucubration*....
* initially I was looking for the analogous of "élucubration" in French, meaning "divagation, absurd theory" (the proper translation might be "rantings" (here) , "ravings" (here) or "wild imagining" (everywhere else here or there)). When I asked Google for a possible translation (here), I got "lucubration" which means "composed by night; that which is produced by meditation in retirement". Well, it was not initially what I intended to say, but since I usually work on my blog during the night, when I got awake by one of the girls, I decided to keep this word.... At least, I learnt something today, appart for the code mentioned above....

Thursday, September 9 2010

Le million ! le milllion !

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 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...

Sunday, August 22 2010

La multiplication des blogs ?

Le transfert de blog suit son petit bonhomme de chemin... Le blog

est le plus complet, et celui qui sert désormais de référence, compte tenu des déconvenues sur le blog de Rennes 1. J'ai malgré tout souhaité maintenir
qui est le blog intial (sur la mise en page, les billets sont plus propres que sur les billets transférés). Ce blog continuera d'exister, mais ne contiendra que les billets de recherche et d'enseignement. Enfin, un nouveau né,
- hébergé chez Hypothèses, portail de carnets de recherche - je mettrai les billets de recherche uniquement.

Wednesday, April 22 2009

Bilan du blog, un an après

Bon visiblement mon blog fait partie de la blogosphère... ce qui apporte une certaine satisfaction, mais aussi une certaine appréhension. Bref, après un an (à peu près) d'existence, je voulais faire un court bilan.

Je trouvais l'idée du blog plus souple pour l'enseignement qu'une page perso classique, permettant de répondre aux questions au fur et à mesure (sans passer par des réponses individuelles par mail), et en mettant en ligne des compléments sur des choses que j'a vu trop rapidement en cours. Je ne sais pas si la solution est meilleure qu'un forum, mais pour ma part, ça me parait plus simple à gérer. Pour la recherche, j'ai longtemps été sceptique, et je crois que je le suis toujours. Le blog est moins clair qu'une page résumant l'ensemble des travaux de recherche (et offrant une certaine vue d'ensemble), même s'il offre la possibilité de compléter les papiers par des bases de données, ou des codes sources. Et parallèlement, je suis ravi par exemple de voir (ici) que Stéphane Loisel rebondi sur mon blog pour écrire un court article passionnant.

En revanche c'est incroyablement coûteux en temps, même si ce n'est pas du temps de perdu, car cela correspond soit au travail qui nous est demandé de faire en temps qu'enseignant chercheur.

Et sinon, pour reprendre un commentaire que j'avais lu sur http://www.mafeco.fr/, et que je partage assez "on a l’impression (en tout cas, c’est mon impression) de participer à ce qui ressemble aux prémices d’une révolution en ce qui concerne la production et la circulations d’idées et de connaissances scientifiques. Je ne sais pas dans quelle mesure les blogs économiques vont se développer en France, mais je suis de plus en plus convaincu que le modèle institutionnel sur lequel est construit la production et circulation d’idées économiques (le système des revues) est voué à évoluer avec internet et les blogs." Cette idée est d'ailleurs énormément développée en ce moment sur la toile, par exemple sur http://blogs.reuters.com/felix-salmon/, qui s'interroge sur le peu d'écoblogueurs (si le mot a un sens) en Allemagne (mais qui est parfois étendu au cas français sur internet, même si "économiste" est souvent réduit à des universitaires qui donnent leur avis sur la politique économique actuelle). En tous les cas, sur l'importance croissante des blogs, je retiendrais l'histoire racontée ici, qui - pour résumer - nous apprend que sur son blog (ici), Hal Turner a obtenu les résultats des stress-tests des banques américaines (qui ne devaient être publiées par les autorités américaines que fin mai) qui annonçait que 16 des 19 banques testées seraient techniquement insolvables et, pour traduire les propos de Hal, "pour le dire franchement, le système bancaire américain serait en train de s'effondrer totalement". Le lendemain matin, lundi, le Dow Jones chute de 3,56%, et les banques s'effondrent, en particulier Bank of America (-24% ), Citigroup (-20% ), JPMorgan (-11%), etc. Bref, pas mal de monde donne beaucoup d'importance à ce blog, même si certaines banques mentionnées ne font pas partie de la liste des banques qui doivent fournir les résultats de leurs stress tests (mentionnée ici). Bref, le trésor américain a démenti les informations officiellement (ici ou ) et donc donné une importance démesurée à ce blog. Bref, quand on fait un blog, on est forcément très narcissique, et on a envie de croire que ce qu'on peut faire va changer le monde...

Friday, February 27 2009

Record de visites sur mon blog !

Bon, si je crois les statistiques de Google Analytics, en moins de 4 mois, j'ai dépassé les 50,000 pages vues et les 15,000 visites (dont plus de 10,000 uniques)... qui aurait pu penser que l'économétrie et l'actuariat intéresseraient autant de monde ! Les mots clés principaux qui font arriver sur ce blog via google sont "chain ladder", "économétrie robuste", "probabilités et statistiques" ou encore "analyse des données"... Bon, pour être (complètement) honnête, les pics ont été observés lorsque j'ai annoncés les mouvements de grèves fin janvier/début février...


Sinon l'évolution via les deux compteurs du site (qui n'ont pas les mêmes décomptes)