These notes are partially based on those of Mathias Drton.

(\(\spadesuit\spadesuit\spadesuit\)) = Contents for graduate students.

Introduction

We have learned a lot about properties of a Markov chain – when we know the transition probability matrix, we can do a lot of analysis. In data analysis, we often do not know the transition probability matrix so we need to estimate it. As long as we have a transition probability matrix, we can use our knowledge about Markov chain to make inference. In this lecture, we will be thinking about how to recover the transition probability matrix and make related statistical inference.

Let \(X_0,\cdots,X_n\) be a set of RVs denoting our data/observations. Now we will assume that these RVs form a Markov chain. Moreover, we assume that each observation \(X_i \in S=\{1,2,\cdots, s\}\). Namely, the state space is finite.

In both Frequentist and Bayesian approaches, the likelihood function plays a key role in inference. Let \(\mathbf{P}=\{p_{ij}\}\) be the transition probability matrix and \(\nu_i = P(X_0=i)\) be the initial probability. The likelihood function can be written as \[L_n(\nu, \mathbf{P}) = \nu_{X_0} \prod_{i=0}^{n-1} p_{X_{i}, X_{i+1}}.\] Generally, we cannot estimate the initial distribution \(\nu\) unless we observed several different Markov chains. So we will make our inference condition on \(X_0\).

Recalled that the joint PMF can be written as \[p(x_0,\cdots, x_1; \mathbf{P}, \nu) = p(x_1,\cdots, x_n; \mathbf{P}|x_0) \times p(x_0; \nu),\] we then use the likelihood function \[L_n(\mathbf{P}) = p(X_1,\cdots, X_n; \mathbf{P}|x_0) = \prod_{i=1}^s \prod_{j=1}^s p_{ij}^{n_{ij}},\] where \(n_{ij} = \sum_{k=0}^{n-1} I(X_k=i, X_{k+1}=j)\) is the number of transition from state \(i\) to state \(j\).

Frequentist Approach

With the likelihood function, we first consider the Frequentist approach to estimate \(\mathbf{P}\). We find the MLE of \(\mathbf{P}\).

The log-likelihood function will be \[\ell_n(\mathbf{P}) = \log L_n(\mathbf{P}) = \sum_{i=1}^s \sum_{j=1}^s n_{ij} \log p_{ij}.\] To find the MLE, note that there are constraints on \(\mathbf{P}\): \[\sum_{j=1}^s p_{ij} = 1\] for every \(i=1,\cdots,s\). These constraints only applies to each row of \(\mathbf{P}\), so we can separate the maximization of \(\ell_b(\mathbf{P})\) into \(s\) different maximizations: \[\begin{align*} \widehat{p}_{i} = {\sf argmax}_{\sum_{j=1}^s p_{ij}=1} \ell_n(p_{i1},\cdots,p_{is}) = {\sf argmax}_{\sum_{j=1}^s p_{ij}=1} \sum_{j=1}^s n_{ij} \log p_{ij}. \end{align*}\] It turns out that this is the MLE of the multinomial distribution so \[\widehat{p}_{ij} = \frac{n_{ij}}{n_{i+}},\] where \(n_{i+} = \sum_{j=1}^s n_{ij}\) is the number of observations from starting at state \(i\).

So now we have the MLE but does this MLE has the usual properties such as statistical consistency and asymptotic normality? We cannot directly apply our conventional statistical knowledge here because each transition from state \(i\) to state \(j\) may not be independent from each other. We need some theories of MLE under the Markov chain structure.

With Proposition 1, we are able to construct a confidence interval of each \(p_{ij}\). Roughly speaking, Proposition 1 shows that the estimating transition probability matrix of a Markov chain behaves like the conventional IID multinomial random variable problem that the sample proportion is a consistent and asymptotic normal estimator with the asymptotic variance being adjusted according to the stationary distribution.

Note that the asymptotic normality implies that \[{\sf Var}(\widehat p_{ij}) \approx \frac{p_{ij}(1-p_{ij})}{{n\pi_i}}.\] This comes from the fact that since the Markov chain will reach stationary distribution, there will be around \(n\pi_i\) numbers of transitions starting from state \(i\).

From the above analysis, we have a simple conclusion:

