Dans ce tutoriel nous allons mettre en pratique ce que nous avons vu lors du cours introductif sur la modélisation statistique des séries temporelles. Nous allons tenter de modéliser l’évolution du nombre mensuel de personne ayant pris l’avion dans le monde. Ces données sont fournies de base dans R

AirPassengers
##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 112 118 132 129 121 135 148 148 136 119 104 118
## 1950 115 126 141 135 125 149 170 170 158 133 114 140
## 1951 145 150 178 163 172 178 199 199 184 162 146 166
## 1952 171 180 193 181 183 218 230 242 209 191 172 194
## 1953 196 196 236 235 229 243 264 272 237 211 180 201
## 1954 204 188 235 227 234 264 302 293 259 229 203 229
## 1955 242 233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417 391 419 461 472 535 622 606 508 461 390 432

On va également utiliser le package astsa (qu’il vous faudra installer si ce n’est pas déjà fait !)

library(astsa)

1. Prise en main des données

plot(AirPassengers)

Une première inspection de la série temporelle nous informe que :

Au fait, la fréquence d’observation est bien mensuelle ?

frequency(AirPassengers)
## [1] 12

Oui c’est bien le cas (car 12 mois dans l’année)

2. Transformation de la série

Il s’agit ici de rendre stationnaire (à une éventuelle saisonnalité près) la série temporelle.

AirPassengers.trans <- diff(log(AirPassengers))
## le log tente de stabiliser la variance
## la differentiation à l'ordre tente de supprimer la tendance linéaire (petite quadratique)
## AirPassenger.trans est une transformation des données initiales que vous déterminerez
par(mfrow = c(1, 2))
plot(AirPassengers.trans)
lines(lowess(AirPassengers.trans), col = "orange")

plot(abs(AirPassengers.trans))
lines(lowess(abs(AirPassengers.trans)), col = "orange")

Ces deux graphiques nous indiquent que :

3. Identification des ordres

a) Composante saisonniere

Pour cela, nous regardons la série, l’ACF et le PACF de la série temporelle transformée.

plot(AirPassengers.trans)

acf.emp <- acf2(AirPassengers.trans)

Sur ce graphique nous pouvons voir un comportement caractéristique d’une composante saisonnière avec des rebonds aux lags 1, 2, 3, …. Une meilleure manière de procéder (non vue en cours car pas d’analyse spectrale) aurait été de faire

spectrum(AirPassengers.trans, log = "no")

Ici il y a des pics marqués tous les ans, confirmant une saisonnalité annuelle.

Regardons donc l’ACF et le PACF de la série stationnarisé selon cette période.

s <- 12
Yt <- diff(AirPassengers.trans, s)
acf.emp2 <- acf2(Yt, max.lag = 5 * 12)

A l’échelle annuelle (donc aux lags 1, 2, 3, …) on remarque que :

  • l’ACF chute brutalement après le lag 1
  • le PACF n’est pas très net du fait du “rebond” au lag 3.

Ceci suggère un MA(1) pur ou un ARMA. Optons pour le MA(1) pur.

b) Composante non saisonnière

On analyse toujours le même graphique que le précédent, on se concentre alors juste sur une saison. Nous constatons que :

  • l’ACF n’a pas vraiment de décrochage sur une annee
  • le PACF n’a pas vraiment de décrochage sur une année

Ceci suggère donc pour la composante non saisonnière un \(ARIMA(p, d, q)\) avec \(p = 1, d = 1, q = 1\).

p <- q <- d <- 1## tailing off

Pour résumer, ceci suggère donc pour la composante saisonnière un \(ARIMA(P,D,Q)\) avec \(P = 0, D = 1, Q = 1\) et pour la composante non saisonnière un \(ARIMA(p, d, q)\) avec \(p = 1, d = 1, q = 1\).

P <- 0; Q <- 1; D <- 1

4. Ajustement de notre premier modèle

A la vue de notre première étude, il semble logique de s’intéresser au modèle \(ARIMA(p,d,q) \times ARIMA(P, D, Q)_s\).

