1 - R et JDemetra+

Désaisonnaliser une série temporelle

Auteur

Alain Quartier-la-Tente

L’objectif de ce TP est d’apprendre à manipuler JDemetra+ sous R à travers le package RJDemetra.

Pour manipuler JDemetra+ sous R il y a actuellement deux façons :

Utiliser le JWSACruncher qui permet, à partir de la console, de mettre à jour un workspace JDemetra+ et d’exporter les résultats sans devoir ouvrir le logiciel. Pour faciliter son utilisation depuis R, le package rjwsacruncher peut être utilisé (voir TP5 - JDemetra+ en production).
Utiliser le package RJDemetra qui permet d’effectuer des désaisonnalisations avec les mêmes algorithmes et paramètres que JDemetra+ et de manipuler des workspaces.

Dans ce TP on utilisera les données du package RJDemetra mais n’hésitez pas à utiliser vos propres séries !

Pour faire de la désaisonnalisation sous R il existe plusieurs packages :

seasonal et x12 qui permettent de faire du X-13ARIMA-SEATS en utilisant les programmes du US Census Bureau
RJDemetra qui est une interface R à JDemetra+ et c’est ce package que l’on va étudier.

RJDemetra est sur le CRAN et se base sur les librairies Java de JDemetra+. Pour l’utiliser il faut avoir Java 8 ou plus. En cas de problème d’installation voir la page : https://github.com/rjdverse/rjdemetra/wiki/Installation-manual.

Le package a aussi un site web (https://rjdverse.github.io/rjdemetra/).

RJDemetra permet :

De faire des modèles RegARIMA, TRAMO-SEATS and X-13-ARIMA comme dans JDemetra+ en définissant sa propre spécification
Manipuler les workspaces de JDemetra+
Importer les modèles workspaces sauvegardés par JDemetra+
Exporter les modèles créés sous RJDemetra

Une nouvelle version de RJDemetra est en cours de développement autour de la version 3.0 de JDemetra+. Elle nécessite toutefois d’avoir une version de Java supérieure à la version 17 et elle n’est pas sur le CRAN. Les fonctionnalités de RJDemetra (et les nouvelles) sont divisées en plusieurs packages disponibles sous https://github.com/rjdverse. Dans cette formation, nous utiliserons les packages suivants RJDemetra, rjd3toolkit, rjd3tramoseats, rjd3providers, rjd3workspace :

packages_to_install <- c("RJDemetra", "remotes", "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")
  )
}

1 RJDemetra v3

1.1 Créer une specification

Les modèles autour de la méthode X-13ARIMA peuvent être estimés avec rjd3x13::x13(), rjd3x13::regarima()) et rjd3x13::x11().

Les spécifications peuvent être crées avec les fonctions rjd3x13::spec_regarima, rjd3x13::spec_x11() ou rjd3x13::spec_x13() et peuvent être modifiées par :

Pour le préajustement : rjd3toolkit::set_arima(), rjd3toolkit::set_automodel(), rjd3toolkit::set_basic(), rjd3toolkit::set_easter(), rjd3toolkit::set_estimate(), rjd3toolkit::set_outlier(), rjd3toolkit::set_tradingdays(), rjd3toolkit::set_transform(), rjd3toolkit::add_outlier(), rjd3toolkit::remove_outlier(), rjd3toolkit::add_ramp(), rjd3toolkit::remove_ramp(), rjd3toolkit::add_usrdefvar() ;
Pour la décomposition : rjd3x13::set_x11() ;
Pour le benchmarking : rjd3toolkit::set_benchmarking().

Exercice

Faire la désaisonnalisation d’une série avec X-13 avec la spécification suivante :

détection automatique du schéma de décomposition, des outliers et du modèle ARIMA ;
une correction des jours ouvrables “working days” et un effet graduel de Pâques.

Faire ensuite un graphique avec la série brute et la série désaisonnalisée.

Indice

Utiliser la spécification RSA4c pour la désaisonnalisation.

Si le modèle créé s’appelle sa_jd3, regarder les valeurs de sa_jd3$result$final et rjd3toolkit::sa_decomposition(sa_jd3).

Solution

sa_jd3 <- rjd3x13::x13(ipi_c_eu[, "FR"], "rsa4")
sa_jd3

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

SARIMA coefficients:
   phi(1)    phi(2)  theta(1) btheta(1) 
  0.05291   0.18672  -0.52138  -0.66132 

Regression model:
             td              lp          easter LS (2008-11-01) LS (2009-01-01) 
         0.6927          2.0903         -2.5476         -9.2744         -7.2838 
AO (2011-05-01) AO (2020-03-01) TC (2020-04-01) 
        13.1870        -21.1492        -35.6481 

 Seasonal filter: S3X3;  Trend filter: H-13 terms
 M-Statistics: q Good (0.355); q-m2 Good (0.389)
 QS test on SA: Good (0.924);  F-test on SA: Good (0.632)

For a more detailed output, use the 'summary()' function.

y <- sa_jd3$result$preadjust$a1
sa <- sa_jd3$result$final$d11final
# ou :
decomp <- rjd3toolkit::sa_decomposition(sa_jd3)
y <- decomp$series
sa <- decomp$sa
# ou on peut directement utiliser les fonctions de rjd3x13 :
plot(sa_jd3, first_date = 2000 #Pour n'afficher le graphique qu'à partir de 200
)

Pour des graphiques ggplot2, on peut également utiliser le package ggdemetra3.

Exercice

Modifier le modèle précédent pour enlever l’effet graduel de Pâques.

Solution

spec_sans_easter_v3 <- 
    sa_jd3$estimation_spec |> 
    rjd3toolkit::set_easter(enabled = FALSE)
sa2_jd3 <- rjd3x13::x13(ipi_c_eu[, "FR"], spec_sans_easter_v3)
sa2_jd3$result$preprocessing

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

