Devoir 3

Le jeu de données Metro Montreal.csv contient le kilométrage mensuel parcouru par l’ensemble des wagons de métro de la Ville de Montréal de janvier 2019 à septembre 2023, ainsi que le kilométrage calculé pour la préparation du budget d’exploitation du métro de 2019 à 2022. On aimerait développer un modèle de décomposition afin d’améliorer les prévisions qui serviront à établir les futurs budgets. Les données de janvier à septembre 2023 serviront de plage d’essai.

rm(list=ls())
setwd("~/USHERBROOKE/ESCUELA/ÉTÉ/MQG813 - Statistiques décisionnelles avancées/Jeux de données-20250506")
source("MQG813_librairie_temporel.R")

## Loading required package: DescTools

## Loading required package: vcd

## Loading required package: grid

## Loading required package: gmodels

## Registered S3 method overwritten by 'gdata':
##   method         from     
##   reorder.factor DescTools

## Loading required package: tidyverse

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: fastDummies
## 
## Loading required package: rio
## 
## Loading required package: car
## 
## Loading required package: carData
## 
## 
## Attaching package: 'car'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## 
## The following object is masked from 'package:purrr':
## 
##     some
## 
## 
## The following object is masked from 'package:DescTools':
## 
##     Recode
## 
## 
## Loading required package: Hmisc
## 
## 
## Attaching package: 'Hmisc'
## 
## 
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## 
## 
## The following objects are masked from 'package:DescTools':
## 
##     %nin%, Label, Mean, Quantile
## 
## 
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## 
## 
## Loading required package: modelr
## 
## Loading required package: pander
## 
## Loading required package: nlme
## 
## 
## Attaching package: 'nlme'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     collapse
## 
## 
## Loading required package: knitr
## 
## Loading required package: ggpubr
## 
## Loading required package: forecast
## 
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo 
## 
## 
## Attaching package: 'forecast'
## 
## 
## The following object is masked from 'package:ggpubr':
## 
##     gghistogram
## 
## 
## The following object is masked from 'package:nlme':
## 
##     getResponse
## 
## 
## The following object is masked from 'package:DescTools':
## 
##     BoxCox
## 
## 
## Loading required package: vars
## 
## Loading required package: MASS
## 
## 
## Attaching package: 'MASS'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## 
## Loading required package: strucchange
## 
## Loading required package: zoo
## 
## 
## Attaching package: 'zoo'
## 
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## 
## Loading required package: sandwich
## 
## 
## Attaching package: 'strucchange'
## 
## 
## The following object is masked from 'package:stringr':
## 
##     boundary
## 
## 
## Loading required package: urca
## 
## Loading required package: lmtest
## 
## 
## Attaching package: 'vars'
## 
## 
## The following object is masked from 'package:DescTools':
## 
##     Phi

#Verifie si une librairie est installee, l'installe au besoin et la charge
list.of.packages <- c("modelr", "rio", "tidyverse", "DescTools", "pastecs", 
                      "psych", "pander", "gmodels", "vcd", "fastDummies", 
                      "questionr", "PerformanceAnalytics", "Hmisc", "FSA", 
                      "car", "robustHD", "hexView", "ggplot2", "dplyr", 
                      "forecast", "latticeExtra", "numbers", "zoo", "lubridate",
                      "stringr", "fastDummies", "descr","tinytex","readr")

new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]

if(length(new.packages)) {install.packages(new.packages)}
lapply(list.of.packages, require, character.only = TRUE)