The analysis of Markov chain is just like the analysis of a contingency table that the empirical transition count matrix is the contingency table of current state (row) versus next state (column).

While a conventional contingency table assumes observations are independent, the weak dependency in the Markov chain does not invalidate the theoretical behavior.

Hypothesis test of the transition probability matrix

Suppose we have observed a long Markov chain \(X_0,X_1,\cdots, X_n\) and we have a hypothesis that the transition probability is a given matrix \(\mathbf{P}\) and the initial distribution \(\nu\).

Question: How do we infer the transition probability matrix \(\mathbf{P}\) and the initial distribution \(\nu\)?

For the initial distribution \(\nu\), unfortunately we cannot test it because we only observed a single Markov chain. No matter how long you have observed the Markov chain, the initial state only occurs once, so there is no way to test its distribution.

On the other hand, the transition probability matrix can be tested from the data. A simple approach to this is exactly how we handle contingency table–the \(\chi^2\)-test–which uses the test statistic \[\begin{equation} T_n = \sum_{i,j: \mathbf{P}_{ij}>0}\frac{ (O_{ij} - E_{ij})^2}{E_{ij}} = \sum_{i,j: \mathbf{P}_{ij}>0}\frac{(n_{ij} - n_{i+} \mathbf{P}_{ij})^2}{n_i \mathbf{P}_{ij}}, \label{eq::chi2} \end{equation}\] where \(O_{ij} = n_{ij}\) is the observed count on the entry \((i,j)\) and \(E_{ij} = n_{i+} \mathbf{P}_{ij}\) is the expected count under the assumed transition probability \(\mathbf{P}\). Under the null hypothesis \(T_n\) will follow from a \(\chi^2\) distribution with a degree of freedom \(d = \sum_{i,j} I(\mathbf{P}_{ij}> 0) -s\), which is the number of non-zero elements in \(\mathbf{P}\) minus \(s\).

Why the degree of freedom is \(d = \sum_{i,j} I(\mathbf{P}_{ij}> 0) -s\)? Each row of \(\mathbf{P}\) is a probability vector so it has to sum into \(1\). Therefore, if it has \(k\) non-zero entries, it has only effectively \(k-1\) free parameters. The summation over all rows leads to \(s\) minus one to the total number of free parameters.

Short but multiple chains

Another common scenario is that instead of having a single long Markov chain, we have many short Markov chains. Suppose we have \(N\) Markov chains \[X^{(1)},\cdots, X^{(N)}\] such that each \[X^{(i)} = (X^{(i)}_0,X^{(i)}_2,\cdots, X^{(i)}_{n_i})\] is a Markov chain from the same transition probability matrix \(\mathbf{P}\) and the initial state \(X^{(1)}_0,\cdots, X^{(N)}_0\sim \nu\), i.e., the initial states are IID from \(\nu\). Here we assume that the number of chains \(N\rightarrow\infty\) but the length of each chain \(n_i\) may not diverge.

Question: How do we infer the transition probability matrix \(\mathbf{P}\) and the initial distribution \(\nu\)?

Unlike the case of single chain, we are able to infer the initial distribution. For the initial distribution, it follows from classical multinomial MLE that we can estimate \[\widehat\nu_j = \frac{1}{N}\sum_{i=1}^n I(X^{(i)}_1 = j).\] And we can construct confidence interval or perform hypothesis test accordingly.

