%% calcBunching
% calculates bunching given parameters for probabilities , data structure, 
% cut off values of z

function [bSim, bStar, b1Tilde, b2Tilde] =     calcBunching(params, data, z1Bar, z0Bar, uz)
   
    % get bunching static amounts
    bStar = calcBFromZ(data.zStar, uz, data);
    if size(z1Bar, 2) > 1
        z1Bar = z1Bar';
    end
    b1Tilde = calcBFromZ(z1Bar, uz, data);
    z1BarLast = z1Bar(data.last);
    b2Tilde = calcBFromZ(z1BarLast, z0Bar, data);
    db2 = b2Tilde - bStar;
    db1 = bStar - b1Tilde;
    
    % adjust for pr, if appropriate 
    if (length(params)>=3)
        if strcmp(data.method, 'nonparametric')

            %piIX = cumsum(params(3:end));
            %cumPr = (exp(piIX) ./ (1+exp(piIX)))';
            mPr = exp(params(3:end)) ./ (1+exp(params(3:end)));
            cumPr = cumprod(mPr)';
            cumPr = cumPr(data.ts);

        else 
            
            % alpha = 3rd parameter, beta = 4th (if there are 4)
            alpha = params(3);
            if length(params) == 3 
                beta = 0;
            else
                beta = params(4);
            end
            cumPr = data.ts;
            for ts = 1:max(data.ts)
                tt = (1:ts);
                ttpr = exp(alpha + beta.*tt)./ ( 1 + exp(alpha + beta.*tt) );
                cumPr(data.ts==ts) = prod(ttpr);
            end
        end       
        
    else
        cumPr = ones(size(b1Tilde));
    end

    % pre-period bunching
    b1Pre = cumPr.*b1Tilde + (1-cumPr).*bStar;        

    % post-period bunching
    bPost =  b2Tilde - (1-cumPr) .* db2  ...
        + cumPr .* (1-cumPr(data.last)).*db1(data.last);        

    bSim = [b1Pre(data.pre) ; bPost(data.pre==0)];
        
end
