TP 6: éléments de correction


On commence par charger les bibliothèques utiles.

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

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



Attaching package: ‘aTSA’


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

    identify



Attaching package: ‘forecast’


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

    forecast



Exercice 1: chronique SNCF


Lecture des données:

In [2]:
sncf <- read.csv(file="sncf.csv",header=TRUE,sep=";")
sncf <- as.vector(t(as.matrix(sncf[,2:13])))
sncf <- ts(sncf,start = c(1963, 1), frequency = 12)
sncf
A Time Series: 18 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
1963175015601820209019102410314028502090185016302420
1964171016001800212021002460320029602190187017702270
1965167016401770219020202610319028602140187017602360
1966181016401860199021102500303029002160194017502330
1967185015901880221021102480288026702100192016702520
1968183417921860213821152485258126392038193617842391
1969179818501981208521202491283427251932208518562553
1970185418232005241822192722291227712153213619102537
1971200818352120230422642175292827382178213720092546
1972208420342152252223182684297127592267215219782723
1973208121122279266122812929308928032296221021352862
1974222322482421271025053021332730442607252521602876
1975248124282596292327953287359831182875275425883266
1976266726682804280629763430370530532764280227073307
1977270625862796297830533463364930952839296628633375
1978282028573306333331413512374431792984295028963611
1979331326442872326733913682393732842849308530433541
1980284829133248325033753640377132593206326931814008

Représentation graphique:

In [3]:
plot(sncf, type="o", pch =20)
grid()

On observe une périodicité annuelle, ce qui semble naturel étant donnée la nature du trafic SNCF.

Regardons le corrélogramme et corrélogramme partiel:

In [4]:
tsdisplay(sncf)

L'aspect "en pont suspendu" de l'ACF est typique d'une chronique non stationnaire qu'il faut dériver saisonnièrement.

Examinons les graphes de la chronique dérivée.

In [5]:
difsncf <- diff(sncf,lag=12)
difsncf2 <- diff(difsncf)

plot(difsncf, type="o", pch=20)
grid()

La chronique dérivée saisonnièrement ne présente plus de saisonnalité flagrante. Néanmoins, son graphe peut faire penser à celui d'une marche aléatoire (qui n'est pas stationnaire et nécessite une dérivation pour le devenir).

In [6]:
plot(difsncf2, type="o", pch=20)
grid()

Effectivement, cette chronique (sncf avec deux dérivations dont une saisonnière) a un graphe d'aspect davantage stationnaire. Bien entendu, ces considérations basées sur les graphiques sont assez empiriques.

Examinons les corrélogrammes de SNCF dérivée saisonnièrement:

In [7]:
arima_sncf <- Arima(sncf,seasonal=c(0,1,0))
tsdisplay(residuals(arima_sncf))

Lorsqu'on observe les corrélations (ACF et PACF) à court-terme (décalages < 10), on constate qu'elles sont assez importantes et ne décroissent pas franchement. C'est signe de non-stationnarité.

Par acquit de conscience, regardons la sortie des tests ADF:

In [8]:
aTSA::adf.test(residuals(arima_sncf))
Augmented Dickey-Fuller Test 
alternative: stationary 
 
Type 1: no drift no trend 
     lag   ADF p.value
[1,]   0 -9.13    0.01
[2,]   1 -6.80    0.01
[3,]   2 -5.16    0.01
[4,]   3 -3.73    0.01
[5,]   4 -3.17    0.01
Type 2: with drift no trend 
     lag    ADF p.value
[1,]   0 -10.53    0.01
[2,]   1  -8.19    0.01
[3,]   2  -6.48    0.01
[4,]   3  -4.91    0.01
[5,]   4  -4.31    0.01
Type 3: with drift and trend 
     lag    ADF p.value
[1,]   0 -11.23    0.01
[2,]   1  -8.94    0.01
[3,]   2  -7.23    0.01
[4,]   3  -5.60    0.01
[5,]   4  -5.00    0.01
---- 
Note: in fact, p.value = 0.01 means p.value <= 0.01 

Tous les tests sont malgré tout en faveur de la stationnarité. Les ACF / PACF précédents ne sont pas typiques d'un AR(p) ou MA(q): cela signifie que si la chronique difsncf est effectivement stationnaire, il faut un modèle "mélange" ARMA(p,q). Je ne suis pas parvenu à trouver un tel modèle. N'oubliez pas que les tests ADF ne nous permettent que de rejeter des hypothèses (type 1 et lag 0, type 1 et lag 1... type 2 et lag 0, type 2 et lag 1...). La non-stationnarité peut s'expliquer par d'autres hypothèses, voir la discussion dans le polycopié.

Nous allons donc dériver dans l'objectif de stationnariser.

In [9]:
arima_sncf <- Arima(sncf,order=c(0,1,0),seasonal=c(0,1,0))
tsdisplay(residuals(arima_sncf))

Nous nous retrouvons dans une situation dans laquelle on sait identifier un modèle.

On commence par ACF et PACF à "court terme" On observe des pics pour lag=1,2 (limite à 2) dans l'ACF et un PACF qui "décroit", signe d'un processus MA(2) (on pourrait essayer MA(1) )

In [10]:
arima_sncf <- Arima(sncf,order=c(0,1,2),seasonal=c(0,1,0))
tsdisplay(residuals(arima_sncf))

Les résidus de ce modèle ne présentent plus de corrélations à court-terme, signe que notre choix semble bon.

On examine les ACF/PACF pour les décalages multiples de 12: on n'a pas franchement de raison de décider entre SAR(1) et SMA(1), on va donc essayer les deux modèles.

