% MCHAIN file mchain.m  Simulation of Markov chains
% Version of 1/2/96  Revised 7/31/97 for version 4.2 and 5.1
% Assumes the procedure chainset has been run
n  = input('Enter the number n of stages   ');
st = input('Enter the initial state  ');
if ~isempty(st)
  s  = MS(st==states);
else
  s = 1;
end
T  = zeros(1,n);           % Trajectory in state numbers
U  = rand(1,n);
for i = 1:n
  T(i) = s;
  s = ((A(s,:) < U(i))&(U(i) <= B(s,:)))*MS';
end
N  = 0:n-1;
tr = [N;states(T)]';
n10 = min(n,11);
TR = tr(1:n10,:);
f  = ones(1,n)/n;
[sn,p] = csort(T,f);
if isempty(PI)
  disp('     State     Frac')
  disp([states; p]')
else
  disp('     State     Frac       PI')
  disp([states; p; PI]')
end
disp('To view the first part of the trajectory of states, call for TR')