INTRODUCTION

Une série temporelle, {yt, t = 1, 2, ··· , T}, est une suite d’observations d’un phénomène, indicées par une date, un temps. Le moment auquel l’observation est faite est généralement une information importante sur le phénomène observé.

Dans le cadre de notre étude statistique,nous serons amené à modéliser le trafic voyageur de la SNCF en 2e classe de 1963 à 1969. Plus précisement, il est question de réaliser une analyse descriptive de cette série temporelle dans un premier temps puis de prédire par la technique statistique de Holt Winters le nombre de voyageurs pour les 12 mois de l’année 1970.

Codes, packages et fonctions R

library(fpp)
library(moments)

PARTIE I- ANALYSE DESCRIPTIVE DE SERIE TEMPORELLE

1- Construction et présentation de la série temporelle

Créons d’abord le vecteur de la serie de données puis faisons sa transformation en serie temporelle.

A- Vecteur de la série de données

a<-c(1750,1560,1820,2090,1910,2410,3140,2850,2090,1850,1630,2420,1710,1600,1800,2120,2100,2460,3200,2960,2190,1870,1770,2270,1670,1640,1770,2190,2020,2610,3190,2860,2140,1870,1760,2360,1810,1640,1860,1990,2110,2500,3030,2900,2160,1940,2750,2330,1850,1590,1880,2210,2110,2480,2880,2670,2100,1920,1670,2520,1834,1792,1860,2138,2115,2485,2581,2639,2038,1936,1784,2391,1798,1850,1981,2085,2120,2491,2834,2725,1932,2058,1856,2553)

B- Tranformation des données en série temporelle (timeseries)

Nous avons des données en série temporelle de 1963 à 1969. La visualisation graphique des données nous fait présumer une série temporelle additive.

serie_t <- ts(a, frequency = 12, start = c(1963, 1))
serie_t  
##       Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
## 1963 1750 1560 1820 2090 1910 2410 3140 2850 2090 1850 1630 2420
## 1964 1710 1600 1800 2120 2100 2460 3200 2960 2190 1870 1770 2270
## 1965 1670 1640 1770 2190 2020 2610 3190 2860 2140 1870 1760 2360
## 1966 1810 1640 1860 1990 2110 2500 3030 2900 2160 1940 2750 2330
## 1967 1850 1590 1880 2210 2110 2480 2880 2670 2100 1920 1670 2520
## 1968 1834 1792 1860 2138 2115 2485 2581 2639 2038 1936 1784 2391
## 1969 1798 1850 1981 2085 2120 2491 2834 2725 1932 2058 1856 2553

C- Ajustement des données en série temporelle

Nous avons choisi d’ajuster les données en prenant en compte le nombre effectif de jours de chaque mois. Le nombre de passagers pouvant être influencé par ce facteur du nombre de jour effectif de chaque mois. Pour cela nous allons calculer une moyenne du nombre de passagers par jour pour chaque mois de la série et utiliser ces données pour nos analyses. C’est une autre forme de transformation des données pour que notre prédiction colle au mieux avec la réalité.

# Moyenne du nombre de passagers par jour pour chaque mois de la série
serie.new <- serie_t / monthdays(serie_t)

2- Visualisation graphique de la série temporelle

La representation sous forme graphique nous permettra d’analyser la distribution des données de la séerie temporelle.

par(mfrow=c(2,1), mar=c(3,3,3,3))
plot(serie_t, xlab="Années",ylab="Passagers",las=1, main=" Graphe de la série brute") 
plot(serie.new, xlab="Années",ylab="Passagers",las=1, main=" Graphe après ajustement de la série", col="blue")

La visualisation après ajustement des données semble légèrement différente de la première visualisation. Elle a pris en compte l’effectivité des jours sur chaque mois.

par(mfrow=c(2,1), mar=c(3,3,3,3))
hist(serie.new,xlab="Passagers",ylab="Fréquence",las=1, main="Graphe en Histogramme")
boxplot(serie.new,ylab="Passagers",las=1, main="Graphe en Boîte à moustâche")

L’histogramme nous montre une plus forte concentration des données entre 55 et 60 passagers journaliers en moyenne sur toute la distribution. La valeur médiane quant à elle est de l’ordre de 68 selon le boxplot.

3- Résumé numérique et indicateurs d’analyse descriptive

Il est question ici des indicateurs de centralité, de dispersion, de forme et de dépendance.

A- Indicateurs de centralité et de dispersion

