Nous nous intéressons dans ce TP à l’étude du nombre de voyageurs sur le réseau SNCF. Les données sont issues du site https://freakonometrics.hypotheses.org.

On commence par charger et observer les données.

sncf=read.table("http://freakonometrics.free.fr/sncf.csv",header=TRUE,sep=";")
train=as.vector(t(as.matrix(sncf[,2:13])))
X=ts(train,start = c(1963, 1), frequency = 12)
plot(X,xlab='Temps',ylab='Trafic',main="Evolution du nombre de voyageurs SNCF")

lag.plot(X,lags=12,layout=c(3,4),do.lines=FALSE)

monthplot(X)

plot(decompose(X))

Les premières observations de la série nous montrent une tendance générale croissante et une saisonnalité marquée sur l’année, avec un pic en décembre et en juillet, ainsi qu’une forte autocorrélation au lag 12.

L’objectif de ce TP est de se ramener par différentiation à un processus stationnaire, que nous modéliserons.

Première différentiation

Le processus n’est clairement pas stationnaire, observons quand même pour mémoire ce que donne une estimation de sa fonction d’autocorrélation.

acf(X)

Nous commençons par différencier une fois.

DeltaX=diff(X)
plot(DeltaX)

acf(DeltaX)

pacf(DeltaX)

Les fonctions d’autocorrélation et d’autocorrélation partielles ressemblent à celles d’un ARMA(\(p\),\(q\)).

Nous allons estimer un modèle ARMA(1,1)

arima(DeltaX,order=c(1,0,1)) -> fit1
fit1
## 
## Call:
## arima(x = DeltaX, order = c(1, 0, 1))
## 
## Coefficients:
##          ar1      ma1  intercept
##       0.4665  -0.9673     6.5497
## s.e.  0.0654   0.0184     1.6558
## 
## sigma^2 estimated as 114900:  log likelihood = -1558.54,  aic = 3125.07
acf(fit1$residuals)

La fonction d’autocorrélation des résidus n’est pas celle d’un bruit blanc.

Essayons d’augmenter l’ordre de la série

arima(DeltaX,order=c(2,0,2)) -> fit2
fit2
## 
## Call:
## arima(x = DeltaX, order = c(2, 0, 2))
## 
## Coefficients:
##          ar1      ar2      ma1     ma2  intercept
##       1.3223  -0.5936  -1.8884  0.9130     5.8821
## s.e.  0.0588   0.0580   0.0416  0.0377     1.9305
## 
## sigma^2 estimated as 94793:  log likelihood = -1538.98,  aic = 3089.97
acf(fit2$residuals)

fit3 <- arima(DeltaX,order=c(2,0,1))
fit3
## 
## Call:
## arima(x = DeltaX, order = c(2, 0, 1))
## 
## Coefficients:
##          ar1     ar2      ma1  intercept
##       0.5388  -0.189  -0.9560     6.4513
## s.e.  0.0690   0.069   0.0183     1.7195
## 
## sigma^2 estimated as 111004:  log likelihood = -1554.89,  aic = 3119.78
acf(fit3$residuals)

fit4 <- arima(DeltaX,order=c(1,0,2))
fit4
## 
## Call:
## arima(x = DeltaX, order = c(1, 0, 2))
## 
## Coefficients:
##          ar1      ma1      ma2  intercept
##       0.2798  -0.7187  -0.2346     6.5067
## s.e.  0.1161   0.1085   0.1003     1.6944
## 
## sigma^2 estimated as 112434:  log likelihood = -1556.23,  aic = 3122.45
acf(fit4$residuals)

Le meilleur modèle semble être le dernier(ARMA(2,2)), mais les résidus obtenus présentent une autocorrélation au lag 12 (1 an), ce qui indique que nous avons manqué quelque chose dans la modélisation.

Nous pouvons sélectionner automatiquement l’ordre

library("forecast")
auto.arima(DeltaX,d=0,seasonal=FALSE)
## Series: DeltaX 
## ARIMA(2,0,2)            with non-zero mean 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2    mean
##       1.3223  -0.5936  -1.8884  0.9130  5.8821
## s.e.  0.0588   0.0580   0.0416  0.0377  1.9305
## 
## sigma^2 estimated as 97050:  log likelihood=-1538.98
## AIC=3089.97   AICc=3090.37   BIC=3110.19

qui nous donne un ARMA(2,2) également.

Nous tentons une autre approche.

Approche 1 : élimination de la saisonnalité

Delta12DeltaX=diff(DeltaX,lag=12)
plot(Delta12DeltaX)

acf(Delta12DeltaX)

pacf(Delta12DeltaX)

Nous avons encore ici des pics au lag 12.

auto.arima(Delta12DeltaX,d=0,D=0,max.p=12,max.q=12,seasonal=FALSE)-> fits1
acf(fits1$residuals)

pacf(fits1$residuals)

Approche 2 : modèle ARMA saisonnier (SARMA)