For the transition probability matrix, it turns out that we can do it by simply stack informations of these \(N\) chains together. Specifically, if we want to estimate \(p_{ij}\), we simply use the sample ratio of transition from state \(i\rightarrow j\) among every transition starting at state \(i\). This is simply because the total likelihood function is the product \[\begin{align*} L_N(\nu, \mathbf{P}) &= \prod_{\ell=1}^N L(\nu, \mathbf{P}|X^{(\ell)}) \\ & = \prod_{\ell=1}^N\nu_{X^{(\ell)}_0} \prod_{k=0}^{n_\ell-1} p_{X^{(\ell)}_{k}, X^{(\ell)}_{k+1}}\\ & = \left(\prod_{\ell=1}^N\nu_{X^{(\ell)}_0} \right) \left(\prod_{\ell=1}^N \prod_{k=0}^{n_\ell-1} p_{X^{(\ell)}_{k}, X^{(\ell)}_{k+1}}\right). \end{align*}\] So the likelihood function for \(\mathbf{P}\) only depends on the second part: \[L_N(\mathbf{P}) = \prod_{\ell=1}^N \prod_{k=0}^{n_\ell-1} p_{X^{(\ell)}_{k}, X^{(\ell)}_{k+1}}.\] Thus, \[\widehat p_{ij} = \frac{N_{ij}}{N_{i+}},\] where \[N_{ij} = \sum_{\ell=1}^N \sum_{k=0}^{n_\ell-1} I(X^{(\ell)}_{k} = i, X^{(\ell)}_{k+1} = j),\qquad N_{i+} = \sum_{\ell=1}^N \sum_{k=0}^{n_\ell-1} I(X^{(\ell)}_{k} = i).\] Namely, \(N_{ij}\) is the total number of observed transition from state \(i\) to state \(j\), regardless of which chain it is from while \(N_i\) is the total number of transition starting from state \(i\). Note that longer chains will naturally provide more information than shorter chains.

The result in Proposition 1 still applies but with one small correction, the asymptotic variance \(p_{ij} (1-p_{ij})n \pi_i\) is replaced by \[{\sf Var}(\widehat p_{ij}) \approx \frac{p_{ij} (1-p_{ij})}{N\cdot \left[\sum_{t=0}^{n-1} \nu \mathbf{P}\right]_i}\] assuming that every short chain has a length \(n\). The quantity \(\left[\sum_{t=0}^{n-1} \nu \mathbf{P}\right]_i\) is the expected number of transition starting from state \(i\). Thus, the quantity \(N\cdot \left[\sum_{t=0}^{n-1} \nu \mathbf{P}\right]_i\) is the rough size of number of transition from state \(i\) when we combine all \(N\) chains together. Note that the same \(\chi^2\) test approach is applicable with replacing \(n_{ij}\) by \(N_{ij}\) and \(n_{i+}\) by \(N_{i+}.\)

Bayesian Approach

When we adopt a Bayesian approach in the Markov chain, we need to put a prior distribution for the transition probability matrix \(\mathbf{P}\). Essentially, we need to put a prior on each \(p_{ij}\) with the constraint that \(\sum_{j=1}^s p_{ij}=1\).

Note that the constraints are on each row of \(\mathbf{P}\) so one possibility is to put a prior on each row of \(\mathbf{P}\). Here we will use the Dirichlet prior.

The Dirichlet distribution is a multivariate distribution over the simplex \(\sum_{i=1}^K x_i = 1\) and \(x_i\geq 0\). Its probability density function is \[p(x_1,\cdots, x_K; \alpha_1,\cdots, \alpha_K) = \frac{1}{B(\alpha)} \prod_{i=1}^K x_i^{\alpha_i-1},\] where \(B(\alpha) = \frac{\prod_{i=1}^K \Gamma(\alpha)}{\Gamma(\sum_{i=1}^K \alpha_i)}\) and \(\alpha = (\alpha_1,\cdots, \alpha_K)\) are the parameters of this distribution. The Dirichlet distribution generates a random vector with length \(K\) and each element of this vector is non-negative and summation of elements is \(1\), meaning that it generates a random probability vector. You can view it as a generalization of the Beta distribution. For \(Z = (Z_1,\cdots, Z_K)\sim{\sf Dirch}(\alpha_1,\cdots, \alpha_K)\), \(\mathbb{E}(Z_i) = \frac{\alpha_i}{\sum_{j=1}^K \alpha_j}\) and the mode of \(Z_i\) is \(\frac{\alpha_i-1}{\sum_{j=1}^K\alpha_j - K}\) so each parameter \(\alpha_i\) determines the relative importance of category (state) \(i\). Because it is a distribution putting probability over \(K\) categories, Dirchlet distribution is very popular in social sciences and linguistics analysis.

