Séries Temporelle Projet 21
Jeremy Hazan
13/05/2020
- Attribution de la série temporelle
- Proposition de modèles et estimation des coefficients
- Modèles
- Proposition d’un premier modele MA(7)
- Proposition d’un deuxième modèle AR(6)
- Proposition d’un troisième modèle AR(5)
- Proposition d’un quatrieme modèle ARMA(3,3)
- Proposition d’un cinquième modèle ARMA(3,2)
- Proposition d’un sixième modèle ARMA(1,3)
- Proposition d’un septième modèle ARMA(1,2)
- Comparaison entre les modèles
- Choix de modèle
- Modèles
Attribution de la série temporelle
serie_21 <- scan(file = "/Users/jeremyhazan/Downloads/serie_21.dat")time_serie_21 <- ts(serie_21,frequency = 12,start=c(1990,1))Représentation de la série
Représentation de la série temporelle
plot.ts(serie_21, lwd=1.5, col="dark green", main="Représentation de la série", xlab = "Temps", ylab = "Valeur")
Affichons la boîte à moustache de la série temporelle
boxplot(time_serie_21, horizontal=T)
On remarque que le premier quartile est aux alentours de 530, le troisième quartile est aux alentours de 590. De plus, l’écart interquartile est de 60. En outre, la médiane est aux alentours de 565.
Affichons le lag plot de la série temporelle
lag.plot(time_serie_21, lags = 9, main = "Lagplot de la série temporelle", do.lines = FALSE, diag.col = "red")
Affichons l’histogramme et le Q-Q Plot de la série
par(mfrow=c(1,2))
hist(time_serie_21, main="Histogramme de la série 21", xlab = "Série 21", ylab = "Fréquence")
qqnorm(time_serie_21)
qqline(time_serie_21, col='red')
Décomposition de la série
Décomposition à l’aide de la fonction
Classiquement, on tente de décomposer les séries temporelles en une somme de trois termes: une tendance, une composante saisonnière et du bruit.C’est pour cette raison que nous allons passer à la représentation de la tendance, de la composante saisonnière et du bruit de la série temporelle
decomposition <- decompose(time_serie_21)
par(mfrow=c(3,1))
plot(decomposition$trend,type="l",main="Tendance",col="dark green", xlab = "Temps", ylab = "");
plot(decomposition$seasonal,type="l",main="Saisonnière",col="blue", xlab = "Temps", ylab = "");
plot(decomposition$random,type="l",main="Bruit",col="red", xlab = "Temps", ylab = "")
ACF / PACF
On trace l’ACF et la PACF empirique de la série
par(mfrow=c(1,2))
acf(decomposition$random, na.action = na.pass, ylab = "ACF")
pacf(decomposition$random, na.action = na.pass, ylab = "ACF Partielle")
Densité Spctrale
Verifions la densité spectrale de la série.
spectrum(time_serie_21, col = "purple")
Construction de la série sans tendance
l = length(time_serie_21)
d=12
q=d/2
M=rep(NA, l) #gérer les valeurs interdites en début de fin de la série M
for (i in (q+1):(l-q)) {
M[i] <- 1/12*(1/2*time_serie_21[i-q]+
sum(time_serie_21[(i-q+1):(i+q-1)])+1/2*time_serie_21[i+q])
}
M <- ts(M,start=c(1990,1),frequency = 12)Calcul d’élimination de la tendance
Première série sans tendance.
X_moinsM <-time_serie_21-Mpar(mfrow=c(3,1))
plot(time_serie_21, main = "Représentation de la série temporelle", xlab = "Temps", ylab = "", col = "pink")
plot(M,type='l', main = "1ère tendance calculée", xlab = "Temps", ylab = "", col = "red")
plot(X_moinsM,type='l', main = "Série corrigée de la tendance", xlab = "Temps", ylab = "", col = "purple")
Calcul de saisonnalité
Premier calcul de saisonnalité.
X_moinsM<- c(X_moinsM,rep(0,11))
mat <- matrix(X_moinsM,ncol = 12, byrow = T)
moy <- apply(mat, 2, mean,na.rm=T)
moy## [1] 1.7244355 0.4676266 0.9911289 -5.1126870 -12.6011212 -13.6908173
## [7] -6.1845172 -0.8402665 2.1639540 12.0369811 13.9057959 7.0604641Estimateur de saisonnalité
es <- numeric(12)
for (i in 1:12) {
es[i] <- moy[i]-mean(moy)
}
B <- rep(es,18)
# On supprime les elements pour avoir 205 valeurs
E <- B[1:(length(B)-11)]
ts_es<- ts(E, start=c(1990,1),frequency = 12)Première série désaisonnalisée
X_des<- serie_21-ts_es
par(mfrow = c(3,1))
plot(time_serie_21, main = "Représentation de la série temporelle", xlab = "Temps", ylab = "", col = "pink");
plot(ts_es,type='l', main = "1ère tendance calculée", xlab = "Temps", ylab = "", col = "light green");
plot(X_des,type="l",col="light blue",main="1ère série désaisonalisée", xlab = "Temps", ylab = "X_des")
Deuxième calcul d’élimination de tendance.
Deuxième calcul d’élimination de tendance
M_2=rep(NA, l)
for (i in (q+1):(l-q)) {
M_2[i] <- 1/12*(1/2*X_des[i-q]+sum(X_des[(i-q+1):(i+q-1)])+1/2*X_des[i+q])
}
M_2 <- ts(M_2,start = c(1990,1),frequency = 12)Deuxième série sans tendance
X_moinsM_2 <- serie_21-M_2
par(mfrow = c(3,1))
plot(time_serie_21, main = "Représentation de la série temporelle", xlab = "Temps", ylab = "", col = "pink");
plot(M_2,type='l', main = "Deuxième tendance estimée", xlab = "Temps", ylab = "", col = "light green");
plot(X_moinsM_2,type="l",main="Deuxième série sans tendance",col="light blue")
Deuxième calcul de saisonnalité.
X_moinsM_3 <- c(X_moinsM_2,rep(0,11))
mat2 <- matrix(X_moinsM_3,ncol =12, byrow = T)
moy2 <- apply(mat2, 2, mean,na.rm=T)
moy2## [1] 1.7244355 0.4676266 0.9911289 -5.1126870 -12.6011212 -13.6908173
## [7] -6.1845172 -0.8402665 2.1639540 12.0369811 13.9057959 7.0604641Estimation de la saisonnalité.
es_2 <- numeric(12)
for (i in 1:12) {
es_2[i] <- moy2[i]-1/12*sum(moy2)
}
D <- rep(es_2,18)
F <- D[1:(length(B)-11)]
ts_es_2<- ts(F, start=c(1990,1),frequency = 12)Vérifions si l’ACF et la PACF correspondent aux ACF et PACF obtenus avec la fonction decompose de R.
par(mfrow=c(2,2))
acf(serie_21-M_2-ts_es_2,na.action = na.pass,main="manuelle") #première décomposition
acf(decomposition$random, na.action = na.pass, main="decompose") #deuxième décomp
pacf(serie_21-M_2-ts_es_2,na.action = na.pass, main="manuelle", ylab = "ACF Partielle")
pacf(decomposition$random, na.action = na.pass, main="decompose", ylab = "ACF Partielle")
Constat / Conclusion
On constate que les ACF et PACF réalisés ci-dessus sont identiques avec ceux réalisés avec la fonction decompos de R.
ACF
D’après l’ACF manuel nous parvenons à déterminer que ses coefficients sont nuls à partir de h=8. En effet, à partir de 8 on se retrouve pratiquement constamment dans l’intervalle de confiance. On peut proposer un modèle MA(7).
PACF
Pour l’ACF partielle on peut déterminer que ses coefficients sont nuls à partir de h=6. On peut donc proposer un AR(6)
Résidus
residu=serie_21-M_2-ts_es_2
residu_test <- ts(residu, start = c(1990,8), end = c(2006,2), frequency = 12)
residu_test## Jan Feb Mar Apr May Jun
## 1990
## 1991 NA 30.1337358 36.5281757 8.6222723 -22.4077937 9.9831089
## 1992 14.1018993 1.5897285 -46.6500425 -24.8202939 8.1586550 7.6511887
## 1993 -37.6179403 -54.4205214 -31.0543957 8.8618945 32.6902969 20.6387263
## 1994 4.8785613 -31.4732278 -34.9322160 -19.3607002 -24.0753567 -6.3016878
## 1995 -17.4080345 -23.8032039 -21.7271540 -2.8699137 2.4489092 1.4724233
## 1996 -20.3071794 -10.4080453 15.6939378 32.9974317 52.2775723 48.6847008
## 1997 20.9365492 38.0668805 49.1773773 4.0977526 -13.8902893 -14.7935524
## 1998 3.1651427 22.6751317 -5.0885672 -29.7514754 -16.7118052 0.2938130
## 1999 6.0671033 11.6733743 45.1979539 55.3795717 11.0783106 -15.9483876
## 2000 6.6865640 9.7583282 44.0370768 20.1227261 -14.5723989 -31.0479421
## 2001 -22.0790915 3.1765282 31.4315698 29.0973797 15.2965775 -3.2727904
## 2002 23.1941469 19.7977880 -18.6206859 -16.5776594 -16.3385464 -50.5133204
## 2003 7.9564257 -24.6870033 -58.1271008 -50.2431020 -38.3965220 -12.2528595
## 2004 -2.7970376 -2.9224557 -13.2841444 -26.7976246 -1.5072054 13.0185917
## 2005 12.0412734 15.3027424 1.7192370 -7.8374187 6.2685597 23.2968394
## 2006 2.2585827 -7.0624886
## Jul Aug Sep Oct Nov Dec
## 1990 NA NA NA NA NA
## 1991 13.0527323 -15.5242780 -11.6014565 -2.0916186 -4.7125406 0.2682683
## 1992 2.5281160 6.4377330 17.7284812 18.0810097 20.5840237 -0.1402006
## 1993 -5.7331625 4.2705912 31.6890042 52.6659492 30.0876294 13.4713140
## 1994 10.6396191 0.4038686 -2.5002023 18.1494318 -9.0559913 -24.0552104
## 1995 18.4718355 30.2323455 22.4271338 -5.2960767 2.8710174 -2.6732311
## 1996 36.4607724 -14.6582803 -51.5219278 -75.8959290 -60.5141626 -14.9249558
## 1997 5.4514804 21.7270526 18.4268527 -5.2967146 -29.6415209 -29.4917725
## 1998 -15.4592175 3.4959313 3.6267419 -14.5646059 -1.4032807 2.0236101
## 1999 -12.1873869 -35.4353581 -40.4526419 -16.9395995 4.8959256 26.3844828
## 2000 -38.6749003 -16.0755011 5.3909836 10.2149805 9.3936663 -15.1806825
## 2001 -9.6022505 -9.6641142 -18.0381447 -14.7815734 -1.3084085 17.6287043
## 2002 -30.3320861 24.0342870 48.3051483 60.0695644 65.7991345 42.6808830
## 2003 0.5824681 -6.2991703 7.2650505 19.4549807 24.4632486 17.7850245
## 2004 13.3214466 -6.4360605 -7.5632176 -0.8217312 -15.8192209 -15.6226140
## 2005 9.6881396 -11.3808467 -33.4738538 -19.9631639 -3.5550660 2.4203663
## 2006y <- window(residu, start=c(1990,8),end=c(2006,8), frequency=12)
length(y)## [1] 193On découpe la série (residu), et on garde les 6 dernières valeurs qui serviront à comparer la performance des modèles proposés.
new_serie=window(y, start=c(1990,8),end=c(2006,2),frequency=12)
vrai=window(residu, start=c(2006,3),c(2006,8),frequency=12)
print(c(length(new_serie),length(vrai)))## [1] 187 6Proposition de modèles et estimation des coefficients
Modèles
Proposition d’un premier modele MA(7)
On commence par proposer un MA(7) puis on essaiera de simplifier afin de diminuer le nombre de paramètres à estimer.
(mod_ARM=arima(new_serie,order=c(0,0,7)))##
## Call:
## arima(x = new_serie, order = c(0, 0, 7))
##
## Coefficients:
## ma1 ma2 ma3 ma4 ma5 ma6 ma7 intercept
## 1.0115 0.1112 -0.3600 -0.5632 -0.6286 -0.4114 -0.1594 0.0228
## s.e. 0.0744 0.1100 0.0977 0.0816 0.1005 0.1045 0.0826 0.1491
##
## sigma^2 estimated as 169.4: log likelihood = -747.84, aic = 1513.69Test de la blancheur des résidus
acf(mod_ARM$residuals, na.action = na.pass)
#pacf(mod_ARIMA$residuals, na.action = na.pass)Tous les coefficients sont à l’intérieur de l’intervalle de confiance, donc on accepte la blancheur des résidus.
Faisons une prédiction puis calculons l’erreur de prédiction.
prediction =predict(mod_ARM,n.ahead=5)
#erreur quadratique moyenne (RMSE)
erreur = sqrt(mean((prediction$pred-vrai)^2))
erreur## [1] 12.00309A l’aide de des intervalles de confiance, nous allons tester la significativité des coefficients.
confint(mod_ARM)## 2.5 % 97.5 %
## ma1 0.8656215 1.157416092
## ma2 -0.1043761 0.326860789
## ma3 -0.5515010 -0.168594856
## ma4 -0.7232238 -0.403194179
## ma5 -0.8255352 -0.431718352
## ma6 -0.6162976 -0.206547159
## ma7 -0.3213922 0.002493277
## intercept -0.2693509 0.3150131710 appartient aux intervalles de confiances de ma1 et ma2
Proposition d’un deuxième modèle AR(6)
(mod_ARM2=arima(new_serie,order=c(6,0,0)))##
## Call:
## arima(x = new_serie, order = c(6, 0, 0))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 intercept
## 1.0944 -0.9246 0.5377 -0.5257 0.1963 -0.1753 0.1702
## s.e. 0.0721 0.1082 0.1240 0.1251 0.1109 0.0741 1.2664
##
## sigma^2 estimated as 187: log likelihood = -755.52, aic = 1527.04Testons la blancheur des résidus
acf(mod_ARM2$residuals,na.action = na.pass)
#pacf(mod_ARM2$residuals,na.action = na.pass)On accepte la blancheur des résidus
prediction2=predict(mod_ARM2,n.ahead=5)
erreur2 = sqrt(mean((prediction2$pred-vrai)^2))
erreur2## [1] 11.69413Testons la significativité des coefficients
confint(mod_ARM2)## 2.5 % 97.5 %
## ar1 0.9532157 1.23565795
## ar2 -1.1365911 -0.71250960
## ar3 0.2945907 0.78084511
## ar4 -0.7709290 -0.28050020
## ar5 -0.0211709 0.41372073
## ar6 -0.3206012 -0.02995028
## intercept -2.3119560 2.65239275Proposition d’un troisième modèle AR(5)
(mod_ARM3=arima(new_serie,order=c(5,0,0)))##
## Call:
## arima(x = new_serie, order = c(5, 0, 0))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 intercept
## 1.0935 -0.8614 0.4626 -0.3753 0.0013 0.2303
## s.e. 0.0732 0.1064 0.1215 0.1094 0.0753 1.5041
##
## sigma^2 estimated as 192.8: log likelihood = -758.27, aic = 1530.54Testons la blancheur des résidus
acf(mod_ARM3$residuals, na.action=na.pass)
On accepte la blancheur des résidus
Prévisions et calcul de l’erreur de prévisions
prediction3=predict(mod_ARM3,n.ahead=5,ci=TRUE)
erreur3 = sqrt(mean((prediction3$pred-vrai)^2))
erreur3## [1] 15.0377Testons la significativité des coefficients
confint(mod_ARM3)## 2.5 % 97.5 %
## ar1 0.9500866 1.2368754
## ar2 -1.0699459 -0.6529023
## ar3 0.2244102 0.7008002
## ar4 -0.5897233 -0.1609311
## ar5 -0.1463245 0.1490172
## intercept -2.7175893 3.1782571Proposition d’un quatrieme modèle ARMA(3,3)
(mod_ARM4=arima(new_serie,order=c(3,0,3)))##
## Call:
## arima(x = new_serie, order = c(3, 0, 3))
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 ma3 intercept
## 0.4427 0.5597 -0.5157 0.5668 -0.9471 -0.6197 0.0032
## s.e. 0.1579 0.2042 0.1172 0.1533 0.0590 0.1119 0.1064
##
## sigma^2 estimated as 166.9: log likelihood = -746.71, aic = 1509.42Testons la blancheur des résidus
acf(mod_ARM4$residuals, na.action=na.pass)
On accepte la blancheur des résidus
Prévisions et calcul de l’erreur de prévisions
prediction4=predict(mod_ARM4,n.ahead=5)
erreur4 = sqrt(mean((prediction4$pred-vrai)^2))
erreur4## [1] 6.582656Test de significativité des coefficients
confint(mod_ARM4)## 2.5 % 97.5 %
## ar1 0.1333286 0.7521227
## ar2 0.1595223 0.9598931
## ar3 -0.7453723 -0.2861204
## ma1 0.2663837 0.8672714
## ma2 -1.0626314 -0.8314767
## ma3 -0.8391005 -0.4003424
## intercept -0.2052597 0.2116634Proposition d’un cinquième modèle ARMA(3,2)
(mod_ARM5=arima(new_serie,order=c(3,0,2)))##
## Call:
## arima(x = new_serie, order = c(3, 0, 2))
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 intercept
## 1.2597 -0.4673 -0.0907 -0.2989 -0.7011 0.0019
## s.e. 0.1059 0.1693 0.1012 0.0746 0.0739 0.0981
##
## sigma^2 estimated as 166.6: log likelihood = -746.61, aic = 1507.22Testons la blancheur des résidus
acf(mod_ARM5$residuals, na.action=na.pass)
On accepte la blancheur des résidus
Prévisions et calcul de l’erreur de prévisions
prediction5=predict(mod_ARM5,n.ahead=5)
erreur5 = sqrt(mean((prediction5$pred-vrai)^2))
erreur5## [1] 6.964171Test de significativité des coefficients
confint(mod_ARM5)## 2.5 % 97.5 %
## ar1 1.0520500 1.4673308
## ar2 -0.7990518 -0.1355758
## ar3 -0.2891767 0.1077128
## ma1 -0.4450673 -0.1527059
## ma2 -0.8459727 -0.5562448
## intercept -0.1903805 0.1941961Proposition d’un sixième modèle ARMA(1,3)
(mod_ARM6=arima(new_serie,order=c(1,0,3)))##
## Call:
## arima(x = new_serie, order = c(1, 0, 3))
##
## Coefficients:
## ar1 ma1 ma2 ma3 intercept
## 0.7014 0.2765 -0.7912 -0.4853 0.0199
## s.e. 0.0596 0.0711 0.0385 0.0640 0.1622
##
## sigma^2 estimated as 187.5: log likelihood = -757.13, aic = 1526.27Testons la blancheur des résidus
acf(mod_ARM6$residuals, na.action=na.pass)
On accepte la blancheur des résidus
Prévisions et calcul de l’erreur de prévisions
prediction6=predict(mod_ARM6,n.ahead=5)
erreur6 = sqrt(mean((prediction6$pred-vrai)^2))
erreur6## [1] 13.30304Test de significativité des coefficients
confint(mod_ARM6)## 2.5 % 97.5 %
## ar1 0.5844682 0.8182696
## ma1 0.1371406 0.4159183
## ma2 -0.8667676 -0.7156865
## ma3 -0.6106971 -0.3599070
## intercept -0.2979944 0.3378018Proposition d’un septième modèle ARMA(1,2)
(mod_ARM7=arima(new_serie,order=c(1,0,2)))##
## Call:
## arima(x = new_serie, order = c(1, 0, 2))
##
## Coefficients:
## ar1 ma1 ma2 intercept
## 0.1731 1.0396 0.2870 0.4674
## s.e. 0.1706 0.1667 0.1399 3.0202
##
## sigma^2 estimated as 217.5: log likelihood = -769.38, aic = 1548.75Testons la blancheur des résidus
acf(mod_ARM7$residuals, na.action=na.pass)
On accepte la blancheur des résidus
Prévisions et calcul de l’erreur de prévisions
prediction7=predict(mod_ARM7,n.ahead=5)
erreur7 = sqrt(mean((prediction7$pred-vrai)^2))
erreur7## [1] 22.11016Test de significativité des coefficients
confint(mod_ARM7)## 2.5 % 97.5 %
## ar1 -0.16131118 0.5075023
## ma1 0.71283225 1.3663292
## ma2 0.01288539 0.5611026
## intercept -5.45218635 6.3868915Comparaison entre les modèles
Modèle <- c("MA(7)","AR(6)","AR(5)","ARMA(3,3)","ARMA(3,2)","ARMA(1,3)","ARMA(1,2)")
AIC <- c(1535.22,1548.39,1552.13,1530.75,1528.2,1571.77,1569.88)
RMSE <- c(18.32,17.36,21.48,12.19,11.23,30.61,30.68)
df=data.frame(Modèle=Modèle,AIC=AIC,RMSE=RMSE)
row.names(df) <- paste("MOD",1:7, sep="")
df## Modèle AIC RMSE
## MOD1 MA(7) 1535.22 18.32
## MOD2 AR(6) 1548.39 17.36
## MOD3 AR(5) 1552.13 21.48
## MOD4 ARMA(3,3) 1530.75 12.19
## MOD5 ARMA(3,2) 1528.20 11.23
## MOD6 ARMA(1,3) 1571.77 30.61
## MOD7 ARMA(1,2) 1569.88 30.68Choix de modèle
En se basant sur le critère de l’AIC, le meilleur modèle est celui qui a le plus faible AIC, donc pour cette série il s’agit du ARMA(3,2).
Affichage du meilleur modèle
Affichons sur un même graphique les données réelles (y) ainsi que les valeurs prédites.
plot.ts(y,type='o',xlab = "Temps", ylab="residu",xlim=c(1990,2006),main='ARMA(3,2)', col='orange')
lines(prediction5$pred,col='green',type='o')
lines(vrai,col='pink',type='o')
legend(x='bottomleft',legend=c('données réelles','prévisions'),col=c('pink','green'),lty=1)
Affichage du pire modèle
Le graphique suivant nous permet d’illustrer le cas ou les prévisions sont les plus éloignées des données réelles.
plot.ts(y,type='o',xlab = "Temps", ylab="residu",xlim=c(1990,2006),main='ARMA(1,2)', col='orange')
lines(prediction7$pred,col='green',type='o')
lines(vrai,col='pink',type='o')
legend(x='bottomleft',legend=c('données réelles','prévisions'),col=c('pink','green'),lty=1)