Introduction

• l’objectif principal de l’analyse d’une série temporelle est la prévision de ses futures réalisations en se basant sur ses valeurs passées

• Une série temporelle \(Y_t\) est communément décomposée en tendance, saisonnalité, bruit:

\[ Y_t = T_t+S_t+\varepsilon_t \]

• la tendance \(T_t\) correspondant à une évolution à long terme de la série, par exemple:

  • tendance linéaire: \(T_t=a+bt\)
  • tendance quadratique: \(T_t=a+bt+ct^2\)
  • tendance logarithmique: \(T_t=\log(t)\)

• la saisonnalité \(S_t\) correspondant à un phénoméne périodique de période identifiée

• l’erreur \(\varepsilon_t\) qui est la partie aléatoire de la série

Le but cette décomposition est de se ramener à un bruit \(\varepsilon_t\) stationnaire. Cette décomposition pouvant être additive, multiplicative \(Y_t = T_t*S_t*\varepsilon_t\) ou des combinaisons des deux:

\[ Y_t = (T_t+S_t)*\varepsilon_t \]

\[ Y_t = (T_t*S_t)+\varepsilon_t... \]

voir Hyndman et al. (2008) Forecasting with exponential smoothing: the state space approach, Springer-Verlag. http://www.exponentialsmoothing.net.

• nous nous intéressons ici aux méthodes de base (modélisation de série chro. linéaires): lissages exponentiels, modéles de régression (régression linéaire, modéles non-paramétriques. . . ), modéles SARIMA

Stationnarité

définition soit un processus aléatoire \((Y_t)_{t \in \textbf{Z}}\), il est dit stationnaire au sens fort (ou strictement) si pour toute fonction f mesurable \(f(Y_1,Y_2,...,Y_t)\) et \(f(Y_{1+h}, Y_{2+h}, ..., Y_{t+h})\) ont la même loi.

Cette notion de stationnarité forte est très difficile à vérifier en pratique. On lui préfère généralement la notion de stationnarité faible qui porte sur les moments d’ordre 1 et 2 du processus.

définition la fonction d’auto-covariance d’un processus \(Y_{t \in \textbf{Z}}\)

\[ \text{cov}(Y_t, Y_{t+h})= \gamma(h) \]

définition la fonction d’auto-corrélation d’un processus \(Y_{t \in \textbf{Z}}\)

\[ \rho(h)= \gamma(h)/\gamma(0) \]

\(\gamma(h)\) et \(\rho(h)\) sont des fonctions symétriques, \(\rho(0)=1\).

définition soit un processus aléatoire \((Y_t)_{t \in \textbf{Z}}\) tel que \(E(Y_t^2) < \infty\), il est dit stationnaire au sens faible (ou d’ordre 2) si son espérance est constante et ses auto-covariances sont stables dans le temps ie:

\[ \forall t \quad E(Y_t)= \mu \]

\[ \forall t \quad, \forall h \quad \text{cov}(Y_t, Y_{t+h})= \gamma(h) \]

On remarque que \(\text{var}(Y_t)=\gamma(0)\) et donc qu’un processus stationnaire faible à une variance constante dans le temps.

En pratique, pour apprécier la stationnarité d’un processus, on commence d’abord par vérifier que sa moyenne et sa variance sont constantes dans le temps.

exercice selon vous quel(s) processus ci-dessous est(sont) stationnaire(s)? Pourquoi?

Voilà quelques examples de processus stationnaires:

• un bruit blanc \(\varepsilon_t\) vérifiant \(E(\varepsilon_t)=\mu\) et \(\text{var}(\varepsilon_t)=\sigma^2\)

preuve on a par définition \(\text{cov}(\varepsilon_t,\varepsilon_{t+h})=0\)

• Le processus gaussien \((Y_t)_{t \in \textbf{Z}}\) tel que \(E(Y_t)=\mu\) et \(\text{cov}(Y_t,Y_{t+h})=\alpha^{|h|}\) (\(|\alpha|<1\)) est faiblement stationnaire. Tout processus gaussien stationnaire faible est stationnaire fort.

• le processus moyenne mobile \(X_t=\varepsilon_t+a_1\varepsilon_{t-1}+a_2\varepsilon_{t-2}+...+a_q\varepsilon_{t-q}\)

preuve

• processus autorégressif d’ordre 1:

\[ Y_t= a Y_{t-1}+ \varepsilon_t \]

en supposant que \(|a|<1\) on a bien \(E(Y_t)=0\) et

\[ \text{var}(Y_t)= \sigma^2(1+a+a^2+...)= \frac{\sigma^2}{1-a^2} \]

pour tout \(h>0\):

\[ \gamma(h)= \sigma^2(a^h+a^{h+2}+...)= \frac{\sigma^2 a^h}{1-a^2} \]

comme de plus \(\gamma(h)=\gamma(-h)\),

\[ \gamma(h)= \frac{\sigma^2 a^{|h|}}{1-a^2} \]

on remarque que pour ce processus \(\rho(h)= a^{|h|}\), donc l’autocorrélation tend vers 0 à une vitesse exponentielle.

En pratique, on ne connait pas explicitement les fonctions d’auto-covariance et d’auto-corrélation. Il est donc nécessaire de les estimer en se basant sur des observations.

définition soit une série d’observations \((y_t)_{t \in (1,...,n)}\), notons \(\bar{y}=\frac{1}{n} \sum_{t=1}^n y_t\), alors la fonction d’auto-covariance empirique vaut, pour tout \(h \in (0,...,n-1)\)

\[ \widehat{\gamma}(h)= \frac{1}{n-h} \sum_{t=h+1}^n (y_t-\bar{y})(y_{t-h}-\bar{y}) \]

définition soit une série d’observations \((y_t)_{t \in (1,...,n)}\), notons \(\bar{y}=\frac{1}{n} \sum_{t=1}^n y_t\), alors la fonction d’auto-corrélation empirique vaut, pour tout \(h \in (0,...,n-1)\)

\[ \widehat{\rho}(h)= \frac{\frac{1}{n-h} \sum_{k=h+1}^n (y_t-\bar{y})(y_{t-h}-\bar{y})}{\frac{1}{n}\sum_{t=1}^n (y_t-\bar{y})^2} \]

Le graphique représentant la fonction d’auto-corrélation empirique est appelé l’auto-corrélogramme.

exemple voilà un exemple de série et son auto-corrélogramme, à votre avis de quel type de série s’agit-il?

Auto-corrélation partielle

Lorsque l’on s’intéresse à caractériser les dépendances d’au moins 3 variables aléatoires, il est nécessaire d’introduire la notion de corrélation partielle. En effet, si l’on considère les variables \(X_1,...,X_k\), \(X_1\) peut être corrélée à \(X_3\) parce que \(X_1\) et \(X_3\) sont toutes deux corrélées à \(X_2\).

définition soit les variables aléatoires \(X_1,...,X_k\), le coefficient de corrélation partielle entre \(X_1\) et \(X_k\) conditionnellement à \(X_2,...,X_{k-1}\) est définie par:

\[ r_{X_2,..,X_{k-1}}(X_1,X_k)= \rho(X_1-P_{M(X_2,..,X_{k-1})}, X_k-P_{M(X_2,..,X_{k-1})}) \]

projection \(P_{M(X_1,..,X_k)}(Y)=X\alpha\) la projection linéaire d’une variable \(Y\) sur \(M(X_1,..,X_k)\) est telle que

\[ E[(Y-X\alpha)^2]\leq E[(Y-X\beta)^2], \forall \beta \]

c’est le vecteur de \(M(X_1,..,X_k)\) le plus proche de \(Y\) au sens de la distance \(d\) définie ci-dessous.

distance on définie la distance entre \(X\) et \(Y\) deux variables aléatoires par \(d(X,Y)=\sqrt{E[(X-Y)^2]}\)

espace engendré l’espace engendré par les variables aléatoires \(X_1, X_2, ...X_k\), noté \(M(X_1,..,X_k)\) est l’ensemble des combinaisons linéaires de ces variables \(M(X_1,..,X_k)=\{\lambda_1X_1+...+\lambda_k X_k, \lambda \in \textbf{R}^k\}\)

définition soit un processus aléatoire \((Y_t)_{t \in \textbf{Z}}\) stationnaire faible centré. La fonction d’auto-corrélation partielle est définie de la manière suivante:

on remarque que \(r(h) = r_{Y_2,...,Y_h} (Y_1,Y_{h+1})= r_{Y_{k+2},...,Y_{k+h}} (Y_{k+1},Y_{k+h+1})\) car le processus est stationnaire.

