fichier <- "../data/data_rte.xlsx"
# # Ou en téléchargeant le fichier depuis internet :
# fichier <- tempfile(fileext = "xlsx")
# url <- "https://aqlt.github.io/formation.2021.rte.cvs/data/data_rte.xlsx"
# download.file(url, fichier)
data_rte <- readxl::read_excel(fichier)
date_deb <- 2006
data_rte <- ts(data_rte[,-1], start = date_deb,
frequency = 12)3 - Qualité du préajustement sous R
Formation - Désaisonnalisation avec JDemetra+ et RJDemetra 
L’objectif de ce TP est d’apprendre à vérifier la qualité du pré-ajustement dans RJDemetra
Si besoin, ci-dessous un exemple de code pour récupérer vos données :
Prenons une spécification par défaut :
library(RJDemetra)
ipi_fr <- ipi_c_eu[, "FR"]
mysa <- x13(ipi_fr)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$regarima)y = regression model + arima (2, 1, 1, 0, 1, 1)
Model: RegARIMA - X13
Estimation span: from 1-1990 to 12-2020
Log-transformation: no
Regression model: no mean, trading days effect(7), leap year effect, Easter effect, outliers(4)
Coefficients:
ARIMA:
Estimate Std. Error T-stat Pr(>|t|)
Phi(1) 0.0003269 0.1077296 0.003 0.9976
Phi(2) 0.1688192 0.0740996 2.278 0.0233 *
Theta(1) -0.5485606 0.1016550 -5.396 1.24e-07 ***
BTheta(1) -0.6660849 0.0422242 -15.775 < 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.014638 *
Tuesday 0.88221 0.22832 3.864 0.000132 ***
Wednesday 1.03996 0.22930 4.535 7.85e-06 ***
Thursday 0.04943 0.22944 0.215 0.829549
Friday 0.91132 0.22988 3.964 8.88e-05 ***
Saturday -1.57769 0.22775 -6.927 1.99e-11 ***
Leap year 2.15403 0.70527 3.054 0.002425 **
Easter [1] -2.37950 0.45391 -5.242 2.71e-07 ***
TC (4-2020) -35.59245 2.17330 -16.377 < 2e-16 ***
AO (3-2020) -20.89026 2.18013 -9.582 < 2e-16 ***
AO (5-2011) 13.49850 1.85694 7.269 2.28e-12 ***
LS (11-2008) -12.54901 1.63554 -7.673 1.60e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.218 on 342 degrees of freedom
Log likelihood = -799.1, aic = 1632, aicc = 1634, bic(corrected for length) = 1.855library(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
Hypothesis:
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
Hypothesis:
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 ' ' 1Concernant 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 .$regarima$residuals.stat$tests :
mysa$regarima$residuals.stat$tests
Normality
Statistic P.value
mean 0.12648 0.8994 ***
skewness -0.01954 0.8799 ***
kurtosis 3.54844 0.0339
Signif. codes: H0 (normality of residuals) is not rejected at
significance levels: 0.1 ***0.05 **
Independence
Statistic P.value
ljung box 55.08622 0.0000
ljung box (residuals at seasonal lags) 3.09960 0.2123 ***
Signif. codes: H0 (independence of residuals) is not rejected at
significance levels: 0.1 ***0.05 **
Linearity
Statistic P.value
ljung box (squared residuals) 34.36237 0.0238
Signif. codes: H0 (no conditional heteroscedasticity of residuals) is not rejected at
significance levels: 0.1 ***0.05 **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, rajouter dans la spécification outlier.usedefcv = FALSE et prendre une valeur de outlier.cv inférieur à 4, qui est la valeur par défaut). 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 <- x13(ipi_fr, x13_spec(mysa, outlier.usedefcv = FALSE,
outlier.cv = 3))
# Bien plus d'outliers sont détectés !
summary(mysa2$regarima)y = regression model + arima (3, 1, 1, 0, 1, 1)
Model: RegARIMA - X13
Estimation span: from 1-1990 to 12-2020
Log-transformation: no
Regression model: no mean, trading days effect(7), leap year effect, Easter effect, outliers(28)
Coefficients:
ARIMA:
Estimate Std. Error T-stat Pr(>|t|)
Phi(1) 0.27920 0.09413 2.966 0.00322 **
Phi(2) 0.36429 0.07687 4.739 3.10e-06 ***
Phi(3) 0.11811 0.07680 1.538 0.12496
Theta(1) -0.63217 0.08106 -7.799 6.84e-14 ***
BTheta(1) -0.34211 0.05369 -6.372 5.72e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Regression model:
Estimate Std. Error T-stat Pr(>|t|)
Monday 0.3870 0.1811 2.137 0.033265 *
Tuesday 0.9721 0.1818 5.348 1.58e-07 ***
Wednesday 0.8970 0.1830 4.903 1.43e-06 ***
Thursday 0.1836 0.1826 1.006 0.315205
Friday 0.9678 0.1846 5.244 2.69e-07 ***
Saturday -1.6723 0.1829 -9.145 < 2e-16 ***
Leap year 1.7207 0.5026 3.423 0.000690 ***
Easter [1] -2.3262 0.3536 -6.578 1.68e-10 ***
TC (4-2020) -34.9803 1.3450 -26.007 < 2e-16 ***
AO (3-2020) -20.5598 1.6847 -12.204 < 2e-16 ***
AO (5-2011) 11.6856 1.3242 8.825 < 2e-16 ***
LS (11-2008) -9.4587 1.0474 -9.030 < 2e-16 ***
LS (1-2009) -7.8963 1.0256 -7.699 1.34e-13 ***
AO (6-2019) -5.7932 1.4256 -4.064 5.94e-05 ***
LS (8-2009) 5.7013 0.8196 6.957 1.66e-11 ***
TC (1-2011) 5.6019 1.0501 5.334 1.70e-07 ***
AO (5-2018) -4.5600 1.3401 -3.403 0.000743 ***
LS (5-2008) -4.7919 0.7824 -6.125 2.38e-09 ***
AO (5-2000) 5.2033 1.3250 3.927 0.000103 ***
AO (6-2003) -5.4817 1.3322 -4.115 4.81e-05 ***
AO (5-1991) -4.5123 1.3544 -3.332 0.000953 ***
LS (5-1994) 3.2027 0.7772 4.121 4.69e-05 ***
TC (12-2009) -4.5311 1.0961 -4.134 4.44e-05 ***
LS (3-1997) 3.5716 0.7567 4.720 3.38e-06 ***
LS (1-1993) -3.4933 0.7611 -4.590 6.15e-06 ***
AO (8-2020) 6.3136 1.6987 3.717 0.000234 ***
TC (11-2000) 5.8864 1.0215 5.762 1.79e-08 ***
TC (8-2015) 3.5207 0.9715 3.624 0.000332 ***
TC (12-1999) 4.4543 1.0307 4.322 2.01e-05 ***
LS (10-1997) 2.7159 0.7569 3.588 0.000379 ***
TC (12-1994) 4.2812 0.9973 4.293 2.27e-05 ***
TC (10-2017) 3.9944 0.9801 4.076 5.65e-05 ***
LS (11-2019) -3.0752 0.9235 -3.330 0.000959 ***
LS (2-2004) 2.1923 0.7651 2.865 0.004409 **
AO (6-2011) -4.8397 1.3372 -3.619 0.000338 ***
TC (11-2011) 4.0079 1.0178 3.938 9.88e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.627 on 317 degrees of freedom
Log likelihood = -685.6, aic = 1455, aicc = 1467, bic(corrected for length) = 1.646La 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 dans RJDemetra : 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 .$regarima$model$spec_rslt) :
mysa <- x13(ipi_fr, x13_spec(mysa),
userdefined = c("diagnostics.fcast-insample-mean",
"diagnostics.fcast-outsample-mean",
"diagnostics.fcast-outsample-variance"))
mysa$regarima$model$spec_rslt Model T.span Log transformation Mean Trading days
1 RegARIMA - X13 from 1-1990 to 12-2020 FALSE FALSE 7
Leap year Easter Outliers
1 TRUE TRUE 4# Pour éviter outputs trop longs, l'affichage est réduit :
mysa$user_definedNames 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`
[1] 0.3057321 0.7599958
attr(,"description")
[1] "T with 340 degrees of freedom"
$`diagnostics.fcast-outsample-mean`
[1] -0.7781656 0.4370126
attr(,"description")
[1] "T with 340 degrees of freedom"
$`diagnostics.fcast-outsample-variance`
[1] 1.1383576 0.3128808
attr(,"description")
[1] "F with 18 degrees of freedom in the nominator and 341 degrees of freedom in the denominator"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 .$regarima$loglik) :
AIC(mysa)[1] 1632.169BIC(mysa)[1] 1698.185# Il y a un peu plus de critères que dans base R : AICc et BICc
mysa$regarima$loglik
logvalue -799.084484
np 17.000000
neffectiveobs 359.000000
aic 1632.168967
aicc 1633.963689
bic 1698.185448
bicc 1.855018
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 ?