Going back to our data analysis problem, by choosing \(K=s\), the Dirchlet distribution can be used as a prior distribution for each row of \(\mathbf{P}\). For each \(p_{i} = (p_{i1},\cdots, p_{is})\), we use the Dirchlet prior with parameter \(\alpha_i = (\alpha_{i1},\cdots, \alpha_{is})\) such that \[p_{i} \sim {\sf Dirch}(p_i;\alpha_i).\] The posterior distribution can then be written as \[\begin{align*} \pi(\mathbf{P}|X_1,\cdots,X_n) &\propto\prod_{i=1}^s\prod_{j=1}^s p_{ij}^{n_{ij}} \prod_{k=1}^s p(p_k;\alpha_k)\\ & \propto\prod_{i=1}^s\prod_{j=1}^s p_{ij}^{n_{ij}} \left(\prod_{k=1}^s p_{k1}^{\alpha_{k1}-1}\times\cdots\times p_{ks}^{\alpha_{ks}-1}\right)\\ & = \prod_{i=1}^s \left(p_{i1}^{n_{i1}+\alpha_{i1}-1}\times\cdots\times p_{is}^{n_{is}+\alpha_{is}-1}\right), \end{align*}\] which is a product of each row.

So the posterior distribution of \(p_i\) given the data is \[\pi(p_i|X_1,\cdots,X_n) = p_{i1}^{n_{i1}+\alpha_{i1}-1}\times\cdots\times p_{is}^{n_{is}+\alpha_{is}-1},\] which is the PDF of \({\sf Dirch}(n_{i1}+\alpha_{i1},\cdots, n_{is}+\alpha_{is})\). This implies that the posterior mean and MAP are \[\widehat{p}_{ij, \pi} = \frac{n_{ij}+\alpha_{ij}}{\sum_{k=1}^s n_{ik}+\alpha_{ik}},\quad \widehat{p}_{ij, MAP} = \frac{n_{ij}+\alpha_{ij}-1}{\sum_{k=1}^s n_{ik}+\alpha_{ik}-K}.\]

The credible interval can be constructed using the level sets of the posterior distribution. In hypothesis test, we can use the Bayes factor and compare it with our prior on the hypotheses to decide if we can reject the null hypothesis or not.

When \(\mathbf{P}= \mathbf{P}(\theta)\)

In some cases, the transition probability matrix \(\mathbf{P}= \mathbf{P}(\theta)\) is determined by the parameter \(\theta \in\mathbb{R}^L\) (namely, each \(p_{ij} = p_{ij}(\theta)\)). And our goal is to estimate \(\theta\) from our data. We will talk about how to make inference in this case using the Frequentist approach (you can use Bayesian as well). Here we first look at a specific example about genetic drift.

This scenario–the parameters of a distribution is controlled by another set of parameters– is very common in real data analysis. When the parametric family is from an exponential family, this model is also known as the curved exponential family. Here we briefly mention a few results about the case of \(\mathbf{P}= \mathbf{P}(\theta)\). Note that in this case, the log-likelihood function is \[\ell_n(\theta) = \sum_{i,j=1}^s n_{ij} \log p_{ij}(\theta)\] and the score equations (score function \(=0\)) are \[\frac{\partial \ell_n(\theta)}{\partial \theta_k} = \sum_{i,j=1}^s \frac{n_{ij}}{p_{ij}(\theta)} \frac{\partial p_{ij}(\theta)}{\partial \theta_k}= 0\] for \(k=1,\cdots, L\) and the Fisher’s information matrix \(I_1(\theta) = \{I_{km}(\theta)\}\) is \[I_{km}(\theta) = \sum_{i,j=1}^s \frac{\pi_i(\theta)}{p_{ij}(\theta)} \frac{\partial p_{ij}(\theta)}{\partial \theta_k} \frac{\partial p_{ij}(\theta)}{\partial \theta_m},\] where \(\pi_i(\theta)\) is the stationary distribution of state \(i\).

Sections I.2 and I.5 in Billingsley (1961, Statistical Inference for Markov Processes)

With Proposition 5, we know that the MLE will be a good estimator of \(\theta_0\) and we can use the MLE with the variance estimator to construct an asymptotically valid confidence interval.

Note that if we are using a Bayesian approach for the problem of \(\mathbf{P}= \mathbf{P}(\theta)\), we need to find a prior distribution of \(\theta\). In the case of Wright-Fisher model with mutation, a possible choice of the prior will be the Beta distribution (for both \(u\) and \(v\)) since it is a conjugate prior fro the binomial distribution.