## Loading required package: pastecs
## 
## Attaching package: 'pastecs'
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
## 
## Loading required package: psych
## 
## Attaching package: 'psych'
## 
## The following object is masked from 'package:Hmisc':
## 
##     describe
## 
## The following object is masked from 'package:car':
## 
##     logit
## 
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
## 
## The following objects are masked from 'package:DescTools':
## 
##     AUC, ICC, SD
## 
## Loading required package: questionr
## 
## Attaching package: 'questionr'
## 
## The following object is masked from 'package:psych':
## 
##     describe
## 
## The following objects are masked from 'package:Hmisc':
## 
##     describe, wtd.mean, wtd.table, wtd.var
## 
## Loading required package: PerformanceAnalytics
## Loading required package: xts
## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################
## 
## Attaching package: 'xts'
## 
## The following objects are masked from 'package:pastecs':
## 
##     first, last
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## 
## Attaching package: 'PerformanceAnalytics'
## 
## The following object is masked from 'package:vcd':
## 
##     Kappa
## 
## The following object is masked from 'package:graphics':
## 
##     legend
## 
## Loading required package: FSA
## Registered S3 methods overwritten by 'FSA':
##   method       from
##   confint.boot car 
##   hist.boot    car 
## ## FSA v0.9.6. See citation('FSA') if used in publication.
## ## Run fishR() for related website and fishR('IFAR') for related book.
## 
## Attaching package: 'FSA'
## 
## The following object is masked from 'package:psych':
## 
##     headtail
## 
## The following object is masked from 'package:car':
## 
##     bootCase
## 
## Loading required package: robustHD
## Loading required package: perry
## Loading required package: parallel
## 
## Attaching package: 'perry'
## 
## The following object is masked from 'package:modelr':
## 
##     mape
## 
## Loading required package: robustbase
## Loading required package: hexView
## Loading required package: latticeExtra
## Loading required package: lattice
## 
## Attaching package: 'latticeExtra'
## 
## The following object is masked from 'package:ggplot2':
## 
##     layer
## 
## The following object is masked from 'package:vcd':
## 
##     rootogram
## 
## Loading required package: numbers
## 
## Attaching package: 'numbers'
## 
## The following object is masked from 'package:PerformanceAnalytics':
## 
##     Omega
## 
## The following object is masked from 'package:psych':
## 
##     omega
## 
## The following objects are masked from 'package:DescTools':
## 
##     GCD, LCM, Primes
## 
## Loading required package: descr
## 
## Attaching package: 'descr'
## 
## The following object is masked from 'package:questionr':
## 
##     freq
## 
## The following object is masked from 'package:gmodels':
## 
##     CrossTable
## 
## Loading required package: tinytex

## [[1]]
## [1] TRUE
## 
## [[2]]
## [1] TRUE
## 
## [[3]]
## [1] TRUE
## 
## [[4]]
## [1] TRUE
## 
## [[5]]
## [1] TRUE
## 
## [[6]]
## [1] TRUE
## 
## [[7]]
## [1] TRUE
## 
## [[8]]
## [1] TRUE
## 
## [[9]]
## [1] TRUE
## 
## [[10]]
## [1] TRUE
## 
## [[11]]
## [1] TRUE
## 
## [[12]]
## [1] TRUE
## 
## [[13]]
## [1] TRUE
## 
## [[14]]
## [1] TRUE
## 
## [[15]]
## [1] TRUE
## 
## [[16]]
## [1] TRUE
## 
## [[17]]
## [1] TRUE
## 
## [[18]]
## [1] TRUE
## 
## [[19]]
## [1] TRUE
## 
## [[20]]
## [1] TRUE
## 
## [[21]]
## [1] TRUE
## 
## [[22]]
## [1] TRUE
## 
## [[23]]
## [1] TRUE
## 
## [[24]]
## [1] TRUE
## 
## [[25]]
## [1] TRUE
## 
## [[26]]
## [1] TRUE
## 
## [[27]]
## [1] TRUE
## 
## [[28]]
## [1] TRUE
## 
## [[29]]
## [1] TRUE

install_formats()
include=F
echo=F

QUESTION 1 (4 points) Présentez et commentez le graphe séquentiel et les graphes saisonniers des kilométrages parcourus entre 2019 et septembre 2023.

commentaire: tendance generale a la baisse entre 2019 et 2023, avec une saisonnalité marquée, les mois d’été sont plus élevés que les mois d’hiver. On observe aussi un creux en avril 2020, probablement dû à un événement exceptionnel (comme la pandémie de COVID-19). On voit des cycles réguliers, ce qui laisse supposer une saisonnalité mensuelle

Le graphe séquentiel ci-dessous illustre une certaine tendance globale à la baisse, avec une valeur anormalement faible en avril 2020, probablement causée par le confinement imposé par le gouvernement en raison de la pandémie.

