function [chain,state]=markov(T,n,s0,V);
%function [chain,state]=markov(T,n,s0,V);
%  chain generates a simulation from a Markov chain of dimension
%  the size of T
%
%  T is transition matrix
%  n is number of periods to simulate
%  s0 is initial state
%  V is the quantity corresponding to each state
%  state is a matrix recording the number of the realized state at time t
%
%r is rows of the transition matrix; c is columns of the transition matrix
[r c]=size(T);
if nargin == 1;%If we only provide the transition matrix as an input of the fct markov, then the program assigns the quantities (1,...,r) to the
  V=[1:r];     %states, the initial state is 1 and it simulates 100 periods
  s0=1;
  n=100;
end;
if nargin == 2;%If we only provide the transition matrix and number of periods, then the program assigns the quantities (1,...,r) 
  V=[1:r];     %to the states, and the initial state is 1
  s0=1;
end;
if nargin == 3;%If we provide everything except the quantities corresponding to each state, then the program assigns the quantities
  V=[1:r];     %(1,...,r)
end;
%If 
if r ~= c;%Error message if the transition matrix is not square
  disp('error using markov function');
  disp('transition matrix must be square');
  return;
end;
%
for k=1:r;%Error message if the matrix is not stochastic(rows do not add up to 1). The progran automatically normalizes to 1
  if sum(T(k,:)) ~= 1;
    disp('error using markov function')
    disp(['row ',num2str(k),' does not sum to one']);
    disp(['normalizing row ',num2str(k),'']);
    T(k,:)=T(k,:)/sum(T(k,:));
  end;
end;
[v1 v2]=size(V);
if v1 ~= 1 |v2 ~=r %Erroe message if the matrix of values does not conform to the transition matrix
  disp('error using markov function');
  disp(['state value vector V must be 1 x ',num2str(r),''])
  if v2 == 1 &v2 == r; %If the values are given in row format, the program will transpose the matrix
    disp('transposing state valuation vector');
    V=V';
  else;
    return;
  end;  
end
if s0 < 1 |s0 > r;%Error message if the initial state is out of the range of the values that the variable can take
  disp(['initial state ',num2str(s0),' is out of range']);
  disp(['initial state defaulting to 1']);
  s0=1;
end;
%
%T
%rand('uniform');
X=rand(n-1,1); %Generates of column of n-1 values drawn from a uniform distrib
s=zeros(r,1);
s(s0)=1;
cum=T*triu(ones(size(T))); %Generates the upper triangular part of a matrix of ones
%This generates the matrix of values for the n simulated periods
for k=1:length(X);
  state(:,k)=s;
  ppi=[0 s'*cum];
  s=((X(k)<=ppi(2:r+1)).*(X(k)>ppi(1:r)))';
end;
chain=V*state;