Etude de la série J de Box et Jenkins


Deux chroniques sont disponibles:

  • $X_t$: débit de gaz alimentant un appareil de chauffage
  • $Y_t$: concentration de CO2 en sortie de l'apparil

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.

In [1]:
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])
Registered S3 method overwritten by 'xts':
  method     from
  as.zoo.xts zoo 

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

Registered S3 methods overwritten by 'forecast':
  method             from    
  fitted.fracdiff    fracdiff
  residuals.fracdiff fracdiff

Loading required package: zoo


Attaching package: 'zoo'


The following objects are masked from 'package:base':

    as.Date, as.Date.numeric


Registered S3 methods overwritten by 'TSA':
  method       from    
  fitted.Arima forecast
  plot.Arima   forecast


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:forecast':

    forecast


The following object is masked from 'package:graphics':

    identify


Warning message:
"package 'dynlm' was built under R version 3.6.3"

1. Représentation graphique


On représente les deux chroniques sur le même graphique:

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


2. Corrélogramme croisé

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.

In [3]:
aTSA::adf.test(X)
Augmented Dickey-Fuller Test 
alternative: stationary 
 
Type 1: no drift no trend 
     lag   ADF p.value
[1,]   0 -2.66    0.01
[2,]   1 -7.66    0.01
[3,]   2 -4.87    0.01
[4,]   3 -4.13    0.01
[5,]   4 -3.77    0.01
[6,]   5 -4.12    0.01
Type 2: with drift no trend 
     lag   ADF p.value
[1,]   0 -2.67  0.0855
[2,]   1 -7.66  0.0100
[3,]   2 -4.88  0.0100
[4,]   3 -4.13  0.0100
[5,]   4 -3.78  0.0100
[6,]   5 -4.13  0.0100
Type 3: with drift and trend 
     lag   ADF p.value
[1,]   0 -2.82   0.229
[2,]   1 -8.14   0.010
[3,]   2 -5.19   0.010
[4,]   3 -4.41   0.010
[5,]   4 -4.05   0.010
[6,]   5 -4.44   0.010
---- 
Note: in fact, p.value = 0.01 means p.value <= 0.01 

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.

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

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

In [6]:
summary(arima_X)
Series: X 
ARIMA(4,0,0) with non-zero mean 

Coefficients:
         ar1      ar2     ar3     ar4     mean
      1.9271  -1.1972  0.0986  0.1215  -0.0593
s.e.  0.0575   0.1259  0.1258  0.0574   0.2123

sigma^2 estimated as 0.03536:  log likelihood=74.79
AIC=-137.59   AICc=-137.3   BIC=-115.45

Training set error measures:
                        ME      RMSE       MAE MPE MAPE      MASE          ACF1
Training set -3.223503e-05 0.1864469 0.1306102 NaN  Inf 0.5134663 -0.0006399262

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.

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

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


3. Premier modèle Arimax (fonction de transfert polynomiale en B)

Nous calculons une régression de $Y$ contre ces variables, en utilisant dynlm de la manière suivante:

In [9]:
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))
Time series regression with "ts" data:
Start = 1, End = 289

Call:
dynlm(formula = Y ~ L(X, -3) + L(X, -4) + L(X, -5) + L(X, -6) + 
    L(X, -7))

Residuals:
   Min     1Q Median     3Q    Max 
-7.184 -2.159 -0.088  2.240  6.871 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 53.33944    0.17758 300.372  < 2e-16 ***
L(X, -3)    -2.77294    0.94067  -2.948  0.00347 ** 
L(X, -4)     3.22131    2.04464   1.575  0.11626    
L(X, -5)    -0.74719    2.33555  -0.320  0.74927    
L(X, -6)    -0.82775    2.04335  -0.405  0.68571    
L(X, -7)     0.07172    0.94027   0.076  0.93925    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.012 on 283 degrees of freedom
Multiple R-squared:  0.1094,	Adjusted R-squared:  0.09369 
F-statistic: 6.954 on 5 and 283 DF,  p-value: 3.821e-06

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.

In [10]:
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...
Call:
arimax(x = Y, order = c(2, 0, 0), xtransf = X, transfer = list(c(0, 7)))