SARIMA coefficients:
   phi(1)    phi(2)  theta(1) btheta(1) 
  0.02527   0.10075  -0.59512  -0.71292 

Regression model:
             td              lp LS (2008-11-01) AO (2011-05-01) AO (2020-03-01) 
         0.7032          2.1641        -12.7560         13.1047        -20.3006 
TC (2020-04-01) 
       -36.4048 

For a more detailed output, use the 'summary()' function.

Exercice

Calculer les p-valeurs associées au modèle Reg-ARIMA de la précédente spécification.

Indice

Récupérer le modèle Reg-ARIMA et utiliser la fonction summary().

Solution

summary(sa2_jd3$result$preprocessing)

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

Coefficients
          Estimate Std. Error  T-stat Pr(>|t|)    
phi(1)     0.02527    0.09998   0.253    0.801    
phi(2)     0.10075    0.07425   1.357    0.176    
theta(1)  -0.59512    0.08863  -6.715 7.65e-11 ***
btheta(1) -0.71292    0.04618 -15.439  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Regression model:
                 Estimate Std. Error  T-stat Pr(>|t|)    
td                0.70321    0.03417  20.582  < 2e-16 ***
lp                2.16412    0.76568   2.826  0.00498 ** 
LS (2008-11-01) -12.75601    1.66740  -7.650 1.99e-13 ***
AO (2011-05-01)  13.10471    1.97806   6.625 1.32e-10 ***
AO (2020-03-01) -20.30056    2.25944  -8.985  < 2e-16 ***
TC (2020-04-01) -36.40477    2.19426 -16.591  < 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: 11
Loglikelihood: -817.6149
Standard error of the regression (ML estimate): 2.330069 
AIC: 1657.23, AICc: 1657.991, BIC: 1699.946

1.2 Créer un workspace

Dans cette partie nous allons créer un workspace depuis R. Pour cela les fonctions qui peuvent être utilisées sont rjd3workspace::jws_new(), rjd3workspace::read_workspace(), rjd3workspace::jws_sap_new(), rjd3workspace::add_sa_item(), rjd3workspace::save_workspace().

Exercice

Créer un workspace qui va contenir une série désaisonnalisée selon 3 spécifications différentes.

Solution

jws <- rjd3workspace::jws_new()
jsap <- rjd3workspace::jws_sap_new(jws, "MP-1")
rjd3workspace::add_sa_item(jsap, "X13 avec Pâques", sa_jd3)
rjd3workspace::add_sa_item(jsap, "X13 sans Pâques", sa2_jd3)
rjd3workspace::add_sa_item(jsap, "TRAMO-SEATS", y, rjd3tramoseats::spec_tramoseats())
rjd3workspace::save_workspace(jws, "ws_v3.xml")

Exercice

Importer le workspace précédent et récupérer :

Le nom du premier multi-processing
Le nombre de modèles dans ce premier multi-processing
L’ensemble des séries brutes
Le 2ème modèle

Solution

# charge tous les modèles du workspace :
jws <- rjd3workspace::.jws_open("ws_v3.xml")
ws <- rjd3workspace::read_workspace(jws)
names(ws$processing)[1]

[1] "MP-1"

length(ws$processing[[1]]) # Nombre de multiprocessing

[1] 3

Autre façon de faire similaire à RJDemetra :

# charge tous les modèles du workspace :
jws <- rjd3workspace::.jws_open("ws_v3.xml")
rjd3workspace::jws_compute(jws)
rjd3workspace::ws_sap_count(jws) # Nombre de multiprocessing

[1] 1

jsap1 <- rjd3workspace::jws_sap(jws, 1) # Le premier multiprocessing
rjd3workspace::sap_name(jsap1)

[1] "MP-1"

rjd3workspace::sap_sai_count(jsap1)

[1] 3

model2 <- rjd3workspace::jsap_sai(jsap1, 2) # On récupère l'objet associé au 2ème modèle
rjd3workspace::read_sai(model2)

$ts
$name
[1] "X13 sans Pâques"

$moniker
$source
[1] ""

$id
[1] "a5e9326a-a1fd-43d8-84ee-3456d5f3e754"

attr(,"class")
[1] "JD3_TSMONIKER"

$metadata
NULL

$data
       Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
