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.