Etude de la serie airmiles

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 :

  • une variable pulse qui prend la valeur 1 en septembre 2001, 0 sinon ;
  • une variable step qui prend la valeur 1 après septembre 2001, 0 avant ;
  • une variable présentant une décroissance géométrique dont il faudra estimer la vitesse de décroissance.


In [1]:
library(forecast)
library(lmtest)
library(TSA)  

options(repr.plot.width=12, repr.plot.height=10) # pour redimensionner les graphiques
Loading required package: zoo


Attaching package: ‘zoo’


The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric



Attaching package: ‘TSA’


The following objects are masked from ‘package:stats’:

    acf, arima


The following object is masked from ‘package:utils’:

    tar


Représentations graphiques, étude qualitative

In [2]:
data(airmiles)
airmiles
A Time Series: 10 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
1996309831743214766338342975359691133647439138772238403956574173849933580773363898423273490136866146
1997288297942971536937179450340020043458549736981956389126403971597132047526344055233212125434447370
1998338092573268844839908269386185153882380541119273421531934317161035309230385705063629155637511830
1999354505763433778542448030408183784046405743527294469189314584104437829297416510473987574939190232
2000361954023723974045664486432629774446760147261820488265574816057939573944429008704184751540727132
2001383065063696933945686100437116644367482646654145495972735016990727077913341815543474980835726959
2002338956713362097242633492395363544089502444060019471833644689936634776695400930253721399742407172
2003371978093509931642910514403702924187512645592190496334394834383837858815425072874023825343698675
2004391801143973643547876012470504394653413051134050543173145239298541816777472056654455365346316602
20054276065741120838520530594815258550047901

On remarque que les données sont disponibles jusque mai 2005.

In [3]:
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.

In [4]:
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.

In [5]:
plot(diff(diff(log(airmiles),12)))

Le graphe a maintenant l'air de celui d'un processus stationnaire, "à l'intervention près".


Modélisation de la chronique initiale

(sans intervention, pour avoir une référence à laquelle comparer)

On travaillera sur log(airmiles).

In [6]:
tsdisplay(log(airmiles))

L'ACF présente une périodicité: cela confirme la nécessité d'une dérivation saisonnière.

In [7]:
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.

In [8]:
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.

In [9]:
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).

In [10]:
arima_air <- Arima(airmiles,order=c(0,1,1),seasonal=c(0,1,1),lambda=0,method="CSS")
summary(arima_air)
Series: airmiles 
ARIMA(0,1,1)(0,1,1)[12] 
Box Cox transformation: lambda= 0 

Coefficients:
          ma1     sma1
      -0.4718  -0.6466
s.e.   0.0903   0.0848

sigma^2 estimated as 0.003116:  part log likelihood=147.68

Training set error measures:
                   ME    RMSE     MAE       MPE     MAPE      MASE        ACF1
Training set 217155.2 1820153 1020236 0.4245917 2.691143 0.3407108 -0.03530456
In [11]:
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.

In [12]:
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.

In [13]:
coeftest(arima_air)
z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.471818   0.090302 -5.2249 1.742e-07 ***
sma1 -0.646567   0.084776 -7.6268 2.407e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Les paramètres sont effectivement significatifs.


Le modèle retenu est donc le suivant (voir les coefficients en sortie de summary plus haut):

$$(1-B)(1-B^{12}) (1+0.37975B)\log(\text{airmiles}_t) = (1-0.655327B^{12}) \varepsilon_t$$

avec une estimation de la variance des résidus à: 0.003244.

Remarques:

  • le modèle n'incorpore pas de constante (c'est normal car il y a deux dérivation)
  • avec la méthode "CSS", summary n'estime pas la vraisemblance, et donc pas plus les critères d'information comme AIC


Modélisation avec intervention

On 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:

In [14]:
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).

In [15]:
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()


Intervention de type pulse

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

In [16]:
fitOLS <- lm(diff(diff(log(airmiles),12))~diff(diff(pulse,12)))
summary(fitOLS)
tsdisplay(residuals(fitOLS),lag.max=35)
Call:
lm(formula = diff(diff(log(airmiles), 12)) ~ diff(diff(pulse, 
    12)))

Residuals:
      Min        1Q    Median        3Q       Max 
-0.176544 -0.025980  0.002795  0.017717  0.225738 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)            0.001448   0.004879   0.297    0.767    
diff(diff(pulse, 12)) -0.245232   0.024397 -10.052   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04879 on 98 degrees of freedom
Multiple R-squared:  0.5076,	Adjusted R-squared:  0.5026 
F-statistic:   101 on 1 and 98 DF,  p-value: < 2.2e-16

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.

In [17]:
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))
Series: airmiles 
Regression with ARIMA(0,1,1)(0,1,0)[12] errors 
Box Cox transformation: lambda= 0 

Coefficients:
          ma1     xreg
      -0.2911  -0.2511
s.e.   0.1026   0.0265

