%% taxReform
%
% Calculate earnings, non-ET tax revenue, and ET revenue, at baseline and 
% for alternative tax systems.

clear
clc


diary(sprintf('../output/tables/tableB6.log'));
fprintf(['Policy reform simulation \n' date '\n']);

% For use in output below
names = {...
    'Base earnings                 ', ...
    'Initially inert share         ', ...
    'Change in earnings            ', ... 
    'Share affected                ', ...
    'Share who adjust              ', ...
    'Change among adjusters        ', ...
    'Percent change among adjusters', ...
};

% Take age 61 earnings distribution as h0.
    zData = importdata('dist6061.csv'); % 1989 distribution
    zGrid = zData.data(:,2);
    zWeight = zData.data(:,1)./sum( zData.data(:,1));
    
% Parameterize, etc.
    tax0 = 0.2721; %average federal MTR + fica, 1989-1990
    tax1 = tax0 + 0.50;
    zStar = 11395;

    data1.t0 = tax0;
    data1.t1 = tax1;
    data1.zStar = zStar;
    
% Only look up to 10k
    zMax = zStar+1e4;
    zGrid = zGrid(zGrid<=zMax);
    zWeight = zWeight(zGrid<=zMax)./sum(zWeight(zGrid<=zMax));
    
    eBase = 0.35;
    eHigh = 0.43;
    eLow = 0.31;
    phi = 278;     
    
%% Simulation 1: Eliminate ET for 62-64 year olds;
    data2 = data1;
    data2.t1 = data2.t0;

    % Case 1: model estimates
        [s1Model, s2Model, sumModel]  = calcStats(eBase,phi, data1, data2, zGrid, zWeight);

    % Case 2: no AC, low e
        [s1Low, s2Low]  = calcStats(eLow,0, data1, data2, zGrid, zWeight);

    % Case 3: no AC, high e
        [s1High, s2High]  = calcStats(eHigh, 0, data1, data2, zGrid, zWeight);

    % Case 4: no AC, actual e
        [s1No, s2No, sumNo]  = calcStats(eBase,0, data1, data2, zGrid, zWeight);

            
    % Output results
        fprintf('A. Eliminate ET for 62-64 year olds \n');
        for rr = 1:7
            fprintf('%s & %7.1f & %7.1f \\\\ \n', names{rr}, sumModel(rr), sumNo(rr));
        end
        fprintf('\n\n');
        

%% Simulation 2: reduce ET brr by 15 percentage points
    data2 = data1;
    data2.t1 = data2.t1-0.17215;

    % Case 1: model estimates
        [s1Model, s2Model, sumModel]  = calcStats(eBase,phi, data1, data2, zGrid, zWeight);

    % Case 2: no AC, low e
        [s1Low, s2Low]  = calcStats(eLow,0, data1, data2, zGrid, zWeight);

    % Case 3: no AC, high e
        [s1High, s2High]  = calcStats(eHigh, 0, data1, data2, zGrid, zWeight);

    % Case 4: no AC, actual e
        [s1No, s2No, sumNo]  = calcStats(eBase,0, data1, data2, zGrid, zWeight);

    % Output results
        fprintf('B: Reduce AET BRR for 62-64 year by 17.22  pppts \n');
        for rr = 1:7
            fprintf('%s & %7.1f & %7.1f \\\\ \n', names{rr}, sumModel(rr), sumNo(rr));
        end
        fprintf('\n\n');

diary off
        