Tag Archives: lissage

Décomposition(s) de Séries Temporelles

La semaine dernière, pour le TP du cours de série temporelle, on avait parlé de décomposition de séries temporelles. On part d’une série temporelle , et on va considérer une décomposition de la forme

par exemple une tendance linéaire, et un cycle annuel,

Considérons la série suivante

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

correspondant à du trafic sur l’autoroute (A7)

>  plot(X,xlim=c(1989,2000))

Pour estimer la tendance, on va faire une régression linéaire, sur le temps, à partir de la base suivante

>  T=time(X)
>  base=data.frame(Y=a7,X1=as.numeric(T))

Le modèle est ici

>  reg1=lm(Y~X1,data=base)

ce qui permet de faire une prévision sur les données passées, mais aussi futures,

>  tendance1=predict(reg1)
>  T=as.numeric(T)
>  lines(T,tendance1,col="red",lwd=3)
>  futur=seq(max(T)+1/12,length=36,by=1/12)
>  tendance2=predict(reg1,
+  newdata=data.frame(X1=futur))
>  lines(futur,tendance2,col="red")

Voilà pour la tendance. Pour le cycle, on va ajuste un modèle de la forme suivante

sur la série des résidus

>  base$Y2=residuals(reg1)

Cette fois, on va définir le mois de la date

> base$mois=round(base$X1-t
+ trunc(base$X1),digit=4)

Et on régresse le résidu sur le mois

>  reg2=lm(Y2~0+as.factor(mois),data=base)

Là encore, on peut faire une régression sur les données passées

>  cycle1=predict(reg2)
>  lines(T,tendance1+cycle1,col="blue",lwd=3)

mais aussi sur les données futures,

>  moisfutur=round(futur-trunc(futur),4)
>  cycle2=predict(reg2,newdata=
+  data.frame(mois=moisfutur))
>  lines(futur,tendance2+cycle2,col="blue")

Cela dit, pour les données futures, on peut aussi tenter de rajouter une région de confiance, en tenant compte de l’erreur d’ajustement,

>  E=residuals(reg2)
>  hist(E)
>  1.96*sqrt(var(E))
[1] 5637.897

On peut alors facilement ajouter notre région de confiance avec

> polygon(c(futur,rev(futur)),
+ c(tendance2+cycle2-1.96*sqrt(var(E)),
+ rev(tendance2+cycle2+1.96*sqrt(var(E)))),
+ col="yellow",border=NA)

Et voilà…. On peut quand même vérifier que le bruit est effectivement non-autocorrélé, et c’est le cas ici

>  plot(acf(E))

Comme on l’a vu en TP, ajuster une droite de régression c’est bien ici, mais peut être qu’on peut envisager un ajustement local. Par exemple prendre la moyenne autour de la date qui nous intéresse, +/- 6 mois (techniquement, si on rajoute 6 mois avant et après, on aura un soucis, alors on fait une petite correction au bord),

ce qui correspond à un lissage par moyenne mobile.

>  n=length(X)
>  F=rep(NA,n)
>  for(t in 7:(n-7)) F[t]=1/12*(sum(X[(t-5):(t+5)])+ (X[t-6]+X[t+6])/2)
>  lines(T,F,lwd=3,col="purple")

On voit ici – par notre ajustement non-paramétrique – que l’ajustement d’une tendance linéaire n’était pas si mauvais que ça

Par contre, pour des séries où la tendance n’est clairement pas linéaire, le lissage par moyenne mobile permet d’avoir une bonne estimation de la tendance de fond

> sncf=read.table(
+ "http://freakonometrics.free.fr/sncf.csv",
+ header=TRUE,sep=";")
> train=as.vector(t(as.matrix(sncf[,2:13])))
> X=ts(train,start = c(1963, 1), frequency = 12)
> plot(X)  
> n=length(X)
> F=rep(NA,n)
> for(t in 7:(n-7)) F[t]=1/12*(sum(X[(t-5):(t+5)])+(X[t-6]+X[t+6])/2)
> lines(as.numeric(time(X)),
+ F,lwd=3,col="purple")

