%
% sal3.m is an example matlab program 
% to compute choice probabilities
%

% ----- save output to file ----- %

clear all
diary sal3.mout
diary off
delete sal3.mout
diary sal3.mout

% ----- display date and time of computation ----- %

%format short
date
time0 = clock;

% ----- set control vars ----- %

iters  = 200;
states = 10;

% ------ set parameters ------ %

beta = 0.95;
c    = 2*ones(states,1);
P    = [0.80 0.20 0 0 0 0 0 0 0 0;
        0.08 0.50 0.42 0 0 0 0 0 0 0;
        0.08 0.21 0.50 0.21 0 0 0 0 0 0;
        0.08 0 0.21 0.50 0.21 0 0 0 0 0;
        0.08 0 0 0.21 0.50 0.21 0 0 0 0;
        0.08 0 0 0 0.21 0.50 0.21 0 0 0;
        0.08 0 0 0 0 0.21 0.50 0.21 0 0;
        0.08 0 0 0 0 0 0.21 0.50 0.21 0;
        0.08 0 0 0 0 0 0 0.21 0.50 0.21;
        0.08 0 0 0 0 0 0 0 0.42 0.50 ] ;

% ----- find limiting distribution ----- %

[v,d]=eig(P') ;

% - find the row of the unit eigenvalue - %

[d1,d2]=find(fix(diag(d)+0.000001))

% - normalize eigenvector - %

mu = abs(v(:,d1))./sum(abs(v(:,d1))) 

% ----- solutions ----- %

w0bar = 4
w1bar = 5

% ----- compute chice probabilities ------ %

for j = (w0bar+1):states ;
 delta0(j,1)=mu(j,1)/sum(mu(5:10,1)) ;
end 

gamma1 = P*delta0 ;

for j = (w1bar+1):states ;
 delta1(j,1)=gamma1(j,1)/sum(gamma1(6:10,1)) ;
end 

gamma2 = eye(10)*delta1;

for j = (w1bar+1):states ;
 delta2(j,1)=gamma2(j,1)/sum(gamma2(6:10,1)) ;
end 

[delta0 delta1 delta2]

comptime = etime(clock, time0)
diary off