Salem Nkunda Nyisingize

Séries temporelles multivariées : quand une série en regarde une autre

Neuf chapitres à considérer les séries « une à la fois », puis la question inévitable : que se passe-t-il quand elles réagissent les unes aux autres ? Corrélations croisées, processus VAR, condition de stationnarité, et le piège du préblanchiment. Avec un laboratoire VAR(1) interactif et du code Python exécutable dans la page.

Pourquoi c'est important

Voici deux séries temporelles. Leur corrélation vaut \(-0{,}02\) : rigoureusement rien. Vous les superposez, vous calculez un coefficient de Pearson, vous concluez qu'elles n'ont aucun rapport, et vous passez à autre chose.

Vous venez de manquer une dépendance de \(-0{,}76\). Elle était là, entière, mais décalée de trois pas de temps. La première série ne réagit pas à la seconde au même instant — elle y réagit trois périodes plus tard. Aucun coefficient de corrélation ne peut voir ça, parce qu'un coefficient de corrélation ne regarde qu'un seul alignement : celui de maintenant.

C'est ce que le cadre vectoriel rend visible. Là où l'analyse univariée traite chaque série « une à la fois », le modèle VAR laisse chacune écouter le passé des autres. Trois taux d'intérêt qui réagissent ensemble aux annonces de la banque centrale, deux capteurs qui subissent le même choc thermique, une variable économique qui en entraîne une autre avec un délai : dans tous ces cas, l'information utile est précisément celle que l'approche série-par-série jette à la poubelle.

Et il y a un piège symétrique, plus vicieux encore. Deux séries lisses et sans aucun lien paraissent corrélées, à presque tous les décalages. Vous verrez ici pourquoi — et le remède, le préblanchiment, tient en quatre lignes. Deux laboratoires vous laissent déformer un modèle VAR à la main, regarder ses racines franchir le cercle unité, et faire s'effondrer des corrélations fantômes d'un clic.

VAR(p) Corrélations croisées Yule–Walker multivarié Préblanchiment Haugh (1976)

Le vocabulaire vectoriel

Un processus stochastique multivarié \(\{\mathbf{X}_t;\, t \in T\}\) est simplement une collection de variables aléatoires vectorielles. Chaque composante \(\{X_{ti}\}\) est une série temporelle univariée ordinaire, et on les empile :

Le vecteur d'observation à l'instant t
\[\mathbf{X}_t = (X_{t1},\, X_{t2},\, \dots,\, X_{tm})'\]

Pour les données de taux d'intérêt américains mensuels de janvier 1991 à avril 2010, \(X_{t1}\) désigne le taux des certificats de dépôt à six mois, \(X_{t2}\) celui des obligations Aaa de Moody's, et \(X_{t3}\) celui des hypothèques conventionnelles à trente ans. Les trois réalisations univariées forment une seule réalisation d'un processus trivarié.

La moyenne devient un vecteur, la covariance une matrice :

Moyenne et covariance
\[\boldsymbol\mu_t = E[\mathbf{X}_t] = (\mu_{t1}, \dots, \mu_{tm})', \qquad \Gamma(t_1,t_2) = E\big[(\mathbf{X}_{t_1} - \boldsymbol\mu_{t_1})(\mathbf{X}_{t_2} - \boldsymbol\mu_{t_2})'\big]\]

L'élément \(\gamma_{ij}(t_1,t_2)\) de cette matrice est la covariance entre \(X_{t_1 i}\) et \(X_{t_2 j}\) : deux séries différentes, à deux instants différents. C'est exactement l'information qu'une analyse univariée jette.

Stationnarité et corrélations croisées

Comme dans le cas univarié, on impose que les moments ne dépendent pas de l'origine du temps. Le processus est stationnaire en covariance si \(E[\mathbf{X}_t] = \boldsymbol\mu\) est constant, si \(\Gamma(0)\) existe et est constante, et si \(\Gamma(t_1,t_2)\) ne dépend que de l'écart \(t_2 - t_1\). On écrit alors :