sigma^2 estimated as 0.002248:  log likelihood=164.88
AIC=-323.75   AICc=-323.5   BIC=-315.94

Training set error measures:
                   ME    RMSE     MAE        MPE     MAPE      MASE        ACF1
Training set 36968.78 1609777 1065378 0.03468143 2.716939 0.3557861 -0.03844535

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

In [18]:
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)
Series: airmiles 
Regression with ARIMA(0,1,1)(0,1,1)[12] errors 
Box Cox transformation: lambda= 0 

Coefficients:
          ma1     sma1     xreg
      -0.3632  -0.5301  -0.2635
s.e.   0.0955   0.1023   0.0311

sigma^2 estimated as 0.001868:  part log likelihood=173.77

Training set error measures:
                   ME    RMSE      MAE       MPE     MAPE      MASE        ACF1
Training set 146326.5 1455987 945333.5 0.3037401 2.399486 0.3156968 -0.05937656

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.

In [19]:
coeftest(arima_airpulse)
z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.363234   0.095490 -3.8039 0.0001424 ***
sma1 -0.530052   0.102317 -5.1805 2.213e-07 ***
xreg -0.263477   0.031076 -8.4784 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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)


Intervention step

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.

In [20]:
fitOLS <- lm(diff(diff(log(airmiles),12))~diff(diff(step,12)))
summary(fitOLS)
tsdisplay(residuals(fitOLS))
Call:
lm(formula = diff(diff(log(airmiles), 12)) ~ diff(diff(step, 
    12)))

Residuals:
      Min        1Q    Median        3Q       Max 
-0.092454 -0.026207  0.001803  0.017326  0.225738 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)           0.001448   0.004535   0.319     0.75    
diff(diff(step, 12)) -0.368984   0.032069 -11.506   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04535 on 98 degrees of freedom
Multiple R-squared:  0.5746,	Adjusted R-squared:  0.5703 
F-statistic: 132.4 on 1 and 98 DF,  p-value: < 2.2e-16

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.

In [21]:
# 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))
Series: airmiles 
Regression with ARIMA(0,1,1)(0,1,0)[12] errors 
Box Cox transformation: lambda= 0 

Coefficients:
          ma1     xreg
      -0.3001  -0.3196
s.e.   0.0989   0.0344

sigma^2 estimated as 0.001898:  part log likelihood=172.46

Training set error measures:
                  ME    RMSE     MAE        MPE     MAPE     MASE        ACF1
Training set 30348.7 1503860 1023016 0.08900405 2.577012 0.341639 -0.04631468

Difficile d'identifier la composante saisonnière par la méthode de Box et Jenkins: on essaie SAR(1) et SMA(1).

In [22]:
# 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)
Series: airmiles 
Regression with ARIMA(0,1,1)(1,1,0)[12] errors 
Box Cox transformation: lambda= 0 

Coefficients:
          ma1     sar1     xreg
      -0.3135  -0.3360  -0.3283
s.e.   0.1181   0.0863   0.0341

sigma^2 estimated as 0.001236:  part log likelihood=188.04

Training set error measures:
                    ME    RMSE      MAE         MPE     MAPE      MASE
Training set -30005.36 1258240 891011.5 -0.09680246 2.215651 0.2975558
                    ACF1
Training set -0.06570495

SAR(1) ne règle pas du tout les corrélations périodiques, qui restent significatives. Cela disqualifie ce modèle.

In [23]:
# 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
Series: airmiles 
Regression with ARIMA(0,1,1)(0,1,1)[12] errors 
Box Cox transformation: lambda= 0 

Coefficients:
          ma1     sma1     xreg
      -0.3721  -0.4744  -0.3179
s.e.   0.0954   0.1014   0.0370

sigma^2 estimated as 0.001638:  part log likelihood=180.33

Training set error measures:
                   ME    RMSE      MAE       MPE     MAPE      MASE        ACF1
Training set 127255.9 1350157 936582.5 0.3358208 2.379255 0.3127744 -0.08314222

Certaines corrélations sont aux limites du seuil de significativité, mais les $p$-valeurs des tests de Portmanteau sont assez grandes.

In [24]:
# 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
Series: airmiles 
Regression with ARIMA(1,1,0)(0,1,1)[12] errors 
Box Cox transformation: lambda= 0 

Coefficients:
          ar1     sma1     xreg
      -0.3914  -0.4839  -0.3323
s.e.   0.0970   0.1012   0.0331

sigma^2 estimated as 0.001612:  part log likelihood=180.63

Training set error measures:
                   ME    RMSE      MAE       MPE     MAPE      MASE        ACF1
Training set 112986.6 1342260 928557.7 0.3045841 2.356858 0.3100945 -0.03845097

Un peu le même constat, mais la variance des résidus est un peu plus faible.

In [25]:
# 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)
Series: airmiles 
Regression with ARIMA(1,1,0)(1,1,0)[12] errors 
Box Cox transformation: lambda= 0 

Coefficients:
          ar1     sar1     xreg
      -0.3170  -0.3401  -0.3382
