Markov Chains

MarkovQ

Source of statement on which code is based: Mårten Blix, 1997, “Rational Expectations in a VAR with Markov Switching”, Sveriges Riksbank Working Paper No. 42, October 1997, p.4.

MarkovQ::usage  =  "MarkovQ[x_MatrixQ] tests whether the supplied \
matrix x is of an appropriate form to be a matrix of Markov transition \
probabilities. Specifically, it tests that the matrix is square, the elements \
are numerical and non-zero and that it returns a column of ones when \
multiplied by a column of ones (ie, that the last transition probability in a \
column equals 1 minus the sum of all the other probabilities in the column so \
that the probabilities are exhaustive)."

MarkovQ::notsq  =  "This matrix is not square."

MarkovQ::negprob  =  "This matrix contains negative or complex \
numbers which cannot be probabilities."

\(MarkovQ[x_?MatrixQ] := 
	
      If[Equal\  @@ \ Dimensions[x]\  (*\ matrix\ is\ square*) , 
		If[
          And\  @@ \ \((NonNegative\  /@ \ Flatten[\ x])\)\  (*\ 
            matrix\ only\ contains\ non - negative\ elements*) , 
		If[
            Table[1, {Length[x]}] . x\  == Table[1. , {Length[x]}], True, \ 
            False], Message[MarkovQ::negprob]], 
        Message[MarkovQ::notsq]]\)

MarkovErgodicProbabilities

Source of statement on which code is based: Mårten Blix, 1997, “Rational Expectations in a VAR with Markov Switching”, Sveriges Riksbank Working Paper No. 42, October 1997, p.4.

For a given matrix of Markov transition probabilities, P, finds the column vector 𝓅 such that 
P.𝓅 = 𝓅, bearing in mind that the elements of 𝓅 must sum to one. Depends on the MarkovQ function as it tests the matrix to ensure it is a proper matrix of Markov transition probabilities.

MarkovErgodicProbabilities::usage  =  
      "For a given matrix of transition probabilities P, gives the vector \
of stationary probabilities, ie the vector p such that P.p = p."

MarkovErgodicProbabilities[x_?MatrixQ] := 
      If[MarkovQ[x], 
        With[{n = Length[x], vars = Table[z[i], {i, 1, Length[x]}]}, 
          Flatten[vars  /.  
              Solve[Join[x, {Table[1, {n}]}] . vars  ==  Join[vars, {1}], 
                vars]]],    
        Message[MarkovErgodicProbabilities::notmark],    
        Message[MarkovErgodicProbabilities::notmark]]

MarkovErgodicProbabilities::notmark = 
      "Can't find ergodic (stationary) probabilities because this is not a \
properly constructed matrix of Markov transition probabilities."

MakeMarkovChain

This is based on code by Pedro J F de Lima for GAUSS, translated into Mathematica.  Depends on the MarkovErgodicProbabilities function.

Needs["Statistics`Common`DistributionsCommon`"]

Needs["Statistics`ContinuousDistributions`"]

?Statistics`Common`DistributionsCommon`*

cumulcounts::notsort = "Vector x should be in ascending order."

cumulcounts::domain  =  
      "Scalar Y isn't within the bounds defined by the elements of X."

\(cumulcounts[x_?VectorQ, y_Real] := 
	
      If[Sort[x] === x\ \ , 
			
        If[y <= Max[x], 
          \(Flatten[Position[x, Min[Select[x, # >= y\ &]]]]
              \)[[1]], 
          Message[cumulcounts::domain]], 
Message[cumulcounts::notsort]]\)

\(makeMarkovChain[n_Integer, \ P_?MatrixQ]\  /; \ MarkovQ[P] := 
	
      Module[{ergodic, cumergodic, \ u, s}, 
		
        ergodic\  = MarkovErgodicProbabilities[P]; 
		
        cumergodic = 
          FoldList[#1\  + \ #2\ &, First[ergodic], Rest[ergodic]]\ ; 
		
        u\  = Table[Random[], {n - 1}]; 
		
        s\  = \ FoldList[
            cumulcounts[
                FoldList[#1\  + \ #2\ &, 
                  First[P[[#1]]], 
                  Rest[P[[#1]]]], #2]
                \ \ &, cumulcounts[cumergodic, Random[]], u]; 

        		{s, ergodic}]\)

Table[Random[], {100000}]; // Timing

testing these functions

test  =  {{0.5, 0.25, 0.25}, {0.25, 0.5, 0.25}, {0.25, 0.25, 0.5}};

test2  = 
       {{0.5, 0.25, 0.25}, {0.25, -0.5, 0.25}, {0.25, 0.25, 0.5}};

test3  =  {{0.5, 0.25, 0.25}, {0.25, 15, 0.25}, {0.25, 0.25, 0.5}};

test4  =  {{0.5, 0.25, 0.25}, {0.25, 0.5, 0.25}};

FoldList[#1  +  #2 &, 
      First[test[[2]]], 
      Rest[test[[2]]]]

makeMarkovChain[200, test] // Timing

MarkovQ[test] // Timing

MarkovQ[test2] // Timing

MarkovQ[test3] // Timing

MarkovQ[test4] // Timing

MarkovErgodicProbabilities[test] // Timing

MarkovErgodicProbabilities[test2]

MarkovErgodicProbabilities[test3]

MarkovErgodicProbabilities[test4]