packages_to_install <- c("RJDemetra", "remotes")
packages <- packages_to_install[! packages_to_install %in% installed.packages()[,"Package"]]
if (length(packages) > 0) {
install.packages(packages)
}
packages_to_install_git <- c("rjd3toolkit", "rjd3x13", "rjd3tramoseats", "rjd3providers", "rjdemetra3")
packages_git <- packages_to_install_git[! packages_to_install_git %in% installed.packages()[,"Package"]]
if (length(packages_git) > 0) {
# # Configurer si besoin le proxy
# proxy <- "proxy_a_definir"
# Sys.setenv(HTTPS_PROXY = proxy)
remotes::install_github(
sprintf("rjdemetra/%s", packages_git),
# option utile dans certaines installations portables de Java :
INSTALL_opts = "--no-multiarch")
}4 - Qualité de la décomposition sous R
Désaisonnaliser une série temporelle
L’objectif de ce TP est d’apprendre à étudier la qualité de la décomposition depuis RJDemetra.
Pour installer tous les packages utiles de ce TP, lancer le programme :
Dans ce TP nous allons voir différentes façon de vérifier la qualité de la décomposition. Tout d’abord, on peut commencer par regarder les statistiques m dont les définitions sont rappelées ci-dessous
| Poids | Description | Si problème | |
|---|---|---|---|
| M1 | 10 | Contribution de l’irrégulier à la variance totale (stationnarisation par différence d’ordre 3). Si trop élevé, difficile d’extraire la saisonnalité. | Points atypiques ou taille des filtres |
| M2 | 11 | Contribution de l’irrégulier à la variance totale (stationnarisation par une droite). | Points atypiques ou taille des filtres |
| M3 | 10 | Mesuré à partir du ratio I/C. Si trop grand on aura du mal à séparer les deux composantes. | Points atypiques ou taille des filtres |
| M4 | 8 | Test autocorrélation sur l’irrégulier (réduire le filtre saisonnier). | Filtre saisonnier plus court |
| M5 | 11 | Mesuré à partir du MCD (nombre de mois nécessaires pour que les variations absolues de la TC l’emporte sur I). | Points atypiques |
| M6 | 10 | Vérifie si la moyenne mobile M3x5 est appropriée (\(1.5 < I/S < 6.5\)). | Prendre filtre plus long |
| M7 | 18 | Permet de voir si la saisonnalité est identifiable (compare part relative de la saisonnalité stable et mobile). | Schéma multiplicatif ? |
| M8 | 7 | Mesure l’évolution de la S de court terme. | Changer filtre saisonnier |
| M9 | 7 | Mesure l’évolution de la S de long terme. | Changer filtre saisonnier |
| M10 | 4 | M8 sur dernières années (\(N-2\) à \(N-5\)). | Changer filtre saisonnier |
| M11 | 4 | M9 sur dernières années (\(N-2\) à \(N-5\)). | Changer filtre saisonnier |
1 RJDemetra
Prenons une spécification par défaut :
library(RJDemetra)
ipi_fr <- ipi_c_eu[, "FR"]
mysa <- x13(ipi_fr, "RSA4")À partir d’un objet "X13", les statistiques m disponibles dans la sous-liste .$decomposition :
mysa$decompositionMonitoring and Quality Assessment Statistics:
M stats
M(1) 0.127
M(2) 0.079
M(3) 1.094
M(4) 0.558
M(5) 1.093
M(6) 0.022
M(7) 0.085
M(8) 0.242
M(9) 0.064
M(10) 0.261
M(11) 0.247
Q 0.355
Q-M2 0.389
Final filters:
Seasonal filter: 3x5
Trend filter: 13 terms Henderson moving average
Exercice
Que signifie ces valeurs des statistiques m plus grandes que 1 ? Est-ce important ? Si oui comment les corriger ?
Indice
Que pensez-vous de la tendance (plot(mysa$final, type_chart = "sa-trend")) ? Quelle est la contribution du cycle à la variance totale (mysa$diagnostics$variance_decomposition) ?
Solution
La tendance est plutôt plate et la contribution du cycle à la variance totale est petite, on peut donc ignorer la statistique m3. La statistique M5 suggère de prendre un filtre saisonnier plus long, par exemple en utilisant le code suivant :
mysa2 <- x13(get_ts(mysa), x13_spec(mysa, x11.seasonalma = "S3X9"))
mysa2$decompositionMonitoring and Quality Assessment Statistics:
M stats
M(1) 0.287
M(2) 0.151
M(3) 1.490
M(4) 1.228
M(5) 1.385
M(6) 0.322
M(7) 0.087
M(8) 0.163
M(9) 0.062
M(10) 0.133
M(11) 0.127
Q 0.541
Q-M2 0.595
Final filters:
Seasonal filter: 3x9
Trend filter: 23-Hendersonmysa2$diagnosticsRelative contribution of the components to the stationary
portion of the variance in the original series,
after the removal of the long term trend
Trend computed by Hodrick-Prescott filter (cycle length = 8.0 years)
Component
Cycle 1.728
Seasonal 51.002
Irregular 1.253
TD & Hol. 2.178
Others 44.904
Total 101.065
Combined test in the entire series
Non parametric tests for stable seasonality
P.value
Kruskall-Wallis test 0.000
Test for the presence of seasonality assuming stability 0.000
Evolutive seasonality test 0.035
Identifiable seasonality present
Residual seasonality tests
P.value
qs test on sa 0.066
qs test on i 0.286
f-test on sa (seasonal dummies) 0.889
f-test on i (seasonal dummies) 0.845
Residual seasonality (entire series) 0.757
Residual seasonality (last 3 years) 0.948
f-test on sa (td) 0.102
f-test on i (td) 0.277Mais en faisant cela on crée de la saisonnalité résiduelle ! C’est donc mieux de rester sur les paramètres par défaut.
Alors que pour changer le filtre saisonnier il suffit d’utiliser le paramètre x11.seasonalma, pour changer la longueur du filtre de Henderson il faut désactiver l’option de recherche automatique de la longueur du filtre (x11.trendAuto = FALSE) et spécifier la longueur dans le paramètre x11.trendma:
new_spec <- x13_spec(mysa, x11.trendma = 15)
new_spec$x11# Colonne trendma inchangée ! x11.mode x11.seasonalComp x11.lsigma x11.usigma x11.trendAuto
Predefined Undefined TRUE 1.5 2.5 TRUE
User_modif <NA> NA NA NA NA
Final Undefined TRUE 1.5 2.5 TRUE
x11.trendma x11.seasonalma x11.fcasts x11.bcasts x11.calendarSigma
Predefined 13 Msr -1 0 None
User_modif 15 <NA> NA NA <NA>
Final 13 Msr -1 0 None
x11.sigmaVector x11.excludeFcasts
Predefined NA FALSE
User_modif NA NA
Final NA FALSEnew_spec <- x13_spec(mysa, x11.trendma = 15, x11.trendAuto = FALSE)
new_spec$x11 x11.mode x11.seasonalComp x11.lsigma x11.usigma x11.trendAuto
Predefined Undefined TRUE 1.5 2.5 TRUE
User_modif <NA> NA NA NA FALSE
Final Undefined TRUE 1.5 2.5 FALSE
x11.trendma x11.seasonalma x11.fcasts x11.bcasts x11.calendarSigma
Predefined 13 Msr -1 0 None
User_modif 15 <NA> NA NA <NA>
Final 15 Msr -1 0 None
x11.sigmaVector x11.excludeFcasts
Predefined NA FALSE
User_modif NA NA
Final NA FALSESur la qualité de la décomposition, la sous liste .$diagnostics contient les contributions des différentes composantes à la variance de la série, le test combiné et les tests sur la saisonnalité et jours ouvrables résiduels :
mysa$diagnosticsRelative contribution of the components to the stationary
portion of the variance in the original series,
after the removal of the long term trend
Trend computed by Hodrick-Prescott filter (cycle length = 8.0 years)
Component
Cycle 1.830
Seasonal 51.089
Irregular 0.927
TD & Hol. 2.179
Others 44.916
Total 100.941
Combined test in the entire series
Non parametric tests for stable seasonality
P.value
Kruskall-Wallis test 0.000
Test for the presence of seasonality assuming stability 0.000
Evolutive seasonality test 0.014
Identifiable seasonality present
Residual seasonality tests
P.value
qs test on sa 0.924
qs test on i 0.643
f-test on sa (seasonal dummies) 0.671
f-test on i (seasonal dummies) 0.453
Residual seasonality (entire series) 0.415
Residual seasonality (last 3 years) 0.954
f-test on sa (td) 0.091
f-test on i (td) 0.333Ces tests sont effectués sur l’ensemble de la série, alors que dans le main result le f-test est effectué sur les 8 dernières années. Il n’est pour l’instant pas possible d’exporter les tests de saisonnalité résiduelle sur les 8 ou 10 dernières années. À partir du packages rjd3toolkit il est en revanche possible de calculer tous les tests à l’exception du f-test1.
Par rapport aux éléments vus en cours, les msr par mois sont exportables en utilisant le paramètre userdefined de x13. Il y a cependant actuellement un bug qui ne permet pas de l’exporter pour le dernier mois :
mysa <- x13(ipi_fr,
userdefined = c("diagnostics.msr-global",
sprintf("diagnostics.msr(%i)", 1:12)))
c(mysa$user_defined)$`diagnostics.msr-global`
[1] 4.224309
$`diagnostics.msr(1)`
[1] 6.918248
$`diagnostics.msr(2)`
[1] 4.679206
$`diagnostics.msr(3)`
[1] 4.351482
$`diagnostics.msr(4)`
[1] 5.808313
$`diagnostics.msr(5)`
[1] 4.446801
$`diagnostics.msr(6)`
[1] 3.175827
$`diagnostics.msr(7)`
[1] 5.102764
$`diagnostics.msr(8)`
[1] 2.756466
$`diagnostics.msr(9)`
[1] 4.216562
$`diagnostics.msr(10)`
[1] 4.750469
$`diagnostics.msr(11)`
[1] 5.088831
$`diagnostics.msr(12)`
NULLPour extraire tous les MSR, préférer la solution suivante :
extract_msr <- function(x, i = 1:12){
jmodel <- suppressWarnings(jx13(get_ts(x), x13_spec(x)))
jres <- jmodel$result@internal$getResults()
jres <- new(Class = "X13_java", internal = jres)
res <- sapply(i, function(i_){
RJDemetra:::result(jres,
sprintf("msr(%i)", i_))
})
names(res) <- sprintf("msr(%i)", i)
res
}
extract_msr(mysa) msr(1) msr(2) msr(3) msr(4) msr(5) msr(6) msr(7) msr(8)
6.918248 4.679206 4.351482 5.808313 4.446801 3.175827 5.102764 2.756466
msr(9) msr(10) msr(11) msr(12)
4.216562 4.750469 5.088831 2.953120 2 rjdemetra3
Prenons une spécification par défaut :
Les statistiques m disponibles dans la sous-liste .$result$mstats :
sa_jd3$result$mstats$m1
[1] 0.1267867
$m2
[1] 0.07923419
$m3
[1] 1.094359
$m4
[1] 0.5580624
$m5
[1] 1.093318
$m6
[1] 0.02176204
$m7
[1] 0.08520675
$m8
[1] 0.2423733
$m9
[1] 0.06360174
$m10
[1] 0.2606264
$m11
[1] 0.247436
$q
[1] 0.3549944
$qm2
[1] 0.3890771
Exercice
Que signifie ces valeurs des statistiques m plus grandes que 1 ? Est-ce important ? Si oui comment les corriger ?
Indice
Que pensez-vous de la tendance (plot(sa_jd3)) ? Quelle est la contribution du cycle à la variance totale (sa_jd3$result$diagnostics$vardecomposition ou rjd3toolkit::diagnostics(sa_jd3)[[1]]) ?
Solution
La tendance est plutôt plate et la contribution du cycle à la variance totale est petite, on peut donc ignorer la statistique m3. La statistique M5 suggère de prendre un filtre saisonnier plus long, par exemple en utilisant le code suivant :
sa2_jd3 <- rjd3x13::x13(
sa_jd3$result$preadjust$a1,
sa_jd3$estimation_spec |>
rjd3x13::set_x11(seasonal.filter = "S3X9"))
sa2_jd3$result$mstats$m1
[1] 0.2865214
$m2
[1] 0.1511407
$m3
[1] 1.490081
$m4
[1] 1.227737
$m5
[1] 1.38513
$m6
[1] 0.3221387
$m7
[1] 0.08654842
$m8
[1] 0.1630884
$m9
[1] 0.06162402
$m10
[1] 0.1329491
$m11
[1] 0.1266407
$q
[1] 0.5406235
$qm2
[1] 0.5948553rjd3toolkit::diagnostics(sa2_jd3)[["residual_tests"]] Statistic P.value
seas.ftest.i 0.8077490 0.63220978
seas.ftest.sa 0.8422258 0.59838207
seas.qstest.i 2.5023511 0.28616820
seas.qstest.sa 5.4440142 0.06574267
td.ftest.i 1.1983675 0.30637684
td.ftest.sa 1.0445406 0.39599185
Description
seas.ftest.i F with 11.0 degrees of freedom in the nominator and 131.0 degrees of freedom in the denominator
seas.ftest.sa F with 11.0 degrees of freedom in the nominator and 131.0 degrees of freedom in the denominator
seas.qstest.i Chi2 with 2.0 degrees of freedom
seas.qstest.sa Chi2 with 2.0 degrees of freedom
td.ftest.i F with 6.0 degrees of freedom in the nominator and 365.0 degrees of freedom in the denominator
td.ftest.sa F with 6.0 degrees of freedom in the nominator and 365.0 degrees of freedom in the denominatorMais en faisant cela on crée de la saisonnalité résiduelle ! C’est donc mieux de rester sur les paramètres par défaut.
Par rapport à RJDemetra, plus d’éléments associés aux MSR sont exportables :
sa_jd3 <- rjd3x13::x13(
ipi_fr, "rsa4",
userdefined = c("decomposition.d9-global-msr",
"decomposition.d9-msr",
"decomposition.d9-msr-table")
)
c(sa_jd3$user_defined)$`decomposition.d9-global-msr`
[1] 4.054405
$`decomposition.d9-msr`
[1] 6.432093 4.778968 5.061714 4.681529 4.177106 3.200465 4.651325 2.617721
[9] 3.885137 4.256683 6.629040 2.577652
$`decomposition.d9-msr-table`
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0.9410370 0.8787867 1.0086908 0.7913817 1.0132386 0.8664120 0.8931384
[2,] 0.1463034 0.1838863 0.1992785 0.1690434 0.2425695 0.2707145 0.1920181
[3,] 6.4320925 4.7789676 5.0617144 4.6815286 4.1771062 3.2004645 4.6513246
[,8] [,9] [,10] [,11] [,12]
[1,] 0.9665417 0.8704061 1.0876705 0.8972216 0.9300347
[2,] 0.3692303 0.2240348 0.2555207 0.1353471 0.3608069
[3,] 2.6177208 3.8851372 4.2566834 6.6290402 2.5776524Notes de bas de page
Entre JDemetra+ 2.x et JDemetra+ 3.x la spécification du f-test a été changée↩︎