%% estimateGMM
%
% Estimate the elasticity and phi via gmm. 
%
% 
% Use data from 1986-1994, on 66-69 year olds, transition from 50% phaseout
% rate to 33%  rate, or from 33% to 0. 
%
% Baseline MTR rates are estimated using TAXSIM on SOI data, see the
% program "mtr030813.do" and the log "computeMTR.log." We use the average
% baseline taxrate for the years underlying each estimator.
%% preliminaries
clear;
clc; 

if exist(sprintf('../output/tables/structuralEstimates%s.log',date),'file'); 
    sprintf('structuralEstimates%s.log',date); 
end
diary(sprintf('../output/tables/structuralEstimates%s.log',date));
fprintf(['Estimate GMM -- many year/sample/tax combinations \n' date '\n']);

optInner = optimset('tolX', 1e-10, 'tolF', 1e-10, ... 
    'MaxFunEvals', 10000, 'MaxIter', 5000, ...
        'Disp','off','LargeScale','off');
optDyn =  optimset('tolX', 1e-10, 'tolF', 1e-10, ...
    'MaxFunEvals', 10000, 'MaxIter', 5000, 'Disp','off');
optOuter = optimset('tolX', 1e-6, 'tolF', 1e-6, 'MaxFunEvals', 10000, ...
        'MaxIter', 5000, 'Disp', 'off','LargeScale','off'); 
optLoud = optimset('tolX', 1e-6, 'tolF', 1e-6, ...
    'MaxFunEvals', 50000, 'MaxIter', 10000, 'Disp', 'Iter-Detailed'); 
optTight = optimset('tolX', 1e-8, 'tolF', 1e-8, 'MaxFunEvals', 10000, ...
        'MaxIter', 5000, 'Disp','off','LargeScale','off');
optBounds = optimset('tolX', 1e-10, 'tolF', 1e-10, ...
    'MaxFunEvals', 10000, 'MaxIter', 5000, 'Disp','off','Jacobian','on');


rng(3);

%% Plot elasticities (Figure 8)
bw = 800;
yStart = 1982; yStop = 1993;
t0 = .2721; t1Pre = t0+0.5 ; t1Post = t0+0.33; % include ficas
[year, bunch, bSE, el, elSE] = calcEl(bw,yStart,yStop,t0,t1Pre,t1Post, optOuter);

fprintf('\n\n')
fprintf('Bunching and elasticities underlying figure 8\n\n')
fprintf(' Year & Bunch & BSE & El & ElSE \n');
for ix = 1:length(year)
    fprintf(' %4.0f & %4.3f & % 4.3f &  %4.3f & % 4.3f \n',...
        year(ix), bunch(ix), bSE(ix), el(ix), elSE(ix));
end

upper = el + 1.96*elSE;
lower = el - 1.96*elSE;

h=figure();
axis();
colorMat = [0.2 0.3 0.4 ; 1 0.5 0.5 ; 1 0.5 0.5];
set(gca, 'colororder', colorMat, 'nextplot', 'replacechildren');
plot(year, el, '-', ...
    year, upper, '--', ...
    year, lower, '--', ...
    'LineWidth', 2);
xlabel('Year')
ylabel('Elasticity')

print(h,sprintf('../output/figures/fig8.pdf'), '-dpdf');



%% point estimate for 1989-1990 transition: sharp change method (TABLE 2)


fprintf('\n\n\n\n\n')
fprintf('Table 2')
fprintf('\n\n\n\n\n')
    tic;
    fprintf('Elasticity and adjustment cost estimates from 1989-1990 change in BRR \n');
    fprintf('           &     (1)        &      (2)       &         (3)     &       (4) \\\\ \\midrule');
    fprintf('           & $\\varepsilon$ &    $\\phi$     &     \\multicolumn{2}{c}{$\\varepsilon |\\phi=0$}    \\\\ \n')
    fprintf(' \\cmidrule(lr){2-2} \n');
    fprintf(' \\cmidrule(lr){3-3} \n');
    fprintf(' \\cmidrule(lr){4-5} \n');
    fprintf('           &                &                &      1990       &      Pre        \\\\ \n')
    fprintf(' \\cmidrule(lr){4-4} \n');
    fprintf(' \\cmidrule(lr){5-5} \n');

    t0Base = 0.2721; %average federal MTR + fica, 1989-1990
    scale = 1000;
    p0 = [0.25 0.15];

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% 66-68 year olds, 1989-1990, s1Main; baseline
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    yearsPre = 1;
    t0 = t0Base;
    t1Pre = 0.5 + t0;
    t1Post = 0.33 + t0;
    bw = 800;
    density = 'flexible';
    name = 'Baseline';
    ages = '72';
    est1990(bw,yearsPre,t0,t1Pre,t1Post,density,optInner,optOuter,p0,name, ages, scale);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% 66-68 year olds, 1989-1990, s1Main, uniform density
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    yearsPre = 1;
    t0 = t0Base ;
    t1Pre = 0.5 + t0;
    t1Post = 0.33 + t0;
    bw = 800;
    name = 'Uniform Density';
    density = 'uniform';
    ages = '72';
    est1990(bw,yearsPre,t0,t1Pre,t1Post,density,optInner,optOuter,p0,name, ages, scale);       

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~percen~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% 66-68 year olds, 1989-1990, s1Main, t1 = .36/.24
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    yearsPre = 1;
    t0 = t0Base ;
    t1Pre = 0.36 + t0;
    t1Post = 0.24 + t0;
    bw = 800;
    density = 'flexible';
    name = 'Benefit Enhancement ';
    ages = '72';
    est1990(bw,yearsPre,t0,t1Pre,t1Post,density,optInner,optOuter,p0,name, ages, scale);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% 66-68 year olds, 1989-1990, s1Main, exclude fica
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    yearsPre = 1;
    t0 = 0.1205;
    t1Pre = 0.5 + t0;
    t1Post = 0.33 + t0;
    bw = 800;
    density = 'flexible';
    name = 'Excluding FICA ';        
    ages = '72';
    est1990(bw,yearsPre,t0,t1Pre,t1Post,density,optInner,optOuter,p0,name, ages, scale);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% 66-68 year olds, 1989-1990, s1Main, bw = 400
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    yearsPre = 1;
    t0 = t0Base ;
    t1Pre = 0.5 + t0;
    t1Post = 0.33 + t0;
    bw = 400;
    density = 'flexible';
    name = 'Bandwidth = \$400';        
    ages = '72';
    est1990(bw,yearsPre,t0,t1Pre,t1Post,density,optInner,optOuter,p0,name, ages, scale);       

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% 66-68 year olds, 1989-1990, s1Main, bw = 1600
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    yearsPre = 1;
    t0 = t0Base ;
    t1Pre = 0.5 + t0;
    t1Post = 0.33 + t0;
    bw = 1600;
    density = 'flexible';
    name = 'Bandwidth = \$1600';        
    ages = '72';
    est1990(bw,yearsPre,t0,t1Pre,t1Post,density,optInner,optOuter,p0,name, ages, scale);       

toc;



%% Elasticities for age 69-70 transition (TABLE 3)

    fprintf('\n\n\n\n\n')
    fprintf('Table 3')
    fprintf('\n\n\n\n\n')

    fprintf('Elasticity and adjustment cost estimates for age 69-70 transition, SC \n');
    fprintf('           &     (1)        &      (2)       &         (3)     \\\\ \\midrule');
    fprintf('           & $\\varepsilon$ &    $\\phi$     & $\\varepsilon |\\phi=0$, Age 69    \\\\ \n')
    fprintf(' \\cmidrule(lr){2-2} \n');
    fprintf(' \\cmidrule(lr){3-3} \n');
    fprintf(' \\cmidrule(lr){4-4} \n');
    
    scale = 1000;
    p0 = [0.25 150/scale];
    t0Base = 0.2783 ; %average federal MTR + fica
    tic;
    
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Baseline
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        ages = 69:71;
        t0 = t0Base;
        t1Pre = 0.33 + t0;
        t1Post = 0 + t0;
        bw = 800;
        density = 'flexible';
        name = 'Baseline';
        agesDensity = '72';
        est6970(bw,ages,t0,t1Pre,t1Post,density,optInner,optOuter,p0,name, agesDensity, scale)

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Uniform
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        ages = 69:71;
        t0 = t0Base;
        t1Pre = 0.33 + t0;
        t1Post = 0 + t0;
        bw = 800;
        density = 'uniform';
        name = 'Uniform Density';
        agesDensity = '72';
        est6970(bw,ages,t0,t1Pre,t1Post,density,optInner,optOuter,p0,name, agesDensity, scale)
        
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Benefit enhancement
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        ages = 69:71;
        t0 = t0Base;
        t1Pre = 0.24 +t0Base;
        t1Post = 0 +t0Base;
        bw = 800;
        density = 'flexible';
        name = 'Benefit Enhancement ';
        agesDensity = '72';
        est6970(bw,ages,t0,t1Pre,t1Post,density,optInner,optOuter,p0,name, agesDensity, scale)
         
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Exclude FICA
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        ages = 69:71;
        t0 =  0.1253; % average federal MTR
        t1Pre = 0.33 + t0;
        t1Post = 0 + t0;
        bw = 800;
        density = 'flexible';
        name = 'Excluding FICA';
        agesDensity = '72';
        est6970(bw,ages,t0,t1Pre,t1Post,density,optInner,optOuter,p0,name, agesDensity, scale)
        
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % BW = 400
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        ages= 69:71;
        t0 = t0Base ;
        t1Pre = 0.33 +t0;
        t1Post = 0 + t0;
        bw = 400;
        density = 'flexible';
        name = 'Bandwidth = \$400';
        agesDensity = '72';
        est6970(bw,ages,t0,t1Pre,t1Post,density,optInner,optOuter,p0,name, agesDensity, scale)
        
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % BW = 1600
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        ages= 69:71;
        t0 = t0Base ;
        t1Pre = 0.33 +t0;
        t1Post = 0 + t0;
        bw = 1600;
        density = 'flexible';
        name = 'Bandwidth = \$1600';
        agesDensity = '72';
        est6970(bw,ages,t0,t1Pre,t1Post,density,optInner,optOuter,p0,name, agesDensity, scale)
        
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % 68-70 year olds
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        ages = 68:70;
        t0 = t0Base ;
        t1Pre = 0.33 + t0;
        t1Post = 0 + t0;
        bw = 800;
        density = 'flexible';
        name = '68-70 year-olds';
        agesDensity = '72';
        est6970(bw,ages,t0,t1Pre,t1Post,density,optInner,optOuter,p0,name, agesDensity, scale)
           
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % 69 and 71 year olds only
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        ages = [69 71];
        t0 = t0Base ;
        t1Pre = 0.33 + t0;
        t1Post = 0 + t0;
        bw = 800;
        density = 'flexible';
        name = '69, 71 year-olds';
        agesDensity = '72';
        est6970(bw,ages,t0,t1Pre,t1Post,density,optInner,optOuter,p0,name, agesDensity, scale)
        
        
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Set up for born Jan-March only 
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
        p0 = [0.25 0.15];
        t0Base = 0.2783 ; %average federal MTR + fica
        % 80's data
        opt80.ages = 69:71;
        opt80.t0 = t0Base;
        opt80.t1Pre = 0.5 + opt80.t0;
        opt80.t1Post = 0 + opt80.t0;
        opt80.agesDensity = '72';
        opt80.opt = optDyn;
        
        % 90's data
        opt90.ages = 69:71;
        opt90.t0 = t0Base;
        opt90.t1Pre = 0.33 + opt90.t0;
        opt90.t1Post = 0 + opt90.t0;
        opt90.agesDensity = '72';
        opt90.opt = optDyn;
        
        bw = 800;
        density = 'flexible';
        
    % Baseline
        scalePhi = 1;
    
        [ests13, estsBS, estPre, estPreBS] = ...
            estBorn13(bw, opt80, opt90, density, optOuter, p0);
        
        % format and print results 
        name = 'Born Jan-March';
        print13(name, ests13, estsBS, scalePhi, estPre, estPreBS);
        
    toc;


