% Optimal_Social_Insurance_and_Health_Inequality.m
% code opimizes social welfare over h1 h2 R and s
% needs subroutines OSIH_fun.m and OSIH_labor_supply_fun.m
% results for Table 1 and 2 are reported in outputvector out


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

% set desired output
fig=3 %(selects Figure to print)
Lauf=1 %(adjusts color of lines in figures)
status_quo=0;  %(set status_quo=1 for status quo policy scenario)
calbrat=0; %(set calbrat=1 to calibrate h2 via the h1/h2-share)

unit_length=45; %(period lenght is 45 years)
labor_guess=[1.01064067319871,0.909922361642239,0.818361010642285,0.735173612489716,0.659645735566390,0.591125220573042,0.529016387049861,0.472774709754942,0.421901931568669,0.375941586502710,0.334474912804406,0.297117141636886,0.263514150880929,0.233339475747535,0.206291667725995,0.182091990779131,0.160482438843887,0.141224052210225,0.124095503200396,0.108891914933515,0.0954238720089830;];
tau_guess=0.155;
R_guess=1;
s_guess=0.187;
h1_guess=2;
N1_guess=1;
% Find the solution via manual iteration over N2_guess (solution fulfils N2=N2_guess)
N2_guess=0.3227;

% Set Parameters:
n = 0:20; %(number of deficits=19 (20 states including no deficits)
sigma=1;
eta=0.58;
epsi=0.14;
z=0.5;
A=36;
rho=0.02;
xi=1;
gamma=0.25;
nu=0.163;

%retirement / policy
T_max=1.78; 
R=1;
tau_g=0.25;
tau=0.155;
s=0.187;
% tau is tau_h in the paper; s is tau_s in the paper; tau_g is tau_w in the paper)

% labor supply
kappa_0=0.5;
chi=0.062;

% health expenditure and aging
alpha_1=1.9;
alpha_1=1.5;
alpha_2=3.8;
beta_1=0.816;
beta_2=0.68;
%beta_2=beta_2*2
b=2.5;

alpha=1/3; r=rho;
omega=(1-alpha)*(alpha/r)^(alpha/(1-alpha));
w=omega.*A;
kappa_0=kappa_0*A^(1-sigma);
delta=0.25;

%**************** Gompertz-Makeham survival prob **********************
alpha_gom=0.11;
R_gom=0.00001;
if Lauf==2
alpha_gom=0.08;
R_gom=0.00018;
end
age=55:100;
S=exp(-R_gom./alpha_gom.*(exp(alpha_gom.*age)-1));
% *********************************************************************


%***** calibrate h1 via the h2-share  **********************+
if calbrat==1
h1=h1_guess;
s=0.187;
tau=0.155;
newguess=ones(1,21);
options=optimset('TolFun',10^(-20),'MaxFunEvals',10^(8),'MaxIter',10^(5),'TolX',10^(-20),'Display','off');
ell= fsolve('OSIH_laborsupply_fun', newguess, options);
a1=alpha_1*exp(-beta_1*h1/A);
g1 = my_poisspdf(n,a1);
h2= tau.*(1+0*s).*sum(g1.*w.*ell)-h1;
hsum=h1+h2;
h2share=h2/(h1+h2);
h1_h2=h1/h2;
L=sum(g1.*ell);
GDP=w.*L./(1-alpha);
health_share=hsum/GDP;
return
end
%************ end calbrate ************************************

if status_quo==0
options = optimset('TolFun',10^(-18),'Display','off');
newguess=[h1_guess;tau_guess;R_guess;s_guess];
lb = [ 0, 0.01, 0.01, 0.01 ];
ub = [ 1000, 0.9, T_max, 0.9 ];
[x,fval] = fmincon(@OSIH_fun,newguess, [], [], [], [], lb, ub );
h1=x(1); tau=x(2); R=x(3); RR=R.*45+20; s=x(4);
end

if status_quo==1
newguess=[h1_guess;R_guess;s_guess];
lb = [0, 0.01, 0.01 ];
ub = [1000 T_max, 0.9 ];
[x,fval] = fmincon(@OSIH_fun,newguess, [], [], [], [], lb, ub );
h1=x(1); R=x(2); RR=R.*45+20
s=0.187
end
newguess=x;

labor_guess=[1.01064067319871,0.909922361642239,0.818361010642285,0.735173612489716,0.659645735566390,0.591125220573042,0.529016387049861,0.472774709754942,0.421901931568669,0.375941586502710,0.334474912804406,0.297117141636886,0.263514150880929,0.233339475747535,0.206291667725995,0.182091990779131,0.160482438843887,0.141224052210225,0.124095503200396,0.108891914933515,0.0954238720089830;];
options=optimset('TolFun',10^(-20),'MaxFunEvals',10^(8),'MaxIter',10^(5),'TolX',10^(-20),'Display','off');
ell= fsolve('OSIH_laborsupply_fun', labor_guess, options);
labor_guess=ell;

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);
h2= tau.*(1+s).*sum(g1.*w.*ell)-h1;
a2=alpha_2*exp(-beta_2*h2/N2_guess/A);

