3 - Qualité du préajustement sous R

Désaisonnaliser une série temporelle

Auteur

Alain Quartier-la-Tente

L’objectif de ce TP est d’apprendre à vérifier la qualité du pré-ajustement.

Pour installer tous les packages utiles de ce TP, lancer le programme :

packages_to_install <- c("rjd3toolkit", "rjd3x13", "rjd3tramoseats", "rjd3providers", "rjd3workspace")

packages <- packages_to_install[! packages_to_install %in% installed.packages()[,"Package"]]
if (length(packages) > 0) {
        install.packages(
      packages, repos = c("https://aqlt.r-universe.dev", "https://cloud.r-project.org")
      )
}

Prenons une spécification par défaut :

library(rjd3toolkit)
library(rjd3x13)
ipi_fr <- RJDemetra::ipi_c_eu[, "FR"]
mysa <- rjd3x13::x13(ipi_fr,spec = "RSA5c")

Comme on l’a vu dans le TP2, les tests de Student peuvent être utilisés pour tester la significativité des coefficients, et on peut également faire des tests de Fisher avec le package car pour voir si l’on peut simplifier les régresseurs jours ouvrables. Voir également le TP2 pour les tests sur la présence de jours ouvrables résiduelle.

summary(mysa$result$preprocessing)

Log-transformation: no 
SARIMA model: (2,1,1) (0,1,1)

Coefficients
            Estimate Std. Error  T-stat Pr(>|t|)    
phi(1)     0.0003291  0.1050495   0.003    0.998    
phi(2)     0.1688151  0.0728072   2.319    0.021 *  
theta(1)  -0.5485603  0.0980879  -5.593 4.58e-08 ***
btheta(1) -0.6660914  0.0513396 -12.974  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Regression model:
                 Estimate Std. Error  T-stat Pr(>|t|)    
monday            0.55932    0.22801   2.453 0.014662 *  
tuesday           0.88221    0.22832   3.864 0.000133 ***
wednesday         1.03996    0.22930   4.535 7.97e-06 ***
thursday          0.04943    0.22944   0.215 0.829552    
friday            0.91132    0.22988   3.964 8.96e-05 ***
saturday         -1.57769    0.22775  -6.927 2.14e-11 ***
lp                2.15402    0.70527   3.054 0.002434 ** 
easter           -2.37950    0.45392  -5.242 2.78e-07 ***
LS (2008-11-01) -12.54899    1.63555  -7.673 1.77e-13 ***
AO (2011-05-01)  13.49848    1.85695   7.269 2.48e-12 ***
AO (2020-03-01) -20.89026    2.18013  -9.582  < 2e-16 ***
TC (2020-04-01) -35.59247    2.17330 -16.377  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Number of observations: 372, Number of effective observations: 359, Number of parameters: 17
Loglikelihood: -799.0845
Standard error of the regression (ML estimate): 2.217552 
AIC: 1632.169, AICc: 1633.964, BIC: 1698.185

library(car)
# On rejette l'hypothèse de nullité globale des coefficients
linearHypothesis(mysa,
                 c("monday","tuesday","wednesday","thursday","friday","saturday"),
                 c(0, 0, 0, 0, 0, 0), test = "F")


Linear hypothesis test:
monday = 0
tuesday = 0
wednesday = 0
thursday = 0
friday = 0
saturday = 0

Model 1: restricted model
Model 2: mysa

  Res.Df Df      F    Pr(>F)    
