Analysez et modélisez des séries temporelles

Stationnarisation de la série

On désigne par $\$\backslash (X\_t\backslash )\$$  la série airpass, et on considère $\$\backslash (Y\_t=\backslash ln\backslash left(X\_t\backslash right)\backslash )\$$ . On travaille en effet sur le logarithme de la série afin de pallier l’accroissement de la saisonnalité. On passe ainsi d’un modèle multiplicatif à un modèle additif.

plot(acf(y,lag.max=36,plot=FALSE),ylim=c(-1,1))

La sortie ACF présente une décroissance lente vers 0, ce qui traduit un problème de non-stationnarité. On effectue donc une différenciation $\$\backslash ((I-B)\backslash )\$$ .

y_dif1=diff(y,lag=1,differences=1)
plot(acf(y_dif1,lag.max=36,plot=FALSE),ylim=c(-1,1))

 La sortie ACF de la série ainsi différenciée présente encore une décroissance lente vers 0 pour les multiples de 12. On effectue cette fois la différenciation $\$\backslash (\backslash left(I-B^\{12\}\backslash right)\backslash )\$$ .

y_dif_1_12=diff(y_dif1,lag=12,differences=1)
plot(acf(y_dif_1_12,lag.max=36,plot=FALSE),ylim=c(-1,1))

La sortiE ACF de la série doublement différenciée semble pouvoir être interprétée comme un autocorrélogramme simple empirique. On identifiera donc un modèle ARMA sur la série  $\$\backslash ((I-B)\backslash left(I-B^\{12\}\backslash right)\backslash ln\backslash left(X\_\{t\}\backslash right)\backslash \; .\backslash )\$$ 

Identification, estimation et validation de modèles

On s’appuie sur les autocorrélogrammes simple et partiels estimés :

y_dif_1_12=diff(y_dif1,lag=12,differences=1)
plot(acf(y_dif_1_12,lag.max=36,plot=FALSE),ylim=c(-1,1))

plot(pacf(y_dif_1_12,lag.max=36,plot=FALSE),ylim=c(-1,1))

Modèle 1

On estime en premier lieu un modèle $\$\backslash (\backslash operatorname\{SARIMA\}(1,1,1)(1,1,1)\_\{12\}\backslash )\$$  au vu des autocorrélogrammes empiriques simples et partiels. Ce modèle s’écrit :

$\$\backslash [\backslash left(I-\backslash varphi\_1\; B\backslash right)\backslash left(I-\backslash varphi\_1^\{\backslash prime\}B^\{12\}\backslash right)(I-B)\backslash left(I-B^\{12\}\backslash right)\backslash ln\backslash left(X\_\{t\}\backslash right)=\backslash left(I+\backslash theta\_1\; B\backslash right)\backslash left(I+\backslash theta\_1^\{\backslash prime\}B^\{12\}\backslash right)\backslash varepsilon\_t\backslash \; .\backslash ]\$$

model1=Arima(y,order=c(1,1,1),list(order=c(1,1,1),period=12),include.mean=FALSE,method="CSS-ML")
summary(model1)
## Series: y
## ARIMA(1,1,1)(1,1,1)[12]
##
## Coefficients:
## ar1 ma1 sar1 sma1
## 0.1666 -0.5615 -0.099 -0.4973
## s.e. 0.2459 0.2115 0.154 0.1360
##
## sigma^2 estimated as 0.00138: log likelihood=245.16
## AIC=-480.31 AICc=-479.83 BIC=-465.93
##
## Training set error measures:
## ME RMSE MAE MPE MAPE
## Training set 0.0006239395 0.03489259 0.02595463 0.01199887 0.4696646
## MASE ACF1
## Training set 0.2144266 -0.01250397

t_stat(model1)
## ar1 ma1 sar1 sma1
## t.stat 0.677738 -2.654214 -0.642984 -3.657670
## p.val 0.497938 0.007949 0.520235 0.000255

