Suppose I have a transition density (matrix) $K$, and I want to calculate the covariance between random variables at two different times.
\begin{align*} Cov_{\mu_0}(X_n,X_{n+k})&=E_{\mu_0}(X_n-E_{\mu_n}X_n)(X_{n+k}-E_{\mu_{n+k}}X_{n+k})\\&=E_{\mu_n}(X_0-E_{\mu_0}X_0)(X_k-E_{\mu_k}X_k)\\&=\int\cdots\int(x_k-E_{\mu_k}X_k)(x_0-E_{\mu_0}X_0)K(x_{k-1},dx_k)K(x_{k-2},dx_{k-1})\cdots K(x_0,dx_1)\mu_n(dx_0)\\&=??? \end{align*}
How can I efficiently calculate this (both in theory and in computer simulation)? I know that
\begin{align*} E_{\mu_0}f(X_k)&=\int\cdots\int f(x_k)K(x_{k-1},dx_k)\cdots K(x_0,dx_1)\mu_0(dx_0)\\&=???(\text{don't know the intermediate steps})\\&=\sum_i\mu^0_iK_{ij}^kf(x_i)\\&=\mu_0'K^{k}f \end{align*}
Hint: use law of iterated expectation to simplify calculation