(**)

mulMat[a_ , b_]:=Module[{all,aRows , aCols , bRows , bCols , cRows , cCols,currentC , currentA , currentB},
aRows = Dimensions[a][[1]];
aCols = Dimensions[a][[2]];
bRows = Dimensions[b][[1]];
bCols = Dimensions[b][[2]];
If[aCols≠ bRows , Throw["The number of columns in matrix a is not the same as the number of columns in matrix b"]];
cRows = aRows;
cCols = bCols;
currentC = Table["□" , {cRows} , {cCols}];
all = {Row[{MatrixForm[a] ,".", MatrixForm[b] , "=" , MatrixForm[currentC]}]};
Do[
Do[
Do[
currentA = 
Table[
If[
i==row&&j == k , 
Style[a[[i , j]] , {Bold,Red}] , 
Style[a[[i,j]] , Blue]] 
, {i , 1 , aRows} , {j , 1 , aCols}];
currentB = 
Table[
If[
i == k &&j==col , 
Style[b[[i , j]] , {Bold,Red}] , 
Style[b[[i,j]] , Blue]] 
, {i , 1 , bRows} , {j , 1 , bCols}];
currentC[[row , col]] = Sum[a[[row , i]] b[[i , col]] , {i , 1 , k}];
AppendTo[all , Row[{MatrixForm[currentA] ,".", MatrixForm[currentB] , "=" , MatrixForm[currentC]}]];
,{k , 1 , aCols}]
,{row , 1 , cRows}]
,{col , 1 , cCols}];
{all , a.b}
]

steps0 =mulMat[{{1 , 2 , 3}} , {{9 , 8 , 7} , {6 , 5 , 4} , {3 , 2 ,1}}];

Manipulate[steps0[[1]][[ii]] , {ii ,1 , Dimensions[steps0[[1]]][[1]] , 1}]

steps1 =mulMat[{{1 , 2 , 3},{4 , 5 , 6},{7 , 8 , 9}} , {{9 , 8 , 7} , {6 , 5 , 4} , {3 , 2 ,1}}];

Manipulate[steps1[[1]][[ii]] , {ii ,1 , Dimensions[steps1[[1]]][[1]] , 1}]

p::usage = "p[i] is the probability of finding the system in state i.";

state = {{p[1] , p[2] , p[3]}};

transition = Table[p[i->j] , {i , 1 , 3} , {j , 1 , 3}];

steps2 = mulMat[state , transition];

Manipulate[steps2[[1]][[ii]] , {ii ,1 , Dimensions[steps2[[1]]][[1]] , 1}]

steps3 = mulMat[transition , transition];

Manipulate[steps3[[1]][[ii]] , {ii ,1 , Dimensions[steps3[[1]]][[1]] , 1}]

(**)

tr::usage = "tr[a , b , p] represents a transition form state a to state b with probability p.";

transition = tr[_ , _ , _Real];

start[tr[a_ , b_ , p_Real]]:=a;
end[tr[a_ , b_ , p_Real]]:=b;
pr[tr[a_ , b_ , p_Real]]:= p;

getNodes[process:{transition...}]:=Union[Join[start/@process , end/@process]];

checkProcess[process:{transition...}]:=
Module[{nodes,checkEscape},
nodes = getNodes[process];
Do[
checkEscape = Total[pr/@Cases[process , tr[nodes[[ii]] , _ , _]]];
If[checkEscape ≠ 1, Throw["Total probability of escaping from node -"<>ToString[nodes[[ii]]]<>"- is not equal to 1.0"]]
,{ii , 1 , Dimensions[nodes][[1]]}];
]