# Résumé statistique de la serie temporelle : indices de centralité et de dispersion
summary(serie.new)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   53.87   59.92   68.15   71.26   82.09  103.23
mean(serie.new) 
## [1] 71.25629
var(serie.new)
## [1] 175.1482
sd(serie.new)
## [1] 13.23435
cvar <- (sd(serie.new)/mean(serie.new))*100
cvar
## [1] 18.57289

Les indices de centralité présentent un nombre de passagers moyen par jour egal à 71. Le plus faible enregistrement de passagers est de 54 pendant que la valeur la plus élevée du nombre de passagers enregistrée est de 103. La variance est de …. ce qui signifie que la dispersion du nombre de passagers par jour est très importante.

La moyenne journalière des passagers étant de 71, l’écart-type de 13 indique que les le nombre de passagers est légèrement dispersé autour de la moyenne. Il varie dans un intervalle fermé entre [71–13; 71+13] soit entre [58;84] passagers. Cette légère dispersion se confirmera effectivement avec le coefficient de variation qui est de 18%.

B- Indicateurs de forme

library(moments)
skewness(serie.new)
## [1] 0.7156052
kurtosis(serie.new) 
## [1] 2.473578

Notre distribution est étalée à droite et donc la majorité des effectifs soit 71% sont concentrées à droite de la moyenne. Le kurtosis étant inférieur à 3, la distribution est platikurtique.

C- Indicateurs de dépendance

# Autocorellation
result_acf=acf(serie.new, lag.max=12,plot = TRUE, ylab="Significativité",las=1, main="Graphe ACF")

print(data.frame(result_acf$lag,result_acf$acf)[1:12,])
##    result_acf.lag result_acf.acf
## 1      0.00000000     1.00000000
## 2      0.08333333     0.48720214
## 3      0.16666667     0.06521557
## 4      0.25000000    -0.21762204
## 5      0.33333333    -0.18896809
## 6      0.41666667    -0.35809519
## 7      0.50000000    -0.42671557
## 8      0.58333333    -0.36313868
## 9      0.66666667    -0.13723798
## 10     0.75000000    -0.16090896
## 11     0.83333333     0.07745376
## 12     0.91666667     0.42957674

Nous pouvons voir à partir de l’échantillon du corrélogramme que l’autocorrélation est particulièrement forte pour les décalages 8 et 9. Il existe bien des correlations non nulles. Si le nombre de mois augmente de 12 jours, le nombre de passagers augmente significativement de 90%. Par contre si le nombre de mois diminue de 6 mois, le nombre de passagers diminue également de 40%. La pointe au décalage 9 indique une forte corrélation entre chaque valeur et la valeur apparaissant 9 mois auparavant. Cette mesure de l’association entre des valeurs de séries actuelles et passées indique que les valeurs au lag 8 et 9 sont les valeurs de séries passées les plus utiles à la prévision de valeurs futures si nous restons dans un décallage de 12 mois. La corrélation positive indique que des valeurs actuelles élevées correspondent à des valeurs élevées au niveau du décalage spécifié. la série est très corrélée avec elle-même pour un décalage de 3. Ainsi, sur la base du graphique, nous voyons qu’il y a présence d’autocorrélation dans cette série temporelle dans les périodes où la ligne dépasse de la bande discontinue

# Autocorellation
result_pacf=pacf(serie.new, max=12,plot = TRUE, ylab="Significativité",las=1, main="Graphe PACF")

print(data.frame(result_pacf$lag,result_pacf$acf)[1:12,])
##    result_pacf.lag result_pacf.acf
## 1       0.08333333      0.48720214
## 2       0.16666667     -0.22573126
## 3       0.25000000     -0.20253640
## 4       0.33333333      0.05266116
## 5       0.41666667     -0.41034950
## 6       0.50000000     -0.22969109
## 7       0.58333333     -0.15029421
## 8       0.66666667     -0.18255683
## 9       0.75000000     -0.50285338
## 10      0.83333333      0.02391978
## 11      0.91666667      0.24529611
## 12      1.00000000      0.34099236

Nous pouvons voir à partir de l’échantillon du corrélogramme que l’autocorrélation aux lag12 dépassent les bornes de significativité positive et présente une forte corellation. Il existe bien des correlations non nulles. Si le nombre de jours augmente de 12 jours, le nombre de passagers augmente significativement de 90%. Par contre si le nombre de jours diminue de 6 jours, le nombre de passagers diminue également de 40%. La corrélation négative indique que des valeurs actuelles élevées correspondent à des valeurs faibles au niveau du décalage spécifié

PARTIE II- PREDICTIONS DU NOMBRE DE PASSAGERS POUR LES 12 PROCHAINS MOIS

1- Choix de la technique Holt-Winters pour la prédiction

