Deux chroniques sont disponibles:
Question: Peut-on établir un modèle permettant d'expliquer la chronique $Y$ à partir de la chronique $X$?
On commence par charger les bibliothèques (à installer au préalable si nécessaire par install.packages("library_name")
) et les données.
library(forecast)
library(lmtest)
library(TSA)
library(aTSA)
library(dynlm)
options(repr.plot.width=12, repr.plot.height=10) # pour redimensionner les graphiques
seriesJ <- read.table(file="seriesJ.txt",sep='\t',header=T)
X <- ts(seriesJ[,1])
Y <- ts(seriesJ[,2])
par(mar=c(5,4,4,4))
plot(X, col="blue", ylab="",type="o",pch=20)
mtext("Input Gas",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) # trace l'axe à droite
mtext("Output CO2",side=4,line=2,col="red")
grid()
Il ne semble pas complètement évident que les évolution de $X_t$ ( input gas ) gouvernent $Y_t$ ( output CO2 )...
Nous allons essayer d'identifier un modèle à fonction de transfert à l'aide du corrélogramme croisé.
Pour ce faire, il faut commencer par identifier un filtre blanchissant $X_t$, en ayant vérifié au préalable que $X_t$ est bien stationnaire.
aTSA::adf.test(X)
Les tests de Dickey-Fuller suggèrent que la série $X_t$ est bien stationnaire: on rejette toujours la stationnarité (sauf dans un cas du type 3, qui n'est de toute façon pas réaliste ici).
Ici, il n'y a donc pas besoin de dériver pour stationnariser les chroniques.
On examine à présent ACF et PACF de la série $X_t$ pour identifier un modèle ARMA.
tsdisplay(X)
Décroissance (certes, peu rapide...) dans l'ACF et corrélations partielles significatives pour les décalages $h=1, 2, 3$ (et 4?) dans le PACF suggèrent un processus AR(3) (ou AR(4)?).
N'oublions pas que le filtre de blanchiment n'intervient que dans l'identification de la fonction de transfert par l'intermédiaire du corrélogramme croisé, et pas dans le modèle final. On peut donc se satisfaire d'un filtre grossier (toutes les hypothèses d'une modélisation par processus ARMA n'ont pas forcément à être satisfaites).
arima_X <- Arima(X,order=c(4,0,0))
tsdisplay(residuals(arima_X))
On constate qu'un modèle AR(4) fournit des corrélations plus faibles qu'un AR(3) (on montre ici AR(4): vous pouvez essayer AR(3)): on choisit donc le modèle AR(4) pour construire le filtre de blanchiment. Cela ne changera rien à l'identification avec le corrélogramme croisé.
Le modèle AR(4) estimé est:
summary(arima_X)
On a donc, avec les notations du cours: $$(1-1.927B+1.197B^2-0.099B^3-0.1215B^4)(X_t+0.059) = \chi_t$$ où $\chi_t$ est un bruit blanc gaussien.
Le filtre blanchissant $X_t$ (c'est-à-dire transformant $X_t$ en $\chi_t$) est donc: $$ (1-1.927B+1.197B^2-0.099B^3-0.1215B^4)$$
Il faut vérifier également que la chronique $Y_t$ filtrée par le filtre AR4 précédent est un processus stationnaire.
Il faut passer model=arima_X
en paramètre à la fonction Arima
pour qu'elle utilise un modèle existant au lieu d'estimer des paramètres.
arima_Yf <- Arima(Y,model=arima_X)
tsdisplay(residuals(arima_Yf))
Effectivement, l'allure du graphe des résidus (c'est-à-dire la chronique $Y_t$ filtrée par le filtre AR4 précédent) ainsi que les corrélogrammes sont cohérents avec une hypothèse de stationnarité.
Nous sommes donc dans les conditions d'application du théorème du cours qui nous dit que les coefficients du développement polynomial de la fonction de transfert sont proportionnels aux coefficients du corrélogramme croisé entre les chroniques $X_t$ et $Y_t$ auxquelles on a appliqué le filtre blanchissant $X_t$.
Avec les notations du cours, $\chi_t$ est residuals(arima_X)
et $\Upsilon_t$ est residuals(arima_Yf)
.
ccf(residuals(arima_X),residuals(arima_Yf))
grid()
On constate que les corrélations croisées significatives ne sont observées que pour des décalages négatifs. Ce constat est cohérent avec l'hypothèse de causalité: $Y_t$ s'explique par les valeurs passées de $X_t$ (attention à la remarque du cours sur les notations utilisées dans ccf
)
Par ailleurs, la première corrélation croisée significative est observée pour un décalage de 3. Le retard est donc $b=3$.
On conclut de ce graphe: $Y_t$ est significativement corrélée avec $X_{t-3}, X_{t-4}, X_{t-5}, X_{t-6}, X_{t-7}$ (de manière limite pour ce dernier).
Nous calculons une régression de $Y$ contre ces variables, en utilisant dynlm
de la manière suivante:
out_lm <- dynlm(Y ~ L(X,-3)+L(X,-4)+L(X,-5)+L(X,-6)+L(X,-7))
summary(out_lm)
tsdisplay(residuals(out_lm))
Le modèle estimé est: $$Y_t = 53.34 - 2.77 X_{t-3} + 3.22 X_{t-4} -0.75 X_{t-5} -0.83 X_{t-6} + 0.07X_{t-7} + R_t$$
Les corrélogrammes nous indiquent que les résidus $R_t$ sont corrélés (on ne peut donc pas interpréter les tests de significativités affichés par dynlm
), et nous permettent d'identifier un modèle AR(2). (En effet, les corrélations partielles pour des décalages 3 et 5 sont à la limite de la significativité mais on cherche un modèle simple.)
On utilise à présent la fonction arimax
, qui nous permet d'indiquer une fonction de transfert et un modèle ARMA sur les résidus du modèle. Attention, dans le paramètre transfer
de cette fonction, le premier coefficient est le degré du dénominateur (0 ici), et le second coefficient le degré du numérateur.
output_XY <- arimax(Y,order=c(2,0,0),transfer=list(c(0,7)),xtransf=X)
summary(output_XY)
tsdisplay(residuals(output_XY))
# pour une raison un peu obscure, tsdiag n'est pas vraiment utilisable...
Le modèle estimé ici s'écrit:
$$Y_t = 53.37 -0.076 X_t + 0.065 X_{t-1} -0.047 X_{t-2} -0.5462 X_{t-3} -0.6586 X_{t-4} -0.852 X_{t-5} -0.486 X_{t-6} -0.366 X_{t-7} + \frac{1}{1-1.540B+0.633B^2} \varepsilon_t $$Les résidus $\varepsilon_t$ du modèle semblent correspondre à un bruit blanc (pas de corrélation visible dans l'ACF), on peut donc à présent interpréter les $Z$-tests:
coeftest(output_XY)
Tous les coefficients sont significatifs, sauf les paramètres devant $X_t$, $X_{t-1}$, et $X_{t-2}$, ce qui est cohérent avec l'allure du corrélogramme croisé discuté précédemment. En effet, nous avions identifié un retard $b=3$.
Nous forçons les paramètres non-significatifs à la valeur 0 par l'intermédiaire du paramètre fixed
: dans la liste, les coefficients sont dans l'ordre de la sortie de summary(output_XY)
.
output_XY <- arimax(Y,order=c(2,0,0),transfer=list(c(0,7)),fixed=c(NA,NA,NA,0,0,0,NA,NA,NA,NA,NA),xtransf=X)
summary(output_XY)
Le modèle est à présent:
$$ Y_t = 53.38 -0.556 X_{t-3} -0.645 X_{t-4} -0.860 X_{t-5} -0.484 X_{t-6} -0.363 X_{t-7} + \frac{1}{1-1.540B+0.630B^2} \varepsilon_t $$Au passage, on peut remarquer que les coefficients sont assez différents de ceux estimés plus haut par dynlm
: passer à un modèle AR(2) sur les résidus n'est pas anecdotique pour la modélisation.
On peut noter la variance ($\sigma^2$) estimée à 0.05823, la log-vraisemblance à -0.82, et le critère d'information AIC à 17.64.
tsdisplay(residuals(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)
Les résidus de ce modèle ne sont pas corrélés temporellement (cf. ACF), leur répartition semble gaussienne (cf. histogramme), on peut donc interprêter les $Z$-tests suivants:
coeftest(output_XY)
Tous les coefficients sont significatifs: le modèle est valide.
(ne pas passer trop de temps sur cette partie: l'objectif est surtout de vous fournir la syntaxe des commandes)
Ce modèle nécessite l'estimation de 5 coefficients dans la fonction de transfert, qui est ici polynomiale. On peut se demander si moins de paramètres ne pourraient pas être utilisés avec une fonction de transfert qui serait le rapport de deux polynomes de degré faible.
L'allure du corrélogramme croisé ne correspond pas vraiment aux cas typiques du cours. En particulier, la fonction de transfert ne peut pas être du type $1/(1-\alpha B)$ car on observerait des coefficients du corrélogramme croisé en "décroissance géométrique".
Dans les cellules suivantes, nous allons donc essayer différentes combinaisons des degrés des numérateurs et dénominateurs.
Attention, ici le décalage est $b=3$: il faut en tenir compte dans les paramètres passés à arimax
. Par exemple, si on souhaite un numérateur de la fonction de transfert de degré 1 avec un retard de 3, le deuxième argument de l'argument transfer
doit être 3+1 = 4, et il faut imposer trois coefficients à 0 dans le modèle (on fait comme si on voulait un numérateur de degré 4, avec les coefficients des termes de degrés 0, 1, et 2 nuls).
# fonction de transfer à numérateur et dénominateur degré 1:
output_XY <- arimax(Y,order=c(2,0,0),transfer=list(c(1,4)),fixed=c(NA,NA,NA,NA,0,0,0,NA,NA),xtransf=X)
summary(output_XY)
tsdisplay(residuals(output_XY)) # des corrélations
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)
coeftest(output_XY)
# sigma^2 estimated as 0.06113: log likelihood = -7.82, aic = 27.64
# Y= 53.351 - (0.4551+0.6320B)/(1-0.6740B) X_{t-3} + 1/(1-1.5177B+0.6267B^2) eps_t
# fonction de transfer à numérateur degré 2 et dénominateur degré 1:
output_XY <- arimax(Y,order=c(2,0,0),transfer=list(c(1,5)),fixed=c(NA,NA,NA,NA,0,0,0,NA,NA,NA),xtransf=X)
summary(output_XY)
tsdisplay(residuals(output_XY)) # des corrélations
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)
coeftest(output_XY)
# sigma^2 estimated as 0.0571: log likelihood = 2.08, aic = 9.83
# Y= 53.362 - (0.53096+0.38013B+0.51801B^2)/(1-0.54903B) X_{t-3} + 1/(1-1.5272B+0.6288B^2) eps_t
# fonction de transfer à numérateur degré 1 et dénominateur degré 2:
output_XY <- arimax(Y,order=c(2,0,0),transfer=list(c(2,4)),fixed=c(NA,NA,NA,NA,NA,0,0,0,NA,NA),xtransf=X)
summary(output_XY)
tsdisplay(residuals(output_XY)) # mieux pour les corrélations
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)
coeftest(output_XY)
# sigma^2 estimated as 0.05849: log likelihood = -1.41, aic = 16.81
# Y= 53.367 - (0.4826+0.3409B)/(1-1.031B+0.2928B^2) X_{t-3} + 1/(1-1.5215B+0.6232B^2) eps_t
# remarque si numérateur de dégré 2 et dénominateur de degré 2, on constate que le coef de degré 2 au dénominateur n'est pas significatif (à faire)
Ces trois modèles sont valides.
Parmi tous les modèles testés, celui de variance minimale (il explique le mieux les données) et de critère AIC minimal (il réalise le meilleur compromis entre complexité et explication des données) est celui donné par:
output_XY=arimax(Y,order=c(2,0,0),transfer=list(c(1,5)),fixed=c(NA,NA,NA,NA,0,0,0,NA,NA,NA),xtransf=X)
Il s'écrit:
$$ Y_t = 53.362 - \frac{0.53096+0.38013B+0.51801B^2}{1-0.54903B} X_{t-3} + \frac{1}{1-1.5272B+0.6288B^2} \varepsilon_t $$avec une variance de $\sigma^2=0.0571$ pour le bruit blanc gaussien $\varepsilon_t$.
C'est le modèle final que nous retenons de notre étude.