s.e.   0.1047   0.0831   0.0284

sigma^2 estimated as 0.001181:  part log likelihood=189.74

Training set error measures:
                  ME    RMSE      MAE         MPE     MAPE      MASE
Training set -8397.6 1235928 864440.4 -0.03169912 2.135571 0.2886823
                    ACF1
Training set -0.04926526

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):

In [26]:
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)
z test of coefficients:

      Estimate Std. Error  z value  Pr(>|z|)    
ar1  -0.391435   0.096953  -4.0374 5.406e-05 ***
sma1 -0.483859   0.101243  -4.7792 1.760e-06 ***
xreg -0.332340   0.033059 -10.0529 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Series: airmiles 
Regression with ARIMA(1,1,0)(0,1,1)[12] errors 
Box Cox transformation: lambda= 0 

Coefficients:
          ar1     sma1     xreg
      -0.3914  -0.4839  -0.3323
s.e.   0.0970   0.1012   0.0331

sigma^2 estimated as 0.001612:  part log likelihood=180.63

Training set error measures:
                   ME    RMSE      MAE       MPE     MAPE      MASE        ACF1
Training set 112986.6 1342260 928557.7 0.3045841 2.356858 0.3100945 -0.03845097

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


Intervention en "décroissance géométrique"

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.

In [27]:
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)
Call:
arimax(x = log(airmiles), order = c(1, 1, 0), seasonal = c(0, 1, 1), method = "CSS", 
    xtransf = pulse, transfer = list(c(1, 0)))

Coefficients:
          ar1     sma1  T1-AR1   T1-MA0
      -0.4994  -0.4238  0.6864  -0.3460
s.e.   0.1149   0.1147  0.0709   0.0267

sigma^2 estimated as 0.001254:  part log likelihood = 192.18

Training set error measures:
                      ME       RMSE        MAE        MPE      MAPE      MASE
Training set 0.003086045 0.03314449 0.02046566 0.01755449 0.1169867 0.2336676
                    ACF1
Training set -0.05127026
z test of coefficients:

        Estimate Std. Error  z value  Pr(>|z|)    
ar1    -0.499397   0.114949  -4.3445 1.396e-05 ***
sma1   -0.423836   0.114672  -3.6961  0.000219 ***
T1-AR1  0.686414   0.070898   9.6817 < 2.2e-16 ***
T1-MA0 -0.345992   0.026750 -12.9344 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ce modèle est valide.

In [28]:
summary(arima_airgeo)
Call:
arimax(x = log(airmiles), order = c(1, 1, 0), seasonal = c(0, 1, 1), method = "CSS", 
    xtransf = pulse, transfer = list(c(1, 0)))

Coefficients:
          ar1     sma1  T1-AR1   T1-MA0
      -0.4994  -0.4238  0.6864  -0.3460
s.e.   0.1149   0.1147  0.0709   0.0267

sigma^2 estimated as 0.001254:  part log likelihood = 192.18

Training set error measures:
                      ME       RMSE        MAE        MPE      MAPE      MASE
Training set 0.003086045 0.03314449 0.02046566 0.01755449 0.1169867 0.2336676
                    ACF1
Training set -0.05127026

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:

  • en septembre 2001: $1-\exp(-0.346) = 29.2\%$
  • en octobre 2001: $1-\exp(-0.346*0.6864) = 21.0\%$
  • en novembre 2001: $1-\exp(-0.346*0.6864^2) = 14.7\%$
  • ...


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


Bonus: prévision avec modèle d'intervention

In [29]:
# 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)
In [30]:
# 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
Series: window(airmiles, end = c(2004, 5)) 
Regression with ARIMA(1,1,0)(0,1,1)[12] errors 
Box Cox transformation: lambda= 0 

Coefficients:
          ar1     sma1     xreg
      -0.4347  -0.6128  -0.3501
s.e.   0.0971   0.1823   0.0293

sigma^2 estimated as 0.0013:  log likelihood=167.75
AIC=-327.51   AICc=-327.02   BIC=-317.6

Training set error measures:
                   ME    RMSE      MAE       MPE     MAPE      MASE        ACF1
Training set 88364.22 1217054 812787.9 0.1611785 2.085826 0.2794969 -0.06683546
In [31]:
# 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
Series: window(airmiles, end = c(2004, 5)) 
Regression with ARIMA(1,1,0)(0,1,1)[12] errors 
Box Cox transformation: lambda= 0 

Coefficients:
          ar1     sma1     xreg
      -0.3537  -0.6377  -0.3379
s.e.   0.1048   0.1506   0.0358

sigma^2 estimated as 0.001621:  log likelihood=157.44
AIC=-306.88   AICc=-306.4   BIC=-296.97

Training set error measures:
                   ME    RMSE      MAE      MPE     MAPE      MASE        ACF1
Training set 130290.2 1322209 916192.5 0.311884 2.377551 0.3150551 -0.05981425
In [ ]: