GEST-S615 Techniques quantitatives de gestion - Travail d’examen - Année académique 2016-2017

A propos des auteurs

Ce travail a été réalisé par :

Michael De Shepper

• Diplôme : ingénieur industriel automaticien (ECAM, 2016)
• Email: michaeldeschepper1992@gmail.com
• Contribution à l’examen: modélisation ARIMA sur R

J’affirme sur l’honneur que j’ai effectué ce travail personnellement.

Thierry Vanderborcht

• Diplôme : ingénieur civil mécanicien (ULB, 2009)
• Email: th.vanderborcht@gmail.com
• Contribution à l’examen: méthode du lissage exponentiel de Winters

J’affirme sur l’honneur que j’ai effectué ce travail personnellement.

Antoine Mercier

• Diplôme : ingénieur industriel mécanicien (ISIB, 2006)
• Email: antoinemercier82@hotmail.com
• Contribution à l’examen: régression linéaire multiple

J’affirme sur l’honneur que j’ai effectué ce travail personnellement.

Introduction

Objectifs poursuivis

Par ce travail, nous espérons pouvoir développer nos compétences en matière d’identification de modèles économiques dans le but d’établir des prévisions sur le plus ou moins long terme. De plus, nous souhaiterions apprendre à utiliser le logiciel R, lequel est aujourd’hui considéré comme étant une référence auprès de la communauté scientifique et universitaire en matière d’analyse statistique. La finalité de notre travail consistera donc en une comparaison des différentes méthodes de prévision et sur un retour d’expérience relatif au logiciel utilisé. Pour la modélisation, nous utiliserons les méthodes suivantes:

• Naïve: rapide et permettant d’avoir une référence
• Décomposition en tendance.saisonnalité.cycle: permet de mieux comprendre la structure de la série
• Lissage exponentiel de Winters: adapté aux séries avec saisonnalité
• Régression multiple
• ARIMA: méthode plus complexe permettant une analyse plus poussée de la série

Présentation des données

Les données choisies pour ce projet sont les suivantes : ventes mensuelles de la marque automobile VW sur une période allant de janvier 2006 à mars 2011 et sur territoire EU30. Ces informations proviennent de l’ACEA (Association des Constructeurs Européens de l’Automobile). Nous avons choisi ce type de données pour les raisons suivantes:

• Ce type de série est analysé et suivi dans toutes les grandes entreprises industrielles. Nous pourrions donc être amenés à réaliser des études similaires à l’avenir.
• Dans le cas extrêmement concurrentiel de l’industrie automobile, les ventes mensuelles sont suivies de près par le Top Management car cette information est directement corrélée avec le chiffre d’affaire et la part de marché.
• Les ventes mensuelles d’automobile peuvent être trouvées relativement facilement.
• le choix de cette série de données est également lié au fait qu’un des membres du groupe a travaillé pour un grand constructeur européen d’automobile.

Le tableau de données et le graphique est respectivement donné à la table 1 et à la figure 1. On peut observer la saisonnalité de la série, ainsi qu’une légère tendance à la baisse.

Tableau de données

Mois	2006	2007	2008	2009	2010	2011
JAN	271303	265907	270829	213583	235886	253903
FEV	226633	229069	253351	222769	218716	239155
MARS	357268	370795	328026	319500	352593	367143
AVR	279664	283026	320752	290961	281579
MAI	318973	314255	291786	305685	273026
JUIN	317925	323487	306408	327456	306833
JUIL	290278	300006	289424	282752	245126
AOUT	221968	224906	201706	203945	186444
SEP	297897	280431	284691	290330	279582
OCT	282953	293365	270500	281739	255350
NOV	296925	280185	228256	262066	256915
DEC	252628	247069	226179	229596	229025

Validation des données

Données réservées: les données réservées pour la préevisions sont les celles des 9 derniers mois de la série. Celle-ci s’étendent donc de juillet 2010 à mars 2011.

Critères de validation: pour valider nos modèles, nous utiliserons le critères RMSE et MAPE. Ceux-ci sont donnés par les formules générales:

\[ RMSE = \sqrt{\frac{1}{n}\displaystyle\sum_{i=1}^{N} (\hat{y_i}-y_i)^2} \] \[ MAPE = \frac{1}{n}\frac{\sum_{i=1}^{N}|\hat{y_i}-y_i|}{|y_i|} \]

avec:

• \(n\) = nombre de données
• \(\hat{y_i}\) = donnée prédite
• \(y_i\) = donnée réservée

Pertinence du choix de critères: Le critère RMSE peut s’utiliser pour des erreurs positives et négatives. Il possèdent l’avantage de garder le même ordre de grandeur que les données. Le critère MAPE ne peut s’utiliser que pour des données strictement positives. Il possède l’avantage qu’il se base sur un pourcentage ce qui est pratique à manipuler. Nous l’utiliserons afin de montrer l’influence du choix de l’indicateur sur la validation du modèle en prédiction.