#Importation des donnees
metro <- import("Metro_Montreal.csv")

#Creation des variables t, annee et trimestre
metro <- metro %>% mutate(Date=as.Date(Date),
                        t = row_number(),
                        mois=factor(month(Date,label = TRUE, abbr = FALSE,)),
                        annee=year(Date))

#plage d’apprentissage (2019 à 2022)
app <- c(1:48)
#(janvier à septembre 2023)
essai<- c(49:57)
########################################
### Graphes sequentiel et saisonnier ###
########################################

#Graphe sequentiel
ggplot(data = metro, aes(x=Date, y = Km_reel)) +
  geom_point() +
  geom_line() +
  geom_smooth(method = "lm", se = FALSE)

## `geom_smooth()` using formula = 'y ~ x'

Graphique saisonnier

Commentaire: Le graphique saisonnier divide les données par mois. Les mois les plus élevés en moyenne sont mars,janvier et mai avec des valeurs autour de 7 500 000 km. Ce sont les mois où l’activité du métro semble la plus intense.

À l’inverse, les mois les plus faibles sont février et avril, avec des moyennes en dessous de 7 000 000 km.

On remarque une forte baisse en avril 2020, ce qui correspond au début de la pandémie de COVID-19. Cette chute affecte fortement la moyenne de ce mois.

Les graphes séquentiels ci-dessous illustrent un certain effet saisonnier. En effet, on remarque qu’en moyenne, le kilométrage parcouru est plus faible en février et en avril. Notons toutefois que la moyenne observée en avril est influencée par la valeur aberrante observée et pourrait être considérée moins représentative de ce qui se produit typiquement durant ce mois. Le mois de mars présente le kilométrage moyen le plus élevé, quoique pas si différente de celles qu’affichent plusieurs autres mois.

metro$moy <- NA
metro[app,] <- metro[app,]  %>% group_by(mois) %>% 
  mutate(moy = mean(Km_reel, na.rm=TRUE)) %>% ungroup()

ggplot(data = metro, aes(x = Date, y = Km_reel, group = mois)) +
  geom_point() +
  geom_line() +
  geom_hline(aes(yintercept = moy, group = mois), color = "red") +
  facet_wrap(~ mois, ncol =12)

## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_hline()`).

  labs(title = "Graphique saisonnier des kilométrages parcourus",
       x = "Mois", y = "Kilométrage réel") +
  theme_minimal()

## NULL

QUESTION 2 (8 points) On tentera d’abord de modéliser l’effet saisonnier à l’aide d’indices propres à chaque mois de l’année. (a) Présentez deux versions des indices saisonniers non ajustés à l’aide des moyennes et des moyennes tronquées des ratios calculés (revenus réels divisés par moyennes mobiles centrées). Suggéreriez-vous d’utiliser les moyennes tronquées ? Justifiez votre réponse.

commentaire: oui, en regardant la table des indices saisonniers par annee et mois on peut observer des valeurs qui sort de la moyenne du mois, comme avril avec le valeur 0.78 et en decembre avec 1.047, on peut dire que les moyennes tronquees sont plus pertinentes car elles ne tiennent pas compte des valeurs extremes qui peuvent fausser l’analyse.

Après avoir calculé les moyennes mobiles et les ratios individuels, voici ci-dessous un résumé des indices saisonniers calculés à partir de moyennes et à partir de moyennes tronquées. On peut remarquer certaines différences non négligeables, notamment pour le mois d’avril, qui affiche une différence de plus de 5 % entre les deux résultats. Cette différence est probablement due à la valeur aberrante, qui joue un grand rôle dans le calcul de la moyenne qui ne compte que trois valeurs. Pour éviter l’effet trop important de potentielles valeurs aberrantes, on suggèrera d’utiliser les moyennes tronquées. Notons toutefois que l’utilisation de moyennes non tronquées se défend aussi, étant donné le peu d’observations à notre disposition.

# Désaisonnalisation
metro <- desaisonnalisation(metro, "Km_reel", "Date", saison = "mois", app = app,
                           moyennes_tronquees = TRUE)
 
# Table des indices saisonniers par année et mois
table_indices <- metro[app, ] %>%
  group_by(annee, mois) %>%
  summarise(valeur = mean(Km_reel_RATIO, na.rm = TRUE)) %>%
  tidyr::pivot_wider(names_from = mois, values_from = valeur) %>%
  ungroup()

## `summarise()` has grouped output by 'annee'. You can override using the
## `.groups` argument.

# Résumé des indices
resume_indices <- metro[app, ] %>%
  group_by(mois) %>%
  summarise(
    moyenne = mean(Km_reel_RATIO, na.rm = TRUE),
    moyenne_tronquee = mean(Trim(Km_reel_RATIO, trim = 1, na.rm = TRUE), na.rm = TRUE),
    difference = moyenne - moyenne_tronquee
  ) %>%
  ungroup()
 resume_indices

## # A tibble: 12 × 4
##    mois      moyenne moyenne_tronquee difference
##    <ord>       <dbl>            <dbl>      <dbl>
##  1 January     1.02             1.01    0.00543 
##  2 February    0.955            0.938   0.0167  
##  3 March       1.05             1.04    0.00141 
##  4 April       0.914            0.969  -0.0549  
##  5 May         1.02             1.01    0.0114  
##  6 June        0.991            0.984   0.00636 
##  7 July        1.00             0.997   0.00420 
##  8 August      1.01             1.01    0.000701
##  9 September   0.988            0.995  -0.00696 
## 10 October     1.02             1.01    0.00948 
## 11 November    1.00             1.00   -0.00173 
## 12 December    1.03             1.03    0.00261

Remarque : Peu importe la réponse donnée en (a), on utilisera les moyennes tronquées (retrait de l’observation la plus élevée et de l’observation la plus faible) pour la suite.

Dans le cas des données sur lesquelles nous travaillons, à quelle statistique équivaut cette moyenne tronquée (avant ajustement) ? commentaire: La moyenne tronquée est équivalente à la moyenne des ratios,où l’on retire la plus haute et la plus basse valeur d’un mois/trimestre/jour donné . Cest nettoyée des extrêmes

Comme chaque mois est représenté à trois reprises sur la plage d’apprentissage, les moyennes tronquées correspondent donc aux médianes de chaque mois.

Présentez les indices saisonniers ajustés et interprétez l’indice le plus faible dans le contexte. C’est avril avec ~3.05% de moins. commentaire:L’indice saisonnier ajusté le plus faible est celui d’avril : 0,9695. Cela signifie que, toutes choses égales par ailleurs, le kilométrage moyen réel du mois d’avril est environ 3 % plus bas que la moyenne annuelle (1-0.9695).

Voici un résumé des indices saisonniers ajustés. L’indice le plus faible correspond à celui du mois de février (0,9386). On pourrait dire qu’en moyenne, le kilométrage effectué en février est 6,14 % inférieur au 1/12 du kilométrage annuel total.

# Table des indices ajustés
table_indices_ajustes <- metro[1:12, c("mois", "Km_reel_INDICES_SAIS")]
table_indices_ajustes

##         mois Km_reel_INDICES_SAIS
## 1    January            1.0116774
## 2   February            0.9386014
## 3      March            1.0450928
## 4      April            0.9695473
## 5        May            1.0080109
## 6       June            0.9844884
## 7       July            0.9969942
## 8     August            1.0110100
## 9  September            0.9951759
## 10   October            1.0128476
## 11  November            1.0051552
## 12  December            1.0255383

Présentez le graphe séquentiel des données réelles et désaisonnalisées. commentaire:La majorité des points oscille autour de 0, ce qui est cohérent avec un modèle ajusté. Cela montre que, dans l’ensemble, le modèle capte bien la tendance globale. Cependant, on observe une valeur aberrante très marquée en avril 2020, avec un résidu standardisé inférieur à –5 écarts-types. Cela indique un événement exceptionnel non capté par la variable temporelle t seule.

    #Graphe sequentiel
    ggplot(data = metro[app, ], aes(x = Date)) +
  geom_point(aes(y = Km_reel, color = "Km_reel")) +
  geom_line(aes(y = Km_reel, color = "Km_reel")) +
  geom_point(aes(y = Km_reel_DESAIS, color = "Km_reel_DESAIS")) +
  geom_line(aes(y = Km_reel_DESAIS, color = "Km_reel_DESAIS"), linetype = "dashed") +
  scale_colour_manual("",
                      breaks = c("Km_reel", "Km_reel_DESAIS"),
                      values = c("black", "red")) +
  labs(title = "Km reel vs Km désaisonnalisée ",
       y = "Kilométrage mensuel")

    #les rouge est plus aleatoire