%% Point estimates, pooling both transitions, static (Table 4)

        scalePhi = 1;

    fprintf('\n\n\n\n\n')
    fprintf('Table 4')
    fprintf('\n\n\n\n\n')


tic;

    % general fields
    t0Base = 0.2721; %average federal MTR + fica, 1989-1990
    scale = 1000;
    
    % Set fields for 1990
    opt90.yMin = 1989;
    opt90.yMax = 1990;
    opt90.t0 = t0Base;
    opt90.t1Pre = opt90.t0 + 0.5;
    opt90.t1Post = opt90.t0 + 0.33;
    opt90.ages = '72';
    opt90.opt = optDyn;
    
    % load the data around age 70
    opt70.ages = 69:71;
    opt70.t0 = t0Base;
    opt70.t1Pre = t0Base + 0.33;
    opt70.t1Post = t0Base;
    opt70.agesD = '72';
    opt70.opt = optDyn;

    % other parameters
    bw = 800;
    p0 = [0.25 0.15];
    method = 'nonparametric';

    % Set up table formatting 
    fprintf('\n\n Point estimates from pooled approach (static) \n'); 
    fprintf('                   &');
    fprintf(' $\\varepsilon$ & $\\phi$ & $\\pi_1$ & $\\pi_2$');
    fprintf(' & $\\pi_3$ & \\pi_4 & \\pi_5 \\\\ \n');
    
    % Baseline: estimate with flexible density
        density = 'flexible';
        name = 'Baseline           ';
        [est, ~, ~, estsBS] = estPooled(bw, opt90, opt70, ...
            density, optOuter, p0, name, scale, method, 0);
        
        % format and print results 
        printPooled(name, est, estsBS, scalePhi);
    
    % estimate with uniform density
        density = 'uniform';
        name = 'Uniform density    ';
        [est, ~, ~, estsBS] = estPooled(bw, opt90, opt70, ...
            density, optOuter, p0, name, scale, method, 0);
        
        % format and print results 
        printPooled(name, est, estsBS, scalePhi);
    
    % Benefit enhancement
        % restore density
        density = 'flexible';        
        name = 'Benefit enhancement';
        
        % Update taxes
        opt90.t0 = t0Base;
        opt90.t1Pre = opt90.t0 + 0.36;
        opt90.t1Post = opt90.t0 + 0.24;
        opt70.t0 = t0Base;
        opt70.t1Pre = t0Base + 0.24;
        opt70.t1Post = t0Base;
        [est, ~, ~, estsBS] = estPooled(bw, opt90, opt70, ...
            density, optOuter, p0, name, scale, method, 0);
        
        % format and print results 
        printPooled(name, est, estsBS, scalePhi);
    
    % Exclude FICA
        name = 'Excluding FICA     ';
        
        % Update taxes
        t0 = 0.1205;
        opt90.t0 = t0;
        opt90.t1Pre = opt90.t0 + 0.5;
        opt90.t1Post = opt90.t0 + 0.33;
        opt70.t0 = t0;
        opt70.t1Pre = t0Base + 0.33;
        opt70.t1Post = t0Base;
        [est, ~, ~, estsBS] = estPooled(bw, opt90, opt70, ...
            density, optOuter, p0, name, scale, method, 0);
        
        % format and print results 
        printPooled(name, est, estsBS, scalePhi);
        
    % BW = 400
        bw = 400;
        name = 'BW = \\$400         ';
        % Change taxes back
        opt90.t0 = t0Base;
        opt90.t1Pre = opt90.t0 + 0.5;
        opt90.t1Post = opt90.t0 + 0.33;
        opt70.t0 = t0Base;
        opt70.t1Pre = t0Base + 0.33;
        opt70.t1Post = t0Base;
        [est, ~, ~, estsBS] = estPooled(bw, opt90, opt70, ...
            density, optOuter, p0, name, scale, method, 0);
        
        % format and print results 
        printPooled(name, est, estsBS, scalePhi);
        
    % BW = 1600
        bw = 1600;
        name = 'BW = \\$1600        ';
        [est, ~, ~, estsBS] = estPooled(bw, opt90, opt70, ...
            density, optOuter, p0, name, scale, method, 0);
        
        % format and print results 
        printPooled(name, est, estsBS, scalePhi);
        
        