1990  92.1  92.3 102.1  93.0  93.3 100.8  92.9  66.7  95.8 105.0  96.7  89.2
1991  92.5  89.2  97.4  93.8  87.5 100.3  93.4  64.3  96.9 103.5  94.0  92.1
1992  90.7  89.0  99.4  93.7  86.1 101.3  90.4  62.9  96.6  98.4  91.9  92.6
1993  82.3  84.0  95.6  88.3  82.2  97.9  85.5  61.3  93.7  93.0  88.3  92.1
1994  83.6  83.7  97.0  88.3  88.3 102.9  87.3  65.9  98.2  98.0  96.8  98.0
1995  91.8  90.1 102.9  90.4  91.6 103.7  90.6  66.8  98.7 101.4  97.2  94.8
1996  92.0  91.1  98.1  94.3  90.5 101.8  96.1  66.3  98.9 105.0  95.0  96.0
1997  91.9  91.3  99.1 102.8  93.2 108.2 100.4  70.5 107.3 114.1  99.6 106.7
1998  98.2  98.7 109.3 103.7  97.6 114.7 106.1  72.1 111.5 112.6 105.6 107.4
1999  97.2  98.3 114.5 104.8  99.9 120.2 105.7  76.1 115.2 115.1 111.1 114.0
2000 103.4 107.5 121.7 105.7 113.1 119.4 108.1  82.0 116.4 121.3 117.2 111.9
2001 110.7 108.9 124.0 109.3 109.8 121.9 112.4  85.5 114.1 123.4 114.2 104.9
2002 108.4 106.7 118.5 113.4 105.6 119.2 113.9  81.4 115.6 121.7 111.0 105.2
2003 106.9 105.4 117.1 112.0 101.5 115.2 111.2  75.7 117.5 122.4 107.8 109.3
2004 104.7 106.7 122.8 112.7 104.5 126.5 111.1  79.7 121.9 118.8 112.2 112.6
2005 107.6 106.3 118.8 113.7 109.7 125.0 106.4  81.7 123.0 115.1 115.5 111.6
2006 108.8 105.9 124.8 108.0 113.1 126.7 108.7  84.1 121.0 121.5 116.6 108.2
2007 111.5 109.6 124.0 111.7 111.7 126.6 116.6  87.1 117.3 127.2 118.0 106.5
2008 113.2 114.4 117.5 120.2 107.6 121.4 115.1  78.8 118.9 118.8  99.5  99.2
2009  90.7  90.1 100.0  93.9  88.3 105.5  97.9  70.3 106.5 104.7  99.0  95.2
2010  90.3  93.1 109.5 100.4  95.5 111.8 100.8  74.5 109.0 105.0 102.7 101.9
2011  99.0 101.6 115.3 101.6 110.1 108.5 101.0  78.3 110.0 106.4 106.3 100.2
2012  99.3  99.9 110.3  99.8  96.1 108.5 103.8  78.8 102.9 107.6 101.9  91.5
2013  96.3  95.7 103.9 103.4  96.2 105.5 105.2  73.4 103.3 109.3  99.0  94.6
2014  96.8  97.1 104.9 102.9  92.4 104.8 103.3  71.5 107.2 107.7  95.5  98.4
2015  94.6  96.1 108.4 102.4  91.0 111.8 101.1  76.1 109.1 107.7 101.8  99.9
2016  95.1  99.6 108.3 103.3  98.6 110.3  95.8  79.7 107.9 103.4 104.7  99.7
2017  98.9  97.1 114.5  98.0 102.7 111.4  99.5  81.4 108.2 113.2 111.0  99.5
2018 101.8  98.6 112.9 103.0  98.7 112.8 106.6  82.6 104.7 116.0 109.6  97.6
2019 103.8 102.0 111.6 107.2 105.2 106.0 109.8  78.8 109.0 116.5 104.0  97.8
2020 101.0 100.1  91.8  66.7  73.7  98.2  97.4  71.7 104.7 106.7 101.6  96.6

attr(,"class")
[1] "JD3_TS"

$domainSpec
Specification

Series
Serie span: All 
Preliminary Check: Yes

Estimate
Model span: All 

Tolerance: 1e-07

Transformation
Function: AUTO
AIC difference: -2
Adjust: NONE

Regression
Calendar regressor: WorkingDays
with Leap Year: Yes
AutoAdjust: TRUE
Test: REMOVE

Easter: No

Pre-specified outliers: 0
Ramps: No

Outliers
Detection span: All 
Outliers type: 
    - AO, critical value : 0 (Auto)
    - LS, critical value : 0 (Auto)
    - TC, critical value : 0 (Auto)
TC rate: 0.7 (Auto)
Method: ADDONE (Auto)

ARIMA
SARIMA model: (0,1,1) (0,1,1)

SARIMA coefficients:
 theta(1) btheta(1) 
        0         0 

Specification X11
Seasonal component: Yes
Length of the Henderson filter: 0
Seasonal filter: FILTER_MSR
Boundaries used for extreme values correction : 
     lower_sigma:  1.5 
     upper_sigma:  2.5
Nb of forecasts: -1
Nb of backcasts: 0
Calendar sigma: NONE

Benchmarking
Is enabled: No

$estimationSpec
Specification

Series
Serie span: All 
Preliminary Check: Yes

Estimate
Model span: All 

Tolerance: 1e-07

Transformation
Function: AUTO
AIC difference: -2
Adjust: NONE

Regression
Calendar regressor: WorkingDays
with Leap Year: Yes
AutoAdjust: TRUE
Test: REMOVE

Easter: No

Pre-specified outliers: 0
Ramps: No

Outliers
Detection span: All 
Outliers type: 
    - AO, critical value : 0 (Auto)
    - LS, critical value : 0 (Auto)
    - TC, critical value : 0 (Auto)
TC rate: 0.7 (Auto)
Method: ADDONE (Auto)

ARIMA
SARIMA model: (0,1,1) (0,1,1)

SARIMA coefficients:
 theta(1) btheta(1) 
        0         0 

Specification X11
Seasonal component: Yes
Length of the Henderson filter: 0
Seasonal filter: FILTER_MSR
Boundaries used for extreme values correction : 
     lower_sigma:  1.5 
     upper_sigma:  2.5
Nb of forecasts: -1
Nb of backcasts: 0
Calendar sigma: NONE

Benchmarking
Is enabled: No

$pointSpec
Specification

Series
Serie span: All 
Preliminary Check: Yes

Estimate
Model span: All 

Tolerance: 1e-07

Transformation
Function: LEVEL
AIC difference: -2
Adjust: NONE

Regression
Calendar regressor: WorkingDays
with Leap Year: Yes
AutoAdjust: FALSE
Test: NO

Easter: No

Pre-specified outliers: 4
    - LS (2008-11-01), coefficient: -12.7560067943903 (ESTIMATED)
    - AO (2011-05-01), coefficient: 13.1047116166881 (ESTIMATED)
    - AO (2020-03-01), coefficient: -20.3005616124476 (ESTIMATED)
    - TC (2020-04-01), coefficient: -36.4047650490134 (ESTIMATED)
Ramps: No

Outliers
Is enabled: No

ARIMA
SARIMA model: (2,1,1) (0,1,1)

SARIMA coefficients:
   phi(1)    phi(2)  theta(1) btheta(1) 
  0.02527   0.10075  -0.59512  -0.71292 

Specification X11
Seasonal component: Yes
Length of the Henderson filter: 0
Seasonal filter: FILTER_MSR
Boundaries used for extreme values correction : 
     lower_sigma:  1.5 
     upper_sigma:  2.5