Pour obtenir l’auto-corrélation partielle d’un processus, voyons les propriétés suivantes.

propriété soit un processus stationnaire faible \((Y_t)_{t \in \textbf{Z}}\) et \(P_{M_{t-1}^k}(Y_t)=b_{k,1} Y_{t-1}+...+b_{k,k} Y_{t-k}\) sa projection sur son passé, alors on a son auto-corrélation partielle d’ordre \(k\) vaut \(r(k)=b_{k,k}\).

algorithme de Durbin l’agorithme de Durbin permet d’obtenir les auto-corrélation partielles \(r(k)=b_{k,k}\) d’un processus processus stationnaire faible \((Y_t)_{t \in \textbf{Z}}\) via la formule récursive suivante:

Comment se ramener à un processus stationnaire?

library(xts)
date1<- strptime("01/01/1900", "%m/%d/%Y")
date2<- strptime("01/01/2000", "%m/%d/%Y")
Date<-seq.POSIXt(date1,date2,by = "year")
n<-length(Date)
t<-c(1:n)
T<-t/20+1
w=2*pi/5
S<-cos(w*t)
eps<-rnorm(n,0,1)
X<-T*S*eps
X<-xts(X,order.by=Date)
T<-xts(T,order.by=Date)
S<-xts(S,order.by=Date)

X<-T+S+eps
MA<-filter(X, filter=array(1/10,dim=10), method = c("convolution"),
             sides = 2, circular = FALSE)

Les processus ARMA

Théorème de représentation de Wold soit \((X_t)_{t \in \textbf{Z}})\) un processus centré et stationnaire du second ordre, alors on a la décomposition suivante:

\[ X_t=\sum_{j=0}^{+\infty} \psi_j \varepsilon_{t-j} +C_t, \;\; t \in \textbf{Z} \]

où

• \(\psi_0=1\) et \(\sum_{j=0}^{+\infty} |\psi_j|< +\infty\)
• \(\varepsilon_t\) est le processus d’innovation de \(X_t\)
• \(C_t\) est un processus déterministe

les moyennes mobiles s’avèrent donc intéressantes pour modéliser des processus stationnaires. La question reste alors l’estimation des coefficients \(\psi\). Nous verrons par la suite que certaines familles de modèles permettent une représentation parcimonieuse de cette moyenne mobile.

définition on appelle opérateur retard \(L\) l’opérateur qui associe à un processus \((X_t)_{t \in \textbf{Z}}\) le processus \((Y_t)_{t \in \textbf{Z}}\) tel que \(Y_t=LX_t=X_{t-1}\). Cet opérateur est:

• linéaire
• inversible: son inverse est \(L^{-1}=F\) tel que \(FX_t=X_{t+1}\) et est appelé l’opérateur avance

De plus on a:

\[ L^nX_t=X_{t-n} \]

et

\[ (\sum_{i=1}^p a_iL^i)X_t=\sum_{i=1}^p a_i X_{t-i} \]

séries en L on peut définir à l’aide de cet opérateur des séries en \(L\) (ou \(F\)). Soit un processus stationnaire \((X_t)_{t \in \textbf{Z}}\) et une suite \((a_i)_{i \in \textbf{Z}}\) de nombres réels absolument sommable \(\sum_{i \in \textbf{Z}} |a_i|<\infty\), alors le processus défini par:

\[ Y_t= \sum_{i \in \textbf{Z}} a_i X_{t-i}=(\sum_{i \in \textbf{Z}} a_i L^i) X_{t} \]

est stationnaire.

Processus autorégressifs moyennes mobiles (ARMA)

Les processus ARMA(p,q) généralise les modèles autorégressifs et moyennes mobiles. Ces modèles sot très utiles en pratique pour modéliser des séries réelles en nécessitant moins de paramètres que les modèles AR ou MA simples.

définition un processus stationnaire \(X_t\) admet une représentation ARMA(p,q) minimale s’il satisfait

\[ X_t +\sum_{i=1}^p \phi_i X_{t-i} = \varepsilon_t+\sum_{j=1}^q \theta_j \varepsilon_{t-j}, \; \; \forall t \in \textbf{Z} \]

soit, avec les notations précédentes