hsum=h1+h2;
h2share=h2/(h1+h2)
GDP=w.*L./(1-alpha);
health_share=hsum/GDP
TT=tau_g*w*L;

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

% individual pension payments
E=0;Old=0;

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,:));
e(:,i)=T_R'.*g(:,i).*ell(i);
sum_e(i)=sum(e(:,i));
old(:,i)=T_R'.*g(:,i);
sum_old(i)=sum(old(:,i));
E=E+sum_e(i);
Old=Old+sum_old(i);
end
x=D/E.*ell;

%Anzahl Young and Old
N2=d*T_R';
N1=g1*R';

% 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); 

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

rest=1-sum(sum(g));
total_welf=sum(sum(welf));
totalpay=sum(sum(pay)); %should be equal to D for balanced budget

%distribution 2nd period
for i=1:1:length(n)
d(i)=sum(g(i,:));
welfare(i)=sum(welf(i,:));                                                                                                                                                                                                                                  
end

mean_n1=g1*n'; mean_n2=d*n';

%Anzahl Young and Old
N2=d*T_R';
N1=g1*R';

LE=20+unit_length.*N1+unit_length.*N2;
popshare=N1/(1+N2);
life_length=20+unit_length.*R+unit_length.*T_R;

unequal_1=GiniCoeff(g1',n');
unequal_2=GiniCoeff(d',n');

% insert donominator from benchmark run
welf_gain=(total_welf/5.0347-1)*100;
out=[round(h1/h2*100)/100 round(tau*1000)/10 round(s*1000)/10 round(RR*10)/10 round(LE*10)/10 round(mean_n1*100)/100  round(mean_n2*100)/100 round(unequal_1*100)/100 round(unequal_2*100)/100 welf_gain]'

if fig==1
%****** survival functions ***********
surv_prob=cumsum(d);
set(gcf,'PaperUnits', 'centimeters','PaperType','a4letter','Units','centimeters','PaperPosition', [00.5 1.5 16 7],'Position',[0.5 1.5 16 7],'Color','white');
set(0,'DefaultTextColor','k')

if Lauf==2
plot(age,S,'b-','LineWidth',1.5)
hold on
plot(life_length,surv_prob,'bo','LineWidth',1.5)
end
if Lauf==1
plot(age,S,'r--','LineWidth',1.5)
hold on
plot(life_length,surv_prob,'ro','LineWidth',1.5)
end
set(gca,'FontName','Times','FontSize',11);
axis([55 105 0 1])
ylabel('survival prob S(age)','Interpreter','LaTex','FontSize',11)
xlabel('age','Interpreter','LaTex','FontSize',11)
end


if fig==2
%*** plot frailty index ***************
set(gcf,'PaperUnits', 'centimeters','PaperType','a4letter','Units','centimeters','PaperPosition', [00.5 1.5 10 7],'Position',[0.5 1.5 10 7],'Color','white');
set(0,'DefaultTextColor','k')
f=n./(n(end)+1);
if Lauf==1
plot(f,g1,'b-','LineWidth',1.5)
end
if Lauf==2
plot(f,d,'r--','LineWidth',1.5)
end
hold on
axis([-0.03 0.8 -0.03 0.45])
xlabel('health deficit index','FontSize',11)
end


if fig==3
   
 %****** plots deficit distribution etc **********
set(gcf,'PaperUnits', 'centimeters','PaperType','a4letter','Units','centimeters','PaperPosition', [00.5 1.5 16 10],'Position',[0.5 1.5 16 10],'Color','white');
set(0,'DefaultTextColor','k')

subplot(2,1,1)
if Lauf==1
   plot(n,g1,'ro--','LineWidth',1.5)
end
if Lauf==2
   plot(n,g1,'bo-','LineWidth',1.5)
end
if Lauf==3
   plot(n,g1,'go-.','LineWidth',1.5,'Color',[0 0.5 0])
end
title('Working Age')
axis([0 10 -Inf 0.4])
set(gca,'FontName','Times','FontSize',11);
%set(gca,'YTick',[0 0.01 0.02 0.03 0.04])
%grid on
hold on


subplot(2,1,2)
if Lauf==1
   plot(n,d,'ro--','LineWidth',1.5)
end
if Lauf==2
   plot(n,d,'bo-','LineWidth',1.5)
end
if Lauf==3
  plot(n,d,'go-.','LineWidth',1.5,'Color',[0 0.5 0])
end
title('Retirement')
axis([0 10 0 0.16])
set(gca,'FontName','Times','FontSize',11);
xlabel('number of deficits','FontSize',11)
%set(gca,'YTick',[0 250 500 750 1000])
%grid on
hold on


end





return

print -depsc  c:\docs\paper\volker_health\volker_health_new_4.eps
eps2pdf('c:\docs\paper\volker_health\volker_health_new_4.eps','C:\Programme\gs\gs9.10\bin\gswin64c.exe')
