%
% sal1.m is an example matlab program 
% to solve a simple Jovanovic model
%

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

clear all
diary sal1.mout
diary off
delete sal1.mout
diary sal1.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))) 

% ----- discretize wage ----- %

w = zeros(states,1);

for j = 1:(states-1);
 w(j+1,1) = w(j,1) + 1;
end;

% ----- value function iteration ----- %

% - initial guess of value function - %

v0  = zeros(states,1);

% - start iterations - %

crit1 = 3;

for j = 1:iters ;
 if crit1 < 0.0001, break, end;

% - compute second period payoffs and policy - %

  payoff11 = 1/(1-beta)*w ;
  payoff12 = c + beta*mu'*v0 ;
  payoff1  = [payoff11  payoff12 ];
  
  [v1,policy1] = max(payoff1') ;
  
% - put matrices in correct dimension - %
  
  v1      = v1';
  policy1 = policy1';
  
% - compute the first period payoffs and policy - %
  
  payoff01 = w + beta*P*v1;
  payoff02 = c + beta*mu'*v0 ;
  payoff0  = [payoff01  payoff02 ];
  
  [newv0,policy0] = max(payoff0') ;
  
% - put matrices back to correct dimension - %
 
  newv0   = newv0';
  policy0 = policy0';
  
 iter = j;

% - compute criterion - %

 crit1 = sqrt((newv0-v0)'*(newv0-v0))/(1+sqrt(v0'*v0)) ;

% - update - %

 v0 = newv0;

end;

% ----- print solutions ----- %

[w,v0,policy0,v1,policy1],  [crit1, iter] 

w0bar = sum(any(policy0-1,2))
w1bar = sum(any(policy1-1,2))

comptime = etime(clock, time0)
diary off