clear, close all

%% Markov Chains
% Consider a chain of "states" that describe a system that evolves in time.
% Let an element of such a chain be x(t_n), where {\bf x} \in R^m
%  be the state of the system at time t_n.  In Macdonald and Ridge (1988), the sociologists  consider  the mobility of working individuals in the UK. They focused on three types of working groups: executives, of which at some given time was given by x_E, skilled workers x_S, and unskilled workers x_U. Time ordering  was generation n. So, at t_n the population was {\bf x}(t_n) = [x_E,x_S,x_U](t_n), their offspring is {\bf x}(t_{n+1}).  They found that
% {\bf x}_{n+1} = T {\bf x}_n,  \quad
% n=0, 1, 2...            (1)
% where the (constant) transition matrix was
% T=\left( \begin{array}{ccc}
% 0.60& 0.29 & 0.16\\
% 0.26 & 0.37 & 0.27 \\
% 0.14& 0.34& 0.57
%  \end{array} \right).
% The transition matrix, representing the probabilities of transitioning,  is interpreted as follows: concerning the first row, there is 0.6 chance that a child of an executive will become an executive. There is a 0,29 chance that a child from a skilled parent becomes an executive,  and there's a 0.16 chance a child from an unskilled parent becomes an executive. We note that the column sum of T is 1:

T=[0.6, 0.29, 0.16;0.26,0.37,0.27;0.14,0.34,0.57]; % the transition matrix
for jj = 1:3
rsum(jj)=sum(T(:,jj));
fprintf('row %d sum = %f \n',jj,rsum(jj))
end

% If the  present worker population distribution at t_n=t_0 = 0, are {\bf x}_0 = [0.3, 0.5, 0.2]^T, we might ask, what the worker population might be after 20 "generations?" We then use (1) to compute:



%nend=input('enter n \n')
nend = 20;
gens=1:20;
x=zeros(3,nend);
x(:,1)=[0.3,0.5,0.2]'; % initial worker population
for jj = 2:nend
    x(:,jj) = T*x(:,jj-1);
    z(:,jj) = x(:,jj)/sum(x(:,jj));
end    
colormap('cool')
bar(gens,x','stacked')
ylabel('number of workers, normalized to workers at t=0')
xlabel('generation')
title('Figure 1')

% The population is compared, as a function of time to the initial population.  Plotted in Figure 1 is the evolution of the worker distribution, with time. We see no growth of the vector {\bf x}(t_n)  due to the structure of the transition matrix T. We can surmize that at least one of the  eigenvalues of T must be  1 and all others between 0 and 1. This is a property of Markovian transition matrices and it is built-in by having row sums equal to 1:

ev = eig(T);
fprintf('The eigenvalues of T: %f, %f, %f. \n', ev(1,1),ev(2,1),ev(3,1))

%
%
% What about the distribution of each worker group? how did this distribution change over time? We plot in Figure (2) the worker distribution, normalized to the total distribution for each generation:

bar(gens,z','grouped')
ylabel('(percentage) distribution of  workers')
xlabel('generation')
title('Figure 2')

% There is a quick readjustment of the relative percentages of workers but overall, the distribution remains fairly unchanging over long periods of time (or generations).
%
% Other things you might try:
% (1)  the command T=rand(3) will generate a random matrix with entries between 0 and 1. Replace this T in the above script and re-run. This T is no longer a row-sum =1 matrix, so the process is no longer a Markov chain (the entries in T are now interpreted as rates, rather than as probabilities).  Does x stay the same size, does it increase? does it decrease? check eig(T) to make some sense of the results.
% (2) we can construct a T with a specifc set of eigenvalues. Let D = diag(lambda1, lambda2, lambda3). The lambdai you specify, and be sure to choose all of the different. Choose these all real, but be sure to choose a set with at least 1 (slightly) larger than 1 in modulus,  both positive and negative. Try all within [0,1) in modulus but a mix of positive and negative ones. To build T: after choosing D, compute U=rand(3), and find U^{-1} (the command UI = U\eye(3)  gives the inverse of U called UI).
% So form T=U*D*UI. The U are going to contain the eigenvectors of T and D the eigenvalues of T. Rerun the scripts above with this T, with the understanding that these are no longer probabilities but instead are just rates (this is no longer a Markov chain generator).