clear all; clc 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. load parameters and data and simulate the model
load('Results\est')
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2. Estimated model parameters
names.Highskilled = 'High skilled';
names.PoorHealth = '10+ GP visits';
names.ExpWorseHealth = 'Expect worse health';
names.Income      = 'Labor income (1,000\pounds)';
names.pensPPP     = 'Has private pension (PPP)';
names.pensEPS     = 'Has employer provided pension (EPS)';
names.cohort      = 'Birth year (minus 1955)';
names.RetAge_cons = 'Constant';
names.RetAge_trend = 'Time trend (minus 25)';
names.RetAge55 = 'Retirement age 55 dummy';
names.RetAge60 = 'Retirement age 60 dummy';
names.RetAge65 = 'Retirement age 65 dummy';
names.Wealth = 'Wealth';

diary off
delete(['Results\Estimates.tex'])
diary(['Results\Estimates.tex'])

fprintf('\\begin{tabular}{clrrrr} \\toprule \n');
fprintf('&& \\multicolumn{2}{c}{Husband} & \\multicolumn{2}{c}{Wife} \\\\ \\hline \n');

fprintf('$\\gamma$ & Joint leisure & $%2.3f$ & ($%2.3f$) & $%2.3f$ & ($%2.3f$) \\\\ \n',par.joint_h,par.SE(1),par.joint_h,par.SE(1));
fprintf('$\\alpha$ & SPA age & \\multicolumn{1}{c}{$-$} & \\multicolumn{1}{c}{$-$} & $%2.3f$ & ($%2.3f$) \\\\ \n',par.SPA_w,par.SE(2));
fprintf('\\multicolumn{6}{l}{} \\\\ \n');

fprintf('\\multicolumn{6}{l}{\\emph{Explanatory variables ($\\beta$)}} \\\\ \n');
num_own = numel(par.ExplVar.h.Own);
num_spouse = numel(par.ExplVar.h.Spouse);
num = num_own+num_spouse;
SE_h = par.SE(3:2+num);
SE_w = par.SE(3+num:2+2*num);
for p=1:num_own
    name = par.ExplVar.h.Own{p};
    name = names.(name);
    fprintf('\\multicolumn{2}{l}{%s} & $%2.3f$ & ($%2.3f$) & $%2.3f$ & ($%2.3f$) \\\\ \n',name,par.beta_h(p),SE_h(p),par.beta_w(p),SE_w(p));
end

for p=num_own+1:num
    name = par.ExplVar.h.Spouse{p-num_own};
    name = names.(name);
    fprintf('\\multicolumn{2}{l}{%s, spouse} & $%2.3f$ & ($%2.3f$) & $%2.3f$ & ($%2.3f$) \\\\ \n',name,par.beta_h(p),SE_h(p),par.beta_w(p),SE_w(p));
end
fprintf('\\multicolumn{6}{l}{} \\\\ \n');

fprintf('\\multicolumn{6}{l}{\\emph{Age variables ($\\delta$)}} \\\\ \n');
num_age = 5;

SE = par.SE(3+2*num:end);
SE_h = SE(1:2:num_age*2);
SE_w = SE(2:2:num_age*2);
for p=1:num_age
    name = par.EstPar{4+p};
    name_p = names.(name);
    fprintf('\\multicolumn{2}{l}{%s} & $%2.3f$ & ($%2.3f$) & $%2.3f$ & ($%2.3f$) \\\\ \n',name_p,par.(name)(1),SE_h(p),par.(name)(2),SE_w(p));
end
fprintf('\\multicolumn{6}{l}{} \\\\ \n');

T       = [par.var_h, 0;
           par.cov  ,  par.var_w]; 
Omega   = T*T';
fprintf('$\\sigma$ & variance & $%2.3f$ && $%2.3f$ &\\\\ \n',Omega(1,1),Omega(2,2));
fprintf('$\\sigma_{hw}$ & covariance & $%2.3f$ & & $%2.3f$ & \\\\ \n',Omega(1,2),Omega(1,2));
% fprintf('\\multicolumn{2}{l}{Objective} & \\multicolumn{4}{c}{$%2.6f$} \\\\ \n',par.fval);
fprintf('\\bottomrule \\end{tabular} \n');