Nous allons utiliser la technique de Holt-Winters. C’est un modele plus évolué qui va au-delà du modele linéaire en prenant en compte la saisonnalité pour la prévision.

Cette technique permet de mettre à jour les prédictions en t+1 sur la base de moyennes pondérées des valeurs passées. Ce lissage exponentiel triple, décompose la variable d’intérêt en deux composants, une tendance (T) et une saisonnalité (S), auxquels on peut ajouter un terme d’erreur (E).

Il existe deux grandes versions de ce modèle : une version additive qui est indiquée lorsque la variance est stable dans le temps, et une version multiplicative plus adéquate lorsque la variance croît/décroît dans le temps. Le lissage expo triple Hotwinter comporte les éléments suivants : (parametre gamma, parametre de la saisonalité c’est-à-dire les moments de pic et de flop et le phénonmene tendanciel).

La fonction Holt-Winters de notre série nous donne les meilleurs parametres alpha, beta et 12 saisonalités. Nous en déduisons que la série est saisonniere.

HoltWinters(serie.new)
## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = serie.new)
## 
## Smoothing parameters:
##  alpha: 0.02817997
##  beta : 0.1704337
##  gamma: 0.3411165
## 
## Coefficients:
##             [,1]
## a    71.52730872
## b    -0.05335708
## s1  -14.80211478
## s2  -11.65712333
## s3  -11.88794986
## s4   -2.30736016
## s5   -4.91712728
## s6   10.19829075
## s7   19.87702026
## s8   16.03135802
## s9   -4.32846114
## s10  -9.05555477
## s11  -9.94123679
## s12   6.83679471

Si les données prédites sont proches de la réalité alors nous pourrons utiliser cela pour faire nos prédictions. Mais pour cela il nous a fallu trouver un système de pondération qui me permet de coller aux valeurs réelles. Plus α , 𝜷, 𝜸 sont petits, plus on tient compte du passé lointain pour la prévision. Plus α , 𝜷, 𝜸 sont grands, plus on tient compte du passé récent pour la prévision.

Utilisons les meilleures valeurs de alpha, beta et gamma produite par l’algorithme Holt-Winters pour prédire les valeurs moyennes des passagers sur les 12 prochains mois avec un intervalle de confiance compris entre 90 et 95%

par(mfrow=c(2,1), mar=c(3,3,3,3))
prediction<- hw(serie.new, h=12, seasonal="add", level=c(90, 95))
prediction
##          Point Forecast    Lo 90     Hi 90    Lo 95     Hi 95
## Jan 1970       56.94778 48.41257  65.48298 46.77745  67.11810
## Feb 1970       59.53035 50.99514  68.06556 49.36002  69.70068
## Mar 1970       60.06542 51.53021  68.60063 49.89509  70.23575
## Apr 1970       70.85835 62.32314  79.39356 60.68802  81.02868
## May 1970       67.50982 58.97461  76.04503 57.33949  77.68015
## Jun 1970       83.41852 74.88331  91.95374 73.24819  93.58886
## Jul 1970       97.28449 88.74928 105.81971 87.11415 107.45483
## Aug 1970       90.72551 82.19029  99.26073 80.55516 100.89585
## Sep 1970       70.79431 62.25908  79.32954 60.62396  80.96466
## Oct 1970       61.32047 52.78524  69.85571 51.15011  71.49083
## Nov 1970       63.22231 54.68707  71.75755 53.05194  73.39268
## Dec 1970       77.04606 68.51081  85.58131 66.87569  87.21644
plot(prediction, xlab="Années",ylab="Passagers",las=1, main=" Graphe de la prédiction moyenne de passagers ")
plot(prediction)
readline(prompt = "Press [enter] to continue") 
## Press [enter] to continue
## [1] ""
lines(fitted(prediction), col="red", lty="longdash")

La courbe de couleur rouge represente les valeurs prédites pour la même période de données.

Ci-dessous les valeurs prédites en terme de moyenne pour la même période couverte par nos données d’origine.

# Moyenne des valeurs prédites pour la meme période du jeu de données
prediction$fitted
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 1963 56.40453 58.98922 59.52201 70.31402 66.96588 82.87079 96.73665 90.17581
## 1964 56.46960 59.05101 59.58289 70.37301 67.02170 82.92841 96.79658 90.23453
## 1965 56.52944 59.10936 59.64334 70.43397 67.08383 82.98922 96.86223 90.30015
## 1966 56.59179 59.17454 59.70845 70.50286 67.14936 83.05846 96.93138 90.36591
## 1967 56.67454 59.26017 59.79712 70.59330 67.24789 83.16075 97.03760 90.47504
## 1968 56.77751 59.36017 59.89964 70.69618 67.34866 83.26052 97.13549 90.56852
## 1969 56.86418 59.44455 59.98282 70.77952 67.43025 83.34031 97.20885 90.64699
##           Sep      Oct      Nov      Dec
## 1963 70.24599 60.76658 62.66949 76.48666
## 1964 70.30450 60.82584 62.72619 76.54902
## 1965 70.36838 60.88802 62.78715 76.61000
## 1966 70.43396 60.95307 62.85159 76.68888
## 1967 70.54296 61.06397 62.97508 76.78897
## 1968 70.63598 61.15509 63.06142 76.87952
## 1969 70.71606 61.23505 63.14044 76.96072