QUESTION 3 (7 points) On tentera maintenant de modéliser l’effet de tendance par un modèle de régression linéaire. (a) Présentez un résumé du modèle de régression linéaire dont la seule variable explicative est la variable incrémentale t. Commentaire : le modele me dit que de base on a 7542964 km au temp 0, T est significativement negatif. 25.4% de la variation du kilométrage désaisonnalisé est expliquée par le modele.

modele1<- lm(Km_reel_DESAIS ~ t, data = metro[app,])
pander(summary(modele1))

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	7542964	79280	95.14	1.805e-54
t	-11147	2817	-3.957	0.0002603

Fitting linear model: Km_reel_DESAIS ~ t
Observations	Residual Std. Error	\(R^2\)	Adjusted \(R^2\)
48	270353	0.254	0.2378

(Intercept) = 7 542 964 : Cela veut dire que, au début de la période (quand t = 0, soit janvier 2019), le kilométrage mensuel désaisonnalisé était d’environ 7,54 millions de kilomètres.

t = –11 147 : Ce coefficient indique que, chaque mois, le kilométrage diminue en moyenne de 11 147 kilomètres. C’est donc une tendance à la baisse.

Cette diminution est statistiquement significative, car la valeur p est très petite (0,00026).

Enfin, le R² = 0,254, ce qui veut dire que le temps explique environ 25 % des variations du kilométrage désaisonnalisé. Le reste s’explique par d’autres facteurs non inclus dans ce modèle.

Présentez le graphe séquentiel des résidus standardisés du modèle de tendance. Comment pouvez-vous expliquer la valeur aberrante ?

Voici ci-desosus le graphe séquentiel des résidus standardisés du modèle estimant l’effet de tendance. La valeur aberrante semble être associée à l’observation du mois d’avril 2020 et pourrait probablement être expliquée par le confinement imposé par le gouvernement en raison de la pandémie.

#Ajout des previsisons et des residus
metro <- ajouter_previsions(metro, modele1,
                           intervalles=0.95, residus=T, 
                           noms_personnalises="TENDANCE")

#Graphe sequentiel des residus standardises
ggplot(data=metro[app,], aes(x=Date, y = ZRES_TENDANCE)) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 0)

(c) Ajoutez une variable explicative appropriée au modèle de tendance qui permettra de tenir compte de cet événement ponctuel. Présentez un résumé de ce nouveau modèle de régression linéaire et interprétez la valeur du R².

Commentaire : la variable BINAIRE a un effet significatif (p < 0.0001) et indique que, toutes choses égales par ailleurs, le kilométrage désaisonnalisé du mois d’avril 2020 est environ 1.43 million de km plus bas que la normale. La variable temporelle t reste significative. Le R² a ameliore beacoup (presque 40%).Cela montre que prendre en compte les événements exceptionnels dans un modèle permet d’améliorer considérablement sa précision et ne pas biaiser l’interprétation ni les prévisions futures.

metro$BINAIRE <- 0
metro[16, "BINAIRE"] <- 1

#Identification du modele
modele_t_binaire <- lm(Km_reel_DESAIS ~ t+ BINAIRE , 
                      data = metro[app,])

#Resume des statistiques du modele
pander(summary(modele_t_binaire))

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	7605252	51647	147.3	4.869e-62
t	-12470	1822	-6.845	1.733e-08
BINAIRE	-1433908	176719	-8.114	2.354e-10

