Freakonometrics

Un billet pour reprendre (et compléter) le code vu hier, en cours de modèles de prévision. On a toujours pour but de prévoir l'évolution du trafic autoroutier, en France, sur l'autoroute A7.

> autoroute=read.table(
+ "http://freakonometrics.blog.free.fr/public/data/autoroute.csv",
+ header=TRUE,sep=";")
> a7=autoroute$a007
> A7=ts(a7,start = c(1989, 9), frequency = 12)

Comme je l'ai dit en cours, les constantes dans les modèles ARIMA, a complexifie inutilement les notations, donc pour faire des choses proches de ce qu'ai pu faire dans le cours, on va centrer la série, et considérer

> A7d=A7-mean(A7)
> plot(A7d,xlim=c(1989,2000))

On avait vu qu'un modèle pouvait convenir, donc on va l'utiliser pour faire de la prévision.

>  fit=arima(A7d,order=c(p=12,d=0,q=0),
+      include.mean =FALSE)
>  summary(fit)
Series: A7d
ARIMA(12,0,0) with zero mean
 
Coefficients:
ar1      ar2      ar3     ar4      ar5     ar6
      0.0981  -0.0353  -0.0375  0.0530  -0.0934  0.0598
s.e.  0.0630   0.0623   0.0600  0.0606   0.0618  0.0568
ar7     ar8     ar9     ar10    ar11    ar12
      -0.0839  0.0278  0.0072  -0.1080  0.1578  0.7887
s.e.   0.0619  0.0629  0.0627   0.0629  0.0637  0.0623
 
sigma^2 estimated as 12823951:  log likelihood=-827.24
AIC=1678.48   AICc=1683.6   BIC=1710.23
 
In-sample error measures:
ME         RMSE          MAE          MPE
257.8363135 3581.0544202 2335.0686684    1.7258390
MAPE         MASE
34.7942409    0.2603139 

La fonction permettant de faire de la prévision est alors

> horizon=48
> library(forecast)
> A7hat=forecast(fit,h=horizon)

On a alors un ensemble de séries temporelles, la première étant la prévision proprement dite , i.e. , où désigne une espérance linéaire (ou en terme technique, une projection sur l'adhérence de l'ensemble des des combinaisons linéaires du passé). On a ensuite les intervalles de confiance estimés, à partir d'une hypothèse de normalité, en utilise un estimateur de ,

> A7hat
Point    Lo 80   Hi 80    Lo 95    Hi 95
Oct 1996    -7276.44 -11865.7 -2687.1 -14295.1  -257.70
Nov 1996   -10453.79 -15065.1 -5842.4 -17506.1 -3401.39
Dec 1996    -6970.67 -11583.5 -2357.8 -14025.3    84.02
Jan 1997   -13404.09 -18021.2 -8786.9 -20465.4 -6342.78

On peut alors représenter ces valeurs graphiquement, comme on l'a fait en cours

>  (TPS=tsp(A7hat$mean))
[1] 1996.750 2000.667   12.000
>  lines(A7hat$mean,col="red",lwd=2)
>  temps=seq(from=TPS[1],to=TPS[2],by=1/TPS[3])
>  lines(temps,A7hat$lower[,1],col="blue",lwd=2)
>  lines(temps,A7hat$lower[,2],col="purple",lwd=2)
>  lines(temps,A7hat$upper[,1],col="blue",lwd=2)
>  lines(temps,A7hat$upper[,2],col="purple",lwd=2)
>  abline(v=TPS[1]-.5/TPS[3],col="green")

En allant un peu plus loin, on peut faire plus joli, en utilisant des aires colorées afin d'illustrer une région de confiance

> plot(A7d,xlim=c(1989,2000))
> polygon(c(temps,rev(temps)),c(A7hat$lower[,2],
+ rev(A7hat$upper[,2])),col="yellow",border=NA)
> polygon(c(temps,rev(temps)),c(A7hat$lower[,1],
+ rev(A7hat$upper[,1])),col="orange",border=NA)
> lines(temps,A7hat$mean,col="red")
> abline(v=TPS[1]-.5/TPS[3],col="green")