In [11]:
arima_sncf <- Arima(sncf,order=c(0,1,2),seasonal=c(0,1,1))
tsdiag(arima_sncf)

Ce modèle (deux dérivations, MA(2), et SMA(1)) fournit des résidus non corrélés, indiscernables d'un bruit blanc gaussien selon les statistiques de Ljung-Box (tests de Portmanteau).

In [12]:
hist(residuals(arima_sncf), freq = F, col = "grey")
curve(dnorm(x, mean = mean(residuals(arima_sncf)), sd = sd(residuals(arima_sncf))), col = 2, add = TRUE)

coeftest(arima_sncf)
z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.669790   0.067904 -9.8637 < 2.2e-16 ***
ma2  -0.188397   0.065424 -2.8796  0.003981 ** 
sma1 -0.551612   0.057147 -9.6525 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

La répartition des résidus a un aspect gaussien.

Rappel: ce qui est important ici est de ne pas rencontrer une distribution franchement non gaussienne, par exemple bimodale, ce qui serait le signe d'un problème de modélisation.

On peut alors interprêter les z-tests de coeftest(construits sous hypothèses résidus = b.b.g.): les coefficients sont tous significatifs.

Rappelons aussi que dans ce modèle où deux dérivations interviennent, il n'y a pas de coefficient constant.

In [13]:
summary(arima_sncf)
Series: sncf 
ARIMA(0,1,2)(0,1,1)[12] 

Coefficients:
          ma1      ma2     sma1
      -0.6698  -0.1884  -0.5516
s.e.   0.0679   0.0654   0.0571

sigma^2 estimated as 15684:  log likelihood=-1269.98
AIC=2547.95   AICc=2548.15   BIC=2561.21

Training set error measures:
                  ME     RMSE      MAE       MPE     MAPE      MASE        ACF1
Training set 6.92879 120.5071 86.19991 0.3002854 3.364295 0.6972554 0.004168261

On note RMSE= 120.5 AIC=2547.95 BIC=2561.21, le modèle s'écrit:

$$ (1-B)(1-B^{12})\text{sncf}_t=(1-0.6698B-0.1884B^2)(1-0.5516B^{12})\varepsilon_t $$

On examine à présent l'hypothèse SAR(1):

In [14]:
arima_sncf <- Arima(sncf,order=c(0,1,2),seasonal=c(1,1,0))
tsdiag(arima_sncf)

Les autocorrélations sont un peu plus importantes (donc les $p$⁻valeurs des Portmanteau plus faibles, mais néanmoins au dessus du risque à 5%).

L'histogramme correspond à ce qu'on attend, les coefficients du modèle sont significatifs:

In [15]:
hist(residuals(arima_sncf), freq = F, col = "grey")
curve(dnorm(x, mean = mean(residuals(arima_sncf)), sd = sd(residuals(arima_sncf))), col = 2, add = TRUE)

coeftest(arima_sncf)
z test of coefficients:

      Estimate Std. Error  z value  Pr(>|z|)    
ma1  -0.699105   0.067057 -10.4256 < 2.2e-16 ***
ma2  -0.160514   0.067398  -2.3816   0.01724 *  
sar1 -0.522773   0.065125  -8.0273 9.966e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
In [16]:
summary(arima_sncf)
Series: sncf 
ARIMA(0,1,2)(1,1,0)[12] 

Coefficients:
          ma1      ma2     sar1
      -0.6991  -0.1605  -0.5228
s.e.   0.0671   0.0674   0.0651

sigma^2 estimated as 16073:  log likelihood=-1272.19
AIC=2552.38   AICc=2552.58   BIC=2565.63

Training set error measures:
                   ME     RMSE      MAE       MPE     MAPE     MASE        ACF1
Training set 4.740897 121.9934 87.32239 0.1887029 3.420704 0.706335 0.006514628

On note: RMSE= 121.99 AIC=2552.38 BIC=2565.63 et on écrit le modèle:

$$ (1+0.52273B^{12})(1-B)(1-B^{12})\text{sncf}_t = (1-0.699105B-0.160514B^2)\varepsilon_t $$

(attention au signe pour la partie SAR)

RMSE et critères d'information sont (un peu) moins bons que pour le premier modèle, auquel on revient pour la prévision. (il faut bien choisir)


Remarque au passage: on aurait pu envisager un modèle du type $(0,1,1)(0,1,1)_{12}$ au moment où on hésitait entre MA(2) et MA(1). Je vous invite à poursuivre l'étude: ce modèle s'avère valide mais explique moins RMSE et AIC/BIC sont moins bons.


Passons à la prévision sur l'année 1981:

In [17]:
arima_sncf <- Arima(sncf,order=c(0,1,2),seasonal=c(0,1,1))

# (attention, selon l'ordre de chargement des packages on peut se retrouver à masquer forecast dans ce cas ajouter "forecast::" devant forecast)
fcast <- forecast::forecast(arima_sncf,h=12)
summary(fcast)
plot(fcast, type="o", pch=20)
grid()
Forecast method: ARIMA(0,1,2)(0,1,1)[12]

Model Information:
Series: sncf 
ARIMA(0,1,2)(0,1,1)[12] 

Coefficients:
          ma1      ma2     sma1
      -0.6698  -0.1884  -0.5516
s.e.   0.0679   0.0654   0.0571

sigma^2 estimated as 15684:  log likelihood=-1269.98
AIC=2547.95   AICc=2548.15   BIC=2561.21

Error measures:
                  ME     RMSE      MAE       MPE     MAPE      MASE        ACF1