Fitting linear model: Km_reel_DESAIS ~ t + BINAIRE (d) Établissez les prévisions et les intervalles de prévision des kilométrages désaisonnalisés. Présentez ensuite le graphe séquentiel des valeurs réelles et prévues des données désaisonnalisées, couvrant les années 2019 à 2022. commentaire:Le graphique montre les valeurs réelles du kilométrage désaisonnalisé en noir, les prévisions du modèle en rouge, et les intervalles de prévision à 95 % en pointillés bleus. Les bandes sont assez étroites, ce qui montre que le modèle fait des prévisions « fiables », sauf pour le point exceptionnel corrigé, qui sort de l’intervalle , ce qui est normal.
Observations	Residual Std. Error	\(R^2\)	Adjusted \(R^2\)
48	174167	0.6971	0.6837

#Ajout des previsisons et des residus
metro <- ajouter_previsions(metro, modele_t_binaire,
                           intervalles=0.95, residus=T, 
                           noms_personnalises="TENDANCE_BINAIRE")

#Graphe sequentiel avec previsions et intervalles de prevision
ggplot(data = metro[app,], aes(x = Date, y = Km_reel_DESAIS)) +
  geom_point() +
  geom_line() +
  geom_line(aes(y = PREV_TENDANCE_BINAIRE), color = "red") +
  geom_line(aes(y = IP_95_INF_TENDANCE_BINAIRE), color = "blue", linetype = "dashed")+
  geom_line(aes(y = IP_95_SUP_TENDANCE_BINAIRE), color = "blue", linetype = "dashed")

QUESTION 4 (4 points) Calculez la série des fluctuations restantes (CI). Si on se contentait des prévisions établies aux questions 2 et 3 pour calculer nos prévisions finales, on souhaiterait donc, à cette étape-ci, que ces fluctuations restantes soient aléatoires et d’effet neutre en moyenne (autour d’une certaine valeur). Vérifiez si chacune de ces deux conditions est plausible.

Bien que la série CI soit centrée autour de 1 (effet neutre), elle n’est pas aléatoire, ce qui montre que les prévisions des questions 2 et 3 ne suffisent pas à tout expliquer. Il faudrait donc ajouter d’autres variables ou tester un modèle plus complexe pour capter la dynamique restante.

À partir des statistiques descriptives ci-dessous, on remarque que la valeur moyenne de la série CI s’élève à 1, ce qui correspond à l’effet neutre attendu. L’écart-type de 2,32 % ainsi que les valeurs minimale (0,959) et maximale (1,06) suggèrent que les fluctuations restantes semblent relativement marginales.

Le graphe séquentiel de la série CI ci-dessous ne semble pas aléatoire, étant donné les longues suites d’observations qui se suivent. Le corrélogramme confirme cette observation. En particulier, l’autocorrélation d’ordre 1 (𝑟1=0,5483) est statistiquement significative au seuil 𝛼=0,05 .

#Création de la serie CI--- ON prend les donnes desaisonalisees que jutilise pour les prev
metro<- metro %>% mutate(CI = Km_reel_DESAIS/PREV_TENDANCE_BINAIRE)

#Graphe sequentiel
ggplot(data=metro[app,], aes(x=Date, y = CI)) +
  geom_point() +
  geom_line() +
  geom_hline(aes(yintercept = 1), color = "red")

# on va autour de 1 a cause de leffet multiplicatif


#Autocorrelations (variable CI)
correlogramme(metro$CI[app])

##    lag           ac          pac
## 1    1  0.548251273  0.548251273
## 2    2  0.326489244  0.037044645
## 3    3  0.465481477  0.390579559
## 4    4  0.365715829 -0.054406127
## 5    5  0.240997433  0.036689263
## 6    6  0.198229157 -0.091755749
## 7    7  0.021413398 -0.243088068
## 8    8 -0.007008421  0.011675467
## 9    9  0.012315635 -0.052420381
## 10  10 -0.137722912 -0.099822846
## 11  11 -0.264098503 -0.161562789
## 12  12 -0.247815834 -0.067250752
## 13  13 -0.278539265 -0.072571425
## 14  14 -0.226186881  0.137775838
## 15  15 -0.322484780 -0.212441149
## 16  16 -0.359414778  0.048539814
## 17  17 -0.310606074 -0.163225707
## 18  18 -0.256620507  0.057372824
## 19  19 -0.229771941  0.005453907
## 20  20 -0.230691728 -0.057015641
## 21  21 -0.246029486 -0.061504114
## 22  22 -0.141915085 -0.008013168
## 23  23 -0.125602594 -0.134979572
## 24  24 -0.169893066 -0.074581000
## 25  25 -0.118958359 -0.023178595

