Introduction

Dans le cadre de ce projet, nous souhaitons modéliser par un processus une série temporelle de données.

Ici, la série concerne les ventes de champagne de janvier 1962 à septembre 1970, avec une observation par mois.

L’ambition est de pouvoir effectuer une prévision sur les 12 dernières observations grâce à une modélisation sur le reste des données.



Importation des données


On importe les données, puis on en affiche un résumé.

data <- read_delim("data.csv", delim = ";", col_types = cols(date = col_date(format = "%d/%m/%Y")), 
    trim_ws = TRUE)
summary(data)
##       date                value      
##  Min.   :1962-01-01   Min.   :  718  
##  1st Qu.:1964-03-01   1st Qu.: 3086  
##  Median :1966-05-01   Median : 4217  
##  Mean   :1966-05-01   Mean   : 4734  
##  3rd Qu.:1968-07-01   3rd Qu.: 5221  
##  Max.   :1970-09-01   Max.   :13916


Vérifions l’hypothèse d’homogéneité de la variance au cours des années avec le test de bartlett

data$date[1:93] %>% 
  format(format = "%Y") %>%
  bartlett.test(data$value[1:93], .)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  data$value[1:93] and .
## Bartlett's K-squared = 15.224, df = 7, p-value = 0.03323

Avec une p-valeur de 0.03, l’hypothèse H0 de l’homogénéité de la variance au cours des années est rejeté.

Nous aurions pu envisager une transformation de Box-cox.


Préparons nos données

  1. On transforme la base en série temporelle avec la fonction ts() de R.

  2. On précise la fréquence de la série (ici mensuelle) et la date de début de la série.

  3. On divise la série en deux parties : une partie d’apprentissage et une partie de test.

data_ts <- ts(data[,2], start = 1962, frequency = 12)

data_ts_train <- window(data_ts, start = c(1962, 1), end = c(1969, 9))
data_ts_test <- window(data_ts, start = c(1969, 10), end = c(1970, 9))


Affichage de la série



Phase d’identification a priori du modèle

Nous avons 105 observations à notre disposition ce qui est peu, nous rechercherons donc à minimiser le nombre de paramètres à estimer afin de s’assurer de la robustesse de ceux-ci.

Un modèle SARIMA est définit comme suit

\[ \Phi_p(B) \nabla^d \quad \Phi_P(B^S)\nabla^D_S \quad X_t = \Theta_q(B) \Theta_Q(B) \mathcal{E}_t \]

Pour \(1 \leq t \leq T\), les \(X_t\) sont nos observations, et nos \(\mathcal{E}_t\) sont des bruits blancs.

Définissons les paramètres des autres composants par étapes.

## Opérateur Nabla

La série temporelle doit être un processus stationnaire si nous voulons la modéliser par un sarima.

Pour observer cela, on décompose la série temporelle en tendance, en saisonnalité et en bruit.

decompose(data_ts_train) %>% 
  autoplot(main = "Décomposition additive de Xt")

on remarque une forte saisonnalité (ici annuelle) et une tendance

acf(data_ts_train, lag.max = 40, main = "Corrélogramme de Xt", ci.type = "ma")

Les autocorrélations sont fortes pour les \(\rho (h)\) dont les h sont des multiples de 12. Il s’agit d’une saisonnalité annuel donc.

Pour traiter cette saisonnalité annuel, on differencie une première fois en saisonnalité par l’application de l’opérateur \(\nabla ^{12} = 1 - B^{12}\) aux données.

data_diff12 <- diff(data_ts_train, 12)

autoplot(data_diff12, main = "Série temporelle revenue champagne (1-B^12)", xlab = "Année", ylab = "Valeur")


On décompose la série temporelle différenciée pour l’analyser.

decompose_diff12 <- decompose(data_diff12)
autoplot(decompose_diff12 , main = "Décomposition du processus Xt(1-B^12)")


Après application de l’opérateur \(\nabla ^{12}\) on observe que cela n’a pas été suffisant pour éliminer complétement la saisonnalité.

De plus, on observe la disparition de la tendance, cela est dû à la présence du facteur \(1 - B\) dans la formule de l’opérateur nabla en saisonnalité.

\[ \nabla ^{12} = 1 - B^{12} = (1 - B)P(B)\]
Vérifions la disparition de la tendance par le test de Mann-Kendall

MannKendall(data_diff12)
## tau = -0.0954, 2-sided pvalue =0.20906

Avec une p-valeur de 0.209, l’hypothèse H0 “il n’y a aucune tendance” n’est pas rejeté. Tandis que l’hypothèse alternative H1 “il y a une tendance” est rejeté.