2- Prédiction du nombre de passagers de Janvier à Décembre de l’année suivante

Nous allons prédire le nombre de voyageurs pour les 12 prochains mois. Mais d’abord évaluons le modele à partir de la prédiction de la moyenne des passagers que nous venons d’effectuer. Nous allons nous baser sur l’analyse des erreurs de prédiction.

A- Evaluation et validation du modele par l’analyse des erreurs de prédiction

Les résidus ou erreurs observées sont définis comme étant les différences entre les valeurs observées et les valeurs estimées par un modèle de régression. L’analyse des résidus a pour objectif de tester la validité de la prédiction.

Il s’agit de vérifier si les résidus du modèle choisi suivent un processus bruit blanc Dans le cas contraire, le modèle n’est pas bon. Il va falloir choisir un autre ou reprendre les étapes précédentes. La satisfaction de la seule condition de bruit blanc suffit pour valider le modèle. Mais si les résidus suivent en plus un processus gaussien (loi normale), le modèle serait encore plus intéressant.

Ces visualisations vont nous permettre de détecter la présence d’une défaillance dans le modèle construit pour la prédiction du nombre de passagers. La valider de notre modèle de prédiction est nécessaire pour s’assurer que le modèle est effectivement capable de prédire avec précision et exactitude les valeurs d’une variable d’intérêt.

# récupérer les résidus 
residus=residuals(prediction)
residus
##               Jan          Feb          Mar          Apr          May
## 1963   0.04708598  -3.27493808  -0.81233399  -0.64734971  -5.35298163
## 1964  -1.30830807  -3.87859888  -1.51837638   0.29365805   0.72023726
## 1965  -2.65847398  -0.53793289  -2.54656760   2.56603124  -1.92254090
## 1966   1.79531093  -0.60311075   0.29155500  -4.16952550   0.91516002
## 1967   3.00287473  -2.47445642   0.84803845   3.07336999   0.81662373
## 1968   2.38378426   2.43293491   0.10036398   0.57048889   0.87714647
## 1969   1.13581689   6.62687792   3.92040346  -1.27951696   0.95684977
##               Jun          Jul          Aug          Sep          Oct
## 1963  -2.53745581   4.55367548   1.75967555  -0.57932589  -1.08915823
## 1964  -0.92841486   6.42922563   5.24934114   2.69550133  -0.50325501
## 1965   4.01077588   6.04099778   1.95791529   0.96494854  -0.56544269
## 1966   0.27487526   0.81055563   3.18247269   1.56604415   1.62757243
## 1967  -0.49407862  -4.13436947  -4.34601194  -0.54296023   0.87151841
## 1968  -0.42718332 -13.87742210  -5.43949101  -2.70264831   1.29652737
## 1969  -0.30697277  -5.78949485  -2.74375990  -6.31605553   5.15204880
##               Nov          Dec
## 1963  -8.33615591   1.57785865
## 1964  -3.72619191  -3.32321137
## 1965  -4.12048633  -0.48097236
## 1966  28.81507591  -1.52758968
## 1967  -7.30841552   4.50135551
## 1968  -3.59475669   0.24951123
## 1969  -1.27377161   5.39411874
#Graphe des résidus : Vérification de la blancheur des données résiduelles 
plot(residus)

# Autocorrelation : Vérifier s'il existe des autocorrélations non nulles aux décalages 1-12, pour les erreurs de prévision dans les échantillons.
acf(residus, lag.max=12, na.action = na.pass)