Training set 6.92879 120.5071 86.19991 0.3002854 3.364295 0.6972554 0.004168261

Forecasts:
         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
Jan 1981       3219.704 3059.210 3380.199 2974.249 3465.160
Feb 1981       3017.417 2848.398 3186.435 2758.925 3275.908
Mar 1981       3315.408 3144.864 3485.952 3054.584 3576.232
Apr 1981       3447.854 3275.798 3619.910 3184.717 3710.991
May 1981       3511.800 3338.245 3685.355 3246.370 3777.229
Jun 1981       3827.575 3652.534 4002.616 3559.873 4095.277
Jul 1981       4021.953 3845.439 4198.468 3751.998 4291.909
Aug 1981       3466.938 3288.962 3644.914 3194.748 3739.129
Sep 1981       3258.280 3078.855 3437.705 2983.873 3532.687
Oct 1981       3346.408 3165.545 3527.271 3069.802 3623.014
Nov 1981       3267.071 3084.782 3449.361 2988.283 3545.859
Dec 1981       3957.400 3773.695 4141.105 3676.447 4238.352

On peut restreindre l'affichage aux dernières années:

In [18]:
plot(fcast, type="o", pch=20, xlim=c(1976,1982))
grid()

On peut également faire une prévision sur l'année 1980, que l'on confronte aux vraies valeurs de cette année:

In [19]:
arima_sncf_b <- Arima(window(sncf,end=c(1979,12)),order=c(0,1,2),seasonal=c(0,1,1))
fcast_b <- forecast::forecast(arima_sncf_b,h=12)
plot(fcast_b, main=" ", type="o", pch=20, xlim=c(1976,1981))
lines(window(sncf,start=c(1979,12)),co=2,type="o")
grid()

On observe que les vraies valeurs sont (à peu de choses près) dans les intervalles de confiance de la prévision (en bleu-gris).



Exercice 2: chronique Busban


Lecture des données:

In [20]:
busban <- read.table(file="busban.txt",sep='\t',header=TRUE)
busban <- ts(busban,start=c(1984,1),frequency=12)
busban
A Time Series: 12 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
19841484149615211514147514311122 8141393151215011558
19851527151915331531148014891147 8401433153615481567
19861566154115821586158615101185 8601459156215821613
19871504152215731598156415461183 8911519161716311658
19881635163016411636163515851222 9131557164316551678
19891639166016821656168216431258 9691588167717181714
19901714173817331715172616501292 9801604167516801729
19911678167416981711168016911290 9611605169717331732
19921687171517281721170216521274 9871668172817921801
19931770178617741771177617011285 9721693179518151829
19941821182918341809180317331301 9831727181818351838
1995179417991826182818071738133310111735183918301871
In [21]:
plot(busban,type="o", pch =20)
grid()

Forte périodicité, tendance à la hausse. Il est logique de voir des valeurs de trafic sur les bus de banlieue plus faibles en été (spécialement en août).

On effectue une décomposition par lissage:

In [22]:
decomp_busban <- decompose(busban)
plot(decomp_busban)

On voit dans les résidus une valeur aberrante en janvier 1987, mais ce n'est pas flagrant. Il y a aussi manifestement des comportements étranges aux étés 1993 et 1994.

Traçons le graphe de la chronique corrigée des variations saisonnières (CVS): c'est la somme de la chronique trend et de la chronique random (autrement dit, la chronique CVS est la chronique initiale à laquelle on a retiré la composante périodique).

In [23]:
plot(busban,type="o", pch =20)
lines(decomp_busban$trend+decomp_busban$random,col=2,type="o")
grid()

Le décrochage en janvier 1987 est bien visible sur la chronique CVS.

A ce stade, on peut s'interroger aussi sur les étés 93 et 94: on y reviendra.


On commence par un modèle SARIMA de la chronique initiale: il servira de point de comparaison pour les modèles d'intervention.

In [24]:
plot(diff(busban,12))  
aTSA::adf.test(diff(busban,12))
Augmented Dickey-Fuller Test 
alternative: stationary 
 
Type 1: no drift no trend 
     lag   ADF p.value
[1,]   0 -4.48  0.0100
[2,]   1 -3.10  0.0100
[3,]   2 -2.46  0.0155
[4,]   3 -2.11  0.0358
[5,]   4 -2.09  0.0382
Type 2: with drift no trend 
     lag   ADF p.value
[1,]   0 -6.20  0.0100
[2,]   1 -4.45  0.0100
[3,]   2 -3.66  0.0100
[4,]   3 -3.21  0.0228
[5,]   4 -3.31  0.0183
Type 3: with drift and trend 
     lag   ADF p.value
[1,]   0 -6.24  0.0100
[2,]   1 -4.50  0.0100
[3,]   2 -3.73  0.0245
[4,]   3 -3.29  0.0761
[5,]   4 -3.42  0.0541
---- 
Note: in fact, p.value = 0.01 means p.value <= 0.01 

Le graphe de la différenciation saisonnière a-t-il un aspect de "marche aléaatoire" ou bien d'un AR avec racine proche de 1? C'est bien entendu discutable.

Les tests ADF suggèrent que la stationnarité: on essaiera un modèle ARIMA avec et sans différantiation supplémentaire.

In [25]:
tsdisplay(busban)
tsdisplay(diff(busban,12))

Comme on s'en doutait, l'ACF de la chronique initiale est périodique: il faut introduire une dérivée saisonnière.

Ensuite, si on n'introduit pas de dérivation, on identifie pour les décalages faibles AR(1) ou AR(2):

In [26]:
arima_bus <- Arima(busban,order=c(1,0,0),seasonal=c(0,1,0),include.mean=T)
tsdisplay(residuals(arima_bus))
arima_bus <- Arima(busban,order=c(2,0,0),seasonal=c(0,1,0),include.mean=T)
tsdisplay(residuals(arima_bus))

C'est bien AR(2) qui s'impose pour obtenir des corrélations court-termes (presque) non significatives.

Pour la composante saisonnière, on teste SAR(1) et SMA(1):

In [27]:
# SAR(1)
arima_bus <- Arima(busban,order=c(2,0,0),seasonal=c(1,1,0),include.mean=T)
tsdisplay(residuals(arima_bus))
tsdiag(arima_bus,gof.lag=15)

# SMA(1)
arima_bus <- Arima(busban,order=c(2,0,0),seasonal=c(0,1,1),include.mean=T)
tsdisplay(residuals(arima_bus))
tsdiag(arima_bus,gof.lag=15)

On ne trouve pas de modèle validant les hypothèses.

diff(busban(12)) ne semble donc finalement pas stationnaire, il faut donc "encore" dériver:

In [28]:
arima_bus <- Arima(busban,order=c(0,1,0),seasonal=c(0,1,0))
tsdisplay(residuals(arima_bus))

On identifie cette fois assez clairement un MA(1):

In [29]:
arima_bus <- Arima(busban,order=c(0,1,1),seasonal=c(0,1,0))
tsdisplay(residuals(arima_bus))

Puis on hésite entre SMA(1) et SAR(1):

In [30]:
arima_bus <- Arima(busban,order=c(0,1,1),seasonal=c(0,1,1))
tsdisplay(residuals(arima_bus))
tsdiag(arima_bus,gof.lag=15)
hist(residuals(arima_bus), freq = F, col = "grey")
curve(dnorm(x, mean = mean(residuals(arima_bus)), sd = sd(residuals(arima_bus))), col = 2, add = TRUE)
coeftest(arima_bus)
summary(arima_bus)
z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.488894   0.092858 -5.2650 1.402e-07 ***
sma1 -0.477789   0.073367 -6.5123 7.401e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Series: busban 
ARIMA(0,1,1)(0,1,1)[12] 

Coefficients:
          ma1     sma1
      -0.4889  -0.4778
s.e.   0.0929   0.0734

sigma^2 estimated as 579.7:  log likelihood=-603.29
AIC=1212.57   AICc=1212.76   BIC=1221.2

Training set error measures:
                     ME     RMSE      MAE         MPE   MAPE      MASE
Training set -0.3526258 22.78819 17.40654 -0.09242704 1.1425 0.5005802
                   ACF1
Training set 0.07304285
In [31]:
arima_bus <- Arima(busban,order=c(0,1,1),seasonal=c(1,1,0))
tsdisplay(residuals(arima_bus))
tsdiag(arima_bus,gof.lag=15)
hist(residuals(arima_bus), freq = F, col = "grey")
curve(dnorm(x, mean = mean(residuals(arima_bus)), sd = sd(residuals(arima_bus))), col = 2, add = TRUE)
coeftest(arima_bus)
summary(arima_bus)
z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.465656   0.093349 -4.9883 6.091e-07 ***
sar1 -0.451991   0.077340 -5.8442 5.091e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Series: busban 
ARIMA(0,1,1)(1,1,0)[12] 

Coefficients:
          ma1     sar1
      -0.4657  -0.4520
s.e.   0.0933   0.0773

sigma^2 estimated as 586.1:  log likelihood=-603.81
AIC=1213.62   AICc=1213.81   BIC=1222.24

Training set error measures:
                    ME     RMSE      MAE         MPE     MAPE      MASE
Training set -0.383825 22.91401 17.50733 -0.08195647 1.142652 0.5034788
                 ACF1
Training set 0.067509

Le modèle en SMA(1) précédent est légèrement meilleur (sigma^2, AIC) on y revient:

In [32]:
arima_bus <- Arima(busban,order=c(0,1,1),seasonal=c(0,1,1))
summary(arima_bus)
plot(residuals(arima_bus))
grid()
residuals(arima_bus)
Series: busban 
ARIMA(0,1,1)(0,1,1)[12] 

Coefficients:
          ma1     sma1
      -0.4889  -0.4778
s.e.   0.0929   0.0734

sigma^2 estimated as 579.7:  log likelihood=-603.29
AIC=1212.57   AICc=1212.76   BIC=1221.2

Training set error measures:
                     ME     RMSE      MAE         MPE   MAPE      MASE
Training set -0.3526258 22.78819 17.40654 -0.09242704 1.1425 0.5005802
                   ACF1
