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)