Auto-Correlation of a discrete Markov chain goal: I am trying to derive a formula for the autocorrelation of a arbitrary discrete Markov chain, that has reached stationarity.
My question is very similar to the ones in question1 and quetsion2. I was not able to extend the answers given there to my more general problem, since I am not that much into probability theory and the other questions make explicit use of some numbers, which I cannot.
Let's assume the state space is $\{1,...,N\}$, we know the transition matrix $T$ and let's call the stationary distribution $\pi$. Then
$$
P(X_n=i|X_{n-1}=j) = T_{ij}\\
P(X_n=i) = \pi(i).
$$
Following the above mentioned posts, we can rewrite the problematic term in the correlation, $\mathbb{E}[X_nX_{n+d}]$, using stationarity and the law of total probability.
$$
\mathbb{E}[X_nX_{n+d}] = \mathbb{E}[X_0X_d] = \sum_{s=1}^N \mathbb{E}[X_0X_d|X_0=s]P(X_0=s) = \sum_{s=1}^N \pi(s)~\mathbb{E}[X_0X_d|X_0=s]
$$
I am not sure how to go on from here with the term $\mathbb{E}[X_0X_d|X_0=s]$. I am lacking a good understanding of the second last equal sign in the answer1 to question1.
Edit: Conif10's answer is pretty and clear. No need to read my confusing attempts.
What I tried so far was not very tricky, but a crude insertion of definitions: First I use the definitions of conditional expectation $\mathbb{E}[X_0X_d|X_0=s]$. I think of the product $X_0X_d =:Y$ as a new random variate $Y$ on the state space $\{1,...,N^2\}$.
$$
\mathbb{E}[X_0X_d|X_0=s] = \sum_{y=1}^{N^2} y P(Y=y|X_0=s)
$$
Then I again use the law of total probability recursivly
$$
P(Y=y|X_0=s) = \sum_{k_1} P(X_dX_0=y|X_0=s|X_1=k_1)P(X_1=k_1|X_0=s) = \\
\sum_{k_1, k_2} P(X_d X_0=y|X_0=s|X_1=k_1|X_2=k_1)P(X_1=k_1|X_0=s)P(X_2=k_2|X_1=k_1) = \\=~~...[\text{recursion}]... ~~=\\
\sum_{k_1,..., k_d} P(X_d X_0=y|X_0=s|X_1=k_1|...|X_d=k_d)P(X_d=k_d|X_{d-a}=k_{d-1})...P(X_1=k_1|X_0=s) = \\
= \sum_{k_1,..., k_d} \delta_{k_d s,y}T_{k_d,k_{d-1}}...T_{k_1s} = \sum_{k_d}\delta_{k_d s,y}T^d_{k_d,s} = \sum_{s'=1}^N\delta_{s' s,y}T^d_{s',s}
$$
Here I evaluated $P(X_d X_0=y|X_0=s|X_1=k_1|...|X_d=k_d)=\delta_{k_d s,y}$, since the value of $X_dX_0$ is fixed to $k_ds$ by the conditions and the deterministic product $k_ds$  can either be equal to $y$ (then $P=1$) or not ($P=0$). I also identified the transition matrix.
This is the part I am very insecure about, especially about the format definition of multiple conditionals like $P(A|B|C)$, which I found nothing about. Additionally I question my use of the law of total probability of the kind $P(A|C) = \sum_B P(A|B|C) P(B|C)$. However, if one puts everything together, the result is
$$
\mathbb{E}[X_0X_d] = \\
=\sum_{s=1}^N \pi(s)~\mathbb{E}[X_0X_d|X_0=s] = \\
=\sum_{s=1}^N \pi(s)~\sum_{y=1}^{N^2} y P(Y=y|X_0=s) =\\
=\sum_{s=1}^N \pi(s)~\sum_{y=1}^{N^2} y \sum_{s'=1}^N\delta_{s' s,y}T^d_{s',s} = \\
=\sum_{s,s'=1}^N \pi(s)~s'sT^d_{s's}
$$
where the (somehow strange) summation over $y$ collapses due to the Kronecker delta.
When I apply this general result to the answer1, I get the exact same. But this is only a indication not a proof for correctness of the above "derivation".
I hope there is a better/shorter/nicer approach.
Can someone either validate my derivation with some notes on the parts where I misuse math or give a better derivation?
 A: Your final result is correct, but I think that the definition of the random variable $Y$ is not necessary and it can lead to confusion.
Using your notation, we start from the result you got using the total law of probability under the stationarity assumption:
$$\mathbb{E}[X_{n+d}X_{n}] = \sum_{s = 1}^N\pi(s)\mathbb{E}[X_{d}X_{0}|X_{0} = s]
$$
Notice that, since the value $X_0 = s$ in the summation is fixed, it is constant with respect to the expectation value. Therefore, we have:
$$\mathbb{E}[X_{n+d}X_{n}] = \sum_{s = 1}^N\pi(s)s\mathbb{E}[X_{d}|X_{0} = s] 
$$
Using the conditional expectation we obtain:
$$\mathbb{E}[X_{n+d}X_{n}] =\sum_{s = 1}^N\pi(s)s\sum_{s' = 1}^Ns'P(X_{d} = s'|X_{0} = s) = \sum_{s,s' = 1}^N\pi(s)s'sT^d_{s's}
$$
which coincides with your final expression.
From this expression, the second last equal sign in answer1 reads:
$$... = \frac{1}{2}(\mathbb{E}[X_{n+d}|X_{n} = 1]-\mathbb{E}[X_{n+d}|X_{n} = -1]) =\frac{1}{2}\left(\sum_{s = \{-1,1\}}sP(X_{n+d} = s|X_{d} = 1) - \sum_{s = \{-1,1\}}sP(X_{n+d} = s|X_{d} = -1)\right) = \frac{1}{2}\left\{1\cdot(1+(2\alpha -1)^d)+(-1)\cdot(1-(2\alpha -1)^d) -[1\cdot(1-(2\alpha -1)^d)+(-1)\cdot(1+(2\alpha -1)^d)]\right\} = (2\alpha -1)^d
$$
In which we are considering all the possible transitions $s\to s'$.
Furthermore, notice that this problem could be also addressed in terms of the joint pdf $\rho(X_{n+d},X_{n})$:
$$\mathbb{E}[X_{n+d}X_{n}] = \sum_{s,s' = 1}^N s's\rho(X_{d} = s', X_{0} = s) = \sum_{s,s' = 1}^N s'sP(X_{d} = s'| X_{0} = s)\pi(X_0 = s) = \sum_{s,s' = 1}^N\pi(s)s'sT^d_{s's}
$$
