Multivariate time series: when one series watches another

Nine chapters spent taking series "one at a time", then the unavoidable question: what happens when they react to each other? Cross-correlations, VAR processes, the stationarity condition, and the prewhitening trap. With an interactive VAR(1) lab and Python code that runs in the page.

📅 August 2026 ⏱️ ~20 min read 🏷️ Time series • VAR • Cross-correlations

← Back to the blog Jump to the VAR lab 🇫🇷 Lire en français

Contents

Why this matters
The vector vocabulary
Stationarity and cross-correlations
Multivariate white noise
The VAR(p) process
Lab: the VAR(1) of Example 10.2
A word on the cross-spectrum
Estimation and order selection
Testing for association: the Haugh procedure
The non-stationary case
Limits and pitfalls
References

Why this matters

Here are two time series. Their correlation is \(-0.02\): strictly nothing. You overlay them, you compute a Pearson coefficient, you conclude they have no relationship, and you move on.

You have just missed a dependence of \(-0.76\). It was there, in full, but shifted by three time steps. The first series does not react to the second at the same instant — it reacts three periods later. No correlation coefficient can see that, because a correlation coefficient examines exactly one alignment: right now.

This is what the vector framework makes visible. Where univariate analysis handles each series "one at a time", the VAR model lets each one listen to the past of the others. Three interest rates responding together to a central bank announcement, two sensors hit by the same thermal shock, one economic variable driving another with a delay: in every one of these cases, the useful information is precisely what the series-by-series approach throws away.

And there is a mirror-image trap, nastier still. Two smooth, entirely unrelated series look correlated at almost every lag. You will see why here — and the remedy, prewhitening, takes four lines. Two labs let you deform a VAR model by hand, watch its roots cross the unit circle, and collapse phantom correlations with a single click.

VAR(p) Cross-correlations Multivariate Yule–Walker Prewhitening Haugh (1976)

The vector vocabulary

A multivariate stochastic process \(\{\mathbf{X}_t;\, t \in T\}\) is simply a collection of vector-valued random variables. Each component \(\{X_{ti}\}\) is an ordinary univariate time series, and we stack them:

The observation vector at time t

\[\mathbf{X}_t = (X_{t1},\, X_{t2},\, \dots,\, X_{tm})\'\]

For the monthly US interest rate data from January 1991 to April 2010, \(X_{t1}\) denotes the six-month certificate of deposit rate, \(X_{t2}\) Moody's seasoned Aaa corporate bond rate, and \(X_{t3}\) the thirty-year conventional mortgage rate. The three univariate realizations form a single realization of a trivariate process.

The mean becomes a vector, the covariance a matrix:

Mean and 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]\]

The entry \(\gamma_{ij}(t_1,t_2)\) of that matrix is the covariance between \(X_{t_1 i}\) and \(X_{t_2 j}\): two different series, at two different times. That is exactly the information a univariate analysis discards.

Stationarity and cross-correlations

As in the univariate case, we require the moments not to depend on the time origin. The process is covariance stationary if \(E[\mathbf{X}_t] = \boldsymbol\mu\) is constant, if \(\Gamma(0)\) exists and is constant, and if \(\Gamma(t_1,t_2)\) depends only on the gap \(t_2 - t_1\). We then write:

Covariance matrix at lag 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]\]

When \(i \neq j\), \(\gamma_{ij}(k)\) is called the cross-covariance at lag \(k\). The diagonal entries of \(\Gamma(0)\) are the component variances. The cross-correlation normalises:

Cross-correlation

\[\rho_{ij}(k) = \frac{\gamma_{ij}(k)}{\sqrt{\gamma_{ii}(0)\,\gamma_{jj}(0)}}\]

Two properties are worth memorising, because they do not look like their univariate counterparts. In the univariate case \(\gamma(k) = \gamma(-k)\): the autocovariance is even. Here that is not true for the cross terms. All we have is:

Cross symmetry