Nb of forecasts: -1
Nb of backcasts: 0
Calendar sigma: NONE

Benchmarking
Is enabled: No

$results
Model: X-13
Log-transformation: no 
SARIMA model: (2,1,1) (0,1,1)

SARIMA coefficients:
   phi(1)    phi(2)  theta(1) btheta(1) 
  0.02527   0.10075  -0.59512  -0.71292 

Regression model:
             td              lp LS (2008-11-01) AO (2011-05-01) AO (2020-03-01) 
         0.7032          2.1641        -12.7560         13.1047        -20.3006 
TC (2020-04-01) 
       -36.4048 

 Seasonal filter: S3X3;  Trend filter: H-13 terms
 M-Statistics: q Good (0.316); q-m2 Good (0.346)
 QS test on SA: Good (1.000);  F-test on SA: Good (0.646)

For a more detailed output, use the 'summary()' function.

# Pour extraire toutes les séries brutes il faut faire une boucle :
all_y <- lapply(
    lapply(
        seq_len(rjd3workspace::sap_sai_count(jsap1)),
        rjd3workspace::.jsap_sai, jsap = jsap1
    ),
    rjd3workspace::get_raw_data
)

1.3 Manipuler les objets Java

La manipulation des objets Java se fait de la même façon qu’avec RJDemetra : .jx13() pour estimer les modèles (nom de la fonction légèrement différente), rjd3toolkit::dictionary() pour connaitre l’ensemble des objets exportables et rjd3toolkit::result() et rjd3toolkit::user_defined() pour exporter des objets.

Exercice

Créer un modèle à partir de la fonction .jx13() et la spécification sans effet graduel de pâques calculée dans la section 1.1.

Solution

myjsa <- rjd3x13::.jx13(ipi_c_eu[, "FR"], spec_sans_easter_v3)
rjd3toolkit::result(myjsa, "sa")

           Jan       Feb       Mar       Apr       May       Jun       Jul
1990  93.24955  95.30940  94.71796  93.81684  93.54929  92.92520  94.21426
1991  93.59169  92.31215  92.64331  92.17822  87.64814  94.71770  92.26714
1992  91.58585  91.95491  92.44861  92.08835  91.15320  90.56855  89.24060
1993  87.99561  87.69727  86.39964  86.75333  87.31361  86.83199  86.72062
1994  89.29811  87.67745  87.97237  89.30396  91.16780  91.47286  90.78785
1995  95.25155  94.23926  93.90860  94.00837  92.19870  91.94726  93.82107
1996  93.17973  92.48812  93.89705  93.13793  91.31516  94.72025  94.11951
1997  93.26078  95.70964  94.51600 101.81259  96.59510  98.42946  98.22737
1998 102.21091 103.17118 101.63558 102.93616 103.62068 102.28550 103.77953
1999 103.79658 102.78743 103.69074 104.14106 106.09004 107.66846 105.67101
2000 110.10497 109.06807 110.12489 109.82704 114.79035 106.77270 110.36351
2001 112.60460 113.23756 114.26231 110.62978 111.90510 111.50689 112.11052
2002 110.38323 110.96526 110.79945 111.78627 107.92165 110.99153 111.32654
2003 108.79624 109.74015 109.32134 109.98806 106.03352 104.23265 108.95547
2004 108.85117 110.81356 110.13728 110.54560 111.16466 112.78890 111.54538
2005 113.93316 110.84172 106.41816 114.18867 113.55450 111.11371 109.26845
2006 112.42833 110.33885 112.86559 111.31529 114.27002 113.01129 111.32242
2007 112.55694 113.84939 115.05181 112.93951 112.66924 115.79095 116.41635
2008 114.19262 115.43079 111.37148 119.20887 111.21601 111.20490 111.96753
2009  94.06595  93.66765  91.59029  92.93345  94.74414  93.45541  94.24009
2010  96.01716  96.45680  98.71937  99.12528 102.35680 100.33373  99.14123
2011 104.73308 104.88254 104.41611 102.25226 114.86761  97.50904 101.41146
2012 102.69646 100.10116 101.67631  99.91395  98.94360 100.23048 101.39058
2013  97.38049  98.67284  97.59268 100.48330  99.43876  99.96091 100.12816
2014  98.00446 100.13410  98.51684  99.69373  98.14253  97.00824  98.38996
2015  98.43075  99.49303  99.43229  99.18875  98.89465 101.71133  96.50599
2016 101.43553 100.46907  96.73895 102.94029 103.57872 100.21169  96.40288
2017 102.83484 101.21980 102.93078 100.60288 104.65915 101.21299 100.09477
2018 103.34498 102.86547 103.87513 103.61900 100.09104 104.87134 104.57061
2019 105.55319 106.23769 105.20469 105.62380 106.14533 100.29916 105.00063
2020 102.89360 103.74645  83.07039  65.24951  79.29577  87.38128  92.35517
           Aug       Sep       Oct       Nov       Dec