toc; 

%% Dynamic estimates (Table 5)
tic;


    fprintf('\n\n\n\n\n')
    fprintf('Table 5')
    fprintf('\n\n\n\n\n')


    % general fields
    scalePhi = 1;
    t0Base = 0.2721; %average federal MTR + fica, 1989-1990
    scale = 1000;
    
    % Set fields for 1990
    opt90.yMin = 1987;
    opt90.yMax = 1992;
    opt90.t0 = t0Base;
    opt90.t1Pre = opt90.t0 + 0.5;
    opt90.t1Post = opt90.t0 + 0.33;
    opt90.ages = '72';
    opt90.opt = optDyn;
    
    % load the data around age 70
    opt70.ages = 67:72;
    opt70.t0 = t0Base;
    opt70.t1Pre = t0Base + 0.33;
    opt70.t1Post = t0Base;
    opt70.agesD = '72';
    opt70.opt = optDyn;

    % other parameters
    bw = 800;
    prs = [0.5 0.1 0.1 0.1 0.1];
    los = log( prs./(1 - prs) ) ;
    p0 = [0.4 0.25 los];
    method = 'nonparametric';

    % Set up table formatting 
    fprintf('\n\n Point estimates from pooled approach \n'); 
    fprintf('                   &');
    fprintf(' $\\varepsilon$ & $\\phi$ & $\\pi_1$ & $\\pi_2$');
    fprintf(' & $\\pi_3$ & \\pi_4 & \\pi_5 \\\\ \n');
    
    % Baseline: estimate with flexible density
        density = 'flexible';
        name = 'Baseline           ';
        [est, ~, ~, estsBS] = estPooled(bw, opt90, opt70, ...
            density, optOuter, p0, name, scale, method, 0);
        
        % format and print results 
        printPooled(name, est, estsBS, scalePhi);
        mainDynamicEsts = est;
        toc;
        
        
    % estimate with uniform density
        density = 'uniform';
        name = 'Uniform density    ';
        [est, ~, ~, estsBS] = estPooled(bw, opt90, opt70, ...
            density, optOuter, p0, name, scale, method, 0);
        
        % format and print results 
        printPooled(name, est, estsBS, scalePhi);
            
    % Benefit enhancement

        prs = [0.9 0.1 0.1 0.1 0.1];
        los = log( prs./(1 - prs) ) ;
        p0 = [0.25 0.15 los];
    
        % restore density
        density = 'flexible';        
        name = 'Benefit enhancement';
        
        % Update taxes
        opt90.t0 = t0Base;
        opt90.t1Pre = opt90.t0 + 0.36;
        opt90.t1Post = opt90.t0 + 0.24;
        opt70.t0 = t0Base;
        opt70.t1Pre = t0Base + 0.24;
        opt70.t1Post = t0Base;
        [est, ~, ~, estsBS] = estPooled(bw, opt90, opt70, ...
            density, optOuter, p0, name, scale, method, 0);
        
        % format and print results 
        printPooled(name, est, estsBS, scalePhi);
    
    % Exclude FICA
        name = 'Excluding FICA     ';
        
        % Update taxes
        t0 = 0.1205;
        opt90.t0 = t0;
        opt90.t1Pre = opt90.t0 + 0.5;
        opt90.t1Post = opt90.t0 + 0.33;
        opt70.t0 = t0;
        opt70.t1Pre = t0Base + 0.33;
        opt70.t1Post = t0Base;
        [est, ~, ~, estsBS] = estPooled(bw, opt90, opt70, ...
            density, optOuter, p0, name, scale, method, 0);
        
        % format and print results 
        printPooled(name, est, estsBS, scalePhi);
        
    % BW = 400
        bw = 400;
        name = 'BW = \\$400         ';
        % Change taxes back
        opt90.t0 = t0Base;
        opt90.t1Pre = opt90.t0 + 0.5;
        opt90.t1Post = opt90.t0 + 0.33;
        opt70.t0 = t0Base;
        opt70.t1Pre = t0Base + 0.33;
        opt70.t1Post = t0Base;
        [est, ~, ~, estsBS] = estPooled(bw, opt90, opt70, ...
            density, optOuter, p0, name, scale, method, 0);
        
        % format and print results 
        printPooled(name, est, estsBS, scalePhi);
        
    % BW = 1600
        bw = 1600;
        name = 'BW = \\$1600        ';
        [est, ~, ~, estsBS] = estPooled(bw, opt90, opt70, ...
            density, optOuter, p0, name, scale, method, 0);
        
        % format and print results 
        printPooled(name, est, estsBS, scalePhi);
