Freakonometrics

To content | To menu | To search

Tag - lissage

Entries feed - Comments feed

Monday, December 3 2012

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

Saturday, March 17 2012

Sondages et prévisions de séries temporelles

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

http://freakonometrics.blog.free.fr/public/perso5/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)

http://freakonometrics.blog.free.fr/public/perso5/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 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 !

Friday, February 24 2012

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.

http://freakonometrics.blog.free.fr/public/perso5/lex03.gif

http://freakonometrics.blog.free.fr/public/perso5/le04.gif désigne un poids attribué à la nouvelle observation dans la fonction de lissage http://freakonometrics.blog.free.fr/public/perso5/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")

http://freakonometrics.blog.free.fr/public/perso5/lissage-exp-1.gif

Sur la figure ci-dessus, non visualise l'impact de http://freakonometrics.blog.free.fr/public/perso5/lisse16.gif sur le lissage. La prévision que l'on fait à la date http://freakonometrics.blog.free.fr/public/perso5/lex14.gif, pour un horizon http://freakonometrics.blog.free.fr/public/perso5/lisse12.gif est alors http://freakonometrics.blog.free.fr/public/perso5/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

http://freakonometrics.blog.free.fr/public/perso5/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.

http://freakonometrics.blog.free.fr/public/perso5/lex50.gif

et

http://freakonometrics.blog.free.fr/public/perso5/HW02.gif