1990  92.38173  93.51234  93.42343  92.18042  90.41279
1991  92.42873  92.22349  92.01173  92.00186  90.59538
1992  93.45078  89.51221  89.56165  89.95835  88.19546
1993  89.35917  86.63394  86.80871  84.00049  87.21539
1994  91.47220  91.08205  91.83094  92.70377  95.22430
1995  92.43935  93.94802  92.67708  93.29762  94.39315
1996  94.63683  93.94307  93.57709  93.75204  93.25425
1997 101.68755  99.71695 102.44668 100.79734 101.82217
1998 103.67043 103.81100 103.12842 104.36147 102.98643
1999 105.51412 107.47894 108.01349 107.32689 110.17856
2000 109.19601 111.03577 111.63442 113.48334 113.71019
2001 112.97886 110.95332 111.23200 110.59473 107.49517
2002 111.59373 109.55119 109.54524 110.06801 106.01718
2003 108.67882 108.40221 110.46920 109.44340 108.13524
2004 110.49903 112.36423 111.93619 108.93019 111.82427
2005 110.17307 113.24026 108.51104 112.00879 113.43584
2006 112.19306 113.67537 112.74018 112.68488 112.70890
2007 114.35913 112.60498 116.41114 113.72338 111.21506
2008 109.87755 109.64124 108.31326  99.94885  99.31443
2009 100.20028  97.57161  97.15487  97.05587  95.50465
2010 100.91566 100.33425 100.28702  98.54035 102.42445
2011 101.41452 101.52787 102.02302 102.59491 103.34800
2012 101.29425  99.61270  98.34668  98.61297  97.28502
2013  97.82035  97.74060 100.19255  98.38971  97.85776
2014  98.04534  99.36208  98.45915  97.19073  99.09353
2015 102.42001 101.51956 100.68665 100.72423 100.52605
2016 101.09905 100.59846  98.34729 100.76466 102.91719
2017 102.74040 103.61311 105.30626 106.79695 105.62276
2018 103.94813 102.71388 105.22599 105.29463 104.37228
2019 102.51548 104.62023 105.52898 102.28202 102.79014
2020  97.79844  97.77412  98.09488 100.12642  99.63962

Pour les révisions, la fonction rjd3x13::x13_revisions() facilite grandement la tâche. L’historique des révisions peut s’exporter à partir de trois paramètres :

data_ids qui permet d’exporter des statistiques ;
ts_ids qui permet d’exporter des estimations d’une composante à une certaine date ;
cmp_ids qui permet l’ensemble des estimations d’une composante à une ensemble de dates.

data_ids <- list(
    # Export du coefficient du premier régresseur jours ouvrables
    list(start = "2005-01-01", id = "regression.td(1)"))
ts_ids <- list(
    # Export de l'historique des estimations de la séries SA de janv 2010 à partir de janvier 2010
    list(period = "2010-01-01", start = "2010-01-01", id = "sa"))
cmp_ids <- list(
    # Export de l'ensemble des tendances estimées entre janv 2010 et dec 2014
    list(start = "2010-01-01", end = "2014-12-01", id = "t"))
rh <- rjd3x13::x13_revisions(
    sa2_jd3$result$preadjust$a1,
    spec = sa2_jd3$result_spec, 
    data_ids, ts_ids, cmp_ids)
plot(rh$data$`regression.td(1)`,
     ylab = "Coef", main = "Coef td")

plot(rh$series$sa,
     ylab = NULL, main = "SA de janv 2010")

ncol(rh$components$t) # nombre de séries exportées

[1] 60

2 RJDemetra v2

Commencez par charger le package RJDemetra :

library(RJDemetra)

2.1 Créer une specification

Dans les prochains exercices, la série utilisée sera ipi_c_eu[, "FR"] qui est l’IPI français. Vous pouvez bien sûr adapter le code pour utiliser vos propres séries. Les fonctions utilisées seront x13(), x13_spec(), regarima_x13, regarima_x13_spec ou regarima. Le détail des spécifications pré-définies par JDemetra+ sont disponibles ici.

Exercice

Faire la désaisonnalisation d’une série avec X-13 avec la spécification suivante :

détection automatique du schéma de décomposition, des outliers et du modèle ARIMA ;
une correction des jours ouvrables “working days” et un effet graduel de Pâques.

Faire ensuite un graphique avec la série brute et la série désaisonnalisée.

Indice

Utiliser la spécification RSA4c pour la désaisonnalisation.

Si le modèle créé s’appelle mysa, regarder les valeurs de mysa$final, mysa$final$series et mysa$final$forecasts.

Solution

library(RJDemetra)
mysa <- x13(ipi_c_eu[, "FR"], spec = "RSA4c")
mysa

RegARIMA
y = regression model + arima (2, 1, 1, 0, 1, 1)
Log-transformation: no
Coefficients:
          Estimate Std. Error
Phi(1)     0.05291      0.108
Phi(2)     0.18672      0.074
Theta(1)  -0.52137      0.103
BTheta(1) -0.66132      0.042

             Estimate Std. Error
Week days      0.6927      0.031
Leap year      2.0903      0.694
Easter [1]    -2.5476      0.442
TC (4-2020)  -35.6481      2.092
AO (3-2020)  -21.1492      2.122
AO (5-2011)   13.1869      1.810
LS (11-2008)  -9.2744      1.758
LS (1-2009)   -7.2838      1.756


Residual standard error: 2.193 on 346 degrees of freedom
Log likelihood = -795.1, aic =  1616 aicc =  1617, bic(corrected for length) = 1.767



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


Final
Last observed values
             y        sa        t           s           i
Jan 2020 101.0 102.87273 103.0457  -1.8727280  -0.1730003
Feb 2020 100.1 103.69025 103.0626  -3.5902540   0.6276448
Mar 2020  91.8  82.69170 103.2654   9.1083000 -20.5736602
Apr 2020  66.7  66.55184 103.6945   0.1481625 -37.1426277
May 2020  73.7  79.28883 104.1379  -5.5888279 -24.8490764
Jun 2020  98.2  87.35362 104.4539  10.8463757 -17.1002284
Jul 2020  97.4  92.26057 104.5518   5.1394323 -12.2912806
Aug 2020  71.7  97.54392 104.3369 -25.8439193  -6.7929846
Sep 2020 104.7  97.75728 103.8361   6.9427184  -6.0788659
Oct 2020 106.7  97.87016 103.1969   8.8298396  -5.3267823
Nov 2020 101.6 100.01475 102.6601   1.5852540  -2.6453446
Dec 2020  96.6  99.61740 102.4081  -3.0173983  -2.7907314

Forecasts:
               y_f     sa_f      t_f         s_f        i_f