Maintenant, si on revient sur notre série initiale, l'autocorrélogramme pourrait laisser penser qu'il y a une racine unité saisonnière,

> plot(acf(A7d,lag=36))

Donc on peut tenter de différencier à l'ordre 12,

> plot(acf(diff(A7d,12),lag=36))

On pourrait avoir l'impression que la série ainsi différenciée est un bruit blanc. On peut faire un test de Box-Pierce pour tester cette hypothèse,

> BP=function(h) Box.test(diff(A7d,12),lag=h,
+ type='Box-Pierce')$p.value
> plot(1:24,Vectorize(BP)(1:24),type='b',
+ ylim=c(0,1),col="red")
> abline(h=.05,lty=2,col="blue")

On pourrait valider ici l'hypothèse de bruit blanc. On va donc ajuster un (pur) modèle ARIMA saisonnier,

> fit2=arima(A7d, order = c(p=0, d=0, q=0),
+ seasonal = list(order = c(0, 1, 0), period = 12))

Si on compare les deux modèles,

> AIC(fit1)
[1] 1511.815
> AIC(fit2)
[1] 1245.374

Le second modèle ayant un critère d'Akaike plus faible, on aurait tendance  à préférer ce modèle au premier. Mais cette méthode de comparaison n'est pas forcément la plus pertinent, si l'objet est de faire une prévision à court terme (disons quelques mois). Pour cela, on va plutôt enlever les dernières observations, ajuster un modèle sur la série restreinte, et comparer les prévisions, avec ce qui a été réellement observé. Disons qu'on va faire du backtesting sur les 9 derniers mois,

> T=length(A7d)
> backtest=9
> subA7d=ts(A7d[-((T-backtest+1):T)],
+ start = c(1989, 9), frequency = 12)
> fit1s=arima(subA7d,order=c(p=12,d=0,q=0),
+ include.mean =FALSE)

pour le premier modèle, et

> fit2s=arima(subA7d, order = c(p=0, d=0, q=0),
+ seasonal = list(order = c(0, 1, 0), period = 12))

pour le second. Ensuite, on peut faire de la prévision, et comparer les prévisions avec les données écartées pour l'estimation du modèle.,

> XA7p1=forecast(fit1s,h=backtest)
> XA7p2=forecast(fit2s,h=backtest)
> A7dobs=A7d[((T-backtest+1):T)]
> (M=cbind(observé=A7dobs,modèle1=as.numeric(XA7p1$mean),
+ modèle2=as.numeric(XA7p2$mean)))
observé     modèle1    modèle2
[1,] -13610.071 -12264.7138 -13558.071
[2,] -10108.071 -12236.6227 -10410.071
[3,]  -8034.071  -9333.8365  -9504.071
[4,]   3899.929   4071.9313   4108.929
[5,]   2009.929    901.3745   2661.929
[6,]   5253.929   7610.9954   5866.929
[7,]  25357.929  29232.2317  30317.929
[8,]  35229.929  29641.5172  31737.929
[9,]   4341.929   5603.8008   4940.929

On peut alors retenir un critère classique, comme l'erreur quadratique, et donc comparer la somme des carrées des erreurs,

> sum((M[,1]-M[,2])^2)
[1] 62677238
> sum((M[,1]-M[,3])^2)
[1] 40253827

Le second modèle semble ici avoir été meilleur, sur les 9 derniers mois. On peut aussi visualiser les erreurs de prévision,

> (TPS=tsp(subA7d))
[1] 1989.667 1995.917   12.000
> temps=seq(TPS[2]+1/TPS[3],by=1/TPS[3],length=backtest)
> rangex=c(TPS[1],TPS[2]+backtest/TPS[3])
> plot(subA7d,type="l",xlim=rangex)
> abline(v=T-11.5,lty=2)
> polygon(c(temps,rev(temps)),c(A7dobs,rev(XA7p1$mean)),
+ col="yellow",border=NA)
> lines(temps,A7dobs)
> lines(temps,XA7p1$mean,col="red")

pour le premier modèle,

et

> polygon(c(temps,rev(temps)),c(A7dobs,rev(XA7p2$mean)),
+ col="yellow",border=NA)
> lines(temps,XA7p2$mean,col="red")

pour le second,

Ce soir, @imparibus proposait sur son blog un billet passionnant sur l'élection présidentielle en France, avec des graphiques superbes, basé sur un lissage temporel des résultats des sondages pour le premier tour de l'élection présidentielle en France (qu'il a fait sous excel, on saluera la performance !). En reprenant les données du site http://www.lemonde.fr/, on peut obtenir la même chose assez simplement (j'ai repris - ou presque - les couleurs utilisées sur le site, François Hollande en rose, Nicolas Sarkozy en bleu, François Bayrou en orange, et Marine Le Pen en noir...). Pour commencer, le code en R pour importer les données et les manipuler est le suivant

sondage=read.table("http://freakonometrics.blog.free.fr/
public/data/sondage1ertour2012.csv",
header=TRUE,sep=";",dec=",")
sondage[36,1]="11/03/12"
sondeur=unique(sondage$SONDEUR)
sondage$date=as.Date(as.character(sondage$DATE),
"%d/%m/%y")

(je fais manuellement une correction de date car je me suis trompé en saisissant les chiffres à la main). A partir de là, on peut reprendre l'idée de faire une régression locale pour trouver une tendance,

CL=brewer.pal(6, "RdBu")
couleur=c(CL[1],CL[6],"orange","grey")
datefin=as.Date("2012/04/12")
k=1
plot(sondage$date,sondage[,k+2],col=couleur[k],
pch=19,cex=.7,xlab="",ylab="",ylim=c(0,40),
xlim=c(min(sondage$date),datefin))
rl=lowess(sondage$date,sondage[,k+2])
lines(rl,lwd=2,col=couleur[k])
for(k in 2:4){
points(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7)
rl=lowess(sondage$date,sondage[,k+2])
lines(rl,lwd=2,col=couleur[k])
}

Le code a l'air long mais il faut définir les couleurs, générer une fenêtre graphique, mettre les régressions locales dedans, etc. En tant que telle, la commande qui permet de faire le lissage est juste

rl=lowess(sondage$date,sondage[,k+2])

Bon, on peut d'ailleurs faire toutes sortes de lissages, avec des splines par exemple (ce que j'ai davantage tendance à utiliser),

D=data.frame(date=min(sondage$date)+0:148)
library(splines)
k=1
plot(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
prl=predict(rs,newdata=D)
prlse=sqrt(prl/100*(1-prl/100)/1000)*100
polygon(c(D$date,rev(D$date)),c(prl+2*prlse,
rev(prl-2*prlse)),col=CL[3],border=NA)
lines(D$date,prl,lwd=2,col=CL[1])
points(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7,)
for(k in 2:4){
points(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7)
}

avec, là encore, la commande suivante pour lisser

rs=lm(sondage[,k+2]~bs(date),data=sondage)

Cette fois le code est un peu plus long, parce que j'ai tracé un intervalle de confiance à 95% (intervalle de confiance classique, en supposant que 1000 personnes ont été interrogées, comme évoqué dans d'anciens billets),

On a ici la courbe suivante

avec l'intervalle de confiance pour François Hollande, mais on peut faire la même chose pour Nicolas Sarkozy,

