Récupération simultanée de la date et de l’heure

require(chron)
x=chron(dates=exemple01$Date.app,times=exemple01$H.deb.app,
        format=c(dates="d/m/Y",times ="h:m:s"))
x[1:2]

## [1] (01/06/06 09:03:51) (02/06/06 17:28:46)

Le résultat sont des dates-heures.

Récupération simultanée de la date et de l’heure comme heure POSIX

L’heure POSIX (POSIX timestamp) est une mesure utilisée principalement dans les systèmes qui respectent la norme POSIX (Portable Operating System Interface). Il s’agit du nombre de secondes écoulées depuis le 1er janvier 1970 à 00:00:00 UTC jusqu’à l’événement à dater, hors secondes intercalaires.

R mesure le temps à l’aide de l’heure POSIX. La fonction strptime() assure la conversion entre représentation comme chaîne de caractères (string) et date-heure POSIX. Si on ne spécifie pas le format d’une date fournie à R, il considère qu’elle est donnée sous la forme internationale %Y-%m-%d, c’est-à-dire : année sur 4 positions–mois sur 2 positions–jour sur 2 positions, les champs étant séparés par des « - ».

Fabriquons le vecteur des chaînes de caractères date-heure puis convertissons-le par strptime() :

dh=paste(exemple01$Date.app,exemple01$H.deb.app)
xp=strptime(dh,"%d/%m/%Y %H:%M:%S")
xp[1:2]

## [1] "2006-06-01 09:03:51 CEST" "2006-06-02 17:28:46 CEST"

Extraction des composantes

Des fonctions permettent l’extraction des composantes d’un objet de type chron ou POSIX :

options(digits=12)
dates(x)[1:4]

## [1] 01/06/06 02/06/06 02/06/06 03/06/06

times(x)[1:4]

## Time in days:
## [1] 13300.3776736 13301.7283102 13301.7853472 13302.6540509

hours(xp);minutes(xp);seconds(xp)

## [1]  9 17 18 15 15

## [1]  3 28 50 41 50

## [1] 51 46 54 50 41

months(xp)[1:3]; days(xp)[1:3]; weekdays(xp)[1:3]

## [1] "juin" "juin" "juin"

## [1] 1 2 2
## 31 Levels: 1 < 2 < 3 < 4 < 5 < 6 < 7 < 8 < 9 < 10 < 11 < 12 < 13 < ... < 31

## [1] "jeudi"    "vendredi" "vendredi"

years(xp)

## [1] 2006 2006 2006 2006 2006
## Levels: 2006

times() extrait l’heure en fraction de jour depuis le 01/01/1970, dates() et times() ne s’appliquent pas à une heure POSIX, hours(), minutes() et seconds() permettent d’extraire des éléments de date.

Extraction naïve des éléments de date

On peut alternativement décomposer l’heure à partir de la chaîne de caractères correspondante, récupérer l’année, le mois… sans considération de structures de date ou d’heure,

heure0=exemple01$H.deb.app
heure=substr(heure0,1,2)
minute=substr(heure0,4,5)
seconde=substr(heure0,7,8)
an=substr(exemple01$Date.app,7,10)
an.num=as.numeric(an)

mais on se prive ainsi de toute structure de date ou de date-heure.

Les structures de séries temporelles dans

Une série temporelle unidimensionnelle est un vecteur à chaque élément duquel est associée une date ou une date-heure, ou si l’on préfère, une série temporelle est une suite de couples (date, valeur), ordonnée sur la date. Une série temporelle multidimensionnelle est une matrice à chaque ligne de laquelle est associée une date ou une date-heure.

La fonction ts() fait partie de stats, chargé au lancement de R. Elle permet de créer des séries temporelles à temps régulier. ?ts fournit la syntaxe et de très utiles exemples.

Fabriquons une série mensuelle supposée commencer en février 2005 à partir d’un vecteur \(x\) obtenu par simulation de 30 valeurs indépendantes d’une normale de moyenne 2, de variance 1, et dont on ne conserve que 3 décimales. L’unité de temps est l’année et il y a 12 observations par unité de temps. Février est la deuxième observation dans une unité de temps, d’où la syntaxe :

set.seed(493)
x=2+round(rnorm(15),3)
x.ts=ts(x,start=c(2005,2),frequency=12)
x.ts

##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 2005       1.545 0.980 1.162 3.208 2.666 2.739 0.413 0.560 1.115 3.298 0.830
## 2006 0.886 1.354 3.141 2.948

Examinons la structure de x.ts :

str(x.ts)

##  Time-Series [1:15] from 2005 to 2006: 1.54 0.98 1.16 3.21 2.67 ...

frequency(x.ts)

## [1] 12

cycle(x.ts)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2005       2   3   4   5   6   7   8   9  10  11  12
## 2006   1   2   3   4

set.seed(8515)
y=10+0.5*round(rnorm(30),3)
xy.ts=ts(cbind(x,y),start=c(2005,2),frequency=12)
str(xy.ts)

##  Time-Series [1:30, 1:2] from 2005 to 2008: 1.54 0.98 1.16 3.21 2.67 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:2] "x" "y"

xy.ts[1:4,]

##          x       y
## [1,] 1.545  9.7715
## [2,] 0.980 11.0025
## [3,] 1.162 10.1980
## [4,] 3.208 10.4780

La définition

Une série temporelle \(\{y_t\}\) est un processus stochastique. Elle est dite strictement stationnaire si la distribution conjointe de \((y_{t_1},\ldots,y_{t_k})\) est identique à celle de \((y_{t_1+t},\ldots,y_{t_k+t})\) quels que soient \(k\) le nombre d’instants considérés, le choix des \(k\) instants eux-mêmes \((t_1,\ldots,t_k)\) et le décalage \(t\). La stationnarité stricte est la propriété d’une série temporelle que la distribution jointe de tout sous-vecteur issu d’elle est invariante quand on translate les instants d’une même quantité.

Question : Ne suffit-il pas de considérer le cas \(k=2\) ? C’est-à-dire : En quoi considérer les triplets de variables aléatoires et non seulement les paires enrichit-il la notion de stationnarité de la série ?

Cette condition est difficile à vérifier. C’est pourquoi on lui préfère une version affaiblie. Une série temporelle \(\{y_t\}\) est faiblement stationnaire si par définition

• \(E[y_t]=\mu\quad\), indépendamment de \(t\),
• \(Cov(y_t,y_{t-l})\) ne dépend que de l’entier \(l\), expression que l’on notera alors \(\gamma_l:=Cov(y_t,y_{t-l})\).

Affirmation : Pour les séries temporelles normalement distribuées, les deux notions de stationnarité coïncident.

Un vecteur gaussien est en effet entièrement déterminé par son espérance \(\mu\) et sa matrice de variance-covariance \(\Sigma\), et celle-ci est entièrement déterminée par la donnée de la covariance pour chaque paire (\(k=2\)) de composantes du vecteur gaussien.

Les moyens de constater une éventuelle stationnarité

Comment apprécier la stationnarité d’une série temporelle ? Si l’on dispose d’une trajectoire de la série \(\{y_t\}\), on peut s’en faire une première idée par l’observation du chronogramme de la trajectoire. Une condition nécessaire de stationnarité est que la moyenne et la variance soient constantes par rapport au temps. Tendance plate, horizontale et fluctuation d’amplitude à peu près constante.