Matrice de covariance au retard k
\[\Gamma(k) = \begin{pmatrix} \gamma_{11}(k) & \cdots & \gamma_{1m}(k) \\ \vdots & \ddots & \vdots \\ \gamma_{m1}(k) & \cdots & \gamma_{mm}(k)\end{pmatrix}, \qquad \gamma_{ij}(k) = E\big[(X_{ti}-\mu_i)(X_{t+k,j}-\mu_j)\big]\]

Quand \(i \neq j\), \(\gamma_{ij}(k)\) s'appelle la covariance croisée au retard \(k\). Les éléments diagonaux de \(\Gamma(0)\) sont les variances des composantes. La corrélation croisée normalise :

Corrélation croisée
\[\rho_{ij}(k) = \frac{\gamma_{ij}(k)}{\sqrt{\gamma_{ii}(0)\,\gamma_{jj}(0)}}\]

Deux propriétés méritent d'être mémorisées, parce qu'elles ne ressemblent pas à leurs analogues univariées. Dans le cas univarié, \(\gamma(k) = \gamma(-k)\) : l'autocovariance est paire. Ici, ce n'est pas vrai pour les termes croisés. On a seulement :

Symétrie croisée
\[\gamma_{ij}(k) = \gamma_{ji}(-k), \qquad \Gamma(-k) = \Gamma'(k)\]

Notez le changement d'ordre des indices. C'est précisément ce qui permet de distinguer « \(X_1\) précède \(X_2\) » de « \(X_2\) précède \(X_1\) » — l'asymétrie porte toute l'information directionnelle. Les valeurs \(\rho_{ij}(0)\) mesurent quant à elles la relation contemporaine, et \(\rho(0)\) porte des 1 sur sa diagonale.

La mise en garde qui fait perdre des heures. Deux définitions incompatibles de la covariance croisée coexistent dans la littérature. Celle utilisée ici, celle de Box et al. (2008), pose \(\gamma_{ij}(k) = E[(X_{ti}-\mu_i)(X_{t+k,j}-\mu_j)]\). Celle de Shumway et Stoffer, de Tsay, de Brockwell et Davis pose \(\gamma^{(a)}_{ij}(k) = E[(X_{t+k,i}-\mu_i)(X_{tj}-\mu_j)]\). Les deux ne sont pas équivalentes : elles inversent le signe du retard. En R, la fonction de base ccf(xi, xj) calcule la seconde ; pour obtenir la première, il faut appeler ccf(xj, xi). Vérifiez toujours quelle convention emploie la source avant d'interpréter un pic à un retard donné — sinon vous conclurez que la cause suit l'effet.

Le bruit blanc multivarié

La série \(\{\mathbf{a}_t\}\) est un bruit blanc multivarié si \(E[\mathbf{a}_t] = \mathbf{0}\) et \(\boldsymbol\rho_a(k) = \mathbf{0}\) pour tout \(k \neq 0\). Rien n'est imposé sur \(\boldsymbol\rho_a(0)\). Autrement dit, les composantes peuvent être corrélées de façon contemporaine : \(\rho_{ij}(0) \neq 0\) est parfaitement permis. On note \(\mathbf{a}_t \sim \text{MWN}(\mathbf{0}, \Gamma_a)\) avec \(\Gamma(0) = \Gamma_a\).

Cela paraît contradictoire au premier abord — du bruit blanc corrélé ? Mais les deux notions portent sur des axes différents. Le bruit blanc interdit toute structure le long de l'axe du temps ; il n'interdit rien en travers, entre composantes au même instant. Deux capteurs qui subissent le même choc thermique à chaque instant produisent exactement cela. L'exemple canonique du chapitre est :

Bruit blanc bivarié à corrélation contemporaine
\[\Gamma_a = \begin{pmatrix} 1 & 1{,}78 \\ 1{,}78 & 4 \end{pmatrix}, \qquad \boldsymbol\rho_a(0) = \begin{pmatrix} 1 & 0{,}89 \\ 0{,}89 & 1 \end{pmatrix}\]

Chaque réalisation prise isolément semble parfaitement non corrélée le long de l'axe du temps ; superposées, elles montent et descendent ensemble. Cette matrice \(\Gamma_a\) est celle qui alimente le VAR(1) du laboratoire ci-dessous.

Le processus VAR(p)

Le modèle central du chapitre. Chaque composante est régressée sur les \(p\) valeurs passées de toutes les composantes, y compris elle-même. En dimension deux et à l'ordre un :

Un VAR(1) bivarié, écrit ligne par ligne
\[X_{t1} = \phi_{11}X_{t-1,1} + \phi_{12}X_{t-1,2} + a_{t1}\] \[X_{t2} = \phi_{21}X_{t-1,1} + \phi_{22}X_{t-1,2} + a_{t2}\]

Ce qui se condense en \(\mathbf{X}_t = \Phi_1 \mathbf{X}_{t-1} + \mathbf{a}_t\). Le cas général :

Définition — processus VAR(p)
\[\mathbf{X}_t = \boldsymbol\alpha + \Phi_1\mathbf{X}_{t-1} + \cdots + \Phi_p\mathbf{X}_{t-p} + \mathbf{a}_t\]

où \(\mathbf{X}_t\) est \(m \times 1\), chaque \(\Phi_k\) est une matrice réelle \(m \times m\), \(\mathbf{a}_t\) est un bruit blanc multivarié de covariance \(\Gamma_a\), et \(\boldsymbol\alpha = (I - \Phi_1 - \cdots - \Phi_p)\boldsymbol\mu\).

En notation opérateur, \(\Phi(B)(\mathbf{X}_t - \boldsymbol\mu) = \mathbf{a}_t\) avec \(\Phi(B) = I - \Phi_1 B - \cdots - \Phi_p B^p\). Le premier vrai résultat concerne la stationnarité, et c'est une condition purement algébrique :

Théorème — condition de stationnarité et de causalité

Le processus est stationnaire et causal si et seulement si les racines de l'équation déterminantale

\[|\Phi(z)| = |I - \Phi_1 z - \cdots - \Phi_p z^p| = 0\]

se situent toutes à l'extérieur du cercle unité.

C'est l'analogue exact de la condition univariée sur \(\phi(z) = 0\). Pour un VAR(1), \(|\Phi(z)| = |I - \Phi_1 z|\), et les racines sont les inverses des valeurs propres de \(\Phi_1\) : la condition revient à exiger que toutes les valeurs propres de \(\Phi_1\) soient de module strictement inférieur à 1. Le laboratoire ci-dessous rend cette condition visible.

Sous stationnarité, le processus admet une représentation linéaire générale \(\mathbf{X}_t - \boldsymbol\mu = \Psi(B)\mathbf{a}_t = \sum_{j \geq 0}\Psi_j \mathbf{a}_{t-j}\) avec \(\Psi(B) = \Phi^{-1}(B)\), et les autocovariances satisfont les équations de Yule–Walker multivariées :

Yule–Walker multivarié
\[\Gamma(k) = \Gamma(k-1)\Phi_1' + \cdots + \Gamma(k-p)\Phi_p', \qquad k = 1,\dots,p\]

Et pour \(k = 0\), la relation \(\Gamma(0) - \Phi_1\Gamma(0)\Phi_1' = \Gamma_a\) — une équation de Lyapunov discrète — détermine la covariance stationnaire. Elle se résout en vectorisant : \(\operatorname{vec}\Gamma(0) = (I - \Phi_1 \otimes \Phi_1)^{-1}\operatorname{vec}\Gamma_a\).

Laboratoire : le VAR(1) de l'exemple 10.2

L'exemple traité dans le chapitre est un VAR(1) bivarié dont les paramètres se reconstituent entièrement à partir des covariances publiées. En appliquant \(\Phi_1' = \Gamma(0)^{-1}\Gamma(1)\) aux matrices du livre, on retrouve :

Les paramètres de l'exemple
\[\Phi_1 = \begin{pmatrix} 0{,}9 & 0{,}6 \\ -0{,}3 & 0{,}7 \end{pmatrix}, \qquad \Gamma_a = \begin{pmatrix} 1 & 1{,}78 \\ 1{,}78 & 4 \end{pmatrix}\]

Ces valeurs redonnent exactement l'équation caractéristique \(1 - 1{,}6z + 0{,}81z^2 = 0\) mentionnée dans le texte, dont les racines sont complexes de module \(1/\sqrt{0{,}81} = 1{,}11\) : à l'extérieur du cercle unité, donc stationnaire, et complexes, d'où le comportement pseudo-cyclique de période environ 13 visible sur les réalisations. Elles redonnent aussi, par l'équation de Lyapunov, la covariance stationnaire du livre :

Covariance et corrélations qui en découlent
\[\Gamma(0) = \begin{pmatrix} 27{,}74 & -0{,}40 \\ -0{,}40 & 13{,}07 \end{pmatrix}, \qquad \boldsymbol\rho(0) = \begin{pmatrix} 1 & -0{,}02 \\ -0{,}02 & 1 \end{pmatrix}, \qquad \boldsymbol\rho(3) = \begin{pmatrix} 0{,}26 & -0{,}76 \\ 0{,}72 & -0{,}06 \end{pmatrix}\]

Voilà le fait le plus instructif de tout le chapitre. La corrélation contemporaine entre les deux séries vaut \(-0{,}02\) : rigoureusement rien. Un analyste qui superpose les deux courbes et calcule un coefficient de Pearson conclut à l'indépendance. Pourtant \(\rho_{12}(3) = -0{,}76\) : décalez l'une de trois périodes, et la relation apparaît, massive. Une relation forte peut être totalement invisible au retard zéro.

Laboratoire 1 — VAR(1) bivarié : paramètres, racines, réalisations

Bougez les \(\phi_{ij}\) et regardez les racines franchir le cercle unité.

0.90 0.60 -0.30 0.70
200
À essayer. Mettez \(\phi_{12} = 0\) et \(\phi_{21} = 0\) : les deux séries deviennent deux AR(1) parfaitement séparés, et les corrélations croisées théoriques s'effondrent vers zéro à tous les retards — sauf au retard 0, où subsiste la corrélation contemporaine du bruit. Puis remontez \(\phi_{12}\) : la structure croisée réapparaît, asymétrique. Enfin poussez \(\phi_{11}\) vers 1,2 : une racine entre dans le cercle unité et le processus explose.

Un mot sur le spectre croisé

Le domaine fréquentiel se généralise de la même manière. Le spectre devient une matrice \(m \times m\) pour chaque fréquence :

Spectre multivarié
\[\mathbf{P}_X(f) = \sum_{k=-\infty}^{\infty} e^{-2\pi i f k}\,\Gamma(k), \qquad |f| \leq 0{,}5\]

Les éléments diagonaux \(P_{jj}(f)\) sont les spectres univariés habituels. Les éléments hors diagonale \(P_{jh}(f)\), appelés spectres croisés, sont complexes. On les décompose en \(P_{jh}(f) = c_{jh}(f) - i\,q_{jh}(f)\), où \(c_{jh}\) est le cospectre et \(q_{jh}\) le spectre de quadrature. En forme polaire \(P_{jh}(f) = \alpha_{jh}(f)e^{i\varphi_{jh}(f)}\), le module est le spectre d'amplitude et \(\varphi_{jh}(f) = \tan^{-1}(-q_{jh}(f)/c_{jh}(f))\) le spectre de phase — qui indique le déphasage entre les deux séries à chaque fréquence. Enfin la cohérence carrée :

Cohérence carrée
\[\kappa^2_{jh}(f) = \frac{|P_{jh}(f)|^2}{P_{jj}(f)P_{hh}(f)} \leq 1\]

Elle s'interprète comme un \(R^2\) fréquence par fréquence : la part de la variabilité de \(X_j\) à la fréquence \(f\) explicable linéairement par \(X_h\) à cette même fréquence. Deux séries peuvent être fortement cohérentes aux basses fréquences et pas du tout aux hautes — un diagnostic que la corrélation croisée seule ne donne pas.

Estimation et choix de l'ordre

Les estimateurs sont les analogues directs du cas univarié. La moyenne est estimée par le vecteur des moyennes empiriques \(\bar{\mathbf{X}}\), et les covariances par :

Covariance croisée empirique
\[\hat\gamma_{ij}(k) = \frac{1}{n}\sum_{t=1}^{n-k}(X_{ti} - \bar X_i)(X_{t+k,j} - \bar X_j)\]

Notez le \(1/n\) et non \(1/(n-k)\) : estimateur biaisé mais de variance moindre, et qui garantit une matrice définie non négative. Pour \(k = 0\), \(\hat\rho_{ij}(0)\) est simplement le coefficient de Pearson usuel calculé sur les \(n\) observations.

Pour les paramètres du VAR, deux voies : les équations de Yule–Walker multivariées, ou les moindres carrés équation par équation. Elles donnent des résultats très proches en pratique. L'ordre \(p\) se choisit par AIC — indispensable, puisqu'un VAR(\(p\)) en dimension \(m\) compte \(m^2 p\) coefficients : un VAR(3) trivarié en compte déjà 27.

# ─── R : base + package vars ───
# la fonction ar() de base calcule les estimations de Yule-Walker pour un VAR ;
# attention, les methodes mle et burg sont univariees uniquement
donnees <- cbind(x1, x2)
ar(donnees, order.max = 8, aic = TRUE, method = 'yw')
ar(donnees, order.max = 8, aic = TRUE, method = 'ols')

library(vars)
ajust <- VAR(donnees, p = 1, type = c("const"), ic = c("AIC"))
serial.test(ajust, lags.pt = 24, type = c("PT.asymptotic"))  # Ljung-Box multivarie
predict(ajust, n.ahead = 25)   # prevision + limites a 95 % par defaut

Le code ci-dessous s'exécute réellement dans votre navigateur (Python compilé en WebAssembly). Modifiez les valeurs et relancez : le premier lancement prend quelques secondes le temps de charger Python.

PYTHON · NUMPY
PYTHON · NUMPY

Tester une association : la procédure de Haugh

Question fréquente : ces deux séries sont-elles liées ? On aimerait tester \(H_0 : \rho_{12}(k) = 0\) en comparant \(\hat\rho_{12}(k)\) à des bandes du type \(\pm 2/\sqrt{n}\). Mais la variance asymptotique de \(\hat\rho_{12}(k)\) vaut :

Variance asymptotique de la corrélation croisée empirique
\[\mathrm{AV}_{12} = \frac{1}{n}\sum_{j=-\infty}^{\infty}\rho_{11}(j)\,\rho_{22}(j)\]

Regardez bien cette formule. Elle fait intervenir les autocorrélations des deux séries. Si les deux sont fortement autocorrélées, la somme est grande, et la variance réelle dépasse de loin \(1/n\). Les bandes \(\pm 2/\sqrt{n}\) sont alors beaucoup trop étroites : on voit des relations partout. La variance ne se réduit à \(1/n\) que si au moins une des deux séries est un bruit blanc.

D'où la procédure recommandée par Haugh (1976), en quatre étapes :

ÉtapeOpération
1Ajuster un ARMA(\(p_i, q_i\)) à chacune des deux séries — les ordres peuvent différer
2Calculer les séries de résidus \(\hat a_{t1}\) et \(\hat a_{t2}\) (c'est le préblanchiment)
3Calculer les corrélations croisées \(\hat\rho^{(a)}_{12}(k)\) entre les deux séries de résidus
4Déclarer \(\hat\rho^{(a)}_{12}(k)\) significative au seuil \(\alpha = 0{,}05\) si \(|\hat\rho^{(a)}_{12}(k)| > 2/\sqrt{n}\)

Les résidus sont du bruit blanc par construction, donc la variance asymptotique redevient \(1/n\) et les bandes usuelles retrouvent leur validité. Comme toujours avec des tests répétés, le seuil s'applique séparément à chaque \(k\) : sur 41 retards testés, deux dépassements fortuits sont attendus sous \(H_0\).

Laboratoire 2 — corrélations croisées avant et après préblanchiment

Deux sources de données au choix, bandes \(\pm 2/\sqrt{n}\) en pointillés.

+3 4 200
L'expérience décisive. Basculez la source sur « AR(1) indépendants ». Les deux séries sont engendrées par deux générateurs sans aucun lien, mais chacune est très autocorrélée. Sans préblanchiment, la corrélogramme croisé montre des pics amples qui traversent les bandes — une relation entièrement fictive. Activez le préblanchiment : les pics s'effondrent dans les bandes. C'est exactement le mécanisme qui, dans le chapitre, conduit à conclure qu'il n'y a pas d'association entre l'activité solaire et la température terrestre : la corrélation croisée brute atteint \(\hat\rho_{12}(6) = 0{,}271\), mais après ajustement d'un AR(9) sur les taches solaires et d'un AR(4) sur les températures, plus aucune corrélation résiduelle ne dépasse les limites \(\pm 2/\sqrt{160} = \pm 0{,}158\).
PYTHON · NUMPY

Le cas non stationnaire

Tout ce qui précède suppose la stationnarité. Quand elle manque, on généralise le modèle ARUMA au cas multivarié en insérant un opérateur diagonal \(\lambda(B)\) :

Modèle VARUMA
\[\Phi(B)\lambda(B)(\mathbf{X}_t - \boldsymbol\mu) = \Theta(B)\mathbf{a}_t, \qquad \lambda(B) = \operatorname{diag}\big(\lambda_1(B), \dots, \lambda_m(B)\big)\]

où chaque \(\lambda_j(z) = 0\) a toutes ses racines sur le cercle unité, tandis que celles de \(|\Phi(z)| = 0\) et \(|\Theta(z)| = 0\) restent à l'extérieur. Le cas qui a reçu le plus d'attention est celui où \(\lambda(B) = \operatorname{diag}\big((1-B)^{d_1}, \dots, (1-B)^{d_m}\big)\) — le pendant multivarié exact du modèle ARIMA(\(p,d,q\)). Si la composante \(j\) est déjà stationnaire, \(d_j = 0\) et aucune différenciation ne lui est appliquée.

Et ensuite ? Différencier chaque composante séparément est parfois la bonne réponse — et parfois une perte d'information sévère. Quand les composantes sont cointégrées, il existe une combinaison linéaire déjà stationnaire, et différencier la détruit. C'est le sujet de l'article suivant, qui reprend exactement ce cadre VAR pour aboutir au modèle à correction d'erreur.

Limites et pièges

Hypothèse ou piègeConséquence → parade
Nombre de paramètres en \(m^2p\)Surajustement rapide dès \(m \geq 4\) → AIC/BIC, VAR bayésien, ou restrictions de nullité
Convention de la covariance croiséeSigne du retard inversé selon la source → vérifier, et préférer ccf(xj, xi) en R pour la convention de ce texte
Lecture d'une CCF bruteCorrélations fantômes entre séries autocorrélées → préblanchir systématiquement
Corrélation contemporaine nulleNe signifie pas indépendance : voir \(\rho_{12}(0) = -0{,}02\) et \(\rho_{12}(3) = -0{,}76\) dans le même modèle
Composantes non stationnairesRégression fallacieuse → tester les racines unitaires, puis cointégration avant de différencier
Bruit blanc à corrélation contemporaine\(\Gamma_a\) non diagonale complique l'interprétation causale → décomposition de Cholesky et ordre des variables, qui n'est pas neutre
VARMA au-delà du VARProblèmes d'identifiabilité sérieux → en pratique, on reste presque toujours sur du VAR pur

Références

  • Woodward, W. A., Gray, H. L. & Elliott, A. C. (2017). Applied Time Series Analysis with R, 2e éd., CRC Press — chapitre 10, la source directe de cet article (définitions, exemple 10.2, procédure de test).
  • Haugh, L. D. (1976). Checking the independence of two covariance-stationary time series: a univariate residual cross-correlation approach, JASA 71 — la procédure de préblanchiment.
  • Brockwell, P. J. & Davis, R. A. (1991, 2002). Time Series: Theory and Methods — la théorie asymptotique des estimateurs multivariés.
  • Reinsel, G. C. (1997). Elements of Multivariate Time Series Analysis, 2e éd., Springer.
  • Box, G. E. P., Jenkins, G. M. & Reinsel, G. C. (2008). Time Series Analysis: Forecasting and Control, 4e éd. — la convention de covariance croisée retenue ici.
  • Lütkepohl, H. (2005). New Introduction to Multiple Time Series Analysis, Springer — la référence exhaustive sur les VAR.
  • Pfaff, B. (2008). VAR, SVAR and SVEC models: implementation within R package vars, Journal of Statistical Software 27(4).

Les simulateurs de cette page sont écrits en JavaScript/Canvas sans dépendance externe ; les cellules Python s'exécutent via Pyodide. Le code est disponible sur GitHub.