%% setData1990
%
% Set up data structures to evaluate the 1989-1990 change 

function [dataB, bunchBSB, dataPre, bunchBSPre] = setData6970(bw, ages,t0,t1Pre,t1Post,density, agesDensity, scale)
    
    % Load the data
    rawBS = importdata(sprintf('ageBS_9099_6773_%g_s1Main.txt',bw)); 
    rawBunch = importdata(sprintf('ageStatsDensPoly_9099_5973_%g_s1Main.txt',bw));    
    rawZStar = importdata('realZStar.csv');
    
    % set up main bunching data
    bunchAge = rawBunch.data(:,1);
    bbIndex = ismember(bunchAge, ages);
    bunch = bw/scale.*rawBunch.data(bbIndex==1,2);
    age = bunchAge(bbIndex==1);
    
    % set up bunchBS matrix
    bsAge = rawBS.data(:,1);
    bsIndex = ismember(bsAge, ages);
    bunchBS = bw/scale.*(rawBS.data(bsIndex,2:end));
    assert(isequal(size(bunchBS,1), size(bunch,1)));    
    
    % Set up the data for estimating both, using pre and post data
    dataB.t1 = t1Pre.*(age<70) + t1Post.*(age>=70);
    dataB.t0 = t0;
    dataB.pre = age<70;
    zStarYear = rawZStar.data(:,1);
    dataB.zStar = mean(rawZStar.data( ((zStarYear>=1990) .* (zStarYear<=1999)==1),3))/scale;
    
    dataB.mu = 10.12; % from tobitClaim log
    dataB.sigma = 1.266; 
    
    dataB.bunch = bunch;
    dataB.rk = 0;
    dataB.density = density;
    bunchBSB = bunchBS;    
    
    % Load flexible distribution, if called for
    dataB.density = density;
    dataB.scale = 1000/scale;
    dataB.bw = bw/scale;
    densityData = importdata('dens72_8399_800_s1Main.csv');
    dataB.income = dataB.scale.*densityData.data(:,2) + dataB.zStar;
    if strcmp(density, 'flexible')
        densityData = importdata(sprintf('dens%s_8399_%g_s1Main.csv', agesDensity, bw));
        dataB.income = dataB.scale.*densityData.data(:,2) + dataB.zStar;
        prob = densityData.data(:,1);
        dataB.dist = cumsum(prob);   
        dataB.dens = prob ./ (bw/scale);
        zStarPoint = logical(densityData.data(:,2) == 0);
        dataB.densAK = prob(zStarPoint)/dataB.bw;
    end    
     
    dataB.taxPre.t0 = t0;
    dataB.taxPre.t1 = t1Pre;
    dataB.taxPre.zStar = dataB.zStar;
    
    dataB.taxPost.t0 = t0;
    dataB.taxPost.t1 = t1Post;
    dataB.taxPost.zStar = dataB.zStar;
    
    % Data for pre period
    dataPre = dataB;
    dataPre.pre = 1;
    dataPre.bunch = bunch(age<70);    
    dataPre.t1 = t1Pre;
    bunchBSPre = bunchBS(age<70,:);
    
end