toc; 


%% Heterogeneity, pooling both transitions, static (table b.3)
tic;


    fprintf('\n\n\n\n\n')
    fprintf('Table B3')
    fprintf('\n\n\n\n\n')

    % general fields
    t0Base = 0.2721; %average federal MTR + fica, 1989-1990
    scale = 1000;
    scalePhi = 1;
    
    % Set fields for 1990
    opt90.yMin = 1989;
    opt90.yMax = 1990;
    opt90.t0 = t0Base;
    opt90.t1Pre = opt90.t0 + 0.5;
    opt90.t1Post = opt90.t0 + 0.33;
    opt90.ages = '72';
    opt90.opt = optDyn;
    
    % load the data around age 70
    opt70.ages = 69:71;
    opt70.t0 = t0Base;
    opt70.t1Pre = t0Base + 0.33;
    opt70.t1Post = t0Base;
    opt70.agesD = '72';
    opt70.opt = optDyn;

    % other parameters
    bw = 800;
    p0 = [0.25 0.15];
    method = 'nonparametric';
    density = 'flexible';
    
    % Set up table formatting 
    fprintf('\n\n Point estimates from pooled approach (static) \n'); 
    fprintf('                   &');
    fprintf(' $\\varepsilon$ & $\\phi$ & $\\pi_1$ & $\\pi_2$');
    fprintf(' & $\\pi_3$ & \\pi_4 & \\pi_5 \\\\ \n');

    names = {'Men', 'Women', 'High Lifetime Income', 'Low Lifetime Income', ...
        'High Lifetime SD', 'Low Lifetime SD' };
    samples = {'s1Male', 's1Female', 's1HighE', 's1LowE', ...
        's1HighV1', 's1LowV1' };

    for ss = 1:6
        sample = samples{ss};
        name = names{ss};
        [est, ~, ~, estsBS] = estPooledH(bw, opt90, opt70, ...
            density, optOuter, p0, name, scale, method, 0, sample);
        printPooled(name, est, estsBS, scalePhi);
        
        
        % report p-values for difference
        
        if (ss ~= 2*floor(ss/2))
           ests1 = est;
           bs1 = estsBS;
        end
        
        if (ss == 2*floor(ss/2))
            diff = ests1 - est;
            diff(2) = diff(2)*1e3;
            diffBS = bs1 - estsBS;
            diffBS(:,2) =  diffBS(:, 2)*1000;

            lCI = percentile(diffBS, .025);
            uCI = percentile(diffBS, .975);


            nPos = sum(diffBS > zeros(size(estsBS)));
            
            fprintf('\n');
            fprintf(' Difference ');
            fprintf('&          %5.2f         ', diff);
            fprintf(' \\\\ \n');
            fprintf('                     ');
            for ii = 1:length(lCI)
                fprintf('& (%5.2f, %5.2f) ', lCI(ii), uCI(ii));
            end
            fprintf(' \\\\ \n');
            fprintf('N Pos');
            for ii = 1:length(lCI)
                fprintf('& %g ', nPos(ii))
            end
            fprintf(' \\\\ \n \n \n');
        end
        
    end
        
