%
% sal2.m is an example matlab program 
% to simulate a simple Jovanovic model
%

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

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

% ----- policy solutions ----- %

w0bar = 4;
w1bar = 5;

% ----- periods and people ----- %

t = 40
n = 500

% ------ draw an initial wage for all agents ----- % 

wage=discreteinvrnd(mu,n,1)-1 ;

% ----- optimal policy conditional on wage ----- %

for j = 1:n 
 if wage(j,1)   >= w0bar
  policy(j,1) = 1 ;
  wageob(j,1) = wage(j,1);
 else 
  policy(j,1) = 2;
  wageob(j,1) = 0;
 end     
end

% ----- draw a wage conditional on previous policy ----- %

for j = 1:n 
 if policy(j,1) == 1 
  temp = markov(P,3,wage(j,1)+1,0:9) ;   
  wage(j,1+1) = temp(2);
 else 
  wage(j,1+1) = discreteinvrnd(mu,1,1)-1 ;
 end     
end

% ----- optimal policy conditional on wage and previous policy ----- %

for j = 1:n 
 if policy(j,1) == 1 & wage(j,1+1) >= w1bar
  policy(j,1+1) = 1 ;
  wageob(j,1+1) = wage(j,1+1);
 elseif policy(j,1) == 2 & wage(j,1+1) >= w0bar 
  policy(j,1+1) = 1 ;
  wageob(j,1+1) = wage(j,1+1);
 else 
  policy(j,1+1) = 2;
  wageob(j,1+1) = 0; 
 end
end

% ----- now automate the above algorithm ----- %

for i = 1:(t-2)

% - draw a wage conditional on previous policy - %

for j = 1:n 
 if policy(j,1+i) == 1 & policy(j,i) == 1
  wage(j,1+1+i) = wage(j,1+i);   
 elseif policy(j,1+i) == 1 & policy(j,i) == 2
  temp=markov(P,3,wage(j,1+i)+1,0:9) ;   
  wage(j,1+1+i) = temp(2);  
 else 
  wage(j,1+1+i) = discreteinvrnd(mu,1,1)-1 ;
 end     
end   
    
% - optimal policy conditional on wage and previous policy - %

for j = 1:n 
 if policy(j,1+i) == 1 & wage(j,1+1+i) >= w1bar
  policy(j,1+1+i) = 1 ;
  wageob(j,1+1+i) = wage(j,1+1+i);
 elseif policy(j,1+i) == 2 & wage(j,1+1+i) >= w0bar 
   policy(j,1+1+i) = 1 ;
   wageob(j,1+1+i) = wage(j,1+1+i);
 else 
  policy(j,1+1+i) = 2;
  wageob(j,1+1+i) = 0; 
 end
end    

end    

% ----- compute labor force participation rate ----- %

lfp = mean((1./policy-1/2)*2) ;

% ----- print a typical simulated agents wage and observed wage ----- %

[wage(14,1:t)' wageob(14,1:t)' policy(14,1:t)']

comptime = etime(clock, time0)
diary off