Box.test.2(model1$residuals,nlag=c(6,12,18,24,30,36),type="Ljung-Box",decim=5)
## Retard p-value
## [1,] 6 0.64051
## [2,] 12 0.81959
## [3,] 18 0.82768
## [4,] 24 0.59646
## [5,] 30 0.75443
## [6,] 36 0.65902

Ce modèle ayant des paramètres non significatifs, on en teste un second.

Modèle 2

Soit le modèle :

$\$\backslash [\backslash left(I-\backslash varphi\_1^\{\backslash prime\}B^\{12\}\backslash right)(I-B)\backslash left(I-B^\{12\}\backslash right)\backslash ln\backslash left(X\_\{t\}\backslash right)=\backslash left(I+\backslash theta\_1\; B\backslash right)\backslash left(I+\backslash theta\_1^\{\backslash prime\}B^\{12\}\backslash right)\backslash varepsilon\_t\backslash \; .\backslash ]\$$

model2=Arima(y,order=c(1,1,1),list(order=c(0,1,1),period=12),include.mean=FALSE,method="CSS-ML")
summary(model2)
## Series: y
## ARIMA(1,1,1)(0,1,1)[12]
##
## Coefficients:
## ar1 ma1 sma1
## 0.1960 -0.5784 -0.5643
## s.e. 0.2475 0.2132 0.0747
##
## sigma^2 estimated as 0.001375: log likelihood=244.95
## AIC=-481.9 AICc=-481.58 BIC=-470.4
##
## Training set error measures:
## ME RMSE MAE MPE MAPE
## Training set 0.0006214569 0.03495868 0.02587023 0.01205357 0.4682205
## MASE ACF1
## Training set 0.2137293 -0.01530519

t_stat(model2)
## ar1 ma1 sma1
## t.stat 0.792074 -2.712668 -7.554412
## p.val 0.428318 0.006674 0.000000

Box.test.2(model2$residuals,nlag=c(6,12,18,24,30,36),type="Ljung-Box",decim=5)
## Retard p-value
## [1,] 6 0.65243
## [2,] 12 0.80854
## [3,] 18 0.82136
## [4,] 24 0.57705
## [5,] 30 0.74768
## [6,] 36 0.67642

Ce modèle ayant des paramètres non significatifs, on en teste un troisième.

Modèle 3

Soit le modèle :

$\$\backslash [(I-B)\backslash left(I-B^\{12\}\backslash right)\backslash ln\backslash left(X\_\{t\}\backslash right)=\backslash left(I+\backslash theta\_1\; B\backslash right)\backslash left(I+\backslash theta\_1^\{\backslash prime\}B^\{12\}\backslash right)\backslash varepsilon\_t\backslash \; .\backslash ]\$$

model3=Arima(y,order=c(0,1,1),list(order=c(0,1,1),period=12),include.mean=FALSE,method="CSS-ML")
summary(model3)

## Series: y 
## ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     sma1
##       -0.4018  -0.5569
## s.e.   0.0896   0.0731
## 
## sigma^2 estimated as 0.001371:  log likelihood=244.7
## AIC=-483.4   AICc=-483.21   BIC=-474.77
## 
## Training set error measures:
##                        ME       RMSE        MAE        MPE      MAPE
## Training set 0.0005730622 0.03504883 0.02626034 0.01098898 0.4752815
##                   MASE       ACF1
## Training set 0.2169522 0.01443892

t_stat(model3)

##              ma1      sma1
## t.stat -4.482494 -7.618978
## p.val   0.000007  0.000000

Box.test.2(model3$residuals,nlag=c(6,12,18,24,30,36),type="Ljung-Box",decim=5)

##      Retard p-value
## [1,]      6 0.51519
## [2,]     12 0.72613
## [3,]     18 0.77822
## [4,]     24 0.50077
## [5,]     30 0.68838
## [6,]     36 0.65352

Les tests de significativité des paramètres et de blancheur du résidu sont validés au niveau 5%.

shapiro.test(model3$residuals)

