%% estimate e, phi, and pr(phi > 0) 
%
% Pools data for 66-68 year olds around 1990 and around age 70, for 90-99
function [ests, estBS, estPre, estPreBS] = ...
    estBorn13(bw, opt80, opt90, density, optOuter, p0)
    
    % Load and set data for 1983-1989
    t0 = opt80.t0;
    t1Pre = opt80.t1Pre;
    t1Post = opt80.t1Post;
    ages = opt80.ages;
    agesDensity = opt80.agesDensity;
    
    scale = 1000;
    [data80, bs80, dataPre80, bsPre80] = setDataBorn13(...
        bw, ages, 8389, t0, t1Pre, t1Post, density, agesDensity, scale);
    data80.wMat = eye(length(ages));
    data80.opt = opt80;
    
    % Load and set the data for 1990-1999
    t0 = opt90.t0;
    t1Pre = opt90.t1Pre;
    t1Post = opt90.t1Post;
    ages = opt90.ages;
    agesDensity = opt90.agesDensity;
    
    [data90, bs90, dataPre90, bsPre90] = setDataBorn13(...
        bw, ages, 9099, t0, t1Pre, t1Post, density, agesDensity, scale);
    data90.wMat = eye(length(ages));
    data90.opt = opt90.opt;
    
    data90.wMat = eye(length(ages));
    data90.method = '';
        
    % clean up some to get dynamics properly set
    data80.t0 = ones(size(data80.t1)).*data80.t0;
    data80.zStar = ones(size(data80.t1)).*data80.zStar;
    data80.last = ages == 69;
    data80.ts = ages';
    data80.ts(ages<70) = ages(ages<70)-64;
    data80.ts(ages>=70) = ages(ages>=70)-69;
    data80.method = '';
    
    data90.t0 = ones(size(data90.t1)).*data90.t0;
    data90.zStar = ones(size(data90.t1)).*data90.zStar;
    data90.last = ages == 69;
    data90.ts = ages';
    data90.ts(ages<70) = ages(ages<70)-64;
    data90.ts(ages>=70) = ages(ages>=70)-69;
    
    % put both into a single structure
    data.data70 = data80;
    data.data90 = data90;
    
    % estimation
    [ests, ~, flag] = fminsearch(@(p) ...
        calcGmmPooled(p, data), p0 , optOuter);
    assert (flag == 1) 
    
    % Return bootstrap distribution of estimates, if requested
    if (nargout>=2)
        
        nBoot = size(bs80, 2);
        estBS = zeros(nBoot, length(p0));
        
        dataBS.data70 = data80;
        dataBS.data90 = data90;
        for bs = 1:nBoot;
           dataBS.data70.bunch = bs80(:, bs);
           dataBS.data90.bunch = bs90(:, bs);
           
           [estsT, ~, flag] = fminsearch(@(p) ...
               calcGmmPooled(p, dataBS), ests, optOuter);
           estBS(bs, :) = estsT;
           if (flag ~= 1)
             fprintf('\n Error: bad flag for %g \n', bs)
           end
        end
    end
    
    % set up data for estimating just elasticity
    dataPre.zStar = kron( [dataPre80.zStar ; dataPre90.zStar], ...
        ones(length(dataPre80.bunch), 1) );
    dataPre.t0 = kron( [dataPre80.t0 ; dataPre90.t0], ...
        ones(length(dataPre80.bunch), 1) );
    dataPre.t1 = kron( [dataPre80.t1 ; dataPre90.t1], ...
        ones(length(dataPre80.bunch), 1) );
    dataPre.bunch = [dataPre80.bunch ; dataPre90.bunch];
    dataPre.pre = ones(length(dataPre.bunch));
    
    dataPre.density = density;
    if strcmp(density, 'flexible')
        dataPre.income = data80.income;
        dataPre.dist = data80.dist;
        dataPre.densAK = data80.densAK;        
    end
    % Estimate e using pre-period data
    wMat = eye(length(dataPre.bunch)); 
    [estPre, objVal, flagT] = fminsearch(...
        @(e) calcGmmNoAC(e, dataPre, wMat), ests(1), optOuter);
    assert (flagT == 1 || objVal<1e-6) 

    % calculate bootstrap standard errors
    if (nargout>=4)

        nBoot = size(bsPre80, 2);
        estPreBS = zeros(nBoot, 1);
        dataPreBS = dataPre;

        for bs = 1:nBoot;
            dataPreBS.bunch = [bsPre80(:, bs) ; bsPre90(:,bs)];

            [estT, ~, flagT] = fminsearch(...
                @(e) calcGmmNoAC(e, dataPreBS, wMat), ests(1), optOuter);
            estPreBS(bs, :) = estT;
           
            if (flagT ~= 1)
                fprintf('\n Error: bad flag for %g \n', bs)
            end
        end
    end
            
end
