% Code for Figure 1 of Brookings Paper
% Flavien MAY 2017

% ancillary files
% data: MR_US1940F.csv
% code: logMRadd_Igrowth.m , survivalcurve.m , life_exp4.m, trunc_data.m

%% Estimates for US 1940 WOMEN with slope for investment 
% data file name is MR_US1940F.csv 

ns=1;                                 % number of simulations
y=1940                                % cohort

filename=strcat(strcat('MR_US',num2str(y),'F.csv')) 
target = csvread(filename,1,0);
target = log(target/100);             % ln per 100 rate       
Tdata =length(target);

% Dataseries goes until age 120 
% but after age 65 (i.e. year 2005) these are just extrapolation
% we truncate to use only the real data

Tdata=min(Tdata,65);
S0=trunc_data(survivalcurve(target),Tdata);
LE0=life_exp4(target,Tdata);

% Parameters for Optimization
nvars=6;                                % number of parameters to optimize upon
options_fmin=optimset('Display','iter','PlotFcns',@optimplotfval,'FunValCheck','on','MaxIter',30*nvars);

weight=ones(Tdata,1); weight(1)=10;     % put extra weight on Infant Rate M

% Initial Guess    
fmin_init=[0.0551    0.0012    0.1489    1.3045    1.7540    1.0224];     % GREAT

disp('Initial Guess')
disp(fmin_init)

% Optimization
criterion_s = @(Y) weight'*(abs(S0-trunc_data(survivalcurve(logMRadd_Igrowth([Y(1:5),1,0,Y(6)])),Tdata)));
[res_US_W1940_fmin_s(:,ns),fval_US_W1940_fmin_s(:,ns)]=fminsearch(criterion_s,fmin_init,options_fmin);

Y=res_US_W1940_fmin_s(:,ns)';

% Simulated LogMR and SurvivalCurve
LMR1=logMRadd_Igrowth([Y(1:5),1,0,Y(6)]);
S1=survivalcurve(LMR1);

LE_US_W1940_fmin_s(1,ns)=LE0;
LE_US_W1940_fmin_s(2,ns)=life_exp4(LMR1,65);                    

disp('fminsearch estimation with growth rate, all in one (7 parameters), survival curve')
disp ('Estimated Parameters')
disp ('  I   ;   delta  ;   sigma_e ; alpha  ;  mu  ;  Growth factor of I');
disp(Y)

figure()
subplot(1,2,1)
plot(1:Tdata,log10(exp(target(1:Tdata))),'-.',1:Tdata,log10(exp(LMR1(1:Tdata))),'LineWidth',2)
title('\fontsize{14}Log Mortality'), xlabel('Age')
subplot(1,2,2)
plot(1:Tdata,S0(1:Tdata),'-.',1:Tdata,S1(1:Tdata),'LineWidth',2)
title('\fontsize{14}Survival Curve'), xlabel('Age')
legend('data','with I growth','Location','NorthEast')

% Life expectancy at until age 65
disp('Life expectancy at until age 65')
disp('Data')
life_exp4(target,65)
disp('Simulation ')
life_exp4(LMR1,65)

% Save result
save('US1940F_series.csv','target','S0','LMR1','S1','-ascii','-tabs')

fid=fopen ('US1940F_est.csv','w');
fprintf ('  I   ;   delta  ;   sigma_e ; alpha  ;  mu  ;  Growth factor of I')
fprintf (fid,'%f\n',res_US_W1940_fmin_s');
fclose (fid);
save('US1940F_est.csv','res_US_W1940_fmin_s','-ascii','-tabs')