toc; 

%% New method with frictional CF earnings

    fprintf('\n\n\n\n\n')
    fprintf('Table B4')
    fprintf('\n\n\n\n\n')

% TABLE H2
        tic;
        fprintf('           &     (1)        &      (2)       &         (3)     &       (4) \\\\ \\midrule');
        fprintf('           & $\\varepsilon$ &    $\\phi$     &     \\multicolumn{2}{c}{$\\varepsilon |\\phi=0$}    \\\\ \n')
        fprintf(' \\cmidrule(lr){2-2} \n');
        fprintf(' \\cmidrule(lr){3-3} \n');
        fprintf(' \\cmidrule(lr){4-5} \n');
        fprintf('           &                &                &      1990         &      Pre      \\\\ \n')
        fprintf(' \\cmidrule(lr){4-4} \n');
        fprintf(' \\cmidrule(lr){5-5} \n');

        p0 = [0.25 .150];
        nDraws = 50;
        t0Base = 0.2721; %average federal MTR + fica, 1989-1990
        scale = 1000; % so that phi and e are on the same scale
        
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % 66-68 year olds, 1989-1990, s1Main
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        yearsPre = 1;
        t0 = t0Base;
        t1Pre = 0.5 + t0;
        t1Post = 0.33 + t0;
        bw = 800;
        name = 'Baseline';
        density = 'flexible';
        ages = '72';
        estimateCFFric(bw,yearsPre, t0, t1Pre, t1Post, density, optInner, optBounds, optOuter, nDraws, p0, name, ages, scale);

        
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % 66-68 year olds, 1989-1990, s1Main, uniform density 
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        yearsPre = 1;
        t0 = t0Base ;
        t1Pre = 0.5 + t0;
        t1Post = 0.33 + t0;
        bw = 800;
        name = 'Uniform Density';
        density = 'uniform';
        ages = '72';
        nDraws = 100;
        estimateCFFric(bw,yearsPre, t0, t1Pre, t1Post, density, optInner, optBounds, optOuter, nDraws, p0, name, ages, scale);    

        
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % 66-68 year olds, 1989-1990, s1Main, benefit enhancement
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        yearsPre = 1;
        t0 = t0Base ;
        t1Pre = 0.36 + t0;
        t1Post = 0.24 + t0;
        bw = 800;
        ages = ''; 
        density = 'flexible';
        name = 'Benefit Enhancement ';
        ages = '72';
        estimateCFFric(bw,yearsPre, t0, t1Pre, t1Post, density, optInner, optBounds, optOuter, nDraws, p0, name, ages, scale);
    
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Exclude FICA
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        yearsPre = 1;    
        t0 = 0.1205;
        t1Pre = 0.5 + t0;
        t1Post = 0.33 + t0;
        bw = 800;
        ages = ''; 
        density = 'flexible';
        name = 'Excluding FICA';
        ages = '72';
        estimateCFFric(bw,yearsPre, t0, t1Pre, t1Post, density, optInner, optBounds, optOuter, nDraws, p0, name, ages, scale);
        
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % BW = 400;
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        yearsPre = 1;    
        t0 = t0Base;
        t1Pre = 0.5 + t0;
        t1Post = 0.33 + t0;
        bw = 500;
        ages = ''; 
        density = 'flexible';
        name = 'Bandwidth = \$400';
        ages = '72';
        estimateCFFric(bw,yearsPre, t0, t1Pre, t1Post, density, optInner, optBounds, optOuter, nDraws, p0, name, ages, scale); 
 toc;

diary off