Prévisions

Méthode de prévision naïve

Nous allons commencer par appliquer la méthode de prévision naïve sur notre série. Cette méthode consiste à reprendre la même valeur mensuelle que le même mois de l’année précédente.

Nous réalisons cette exercice sur nos données réservées, à savoir les 9 derniers mois de la série [juillet-2010 ; mars-2011].

Nous obtenons le résultat ci-dessous :

Les résultats concernant nos critères sont les suivants:

• RMSE : 19762
• MAPE : 6.75%

Comme on pouvait s’y attendre, le résultat est relativement bon car notre série ne présente pas une tendance forte.

Dans les paragraphes suivants nous allons raffiner les méthodes utilisées, et nous verrons si nous obtenons une meilleure prévision que par la méthode naïve.

Nous commencerons par décomposer notre série en une « tendance », une « saisonnalité » et « cycle » dans le cadre d’un modèle multiplicatif. Nous savons que la plupart des phénomènes économiques sont de nature multiplicative (comme les ventes de champagne vu au cours), mais afin de confirmer nous avons réalisé l’exercice par la méthode additive et nous avons constaté que les résultats étaient de l’ordre de 50% moins bons.

Calcul de la tendance

Notre série présentant un caractère saisonnier, il est donc nécessaire de calculer la tendance sur base annuelle afin de s’affranchir de cette saisonnalité. Nous avons calculé – par la méthode des moindres carrés – la droite de tendance sur base des ventes mensuelles moyennes de chaque année complète (2006 à 2010).

Les résultats sont les suivants:

• Constante : 293391.32
• Coefficient de pente : -6406.66
• R² : 0.9396

Méthode de prévision naïve améliorée

Connaissant la tendance, nous pouvons améliorer la prévision naïve en y ajoutant cette nouvelle donnée. Dans la suite du rapport, nous appellerons cette méthode la « prévision naïve améliorée ».

Nous obtenons le résultat suivant :

• RMSE : 19171
• MAPE : 6.41%

Soit une légère amélioration par rapport à la prévision naïve « simple ». Passons maintenant à la décomposition saisonnière de notre série de données.

Décomposition saisonnière

Nous avons observé un caractère saisonnier de type « annuel » sur notre série de données. Nous allons maintenant calculer les facteurs saisonniers. Pour ce faire, nous calculons d’abord la moyenne mobile d’ordre 12 (CMA(12)), ensuite le rapport entre la donnée (vente mensuelle) et cette moyenne mobile. On calcule ensuite la moyenne de ce rapport.

Remarque: nous n’avons pas besoin d’utiliser des moyennes élaguées car il n’y a pas de rapport extrême dans nos données.

Les résultats sont repris dans le tableau ci-dessous où la ligne:

• « Provisoire » est la moyenne.
• « Définitif » est la normalisation en base 1 de la moyenne.

Décomposition saisonnière

Année	JAN	FEV	MARS	AVR	MAI	JUIN	JUIL	AOUT	SEP	OCT	NOV	DEC
2006							1.021	0.781	1.046	0.991	1.040	0.885
2007	0.929	0.799	1.296	0.990	1.101	1.137	1.054	0.787	0.984	1.030	0.982	0.871
2008	0.958	0.901	1.170	1.147	1.055	1.120	1.071	0.757	1.074	1.027	0.869	0.856
2009	0.807	0.842	1.206	1.095	1.143	1.217	1.047	0.753	1.067	1.032	0.966	0.853
2010	0.885	0.827
2011
Provisoire	0.895	0.842	1.224	1.078	1.099	1.158	1.048	0.769	1.043	1.020	0.964	0.866
Définitif	0.906	0.834	1.205	1.095	1.100	1.136	1.050	0.768	1.056	1.028	0.973	0.863

Méthode de prévision « Tendance * Saisonnalité * Cycle »

Pour rappel, nous avons choisi un cycle modèle multiplicatif. Pour la valeur du cycle, nous avons cherché l’optimum (sur base du critère MAPE) entre la valeur minimum et la valeur maximum calculée pour notre série. Nous avons remarqué que c’est la valeur maximum du facteur cycle qui minimise le critère MAPE. De plus, nous avons calculé un intervalle de prévision de 80% (basé sur les quantiles 0.10 et 0.90) autour de la prévision. Cette représentation permet de mieux se rendre compte visuellement de la précision de la prévision par rapport aux données réservée. Les résultats sont repris ci-dessous:

Equation:

\[ \hat{y_t} = T_t \times S_t \times C \]

Remarque: l’effet du facteur cycle n’est pas négligeable dans notre approche. En effet, en prenant la valeur minimum du facteur cycle, le critère MAPE passe à 7.96%.

Méthode du lissage exponentiel de Winters

Nous passons maintenant à un niveau de complexité supérieur avec la méthode du lissage exponentiel. Nous utilisons le lissage exponentiel de Winters afin de prendre directement en compte la saisonnalité de notre série de donnée. Remarque: malgré l’utilisation du mot « lissage », il s’agit bien d’une méthode de prévision.

Equations:

Si \(s=12\):

\[ \hat{y_t}(h) = (S_t + hT_t)I_{t+h-12} \text{ si } h \leq 12\] Equations de mise à jour pour niveau, pente de tendance et coefficient saisonnier:

\[ S_t = \alpha (y_t/I_{t-12}) + (1 - \alpha)(S_{t-1} + T_{t-1}) \] \[ T_t = \gamma (S_t - S_{t-1}) + (1 - \gamma)T_{t-1} \text{ inchangée} \] \[ I_t = \delta (y_t / S_t) + (1 - \delta) I_{t-12} \]

Et où les valeurs initiales sont:

\[ S_0 = \text{moyenne}(y_{t-12};y_{t-1})\text{ soit la moyenne des ventes sur les 12 mois précédents} \] \[ T_0 = 0 \] \[ I_0 = y_{t-1}/S_0 \]

Dans le cadre de cette méthode réalisée sous Excel, il est nécessaire d’utiliser le Solver afin de calculer les valeurs optimales de \(\alpha\), \(\gamma\) et \(\delta\) qui minimisent le critère choisi. Le résultat dépend donc du critère que l’on choisit de minimiser ! Attention : il faut réaliser l’optimisation sur base de la série, et non sur les données réservées.

Remarque: Le Solver d’Excel permet de balayer différentes valeurs pour ces 3 paramètres tout en fixant des contraintes sur ceux-ci.

La méthode du lissage exponentiel de Winters s’applique donc bien à notre type de série de données.

Régression multiple

Rappels théoriques

La forme générale d’une régression multiple est:

\[ y_t = \beta_0 + \beta_1x_{1,i} + \beta_2x_{2,i} + \beta_3x_{3,i} + \text{...} + \beta_{k}x_{k,i} + e_i \text{,}\] où \(y_i\) est la variable à prévoir (au sens de la prévision) et \(x_{1,i}\text{,...,}x_{k,i}\) sont les \(k\) variables explicatives. Chaque variable explicative doit être numérique. Les coefficients \(\beta_1\text{,...,}\beta_k\) mesurent l’effet de chaque variable explicative après avoir tenu compte de l’effet de toutes les autres variables explicatives dans le modèle. Par conséquent, les coefficients mesurent l’effet marginal des variables explicatives.

Comme pour une régression linéaire simple, quand on établit des prévisions, les hypothèses suivantes concernant les erreurs \((e_1\text{,...,}e_N)\) doivent être respectées:

• les erreurs ont une moyenne nulle;
• les erreurs ne sont pas correlées les unes avec les autres;
• les erreurs ne sont pas correlées avec aucunes des variables explicatives \(x_{j,i}\)

Il peut être également utile d’avoir une distribution normale des erreurs avec une variance constante pour pouvoir produire des intervalles de prédiction, mais cela n’est pas nécessaire dans le cadre des prévisions.

Estimation du modèle

Les valeurs des coefficients \(\beta_0\text{,...,}\beta_k\) sont obtenues en trouvant le minimum de la somme des erreurs au carré (estimation des “moindres carrés”). Autrement dit, on trouve les valeurs \(\beta_0\text{,...,}\beta_k\) qui minimisent:

\[ \displaystyle\sum_{i=1}^{N} e_i^2 = \displaystyle\sum_{i=1}^{N} (y_i - \beta_0 - \beta_1x_{1,i} - \text{...} - \beta_{k}x_{k,i})^2 \]

En pratique, le calcul est de nos jours toujours réalisé sur ordinateur avec l’utilisation d’outils informatiques tels que Excel ou R dans le cas présent. Trouver la meilleure estimation des coefficients, c’est “fitter”(terme anglais) le modèle aux données.

Quand on fait référence aux coefficients estimés, on utilise la notation \(\hat{\beta_0} \text{,...,} {\hat\beta_k}\).

Valeurs fittées, valeurs de prévision et résidus

Les prévisions de \(y\) peuvent être calculées en ignorant l’erreur dans l’équation de régression. On a donc:

\[ \hat{y} = \hat{\beta_0} + \hat{\beta_1}x_1 + \hat{\beta_2}x_2 + \text{...} + \hat{\beta_k}x_k \] Venir brancher les valeurs de \(x_1 \text{,...,} x_k\) dans la partie droite de cette équation donne une prédiction de \(y\) pour cette combinaison de variables explicatives.

