
function r=imro3tau(tau)

% function is to obtain tau
% output is gov't deficit

global sigma rho hhat beta pee peu pi0 pi1 theta maxast incast nsim

%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-4 iter<2000]
    % 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;
end

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;

%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-4 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

% obtain government balance
sf=sum(f);
r=(sf(2)+sf(4))*(1-tau)*theta-sf(1)*tau;