\[\gamma_{ij}(k) = \gamma_{ji}(-k), \qquad \Gamma(-k) = \Gamma\'(k)\]

Note the swap in the index order. That is precisely what makes it possible to tell "\(X_1\) leads \(X_2\)" from "\(X_2\) leads \(X_1\)" — the asymmetry carries all the directional information. The values \(\rho_{ij}(0)\), meanwhile, measure the contemporaneous relationship, and \(\rho(0)\) has ones along its diagonal.

The warning that costs people hours. Two incompatible definitions of cross-covariance coexist in the literature. The one used here, from Box et al. (2008), sets \(\gamma_{ij}(k) = E[(X_{ti}-\mu_i)(X_{t+k,j}-\mu_j)]\). The one used by Shumway and Stoffer, by Tsay, and by Brockwell and Davis sets \(\gamma^{(a)}_{ij}(k) = E[(X_{t+k,i}-\mu_i)(X_{tj}-\mu_j)]\). The two are not equivalent: they flip the sign of the lag. In R, the base function ccf(xi, xj) computes the second; to get the first you must call ccf(xj, xi). Always check which convention a source uses before interpreting a peak at a given lag — otherwise you will conclude that the cause follows the effect.

Multivariate white noise

The series \(\{\mathbf{a}_t\}\) is multivariate white noise if \(E[\mathbf{a}_t] = \mathbf{0}\) and \(\boldsymbol\rho_a(k) = \mathbf{0}\) for every \(k \neq 0\). Nothing is imposed on \(\boldsymbol\rho_a(0)\). In other words, the components may be contemporaneously correlated: \(\rho_{ij}(0) \neq 0\) is perfectly allowed. We write \(\mathbf{a}_t \sim \text{MWN}(\mathbf{0}, \Gamma_a)\) with \(\Gamma(0) = \Gamma_a\).

This sounds contradictory at first — correlated white noise? But the two notions live on different axes. White noise forbids any structure along the time axis; it forbids nothing across, between components at the same instant. Two sensors experiencing the same thermal shock at every instant produce exactly this. The canonical example from the chapter is:

Bivariate white noise with contemporaneous correlation

\[\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}\]

Each realization taken on its own looks perfectly uncorrelated along the time axis; overlaid, they rise and fall together. This \(\Gamma_a\) is the one feeding the VAR(1) in the lab below.

The VAR(p) process

The central model of the chapter. Each component is regressed on the \(p\) past values of all components, including itself. In two dimensions and at order one:

A bivariate VAR(1), written row by row

\[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}\]

which condenses into \(\mathbf{X}_t = \Phi_1 \mathbf{X}_{t-1} + \mathbf{a}_t\). The general case:

Definition — VAR(p) process

\[\mathbf{X}_t = \boldsymbol\alpha + \Phi_1\mathbf{X}_{t-1} + \cdots + \Phi_p\mathbf{X}_{t-p} + \mathbf{a}_t\]

where \(\mathbf{X}_t\) is \(m \times 1\), each \(\Phi_k\) is a real \(m \times m\) matrix, \(\mathbf{a}_t\) is multivariate white noise with covariance \(\Gamma_a\), and \(\boldsymbol\alpha = (I - \Phi_1 - \cdots - \Phi_p)\boldsymbol\mu\).

In operator notation, \(\Phi(B)(\mathbf{X}_t - \boldsymbol\mu) = \mathbf{a}_t\) with \(\Phi(B) = I - \Phi_1 B - \cdots - \Phi_p B^p\). The first real result concerns stationarity, and it is a purely algebraic condition:

Theorem — stationarity and causality condition

The process is stationary and causal if and only if the roots of the determinantal equation

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

all lie outside the unit circle.

This is the exact analogue of the univariate condition on \(\phi(z) = 0\). For a VAR(1), \(|\Phi(z)| = |I - \Phi_1 z|\), and the roots are the reciprocals of the eigenvalues of \(\Phi_1\): the condition amounts to requiring that every eigenvalue of \(\Phi_1\) have modulus strictly less than 1. The lab below makes that condition visible.

Under stationarity the process admits a general linear representation \(\mathbf{X}_t - \boldsymbol\mu = \Psi(B)\mathbf{a}_t = \sum_{j \geq 0}\Psi_j \mathbf{a}_{t-j}\) with \(\Psi(B) = \Phi^{-1}(B)\), and the autocovariances satisfy the multivariate Yule–Walker equations:

Multivariate Yule–Walker

