4 - Qualité de la décomposition sous R

Désaisonnaliser une série temporelle

Auteur

Alain Quartier-la-Tente

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 :

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

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$decomposition
Monitoring 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 ?

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"))
mysa2$decomposition
Monitoring 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-Henderson
mysa2$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         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.277

Mais 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             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         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.333

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

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 ?

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

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.5948553
rjd3toolkit::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 denominator

Mais 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.5776524

Notes de bas de page

  1. Entre JDemetra+ 2.x et JDemetra+ 3.x la spécification du f-test a été changée↩︎