Etude de la chronique ozone¶
Il s'agit d'une mesure de la pollution moyenne mensuelle dans la ville de Los Angeles de 1955 à 1972.
Une loi imposant le pot d'échappement catalytique a été appliquée à partir de janvier 1960.
Question: cette loi a-t-elle un impact significatif ? Si oui, pour quelle diminution de la pollution ?
Nous allons représenter graphiquement les données, puis modéliser la chronique ozone sans tenir compte de la loi (modélisation de référence), puis enfin incorporer la loi dans le modèle à l'aide d'une intervention. L'objectif est de déterminer quel est le modèle qui explique le mieux les données.
1. Représentations graphiques¶
On commence par charger les bibliothèques utiles et la chronique.
library(forecast)
library(lmtest)
library(aTSA)
options(repr.plot.width=12, repr.plot.height=10) # pour redimensionner les graphiques
ozone <- read.table("ozone.txt",header=T)
ozone <- ts(ozone,start=c(1955,1),frequency=12)
plot(ozone, type="o", pch =20, main="graphe de la chronique Ozone")
grid()
La chronique n'est pas stationnaire, on va être amené à la différencier (ou dériver). Nous n'introduirons pas de transformation de la chronique comme un passage au logarithme par exemple.
Pour commencer, examinons le graphique de la différenciation saisonnière:
plot(diff(ozone,12), ylab="(1-B^12)Ozone", main="graphe de la chronique (1-B^12)Ozone")
grid()
Le graphe ressemble à celui d'une marche aléatoire: il faudra sans doute encore différencier pour arriver à une chronique stationnaire. Néanmoins un processus AR(1) pourrait fournir une représentation similaire.
C'est d'ailleurs ce que semble suggérer les tests augmentés de Dickey-Fuller: dans tous les cas on peut rejeter l'hypothèse de non-stationnarité de la chronique $(1-B^{12})\text{Ozone}_t$ (et donc accepter sa stationnarité).
adf.test(diff(ozone,12))
plot(diff(diff(ozone,12)), ylab="(1-B)(1-B^12)Ozone", main="graphe de la chronique (1-B)(1-B^12)Ozone")
grid()
Cette chronique différenciée deux fois a l'air stationnaire. Dans la méthode de Box et Jenkins il faut s'attendre à devoir appliquer une différenciation saisonnière, et peut-être aussi une différenciation.
Rappelons que l'ordre des différenciations est indifférent, car $(1-B)(1-B^{12}) = (1-B^{12})(1-B)$.
Examinons une décomposition par lissage de la chronique ozone:
plot(decompose(ozone,type="additive"))
La rupture de tendance en 1960 n'est pas claire: il n'est pas évident que le modèle d'intervention sera le plus adapté. La tendance est globalement à une baisse de la pollution à l'ozone au cours du temps.
2. Modèle sans intervention¶
(modélisation de référence)
On commence par examiner les corrélogrammes de la chronique initiale:
tsdisplay(ozone)
Le corrélogramme (ACF) est périodique: il faut utiliser une différenciation saisonnière, comme le suggéraient nos graphes initiaux.
tsdisplay(diff(ozone,12))
L'ACF présente une décroissance relativement lente: on peut dériver. On peut aussi préférer l'identification alternative: décroissance des corrélations (ACF), corrélation significative pour un décalage de 1 dans le PACF, d'où l'identification d'un processus AR(1) sur les décalages "de court-terme".
2a. On considère $(1-B)(1-B^{12})\text{Ozone}_t$ stationnaire¶
On commence par le modèle avec différenciation. Remarquons que les tests de Dickey-Fuller nous indiquaient la stationnarité de $(1-B^{12})\text{Ozone}_t$: une dérivation supplémentaire peut être jugée inutile. Nous étudions tout de même ce cas pour des raisons pédagogiques, mais on peut aussi directement passer au paragraphe 2b.
tsdisplay(diff(diff(ozone,12)))
On remarque que l'ACF pour un décalage de 1 est négatif (presque -0.4), peut-être cette dérivation ne s'imposait-elle pas? (comme le suggéraient les graphiques initiaux et les tests de Dickey-Fuller)
On identifie pour la composante court-terme un processus MA(1) (MA(2)?). En effet, il y a globalement décroissance des corrélations partielles, et des corrélations significatives pour des décalages de 1 (2?). Rappelons qu'on n'inclut pas de constante dans le modèle lorsqu'on applique deux dérivations.
arima_ozone <- Arima(ozone,order=c(0,1,1),seasonal=c(0,1,0))
# remarque: la moyenne est nulle par défaut
tsdisplay(residuals(arima_ozone))
Il n'y a plus de corrélations à court-terme. Il n'est pas nécessaire d'envisager un processus MA(2). On identifie maintenant un processus pour la composante saisonnière, en examinant les décalages multiples de 12. Remarquons que les corrélations périodiques sont plus lisibles une fois les corrélations court-termes traitées. Comme l'ACF est significatif pour un décalage de 12, non significatif pour un décalage de 24, et que le PACF indique des corrélations partielles significatives et décroissantes entre les décalages 12 et 24, on identifie un processus SMA(1) pour la composante saisonnière.
arima_ozone <- Arima(ozone,order=c(0,1,1),seasonal=c(0,1,1))
tsdisplay(residuals(arima_ozone))
tsdiag(arima_ozone,gof.lag=15) # il faut examiner les décalages au moins jusque la valeur de la période
Il n'y a plus de corrélations significatives, on voit que les $p$-valeurs des tests de Portmanteau (Ljung-Box statistic) sont tout juste au dessus du risque de 5%: on ne peut pas rejeter l'hypothèse de résidus formant un bruit blanc gaussien.
hist(residuals(arima_ozone), freq = F, col = "grey")
curve(dnorm(x, mean = mean(residuals(arima_ozone)), sd = sd(residuals(arima_ozone))), col = 2, add = TRUE)
L'histogramme des résidus, superposé à la courbe gaussienne dont la moyenne et l'écart-type sont calculés à partir des données, montre une légère asymétrie dans la distribution des résidus.
Maintenant que l'hypothèse de bruit blanc gaussien pour les résidus du modèle est validée, on peut tester la significativité des coefficients:
coeftest(arima_ozone)
Les coefficients sont tous significativement non nuls: dans tous les cas la $p$-valeur est très faible, on peut donc rejeter l'hypothèse nulle ${\cal H}_0$ (= "le coefficient est nul").
summary(arima_ozone)
Le modèle final s'écrit:
$$ (1-B)(1-B^{12}) Ozone_t = (1-0.7674B)(1-0.7579B^{12}) \varepsilon_t$$avec un écart-type des résidus estimé à 8.115868 (RMSE).
2b. On considère $(1-B^{12})\;\text{Ozone}_t$ stationnaire¶
Si on décide de modéliser sans dériver initialement (on utilise tout de même une différenciation saisonnière) mais en utilisant un processus AR(1) pour modéliser les corrélations court-termes:
arima_ozone <- Arima(ozone,order=c(1,0,0),seasonal=c(0,1,0),include.constant=T)
tsdisplay(residuals(arima_ozone))
Remarquons qu'on inclut la moyenne dans le modèle (include.constant=T), ce qui n'est pas le comportement par défaut de Arima. On verra par la suite si la moyenne est significative.
Les corrélations à court-terme ne sont plus significatives.
Si on examine les corrélations pour les décalages multiple de 12, on identifie un processus SMA(1) (car décroissance pour PACF(12)-PACF(24), et ACF(12) significatif / ACF(24) non-significatif).
arima_ozone <- Arima(ozone,order=c(1,0,0),seasonal=c(0,1,1),include.constant=T)
tsdisplay(residuals(arima_ozone))
tsdiag(arima_ozone,gof.lag=15)
Les tests de Portmanteau ne permettent pas de rejeter l'hypothèse de résidus issus d'un bruit blanc gaussien.
hist(residuals(arima_ozone), freq = F, col = "grey")
curve(dnorm(x, mean = mean(residuals(arima_ozone)), sd = sd(residuals(arima_ozone))), col = 2, add = TRUE)
On observe toujours une légère asymétrie.
coeftest(arima_ozone)
summary(arima_ozone)
Tous les coefficients sont significatifs. Le modèle identifié est:
$$ (1-0.364B)((1-B^{12}) \,\text{Ozone}_t + 0.1254) = (1-0.7332B^{12}) \varepsilon_t$$avec un écart-type des résidus estimé à 7.94899.
Si on compare les deux modèles identifiés, on voit que le second modèle induit une composante aléatoire d'amplitude plus faible (RMSE plus petit). Il explique donc mieux les données. Il a également une valeur du critère AIC (légèrement) plus faible, et réalise donc un meilleur compromis entre explication des données et complexité du modèle.
On retient donc ce dernier modèle, de type $\text{SARIMA}(1,0,0)(0,1,1)_{12}$, pour modéliser les données sans tenir compte de la mise en oeuvre de la loi.
3. Modèle avec intervention¶
On incorpore à présent une fonction d'intervention "en échelon".
On commence par créer une nouvelle série chronologique, appelée loi, qui vaut 1 à partir de janvier 1960, et 0 avant.
loi <- ts(as.numeric(time(ozone)>=1960),start=c(1955,1),frequency=12)
plot(ozone, type="o", pch =20, main="Chronique Ozone et intervention")
lines(loi*80,col=2)
grid()
On commence par estimer les coefficients $\mu$ et $\alpha$ d'un modèle du type: $$\text{Ozone}_t = \mu + \alpha \;\text{loi}_t + R_t$$ par régression au sens des moindres carrés.
fitOLS <- lm(ozone~loi);
summary(fitOLS)
tsdisplay(residuals(fitOLS))
Les résidus $R_t$ ne forment pas un bruit blanc gaussien (on observe une corrélation temporelles), ce à quoi on pouvait se douter vue la composante saisonnière. Rappelons que les différents tests statistiques définis dans le cadre de la régression sont basés sur l'hypothèse de résidus indépendants et identiquement distribués selon une loi normale.
Ici, nous allons plutôt déterminer les coefficients d'un modèle: $$(1-B^{12})\;\text{Ozone}_t = \mu + \alpha\; (1-B^{12})\;\text{loi}_t + R_t$$
fitOLS <- lm(diff(ozone,12)~diff(loi,12))
summary(fitOLS)
tsdisplay(residuals(fitOLS))
Notons que le cofficient $\mu$ est déclaré non significatif par la fonction lm. Néanmoins ce test n'a de sens que dans les conditions rappelées précédemment, qui ne sont pas remplies ici. De manière générale, on ne peut pas juger de la pertinence de notre modèle d'intervention tant qu'on n'a pas une modélisation fondée sur des hypothèses validées.
Les corrélogrammes montrent que les résidus de ce modèle sont corrélés. Néanmoins, ils sont sans doute à présent stationnaires. Nous allons essayer d'identifier un modèle où les résidus sont modélisés par un processus SARMA. Le couple ACF/PACF suggère un processus AR(1) ou MA(1) pour la composante court-terme.
3a. Alternative 1¶
Commençons par un processus AR(1). On va préciser à la fonction Arima que l'on cherche une régression de ozone contre loi avec des résidus formant un processus SARIMA, de la manière suivante:
arima_ozloi <- Arima(ozone,order=c(1,0,0),seasonal=c(0,1,0),xreg=loi, include.constant=T)
summary(arima_ozloi)
tsdisplay(residuals(arima_ozloi))
Dans ce modèle, les résidus n'ont plus de corrélation court-terme. Restent des corrélations saisonnières, pour lesquelles on identifie un processus SMA(1):
arima_ozloi <- Arima(ozone,order=c(1,0,0),seasonal=c(0,1,1),xreg=loi, include.constant=T)
summary(arima_ozloi)
tsdisplay(residuals(arima_ozloi))
tsdiag(arima_ozloi,gof.lag=15)
Il n'y a plus de corrélations significatives dans les résidus du modèle, ce qui est confirmé par les tests de Portmanteau. Remarquons qu'ici, il n'y a pas d'hésitation: les $p$-valeurs de ce test sont au dessus du risque.
hist(residuals(arima_ozloi), freq = F, col = "grey")
curve(dnorm(x, mean = mean(residuals(arima_ozloi)), sd = sd(residuals(arima_ozloi))), col = 2, add = TRUE)
On peut toujours noter une légère asymétrie.
coeftest(arima_ozloi)
Les coefficients sont significatifs, le modèle est validé. Les intervalles de confiance sont donnés par:
coefci(arima_ozloi)
summary(arima_ozloi)
Le modèle avec intervention s'écrit donc:
$$ (1-B^{12})\, \text{Ozone}_t = -0.074 -9.1819 (1-B^{12}) \;\text{loi}_t + \frac{1-0.7511B^{12}}{1-0.3012B} \varepsilon_t $$Autrement dit: $$ \text{Ozone}_t = -9.1819 \;\text{loi}_t + R_t$$ où $$ (1-0.3012B)((1-B^{12})\,R_t + 0.074) = (1-0.7511B^{12})\; \varepsilon_t $$
L'écart-type des résidus $\varepsilon_t$ est estimé à 7.749069.
Remarquons qu'on n'écrit pas les polynomes symboliques $1-B$ ou $1-B^{12}$ au dénominateur d'une fraction rationnelle en $B$ car ils ne sont pas inversibles.
3b. Alternative 2¶
Pour la modélisation du processus court-terme, on pouvait hésiter avec un processus MA(1):
arima_ozloi <- Arima(ozone,order=c(0,0,1),seasonal=c(0,1,0),xreg=loi, include.constant=T)
summary(arima_ozloi)
tsdisplay(residuals(arima_ozloi))
... pour lequel on identifie un processus SMA(1) dans la partie saisonnière:
arima_ozloi <- Arima(ozone,order=c(0,0,1),seasonal=c(0,1,1),xreg=loi, include.constant=T)
summary(arima_ozloi)
tsdisplay(residuals(arima_ozloi))
tsdiag(arima_ozloi,gof.lag=15)
hist(residuals(arima_ozloi), freq = F, col = "grey")
curve(dnorm(x, mean = mean(residuals(arima_ozloi)), sd = sd(residuals(arima_ozloi))), col = 2, add = TRUE)
coeftest(arima_ozloi)
Ce modèle (de type $\text{SARIMA}(0,0,1)(0,1,1)_{12}$) est aussi valide que le précédent. Il s'écrit:
$$ (1-B^{12}) \,\text{Ozone}_t = - 0.0748 -9.2251 (1-B^{12}) \;\text{loi}_t + (1+0.2828B)(1-0.7363B^{12}) \varepsilon_t $$Autrement dit: $$ \text{Ozone}_t = - 9.2251 \;\text{loi}_t + R_t$$ où $$ (1-B^{12})\;R_t = - 0.0748+ (1+0.2828B)(1-0.7363B^{12})\; \varepsilon_t $$
Néanmois, il présente des valeurs de RMSE et AIC légèrement plus grandes: on préfèrera le premier modèle de type $\text{SARIMA}(1,0,0)(0,1,1)_{12}$.
4. Conclusion sur le modèle retenu¶
Le modèle retenu finalement s'écrit: $$ \text{Ozone}_t = -9.1819 \;\text{loi}_t + R_t$$ où $$ (1-0.3012B)((1-B^{12})\,R_t + 0.074) = (1-0.7511B^{12})\; \varepsilon_t $$
Le drift négatif modélise la décroissance globale du taux d'ozone au cours du temps, ce qui est visible dans le graphique de la chronique. Néanmoins, nous superposons à cette décroissance intrinsèque à la chronique une intervention extrinsèque par l'intermédiaire de la variable loi. Le coefficient négatif ($-9.1819$) montre que l'intervention se traduit par un décrochage "vers le bas" d'amplitude $9.1819$, d'effet permanent à partir de janvier 1960. Cette valeur est à mettre en relation avec l'amplitude des fluctuations aléatoires et avec la gamme des valeurs prises par la chronique (entre 20 et 80 unités).
Ce dernier modèle (avec intervention) a une RMSE plus faible que les modèles qui ne tiennent pas compte de la loi. Les critères d'information (AIC/AICc/BIC) sont également plus faibles, signe que ce modèle réalise un meilleur compromis entre explication des données et complexité du modèle.
On peut donc conclure: notre modélisation de l'influence de la loi comme un échelon explique mieux les données qu'une modélisation n'en tenant pas compte. Selon notre modèle, l'influence de la loi a permis de diminuer significativement (au sens statistique) la pollution à l'ozone. On peut même quantifier cette influence à une diminution de $9.1819$ unités (non précisées ici), avec un intervalle de confiance de l'estimation au risque de 5% qui est $[-14.5923581 ; -3.7714197]$.