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?
 A: 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.
A: 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*} 
