3 - Qualité du préajustement sous R

Désaisonnaliser une série temporelle

Auteur

Alain Quartier-la-Tente

L’objectif de ce TP est d’apprendre à vérifier la qualité du pré-ajustement dans 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")
}

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.855

library(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 ' ' 1

Concernant 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.646

La 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_defined

Names 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.169

BIC(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 ?