%% setData1990
%
% Set up data structures to evaluate the 1989-1990 change, assuming
% pre-period bunching is attenuated by frictions

function [data, bunchBSB, dataPre, bunchBSPre, dataPost, bunchBSPost] = setDataCFFric(bw, yearsPre, t0,t1Pre,t1Post, density, ages, nDraws, scale)
    
    % Load the data
    rawBS = importdata(sprintf('yearBS_8594_6768_%g_s1Main.txt',bw));
    rawBunch = importdata(sprintf('yearStatsDensPoly_8594_6768_%g_s1Main.txt',bw));
    rawZStar = importdata('realZStar.csv');
    
    bunchYear = rawBunch.data(:,1);
    
    % align bunch and bunch bs matrice
    bbIndex =logical( (bunchYear>=1990-yearsPre).*(bunchYear<=1990) );
    bunch = bw.*rawBunch.data(bbIndex,2)/scale;
    bsYear = rawBS.data(:,1);
    bsIndex = logical( (bsYear>=1990-yearsPre).*(bsYear<=1990) );
    bunchBSB = bw.*(rawBS.data(bsIndex,2:end))/scale;
    assert(isequal(size(bunchBSB,1), size(bunch,1)));
    
    zStar = rawZStar.data(:,3)/scale;
    zStarYear = rawZStar.data(:,1);
    
    % Calculate density, if called for 
    if strcmp(density, 'flexible')
        data.scale = 1000/scale;
        data.bw = bw/scale;
        densityData = importdata(sprintf('dens%s_8389_%g_s1Main.csv', ages, bw));
        prob = densityData.data(:,1);
        data.dist = cumsum(prob);       
        data.income = data.scale.*densityData.data(:,2);
        zStarPoint = logical(abs(data.income-zStar(zStarYear==1989))<data.bw/2);
        data.densAK = prob(zStarPoint)/data.bw;
        data.zDraws = rand(nDraws,1);
    end    
    
    % Set up the data for estimating both, using pre and post data
    data.tax0.zStar = 1e4;
    data.tax0.t0 = t0;
    data.tax0.t1 = t0;
    
    data.tax1.zStar = zStar(logical( (zStarYear<1990).*(zStarYear>=1990-yearsPre) ) );
    data.tax1.t0 = t0.*ones(yearsPre,1);
    data.tax1.t1 = t1Pre.*ones(yearsPre,1);
    
    data.tax2.zStar = zStar(zStarYear==1990);
    data.tax2.t0 = t0;
    data.tax2.t1 = t1Post;
    
    data.bunch = bunch;
    
    % Data for pre period
    dataPre = data;
    dataPre.t1 = t1Pre.*ones(yearsPre,1) ;
    dataPre.t0 =  t0.*ones(yearsPre,1);
    dataPre.pre = ones(yearsPre,1);
    dataPre.zStar = zStar( (zStarYear<1990).*(zStarYear>=1990-yearsPre)==1 );
    dataPre.mu = 10.12; % from tobitClaim log
    dataPre.sigma = 1.266; 
    
    dataPre.bunch = bunch(1:end-1);
    dataPre.rk = zeros(size(dataPre.bunch));
    dataPre.density = density;
    bunchBSPre = bunchBSB(1:end-1,:);
    
    % Data for post period
    dataPost = data;
    dataPost.t1 = t1Post;
    dataPost.t0 =  t0;
    dataPost.pre = 1; % set to one because we set phi = 0 when using this sample 
    dataPost.zStar = zStar(zStarYear==1990);
    
    dataPost.mu = 10.12; % from tobitClaim log
    dataPost.sigma = 1.266; 
    
    dataPost.bunch = bunch(end);
    dataPost.rk = 0;
    dataPost.density = density;
    
    bunchBSPost = bunchBSB(end,:);
    
end