##
## Shapiro-Wilk normality test
##
## data: model3$residuals
## W = 0.98637, p-value = 0.1674

Le test de normalité est également validé pour ce modèle.

Prévision à l’aide du modèle retenu (3) de l’année 1961

pred_model3=forecast(model3,h=12,level=95)
pred=exp(pred_model3$mean)
pred_l=ts(exp(pred_model3$lower),start=c(1961,1),frequency=12)
pred_u=ts(exp(pred_model3$upper),start=c(1961,1),frequency=12)
ts.plot(x,pred,pred_l,pred_u,xlab="t",ylab="Airpass",col=c(1,2,3,3),lty=c(1,1,2,2),lwd=c(1,3,2,2))

ts.plot(window(x,start=c(1960,1)),pred,pred_l,pred_u,xlab="t",ylab="Airpass",col=c(1,2,3,3),lty=c(1,1,2,2),lwd=c(1,3,2,2))

Analyse a posteriori

On tronque la série de l’année 1960, qu’on cherche ensuite à prévoir à partir de l’historique 1949-1959.

x_tronc=window(x,end=c(1959,12))
y_tronc=log(x_tronc)
x_a_prevoir=window(x,start=c(1960,1))

On vérifie que le modèle 3 sur la série tronquée est toujours validé.

model3tronc=Arima(y_tronc,order=c(0,1,1),list(order=c(0,1,1),period=12),include.mean=FALSE,method="CSS-ML")
summary(model3tronc)

## Series: y_tronc 
## ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
## ma1 sma1
## -0.3484 -0.5623
## s.e. 0.0943 0.0774
## 
## sigma^2 estimated as 0.001338: log likelihood=223.63
## AIC=-441.26 AICc=-441.05 BIC=-432.92
## 
## Training set error measures:
## ME RMSE MAE MPE MAPE
## Training set 0.00104934 0.03443221 0.02590904 0.01899277 0.4738142
## MASE ACF1
## Training set 0.2113963 0.004637637

t_stat(model3tronc)

## ma1 sma1
## t.stat -3.695894 -7.262873
## p.val 0.000219 0.000000

Box.test.2(model3tronc$residuals,nlag=c(6,12,18,24,30,36),type="Ljung-Box",decim=5)

## Retard p-value
## [1,] 6 0.52539
## [2,] 12 0.85631
## [3,] 18 0.87341
## [4,] 24 0.78327
## [5,] 30 0.90181
## [6,] 36 0.84635

shapiro.test(model3tronc$residuals)

## 
## Shapiro-Wilk normality test
## 
## data: model3tronc$residuals
## W = 0.988, p-value = 0.3065

On constate que la réalisation 1960 est dans l’intervalle de prévision à 95% (basé sur les données antérieures à 1959).

pred_model3tronc=forecast(model3tronc,h=12,level=95)
pred_tronc=exp(pred_model3tronc$mean)
pred_l_tronc=ts(exp(pred_model3tronc$lower),start=c(1960,1),frequency=12)
pred_u_tronc=ts(exp(pred_model3tronc$upper),start=c(1960,1),frequency=12)
ts.plot(x_a_prevoir,pred_tronc,pred_l_tronc,pred_u_tronc,xlab="t",ylab="Airpass",col=c(1,2,3,3),lty=c(1,1,2,2),lwd=c(3,3,2,2))
legend("topleft",legend=c("X","X_prev"),col=c(1,2,3,3),lty=c(1,1),lwd=c(3,3))
legend("topright",legend=c("int95%_inf","int95%_sup"),col=c(3,3),lty=c(2,2),lwd=c(2,2))

On calcule les RMSE et MAPE.

rmse=sqrt(mean((x_a_prevoir-pred_tronc)^2))
rmse

## [1] 18.59359

mape=mean(abs(1-pred_tronc/x_a_prevoir))*100
mape

## [1] 2.904473

L’interprétation des critères d’erreur dépend de la série et de la qualité de prévision exigée. Dans le cas présent, un MAPE de 2.9% semble satisfaisant a priori.