clear all
clc

trans = [0.5 0.5 0.0 0.0;
        0.1 0.9 0.0 0.0;
        0.0 0.0 0.3 0.7;
        0.0 0.0 0.2 0.8];
    
 n=100;
 state = zeros(4,n);
 path=zeros(n);
 
 for i=1:4;
     rand1 = rand/2;
     rand2 = rand/2;
     state(:,1)=[rand1 0.5-rand1 rand2 0.5-rand2];
     for j=1:n;
        x = rand;
        y = find(x<cumsum(state(:,j)));
        path(j,i) = min(y);
        state(:,j+1)=trans(path(j,i),:)';
         for k=1:4;
             count(j,i,k)=length(find(path(1:j,i)==k))/j;
         end;
     end;
 end;

 for i=1:4;
     figure(i);
     stairs(path(:,i));
     xlabel('time')
     ylabel('X(i) state')
     title(['Sample Path #' int2str(i)]);
     saveas(figure(i),['Prob1 Sample path ' int2str(i) '.png'])
 end;

 x1 = [1:1:100];
 for i=1:4;
     figure(i+4);
     hold on
     plot(x1,count(:,i,1));
     plot(x1,count(:,i,2));
     plot(x1,count(:,i,3));
     plot(x1,count(:,i,4));
     xlabel('i');
     ylabel('Relative frequency');
     title(['sample path relative freq' int2str(i)]);
     saveas(figure(i+4),['Prob1 sample path rel freq ' int2str(i) '.png'])
     hold off
 end;
 
 for i=1:4;
 [v,d] = eig(trans');
 v1_eigen = v(1:4,1)/sum(v(1:4,1))
 v2_Rel_Freq = [count(n,i,1), count(n,i,2), count(n,i,3), count(n,i,4)]
 end;
 
 vectoreig1=v(1:4,2)/sum(v(1:4,2))
 vectoreig2=v(1:4,4)/sum(v(1:4,4))
 
 