Training set 0.07304285
A Time Series: 12 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
1984 0.85678758 0.39246220 0.27805560 0.20107692 0.12512701 0.06363781 -0.23316702 -0.49321389 0.10944710 0.21166481 0.18186965 -0.51893952
1985 -4.81217143-16.22190169-16.67070010 -3.43967429-12.48574479 41.71188275 -9.38978891 -3.68926720 10.82728437 -9.14490035 16.28032515-26.32955611
1986 18.54660489-15.54189379 14.71290486 14.93882505 52.66695461-37.31471723-14.16476262-24.17439982 -0.60871258 -6.39408534 13.47156597 3.86111985
1987-90.77789831-12.30402523 14.31196150 31.51468081 2.69941692 29.67802197-21.41616235 14.34262199 41.09990569 12.27679658 7.76099795 -1.45658577
1988 41.57767808 12.25141881-24.31485942-30.20020256 12.16500714-12.54626670-23.20147191-16.53699257 24.09304646 -3.92206618 -3.07973008 -7.99609253
1989 0.22440117 22.27682121 7.42605595-26.10412639 27.08137358 15.41258621-22.60386477 6.46590501 -6.47843037 -7.66010137 24.69348299-18.01809493
1990 32.14896692 29.30378627-14.32657461-13.20731606 -2.42496159-37.14574379 -5.55767150-17.34897736 -8.08642271-24.09842001-34.19602144 21.90612269
1991-20.72453021-31.63869615 -0.15724030 27.95909759-26.40411137 56.90985202 -9.15553670-28.46635535 6.27161077 14.44147108 27.35131197-18.17440629
1992-17.90316697 12.97146292 2.65683310 -5.30594863 -9.74030814-32.40371420-10.50969134 25.39963835 59.06340889 2.31078285 38.82429972 13.90870209
1993 16.49125906 6.44201808-23.61108193-10.69904734 15.35488616-30.69993841-50.46129644-36.07966776 44.64753996 51.13548775 -0.99012147 2.09245878
1994 28.65340723 5.23425390 6.77310478-18.28527673-10.10400456-18.19459945-41.83393110-30.90366867 37.65124474 21.41030584 -4.95033913-12.18915763
1995-28.75770327-21.25171181 13.62360989 23.34188231 -4.14463676 -7.35929203 7.66433526 -5.24646182 2.64312643 15.72692454-25.67761034 20.77885328

On remarque qu'un résidu semble aberrant en janvier 87

Contrairement aux étés 93 et 94 (qui ne présentent pour leur part pas de résidu aberrant), le "décrochage" en janvier 1987 ne s'explique pas par la dynamique propre de la chronique modélisée par ce processus SARIMA.

Ce modèle s'écrit:

$$ (1-B)(1-B^{12}) \text{busban}_t = (1-0.4889 B)(1 -0.4778 B^{12}) \varepsilon_t$$


On passe au modèle avec intervention.

Avec pulse l'intervention n'aura qu'un effet très limité sur la modélisation: on peut directement tester le modèle "sans intervention" obtenu par la méthode de Box et Jenkins.

In [33]:
pulse <- ts(as.numeric(time(busban)==1987),start=c(1984,1),frequency=12)
plot(busban, type="o", pch =20)
lines(pulse*1800,col=3)
grid()
In [34]:
arima_buspulse <- Arima(busban,order=c(0,1,1),seasonal=c(0,1,1),xreg=pulse)
summary(arima_buspulse)
tsdisplay(residuals(arima_buspulse))
tsdiag(arima_buspulse,gof.lag=15)
hist(residuals(arima_buspulse), freq = F, col = "grey")
curve(dnorm(x, mean = mean(residuals(arima_buspulse)), sd = sd(residuals(arima_buspulse))), col = 2, add = TRUE)
coeftest(arima_buspulse)
Series: busban 
Regression with ARIMA(0,1,1)(0,1,1)[12] errors 

Coefficients:
          ma1    sma1      xreg
      -0.4998  -0.434  -64.4026
s.e.   0.0917   0.077   16.9059

sigma^2 estimated as 528.4:  log likelihood=-596.41
AIC=1200.82   AICc=1201.13   BIC=1212.32

Training set error measures:
                     ME     RMSE      MAE         MPE     MAPE      MASE
Training set -0.3752529 21.67195 17.06438 -0.09033569 1.119376 0.4907402
                   ACF1
Training set 0.06520449
z test of coefficients:

       Estimate Std. Error z value  Pr(>|z|)    
ma1   -0.499806   0.091662 -5.4527 4.961e-08 ***
sma1  -0.433989   0.077034 -5.6337 1.764e-08 ***
xreg -64.402601  16.905914 -3.8095 0.0001393 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Toutes les hypothèses sont validées.

Au passage: il n'y a plus le résidu très négatif en janvier 1987 que l'on voyait dans l'histogramme du modèle sans intervention.

In [35]:
summary(arima_buspulse)
Series: busban 
Regression with ARIMA(0,1,1)(0,1,1)[12] errors 

Coefficients:
          ma1    sma1      xreg
      -0.4998  -0.434  -64.4026
s.e.   0.0917   0.077   16.9059

sigma^2 estimated as 528.4:  log likelihood=-596.41
AIC=1200.82   AICc=1201.13   BIC=1212.32

Training set error measures:
                     ME     RMSE      MAE         MPE     MAPE      MASE
Training set -0.3752529 21.67195 17.06438 -0.09033569 1.119376 0.4907402
                   ACF1
Training set 0.06520449

Le modèle identifié s'écrit:

$$ \text{busban}_t = -64.4 \text{pulse}_t + u_t $$

avec
$$(1-B)(1-B^{12}) u_t = (1-0.4998B)(1-0.434B^{12}) \varepsilon_t$$

RMSE et critères d'information sont plus petits que dans le modèle sans intervention.