1    348                        
2    342  6 83.415 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# On pourrait rassembler les jours de la semaine :
linearHypothesis(mysa,
                 c("monday = tuesday","tuesday = wednesday","
                   wednesday = thursday", "thursday = friday"), test = "F")


Linear hypothesis test:
monday - tuesday = 0
tuesday - wednesday = 0
wednesday - thursday = 0
thursday - friday = 0

Model 1: restricted model
Model 2: mysa

  Res.Df Df      F  Pr(>F)  
1    346                    
2    342  4 2.1504 0.07429 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Concernant la qualité du modèle RegARIMA, on peut citer trois tests :

Le test d’indépendance des résidus
Le test d’homoscédasticité des résidus
Le test de normalité des résidus

Ces trois tests, également disponibles par des fonctions spécifiques sous R (la commande residuals(mysa) permet de récupérer les résidus du modèle), sont également disponibles dans le sous objet .$result$preprocessing$diagnostics :

mysa$result$preprocessing$diagnostics

$mean
Value: 0.1264845 
P-Value: 0.8994 

$skewness
Value: -0.01953818 
P-Value: 0.8799 

$kurtosis
Value: 3.548441 
P-Value: 0.0339 

$doornikhansen
Value: 5.453171 
P-Value: 0.0654 

$lb
Value: 55.08611 
P-Value: 0.0000 

$bp
Value: 52.42581 
P-Value: 0.0001 

$seaslb
Value: 3.099964 
P-Value: 0.2123 

$seasbp
Value: 2.979744 
P-Value: 0.2254 

$lb2
Value: 34.36227 
P-Value: 0.0238 

$bp2
Value: 33.25964 
P-Value: 0.0316 

$nruns
Value: 0.2895201 
P-Value: 0.7722 

$lruns
Value: 14.58226 
P-Value: 1.0000 

$nudruns
Value: -0.8784381 
P-Value: 0.3797 

$ludruns
Value: 6.191513 
P-Value: 1.0000

L’hétéroscédasticité et la non-normalité proviennent souvent de la présence de points atypiques non corrigés (pour jouer sur le seuil de détection, prendre une valeur de critical.value inférieur à 4, qui est la valeur par défaut, dans la fonction set_outliers). Changer le schéma de décomposition peut aussi aider (transform.function = "None" pour un modèle additif ou transform.function = "Log" pour un modèle multiplicatif) :

mysa2 <- rjd3x13::x13(
  ipi_fr, 
  rjd3x13::spec_x13("RSA5c") |> 
    set_outlier(critical.value = 3)
  )
# Bien plus d'outliers sont détectés !
summary(mysa2$result$preprocessing)

Log-transformation: no 
SARIMA model: (3,1,1) (0,1,1)

Coefficients
          Estimate Std. Error T-stat Pr(>|t|)    
phi(1)     0.27046    0.10353  2.612  0.00941 ** 
phi(2)     0.36526    0.08059  4.532 8.26e-06 ***
phi(3)     0.11235    0.07978  1.408  0.16003    
theta(1)  -0.58596    0.09184 -6.380 6.21e-10 ***
btheta(1) -0.40749    0.07589 -5.370 1.52e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Regression model:
                Estimate Std. Error  T-stat Pr(>|t|)    
monday            0.3888     0.1853   2.098 0.036688 *  
tuesday           0.9865     0.1858   5.309 2.08e-07 ***
wednesday         0.9400     0.1872   5.021 8.59e-07 ***
thursday          0.1620     0.1868   0.867 0.386634    
friday            0.9605     0.1891   5.080 6.45e-07 ***
saturday         -1.6873     0.1865  -9.045  < 2e-16 ***
lp                1.9863     0.5177   3.837 0.000150 ***
easter           -2.3026     0.3610  -6.378 6.31e-10 ***
AO (1991-05-01)  -4.4831     1.4122  -3.174 0.001648 ** 
LS (1993-01-01)  -3.3976     0.8510  -3.993 8.12e-05 ***
LS (1994-05-01)   3.1013     0.8612   3.601 0.000367 ***
TC (1994-12-01)   4.4390     1.0772   4.121 4.82e-05 ***
LS (1997-03-01)   3.4330     0.8465   4.056 6.29e-05 ***
LS (1997-10-01)   2.6811     0.8469   3.166 0.001696 ** 
TC (1999-12-01)   4.2817     1.1080   3.864 0.000135 ***
AO (2000-05-01)   5.4149     1.3765   3.934 0.000103 ***
TC (2000-11-01)   5.6357     1.0978   5.133 4.96e-07 ***
AO (2003-06-01)  -5.4691     1.3722  -3.986 8.35e-05 ***
LS (2008-05-01)  -4.6136     0.8664  -5.325 1.91e-07 ***
LS (2008-11-01)  -9.6863     1.1118  -8.712  < 2e-16 ***
LS (2009-01-01)  -7.5793     1.0958  -6.917 2.55e-11 ***
LS (2009-08-01)   5.5337     0.8948   6.184 1.91e-09 ***
TC (2009-12-01)  -4.7614     1.1491  -4.144 4.39e-05 ***
TC (2011-01-01)   4.7682     1.0969   4.347 1.86e-05 ***
AO (2011-05-01)  11.6597     1.3738   8.487 8.17e-16 ***
AO (2011-06-01)  -5.2043     1.3823  -3.765 0.000198 ***
TC (2012-10-01)  -3.0496     1.0626  -2.870 0.004381 ** 
TC (2015-08-01)   3.5699     1.0582   3.374 0.000833 ***
TC (2017-10-01)   4.0720     1.0691   3.809 0.000167 ***
AO (2018-05-01)  -4.3495     1.3934  -3.121 0.001965 ** 
AO (2019-06-01)  -5.7232     1.4681  -3.898 0.000118 ***
LS (2019-11-01)  -3.0612     0.9948  -3.077 0.002271 ** 
AO (2020-03-01) -20.6326     1.6975 -12.155  < 2e-16 ***
TC (2020-04-01) -34.9562     1.3959 -25.042  < 2e-16 ***
AO (2020-08-01)   6.0365     1.7042   3.542 0.000456 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Number of observations: 372, Number of effective observations: 359, Number of parameters: 41
Loglikelihood: -692.0561
Standard error of the regression (ML estimate): 1.655325 
AIC: 1466.112, AICc: 1476.977, BIC: 1625.328

La qualité des prévisions peut également être vérifiée à travers plusieurs tests :

Est-ce que la moyenne des erreurs prévisions in sample (i.e. : modèle estimé sur toute la période) et la moyenne des prévisions out of sample (i.e. : modèle estimé de manière dynamique en ajoutant une a à une les nouvelles données) sont nulles ? Ces tests sont sensibles à la non-normalité des résidus
Est-ce que les variances des erreurs de prévision in sample et out of sample sont les mêmes ? Ce test est sensible à la non-normalité des résidus
Est-ce qu’il y a “trop” d’outliers ? Dans JDemetra+, on considère par défaut qu’il y a trop d’outliers si la proportion d’outliers par rapport aux nombres d’observations est supérieure à 5 %.

Les trois premiers tests ne sont pas par défaut exportés : il faut les rajouter à la main avec le paramètre userdefined. Ils seront alors disponibles dans la sous-liste .$user_defined. Concernant la proportion d’outliers, elle peut être calculée à la main à partir du nombre d’outliers (par exemple disponible dans .$result_spec$regarima$regression$outliers) :

mysa <- x13(
  ipi_fr, 
  mysa$estimation_spec,
  userdefined = c("diagnostics.fcast-insample-mean",
                  "diagnostics.fcast-outsample-mean",
                  "diagnostics.fcast-outsample-variance")
)

# Pour éviter outputs trop longs, l'affichage est réduit :
mysa$user_defined

Names of additional variables (3):
diagnostics.fcast-insample-mean, diagnostics.fcast-outsample-mean, diagnostics.fcast-outsample-variance

# Pour supprimer cela, vous pouvez par exemple utiliser le code suivant :
c(mysa$user_defined)

$`diagnostics.fcast-insample-mean`
Value: 0.3057319 
P-Value: 0.7600 

$`diagnostics.fcast-outsample-mean`
Value: -0.7781645 
P-Value: 0.4370 

$`diagnostics.fcast-outsample-variance`
Value: 1.138358 
P-Value: 0.3129

Vous pouvez bien sûr utiliser votre tests préféré à partir de ceux disponibles sous R (autre test de normalité…).

Pour comparer différents modèles, vous pouvez également utiliser les critères d’information (mais il faut que les modèles ARIMA aient les mêmes ordres de différenciation !). Vous pouvez pour cela utiliser les fonctions de bases de R (AIC(), BIC()…) ou prendre ceux de JDemetra+ (affichés lors du summary(), qu’on peut également retrouver par la commande .$result$preprocessing$estimation$likelihood) :

AIC(mysa)

[1] 1632.169

BIC(mysa)

[1] 1698.185

# Il y a un peu plus de critères que dans base R :
mysa$result$preprocessing$estimation$likelihood

Number of observations: 372 
Number of effective observations: 359 
Number of parameters: 17 

Loglikelihood: -799.0845 
Standard error of the regression (ML estimate): 2.217552 
AIC: 1632.169 
AICC: 1633.964 
BIC: 1698.185

Exercice

Prenez une série et étudier la qualité du modèle RegARIMA. Essayer de changer quelques paramètres : est-ce que le nouveau modèle vous parait meilleur ou moins bien que l’ancien ?