2 - Correction des jours ouvrables

Désaisonnaliser une série temporelle

Auteur

Alain Quartier-la-Tente

L’objectif de ce TP est d’apprendre à créer des régresseurs jours ouvrables personnalisés, à les utiliser dans RJDemetra/JDemetra+ et à vérifier la qualité de la correction.

Pour installer tous les packages utiles de ce TP, lancer le programme :

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

Pour créer des régresseurs jours ouvrables personnalisés, deux solutions :

  1. Le faire depuis JDemetra+, en créant un calendrier personnalisé puis en exportant les régresseurs. Voir par exemple la documentation de JDemetra+ ici et là.

  2. Créer le calendrier depuis R à l’aide du package rjd3toolkit.

Dans ce TP, nous nous intéresserons uniquement à la seconde option. En effet, le package R est plus flexible et permet de créer des régresseurs moins classiques que les working days et trading days.

1 Création de régresseurs JO avec rjd3toolkit

1.1 Création d’un calendrier

Par défaut, les régresseurs jours ouvrables de JDemetra+ ne prennent pas en compte les spécificité calendaires des pays : on ne prend pas en compte les jours fériés. Pour les prendre en compte, il faut créer son propre calendrier où l’on supposera qu’un jour férié de plus dans le mois a le même effet qu’un dimanche.

library(rjd3toolkit)

Trois fonctions peuvent être utilisées pour ajouter des jours fériés :

  1. fixed_day() pour ajouter un jour férié qui tombe à date fixe. Par exemple, pour ajouter le 8 mai :
fixed_day(month =  5, day = 8)
  1. easter_day() pour ajouter un jour férié dont le jour dépend de Pâques : le paramètre offset permet de spécifier le nombre de jours avant (si négatif) ou après Pâques (si positif). Par exemple, pour ajouter la Pentecôte qui a lieu 60 jours après Pâques :
easter_day(offset = 60)
  1. single_day() pour ajouter un jour ferié qui n’a eu lieu qu’une seule fois.
single_day("1993-01-15")
  1. fixed_week_day() qui permet d’ajouter des jours fériés qui apparaissent certaines semaines de certains mois (par exemple le premier lundi du mois de septembre aux USA)
fixed_week_day(9, 1, 1) # first Monday(1) of September.
  1. special_day() qui permet d’ajouter des jours fériés par rapport à des dates déjà connues dans JDemetra+ (voir tableau ci-dessous). Comme pour la fonction easter_day(), le paramètre offset permet de spécifier la position du jour voulu par rapport rapport à la fête pré-spécifié (par défaut offset = 0, le jour férié coïncide avec le jour pré-spécifié). Par exemple, pour ajouter le nouvel an :
special_day("NEWYEAR")
Jours pré-spécifiés
Event Définition
NEWYEAR Fête fixe, 1er janvier.
SHROVEMONDAY Fête mobile, lundi avant le mecredi des cendres (48 jours avant pâques).
SHROVETUESDAY Fête mobile, mardi avant le mecredi des cendres (47 jours avant pâques).
ASHWEDNESDAY Fête mobile, 46 jours avant Pâques.
EASTER Fête mobile, Pâques, varie entre le 22 mars et le 25 avril.
MAUNDYTHURSDAY Fête mobile, le jeudi avant Pâques.
GOODFRIDAY Fête mobile, le vendredi avant Pâques.
EASTERMONDAY Fête mobile, le lendemain de Pâques.
ASCENSION Fête mobile, célébrée un jeudi, 40 jours après Pâques.
PENTECOST Fête mobile, 50 jours après Pâques.
CORPUSCHRISTI Fête mobile, 60 jours après Pâques.
WHITMONDAY Fête mobile, le jour après la Pentecôte.
MAYDAY Fête fixe, 1er mai.
ASSUMPTION Fête fixe, 15 août.
HALLOWEEN Fête fixe, 31 octobre.
ALLSAINTSDAY Fête fixe, 1er novembre.
ARMISTICE Fête fixe, 11 novembre.
CHRISTMAS Fête fixe, 25 décembre.
Exercice

Créer un calendrier qui contient tous les jours fériés de votre pays.

Deux exemples :

  1. Calendrier associé à la France :
FR <- list(
    special_day("NEWYEAR"),
    special_day("EASTERMONDAY"), # Lundi de Pâques
    special_day("MAYDAY"), # 1er mai
    special_day("ASCENSION"), # Jour de l'Ascension
    fixed_day(5, 8),
    special_day("WHITMONDAY"), # Lundi de Pentecôte
    fixed_day(7, 14),
    special_day("ASSUMPTION"), # Assomption
    special_day("ALLSAINTSDAY"), # Toussaint
    special_day("ARMISTICE")
)
CAL <- national_calendar(FR)
  1. Calendrier associé à la Macronia, la difficulté étant qu’il faut ajouter à la main des jours associés aux fêtes musulmanes