Box.test(residus, lag=12, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  residus
## X-squared = 11.751, df = 12, p-value = 0.4659

Selon la valeur de p-value nous avons les hypothèses suivantes —> H0 : la série est un bruit blanc; H1 : la série n’est pas un bruit blanc. La p-value étant > 0.05, on ne peut rejetter H0, la série est un bruit blanc. Vérifions cela avec la representation graphique de la blancheur des données.

Essayons à nouveau de déterminer la moyenne des résidus pour en avoir la certitude sur la blancheur des données.

# Moyenne des résidus
mean(prediction$residuals)
## [1] 0.02804928

Puisque la moyenne des residus est proche de 0, cela montre que nous avons un processus bruit blanc.

Etant conforté dans cette hypothèse, nous pouvons continuer en vérifiant si en plus d’être un bruit blanc, les erreurs suivent un processus gaussien de moyenne nulle (Bruit blanc gaussien centré) avec le test de shapiro.

Les hypothèses sont les suivantes —-> H0 : la série suit une loi normale et H1 : la série ne suit pas une loi normale.

shapiro.test(residus)
## 
##  Shapiro-Wilk normality test
## 
## data:  residus
## W = 0.79483, p-value = 1.825e-09

Puisque la p-value est < 0.05, on rejette H0, la série ne suit pas une loi normale.

B- Détermination du nombre de passagers sur les 12 prochains mois

Utilisons à présent la moyenne des résidus obtenue pour pouvoir déterminer les prédictions des passagers de 12 prochains mois selon le nombre effectif de jours de chaque mois.

data.frame(prediction)
##          Point.Forecast    Lo.90     Hi.90    Lo.95     Hi.95
## Jan 1970       56.94778 48.41257  65.48298 46.77745  67.11810
## Feb 1970       59.53035 50.99514  68.06556 49.36002  69.70068
## Mar 1970       60.06542 51.53021  68.60063 49.89509  70.23575
## Apr 1970       70.85835 62.32314  79.39356 60.68802  81.02868
## May 1970       67.50982 58.97461  76.04503 57.33949  77.68015
## Jun 1970       83.41852 74.88331  91.95374 73.24819  93.58886
## Jul 1970       97.28449 88.74928 105.81971 87.11415 107.45483
## Aug 1970       90.72551 82.19029  99.26073 80.55516 100.89585
## Sep 1970       70.79431 62.25908  79.32954 60.62396  80.96466
## Oct 1970       61.32047 52.78524  69.85571 51.15011  71.49083
## Nov 1970       63.22231 54.68707  71.75755 53.05194  73.39268
## Dec 1970       77.04606 68.51081  85.58131 66.87569  87.21644
write.csv(prediction,"D:/Insseds/Export Rxls/datadef.csv", row.names=TRUE, sep=";" ,dec=",")
setwd("D:/Insseds/Dataframe") #Donner le chemin d'acces au jeu de données
data_def <- read.table("data_def.csv",header=TRUE,sep=",",check.names=FALSE, row.names=1,stringsAsFactors = TRUE) # Lire le jeu de données 
library(moments)
library(DMwR2)
data_def$Lo.80=NULL
data_def$Hi.80=NULL
data_def
##          Point.Forecast    Lo.95     Hi.95
## Jan 1970       56.94778 46.77745  67.11810
## Feb 1970       59.53035 49.36002  69.70068
## Mar 1970       60.06542 49.89509  70.23575
## Apr 1970       70.85835 60.68802  81.02868
## May 1970       67.50982 57.33949  77.68015
## Jun 1970       83.41852 73.24819  93.58886
## Jul 1970       97.28449 87.11415 107.45483
## Aug 1970       90.72551 80.55516 100.89585
## Sep 1970       70.79431 60.62396  80.96466
## Oct 1970       61.32047 51.15011  71.49083
## Nov 1970       63.22231 53.05194  73.39268
## Dec 1970       77.04606 66.87569  87.21644

Rappelons que nous avons au départ procédé à une transformation des données les en ajustant au nombre de jours effectifs de chaque mois. Nous allons donc multiplier le résultat de cette prédiction par le vecteur du nombre de jours dans chaque mois pour obtenir le nombre de passagers prédit pour les 12 prochains mois.

data_def* c(31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)
##          Point.Forecast    Lo.95    Hi.95
## Jan 1970       1765.381 1450.101 2080.661
## Feb 1970       1666.850 1382.081 1951.619
## Mar 1970       1862.028 1546.748 2177.308
## Apr 1970       2125.751 1820.641 2430.860
## May 1970       2092.804 1777.524 2408.085
## Jun 1970       2502.556 2197.446 2807.666
## Jul 1970       3015.819 2700.539 3331.100
## Aug 1970       2812.491 2497.210 3127.771
## Sep 1970       2123.829 1818.719 2428.940
## Oct 1970       1900.935 1585.654 2216.216
## Nov 1970       1896.669 1591.558 2201.780
## Dec 1970       2388.428 2073.146 2703.710

CONCLUSION

Le modèle est acceptable car les résidus sont un bruit blanc mais peut être amélioré car ce n’est pas un blanc gaussien et centré.