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


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

    % Bootstrap standard errors
    nBoot = size(bunchBS,2);
    estV = zeros(2,nBoot);
    dataBS = data;
    jacDet = zeros(1, nBoot);
    for bs = 1:nBoot
        
        % Find point estimates
        dataBS.bunch = ( bunchBS(:,bs)); 
        [est, ~, flag] = fminsearch( ...
            @(p) calcGmmInert([p 0], dataBS, wMat), estEP, optOuter);
        estV(:, bs) = est;
        assert (flag == 1)        
        
        
        % Calculate determinant         
        [zb1, zb0, zOpt, flag] = calcCutoffs(est, dataBS); 
        
        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 = 2*sum(estV<1e-6,2)./nBoot;
end