\[\Gamma(k) = \Gamma(k-1)\Phi_1\' + \cdots + \Gamma(k-p)\Phi_p\', \qquad k = 1,\dots,p\]

And for \(k = 0\), the relation \(\Gamma(0) - \Phi_1\Gamma(0)\Phi_1\' = \Gamma_a\) — a discrete Lyapunov equation — determines the stationary covariance. It is solved by vectorising: \(\operatorname{vec}\Gamma(0) = (I - \Phi_1 \otimes \Phi_1)^{-1}\operatorname{vec}\Gamma_a\).

Lab: the VAR(1) of Example 10.2

The example worked through in the chapter is a bivariate VAR(1) whose parameters can be fully reconstructed from the published covariances. Applying \(\Phi_1\' = \Gamma(0)^{-1}\Gamma(1)\) to the book's matrices recovers:

The parameters of the example

\[\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}\]

These values reproduce exactly the characteristic equation \(1 - 1.6z + 0.81z^2 = 0\) mentioned in the text, whose roots are complex with modulus \(1/\sqrt{0.81} = 1.11\): outside the unit circle, hence stationary, and complex, hence the pseudo-cyclic behaviour of period roughly 13 visible in the realizations. Through the Lyapunov equation they also reproduce the book's stationary covariance:

The resulting covariance and correlations

\[\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}\]

Here is the single most instructive fact in the whole chapter. The contemporaneous correlation between the two series is \(-0.02\): strictly nothing. An analyst who overlays the two curves and computes a Pearson coefficient concludes independence. And yet \(\rho_{12}(3) = -0.76\): shift one of them by three periods and the relationship appears, and it is massive. A strong relationship can be completely invisible at lag zero.

φ₁₁ 0.90 φ₁₂ 0.60 φ₂₁ -0.30 φ₂₂ 0.70

n 200

Things to try. Set \(\phi_{12} = 0\) and \(\phi_{21} = 0\): the two series become two perfectly separate AR(1) processes, and the theoretical cross-correlations collapse towards zero at every lag — except at lag 0, where the contemporaneous correlation of the noise remains. Then raise \(\phi_{12}\) again: the cross structure reappears, asymmetric. Finally push \(\phi_{11}\) towards 1.2: a root enters the unit circle and the process explodes.

A word on the cross-spectrum

The frequency domain generalises in the same way. The spectrum becomes an \(m \times m\) matrix at each frequency:

Multivariate spectrum

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

The diagonal entries \(P_{jj}(f)\) are the usual univariate spectra. The off-diagonal entries \(P_{jh}(f)\), called cross-spectra, are complex valued. They decompose as \(P_{jh}(f) = c_{jh}(f) - i\,q_{jh}(f)\), where \(c_{jh}\) is the cospectrum and \(q_{jh}\) the quadrature spectrum. In polar form \(P_{jh}(f) = \alpha_{jh}(f)e^{i\varphi_{jh}(f)}\), the modulus is the amplitude spectrum and \(\varphi_{jh}(f) = \tan^{-1}(-q_{jh}(f)/c_{jh}(f))\) the phase spectrum — which gives the phase shift between the two series at each frequency. Finally the squared coherency:

Squared coherency

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

It reads as an \(R^2\) frequency by frequency: the share of the variability of \(X_j\) at frequency \(f\) linearly explainable by \(X_h\) at that same frequency. Two series can be highly coherent at low frequencies and not at all at high ones — a diagnostic the cross-correlation alone does not give.

Estimation and order selection

The estimators are direct analogues of the univariate case. The mean is estimated by the vector of sample means \(\bar{\mathbf{X}}\), and the covariances by:

Sample cross-covariance

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

Note the \(1/n\) rather than \(1/(n-k)\): a biased estimator, but with smaller variance, and one that guarantees a non-negative definite matrix. For \(k = 0\), \(\hat\rho_{ij}(0)\) is simply the usual Pearson coefficient computed on the \(n\) observations.

For the VAR parameters there are two routes: the multivariate Yule–Walker equations, or equation-by-equation least squares. They give very similar results in practice. The order \(p\) is chosen by AIC — indispensable, since a VAR(\(p\)) in dimension \(m\) has \(m^2 p\) coefficients: a trivariate VAR(3) already has 27.

# ─── R: base + the vars package ───
# base ar() computes Yule-Walker estimates for a VAR;
# note that the mle and burg methods are univariate only
data <- cbind(x1, x2)
ar(data, order.max = 8, aic = TRUE, method = 'yw')
ar(data, order.max = 8, aic = TRUE, method = 'ols')

library(vars)
fit <- VAR(data, p = 1, type = c("const"), ic = c("AIC"))
serial.test(fit, lags.pt = 24, type = c("PT.asymptotic"))  # multivariate Ljung-Box
predict(fit, n.ahead = 25)   # forecasts + 95% limits by default

The code below really runs in your browser (Python compiled to WebAssembly). Edit the values and re-run: the first launch takes a few seconds while Python loads.

PYTHON · NUMPY

# Full reconstruction of Example 10.2 from the chapter.
# Starting from Phi_1 and Gamma_a, we recover EVERY published quantity.
import numpy as np

PHI = np.array([[ 0.9, 0.6],
                [-0.3, 0.7]])
GA  = np.array([[1.00, 1.78],
                [1.78, 4.00]])

# 1. characteristic equation |I - PHI z| = 0 for a bivariate VAR(1)
#    = 1 - tr(PHI) z + det(PHI) z^2
tr, det = np.trace(PHI), np.linalg.det(PHI)
print(f"characteristic equation: 1 - {tr:.2f} z + {det:.2f} z^2 = 0")
racines = np.roots([det, -tr, 1])
print(f"roots: {np.round(racines, 3)}   moduli: {np.round(abs(racines), 3)}")
print("stationary:", np.all(abs(racines) > 1), " (roots outside the unit circle)")
print(f"pseudo-cycle period: {2*np.pi/abs(np.angle(racines[0])):.1f} periods\n")

# 2. Gamma(0) from the discrete Lyapunov equation: G - PHI G PHI' = GA
G0 = gamma_zero(PHI, GA)
print("Gamma(0) =\n", np.round(G0, 2), "   (book: 27.74, -0.40 / -0.40, 13.07)\n")

# 3. multivariate Yule-Walker: Gamma(k) = Gamma(k-1) PHI'
G = [G0]
for k in range(1, 6):
    G.append(G[-1] @ PHI.T)
print("Gamma(1) =\n", np.round(G[1], 2), "   (book: 24.73, -8.60 / 7.48, 9.27)\n")

D = np.diag(1/np.sqrt(np.diag(G0)))
for k in range(4):
    print(f"rho({k}) =", np.round(D @ G[k] @ D, 2).tolist())

plot_var(*simule_var1(PHI, GA, n=75, seed=3), "A VAR(1) realization, n = 75")

PYTHON · NUMPY

Testing for association: the Haugh procedure

A frequent question: are these two series related? We would like to test \(H_0 : \rho_{12}(k) = 0\) by comparing \(\hat\rho_{12}(k)\) with bands of the form \(\pm 2/\sqrt{n}\). But the asymptotic variance of \(\hat\rho_{12}(k)\) is:

Asymptotic variance of the sample cross-correlation

\[\mathrm{AV}_{12} = \frac{1}{n}\sum_{j=-\infty}^{\infty}\rho_{11}(j)\,\rho_{22}(j)\]

Look closely at that formula. It involves the autocorrelations of both series. If both are strongly autocorrelated the sum is large, and the true variance far exceeds \(1/n\). The bands \(\pm 2/\sqrt{n}\) are then far too narrow: you see relationships everywhere. The variance reduces to \(1/n\) only if at least one of the two series is white noise.

Hence the procedure recommended by Haugh (1976), in four steps:

Step	Operation
1	Fit an ARMA(\(p_i, q_i\)) to each of the two series — the orders may differ
2	Compute the residual series \(\hat a_{t1}\) and \(\hat a_{t2}\) (this is the prewhitening)
3	Compute the cross-correlations \(\hat\rho^{(a)}_{12}(k)\) between the two residual series
4	Declare \(\hat\rho^{(a)}_{12}(k)\) significant at level \(\alpha = 0.05\) if \(\|\hat\rho^{(a)}_{12}(k)\| > 2/\sqrt{n}\)

The residuals are white noise by construction, so the asymptotic variance returns to \(1/n\) and the usual bands regain their validity. As always with repeated tests, the level applies separately to each \(k\): across 41 lags tested, two chance exceedances are expected under \(H_0\).

shift k +3 AR order 4 n 200

The decisive experiment. Switch the source to "independent AR(1)". The two series are generated by two unrelated generators, but each is strongly autocorrelated. Without prewhitening, the cross-correlogram shows wide peaks crossing the bands — an entirely fictitious relationship. Turn prewhitening on: the peaks collapse back inside the bands. This is exactly the mechanism that, in the chapter, leads to the conclusion that there is no association between solar activity and the earth's temperature: the raw cross-correlation reaches \(\hat\rho_{12}(6) = 0.271\), but after fitting an AR(9) to the sunspot numbers and an AR(4) to the temperatures, no residual correlation exceeds the limits \(\pm 2/\sqrt{160} = \pm 0.158\).

PYTHON · NUMPY

# The Haugh (1976) procedure on two STRICTLY independent series that are
# nonetheless strongly autocorrelated. How many false positives without prewhitening?
import numpy as np

def two_independent_ar(n, phi=0.95, seed=0):
    rng = np.random.default_rng(seed)
    x1 = np.zeros(n); x2 = np.zeros(n)
    e1 = rng.normal(0, 1, n); e2 = rng.normal(0, 1, n)   # no link whatsoever between e1 and e2
    for t in range(1, n):
        x1[t] = phi*x1[t-1] + e1[t]
        x2[t] = phi*x2[t-1] + e2[t]
    return x1, x2

n, K = 200, 20
x1, x2 = two_independent_ar(n, phi=0.95, seed=4)
limit = 2/np.sqrt(n)

k, r_raw = ccf(x1, x2, K)                  # without prewhitening
a1 = residus_ar(x1, p=4)                   # Haugh steps 1 and 2
a2 = residus_ar(x2, p=4)
_, r_white = ccf(a1, a2, K)                # steps 3 and 4

print(f"limits: +/- {limit:.3f}   ({2*K+1} lags tested)")
print(f"exceedances WITHOUT prewhitening: {int(np.sum(abs(r_raw) > limit)):3d}")
print(f"exceedances WITH prewhitening   : {int(np.sum(abs(r_white) > limit)):3d}"
      f"   (expected under H0: ~{0.05*(2*K+1):.1f})")
print(f"\nmax |rho| raw        : {abs(r_raw).max():.3f}")
print(f"max |rho| prewhitened: {abs(r_white).max():.3f}")

plot_ccf2(k, r_raw, r_white, limit,
          "Cross-correlations: raw (a) then after prewhitening (b)")

The non-stationary case

Everything above assumes stationarity. When it fails, the ARUMA model generalises to the multivariate case by inserting a diagonal operator \(\lambda(B)\):

VARUMA model

\[\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)\]

where each \(\lambda_j(z) = 0\) has all its roots on the unit circle, while those of \(|\Phi(z)| = 0\) and \(|\Theta(z)| = 0\) remain outside. The case that has received the most attention is \(\lambda(B) = \operatorname{diag}\big((1-B)^{d_1}, \dots, (1-B)^{d_m}\big)\) — the exact multivariate counterpart of the ARIMA(\(p,d,q\)) model. If component \(j\) is already stationary, \(d_j = 0\) and no differencing is applied to it.

What comes next? Differencing each component separately is sometimes the right answer — and sometimes a severe loss of information. When the components are cointegrated, a linear combination of them is already stationary, and differencing destroys it. That is the subject of the next article, which picks up this exact VAR framework and carries it to the error correction model.

Cointegration →

Limits and pitfalls

Assumption or pitfall	Consequence → remedy
Parameter count in \(m^2p\)	Rapid overfitting once \(m \geq 4\) → AIC/BIC, Bayesian VAR, or zero restrictions
Cross-covariance convention	Lag sign flipped depending on the source → check it, and prefer `ccf(xj, xi)` in R for the convention used here
Reading a raw CCF	Phantom correlations between autocorrelated series → prewhiten systematically
Zero contemporaneous correlation	Does not mean independence: see \(\rho_{12}(0) = -0.02\) alongside \(\rho_{12}(3) = -0.76\) in the same model
Non-stationary components	Spurious regression → test for unit roots, then for cointegration, before differencing
White noise with contemporaneous correlation	A non-diagonal \(\Gamma_a\) complicates causal interpretation → Cholesky decomposition and variable ordering, which is not neutral
VARMA beyond VAR	Serious identifiability problems → in practice one almost always stays with a pure VAR

References

Woodward, W. A., Gray, H. L. & Elliott, A. C. (2017). Applied Time Series Analysis with R, 2nd ed., CRC Press — chapter 10, the direct source of this article (definitions, Example 10.2, testing procedure).
Haugh, L. D. (1976). Checking the independence of two covariance-stationary time series: a univariate residual cross-correlation approach, JASA 71 — the prewhitening procedure.
Brockwell, P. J. & Davis, R. A. (1991, 2002). Time Series: Theory and Methods — the asymptotic theory of multivariate estimators.
Reinsel, G. C. (1997). Elements of Multivariate Time Series Analysis, 2nd ed., Springer.
Box, G. E. P., Jenkins, G. M. & Reinsel, G. C. (2008). Time Series Analysis: Forecasting and Control, 4th ed. — the cross-covariance convention adopted here.
Lütkepohl, H. (2005). New Introduction to Multiple Time Series Analysis, Springer — the comprehensive reference on VAR models.
Pfaff, B. (2008). VAR, SVAR and SVEC models: implementation within R package vars, Journal of Statistical Software 27(4).

The simulators on this page are written in plain JavaScript/Canvas with no external dependency; the Python cells run through Pyodide. The source is available on GitHub.