On continuera en cours cette semaine autour du lissage, par moyenne mobile, et avec du lissage exponentiel ensuite. A suivre donc.

Elections Sans Bipartisme

Hier, sur Twitter, @JF_Godbout partageait un joli graphique relatif aux élections québécoises, avec les nombres de votes obtenus (ici en pourcentage des votes totaux) et le pourcentage de sièges que cela donne,

Il faut dire qu’hier, c’était jour d’élection. Ce sont des élections à un tour, avec plusieurs partis (disons 4 grands partis si on se limite aux élections les plus récentes). Jusqu’à lundi, le PQ (Parti Québécois) avait une majorité (relative), et depuis lundi, le PLQ (Parti Libéral) est revenu au pouvoir. Et il y a à côté de ces deux grands partis d’autres partis, que l’on va regrouper. @JF_Godbout m’ayant envoyé ses données, j’ai pu reproduire les graphiques.

La grande différence ici, entre nos approches est – je pense – la méthode de lissage: alors que  @JF_Godbout faisait des régressions locales, j’ai opté pour des splines.

reg=glm(Seat~bs(Election),data=d,family=binomial)
yp=predict(reg,newdata=data.frame(Election=1970:2015),type="response")

J’utilise des splines parce que ça se généralise facilement à d’autres lois, et qu’on peut aussi les utiliser quand on régresse sur plusieurs variables continues. Si on rajoute les résultats de l’élection de lundi, on obtient le graphique suivant

Ce qui me gêne dans cette modélisation, c’est qu’on considère nos trois dernières régressions comme indépendantes…. alors que ce n’est pas le cas: on a trois modalités possibles, et la somme doit forcément faire 1. D’ailleurs, quand on somme les prédictions, on dépasse 1.

Avec deux partis, notre régression logistique nous garantie d’avoir une somme à un, mais plus ici. A moins d’utiliser un modèle multinomial. Si on oublie le lissage, l’idée est d’écrire

qui est tout simplement l’extension du modèle logistique. La difficulté technique est que, sous R, quand on fait une régression logistique, la variable d’intérêt doit être une variable prenant trois facteurs, et non pas un vecteur de probabilités. Au lieu d’utiliser des taux, ou des pourcentages, je devrais utiliser des nombres d’observations de chaque modalités. On va donc crééer une fausse base

> vote=data.frame(
+ plq=plq$Vote,
+ pq=pq$Vote,
+ autre=autre$Vote,
+ Election=pq$Election)
> vote
      plq     pq  autre Election
1  0.5465 0.3022 0.1513     1972
2  0.3378 0.4137 0.2485     1976
3  0.4608 0.4926 0.0466     1981
4  0.5599 0.3869 0.0532     1985
5  0.4995 0.4016 0.0989     1989
6  0.4440 0.4475 0.1085     1994
7  0.4355 0.4287 0.1358     1998
8  0.4599 0.3324 0.2077     2003
9  0.3308 0.2835 0.3857     2007
10 0.4208 0.3517 0.2275     2008
11 0.3120 0.3195 0.3685     2012
12 0.4140 0.2538 0.3312     2014
> M=t(round(vote[,1:3]*1000))
> M
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
plq    546  338  461  560  500  444  436  460  331   421   312   414
pq     302  414  493  387  402  448  429  332  284   352   320   254
autre  151  249   47   53   99  108  136  208  386   227   368   331

On réplique alors autant de fois le nom du parti qu’on a de (pseudo) votes

> Y=X=NULL
> for(i in 1:3){
+ for(j in 1:ncol(M)){
+ Y=c(Y,rep(c("PLQ","PQ","Autre")[i],M[i,j]))
+ X=c(X,rep(pq$Election[j],M[i,j]))
+ }}
> basev=data.frame(Y=as.factor(Y),X=X)