#on identifie autocorrelation forte, on essaye de corriger

#Identification du modele
modele_CI <- lm(CI ~ lag(CI, 1), 
                data = metro[app,])

#Resume des statistiques du modele
pander(summary(modele_CI))

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	0.4365	0.1244	3.51	0.001031
lag(CI, 1)	0.5642	0.1244	4.535	4.252e-05

Fitting linear model: CI ~ lag(CI, 1)
Observations	Residual Std. Error	\(R^2\)	Adjusted \(R^2\)
47	0.0195	0.3137	0.2984

#Ajout des previsisons et des residus
#metro <- ajouter_previsions(metro, modele_CI, plage_connue=app,
#                           intervalles=0.95, residus=T, 
 #                          noms_personnalises="CI")

#Graphe sequentiel des residus standardises
#ggplot(data=metro, aes(x=Date, y = ZRES_CI)) +
 # geom_point() +
  #geom_line() +
  #geom_hline(yintercept = 0)

#Autocorrelations des residus
#correlogramme(metro$ZRES_CI[app])
#Il reste encore dautocorrelation maintenant  dans le 3

#Identification du modele
#modele_CI_2 <- lm(CI ~ lag(CI, 1) +lag(CI,3),
#                data = metro[app,])




#Ajout des previsisons et des residus
#metro <- ajouter_previsions(metro, modele_CI_2, plage_connue=app,
 #                          intervalles=0.95, residus=T, 
  #                         noms_personnalises="CI_2")

#Graphe sequentiel des residus standardises
#ggplot(data=metro, aes(x=Date, y = ZRES_CI_2)) +
#  geom_point() +
 # geom_line() +
#  geom_hline(yintercept = 0)

#Autocorrelations des residus
#correlogramme(metro$ZRES_CI_2[app])


#Histogramme des residus
#ggplot(metro[app,], aes(x=ZRES_CI_2)) + geom_histogram()

#Test de normalite des residus (Shapiro-Wilk)
#pander(shapiro.test(metro$ZRES_CI_2[app]))

#Graphe sequentiel avec previsions et intervalles de prevision (plage d'apprentissage)---> pas demande dans cette question
#ggplot(data = metro[app,], aes(x = Date, y = CI)) +
 # geom_point() +
  #geom_line() +
  #geom_line(aes(y = PREV_CI_2), color = "red") +
  #geom_line(aes(y = IP_95_INF_CI_2), color = "blue", linetype = "dashed")+
  #geom_line(aes(y = IP_95_SUP_CI_2), color = "blue", linetype = "dashed")