fit <- sarima(log(AirPassengers), p, d, q, P, D, Q, s)
## initial  value -3.085211 
## iter   2 value -3.225399
## iter   3 value -3.276697
## iter   4 value -3.276902
## iter   5 value -3.282134
## iter   6 value -3.282524
## iter   7 value -3.282990
## iter   8 value -3.286319
## iter   9 value -3.286413
## iter  10 value -3.288141
## iter  11 value -3.288262
## iter  12 value -3.288394
## iter  13 value -3.288768
## iter  14 value -3.288969
## iter  15 value -3.289089
## iter  16 value -3.289094
## iter  17 value -3.289100
## iter  17 value -3.289100
## iter  17 value -3.289100
## final  value -3.289100 
## converged
## initial  value -3.288388 
## iter   2 value -3.288459
## iter   3 value -3.288530
## iter   4 value -3.288649
## iter   5 value -3.288753
## iter   6 value -3.288781
## iter   7 value -3.288784
## iter   7 value -3.288784
## iter   7 value -3.288784
## final  value -3.288784 
## converged

fit <- sarima(log(AirPassengers), p, d, q, P, D, Q, s)
## initial  value -3.085211 
## iter   2 value -3.225399
## iter   3 value -3.276697
## iter   4 value -3.276902
## iter   5 value -3.282134
## iter   6 value -3.282524
## iter   7 value -3.282990
## iter   8 value -3.286319
## iter   9 value -3.286413
## iter  10 value -3.288141
## iter  11 value -3.288262
## iter  12 value -3.288394
## iter  13 value -3.288768
## iter  14 value -3.288969
## iter  15 value -3.289089
## iter  16 value -3.289094
## iter  17 value -3.289100
## iter  17 value -3.289100
## iter  17 value -3.289100
## final  value -3.289100 
## converged
## initial  value -3.288388 
## iter   2 value -3.288459
## iter   3 value -3.288530
## iter   4 value -3.288649
## iter   5 value -3.288753
## iter   6 value -3.288781
## iter   7 value -3.288784
## iter   7 value -3.288784
## iter   7 value -3.288784
## final  value -3.288784 
## converged

Les graphiques produits automatiquement par la librairie astsa lors de l’ajustement nous indiquent que :

Aussi nous pouvons dire que le modèle correct.

Voyons maintenant plus en détail notre modèle ajusté.

fit$ttable
##      Estimate     SE t.value p.value
## ar1    0.1960 0.2475  0.7921  0.4298
## ma1   -0.5784 0.2132 -2.7127  0.0076
## sma1  -0.5643 0.0747 -7.5544  0.0000

Ce tableau nous indique que le coeff ar1 n’est pas stat. significatif, aussi il semblerait pertinent de prendre un ARMA(0, 1, 1) pour la partie mensuelle. On le fait pas je suis à la bourre pour ce cours.

5. Prévisions

Une fois arrivés à cette étape, nous avons clairement fait le plus dur et nous pouvons apprécier à sa juste valeur notre travail éprouvant (en ferais-je trop ?) en visualisant nos prévisions (ici sur 2 ans) :

n.ahead <- 2 * 12
prevision <- sarima.for(log(AirPassengers), n.ahead = n.ahead, p, d, q, P, D, Q, s)

n.ahead <- 2 * 12
prevision <- sarima.for(log(AirPassengers), n.ahead = n.ahead, 0, 1, 1, 0, 1, 1, 12)

Puisque mon chef est pénible et ne sais pas calculer l’exponentielle de tête, je vais être sympa et revenir à l’échelle initiale

plot(ts.union(AirPassengers, exp(prevision$pred),
              exp(prevision$pred + qnorm(0.95) * prevision$se),
              exp(prevision$pred - qnorm(0.95) * prevision$se)),
     plot.type = "single", col = c("black", "red", "grey", "grey"),
     ylab = "Nombre de passagers")

Et voilà notre modélisation est finie et nous pouvons prendre un thé // café ;-)