diary off
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3. model-fit

seed = 9210; % seed used to draw the normal errors for both spouses
S    = 2000;%5000;  % number of simulated shocks - used to evaluate the integral

data = fun.load_data(seed,S);
data = fun.explanatory_vars(data,par.ExplVar,par.ExplVar_excl,par.EstPar);

par.age_h = par.RetAgeDum_fun(par.age_min:par.age_max,1,par); 
par.age_w = par.RetAgeDum_fun(par.age_min:par.age_max,2,par);
[ExpRetAge.h,ExpRetAge.w]   = fun.calc_ret_combi_mex(data,par);
  

%% Histograms
members = {'h','w'};
FontSize = 15;
for m=1:numel(members)
figure(m)
set(gca,'FontSize',FontSize)
x = 50:1:70;
[h1]=hist(ExpRetAge.(members{m})(:),x);
[h2]=hist(data.(members{m}).ExpRetAge(:),x);
b=bar(x,[h2./sum(h2);h1./sum(h1)]');
%b=bar(x,[h2./sum(h2);h1./sum(h1)]','FaceColor','flat');
%b(1).FaceColor = .0*ones(1,3);
%b(2).FaceColor = .9*ones(1,3);
xlim([49 71])
if m==1
    legend('Data','Model','Location','NorthWest')
    title('Men')
else
    title('Women')
end
xlabel('Planned retirement age')
ylabel('share')
set(gca,'FontSize',FontSize)
print('-depsc',['Results\Fit' members{m}])
end


%%
ret_year_h_sim = ExpRetAge.h - data.h.Age;
ret_year_w_sim = ExpRetAge.w - data.w.Age;
%ret_diff_sim = mean(ret_year_h_sim - ret_year_w_sim,2);
ret_diff_sim = ret_year_h_sim - ret_year_w_sim;

ret_year_h = data.h.ExpRetAge - data.h.Age;
ret_year_w = data.w.ExpRetAge - data.w.Age;
ret_diff_data  = ret_year_h-ret_year_w;
ret_diff_data(isnan(ret_diff_sim(:,1))) = NaN; % add NaNs where they are in the model due to missing information on X's

figure(1)
histogram(ret_diff_data,-30.5:30.5,'Normalization','probability')
hold on
histogram(ret_diff_sim(:),-30.5:30.5,'Normalization','probability')
hold off
set(gca,'FontSize',FontSize)
xlabel('Planned retirement year difference')
ylabel('share')
legend('Data','Model')
set(gca,'FontSize',FontSize)
print('-depsc','Results\Fit_joint')

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 5. Sensitivity measures for gamma
pari = 1; % gamma is the first parameter

[num_par,num_mom] = size(par.sens.M1);

par.sens.M1s = abs(par.sens.M1e)./max(abs(par.sens.M1e),[],2);
par.sens.M2s = abs(par.sens.M2e)./max(abs(par.sens.M2e),[],2);
par.sens.M3s = abs(par.sens.M3e)./max(abs(par.sens.M3e),[],2);
par.sens.M4s = abs(par.sens.M4e)./max(abs(par.sens.M4e),[],2);
par.sens.M5s = abs(par.sens.M5e)./max(abs(par.sens.M5e),[],2);
par.sens.M6s = abs(par.sens.M6e)./max(abs(par.sens.M6e),[],2);

vec_measures ={{'M1','M2','M3','M4','M5','M6'} , {'M1e','M2e','M3e','M4e','M5e','M6e'}, {'M1s','M2s','M3s','M4s','M5s','M6s'} };
vec_measures_name = { {'M_1','M_2','M_3','M_4','M_5','M_6'} , {'\mathcal{E}_1','\mathcal{E}_2','\mathcal{E}_3','\mathcal{E}_4','\mathcal{E}_5','\mathcal{E}_6'},...
    {'\overline{M}_1','\overline{M}_2','\overline{M}_3','\overline{M}_4','\overline{M}_5','\overline{M}_6'} };
vec_name = {['Results\sens_raw.tex'],['Results\sens_e.tex'],['Results\sens_s.tex']};


measures        = vec_measures{2};
measures_name   = vec_measures_name{2};

name = ['Results\sens_e1.tex'];
reg_str = {'Constant','High skilled, husband','10+ GP visits, husband','Expect worse health, husband','Labor income, husband','Has private pension, husband','Has employer provided pension, husband','Birth year (minus 1955), husband',...
        'High skilled, wife','10+ GP visits, wife','Expect worse health, wife','Labor income, wife','Has private pension, wife','Has employer provided pension, wife','Birth year, wife','Birth year, wife in 1951--1955','Birth year, wife later than 1955'};
first = {'Regression, husband','Regression, wife'};

diary off
delete(name)
diary(name)

fprintf('\\begin{tabular}{ll*{%d}{r}} \n \\toprule \\multicolumn{2}{l}{Moment} ',numel(measures));
for m=1:numel(measures), fprintf('& \\multicolumn{1}{c}{$%s$} ',measures_name{m}); end 
fprintf('\\\\ \\cmidrule(lr){3-%d} \n',numel(measures)+2);

for i=1:numel(first)
    fprintf('\\multicolumn{2}{l}{\\emph{%s}} & \\multicolumn{%d}{c}{} \\\\ \n',first{i},numel(measures))
    
    for j=1:numel(reg_str)
        id = j+(i-1)*numel(reg_str);
        fprintf('%d & %s ',id , reg_str{j});
        for m=1:numel(measures), fprintf('& $%2.3f$ ',par.sens.(measures{m})(pari,id)); end
        fprintf('\\\\ \n');
    end
end
fprintf('\\bottomrule \\end{tabular}\n');
diary off

id_now = j+(i-1)*numel(reg_str);

name = ['Results\sens_e2.tex'];
reg_str = {'Share at ages 50--54','Share at age 55','Share at ages 56--59','Share at age 60','Share at ages 61--64','Share at age 65'};
first = {'Retirement age, husband','Retirement age, wife'};

diary off
delete(name)
diary(name)

fprintf('\\begin{tabular}{ll*{%d}{r}} \n \\toprule \\multicolumn{2}{l}{Moment} ',numel(measures));
for m=1:numel(measures), fprintf('& \\multicolumn{1}{c}{$%s$} ',measures_name{m}); end 
fprintf('\\\\ \\cmidrule(lr){3-%d} \n',numel(measures)+2);

for i=1:numel(first)
    fprintf('\\multicolumn{2}{l}{\\emph{%s}} & \\multicolumn{%d}{c}{} \\\\ \n',first{i},numel(measures))
    
    for j=1:numel(reg_str)
        id = id_now+j+(i-1)*numel(reg_str);
        fprintf('%d & %s ',id , reg_str{j});
        for m=1:numel(measures), fprintf('& $%2.3f$ ',par.sens.(measures{m})(pari,id)); end
        fprintf('\\\\ \n');
    end
end

id_now = id_now+j+(i-1)*numel(reg_str);
reg_str = {'var$(e_h)$','var$(e_w)$','cov$(e_h,e_w)$','diff [-2,-1]','diff [1,2]','Joint retirement'};
fprintf('\\multicolumn{2}{l}{\\emph{Simultaneous retirement}} & \\multicolumn{%d}{c}{} \\\\ \n',numel(measures))
for j=1:numel(reg_str)
    id = id_now+j;
    fprintf('%d & %s ',id , reg_str{j});
    for m=1:numel(measures), fprintf('& $%2.3f$ ',par.sens.(measures{m})(pari,id)); end
    fprintf('\\\\ \n');
end

fprintf('\\bottomrule \\end{tabular}\n');
diary off

%%
t = par.joint_h/par.SE(1)
p = 2*(1-normcdf(t))