Un processus \((X_t)_{t\in\mathbb Z}\) est un \(SARMA(p,q)(P,Q)_s\) s’il est stationnaire et qu’il existe \(c\in\mathbb R\), \((\varphi_1,...,\varphi_p)\in\mathbb R^p\) et \((\theta_1,...,\theta_q)\in\mathbb R^q\) tels que \[ X_t=c+\varphi_1X_{t-1}+...+\varphi_pX_{t-p}+b_t+\theta_1b_{t-1}+...+\theta_q b_{t-q} \] où \((b_t)_{t\in\mathbb Z}\) vérifie une équation de type ARMA(\(P\),\(Q\)) avec une période de \(s\), c’est-à-dire qu’il existe \(c^{(s)}\in\mathbb R\), \((\varphi_1^{(s)},...,\varphi_P^{(s)})\in\mathbb R^P\) et \((\theta_1^{(s)},...,\theta_Q^{(s)})\in\mathbb R^Q\) tel que \[ b_t=c^{(s)}+\varphi_1^{(s)}X_{t-s}+...+\varphi_P^{(s)}X_{t-sP}+Z_t+\theta_1^{(s)}Z_{t-s}+...+\theta_Q^{(s)} Z_{t-sQ}, \] où \((Z_t)_{t\in\mathbb Z}\) est un bruit blanc.

Cette modélisation semble ici raisonnable car les résidus de l’estimation d’un ARMA(\(p\),\(q\)) sur le processus \(\Delta X\) semblent être autocorrélés au lag 12 (1 an).

fitsarma1 <- arima(DeltaX,order=c(1,0,1),seasonal=list(order=c(1,0,1),period=12))
fitsarma1
## 
## Call:
## arima(x = DeltaX, order = c(1, 0, 1), seasonal = list(order = c(1, 0, 1), period = 12))
## 
## Coefficients:
##          ar1      ma1    sar1     sma1  intercept
##       0.2596  -0.9089  0.9879  -0.5514     5.3502
## s.e.  0.0784   0.0314  0.0057   0.0599    12.2772
## 
## sigma^2 estimated as 15362:  log likelihood = -1356.9,  aic = 2725.8
acf(fitsarma1$residuals)

pacf(fitsarma1$residuals)

Effectivement cela semble mieux. Essayons d’obtenir un meilleur modèle :

fitsarma2 <- auto.arima(DeltaX,d=0,D=0)
fitsarma2
## Series: DeltaX 
## ARIMA(0,0,1)(0,0,1)[12] with non-zero mean 
## 
## Coefficients:
##           ma1    sma1    mean
##       -0.9427  0.7778  6.4107
## s.e.   0.0208  0.0456  1.8510
## 
## sigma^2 estimated as 64474:  log likelihood=-1500.36
## AIC=3008.73   AICc=3008.92   BIC=3022.21

La fonction auto.arima sélectionne un modèle dont le critère AIC (ainsi que le critère BIC) est plus grand.

BICsarma1 = -2*fitsarma1$loglik+6*log(length(DeltaX))
BICsarma1
## [1] 2746.024

La fonction n’a donc pas permis de détecter le meilleur modèle. La fonction est en effet conçue pour chercher rapidement le meilleur modèle pour des séries temporelles multivariées. Elle ne calcule pas les coefficients de tous les modèles possibles et fait son choix à partir d’approximations des critères d’information. Pour une recherche exaustive avec un caclul exact des critères d’information

fitsarma3 <- auto.arima(DeltaX,d=0,D=0,stepwise=FALSE,approximation=FALSE)
fitsarma3
## Series: DeltaX 
## ARIMA(0,0,3)(0,0,2)[12] with non-zero mean 
## 
## Coefficients:
##           ma1      ma2      ma3    sma1    sma2    mean
##       -0.5146  -0.3035  -0.1345  0.7561  0.6786  6.1108
## s.e.   0.0730   0.0682   0.0704  0.0617  0.0635  1.6772
## 
## sigma^2 estimated as 37106:  log likelihood=-1442.31
## AIC=2898.62   AICc=2899.16   BIC=2922.22

Le modèle sélectionné est meilleur, mais les valeurs des critères BIC et AIC sont aussi plus élevées que celles du premier modèle, ce qui est étrange.

Approche 3 : modélisation d’un processus de type SARMA sur la série temporelle différenciée deux fois

Reprenons la série Delta12DeltaX.

acf(Delta12DeltaX)

pacf(Delta12DeltaX)

Essayons de modéliser cette série comme un processus SARMA.

fitsarma5 <- arima(Delta12DeltaX,order=c(1,0,1),seasonal=list(order=c(1,0,1),period=12))
fitsarma5
## 
## Call:
## arima(x = Delta12DeltaX, order = c(1, 0, 1), seasonal = list(order = c(1, 0, 
##     1), period = 12))
## 
## Coefficients:
##          ar1      ma1     sar1     sma1  intercept
##       0.2262  -0.8956  -0.2351  -0.3952     0.5503
## s.e.  0.0836   0.0396   0.1306   0.1217     0.6309
## 
## sigma^2 estimated as 15103:  log likelihood = -1267.95,  aic = 2547.9
acf(fitsarma5$residuals)

pacf(fitsarma5$residuals)

Nous obtenons bien ici également un bruit blanc.