%% This file generates graphs and tables from the structural model
% All inputs are provided, and can otherwise be computed by running
% main.m with the appropriate choices.

% For BPP-coefficient related statistics, see bpp_coefficients.R

clear 

% Load all relevant outputs from the several different runs of the model
NL=load('Outputs/case1');
Ca=load('Outputs/case2');
D1=load('Outputs/case3');
KO=load('Outputs/case4');

Ca_cb=load('Outputs/case2_cb');
D1_cb=load('Outputs/case3_cb');
KO_cb=load('Outputs/case4_cb');

Robust_tr=load('Outputs/case5');
Robust_per=load('Outputs/case6');
Robust_can=load('Outputs/case7');

% Load the data for the variance of log earnings and variance of log
% consumption that have been computed in the "Earnings" folder
dataearn=csvread('../Earnings/Outputs/earnvarprofile_age_kaplan.csv',1,1);
earnvari=smooth(dataearn(2:37,2),5);
earnvark=[earnvari(3) earnvari(8) earnvari(13) earnvari(18)...
    earnvari(23) earnvari(28) earnvari(33)]; %5y averages to compare

datakar=csvread('../Earnings/Outputs/data_varlogcons_ka_reg.csv',1,0);
ages=datakar(:,2);
datanormkar= datakar(:,1)-datakar(1,1);

% Figure 6

figure
plot(ages,earnvark,'k-',ages,datakar,'k--')
legend('Earnings','Consumption','Location','NorthWest')
xlabel('Age')
ylabel('Variance')
axis([25 60 0.20 0.55])
set(gca,'FontSize',16)


% Figure 7

figure    
plot(ages,datanormkar,'k:',...
    ages,Ca.lvarconsgroup-Ca.lvarconsgroup(1),'k--',...
    ages,NL.lvarconsgroup-NL.lvarconsgroup(1),'k');
legend('CEX Data','Canonical','NL process','Location','NorthWest')
xlabel('Age')
ylabel('Var. of log consumption (age t-age 25)')
axis([25 60 -0.04 0.1])
set(gca,'FontSize',14)


% Figure 8

figure    
plot(ages,datanormkar(:,1),'k:',...
    ages,Ca.lvarconsgroup-Ca.lvarconsgroup(1),'k--',...
    ages,KO.lvarconsgroup-KO.lvarconsgroup(1),'ko-',...
    ages,D1.lvarconsgroup-D1.lvarconsgroup(1),'kx-',...
    ages,NL.lvarconsgroup-NL.lvarconsgroup(1),'k-');
legend('CEX Data','Canonical',...
    'Normal, age-dependent','Non-normal, age-dependent',...
    'NL process',...
    'Location','NorthWest')
xlabel('Age')
ylabel('Var. of log consumption (age t-age 25)')
axis([25 60 -0.05 0.15])
set(gca,'FontSize',14)


% Figures 9a and 9b

figure    
plot(25:60,Ca.savbyage,'k--',...
    25:60,KO.savbyage,'ko-',...
    25:60,D1.savbyage,'kx-',...
    25:60,NL.savbyage,'k-');
legend('Canonical',...
    'Normal, age-dependent','Non-normal, age-dependent',...
    'NL process',...
    'Location','Best')
xlabel('Age')
ylabel('Average wealth')
set(gca,'FontSize',14)


figure    
plot(25:60,Ca.varage,'k--',...
    25:60,KO.varage,'ko-',...
    25:60,D1.varage,'kx-',...
    25:60,NL.varage,'k-');
legend('Canonical',...
    'Normal, age-dependent','Non-normal, age-dependent',...
    'NL process',...
    'Location','NorthWest')
xlabel('Age')
ylabel('Variance of wealth')
set(gca,'FontSize',14)


% Figure D.6.

figure    
plot(ages,datanormkar,'k:',...
    ages,NL.lvarconsgroup-NL.lvarconsgroup(1),'k',...
    ages,Robust_tr.lvarconsgroup-Robust_tr.lvarconsgroup(1),'kd-');
legend('Data',...
    '8 gridpoints','16 gridpoints','Location','NorthWest')
xlabel('age')
ylabel('calib var of log consumption (var in t - var at 25)')

% Figure D.7.

figure    
plot(ages,datanormkar,'k:',...
    ages,NL.lvarconsgroup-NL.lvarconsgroup(1),'k',...
    ages,Robust_per.lvarconsgroup-Robust_per.lvarconsgroup(1),'kd-');
legend('Data',...
    '18 gridpoints','36 gridpoints','Location','NorthWest')
xlabel('age')
ylabel('calib var of log consumption (var in t - var at 25)')

% Figure D.8.

figure    
plot(ages,datanormkar,'k:',...
    ages,Ca.lvarconsgroup-Ca.lvarconsgroup(1),'k--',...
    ages,NL.lvarconsgroup-NL.lvarconsgroup(1),'k',...
    ages,Robust_can.lvarconsgroup-Robust_can.lvarconsgroup(1),'kd-');
legend('Data','Canonical',...
    'NL process','Canonical (alternative discret.)','Location','NorthWest')
xlabel('age')
ylabel('calib var of log consumption (var in t - var at 25)')

% Display wealth distributions

disp('Wealth distribution, canonical process')
disp(Ca.percentable)
disp('Wealth distribution, NL process')
disp(NL.percentable)

% Display consumption equivalents (Table 4)

disp(NL.xi)
disp(Ca_cb.xi)
disp(KO_cb.xi)
disp(D1_cb.xi)