k=2
plot(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
prl=predict(rs,newdata=D)
prlse=sqrt(prl/100*(1-prl/100)/1000)*100
polygon(c(D$date,rev(D$date)),c(prl+2*prlse,
rev(prl-2*prlse)),col=CL[4],border=NA)
lines(D$date,prl,lwd=2,col=CL[6])
points(sondage$date,sondage[,k+2],col=
couleur[k],pch=19,cex=.7,)
for(k in c(1,3:4)){
points(sondage$date,sondage[,k+2],col=
couleur[k],pch=19,cex=.7)
}

Mais comme auparavant, on peut aussi visualiser les "tendances" des quatre candidats en tête,

Voilà pour le début. En fait, idéalement, on voudrait faire un peu de prévision... Le hic est que les code usuel pour faire de la prévision (avec des ARIMA, du lissage exponentiel ou tout autre modèle classique) nécessitent de travailler avec des séries temporelles, telles qu'elles sont classiquement définies, c'est à dire "régulièrement espacées dans le temps".

Je ne vais pas commencer à réclamer plus de sondages, loin de moins cette idée, mais pour mon modèle, j'avoue que j'aurais préféré avoir plus de points. Beaucoup plus de points. Un sondage par jour aurait été idéal en fait...

Qu'à cela ne tienne, on va simuler des sondages. L'idée est simple: on a une tendance qui nous donne une probabilité jour par jour (c'est exactement ce qui a été calculé pour tracer les courbes), via

D=data.frame(date=min(sondage$date)+0:148)
prl=predict(rs,newdata=D)

On va ensuite utiliser une hypothèse de normalité multivariée sur nos 4 candidats, auquel j'en ai en rajouté un cinquième fictif (mais ça ne servait à rien) en utilisant l'analogue multivariée de l'expression suivante (évoqué dans le premier billet de l'année sur les élections)

On peut alors générer, jour après jour, des sondages, en simulant des vecteurs Gaussiens. Je vais les générer indépendants les uns des autres, parce que c'est plus simple, et que cela ne me semble pas être une trop grosse hypothèse. Pour générer un ensemble de sondages, on utilise la fonction suivante

library(mnormt)
simulation=function(S){
proba=c(S,100-sum(S))/100
variance= -proba%*%t(proba)
diag(variance)=proba*(1-proba)
variance=variance/1000
return(rmnorm(1,proba[1:4],variance[1:4,1:4]))}
simulsondages=function(M){
Z=rep(NA,ncol(M))
for(i in 1:nrow(M)){
Z=rbind(Z,simulation(M[i,]))
}
return(Z[-1,])}
prediction4=matrix(NA,nrow(D),4)
for(k in 1:4){
rs=lm(sondage[,k+2]~bs(date),data=sondage)
prediction4[,k]=predict(rs,newdata=D)
}

Par exemple, si on la fait tourner une fois, on obtient le graphique suivant

set.seed(1)
S4=simulsondages(prediction4)
S100=100*S4
k=1
plot(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
for(k in 2:4){
points(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7)
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
}

(les courbes lissées sont celles obtenues sur les vrais sondages). Là on peut être heureux, parce qu'on a des vraies séries temporelles. On peut alors faire de la prévision, par exemple en faisant du lissage exponentiel (optimisé),

library(forecast)
k=1
X=S100[,k]
ETS=ets(X)
F=forecast(ETS,h=60)
fdate=max(D$date)+1:60
k=1
plot(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
for(k in 2:4){
points(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7)
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),lwd=2,col=couleur[k])
}
polygon(c(fdate,rev(fdate)),c(as.numeric(F$lower[,2]),
rev(as.numeric(F$upper[,2]))),col=CL[3],border=NA)
polygon(c(fdate,rev(fdate)),
c(as.numeric(F$lower[,1]),rev(as.numeric(F$upper[,1]))),
col=CL[2],border=NA)
lines(fdate,as.numeric(F$mean),lwd=2,col=CL[1])

Là encore le code peut paraître long, mais c'est surtout la partie associée au graphique qui prend de la place (par soucis purement esthétique, on passe du temps sur les codes graphiques depuis le début). Le code qui modélise la série, et qui la projette, est ici donné par les deux lignes suivantes

ETS=ets(X)
F=forecast(ETS,h=60)

Pour François Hollande, sur la simulation des sondages passés que l'on vient d'effectuer, on obtient la prévision suivante

On peut bien entendu faire une projection similaire pour Nicolas Sarkozy,

k=2
X=S100[,k]
ETS = ets(X)
F=forecast(ETS,h=60)
k=1
plot(D$date,S100[,k],col=couleur[k],pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
for(k in 2:4){
points(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7)
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
}
polygon(c(fdate,rev(fdate)),c(as.numeric(F$lower[,2]),
rev(as.numeric(F$upper[,2]))),col=CL[4],border=NA)
polygon(c(fdate,rev(fdate)),
c(as.numeric(F$lower[,1]),rev(as.numeric(F$upper[,1]))),
col=CL[5],border=NA)
lines(fdate,as.numeric(F$mean),lwd=2,col=CL[6])

On notera que je ne me suis pas trop fatigué sur la projection: autant sur la génération de sondages passés on a tenu compte de corrélation (i.e. si un candidat a un score élevé, ça se fera au détriment des autres, d'où la corrélation négative utilisée dans la simulation), autant ici, les projections sont faites de manière complétement indépendantes. En se fatiguant un peu plus, on pourrait se lancer dans un modèle vectoriel gaussien. Ou mieux, on pourrait travailler sur des processus de Dirichlet (surtout que R a un package dédié à ce genre de modèle) parce qu'on travaille depuis le début sur des taux, comme évoqué auparavant. Mais commençons par faire simple pour un premier billet rapide sur ce sujet.

On peut ensuite s'amuser à générer plusieurs jeux de sondages, et de lancer des prévisions dessus,

set.seed(1)
for(sim in 1:20){
Ssim=simulsondages(prediction4)
F1=forecast(ets(100*Ssim[,1]),h=60)
F2=forecast(ets(100*Ssim[,2]),h=60)
lines(fdate,as.numeric(F1$mean),col=CL[1])
lines(fdate,as.numeric(F2$mean),col=CL[6])}

Une fois qu'on a fait tout ça, on a presque fini. On peut regarder au 22 avril qui est en tête (voire, par scénario, calculer la probabilité qu'un des deux candidats soit en tête), c'est à dire au soir du 1er tour

set.seed(1)
VICTOIRE=rep(NA,1000)
for(sim in 1:1000){
Ssim=simulsondages(prediction4)
F1=forecast(ets(100*Ssim[,1]),h=60)
F2=forecast(ets(100*Ssim[,2]),h=60)
VICTOIRE[sim]=(F1$mean[30]>F2$mean[30])}
mean(VICTOIRE)

Et voilà. Ah oui, je n'ai pas laissé le résultat. Tout d'abord parce que je suis un statisticien, pas un prédicateur. Mais aussi parce que si ce genre de pronostic amuse des gens... ils n'ont qu'à se mettre à R pour faire tourner les codes !

I recently received an email about forecasting and rules of thumb. "Dans la profession [...] se transmet une règle empirique qui voudrait que l'on prenne un historique du double de l'horizon de prévision : 20 ans de données pour une prévision à 10 ans, etc... Je souhaite savoir si cette règle n'aurait pas, par hasard, un fondement théorique quitte à ce que le rapport ne soit pas de 2 pour 1, mais de 3 pour 1 ou de 1 pour 1 par exemple." To summarize briefly, the rule is to consider a 2-1 ratio for the period of observation vs. forecast horizon. And the interesting question is if there are justifications for such a rule...

• It depends on the series

library(forecast)
X = AirPassengers
ETS = ets(X)
plot(forecast(ETS,h=length(X)/2))

• It depends on the model

50 days ago (here), I was supposed to be very optimistic about the probability that I could reach a million viewed pages on that blog (over a bit more than two years). Unfortunately, the wind has changed and today, the probability is quite low...

 base=read.table("millionb.csv",sep=";",header=TRUE)
X1=cumsum(base$nombre)
base=read.table("million2b.csv",sep=";",header=TRUE)
X2=cumsum(base$nombre)
X=X1+X2
 D=as.Date(as.character(base$date),"%m/%d/%Y")
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 1,7),method="CSS") 
 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]))
}
plot( D[length(D)-length(P)]+1:220,c(P,rep(NA,220-length(P))),
ylab="Probability to reach 1,000,000",xlab="",
type="l",col="red",ylim=c(0,1))