makeGraph[process:{transition...}]:=Module[{checkEscape},
checkProcess[process];
Graph[Labeled[start[#] \[DirectedEdge]end[#] , pr[#]]&/@process,VertexLabels->"Name"]
];

getMatrix[process:{transition...} , OptionsPattern[{"rightMatrix"->False}]]:=
Module[{nodes , pm},
checkProcess[process];
nodes = getNodes[process];
pm = Table[
Total[pr/@Cases[process , tr[nodes[[ii]] , nodes[[jj]] , _Real]]]
,{ii , 1 , Dimensions[nodes][[1]]} , {jj , 1 , Dimensions[nodes][[1]]}];
If[OptionValue["rightMatrix"] , Transpose[pm],pm]
]

weather = {tr["Asunny" , "Brainy" , 0.1] , tr["Asunny" , "Asunny" , 0.9] , tr["Brainy" , "Asunny" , 0.5] , tr["Brainy" , "Brainy" , 0.5]};

makeGraph[weather]

weatherMatrix = getMatrix[weather];

weatherMatrix//MatrixForm

weatherEigen = Eigensystem[Transpose[weatherMatrix]];

weatherEigen[[2,1]]/Total[weatherEigen[[2,1]]]

%.weatherMatrix

weatherEigen[[2,2]]

stockMarket = {tr["bull" , "bull" , 0.9] , tr["bear" , "bear" , 0.8] , tr["stagnant" , "stagnant" , 0.5] , tr["bull" , "bear" , 0.075] , tr["bull" , "stagnant" , 0.025] , tr["bear" , "bull" , 0.15] , tr["bear" , "stagnant" , 0.05] , tr["stagnant" , "bull" , 0.25],tr["stagnant" , "bear" , 0.25]};

makeGraph[stockMarket]

stockMarket//getNodes

stockMarketMatrix = getMatrix[stockMarket];

stockMarketMatrix//MatrixForm

stockMarketEigen = Eigensystem[stockMarketMatrix//Transpose];

stockMarketEigen[[2,1]]/Total[stockMarketEigen[[2,1]]]

%.stockMarketMatrix

getNodes[stockMarket]

evolveState[n_ , initial_List][process:{transition...}]:=Module[{t , s , evolution},
t = getMatrix[process];
If[Dimensions[initial][[1]] ≠ Dimensions[t][[1]], Throw["Initial state size does not match probability matrix size."]];
s = initial / Total[initial];
evolution = NestList[#.t& , s , n];
Graphics[
{
Table[{GrayLevel[1.0-evolution[[ii , jj]]] , Disk[{jj , -ii} , 0.2]} , {ii , 1 ,Dimensions[evolution][[1]]} , {jj , 1 , Dimensions[evolution][[2]]}],
Table[Table[{GrayLevel[1.0-t[[jj , kk]]] , Arrow[{{jj , -ii} , {kk , -(ii+1)}}]} , {jj , 1 , Dimensions[evolution][[2]]} , {kk , 1 , Dimensions[evolution][[2]]}] , {ii , 1 ,Dimensions[evolution][[1]]-1}],
{Black,Text[#[[2]],{#[[1]] , -1+1}]&/@Partition[Riffle[Range[Dimensions[s][[1]]] , s],2]},
{Black,Text[#[[2]],{#[[1]] , -n-1-1}]&/@Partition[Riffle[Range[Dimensions[s][[1]]] , evolution[[-1]]],2]}
}
]
];

drawTrajectories[traj_]:=Module[{},
Table[Graphics[Line[Table[{traj[[ii]][[jj]],-jj} , {jj , 1 , Dimensions[traj[[ii]]][[1]]}]]] , {ii , 1 , Dimensions[traj][[1]]}]
];

generateTrajectories[n_ ,nt_, initial_List][process:{transition...}]:=Module[{t , s , evolution , all , repr , tab},
t = getMatrix[process];
If[Dimensions[initial][[1]] ≠ Dimensions[t][[1]], Throw["Initial state size does not match probability matrix size."]];
s = initial / Total[initial];
evolution = NestList[#.t& , s , n];
all = Table[RandomChoice[#-> Range[Dimensions[s][[1]]]]&/@evolution , {nt}];
tab=Table[Table[N[Count[all[[1;;tt , tm]] ,st]/tt] , {tm , 1 , Dimensions[evolution][[1]]} , {st , 1 , Dimensions[evolution][[2]]}] , {tt , 1 , nt}];
Table[ 
Graphics[
{
Table[{GrayLevel[1.0-tab[[tt]][[ii , jj]]] , Disk[{jj , -ii} , 0.2]} , {ii , 1 ,Dimensions[tab[[tt]]][[1]]} , {jj , 1 , Dimensions[tab[[tt]]][[2]]}],
{Black,Text[#[[2]],{#[[1]] , -1+1}]&/@Partition[Riffle[Range[Dimensions[s][[1]]] , s],2]},
{Black,Text[#[[2]],{#[[1]] , -n-1-1}]&/@Partition[Riffle[Range[Dimensions[s][[1]]] , tab[[tt]][[-1]]],2]},
{Red , Line[Table[{all[[tt]][[jj]],-jj} , {jj , 1 , Dimensions[all[[tt]]][[1]]}]]}
}
]
, {tt , 1 , nt}]
];

{0.3124999999999998`,0.6250000000000001`,0.06250000000000007`}

getNodes[stockMarket]

Manipulate[evolveState[8 , {pBear , pBull , pStagnant}][stockMarket] , {{pBear , 1} , 0 , 1} , {{pBull , 0} , 0 , 1} , {{pStagnant , 0} , 0 , 1}]

trajs = generateTrajectories[8 , 1000,{1 , 0 , 0}][stockMarket];

Manipulate[trajs[[ii]] , {ii , 1 , Dimensions[trajs][[1]] , 1}]