% This script calculates and plots life cycle earning dynamics and generates Figure A1

%% Setup

%Education groups
edus = {'COL','HS'};

%Model age periods
startage = 21;
endage = 65;
periodlenght = 3;

agesModel = startage:periodlenght:endage;

%Data age
agesHubbard = 20:65;

%Initialize variables
earningsLifecycle = zeros(2,1);
earningsLifecycle3y = earningsLifecycle;

FigureA1 = figure('Visible', 'off');


%% Calculate and Plot Lifecycle Earning Patterns

for i = 1:length(edus)

    edu = edus{i};

    %age profile from Hubbard, Skinner, Zeldes (1994): Data Table A.2
    % college = 16+ years of education
    if strcmp(edu,'COL')
        const  = 72270;
        b_age  = -5579;
        b_age2 = 19200 / 100;
        b_age3 = - 18076 / 10000;

        %use Hubbard data Table A.2, high school = 12-15 years
    elseif strcmp(edu,'HS')
        const  = -11833;
        b_age  = 1421;
        b_age2 = -137 / 100;
        b_age3 = - 2186 / 10000;

    else
        warning(['Not a correct education case: edu = ',edu])
    end

    %calculate annual polynomial of average lifecycle earnings from Hubbard et al
    earningsLifecycle(i,1:length(agesHubbard)) = const + b_age.*agesHubbard + b_age2.*agesHubbard.^2 + b_age3.*agesHubbard.^3;

end

% convert annual life-cycle to 3y model period
bb = 1;
for aa = 2:periodlenght:length(earningsLifecycle)
    earningsLifecycle3y(:,bb) = mean(earningsLifecycle(:,aa-1:aa+1),2);
    bb = bb+1;
end

avg_earn = shareCollege.*mean(earningsLifecycle3y(1,:)) + (1-shareCollege).*mean(earningsLifecycle3y(2,:));

earningsMultipliers = earningsLifecycle3y ./ avg_earn;

figure(FigureA1)
plot(agesModel,earningsMultipliers,'x--')
legend(edus)
xlabel('Age')
ylabel('Income Multiplier')


saveas(FigureA1,fullfile(outFolder,'FigureA1.png'))
close all

disp('Success!')