Projet Séries Temporelles
Samir Aidoudi
04/02/2022
Chargement de Libraries usuelles et installation.
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.4 v purrr 0.3.4
## v tibble 3.1.2 v dplyr 1.0.6
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 1.4.0 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layout##
## Attaching package: 'xts'## The following objects are masked from 'package:dplyr':
##
## first, last#1 Analyse transverse et transformation des données
#1.1 Import des donnes et Data quality
Chargement des données des passagers de la SNCF de 1963 à 1980.
sncf=read.table("http://freakonometrics.free.fr/sncf.csv",header=TRUE,sep=";")On regarde les données
View(sncf)
summary(sncf)## ANNEE JANVIER FEVRIER MARS AVRIL
## Min. :1963 Min. :1670 Min. :1560 Min. :1770 Min. :1990
## 1st Qu.:1967 1st Qu.:1816 1st Qu.:1678 1st Qu.:1865 1st Qu.:2151
## Median :1972 Median :2044 Median :1942 Median :2136 Median :2470
## Mean :1972 Mean :2195 Mean :2101 Mean :2309 Mean :2555
## 3rd Qu.:1976 3rd Qu.:2620 3rd Qu.:2546 3rd Qu.:2746 3rd Qu.:2894
## Max. :1980 Max. :3313 Max. :2913 Max. :3306 Max. :3333
## MAI JUIN JUILLET AOUT SEPTEMBRE
## Min. :1910 Min. :2175 Min. :2581 Min. :2639 Min. :1932
## 1st Qu.:2111 1st Qu.:2486 1st Qu.:2939 1st Qu.:2762 1st Qu.:2143
## Median :2272 Median :2703 Median :3165 Median :2880 Median :2228
## Mean :2489 Mean :2888 Mean :3249 Mean :2928 Mean :2426
## 3rd Qu.:2931 3rd Qu.:3394 3rd Qu.:3636 3rd Qu.:3084 3rd Qu.:2820
## Max. :3391 Max. :3682 Max. :3937 Max. :3284 Max. :3206
## OCTOBRE NOVEMBRE DECEMBRE
## Min. :1850 Min. :1630 Min. :2270
## 1st Qu.:1937 1st Qu.:1774 1st Qu.:2445
## Median :2144 Median :1994 Median :2638
## Mean :2359 Mean :2205 Mean :2861
## 3rd Qu.:2790 3rd Qu.:2677 3rd Qu.:3297
## Max. :3269 Max. :3181 Max. :4008Pas de données maquantes ni incoherente nous allons confirmer cela avec les graphes.
Transformation en serie temporelle
train=as.vector(t(as.matrix(sncf[,2:13])))
ts_train=ts(train,start = c(1963, 1), frequency = 12)Voici des informations usuelle concernant la série.
ts_info(ts_train)## The ts_train series is a ts object with 1 variable and 216 observations
## Frequency: 12
## Start time: 1963 1
## End time: 1980 121.2 Graphiques usuels
month plot
Le nombre de passager augmente durant les vacances scolaires.
Observons le Lag plot pour voir l’influence du passé sur la série.
On observe un Lag de 12 mois.
Regardons ce Lag au fil des 3 premieres années.
La série semble avoir à la fois une tendance et une composante saisonnière. Confirmons cela avec un autre graphique.
1.3 Determination du type de Serie
Etude du cycle de la série.
On voit bien que les périodes de vacances scolaires sont les mois ou la fréquentation est le plus élevés.
La série est elle Additive ou Multiplicative ?
Tout d’abord reprenons le graphe.
L’amplitude de la composante saisonnière augmente au fil des années, si la série était additive cette composante serait stationnaire (indépendamment de l’effet de la tendance) La série est donc Multiplicative.
Décomposons la série 
Etudions la saisonalité de la série.
Nous pouvons tester la stationnarité de la serie avec un test de Augmented Dickey-Fuller .
adf.test(ts_train, alternative ="stationary", k=12)##
## Augmented Dickey-Fuller Test
##
## data: ts_train
## Dickey-Fuller = -1.3766, Lag order = 12, p-value = 0.8375
## alternative hypothesis: stationaryEtant donnée la p-valeur nous pouvons conclure à la non stationnarité de la série.
Observons la saisonalité à l’aide de graphiques.
Ces graphiques mettent en évidence la saisonalité de la série.
Etude de la variance et de la moyenne.
Pour cela nous appliquons un test de Box-Pierce.
Box.test(ts_train,lag=20,type="Box-Pierce")##
## Box-Pierce test
##
## data: ts_train
## X-squared = 1089.9, df = 20, p-value < 2.2e-16La probabilité que la série soit un bruit blanc est évidement nulle Cela conforte notre analyse graphique que l’espérance et la variance sont non constante On peut également utiliser un graphique.
## integer(0)Grace au graphique on voit que la variance et la moyenne augmentent également, elles ne sont donc pas stationnaires (comme vu précédemment). #Ainsi afin d’étudier notre série nous devons transformer la série afin de rendre la variance et la moyenne constante.
2 Modélisation
2.1 Analyse des Corrélations
On affiche le graphique ACF pour rechercher le coéfficient AR
acf(diff(diff(log(ts_train),lag=12,difference=1),lag=1,difference=1))
Ainsi en appliquant le critere vu en cours on choisit p=1.
On affiche le graphique PACF pour rechercher le coéfficient MA
pacf(diff(diff(log(ts_train),lag=12,difference=1),lag=1,difference=1))
Ainsi en appliquant le critere vu en cours on choisit q=1.
2.2 Choix du Modèle Manuellement
md1=arima(diff(diff(log(ts_train),lag=1,difference=1),lag=1,difference=1),order=c(1,0,1))
t_stat(md1)## ar1 ma1 intercept
## t.stat -2.424902 -85.62507 0.064286
## p.val 0.015312 0.00000 0.948743Les valeurs des test ne confirme qu’il faut changer l’ordre.
Choix avec d’autres coefficient. Apres avoir chercher de façon empirique, voici le meilleur modèle ARMA en manuel. Ce modèle à le meilleur compromis entre la simplicité et les résultats au test.
md2=arima(diff(diff(log(ts_train),lag=1,difference=1),lag=1,difference=1),order=c(3,0,1))On peut regarder la matrice de correlations des estimateurs.
cor.arma(md2)## ar1 ar2 ar3 ma1 intercept
## ar1 1.000000e+00 2.002232e-01 1.852890e-01 -3.269522e-07 1.516138e-02
## ar2 2.002232e-01 1.000000e+00 2.003658e-01 -2.725740e-07 7.311459e-03
## ar3 1.852890e-01 2.003658e-01 1.000000e+00 -2.114444e-07 1.902618e-03
## ma1 -3.269522e-07 -2.725740e-07 -2.114444e-07 1.000000e+00 -2.073271e-06
## intercept 1.516138e-02 7.311459e-03 1.902618e-03 -2.073271e-06 1.000000e+00Tous les coéfficients sont bien inferieurs à 0.9
t_stat(md2)## ar1 ar2 ar3 ma1 intercept
## t.stat -3.586708 -3.414976 -3.599056 -84.42742 0.056132
## p.val 0.000335 0.000638 0.000319 0.00000 0.955237Les coéfficients ont une p.val faible.
summary(md2)##
## Call:
## arima(x = diff(diff(log(ts_train), lag = 1, difference = 1), lag = 1, difference = 1),
## order = c(3, 0, 1))
##
## Coefficients:
## ar1 ar2 ar3 ma1 intercept
## -0.2391 -0.2281 -0.2396 -1.0000 0e+00
## s.e. 0.0667 0.0668 0.0666 0.0118 1e-04
##
## sigma^2 estimated as 0.02202: log likelihood = 101.26, aic = -190.52
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.003190166 0.1484081 0.1247135 8.572188 198.7046 0.3883776
## ACF1
## Training set 0.01059991AIC =-190
checkresiduals(md2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(3,0,1) with non-zero mean
## Q* = 378.64, df = 19, p-value < 2.2e-16
##
## Model df: 5. Total lags used: 24Les résidus semblent se comporter comme un bruit blanc Le test de Ljung-Box test a une très faible une p-value.
2.3 Modele Auto Arima
On peut sélectionner un modèle automatique de 2 façons : Premièrement avec la fonction armaselect
armaselect(diff(diff(log(ts_train),lag=12,difference=1),lag=1,
difference=1))## p q sbc
## [1,] 15 0 -1150.815
## [2,] 12 0 -1150.277
## [3,] 13 0 -1149.574
## [4,] 0 13 -1145.498
## [5,] 14 0 -1144.649
## [6,] 12 1 -1143.524
## [7,] 2 12 -1142.869
## [8,] 2 13 -1142.462
## [9,] 1 13 -1141.332
## [10,] 0 14 -1141.305Ou avec la fonction auto.arima
mda=auto.arima(diff(diff(log(ts_train),lag=12,difference=1),lag=1,
difference=1))
summary(mda)## Series: diff(diff(log(ts_train), lag = 12, difference = 1), lag = 1, difference = 1)
## ARIMA(1,0,3)(0,0,2)[12] with zero mean
##
## Coefficients:
## ar1 ma1 ma2 ma3 sma1 sma2
## -0.8144 0.1049 -0.6525 -0.2062 -0.5937 0.1419
## s.e. 0.1281 0.1404 0.1120 0.0754 0.0758 0.0744
##
## sigma^2 = 0.002374: log likelihood = 325.48
## AIC=-636.96 AICc=-636.38 BIC=-613.77
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.001728101 0.04800068 0.03713941 -3.935493 296.2438 0.4142767
## ACF1
## Training set -0.008682619#AIC -637
On observe les résidus
checkresiduals(mda)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,3)(0,0,2)[12] with zero mean
## Q* = 22.242, df = 18, p-value = 0.2214
##
## Model df: 6. Total lags used: 24Le modèle semble être cohérent. Les résidus semblent se comporter comme un bruit blanc gaussien. Cependant le test de Ljung-Box test à une p-value de 0.22 essayons d’améliorer ce résultat. Pour se faire essayons de modéliser la série à l’aide d’un modèle SARMA.
2.4 Modèle Auto SARMA
Nous recherchons par la suite un modèle SARIMA en mode automatique. #Pour cela on conserve la saisonnalité et on applique la fonction Auto.arima .
mdsa=auto.arima(diff(log(ts_train),lag=1, difference=1))
summary(mdsa) ## Series: diff(log(ts_train), lag = 1, difference = 1)
## ARIMA(1,0,3)(0,1,2)[12]
##
## Coefficients:
## ar1 ma1 ma2 ma3 sma1 sma2
## -0.8144 0.1049 -0.6525 -0.2062 -0.5937 0.1419
## s.e. 0.1281 0.1404 0.1120 0.0754 0.0758 0.0744
##
## sigma^2 = 0.002374: log likelihood = 325.48
## AIC=-636.96 AICc=-636.38 BIC=-613.77
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.001631537 0.04664192 0.03507763 -89.9094 210.333 0.6466543
## ACF1
## Training set -0.008773206Nous obtenons un AIC de -637.
checkresiduals(mdsa) 
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,3)(0,1,2)[12]
## Q* = 23.439, df = 18, p-value = 0.1743
##
## Model df: 6. Total lags used: 24Le modèle semble être cohérent. Les résidus semblent se comporter comme un bruit blanc gaussien Cependant le test de Ljung-Box test à une p-value de 0.17.
2.5 Conclusion Choix de modèles.
Le meilleur modèle est l’ARMA md2 (3,0,1). Il s’agit du meilleur modèle ceci en se basant sur 2 critères essentiels la parcimonie et les résultats au test.
3 Prédiction
3.1 Prédiction avec le modèle ARMA
fc2=forecast(md2, h = 12)Prévison sur 1 an avec le meilleur modele ARMA.
Nous observons la prévision sur 1 an avec un seuil à 95% et 80%.
3.2 Prédiction avec le modèle SARMA
fcsa=forecast(mdsa, h = 12)Prévison sur 1 an avec le modele SARMA.
Nous observons la prévision sur 1 an avec un seuil à 95% et 80%.
4 Conclusion
Cette étude fu très intéressante nous avons dû retravailler la série notamment en la transformant afin de supprimer les effets de saisonnalité.
Nous obtenons un Modèle ARMA intéressant qui nous permet de faire de bonnes prédictions sur 1 an.