La chronique airmiles (dans la bibliothèque TSA), donne le nombre de passagers transportés sur les lignes aériennes aux USA de janvier 1996 à mai 2007.
Proposez des modèles d’intervention basés sur :
library(forecast)
library(lmtest)
library(TSA)
options(repr.plot.width=12, repr.plot.height=10) # pour redimensionner les graphiques
data(airmiles)
airmiles
On remarque que les données sont disponibles jusque mai 2005.
plot(airmiles, type="o", pch =20)
grid()
plot(log(airmiles), type="o", pch =20)
grid()
On identifie plutôt un modèle multiplicatif: l'amplitude des fluctuations périodiques étant (un peu) plus grande lorsque la tendance est plus élevée (cf séance 1).
Lorsqu'on examine le logarithme de la chronique, les fluctuations périodiques semblent plus régulières. On travaillera donc sur le logarithme de la chronique. Notons qu'il s'agit d'un choix de modélisation discutable qui pourrait être remis en question.
Etant donnée la saisonnalité, nous allons opérer une dérivée saisonnière.
plot(diff(log(airmiles),12))
grid()
L'aspect du graphe est plutôt celui d'une "marche aléatoire" que celui d'un processus stationnaire.
Examinons le graphe de cette chronique encore dérivée.
plot(diff(diff(log(airmiles),12)))
Le graphe a maintenant l'air de celui d'un processus stationnaire, "à l'intervention près".
(sans intervention, pour avoir une référence à laquelle comparer)
On travaillera sur log(airmiles)
.
tsdisplay(log(airmiles))
L'ACF présente une périodicité: cela confirme la nécessité d'une dérivation saisonnière.
tsdisplay(diff(log(airmiles),12))
L'ACF de la chronique dérivée saisonnièrement présente une décroissance lente: elle n'est donc toujours pas stationnaire, on dérive à nouveau la chronique.
tsdisplay(diff(diff(log(airmiles),12)))
On regarde les corrélation "à court-terme" (décalages < 6). L'explication la plus simple pourrait être: décroissance dans le PACF et décalage 1 significatif dans l'ACF, ce qui indiquerait un processus MA(1) à court-terme. C'est bien entendu assez discutable (nous avons examiné l'hypothèse AR(1) dans la correction de cette chronique lors de la séance 4). Ne perdons pas de vue que la modélisation sera nécessairement imparfaite sans intervention, puisque la modélisation de la chronique nécessite, justement, une intervention.
Nous allons donc examiner les résidus d'un tel modèle. Rappelons que dans les arguments d'Arima
, order=c(0,1,1)
indique une dérivation et un processus MA(1), seasonal=c(0,1,0)
une dérivation saisonnière, et lambda=0
un passage préalable au logarithme. Ici on passe aussi l'option method="CSS"
pour éviter des problèmes de convergence de l'algorithme d'optimisation utilisé pour l'estimation par maximum de vraisemblance.
arima_air <- Arima(airmiles,order=c(0,1,1),seasonal=c(0,1,0),lambda=0,method="CSS")
tsdisplay(residuals(arima_air))
Les résidus de ce modèle n'ont pas de corrélations pour des décalages faibles: les corrélations que l'on observait dans la chronique $(1-B)(1-B^{12})\text{Airmiles}_t$ sont donc bien modélisées par un processus MA(1).
Restent les corrélations saisonnières. L'explication la plus simple est: corrélation significative pour un décalage de 12 (dans l'ACF), corrélations partielles en décroissance entre 12 et 24 (dans le PACF), donc plutôt SMA(1).
arima_air <- Arima(airmiles,order=c(0,1,1),seasonal=c(0,1,1),lambda=0,method="CSS")
summary(arima_air)
tsdiag(arima_air,gof.lag=15)
Remarque : d'après help(tsdiag)
: gof.lag
est " the maximum number of lags for a Portmanteau goodness-of-fit test ". On le prend plus grand que la période pour vérifier l'absence de corrélations périodiques résiduelles.
On voit un résidu très négatif en septembre 2001, dont l'effet n'est pas du tout pris en compte dans ce modèle. Il n'y a pas de corrélation significative dans les résidus, les tests de Portmanteau (Ljung-Box) ne permettent pas de rejeter l'hypothèse d'un bruit blanc gaussien.
hist(residuals(arima_air), freq = F, col = "grey")
curve(dnorm(x, mean = mean(residuals(arima_air)), sd = sd(residuals(arima_air))), col = 2, add = TRUE)
Difficile de juger de la gaussianité des résidus sur cet histogramme. Néanmoins, on note un résidu très négatif, qui correspond à septembre 2001.
Maintenant qu'on a validé l'hypothèse de résidus formant un bruit blanc gaussien, il est légitime d'examiner les z-tests de significativité des paramètres.
coeftest(arima_air)
Les paramètres sont effectivement significatifs.
Le modèle retenu est donc le suivant (voir les coefficients en sortie de summary
plus haut):
avec une estimation de la variance des résidus à : 0.003244.
Remarques:
summary
n'estime pas la vraisemblance, et donc pas plus les critères d'information comme AICOn commence par examiner la décomposition par lissage. On voit une rupture de tendance en septembre 2001, qui a aussi un effet dans les résidus:
decomp_air <- decompose(airmiles,type="multiplicative")
plot(decomp_air)
Ceci justifie le modèle d'intervention. On commence par définir deux interventions potentielles: step (1 après septembre 2001, et 0 avant) et pulse (1 en septembre 2001, 0 sinon).
step <- ts(as.numeric(seq(airmiles)>=69),start=c(1996,1),frequency=12)
pulse <- ts(as.numeric(seq(airmiles)==69),start=c(1996,1),frequency=12)
plot(log(airmiles), type="o", pch =20)
lines(step*17.7,col=2)
lines(pulse*17.7,col=3) # 17.7 est une valeur ad-hoc pour que l'intervention apparaisse sur le graphe
grid()
On commence par une régression linéaire simple de la chronique $(1-B)(1-B^{12})\text{log}(\text{airmiles}_t)$ contre la chronique (déterministe) pulse. On dérive la chronique car on sait qu'il s'agit d'une opération nécessaire pour la rendre stationnaire (rappelons que le cas non stationnaire est source de régression fallacieuses).
fitOLS <- lm(diff(diff(log(airmiles),12))~diff(diff(pulse,12)))
summary(fitOLS)
tsdisplay(residuals(fitOLS),lag.max=35)
Le modèle estimé à ce stade est: $$(1-B)(1-B^{12})\text{log}(\text{airmiles}_t) = 0.00145 -0.245232 (1-B)(1-B^{12})\text{pulse}_t + R_t$$
Les résidus R_t ont l'air stationnaires (peu de corrélations significatives, graphe d'aspect différent d'une marche aléatoire, mais ne forment pas un bruit blanc (des corrélations sont significatives). Rappelons qu'on ne peut donc pas interpréter les tests de significativité des coefficients.
Remarquons tout de même la très faible valeur de intercept , il est bien légitime de chercher des modèles sans constantes dès lors qu'on a dérivé deux fois.
On va donc chercher à identifier un modèle ARMA sur ces résidus. En examinant les corrélation court-terme, on peut légitimement hésiter en AR(1) et MA(1). On utilise Arima
, avec l'argument xreg=pulse
pour préciser l'intervention.
arima_airpulse <- Arima(airmiles,order=c(0,1,1),seasonal=c(0,1,0),xreg=pulse,lambda=0)
summary(arima_airpulse)
tsdisplay(residuals(arima_airpulse))
En examinant les corrélations pour les décalages multiples de 12, on peut pencher pour SMA(1), mais on pourrait également essayer SAR(1) (que je ne fais pas ici).
arima_airpulse <- Arima(airmiles,order=c(0,1,1),seasonal=c(0,1,1),xreg=pulse,lambda=0,method="CSS")
summary(arima_airpulse)
tsdisplay(residuals(arima_airpulse))
tsdiag(arima_airpulse)
Les résidus de ce modèle ne sont pas distinguables d'un bruit blanc (pas de corrélation significative, grandes $p$-valeurs pour les tests de Portmanteau). Néanmoins, le graphe des résidus nous montre que le 11-septembre n'a pas un effet limité à septembre 2001: octobre et novembre présentent également des résidus fortement négatifs.
coeftest(arima_airpulse)
Les coefficients du modèle sont significatifs.
Le modèle obtenu s'écrit:
$$ \log(\text{airmiles}_t) = -0.26347 \text{pulse}_t + u_t $$où $$(1-B^{12})(1-B)u_t = (1-0.3632B)(1-0.5301B^{12})\varepsilon_t$$ les résidus $\varepsilon_t$ ayant une variance de 0.001868.
Autrement écrit: $$\text{airmiles}_t = \exp(-0.26347 \text{pulse}_t) * \exp(u_t)$$
Le 11-septembre se traduit donc dans ce modèle par une perte de $1-\exp(-0.26347) = 23.16\%$ de passagers en septembre 2001 (et retour "à la normale" juste après)
La remarque sur les résidus du modèle précédent nous invite à introduire une intervention d'effet plus long que "pulse".
Nous effectuons la même démarche pour l'intervention step.
fitOLS <- lm(diff(diff(log(airmiles),12))~diff(diff(step,12)))
summary(fitOLS)
tsdisplay(residuals(fitOLS))
Les résidus semblent stationnaires mais sont corrélés. On hésite entre AR(1) et MA(1), et on va tester ces deux hypothèses.
# MA(1) ?
arima_airstep <- Arima(airmiles,order=c(0,1,1),seasonal=c(0,1,0),xreg=step,lambda=0,method="CSS")
summary(arima_airstep)
tsdisplay(residuals(arima_airstep))
Difficile d'identifier la composante saisonnière par la méthode de Box et Jenkins: on essaie SAR(1) et SMA(1).
# avec SAR(1)
arima_airstep <- Arima(airmiles,order=c(0,1,1),seasonal=c(1,1,0),xreg=step,lambda=0,method="CSS")
summary(arima_airstep)
tsdisplay(residuals(arima_airstep))
tsdiag(arima_airstep,gof.lag=15)
SAR(1) ne règle pas du tout les corrélations périodiques, qui restent significatives. Cela disqualifie ce modèle.
# avec SMA(1) en plus de MA(1)
arima_airstep <- Arima(airmiles,order=c(0,1,1),seasonal=c(0,1,1),xreg=step,lambda=0,method="CSS")
summary(arima_airstep)
tsdisplay(residuals(arima_airstep))
tsdiag(arima_airstep,gof.lag=15)
# sigma^2 estimated as 0.001638: part log likelihood=180.33
Certaines corrélations sont aux limites du seuil de significativité, mais les $p$-valeurs des tests de Portmanteau sont assez grandes.
# avec AR(1)/SMA(1)
arima_airstep <- Arima(airmiles,order=c(1,1,0),seasonal=c(0,1,1),xreg=step,lambda=0,method="CSS")
summary(arima_airstep)
tsdisplay(residuals(arima_airstep))
tsdiag(arima_airstep,gof.lag=15)
# sigma^2 estimated as 0.001612: part log likelihood=180.63
Un peu le même constat, mais la variance des résidus est un peu plus faible.
# avec AR(1)/SAR(1)
arima_airstep <- Arima(airmiles,order=c(1,1,0),seasonal=c(1,1,0),xreg=step,lambda=0,method="CSS")
summary(arima_airstep)
tsdisplay(residuals(arima_airstep))
tsdiag(arima_airstep,gof.lag=15)
On note des corrélations significatives, et un rejet de l'hypothèse de bruit blanc gaussien par les tests de Portmanteau.
On garde donc le modèle suivant (il faut bien en choisir un):
arima_airstep <- Arima(airmiles,order=c(1,1,0),seasonal=c(0,1,1),xreg=step,lambda=0,method="CSS")
hist(residuals(arima_airstep), freq = F, col = "grey")
curve(dnorm(x, mean = mean(residuals(arima_airstep)), sd = sd(residuals(arima_airstep))), col = 2, add = TRUE)
coeftest(arima_airstep)
summary(arima_airstep)
Les coefficients du modèle sont significatifs. Remarquons que cette fois, les résidus présentent quelques grandes valeurs inattendues (par rapport à une répartition gaussienne): ce sont les valeurs positives relativement grandes observées dans le graphe des résidus juste après septembre 2001. Avec une intervention "step", on a donc tendance à "surcorriger" l'effet du 11-septembre dans les mois qui suivent.
Le modèle s'écrit: $$\log(\text{airmiles}_t) = -0.3323 \text{step}_t + u_t$$ où $$(1-B)(1-B^{12})u_t = \frac{1-0.4839B^12}{1+0.3914B}\varepsilon_t$$ et la variance des résidus est estimée à 0.001612.
Notons que la variance est plus faible qu'avec "pulse", ce dernier modèle "colle" donc mieux aux données.
On peut réécrire le modèle en:
$$\log(\text{airmiles}_t) = \exp(-0.3323 \text{step}_t) * exp(u_t)$$le 11 septembre se traduit donc dans ce modèle par une perte de $1-\exp(-0.3323) = 28.3\%$ de passagers à partir de septembre 2001 (l'effet est permanent avec une intervention "step").
La remarque précédente sur la "surcorrection" de l'effet 11-septembre avec un intervention "step" nous amènent à chercher une introduire une intervention dont l'effet s'atténue au cours du temps.
On a vu dans le cours qu'une possibilité est d'envisager une intervention en $$ \frac{1}{1-\lambda B} \text{pulse}_t$$
Il faut utiliser la fonction arimax
qui permet une fonction de transfert avec un dénominateur. Ici, le dénominateur est de degré 1 et le numérateur de degré 0, d'où l'option transfer=list(c(1,0))
. Comme arimax
ne dispose pas de l'option lambda=0Ì€
, il faut travailler sur la chronique ̀log(airmiles)
.
Par ailleurs, nous allons tenter le même modèle ARIMA sur les résidus que pour l'intervention "step". En toute rigueur, il faudrait examiner les ACF/PACF des résidus du modèle:
$$ \log(\text{airmiles}_t) = \alpha_0 + \alpha_1 \frac{1}{1-\lambda B} \text{pulse}_t + R_t$$
dont les paramètres seraient estimés par moindres carrés des résidus par exemple, puis identifier le modèle sur $R_t$ par la méthode de Box et Jenkins. Malheureusement, arimax
ne le permet pas.
arima_airgeo <- arimax(log(airmiles),order=c(1,1,0),seasonal=c(0,1,1),transfer=list(c(1,0)),xtransf=pulse,method="CSS")
summary(arima_airgeo)
tsdisplay(residuals(arima_airgeo))
tsdiag(arima_airgeo,gof.lag=25)
hist(residuals(arima_airgeo), freq = F, col = "grey")
curve(dnorm(x, mean = mean(residuals(arima_airgeo)), sd = sd(residuals(arima_airstep))), col = 2, add = TRUE)
coeftest(arima_airgeo)
Ce modèle est valide.
summary(arima_airgeo)
Il s'écrit:
$$ \log(\text{airmiles}_t) = \frac{-0.3460}{1-0.6864B} \text{pulse}_t + u_t$$où $$(1-B)(1-B^{12})u_t = \frac{1-0.4238B^{12}}{1+0.4994B}\varepsilon_t$$ La variance des résidus est estimée à 0.001254. Il s'agit de la plus petite variance parmi tous les modèles: ce dernier modèle explique le mieux les données
On peut donc écrire: $$ \text{airmiles}_t = \exp\left(\frac{-0.3460}{1-0.6864B} \text{pulse}_t\right) * \exp(u_t) $$ soit $$ \text{airmiles}_t = \exp\left( -0.3460*\text{pulse}_t - 0.3460*0.6864 \text{pulse}_{t-1} - 0.3460*0.6864^2 \text{pulse}_{t-2}-...\right) * \exp(u_t) $$
Le 11-septembre se traduit donc dans ce modèle par une perte de passagers par rapport à ce qui serait attendu sans intervention 2001:
L'effet 11-septembre s'est atténué à 99\% (son effet a quasiment disparu dans le modèle) lorsque $-0.346*0.6864^k= -0.346*0.01$, soit après $k=\log(0.01)/\log(0.6864) = 12.23$ mois (selon ce modèle particulier).
# représentation graphique de l'intervention à décroissance géométrique dont on a estimé les paramètres juste avant:
geo <- filter(pulse,filter=0.6864,method="rec",sides=1)
plot(log(airmiles))
lines(log(airmiles[68])-0.3460*geo,col=2)
# voici une manière de revenir à Arima pour pouvoir faire des prévisions
# (avec la fonction d'intervention précédemment calculée)
# on crée le modèle jusque mai 2004, puis on prédit sur un horizon de 12.
output_airgeo <- Arima(window(airmiles,end=c(2004,5)),order=c(1,1,0),seasonal=c(0,1,1),xreg=window(geo,end=c(2004,5)),lambda=0)
summary(output_airgeo) # on remarque que l'estimation des paramètres diffère légèrement
# prévision à partir de juin 2004, sur un horizon de 12.
fcast_airgeo <- forecast(output_airgeo,h=12,xreg=window(geo,start=c(2004,6)),lambda=0,biasadj=T)
# représentation graphique:
plot(fcast_airgeo, main=" ", type="o", pch=20) # les valeurs de la chronique jusque mai 2004, puis la prédiction
grid()
lines(window(airmiles,start=c(2004,6)),col=2,type="b") # en rouge, les "vraies" valeurs
# par comparaison: cas de la prévision avec step
output_airstep <- Arima(window(airmiles,end=c(2004,5)),order=c(1,1,0),seasonal=c(0,1,1),xreg=window(step,end=c(2004,5)),lambda=0)
summary(output_airstep)
# prévision à partir de juin 2004, sur un horizon de 12.
fcast_airstep <- forecast(output_airstep,h=12,xreg=window(step,start=c(2004,6)),lambda=0,biasadj=T)
# représentation graphique:
plot(fcast_airstep, main=" ", type="o", pch=20) # les valeurs de la chronique jusque mai 2004, puis la prédiction
grid()
lines(window(airmiles,start=c(2004,6)),col=2,type="b") # en rouge, les "vraies" valeurs
# on constate que la prévision est plus éloignée des vraies valeurs