Jan 2021  94.29728 101.0937 102.3876  -6.7963909 -1.2939310
Feb 2021  97.89298 101.6869 102.4525  -3.7939488 -0.7655396
Mar 2021 113.65190 102.1478 102.4593  11.5041248 -0.3115701
Apr 2021 102.34532 102.1806 102.3419   0.1647274 -0.1612727
May 2021  96.14552 101.6436 102.1712  -5.4980759 -0.5276194
Jun 2021 112.15785 101.2156 102.0241  10.9422549 -0.8085520
Jul 2021 104.38495 101.5439 101.9635   2.8410057 -0.4195241
Aug 2021  79.02901 102.3820 102.0412 -23.3530134  0.3408366
Sep 2021 109.40288 102.3705 102.1977   7.0323700  0.1728109
Oct 2021 108.22472 101.8558 102.3655   6.3689511 -0.5096780
Nov 2021 106.22014 102.4265 102.5266   3.7936049 -0.1000754
Dec 2021  99.67523 102.9462 102.6812  -3.2709156  0.2649189


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


Additional output variables

y <- mysa$final$series[,"y"]
# De façon équivalente :
y <- get_ts(mysa)
sa <- mysa$final$series[,"sa"]
plot(y)
lines(sa, col = "red")

# ou on peut directement utiliser les fonctions de RJDemetra :
plot(mysa, first_date = 2000, #Pour n'afficher le graphique qu'à partir de 200
     type_chart = "sa-trend" # Pour faire le graphique avec y, sa et tendance
)

Pour des graphiques ggplot2, on peut également utiliser le package ggdemetra :

library(ggdemetra)
# y <- ggdemetra::raw(mysa)
# sa <- ggdemetra::seasonaladj(mysa)
p_sa <- 
    ggplot(data = ggdemetra::ts2df(y), 
           mapping = aes(x = date, y = y)) +
    geom_line(color = "#F0B400") +
    labs(title = "Y, Sa, Trend",
         x = NULL, y = NULL) +
    geom_sa(component = "y_f", linetype = 2,
            spec = x13_spec(mysa), frequency = 12,
            color = "#F0B400") + 
    geom_sa(component = "sa", color = "#155692") +
    geom_sa(component = "sa_f", color = "#155692", linetype = 2)+ 
    geom_sa(component = "t", color = "#1E6C0B") +
    geom_sa(component = "t_f", color = "#1E6C0B", linetype = 2) + 
    theme_bw()
p_sa

Exercice

Modifier le modèle précédent pour enlever l’effet graduel de Pâques.

Solution

spec_sans_easter <- x13_spec(mysa,
                             easter.enabled = FALSE)
mysa2 <- x13(ipi_c_eu[, "FR"], spec_sans_easter)
mysa2$regarima

y = regression model + arima (2, 1, 1, 0, 1, 1)
Log-transformation: yes
Coefficients:
          Estimate Std. Error
Phi(1)     0.05032      0.118
Phi(2)     0.09575      0.082
Theta(1)  -0.55559      0.110
BTheta(1) -0.73033      0.039

              Estimate Std. Error
Week days     0.007179      0.000
AO (5-2011)   0.124297      0.018
LS (11-2008) -0.086926      0.017
LS (1-2009)  -0.071005      0.017


Residual standard error: 0.02133 on 338 degrees of freedom
Log likelihood =   838, aic =  1548 aicc =  1548, bic(corrected for length) = -7.561

Exercice

Calculer les p-valeurs associées au modèle Reg-ARIMA de la précédente spécification.

Indice

Récupérer le modèle Reg-ARIMA et utiliser la fonction summary().

Solution

reg_sum <- summary(mysa2$regarima)
reg_sum

y = regression model + arima (2, 1, 1, 0, 1, 1)

Model: RegARIMA - X13
Estimation span: from 1-1990 to 12-2019
Log-transformation: yes
Regression model: no mean, trading days effect(1), no leap year effect, no Easter effect, outliers(3)

Coefficients:
ARIMA: 
          Estimate Std. Error  T-stat Pr(>|t|)    
Phi(1)     0.05032    0.11801   0.426    0.670    
Phi(2)     0.09575    0.08224   1.164    0.245    
Theta(1)  -0.55559    0.10956  -5.071 6.45e-07 ***
BTheta(1) -0.73033    0.03938 -18.545  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Regression model: 
               Estimate Std. Error T-stat Pr(>|t|)    
Week days     0.0071788  0.0003171 22.637  < 2e-16 ***
AO (5-2011)   0.1242970  0.0180390  6.890 2.63e-11 ***
LS (11-2008) -0.0869263  0.0171661 -5.064 6.69e-07 ***
LS (1-2009)  -0.0710049  0.0171449 -4.141 4.34e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.02133 on 338 degrees of freedom
Log likelihood =   838, aic =  1548, aicc =  1548, bic(corrected for length) = -7.561

On peut récupérer ces valeurs en exploitant l’objet reg_sum :

reg_sum$coefficients

$arima
             Estimate Std. Error     T-stat     Pr(>|t|)
Phi(1)     0.05031546 0.11801332   0.426354 6.701142e-01
Phi(2)     0.09575287 0.08224054   1.164303 2.451013e-01
Theta(1)  -0.55559277 0.10955864  -5.071191 6.449100e-07
BTheta(1) -0.73033107 0.03938071 -18.545401 0.000000e+00

$regression
                 Estimate   Std. Error    T-stat     Pr(>|t|)
Week days     0.007178836 0.0003171301 22.636876 0.000000e+00
AO (5-2011)   0.124296961 0.0180390210  6.890449 2.628564e-11
LS (11-2008) -0.086926347 0.0171661398 -5.063826 6.685537e-07
LS (1-2009)  -0.071004892 0.0171448840 -4.141462 4.338403e-05

$fixed_out
NULL

$fixed_var
NULL

2.2 Créer un workspace