Coefficients:
         ar1      ar2  intercept   T1-MA0  T1-MA1   T1-MA2   T1-MA3   T1-MA4
      1.5395  -0.6332    53.3713  -0.0758  0.0646  -0.0474  -0.5462  -0.6586
s.e.  0.0471   0.0512     0.1509   0.0761  0.0823   0.0822   0.0822   0.0822
       T1-MA5   T1-MA6   T1-MA7
      -0.8521  -0.4864  -0.3661
s.e.   0.0823   0.0824   0.0768

sigma^2 estimated as 0.05795:  log likelihood = -0.11,  aic = 22.22

Training set error measures:
                       ME      RMSE       MAE          MPE      MAPE      MASE
Training set 0.0007817914 0.2407192 0.1750768 -0.002986802 0.3260501 0.2921248
                   ACF1
Training set 0.03195547

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:

In [11]:
coeftest(output_XY)
z test of coefficients:

           Estimate Std. Error  z value  Pr(>|z|)    
ar1        1.539480   0.047117  32.6736 < 2.2e-16 ***
ar2       -0.633190   0.051187 -12.3701 < 2.2e-16 ***
intercept 53.371263   0.150906 353.6727 < 2.2e-16 ***
T1-MA0    -0.075778   0.076122  -0.9955    0.3195    
T1-MA1     0.064639   0.082299   0.7854    0.4322    
T1-MA2    -0.047378   0.082201  -0.5764    0.5644    
T1-MA3    -0.546225   0.082185  -6.6463 3.006e-11 ***
T1-MA4    -0.658598   0.082205  -8.0117 1.131e-15 ***
T1-MA5    -0.852063   0.082274 -10.3564 < 2.2e-16 ***
T1-MA6    -0.486423   0.082422  -5.9016 3.600e-09 ***
T1-MA7    -0.366137   0.076802  -4.7673 1.867e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

In [12]:
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)
Call:
arimax(x = Y, order = c(2, 0, 0), fixed = c(NA, NA, NA, 0, 0, 0, NA, NA, NA, 
    NA, NA), xtransf = X, transfer = list(c(0, 7)))

Coefficients:
         ar1      ar2  intercept  T1-MA0  T1-MA1  T1-MA2   T1-MA3   T1-MA4
      1.5379  -0.6291    53.3764       0       0       0  -0.5558  -0.6445
s.e.  0.0470   0.0509     0.1552       0       0       0   0.0778   0.0812
       T1-MA5   T1-MA6   T1-MA7
      -0.8598  -0.4836  -0.3633
s.e.   0.0810   0.0810   0.0773

sigma^2 estimated as 0.05823:  log likelihood = -0.82,  aic = 17.64

Training set error measures:
                       ME      RMSE       MAE          MPE      MAPE      MASE
Training set 0.0008444532 0.2413084 0.1749322 -0.003005804 0.3258383 0.2918835
                   ACF1
Training set 0.03481904

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.

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

In [14]:
coeftest(output_XY)
z test of coefficients:

           Estimate Std. Error  z value  Pr(>|z|)    
ar1        1.537908   0.047025  32.7041 < 2.2e-16 ***
ar2       -0.629109   0.050904 -12.3588 < 2.2e-16 ***
intercept 53.376447   0.155169 343.9893 < 2.2e-16 ***
T1-MA3    -0.555777   0.077798  -7.1439 9.074e-13 ***
T1-MA4    -0.644471   0.081157  -7.9410 2.005e-15 ***
T1-MA5    -0.859807   0.081014 -10.6131 < 2.2e-16 ***
T1-MA6    -0.483575   0.081028  -5.9680 2.402e-09 ***
T1-MA7    -0.363310   0.077342  -4.6974 2.634e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Tous les coefficients sont significatifs: le modèle est valide.


4. Modèles à fonction de transfert quotient de deux polynomes en B

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

In [15]:
# 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
Call:
arimax(x = Y, order = c(2, 0, 0), fixed = c(NA, NA, NA, NA, 0, 0, 0, NA, NA), 
    xtransf = X, transfer = list(c(1, 4)))

Coefficients:
         ar1      ar2  intercept  T1-AR1  T1-MA0  T1-MA1  T1-MA2   T1-MA3
      1.5177  -0.6267    53.3509   0.674       0       0       0  -0.4551
s.e.  0.0474   0.0498     0.1323   0.022       0       0       0   0.0740
       T1-MA4
      -0.6320