Selon ce modèle, on peut quantifier la perte de passager en janvier 1987 à 64400 passagers par jour (l'unité de la chronique est le millier de passagers par jour pendant le mois considéré).

L'intervalle de confiance de cette estimation est fournie par:

In [36]:
coefci(arima_buspulse)
A matrix: 3 × 2 of type dbl
2.5 %97.5 %
ma1 -0.6794606 -0.3201522
sma1 -0.5849737 -0.2830047
xreg-97.5375831-31.2676194


Avec une intervention décroissante au cours du temps:

In [37]:
arima_busgeo <- arimax(busban,order=c(0,1,1),seasonal=c(0,1,1),transfer=list(c(1,0)),xtransf=pulse)
summary(arima_busgeo)
tsdisplay(residuals(arima_busgeo))
hist(residuals(arima_busgeo), freq = F, col = "grey")
curve(dnorm(x, mean = mean(residuals(arima_busgeo)), sd = sd(residuals(arima_busgeo))), col = 2, add = TRUE)
coeftest(arima_busgeo)
Call:
arimax(x = busban, order = c(0, 1, 1), seasonal = c(0, 1, 1), xtransf = pulse, 
    transfer = list(c(1, 0)))

Coefficients:
          ma1     sma1  T1-AR1    T1-MA0
      -0.5030  -0.4391  0.0052  -64.8518
s.e.   0.0915   0.0764  0.0211   16.9798

sigma^2 estimated as 515.4:  log likelihood = -596.36,  aic = 1200.71

Training set error measures:
                     ME     RMSE      MAE         MPE     MAPE      MASE
Training set -0.3774758 21.65724 17.04685 -0.09110128 1.118472 0.1255317
                   ACF1
Training set 0.06841138
z test of coefficients:

          Estimate  Std. Error z value  Pr(>|z|)    
ma1     -0.5030052   0.0914796 -5.4986 3.829e-08 ***
sma1    -0.4391477   0.0764366 -5.7453 9.178e-09 ***
T1-AR1   0.0052122   0.0210533  0.2476 0.8044668    
T1-MA0 -64.8517867  16.9798427 -3.8193 0.0001338 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ce modèle s'écrit:

$$ \text{busban}_t = \frac{-64.85}{1-0.0052B} \text{pulse}_t + u_t $$

avec
$$(1-B)(1-B^{12}) u_t = (1-0.503B)(1-0.439B^{12}) \varepsilon_t $$

Néanmoins, le coefficient au dénominateur de l'intervention n'est pas significatif.

Conclusion: il n'est donc pas pertinent de chercher une intervention en $ 1/(1-\delta B) \text{pulse}_t $ (puisque $\delta$ n'est pas significatif)

L'effet de la grève est bien ponctuel. Cette conclusion est cohérente avec le graphe des résidus du modèle sans intervention dans lequel on ne voyait qu'un seul résidu aberrant en janvier 1987, la situation semblant être revenue à la normale dès février 1987.


Voyons le côté positif du confinement: vos successeurs auront de beaux exercices avec modèle d'intervention!



Exercice 3: SOI et recruit


On charge la bibliothèque dynlm, les données, on les met en forme, et on fait un graphique (sur lequel il est difficile de voir grand-chose...):

In [38]:
library(dynlm)

soirecrut <- read.table("soirecruit.txt",sep='\t',header=T)

X=ts(soirecrut[,1],start=c(1950,1),frequency=12)  # soi
Y=ts(soirecrut[,2],start=c(1950,1),frequency=12)  # recruit

par(mar=c(5,4,4,4))  
plot(X, col="blue", ylab="",type="o",pch=20)
mtext("SOI",side=2,line=2,col="blue")
par(new=T)   # pour superposer le prochain graphe
plot(Y, col="red", type="o", xaxt="n", yaxt="n", xlab="", ylab="", pch=20)
axis(4)  # tracer l'axe à droite
mtext("Recruit",side=4,line=2,col="red")
grid()

On cherche à expliquer $Y$ en fonction de $X$. On va commencer par chercher un filtre blanchissant $X$, que l'on appliquera à $Y$. Ensuite on pourra examiner le corrélogramme croisé de ces chronique filtrées afin d'identifier la fonction de transfert liant $Y$ à $X$.

In [39]:
tsdisplay(X)

Belle périodicité de l'ACF: il y a saisonnalité sur une période de 12 mois.

In [40]:
tsdisplay(diff(X,12))

La décroissance lente de l'ACF nous indique une chronique non stationnaire que l'on dérive donc.

In [41]:
tsdisplay(diff(diff(X,12)))

Après ces deux dérivations, la chronique semble bien stationnaire: on peut même identifier aisément une combinaison MA(1) + SMA(1).

In [42]:
dX <- diff(diff(X,12))
dY <- diff(diff(Y,12))
aTSA::adf.test(dX)
aTSA::adf.test(dY)
Augmented Dickey-Fuller Test 
alternative: stationary 
 
Type 1: no drift no trend 
     lag   ADF p.value
[1,]   0 -32.8    0.01
[2,]   1 -25.0    0.01
[3,]   2 -18.8    0.01
[4,]   3 -14.5    0.01
[5,]   4 -12.2    0.01
[6,]   5 -10.2    0.01
Type 2: with drift no trend 
     lag   ADF p.value
[1,]   0 -32.8    0.01
[2,]   1 -24.9    0.01
[3,]   2 -18.8    0.01
[4,]   3 -14.5    0.01
[5,]   4 -12.2    0.01
[6,]   5 -10.2    0.01
Type 3: with drift and trend 
     lag   ADF p.value
[1,]   0 -32.7    0.01
[2,]   1 -24.9    0.01
[3,]   2 -18.8    0.01
[4,]   3 -14.4    0.01
[5,]   4 -12.2    0.01
[6,]   5 -10.2    0.01
---- 
Note: in fact, p.value = 0.01 means p.value <= 0.01 
Augmented Dickey-Fuller Test 
alternative: stationary 
 
Type 1: no drift no trend 
     lag    ADF p.value
