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