%% Elasticities
%
% Returns
%
%   year, year
%   bunch, bunching amount
%   bSE, bunching SE
%   el, elasticity implied by bunching amount
%   elSE, delta-method SE. 
%
% Input: year range 

function [year, bunch,bSE, el, elSE, elFlex] = calcEl(bw,yStart,yStop,t0,t1Pre,t1Post ,opt)

    %yearStatsDensPoly_7803_6668_800_s1Main
    rawBunch = importdata(sprintf('yearStatsDensPoly_7802_6668_%g_s1Main.txt',bw));
    rawZStar = importdata('realZStar.csv');
        
    bunchYear = rawBunch.data(:,1);
    
    % align bunch and bunch bs matrices
    bbIndex =(bunchYear>=yStart).*(bunchYear<=yStop);
    bunch = bw/1e3.*rawBunch.data(bbIndex==1,2);
    year = bunchYear(bbIndex==1);
    bSE = bw/1e3.*(rawBunch.data(bbIndex==1,3));
    zStarPost = rawZStar.data(:,3)/1e3;
    zStarYear = rawZStar.data(:,1);

    % Set zStar and t1; t0 fixed as input.
    zStar = zStarPost( (zStarYear>=yStart).*(zStarYear<=yStop)==1 );
    t1 = t1Pre.*(bunchYear(bbIndex==1)< 1990) + t1Post.*(bunchYear(bbIndex==1)>=1990);
    
    % Calculate elasticity
    el = log( bunch./zStar + 1) ./ log( (1-t0)./(1-t1) );
    
    % delta method standard errors
    dedb = 1 ./( (bunch + zStar).*log( (1-t0)./(1-t1) ) );
    elSE = bSE.*dedb;
    
    % calcalate flexible elasticity, if called for
    if nargout >= 5
       
        % get the density data
        data.density = 'flexible';
        data.t0 = t0;
        densityData = importdata(sprintf('dens%s_8399_%g_s1Main.csv', '72', bw));
        prob = densityData.data(:,1);
        %prob = linspace(0, 1, length(data.dist)); 
        data.dist = cumsum(prob);      
        zStarPoint = logical(densityData.data(:,2) == 0);
        data.densAK = prob(zStarPoint)/(bw/1000);    
        
        % loop over years
        elFlex = yStart:yStop;        
        for yy = 1:length(elFlex)
            
            % update year-specific parameters
            data.t1 = t1(yy);
            data.zStar = zStar(yy);
            data.income = densityData.data(:,2) + data.zStar;
            data.bunch = bunch(yy);
            
            % Estimate 
            [est, ~, flag] = fzero(@(e0) calcErrDirect(e0, data),[0.001 5], opt);
            assert (flag == 1)
            elFlex(yy) = est;

        end
    end
    
end