\[ \Phi(L) X_t= \Theta(L)\varepsilon_t, \; \; \forall t \in \textbf{Z} \]

avec les conditions suivantes:

1. \(\phi_p \neq 0\) et \(\theta_q \neq 0\)
2. \(\Phi\) et \(\Theta\) n’ont pas de racines communes et leurs racines sont de modules \(>1\)
3. \(\varepsilon_t\) est un BB de variance \(\sigma^2\)

propriété si \(X_t\) est un processus stationnaire de représentation ARMA(p,q) minimale:

\[ \Phi(L) X_t= \Theta(L)\varepsilon_t, \; \; \forall t \in \textbf{Z} \]

alors:

1. \(X\) admet la repésentation \(MA(\infty)\):

\[ X_t=\frac{\Theta(L)}{\Phi(L)} \varepsilon_t= \sum_{j=0}^{+\infty} h_j \varepsilon_{t-j}, \;\; h_0=1 \]

2. \(X\) admet la repésentation \(AR(\infty)\):

\[ \frac{\Phi(L)}{\Theta(L)} X_t= \sum_{j=0}^{+\infty} \pi_j X_{t-j}= \varepsilon_t, \;\; \pi_0=1 \]

3. l’innovation de \(X\) est \(\varepsilon\)

détermination des ordres

l’intérêt des modèles ARMA tient dans leur parcimonie. En effet, pour un ajustement aux données de la même qualité, un ARMA nécessite moins de paramètres qu’un AR ou un MA pur (on le comprend en regardant la décomposition AR ou MA infinie). En revanche, le choix des ordre est plus complexe que pour les modèles AR ou MA purs.

Identification et estimation de modèles ARMA

En pratique, lorsque l’on doit ajuster un modèle AR, MA ou ARMA à des données réelles la première question qui se pose est celle du choix des ordres p et q du modèle ARMA (on considère que les AR et MA sont un cas particulier d’ARMA avec respectivement \(q=0\) et \(p=0\)). Pour choisir ces ordres, nous pouvons exploiter les résultats suivants:

Ses résultats sont facilement exploitables dans le cas dun AR ou MA pur. Dans le cas d’un ARMA il existe un moyen de déterminer les ordres \(p\) et \(q\) à l’aide de la méthode du coin basée sur certains déterminant de matrices de corrélations que nous ne développons pas ici. Une autre approche est de considérer un ensemble crédible de modèles ARMA(p,q) puis par de sélectionner un candidat par une méthode de sélection de modèle basée sur des critères d’information de type BIC ou AIC.

Une fois l’ordre et le type de modèle choisis, pour estimer les coefficients des modèles plusieurs approches sont possibles.

1. Equations de Yule-Walker. Nous avons obtenus précédemment que dans le cas d’un processus AR, les coefficients vérifient:

2. Moindre carrés conditionnels. si \(\varepsilon_t\) est gaussien, le vecteur des observations \((X_1,...,X_n)\) est gaussien et on peut alors calculer sa vraisemblance. Une approche simplifiée est de calculer la densité de \((X_{p+1},...,X_n)\) conditionnellement à \((X_1,...,X_p)\) qui est gaussienne. Sachant que \(X_t/X_{t-1},...,X_{t-p}\) est une gaussienne d’espérance \(\sum_{j=1}^p \phi_j X_{t-j}\) et de variance \(\sigma^2\) on a, en utilisant la formule des conditionnements successifs:

\[ f(X_{p+1},...,X_n/X_{1},...,X_p)= \prod_{t=p+1}^{n} f(X_t/X_{t-1},...,X_{t-p}) \]

\[ =\frac{1}{(2\pi \sigma^2)^{\frac{n-p}{2}}} \text{exp}(-\frac{\sum_{t=p+1}^{n} (X_t-\sum_{j=1}^{p}\phi_j X_{t-j} )^2}{2\sigma^2}) \]

maximiser cette vraisemblance conditionnelle en \(\phi\) revient à, en passant au log, minimiser:

\[ (\phi_1,...,\phi_p,\sigma^2) \rightarrow \frac{\sum_{t=p+1}^{n} (X_t-\sum_{j=1}^{p}\phi_j X_{t-j} )^2}{\sigma^2} + (n-p)\ln(\sigma^2) \]

à \(\sigma^2\) fiwé cela correspond au critère des moindres carrés.

Revenons sur le choix de l’ordre d’un ARMA(p,q). Nous avons vu précédemment que l’étude des autocorrélations et autocorrélations partielles peuvent permettre de préselectionner un certains nombre de modèles plausibles. On peut ensuite, une fois les paramètres de ces modèles estimés sélectionner celui qui minimise les critères suivants:

• AIC Akaike Information Criterion, adapté au problème de la prévision, défini par:

\[ \text{AIC}(\phi,\theta,\sigma^2)= -2\log(L(\theta,\phi,\sigma^2)) +2k \]

où \(L\) est la vraisemblance, \(k\) est le ombre de paramètres dans le modèle donc \(p+q\) ou \(p+q+1\) si la constante est dans le modèle.

• BIC Bayesian Information Criterion, adapté au problème de la prévision, défini par:

\[ \text{BIC}(\phi,\theta,\sigma^2)= -2\log(L(\theta,\phi,\sigma^2)) +nk \]

Voilà ci-dessous un exemple obtenu pour un processus AR(p) simulé ainsi:

library(polynom)
n<-500
sigma<-1
set.seed(100)
eps<-rnorm(n,0,sd=sigma)
pol<-poly.calc(x=c(1.5,2,1.6))
pol<--pol/pol[1]
pol

## -1 + 1.791667*x - 1.0625*x^2 + 0.2083333*x^3

x1<-arima.sim(n = n, list(ar = pol[2:length(pol)],innov=eps))
plot(x1)

par(mfrow=c(1,2))
acf(x1)
pacf(x1)

L’estimation d’un modèle AR(p) par maximum de vraisemblance sur cette série se fait via la commande:

 x1.model<-arima(x1, order = c(p,0,0), method = c("ML"), SSinit = c("Rossignol2011"),
                 optim.method = "BFGS", include.mean = F)

On effectue ce calcule pour \(p=1,...,20\) et on obtient:

Les processus ARIMA

Nous avons vu dans un chapitre antérieur comment corriger une série non-stationnaire de composantes déterministes telles qu’une tendance ou une saisonnalité. Entre autre, nous avons étudié l’opérateur de différenciation (rappel):

Notons \(\Delta\) l’opérateur de différenciation: \(\Delta X_t=X_t-X_{t-1}\). L’opérateur de différenciation d’ordre \(k\) correspondant est: \(\Delta^k X_t=\Delta(\Delta^{k-1} X_t)\)

Propriété soit un processus \(y\) admettant une tendance polynomiale d’ordre \(k\):

\[ y_t= \sum_{j=0}^{k} a_j t^j +\varepsilon_t \]

alors le processus \(\Delta y\) admet une tendance polynomiale d’ordre \(k-1\).

Avec les notations des modèles ARMA, on peut remarquer que \(\Delta X_t=(1-L)X_t\) et plus généralement \(\Delta^d X_t=(1-L)^d X_t\). Ainsi un processus ARIMA est défini ainsi:

définition un processus stationnaire \(X_t\) admet une représentation ARIMA(p,d,q) minimale s’il satisfait

\[ \Phi(L) (1-L)^d X_t= \Theta(L)\varepsilon_t, \; \; \forall t \in \textbf{Z} \]

avec les conditions suivantes:

1. \(\phi_p \neq 0\) et \(\theta_q \neq 0\)
2. \(\Phi\) et \(\Theta\), polynômes de degrés resp. \(p\) et \(q\), n’ont pas de racines communes et leurs racines sont de modules \(>1\)
3. \(\varepsilon_t\) est un BB de variance \(\sigma^2\)

Un processus ARIMA(p,d,q) convient pour modéliser une série temporelle comprenant une tendance polynômiale de degrés \(d\), l’opérateur \((1-L)^d\) permettant de transformer un polynôme de degré \(d\) en une constante.

Pour estimer les paramètres d’un modèle ARIMA, on procède de même que pour un ARMA sur le processus différencié \((1-L)^d X_t\).

Les processus SARIMA

Certaines série chronologiques présentent des saisonnalités. Nous avons déjà abordé le cas de séries dont les composantes saisonnières sont déterministes et nous avons procédé par des méthodes de régressions (paramétriques, non-paramétriques) ou par différenciation pour les désaisonnaliser et ensuite leur appliquer une modélisation ARMA. Une autre démarche consiste à intégrer, dans un modèle ARIMA, des décalages multiples de la saisonnalité (par exemple retard de 12 pour des séries mensuelles possédant une saisonnalité annuelle). En théorie, si on choisit des ordres \(p\) et \(q\) suffisamment important, ces retards sont naturellement intégrés au modèle ARIMA(p,d,q), l’inconvénient étant qu’on rajoute un grand nombre de retards intermédiaires potentiellement non-significatifs.

Par exemple, la série suivante:

set.seed(100)
n<-500
sigma<-1
eps<-rnorm(n,0,sd=sigma)
x1<-arima.sim(n = n, list(ar =c(.6,rep(0,10),.5,-.30),innov=eps))

Où les lag \(12\) et \(13\) sont significatifs alors que les lag intermédaires ne sont pas significativement différents de 0. L’autocorrélogramme indiquant clairement une saisonnalité (pas immédiat en se basant sur la trajectoire).

Box et Jenkins propose des modèles multiplicatifs SARIMA définis ainsi.

définition un processus stationnaire \(X_t\) admet une représentation SARIMA[(p,d,q);(P,D,Q);s] minimale s’il satisfait

\[ (1-L)^d \Phi_p(L) (1-L^s)^D \Phi_P(L^s) X_t= \Theta_q(L)\Theta_Q(L^s)\varepsilon_t, \; \; \forall t \in \textbf{Z} \]

avec des conditions similaires que pour les ARIMA.

\(s\) correspond à la période du processus SARIMA qu’on peut identifier en regardant l’autocorrélogramme (cf exemple précédant) ou la densité spectrale. Les entiers \(d\) et \(D\) sont choisis de sorte que la série différenciée: \((1-L)^d (1-L^s)^D\) soit stationnaire (en pratique regarder la trajectoire, les autocorrélations). Les ordres \(p\) et \(q\) s’obtiennent comme pour les modèles ARMA(p,q) classiques (autocorrélation partielle et simple), les ordres \(P\) et \(Q\) en regardant les ordre multiples de \(s\) de l’autocorrélogramme. En pratique, on commencera par différencier en \(1-B^s\) (choisir \(D\)) avant de choisir \(d\) car \(1-B^s=(1-B)(1+B+B^2+...+B^{k-1})\).

L’approche de Box-Jenkins

L’approche de box jenkins, du nom des statisticiens George Box et Gwilym Jenkins, est une méthode empirique de choix et construction de modèle SARIMA mise au point dans les années 1970. Pour simplifier les choses nous nous intéressons ci-dessous aux processus ARIMA.

• le premier paramètre à choisir est le degré de différentiation \(d\). Plusieurs moyens sont possibles pour détecter une non-stationnarité. La représentation temporelle de la série peut faire apparaitre des tendances polynomiale. On peut aussi calculer les autocorrélation empiriques et analyser la vitesse de décroissance vers 0 de cette fonction, si cette décroissance est lente (plus lente qu’exponentiel) on peut suspecter une non-stationnarité. Notons qu’en pratique le cas \(d>2\) est rarement rencontré. En effet, il est dangereux de sur-différencier un processus, cela peut conduire à une non-inversibilité des ARMA associés.

Par exemple dans le cas du processus \(X_t\) suivant:

• les ordres \(p\) et \(q\) des modèles AR et MA sont ensuites obtenus en regardant les autocorrélations et autocorrélations partielles. Bien souvent, ces ordres n’apparaissent pas de manière évidente. On peut dans ce cas obtenir des bornes supérieures pour \(p\) et \(q\) puis sélectionner un modèle en minimisant un critère pénalisé de type AIC ou BIC.

• l’estimation des paramètres se fait ensuite par les méthodes précédemment citées (maximum de vraisemblance, moindres carrés conditionnels).

• d’autres approches que la différenciation sont possibles pour “stationnariser” un processus. Par exemple, il est clair que si un processus possède une tendance non-polynomiale, par exemple exponentielle ou logarithmique, l’opérateur \((1-L)^d\) ne suffira pas à le rendre stationnaire. On peut alors soit modéliser par régression la partie déterministe du processus soit appliquer des transformations à la série. Une transformation courante est le passage au log de la série (pour gérer des phénomènes de variance multiplicative par exemple).

Diagnostics et tests