Testing the Independence versus Markov Dependence

When we observe a sequence of RVs \(X_0,\cdots, X_n\). Sometimes we are not sure if they are IID or they form a Markov chain. We can do a hypothesis test to examine which assumption is more plausible. Determining if the sequence is IID is very important because if they are IID, we have a lot more statistical tools to analyze the underlying population.

In this case, the null and the alternative hypotheses are \[\begin{align*} H_0:X_1,\cdots, X_n\mid X_0 &\mbox{ are IID.}\\ H_1:X_1,\cdots, X_n\mid X_0 &\mbox{ form a Markov chain with a t.p.m. of non-identical rows.} \end{align*}\] Under \(H_0\), the observations are like IID from a PMF \(p_i = P(X_1= i)\) whereas under \(H_1\), the underlying probability model is described by the transition probability matrix \(\mathbf{P}= \{p_{ij}\}\).

An interesting fact is – even under \(H_0\), we can still write the model using the transition probability matrix but with the requirement that every row is the same. Namely, we can use the quantity \(\mathbf{P}= \{p_{ij}\}\), where \(p_{ij} = P(X_n=j\mid X_{n-1}=i)\) to describe the case of both \(H_0\) and \(H_1\). In this scenario, \[\begin{equation} H_0: p_{ij} = p_{j} \mbox{ for all $i,j\in S$}. \label{eq::H0} \end{equation}\] Note that identifying a common quantity that both \(H_0\) and \(H_1\) is very important in hypothesis testing because we can then transform the problem of testing a ‘statement/property’ into testing some parameters of interest.

There are multiple ways to test equation \(\eqref{eq::H0}\). We use likelihood ratio test (LRT) and the Bayes factor in this case.

LRT. Recall that the likelihood ratio test would use the MLE under \(H_0\) and \(H_0\cup H_1\) to perform the test. Under \(H_0\), \(\widehat{p}_{ij, H_0} =\widehat{p}_j= \frac{n_{+j}}{n}\) while under \(H_1\), \(\widehat{p}_{ij} = \frac{n_{ij}}{n_{i+}}\). Thus, the test statistic is \[\begin{align*} T_n &= 2 (\ell_{n}(\widehat{\mathbf{P}}_{MLE})-\ell_{n}(\widehat{\mathbf{P}}_{MLE, H_0})) \\ &= 2 \sum_{i,j=1}^s n_{ij} \log \left(\frac{\widehat{p}_{ij}}{\widehat{p}_j}\right)\\ & = 2 \sum_{i,j=1}^s n_{ij} \log \left(\frac{n_{ij}\cdot n}{n_{i+}\cdot n_{+j}}\right). \end{align*}\]

To compute the p-value using LRT, we use the \(\chi^2\) distribution as a reference distribution of \(T_n\). Now here comes the question: what is the degree of freedom in this case? First, we think about \(H_1\). There are totally \(s^2\) parameters in \(\mathbf{P}\). Each row has to sum to \(1\) so there are \(s\) constraints. Thus, under \(H_1\), there are \(s^2-s\) degrees of freedom. What about \(H_0\)? Under \(H_0\), all rows are identical so there are at most \(s\) parameters. But all these \(s\) numbers have to sum to \(1\), leading to one constrain so there are totally \(s-1\) degrees of freedom. Thus, the remaining degrees of freedom is \((s^2-s) - (s-1) = (s-1)^2\). So \(T_n \sim \chi^2_{(s-1)^2}\).

\(\chi^2\)-test from contingency table. Alternatively, you may use the \(\chi^2\) independence test from contingency table, which uses the test statistic: \[W_n =\sum_{i,j=1}^s \frac{(n_{ij} - n_{i+}\cdot n_{+j}/n)^2}{n_{i+}\cdot n_{+j}/n}\] because under \(H_0\), the expected number of observed transition from \(i\) to \(j\) is \(E_{ij} = n_{i+} \cdot \widehat p_{ij} = n_{i+} \cdot n_{+j}/n\) while the observed count \(O_{ij} = n_{ij}\). You can show that \(W_n\) is asymptotically the same as \(T_n\) and follows from \(\chi^2_{(s-1)^2}\).

