Let $X_1, \dots, X_m$ be a $n$-state ergodic stationary Markov chain with transition probability matrix $P$ and stationary distribution $\pi$. For $(i,j) \in [n]^2$, let $T_{ij}$ the random variable corresponding to the number of observed transitions from $i$ to $j$. I am looking for a closed form expression for the second moment $\mathbf{E}_{\pi}\left[ (T_{ij})^2 \right]$ when the chain is started from its stationary distribution.
My goal is to compute, or at least get a fair control over the variance of $T_{ij}$. I can assume reversibility of the chain if really needed.
Let $U_k$ be the indicator of the $k$-th transition is from $i$ to $j$: $\{X_{k-1} = i, X_k = j\}$. Then $$ T_{ij} = \sum_{k=1}^m U_k$$
$$ \begin{align} E[T_{ij}] &= \sum_{k=1}^m E[U_k] \\ &= \sum_{k=1}^m \Pr\{X_{k-1} = i, X_k = j\} \\ &= \sum_{k=1}^m \pi P^{k-1}e_ip_{i,j} \end{align}$$
where $\pi$ is a row vector containing the pmf of $X_0$, $e_i$ is a column vector with all the entries being $0$ except the $i$-th entry being $1$, $p_{i,j}$ is the $(i,j)$-th entry of the transition matrix $P$, the $1$-step transition probability from state $i$ to state $j$.
To compute the variance of $T_{ij}$, recall that $U_k$ is an indicator and thus $$ Var[U_k] = \pi P^{k-1}e_ip_{i,j}(1 - \pi P^{k-1}e_ip_{i,j})$$
And for $k < l$, $$ E[U_kU_l] = \Pr\{X_{k-1} = i, X_k = j, X_{l-1} = i, X_{l} = j\} $$
Obviously when $l = k+1$, the above probability is $0$. When $l > k + 1$,
$$ \begin{align} E[U_kU_l] &= \Pr\{X_{k-1} = i, X_k = j, X_{l-1} = i, X_{l} = j\} \\ &= \pi P^{k-1}e_ip_{i,j} p^{(l-1-k)}_{j,i}p_{i,j} \end{align}$$
where $p^{(l-1-k)}_{j,i} = e_i^T P^{l-1-k}e_j$ is the $(j, i)$-th entry of $P^{l-1-k}$, the $(l-1-k)$-step transition probability from state $j$ to state $i$.
So you can compute the covariance, and thus the variance accordingly.