Une fois le modèle estimé, il est important de valider ou non nos choix de modélisation. Pour cela différents tests et diagnostics existent.

Lissage exponentiel simple

Un algorithme de base pour la prévision de séries temporelles univariées est le lissage exponentiel. On peut voir le lissage exponentiel comme une méthode de prévision mais également, comme son nom l’indique, comme une technique de lissage de données.

définition soit une série temporelle \(y_t\). On appelle lissage exponentiel simple de paramètre \(\alpha \in [0,1]\) de cette série le processus \(\widehat{y}_{t}\) définie ainsi:

\[ \widehat{y}_{t+1/t}= \alpha y_t +(1-\alpha) \widehat{y}_{t/t-1} \]

on a donc:

\[ \widehat{y}_{t+1/t}=\sum_{i=0}^{t-1} \alpha (1-\alpha)^i y_{t-i} \]

la prévision de l’instant \(t+1\) est donc une somme pondérée des valeurs passées de la série, les poids décroissant exponentiellement dans le passé. La mémoire de la prévision dépend de \(\alpha\). Plus \(\alpha\) est proche de 1 plus les observations récentes influent sur la prévision, à l’inverse un \(\alpha\) proche de 0 conduit à une prévision très stable prenant en compte un passé lointain.

Une autre façon d’écrire le lissage exponentiel (error correction form):

\[ \widehat{y}_{t+1/t}= \widehat{y}_{t/t-1} + \alpha (y_t -\widehat{y}_{t/t-1} ) \]

\(\widehat{y}_{t+1/t}\) est une prévision à horizon 1. Il est parfois nécessaire d’effectuer une prévision à un horizon \(h\) quelconque. On notera par la suite \(\widehat{y}_{t+h/t}\) la prévision de \(y_{t+h}\) conditionellement à \((y_1,...,y_t)\). Pour le lissage exponentiel simple cette prévision est tout simplement: \(\widehat{y}_{t+h/t}=\widehat{y}_{t+1}\) car on approxime le futur de la série à une constante.

le prédicteur obtenu par lissage exponentiel simple peut être vu comme celui qui minimise un problème de moindres carrés pondérés.

En effet, si on note: \[ \widehat{a}=\text{argmin}_a \sum_{i=0}^{t-1} (1-\alpha)^i(y_{t-i}-a)^2 \]

on a \(\widehat{a}= \sum_{i=0}^{t-1} (1-\alpha)^i y_{t-i} / \sum_{i=0}^{t-1} (1-\alpha)^i\),

soit \(\widehat{a}= \frac{\alpha}{1-(1-\alpha)^t} \sum_{i=0}^{t-1} (1-\alpha)^i y_{t-i}\)

ce qui, si \(t\) est suffisamment grand tend vers l’estimateur de lissage exponentiel simple définie précédemment.

Si on défnie \(y_t\) comme \((0,...,0, y_1,...,y_t)\) ie valant \(0\) pour les indices temporelles inférieurs à 1 alors

\[ \widehat{y}_{t+1}=\text{argmin}_a \sum_{i=0}^{\infty} (1-\alpha)^i(y_{t-i}-a)^2 \]

Voilà un exemple de lissage exponentiel simple sur des données simulées de la façon suivante:

set.seed(150)
n<-100
t<-c(1:n)
eps<-rnorm(n,0,1/2)
T<-log(t)+pmax(t-50,0)/10-pmax(t-70,0)/5
X<-T+eps

Lissage exponentiel double (ou de Holt)

Holt (1957) a étendu le lissage exponentiel simple au cas du lissage exponentielle linéaire. L’idée est d’ajuster une droite au lieu d’une constante dans l’approximation locale de la série.

définition soit une série temporelle \(y_t\). On appelle lissage exponentiel double (ou de Holt) de paramètre \(\alpha \in [0,1]\) de cette série le processus \(\widehat{y}_{t}\) définie ainsi:

\[ \widehat{y}_{t+h/t}=l_t+b_t h \]

avec:

