Freakonometrics

Pour la fin de cours , quelques transparents sur l'identification et l'estimation des modèles SARIMA, quelques compléments sur les tests (racines unités et non-stationnarité, ainsi que saisonnalité), et enfin, quelques pistes pour construire des prédictions (avec une quantification de l'incertitude). Les transparents sont en ligne ici. Maintenant que les transparents sont finis (et en ligne) les prochains billets seront orientés autour de la modélisation et des aspects computationnels.

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,