%This program hopefully solves the model in Hansen & Imrohoroglu (JPE 1992)

clear all;
time=cputime;
global sigma rho hhat beta pee peu pi0 pi1 theta maxast incast nsim

%gtheta=0;gtau=0;gdeficit=0;gemploy=0;gutil=0;

sigma=0.67;   %Blend consumption-leisure
rho=2.5;      %Risk aversion parameter.
hhat=0.45;    %Working hours
beta=.995;    %Discount factor.
pee=.9681;    %Probability of receiving offer when offered last period.
peu=.5;       %Probability of receiving offer when not offered last period.
pi0=.1;        %Prob. not getting caught if not worked last period
pi1=0;        %Prob. not getting caught if worked last period (quitter)
maxast=4;   %Maximum number of assets.
incast=.1;  %Minimum number of assets.
nsim=50000;   %Number of simulations at end.

for theta=0.35:0.05:0.35;
disp('theta');disp(theta);
% Call routine that finds tau
%opt=optimset('TolX',0.01);
tau=fzero('imro3tau',[.0 .5]);

%Tabulate return function, i.e.,
% create matrices (a')x(a) with possible utility payoffs
% (Note that -10000 stands for minus infinity)
nasts=round(maxast/incast+1);
util1=-10000*ones(nasts,nasts);  %Utility when employed.
util2=-10000*ones(nasts,nasts);  %Utility when getting benefits
util3=-10000*ones(nasts,nasts);  %Utility when getting nothing

for a=0:incast:maxast
   i=round(a/incast+1);
   apr=0;
   while [apr<a+1-tau-.0001 apr<maxast+.0001]
      j=round(apr/incast+1);
      util1(j,i)=((((1-tau)+a-apr)^(1-sigma)*(1-hhat)^(sigma))^(1-rho)-1)/(1-rho);
      apr=apr+incast;
   end
   apr=0;
   while [apr<a+(1-tau)*theta-.0001 apr<maxast+.0001]
      j=round(apr/incast+1);
      util2(j,i)=(((1-tau)*theta+a-apr)^((1-sigma)*(1-rho))-1)/(1-rho);
      apr=apr+incast;
   end
   apr=0;
   while [apr<a-.0001 apr<maxast+.0001]
      j=round(apr/incast+1);
      util3(j,i)=((a-apr)^((1-sigma)*(1-rho))-1)/(1-rho);
      apr=apr+incast;
   end
end 

%Iterate on Belman's equation.
v=zeros(3,nasts);
% the first dimension is combinations of s'={e,u} e eta'={0,1}, that is
% 1=e0, 2=e1, 3=u
test=10;
iter=0; 
while [max(max(test))>1E-6 iter<4000]
    % assets decisions, given eta', miu, and s
    re=util1+beta*(pee*v(2,:)+(1-pee)*v(3,:))'*ones(1,nasts); % employed
    ru=util2+beta*(peu*v(1,:)+(1-peu)*v(3,:))'*ones(1,nasts); % unemployed, didn't get offer (with benefits)
    ru0=util3+beta*(pee*v(1,:)+(1-pee)*v(3,:))'*ones(1,nasts); % unemployed, got offer, without benefits (miu=0)
    ru1=util2+beta*(pee*v(1,:)+(1-pee)*v(3,:))'*ones(1,nasts); % unemployed, got offer, with benefits (miu=1)
    [tve,aae]=max(re);
    [tvu,aau]=max(ru);
    [tvu0,aau0]=max(ru0);
    [tvu1,aau1]=max(ru1);
    % eta' decision (only if got offer), given eta
    r0(1,:)=pi0*tvu1+(1-pi0)*tvu0;   %decided not to work, hadn't worked before
    r0(2,:)=tve;                     %decided to work, hadn't worked before
    r1(1,:)=pi1*tvu1+(1-pi1)*tvu0;   %decided not to work, worked before
    r1(2,:)=tve;                      %decided to work, worked before
    [tv0,mm0]=max(r0);
    [tv1,mm1]=max(r1);
    tv=[tv0;tv1;tvu];                %note order is important
    test=abs(tv-v);
    v=tv;
    iter=iter+1;
    % disp(max(max(test)));
end
disp('Dynamic Program took');;disp(iter);;disp('iterations');

decism=[mm0-ones(1,nasts);mm1-ones(1,nasts)]; % choice to work (1) or not (0) , contingent on worked last period eta={0,1}
decisa=[aae;aau;aau0;aau1]; % assets next period contingent on {employed(working), no offer, without benefits, with benefits}
decisar=(decisa-1)*incast;

x=0:incast:maxast;
figure(1);plot(x,decisar);
%figure(2);
%subplot(2,2,1);plot(x,decisar(1,:));title('employed');
%subplot(2,2,2);plot(x,decisar(2,:));title('no offer');
%subplot(2,2,3);plot(x,decisar(3,:));title('w/o benefits');
%subplot(2,2,4);plot(x,decisar(4,:));title('w/ benefits');

%Compute invariant distribution.
da1=zeros(nasts,nasts); da2=da1;da3=da1;da4=da1;
for i=1:nasts
   da1(decisa(1,i),i)=1;
   da2(decisa(2,i),i)=1;
   da3(decisa(3,i),i)=1;
   da4(decisa(4,i),i)=1;
end
dm0w=diag(decism(1,:)); %decides to work if was NOT working before
dm0n=diag(ones(1,nasts)-decism(1,:)); %decides not to work when NOT before
dm1w=diag(decism(2,:)); %decides to work if WAS working before
dm1n=diag(ones(1,nasts)-decism(2,:)); % not to work when WAs before

f=ones(nasts,4)/(4*nasts); %Initial guess.
diter=0;
dtest=10;
while [max(max(dtest))>1E-6 diter<500]
    fpr=[(dm1w*pee*da1*f(:,1)+dm0w*(peu*da2*f(:,2)+pee*da3*f(:,3)+pee*da4*f(:,4))) ...                %working
        ((1-pee)*da1*f(:,1)+(1-peu)*da2*f(:,2)+(1-pee)*da3*f(:,3)+(1-pee)*da4*f(:,4)) ...             %no offer
        ((1-pi1)*dm1n*pee*da1*f(:,1)+(1-pi0)*dm0n*(peu*da2*f(:,2)+pee*da3*f(:,3)+pee*da4*f(:,4))) ...  %with offer NO benefits
           (pi1*dm1n*pee*da1*f(:,1)    +pi0*dm0n*(peu*da2*f(:,2)+pee*da3*f(:,3)+pee*da4*f(:,4)))];     %with offer with benefits
    dtest=abs(fpr-f);
    f=fpr;
    diter=diter+1;
end
%x=0:incast:maxast;
%figure(3);plot(x,f);
%figure(4);
%subplot(2,2,1);plot(x,f(:,1));title('employed');
%subplot(2,2,2);plot(x,f(:,2));title('no offer');
%subplot(2,2,3);plot(x,f(:,3));title('w/o benefits');
%subplot(2,2,4);plot(x,f(:,4));title('w/ benefits');


%Compute average asset holdings using invariant distribution.
fa=f(:,1)+f(:,2)+f(:,3)+f(:,4);
assets=0:incast:incast*(nasts-1);
massets=sum(assets'.*fa); 
disp('Mean Asset Holdings Using Invariant Distribution');disp(massets);

%Compute average asset holding by simulating decision rules
%     Simulate nsim periods, and throw away the first 100.
a=massets;s=1;eta=0;  % initial states
rand('state',123456); % resets uniform distrib. drawer
state=rand(nsim,2); % employment and moral hazard
suma=zeros(nsim,1);
for i=1:nsim
   ia=round(a/incast+1);
   if s==1;                % offer last period
       probs=pee;
   else;                   % NO offer last period
       probs=peu;
   end;
   %  remember: decisar is assets', contingent on {employed(working), no offer, without benefits, with benefits}
   %  remember: decism is choice to work (1) or not (0) , contingent on worked last period eta={0,1}
   if state(i,1)<probs; % get offer this period
       s=1; 
       if eta==0;  %haven't worked last period
           if decism(1,ia)==1; %decides to work
               eta=1;              
               apr=decisar(1,ia);
           else;               %decides not to work
               eta=0;
               if state(i,2)<pi0; % got benefits
                   apr=decisar(4,ia);
               else;              % didn't get benefits
                   apr=decisar(3,ia);
               end;
           end;
       end;
       if eta==1;  %have worked last period
           if decism(2,ia)==1; %decides to work
               eta=1;              
               apr=decisar(1,ia);
           else;               %decides not to work
               eta=0;
               if state(i,2)<pi1; % got benefits
                   apr=decisar(4,ia);
               else;              % didn't get benefits
                   apr=decisar(3,ia);
               end;
           end;
       end;             
   else;                % no offer this period
       s=2; 
       eta=0; 
       apr=decisar(2,ia);
   end;
   suma(i,1)=apr;
   a=apr;
end
masset=mean(suma(101:nsim));
disp('Mean Asset Holding from Simulation');;disp(masset);

% obtain government balance
% remember f (working, not offer, no benefits, benefits)
sf=sum(f);
deficit=(sf(2)+sf(4))*(1-tau)*theta-sf(1)*tau;

% obtain average utility
uuu=(1-beta)*(tve*f(:,1)+tvu*f(:,2)+tvu0*f(:,3)+tvu1*f(:,4));

time=cputime-time;
disp('time elapsed:');disp(time);disp('seconds');

%gtheta=[gtheta theta];
%gtau=[gtau tau];
%gdeficit=[gdeficit deficit];
%gemploy=[gemploy sf(1)];
%gutil=[gutil uuu];

end;

%ngtheta=size(gtheta);
%gtheta=gtheta(2:ngtheta(2));
%gtau=gtau(2:ngtheta(2));
%gdeficit=gdeficit(2:ngtheta(2));
%gemploy=gemploy(2:ngtheta(2));
%gutil=gutil(2:ngtheta(2));
%figure(5);plot(gtheta,gutil);

