function f = OSIH_fun(newguess)
% function file for Optimal_Social_Insurance_and_Health_Inequality.m
% computes the (constrained) optimal policy

global a1 sigma nu n r T_max beta_1 beta_2 rho delta chi alpha gamma A alpha_1 alpha_2 b xi tau s h1 w eta kappa_0 labor_guess z tau_g N1_guess N2_guess status_quo

if status_quo==0
h1=newguess(1);
tau=newguess(2);
R=newguess(3);
s=newguess(4);
end

% set constraints for constrained optimization
if status_quo==1
h1=newguess(1);
tau=0.155;
R=newguess(2);
s=0.187;
s=newguess(3);
end

omega=(1-alpha)*(alpha/r)^(alpha/(1-alpha));
w=omega.*A;

options=optimset('TolFun',10^(-20),'MaxFunEvals',10^(8),'MaxIter',10^(5),'TolX',10^(-20));
ell= fsolve('OSIH_laborsupply_fun', labor_guess, options);

a1=alpha_1*exp(-beta_1*h1/N1_guess/A);
g1 = my_poisspdf(n,a1);

L=sum(g1.*ell);
T=T_max.*exp(-chi.*n);
T_R=max(T-R,eps);
R=min(R,T);
TT=tau_g*w*L;
h2= (tau.*(1+s).*sum(g1.*w.*ell)-h1);
a2=alpha_2*exp(-beta_2*h2/N2_guess/A);
hsum=h1+h2;

for i=1:1:length(n)
lambda=a2+b*n(i);
g2(:,i) = my_poisspdf(n,lambda);
g(:,i)=g1(i)*g2(:,i);
d(i)=sum(g(i,:));
end

N1=g1*R';
N2=d*T_R';

% pension budget
D=sum((1+r).*R.*s.*g1.*w.*ell);

% individual pension payments
for i=1:1:length(n)
    if i==1
        E=0;
    end
    
lambda=a2+b*n(i);   
g2(:,i) = my_poisspdf(n,lambda);
g(:,i)=g1(i)*g2(:,i);
e(:,i)=T_R'.*g(:,i).*ell(i);
sum_e(i)=sum(e(:,i));
E=E+sum_e(i);
end

if E==0
    x=eps.*ell;
else
x=D/E.*ell;
end

% consumption / utility period 1
c1=(1-tau_g-tau-s).*w.*ell +z*w*ell(3)+TT/N1;
if sigma==1
u1=log(c1)-kappa_0.*exp(delta.*n).*ell.^(1+1./eta)./(1+1./eta);
else
u1=(c1.^(1-sigma)-1)./(1-sigma)-kappa_0.*exp(delta.*n).*(ell.^(1+1./eta)-1)./(1+1./eta);
end

% utility period 2 and welfare
c2=zeros(1,length(n));

for i=1:1:length(n)
lambda=a2+b*n(i);
g2(:,i) = my_poisspdf(n,lambda);
g(:,i)=g1(i)*g2(:,i);
welf1(:,i)=(1-exp(-rho*R(i)))/rho.*u1(i).*g(:,i);
c2(i,:)=(1-tau)*x(i); %+z*w*ell(3);
if sigma==1
u2(i,:)=log(c2(i,:));
else
u2(i,:)=((c2(i,:)).^(1-sigma)-1)./(1-sigma);
end

pay(i,:)=(x(i)*T_R)'.*g(:,i);
fack=max(0, ( exp(-rho.*R(i)*ones(1,length(n)) )-exp(-rho.*T) )'./rho);
welf2(:,i)=( fack.*u2(i,:)'-nu.*(1+xi.*n(i)).*R(i).^(1+1/gamma) ).*g(:,i);
welf(:,i)=welf1(:,i)+welf2(:,i);

end

total_welf=sum(sum(welf));
f=-total_welf;