QUESTION 5 (7 points) Afin de conserver un mod`ele simple, on d´ecide de ne pas mod´eliser les fluctuations restantes. On d´esire maintenant ´etablir les pr´evisions finales en multipliant les estimations des effets saisonniers et les pr´evisions de l’effet de tendance.

ici prev restant =ci, on va pas les ajouter

#Calcul des previsions
metro <- metro %>% mutate(KM_PREV = Km_reel_INDICES_SAIS*PREV_TENDANCE_BINAIRE)
# on utilise les prev de ci

Présentez le graphe séquentiel des kilométrages réels, prévus ainsi que budgétés pour les années 2019 à 2022.

#graphique sequentiel des km reel, prevus et budgetes
ggplot(data = metro[app,], aes(x = Date)) +
  geom_point(aes(y = Km_reel, color = "Km_reel")) +
  geom_line(aes(y = Km_reel, color = "Km_reel")) +
  geom_point(aes(y = KM_PREV, color = "KM_PREV")) +
  geom_line(aes(y = KM_PREV, color = "KM_PREV"), linetype = "dashed") +
  geom_point(aes(y = Km_budget, color = "Km_budget")) +
  geom_line(aes(y = Km_budget, color = "Km_budget"), linetype = "dotted") +
  scale_colour_manual("",
                      breaks = c("Km_reel", "KM_PREV", "Km_budget"),
                      values = c("black", "red", "blue")) +
  labs(title = "Kilométrages réels, prévus et budgétés",
       y = "Kilométrage mensuel")

Le graphique compare :les kilométrages réels (en noir), les kilométrages prévus par le modèle (en rouge) et les kilométrages budgétés (en bleu). Entre 2020 et 2022, les prévisions surestiment légèrement les valeurs réelles, mais restent inférieures aux valeurs budgétées, surtout après la “pandémie”.Mais u debut 2019 et fin 2022, les prévisions sont très proches des réels, ce qui montre une bonne qualité d’estimation en début et fin de période. Globalment les prévisions du modèle sont plus proches des valeurs réelles que les valeurs budgétées, ce qui montre une meilleure performance du modèle par rapport à la planification initiale.

Quel mois présente le taux d’erreur absolu moyen le plus faible (donc le meilleur, le plus bas) ? Justifiez votre réponse à l’aide des graphes séquentiels saisonniers des taux d’erreurs absolus sur la plage d’apprentissage.

Le mois qui présente le taux d’erreur absolu moyen le plus faible est juillet. Sur le graphique, la ligne rouge moyenne pour juillet est la plus basse parmi tous les mois, ce qui indique que les erreurs de prévision sont très faibles durant cette période.

#Creation des variables residu, residu_abs et residu_taux.. on avait pas les residu finaux
metro <- metro %>% mutate(residu = Km_reel-KM_PREV,
                        residu_abs = abs(residu),
                        residu_taux = residu_abs/Km_reel)

 #Creation des variables correspondant aux moyennes mensuelles 
#des residus absolus sur la plage d'app
metro$moy_res_abs <- NA
metro[app,] <- metro[app,] %>% group_by(mois) %>% mutate(moy_res_abs = mean(residu_abs, na.rm=TRUE)) %>% ungroup()

#Graphes saisonniers des residus absolus
ggplot(data=metro[app,], aes(x=Date, y = residu_abs)) +
  geom_point() +
  geom_line() +
  geom_hline(aes(yintercept = moy_res_abs, group = mois), color = "red") +
  facet_wrap(~ mois, ncol = 12)

(c) Présentez les indices de performance du modèle de décomposition sur la plage d’apprentissage et sur la plage d’essai. Comparez les résultats de la plage d’apprentissage aux indices de performance des kilométrages budgétés.

Sur la plage d’apprentissage, le modèle de décomposition performe mieux que les données budgétées. En effet, l’erreur absolue moyenne s’élève à 132 489 km pour le modèle de décomposition, alors qu’il s’élève à 228 654 km pour les kilométrages budgétés. Il en est de même pour le MAPE (1,80 % vs 3,29 %)

#Indices de performance
indices_performance("Km_reel", "KM_PREV", metro,app=app)

##                N      MAD     MAPE
## Apprentissage 48 132488.9 0.018034
## Essai          9 197571.0 0.028130

indices_performance("Km_reel", "Km_budget", metro[app,])

##                     N      MAD     MAPE
## Valeur de l'indice 48 228653.6 0.032898

Malgré les différences importantes présentées en (b), pourquoi les résultats présentés ne démontrent-ils pas que le modèle de décomposition développé dans ce devoir performerait mieux que le modèle actuellement utilisé par la Société de transport de Montréal pour établir un futur budget annuel ?

Même si le modèle donne de bons résultats entre 2019 et 2023, on ne peut pas dire avec certitude qu’il ferait mieux que celui de la STM. Il a été testé sur une période courte, et on ne connaît pas tous les éléments pris en compte dans la planification réelle. Il n’y a pas eu de comparaison complète avec le modèle STM, et on ne dispose pas des détails sur le modèle actuellement utilisé. Il faudrait donc plus de tests sur plusieurs années pour vraiment comparer

Les résultats présentés en (c) ne démontrent pas que le modèle de décomposition performe mieux, car ledit modèle a été créé à partir des données mêmes sur lesquelles on vient d’établir des prévisions. Lors de la création du budget, les administrateurs n’avaient évidemment pas accès aux données réelles. Il faudrait donc comparer les résultats sur la plage d’essai si on souhaite vérifier si le modèle de décomposition se généralise bien sur des données qui lui sont inconnues..

Devoir 3

Gaby

2025-05-28