Dans ce cas, la prédiction est http://freakonometrics.blog.free.fr/public/perso5/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 http://freakonometrics.blog.free.fr/public/perso5/lisse15.gif (http://freakonometrics.blog.free.fr/public/perso5/lisse16.gif étant ici fixé),
http://freakonometrics.blog.free.fr/public/perso5/lissage-exp-2.gif
(le lissage simple - avec le même poids http://freakonometrics.blog.free.fr/public/perso5/lisse16.gif - apparaissant ici en trait clair).

Monday, February 20 2012

Lissage exponentiel et séries temporelles

L'idée du lissage exponentiel est de supposer que http://freakonometrics.blog.free.fr/public/perso5/lex01.gif, où on suppose que le bruit http://freakonometrics.blog.free.fr/public/perso5/lex02.gif est ici un bruit indépendant (imprévisible). On va alors supposer que

http://freakonometrics.blog.free.fr/public/perso5/lex03.gif

http://freakonometrics.blog.free.fr/public/perso5/le04.gif désigne un poids attribué à la nouvelle observation dans la fonction de lissage http://freakonometrics.blog.free.fr/public/perso5/lex05.gif.

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

http://freakonometrics.blog.free.fr/public/perso5/lex08.gif

i.e., en itérant

http://freakonometrics.blog.free.fr/public/perso5/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.

  • http://freakonometrics.blog.free.fr/public/perso5/lex05.gif sera interprété comme un niveau
  • http://freakonometrics.blog.free.fr/public/perso5/lex11.gif correspondant à la variation du processus
  • http://freakonometrics.blog.free.fr/public/perso5/lex12.gif qui correspond à un cycle (éventuel)

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

http://freakonometrics.blog.free.fr/public/perso5/lex50.gif

et

http://freakonometrics.blog.free.fr/public/perso5/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 http://freakonometrics.blog.free.fr/public/perso5/lex16.gif, à la date http://freakonometrics.blog.free.fr/public/perso5/lex14.gif, pour un horizon http://freakonometrics.blog.free.fr/public/perso5/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)

De la moyenne mobile aux notions de plus proches voisins (et de lissage)

En séries temporelles, il est naturelle de parler de moyennes mobiles. En fait, on peut adapter cette idée sur d'autres types de données, comme en économétrie (sur données individuelles). Ça permet assez facilement de construire un estimateur non-paramétrique (ou presque) de l'espérance conditionnelle. L'idée est relativement simple (sur le principe). La prévision consiste à approcher http://freakonometrics.blog.free.fr/public/perso5/lissage-03.gif. Dans le modèle linéaire, on supposer que l'on projeter sur l'espace des combinaisons linéaires, mais l'hypothèse était probablement un peu forte. Surtout si on cherche un ajustement local. Considérons alors http://freakonometrics.blog.free.fr/public/perso5/lissage-01.gif un voisinage de http://freakonometrics.blog.free.fr/public/perso5/lissage-08.gif, de telle sorte que

http://freakonometrics.blog.free.fr/public/perso5/lissage-02.gif

Notons

http://freakonometrics.blog.free.fr/public/perso5/lissage-04.gif

le nombre de points au voisinage de http://freakonometrics.blog.free.fr/public/perso5/lissage-08.gif, et

http://freakonometrics.blog.free.fr/public/perso5/lissage-05.gif

les indices des observations http://freakonometrics.blog.free.fr/public/perso5/lissage-12.gif au voisinage de http://freakonometrics.blog.free.fr/public/perso5/lissage-08.gif. On va alors considérer comme prévision

http://freakonometrics.blog.free.fr/public/perso5/lissage-09.gif

C'est ce que fait la fonction suivante,
> plot(cars)
> k=10
> X=cars$speed
> Y=cars$dist
> meank=function(x){
+ D=abs(x-X)
+ i=which(rank(D)<=k)
+ Yi=Y[i]
+ return(mean(Yi))
+ }
> MeanK=Vectorize(meank)
> u=seq(5,25,by=.2)
> v=MeanK(u)
> lines(u,v,col="red",type="s")
On peut aussi visualiser la construction de cet estimateur
> k=10
> dessin=function(x){
+ plot(cars)
+ u=seq(5,25,by=.2)
+ v=MeanK(u)
+ D=abs(x-X)
+ i=which(rank(D)<=k)
+ polygon(c(min(X[i])-.5,min(X[i])-.5,max(X[i])+
+ .5,max(X[i])+.5),c(0,120,120,0), + col="light blue",border=NA) + abline(v=x,col="red",lty=2) + points(X[i],Y[i],pch=19,col="blue") + points(x,mean(Y[i]),pch=19,col="red") + lines(u,v,col="red",type="s") + } > dessin(15)

http://freakonometrics.blog.free.fr/public/perso5/lissage-cars-1.gif
Le choix de la distance n'a pas d'impact pour la régression simple (car le choix de la distance ne changera pas l'ordre: le point le plus loin sera toujours le plus loin, quelle que soit la distance considérée). Mais si on régresse deux deux variables continues, le choix de la distance ne sera pas neutre. C'est - à peu de choses près - cette fonction qui est utilisée quand on a fait un lissage par moyenne mobile (afin d'enlever une composante saisonnière), en notant que les choses sont plus simples sur les séries temporelles, les observations étant régulièrement espacées,

http://freakonometrics.blog.free.fr/public/perso5/lissage-MA.gif

En fait, au lieu de prendre la moyenne sur ces points, on peut aussi prendre une moyenne pondérée: on mettra un poids http://freakonometrics.blog.free.fr/public/perso5/lissage-11.gif d'autant plus important que l'observation http://freakonometrics.blog.free.fr/public/perso5/lissage-12.gif est proche de http://freakonometrics.blog.free.fr/public/perso5/lissage-08.gif,

http://freakonometrics.blog.free.fr/public/perso5/lissage-10.gif

C'est l'estimateur proposé par Nadaraya-Watson (détaillé e.g. dans les notes de cours de Bruce Hansen).

> library(np)
> reg=npreg(dist~speed,data=cars)
> plot(reg,plot.errors.method="asymptotic",
+      plot.errors.style="band",col=c("red","black"))
> points(cars)

Une autre piste pourrait être de noter que faire une moyenne est probablement trop simple, et que faire un ajustement linéaire local pourrait donner de bien meilleurs résultats

> lissagek=function(x){
+ D=abs(x-X)
+ i=which(rank(D)<=k)
+ reg=lm(dist~speed,data=cars[i,])
+ Yi=predict(reg,newdata=data.frame(speed=x))
+ return(mean(Yi))
+ }

http://freakonometrics.blog.free.fr/public/perso5/lissage-cars-2b.gif

Mais on ne va pas ouvrir davantage la boite de Pandore des méthodes de lissage non-paramétriques, ce n'est pas (explicitement) l'objet du cours...

Régression et lissage par moyenne mobile

Notons http://freakonometrics.blog.free.fr/public/perso5/TS01.gif une série temporelle, observée jusqu'à la date http://freakonometrics.blog.free.fr/public/perso5/TS02.gif. 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)

http://freakonometrics.blog.free.fr/public/perso5/TS04.gif

> 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 http://freakonometrics.blog.free.fr/public/perso5/TS06.gif 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 (http://freakonometrics.blog.free.fr/public/perso5/TS07.gif=12), i.e.
http://freakonometrics.blog.free.fr/public/perso5/TS08.gif
Les coefficients http://freakonometrics.blog.free.fr/public/perso5/TS09.gif 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 http://freakonometrics.blog.free.fr/public/perso5/TS10.gif 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

http://freakonometrics.blog.free.fr/public/perso5/TS30.gif

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 http://freakonometrics.blog.free.fr/public/perso5/TS01.gif est d'effectuer un lissage par une moyenne mobile (ou glissante) dont la période est précisément http://freakonometrics.blog.free.fr/public/perso5/TS07.gif. Si http://freakonometrics.blog.free.fr/public/perso5/TS07.gif est pair (ce qui sera souvent le cas)
http://freakonometrics.blog.free.fr/public/perso5/TS13.gif
Si la série était parfaitement cyclique (http://freakonometrics.blog.free.fr/public/perso5/TS14.gif pour tout http://freakonometrics.blog.free.fr/public/perso5/TS16.gif) alors http://freakonometrics.blog.free.fr/public/perso5/TS17.gif 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. http://freakonometrics.blog.free.fr/public/perso5/TS32.gif
On utilise alors la méthode Buys-Balhot pour extraire la saisonnalité

http://freakonometrics.blog.free.fr/public/perso5/TS-18.gif

On a alors

http://freakonometrics.blog.free.fr/public/perso5/TS31.gif

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 http://freakonometrics.blog.free.fr/public/perso5/TS-20.gif,

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