Bayes Factor. To use the Bayes factor, we need to put priors on the parameters \(p_i = (p_{i1},\cdots, p_{is})\). As we have mentioned, the Dirichlet piror seems to be a good choice. Under \(H_1\), although we can use different priors on different rows, here we use the same prior for every row. Moreover, we would use the same prior for both \(H_0\) and \(H_1\). This choice reflects the fact that we do not have any different believes of \(\mathbf{P}\) under the two hypotheses. So our prior is \[p_i\sim {\sf Dirichlet}(\alpha_1,\cdots, \alpha_s)\] for some given hyperparameters \(\alpha_1,\cdots, \alpha_s\). The Bayes factor is \[\begin{align*} {\sf BF} &= \frac{\int p({\sf Data}|\mathbf{P}) \pi(\mathbf{P}|H_0)d\mathbf{P}}{\int p({\sf Data}|\mathbf{P}) \pi(\mathbf{P}|H_1)d\mathbf{P}}\\ & = \frac{\frac{1}{B(\alpha)}\int \prod_{j=1}^s p_{j}^{n_{+j}} p_j^{\alpha_j-1}dp_j}{\frac{1}{B^s(\alpha)}\int \prod_{i,j=1}^s p_{ij}^{n_{ij}}p_{ij}^{\alpha_j-1}dp_{ij}}\\ & = \frac{\frac{1}{B(\alpha)} \prod_{j=1}^s\int p_{j}^{n_{+j}+\alpha_j-1} dp_j}{\frac{1}{B^s(\alpha)}\prod_{i,j=1}^s\int p_{ij}^{n_{ij}-{\alpha_j}-1}dp_{ij}}\\\ & = \frac{B(n_{+1}+\alpha_1,\cdots,n_{+s}+\alpha_s)/B(\alpha_1,\cdots, \alpha_s)}{\prod_{i=1}^s\left[B(n_{i1}+\alpha_1,\cdots,n_{is}+\alpha_s)/B(\alpha_1,\cdots, \alpha_s)\right]}. \end{align*}\]

Asymptotic normality of the MLE (\(\spadesuit\spadesuit\spadesuit\))

To derive the theory of MLE for Markov chain, we first introduce the concept of snake chain and state two useful lemmas. Let \(\{X_n\}\) be a homogeneous Markov chain with a state space \(S\) and a transition probability matrix \(\mathbf{P}\). Define \(Y_n = (X_n, X_{n+1})\). Then \(\{Y_n\}\) is called the snake chain and it is also a homogeneous Markov chain with a state space \(S_Y = \{(i_0,i_1)\in S^2: p_{i_0,i_1}>0\}\).

Going back to our data analysis problem, the snake chain provides a useful framework for analyzing the quantity \[n_{ij} = \sum_{k=0}^{n-1}I(X_{k}= i, X_{k+1}=j) = \sum_{k=0}^{n-1}I(Y_{k} = (i,j))\] can be written as summation of indicator function with snake chains.

Another problem we need to address in our MLE is that the ‘sample size’ for estimating \(p_{ij}\) is \(n_{i+}\) (see the denominator), which is also random. So the conventional central limit theorem cannot be applied. Here we introduce another useful lemma (sometimes it is called Anscombe Lemma) for handling this case.

See Anscombe (1952, Proceedings of the Cambridge Philosophical Society) or Theorem 11.6.1 in Shorack (2000, Probability for Statisticians).

Note that there are modified version of this lemma and is often related to the Poissonization techniques. It is particularly useful when we are trying to derive the central limit theory of an integrated error of a function estimation¹.

Proof of Proposition 1

With these two lemmas, we can show that the MLE also works in Markov chain.

Before proceeding to the proof, we first compare this asymptotic normality to the one for multinomial distribution. In multinomial case, we have \[\frac{\widehat{p}_{i}- p_{i}}{\sqrt{p_{i}(1-p_{i})/{n}}}\overset{D}{\rightarrow}N(0,1).\] Here you see that the difference is that the sample size changes from \(n\) to \(n\pi_i\). So we can say that \(n\pi_i\) behaves like the effective sample size. This quantity makes sense because we know that the stationary distribution characterizes the long run proportion of each state. The transition probability \(p_{ij}\) is estimated by the proportion of transiting from state \(i\) to state \(j\) over the the total number of transiting from state \(i\). So the long run proportion of state \(i\) determines how much transition from state \(i\) we will observe.