Nous avons stationnariser la série et n’appliquerons donc pas un second opérateur nabla sur les données.

Nous avons pu éliminer la tendance, ce qui est positif pour notre modèle. Mais cela induit une augmentation de l’hétérogénéité de la variance qui est elle négative pour notre modélisation.

Quantifions l’augmentation en effectuant à nouveau un test d’homogénéité de la variance entre les années.

df_diff12 <- data.frame(Y=data_diff12, date=round(time(data_diff12)))
bartlett.test(df_diff12$value, df_diff12$date)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  df_diff12$value and df_diff12$date
## Bartlett's K-squared = 22.826, df = 7, p-value = 0.001828

La dégradation est significative, nous sommes passés d’une p-valeur de 0.03 à 0.001. C’est à dire, d’une probabilité d’homogénéité de la variance entre les années de 3% à 0.1%.

## Estimation paramètres p, q, P, Q
Après estimation de nos paramètres d & D, estimons p et q (P et Q pour la saisonnalité).

Pour trouver la valeur du paramètre q, on regarde le corrélograme (autocorrélation) de nos données puisque \(\forall h > q ,\rho (h) = 0\) , ce qui correspond à \(\hat{\rho} (h) \in [interval de Bartlett]\) . Puisque nous ne pouvons que estimer les \(\rho (h)\). A noter que l’interval de Bartlett n’est pas constant.

acf(data_diff12, lag.max = 75, main = "Corrélogramme Xt(1-B^12)",ci.type = "ma")

D’après le corrélograme de \((1 - B^{12})Xt\) on remarque que le processus présente une autocorrélation fortement significative en \(\hat{\rho} (12)\)

Par la corrélation significative de \(\hat{\rho} (12)\), nous pourions implémenter un \(MA_{12}(1)\) dans notre modèle SARIMA. Un \(MA(12)\) n’était pas envisageable car cela induirait l’estimation de 12 paramètres, ce qui est impossible compte tenu du nombre d’observations.

Voyons s’il est plus judicieux d’implémenter une partie autoregressive.

Pour trouver la valeur du paramètre p, on regarde la fonction d’autocorrélation partielle estimée de notre série différencié puisque \(\forall h > p ,r(h) = 0\) , ce qui correspond à \(\hat{r} (h) \in [interval de Quenouille]\) . Puisque nous ne pouvons que estimer les \(r(h)\) 

pacf(data_diff12, lag.max = 75, main = "Fonction d'autocorrélation partielle Xt(1-B^12)") %>% 
  autoplot()

D’après la fonction d’autocorrélation partielle de \((1 - B^{12})Xt\) on remarque que le processus présente une autocorrélation partielle fortement significative en \(\hat{r} (12)\)

Par la corrélation significative de \(\hat{r} (12)\), nous pourions implémenter un \(AR_{12}(1)\) dans notre modèle SARIMA. Un \(AR(12)\) n’était pas envisageable car cela induirait l’estimation de 12 paramètres, ce qui est impossible compte tenu du nombre d’observations (estimations peu robustes).
Nous sommes face à un dilemne, soit nous implémentons un terme \(MA_{12}(1)\) soit un terme \(AR_{12}(1)\), sachant que les valeurs sont toute deux significatives. Pour trancher, regardons les modèles associés :

Implémentation d’un terme \(MA_{12}(1) => SARIMA_{12}[(0,0,0),(0,1,1)]\) \[(1 - B^{12}) Xt = (1- \theta_{12}B^{12})\mathcal{E}_t \]

Implémentation d’un terme \(AR_{12}(1) => SARIMA_{12}[(0,0,0),(1,1,0)]\) \[(1 - B^{12})(1 - \varphi_{12}B^{12})Xt = \mathcal{E}_t \]

Par lecture des 2 modèles possibles, on remarque que sur l’écriture du modèle \(SARIMA_{12}[(0,0,0),(0,1,1)]\) les prédictions ne sont pas faites à partir des données précédemment observés, nous préferons alors implémenter le terme \(AR_{12}(1)\) dans notre modèle.
Création de notre premier modèle

model1 <- arima(data_ts_train, order = c(0,0,0), seasonal = list(order = c(0,1,1), period = 12))
model1
## 
## Call:
## arima(x = data_ts_train, order = c(0, 0, 0), seasonal = list(order = c(0, 1, 
##     1), period = 12))
## 
## Coefficients:
##          sma1
##       -0.1934
## s.e.   0.0956
## 
## sigma^2 estimated as 964549:  log likelihood = -673.23,  aic = 1350.46


Vérifions que la série des résidus \(\hat{\mathcal{E}_t}\) est cohérente avec l’hypothèse selon laquelle \(\mathcal{E}_t\) est un bruit blanc grâce au test de Ljung-Box, il s’agit d’un équivalent du test Portmanteau (un khi²)

Box.test(model1$residuals, lag = 12, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  model1$residuals
## X-squared = 12.831, df = 12, p-value = 0.3815

Avec une p-valeur non significative de 0.38, l’hypothèse H0 “les résidus sont non corrélés” n’est pas rejeté, tandis que l’hypothèse H1 “les résidus sont corrélés” est rejeté. Le test est peu significatif, recherchons un terme simple à rajouter à notre modèle. 

acf(model1$residuals, lag.max = 12, main = "Corrélogramme (I-B^12)(1+0.2B^12)Xt",ci.type = "ma")

pacf(model1$residuals, lag.max = 12, main = "Fonction d'autocorrélation partielle (I-B^12)(1+0.2B^12)Xt") %>% 
  autoplot()

Aucune valeur n’est significative sur la fonction d’autocorrélation partielle, tandis que \(\hat{\rho} (1) = 0.185\) est presque significative sachant l’interval de Bartlett à 0.2251079 pour cette valeur.On accepte alors l’entré de ce terme dans notre modèle.

# Estimation des paramètres du modèle

Notre modèle SARIMA est définit tel que \[SARIMA_{12}[(0,0,1),(1,1,0)]\] \[(1 - B^{12})(1 - \varphi_{12}B^{12})Xt = (1- \theta_{1}B^{1})\mathcal{E}_t \]

model2 <- arima(data_ts_train, order = c(0,0,1), seasonal = list(order = c(1,1,0), period = 12))
model2
## 
## Call:
## arima(x = data_ts_train, order = c(0, 0, 1), seasonal = list(order = c(1, 1, 
##     0), period = 12))
## 
## Coefficients:
##          ma1     sar1
##       0.3578  -0.3361
## s.e.  0.1115   0.1064
## 
## sigma^2 estimated as 841116:  log likelihood = -668.24,  aic = 1342.48

Le paramètre de l’AR(1) de saisonnalité 12 est estimé à -0.34 Et le paramètre du MA(1) est estimé à 0.36. Avec des écart-types respectifs de 0.11 et 0.10

On remplace les paramètres estimés \[(1 - B^{12})(1 + 0.34B^{12})Xt = (1- 0.36B^{1})\mathcal{E}_t \]

# Tests statistiques
## Significativité des coefficients

On test la significativité des coffecients de notre ARIMA. \[ \left\{ \begin{array}{ll} coefficient = 0 & \mbox{(H0)} \\ coefficient \ne 0 & \mbox{(H1)} \end{array} \right. \]

require(lmtest)
## Le chargement a nécessité le package : lmtest
coeftest(model2)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)   
## ma1   0.35775    0.11152  3.2080 0.001336 **
## sar1 -0.33607    0.10643 -3.1576 0.001591 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Le test est hautement signficatif, on rejete H0 “les coefficients sont nulles” et on ne rejete pas H1 “les coefficients sont non nulles”.
## Test des résidus

Vérifions que la série des résidus \(\hat{\mathcal{E}_t}\) est cohérente avec l’hypothèse selon laquelle \(\mathcal{E}_t\) est un bruit blanc grâce au test de Ljung-Box, il s’agit d’un équivalent du test Portmanteau (un khi²)

Box.test(model2$residuals, lag = 12, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  model2$residuals
## X-squared = 5.7129, df = 12, p-value = 0.9299

Avec une p-valeur non significative de 0.93, l’hypothèse H0 “les résidus sont non corrélés (indépendants)” n’est pas rejeté, tandis que l’hypothèse H1 “les résidus sont corrélés” est rejeté. Le test est hautement significatif.

Vérification visuel

autoplot(model2$residuals, main = "Résidus du modèle SARIMA12 (0,0,1)(1,1,0)", xlab = "Année", ylab = "Valeur")

acf(model2$residuals, lag.max = 24, main = "ACF des résidus du modèle SARIMA12 (0,0,1)(1,1,0)", ci.type = "ma") %>% 
  autoplot()

pacf(model2$residuals, main = "PACF des résidus du modèle SARIMA12 (0,0,1)(1,1,0)") %>% 
  autoplot()

On constate visuellement que aucune corrélation ou autocorrélation partielle des résidus n’est présente. De ce fait on peut conclure que le modèle SARIMA12 (0,0,1)(0,1,1) est correct pour notre série temporelle.
## Prédiction sur la série temporelle de test

Nous avons entrainé notre modèle sur l’ensemble des données excepté la dernière année, prédisons cette dernière avec notre modèle.

pred = predict(model2, 12)
data$predict = c(rep(NA,93) , pred$pred )
names(data) = c("date","Xt","prédiction du modèle SARIMA")

Affichage de la prédiction

affichage_pred <- function(data){

don=xts( x=data[,-1], order.by=data$date)

p <- dygraph(don, main = "Affichage de la présivion par modèle SARIMA",ylab = "Revenue vente de champagne") %>%
  dyOptions(labelsUTC = TRUE, fillGraph=TRUE, fillAlpha=0.1) %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 1)
 p}  

affichage_pred(data)



# Choix d’un modèle
Comparons notre modèle avec d’autres modèles qui auraient pu être choisit.

Nous avons créé 3 modèles Sarima suplémentaire, le premier est le choix d’un MA de saisonnalité 12 à la place de l’AR de saisonnalité 12 que nous avions préféré.

Afin d’observer l’impact de la différenciation, pour les 2 autres modèles, l’un sera forcé à n’avoir aucune différenciation tandis que l’autre sera forcé à avoir une différenciation supplémentaire, les estimations des p, q, P, Q ont eux été réalisé comme notre modèle originel.

On obtient alors les modèles suivants :

Sarima AR de saisonnalité 12

\[SARIMA_{12}[(0,0,1),(1,1,0)]\]

Sarima MA de saisonnalité 12

\[SARIMA_{12}[(0,0,1),(0,1,1)]\]

Sarima sans différenciation

\[SARIMA_{12}[(1,0,0),(0,0,4)]\]

Sarima avec différenciations supplémentaires

\[SARIMA_{12}[(2,1,0),(0,3,4)]\]

Affichage des prédictions réalisés par les 4 Sarima

affichage_pred(data[,1:6])


Par évaluation selon le critère d’information (Akaïde), les scores suivants sont obtenus

Modèle Sarima AIC
Sarima AR 1342.48
Sarima MA 1344.16
Sarima sans différenciation 1589.39
Sarima avec différenciations supplémentaires 1032.18

Selon ce critère, on retient le modèle ayant obtenu le minimum de cette quantité. On devrait donc retenir le modèle avec de nombreuses différenciations. Confirmons ce choix en comparant les prédictions sur la dernière année pour ces modèles.

Pour la vérification de nos résultats on utilise les mesures suivantes :

  • MAPE : L’erreur absolue moyenne en pourcentage

  • MAE : L’erreur absolue moyenne

  • MSE : Le carré moyen des erreurs

  • RMSE : L’erreur quadratique moyenne


De plus, pour mettre en perspective les scores obtenus par nos Sarima, nous ajoutons les modèles suivants :

  • meanf : prend la valeur moyenne de toute les observations pour toutes les prédictions,
  • naive : prend la dernière observation pour toutes les prédictions,
  • drift : prend la première et la dernière observations et trace une lignes entre les deux, on utilise la courbe pour les prédictions,
  • snaive : prend la dernière valeur de la saison précédente comme prédiction (ex : sept 1962 = sep 1962 + erreur)
  • Prophet : modèle additif de prédiction temporelle
  • Buys-Ballot : modèle par décomposition de la tendance et de la saisonnalité
    Tableau des résultats.
Sarima AR Sarima MA Sarima sans différenciation Sarima avec différenciation supplémentaire Moyenne Dernière Observation Drift Dernière Saisonnalité Prophet Buys-Ballot
MAPE 0.0922 0.0727 0.1609 0.1863 0.1687 0.4502 0.3791 0.5516 0.6084 0.0681
MAE 427.5 318.8 718.0 929.2 543.6 1464.8 1888.2 2199.5 2387.0 304
MSE 237629.2 139554 639349.4 1239383 489467.1 4792691 9296286 8189765 8760018 118010.3
RMSE 487.47 373.56 799.59 1113.27 699.61 2189.22 3048.98 2861.77 2959.73 343.52


On remarque que nos modèles Sarima ont tous été meilleurs que le modèle Prophet, cela peux s’expliquer par la complexité trop élévé du modèle Prophet. Tandis que le modèle ayant obtenu les meilleurs résultats est le modèle de Buys-Ballot.

Conclusion

Nous sommes très satisfait de notre travail et de nos nouvelles capacités pour trouver et construire des modèles de séries chronologiques