La covariance \(\gamma_l:=Cov(y_t,y_{t-l})\) est appelée autocovariance d’ordre \(l\) (lag-\(l\) autocovariance). Vue comme fonction de \(l\), elle est appelée la fonction d’autocovariance (\(l\in\mathbb{Z}\)) . Cette fonction vérifie :

• \(\gamma_0=Var(y_t)\geq 0\),
• \(|\gamma_l|\leq \gamma_0\),
• \(-\gamma_l=\gamma_{-l}\quad\),fonction impaire, que par conséquent l’on ne représente que sur \(\mathbb{N}\).

La variance d’une combinaison linéaire de \(n\) v.a. \(y_{t_1},\ldots,y_{t_k}\) est positive. Une condition qui exprime cela est :

\[ \forall n\mbox{ et }\forall \mbox{ vecteur } a=(a_1,\ldots,a_n), \sum_{i,j=1}^n a_i\cdot\gamma_{i-j}\cdot a_j\geq0 \]

On dit que \(\gamma_l\) est une fonction de type positif.

On définit le coefficient d’autocorrélation d’ordre \(l\) par

\[ \begin{eqnarray*} \rho_l &:=& \frac{Cov(y_t,y_{t-l})}{\sqrt{Var(y_t)\cdot Var (y_{t-l})}} \\ &=& \frac{Cov(y_t,y_{t-l})}{Var(y_t)} \\ &=& \frac{\gamma_l}{\gamma_0} \end{eqnarray*} \]

\(l\mapsto \rho_l\) est la fonction d’autocorrélation théorique de la série \(\{y_t\}\). Son graphe s’appelle le corrélogramme. On a que \(\rho_0=1\) et \(-1\leq\rho_l\leq 1\).

La fonction d’autocovariance empirique \(l\mapsto \hat{\gamma}_l\) est donnée par :

\[ \hat{\gamma}_l:=\frac{1}{T}\cdot\left(\sum_{t=l-1}^T (y_t-\bar{y})\cdot(y_{t-l}-\bar{y})\right) \]

Quant à la fonction d’autocorrélation empirique \(l\mapsto \hat{\rho}_l\), elle est donnée par :

\[ \hat{\rho}_l:=\frac{\hat{\gamma}_l}{\sum_{t=1}^T (y_t-\bar{y})^2}\,,\quad 0\leq l\leq T-1. \]

Définition : un bruit blanc \(\{z_t\}\) est une suite de v.a. non corrélées (mais pas nécessairement indépendantes), de moyenne nulle et de variance constante \(\sigma_z^2\). C’est donc une série faiblement stationnaire. Un bruit blanc gaussien, où en plus \(z_t\sim\mathcal{N}(0,\sigma_z^2)\), est noté \(BB(0,\sigma_z^2)\). Il est strictement stationnaire parce que gaussien et faiblement stationnaire.

Affirmation : Si \(y_t\), \(t=1,\ldots,T\), est une observation d’une suite de v.a. i.i.d. de moment d’ordre 2 fini, i.e. \(E[y_t^2]<\infty\), alors \(\hat{\rho}_l\sim\mathcal{N}(0,\frac{1}{T})\).

En s’appuyant sur cette affirmation, on peut tracer des intervalles autour de 0 qui doivent contenir les \(\hat{\rho}_l\). Cependant il n’est pas facile de vérifier qu’une série est formée de v.a. i.i.d.. On peut par contre tester la nullité des coefficients d’autocorrélation :

Théorème du portemanteau (test de blancheur) : Soit la série observée \(y_t\) , \(t=1,\ldots,T\). La statistique suivante, dépendant d’un décalage \(h\) qu’il s’agit de faire varier, en répétant le test pour plusieurs valeurs de \(h\), permet de tester l’hypothèse

\(H_0 : \rho_1=\rho_2=\ldots=\rho_h=0\)

contre l’hypothèse

\(H_1\) : au moins un des \(\rho_j\) est non nul.

Voici cette statistique (statistique de Box-Pierce) :

