%% calcStats
% 
% Calculates counterfactual statistics under new tax policy

function [stats1, stats2, statsSum] = calcStats(e,phi, data1, data2, zGrid, zWeight)

% Optimal earnings under t1
    deltaZStar1 = data1.zStar.*( (1-data1.t0)/(1-data1.t1) )^e;
    z1Opt = zGrid.*(zGrid<data1.zStar) + ...
        data1.zStar.*(zGrid>=data1.zStar).*(zGrid<=deltaZStar1) + ...
        zGrid.*(zGrid>deltaZStar1).*( (1-data1.t1)/(1-data1.t0) )^e;

% Frictional earnings under t1
    uOptGrid = uOpt(zGrid, e, data1);    
    nGrid = zGrid./( (1-data1.t0).^e );
    uDiffGrid = calcUDiff(uOptGrid, zGrid, nGrid, e, phi, data1);
    z1 = z1Opt.*(uDiffGrid>=phi) + zGrid.*(uDiffGrid<phi);

% Optimal earnings under t2
    z2Linear = nGrid.*(1-data2.t0).^e;
    deltaZStar2 = data2.zStar.*( (1-data2.t0)/(1-data2.t1) )^e;
    z2Opt = z2Linear.*(zGrid<data2.zStar) + ...
        data2.zStar.*(z2Linear>=data2.zStar).*(z2Linear<=deltaZStar2) + ...
        z2Linear.*(z2Linear>deltaZStar2).*( (1-data2.t1)/(1-data2.t0) )^e;

% Frictional earnings under t2
    uOpt2Grid = uQL(z2Opt, nGrid, e, data2);
    uDiff2Grid = calcUDiff(uOpt2Grid, z1, nGrid, e, phi, data2);       
    z2 = z1.*(uDiff2Grid<0) + z2Opt.*(uDiff2Grid>=0);
    
% Calculate revenue under baseline and new taxes
    taxRevenue1Grid = z1.*data1.t0;
    ssRevenue1Grid = (z1-data1.zStar).*(z1>data1.zStar).*(data1.t1-data1.t0);
    taxRevenue2Grid = z2.*data2.t0;
    ssRevenue2Grid = (z2-data1.zStar).*(z2>data2.zStar).*(data2.t1-data2.t0);
    
% Integrate up to get revenue, etc.
    earnings1 = zWeight'*z1;
    taxRevenue1 = zWeight'*taxRevenue1Grid;
    ssRevenue1 = zWeight'*ssRevenue1Grid;
    stats1 = [earnings1 ; taxRevenue1; ssRevenue1];
    
    earnings2 = zWeight'*z2;
    taxRevenue2 = zWeight'*taxRevenue2Grid;
    ssRevenue2 = zWeight'*ssRevenue2Grid;
    stats2 = [earnings2 ; taxRevenue2 ; ssRevenue2];
    
    inert1 = zWeight'*(abs(z1-z1Opt)>1e-6);
    
    changeAbove = data1.t1~=data2.t1;
    changeBelow = data1.t0~=data2.t0;
    shareAbove = zWeight'*(zGrid>data1.zStar);
    shareBelow = zWeight'*(zGrid<data1.zStar);
    shareAffected = shareAbove*changeAbove + shareBelow*changeBelow;
    adjusted = logical(abs(z2-z1)>1e-6);
    shareAdjusting = zWeight'*adjusted;
    
    if shareAdjusting > 0
        deltaCop = (earnings2-earnings1)/shareAdjusting;
        percentCop = (zWeight(adjusted)'*(z2(adjusted)./z1(adjusted)-1))/shareAdjusting;
    else
        deltaCop = 0;
        percentCop = 0;
    end
    
    statsSum = [earnings1 ; 100*inert1 ; earnings2-earnings1 ; ...
        100*shareAffected ; 100*shareAdjusting ; deltaCop; 100*percentCop];
