Covariance of unconditional density

#1
I do have a minor problem with a simple covariance matrix.

Let x be a first order Markov chain of length n defined by the transition matrix P, and let \(y|x \sim N(\mu_{y|x}, \Sigma_{y|z} ) \). Find E(y) and Cov(y). Note that we assume \( \Sigma_{y|x}\) to be dense, and defined by a spatial correlation function \(\rho(h)\) where h = |i-j| for i,j = 1,.., n.

The mean and variance are easy to find,
\(\mu_i = E(y_i) = \sum_{x} \mu_{y|x}p_s(x)\)
\(Var(y_i) = \sum_{x}\sigma_{y|x}^2 p_s(x) + \sum_{x}(\mu_{y|x}-\mu_i)^2p_s(x)\)
where \(p_s(x)\) is the stationary distribution of x.

However, I get stuck when calculating the covariances using the total law of covariance. I get something like
\(\Sigma_y = Cov(y) = E(Cov(y|x)) + Cov(E(y|x)) = E(\Sigma_{y|x}) + Cov(\mu_{y|x})\).

From intuition I guess each element should be something like (j > i), but I can't prove it
\( [\Sigma_y ]_{ij} = \sum_{x_i} \sum_{x_j} (\sigma_{x_i} \sigma_{x_j} \rho(|i-j|) + (\mu_{x_i}-\mu_i)(\mu_{x_j}-\mu_i))p(x_j|x_i)\)