C’est un peu artificiel, j’en conviens, mais ça marche… On peut alors faire la régression

> library(nnet)
> regv=multinom(Y~bs(X),data=basev)
# weights:  15 (8 variable)
initial  value 13186.643301 
iter  10 value 12241.931805
final  value 12221.591172 
converged

et des prévisions

> ypv=predict(regv,newdata=
+ data.frame(X=1973:2014),type="probs")
> ypv
        Autre       PLQ        PQ
1  0.18426281 0.4705415 0.3451957
2  0.16051060 0.4767408 0.3627486
3  0.14149153 0.4805514 0.3779571
4  0.12631826 0.4826100 0.3910718
5  0.11426236 0.4834165 0.4023212
6  0.10474046 0.4833543 0.4119053
7  0.09729304 0.4827117 0.4199953
8  0.09156246 0.4817016 0.4267360
9  0.08727323 0.4804776 0.4322491
10 0.08421535 0.4791477 0.4366370
11 0.08223095 0.4777835 0.4399855
12 0.08120364 0.4764292 0.4423672
13 0.08105030 0.4751066 0.4438431
14 0.08171457 0.4738203 0.4444651

Et on notera que cette fois, nos prévisions somment à un.

Cela dit, visuellement, on obtient des choses proches de ce qu’on avait auparavant,

mais intellectuellement, je trouve plus logique d’utiliser autre chose que des régressions logistiques quand on est dans une élection où plusieurs partis sont en compétition.

Introduction aux modèles linéaires généralisés

J’ai un peu d’avance dans le cours. Je vais mettre en ligne les transparents pour la semaine prochaine (normalement), où nous aborderons la classe des modèles linéaires généralisés. Les transparents sont en ligne ici.

Je n’ai pas mis de section sur lesGeneralized Additive Models, on se contentera de la section sur le lissage évoquée à la fin des transparents sur la modélisation de la fréquence. Afin de légitimer les méthodes de lissage (sur l’âge de l’assuré en particulier), je renvoie vers un graphique produit il y a plusieurs années par un cabinet de conseil, qui notait que la forme de la fonction de lissage, liant l’âge à la fréquence est identique, dans tous les pays,

https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2013/02/assurance4.jpgMais je pense que je ferais un billet dédié au lissage, dans la problématique de la tarification en assurance IARD.

Méthodes de lissage en assurance

Dans le cours de modèles de prévision, j’avais abordé les méthodes de régression locale rapidement, en finissant la section sur les données individuelles. On verra plus d’outils la session prochaine dans le cours actuariat IARD, a.k.a. ACT2040 (que je donnerai cet hiver). En attendant, je mets en ligne les transparents du cours que vient de donner Julien Tomas à l’Institut de science financière et d’assurances (à Lyon, en France), dans un contexte d’assurance dépendance.

(et que Julien m’a autorisé à diffuser).

Régression: données individuelles ou temporelles ?

Avant de parler davantage du premier devoir, un mot pour expliquer pourquoi dans la première partie il est important d’éviter d’utiliser des données temporelles (le cours présente des outils pour modéliser des données individuelles). Considérons un exemple simple (évoqué voilà plusieurs mois ici même), où on cherche à prédire un tour de poitrine conditionnellement à un tour de taille (pour une personne de sexe féminin, cf. le mannequin ci-contre). L’idée intuitive est qu’il doit exister une relation positive, possiblement linéaire entre les deux mensurations. A partir de données, on devrait pouvoir quantifier cette relation. On dispose pour cela de données, en ligne sur

> mensurations<-read.table("
+ http://freakonometrics.free.fr/mensurations.csv",
+ header=TRUE,sep=";")
> head(mensurations[,-1])
POITRINE TAILLE HANCHE HAUTEUR POIDS      AGE
1       NA     NA     NA      NA    NA       NA
2       94     58     89     165    51 22.24230
3       NA     NA     NA      NA    NA       NA
4       94     58     89     165    51 22.40383
5       91     58     97     165    55       NA
6       NA     NA     NA      NA    NA       NA