jours_macronia <- list(
  special_day("NEWYEAR"),
    special_day("EASTERMONDAY"), # Lundi de Pâques
  fixed_day(4, 4), # Jour de l'indépendance de la Macronia
  special_day("MAYDAY"), # 1er mai
  special_day("ASCENSION"), # Jour de l'Ascension
    special_day("WHITMONDAY"), # Lundi de Pentecôte
  special_day("ASSUMPTION"), # Assomption de Marie
    special_day("ALLSAINTSDAY"), # Toussaint
  special_day("CHRISTMAS") # Noël
)
# # Manque Début ramadan et jours décrétés
# # On récupère ces jours construisant un fichier Excel
# jours_mobiles <- readxl::read_excel("../data/DateFetesMusulmanes_Macronia.xlsx") |>
#   as.data.frame()
# jours_mobiles <- jours_mobiles[,-1, drop = FALSE] # on enlève l'année

jours_macronia <- c(
  jours_macronia,
  list(
    # DEBUT_RAMADAN,
single_day("2000-11-28"),
single_day("2001-11-17"),
single_day("2002-11-07"),
single_day("2003-10-27"),
single_day("2004-10-15"),
single_day("2005-10-05"),
single_day("2006-09-24"),
single_day("2007-09-14"),
single_day("2008-09-02"),
single_day("2009-08-22"),
single_day("2010-08-12"),
single_day("2011-08-01"),
single_day("2012-07-20"),
single_day("2013-07-10"),
single_day("2014-06-29"),
single_day("2015-07-18"),
single_day("2016-06-07"),
single_day("2017-05-27"),
single_day("2018-05-17"),
single_day("2019-05-06"),
single_day("2020-04-24"),
single_day("2021-04-14"),
single_day("2022-04-03"),
single_day("2023-03-23"),
single_day("2024-03-12")
  ))

CAL <- national_calendar(jours_macronia)

1.2 Création de régresseurs JO

Le modèle général de correction de jours ouvrables peut s’écrire de la façon suivante : \[ X_t = \sum_{i=1}^{7} \alpha_i N_{it} + \varepsilon_t \] Avec :

  • \(N_{it}\) le nombre de jours de lundis (\(i=1\)), …, dimanches et jours fériés (\(i=7\))

  • \(\alpha_i\) l’effet d’un jour de type \(i\)

Pour éviter les problèmes de multi-colinéarité, on réécrit le modèle en utilisant une modalité de référence (ici dimanche). On désaisonnalise également les régresseurs en enlevant la moyenne de long-terme : \[X_t = \sum_{i=1}^{6} \beta_i (N_{it} - N_{7t}) + \bar{\alpha} \underbrace{(N_t - \bar{N}_t)}_{LY_t} + \varepsilon_t\] Ce modèle peut être simplifié si en faisant des hypothèses sur les effets des jours ouvrés :

  • L’hypothèse working days correspond au cas où l’on suppose que tous les jours de la semaine (lundi à vendredi) ont le même effet (\(\alpha_1=\dots=\alpha_5\)), les samedis et les dimanches (et jours fériés) ont le même effet (\(\alpha_6=\alpha_7\)) et sont utilisés en tant que variable de contraste.

  • L’hypothèse trading days correspond au cas où l’on suppose que tous les jours ont un effet différent et les dimanches (et jours fériés) sont utilisés en tant que variable de constrate.

Sous JDemetra+ on ne peut utiliser que ces deux hypothèses mais rjd3toolkit permet de construire d’autres types de JO.

De manière plus générale, lorsque l’on utilise une variable de contraste, les régresseurs \(CJO_{t,i}\) associé au groupe \(i\) est calculé de la façon suivante : \[ CJO_{t,i} = \underbrace{\sum_{j\in\text{groupe }i}N_{jt}}_{ \text{nb de jours du groupe }i } - \frac{\sum_{j\in\text{groupe }i}1}{\sum_{j\in\text{groupe }0}1} \times \underbrace{\sum_{j\in\text{groupe }0}N_{jt}}_{ \text{nb de jours du groupe contraste} } \] Dans le cas working days, il y a 2 jours dans le groupe contraste (samedi et dimanche, \(\sum_{j\in\text{groupe }0}1=2\)) et 5 jours dans le groupe 1 (lundi à vendredi, \(\sum_{j\in\text{groupe }1}1=5\)). Au mois \(t\), le régresseurs JO type de jours est donc égal au nombre de jours de la semaine dans le mois, mois \(5/2\times\) nombre de jours de week-end.

Les régresseurs JO peuvent être créés à partir de 2 fonctions : htd() qui permet de les créer à partir d’un calendrier spécifique et td(). Dans ces fonctions, le paramètre le plus important est groups pour permet de faire des hypothèses sur les jours. C’est un vecteur de longueur 7 (le nombre de jours de la semaine) dont chaque élément indique à quel groupe le jour de la semaine associé correspond. La variable de contraste est associé au groupe 0.
Par exemple, groups = c(1,2,3,4,5,6,0) correspond au trading days et groups = c(1,1,1,1,1,0,0) correspond au working days.

Par exemple :

groups <- c(1, 2, 3, 4, 5, 6, 0)
frequency <- 12
start <- c(2000,1)
wkd <- calendar_td(CAL, frequency = frequency, start = start, length = 12*35,
                   groups = groups)
wkd <- ts(wkd, start = start, frequency = frequency)
Exercice

Comparer le régresseurs JO working days créé avec personnalisé et celui sans hypothèse sur les jours fériés (fonction td()).

Les régresseurs sont bien différents :

groups <- c(1, 1, 1, 1, 1, 0, 0)
frequency <- 12
start <- c(2000,1)
wkd <- calendar_td(CAL, frequency = frequency, start = start, length = 12*35,
                   groups = groups)
wkd_def <- td(frequency = frequency, start = start, length = 12*35,
              groups = groups)
round(wkd - wkd_def,1)
      Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
2000  2.5  0.0  0.7 -1.7  4.5 -5.5  0.0 -1.0  0.0  0.0 -4.5 -1.0
2001 -1.0  0.0  0.7 -1.7  1.0 -2.0  0.0 -1.0  0.0  0.0 -1.0 -1.0
2002 -1.0  0.0  0.7 -1.7 -2.5  1.5  0.0 -1.0  0.0  0.0 -4.5 -1.0
2003 -1.0  0.0  0.7 -1.7  1.0 -2.0  0.0 -1.0  0.0 -3.5  2.5 -1.0
2004 -1.0  0.0  0.7  1.8  1.0  1.5  0.0  2.5  0.0 -3.5 -1.0  2.5
2005  2.5  0.0 -2.8  1.8  1.0  1.5  0.0 -1.0  0.0 -3.5 -1.0  2.5
2006  2.5  0.0  0.7 -1.7  1.0 -2.0  0.0 -1.0  0.0  0.0 -1.0 -1.0
2007 -1.0  0.0  0.7 -1.7 -2.5  1.5  0.0 -1.0 -3.5  0.0 -1.0 -1.0
2008 -1.0  0.0 -2.8  1.8  1.0  1.5  0.0 -1.0 -3.5  0.0  2.5 -1.0
2009 -1.0  0.0  0.7  1.8  1.0 -2.0  0.0  2.5  0.0  0.0  2.5 -1.0
2010 -1.0  0.0  0.7  1.8  1.0  1.5  0.0 -1.0  0.0  0.0 -1.0  2.5
2011  2.5  0.0  0.7 -1.7  8.0 -5.5  0.0 -4.5  0.0  0.0 -1.0  2.5
2012  2.5  0.0  0.7 -1.7 -2.5  1.5 -3.5 -1.0  0.0  0.0 -1.0 -1.0
2013 -1.0  0.0  0.7 -1.7 -2.5  1.5 -3.5 -1.0  0.0  0.0 -1.0 -1.0
2014 -1.0  0.0  0.7 -1.7  1.0 -2.0  0.0 -1.0  0.0  0.0  2.5 -1.0
2015 -1.0  0.0  0.7  1.8 -2.5  1.5  0.0  2.5  0.0  0.0  2.5 -1.0
2016 -1.0  0.0 -2.8  1.8  1.0 -2.0  0.0 -1.0  0.0  0.0 -1.0  2.5
2017  2.5  0.0  0.7 -1.7  1.0 -2.0  0.0 -1.0  0.0  0.0 -1.0 -1.0
2018 -1.0  0.0  0.7 -1.7 -6.0  1.5  0.0 -1.0  0.0  0.0 -1.0 -1.0
2019 -1.0  0.0  0.7 -1.7 -2.5 -2.0  0.0 -1.0  0.0  0.0 -1.0 -1.0
2020 -1.0  0.0  0.7 -1.7  1.0 -2.0  0.0  2.5  0.0  0.0  2.5 -1.0
2021 -1.0  0.0  0.7 -1.7  1.0  1.5  0.0  2.5  0.0  0.0 -1.0  2.5
2022  2.5  0.0  0.7 -1.7  4.5 -2.0  0.0 -1.0  0.0  0.0 -1.0  2.5
2023  2.5  0.0 -2.8 -1.7 -2.5  1.5  0.0 -1.0  0.0  0.0 -1.0 -1.0
2024 -1.0  0.0 -2.8 -1.7 -2.5  1.5  0.0 -1.0  0.0  0.0 -1.0 -1.0
2025 -1.0  0.0  0.7 -1.7  1.0 -2.0  0.0 -1.0  0.0  0.0  2.5 -1.0
2026 -1.0  0.0  0.7  1.8 -2.5  1.5  0.0  2.5  0.0  0.0  2.5 -1.0
2027 -1.0  0.0 -2.8  5.3  1.0  1.5  0.0  2.5  0.0  0.0 -1.0  2.5
2028  2.5  0.0  0.7 -1.7  1.0 -2.0  0.0 -1.0  0.0  0.0 -1.0 -1.0
2029 -1.0  0.0  0.7 -1.7 -2.5  1.5  0.0 -1.0  0.0  0.0 -1.0 -1.0
2030 -1.0  0.0  0.7 -1.7  1.0 -2.0  0.0 -1.0  0.0  0.0 -1.0 -1.0
2031 -1.0  0.0  0.7 -1.7  1.0 -2.0  0.0 -1.0  0.0  0.0  2.5 -1.0
2032 -1.0  0.0 -2.8  5.3  1.0  1.5  0.0  2.5  0.0  0.0 -1.0  2.5
2033  2.5  0.0  0.7 -1.7  4.5 -2.0  0.0 -1.0  0.0  0.0 -1.0  2.5
2034  2.5  0.0  0.7 -1.7 -2.5  1.5  0.0 -1.0  0.0  0.0 -1.0 -1.0

1.3 Régresseur leap year

Le régresseur année bissextile (leap year), \(LY_t\) doit être créé à la main. Il est égal à la différence entre le nombre de jours dans le mois \(t\) et le nombre de jours moyens dans le mois \(t\), \(\bar N_t\). Tous les mois ont le même nombre de jours, sauf le mois de février qui est de 29 jours tous les 4 ans. \(\bar N_t\) est donc égal à 30 ou 31 si le mois considéré n’est pas un mois de février (et donc \(N_t - \bar N_t=0\)) à 28,25 en février1. \[ LY_{t} = \begin{cases} 0,75 & \mbox{si } t \mbox{ est un mois de février bissextil } \\ -0,25 & \mbox{si } t \mbox{ est un mois de février non bissextil } \\ 0 & \mbox{sinon} \end{cases} \]

Exercice

Créer une fonction leap_year qui permet de générer le régresseur leap year.

leap_year <- function(start = 1990, end = 2030, frequency = 12){
    ly <- ts(0, start = start, end = end, frequency = 12)
    mois_feb <- cycle(ly) == 2
    annees <- trunc(round(time(ly), 3)) # arrondi car parfois des pbs avec fonction time
    # On utilise la définition exacte
    is_ly <- (annees %% 400 == 0) |
        ((annees %% 4 == 0) & (annees %% 100 != 0))
    ly[mois_feb] <- 28 - 28.2425
    ly[mois_feb & is_ly] <- 29 - 28.2425
    # on change si besoin la fréquence
    stats::aggregate(ly, nfrequency = frequency) 
}
leap_year(frequency = 12)
         Jan     Feb     Mar     Apr     May     Jun     Jul     Aug     Sep
1990  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
1991  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
1992  0.0000  0.7575  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
1993  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
1994  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
1995  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
1996  0.0000  0.7575  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
1997  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
1998  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
1999  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2000  0.0000  0.7575  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2001  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2002  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2003  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2004  0.0000  0.7575  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2005  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2006  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2007  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2008  0.0000  0.7575  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2009  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2010  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2011  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2012  0.0000  0.7575  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2013  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2014  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2015  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2016  0.0000  0.7575  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2017  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2018  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2019  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2020  0.0000  0.7575  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2021  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2022  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2023  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2024  0.0000  0.7575  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2025  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2026  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2027  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2028  0.0000  0.7575  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2029  0.0000 -0.2425  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
2030  0.0000                                                                
         Oct     Nov     Dec
1990  0.0000  0.0000  0.0000
1991  0.0000  0.0000  0.0000
1992  0.0000  0.0000  0.0000
1993  0.0000  0.0000  0.0000
1994  0.0000  0.0000  0.0000
1995  0.0000  0.0000  0.0000
1996  0.0000  0.0000  0.0000
1997  0.0000  0.0000  0.0000
1998  0.0000  0.0000  0.0000
1999  0.0000  0.0000  0.0000
2000  0.0000  0.0000  0.0000
2001  0.0000  0.0000  0.0000
2002  0.0000  0.0000  0.0000
2003  0.0000  0.0000  0.0000
2004  0.0000  0.0000  0.0000
2005  0.0000  0.0000  0.0000
2006  0.0000  0.0000  0.0000
2007  0.0000  0.0000  0.0000
2008  0.0000  0.0000  0.0000
2009  0.0000  0.0000  0.0000
2010  0.0000  0.0000  0.0000
2011  0.0000  0.0000  0.0000
2012  0.0000  0.0000  0.0000
2013  0.0000  0.0000  0.0000
2014  0.0000  0.0000  0.0000
2015  0.0000  0.0000  0.0000
2016  0.0000  0.0000  0.0000
2017  0.0000  0.0000  0.0000
2018  0.0000  0.0000  0.0000
2019  0.0000  0.0000  0.0000
2020  0.0000  0.0000  0.0000
2021  0.0000  0.0000  0.0000
2022  0.0000  0.0000  0.0000
2023  0.0000  0.0000  0.0000
2024  0.0000  0.0000  0.0000
2025  0.0000  0.0000  0.0000
2026  0.0000  0.0000  0.0000
2027  0.0000  0.0000  0.0000
2028  0.0000  0.0000  0.0000
2029  0.0000  0.0000  0.0000
2030                        
# ou rjd3toolkit::lp_variable()

On peut également uiliser la fonction rjd3toolkit::ts_adjust() pour préajuster de l’effet année bissextile.

1.4 Exercice bilan

Exercice

Créer un objet regresseurs_JO qui contiendra tous les jeux de régresseurs plausibles. Par exemple :

  • le régresseur leap year

  • le jeu de régresseur trading days (TD7, lundi à samedi, dimanche = contraste)

  • le jeu de régresseur working days (TD2, lundi =… = vendredi, samedi=dimanche=contraste)

  • le jeu TD3 : lundi = … = vendredi, samedi et dimanche = contraste

frequency <- 12

gen_calendrier <- function(cal, frequency, start = c(1990, 1), end = c(2030, 1)) {
    length = (end[1] - start[1]) * frequency + end[2] - start[2]
    ly <- rjd3toolkit::lp_variable(frequency = frequency, start = start,
                                   length = length)
    # N'hésitez pas à ajouter les votres !
    TD7 <- calendar_td(cal, frequency = frequency, start = start, length = length,
                       groups = c(1, 2, 3, 4, 5, 6, 0))
    TD4 <- calendar_td(cal, frequency = frequency, start = start, length = length,
                       groups = c(1, 1, 1, 1, 2, 3, 0))
    TD3 <- calendar_td(cal, frequency = frequency, start = start, length = length,
                       groups = c(1, 1, 1, 1, 1, 2, 0))
    TD3c <- calendar_td(cal, frequency = frequency, start = start, length = length,
                        groups = c(1, 1, 1, 1, 2, 2, 0))
    TD2 <- calendar_td(cal, frequency = frequency, start = start, length = length,
                       groups = c(1, 1, 1, 1, 1, 0, 0))
    TD2c <- calendar_td(cal, frequency = frequency, start = start, length = length,
                        groups = c(1, 1, 1, 1, 1, 1, 0))
    
    reg_jo <- ts(cbind(TD2, TD2c, TD3, TD3c, TD4, TD7),
                 start = start, frequency = frequency)
    reg_jo <- ts.intersect(reg_jo,
                           ly)
    colnames(reg_jo) <- c(
        "TD2_semaine",
        "TD2c_lundi_samedi",
        sprintf("TD3_%s", c("semaine", "samedi")),
        sprintf("TD3c_%s", c("lundi_jeudi", "vendredi_samedi")),
        sprintf("TD4_%s", c("lundi_jeudi", "vendredi", "samedi")),
        sprintf("TD7_%s", c("lundi", "mardi", "mercredi", "jeudi", "vendredi", "samedi")),
        "leap_year")
    reg_jo
}

regresseurs_JO_mens <- gen_calendrier(CAL, frequency = 12)
regresseurs_JO_trim <- gen_calendrier(CAL, frequency = 4)
# Si l'on utilise des séries mensuelles :
regresseurs_JO <- regresseurs_JO_mens

Si l’on veut jongler entre R et JDemetra+, il faut ajouter les nouvelles variables dans le dictionnaire de variables de JDemetra+, voir le TP associé. Une solution est de créer un workspace vide depuis R et de l’utiliser pour charger vos données. Ci-dessous un code qui vous permet d’exporter ce workspace :

# On va créer un nouveau workspace où l'on va ajouter les nouvelles variables
library(RJDemetra)
library(rJava)
complete_variables <- function(liste_var, workspace){
    if(!is.mts(liste_var))
        stop("liste_var doit être de type mts")
    context_dictionary <- .jcall(workspace,"Lec/tstoolkit/algorithm/ProcessingContext;", "getContext")
    ts_variable_managers <- context_dictionary$getTsVariableManagers()
    ts_variables <- .jnew("ec/tstoolkit/timeseries/regression/TsVariables")
    jd_r_variables <- ts_variable_managers$get("r")
    if (is.null(jd_r_variables)) {
        ts_variable_managers$set("r",
                                 .jnew("ec/tstoolkit/timeseries/regression/TsVariables"))
        jd_r_variables <- ts_variable_managers$get("r")
    }
    jd_var_names <- jd_r_variables$getNames()
    
    model_var_names <- colnames(liste_var)
    
    for (i in seq_along(model_var_names)) {
        name <- model_var_names[i]
        dictionary_var <- jd_r_variables$get(name)
        tsvar <- .jnew("ec/tstoolkit/timeseries/regression/TsVariable",
                       name, RJDemetra:::ts_r2jd(liste_var[, i]))
        if (is.null(dictionary_var)) {
            jd_r_variables$set(name, tsvar)
        } else {
            warning(sprintf("La variable %s existe déjà", name))
        }
    }
}
# On renmme les variables pour éviter des conflits sur les noms
colnames(regresseurs_JO_trim) <- paste0(colnames(regresseurs_JO_trim), "_trim")
colnames(regresseurs_JO_mens) <- paste0(colnames(regresseurs_JO_mens), "_mens")
# Création d'un nouveaux
jws <- new_workspace()
# regresseurs_JO est l'objet mts qui contient tous vos régresseurs
# Il doit donc déjà être créé (voir code ci-dessus) !
complete_variables(regresseurs_JO_mens, jws)
complete_variables(regresseurs_JO_trim, jws)
save_workspace(jws,"wk_CJO.xml")

1.5 Effet graduel de Pâques

Prenons l’exemple de la vente de chocolats. Il est assez commun d’offrir des chocolats à Pâques : il y a donc une hausse des ventes autour du lundi de Pâques. Toutefois, ces ventes ne se font pas le jour de Pâques mais plusieurs jours avant, et plus on se rapproche du jour J, plus ces ventes sont importantes. C’est ce que l’on appel l’effet graduel de Pâques. Sous JDemetra+ on peut définir le nombre de jours avant Pâques pour lequel on considère qu’il y a un effet (easter_day.duration, entre 1 et 20) ou laisser ce choix à JDemetra+.

Exercice

Serait-il pertinent de considérer un effet graduel de Noël dans le modèle Reg-ARIMA ?

Non car l’effet graduel de Noël est en fait saisonnier car c’est un jour fixe ! Pour Pâques, comme c’est une fête mobile, les jours précédents peuvent être dans des mois différents en fonction de l’année considérée. Je ne suis pas entré dans les détails mais le régresseur utilisé pour la correction de l’effet graduel de Pâques est désaisonnalisé pour ne prendre en compte que l’effet voulu

Le régresseur associé à l’effet graduel de Pâques peut être généré en utilisant la fonction rjd3toolkit::easter_variable().

2 Utilisation des régresseurs dans RJDemetra

Dans RJDemetra, pour utiliser nos régresseurs jours ouvrables personnalisés, il faut créer sa propre spécification (fonctions x13_spec() ou regarima_spec_x13()) en utilisant l’option usrdef.varEnabled = TRUE, en spécifiant les régresseurs dans usrdef.var et indiquant que les régresseurs sont des régresseurs calendaires avec l’option usrdef.varType = "Calendar". Par exemple :

library(RJDemetra)
ipi_fr <- ipi_c_eu[, "FR"]
# On arrête la série en décembre 2019 pour éviter les changements de résultats
# liés aux futures actualisation des données de RJDemetra
ipi_fr <- window(ipi_fr, end = c(2019, 12))
# on garde le jeu reg6
wkd <- regresseurs_JO[,c(grep("TD7", colnames(regresseurs_JO), value = TRUE),
                         "leap_year")]
# Pour simplifier l'output, on enlève le "TD7_"
# mais ce n'est pas obligatoire
colnames(wkd) <- gsub("TD7_", "", colnames(wkd))
myspec1 <- regarima_spec_x13(spec = "RG5c",
                             usrdef.varEnabled = TRUE,
                             usrdef.var = wkd,
                             usrdef.varType = "Calendar",
                             easter.enabled = FALSE)
myreg1 <- regarima(ipi_fr, myspec1)
summary(myreg1)
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(7), no leap year effect, no Easter effect, outliers(2)

Coefficients:
ARIMA: 
          Estimate Std. Error  T-stat Pr(>|t|)    
Phi(1)     0.28608    0.13444   2.128   0.0340 *  
Phi(2)     0.24312    0.08163   2.978   0.0031 ** 
Theta(1)  -0.31214    0.13573  -2.300   0.0221 *  
BTheta(1) -0.63321    0.04438 -14.269   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Regression model: 
              Estimate Std. Error  T-stat Pr(>|t|)    
lundi         0.004817   0.001512   3.185 0.001578 ** 
mardi         0.009870   0.001730   5.705 2.50e-08 ***
mercredi      0.008431   0.001744   4.834 2.01e-06 ***
jeudi         0.006585   0.001777   3.705 0.000246 ***
vendredi      0.004267   0.001740   2.452 0.014696 *  
samedi       -0.015883   0.001588 -10.000  < 2e-16 ***
leap_year     0.021423   0.006296   3.403 0.000746 ***
LS (11-2008) -0.092121   0.015309  -6.017 4.50e-09 ***
LS (1-2009)  -0.067697   0.015325  -4.417 1.34e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.01943 on 333 degrees of freedom
Log likelihood = 871.8, aic =  1490, aicc =  1491, bic(corrected for length) = -7.663
Pour information

Pour ne pas repartir de zéro et avoir des modèles, il y a sous JDemetra+ 13 spécifications pré-définies décrites ici ou dans l’aide des fonction RJDemetra associés (?regarima, ?x13 ou ?tramoseats).

Pour faire des tests multiples sur les régresseurs jours ouvrables, on peut utiliser la fonction car::linearHypothesis(). Dans le modèle précédent, il parait clair que les régresseurs jours ouvrables sont significatifs. Toutefois, on peut se demander, si par parcimonie on peut simplifier le modèle en regroupant les jours de la semaine :

library(car)
linearHypothesis(myreg1,
                 c("lundi", "mardi", "mercredi", "jeudi", "vendredi", "samedi"),
                 c(0, 0, 0, 0, 0, 0), test = "F")

Linear hypothesis test:
lundi = 0
mardi = 0
mercredi = 0
jeudi = 0
vendredi = 0
samedi = 0

Model 1: restricted model
Model 2: myreg1

  Res.Df Df     F    Pr(>F)    
1    339                       
2    333  6 116.8 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
linearHypothesis(myreg1,
                 c("lundi = mardi","mardi = mercredi","mercredi = jeudi","jeudi = vendredi"), test = "F")

Linear hypothesis test:
lundi - mardi = 0
mardi - mercredi = 0
mercredi - jeudi = 0
jeudi - vendredi = 0

Model 1: restricted model
Model 2: myreg1

  Res.Df Df      F  Pr(>F)  
1    337                    
2    333  4 2.1578 0.07349 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Exercice

Essayez maintenant sur vos propres séries. Tester de laisser ou non l’effet graduel de Pâques (easter.enabled = TRUE).

3 Test de la présence de jours ouvrables résiduels

Un point important lorsque le fait de la correction de jours ouvrables est de tester s’il reste un effet jour ouvrable après la correction. La fonction rjd3toolkit::td_f() peut aider à le faire.

Généralement ce test est effectué après la décomposition, sur la composante désaisonnalisée ou sur l’irrégulier. Plutôt que la fonction regarima on va utiliser la fonction x13 qui effectue la décomposition sur la série linéarisée. Ces tests sont disponibles dans le sous-objet .$diagnostics ("f-test on sa (td)" et "f-test on i (td)") :

myspec1_sa <- x13_spec(spec = "RSA5c",
                       usrdef.varEnabled = TRUE,
                       usrdef.var = wkd,
                       usrdef.varType = "Calendar",
                       easter.enabled = FALSE)
mysa <- x13(ipi_fr, myspec1_sa)
# On retrouve d'ailleurs la partie regarima
# summary(mysa$regarima)
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.084
 Seasonal     60.997
 Irregular     0.780
 TD & Hol.     2.580
 Others       34.241
 Total       100.682

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

Residual seasonality tests
                                      P.value
 qs test on sa                          0.078
 qs test on i                           0.054
 f-test on sa (seasonal dummies)        0.998
 f-test on i (seasonal dummies)         0.988
 Residual seasonality (entire series)   1.000
 Residual seasonality (last 3 years)    0.901
 f-test on sa (td)                      0.068
 f-test on i (td)                       0.229

Sous JDemetra+, les tests affichés portent sur les 8 dernières années et dans RJDemetra sur la série entière ! Pour reproduire les résultats de JDemetra+, utiliser la fonction rjd3toolkit::td_f(). Pour le test, six spécifications différentes sont possibles :

  1. Par défaut sous JDemetra+ 2.x.y et model = "R100" sous rjd3toolkit \[ y_t=c + \alpha y_{t-1} + \sum_{i=1}^{6} \beta_i (N_{it} - N_{7t}) + \varepsilon_t \]

  2. model = "D1" \[ \Delta y_t - \overline{\Delta y} =\sum_{i=1}^{6} \beta_i \Delta(N_{it} - N_{7t}) + \varepsilon_t \]

  3. model = "DY" \[ \Delta_s y_t - \overline{\Delta_s y} =\sum_{i=1}^{6} \beta_i \Delta_s(N_{it} - N_{7t}) + \varepsilon_t \]

  4. model = "D1DY" \[ \Delta_s \Delta y_t - \overline{\Delta_s \Delta y} =\sum_{i=1}^{6} \beta_i \Delta_s\Delta(N_{it} - N_{7t}) + \varepsilon_t \]

  5. model = "AIRLINE" \[ y_t =\sum_{i=1}^{6} \beta_i (N_{it} - N_{7t}) + \varepsilon_t\text{ avec }\varepsilon_t\sim ARIMA(0,1,1)(0,1,1) \]

  6. model = "R011" (par défaut sous JDemetra+ 3.x.y) \[ y_t =\sum_{i=1}^{6} \beta_i (N_{it} - N_{7t}) + \varepsilon_t\text{ avec }\varepsilon_t\sim ARIMA(0,1,1) \]

  7. model = "WN" \[ y_t - \bar y =\sum_{i=1}^{6} \beta_i (N_{it} - N_{7t}) + \varepsilon_t \]

avec \(y_t\) pris en logarithme si le schéma est multiplicatif. Dans tous les cas \((H_0):\beta_1=\dots = \beta_6=0\) et les régresseurs utilisés ne prennent pas en compte le calendrier personnalisé que l’on a créé !

library(rjd3toolkit)
sa <- mysa$final$series[,"sa"]
i <- mysa$final$series[,"i"]
if (mysa$regarima$model$spec_rslt[,"Log transformation"]) {
    sa <- log(sa)
    i <- log(i)
}
rjd3toolkit::td_f(sa, nyears = 8, model = "R100")
Value: 6.203392 
P-Value: 0.0000 
rjd3toolkit::td_f(i, nyears = 8, model = "R100")
Value: 3.759636 
P-Value: 0.0022 
# Résultats différents sur l'ensemble de la série
rjd3toolkit::td_f(sa, nyears = 0, model = "R100")
Value: 3.128356 
P-Value: 0.0053 
rjd3toolkit::td_f(i, nyears = 0, model = "R100")
Value: 1.161422 
P-Value: 0.3265

En réalité les tests affichés dans JDemetra+ et RJDemetra ne sont pas effectuées sur les composantes “finales” de l’irrégulier et la série désaisonnalisée, mais sur les composantes issues de X-11 avant ajout des points atypiques du préajustement. Pour X-13-ARIMA ces séries ne sont pas exportées par défaut (et ne sont pas exportables depuis JDemetra+). Pour les extraires, nous allons réestimer le modèle2 :

jmodel <- suppressWarnings(jx13(get_ts(mysa), x13_spec(mysa)))

sa <- get_indicators(jmodel, "decomposition.sa_cmp")[[1]]
i <- get_indicators(jmodel, "decomposition.i_cmp")[[1]]
# Pour mettre tous les résultats sous forme de matrice :
t(simplify2array(
    list(
        rjd3toolkit::td_f(sa, nyears = 0, model = "R100"),
        rjd3toolkit::td_f(sa, nyears = 8, model = "R100"),
        rjd3toolkit::td_f(i, nyears = 0, model = "R100"),
        rjd3toolkit::td_f(i, nyears = 8, model = "R100"))
))
     value    pvalue      
[1,] 2.646733 0.01594339  
[2,] 6.207867 1.830266e-05
[3,] 1.166123 0.3238788   
[4,] 3.728232 0.002376407

En utilisant la fonction rjd3toolkit::sarima_estimate() et le package car, vous pouvez aussi construire vous-même le test3 :

car::linearHypothesis(
    rjd3toolkit::sarima_estimate(
        sa,
        order = c(0, 1, 1),
        seasonal = c(0, 0, 0),
        mean = FALSE,
        xreg = rjd3toolkit::td(s = sa)
    ),
    c("group_1 = 0", "group_2 = 0", "group_3 = 0", 
      "group_4 = 0", "group_5 = 0", "group_6 = 0"),
    test = "F"
)

Linear hypothesis test:
group_1 = 0
group_2 = 0
group_3 = 0
group_4 = 0
group_5 = 0
group_6 = 0

Model 1: restricted model
Model 2: rjd3toolkit::sarima_estimate(sa, order = c(0, 1, 1), seasonal = c(0, 
    0, 0), mean = FALSE, xreg = rjd3toolkit::td(s = sa))

  Res.Df Df     F Pr(>F)
1    357                
2    351  6 1.473 0.1864

Une autre solution est de passer par les objets jSA et d’exporter les indicateurs "diagnostics.td-i-last" et "diagnostics.td-sa-last", ou bien de rajouter ces indicateurs dans le paramètre userdefined de x13() :

jsa <- jx13(get_ts(mysa), x13_spec(mysa))
get_indicators(jsa,"diagnostics.td-sa-last")
$`diagnostics.td-sa-last`
[1] 5.1596190676 0.0001381414
attr(,"description")
[1] "F with 6 degrees of freedom in the nominator and 88 degrees of freedom in the denominator"
get_indicators(jsa,"diagnostics.td-i-last")
$`diagnostics.td-i-last`
[1] 3.306097404 0.005557358
attr(,"description")
[1] "F with 6 degrees of freedom in the nominator and 88 degrees of freedom in the denominator"
mysa <- x13(ipi_fr, myspec1_sa, userdefined = c("diagnostics.td-sa-last", "diagnostics.td-i-last"))
mysa$user_defined$`diagnostics.td-sa-last`
[1] 5.1596190676 0.0001381414
attr(,"description")
[1] "F with 6 degrees of freedom in the nominator and 88 degrees of freedom in the denominator"
mysa$user_defined$`diagnostics.td-i-last`
[1] 3.306097404 0.005557358
attr(,"description")
[1] "F with 6 degrees of freedom in the nominator and 88 degrees of freedom in the denominator"
# # ou :
# t(simplify2array(
#   mysa$user_defined
# ))

Privilégier plutôt ces deux dernières solutions, cela évite se regarder le schéma de décomposition (mais on ne peut pas personnaliser le test).

Notes de bas de page

  1. En réalité, la vraie valeur est 28,2425. En effet, une année bissextile est une année divisible par 4 mais pas par 100, sauf si elle est divisible par 400 : 1900 n’était pas une année bissextile mais 2000 l’était !↩︎

  2. Vous remarquerez que le résultat est différent de celui affiché dans JDemetra+ : cela vient de la méthode d’estimation du test utilisée !↩︎

  3. Vous pouvez également utiliser le code vu dans la section 2 pour estimer un modèle automatique.↩︎