s.e.   0.0936

sigma^2 estimated as 0.06113:  log likelihood = -7.82,  aic = 27.64

Training set error measures:
                        ME      RMSE       MAE          MPE     MAPE     MASE
Training set -0.0001853234 0.2472487 0.1726197 -0.001528223 0.321501 0.288025
                   ACF1
Training set 0.02605155
z test of coefficients:

           Estimate Std. Error z value  Pr(>|z|)    
ar1        1.517672   0.047424  32.002 < 2.2e-16 ***
ar2       -0.626729   0.049769 -12.593 < 2.2e-16 ***
intercept 53.350931   0.132294 403.274 < 2.2e-16 ***
T1-AR1     0.673994   0.022023  30.604 < 2.2e-16 ***
T1-MA3    -0.455060   0.074042  -6.146 7.946e-10 ***
T1-MA4    -0.631985   0.093600  -6.752 1.459e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
In [16]:
# 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
Call:
arimax(x = Y, order = c(2, 0, 0), fixed = c(NA, NA, NA, NA, 0, 0, 0, NA, NA, 
    NA), xtransf = X, transfer = list(c(1, 5)))

Coefficients:
         ar1      ar2  intercept  T1-AR1  T1-MA0  T1-MA1  T1-MA2   T1-MA3
      1.5272  -0.6288    53.3618  0.5490       0       0       0  -0.5310
s.e.  0.0467   0.0495     0.1375  0.0392       0       0       0   0.0738
       T1-MA4   T1-MA5
      -0.3801  -0.5180
s.e.   0.1017   0.1086

sigma^2 estimated as 0.0571:  log likelihood = 2.08,  aic = 9.83

Training set error measures:
                       ME      RMSE       MAE          MPE      MAPE      MASE
Training set 0.0001700879 0.2389594 0.1681788 -0.001732428 0.3130213 0.2806151
                   ACF1
Training set 0.02877323
z test of coefficients:

           Estimate Std. Error  z value  Pr(>|z|)    
ar1        1.527181   0.046723  32.6859 < 2.2e-16 ***
ar2       -0.628841   0.049471 -12.7114 < 2.2e-16 ***
intercept 53.361773   0.137503 388.0769 < 2.2e-16 ***
T1-AR1     0.549027   0.039191  14.0089 < 2.2e-16 ***
T1-MA3    -0.530964   0.073814  -7.1933 6.325e-13 ***
T1-MA4    -0.380125   0.101704  -3.7376 0.0001858 ***
T1-MA5    -0.518006   0.108562  -4.7715 1.829e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
In [17]:
# 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)
Call:
arimax(x = Y, order = c(2, 0, 0), fixed = c(NA, NA, NA, NA, NA, 0, 0, 0, NA, 
    NA), xtransf = X, transfer = list(c(2, 4)))

Coefficients:
         ar1      ar2  intercept  T1-AR1   T1-AR2  T1-MA0  T1-MA1  T1-MA2
      1.5215  -0.6232    53.3667  1.0307  -0.2928       0       0       0
s.e.  0.0471   0.0499     0.1391  0.0892   0.0683       0       0       0
       T1-MA3   T1-MA4
      -0.4827  -0.3409
s.e.   0.0701   0.1293

sigma^2 estimated as 0.05849:  log likelihood = -1.41,  aic = 16.81

Training set error measures:
                        ME      RMSE       MAE          MPE     MAPE      MASE
Training set -0.0001953837 0.2418558 0.1699317 -0.002523852 0.316484 0.2835399
                  ACF1
Training set 0.0284708
z test of coefficients:

           Estimate Std. Error  z value  Pr(>|z|)    
ar1        1.521549   0.047082  32.3171 < 2.2e-16 ***
ar2       -0.623208   0.049922 -12.4836 < 2.2e-16 ***
intercept 53.366685   0.139140 383.5458 < 2.2e-16 ***
T1-AR1     1.030746   0.089230  11.5516 < 2.2e-16 ***
T1-AR2    -0.292786   0.068266  -4.2889 1.795e-05 ***
T1-MA3    -0.482654   0.070082  -6.8870 5.697e-12 ***
T1-MA4    -0.340854   0.129314  -2.6359  0.008392 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ces trois modèles sont valides.


5. Conclusion

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.

In [ ]: