4 - Qualité de la décomposition sous R

Formation - Désaisonnalisation avec JDemetra+ et RJDemetra

Auteur

Alain Quartier-la-Tente

L’objectif de ce TP est d’apprendre à étudier la qualité de la décomposition depuis RJDemetra

Si besoin, ci-dessous un exemple de code pour récupérer vos données :

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)

Prenons une spécification par défaut :

library(RJDemetra)
ipi_fr <- ipi_c_eu[, "FR"]
mysa <- x13(ipi_fr)

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

À partir d’un objet "X13", les statistiques m disponibles dans la sous-liste .$decomposition :

mysa$decomposition
Monitoring and Quality Assessment Statistics:
      M stats
M(1)    0.163
M(2)    0.089
M(3)    1.181
M(4)    0.558
M(5)    1.020
M(6)    0.090
M(7)    0.083
M(8)    0.244
M(9)    0.062
M(10)   0.272
M(11)   0.256
Q       0.368
Q-M2    0.402

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 ?

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

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

Mais en faisant cela on crée de la saisonnalité résiduelle ! (mysa2$diagnostics), 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             FALSE
new_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             FALSE

Sur 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$diagnostics
Relative 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         2.251
 Seasonal     59.750
 Irregular     1.067
 TD & Hol.     2.610
 Others       33.718
 Total        99.395

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.034
 
 Identifiable seasonality present

Residual seasonality tests
                                      P.value
 qs test on sa                          0.985
 qs test on i                           0.865
 f-test on sa (seasonal dummies)        0.958
 f-test on i (seasonal dummies)         0.893
 Residual seasonality (entire series)   0.876
 Residual seasonality (last 3 years)    0.906
 f-test on sa (td)                      0.987
 f-test on i (td)                       0.993

Ces 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 rjd3sa il est en revanche possible de calculer la tous les tests à l’exception du f-test.

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)`
NULL

Pour 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