\[ \begin{eqnarray*} Q(h) &=& \sum_{j=1}^h \left(\frac{\hat{\rho}_j -0}{1/\sqrt{T}}\right)^2 && \text{puisque la variance des }\hat{\rho}_j \text{vaut }1/T \text{ (selon l'affirm. ci-dessus)} \\ &=& T\cdot\sum_{j=1}^h \hat{\rho}_j^2\quad\sim\chi^2_h \\ \end{eqnarray*} \]

Remarque : Quand ce test est appliqué sur les résidus d’un ajustement estimant \(m\) paramètres, par exemple dans le cadre d’une régression multilinéaire, la loi approchée de la statistique ci-dessus, à supposer que l’hypothèse \(H_0\) soit vérifiée, est \(\chi^2_{h-m}\).
Voir: test de Durbin-Watson, un test d’absence d’autocorrélation d’ordre 1 sur le résidu d’une régression linéaire

Remarque : Variante du test pour de petits échantillons : la statistique de Ljung-Box \(Q^*(h)=T\cdot(T+2)\cdot\sum_{k=1}^h \hat{\rho}_k^2/(T-k)\), qui a une distribution mieux approchée par un \(\chi^2\) que la statistique de Box-Pierce.

require(caschrono)
set.seed(123)
y1=arima.sim(n=100,list(ar=-.7),sd=sqrt(4))
y2=arima.sim(n=100,list(ar=c(rep(0,11),-.7)),sd=sqrt(4))
retard=c(3,6,9,12)
a1=Box.test.2(y1,nlag=retard,type="Ljung-Box",decim=2)
a2=Box.test.2(y2,nlag=retard,type="Ljung-Box",decim=2)
a12=cbind(a1,a2[,2])
colnames(a12)=c("Retard","p-val. y1","p-val.y2")
a12

##      Retard p-val. y1 p-val.y2
## [1,]      3         0     0.66
## [2,]      6         0     0.49
## [3,]      9         0     0.43
## [4,]     12         0     0.00

data(mtcars)
head(mtcars)

##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

model <- lm(mpg ~ disp+wt, data=mtcars)

library(car)

## Le chargement a nécessité le package : carData

durbinWatsonTest(model)

##  lag Autocorrelation D-W Statistic p-value
##    1        0.341622      1.276569   0.036
##  Alternative hypothesis: rho != 0

L’opérateur retard (the lag operator)

L’opérateur retard est défini comme suit : \[ By_t:=y_{t-1}. \] On a par exemple que \(B^2y_t=B(By_t)=y_{t-2}\). On peut construire, à l’aide de l’opérateur retard, l’opérateur de différence saisonnière : comparer janvier avec janvier, février avec février, etc. Il suffit de définir

\[ \Delta_{12}:=1-B^{12}. \]

On a donc \(\Delta_{12}\,y_t=y_t-y_{t-12}\).

Il est possible de calculer algébriquement avec l’opérateur retard, notamment et en parfaite analogie avec l’opérateur multiplication (considéré alors comme un autre opérateur commutatif) :

\[ \begin{array}{ll} (Id-B)(Id+B+B^2+B^3+B^4) \\ \mbox{pas de multiplication ici : la seconde parenthèse est l'argument de la première} \\ = Id-B^5 \\ \mbox{et, si }y_t\mbox{ est tel que la répétition du lag à l'infini }B^\infty\,y_t \mbox{ donne zéro,} \\ (Id-B)\left(\sum_{k=0}^\infty B^k\right)\\ = Id\,, \\ \mbox{si bien que} \\ \sum_{k=0}^\infty B^k \\ =\frac{Id}{Id-B} \,, \\ \mbox{cette formule fractionnaire ne faisant sens que si elle signifie effectivement une série} \\ \mbox{de puissances, la fraction ne signifiant dans ce contexte pas une division} \\ \mbox{mais bien une opération réciproque. Nous ne faisons en effet qu'emprunter} \\ \mbox{à l'écriture multiplicative, mais ne nous y trompons pas !} \end{array} \]

A l’aide de l’opérateur retard et en écrivant l’opérateur identité \(Id\) par 1 pour rester dans l’analogie avec la multiplication, on obtient à nouveau le résultat précédent à propos du modèle AR(1) :

\[ \begin{eqnarray*} (1-\phi\cdot B)y_t &=& c+z_t \\ \mbox{et, en supposant que }|\phi|<1\,, \\ y_t &=& \underbrace{\frac{c}{1-\phi\cdot B}}_{=\frac{c}{1-\phi}=E[y_t]}+\frac{1}{1-\phi\cdot B}\cdot z_t \\ \mbox{Il est clair alors que } \\ \frac{1}{1-\phi\cdot B} &=& 1+\phi\cdot B+\phi^2\cdot B^2+\ldots \end{eqnarray*} \]

A quel genre de processus aléatoire avons-nous à faire ?

Nous avons obtenu précédemment que

\[ E[y_t]=\mu=\frac{c}{1-\phi}. \]

Nous pouvons comprendre rien qu’en regardant cette expression de la moyenne que si \(\phi=1\) le processus n’est pas stationnaire -On parle alors de marche aléatoire- et que si \(|\phi|>1\), le processus est explosif.

Processus autorégressifs d’ordres supérieurs

Processus autorégressif d’ordre p, noté AR(p) :

\[ y_t=c+\phi_1\cdot y_{t-1}+\phi_2\cdot y_{t-2}+\ldots+\phi_p\cdot y_{t-p}+z_t\,, \]

où bien évidemment \(\phi_p\neq0\). Avec l’opérateur retard, on peut écrire cette autorégression comme

\[ \Phi(B)y_t=c+z_t\,, \]

avec \(\Phi(B):=1-\phi_1\cdot B-\phi_2\cdot B^2-\ldots-\phi_p\cdot B^p\), et où \(\Phi\) est appelé l’opérateur d’autorégression, ici d’ordre \(p\).

Etudions la question de la stationnarité au travers du cas \(p=2\).

Appelons \(s_1\) et \(s_2\) les racines (\(\in\mathbb{C}\)) de

\[ 1-\phi_1\cdot z-\phi_2\cdot z^2=0. \]

On a que

\[ \begin{eqnarray*} && 1-\phi_1\cdot z-\phi_2\cdot z^2 \\ &=& (1-z/s_1)(1-z/s_2) \end{eqnarray*} \]

et on voit que le développement en série de l’inverse de ce polynôme est possible si ses racines sont en norme strictement supérieures à 1. Dans ce cas, \(\{y_t\}\) est stationnaire, de moyenne \(\mu=\frac{1}{1-\phi_1-\phi_2}\).

On le comprend en considérant le calcul suivant :

\[\begin{eqnarray*} && \frac{1}{1-\phi_1\cdot z-\phi_2\cdot z^2} \\ &=& \frac{1}{1-z/s_1}\cdot \frac{1}{1-z/s_2} \\ &=& \left(1+\frac{z}{s_1}+\frac{z^2}{s_1^2}+\cdots\right)\cdot \left(1+\frac{z}{s_2}+\frac{z^2}{s_2^2}+\cdots\right) \\ &=& 1+z\cdot\left(\frac{1}{s_1}+\frac{1}{s_2}\right)+z^2\cdot \left(\frac{1}{s_1^2}+\frac{1}{s_2^2}+\frac{1}{s_1\cdot s_2}\right)+\cdots \end{eqnarray*}\]

Pour \(p>2\), il en va de même :

Affirmation : \(\{y_t\}\) se laisse écrire comme une série linéaire si les racines du polynôme sont toutes strictement plus grandes que 1 en norme. Dans ce cas, on a que

\[ \begin{array}{ll} E[y_t]=\mu=\frac{c}{1-\phi_1-\phi_2-\ldots-\phi_p}\,, \\ \mbox{ainsi que} \\ y_t=\mu+\frac{1}{1-\phi_1\cdot B-\phi_2\cdot B^2-\ldots-\phi_p\cdot B^p}\cdot z_t. \end{array} \]

Fonction d’autocorrélation d’un AR

Partant de la représentation en MA(\(\infty\)) d’un AR(1), on obtient :

\[\begin{eqnarray*} Var(y_t) &=& Var(1+\phi\cdot B + \phi^2\cdot B^2+\cdots)z_t \\ &=& Var(z_t)+\underbrace{Var(\phi\cdot B z_t)}_{\phi^2\cdot Var(B z_t)} + \underbrace{Var(\phi^2\cdot B^2 z_t)}_{\phi^4\cdot Var(B^2 z_t)}+\cdots \\ &=& \sigma_z^2\cdot(1+\phi^2+\phi^4+\cdots) \\ &=& \frac{\sigma_z^2}{1-\sigma_z^2} \end{eqnarray*}\]

La fonction d’autocorrélation est donc

\[\begin{eqnarray*} \rho_l &=& \frac{Cov(y_t,y_{t-l})}{Var(y_t)} \\ &=& \frac{Cov(y_{t},\,\phi^l\cdot y_{t})}{Var(y_t)} \\ &=& \phi^l \end{eqnarray*}\]

En effet :

\[ \begin{eqnarray*} Cov(y_t,y_{t-l}) &=& E\left[(y_t-E[y_t])\cdot(y_{t-l}-E[y_{t-l}])\right] \\ &=& E\left[\left(\sum_{j=0}^\infty\phi^j\cdot z_{t-j}\right)\cdot \left(\sum_{k=0}^\infty\phi^k\cdot z_{t-l-k}\right) \right] \\ &=& E\left[\left(\sum_{j=0}\sum_{k=0}\phi^{j+k}\cdot z_{t-j}\cdot z_{t-l-k} \right) \right] \\ &=& \sum_{j=0}\sum_{k=0\\k+l\neq j}\underbrace{E[z_{t-j}\cdot z_{t-l-k}]}_{\underbrace{E[z_{t-j}]\cdot E[z_{t-l-k}]}_{=0\cdot 0}} + \sum_{k=0}^\infty\phi^{2k+l}\underbrace{E[z_{t-l-k}^2]}_{Var(z_{t-l-k})+\underbrace{E[z_{t-l-k}]}_{=0}^2} \\ &=& \sum_{k=0}^\infty\phi^{2k+l}Var(z_{t-l-k}) \\ &=& \phi^l\cdot\sum_{k=0}^\infty\phi^{2k}\cdot\sigma_z^2 \end{eqnarray*} \]

Affirmation : Cette fonction décroît exponentiellement vers 0, en oscillant si \(\phi<0\).

Affirmation : La fonction d’autocovariance d’un AR(\(p\)) obéit à : \[ \begin{align*} \gamma_0 &= \phi_1\gamma_1 + \phi_2\gamma_2 + \cdots + \phi_p\gamma_p + \sigma_z^2, \\ \gamma_l &= \phi_1\gamma_{l-1} + \phi_2\gamma_{l-2} + \cdots + \phi_p\gamma_{l-p}, \quad l \geq 1. \end{align*} \]

et la fonction d’autocorrélation d’un AR(\(p\)) obéit à :

\[ \begin{equation*} \rho_l = \phi_1\rho_{l-1} + \phi_2\rho_{l-2} + \cdots + \phi_p\rho_{l-p}, \quad l \geq 1. \end{equation*} \]

Preuve : Il faut considérer un AR(\(p\)) qu’on supposera stationnaire :

\[ y_t=c+\phi_1\cdot y_{t-1}+\phi_2\cdot y_{t-2}+\ldots+\phi_p\cdot y_{t-p}+z_t. \]

Dans ce cas, on peut prendre l’espérance de part et d’autre de l’équation :

\[\begin{eqnarray*} \mu=c+\phi_1\cdot\mu+\phi_2\cdot\mu+\ldots+\phi_p\cdot\mu \\ \end{eqnarray*}\]

En soustrayant cette équation-ci à cette équation-là, on obtient

\[ \begin{eqnarray*} && y_t &=& c+\phi_1\cdot y_{t-1}+\phi_2\cdot y_{t-2}+\ldots+\phi_p\cdot y_{t-p}+z_t \\ \Leftrightarrow\qquad && y_t-\mu &=& c+\phi_1\cdot y_{t-1}+\phi_2\cdot y_{t-2}+\ldots+\phi_p\cdot y_{t-p}+z_t \\ && && -c-\phi_1\cdot\mu-\phi_2\cdot\mu-\ldots-\phi_p\cdot\mu\\ \Leftrightarrow\qquad && y_t-\mu &=& \phi_1\cdot(y_{t-1}-\mu)+\phi_2\cdot(y_{t-2}-\mu)+\ldots+\phi_p\cdot(y_{t-p}-\mu)+z_t.\\ \end{eqnarray*} \]

Si maintenant nous multiplions de chaque côté de cette dernière équation par \((y_{t-l}-\mu)\) puis prenons l’espérance de chaque côté, nous obtenons les autocovariances :

\[\begin{eqnarray*} \Leftrightarrow\qquad && y_t-\mu &=& \phi_1\cdot(y_{t-1}-\mu)+\phi_2\cdot(y_{t-2}-\mu)+\ldots+\phi_p\cdot(y_{t-p}-\mu)+z_t\\ \Leftrightarrow\qquad && (y_t-\mu)\cdot(y_{t-l}-\mu) &=& \phi_1\cdot(y_{t-1}-\mu)\cdot(y_{t-l}-\mu)+\ldots+\phi_p\cdot(y_{t-p}-\mu)\cdot(y_{t-l}-\mu)+z_t\\ \\ \Leftrightarrow\qquad && \gamma_l &=& \begin{cases} \phi_1\cdot\gamma_{1}+\phi_2\cdot\gamma_{2}+\ldots+\phi_p\cdot\gamma_{p}+\sigma_z^2 & \text{si } l=0 \\ \phi_1\cdot\gamma_{l-1}+\phi_2\cdot\gamma_{l-2}+\ldots+\phi_p\cdot\gamma_{l-p} & \text{si } l < 1,2,\ldots \end{cases}\\ \end{eqnarray*}\]

Ceci démontre l’affirmation.

On appelle équations de Yule-Walker la présentation matricielle des \(p+1\) premières équations ci-dessus, i.e. pour les valeurs de \(l\) allant de 0 à \(p\).

\[ \begin{equation} \begin{cases} \Gamma_p \cdot\phi = \gamma_p \\ \sigma_z^2 = \gamma_0 - \phi^T\cdot \gamma_p \end{cases} \end{equation} \]

où

\[ \Gamma_p = \begin{pmatrix} \gamma_{1-1} & \gamma_{2-1} & \dots & \gamma_{p-1} \\ \gamma_{1-2} & \gamma_{2-2} & \dots & \gamma_{p-2} \\ \vdots & \vdots & \ddots & \vdots \\ \gamma_{1-p} & \gamma_{2-p} & \dots & \gamma_{p-p} \end{pmatrix} =\begin{pmatrix} \gamma_{0} & \gamma_{1} & \dots & \gamma_{p-1} \\ \gamma_{-2} & \gamma_{0} & \dots & \gamma_{p-2} \\ \vdots & \vdots & \ddots & \vdots \\ \gamma_{1-p} & \gamma_{2-p} & \dots & \gamma_{0} \end{pmatrix}, \]

\[ \phi = \begin{pmatrix} \phi_1 \\ \phi_2 \\ \vdots \\ \phi_p \end{pmatrix}, \]

et

\[ \gamma_p = \begin{pmatrix} \gamma_1 \\ \gamma_2 \\ \vdots \\ \gamma_p \end{pmatrix}. \]

Notons que si l’on connaît la matrice des autocovariances, ces équations forment un système linéaire dont la solution est \(\phi\) et \(\sigma_z^2\).

Fonction d’autocorrélation d’un MA

Considérée comme un MA(1), la série \(\{y_t\}\) s’écrit comme \(y_t=z_t+\theta\cdot z_{t-1}\).

Sa variance est

\[ Var(y_t)=(1+\theta^2)\cdot\sigma_z^2 \]

\[ \rho(l) = \begin{cases} 1 & \text{si }l=0 \\ \frac{\theta}{1+\theta^2} & \text{si }l=1 \\ 0 & \text{si }l\gt 1 \end{cases} \]

La covariance est

\[ \begin{eqnarray*} Cov(y_t,y_{t-1}) &=& Cov(z_t+\theta\cdot z_{t-1},\, z_{t-1}+\theta\cdot z_{t-2}) \\ &=& \theta\cdot\sigma_z^2 \end{eqnarray*} \]

La fonction réelle \(x \mapsto \frac{x}{1+x^2}\) prend ses valeurs entre \(-0,5\) et \(+0,5\). Par conséquent il n’est pas possible de décrire des phénomènes à forte autocorrélation à l’aide d’un processus MA(1).

Quant à la fonction d’autocovariance (ACF) d’un MA(\(q\)) :

\[ \gamma_h = \begin{cases} (1 + \theta^2_1 + \theta^2_2 + \dots + \theta^2_q)\cdot\sigma_z^2 & \text{if } h = 0, \\ (\theta_h + \theta_{h+1}\theta_1 + \dots + \theta_q\theta_{q-h})\cdot\sigma_z^2 & \text{if } 1 \leq h \leq q, \\ 0 & \text{if } h > q. \end{cases} \]

En général, la fonction d’autocorrélation d’un processus MA(\(q\)) est nulle à partir de l’ordre \(q+1\). Si l’on observe une trajectoire d’un MA(\(q\)), on peut donc s’attendre à ce que la fonction d’autocorrélation empirique de la série ne soit pas significativement différente de 0 au-delà de l’ordre \(q\). Et si une fonction d’autocorrélation empirique est plus ou moins nulle à partir d’un certain ordre \(q+1\), on peut penser qu’il s’agit de la fonction d’autocorrélation empirique d’une série MA(\(q\)).

Voir : Formule de Bartlett

Prévision basée sur l’espérance conditionnelle

Nous suivons ici Hamilton, 1994, pp.72

Supposons que nous souhaitions prévoir la valeur d’une variable \(y_{t+1}\) sur la base d’un ensemble de variables \(M_t\), observées à la date \(t\). Par exemple, nous pourrions vouloir prévoir \(y_{t+1}\) sur la base de ses \(m\) valeurs les plus récentes. Dans ce cas, \(M\) se compose d’une constante et des variables \(y_{t},\, y_{t-1},\, \ldots,\, y_{t-m+1}\).

Soit \(y^*_{t+1\,|t}\) une prévision de \(y_{t+1}\) basée sur \(M_t\). Pour évaluer l’utilité de cette prévision, nous devons spécifier une fonction de perte. Nous choisissons une fonction de perte quadratique. Cela signifie que nous choisissons la prévision \(y^*_{t+1\,|t}\) de manière à minimiser

\[ E\left[\left(y_{t+1}-y^*_{t+1\,|t}\right)^2\right]. \]

Cette expression est connue sous le nom d’erreur quadratique moyenne associée à la prévision \(y^*_{t+1\,|t}\), et notée \(MSE(y^*_{t+1\,|t})\).

La prévision présentant l’erreur quadratique moyenne la plus faible est l’espérance de \(y_{t+1}\) conditionnellement à \(M_t\) :

\[ y^*_{t+1\,|t}=E\left[y_{t+1}\,|M_t\right]. \]

Pour vérifier cette affirmation, considérons que l’on base \(y^*_{t+1\,|t}\) sur une fonction \(g(M_t)\) autre que l’espérance conditionnelle,

\[ y^*_{t+1\,|t}=g(M_t). \]

Pour cette règle de prévision candidate, le \(MSE\) serait le suivant

\[ \begin{eqnarray*} E\left[\left(y_{t+1}-g(M_t)\right)^2\right] &=& E\left[\left(y_{t+1}-E\left[y_{t+1}\,|M_t\right]+E\left[y_{t+1}\,|M_t\right]-g(M_t)\right)^2\right] \\ &=& \quad\; E\left[\left(y_{t+1}-E\left[y_{t+1}\,|M_t\right]\right)^2\right] \\ &+& 2\cdot E\left[\left(y_{t+1}-E\left[y_{t+1}\,|M_t\right]\right)\cdot\left(E\left[y_{t+1}\,|M_t\right]-g(M_t)\right)\right] \\ &+& \quad\; E\left[\left(E\left[y_{t+1}\,|M_t\right]-g(M_t)\right)^2\right] \\ \end{eqnarray*} \]

Ecrivons le terme central du côté droit de cette équation de la façon suivante :

\[ 2\cdot E\left[\eta_{t+1}\right], \]

où

\[ \eta_{t+1}:=\left(y_{t+1}-E\left[y_{t+1}\,|M_t\right]\right)\cdot\left(E\left[y_{t+1}\,|M_t\right]-g(M_t)\right)\,. \]

Considérons d’abord l’espérance de \(\eta_{t+1}\) conditionnellement à \(M_t\). Sous condition de \(M_t\), les termes \(E\left[y_{t+1}\,|M_t\right]\) et \(g(M_t)\) sont des constantes connues et peuvent être factorisés hors de cette espérance :

\[ \begin{eqnarray*} E\left[\eta_{t+1}\,|M_t\right] &=& \left(E[y_{t+1}\,|M_t]-g(M_t)\right)\cdot\left(E\left[y_{t+1}-E[y_{t+1}\,|M_t]\;|M_t\right]\right) \\ &=& \left(E[y_{t+1}\,|M_t]-g(M_t)\right)\cdot 0 \\ &=& 0. \end{eqnarray*} \]

Par une application directe de la loi des espérances itérées, il s’ensuit que

\[ E\left[\eta_{t+1}\right]=E_{M_t}\left(\eta_{t+1}\,|M_t\right)=0. \]

Si l’on substitue ce résultat à l’expression ci-dessus de l’erreur quadratique moyenne, on obtient

\[ \begin{eqnarray*} E\left[\left(y_{t+1}-g(M_t)\right)^2\right] &=& E\left[\left(y_{t+1}-E\left[y_{t+1}\,|M_t\right]+E\left[y_{t+1}\,|M_t\right]-g(M_t)\right)^2\right] \\ &=& E\left[\left(y_{t+1}-E\left[y_{t+1}\,|M_t\right]\right)^2\right] \\ &+& E\left[\left(E\left[y_{t+1}\,|M_t\right]-g(M_t)\right)^2\right] \\ \end{eqnarray*} \]

Le deuxième terme du côté droit ne peut pas être rendu plus petit que zéro, et le premier terme ne dépend pas de \(g(M_t)\). La fonction \(g(M_t)\) qui rend l’erreur quadratique moyenne aussi petite que possible est la fonction qui met le second terme à zéro :

\[ E[y_{t+1}\,|M_t]=g(M_t). \]

Ainsi, la prévision \(g(M_t)\) qui minimise l’erreur quadratique moyenne est l’espérance conditionnelle \(E[y_{t+1}\,|M_t]\), comme on l’a dit.

L’erreur quadratique moyenne de cette prévision optimale est la suivante

\[ E\left[\left(y_{t+1}-g(M_t)\right)^2\right] = E\left[\left(y_{t+1}-E\left[y_{t+1}\,|M_t\right]\right)^2\right] \]

Nous allons maintenant restreindre la classe des prévisions considérées en exigeant que la prévision \(y^*_{t+1\,|t}\) soit une fonction linéaire de \(M_t\) :

\[ \begin{eqnarray*} y^*_{t+1\,|t} &=& \alpha^T\cdot M_t \\ &=& (1,\,\alpha_1,\,\ldots\,,\alpha_t)\cdot (\text{const},\, y_1,\,\ldots,\,y_t)^T \\ &=& \text{const} +\sum_{i=1}^t\alpha_i\cdot y_i. \end{eqnarray*} \]

Supposons que nous trouvions une valeur pour \(\alpha\) telle que l’erreur de prévision \(y_{t+1}-\alpha^T\cdot M_t\) ne soit pas corrélée avec \(M_t\) :

\[ E\left[\left(y_{t+1}-\alpha^T\cdot M_t\right)\cdot M_t^T\right]=0^T. \]

Si cette relation se vérifie, la prévision \(\alpha^T\cdot M_t\) est appelée projection linéaire de \(y_{t+1}\) sur \(M_t\).

La projection linéaire produit l’erreur quadratique moyenne la plus faible parmi la classe des règles de prévision linéaires. C’est le BLP (best linear predictor).

Affirmation : Soit une trajectoire d’une série stationnaire de moyenne \(\mu\), trajectoire que nous noterons \(y_1,\,\ldots,\,y_T\) en invitant le lecteur à ne pas confondre la variable aléatoire de l’une de ses réalisations possibles. Alors le BLP vérifie la relation

\[ y^*_{t+1\,|t}=\mu+\sum_{i=1}^T \alpha_i\cdot(y_{T+1-i}-\mu). \]

Ceci signifie que l’espérance du BLP est \(\mu\), et que par conséquent le BLP d’une série stationnaire est un estimateur sans biais de \(y_{t+1}\).

Remarque : Si \(y_t\) est gaussien, le BLP est aussi le meilleur prédicteur tout court.

Fonction d’autocorrélation partielle (PACF)

Considérons \(\{y_t\}\) centrée et ses régressions linéaires sur son passé résumées à la dernière, respectivement aux deux, trois dernières, etc., observations :

\[ \begin{eqnarray*} y_t &=& \phi_{1,1}\cdot y_{t-1}+u_{1,t} \\ y_t &=& \phi_{1,2}\cdot y_{t-1}+\phi_{2,2}\cdot y_{t-2}+u_{2,t} \\ y_t &=& \phi_{1,3}\cdot y_{t-1}+\phi_{2,3}\cdot y_{t-2}+\phi_{3,3}\cdot y_{t-3}+u_{3,t} \\ &\vdots& \end{eqnarray*} \]

Par exemple, \(\phi_{1,2}\cdot y_{t-1}+\phi_{2,2}\cdot y_{t-2}\) désigne le BLP de \(y_t\) connaissant \(y_{t-1}\) et \(y_{t-2}\) et s’obtient en résolvant l’équation

\[ \begin{pmatrix} \gamma_0 & \gamma_1 \\ \gamma_1 & \gamma_0 \\ \end{pmatrix} \cdot \begin{pmatrix} \phi_{1,2} \\ \phi_{2,2} \\ \end{pmatrix} = \begin{pmatrix} \gamma_1 \\ \gamma_2 \\ \end{pmatrix} . \]

Considérons les \(\phi_{k,k}\), \(k=1,2,\ldots\). Ils admettent la même interprétation que les coefficients d’une régression linéaire classique : \(\phi_{k,k}\) représente l’apport d’explication de \(y_{t-k}\) à \(y_t\), ceteris paribus, i.e. étant donné que l’on régresse également sur \(y_{t-1},\,\ldots,\,y_{t-k+1}\). Vus comme fonction de \(k\), les \(\phi_{k,k}\) définissent la fonction d’autocorrélation \(k\mapsto \phi_{k,k}\) partielle (PACF).

Supposons en particulier que \(y_t\) soit autorégressif AR(3), pour fixer les idées. Alors il est clair que \(y_{t-4}\) n’apporte rien de plus que \(y_t-1\), \(y_t-2\) et \(y_t-3\) quant il s’agit d’expliquer \(y_t\), et on montre qu’en effet pour un AR(3), \(\phi_{k,k}=0\) pour \(k>3\). De façon générale :

Affirmation : La PACF d’un AR(\(p\)) est nulle à partir de l’ordre \(p+1\).

Affirmation : La PACF d’un processus qui en plus a une composante moyenne mobile (un ARMA, donc) a un décroissance exponentielle ou sinusoïdale amortie exponentiellement.

Calcul de la PACF : Voir l’algorithme de Durbin-Levinson, qui calcule la PACF itérativement. L’algorithme fournit alors l’estimateur \(\hat{\phi}_{k,k}\).

D’un point de vue pratique, on pensera qu’une série temporelle suit un modèle AR(\(p\)) si les estimateurs \(\hat{\phi}_{k,k}\approx 0\) pour \(k\gt p\).

\[ \begin{eqnarray} \hat{\phi}_{k,k} &\longrightarrow& \phi_{k,k} \qquad &\text{pour }T\rightarrow\infty&\\ \hat{\phi}_{l,l} &\longrightarrow& 0 &\text{pour }T\rightarrow\infty\text{ et pour tout }l\gt p& \\ Var\left(\hat{\phi}_{l,l}\right) &\approx& 1/T &\forall l\gt p.& \end{eqnarray} \]

Si donc la PACF empirique d’une série n’est plus significativement différente de 0 à partir d’un certain ordre \(k\), on essayera de lui ajuster un modèle AR(\(k-1\)).

Prévision d’un modèle AR(\(p\))

Prévision d’un modèle MA(\(q\))

Nous ne le faisons ici que pour l’horizon 1.

\[ y_{t+1}=\mu+z_{t+1}+\theta_1\cdot z_t+\ldots+\theta_q\cdot z_{t+1-q} \]

et donc la prédiction BLP est

\[ \begin{eqnarray*} y_{t+1\,|t} &=& E\left[y_{t+1}\,|\,y_1,\,y_2,\,\ldots,\,y_t\right] \\ &=& \mu+\theta_1\cdot z_t+\theta_2\cdot z_{t-1}+\ldots+\theta_q\cdot z_{t-q} \end{eqnarray*} \]

et l’erreur de prédiction est \(e_t(1)=z_{t+1}\) (comme ci-dessus).

Notons juste que pour un horizon \(h>q\), la prévision est \(\mu\) (sans preuve ici).

fonction de vraisemblance d’un processus gaussien AR(1)

La série temporelle est

\[ Y_t=c+\phi\cdot Y_{t-1}+Z_t\,, \qquad Z_t\sim BB(0,\sigma^2). \]

On sait qu’alors \(Y_t\) suit une loi normale avec \(E[Y_1]=\frac{c}{1-\phi}\) et \(Var(Y_1)=\frac{\sigma^2}{1-\phi^2}\)

La fonction de densité de probabilité (fdp) de \(Y_1\) est par conséquent égale à

\[ f_{Y_1}(y_1|(c,\phi,\sigma^2))=\frac{1}{\sqrt{2\pi\frac{\sigma^2}{1-\phi^2}}}\cdot\exp\left(-\frac{(y_1-\frac{c}{1-\phi})^2}{2\frac{\sigma^2}{1-\phi^2}}\right). \]

Sachant que \(Y_1=y_1\), \(Y_2\) est de loi normale de moyenne \(c+\phi\cdot y_1\) et de variance \(\sigma^2\), et donc la densité conditionnelle associée est

\[ f_{Y_2\,|Y_1=y_1}(y_2|(c,\phi,\sigma^2))=\frac{1}{\sqrt{2\pi\sigma^2}}\cdot\exp\left(-\frac{(y_2-c-\phi\cdot y_1)^2}{2\sigma^2}\right). \]

On en déduit la f.d.p. conjointe au couple \((Y_1,Y_2)\) :

\[ f_{Y_1,Y_2}(y1,y2|(c,\phi,\sigma^2))=f_{Y_1}(y_1|(c,\phi,\sigma^2))\cdot f_{Y_2\,|Y_1=y_1}(y_2|(c,\phi,\sigma^2)) \]

On observe par ailleurs que \(Y_t\) ne dépend explicitement que de \(y_{t-1}\) :

\[ \begin{eqnarray*} f_{Y_t\,|Y_{t-1}=y_{t-1},\cdots,Y_1=y_1}(y_t|(c,\phi,\sigma^2)) &=& f_{Y_t\,|Y_{t-1}=y_{t-1}}(y_t|(c,\phi,\sigma^2)) \\ &=& \frac{1}{\sqrt{2\pi\sigma^2}}\cdot\exp\left(-\frac{(y_t-c-\phi\cdot y_{t-1})^2}{2\sigma^2}\right). \end{eqnarray*} \]

La f.d.p. conjointe des observations est par conséquent :

\[ f_{Y_1,\cdots, Y_T}(y_1,\cdots,y_T|(c,\phi,\sigma^2))=f_{Y_1}(y_1|(c,\phi,\sigma^2))\cdot\prod_{t=2}^T f_{Y_t\,|Y_{t-1}=y_{t-1}}(y_t|(c,\phi,\sigma^2))\;. \]

Ceci nous permet d’obtenir la fonction de \(\log-\)vraisemblance :

\[ \log\mathcal{L}(c,\phi,\sigma^2)=-\frac{1}{2}\ln\left(2\pi\frac{\sigma^2}{1-\phi^2}\right)-\frac{1}{2}\frac{(y_1-\frac{c}{1-\phi})^2}{\frac{\sigma^2}{1-\phi^2}}-\frac{T-1}{2}\ln(2\pi\sigma^2)-\frac{1}{2}\sum_{t=2}^T\frac{(y_t-c-\phi y_{t-1})^2}{\sigma^2}\;. \]

Si l’on travaille conditionnellement à la première valeur \(y_1\), la \(\log-\)vraisemblance se simplifie en la \(\log-\)vraisemblance conditionnelle :

\[ \log\mathcal{L}(c,\phi,\sigma^2)=-\frac{T-1}{2}\ln(2\pi\sigma^2)-\sum_{t=2}^T\frac{(y_t-c-\phi y_{t-1})^2}{\sigma^2}\;. \]

Sa maximisation donne l’estimateur du maximum de vraisemblance conditionnelle. Les estimateurs \(\hat{c}\), \(\hat{\phi}\) et \(\hat{\sigma}^2\), résultats de cette maximisation, sont obtenus à l’aide de la fonction arima().

Un raisonnement analogue permet de déterminer la fonction de vraisemblance d’un processus gaussien MA(1). Dans ce cas encore, la fonction arima() fournit les estimateurs.

Remarque : Les résidus \(u_t:=y_t-\hat{y}_t\) font sens ici comme en régression linéaire et si le modèle -disons AR(\(p\))- est correct, alors les résidus sont un bruit blanc et passent le test du portemanteau pour tous les retards raisonnables, i.e. entre 1 et environ un quart de la longueur \(T\) de la série des observations.

Modèle connu, paramètres à déterminer

On simule une série AR(1) vérifiant \(y_t=-0.7\cdot y_{t-1}+z_t\) avec \(z_t\sim BB(0,4)\). On commence par tester la blancheur de cette série, évidemment non blanche par construction.

La simulation :

require(cashchrono)

## Le chargement a nécessité le package : cashchrono

## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : aucun package nommé 'cashchrono' n'est trouvé

set.seed(123) # afin de pouvoir reproduire les résultats aléatoires
y=arima.sim(n=100,list(ar=-0.7),sd=sqrt(4))

Le test de blancheur :

ret=c(3,6,9,12)
# les différentes valeurs de retard choisies pour le test du portemanteau
Box.test.2(y,nlag=ret,type="Ljung-Box",decim=2)

##      Retard p-value
## [1,]      3       0
## [2,]      6       0
## [3,]      9       0
## [4,]     12       0

On en conclut que l’hypothèse de blancheur \(H_0\) doit être rejetée absolument car la \(p\)-value vaut 0 pour tous les retards.

Nous feignons ensuite d’oublier comment nous avons obtenu cette série et la considérons comme une série dont nous ne savons rien sinon que c’est un AR(1). Nous allons donc essayer de retrouver les justes paramètres. Nous devrions sans mystère retrouver les paramètres ayant servi à la simulation de notre série.

plot.ts(y)

Cette série est de moyenne nulle; ceci, il vaut mieux le préciser à la fonction arima() par l’option include.mean=FALSE. Par ailleurs les arguments de arima() sont \(p\), \(d\), \(q\), c’est-à-dire l’ordre d’autorégression, l’ordre de différenciation et l’ordre de moyenne mobile. On précise à cette fonction le modèle, ici AR(1), et celle-ci fait le fitting des paramètres du modèle donné.

my=arima(y,order=c(1,0,0), include.mean=FALSE)
my

## 
## Call:
## arima(x = y, order = c(1, 0, 0), include.mean = FALSE)
## 
## Coefficients:
##           ar1
##       -0.6324
## s.e.   0.0776
## 
## sigma^2 estimated as 3.05:  log likelihood = -197.91,  aic = 399.82

La moyenne est ici absente, comme nous l’avons demandé.

La commande summary(my) donne quelques résultats complémentaires : les valeurs ajustées, les résidus, la matrice de covariance des estimateurs. Quant à str(my), cette commande nous renseigne sur l’emplacement de ces résultats et la façon de les récupérer. On peut calculer des t-statistiques approchées en effectuant le quotient, paramètre estimé / estimation de l’écart type de l’estimateur. La fonction t_stat() de caschrono effectue ces quotients et donne les \(p\)-values approchées. Mais ces \(p\)-values n’ont de sens que si le résidu peut être considéré comme un bruit blanc. Il faut donc, avant de commenter les \(p\)-values des paramètres, tester la blancheur du résidu.

residuals(my)

## Time Series:
## Start = 1 
## End = 100 
## Frequency = 1 
##   [1]  2.577502263 -1.170373130 -1.914415109 -0.446538248 -2.015146304
##   [6] -1.344962301 -1.230571789 -3.302591977  1.853907812  0.068718804
##  [11] -2.130376990  2.559276273  0.647372382 -0.503873517  1.769729963
##  [16]  1.649691052  1.599120256  1.297105614  1.070915205 -0.172819440
##  [21] -0.569263192 -0.749466729 -1.346041022 -0.352333404 -2.547151702
##  [26]  4.520331562  1.995181533 -2.114905790 -0.745943874 -0.920754794
##  [31]  1.614191135 -0.310102238  0.618255938 -0.169452738 -0.003232489
##  [36]  2.685240919 -0.600079954  3.167421827 -3.396533141  1.587799536
##  [41] -0.124279393  0.675540645  0.559542737 -0.916124826 -0.660511051
##  [46] -1.996263552 -2.034583275  0.675568105  0.807452075  0.107727844
##  [51]  1.836169935  3.981416854 -1.175923534 -4.416291444  2.182037592
##  [56] -1.673679739 -1.101502015  1.951939274 -0.638669047 -2.354573478
##  [61]  0.466735376 -0.375168620  0.098464180  0.708926935 -0.750231783
##  [66]  1.345071363 -0.567457946  0.781893353  2.066020202  0.811528993
##  [71] -0.669476967  2.353981673  1.792335236  1.098831850  0.401942732
##  [76] -1.235202799  2.791715053 -1.433644714  4.618955236  2.598687126
##  [81] -0.351898616 -2.104647679 -1.245868773  0.487289755 -0.509557079
##  [86] -0.650433211 -1.887536872  0.027528026 -1.646033653 -3.176477739
##  [91] -0.646680868  1.809725292 -1.255072467  1.366729039 -3.423468055
##  [96]  0.238860856  0.801332062  0.698367000  0.103421214 -1.220138158

Pour accéder aux analyses :

require(broom)

## Le chargement a nécessité le package : broom

# Le package broom permet de mettre l'output de summary()
# sous forme de tableau data.frame, facile à manipuler.
library(broom)
tidy(my)

## # A tibble: 1 × 3
##   term  estimate std.error
##   <chr>    <dbl>     <dbl>
## 1 ar1     -0.632    0.0776

Quant à voir comment est structuré le data.frame ainsi produit :

str(my)

## List of 14
##  $ coef     : Named num -0.632
##   ..- attr(*, "names")= chr "ar1"
##  $ sigma2   : num 3.05
##  $ var.coef : num [1, 1] 0.00601
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "ar1"
##   .. ..$ : chr "ar1"
##  $ mask     : logi TRUE
##  $ loglik   : num -198
##  $ aic      : num 400
##  $ arma     : int [1:7] 1 0 0 0 1 0 0
##  $ residuals: Time-Series [1:100] from 1 to 100: 2.578 -1.17 -1.914 -0.447 -2.015 ...
##  $ call     : language arima(x = y, order = c(1, 0, 0), include.mean = FALSE)
##  $ series   : chr "y"
##  $ code     : int 0
##  $ n.cond   : int 0
##  $ nobs     : int 100
##  $ model    :List of 10
##   ..$ phi  : num -0.632
##   ..$ theta: num(0) 
##   ..$ Delta: num(0) 
##   ..$ Z    : num 1
##   ..$ a    : num -0.646
##   ..$ P    : num [1, 1] 0
##   ..$ T    : num [1, 1] -0.632
##   ..$ V    : num [1, 1] 1
##   ..$ h    : num 0
##   ..$ Pn   : num [1, 1] 1
##  - attr(*, "class")= chr "Arima"

On poursuit en testant la blancheur des résidus du modèle AR(1) obtenu.

ret=c(3,6,9,12)
Box.test.2(residuals(my),nlag=ret,type="Ljung-Box",decim=4,fitdf=1)

##      Retard p-value
## [1,]      3  0.8164
## [2,]      6  0.9647
## [3,]      9  0.9889
## [4,]     12  0.9812

La commande supplémentaire fitdf=1 indique le nombre de degrés de liberté devant être soustraits, ici 1 en raison du fait que la série temporelle est formée de résidus résultant d’un ajustement à un paramètre.

Clairement, on accepte \(H_0\) et concluons que la série des résidus est blanche. Ceci démontre la qualité du fit AR(1) : le modèle que la fonction arima() propose explique l’autocorrélation de la série \(\{y\}\) puisque la série des résidus ne montre plus aucune autocorrélation non expliquée.

D’autre part et comme nous l’avons dit plus haut, le fait que la série des résidus est blanche autorise à calculer la \(p\)-value pour chaque paramètre estimé, puisqu’alors elle fait bien sens.

Dans notre cas AR(1) à moyenne nulle, il n’y a qu’un seul paramètre et son estimation est obtenue par

my$coef

##        ar1 
## -0.6324448

et sa variance obtenue par

my$var.coef

##             ar1
## ar1 0.006014563

La statistique de Student est le quotient

(my$coef)/(my$var.coef)^0.5

##           ar1
## ar1 -8.154936

La \(p\)-value est la probabilité qu’une v.a. de loi de Student, i.e. plus ou moins une loi \(\mathcal{N}(0,1)\), dépasse ce quotient en valeur absolue. Comme cette probabilité est très faible, on refuse l’hypothèse que le coefficient d’autorégression soit non significatif.

On peut aussi utiliser t_stat() qui fait exactement ce travail :

t_stat(my)

##              ar1
## t.stat -8.154936
## p.val   0.000000

Modèle et paramètres inconnus

Le modèle est supposé cette fois si général, à savoir un ARMA(\(p\),\(q\)), que nous pouvons presque le considérer comme inconnu (si ce n’est qu’il est stationnaire). En effet, même un modèle saisonnier SARMA peut-être remplacé par un ARMA, au paramétrage certes plus compliqué mais tout de même. Savoir qu’une série observée suit un modèle ARMA ne nous avance donc pas beaucoup quant à la détermination exacte de sa modélisation, paramètres y compris.

Nous nous donnons une trajectoire \(\{y_t\}_{t=1:T}\), et nous jugeons, après visualisation de son chronogramme, qu’il s’agit d’une série stationnaire. Nous voulons par conséquent lui ajuster un modèle ARMA(\(p\),\(q\)).

La première étape consiste à choisir les ordres \(p\) et \(q\). C’est ce qu’on appelle l’étape d’identification.

Vient ensuite l’étape d’estimation. Il s’agit alors d’estimer le modèle pour confirmer ses choix et éventuellement les améliorer. Le choix des ordres \(p\) et \(q\) est en effet rarement unique, et pour chaque choix, on doit estimer le modèle et en tester la qualité. On mesure la qualité d’un modèle à la blancheur de son résidu. Si le résidu n’est pas blanc, c’est que le modèle ne capte pas toute la dynamique du phénomène observé.

Enfin et si le modèle retenu est satisfaisant, on peut l’utiliser pour prédire la série, étape de prédiction.

La méthode MINIC

MINIC est un acronyme pour Minimum Information Criterion. Cette méthode aide à découvrir les ordres \(p\) et \(q\) pour une série stationnaire susceptible de recevoir une modélisation ARMA.

On dispose d’une série \(y_t\), \(t = 1, \dots, n\), centrée, observation d’une série de moyenne nulle, stationnaire et inversible suivant un modèle ARMA, d’ordres \(p\) et \(q\) inconnus :

\[y_t = \frac{1 + \theta_1 B + \dots + \theta_q B^q}{1 - \phi_1 B - \dots - \phi_p B^p} z_t\]

où \(z_t \sim BBN(0, \sigma_z^2)\).

On a noté qu’un modèle ARMA inversible peut être représenté comme un AR(\(\infty\)) ; en choisissant un ordre d’autorégression suffisamment grand, on peut approcher de façon satisfaisante le modèle ARMA de la série par un modèle AR d’ordre fini. Le résidu de l’ajustement de la série à un tel modèle est alors proche du bruit blanc sous-jacent à la série, même si les paramètres trop nombreux sont mal estimés. Par ailleurs, l’ajustement d’un modèle autorégressif à une série peut se faire rapidement par différentes méthodes, MCO et Yule-Walker notamment. On obtient ainsi un substitut de résidu, proche du résidu qu’on aurait obtenu en régressant avec le modèle correct. Une fois un tel résidu disponible, il est facile d’ajuster un modèle ARMA en considérant les résidus retardés comme des covariables.

Par exemple, supposons qu’on ait obtenu des résidus \(z_t\) après une régression AR d’ordre élevé. On veut ajuster un ARMA(2, 3) à la série. On peut effectuer l’ajustement linéaire donc rapide :

\[y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + \theta_1 z_{t-1} + \theta_2 z_{t-2} + \theta_3 z_{t-3} + z_t, \quad t = 4, \dots, T,\]

pour obtenir des estimations des \(\phi\) et \(\theta\).

On voit que les ajustements d’ARMA de différents ordres \(p\) et \(q\) sont ainsi facilement effectués ; ils doivent ensuite être comparés. Comme ces modèles ne sont pas emboîtés, on doit choisir d’après un critère d’information. On retient le modèle correspondant au minimum de ce critère.

Calcul de la PACF

Le logiciel calcule la PACF sur la base de l’algorithme de Durbin-Levinson. Il s’agit d’un calcul itératif de la PACF à partir de la fonction d’autocorrélation (ACF) dont le graphe est le corrélogramme. On emploie le package FitAR.

On simule la série \(\{y\}\) suivante :

\[ \begin{array}{ll} y_t&:=&4+\frac{(1+0.6\cdot B^2)z_t}{1+0.8\cdot B-0.2\cdot B^4}\,, \\ z_t&\sim& \mathcal{N}(0,1.5) \end{array} \]

set.seed(281)
y=4+arima.sim(n=200,list(ar=c(0,-0.8,0,0,0.2),ma=c(0,0.6)),
sd = sqrt(1.5),n.start=50)
plot.ts(y)

Le corrélogramme est

acf(y)

Les deux lignes horizontales indiquent l’intervalle de confiance à 95% autour de 0.

Lire attentivement la description des options des fonctions acf() et pacf() dans l’aide, pour savoir comment elles gèrent par défaut les éventuelles valeurs manquantes.

La PACF est :

pacf(y)

Les ACF et PACF théoriques sont :

acf.th=ARMAacf(ar=c(0,-0.8,0,0,0.2),ma=c(0,0.6),lag.max=20,pacf=FALSE)