1.1 Markov chains on a discrete state spaces

Formally we can think of a stochastic process as a probability measure on S^\infty (the joint distribution of all X_n) but it is not convenient to put your hands on this object directly. A famous result in measure theory called Kolmogorov extension theorem say that is enough to specify all the joint distribution of X_{n_1}, X_{n_2}, \cdots, X_{n_k} for all choices of k and n_1 < n_2 \cdots < n_k.

To start, in this chapter, we assume that S is discrete (either finite or countable) and without loss of generality we can write S=\{1, 2, \cdots, N\} with N finite or not (by relabeling the state as needed). In this context we only need to specify the finite dimensional distributions \begin{aligned} &P \left( X_0= i_0, X_1=i_1, \cdots, X_n = i_n \right) =P \left( X_0=i_0\right) P\left( X_1=i_1 \vert X_0=i_0 \right)P\left( X_1=i_2 \vert X_1=i_1, X_0=i_0 \right) \cdots \\ &\cdots P\left( X_n=i_n \vert X_{n-1}=i_{n-1}, \cdots, X_0=i_0 \right) \quad \end{aligned} which we have written in terms of conditional probabilities using the product rule.

1.2 Matrix notation

For time homogeneous Markov chains it is natural to use the (possibly infinite) vector/matrix notation:

• Initial distribution \mu(i) for i \in S where \mu(i)=P(X_0=i). We have \mu(i)\ge 0 \quad \textrm{ and } \sum_{i \in S} \mu(i) =1 and \mu(i) is called a probability vector .

• Transition matrix P(i,j) for i,j \in S where P(i,j)=P(X_n=j|X_{n-1}=i). We have P(i,j) \ge 0 \quad \textrm{ and } \sum_{j \in S} P(i,j)=1 and P(i,j) is called a the transition probability matrix

• All quantities of interest for the Markov chain can be computed using these two objects. For example we have P\left\{ X_0 =i_0, X_1=i_1, X_2=i_2\right\} \,=\, \mu(i_0) P(i_0,i_1) P(i_1, i_2) \,. or P\left\{ X_2=i \right\} \,=\, \sum_{i_0\in S} \sum_{i_1 \in S} \mu(i_0) P(i_0,i_1) P(i_1, i) and so on.

1.4 Memoryless property

Markov chains forgets their past. For example if we observe that X_m=i for some i \in S then the Markov chain starts anew at i: conditional on the event \{X_m=i\}, the stochastic process Y_n=X_{m+n}, n=0,1,2,\cdots is a Markov chain with transition matrix P and initial condition i.
Indeed we have \begin{aligned} &P( X_{m+1} = i_{m+1} \cdots X_{m+n}=i_{m+n}| X_m =i ) = \frac{P (X_m=i, X_{m+1} = i_{m+1} \cdots X_{m+n}=i_{m+n})}{P(X_m=i)}\\ &= \frac{\mu P^m(i) P(i,i_{m+1}) \cdots P(i_{m+n-1},i_{m+n})}{\mu P^m(i)} = P(i,i_{m+1}) \cdots P(i_{m+n-1},i_{m}) \end{aligned}

Actually a stronger statement holds and shows that the Markov chain after time n is independent of the past!

1.5 Stationary and limiting distributions

Basic question: For Markov chain understand the distribution of X_n for large n, for example we may want to know whether the limit \lim_{n\to \infty} P( X_n=i) \,=\, \lim_{n \to \infty} \mu P^n(i) exists or not, whether it depends on the choice of initial distribution \mu and how to compute it.

Remark: Limiting distributions are always stationary distributions: If \lim_{n\to \infty} \mu P^n = \pi then \pi P \,=\, (\lim_{n \to \infty} \mu P^n ) P \,=\, \lim_{n \to \infty} \mu P^{n+1} \,=\, \lim_{n \to \infty} \mu P^{n} \,=\, \pi \,.

Remark: If \pi is stationary and is the initial condition then X_n is a sequence of identically distributed random variables.

1.6 Examples

1.7 Simulation of Markov chains

• It is relatively easy to compute the distribution of X_n by matrix multiplication, that is \mu P^n if P is not too large.

• It is also not difficult to generate the paths X_0, X_1, X_2, \cdots of the Markov chain.

  • Set X_0=i (or generate the random variable X_0 with distribution \mu)

  • Generate the RV Z with probability distribution (P(i,1), \cdots P(i,n)) and set X_1=Z.

  • and so on.

• To generate a discrete RV with finite state space do

import numpy as np

# Create an array to describe the state space
elements = np.array([1, 2, 3, 4, 5])

# Define nonuniform probabilities for each element
probabilities = np.array([0.1, 0.2, 0.3, 0.2, 0.2])

# Use random.choice with nonuniform probabilities
random_element = np.random.choice(elements, p=probabilities)

2.1 Existence of stationary distributions

Stationary distributions always exist for finite state Markov chains. This will not be the case if the state space is countable.

2.2 Uniqueness of stationary distributions

It is easy to build examples of Markov chain with multiple stationary distributions, for example P \,=\, \left( \begin{array}{cccc} \frac{1}{2} & \frac{1}{2} & 0 & 0 \\ \frac13 & \frac23 & 0 & 0 \\ 0 & 0 & \frac16 & \frac56 \\ 0 & 0 & \frac45 & \frac15 \end{array} \right) \quad \quad \pi_1=\left( \frac35 , \frac 25 ,0,0\right) \textrm{ and } \pi_2= \left( 0, 0 ,\frac{25}{49}, \frac{24}{49}\right) \quad \textrm{ are stationary} If the Markov chain starts in state 1 or 2 it will only visits states 1 and 2 in the future and so we can build a stationary distribution (by Theorem 2.1 using an initial \mu_0 which vanishes on 3 and 4 ) which will vanish on the state 3 and 4.

In an irreducible Markov chain, starting with any state the Markov chain will eventually visit any other state.

2.3 Convergence to stationary distribution

2.4 The period of a Markov chain

As we have seen before Markov chain can exhibit periodic behavior and circle around various part of the state space.

For example for random walks with periodic boundary conditions the period is 2 is the number of sites is even and 1 if the number of sites is odd.

2.5 Hitting times, return times adn the strong Markov property

We now have good understanding the asymptotic behavior of P^n(i,j) or \mu P^n as n \to \infty. When simulating or observing a Markov chain we usually observe a single sequence X_0, X_1, X_2, \cdots and our next goal is to derive a law large numbers for this sequence of random variables. To do this it will be useful to track how often the Markov chain visits some state.

Hitting time and return times are closely related and are distinct only from the way they consider what happen at time 0 (both will be useful later on).

2.6 Strong law of large numbers for Markov chains

Irreducibility means that starting from i the Markov chain will reach j and thus for any i we have P\{\tau(j) <\infty |X_0=i\} >0 If the state space is finite then we have something much stronger. Not only we have P\{\tau(j) <\infty |X_0=i\} =1, the expectation of \tau(j) is actually finite.

2.7 Exercises

3.2 Communication diagram

To figure out the class structure it is convenient to build a directed graph where the vertices are the state and each directed edges corespond to a pair (i,j) with P(i,j)>0.

For example

3.3 Decomposition of reducible Markov chains

Markov chain with closed classes R_1, \cdots R_L and transient classes denoted by T_1, \cdots T_K (set T= T_1 \cup \cdots T_K). After reordering the states the transition matrix has the block form P\,=\, \begin{matrix} R_1 \\ R_2 \\ R_3 \\ \vdots \\ R_L \\ T \end{matrix} \begin{pmatrix} P_1 & & & & &\\ & P_2 & & 0 & &\\ & & P_3 & & & \\ & 0 & & \ddots & &\\ & & & &P_L &\\ & & S & & & Q \end{pmatrix} \tag{3.1}

where P_l gives the transition probabilities within the class R_l, Q the transition within the transient classes and S\not= 0 the transition from the transient classes into the closed classes.

It is easy to see that P^n has the form P^n\,=\, \begin{matrix} R_1 \\ R_2 \\ R_3 \\ \vdots \\ R_L \\ T \end{matrix} \begin{pmatrix} P_1^n & & & & &\\ & P_2^n & & 0 & &\\ & & P_3^n & & & \\ & 0 & & \ddots & &\\ & & & &P_L^n &\\ & & S_n & & & Q^n \end{pmatrix} for some matrix S_n.

3.4 Absorption time

• Starting in a closed class R_l the Markov chain stays in R_l forever and its behavior is dictated entirely by the irreducible matrix P_l with state space R_l.

• Starting from a transient class the Markov chain will eventually exit the class T and is absorbed into some closed class. One basic question is: What is the expected time until absorption?

• We define the absorption time \tau_{\textrm{abs}} = \min\{ n \ge 0 \,;\, X_n \notin T \}

3.5 Hitting time in irreducible Markov chain

• Suppose X_n is irreducible and we are interested in computing the expected time to reach one state from another state that is E[ \tau(i)| X_0=j] \textrm{ for } j \not= i.

• If i=j then from Theorem 2.6 that E[ \tau(i)| X_0=i]=\pi(i)^{-1} where \pi is the stationary distribution.

• If i \not= j we use the following idea. Relabel the states so that the first state is i and turn i into an absorbing state
  P = \begin{pmatrix} P(i,i) & R \\ S & Q \end{pmatrix} \quad \quad \longrightarrow \quad \quad \widetilde{P} = \begin{pmatrix} 1 & 0 \\ S & Q \end{pmatrix} \tag{3.2} The fact that P is irreducible implies that S\setminus \{i\} is a transient class for the modified transition matrix \widetilde{P}. The return time \tau(i) starting from i\not= j for the Markov chain with transition P is exactly the same as the absorption time for the Markov chain with transtion \widetilde{P}. Indeed the hitting time does not depend on the submatrix R. Therefore from Theorem 3.1 we obtain immediately