Dans cette partie nous allons créer un workspace depuis R. Pour cela les fonctions qui peuvent être utilisées sont new_workspace(), load_workspace(), new_multiprocessing(), add_sa_item(), save_workspace(), compute(), get_object(), get_name(), get_ts() ou count().

Exercice

Créer un workspace qui va contenir une série désaisonnalisée selon 3 spécifications différentes.

Solution

jws <- new_workspace()
new_multiprocessing(jws, "MP-1")
add_sa_item(jws, "MP-1", mysa, "X13 avec Pâques")
add_sa_item(jws, "MP-1", mysa2, "X13 sans Pâques")
add_sa_item(jws, "MP-1", tramoseats(ipi_c_eu[, "FR"]), "TRAMO-SEATS")
save_workspace(jws, "mon_premier_workspace.xml")

Exercice

Importer le workspace précédent et récupérer :

Le nom du premier multi-processing
Le nombre de modèles dans ce premier multi-processing
L’ensemble des séries brutes
Le 2ème modèle

Solution

jws <- load_workspace("mon_premier_workspace.xml")
compute(jws)
count(jws) # Nombre de multiprocessing

[1] 1

jmp1 <- get_object(jws, 1) # Le premier multiprocessing
get_name(jmp1)

[1] "MP-1"

count(jmp1)

[1] 3

all_y <- get_ts(jmp1) # toutes les séries brutes
model2 <- get_object(jmp1, 2) # On récupère l'objet associé au 2ème modèle
get_model(model2, jws)

RegARIMA
y = regression model + arima (2, 1, 1, 0, 1, 1)
Log-transformation: yes
Coefficients:
          Estimate Std. Error
Phi(1)     0.05032      0.118
Phi(2)     0.09575      0.082
Theta(1)  -0.55559      0.110
BTheta(1) -0.73033      0.039

              Estimate Std. Error
Week days     0.007179      0.000
AO (5-2011)   0.124297      0.018
LS (11-2008) -0.086926      0.017
LS (1-2009)  -0.071005      0.017


Residual standard error: 0.02133 on 338 degrees of freedom
Log likelihood =   838, aic =  1548 aicc =  1548, bic(corrected for length) = -7.561



Decomposition
Monitoring and Quality Assessment Statistics:
      M stats
M(1)    0.080
M(2)    0.044
M(3)    0.925
M(4)    0.324
M(5)    1.018
M(6)    0.172
M(7)    0.075
M(8)    0.207
M(9)    0.069
M(10)   0.229
M(11)   0.207
Q       0.311
Q-M2    0.344

Final filters: 
Seasonal filter:  3x5
Trend filter:  13 terms Henderson moving average


Final
Last observed values
             y       sa        t         s         i
Jan 2019 103.9 104.9497 104.8872 0.9899977 1.0005966
Feb 2019 101.9 106.2130 105.1277 0.9593928 1.0103234
Mar 2019 111.0 104.7070 105.2207 1.0601007 0.9951182
Apr 2019 107.4 105.0688 105.1200 1.0221871 0.9995130
May 2019 105.5 108.7078 104.9337 0.9704918 1.0359658
Jun 2019 105.8 101.5038 104.6881 1.0423256 0.9695826
Jul 2019 110.1 105.4918 104.3766 1.0436828 1.0106850
Aug 2019  78.7 102.6361 104.0097 0.7667871 0.9867928
Sep 2019 108.5 104.0439 103.5870 1.0428287 1.0044106
Oct 2019 116.8 104.5857 103.1442 1.1167879 1.0139754
Nov 2019 103.8 101.7786 102.7722 1.0198610 0.9903319
Dec 2019  97.7 101.5883 102.5735 0.9617252 0.9903948

Forecasts:
               y_f     sa_f      t_f       s_f       i_f
Jan 2020 101.86190 102.9804 102.5981 0.9891385 1.0037265
Feb 2020 100.79119 103.1440 102.8002 0.9771887 1.0033443
Mar 2020 111.43956 102.5781 103.0224 1.0863878 0.9956869
Apr 2020 105.54428 103.3790 103.2028 1.0209447 1.0017079
May 2020  96.34392 104.3346 103.2155 0.9234131 1.0108422
Jun 2020 112.45079 102.5347 103.0496 1.0967099 0.9950029
Jul 2020 106.92345 102.3231 102.8558 1.0449593 0.9948203
Aug 2020  76.61765 102.5744 102.7208 0.7469472 0.9985747
Sep 2020 110.47115 103.3476 102.6644 1.0689277 1.0066554
Oct 2020 111.53872 102.2010 102.6987 1.0913668 0.9951530
Nov 2020 104.52886 102.5479 102.7719 1.0193176 0.9978206
Dec 2020 101.73775 103.4319 102.9459 0.9836210 1.0047204


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.984
 Seasonal     62.978
 Irregular     0.998
 TD & Hol.     2.273
 Others       33.676
 Total       101.909

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

Residual seasonality tests
                                      P.value
 qs test on sa                          1.000
 qs test on i                           0.965
 f-test on sa (seasonal dummies)        0.625
 f-test on i (seasonal dummies)         0.474
 Residual seasonality (entire series)   0.779
 Residual seasonality (last 3 years)    0.878
 f-test on sa (td)                      0.076
 f-test on i (td)                       0.385


Additional output variables

2.3 Manipuler les objets Java

L’objectif de cette partie est de manipuler la fonction jx13() pour gagner en temps de calcul.

Exercice

Créer un modèle à partir de la fonction jx13() et la spécification sans effet graduel de pâques crée dans la section 2.1.

Solution

myjsa <- jx13(ipi_c_eu[, "FR"], spec_sans_easter)
get_indicators(myjsa, "sa")

$sa
           Jan       Feb       Mar       Apr       May       Jun       Jul