Si on regarde rapidement la corrélation entre la taille et la poitrine, cette dernière est faible,

> Z=mensurations[,c("TAILLE","POITRINE")]
> I=apply(Z,1,function(x) sum(is.na(x)))
> cor(Z[I==0,])
TAILLE   POITRINE
TAILLE   1.00000000 0.09700002
POITRINE 0.09700002 1.00000000

… mais significativement non nulle,

> cor.test(Z[,1],Z[,2],,method = "pearson")

Pearsons product-moment correlation

data:  Z[, 1] and Z[, 2]
t = 2.4559, df = 635, p-value = 0.01432
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.0194634 0.1733769
sample estimates:
cor
0.09700002

Avec une https://latex.codecogs.com/gif.latex?p-value de l’ordre de 1.5%, on va rejeter l’hypothèse https://latex.codecogs.com/gif.latex?H_0 de nullité de la corrélation (c’est ici un test basé sur la transformation de Fisher). On peut aussi faire une régression simple,

> reg=lm(POITRINE~TAILLE,data=mensurations)
> summary(reg)

Call:
lm(formula = POITRINE ~ TAILLE, data = mensurations)

Residuals:
Min      1Q  Median      3Q     Max
-9.7253 -3.4424 -0.2091  1.5576 13.8579

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 82.67798    2.82953  29.220   <2e-16 ***
TAILLE       0.11663    0.04749   2.456   0.0143 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.133 on 635 degrees of freedom
(22 observations deleted due to missingness)
Multiple R-squared: 0.009409,	Adjusted R-squared: 0.007849
F-statistic: 6.031 on 1 and 635 DF,  p-value: 0.01432

