Freakonometrics

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 !

L'idée du lissage exponentiel est de supposer que , où on suppose que le bruit est ici un bruit indépendant (imprévisible). On va alors supposer que

où désigne un poids attribué à la nouvelle observation dans la fonction de lissage .

On parle de lissage pour moyenne mobile à poids exponentiel (ewma, exponentially weighted moving average, évoqué ici ou là sur ce blog). En effet

i.e., en itérant

(même s'il faut faire attention à cette décomposition qui, en pratique, s’arrête avec la première observation).
La méthode de Holt-Winters propose de généraliser ce lissage, en introduisant deux (voire trois) composantes.

• sera interprété comme un niveau
• correspondant à la variation du processus
• qui correspond à un cycle (éventuel)

On suppose désormais, dans cette méthode (dite de Holt-Winters), que

et

D'un point de vue plus pratique, la programmation se fait comme suit,

> HW=HoltWinters(X,alpha=.2,beta=0)
> plot(HW)

> P=predict(HW,24,prediction.interval=TRUE)
> plot(HW,xlim=range(c(time(X),time(P))))
> polygon(c(time(P),rev(time(P))),c(P[,2],rev(P[,3])),
+ col="yellow",border=NA)
> lines(P[,1],col="red",lwd=3)

Notons une série temporelle, observée jusqu'à la date . La première peut être de penser que l'on a une tendance linéaire. On peut alors chercher à résoudre le problème (classique)

> autoroute=read.table(
+ "http://freakonometrics.blog.free.fr/public/data/autoroute.csv",
+ header=TRUE,sep=";")
> a7=autoroute$a007
> X=ts(a7,start = c(1989, 9), frequency = 12)
> T=time(X)
> S=cycle(X)
> B=data.frame(x=as.vector(X),T=as.vector(T),S=as.vector(S))
> regT=lm(x~T,data=B)
> plot(X)
> abline(regT,col="red")

On posera alors la série dont la tendance linéaire a été enlevée. On peut aussi visualiser la tendance sous cette forme

> X1=predict(regT)
> B$X1=X1
> Y=B$X1
> plot(X,xlab="",ylab="")
> YU=apply(cbind(X,Y),1,max)
> YL=apply(cbind(X,Y),1,min)
> i=which(is.na(Y)==FALSE)
> polygon(c(T[i],rev(T[i])),c(Y[i],rev(YU[i])),col="blue")
> polygon(c(T[i],rev(T[i])),c(Y[i],rev(YL[i])),col="red")
> lines(B$T,Y,lwd=3)

On peut alors étudier la composante cyclique, par la méthode de Buys-Balhot: si la série est mensuelle, on va chercher un cycle mensuel (=12), i.e.

Les coefficients sont estimés par moindres carrés. Cette écriture un peu prétentieuse veut juste dire que l'on va régresser la série sur le mois (pris en tant que facteur)

> B$res1=X-X1
> regS=lm(res1~0+as.factor(S),data=B)
> B$X2=predict(regS)
> plot(B$S,B$res1,xlab="saisonnalité")
> lines(B$S[1:4],B$X2[1:4],col="red",lwd=2)
> lines(B$S[5:13],B$X2[5:13],col="red",lwd=2)

(on a coupé la tendance en deux - sur le graphique - car les données ne commençaient pas en janvier mais en septembre). On a alors

où le terme bleu est la tendance (ici linéaire), le terme rouge est le cycle, et le terme vert est la composante (supposée) aléatoire.

> plot(X)
> lines(B$T,B$X1+B$X2,col="red")

ou encore

> plot(X,xlab="",ylab="")
> YU=apply(cbind(X,Y),1,max)
> YL=apply(cbind(X,Y),1,min)
> i=which(is.na(Y)==FALSE)
> polygon(c(T[i],rev(T[i])),c(Y[i],rev(YU[i])),col="blue")
> polygon(c(T[i],rev(T[i])),c(Y[i],rev(YL[i])),col="red")
> lines(B$T,Y,lwd=3)

On notera qu'il est très facile de faire de la prévision à l'aide de ce modèle: il suffit de prolonger la tendance (linéaire), de rajouter la composante cyclique, et l'aléa sera la variance du bruit résiduel,

> n=length(T)
> T0=T[(n-23):n]+2
> plot(X,xlim=range(c(T,T0)))
> X1p=predict(regT,newdata=data.frame(T=T0))
> Yp=X1p+B$X2[(n-23):n]
> se=sd(X-Y)
> YpU=Yp+1.96*se
> YpL=Yp-1.96*se
> polygon(c(T0,rev(T0)),c(YpU,rev(YpL)),col="yellow",border=NA)
> lines(T0,Yp,col="red",lwd=3)
> lines(B$T,Y,lwd=3,col="red")

avec ici un prévision sur un horizon de 24 mois,

Une autre manière de prendre en compte la saisonnalité de la série est d'effectuer un lissage par une moyenne mobile (ou glissante) dont la période est précisément . Si est pair (ce qui sera souvent le cas)

Si la série était parfaitement cyclique ( pour tout ) alors serait plate. Cette série correspond à un lissage sur une année glissante.

> MA=function(x){(sum(x[2:12])+(x[1]+x[13])/2)/12}
> Y=rep(NA,length(X))
> for(t in 7:(length(X)-7)){
+ Y[t]=MA(X[(t-6):(t+6)])
+ }
> plot(X,xlab="",ylab="")
> YU=apply(cbind(X,Y),1,max)
> YL=apply(cbind(X,Y),1,min)
> i=which(is.na(Y)==FALSE)
> polygon(c(T[i],rev(T[i])),c(Y[i],rev(YU[i])),col="blue")
> polygon(c(T[i],rev(T[i])),c(Y[i],rev(YL[i])),col="red")
> lines(B$T,Y,lwd=3)

On peut alors reprendre la modélisation précédente, en remplaçant la tendance linéaire par cette série lissée, i.e.
On utilise alors la méthode Buys-Balhot pour extraire la saisonnalité

On a alors

où le terme bleu est la tendance (qui n'est plus - parfaitement - linéaire), le terme rouge est le cycle, et le terme vert est la composante (supposée) aléatoire (sous en tous les cas imprévisible... mais on passera du temps par la suite sur ce dernier terme).

> B$res2=X-Y
> regS2=lm(res2~0+as.factor(S),data=B)
> B$X3=NA
> B$X3[i]=predict(regS2)
> plot(X,xlab="",ylab="")
> YU=apply(cbind(X,Y+B$X3),1,max)
> YL=apply(cbind(X,Y+B$X3),1,min)
> i=which(is.na(Y)==FALSE)
> polygon(c(T[i],rev(T[i])),c(Y[i]+B$X3[i],rev(YU[i])),col="blue")
> polygon(c(T[i],rev(T[i])),c(Y[i]+B$X3[i],rev(YL[i])),col="red")
> lines(B$T,Y+B$X3,lwd=3)

Cette décomposition, aussi simpliste qu'elle paraisse, est celle qui est obtenue sous R

> bruit=X-(Y+B$X3)
> D=decompose(X)
> plot(D,col="blue")

On notera que les résidus obtenus ici coïncident avec ,

> plot(bruit)
> lines(D$random,col="red")

Il est d'ailleurs possible d'utiliser une décomposition non plus additive (comme on vient de le faire) mais multiplicative...