1990  93.27297  96.29529  94.47992  93.73334  93.94507  92.83948  94.28898
1991  93.61809  93.12145  92.56800  92.23966  88.02740  94.55282  92.48219
1992  91.62891  91.43164  92.38484  92.19526  91.08704  90.59578  89.51814
1993  87.30634  87.99736  86.85921  86.97894  86.99296  87.29234  86.76534
1994  88.69622  87.77278  88.37550  89.30617  91.27332  91.48108  90.65425
1995  95.09916  94.52731  93.92731  93.90941  92.37904  92.04264  93.77960
1996  93.06670  91.54362  94.26169  93.30918  91.32878  94.96354  94.26156
1997  93.06742  95.67044  95.13059 101.85852  96.43204  98.44794  98.23010
1998 102.13804 103.25275 102.01932 102.89690 103.62209 101.81934 103.60992
1999 103.71964 102.64444 103.76759 104.04665 106.13300 106.82584 105.69365
2000 110.38789 107.38428 109.60527 110.17270 114.61835 106.35539 110.67086
2001 112.45025 113.34404 113.88119 110.73914 111.61677 111.49774 112.13527
2002 110.26696 110.99285 111.13280 111.58735 107.54579 111.87216 110.92070
2003 108.72641 109.77205 109.75157 109.85823 105.81090 105.36978 108.56247
2004 109.12622 109.43552 109.52322 110.42603 111.47708 112.70525 111.38248
2005 114.72877 111.07384 106.23729 114.40231 113.83685 111.14636 109.33782
2006 112.84431 110.57123 111.88826 111.79247 114.34733 112.62975 111.48879
2007 112.53191 114.29580 114.23700 113.06787 112.93424 115.43386 116.32302
2008 114.06045 113.94480 111.07963 118.83142 111.92185 110.96502 111.54674
2009  93.60983  93.32624  92.17978  92.82197  94.61843  94.29186  94.42552
2010  95.51379  96.12075  98.46494  99.01596 102.81820 100.22951  99.30088
2011 104.88512 104.72507 103.63569 102.25128 116.10802  97.55663 101.60353
2012 102.83064  98.48657 101.60323  99.95552  99.30289 100.25687 101.48753
2013  97.45647  98.45690  97.93752 100.35136  99.79029 100.14504 100.08302
2014  97.89835 100.00657  98.78845  99.56271  98.29814  97.26484  98.38400
2015  98.27255  99.53811  99.50288  99.07983  99.28072 101.49169  96.66780
2016 101.32397  99.22467  96.95194 102.97678 105.03035 100.37169  96.68839
2017 102.72859 101.41668 102.63069 100.60801 106.58076 101.54655 100.67684
2018 103.19713 103.14126 103.70513 103.35321 102.04279 105.47625 105.16178
2019 104.94974 106.21301 104.70704 105.06883 108.70777 101.50380 105.49182
           Aug       Sep       Oct       Nov       Dec
1990  92.81532  93.42471  93.21965  92.25071  90.06120
1991  91.72254  92.21808  91.96428  92.02285  90.46195
1992  91.95223  89.69252  89.81277  90.02832  88.33305
1993  87.30196  87.03081  87.21300  84.45288  87.46656
1994  91.37247  91.19741  91.98386  92.78677  95.17153
1995  92.38821  93.99473  92.85706  93.36671  94.38969
1996  93.75684  94.14529  93.72433  93.79062  93.43034
1997 101.81145  99.64928 101.86478 100.90774 101.70138
1998 103.62658 103.60456 103.09096 104.46065 102.99964
1999 106.07083 107.12621 108.29389 107.17738 110.10121
2000 111.01860 110.92744 111.42763 113.10692 114.51591
2001 115.50617 111.31601 110.65583 110.20538 108.13496
2002 112.75124 109.55689 109.18036 109.90888 106.38982
2003 107.51695 108.10694 110.01758 109.50133 108.16481
2004 110.31329 111.75188 112.44658 108.46909 111.73847
2005 110.00423 112.59573 109.23754 111.65979 113.81029
2006 112.62602 113.48750 112.76155 112.55067 113.59525
2007 115.82841 112.91803 115.49734 113.69439 112.28630
2008 109.36385 109.09852 108.03726 100.63716  99.86454
2009  96.91866  97.94748  97.85917  97.60414  96.07651
2010  99.54513 100.43464 100.75477  98.83302 102.97020
2011 101.60617 101.45456 102.21445 102.59847 103.85017
2012 101.83135  99.96547  98.23490  98.64730  97.23489
2013  96.84385  97.95614  99.86479  98.42704  97.79728
2014  96.12996  99.28707  98.34888  97.27581  99.01780
2015 102.05354 101.35873 100.66252 100.76102 100.45864
2016 101.62022 100.64268  98.76241 100.83623 102.82647
2017 103.73341 103.82246 105.29816 106.59767 105.53516
2018 105.03058 103.00706 104.49092 104.89205 103.66170
2019 102.63605 104.04393 104.58566 101.77858 101.58827

Exercice

Toujours avec la même spécification, extraire les révisions de la séries désaisonnalisée du point de janvier 2005 (i.e. : série désaisonnalisée lorsqu’on a les données jusqu’en janvier 2005, puis jusqu’en février 2005, etc.).

Indice

Récupérer l’ensemble des dates de fin d’estimation avec la fonction window(time(ipi_c_eu[, "FR"]), start = 2005).

Solution

dates <- as.numeric(window(time(ipi_c_eu[, "FR"]), start = 2005))
estimations <- sapply(dates, function(last_date_estimation){
    myjsa <- jx13(window(ipi_c_eu[, "FR"], end = last_date_estimation), spec_sans_easter)
    sa <- get_indicators(myjsa, "sa")$sa
    window(sa, start = 2005, end = 2005) # Pour ne récupérer que la valeur en 2005
})
estimations <- ts(estimations, start = 2005, frequency = 12)
plot(estimations)