%% estimateCFFric
%
% Estimates e and phi assuming cf earnings affected by frictions

function estimateCFFric(bw, yearsPre, t0, t1Pre, t1Post, density, ...
    optInner, optBounds, optOuter, nDraws, p0, name, ages, scale)

    % Set up data    
    [data, bunchBSB, dataPre, bunchBSPre, dataPost, bunchBSPost] = ...
        setDataCFFric(bw, yearsPre, t0,t1Pre,t1Post, density, ages, nDraws, scale);
    wMat = eye(yearsPre+1);
        
    [estEP, objVal, flag] = fminsearch( ...
        @(ep) calcCFFricGMM(ep, wMat, data, optBounds, optInner, nDraws, density), ...
        p0, optOuter);
    assert (flag == 1);
    
    objValOM = floor(log(objVal)/log(10));
    objValFinal = objVal;
    flagFinal = flag;
    
    % Calculate bootstrap standard errors 
    nBoot = size(bunchBSB,2);
    
    estV = zeros(2,nBoot);
    badSoln = zeros(1,nBoot);
    dataBS = data;
    dataPreBS = dataPre;
    for bs = 1:nBoot
        dataBS.bunch = ( bunchBSB(:,bs)); 
        dataPreBS.bunch = bunchBSB(1:yearsPre,bs);
        
        [estV(:,bs), ~, flag] = fminsearch(@(ep) calcCFFricGMM(ep, ...
            wMat, dataBS, optBounds, optInner, nDraws, density), ...
            estEP, optOuter);
        assert(flag==1);
%         fprintf('.');
%         if bs == 50*floor(bs/50)
%             fprintf('\n')
%         end
    end
    
    
    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-3,2)./nBoot;
    if p(1)<=0.01 
        eStars = '***';
    elseif p(1)<=0.05
        eStars = '** ';
    elseif p(1)<=0.1
        eStars = '*  ';
    else
        eStars = '   ';
    end
    if p(2) <=0.01
        phiStars = '***';
    elseif p(2)<=0.05
        phiStars = '** ';
    elseif p(2)<=0.1
        phiStars = '*  ';
    else
        phiStars = '   ';
    end
    
    % Estimate e using pre-period data
    wMat = eye(yearsPre);
    [estElPre,objVal,flag] = fminunc(@(el) ...
        calcGmmNoAC(el, dataPre, wMat), 0.1,optOuter);
    assert (objVal < 1e-6 || flag==1) 

    % Bootstrap SEs
    [~, ~, ePreLB, ePreUB, estEPre_p] = calcBSSE_e(estElPre, ...
        dataPre,wMat,bunchBSPre, optOuter);

    % Stars
    if estEPre_p <=0.01 
        ePreStars = '***';
    elseif estEPre_p <=0.05
        ePreStars = '** ';
    elseif estEPre_p <=0.1
        ePreStars = '*  ';
    else
        ePreStars = '   ';
    end

    % Estimate e using post-period data
    wMat = 1;
    [estElPost,objVal,flag] = fminunc(@(el) calcGmmNoAC(el, dataPost, wMat), 0.1,optOuter);
    assert (objVal < 1e-6 || flag==1) 

    % Bootstrap SEs
    [~, ~, ePostLB, ePostUB, estEPost_p] = calcBSSE_e(estElPost, ...
        dataPost,wMat,bunchBSPost, optOuter);

    % Stars
    if estEPost_p <=0.01 
        ePostStars = '***';
    elseif estEPost_p <=0.05
        ePostStars = '** ';
    elseif estEPost_p <=0.1
        ePostStars = '*  ';
    else
        ePostStars = '   ';
    end
            
    % Output results 
    fprintf(' %s: &     %4.2f      &     %6.0f      ',name, estEP(1), scale*estEP(2))
    fprintf(' &     %4.2f       &     %4.2f       \\\\ \n', estElPost, estElPre);

    fprintf('           & [%4.2f, %4.2f]%s & [%6.0f, %6.0f]%s ', lb(1), ub(1), eStars, scale*lb(2), scale*ub(2), phiStars);
    fprintf(' & [%4.2f, %4.2f]%s  & [%4.2f, %4.2f]%s  \\\\ \n', ePostLB, ePostUB, ePostStars, ePreLB, ePreUB, ePreStars);
    %fprintf('\n Obj val (order of magnitude) & %5.1f \n', objValOM);
    %fprintf('Long e  : %10.5f \n', estEP(1));
    %fprintf('Long phi: %10.5f \n', estEP(2));
    %fprintf('Flag: %g\n', flagFinal);
    %fprintf('Obj val final: %g \n', objValFinal);
    if (sum(badSoln))
        fprintf('%g bad solutions \n', sum(badSoln))
    end
end