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

function [dataB bVarB bunchBSB dataPre bVarPre bunchBSPre dataPost bVarPost bunchBSPost] = ...
    setData1990(bw, yearsPre,t0,t1Pre,t1Post,density, ages, scale)
    

    % Load the data
    rawBS = importdata(sprintf('yearBS_7802_6668_%g_s1Main.txt',bw));
    rawBunch = importdata(sprintf('yearStatsDensPoly_7802_6668_%g_s1Main.txt',bw));
    rawZStar = importdata('realZStar.csv');
    
    bunchYear = rawBunch.data(:,1);
    assert (isequal(rawBS.data(:,1), bunchYear));
    
    % align bunch and bunch bs matrices
    bbIndex = logical((bunchYear>=1985).*(bunchYear<=1994));
    bunch = bw/scale.*rawBunch.data(bbIndex==1,2);
    year = bunchYear(bbIndex==1);
    bunchVar = (rawBunch.data(bbIndex==1,3));
    
    bunchBS = bw/scale .* (rawBS.data(bbIndex,2:end));
    assert(isequal(size(bunchBS,1), size(bunch,1)));
    
    zStarPost = rawZStar.data(:,3)/scale;
    zStarYear = rawZStar.data(:,1);

    % Set up the data for estimating both, using pre and post data
    dataB.t1 = [t1Pre.*ones(yearsPre,1) ; t1Post];
    dataB.t0 = t0.*ones(yearsPre+1,1);
    dataB.pre = [ones(yearsPre,1) ; 0];
    dataB.last = logical([zeros(yearsPre-1, 1) ; 1 ; 0]);
    dataB.zStar = zStarPost( (zStarYear<=1990).*(zStarYear>=1990-yearsPre)==1 );
    
    dataB.mu = 10.12; % from tobitClaim log
    dataB.sigma = 1.266; 
    
    sample = ( (year<=1990).* (year>=1990-yearsPre));
    dataB.bunch = bunch(sample==1);
    dataB.rk = zeros(yearsPre+1,1);
    dataB.density = density;
    
    bVarB = bunchVar(sample==1);
    bunchBSB = bunchBS(sample==1,:);
        
    % Strucutre available in both periods with pre-period tax info
    dataB.tPre.t1 = t1Pre;
    dataB.tPre.t0 = t0;
    dataB.tPre.zStar = dataB.zStar(1:end-1);
    
    dataB.tPost.t1 = t1Post;
    dataB.tPost.t0 = t0;
    dataB.tPost.zStar = dataB.zStar(end);
        
    % Calculate density, if called for 
    if strcmp(density, 'flexible')
        dataB.scale = 1000/scale;
        dataB.bw = bw/scale;
        densityData = importdata(sprintf('dens%s_8399_%g_s1Main.csv', ages, bw));
        prob = densityData.data(:,1);
        dataB.dens = prob ./ (bw/scale);
        dataB.dist = cumsum(prob);   
        dataB.income = dataB.scale.*densityData.data(:,2) + dataB.zStar(1);
        zStarPoint = logical(densityData.data(:,2) == 0);
        dataB.densAK = prob(zStarPoint)/dataB.bw;
    end    
   
    % Data for pre period
    dataPre = dataB;
    dataPre.t1 = t1Pre.*ones(yearsPre,1) ;
    dataPre.t0 =  t0.*ones(yearsPre,1);
    dataPre.pre = ones(yearsPre,1);
    dataPre.zStar = zStarPost( (zStarYear<1990).*(zStarYear>=1990-yearsPre)==1 );
    dataPre.mu = 10.12; % from tobitClaim log
    dataPre.sigma = 1.266; 
    
    sample = ( (year<1990).* (year>=1990-yearsPre));
    dataPre.bunch = bunch(sample==1);
    dataPre.rk = zeros(sum(sample),1);
    dataPre.density = density;
    
    bVarPre = bunchVar(sample==1);
    bunchBSPre = bunchBS(sample==1,:);
    
    % Data for post period
    dataPost = dataB;
    dataPost.t1 = t1Post;
    dataPost.t0 =  t0;
    dataPost.pre = 1; % set to one because we set phi = 0 when using this sample 
    dataPost.zStar = zStarPost(zStarYear==1990);
    
    dataPost.mu = 10.12; % from tobitClaim log
    dataPost.sigma = 1.266; 
    
    sample = (year==1990);
    dataPost.bunch = bunch(sample==1);
    dataPost.rk = 0;
    dataPost.density = density;
    
    bVarPost = bunchVar(sample==1);
    bunchBSPost = bunchBS(sample==1,:);
    
end