3.6 Examples

Example(continued): Random walk with absorbing boundary conditions on \{0,1,\cdots,5\}. We get
Q\,=\, \begin{pmatrix} 0 & 1/2 & 0 & 0 \\ 1/2 & 0 & 1/2 & 0 \\ 0 & 1/2 & 0 & 1/2\\ 0 & 0 & 1/2 & 0 \end{pmatrix} \quad \implies \quad M=(I-Q)^{-1} \,=\, \begin{pmatrix} 1.6 & 1.2 & 0.8 & 0.4 \\ 1.2 & 2.4 & 1.6 & 0.8 \\ 0.8 & 1.6 & 2.4 & 1.2\\ 0.4 & .8 & 1.2 & 1.6 \end{pmatrix} and thus the expected times until absorption are 4 for states 1 and 4 and 6 for states 2 and 3.

Example: Random walk with reflecting boundary conditions (see Example 1.5) with N=5 and we compute E[ \tau(0) | X_0=i]. The stationary distribution is \pi = \left(\frac{1}{10}, \frac{2}{10}, \frac{2}{10}, \frac{2}{10}, \frac{2}{10}, \frac{1}{10}\right) so E[ \tau(0) | X_0=0]=10. For i\not=0 we delete the first row and column from P and have Q\,=\, \left( \begin{array}{ccccc} 0 & 1/2 & 0 & 0 &0 \\ 1/2 & 0 & 1/2 & 0 &0 \\ 0 & 1/2 & 0 & 1/2 &0 \\ 0 & 0 & 1/2 & 0 &1/2 \\ 0 & 0 & 0 & 1 &0 \end{array} \right) \,, \quad M=(I-Q)^{-1} \,=\, \left( \begin{array}{ccccc} 2 & 2 & 2 & 2 & 1\\ 2 & 4 & 4 & 4 &2 \\ 2 & 4 & 6 & 6& 3\\ 2 & 4 & 6 & 8 & 4 \\ 2 & 4 & 6 & 8 & 5 \end{array} \right) and so the expected return times to 1 are 10, 9, 16, 21, 24, 25 respectively.

3.7 Absorption probabilities

• If X_0 = i \in T belongs to some transient class and if there are more than one closed classes, say R_1, R_2, \cdots, R_L then the Markov may be absorbed in distinct closed class and so we wisht to compute the probabilities P\{ X_n \,\,\textrm{ reaches class } R_l \,|\, X_0=i \}

• Without loss of generality we can assume that each closed class is an absorbing state r_1, \cdots r_L (we can always collapse a closed class into a absorbing state since it does not matter which state in the closed class we first visit). We denote the transient states by t_1, \cdots, t_M and upon reordering the state the transition matrix has the form P\,=\, \begin{matrix} r_1 \\ \vdots \\ r_L \\ t_1 \\ \vdots \\ t_M \end{matrix} \begin{pmatrix} & & & & & \\ & I & & &0& \\ & & & & & \\ & & & & & \\ & S & & &Q & \\ & & & & & \\ \end{pmatrix} \tag{3.3}

• We wish to compute A( t_i, r_j) = P\{ X_n \, \textrm{ eventually reaches } r_j \,|\, X_0= t_i\}\,. and it will be convenient to set A(r_l, r_l)=1 and A(r_k, r_l)=0 if k\not= l.

3.8 Tennis game

We model tennis as a Markov chain: For each point player A will win with probability p and player B will be win with probability q=1-p and we assume all points are independent. Winning a game leads to the Markov chain with following communication diagram. This problem illustrates the power of compunding. For example my fellow countryman Roger Federer won 54% of all points played during his career. But he won 80% all of his game.

3.9 Gambler’s ruin

• A gambler is coming to the casino to play the game of dice called craps. At this game the standard bet (pass line bet) has even money odds and the probability to win is \frac{244}{495}=0.49292929..., one of the casino games with best odds.

• Suppose that the player start with a fortune of j and bets 1 in every game. Their ultimate goal is to reach a fortune of N and, if they do, they would stop. Of course they could lose all their money before reaching their goal.

• The game is described by a random walk on \{0,1, \cdots, N\} with absorbing boundary conditions and we want to compute \alpha_j = P(X_n \textrm{ reaches } N \textrm{ before reaching } 0| X_0=j)

• Conditioning on the first step gives the second order homogeneous linear difference equation \alpha_j = p \alpha_{j+1} + (1-p) \alpha_{j-1} \quad \textrm{ with boundary conditions } \alpha_0=0\,, \alpha_N=1

