# Covariance of two time series driven by a restricted VAR(1) model

Suppose that I have two time series $X_n$ and $Y_n$ where:

$$X_n = \rho_x X_{n-1} + \epsilon_n$$ and $$Y_n = \rho_y Y_{n-1} + \rho_{xy}X_n +z_n$$

Here, $z_n,\epsilon_n$ are independent random variables, and the $\rho$'s can be thought of as the correlations. Now, I am trying to find a recursive relationship for $Cov(X_n, Y_n)$, and then take the limit. However, I am finding too many terms to deal with. Is there a simpler way than to go straight by the definition?

What I have so far is to start with $Cov(X_1,Y_1) = Cov(\rho_x X_0 + \epsilon_1, \rho_y Y_{0} + \rho_{xy}X_1 +z_1) = \rho_{xy}$ due to $X_0, Y_0$ being constants. Then, $Cov(X_2,Y_2) = \rho_x \rho_y Cov(X_1,Y_1)+ \rho_x^2 \rho_{xy}+ \rho_{xy}$.

Then, $Cov(X_n,Y_n) = \rho_x \rho_y Cov(X_{n-1},Y_{n-1})+ \rho_x^2 \rho_{xy}+ \rho_{xy}$.

However, this still gives me a mess I cannot rewrite in terms of any known summation formula. Did I do it correctly?

• This is almost a VAR(1) model except for the fact that $X_n$ enters the second equation contemporaneously rather than with a lag. But you could substitute $X_n$ in the second equation by its expression from the first equation to get a restricted VAR(1) model (or something very similar to it). As currently formulated, this is kind of a structural VAR(1) model. Perhaps this insight could be helpful (e.g. refer to literature on VAR and structural VAR models). Also, consider posting your intermediate result after which you feel it becomes too complicated. – Richard Hardy Feb 9 '16 at 8:36

One direction to go in might be something like: $$X_n = \rho_x X_{n-1} + \epsilon_n$$ and $$Y_n = \rho_y Y_{n-1} + \rho_{xy}X_n +v_n$$

Substitute for $$X_n$$ in the second equation: $$Y_n = \rho_y Y_{n-1} + \rho_{xy} \rho_x X_{n-1} + \rho_{xy} \epsilon_n +v_n$$

Let $$Z_n = \left[ \begin{array}{c} X_n \\ Y_n \end{array} \right]$$. As suggested by Richard, we can write the above equations as a single VAR(1) process $$Z_n = A Z_{n-1} + B U_n$$.

$$\underbrace{\left[ \begin{array}{c} X_n \\ Y_n \end{array} \right] }_{Z_n} = \underbrace{\left[ \begin{array}{cc} \rho_x & 0 \\ \rho_{xy} \rho_x&\rho_y \end{array} \right]}_{A} \underbrace{\left[ \begin{array}{cc} X_{n-1} \\ Y_{n-1} \end{array} \right]}_{Z_{n-1}} + \underbrace{\left[ \begin{array}{cc} 1 & 0 \\ \rho_{xy} & 1 \end{array} \right]}_B \underbrace{\left[ \begin{array}{c} \epsilon_n \\ v_n \end{array} \right]}_{U_n}$$

The process is is mean zero so we can write the covariance as:

$$E\left[ Z_n Z'_n \right] = E\left[ \left( AZ_{n-1} + B U_{n-1}\right) \left( AZ_{n-1} + B U_{n-1}\right) ' \right]$$ $$= A E\left[ Z_{n-1} Z_{n-1}'\right] A' + BB'$$ If the process is stationary, $$E\left[ Z_n Z'_n \right] = E\left[ Z_{n-1} Z'_{n-1} \right] = \Sigma$$. Under some technical conditions $$\Sigma$$ is the solution to: $$\Sigma = A \Sigma A' + BB'$$ This is a linear system of equations. Perhaps more conveniently for solving numerically with standard software, you can write $$\Sigma$$ as a vector. The resulting system is: $$vec(\Sigma) = (I - A \otimes A)^{-1} vec(BB')$$ Where $$\otimes$$ is the Kronecker product and vec is the vec operator. Basically, $$\Sigma$$ is the solution to a system of 4 equations.

After taking Richard's hint into account, you can appeal to the general result of an autocovariance function of a stable $VAR(1)$

$$y_t=c+\Phi y_{t-1}+\epsilon_t,$$

with $\Omega$ the variance-covariance matrix of $\epsilon_t$: $$\Gamma_j=\sum_{i=0}^\infty\Phi^{j+i}\Omega(\Phi^{i})^\top$$ This follows from writing the $VAR(1)$ as a vector moving average $$y_t=\mu+\sum_{j=0}^\infty\Phi^j\epsilon_{t-j},\quad\text{where}\quad \mu=(I-\Phi)^{-1}c,$$ and hence \begin{eqnarray*} \Gamma_j &=& E(y_t - \mu)(y_{t-j}-\mu)^\top\\ &=& E\left(\sum_{i=0}^{\infty}\Phi^i\epsilon_{t-j}\right)\left(\sum_{k=0}^{\infty}\Phi^k\epsilon_{t-k-j}\right)^\top\\ &=& E\left(\sum_{i=j}^{\infty}\Phi^i\epsilon_{t-i}\right)\left(\sum_{k=0}^{\infty}\Phi^k\epsilon_{t-k-j}\right)^\top\\ &=& E\left(\sum_{i=j}^{\infty}\Phi^i\epsilon_{t-i}\right)\left(\sum_{k=0}^{\infty}\epsilon_{t-k-j}{\Phi^k} ^\top\right)\\ &=& \sum_{i=j}^{\infty}\sum_{k=0}^{\infty}\Phi^iE\left(\epsilon_{t-i}\epsilon_{t-k-j}^\top\right){\Phi^k}^\top\\ &=& \Phi^j\Omega(\Phi^0)^\top + \Phi^{j+1}\Omega\Phi^\top + \Phi^{j+2}\Omega(\Phi^2)^\top + \ldots\\ &=& \sum_{i=0}^{\infty}\Phi^{j+i} \Omega {\Phi^i}^\top \end{eqnarray*}