Freakonometrics

To content | To menu | To search

Tag - SARIMA

Entries feed - Comments feed

Wednesday, December 5 2012

Syntaxe pour les SARIMA

La forme générale (enfin, la plus simple) des modèles SARIMA est la suivante

http://latex.codecogs.com/gif.latex?(1-L)^d(1-L^s)^{d_s}\Phi(L)\Phi_s(L^s)%20(X_t-m)%20=%20\Theta(L)\Theta_s(L^s)\varepsilon_t

Afin de comprendre comment les écrire sous R, commençons par un processus simple, autorégressive, de la forme
http://latex.codecogs.com/gif.latex?(1-\phi_1%20L)X_t%20=%20\varepsilon_t

La syntaxe est ici

> arima(X, order = c(p=1, d=0, q=0))

Supposons que l'on souhaite rajouter une composante moyenne mobile, i.e.

http://latex.codecogs.com/gif.latex?(1-\phi_1%20L)X_t%20=%20(1-\theta_1%20L)\varepsilon_t

La syntaxe devient

> arima(X, order = c(p=1, d=0, q=1))

Si on suppose maintenant que pour la composante autorégressive, on a une racine unité et le processus s'écrit

http://latex.codecogs.com/gif.latex?(1-%20L)X_t%20=%20(1-\theta_1%20L)\varepsilon_t

L’estimation de ce modèle se fait avec la commande

> arima(X, order = c(p=0, d=1, q=1))

Si pour finir, on veut calibrer un modèle de la forme

http://latex.codecogs.com/gif.latex?(1-%20L)(1-\phi_1%20L-\phi_2%20L^2)X_t%20=%20(1-\theta_1%20L)\varepsilon_t

on utilise la commande

> arima(X, order = c(p=2, d=1, q=1))

Supposons maintenant que notre série http://latex.codecogs.com/gif.latex?(X_t) ait été obtenu par différenciation saisonnière d'une série http://latex.codecogs.com/gif.latex?(Y_t), au sens où http://latex.codecogs.com/gif.latex?X_t=(1-L^{12})Y_t=\Delta^{12}%20Y_t, alors on devrait écrire

http://latex.codecogs.com/gif.latex?(1-%20L)(1-\phi_1%20L-\phi_2%20L^2)(1-L^{12})Y_t=%20(1-\theta_1%20L)\varepsilon_t

Pour estimer un tel modèle, la syntaxe est alors

> arima(X, order = c(p=2, d=1, q=1),
+ seasonal = list(order = c(0, 1, 0), period = 12))

En particulier, on peut avoir deux modèles assez proches pour modéliser une série cyclique: un bruit blanc avec une intégration saisonnière, i.e.

http://latex.codecogs.com/gif.latex?(1-L^{12})%20X_t%20=%20\varepsilon_t

dont la syntaxe serait

> arima(X, order = c(p=0, d=0, q=0),
+ seasonal = list(order = c(0, 1, 0), period = 12))

et un simple processus autorégressif à l'ordre 12. Mais là encore, soit je veux imposer que seule la dernière composante soit non nulle i.e.

http://latex.codecogs.com/gif.latex?(1-\phi_{12}%20L^{12})%20X_t%20=%20\varepsilon_t

ce qui s'écrit

> arima(X, order = c(p=0, d=0, q=0),
+ seasonal = list(order = c(1, 0, 0), period = 12))

soit je prends un processus plus général, de la forme

http://latex.codecogs.com/gif.latex?(1-\phi_1%20L-\cdots%20\phi_{12}%20L^{12})%20X_t%20=%20\varepsilon_t

dont la syntaxe devient

> arima(X, order = c(p=12, d=0, q=0))

En fait, on peut introduire un polynôme autorégressive avec uniquement des retards à l'ordre 12, en plus de la différentiation saisonnière, i.e.

http://latex.codecogs.com/gif.latex?(1-L)(1-L^{12})(1-\phi_1%20L)(1-\phi_{12}%20L^{12})%20X_t%20=%20(1-\theta_1%20L)\varepsilon_t

Ce modèle s'écrire alors sous la forme

> arima(X, order = c(p=1, d=1, q=1),
+ seasonal = list(order = c(1, 1, 0), period = 12))

On peut bien sûr aller plus loin, en autorisant non seulement une composante autorégressive saisonnière, mais pourquoi pas, une composante moyenne mobile saisonnière, i.e.

http://latex.codecogs.com/gif.latex?(1-L)(1-L^{12})(1-\phi_1%20L)(1-\phi_{12}%20L^{12})%20X_t%20=%20(1-\theta_1%20L)\varepsilon_t(1-\theta_{12}L^{12})\varepsilon_t

qui s'écrit, sous R,

> arima(X, order = c(p=1, d=1, q=1),
+ seasonal = list(order = c(1, 1, 1), period = 12))

On a vu la forme la plus générale des SARIMA. Enfin, comme je le disais en préambule, c'est aussi la plus simple. Car si je suppose qu'une autre cycle se superpose au cycle annuel, je peux le faire. Mais on va peut être essayer d'éviter de trop compliquer les notations...

Friday, November 30 2012

Estimation et prévision pour des séries temporelles

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.

Wednesday, November 21 2012

Introduction aux modèles SARIMA

Quelques transparents en plus, qui devraient correspondre aux deux prochains cours de séries temporelles, sur les processus autorégressifs (AR) et moyennes mobiles (MA), les ARMA, les ARIMA (intégrés) et les SARIMA (saisonniers). J'ai mis des notes sur les tests de racine unités, je rajouterais quelques transparents la semaine prochaine sur les tests de saisonnalité, et quelques exemples pratiques de prévision. Les transparents sont en ligne ici,

Wednesday, September 19 2012

Prévision, choix de modèle(s) et backtesting

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 http://latex.codecogs.com/gif.latex?X_t-\overline%20X

> A7d=A7-mean(A7)
> plot(A7d,xlim=c(1989,2000))

On avait vu qu'un modèle http://latex.codecogs.com/gif.latex?AR(12) 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 http://latex.codecogs.com/gif.latex?{}_T\widehat{X}_{T+h}, i.e. http://latex.codecogs.com/gif.latex?EL(X_{T+h}|\underline{X}_T), où http://latex.codecogs.com/gif.latex?EL(\cdot) 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 http://latex.codecogs.com/gif.latex?{}_T\text{mse}_{T+h}=\mathbb{E}([{}_T\widehat{X}_{T+h}-X_{T+h}]^2|\underline{X}_T),

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