(on notera que le test sur la corrélation, le test de significativité de la – seule – variable explicative, et le test de significativité globale donnent ici exactement la même https://latex.codecogs.com/gif.latex?p-value de l’ordre de 1.5%). Là encore, on voit que le tour de taille a un effet significativement non nul… mais marginal. Graphiquement, la prédiction avec le tour de taille est le trait rouge,

La valeur moyenne est en trait pointillé. On peut voir que la différence est infime,

> base=data.frame(TAILLE=seq(50,70,by=5))
> base$Pred.POITRINE=predict(reg,newdata=base)
> base$Moy.POITRINE=mean(POITRINE,na.rm=TRUE)
> base
TAILLE Pred.POITRINE Moy.POITRINE
1     50      88.50935     89.61538
2     55      89.09248     89.61538
3     60      89.67562     89.61538
4     65      90.25876     89.61538
5     70      90.84189     89.61538

Le soucis est que cette base est en fait une base sur données individuelles. En fait, il s’agit des mensurations des demoiselles qui ont posé comme playmate dans la version U.S. du magasine Playboy. Si on se restreint sur des sous-période, par exemple les années 70, on a des prédictions relativement différentes,

> submensurations=mensurations[abs(DATE-1975)<=5,]
> subreg=lm(POITRINE~TAILLE,data=submensurations)
> summary(subreg)

Call:
lm(formula = POITRINE ~ TAILLE, data = submensurations)

Residuals:
Min     1Q Median     3Q    Max
-8.655 -1.871 -0.655  1.750 13.155

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  66.1485     7.3839   8.958 6.12e-15 ***
TAILLE        0.4053     0.1230   3.296   0.0013 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.574 on 117 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.08496,	Adjusted R-squared: 0.07714
F-statistic: 10.86 on 1 and 117 DF,  p-value: 0.001299

> base$Pred70s.POITRINE=predict(subreg,newdata=base)
> base
TAILLE Pred.POITRINE Moy.POITRINE Pred70s.POITRINE
1     50      88.50935     89.61538         86.41274
2     55      89.09248     89.61538         88.43917
3     60      89.67562     89.61538         90.46559
4     65      90.25876     89.61538         92.49201
5     70      90.84189     89.61538         94.51844

que l’on peut visualiser ci-dessous,

L’explication est que la morphologie des femmes (ou en tous les cas de celles qui ont posé nues dans Playboy) a beaucoup changé en 60 ans, comme évoquéprécédemment. On peut le voir en regardant l’évolution de la pente, par date (https://latex.codecogs.com/gif.latex?\pm5ans)

(le trait horizontal est la pente obtenue sur l’ensemble de la période). On peut d’ailleurs intégrer l’année dans notre modèle, en rajoutant le temps comme variable explicative dans la régression,

> reg=lm(POITRINE~TAILLE+DATE,data=mensurations)
> summary(reg)

Call:
lm(formula = POITRINE ~ TAILLE + DATE, data = mensurations)

Residuals:
Min      1Q  Median      3Q     Max
-9.7304 -2.2184 -0.5405  2.0742 12.1041

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 368.296266  18.008455  20.451  < 2e-16 ***
TAILLE        0.323319   0.042140   7.673 6.39e-14 ***
DATE         -0.150301   0.009393 -16.002  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.491 on 634 degrees of freedom
(22 observations deleted due to missingness)
Multiple R-squared: 0.2944,	Adjusted R-squared: 0.2922
F-statistic: 132.3 on 2 and 634 DF,  p-value: < 2.2e-16

Visuellement, on aurait la prédiction suivante,

Si on se place autour de 1970,

> submensurations=mensurations[abs(DATE-1970)<=5,]
> subreg=lm(POITRINE~TAILLE,data=submensurations)
> predict(subreg,newdata=data.frame(TAILLE=c(55,70)))
1        2
89.14361 95.75476

alors qu’autour de 1990,

> submensurations=mensurations[abs(DATE-1990)<=5,]
> subreg=lm(POITRINE~TAILLE,data=submensurations)
> predict(subreg,newdata=data.frame(TAILLE=c(55,70)))
1        2
87.69853 92.70228

Notons qu’on fait ici des régressions locales (dans ces deux derniers calculs), ce qui n’est pas ce qui est représenté visuellement… Si on veut une meilleure modélisation locale, il faudrait quitter le monde des modèles linéaires pour aller un peu plus loin. On verra en fin de cours si on a le temps d’en parler

Cela dit, globalement, on avait un effet positif, mais il aurait été tout à fait possible d’avoir un effet (global négatif). C’est ce qui arrive dans le paradoxe de Simpson.

Sondages et prévisions de séries temporelles

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),

https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic10.gif

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)

https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/ic-sondage01.gif

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 unpré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 !

Sur le lissage exponentiel

Sur le lissage exponentiel, je pourrais renvoyer vers Gardner (1985), et la version remise au goût du jour, Gardner (2005). J’avais pris comme position, dans le cours, de présenter rapidement le lissage exponentiel, en notant qu’on les évoquerais à nouveau plus tard comme cas particulier des modèles ARIMA. Comme le note la dernière version, Gardner (2005), ce n’est pas forcément l’unique manière de voir, et le cours pourrait être orienté complètement autour des notions de lissage exponentiel (c’est d’ailleurs le point du vue de livre Hyndman, Koehler, Ord & Snyder (2008)), “when Gardner (2005) appeared, many believed that exponential smoothing should be disregarded because it was either a special case of ARIMA modeling or an ad hoc procedure with no statistical rationale. As McKenzie (1985) observed, this opinion was expressed in numerous references to my paper. Since 1985, the special case argument has been turned on its head, and today we know that exponential smoothing methods are optimal for a very general class of state-space models that is in fact broader than the ARIMA class.”

N’étant toujours pas convaincu, j’évoquerais à nouveau le lissage exponentiel quand on présentera les modèles ARIMA. En attendant, un peu de code pour mieux comprendre ce qui est fait quand on fait du lissage exponentiel.

