Markov Chains

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2006

Introduction

Suppose that a system has n possible states and that the probability vector {(Subscript[x, i])} contains a list of probabilities for the occurrence of these states at a particular time.  At a later time, or the next computation cycle, the probability vector becomes y = T . x where T is the transition matrix consisting of the probabilities (Subscript[T, i, j]) for a transition from state j to state i during the specified interval.  If the transition matrix is independent of time, the evolution of the system is described as a Markov chain consisting of the sequence of probability vectors (Subscript[y, 1]), (Subscript[y, 2]), ⋯, (Subscript[y, n]) where (Subscript[y, n]) = (T^n) x.  All elements of each probability vector and of the transition matrix must be nonnegative real numbers and the sum of the elements in any probability vector and in each column of the transition matrix must be unity.

Suppose that for a given T there exists a probability vector x such that x = T . x.  We then call x an equilibrium vector of T.  Every transition matrix has at least one equilibrium vector.  We will use two important theorems about equilibrium vectors.

Theorem 1: If T is the transition matrix for a Markov chain, then λ=1 is an eigenvalue of T.  If x is an eigenvector with λ=1 and nonnegative components, then x/(∑_i) (Subscript[x, i]) is an equilibrium vector of T.

Theorem 2: If T is the transition matrix for a Markov chain and there exists some power of T for which every element is nonzero, then T has one unique equilibrium vector.

These theorems provide a simple algorithm for construction of the equilibrium vector known as the power method.  If we assume that there exists a unique equilibrium vector, we can construct that vector by choosing a random probability vector x and evaluating the limit of (T^n) x as x -> ∞.  Otherwise, if we evaluate the limit of (T^n) directly its transpose provides a list of equilibrium vectors.  However, this list is usually highly degenerate.

Initialization

ClearAll["Global`*"];
Off[General::spell,  General::spell1]

Implementation

Function for construction of equilibrium vectors

The function EquilibriumVector implements two closely related algorithms for the construction of equilibrium vectors for a Markovian transition matrix.  If we assume that the equilibrium vector is unique, a random probability vector x is constructed and the limit of (T^n) x is evaluated using the rms deviation as a measure of convergence.  Otherwise, the limit of (T^n) is evaluated and its transpose is returned as a list of equilibrium vectors which might contain many duplications.

\(\(\(ClearAll[EquilibriumVector];\)\(
\)
    \)\)
\(\(\(EquilibriumVector::usage = "EquilibriumVector[T,options] constructs the equilibrium vector(s) for a Markov transition matrix T.  The equilibrium vector is assumed to be unique by default, but with option Unique->False the transpose of the limiting power of T is returned.  Additional options include the maximum number of iterations, MaxIterations, and the convergence criterion, Tolerance.";\)\(
\)
    \)\)
EquilibriumVector::unsquare = "object with dimensions `1` is not \
square";
EquilibriumVector::badvalues = "illegal values in transition matrix\
";
EquilibriumVector::badcolumn = "unnormalized column";
EquilibriumVector::fail = "failed";
\(\(\(EquilibriumVector::iterations = "Warning: failed to converge in `1` iterations";\)\(
\)
    \)\)
\(\(\(Options[
          EquilibriumVector] = {Tolerance -> 10\ $MachineEpsilon, 
          MaxIterations -> $RecursionLimit, 
          Unique -> True};\)\(
\)
    \)\)
\(EquilibriumVector[T_?MatrixQ, opts___] := 
      With[{defaults = Options[EquilibriumVector]}, 
        Module[{n = 0, m, error = Infinity, U, V, goof = False, tolerance, 
            nmax, result}, 
          m = Dimensions[T]; 

 (*\ 
            Verify\ that\ T\ is\ a\ proper\ transition\ \(\(matrix\)\(.\)\)\ *) 


          If[Not[MatchQ[m, {k_, k_}]], \ 
            Message[EquilibriumVector::unsquare, m]; goof = True]; 

          If[Count[Flatten[T], \((x_ /; x < 0 || x > 1)\)] > 0, 
            Message[EquilibriumVector::badvalues]; goof = True]; 

          If[Count[\((1 - Plus @@ # &)\)\  /@ Transpose[T], 
                x_ /; x > tolerance] > 0, 
            Message[EquilibriumVector::badcolumn]; goof = True]; 

 (*\ 
            Initialize\ iteration\ parameters\ and\ \(\(options\)\(.\)\)*) 

tolerance = \(Tolerance /. {opts}\) /. 
              
              defaults; 
nmax = \(MaxIterations /. {opts}\) /. defaults; 
unique = \(Unique /. {opts}\) /. 
              defaults; 
m = 
            m[[1]]; 
Which[goof, 
            Message[EquilibriumVector::fail], 

 (*\ 
              Construct\ dominant\ eigenvector, \ assumed\ to\ be\ unique, \ 
              from\ random\ initial\ \(\(state\)\(.\)\)\ *) 

unique, 
V = 
              Table[Random[Real, {0, 1}], {m}]; 
V = 
              V/Plus @@ V; 
While[
              error ≥ tolerance && n ≤ nmax, U = V; 
              V = T . V; error = Sqrt[Plus @@ \((\((V - U)\)\^2)\)/m]; 
              n = n + 1]; \ 
            result = V, 
				
 (*\ 
              Evaluate\ limit\ of\ T\^n\ *) 

True, 
            U = T; 
While[
              error ≥ tolerance && n ≤ nmax, V = U; 
              U = T . U; error = Sqrt[Plus @@ \((Flatten[V - U]\^2)\)]/m; 
              n = n + 1]; 
result = Transpose[V]]; 

 (*\ 
            Return\ \(\(result\)\(.\)\)\ *) 

If[n ≥ nmax, 
            Message[EquilibriumVector::iterations, n]]; result]]\)

?EquilibriumVector

Test cases

The following function constructs a random transition matrix.

randomT[n_] := 
      Transpose[
        Map[#/Plus @@ # &,  
          Table[Random[Real, {0, 1}], {n}, {n}]]]

testMatrix = randomT[100];

Below we compare equilibrium vectors constructed using both versions of our algorithm with the result of the built-in Eigensystem function.  Note that we must normalize the latter result.

Timing[result1 = EquilibriumVector[testMatrix];]

Timing[result2 = 
          EquilibriumVector[testMatrix, Unique -> False];]

Sqrt[(Plus @@ (((result1 - First[result2]))^2))/
        Length[result1]]

Timing[result3 = (#/Plus @@ # &)@Rest[
                Eigensystem[testMatrix]][[1, 
              1]];]

Sqrt[(Plus @@ (((result1 - result3))^2))/
        Length[result1]]

Thus, we find good numerical agreement between the three methods.  The built-in function is the most efficient, producing the entire eigensystem in about the same time as a single equilibrium vector by the power method.  Evaluating the limiting power of the transition matrix directly is rather slow.

Evaluate all equilibrium vectors of (({{1/3}, {1/3}, {1/3}})({{0}, {0}, {1}})({{1/2}, {1/2}, {0}})).

Evaluate all equilibrium vectors of (({{1/3}, {2/3}, {0}})({{0}, {0}, {1}})({{1/2}, {1/2}, {0}})) and compare with the previous exercise.  Justify the claim that if a transition matrix can be arranged in block diagonal form by renumbering the state space, the number of equilibrium vectors will be at least as large as the number of blocks.