Lorsque ce calcul est fait en utilisant des valeurs de variables explicatives provenant des données qui ont été utilisées pour estimer les modèle, on appelle les valeurs résultantes de \(\hat{y}\) les “valeurs fittées”. Celles-ci sont des “prédictions” des données utilisées pour esimer le modèle. Elles ne représentent pas des prévisions authentiques puisque la valeur réelle de \(y\) pour ce jeu de variables explicatives a été utilisé pour estimer le modèle, par conséquent la valeur de \(\hat{y}\) est affectée par la valeur réelle de \(y\).

Quand les valeurs de \(x_1 \text{,...,} x_k\) sont de nouvelles valeurs (autrement dit, ne faisant pas partie des données utilisées pour estimer le modèle), la valeur résultante de \(\hat{y}\) est une prévision authentique. C’est pour cette raison que, dans le cadre de cet exercice, nous avons réservé les neuf derniers mois de la série VW, afin de pouvoir évaluer de manière authentique la précision de la prévision en la comparant aux données réservées par l’intermédiaire du calcul des valeurs de RMSE et de MAPE.

Par ailleurs, lors de l’évaluation de la précision de la prévision sur base de données réservées, on doit vérifier, comme pour les régressions simples, que les “résidus”, tels que définis ci-après, ont une moyenne nulle et ne sont corrélés avec aucune des variables explicatives. Les “résidus” sont définis comme étant la différence entre les \(y\) observations et les valeurs fittées:

\[ e_i = y_i - \hat{y_i} = y_i - \beta_0 - \beta_1x_{1,i} - \text{...} - \beta_{k}x_{k,i} \]

Application de la régression multiple aux séries temporelles

Le principe de la méthode de régéression multiple peut être facilement étendu aux séries temporelles, en particulier celles qui présentent des composantes remarquables telles qu’une composante de tendance et une composante de saisonalité, comme c’est le cas pour notre série VW.

Lorsque ces deux composantes sont présentes, nous pouvons tirer deux avantages en les modélisant directement. D’une part, il est possible de faire des prévisions au delà de la période qui suit la dernière observation, d’autre part il est possible de comprendre et d’interpréter les composants aux-mêmes et ainsi aboutir à une compréhension plus riche de la série elle-même.

Modélisation de la tendance

Quand une série présente une tendance linéaire, une régression linéaire de \(y_t\) sur le temps peut modéliser la tendance:

\[ y_t = \beta_0 + \beta_1t \] Une caractéristique séduisante d’un modèle basé sur une régression est que les coefficients de \(t\) peuvent être interpretés directement comme le changement de \(y\) (ventes mensuelle de VW) par unité de temps (dans notre cas, des mois).

Modélisation de la saisonnalité

La version la plus simple d’une composante saisonnière est celle qui ajoute une valeur différente à la série (en complément à la tendance) pour chaque saison.

Ce comportement peut être modélisé en introduisant des variables indicatrices pour chaque saison. Par exemple, nous pouvons définir nos variables indicatrices comme étant:

• Saison1 = 1 si nous sommes au mois de Janvier, sinon Saison1 = 0
• Saison2 = 1 si nous sommes au mois de Février, sinon Saison2 = 0
• etc.

Cette idée est une extension des modèles de régression linéaires impliquant des variables indicatrices pour tenir compte de variables catégoriques, par conséquent il est important de rappeler que si l’on a un niveau \(k\) de saisons, nous ne devons introduire que \(k - 1\) variables indicatrices, autrement cela générerait de la colinéarité.

Le choix de la variable (le mois) qui sera ommis n’a pas d’importance. Le modèle va estimer le niveau de vente avec le mois ommis (niveau de base) et les coefficients pour chaque variable indicatrice va estimer le delta vers les haut ou vers le bas dans la série par rapport à ce niveau de base.

Dans notre cas, nous choisirons donc 11 variables indicatrices et nous ommetrons le mois de Janvier. On pourra donc dire que, par exemple, le coefficient associé au mois de Février est une mesure de la différence entre Février par rapport à Janvier sur la variable prédite (la vente mensuelle).

Modélisation de la série VW

Nous pouvons maintenant modéliser la série VW en utilisant un modèle de régression linéraire avec tendance et variables indicatrices mensuelles:

\[ y_t = \beta_0 + \beta_1t + \beta_2d_{2,t} + \beta_3d_{3,t} + \text{...} + \beta_{12}d_{12,t} + e_t \] où \(d_{i,t}=1\) si \(t\) est à la saison \(i\) et \(d_{i,t}=0\) dans les autres cas.

Interprétation des résultats

Le résultat donné par le logiciel R est le suivant:

## 
## Call:
## tslm(formula = VW_training_set ~ trend + season)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -35748  -8871    848   8319  29556 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 263394.8     7505.0  35.096  < 2e-16 ***
## trend         -475.7      132.3  -3.595 0.000862 ***
## season2     -20918.3     9528.0  -2.195 0.033847 *  
## season3      95086.3     9530.7   9.977 1.57e-12 ***
## season4      41122.0     9535.3   4.313 9.90e-05 ***
## season5      51146.3     9541.8   5.360 3.50e-06 ***
## season6      67298.8     9550.0   7.047 1.41e-08 ***
## season7      39113.4    10105.0   3.871 0.000382 ***
## season8     -37894.6    10105.9  -3.750 0.000547 ***
## season9      37787.1    10108.4   3.738 0.000567 ***
## season10     32064.8    10112.8   3.171 0.002876 ** 
## season11     17259.3    10118.8   1.706 0.095637 .  
## season12    -10255.0    10126.6  -1.013 0.317155    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15060 on 41 degrees of freedom
## Multiple R-squared:  0.8912, Adjusted R-squared:  0.8594 
## F-statistic: 27.99 on 12 and 41 DF,  p-value: 6.1e-16

La première colonne donne les estimmations pour chaque coefficient \(\beta\) et la seconde colonne donne son “erreur standard”, c’est-à-dire l’écart-type qui aurait été obtenu en estimant les coefficients \(\beta\) de manière répétitive sur des jeux de données similaires. L’erreur standard fournit une mesure de l’incertitude sur les coefficients \(\beta\) estimés.

Dans le cadre des prévisions, les deux dernières colonnes n’ont pas d’intérêt. La “t value” est le ratio entre un coefficient \(\beta\) et son erreur standard et la dernière colonne donne la “p value”: la probabilité que le coefficient \(\beta\) estimé soit aussi grande si il n’y avait pas de réelle relation entre la vente mensuelle et la variable explicative. Elles sont utiles quand on étudie l’effet de chaque variable explicative, mais pas particulièrement utile dans le cadre des prévisions.

Dans notre cas, on peut donc dire qu’il y a une tendance négative de 475.7 véhicules vendus par mois. Par ailleurs, en moyenne, au mois de Février, il y a une diminution des ventes de 20918.3 véhicules par rapport au mois de Janvier. Au mois de Mars, il y a une augmentation des ventes de 95086.3 véhicules par rapport au mois de Janvier, et ainsi de suite.

Par ailleurs, la valeur “Multiple R-Quared” indique que 89.1% de la variation des ventes mensuelles peut être expliquée par ce modèle.

La figure suivante montre les valeurs réelles (sans les données réservées) en comparaison avec le modèle:

Diagnostic des résidus

Afin de vérifier qu’il n’y a pas d’autocorrélation dans les résidus, commençons par afficher un graphique des résidus:

On constate qu’on a bien un comportement similaire à du bruit blanc (moyenne nulle, pas de “pattern” évident).

Analysons ensuite un corrélogramme des résidus:

En principe, on peut considérer qu’il n’y a pas d’autocorrélation dans les résidus si 95% des pics du corrélogramme restent à l’intérieur de l’intervalle \(\pm 2 / \sqrt{T}\), où \(T\) est la longueur de la série temporelle. Ces bornes sont généralement représentées par des lignes interrompues bleues dans les corrélogrammes.

On constate que le modèle a bien capturé le comportement de la série, puisqu’il ne subsiste pas d’autocorrélation dans les résidus.

Un autre test qui est généralement utilisé dans le cas de modèles de régressions multiples est le test de Durbin-Watson. Il est utilisé pour tester l’hypothèse qu’il n’y a pas d’autocorrélation au retard 1 dans les résidus. S’il n’y a pas d’autocorrélation, le distribution de Durbin-Watson sera symmétrique autour de la valeur 2. Un valeur de p-value faible indique que l’autocorrélation dans les résidus est significative. Le programme R nous retourne le résultat suivant:

## 
##  Durbin-Watson test
## 
## data:  fit_VW
## DW = 2.0138, p-value = 0.988
## alternative hypothesis: true autocorrelation is not 0

D’après les résultats de ce test, on peut considérer que l’autocorrélation dans les résidus n’est pas significative.

Prévisions

La figure suivante indique les prévisions sur un horizon de 9 mois du base du modèle de régression multiple avec les intervalles de prediction de 80% et 95%:

On peut également comparer plus spécifiquement les données réservées avec les prévisions sur l’horizon de 9 mois:

Evaluation de la précision du modèle

La précision de ce modèle est évaluée sur les données réservées et le programme R nous renvoie le résultat suivant:

##           ME         RMSE          MAE          MPE         MAPE 
## 4.717569e+03 2.085133e+04 1.719536e+04 1.161472e+00 6.560648e+00 
##         MASE         ACF1    Theil's U 
## 9.933524e-01 4.407019e-01 2.822771e-01

Méthode ARIMA

Rappels théoriques

Les modèles ARIMA sont une large classe de modèles de regression de séries temporelles.

Contrairement aux modèles de regression multiples étudiés plus haut où l’on régresse linéairement la variable dépendante sur des variables explicatives avec l’hypothèse crutiale que l’erreur est du bruit blanc (en d’autres termes les erreurs sont indépendantes, normales et homoscédastiques), avec les séries temporelles on peut également régresser “aujourd’hui” (variable dépendante) sur “hier” (variable explicative), c’est l’autoregression (AR):

\[ X_t = \phi X_{t-1} + \epsilon_t \] où \(\epsilon_t\) est là encore du bruit blanc (white noise, WN) avec une moyenne nulle.. Ainsi, \(X_t\) peut être expliqué par \(p\) valeurs passées \(X_{t-1},X_{t-2}\text{,...,}X_{t-p}\).

Le modèle de moyenne mobile (Moving Average ou MA), quant à lui, part de l’idée de régresser “aujourd’hui” à partir du “bruit d’hier”:

\[ X_t = \epsilon_t + \theta \epsilon_{t-1} \] où \(\epsilon_t\) est là encore du bruit blanc (white noise, WN) avec une moyenne nulle.

On peut par ailleurs modèliser l’erreur du modèle AR par un modèle MA:

\[ \epsilon_t = W_t + \theta W_{t-1} \]

où \(W_t\) est du bruit blanc. Dans ce cas, l’erreur au temps \(t\) est corrélée avec l’erreur au temps \(t-1\) puisqu’elles ont toutes les deux un \(W_{t-1}\).

Ainsi, on définit le modèle ARMA en combinant la première et la troisième équation:

\[ X_t = \phi X_{t-1} + W_t + \theta W_{t-1} \] En d’autres termes, ce modèle ARMA est constitué d’autogression et d’erreurs autocorrélées.

Une série est stationnaire lorsque:

• sa moyenne est constante au fil du temps
• la structure de corrélation reste constante au fil du temps

Il n’est valide de modéliser une série temporelle par un modèle ARMA que lorsque cete série est stationnaire. En effet, il a été prouvé d’une part, que toute série stationnaire peut être représentée par une combinaison linéaire de bruit blanc:

\[ X_t = W_t + a_1 W_{t-1} + a_2 W_{t-2} + \text{...} \]

D’autre part, il est possible de montrer que tout modèle ARMA est une combinaison de bruit blanc. Par conséquent, les modèles ARMA sont tout à fait indiqués pour capturer la dynamique d’une série stationnaire.

Enfin, une série affiche un comportement de type ARIMA (integrated ARMA) si ses données différenciées affichent un comportement de type ARMA.

Pour des raisons pratiques, on représente un modèle ARIMA par la notation \(\text{ARIMA(p, d, q)}\) où:

• \(p\) est l’ordre AR ordinaire
• \(d\) est le nombre de différences ordinaires, ou ordre d’intégration
• \(q\) est l’ordre MA ordinaire

Comment traiter les séries présentant de la saisonnalité? Considérons une série purement saisonière de prériode S = 12, elle pourrait être modèlisée par un modèle autoregressif saisonnier \(\text{SAR(P = 1)}_{S = 12}\):

\[ X_t = \phi X_{t-12} + W_t \] Cependant, les séries purement saisonières sont rares, et on doit généralement modéliser la partie saisonière avec la partie non saisonière. Par extension des modèles ARIMA, ces modèles sont donc représentés par un modèle \(\text{SARIMA(p, d, q)} \times \text{(P, D, Q)} _S\). Les lettres \(\text{(p, d, q)}\) représentent les composants ordinaires, et les lettres \(\text{(P, D, Q)} _S\) représentent les composants saisonniers.

Lorsque l’on utilise un programme tel que R, le problème revient à déterminer les ordres \(\text{p, d, q, P, D, Q}\) et \(\text{S}\) à partir de l’analyse du comportement de la série elle-même et de l’analyse des corrélogrammes, et le programme se charge de déterminer la valeurs des coefficients.

Modélisation de la série VW

Différenciation

Puisque nous savons déjà que la série VW est saisonière de période 12, nous pouvons appliquer une différence saisonnière de période 12 de la série (\(\nabla_{12} \text{VW}\)) afin de filtrer la saisonnalité et vérifier si la série résultante est stationnaire:

On peut remarquer que puisque l’on a différencié sur une période 12, 12 observations ont été retirées de la série différenciée (les 12 premiers mois), ce qui est bien visible sur ce graphique. Par ailleurs, la série diférenciée semble stationnaire et sa moyenne considérée comme nulle, ce qui nous indique que d = 0 et D = 1.

Au passage, on peut également examiner le corrélogramme de la série initiale pour visualiser la persistance de l’autocorrélation aux retards 1S, 2S, 3S etc.:

En examinant le corrélogramme de la série \(\nabla_{12} \text{VW}\) on peut constater que l’autocorrélation de retard 12 est significative et négative, ce qui est souvent le cas avec les séries présentant de la saisonnalité. Pour la compenser, on aura besoin d’un polynôme saisonier AR ou MA.

Détermination des ordres du modèle ARIMA

Modèle manuel

Le polynôme que nous cherchons est donc de type SARIMA(p,0,q)(P,1,Q)[12] pour lequel p, q, P, Q sont encore à déterminer. Après quelques essais, il a été nécessaire de limiter la taille des coefficients à \(p_{max}=2,q_{max}=2,P_{max}=1,Q_{max}=1\). En effet, des coefficients supérieurs provoquent des non stationnarités ou erreurs de convergence. Avec cette limitation, il reste 36 possibilités de modèle ARIMA à tester. Pour éviter de devoir essayer toutes les solutions une par une, une boucle de programmation a été créée dans R. Son but est de tester toute les possibilités et de calculer pour chacune d’entre elles le critère RSME. Elle permet également de garder en mémoire la configuration pour laquelle le critère est le plus bas. La table ci-dessous reprend les critères RMSE que la boucle a gardé en mémoire.

Critère RMSE pour plusieurs combinaisons d ordres

p	d	q	P	D	Q	RMSE
0	0	0	1	1	1	20758.11
0	0	1	0	1	0	19374.25
0	0	2	0	1	0	16105.61
1	0	2	0	1	0	14563.52

Le meilleur modèle SARIMA est donc celui pour lequel le critère RMSE entre les données prédites et les données réservées est le plus petit. On peut voir que le modèle SARIMA(1,0,2)(0,1,0)[12] a le critère le plus petit. C’est donc celui-ci que nous retiendrons. Le programme R nous retourne

Le résultat donné par le logiciel R pour ce modèle en particulier est le suivant:

##           Estimate       SE t.value p.value
## ar1         0.6127   0.2054  2.9827  0.0044
## ma1        -0.6077   0.2175 -2.7939  0.0074
## ma2         0.4090   0.1739  2.3519  0.0227
## constant -475.0409 498.5834 -0.9528  0.3453

Pour ce modèle, le polynôme Arima est le suivant:

\[ (1-0.6127B)(1-B^{12})y_t=(1 -0.6077B)e_t+(1+0.4090B^2)e_t \]

Diagnostic des résidus

Les principes de l’analyse des résidus sont les mêmes que pour les régressions, c’est-à-dire que l’on doit s’assurer que les résidus sont du bruit blanc. Si ce n’est pas le cas, alors nous n’avons pas trouvé le meilleur modèle.

R fournit les graphiques suivants d’analyse des résidus:

• Le graphique des résidus doit être inspecté pour repérer des ‘patterns’ éventuels. Il est difficile de dire si l’on a du bruit blanc à partir d’un graphe, en revenche il est facile de repérer que ce n’est pas du bruit blanc, par exemple si il y a des patterns évidents dans les rédisus. Dans notre cas, d’après ce graphe on peut conclure que les résidus sont probablement du bruit blanc.

• Le corrélogramme des résidus peut être également utilisé pour estimer la “blancheur” si 95% des valeurs se trouvent à l’interieur des lignes interrompues bleues, ce qui le cas dans ci-dessous.

• Le “Q-Q plot” estime la normalité. Ici, on peut voir que l’hypothèse de normalité est raisonnable.

• Enfin, nous disposons d’un “Q-statistique” qui teste la blancheur dans les résidus, basé sur le test de Ljung-Box. A partir du moment où la majorité des points se trouvent au dessus de la ligne interrompue bleue (ce qui est le cas ici), on peut supposer que le bruit est blanc.

Le modèle est donc satisfaisant et peut donc être raisonnablement considéré comme étant le plus satisfaisant des modèles Arima.

Prévisions

La figure suivante indique les prévisions sur un horizon de 9 mois du base du modèle ARIMA(1,0,2)(0,1,0)[12] avec les intervalles de prediction de 80% et 95%:

La figure suivante montre la superposition des prévisions pour les 9 derniers mois avec les données réservées:

La précision de ce modèle est évaluée sur les données réservées et le programme R nous renvoie le résultat suivant:

##           ME         RMSE          MAE          MPE         MAPE 
## 3.027768e+03 1.456352e+04 1.162627e+04 1.059468e+00 4.502145e+00 
##         MASE         ACF1    Theil's U 
## 6.716339e-01 3.541786e-01 1.958075e-01

Modèle automatique

R possède une fonction (“auto.arima”) automatisant la recherche du moins mauvais modèle ARIMA. Cette fonction teste la stationnarité de la série grâce à un test KPSS, lequel permet de choisir l’ordre de différentiation, laquelle peut ensuite être appliquée à la série initiale. Un processus itératif permet ensuite de déterminer le moins mauvais modèle en minimisant le critère d’Akaike (\(AIC_c\)). Le résultat de la modélisation est un ARIMA(0,0,0)(1,1,0)[12] avec dérive. Son équation est la suivante :

\[ (1+0.4397)(1-B^{12})y_t=-548.1784 \] La figure suivante montre le corrélogramme des résidu de ce modèle.

La précision de ce modèle est évaluée sur les données réservées et le programme R nous renvoie le résultat suivant:

##           ME         RMSE          MAE          MPE         MAPE 
## 8.514228e+03 2.395220e+04 2.024661e+04 2.780155e+00 7.746277e+00 
##         MASE         ACF1    Theil's U 
## 1.169619e+00 4.152390e-01 3.369245e-01

D’après les valeurs de RMSE, on peut voir que le modèle ARIMA obtenu avec la fonction automatique donne un moins bon résultat que celui obtenu “manuellement”. C’est donc ce dernier que nous retiendrons.

Conclusions

Comparaison des méthodes utilisées

Les principaux résultats sont repris dans le tableau ci-dessous:

Résultats

Méthode/Critère	Logiciel	RMSE	MAPE
Naïve	Excel	19762	6.75%
Naïve améliorée	Excel	19171	6.41%
TSC	Excel	23597	6.74%
Lissage exponentiel de Winters	Excel	19123 / 19217	5.17% / 5.13 %
Regression multiple	R	20851	6.56%
ARIMA (1,0,2)(0,1,0)[12]	R	14563	4.50%
ARIMA (0,0,0)(1,1,0)[12]	R	23952	7.74%

On peut donc conclure d’après ce tableau de synthèse que la meilleure methode pour prédire les ventes mensuelles de VW est le modèle ARIMA(1,0,2)(0,1,0)[12].

Comparaison des logiciels utilisés

Les avantages et inconvénients des différents logiciels utilisés sont repris ci-dessous:

Microsoft Excel

Avantages

• Outil généraliste de base largement répandu dans toutes les entreprises.
• Rapidité de prise en main de l’outil.
• Rapidité d’obtention des résultats pour une série de données.

Inconvénients

• Outil généraliste qui a ses limites car il n’est pas spécifique à notre application.
• Risque lié à l’utilisation en mode « boîte noire » de certaines fonctions préinstallées.
• Les feuilles Excel doivent potentiellement être adaptées lorsque l’on change le set de données.

R

En comparaison avec Excel, R est un logiciel beaucoup moins intuitif. En effet, celui-ci est entièrement en ligne de code et nécessite de connaitre les bases de la programmation. De plus, son utilisation nécessite l’apprentissage de commandes très spécifiques, lesquelles peuvent demander beaucoup de temps et de pratique. Une fois ces barrières dépassées, celui-ci permet d’avoir une combinaison d’outil statistiques très complet et très rapide à utiliser. La prévision par des modèles complexes tels que Arima, Holtz Winter,… ne nécessitent en effet que peu de lignes de commande. Cet avantage est également un désavantage, car il est facile d’obtenir un résultat sans devoir connaître exactement le fonctionnement de l’outil utilisé (processus sous jacent).Il est donc nécéssaire de passer beaucoup de temps dans la documentation des fonctions afin de paraméter correctement les options, lesquelles peuvent fortement influencer le résultat.

Un autre grand avantage de ce programme est sa facilité d’adaptation pour traiter une autre série. En comparaison, il est beaucoup plus fastidieux d’adapter un tableur pour une nouvelle série sans devoir tout modifier (réintroduire chaque cellule, changer les formules, adapter les graphiques,…). Nous saluerons au passage la grande communauté présente autour de ce logiciel. Grâce à elle, il est très facile de trouver réponse à ses questions, que ce soit par le biais de forum ou de tutoriaux qui sont constamment mis à jours.

Pour terminer, on signalera que ce document a été entièrement rédigé dans R Studio avec l’application R Markdown, ce qui permet d’intégrer le code complet (dont la majeure partie a été évidemment masquée pour une meilleure visibilité) avec le texte et sa mise en forme. Cela permet une reproductibilité remarquable du travail réalisé, et permet de générer le résultat dans les format html, pdf et MS Word. R Studio permet par ailleurs de travailler de manière collaborative via l’intégration avec la célèbre plateforme Github.com et permet également la publication automatique des documents sur RPubs.com. On peut donc retrouver ce travail également à l’adresse http://rpubs.com/antoinemercier82/240516.