%% f_BSSE
%
% Calculate bootstrap standard errors for (e, p)


function [se, lb, ub, p, estV, jacDet] = calcBSSE6970(estEP, data,wMat,bunchBS, optInner, optOuter)

    % Bootstrap standard errors
    nBoot = size(bunchBS,2);
    estV = zeros(2,nBoot);
    jacDet = zeros(1, nBoot);
    dataBS = data;
    
    for bs = 1:nBoot
        
        dataBS.bunch = ( bunchBS(:,bs)); 
        [est, obj, flag] = fminsearch(@(p) ...
            calcGmm6970([p 0], dataBS, wMat, optInner), estEP, optOuter);
        assert ((flag == 1) || obj < 1e-8)
        estV(:,bs) = est;
        
        % Calculate z cutoffs for determinant
            % pull out parameters
            el = est(1);
            phi = est(2);

            % Calculate zLower
            z0 = data.zStar+phi./(data.taxPre.t1-data.taxPre.t0 +0.00001);
            zb1 =fsolve(@(zOld) uOpt(zOld,el,data.taxPre) ...
                - uFric(zOld,el,data.taxPre) - phi, z0, optInner);

            % Calculate zUpperPre, zUpperPost
            zUpperPre = data.zStar.*( (1-data.t0) ./ (1-data.t1) ).^el;    
            z0 = zUpperPre(data.pre==1); 
            [zUpperPost,~,flag] = fsolve(@(zOld) uOpt(zOld,el,data.taxPost) ...
                - uAK(zOld, el, data.taxPost) - phi, z0, optInner);
            if ~(flag==1)
                z0 = 1/2 * (zb1 + zUpperPre(data.pre==1) );
                [zUpperPost,~,flag] = fsolve(@(zOld) uOpt(zOld,el,data.taxPost) ...
                    - uAK(zOld, el, data.taxPost) - phi, z0, optInner);   
            end
            assert (flag==1)

            % Can't have more bunching post than pre
            tt = zUpperPre(data.pre);
            zUpperPost = min(zUpperPost, tt(end));

            % Get sizes right for combine pre/psot 
            zUpperPost = [zUpperPost ; repmat(zUpperPost(end),sum(data.pre==0),1)];

            zOpt = zUpperPre;
            zb0 = zUpperPost;

        % Calculate determinant    
        if strcmp(data.density, 'flexible')
            hL = interp1(data.income, data.dens, zb1(1));
            hU1 = interp1(data.income, data.dens, zOpt(1));
            hU2 = interp1(data.income, data.dens, zb0(1));
        elseif strcmp(data.density, 'uniform');
            hL = 1;
            hU1 = 1;
            hU2 = 1;
        end
        
        jacDet(bs) = calcDbDiff(est(1), hL, hU1, hU2, ...
            dataBS.t0(1), dataBS.t1(1), dataBS.t1(2), ...
            zb0(1),zb1(1),dataBS.zStar(1));
        
    end
    
    se = [std(estV(1,:)) ; std(estV(2,:))];
    lb = [percentile(estV(1,:)',0.025) ; percentile(estV(2,:)',0.025)];
    ub = [percentile(estV(1,:)',0.975) ; percentile(estV(2,:)',0.975)];
    p = sum(estV<1e-6,2)./nBoot;
    
end