P OLYTECH ’L ILLE D ÉPARTEMENT G.I.S. 2007-2008
Introduction aux séries temporelles Julien JACQUES
Table des matières 1 Introduction et premières définitions 1.1 Tendances et composantes saisonnières . 1.2 Indices descriptifs d’une série temporelle 1.2.1 Indices de tendances centrales . . 1.2.2 Indices de dispersion . . . . . . . 1.2.3 Indices de dépendance . . . . . . 1.3 Mise en oeuvre sous R . . . . . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
3 5 5 5 6 6 8
2 Lissages exponentiels 2.1 Lissage exponentiel simple . . . . . . . . . 2.2 Lissage exponentiel double . . . . . . . . . 2.3 Méthode de Holt-Winters . . . . . . . . . . 2.3.1 Méthode non saisonnière . . . . . . 2.3.2 Méthode saisonnière additive . . . . 2.3.3 Méthode saisonnière multiplicative 2.4 Mise en oeuvre sous R . . . . . . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
9 9 9 10 10 10 10 11
3 Estimation et élimination de la tendance et de la saisonnalité 3.1 Processus stationnaire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Une estimation paramétrique de la tendance . . . . . . . . . . . . . . . . . . 3.3 Estimation non paramétrique : moyenne mobile . . . . . . . . . . . . . . . . 3.4 Elimination de la tendance et de la saisonnalité par la méthode des différences 3.5 Mise en oeuvre sous R . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
16 16 16 17 20 20
4 Modélisation des séries stationnaires 4.1 Auto-corrélation partielle . . . . . . . . . . . . . . . . 4.2 Quelques exemples simples de processus stationnaires. 4.2.1 Bruit blanc . . . . . . . . . . . . . . . . . . . 4.2.2 Moyenne mobile . . . . . . . . . . . . . . . . 4.2.3 Processus auto-régressif d’ordre 1 : AR1 . . . 4.3 Les processus linéaires stationnaires . . . . . . . . . . 4.3.1 Les processus moyenne mobile . . . . . . . . 4.3.2 Les processus auto-régressif . . . . . . . . . . 4.3.3 Les processus mixte ARM Ap,q . . . . . . . . 4.3.4 Estimation, prévision . . . . . . . . . . . . . . 4.3.5 Illustrations . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
22 22 22 22 22 23 23 24 24 25 25 26
5 Processus non stationnaire : ARIMA et SARIMA 5.1 Mise en oeuvre sous R : processus ARMA, ARIMA, SARIMA . . . . . . . . . . . . . . . . . . . . .
35 36
6 Processus ARCH et GARCH 6.1 Définitions . . . . . . . . . . . . 6.2 Rappels de probabilités . . . . . 6.3 Propriétés des processus ARCH 6.4 Processus GARCH et propriétés 6.5 Mise en oeuvre sous R . . . . .
37 38 38 38 39 39
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . . . . . . . .
. . . . .
. . . . . . . . . . .
. . . . .
. . . . . . . . . . .
. . . . .
. . . . . . . . . . .
. . . . .
. . . . . . . . . . .
. . . . .
. . . . . . . . . . .
. . . . .
. . . . . . . . . . .
. . . . .
. . . . . . . . . . .
. . . . .
. . . . . . . . . . .
. . . . .
. . . . . . . . . . .
. . . . .
. . . . . . . . . . .
. . . . .
. . . . . . . . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
Les différentes démonstrations nécessaires à ce cours sur les séries temporelles ne figurent pas dans ce document, et sont généralement données en exercice. Elles seront néanmoins faites en cours.
1 Introduction et premières définitions Une série temporelle (ou série chronologique) est une suite finie (xt )1≤t≤n , où t représente le temps (en minute, jour, année...).
9000 7000
8000
USAccDeaths
10000
11000
Voici quelques exemples de séries temporelles : – Ex. 1 : Nombre de morts accidentelles aux Etats-Unis de 1973 à 1978 (figure 1).
1973
1974
1975
1976
1977
1978
1979
Time
F IG . 1 – Nombre de morts accidentelles aux Etats-Unis de 1973 à 1978
400 300 200 100
AirPassengers
500
600
– Ex. 2 : Nombre de passagers par mois (en milliers) dans les transports aériens, de 1949 à 1960 (figure 2).
1950
1952
1954
1956
1958
1960
Time
F IG . 2 – Nombre de passagers (en milliers) dans les transports aériens
3
100 0
50
sunspot.year
150
– Ex. 3 : Nombre annuel de tâches solaires observées à la surface du soleil de 1700 à 1980 (figure 3).
1700
1750
1800
1850
1900
1950
Time
F IG . 3 – Nombre annuel de tâches solaires – Ex. 4 : Taille de la population française (en milliers) de 1985 à 2005 (figure 4).
F IG . 4 – Population française de 1985 à 2005
3500 3000 2500 1500
2000
EuStockMarkets[, 3]
4000
– Ex. 5 : Valeurs de clôtures journalières du CAC40 de 1991 à 1998 (figure 5).
1992
1993
1994
1995
1996
1997
1998
Time
F IG . 5 – Valeurs de clôtures journalières du CAC40 de 1991 à 1998 Excepté l’exemple 4, ces données sont disponibles dans le logiciel R sous les noms : EuStockMarkets, USAccDeaths, AirPassengers et sunspot.year. Exercice : repérer les tendances (croissance, décroissance, linéaire, quadratique...) et saisonnailités (périodicités) de 4
chacune de ces séries. Un des objectifs principaux de l’étude d’une série temporelle est la prévision des réalisations futures, très souvent pour des raisons économiques (prévoir l’évolution de la vente d’un produit pour ajuster au mieux les moyens de production, prévoir l’évolution d’un marché financier ...). genre de Bien entendu, aucun modèle ne correspond exactement à la réalité, et il est impossible de prévoir parfaitement le devenir d’une série temporelle. Lorsque cela sera possible, nous donnerons des intervalles de prévisions, afin de pouvoir apporter une information quant à la précision de la prévision. Pour ce faire, il existe un large choix de modèle utilisable : – les modèles de régression, comme par exemple : xt = α1 t2 + α2 t + α3 + ǫt ,
t = 1, . . . , n.
Une fois les coefficients de ce modèle estimés, la prévision de xt+1 sera x ˆt+1 = αˆ1 (t + 1)2 + α ˆ 2 (t + 1) + α ˆ3 . – les lissages exponentiels qui sont très simples à mettre en oeuvre, et qui feront l’objet d’un chapitre suivant, – les modèles de type ARMA, qui consistent à enlever de la série les tendances et saisonnalités (ou prédiodicités) évidentes et à modéliser le résidu restant. Ces méthodes sont plus sophistiquées et plus lourdes numériquement (temps de calcul) que les précédentes, mais également plus performantes. Parmi les 5 exemples précédents, celui relatif au nombre de passagers dans les transports aériens (figure 2) est une série assez typique de ce que l’on rencontre en économétrie, et elle donne lieu à de bonnes prévisions pour toutes les méthodes classiques. Au contraire, l’évolution des marchés boursiers (figure 5) est beaucoup plus difficile à prévoir.
1.1 Tendances et composantes saisonnières On parle de tendance lorsque la série (xt )1≤t≤n peut s’écrire, à une erreur d’ajustement ǫt près, comme une combinaison linéaire de m fonctions, indépendantes du temps, choisies a priori (par exemple fonction puissance, exponentielle, logarithmique...) : xt =
m X
1 ≤ t ≤ n.
αj fj (t) + ǫt
j=1
Lorsque xt = αt + β + ǫt la tendance est linéaire (m = 1 et f (t) = αt + β). Une tendance polynomiale se traduira par xt = α1 tp + αp−1 tp−1 + . . . + αp+1 + ǫt . Dans l’exemple 5, la tendance semble être exponentielle.
On parle de composante périodique lorsque la série (xt )1≤t≤n peut se décomposer en : 1 ≤ t ≤ n,
xt = st + ǫt
où st est périodique, c’est-à-dire st+T = st , avec T la période (supposée entière). Lorsque la période est de 6 mois ou 1 an, on parle généralement de composante saisonnière.
Enfin, il est fréquent qu’une série comporte à la fois une tendance et une composante périodique (cf. exemple 2).
1.2 Indices descriptifs d’une série temporelle 1.2.1 Indices de tendances centrales Nous utilisons comme indicateur de la tendance centrale la moyenne : n
x¯n =
1X xt . n t=1 5
1.2.2 Indices de dispersion Nous utilisons comme indicateur de dispersion la variance empirique (et sa racine carrée, l’écart-type empirique) : n
σ ˆn (0) =
1X (xt − x ¯n )2 . n t=1
1.2.3 Indices de dépendance Ces notions, plus spécifiques à l’étude de série temporelle, renseignent sur la dépendance entre les données xt . Auto-covariance L’auto-covariance empirique d’ordre 1 renseigne sur la dépendance entre deux données successives : σ ˆn (1) =
n−1 1 X (xt − x ¯n )(xt+1 − x ¯n ), n − 1 t=1
l’auto-covariance empirique d’ordre 2 renseigne sur la dépendance entre deux données écartées de deux pas de temps : σ ˆn (2) =
n−2 1 X (xt − x ¯n )(xt+2 − x ¯n ), n − 2 t=1
et ainsi de suite. Pour des raisons de bon sens statistique, nous ne considèrerons les covariances empiriques que jusqu’à un ordre h pas trop grand. On appelle fonction d’auto-covariance (empirique) la fonction qui à j associe σ ˆn (j). Auto-corrélation Les auto-corrélations empiriques sont les quotients des covariances empiriques par la variance empirique : ρˆn (j) =
σ ˆn (j) σ ˆn (0)
j = 0, . . . , h.
Ce sont les auto-corrélations empiriques que nous utiliserons pour caractériser la dépendance entre les variables. On appelle fonction d’auto-corrélation (empirique) la fonction qui à j associe ρˆn (j).
La représentation graphique des nuages de points (xt , xt+1 ), pour t = 1, . . . , n − 1, est une bonne illustration de la valeur de l’auto-corrélation d’ordre 1 ρˆn (1) : plus le nuage sera arrondi, plus ρˆn (1) sera proche de 0, et plus le nuage sera allongé, plus ρˆn (1) sera proche de 1. Ceci est vrai pour toutes les auto-corrélations ρˆn (k), 1 ≤ k ≤ n. Proposition 1 : pour une série tendance linéaire pure xt = at + b, t = 1, . . . , n, on a : ρˆn (h) −−−−→ 1 n→∞
La preuve est un bon exercice. Proposition 2 : pour une série périodique pure xt = a cos 2tπ T , t = 1, . . . , n, on a : cos
2hπ
ρˆn (h) −−−−T−→ 1 n→∞
Ces deux propositions montrent comment repérer grâce à l’auto-corrélation qu’une série temporelle admet une tendance (l’auto-corrélation tend vers 1) ou une saisonnalité (la saisonnalité se voit sur l’auto-corrélation).
6
11000
11000
7000
9000
11000
11000
x[7:length(x)]
11000
11000
9000
x[1:(length(x) − 3)]
7000
7000
9000
x[6:length(x)]
11000 9000
9000
x[1:(length(x) − 2)]
7000
x[5:length(x)]
x[1:(length(x) − 1)]
7000
11000 7000
7000
9000
9000
9000
x[4:length(x)]
11000 7000
9000
x[3:length(x)]
11000 9000 7000
x[2:length(x)]
7000
7000
x[1:(length(x) − 4)]
9000
11000
7000
x[1:(length(x) − 5)]
9000
11000
x[1:(length(x) − 6)]
7000
9000
11000
x[1:(length(x) − 7)]
0.4
0.8 7000
−0.4 0.0
ACF
9000
x[9:length(x)]
9000 7000
x[8:length(x)]
Series x
7000
9000
11000
x[1:(length(x) − 8)]
0.0
0.5
1.0
1.5
Lag
F IG . 6 – Nuages de points (xt , xt+k ) pour k = 1, . . . , 8 et auto-corrélation pour la série temporelle du nombre de morts accidentelles aux Etats-Unis de 1973 à 1978
7
1.3 Mise en oeuvre sous R Quelques fonctions R utiles à l’étude des séries temporelles : – Lire un fichier de données en sautant les k premières lignes : data=scan(file=”donnee.dat”,skip=k). – Créer un objet de type série temporelle : serie <- ts (data,start,end,frequency). data contient le vecteur des données (un fichier contenant les données peut être mentionné en remplaçant data par file=”donnees.dat”), start et end mentionne les dates de début et de fin de la série, et enfin frequency mentionne le nombre de données par unité de temps (par exemple, si les dates de début et de fin sont des années, et que les données sont mensuelles, il faudra indiquer frequency=12). – Représenter graphiquement un objet de type série temporelle : plot.ts(serie) – La fonction acf(x, lag.max = 10, type = c("correlation", "covariance"), plot = TRUE) calcule (et trace si l’option plot est à TRUE) les lag.max premières auto-corrélations et auto-covariances.
8
2 Lissages exponentiels Les méthodes de lissages exponentiels constituent un outil permettant de réaliser des prévisions à partir de l’observation d’une série temporelle. Ces méthodes étant relativement basiques et simples de mise en oeuvre, elle sont souvent utilisées dans l’industrie, notamment lorsque le nombre de prévisions à réaliser est important (par exemple, prévisions des ventes de centaines de produits dans une grande surface). Nous présentons trois types de lissage exponentiel : – le lissage exponentiel simple qui consiste à ajuster localement à la série temporelle une constante, – le lissage exponentiel double qui ajuste quant à lui une droite, – le lissage exponentiel de Holt-Winters qui considère des fonctions plus complexes (polynomiales, périodiques...).
2.1 Lissage exponentiel simple Disposant d’une série temporelle x1 , . . . , xn , l’objectif du lissage exponentiel est d’estimer la valeur xn+h non encore observée. Nous noterons x ˆn,h cette prévision. Etant donné une constante de lissage 0 < α < 1, on définit la prévision par lissage exponentiel simple : x ˆn,h = (1 − α)
n−1 X
αj xn−j .
(1)
j=0
La prévision est une moyenne de toutes les observations passées, pondérée de sorte à ce que plus l’observation soit ancienne moins elle ait d’importance. Une constante de lissage α proche de 1 (≥ 0.7) donne une importance significative aux observations éloignées, tandis qu’un α proche de 0 (≤ 0.3) tend à négliger ces observations éloignées. Remarque : la prévision x ˆn,h ne dépend pas de h ! Formules récursives de mise à jour La définition (1) vérifiant la formule réccursive suivante x ˆn,h = (1 − α)xn + αˆ xn−1,h , la prévision xˆn,h peut être obtenue immédiatement à partir de la connaissance de : 1- la prévision x ˆn−1,h basée sur les n − 1-èmes premières observations, 2- l’observation xn . L’utilisation de cette récurrence permet de réaliser des algorithmes très rapides d’estimation de la prévision par lissage exponentiel (en initialisant à xˆ1,h = x1 ).
2.2 Lissage exponentiel double On ajuste au voisinage de l’instant n une droite d’équation yt = a1 + a2 (t − n). La prévision par lissage exponentiel double est : x ˆn,h = a ˆ1 + a ˆ2 h
avec
a ˆ1 = 2L1 (n) − L2 (n)
et
a ˆ2 =
1−α (L1 (n) − L2 (n)) α
Pn−1 Pn−1 où L1 (n) = (1 − α) j=0 αj xn−j et L2 (n) = (1 − α) j=0 αj L1 (n − j) sont deux lissages exponentiels simples successifs. Remarque : comme pour le lissage exponentiel simple, l’estimateur de la prévision est la meilleure approximation au sens des moindres carrés pondérés. Formules récursives de mise à jour a ˆ1 (n) = a ˆ2 (n) =
a ˆ1 (n − 1) + a ˆ2 (n − 1) + (1 − α2 )(xn − x ˆn−1,1 ), 2 a ˆ2 (n − 1) + (1 − α )(xn − x ˆn−1,1 ),
avec a ˆ1 (0) = x1 et a ˆ2 (0) = x2 − x1 comme valeurs initiales. 9
2.3 Méthode de Holt-Winters 2.3.1 Méthode non saisonnière Comme la méthode de lissage exponentiel doubles celle de Holt-Winters non saisonnière revient à estimer au voisinage de l’instant n une droite yt = a1 + a2 (t − n). La prévision prend la forme x ˆn,h = a ˆ1 (n) + hˆ a2 (n). La variante par rapport à la méthode de lissage exponentiel double est au niveau des formules de mise à jour dans l’estimation des paramètres a1 et a2 . Soit deux constantes de lissages 0 < β < 1 et 0 < γ < 1. Les formules de mise à jour sont : a ˆ1 (n) = a ˆ2 (n) =
βxn + (1 − β)[ˆ a1 (n − 1) + a ˆ2 (n − 1)], γ[ˆ a1 (n) − a ˆ1 (n − 1)] + (1 − γ)ˆ a2 (n − 1).
On remarquera que les formules de mise à jour du lissage exponentiel double sont un cas particulier de ces dernières, avec β = 1 − α2 et γ = 1−α 1+α . Remarque : l’introduction de deux constantes rend la méthode plus souple que le lissage exponentiel double. 2.3.2 Méthode saisonnière additive On cherche maintenant à ajuster au voisinage de l’instant n une droite d’équation yt = a1 + a2 (t − n) + st , où st est une composante périodique de période T . Les formules récursives de mise à jour sont : a ˆ1 (n) =
β(xn − sˆn−T ) + (1 − β)[ˆ a1 (n − 1) + a ˆ2 (n − 1)],
a ˆ2 (n) = sˆn =
γ[ˆ a1 (n) − a ˆ1 (n − 1)] + (1 − γ)ˆ a2 (n − 1), δ[xn − a ˆ1 (n)] + (1 − δ)ˆ sn−T .
Les prévisions sont de la forme : x ˆn,h x ˆn,h
= a ˆ1 (n) + hˆ a2 (n) + sˆn+h−T 1 ≤ h ≤ T, = a ˆ1 (n) + hˆ a2 (n) + sˆn+h−2T T + 1 ≤ h ≤ 2T.
et ainsi de suite pour h ≥ 2T . Se référer à Gouriéroux et Monfort 1983 [2] pour les valeurs d’initialisation. 2.3.3 Méthode saisonnière multiplicative On ajuste au voisinage de l’instant n une droite d’équation yt = [a1 + a2 (t − n)] × st , où st est une composante périodique de période T . Les formules récursives de mise à jour sont : a ˆ1 (n) = a ˆ2 (n) = sˆn
=
xn + (1 − β)[ˆ a1 (n − 1) + a ˆ2 (n − 1)], sˆn−T γ[ˆ a1 (n) − a ˆ1 (n − 1)] + (1 − γ)ˆ a2 (n − 1), xn + (1 − δ)ˆ sn−T . δ a ˆ1 (n) β
10
Les prévisions sont de la forme : xˆn,h xˆn,h
= [ˆ a1 (n) + hˆ a2 (n)]ˆ sn+h−T 1 ≤ h ≤ T, = [ˆ a1 (n) + hˆ a2 (n)]ˆ sn+h−2T T + 1 ≤ h ≤ 2T.
Se référer également à [2] pour les valeurs d’initialisation.
2.4 Mise en oeuvre sous R Les méthodes de lissages exponentiels sont disponibles sous R, grâce à la fonction HoltWinters. Pour une série temporelle x, cette procédure permet : – un lissage exponentiel simple : xlisse <- HoltWinters(x,alpha=α,beta=0,gamma=0), – un lissage de Holt-Winters sans composante saisonnière : xlisse <- HoltWinters(x,alpha=α,beta=β,gamma=0), – un lissage Holt-Winters additif : xlisse <- HoltWinters(x,alpha=α,beta=β,gamma=γ,seasonal=”add”), – un lissage Holt-Winters multiplicatif : xlisse <- HoltWinters(x,alpha=α,beta=β,gamma=γ,seasonal=”mul”). A noter que pour un lissage de Holt-Winters avec composante saisonnière la série temporelle x doit être un objet de type série temporelle, défini avec la fonction ts en précisant la saisonnalité. L’affichage et la visualisation des résultats peuvent être réalisés à l’aide des commandes : – summary(xlisse) : description de l’objet xlisse obtenu précédemment par la procédure HoltWinters, – plot(xlisse) : représentation des valeurs observées et des valeurs lissées, – plot(xlisse$fitted[,1]) : représentation des paramètres de la méthode de lissage exponentiel. Les prévisions sont réalisées à l’aide de la fonction predict : h<-5 ; p<-predict(xlisse,n.ahead=h). Un intervalle de confiance (dont le fondement théorique n’a pas été étudié dans ce cours) peut être obtenu en validant (à TRUE) l’option prediction.interval. Attention les paramètres du lissage exponentiel de ce cours (α, β, γ, δ) ne correspondent pas à leur homonyme implémenté dans la fonction HoltWinters de R. Voici un tableau présentant la correspondance des paramètres : méthode lissage exponentiel simple
paramètres du cours α
paramètres de R 1−α (β = γ = 0) Holt-Winters double β α γ β (γ = 0) Holt-Winters avec composante saisonnière β α γ β δ γ Ainsi, l’interprétation précédente du cours relativement au paramètre α du lissage exponentiel simple se trouve inversée pour le paramètre α de R.
11
F IG . 7 – Lissage et prévision par lissage exponentiel double d’un bruit blanc gaussien
12
F IG . 8 – Lissage et prévision par lissage exponentiel double de la série X(t) = 0.5t + 2ǫt avec ǫt ∼ N (0, 1)
13
F IG . 9 – Lissage et prévision par lissage exponentiel double de la série X(t) = 0.5t+ǫt +3 cos t π6 avec ǫt ∼ N (0, 1)
14
30 0
10
20
bb
40
50
Lissage Exponentiel Simple
2
4
6
8
10
8
10
8
10
Time
30 0
10
20
bb
40
50
Lissage Exponentiel Double, alpha=0.5
2
4
6 Time
30 0
10
20
bb
40
50
HoltWinters Seasonal
2
4
6 Time
F IG . 10 – Lissage et prévision par lissage exponentiel simple, double, et Holt-Winters avec composante saisonnière de la série X(t) = 0.5t + ǫt + 3 cos t π6 avec ǫt ∼ N (0, 1) 15
3 Estimation et élimination de la tendance et de la saisonnalité 3.1 Processus stationnaire On dit qu’une suite de variables aléatoires (X1 , . . . , Xn , . . .) est stationnaire si les espérances sont constantes E[Xn ] := µ
∀n,
et si les covariances sont stables par translation dans le temps, c’est-à-dire, pour tout h Cov(Xn , Xn+h ) := σ(h)
∀n.
Une telle suite est appelée processus stationnaire. On appelle fonction d’auto-covariance du processus stationnaire la suite σ(h), et fonction d’auto-corrélation du processus stationnaire la suite ρ(h) := σ(h) σ(0) . Propriété : σ(h) = σ(−h). Une série temporelle x1 , . . . , xn , comme nous l’avons présentée précédemment, est la réalisation des n premiers termes d’un processus stochastique (Xt )t . C’est ce processus que l’on cherche désormais à modéliser. Pour cela, la démarche suivante doit être adoptée : – représenter graphiquement la série afin de repérer les tendances et saisonnalités, – estimer et supprimer les tendances et saisonnalités (partie déterministe du processus stochastique), – choisir un modèle pour les résidus (partie aléatoire du processus stochastique) et l’estimer, – prédire les réalisations futures à l’aide de ce modèle. L’objectif de cette section est de donner quelques méthodes pour estimer et supprimer les tendances et saisonnalités. La fin de ce cours sera concentré sur la modélisation de processus stationnaires.
3.2 Une estimation paramétrique de la tendance Nous supposons que la série temporelle étudiée soit la réalisation d’un processus stochastique composé d’une tendance déterministe mt et d’une partie aléatoire ǫt (supposée de moyenne nulle) : X t = mt + ǫ t . Une méthode simple consiste à supposer que cette tendance est linéaire : mt = a + bt, et d’estimer les paramètres a et b par moindres carrés. Ainsi, si on observe la série x1 , . . . , xn , il faut trouver a et b qui minimisent la quantité : n X
(xt − a − bt)2
t=1
Les solutions de ce problème sont : a ˆ
=
ˆb =
! 2n + 1 − txt + n¯ x , 3 t=1 ! n X n+1 12 txt − n¯ x . n(n2 − 1) t=1 2 6 n(n − 1)
n X
Exercice : écrire le problème sous forme matricielle... et le résoudre. L’hypothèse de linéairité de la tendance convient très bien à certaines séries temporelles : par exemple, celles des figures 1, 2 (et encore...), 4... Mais ce n’est pas le cas de toutes les séries : voir par exemple celle représentant le cours 16
du CAC40, figure 5. Il est alors possible de supposer que la tendance soit de forme polynomiale : mt = a + bt + ct2 et d’estimer les paramètres a, b et c par moindres carrés. Mais il est parfois difficile d’estimer le degré du polynôme, et lorsque le degré est trop important, le nombre de paramètres à estimer devient grand et les calculs fastidieux. Dans cette situation, on a recourt à une méthode d’estimation non paramétrique.
3.3 Estimation non paramétrique : moyenne mobile Tendance Supposons que la tendance mt soit linéaire dans un petit intervalle [t − q, t + q] autour de t. Dans ce cas, un bon estimateur de la tendance est : q X 1 xt+k . m ˆt = 2q + 1 k=−q
On peut donc estimer la tendance à chaque temps t en calculant la moyenne sur les observations étant dans une fenêtre de largeur 2q + 1 autour de t : c’est ce que l’on appelle une estimation par moyenne mobile. Pour éviter les problèmes de bord, on suppose que xt = x1 si t < 1 et xt = xn si t > n. Tendance et saisonnalité Supposons désormais que le processus ne comporte par uniquement une tendance, mais également une saisonnalité : X t = m t + st + ǫ t , avec st une fonction T -périodique. Le principe d’estimation est (en simplifiant légèrement) le suivant : on estime la tendance moyenne sur une période, puis on estime la composante saisonnière en moyennant sur toutes les périodes les écarts à la tendance moyenne de la période. Remarque : il est également possible de considérer une composante saisonnière multiplicative. Application à la série du nombre de morts accidentelles aux Etats-Unis La figure 11 représente la série du nombre de morts accidentelles aux Etats-Unis de 1973 à 1978 (déjà présentée en introduction), la figure 12 représente la tendance (estimée par moyenne mobile) et la figure 13 représente la composante saisonnière. La figure 14 représente la série après élimination de la tendance et la figure 15 représente la série après élimination de la tendance et de la composante saisonnière. La série 15 ainsi obtenue est une série (supposée) stationnaire, sur laquelle nous chercherons plus tard à ajuster un modèle.
17
11000 10000 9000 7000
8000
USAccDeaths
1973
1974
1975
1976
1977
1978
1979
Time
8400 8600 8800 9000 9200 9400 9600
m$trend
F IG . 11 – Série USAccDeaths : nombre de morts accidentelles aux Etats-Unis de 1973 à 1978
1973
1974
1975
1976
1977
1978
1979
1978
1979
Time
0 −1500 −1000 −500
m$seasonal
500
1000 1500
F IG . 12 – Tendance de la série USAccDeaths
1973
1974
1975
1976
1977
Time
F IG . 13 – Composante saisonnière de la série USAccDeaths
18
2000 1000 0 −2000
−1000
USAccDeaths − m$trend
1973
1974
1975
1976
1977
1978
1979
Time
200 0 −200 −400
USAccDeaths − m$trend − m$seasonal
400
F IG . 14 – Série USAccDeaths après élimination de la tendance
1973
1974
1975
1976
1977
1978
1979
Time
F IG . 15 – Série USAccDeaths après élimination de la tendance et de la composante saisonnière
19
3.4 Elimination de la tendance et de la saisonnalité par la méthode des différences Cette méthode permet de supprimer les tendance et saisonnalité d’une série temporelle sans les estimer. Soit ∆T l’opérateur qui fait passer de (Xt ) à (Xt − Xt−T ) : ∆T Xt = Xt − Xt−T . On note ∆ l’opérateur ∆1 , et ∆kT l’opérateur ∆T ◦ . . . ◦ ∆T . | {z } k fois Supposons que le processus admette une tendance polynomiale d’ordre k : Xt =
k X
aj tj +ǫt .
j=0
| {z } mt
En appliquant ∆ à mt , on obtient un polynôme de degré k − 1. Ainsi, en appliquant k fois ∆, on élimine la tendance. Remarque : il est important de remarquer que si l’on applique ∆t quelque soit t, le résultat est le même quant à l’élimination de la tendance. Comme en pratique il n’est pas évident de connaître le degré k, on appliquera l’opérateur ∆ jusqu’à ce que la moyenne du processus soit nulle (k sera généralement 1, 2 ou 3). Supposons maintenant que le processus admette une tendance mt et une saisonnalité, de période T : X t = m t + st + ǫ t . Dans ce cas, ∆T Xt = (mt − mt−T ) + (ǫt − ǫt−T ) est un processus désaisonnalisé. De plus, si le processus admet une tendance linéaire, elle est également supprimé (cf. remarque précédente). Si la tendance est plus que linéaire, il suffit d’appliquer la procédure précédente pour finir de supprimer la tendance, et obtenir ainsi un processus que l’on supposera stationnaire. La figure 16 illustre l’élimination de la tendance linéaire et de la saisonnalité de la série xt = ǫt ∼ N (0, 1).
t 2
+ 3 cos tπ 6 ǫt avec
3.5 Mise en oeuvre sous R La fonction decompose permet d’extraire d’une série temporelle (via la méthode de la moyenne mobile) : serie_decomp<-decompose(serie,type=c(”additive”,”mutliplicative”)) – la composante saisonnière : serie_decomp$seasonal, que l’on suppose additive ou multiplicative dans l’option type, – la tendance : serie_decomp$trend, – le partie aléatoire stationnaire de la série : serie_decomp$random. La fonction diff.ts(serie,lag=T,difference=k) permet d’appliquer l’opérateur de différeciation ∆kT .
20
30 0
10
bb
50
Serie avec tendance lineaire et saisonnalité (période 12)
0
20
40
60
80
100
Time
2 0 −4
−2
bb_diff1
4
Serie avec la tendance eliminée par la méthode des différences
0
20
40
60
80
100
Time
6 4 2
bb_diff12
8
Serie avec la saisonnalité (et la tendance) eliminées par la méthode des différences
0
20
40
60
80
Time
F IG . 16 – Elimination de la tendances et saisonnalité par la méthode des différences (figure du haut : série xt = tπ t 2 + 3 cos 6 ǫt avec ǫt ∼ N (0, 1) ; figure du milieu : série xt − xt−1 ; figure du bas : série xt − xt−12
21
4 Modélisation des séries stationnaires Nous présentons dans cette section comment modéliser une série, qui une fois tendances et saisonnalités supprimées, est stationnaire. A noter que le seul fait de supprimer la tendance et la saisonnalité ne rend pas la série nécessairement stationnaire, puisque cela n’affecte pas la variance et l’auto-covariance, qui dovient être constantes pour un processus stationnaire.
4.1 Auto-corrélation partielle Le coefficient de corrélation partielle entre les deux variables X1 et Xn d’un processus stochatistique (Xt )t est le coefficient de corrélation entre les deux variables auxquelles on a retranché leurs meilleurs explications en terme de X2 , . . . , Xn−1 : rX2 ,...,Xn−1 (X1 , XN ) = corr(X1 − PX2 ,...,Xn−1 (X1 ), Xn − PX2 ,...,Xn−1 (Xn )), où corr est le coefficient de corrélation classique (quotient de la covariance par le produit des écarts-types), et où PX2 ,...,Xn−1 (X1 ) est la projection1 de la variable X1 dans l’espace vectoriel engendré par les variables X2 , . . . , Xn−1 . Ce coefficient exprime la dépendance entre les variables X1 et Xn qui n’est pas due aux autres variables X2 , . . . , Xn−1 . La fonction d’auto-corrélation partielle r(h) d’un processus stationnaire est défini par : r(h) r(h)
= rX2 ,...,Xh (X1 , Xh+1 ) = r(−h) ∀h 6= 0
r(1)
= ρ(1)
∀h ≥ 2
L’algorithme de Durbin-Watson, que nous ne présentons pas ici, permet d’estimer les auto-corrélations partielles d’un processus stationnaire. Dans le logiciel R, la fonction pacf permet ces estimations.
4.2 Quelques exemples simples de processus stationnaires. 4.2.1 Bruit blanc Un bruit blanc est une suite de variables aléatoires indépendantes, d’espérance et de variance constantes. Si l’espérance est nulle, le bruit blanc est centré. Ces processus ne peuvent être employés en prévision (puisque chaque réalisation du processus est indépendante des réalisations passées), mais sert de base à la construction de beaucoup d’autre processus. Exercice : quelle est la fonction d’auto-covariance d’un tel processus ? 4.2.2 Moyenne mobile On appelle moyenne mobile (Moving Average) un processus de la forme Xt = ǫt + a1 ǫt−1 + . . . + aq ǫt−q , où les ǫj pour t − q ≤ j ≤ t sont des bruits blancs centrés de variance σ 2 , et où l’entier non nul q est l’ordre de la moyenne mobile. Ces processus, que l’on note M Aq , seront étudiés plus en détail plus loi. Exercice : quelle est la fonction d’auto-covariance d’un tel processus ? 1 la projection P X2 ,...,Xn−1 (X1 ) = γ2 X2 + . . . + γn−1 Xn−1 est caractérisée par le fait qu’elle est orthogonale à toutes les variables X2 , . . . , Xn−1 .
22
4.2.3 Processus auto-régressif d’ordre 1 : AR1 Un processus auto-régressif d’ordre 1 est un processus stationnaire centré qui vérifie les équations de récurrence : Xt = aXt−1 + ǫt , où |a| < 1 et ǫt un bruit blanc centré de variance σ 2 . Propriétés : E[Xt ] = σ(h)
0 σ2
=
ah 1 − a2
∀h ≥ 0
Ainsi, l’auto-covariance (et l’auto-corrélation) tendent vers 0 lorsque h tend vers l’infini.
4.3 Les processus linéaires stationnaires Un processus discret linéaire (Xt ) est obtenu à partir d’une suite d’aléas indépendants ǫt (chocs aléatoires) sous la forme d’une combinaison linéaire infinie des chocs passés : Xt =
∞ X
t∈Z
bj ǫt−j
(2)
j=0
Remarque : en comparant avec la définition précédente d’un processus moyenne mobile, vous comprendrez pourquoi on parle parfois de processus moyenne mobile infinie pour un processus linéaire. P Pour que cette somme ait un sens, nous supposons que la série des bj est absolument convergente : ∞ j=0 |bj | < ∞ et que b0 = 1. De plus, dans le but de faire de la prévision, nous supposons que ce processus est stationnaire. Ainsi, les ǫt ont même espérance et même variance. Le processus (ǫt ) est un bruit blanc. Nous supposons désormais, sans perte de généralité, que le processus est centré, et donc que les ǫt sont d’espérance nulle. Nous notons σ 2 la variance des ǫt . Interprétation des ǫt
L’équation (2) peut également d’écrire : ǫt = Xt −
∞ X
bj ǫt−j
j=1
et, si on remplace les ǫt−j par leur valeurs, on obtient ǫt = Xt −
∞ X
aj Xt−j
j=1
sous réserve que les coefficients aj (qui sont fonction des bj ) vérifient processus est inversible). Le modèle linéaire stationnaire discret peut ainsi s’écrire : Xt = ǫt +
∞ X
aj Xt−j .
P∞
j=0
|aj | < ∞ (auquel cas nous dirons que le
(3)
j=1
Xt est alors la somme P∞ d’un choc aléatoire à l’instant t (ǫt ), indépendant de l’historique, et d’une fonction linéaire de son passé xˆt = j=1 aj Xt−j , qui peut être vue comme la prédiction de Xt à partir des observations passées. ǫt est l’innovation contenue dans le processus au temps t. On dit que (ǫt ) est le processus d’innovation. 23
Remarque : les expressions (2) et (3) d’un même processus linéaire stationnaire sont appelées expression en somme infinie et expression auto-régressive. Le nombre de paramètres d’un processus linéaire stationnaire étant infini, nous considérons les modèles tronqués suivants : – le processus moyenne mobile d’ordre q, noté M Aq : Xt =
q X
bj ǫt−j
(4)
j=0
– le processus auto-régressif d’ordre p, noté ARp : Xt = ǫt +
p X
aj Xt−j
(5)
j=1
Remarque : t ∈ Z. 4.3.1 Les processus moyenne mobile On peut montrer que le processus moyenne mobile d’ordre q a pour auto-covariance : Pq−h σ 2 k=0 bk bk+h ∀h ≤ q σ(h) = 0 ∀h > q Ainsi, nous reconnaissons un processus moyenne mobile et identifions son ordre lorsque son auto-covariance (et par conséquent son auto-corrélation) est nulle à partir d’un certain ordre. Quant à l’auto-corrélation partielle, elle tend vers 0 lorsque h grandit : lim r(h) = 0
h→∞
4.3.2 Les processus auto-régressif Le processus AR1 Xt = ǫt + aXt−1 On a vu précédemment que l’auto-covariance d’un AR1 est : σ(h) = σ 2
ah . 1 − a2
On en déduit l’auto-corrélation : ρ(h) = ah . En remplacant dans l’expresion Xt = ǫt + aXt−1 d’un AR1 , Xt−1 par son expression, et ainsi de suite, on obtient l’expression en somme infinie d’un AR1 : Xt =
∞ X
aj ǫt−j .
j=0
On déduit alors de l’existance de ce processus la condition sur a suivante : |a| < 1. Ainsi l’auto-covariance et l’auto-corrélation d’un AR1 tendent vers 0 lorsque h tend vers l’infini.
24
Les processus ARp
Xt = ǫt +
p X
aj Xt−j avec ap 6= 0.
j=1
La variance d’un tel processus est :
σ(0) = σ 2 +
p X
aj σ(j)
j=1
et l’auto-covariance, pour tout h > 0 : σ(h) =
p X
aj σ(h − j).
j=1
On peut montrer que l’auto-covariance et l’auto-corrélation tendent vers 0 lors h tend vers l’infini. Quant à l’auto-corrélation partielle, elle vérifie : ∀h > p,
r(h)
=
0
r(p)
=
ap .
4.3.3 Les processus mixte ARM Ap,q Le nombre de paramètre d’un processus moyenne mobile ou d’un auto-régressif peut parfois être réduit en considérant un processus mixte : auto-régressif moyenne mobile d’ordres p et q : Xt −
p X
ak Xt−k = ǫt −
q X
bj ǫt−j .
j=1
k=1
La stationnarité d’un ARM Ap,q est assurée lorsque toutes les racines du polynôme A(z) = 1 + a1 z + . . . + ap z p sont de module strictement supérieur à 1. Ce polynôme forme avec B(z) = 1 + b1 z + . . . + bq z p les deux polynômes caractéristiques du processus. On supposera également que les polynômes A et B n’ont pas de racine commune, afin de s’assurer qu’il n’y a pas de représentation plus courte. Exemple : soit l’ARM A1,1 défini par Xn + aXn−1 = ǫn + aǫn−1 . Quelles sont les racines des deux polynômes caractéristiques ? Existe-t-il une écriture plus simple de ce processus ? La bonne qualité de modélisation des processus ARMA est due à la propriété suivante : Proposition : pour une large classe de suite d’auto-covariances (σh )h et pour tout entier K positif, il existe un processus ARMA tel que ses K premières auto-covariances coïncident avec σ1 ≤ . . . ≤ σK . Exercice : quelle est la fonction d’auto-covariance d’une ARM A1,1 ? La talbeau 4.3.3 récapitule les principales propriétés des processus M Aq , ARp et ARM Ap,q . modèle auto-covariance
M Aq σ(h) = 0 ∀h > q
auto-corrélation
ρ(h) = 0 ∀h > q
auto-corrélation partielle
lim r(h) = 0
h→∞
ARp lim σ(h) = 0
ARM Ap,q ∀h > q, lim σ(h) = 0
lim ρ(h) = 0
∀h > q, lim ρ(h) = 0
h→∞
h→∞
h→∞
h→∞
r(h) = 0 ∀h > p et r(p) = ap
TAB . 1 – Récapitulatif des propriétés des processus M Aq , ARp et ARM Ap,q .
4.3.4 Estimation, prévision Cette partie n’est pas détaillée dans ce cours. Nous donnons juste quelques pistes dans le cadre d’un modèle mixte ARM Ap,q . 25
Estimation et choix de modèle Exemple (en cours) : cas d’un AR2 . Dans le cas général, l’estimation des paramètres des modèles ARMA peut être faite par maximum de vraisemblance. L’expression de la vraisemblance étant généralement trop complexe pour que l’on puisse obtenir un maximum explicite, des algorithmes numériques sont utilisés. Afin de choisir entre plusieurs modèles possibles (notamment pour le choix des ordres p et q d’un modèle ARMA), des critères de choix de modèles (AIC, BIC, ...) peuvent être employés. Plus p et q seront grands, plus le nombre de paramètres du modèle ARMA sera important, ce qui permettera de s’ajuster au mieux aux données. Mais en contrepartie, l’estimation de ces paramètres demandera un plus grand nombre de données. Pour un historique de données fixé, plus le nombre de paramètres sera grand moins bonne sera l’estimation. Les critères de type BIC et AIC ont pour but de déterminer le modèle qui effectue le meilleur compromis entre pouvoir prédictif et complexité. Prévision L’objectif est de prévoir la valeur que va prendre la variable aléatoire Xn+h , h étant appelé l’horizon de la prévision, ayant observé la réalisation des variables aléatoire X1 , . . . , Xn . Dans le cadre d’une modélisation ARMA, on choisit d’estimer Xn+h par une combinaison linéaire des Xj précédents (1 ≤ j ≤ n) : ˆ n,h = c1,h X1 + . . . + cn,h Xn . X Les coefficients cj,h sont estimés de sorte à ce qu’il minimise : E[(Xn+h − c1,h X1 − . . . − cn,h Xn )2 ]. ˆ n,h ainsi défini n’est autre que la projection de Xn+h sur l’espace vectoriel engendré par les variables L’estimateur X X 1 , . . . , Xn . Exemple (en cours) : on montre que l’erreur de prévision au rang 1 pour un ARp n’est autre que le bruit d’innovation. Le calcul de la projection dans le cas général n’est pas détaillé dans ce cours, mais voici néanmoins un résultat important sur le comportement de la prévision : la variance de l’erreur de prevision à horizon h croît depuis la variance du bruit d’innovation jusqu’à la variance du processus lui-même. Intervalle de confiance Si on suppose que le bruit d’innovation ǫt est gaussien, les variables aléatoires Xt sont elles aussi gaussiennes, tout comme l’erreur de prédiction. Ainsi, il sera possible de construire des intervalles de confiances sur la prédiction, ce que permettent les fonctions du logiciel R décrites plus loin. 4.3.5 Illustrations Voici quelques simulations de processus ainsi que leurs auto-corrélation et auto-corrélation partielle empiriques : – AR1 avec coefficient positif puis négatif : figures 17 et 18, – AR2 : figures 19 et 20, – M A1 avec coefficient positif puis négatif : figures 21 et 22, – M A3 : figure 23, – ARM A2,2 : figure 24.
26
4 2 0 −4
−2
ar1
0
200
400
600
800
1000
Time
0.4 0.0
ACF
0.8
Series ar1
0
5
10
15
20
25
30
25
30
Lag
0.4 0.0
Partial ACF
0.8
Series ar1
0
5
10
15
20
Lag
F IG . 17 – Simulation d’un AR1 : Xt = 0.8Xt−1 + ǫt , auto-corrélation et auto-corrélation partielle.
27
4 2 0 −4
ar1
0
200
400
600
800
1000
Time
0.0 −0.5
ACF
0.5
1.0
Series ar1
0
5
10
15
20
25
30
25
30
Lag
−0.4 −0.8
Partial ACF
0.0
Series ar1
0
5
10
15
20
Lag
F IG . 18 – Simulation d’un AR1 : Xt = −0.8Xt−1 + ǫt , auto-corrélation et auto-corrélation partielle.
28
6 4 2 −2 −6
ar2
0
200
400
600
800
1000
Time
0.6 0.2 −0.2
ACF
1.0
Series ar2
0
5
10
15
20
25
30
25
30
Lag
0.6 0.2 −0.2
Partial ACF
Series ar2
0
5
10
15
20
Lag
F IG . 19 – Simulation d’un AR2 : Xt = 0.9Xt−2 + ǫt , auto-corrélation et auto-corrélation partielle.
29
5 0 −5
ar2
0
200
400
600
800
1000
Time
0.0 −0.5
ACF
0.5
1.0
Series ar2
0
5
10
15
20
25
30
25
30
Lag
0.0 −0.4 −0.8
Partial ACF
Series ar2
0
5
10
15
20
Lag
F IG . 20 – Simulation d’un AR2 : Xt = −0.5Xt−1 − 0.9Xt−2 + ǫt , auto-corrélation et auto-corrélation partielle.
30
4 2 0 −4
−2
ma1
0
200
400
600
800
1000
Time
0.5 0.0 −0.5
ACF
1.0
Series ma1
0
5
10
15
20
25
30
25
30
Lag
−0.1 −0.3 −0.5
Partial ACF
Series ma1
0
5
10
15
20
Lag
F IG . 21 – Simulation d’un M A1 : Xt = ǫt − 0.8ǫt−1 , auto-corrélation et auto-corrélation partielle.
31
4 2 0 −4
−2
ma1
0
200
400
600
800
1000
Time
0.4 0.0
ACF
0.8
Series ma1
0
5
10
15
20
25
30
25
30
Lag
0.2 −0.2
Partial ACF
0.4
Series ma1
0
5
10
15
20
Lag
F IG . 22 – Simulation d’un M A1 : Xt = ǫt + 0.8ǫt−1 , auto-corrélation et auto-corrélation partielle.
32
10 5 0 −10
ma3
0
200
400
600
800
1000
Time
0.4 0.0
ACF
0.8
Series ma3
0
5
10
15
20
25
30
25
30
Lag
0.4 0.2 −0.2
Partial ACF
Series ma3
0
5
10
15
20
Lag
F IG . 23 – Simulation d’un M A3 , auto-corrélation et auto-corrélation partielle.
33
5 0 −5
arma22
0
200
400
600
800
1000
Time
0.0 −0.5
ACF
0.5
1.0
Series arma22
0
5
10
15
20
25
30
25
30
Lag
0.0 −0.4 −0.8
Partial ACF
0.4
Series arma22
0
5
10
15
20
Lag
F IG . 24 – Simulation d’un ARM A2,2 , auto-corrélation et auto-corrélation partielle.
34
5 Processus non stationnaire : ARIMA et SARIMA Les modèles ARMA sont destinés à modéliser des processus stationnaires. En pratique, les séries temporelles sont généralement non stationnaires, et un pré-traitement est nécessaire pour supprimer les tendances et saisonnalités. Une fois la série « stationnarisée »anlaysée, et les valeurs futures prédites, il est ensuite nécessaire de revenir à la série initiale. Exercice : écrire le retour à la série initiale dans le cas de la modélisation d’une série différenciée par ∆. Les processus ARIMA et SARIMA sont la généralisation des modèles ARMA pour des processus non stationnaire, admettant une tendance (ARIMA), ou encore une tendance et une saisonnalité (SARIMA). Soit ∆ l’opération de différenciation précédemment défini (section 3), qui associe au processus (Xt )t∈N le processus (Xt − Xt−1 )t∈N ). Nous rappelons que le processus obtenu en différenciant deux fois de suite, (Xt − 2Xt−1 + Xt−2 )t∈N , est noté ∆2 . Et ainsi de suite. Un processus (Xt ) est un processus ARIM Ap,d,q si le processus Yt = ∆d Xt est un processus ARM Ap,q . Les processus ARIM Ap,d,q sont bien adaptés aux séries temporelles présentant une tendance polynômiale de degré d − 1. Soit ∆T l’opérateur précédemment défini (section 3), qui fait passer de (Xt ) à (Xt − Xt−T ). Un processus (Xt ) est un processus SARIM Ap,d,q,T si le processus Yt = ∆T ◦ ∆d Xt est un processus ARM Ap,q . Les processus SARIM Ap,d,q,T sont bien adaptés aux séries temporelles présentant une période de longueur T et une tendance polynômiale de degré d. Remarque : en réalité, les processus SARIMA sont plus complexes et comportent d’autres paramètres (pour permettre des périodes de longueurs différentes notamment). Nous ne considèrerons que cette version SARIM Ap,d,q,T dans ce cours, se référer à [2] pour des compléments sur ces processus.
35
5.1 Mise en oeuvre sous R : processus ARMA, ARIMA, SARIMA La fonction arima.sim(modele,n) permet de simuler un processus ARM Ap,q défini par Xt −
p X
ak Xt−k = ǫt −
q X
bj ǫt−j .
j=1
k=1
Les paramètres ak et bj du processus sont précisés dans le paramètre modele de la fonction : modele<-list(ar=c(a1, . . . , ap ),ma=c(b1 , . . . , bq )). Pour simuler un modèle ARIM Ap,d,q il faut ajouter le composant order=c(p,d,q) dans le paramètre modele de la fonction arima.sim. La fonction ar permet d’estimer les paramètres d’un processus ARp : out<-ar(data,aic=TRUE,order.max=NULL) L’ordre p du processus auto-régressif est choisi (inférieur à order.max) à l’aide du critère AIC (si l’option aic est validée). La fonction arima permet d’estimer les paramètres : – d’un ARM Ap,q : out<-arima(serie,order=c(p,0,q)) – d’un ARIM Ap,d,q : out<-arima(serie,order=c(p,d,q)) – d’un SARIM Ap,d,q,T : out<-arima(serie,order=c(p,d,q),seasonal=list(order=c(P,D,Q),period=T) Les paramètres P, D, Q du modèle SARIMA ne sont pas abordés dans ce cours, nous leur donnerons par défaut la valeur des paramètres p, d, q. Parmi les sorties de cette fonction, on peut obtenir : – out$coef : estimation des coefficients, – out$aic : valeur du critère AIC, – out$resid : estimation des résidus. La fonction p=predict(out,h) permet d’effectuer une prévision à l’horizon h. Parmi les sorties de cette fonction, p$pred contient les prévisions, et p$se contient l’écart-type de l’erreur de prévision. Il n’existe pas de fonction prédéfinie pour calculer un intervalle de confiance sur les prévisions, mais cela peut être fait manuellement grâce à ces deux sorties de la fonction predict.
36
6 Processus ARCH et GARCH
−0.05 −0.15
−0.10
x
0.00
0.05
Les séries précédemment étudiées étaient supposées stationnaires. Si besoin, tendances et saisonnalités étaient supprimées pour obtenir une série résiduelle stationnaire. Néanmoins, toutes les séries résiduelles obtenues de la sorte ne sont pas nécessairement stationnaires. C’est le cas par exemple de la série représentée par la figure 25, qui contient les évolutions journalières de la bourse des valeurs de New-York (NYSE) du 19 octobre 1984 au 31 décembre 1991.
0
500
1000
1500
2000
Index
F IG . 25 – Evolution journalière de la bourse des valeurs de New-York (1984-1991)
−2
−1
0
x
1
2
3
La figure 26 représente la simulation d’un processus ARCH2 .
0
200
400
600
800
1000
Time
F IG . 26 – Simulation d’un ARCH2 Comme on peut le voir, la moyenne semble constante alors que la variance change au cours du temps (on qualifie ce comportement d’« hétéroscédastique »). De plus, les moments de grande variabilité semblent regroupés. Les modèles de type ARIMA qui supposent un comportement « homoscédastique »(variance constante), ne sont pas adaptés à ce type de série. Nous présentons dans cette section des modèles adaptés à ce type de série : les processus ARCH (AutoRegressive Conditionally Heteroscedastic) introduits pas Engle vers 1982, ainsi que leur généralisation, les processus GARCH.
37
6.1 Définitions Un processus ARCHp est défini par : Xt = ǫt |Xt−1 , Xt−2 , . . . ∼ σt2
=
ǫt N (0, σt2 ) 2 2 α0 + α1 Xt−1 + . . . + αp Xt−p
La variance condionnelle dépend du temps : si les valeurs précédentes sont grandes (en valeur absolue), la variance sera grande, et inversement. Ainsi, si on observe un choc dans la série (valeur anormalement grande), elle sera suivi d’une période de haute volatilité, dont la durée dépend de l’ordre p du modèle ARCH. Exercice : cela correspond-t-il à la série NYSE (figure 25) ?
6.2 Rappels de probabilités Soit un couple de variables aléatoires continues (X, Y ), de dentsité f (., .). Soit fY (.) la densité marginale de Y . La densité conditionnelle de X sachant que Y = y est définie par : fX|Y (x, y) =
f (x, y) fY (y)
si fY (y) > 0
L’espérance conditionelle est donc l’espérance par rapport à cette loi : Z E[X|Y = y] = xfX|Y (x, y)dx R
La variance conditionnelle est 2
σt2 = V(Xt |It−1 ) = E[X 2 |Y = y] − E[X|Y = y]
Soit Is−1 = {Xs−1 , Xs−2 , . . .} l’ensemble des observations précédant l’instant s. On a : E[Xt |Is ] = E[Xt |Ir ] =
Xt ∀t ≤ s E[E[Xt |Is ]|Ir ]
∀r ≤ s
et en particulier E[Xt ] = E[E[Xt |Is ]].
6.3 Propriétés des processus ARCH Si Xt est un processus ARCH, alors : E[Xt ] =
0
E[Xt |It−1 ] =
0
V(Xt ) = V(Xt |It−1 ) = Cov(Xt , Xt+h ) = σh = Cov(Xt , Xt+h |It−1 ) =
1−
α P0p
i=1
αi
si
p X
αi < 1
i=1
2 2 α0 + α1 Xt−1 + . . . + αp Xt−p 0 ∀h > 0
0
Remarque : un processus ARCH est conditionnellement hétéroscédastique mais inconditionnellement homoscédastique ! Pp Condition suffisante de stationnarité : i=1 αi < 1. On peut montrer également, ce qui peut être intéressant pour détecter un ARCH en pratique, que la distribution d’un processus ARCH a – un skewness nul (moment centré d’ordre 3) : la distribution est donc symétrique, – un kurtosis (moment centré d’ordre 4) supérieur à 3 : la distribution est donc plus applatie qu’une gaussienne. 38
6.4 Processus GARCH et propriétés Un processus GARCHp,q est défini par : Xt
=
ǫt |Xt−1 , Xt−2 , . . . ∼ σt2
=
ǫt N (0, σt2 ) 2 2 2 2 α0 + α1 Xt−1 + . . . + αp Xt−p + β1 σt−1 + . . . + βq σt−q
avec α0 > 0, αi ≥ 0 pour i = 1, . . . , p et βj ≥ 0 pour j = 1, . . . , q. Un processus GARCH peut être vu comme un processus ARCH d’ordre infini, et peut ainsi représenter formellement de façon plus parcimonieuse un processus ARCH comprenant un nombre élevé de paramètres. Remarque évidente : un GARCHp,0 est un ARCHp . Si Xt est un processus GARCHp,q , alors : E[Xt ] =
0
E[Xt |It−1 ] = Cov(Xt , Xt+h ) = σh =
0 0
Cov(Xt , Xt+h |It−1 ) =
0
∀h > 0
Propriété : soit Xt un processus GARCHp,q , et soit m = max(p, q). Le processus Xt2 admet une représentation ARM A(m, q). Ainsi, pour identifier un GARCHp,q , on identifiera tout d’abord le processus ARM A(m, q) qui modélise Xt2 . Pour identifier p dans le cas où m = q (p ≤ q), il faut effectuer des tests de significativité des coefficients αq , . . . , α1 du processus ARM A(m, q).
6.5 Mise en oeuvre sous R Pour utiliser les fonctions spécifiques à l’étude des modèles ARCH et GRACH, il faut avant tout charger le package tseries à l’aide de la commande library(tseries). La fonction garch permet d’estimer un GARCHp,q : serie<-garch(data,order=c(q,p)) Parmi les sorties de cette fonction : coef, residuals, fitted.values.
Références [1] Bosq D. et Lecoutre J-P. Analyse et prévision des séries chronologiques, Masson, Paris, 1992. [2] Gouriéroux C. et Montfort A. Cours de séries temporelles, Economica, Paris, 1983. [3] Shumway R.H. et Stoffer D.S. Times Series Analysis and Its Applications, With R Example, Springer, 2006. Ce cours est en partie inspiré du cours de M.-Cl. Viano disponible à l’adresse suivante : http ://math.univ-lille1.fr/ viano/economcours.pdf
39