[1,]   0 -14.18    0.01
[2,]   1 -11.18    0.01
[3,]   2  -9.68    0.01
[4,]   3  -9.45    0.01
[5,]   4  -8.33    0.01
[6,]   5  -7.31    0.01
Type 2: with drift no trend 
     lag    ADF p.value
[1,]   0 -14.17    0.01
[2,]   1 -11.17    0.01
[3,]   2  -9.67    0.01
[4,]   3  -9.44    0.01
[5,]   4  -8.32    0.01
[6,]   5  -7.30    0.01
Type 3: with drift and trend 
     lag    ADF p.value
[1,]   0 -14.16    0.01
[2,]   1 -11.16    0.01
[3,]   2  -9.67    0.01
[4,]   3  -9.43    0.01
[5,]   4  -8.31    0.01
[6,]   5  -7.30    0.01
---- 
Note: in fact, p.value = 0.01 means p.value <= 0.01 

Après les deux dérivations, les deux chroniques semblent effectivement stationnaires si on se fie aux tests ADF.

On cherche le filtre de blanchiment de X (qu'il faut différentier pour stationnariser comme discuté précédemment): on teste MA(1)+SMA(1)/

In [43]:
arima_X <- Arima(X,order=c(0,1,1),seasonal=c(0,1,1))
tsdiag(arima_X)
coeftest(arima_X)
z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.674462   0.038011 -17.744 < 2.2e-16 ***
sma1 -0.999978   0.034817 -28.721 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ce modèle convient en effet. Rappelons qu'on ne cherche pas un modèle parfait de $X$ (on voit ici que les tests de Portmanteau sont un peu limite): le but est de se mettre dans les conditions du théorème permettant d'interprêter le corrélogramme croisé. Ce filtre n'interviendra pas dans le modèle final.

On applique le filtre blanchissant X à Y, et on calcule les corrélations croisées (c'est le ccf entre chi et upsilon avec les notations du cours).

In [44]:
arima_Yf <- Arima(Y,model=arima_X)
ccf(residuals(arima_X),residuals(arima_Yf),lag.max=50)

Conclusion que l'on tire de ce graphe:

Cela a un sens de chercher un modèle à fonction de transfert: recruit semble corrélé aux valeurs passées de SOI. Ici on a globalement une décroissance géométrique à partir du décalage 5, donc une fonction de transfert avec dénominateur de degré 1 (on pourrait essayer degré 2 au dénominateur, pour tenir compte des corrélations positives donnant un aspect vaguement "sinusoidal")

Les autres décalages sont (assez) proches des bornes des intervalles de confiance.

Au passage, si vous faites ̀ccf(X,Y,lag.max=50) vous allez voir de nombreuses corrélations fallacieuses.


On effectue une régression dynamique entre chroniques stationnaires: on régresse dY contre dX pour les décalages compris entre 5 et 28.

In [45]:
out_lm <- dynlm(dY ~ L(dX,-(5:28)))
summary(out_lm)
tsdisplay(residuals(out_lm))
Time series regression with "ts" data:
Start = 1951(2), End = 1985(5)

Call:
dynlm(formula = dY ~ L(dX, -(5:28)))

Residuals:
    Min      1Q  Median      3Q     Max 
-57.153  -7.824   1.434   8.193  31.856 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -0.02759    0.65209  -0.042 0.966268    
L(dX, -(5:28))-5   2.27563    2.12755   1.070 0.285467    
L(dX, -(5:28))-6   8.29453    2.47846   3.347 0.000898 ***
L(dX, -(5:28))-7  14.67181    2.63782   5.562 4.98e-08 ***
L(dX, -(5:28))-8   2.83835    2.69205   1.054 0.292382    
L(dX, -(5:28))-9   2.66997    2.71138   0.985 0.325373    
L(dX, -(5:28))-10 -1.24901    2.70842  -0.461 0.644944    
L(dX, -(5:28))-11 -2.74481    2.69176  -1.020 0.308502    
L(dX, -(5:28))-12 -1.80479    2.67812  -0.674 0.500774    
L(dX, -(5:28))-13  0.14053    2.66519   0.053 0.957975    
L(dX, -(5:28))-14  2.03605    2.63735   0.772 0.440582    
L(dX, -(5:28))-15 -2.04160    2.57487  -0.793 0.428326    
L(dX, -(5:28))-16 -1.22056    2.51424  -0.485 0.627625    
L(dX, -(5:28))-17  0.66499    2.49310   0.267 0.789817    
L(dX, -(5:28))-18  6.38920    2.56657   2.489 0.013215 *  
L(dX, -(5:28))-19 11.88751    2.62930   4.521 8.19e-06 ***
L(dX, -(5:28))-20  6.86317    2.67613   2.565 0.010706 *  
L(dX, -(5:28))-21  6.08128    2.70514   2.248 0.025136 *  
L(dX, -(5:28))-22 -1.87056    2.71383  -0.689 0.491068    
L(dX, -(5:28))-23 -1.69521    2.70864  -0.626 0.531780    
L(dX, -(5:28))-24  0.65945    2.70027   0.244 0.807192    
L(dX, -(5:28))-25  5.44505    2.69574   2.020 0.044085 *  
L(dX, -(5:28))-26  4.99535    2.65641   1.880 0.060792 .  
L(dX, -(5:28))-27  4.12098    2.48832   1.656 0.098506 .  
L(dX, -(5:28))-28  3.34263    2.13081   1.569 0.117531    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 13.23 on 387 degrees of freedom
Multiple R-squared:  0.1592,	Adjusted R-squared:  0.1071 
F-statistic: 3.054 on 24 and 387 DF,  p-value: 3.097e-06

De nombreux paramètres semblent non significatifs, néanmoins on ne peut pas interpréter ces tests tant qu'on ne s'est pas assurés que les résidus forment un b.b.g.

Ce n'est manifestement pas le cas ici: on identifierait plutôt AR(1)+SMA(1) sur les résidus.


Néanmoins, plutôt qu'estimer 22 coefficients, on voudrait fixer une fonction de transfert en $\frac{B^5}{1-\alpha B}$ ($B^5$ pour le retard), comme évoqué dans la discussion du ccf.

In [46]:
output_XY <- arimax(Y,order=c(1,1,0),seasonal=c(0,1,1),transfer=list(c(1,5)),fixed=c(NA,NA,NA,0,0,0,0,0,NA),xtransf=X,method="ML")
summary(output_XY)
Call:
arimax(x = Y, order = c(1, 1, 0), seasonal = c(0, 1, 1), fixed = c(NA, NA, NA, 
    0, 0, 0, 0, 0, NA), method = "ML", xtransf = X, transfer = list(c(1, 5)))

Coefficients:
         ar1     sma1  T1-AR1  T1-MA0  T1-MA1  T1-MA2  T1-MA3  T1-MA4    T1-MA5
      0.3296  -0.9985  0.8378       0       0       0       0       0  -19.7263
s.e.  0.0453   0.0319  0.0360       0       0       0       0       0    1.3493

sigma^2 estimated as 54.27:  log likelihood = -1507.37,  aic = 3022.74

Training set error measures:
                       ME     RMSE      MAE       MPE     MAPE      MASE
Training set -0.005970017 7.259283 4.659853 0.6757152 11.22912 0.5825106
                   ACF1
Training set 0.02807246

Remarque importante : il y a un souci d'estimation, avec le coef sma1 très proche de -1 !!! (et qui compense donc la dérivation saisonnière)

(sans l'option method="ML", l'optimisation ne converge pas du tout.

On touche là une limite de l'algorithme d'optimisation de R (fonction optim) qui a l'air assez critiquée...

(on peut très bien évaluer ce modèle sous SAS par exemple)


On continue tout de même pour la beauté de l'exercice pédagogique:

In [47]:
output_XY <- arimax(Y,order=c(1,1,0),seasonal=c(0,1,1),transfer=list(c(1,5)),fixed=c(NA,NA,NA,0,0,0,0,0,NA),xtransf=X,method="ML")
summary(output_XY)

tsdiag(output_XY) 

hist(residuals(output_XY), freq = F, col = "grey")
curve(dnorm(x, mean = mean(residuals(output_XY),na.rm=T), sd = sd(residuals(output_XY),na.rm=T)), col = 2, add = TRUE)
Call:
arimax(x = Y, order = c(1, 1, 0), seasonal = c(0, 1, 1), fixed = c(NA, NA, NA, 
    0, 0, 0, 0, 0, NA), method = "ML", xtransf = X, transfer = list(c(1, 5)))

Coefficients:
         ar1     sma1  T1-AR1  T1-MA0  T1-MA1  T1-MA2  T1-MA3  T1-MA4    T1-MA5
      0.3296  -0.9985  0.8378       0       0       0       0       0  -19.7263
s.e.  0.0453   0.0319  0.0360       0       0       0       0       0    1.3493

sigma^2 estimated as 54.27:  log likelihood = -1507.37,  aic = 3022.74

Training set error measures:
                       ME     RMSE      MAE       MPE     MAPE      MASE
Training set -0.005970017 7.259283 4.659853 0.6757152 11.22912 0.5825106
                   ACF1
Training set 0.02807246

Il reste une corrélation significative pour un décalage de 4, et on voit des valeurs négatives aberrantes, mais il y en a peu.

Remember: "all models are wrong..."

In [48]:
summary(output_XY)
Call:
arimax(x = Y, order = c(1, 1, 0), seasonal = c(0, 1, 1), fixed = c(NA, NA, NA, 
    0, 0, 0, 0, 0, NA), method = "ML", xtransf = X, transfer = list(c(1, 5)))

Coefficients:
         ar1     sma1  T1-AR1  T1-MA0  T1-MA1  T1-MA2  T1-MA3  T1-MA4    T1-MA5
      0.3296  -0.9985  0.8378       0       0       0       0       0  -19.7263
s.e.  0.0453   0.0319  0.0360       0       0       0       0       0    1.3493

sigma^2 estimated as 54.27:  log likelihood = -1507.37,  aic = 3022.74

Training set error measures:
                       ME     RMSE      MAE       MPE     MAPE      MASE
Training set -0.005970017 7.259283 4.659853 0.6757152 11.22912 0.5825106
                   ACF1
Training set 0.02807246

Le modèle s'écrit:

$$ \text{Recruit}_t = \frac{-19.7263}{1-0.8378B} \text{SOI}_{t-5}+ u_t $$

avec $$ (1-B)(1-B^{12})u_t = \frac{1-0.99847B^{12}}{1-0.3296B} \varepsilon_t $$

On suspecte un pb d'estimation vu que la dérivation saisonnière compense quasiment le SMA(1)...

En tout cas, on a mis en évidence un retard de 5 mois dans l'influence de SOI sur RECRUIT.


Remarque: sous SAS, on obtenait le modèle: $$ \text{Recruit}_t = \frac{-16.664}{1-0.453B} \text{SOI}_{t-5} + u_t $$ avec $$(1-B)(1-B^{12})u_t = \frac{1-0.889B^{12}}{1-0.401B} \varepsilon_t$$ qui n'est finalement pas tellement éloigné.

In [ ]: