library(lmtest)
library(TSA)
library(aTSA)
library(lmtest)
library(forecast)
options(repr.plot.width=12, repr.plot.height=10) # pour redimensionner les graphiques
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
Représentation graphique:
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:
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.
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).
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:
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:
aTSA::adf.test(residuals(arima_sncf))
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.
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) )
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.
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).
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)
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.
summary(arima_sncf)
On note RMSE= 120.5 AIC=2547.95 BIC=2561.21
, le modèle s'écrit:
On examine à présent l'hypothèse SAR(1):
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:
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)
summary(arima_sncf)
On note: RMSE= 121.99 AIC=2552.38 BIC=2565.63
et on écrit le modèle:
(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:
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()
On peut restreindre l'affichage aux dernières années:
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:
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).
busban <- read.table(file="busban.txt",sep='\t',header=TRUE)
busban <- ts(busban,start=c(1984,1),frequency=12)
busban
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:
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).
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.
plot(diff(busban,12))
aTSA::adf.test(diff(busban,12))
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.
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):
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):
# 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:
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):
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):
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)
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)
Le modèle en SMA(1) précédent est légèrement meilleur (sigma^2, AIC) on y revient:
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)
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.
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()
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)
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.
summary(arima_buspulse)
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:
coefci(arima_buspulse)
Avec une intervention décroissante au cours du temps:
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)
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!
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...):
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$.
tsdisplay(X)
Belle périodicité de l'ACF: il y a saisonnalité sur une période de 12 mois.
tsdisplay(diff(X,12))
La décroissance lente de l'ACF nous indique une chronique non stationnaire que l'on dérive donc.
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).
dX <- diff(diff(X,12))
dY <- diff(diff(Y,12))
aTSA::adf.test(dX)
aTSA::adf.test(dY)
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)/
arima_X <- Arima(X,order=c(0,1,1),seasonal=c(0,1,1))
tsdiag(arima_X)
coeftest(arima_X)
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).
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.
out_lm <- dynlm(dY ~ L(dX,-(5:28)))
summary(out_lm)
tsdisplay(residuals(out_lm))
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.
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)
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:
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)
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..."
summary(output_XY)
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é.