\[ \left \{ \begin{array}{l @{ = } l} l_t & l_{t-1}+b_{t-1}+(1-(1-\alpha)^2) (y_t-\widehat{y}_{t/t-1})\\ b_t & b_{t-1}+\alpha^2 (y_t-\widehat{y}_{t/t-1})\\ \end{array} \right. \]

\(l_t\) et \(b_t\) minimise à chaque instant:

\[ (\widehat{l},\widehat{b})=\text{argmin}_{l,b} \sum_{i=0}^{\infty} (1-\alpha)^i(y_{t-i}-(l+bi))^2 \]

un exemple de prévisions à horizon \(1,2,...,20\):

Lissage exponentiel de Holt-Winters

Cette approche est une généralisation du lissage double, qui permet entre autre de proposer les modèles suivants:

• tendance linéaire locale
• tendance linéaire locale + saisonnalité (modèle additif)
• tendance linéaire locale * saisonnalité (modèle multiplicatif)

Dans ce cas 2 paramètres de lissage entrent en jeu et on ajuste au voisinage de \(t\) un fonction linéaire \(l_t+hb_t\), \(h\) étant l’horizon de prévision.

définition soit une série temporelle \(y_t\). On appelle lissage exponentiel double de Holt-Winters de paramètres \(\alpha \in [0,1]\) et \(\beta \in [0,1]\) de cette série le processus \(\widehat{y}_{t}\) définie ainsi:

\[ \widehat{y}_{t+h}=l_t + h b_t \]

avec

\[ \left \{ \begin{array}{l @{ = } l} l_t & \alpha y_t +(1-\alpha)(l_{t-1}+b_{t-1})\\ b_t & \beta (l_t-l_{t-1})+(1-\beta)b_{t-1}\\ \end{array} \right. \]

la encore \(l_t\) est une estimation du niveau de la série, \(b_t\) de sa pente (localement en temps).

Pour ajuster une composante saisonnière, on considère le modèle (localement au voisinage de \(t\)):

\[ \widehat{y}_{t+h}=l_t + h b_t +s_t \]

ou \(s_t\) est une composante périodique de période \(T\).

définition soit une série temporelle \(y_t\). On appelle lissage exponentiel double de Holt-Winters saisonnier de paramètres \(\alpha \in [0,1]\), \(\beta \in [0,1]\), \(\delta \in [0,1]\) de cette série le processus \(\widehat{y}_{t}\) définie ainsi:

\[ \widehat{y}_{t+h}=l_t + h b_t+s_t \]

avec

\[ \left \{ \begin{array}{l @{ = } l} l_t & \alpha (y_t-s_{t-T}) +(1-\alpha)(l_{t-1}+b_{t-1})\\ b_t & \beta (l_t-l_{t-1})+(1-\beta)b_{t-1}\\ s_t & \delta (y_t- l_t) +(1-\delta)s_{t-T}\\ \end{array} \right. \]

Reprenons les données simulées précédamment en ajoutant une composante périodique de période \(T=10\).

n<-100
t<-c(1:n)
eps<-rnorm(n,0,1/2)
w<-2*pi/10
S<-cos(w*t)
T<-log(t)+pmax(t-50,0)/10-pmax(t-70,0)/5
X<-T+S+eps

Voilà ce qu’on obtient en appliquant un lissage exponentiel saisonnier de paramètres \(\alpha =0.2\), \(\beta =0.2\), \(\delta =0.2\):

Il existe une variante multiplicative de ce lissage saisonnier.

définition soit une série temporelle \(y_t\). On appelle lissage exponentiel double de Holt-Winters saisonnier multiplicatif de paramètres \(\alpha \in [0,1]\), \(\beta \in [0,1]\), \(\delta \in [0,1]\) de cette série le processus \(\widehat{y}_{t}\) définie ainsi:

\[ \widehat{y}_{t+h}=(l_t + h b_t)s_t \]

\[ \left \{ \begin{array}{l @{ = } l} l_t & \alpha \frac{y_t}{s_{t-T}} +(1-\alpha)(l_{t-1}+b_{t-1})\\ b_t & \beta (l_t-l_{t-1})+(1-\beta)b_{t-1}\\ s_t & \delta \frac{y_t}{l_t} +(1-\delta)s_{t-T}\\ \end{array} \right. \]

HoltWinters(x, alpha = NULL, beta = NULL, gamma = NULL,
            seasonal = c("additive", "multiplicative"),
            start.periods = 2, l.start = NULL, b.start = NULL,
            s.start = NULL,
            optim.start = c(alpha = 0.3, beta = 0.1, gamma = 0.1),
            optim.control = list())