Consistency:
We construct a snake chain \(Y_n = (X_{n+1},X_{n})\) with a state space \(S_Y = \{(i_0,i_1)\in S^2: p_{i_0i_1}>0\}\). Using the snake chain lemma along with the fact that \(\{X_n\}\) is irreducible, \(\{Y_n\}\) is also irreducible with a transition probability matrix \(\mathbf{P}_Y\) with element \[P(Y_n = (i_n,j_n)\mid Y_{n-1} = (i_{n-1}, j_{n-1})) = p_{i_nj_n} I(i_n = j_{n-1})\] and stationary distribution \(\nu(i,j) = \pi_i p_{ij}\).

Then we apply the Ergodic theorem to \(\{Y_n\}\) and \(\{X_n\}\), leading to \[\begin{align*} \frac{n_{ij}}{n} &= \frac{1}{n}\sum_{k=0}^{n-1} I(Y_k = (i,j)) \overset{a.s.}{\rightarrow} \pi_i p_{ij}\\ \frac{n_{i+}}{n} &= \frac{1}{n}\sum_{k=0}^{n-1} I(X_k = i) \overset{a.s.}{\rightarrow} \pi_i. \end{align*}\] Thus, \[\widehat{p}_{ij} = \frac{n_{ij}/n}{n_{i+}/n} \overset{a.s.}{\rightarrow} \frac{\pi_1 p_{ij}}{\pi_i} = p_{ij}.\]

Asymptotic normality:
We will use Anscombe’s lemma to prove it. Here the trick is to identify the random variable \(Z_m\) and \(W_n\). When we focus on the state \(i\), you can easily see that a natural choice of \(W_n\) is \(W_n = n_{i+}\) since it serves as an effective sample size.

How about \(Z_m\)? We want to construct \(S_{W_n} = S_{n_{i+}} = Z_1+\cdots+Z_{n_{i+}}\) to be related to the denominator \(\widehat{p}_{ij}- p_{ij} \propto n_{ij} - n_{i+}p_{ij}\). So it turns out that a possible choice is \[Z_m = \begin{cases} 1-p_{ij}\quad&\mbox{if $X_{\tau_m +1}= j$}\\ -p_{ij} \quad&\mbox{if $X_{\tau_m+1}\neq j$} \end{cases},\] where \(\tau_m\) is the \(m\)-th visit to state \(i\). Since \(X_n\) is a homogeneous Markov chain, the random variables \(Z_1,Z_2,\cdots\) are IID (recall regenerative cycle lemma). This choice leads to \[\begin{align*} \mathbb{E}(Z_m) &= (1-p_{ij})p_{ij} - p_{ij} (1-p_{ij}) = 0\\ \mathbb{E}(Z_m^2)& = (1-p_{ij})^2 p_{ij}+ p_{ij}^2 (1-p_{ij})\\ & = p_{ij}(1-p_{ij}). \end{align*}\] Moreover, \[S_{W_n} = Z_1+\cdots+Z_{n_i+} = n_{ij}(1-p_{ij}) + (n_{i+}-n_{ij}) (-p_{ij}) = n_{ij} - n_{i+}p_{ij}.\] Thus, Ancombe’s lemma implies that \[\begin{align*} \frac{S_{n_{i+}}}{\sqrt{p_{ij}(1-p_{ij})n_{i+}}} &= \frac{n_{ij}-n_{i+}p_{ij}}{\sqrt{p_{ij}(1-p_{ij}n_{i+})}}\\ & = \frac{\widehat{p}_{ij}- p_{ij}}{\sqrt{p_{ij}(1-p_{ij})/n_{i+}}} \overset{D}{\rightarrow} N(0,1). \end{align*}\]

Finally, observing that \(\frac{n_{i+}}{n\pi_i}\overset{a.s.}{\rightarrow}1\) we then conclude the result.