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 <- CAL_FR <- 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_mobiles <- lapply(jours_mobiles, function(x){
#   na.omit(as.character(format(x, "%Y-%m-%d")))
#   })
# cat(unlist(lapply(names(jours_mobiles), function(day){
#   c(
#     sprintf("# %s", day),
#     sprintf('single_day("%s")', jours_mobiles[[day]])
#   )
# })),
#   sep= ",\n")

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é, ou créer le calendrier depuis JDemetra+. 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 deux groupes de variables en fonction de la fréquence
ctxt <- rjd3toolkit::modelling_context(
    # on appelle "CAL" le calendrier
    calendars = list(CAL = CAL), 
    # on crée un groupe de variables "cjo_mens" contenant les régresseurs mensuels
    # et un groupe de variables "cjo_trim" contenant les régresseurs trimestriels
    variables = list(cjo_mens = regresseurs_JO_mens,
                     cjo_trim = regresseurs_JO_trim)
) 

jws <- rjd3workspace::jws_new(ctxt)
# On peut également ajouter les calendriers et les variables avec les fonctions :
# rjd3workspace:::add_variables()
# rjd3workspace:::add_calendar()
# Pour modifier un workspace existant :
# rjd3workspace::set_context()
rjd3workspace::save_workspace(jws, "wk_CJO_v3.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 du calendrier personnalisé

Dans la V3 de RJDemetra, pour utiliser notre calendrier personnalisé, il faut :

  1. créer sa propre spécification (fonctions rjd3x13::x13_spec() ou rjd3x13::regarima_spec)

  2. ajouter les régresseurs dans le contexte grâce au paramètre context de rjd3x13::x13() et avec la fonction rjd3toolkit::modelling_context()

  3. modifier la spécification avec la fonction rjd3toolkit::set_tradingdays().

library(rjd3x13)
ipi_fr <- RJDemetra::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))
ctxt <- rjd3toolkit::modelling_context(
    # on crée un groupe de variables "cjo" contenant les régresseurs
    variables = list(cjo = regresseurs_JO)
) 

spec1_jd3 <- rjd3x13::spec_regarima("rg5c") |> 
    rjd3toolkit::set_tradingdays(
        option = "UserDefined",
        uservariable = paste0("cjo.", c(grep("TD7", colnames(regresseurs_JO), value = TRUE),
                                        "leap_year"))
    )
reg1_jd3 <- rjd3x13::regarima(ipi_fr, spec1_jd3, ctxt)
summary(reg1_jd3)
Method: RegARIMA 
Log-transformation: yes 
SARIMA model: (2,1,1) (0,1,1)

Coefficients
          Estimate Std. Error  T-stat Pr(>|t|)    
phi(1)     0.26571    0.12732   2.087  0.03765 *  
phi(2)     0.24252    0.07797   3.110  0.00203 ** 
theta(1)  -0.31346    0.12788  -2.451  0.01475 *  
btheta(1) -0.62456    0.05716 -10.926  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Regression model:
                  Estimate Std. Error  T-stat Pr(>|t|)    
cjo.TD7_lundi     0.003563   0.001535   2.322 0.020849 *  
cjo.TD7_mardi     0.009884   0.001696   5.827 1.33e-08 ***
cjo.TD7_mercredi  0.008337   0.001710   4.877 1.68e-06 ***
cjo.TD7_jeudi     0.007173   0.001752   4.093 5.35e-05 ***
cjo.TD7_vendredi  0.004172   0.001707   2.444 0.015026 *  
cjo.TD7_samedi   -0.016108   0.001557 -10.348  < 2e-16 ***
cjo.leap_year     0.021401   0.006163   3.472 0.000584 ***
easter           -0.013402   0.004103  -3.266 0.001204 ** 
LS (2008-11-01)  -0.091168   0.015194  -6.000 5.15e-09 ***
LS (2009-01-01)  -0.069050   0.015216  -4.538 7.95e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Number of observations: 360, Number of effective observations: 347, Number of parameters: 15
Loglikelihood: 877.1828, Adjusted loglikelihood: -725.5375
Standard error of the regression (ML estimate): 0.01913804 
AIC: 1481.075, AICc: 1482.525, BIC: 1538.815

Une autre façon de faire est d’ajouter le calendrier du personnalisé comme calendrier par défaut, ce qui permet de garder les options de JDemetra+ comme le préajustement de l’effet leap-year ou la sélection automatique des jeux de régresseurs.

ctxt <- rjd3toolkit::modelling_context(
    # on appelle "FR" le calendrier francais
  # Remplacer ici par votre calendrier
    calendars = list(FR = CAL_FR), 
    # on crée un groupe de variables "cjo" contenant les régresseurs
    variables = list(cjo = regresseurs_JO)
) 

spec2_jd3 <- rjd3x13::spec_regarima("rg5c") |> 
    rjd3toolkit::set_tradingdays(
        calendar.name = "FR"
    )
reg2_jd3 <- rjd3x13::regarima(ipi_fr, spec2_jd3, ctxt)
summary(reg2_jd3)
Method: RegARIMA 
Log-transformation: yes 
SARIMA model: (0,1,1) (0,1,1)

Coefficients
          Estimate Std. Error T-stat Pr(>|t|)    
theta(1)  -0.49888    0.04590 -10.87   <2e-16 ***
btheta(1) -0.63802    0.04539 -14.06   <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.002512   0.001290   1.947 0.052342 .  
tuesday          0.009517   0.001380   6.896 2.68e-11 ***
wednesday        0.005868   0.001404   4.181 3.71e-05 ***
thursday         0.012120   0.001464   8.276 3.05e-15 ***
friday           0.005087   0.001428   3.562 0.000421 ***
saturday        -0.014367   0.001290 -11.140  < 2e-16 ***
easter          -0.009844   0.003566  -2.760 0.006092 ** 
LS (2008-11-01) -0.080722   0.014104  -5.723 2.32e-08 ***
LS (2009-01-01) -0.069788   0.014117  -4.944 1.21e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Number of observations: 360, Number of effective observations: 347, Number of parameters: 12
Loglikelihood: 916.6591, Adjusted loglikelihood: -686.0612
Standard error of the regression (ML estimate): 0.0170764 
AIC: 1396.122, AICc: 1397.056, BIC: 1442.314

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(reg2_jd3,
                 c("monday", "tuesday", "wednesday", "thursday", "friday", "saturday"),
                 c(0, 0, 0, 0, 0, 0), test = "F")

Linear hypothesis test:
monday = 0
tuesday = 0
wednesday = 0
thursday = 0
friday = 0
saturday = 0

Model 1: restricted model
Model 2: reg2_jd3

  Res.Df Df      F    Pr(>F)    
1    341                        
2    335  6 193.04 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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 .$result$diagnostics ("td.ftest.sa" et "td.ftest.i") :

myspec1_sa <- rjd3x13::spec_x13("rsa5c") |> 
    rjd3toolkit::set_tradingdays(
        calendar.name = "FR"
    )
mysa <- rjd3x13::x13(ipi_fr, myspec1_sa, ctxt)
# On retrouve d'ailleurs la partie regarima
# summary(mysa$result$preprocessing)
mysa$result$diagnostics
$vardecomposition
$vardecomposition$cycle
[1] 0.08658561

$vardecomposition$seasonal
[1] 0.5463659

$vardecomposition$irregular
[1] 0.002511745

$vardecomposition$calendar
[1] 0.02553243

$vardecomposition$others
[1] 0.2807974

$vardecomposition$total
[1] 0.9417932


$seas.ftest.i
Value: 0.5949302 
P-Value: 0.8303 

$seas.ftest.sa
Value: 0.5909521 
P-Value: 0.8335 

$seas.qstest.i
Value: 0 
P-Value: 1.0000 

$seas.qstest.sa
Value: 0 
P-Value: 1.0000 

$td.ftest.i
Value: 3.240159 
P-Value: 0.0041 

$td.ftest.sa
Value: 3.035765 
P-Value: 0.0066 
# Pour afficher sous forme de tableau :
round(t(sapply(
    mysa$result$diagnostics[c("td.ftest.sa", "td.ftest.i")],
    unlist
)), 4)
             value pvalue
td.ftest.sa 3.0358 0.0066
td.ftest.i  3.2402 0.0041

Sous JDemetra+ 2.x.y, les tests affichés portent sur les 8 dernières années mais sur JDemetra+ 3.x.y il est calculé sur l’ensemble de la série. 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+ 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$result$decomposition$d11
i <- mysa$result$decomposition$d13
# On restreint la série car la composante contient les séries prolongées
sa <- window(sa, start = start(ipi_fr), end = end(ipi_fr))
i <- window(i, start = start(ipi_fr), end = end(ipi_fr))
if (mysa$result$preprocessing$description$log) {
    sa <- log(sa)
    i <- log(i)
}
rjd3toolkit::td_f(sa, nyears = 0, model = "R011")
Value: 3.035765 
P-Value: 0.0066 
rjd3toolkit::td_f(i, nyears = 0, model = "R011")
Value: 3.240159 
P-Value: 0.0041 
# Résultats sur les 8 dernières années
rjd3toolkit::td_f(sa, nyears = 0, model = "R011")
Value: 3.035765 
P-Value: 0.0066 
rjd3toolkit::td_f(i, nyears = 0, model = "R011")
Value: 3.240159 
P-Value: 0.0041 
# Pour mettre tous les résultats sous forme de matrice :
round(t(sapply(
  list(
    rjd3toolkit::td_f(sa, nyears = 0, model = "R011"),
    rjd3toolkit::td_f(sa, nyears = 8, model = "R011"),
    rjd3toolkit::td_f(i, nyears = 0, model = "R011"),
    rjd3toolkit::td_f(i, nyears = 8, model = "R011")
  ),
  unlist
)),4)
      value pvalue
[1,] 3.0358 0.0066
[2,] 2.3481 0.0377
[3,] 3.2402 0.0041
[4,] 2.3973 0.0342

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

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 3.0358 0.006596 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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 pouvez également utiliser le code vu dans la section 2 pour estimer un modèle automatique.↩︎