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.
# Construction 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)
# Affichage de la serie de données par mois de 1963 à 1969
serie_t <- ts(a, frequency = 12, start = c(1963, 1))
serie_t # Le graphique nous présente une série additive.
## 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
# Structure de la serie temporelle
str(serie_t)
## Time-Series [1:84] from 1963 to 1970: 1750 1560 1820 2090 1910 2410 3140 2850 2090 1850 ...
# Structure de la serie temporelle
attributes(serie_t)
## $tsp
## [1] 1963.000 1969.917 12.000
##
## $class
## [1] "ts"
# Résumé statistique de la serie temporelle
summary(serie_t)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1560 1850 2095 2172 2481 3200
# Visualisation graphique en courbe de la serie temporelle
plot.ts(serie_t,xlab="Années",ylab="Passagers",las=1, main="Tendance de la série temporelle")
# Visualisation graphique en histogramme de la serie temporelle
hist(serie_t,xlab="Passagers",ylab="Fréquence",las=1, main="Fréquence de la série temporelle") #Les valeurs entre 1500 et 2000 sont les plus representatives dans la série de données.
# Visualisation graphique sur une boîte à moustâche de la serie temporelle
boxplot(serie_t,ylab="Passagers",las=1, main="Boîte à moustâche de la série temporelle") #Les valeurs entre 1500 et 2000 sont les plus representatives dans la série de données.
# Moyenne
mean(serie_t) #Les valeurs entre 1500 et 2000 sont les plus representatives dans la série de données.
## [1] 2171.738
# Variance
var(serie_t) #Les valeurs entre 1500 et 2000 sont les plus representatives dans la série de données.
## [1] 178575.5
# Ecart-Type
sd(serie_t) #Les valeurs entre 1500 et 2000 sont les plus representatives dans la série de données.
## [1] 422.5819
library(moments)
skewness(serie_t)
## [1] 0.7182608
kurtosis(serie_t)
## [1] 2.585744
# Autocorellation
acf(serie_t,lag.max=12,plot = TRUE)
acf(serie_t,lag.max=12,plot = FALSE)
##
## Autocorrelations of series 'serie_t', by lag
##
## 0.0000 0.0833 0.1667 0.2500 0.3333 0.4167 0.5000 0.5833 0.6667 0.7500 0.8333
## 1.000 0.495 0.031 -0.199 -0.178 -0.351 -0.458 -0.349 -0.129 -0.144 0.038
## 0.9167 1.0000
## 0.435 0.762
# Autocorellation partielle
pacf(serie_t,lag.max=12,plot = TRUE)
pacf(serie_t,lag.max=12,plot = FALSE)
##
## Partial autocorrelations of series 'serie_t', by lag
##
## 0.0833 0.1667 0.2500 0.3333 0.4167 0.5000 0.5833 0.6667 0.7500 0.8333 0.9167
## 0.495 -0.285 -0.112 0.006 -0.416 -0.235 -0.125 -0.210 -0.483 -0.036 0.245
## 1.0000
## 0.366
Nous allons utiliser la technique de Holt-Winters. 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_t)
## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = serie_t)
##
## Smoothing parameters:
## alpha: 0.01771231
## beta : 0.3099923
## gamma: 0.3346422
##
## Coefficients:
## [,1]
## a 2161.165644
## b -1.650786
## s1 -420.947398
## s2 -485.242941
## s3 -327.845171
## s4 -99.715999
## s5 -109.920570
## s6 276.199906
## s7 661.565776
## s8 539.981819
## s9 -159.893503
## s10 -239.975485
## s11 -328.243309
## s12 254.439558
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
xlisse <- HoltWinters(serie_t,alpha=0.01771231,beta=0.3099923,gamma=0.3346422,seasonal="add")
xlisse
## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = serie_t, alpha = 0.01771231, beta = 0.3099923, gamma = 0.3346422, seasonal = "add")
##
## Smoothing parameters:
## alpha: 0.01771231
## beta : 0.3099923
## gamma: 0.3346422
##
## Coefficients:
## [,1]
## a 2161.165663
## b -1.650788
## s1 -420.947417
## s2 -485.242966
## s3 -327.845193
## s4 -99.716016
## s5 -109.920588
## s6 276.199888
## s7 661.565771
## s8 539.981807
## s9 -159.893516
## s10 -239.975508
## s11 -328.243329
## s12 254.439535
HoltWinters nous a fait les prévisions pour la même période couverte par nos données d’origine.
xlisse$fitted
## xhat level trend season
## Jan 1964 1678.069 2111.439 5.6920163 -439.06250
## Feb 1964 1567.418 2117.697 5.8673397 -556.14583
## Mar 1964 1765.292 2124.141 6.0462343 -364.89583
## Apr 1964 2087.143 2130.802 6.2368059 -49.89583
## May 1964 2067.476 2137.621 6.4172110 -76.56250
## Jun 1964 2435.064 2144.614 6.5957904 283.85417
## Jul 1964 3174.322 2151.652 6.7327036 1015.93750
## Aug 1964 2891.651 2158.839 6.8736926 725.93750
## Sep 1964 2139.277 2166.924 7.2489770 -34.89583
## Oct 1964 1907.286 2175.071 7.5274813 -275.31250
## Nov 1964 1684.782 2181.938 7.3227550 -504.47917
## Dec 1964 2474.082 2190.770 7.7906607 275.52083
## Jan 1965 1773.050 2194.946 6.6701125 -428.56627
## Feb 1965 1660.460 2199.791 6.1042972 -545.43581
## Mar 1965 1858.038 2205.533 5.9919599 -353.48674
## Apr 1965 2176.379 2209.966 5.5085699 -39.09537
## May 1965 2155.428 2215.715 5.5833593 -65.87134
## Jun 1965 2515.791 2218.900 4.8397703 292.05086
## Jul 1965 3255.144 2225.409 5.3570431 1024.37820
## Aug 1965 2983.016 2229.612 4.9993595 748.40495
## Sep 1965 2218.534 2232.432 4.3239177 -18.22240
## Oct 1965 1951.689 2235.365 3.8927142 -287.56901
## Nov 1965 1764.788 2237.811 3.4441871 -476.46668
## Dec 1965 2453.024 2241.170 3.4178953 208.43607
## Jan 1966 1783.407 2242.941 2.9071289 -462.44039
## Feb 1966 1697.211 2246.319 3.0531411 -552.16119
## Mar 1966 1868.671 2248.358 2.7390160 -382.42623
## Apr 1966 2219.017 2250.944 2.6914049 -34.61790
## May 1966 2140.624 2249.579 1.4339436 -110.38838
## Jun 1966 2574.755 2250.470 1.2657947 323.01885
## Jul 1966 3254.232 2250.412 0.8553387 1002.96446
## Aug 1966 2954.888 2247.296 -0.3758469 707.96774
## Sep 1966 2201.233 2245.948 -0.6772175 -44.03760
## Oct 1966 1929.215 2244.540 -0.9036143 -314.42134
## Nov 1966 1764.942 2243.828 -0.8443986 -478.04071
## Dec 1966 2442.853 2260.431 4.5642366 177.85761
## Jan 1967 1813.242 2262.996 3.9445986 -453.69896
## Feb 1967 1700.771 2267.592 4.1464259 -570.96717
## Mar 1967 1888.038 2269.776 3.5382171 -385.27660
## Apr 1967 2166.767 2273.172 3.4940836 -109.89934
## May 1967 2160.708 2277.432 3.7314628 -120.45509
## Jun 1967 2582.164 2280.265 3.4530395 298.44575
## Jul 1967 3214.057 2281.909 2.8920893 929.25609
## Aug 1967 2969.867 2278.884 1.0578896 689.92535
## Sep 1967 2216.450 2274.630 -0.5885849 -57.59148
## Oct 1967 1959.875 2271.979 -1.2279768 -310.87623
## Nov 1967 2114.360 2270.045 -1.4469181 -154.23762
## Dec 1967 2397.602 2260.727 -3.8867592 140.76125
## Jan 1968 1814.178 2259.009 -3.2147107 -441.61600
## Feb 1968 1645.660 2256.145 -3.1058741 -607.37930
## Mar 1968 1865.410 2255.631 -2.3023673 -387.91878
## Apr 1968 2155.213 2253.233 -2.3320721 -95.68797
## May 1968 2111.046 2250.596 -2.4265829 -137.12368
## Jun 1968 2510.698 2248.239 -2.4048715 264.86290
## Jul 1968 3062.280 2245.379 -2.5459683 819.44659
## Aug 1968 2820.475 2234.309 -5.1885231 591.35454
## Sep 1968 2123.851 2225.906 -6.1849438 -95.87047
## Oct 1968 1887.560 2218.200 -6.6563222 -323.98376
## Nov 1968 1905.706 2212.402 -6.3903559 -300.30554
## Dec 1968 2377.793 2203.856 -7.0586059 180.99531
## Jan 1969 1754.945 2197.031 -6.9860894 -435.10020
## Feb 1969 1624.783 2190.808 -6.7496889 -559.27511
## Mar 1969 1792.837 2188.047 -5.5130952 -389.69714
## Apr 1969 2080.041 2185.867 -4.4799530 -101.34612
## May 1969 2041.198 2181.475 -4.4527246 -135.82387
## Jun 1969 2430.814 2178.418 -4.0200499 256.41574
## Jul 1969 2833.017 2175.464 -3.6895860 661.24265
## Aug 1969 2699.809 2171.792 -3.6841888 531.70102
## Sep 1969 2040.917 2168.554 -3.5458710 -124.09085
## Oct 1969 1850.874 2163.079 -4.1438998 -308.06095
## Nov 1969 1819.285 2162.604 -3.0066367 -340.31220
## Dec 1969 2342.779 2160.247 -2.8050447 185.33672
# Faisons la représentation graphique de notre lissage exponentiel
plot(xlisse)
Nous allons prédire le nombre de voyageurs sans intervalle de confiance et avec intervalle de confiance. Ceci permettra de mieux nous ajuster en terme de prise de décisions. L’algorithme nous a permis de calculer les valeurs optimales des parametres.
xlisse <- HoltWinters(serie_t,alpha=0.01771231,beta=0.3099923,gamma=0.3346422,seasonal="add")
serie_t.prev<-predict(xlisse,n.ahead=12,prediction.interval=FALSE)
serie_t.prev
## Jan Feb Mar Apr May Jun Jul Aug
## 1970 1738.567 1672.621 1828.368 2054.846 2042.991 2427.461 2811.176 2687.941
## Sep Oct Nov Dec
## 1970 1986.415 1904.682 1814.764 2395.796
xlisse <- HoltWinters(serie_t,alpha=0.01771231,beta=0.3099923,gamma=0.3346422,seasonal="add")
serie_t.prev<-predict(xlisse,n.ahead=12,prediction.interval=TRUE)
serie_t.prev
## fit upr lwr
## Jan 1970 1738.567 2083.281 1393.854
## Feb 1970 1672.621 2017.427 1327.815
## Mar 1970 1828.368 2173.316 1483.420
## Apr 1970 2054.846 2399.996 1709.697
## May 1970 2042.991 2388.411 1697.571
## Jun 1970 2427.461 2773.231 2081.690
## Jul 1970 2811.176 3157.387 2464.965
## Aug 1970 2687.941 3034.693 2341.189
## Sep 1970 1986.415 2333.817 1639.013
## Oct 1970 1904.682 2252.854 1556.510
## Nov 1970 1814.764 2163.834 1465.693
## Dec 1970 2395.796 2745.903 2045.688
# Représentation graphique de la prévision
par(mfrow=c(1,2), mar=c(3,3,3,3))
ts.plot(serie_t, serie_t.prev, col=c("black","red","green"))
library("forecast")
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
plot(forecast(xlisse))
par(mfrow=c(1,1) , mar=c(0,0,0,0))
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 normal), le modèle serait encore plus intéressant.
Nous allons nous interesser à l’analyse des erreurs de prévision dans l’échantillon de données.
# faire le lissage exponentielle
xlisse <- HoltWinters(serie_t)
# récupérer les résidus
residus=residuals(xlisse)
residus
## Jan Feb Mar Apr May
## 1964 31.9310897 32.5815098 34.7081814 32.8566131 32.5241019
## 1965 -103.0501318 -20.4596420 -88.0382830 13.6211473 -135.4275057
## 1966 26.5927513 -57.2106051 -8.6712478 -229.0173989 -30.6244264
## 1967 36.7581691 -110.7711115 -8.0378878 43.2331339 -50.7083350
## 1968 19.8220788 146.3401552 -5.4100165 -17.2129344 3.9542472
## 1969 43.0548925 225.2168951 188.1629093 4.9590424 78.8016792
## Jun Jul Aug Sep Oct
## 1964 24.9355674 25.6778635 68.3493562 50.7230867 -37.2861517
## 1965 94.2092571 -65.1437868 -123.0160555 -78.5337220 -81.6888035
## 1966 -74.7550262 -224.2318904 -54.8876538 -41.2329033 10.7847831
## 1967 -102.1640616 -334.0569227 -299.8671485 -116.4503981 -39.8750646
## 1968 -25.6975024 -481.2800577 -181.4748967 -85.8506130 48.4396175
## 1969 60.1863588 0.9830088 25.1914038 -108.9170565 207.1261049
## Nov Dec
## 1964 85.2181836 -204.0818890
## 1965 -4.7884300 -93.0242552
## 1966 985.0575202 -112.8526770
## 1967 -444.3605065 122.3980770
## 1968 -121.7062373 13.2072267
## 1969 36.7153074 210.2211246
#Graphe des résidus : Vérification de la blancheur des résidus
plot(residus)
# Autocorrelation
acf(residus, lag.max=12, na.action = na.pass)
Box.test(residus, lag=20, type="Ljung-Box")
##
## Box-Ljung test
##
## data: residus
## X-squared = 20.89, df = 20, p-value = 0.4036
hist(residus)
shapiro.test(residus)
##
## Shapiro-Wilk normality test
##
## data: residus
## W = 0.77549, p-value = 3.97e-09
# La moyenne des erreurs de prédiction
mean(residus, na.rm=TRUE)
## [1] -18.5739