Commençons par le lissage exponentiel simple, i.e.

https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex03.gif

https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/le04.gif désigne un poids attribué à la nouvelle observation dans la fonction de lissage https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex05.gif. Le code pour lisser une série est alors assez simple,

> library(datasets)
> X=as.numeric(Nile)
> Lissage=function(a){
+  T=length(X)
+  L=rep(NA,T)
+  L[1]=X[1]
+  for(t in 2:T){L[t]=a*X[t]+(1-a)*L[t-1]}
+  return(L)
+ }
> plot(X,type="b",cex=.6)
> lines(Lissage(.2),col="red")

https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lissage-exp-1.gif

Sur la figure ci-dessus, non visualise l’impact de https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse16.gif sur le lissage. La prévision que l’on fait à la date https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex14.gif, pour un horizon https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse12.gif est alors https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse11.gif. Il est alors possible de voir le poids comme un paramètre et on va alors chercher le poids optimal. La stratégie classique est de minimiser l’erreur de prédiction commise à un horizon de 1

https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse14.gif

> V=function(a){
+  T=length(X)
+  L=erreur=rep(NA,T)
+  erreur[1]=0
+  L[1]=X[1]
+  for(t in 2:T){
+  L[t]=a*X[t]+(1-a)*L[t-1]
+  erreur[t]=X[t]-L[t-1] }
+ return(sum(erreur^2))
+ }

Ici, on obtient comme poids optimal

> A=seq(0,1,by=.02)
> Ax=Vectorize(V)(A)
> plot(A,Ax,ylim=c(min(Ax),min(Ax)*1.05))
> optimize(V,c(0,.5))$minimum
[1] 0.246581

On notera que c’est ce que suggère la fonction de R,

> hw=HoltWinters(X,beta=FALSE,gamma=FALSE,l.start=X[1])
> hw
Holt-Winters exponential smoothing without trend an seasonal comp.

Call:
HoltWinters(x = X, beta = FALSE, gamma = FALSE, l.start = X[1])

Smoothing parameters:
alpha:  0.2465579
beta :  FALSE
gamma:  FALSE

Coefficients:
[,1]
a 805.0389

> plot(hw)
> points(2:(length(X)+1),Vectorize(Lissage)(.2465),col="blue")

Dans le cas du lissage exponentiel double, i.e.

https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex50.gif

et

https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/HW02.gif

Dans ce cas, la prédiction est //f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse13.gif" alt="https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse13.gif" />. Le code pour faire un lissage double est là encore assez simple,

> Lissage=function(a,b){
+  T=length(X)
+  L=B=rep(NA,T)
+  L[1]=X[1]; B[1]=0
+  for(t in 2:T){
+  L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1])
+  B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] }
+ return(L)
+ }

Sur la figure suivante, on visualise l’évolution du lissage en fonction de https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse15.gif (https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse16.gif étant ici fixé),
https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lissage-exp-2.gif
(le lissage simple – avec le même poids https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse16.gif – apparaissant ici en trait clair).

Lissage exponentiel et séries temporelles

L’idée du lissage exponentiel est de supposer que https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex01.gif, où on suppose que le bruit https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex02.gif est ici un bruit indépendant (imprévisible). On va alors supposer que

https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex03.gif

où https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/le04.gif désigne un poids attribué à la nouvelle observation dans la fonction de lissage https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex05.gif.

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

https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex08.gif

i.e., en itérant

https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex09.gif

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

  • https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex05.gif sera interprété comme un niveau
  • https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex11.gif correspondant à la variation du processus
  • https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex12.gif qui correspond à un cycle (éventuel)

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

https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex50.gif

et

https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/HW02.gif

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

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

Cette trois composantes seront utilisées pour faire une prévision de la série: la prédiction sera de la forme https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex16.gif, à la date https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex14.gif, pour un horizon https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex17.gif (avec les notations du précédant billet).

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