• To solve this second order homogeneous linear difference equations we look for solutions of the form \alpha_j=x^j which leads to the equation px^2 - x + (1-p) =0 \quad \implies x_1=1\,, x_2= \frac{1-p}{p} \tag{3.4} This leads to the system of equation (for p\not= \frac12) \alpha_j= C_1 + C_2 \left(\frac{1-p}{p} \right)^j \quad \textrm{ with } \quad \left\{ \begin{array}{r} C_1+C_2=0 \\ C_1 + C_2\left(\frac{1-p}{p} \right)^N=1 \end{array} \right.

4.1 Examples

We introduce some basic examples of Markov chains on a countable state space

4.2 Transience/recurrence dichotomy

In Markov chains on a countable state space a new phenomenon occurs compared to finite state space. Suppose the the Markov chain is irreducible (or starts in a closed class), from a state i the Markov chain will return to i with positive probaility but it is also possible that the Markov chain does not return to i and “wander away to infinity”.

We introduce the corresponding definitions of transience and recurrence of a state. Recall that the return time to state i is given by
\tau(i) = \min \{ n \ge 1 \,;\, X_n =i \}\,. \quad (\text{return time})

4.3 Transience/recurrence for the simple random walk.

We analyze the recurence and transience propeties for the simple random walk on \mathbb{Z}^d (see Example 4.5).
As we will see this depends on the dimension. To prove recurrence transience here we compute/estimate directly \sum_{n=0}^\infty P^n(0,0) = \sum_{n=0}^\infty P^{2n}(0,0) since X_n is periodic with period 2.

d=1: To return to 0 in 2n steps the Markov chain must take exactly n steps to the left and n steps to the right and thus we have P^{2n}(0,0) = {2n \choose n} \frac{1}{2^{2n}} By Stirling’s formula we have n! \sim \sqrt{2 \pi n}e^{-n} n^n where a_n \sim b_n means that \lim a_n/b_n =1. Thus we have {2n \choose n} \frac{1}{2^{2n}} \sim \frac{1}{2^{2n}} \frac{ \sqrt{2 \pi 2n} e^{-2n} (2n)^{2n}} { 2 \pi n e^{-2n} n^{2n}} = \frac{1}{\sqrt{\pi n}}\,.
Recalling that if a_n \sim b_n then \sum a_n converges if and only if \sum b_n converges we see that the random walk in d=1 is recurrent.

4.4 Transience/recurrence for the succcess run chain

Continuing with Example 4.4 we consider the return time to state 0, \tau(0) whose pdf we can compute explicitly since to return to 0 the only possible paths are 0 \rightarrow 0, 0 \rightarrow 1 \rightarrow 0, 0 \rightarrow 1 \rightarrow 2 \rightarrow 0, and so on. We set u_n \equiv p_0p_1 \cdots p_{n-1} and find P( \tau(0) = k|X_0=0) = p_0 p_1 p_{k-2}q_{k-1}= p_0 p_1 p_{k-2}(1-p_{k-1}) = u_{k-1} - u_k\,. Therefore P( \tau(0) \le n |X_0=0) = \sum_{k=1}^n P( \tau(0) = k|X_0=0) = (1- u_0) + (u_0-u_1) + \cdots + (u_{n-1} - u_n) = 1 - u_n and P( \tau(0) < \infty|X_0=0) = 1 - \lim_{n\to \infty} u_n and we obtain that \textrm{The success run chain is recurrent if and only if } \lim_{n \to \infty} u_n = \lim_{n \to \infty} p_0 \cdots p_{n-1} =0 To have a better handle on this criterion we need a little result from analysis about infinite products.

4.5 Another criterion for transience

We add one more method to establish transience. For this pick reference state j and consider the hitting time to state j, \sigma(j) (not the return time but we will play with both.) \alpha(i)\,=\, P \left\{ \sigma(j)< \infty | X_0= i \right\} By definition we have \alpha(j)=1 since \sigma(j)=0 if X_{0} = j. If the chain is transient we must have \alpha(i) < 1 \quad \textrm { for } i \not= j \,. since by Theorem 4.1 for i \not=j we have \alpha(i)= P \left\{ \tau(j) < \infty |X_0= i \right\} < 1.

Let us derive an equation for \alpha(i) by conditioning of the first step. For i \not=j \begin{aligned} \alpha(i) &= P \left( \sigma(j) < \infty | X_0= i \right) = P \left( \tau(j) < \infty | X_0= i \right) \\ &= \sum_{k \in S} P \left( \tau(j) < \infty | X_1=k \right) P(i,k) = \sum_{k \in S} P \left( \sigma(j)< \infty | X_0=k \right) P(i,k) \\ &= \sum_{k \in S} P(i,k) \alpha(k) \end{aligned} and thus \alpha(i) satisfies the equation P\alpha(i) = \alpha(i) \,, \quad i \not= j \,.

4.6 Transience/recurrence for the RW on \{0,1,2,\cdots \}

Continuing with Example 4.1 we use Theorem 4.2. We pick 0 as the reference state and for j \not=0 solve the equation P\alpha(j) = P(j, j-1) \alpha(j-1) + P(j, j+1) \alpha(j+1) = (1-p) \alpha(j-1) + p \alpha(j+1) = \alpha(j) whose solution is (like for the Gambler’s ruin problem) \alpha(j) = \left\{ \begin{array}{cl} C_1 + C_2 \left( \frac{1-p}{p}\right)^j & \textrm{ if } p \not= \frac{1}{2} \\ C_1 + C_2 j & \textrm{ if } p = \frac{1}{2} \end{array} \right. \,. Using that \alpha(0)=0 we find \alpha(j) = \left\{ \begin{array}{cl} (1-C_2) + C_2 \left( \frac{1-p}{p}\right)^j & \textrm{ if } p \not= \frac{1}{2} \\ (1-C_2) + C_2 j & \textrm{ if } p = \frac{1}{2} \end{array} \right. \,. and we see the condition 0 < \alpha(i) < 1 is possible only if (1-p)/p < 1 (that is p > 1/2) and by choosing C_2=1. Thus we conclude \textrm{ The random walk on } \{0,1,2, \cdots\} \textrm{ is } \left\{ \begin{array}{l} \textrm{transient for } p> \frac{1}{2} \\ \textrm{recurrent for } p \le \frac{1}{2} \end{array} \right. \,.

5.1 Positive recurrence versus null recurrence

A finite state irreducible Markov chain is always recurrent, P(\tau(i) < \infty |X_0=j)=1 and we have proved Kac’s formula for the invariant measure \pi(i) = E[ \tau(i)|X_0=i]^{-1}, that is the random variable \tau(i) has finite expectation.
For a countable state space it is possible for a Markov chain to be recurrent but also that \tau(i) does not have finite expectation. This motivates the following definitions.

We first investigate the relation between recurrence and existence of invariant measures. We first show that if one state j is positive recurrent then there exists a stationary distribution. The basic idea is to decompose any path of the Markov chain into successive visits to the state j. To build up our intuition if a stationary distribution were to exists it should measure the amount of time spent in state i and to measure this we introduce \mu(i) = E\left[ \sum_{n=0}^{\tau(j)-1} {\bf 1}_{\{X_n=i\}}| X_0=j \right] \quad = \textrm{number of visits to $i$ between two successive visits to $j$}.

Note that \mu(j)=1 since X_0=j, and if j is positive recurrent \sum_{i} \mu(i) = \sum_{i\in S} E\left[ \sum_{n=0}^{\tau(j)-1} {\bf 1}_{\{X_n=i\}}| X_0=j \right] \,=\, E\left[ \tau(j) \,|\, X_0=j \right] < \infty

5.2 Stationarity and irreducibility implies positive recurrence

5.3 Ergodic theorem for countable Markov chains

5.4 Convergence to equilibrium via coupling

We continue our theoretical consideration by proving that if the chain is aperiodic then the distribution of X_n converges to \pi(j).

5.5 Examples

Positive recurrence for the random walk on \{0,1,2, \cdots\}. Continuing with Example 4.1, we use Theorem 5.3 to establish positive recurrence by computing the stationary distribution. We find the equations \pi(0) (1-p) + \pi(1)(1-p) = \pi(0) \quad \textrm{ and } \quad \pi(j+1)(1-p) + \pi(j-1)p = \pi(j) \quad \textrm{ for }j \ge 1 The second equation has the general solution \pi(n)=C_1 + C_2 \left(\frac{p}{1-p}\right)^n \,\,\,(\textrm{if }p\not=\frac12) \quad \textrm{ and } \quad \pi(n)=C_1 + C_2n\,\,\, (\textrm{ if } p=\frac12) and the first equation gives \pi(1)=\frac{p}{1-p}\pi(0). For p=\frac{1}{2} we cannot find a normalized solution. For p < \frac{1}{2} we can choose C_1=0 and \pi(n) = \left(\frac{p}{1-p}\right)^n \pi(0) can be normalized toi find \pi(n) = \frac{1-2p}{1-p} \left(\frac{p}{1-p}\right)^n. If p> \frac{1}{2} we already know the Markov chain is transient and thus \textrm{ The random walk on } \{0,1,2, \cdots\} \textrm{ is } \left\{ \begin{array}{l} \textrm{transient for } p> \frac{1}{2} \\ \textrm{null recurrent for } p = \frac{1}{2} \\ \textrm{positive recurrent for } p < \frac{1}{2} \\ \end{array} \right. \,.

6.1 Moment generating function

6.2 Branching process

6.3 Examples

• If p(0)=\frac14, p(1)=\frac14, p(2)=\frac12 \implies \phi(s)=\frac14 + \frac14 s + \frac12 s^2 then \mu=\frac{5}{4}>1 and solving \phi(a)=a gives a=1,\frac12 so the extinction probability is 1/2.

• If p(0)=\frac14, p(1)=\frac12, p(2)=\frac11 \implies \phi(s)=\frac14 + \frac12 s + \frac14 s^2 then \mu=1 and solving \phi(a)=a gives a=1 so the extinction probability is 1.

• If p(0)=\frac12, p(1)=\frac14, p(2)=\frac14 \implies \phi(s)=\frac12 + \frac14 s + \frac14 s^2 then \mu=\frac{3}{4}<1 and solving \phi(a)=a gives a=1,2 so the extinction probability is 1.

6.4 Repair shop example

Recall the Markov chain in Example 4.3. To establish if the system is postive recurrent we try to solve \pi P=\pi and find the system of equations \begin{aligned} \pi(0)&= \pi(0)a_0 + \pi(1)a_0 \\ \pi(1)&= \pi(0)a_1 + \pi(1) a_1 + \pi(2)a_0 \\ \pi(2)& = \pi(0)a_2 + \pi(1)a_2 +\pi(2) a_1 + \pi(3)a_0 \\ & \vdots \\ \pi(n)& = \pi(0)a_n + \sum_{j=1}^{n+1} \pi(j) a_{n+1-j} \end{aligned}

We solve this using the generating functions \psi(s)=\sum_{n=0}^\infty s^n \pi(n) and \phi(s)=\sum_{k=0}s^k a_k: \begin{aligned} \psi(s) = \sum_{n=0}^\infty s^n\pi(n) &= \pi(0) \sum_{n=0}^\infty s^n a_n + \sum_{n=0}^\infty s^n \sum_{j=1}^{n+1} \pi(j) a_{n+1-j}\\ & = \pi(0) \sum_{n=0}^\infty s^n a_n + s^{-1} \sum_{j=1}^\infty \pi(j) s^j \sum_{n=j-1}^\infty s^{n+1-j} a_{n+1-j} \\ &=\pi(0) \phi(s) + s^{-1}(\psi(s)-\pi(0)) \phi(s) \end{aligned}

7.1 Balance equation

Consider a Markov chain with transition probabilities P(i,j) and stationary distribution \pi(i). We can rewrite the equation for stationarity, \pi(i)=\sum_{j} \pi(j) P(j,i) as
\sum_{j} \pi(i) P(i,j) \,=\, \sum_{j} \pi(j) P(j,i) \,. \tag{7.1} which we are going to interpret as a balance equation.

We introduce the (stationary) probability current from i to j as J(i,j) \equiv \pi(i) P(i,j) \tag{7.2} and Equation 7.1 can be rewritten as \underbrace{\sum_{ j} J(i,j)}_{\textrm{total current out of state } i} \,=\, \underbrace{\sum_j J(j,i) }_{\textrm{total current into of state i}} \, \quad \quad \textrm{ balance equation } \tag{7.3} i.e., to be stationary the total probability current out of i must be equal to the total probability current into i.

7.2 Detailed balance

A stronger condition for stationarity can be expressed in terms of the balance between the currents J(i,j). A Markov chain X_n satisfies detailed balance if there exists \pi(i)\ge 0 with \sum_{i} \pi(i)= 1 such that for all i,j we have \pi(i) P(i,j) \,=\, \pi(j) P(j,i) \quad \textrm{ or equivalently } \quad J(i,j) = J(j,i) \tag{7.4} This means that for every pair i,j the probability currents J(i,j) and J(j,i) balance each other.
Clearly Equation 7.4 is a stronger condition than Equation 7.3 and thus we have

But it is easy to see that detailed balance is a stronger condition than stationarity. The property of detailed balance is often also called (time)-reversibility since we have the following results which states that the probability of any sequence is the same as the probability of the time reversed sequence.

7.3 Examples

8.1 MCMC

Suppose you are given a certain probability distribution \pi on a set S and you goal is to generate a sample from this distribution. The Monte-Carlo Markov chain (MCMC) method consists in constructing an irreducible Markov chain X_n whose stationary distribution is \pi. Then to generate \pi one simply runs the Markov chains X_n long enough such that it is close to its equilibrium distribution. It turns out that using the detailed balance condition is a very useful tool to construct the Markov chain in this manner.

A-priori this method might seems an unduly complicated way to sample from \pi. Indeed why not simply simulate from \pi directly?

To dispel this impression we will consider some concrete examples but for example we will see that often one want to generate a uniform distribution on a set S whose cardinality |S| might be very difficult.

A large class of model are so-called “energy models” which are described by an explicit function f:S \to \mathbb{R} interpreted as the energy (or weight) of the state. Then one is interested in the stationary distribution \pi(i) = Z^{-1} e^{-f(i)} \quad \textrm { where } Z= \sum_{i}e^{-f(i)} \textrm{ is a normalization constant} As we will see computing the normalization constant Z can be very difficult and so one cannot directly simulate from \pi!

The measure \pi assign biggest probability to the states i where f(i) is the smallest, that is they have the smallest energy.
Often (e.g. in economics) the measure is written with a different sign convention \pi(i) = Z^{-1}e^{f(i)} which now favor the states with the highest values of f.

8.2 Proper q-coloring of a graph

For a graph G=(E,V) a proper q-coloring consists of assigning to each vertex v of the graph one of q colors subject to the constraint that if 2 vertices are linked by an edge they should have different colors. The set S' of all proper q-coloring which is a subset of S \,=\, \left\{ 1, \cdots, q\right\}^V. We denote the elements of S by \sigma = \{ \sigma(v)\}_{v \in V} with \sigma(v) \in \{1, \cdots, q\}. The uniform distribution on all such proper coloring is \pi(\sigma) = 1/|S'| for all \sigma \in S'. Even for moderately complicated graph it can be very difficult to compute to |S'|.

Consider the following transition rules:
1. Choose a vertex v of at random and choose a color a at random. 2. Set \sigma'(v) = a and \sigma'(w)=\sigma(w) for w \not=v. 3. If \sigma' is a proper q-coloring then set X_{n+1}=\sigma'. Otherwise set X_{n}=\sigma.

The transition probabilities are then given by \begin{aligned} P(\sigma, \sigma') =\left\{ \begin{array}{cl} \frac{1}{q |V|} & \textrm{ if } \sigma \textrm{ and } \sigma' \textrm{ differ at exactly one vertex} \\ 0 & \textrm{ if }\sigma \textrm{ and } \sigma' \textrm{ differ at more than one vertex} \\ 1 - \sum_{\sigma' \not = \sigma} P(\sigma, \sigma') & \textrm{ if } \sigma'=\sigma \end{array} \right. \end{aligned} Note that P(\sigma, \sigma) is not known explicitly either but is also not used to run the algorithm. The cardinality |S'| is also not needed.

In order to check that the uniform distribution is stationary for this Markov chain it note that P(\sigma, \sigma') is symmetric matrix. If one can change \sigma into \sigma' by changing one color then one can do the reverse transformation too.

8.5 Uniform distribution on a graph

Suppose we are given a graph G=(E,V). Often in application to computer networks or social networks the graph is not fully known. Only local information is available, given a vertetx one knows the vertices one is connected. Nonetheless you can run a random walk on the graph by choosing a edge at random and moving to the corresponding vertex, i.e.  Q(v,w) = \frac{1}{\textrm{ deg}(v)} \textrm{ if } v \sim w and the stationary distribution is \pi(v) \propto deg(v).

Suppose you wish to generate a uniform distribution on the graph. Then we can use the Metropolis-Hasting to generate a Markov chain with a uniform stationary distribution. We have \alpha(v,w) = \min\left\{ 1, \frac{Q(w,v)}{Q(v,w)}\right\} =\min\left\{ 1, \frac{\textrm{ deg}(v)}{\textrm{ deg}(w)}\right\} \textrm{ if } v \sim w and thus the transition matrix P(v,w) = Q(v,w) \alpha(v,w) = \frac{1}{\textrm{ deg}(v)} \min\left\{ 1, \frac{\textrm{ deg}(v)}{\textrm{ deg}(w)}\right\} = \min\left\{ \frac{1}{\textrm{ deg}(v)}, \frac{1}{\textrm{ deg}(w)} \right\} generates a uniform distribution on a graph.

8.7 Glauber algorithm

Another algorithm which is widely used for Monte-Carlo Markov chain is the Glauber algorithm which appear in the literature under a variety of other names such as Gibbs sampler in statistical applications, logit rule in economics and social sciences, heat bath in physics, and undoubtedly under various other names.

The Glauber algorithm is not quite as general as the Metropolis algorithm. Indeed we assume that the state space S has the following structure S \subset \Omega^V where both \Omega and V are finite sets. For example S \subset \{0,1\}^m in the case of the knapsack problem or S \subset \{ 1, \cdots, q\}^V for the case of the proper q-coloring of a graph. We denote by \sigma \,=\, \{ \sigma(v) \}_{v \in V} \,, \quad \sigma(v) \in \Omega \,. the elements of S.

It is useful to introduce the notation \sigma_{-v} \,=\, \{ \sigma(w) \}_{w \in V, w\not =v} and we write \sigma \,=\, ( \sigma_{-v}, \sigma(v) ) \,. to single out the v entry of the vector \sigma.

8.8 Ising model on a graph

Let G=(E,V) be a graph and let S= \{-1,1\}^V. That is to each vertex assign the value \pm 1, you can think of a magnet at each vertex pointing either upward (+1) or downward (-1). To each \sigma \in S we assign an “energy” H(\sigma) given by H(\sigma) \,=\, - \sum_{ e=(v,w) \in E} \sigma(v) \sigma(w) \,. The energy \sigma is minimal if \sigma(v) \sigma(w) =1 i.e., if the magnets at v and w are aligned. Let us consider the probability distribution \pi_\beta (\sigma) \,=\, \frac{e^{-\beta H(\sigma)} }{Z_\beta} \,, \quad Z_\beta \,=\, \sum_{\sigma} e^{-\beta H(\sigma)} \,. The distribution \pi_\beta is concentrated around the minima of H(\sigma). To describe the Glauber dynamics note that H ( \sigma_{-v}, 1) - H( \sigma_{-v}, -1) \,=\, -2 \sum_{ w\,;\, w \sim v} \sigma(w) and this can be computed simply by looking at the vertices connected to v and not at all the graph. So the transition probabilities for the Glauber algorithm are given by picking a vertex at random and then updating with probabilities \frac{ \pi( \sigma_{-v}, \pm 1)} { \pi( \sigma_{-v}, 1) + \pi( \sigma_{-v}, -1) } \,=\, \frac{1}{ 1 + e^{\pm \beta \left[ H ( \sigma_{-v}, 1) - H( \sigma_{-v}, -1)\right]}} \,=\, \frac{1}{1 + e^{\mp 2 \beta \sum_{ w\,;\, w \sim v} \sigma(w)}} \,.

8.9 Simulated annealing

8.10 Parallel tempering

9.1 Total variation norm

• Given two probability measure \mu and \nu on S the total variation distance between \mu and \nu is given by \|\mu - \nu\|_{\textrm{TV}}= \sup_{A \subset \Omega}|\mu(A) - \nu(A)| that is the largest distance between the measures of sets.

• The supremum over all set is not convenient to compute but we have the formula

9.2 Total variation and coupling

9.3 Coupling of Markov chains

• A coupling of a Markov chain is a stochastic process (X_n,Y_n) with state space S\times S such that both X_n and Y_n are Markov chains with transition matrix P but with possibly different initial conditions.

• A Markovian coupling is a coupling such that the joint (X_n,Y_n) is itself a Markov chain with some transition matrix Q((i,j), (k,l)) which must then satisfy \sum_{k} Q((i,j), (k,l)) = P(j,l) \quad \textrm{ and } \sum_{l} Q((i,j), (k,l)) = P(i,k) \,.

• Recall that Theorem 2.3 where we proved the convergence of aperiodic irreducible Markov chains to their stationary distribution we used an independent coupling. We will expand upon this idea by considering other more interesting couplings.

• The coupling time is defined by \sigma= \inf \{ n \ge 0 \,:\, X_n=Y_n \} that is, it is the first time that the Markov chain visit the same state.

• After the coupling time \sigma, X_n and Y_n have the same distribution (by the strong Markov property) so we can always modified a coupling so that, after time \sigma, X_n and Y_n are moving together, i.e. X_n=Y_n for all n \ge \sigma. We will always assume this to be true in what follows.

9.4 Speed of convergence via coupling

We can use this result to bound the distance to the stationary measure

9.5 Mixing time for the sucess run chain

We start with a countable state space example, the success run chain which is very easy to analyze.

Parenthetically, it should be noted that the supremum over i in d(n) is often not well suited for countable state space. It may often happen that the number of steps it take to be close to the stationary distribution may depend on where you start.

We consider a special case of Example 4.4 with constant succes probability. P(n,0)=(1-p)\,, \quad P(n,n+1)= p\,, \quad n=0,1,2,\cdots

Suppose X_0=i and Y_0=j (with i\not=j) we couple the two chains by moving them together.

• Pick a random number U, if U \le (1-p) then set X_1=Y_1=0 the coupling time \sigma=1.

• If U \ge 1-p then move X_1=i+1, Y_1=j+1.

Clearly we have P( \sigma >n| X_0=i, Y_0=j) = p^n and thus we find for the mixing time
d(n)=p^n < \varepsilon \iff n \ge \frac{\ln(p)}{\ln \varepsilon}

9.6 Mixing time for the random walk on the hypercube

Recall that for the random walk on the hypercube (see Example 1.4) a state is a d-bit \sigma=(\sigma_1, \cdots, \sigma_d) with \sigma_j\in \{0,1\}.

The Markov chain is periodic with period 2 so to make it aperiodic we consider its “lazy” version in which we do not move with probability 1/2, that is we consider instead the transition matrix \frac{P+I}{2} which makes it periodic. The Markov chain moves then as follows, we pick a coordinate \sigma_j at random and then replace by a random bit 0 or 1.

To couple the Markov chain we simply move move them together: If X_n=\sigma and Y_n=\sigma' we pick a coordinate j \in \{1,\cdots, d\} at random and then replace both \sigma_j and \sigma'_j by the same random bit. This is a coupling and after a move the chosen coordinates coincide.

Under this coupling X_n and Y_n will get closer to each other if we select a j such that \sigma_j \not= \sigma'_j and we will couple when all the coordinates have been selected. The distribution of the coupling time is thus exactly the same as the time \tau need to collect d toys in the coupon collector problem.

We now get a bound on the tail as follows. We denote by A_i the events that coordinate i has not been selected after m steps. We have, using the inequality (1-x)\le e^{-x} \begin{aligned} P\{\sigma > n\} = P\left(\bigcup_{i=1}^d A_i\right) \le \sum_{i=1}^d P(A_i) = \sum_{i=1}^d \left(1- \frac{1}{d}\right)^n \le d e^{-\frac{n}{d}} \end{aligned} So if we pick n = d\ln(d) + c d we find P\{\sigma > d\ln(d) + c d\} \le d e^{-\frac{d \ln(d) + cd}{d}}= e^{-c} \implies t_{\textrm{mix}}(\varepsilon) \le d\ln(d) + \ln(\varepsilon^{-1})d\,.