function [mutildeH,sigmaPerm,deltaBoth,eta,Nobs,growthall]=panelmuH(Year1,IncPctileCutoff,MuTildePctile,DeltaDrop,deltaCD,pbar,PermTrans,IncomeType);

%function [mutildeH,sigmaPerm,deltaBoth,eta,Nobs,growthall]=panelmuH(Year1,IncPctileCutoff,MuTildePctile,DeltaDrop,deltaCD,pbar,PermTrans,IncomeType);
%
%  Estimate muH, sigmaH, and delta in the IRS public use panel data using the
%  cross-section of growth rates of incomes
%    Basic method: throw away the bottom growth rates (DeltaDrop) and use the 
%        upper part of the distribution to recover muH and sigmaH
%    
%    Year1 = first year (second year is Year1+1)    e.g. 1989 for 1989-1990
%    IncPctileCutoff = Look only at initial incomes above this percentile
%    MuTildePctile   = Percentile of the growth rate distribution to estimate mutilde
%    DeltaDrop = 40 ==> drop growth rates that lead income to fall by more than 40%
%    PermTrans = 1/2 = Fraction of var(growth) that is the perm random walk
%                      -- use this fraction in computing eta
% IncomeType: 0=Wage+Salaries, 1=Wages/Salaries+EntrepreneurialIncome,  2=Entrepreneurial only
%
%  09/23/16: Update includebusiness to Piketty-Saez "entrepreneurial income"
%       “data11+data17+data21+data75+data79” = Wages/Salaries + Entrepreneurial Income


if ~exist('IncomeType'); IncomeType=0; end;
years=[Year1 Year1+1]
cpiSaez % loads the CPI data and SaezMeanIncome for normalization

wagevar=11;  % Wage and salaries
entrevar1=17;%   'data17 business net income (schedule C)'	  
entrevar2=21;%   'data21 farm income'                             
entrevar3=75;%   'data75 net partnership income or loss '         
entrevar4=79;%   'data79 small business corp total income or loss'
entrevars=[entrevar1 entrevar2 entrevar3 entrevar4];
idnum=104  % U of M panel file ID number

NumtoShow=100
step=1
recs=(1:step:step*NumtoShow)';


% ----------------
%   First Year
% ----------------
disp ' '; 
fname=['p' num2str(years(1))];
eval(['load /home/data/soi/taxsim/txt/' fname '.txt']);
eval(['data1=' fname '; clear ' fname]); 

% Data before cleaning
entreincome1=sum(data1(:,entrevars)')';
fmt='%4.0f %8.0f %8.0f %15.0f %15.0f %10.0f';
cshow(' ',[recs data1(recs,[1 2 wagevar]) entreincome1(recs)  data1(recs,idnum)],fmt,'Line Weight Married Wage Business PanelID');


% Clean out if not married filing jointly
married=data1(:,2); % married=2 for married-joint filing

disp ' ';
if IncomeType==0;
    disp 'Income Measure = Wages and Salaries';
    income1=data1(married==2,wagevar);
elseif IncomeType==1;
    disp 'Income Measure = Wages/Salaries + Piketty-Saez Entrepreneurial Income';
    disp '      data11+data17+data21+data75+data79 = Wages/Salaries + Entrepreneurial Income'
    disp 'data17 business net income (schedule C)'	
    disp 'data21 farm income'
    disp 'data75 net partnership income or loss '
    disp 'data79 small business corp total income or loss'
    income1=data1(married==2,wagevar)+entreincome1(married==2);
elseif IncomeType==2;
    disp 'Income Measure = Piketty-Saez Entrepreneurial Income';
    disp '      data17+data21+data75+data79 = Entrepreneurial Income'
    disp 'data17 business net income (schedule C)'	
    disp 'data21 farm income'
    disp 'data75 net partnership income or loss '
    disp 'data79 small business corp total income or loss'
    income1=entreincome1(married==2);
end;

id1=data1(married==2,idnum);
weight1=data1(married==2,1); % sample weight -- this many "people" with this observation
fprintf('There are %12.0f kept observations in total\n',length(income1));



% ----------------
%   Second Year
% ----------------
disp ' '; 
fname=['p' num2str(years(2))];
eval(['load /home/data/soi/taxsim/txt/' fname '.txt']);
eval(['data2=' fname '; clear ' fname]); 

% Data before cleaning
entreincome2=sum(data2(:,entrevars)')';
cshow(' ',[recs data2(recs,[1 2 wagevar]) entreincome2(recs)  data2(recs,idnum)],fmt,'Line Weight Married Wage Business PanelID');

% Clean out if not married filing jointly
married=data2(:,2); % married=2 for married-joint filing

disp ' ';
if IncomeType==0;
    disp 'Income Measure = Wages and Salaries';
    income2=data2(married==2,wagevar);
elseif IncomeType==1;
    disp 'Income Measure = Wages/Salaries + Piketty-Saez Entrepreneurial Income';
    disp '      data11+data17+data21+data75+data79 = Wages/Salaries + Entrepreneurial Income'
    income2=data2(married==2,wagevar)+entreincome2(married==2);
elseif IncomeType==2;
    disp 'Income Measure = Piketty-Saez Entrepreneurial Income';
    disp '      data17+data21+data75+data79 = Entrepreneurial Income'
    income2=entreincome2(married==2);
end;

id2=data2(married==2,idnum);
weight2=data2(married==2,1); % sample weight -- this many "people" with this observation
fprintf('There are %12.0f kept observations in total\n',length(income2));


% CPI correction
income1=income1/cpi(years(1)-cpiyrs(1)+1)*100;  % cpiyears(1)=1913   2012 dollars
income2=income2/cpi(years(2)-cpiyrs(1)+1)*100;  % cpiyears(1)=1913  2012 dollars

% Normalize by aggregate -- a bit noisy because of small sample sizes...
%income1=income1/mean(income1);
%income2=income2/mean(income2);

% Normalize by aggregate SaezMeanIncome -- AGI excluding capital gains
income1=income1/SaezMeanIncome(years(1)-cpiyrs(1)+1);  % cpiyears(1)=1913
income2=income2/SaezMeanIncome(years(2)-cpiyrs(1)+1); 

% Replace zeros with $1 for taking logs. We do not use those obs except to count
income1(income1==0)=1; 
income2(income2==0)=1; 


% ------------------------------
% Sort data in Year 1
% ------------------------------

% Now sort the data from first year
[ysort1,indx]=sort(income1);
weight1=weight1(indx); % sorted weight
id1=id1(indx);
N1=length(ysort1);

% Find the IncPctileCutoff percentile
if length(IncPctileCutoff)==1;
    nstart=floor(IncPctileCutoff*N1/100);
    disp ' ';
    fprintf('The %5.2f percentile (normalized) income is %8.3f\n',[IncPctileCutoff ysort1(nstart)]);
    fprintf('and there are %7.0f observations above that\n',N1-nstart);
    nend=N1;
else; % IncPctileCutoff=[90 95]
    nstart=floor(IncPctileCutoff(1)*N1/100);
    nend=ceil(IncPctileCutoff(2)*N1/100);
    disp ' ';
    fprintf('The %5.2f percentile (normalized) income is %8.3f\n',[IncPctileCutoff(1) ysort1(nstart)]);
    fprintf('The %5.2f percentile (normalized) income is %8.3f\n',[IncPctileCutoff(2) ysort1(nend)]);
    fprintf('and this range includes %7.0f observations\n',nend-nstart);
    
end;
% -------------------------------
% Match and compute income growth
% -------------------------------

disp ' ';
fmt='%7.0f %7.0f %8.0f %8.0f %15.3f %15.3f %8.2f';
tle='ID1 ID2 Weight1 Weight2 Income1 Income2 Growth';
cshow(' ',[],fmt,tle);
growth=zeros(N1,1)*NaN;
for i=nstart:nend;
    indx=find(id2==id1(i));
    indx1=find(id1==id1(i)); % Check for unique ID in both years!
    if ~isempty(indx) & length(indx)==1 & length(indx1)==1;
        if income2(indx)>0;
            growth(i)=100*log(income2(indx)/ysort1(i));
        else; % Percent change if income2 is negative to help with DeltaDrop
            growth(i)=100*(income2(indx)/ysort1(i)-1);  
        end;
        if i>nend-NumtoShow;
            cshow(' ',[id1(i) id2(indx) weight1(i) weight2(indx) ysort1(i) income2(indx) growth(i)],fmt);
        end;
    end;
end;

% Estimating delta based on fraction of large negative growth rates
keepit=~isnan(growth);
growthall=growth(keepit)/100;
tooneg=(growthall<log(1-DeltaDrop/100)); % i.e. if income falls by more than 60% ==> DeltaDrop=60
deltaBoth=sum(tooneg)/(sum(keepit)); % Our point estimate for delta+deltabar
Nobs=length(growthall);

% Note well: the panel sample has no sampling weights (the weight variable
% is the same for all observations, so no need to weight. (This contrasts
% with the cross-section data).

g=growthall(~tooneg);
disp ' ';
fprintf('Dropping %4.0f observations where incomes fall by more than %4.1f percent\n',[sum(tooneg) DeltaDrop]);
disp ' ';


mutildeH=getprctile(g,MuTildePctile);
GrowthPercentile=100-2*(100-MuTildePctile); % E.g. mutilde is 95th, then use 90th and above for std
i=floor(GrowthPercentile/100*length(g)); % All entries greater than this.
sigmaH=std(g(i:end));
sigmaPerm=sqrt(PermTrans)*sigmaH;



% Compute eta
deltabar=deltaBoth-deltaCD;
if deltabar<0.01; 
    disp 'PROBLEM: deltabar<0.01 in panelmuH.m. Setting deltabar=.01...';
    deltabar=.01
end;
mm=mutildeH./sigmaPerm.^2;
squig=-mm+sqrt(mm.^2+2*(deltabar+pbar)./sigmaPerm.^2);
eta=1./squig;

tle='  %4.0f-%4.0f  NObs=%5.0f   mutildeH=%5.2f   StdevH=%5.2f  deltaBoth=%6.4f  eta=%8.4f\n';
fprintf(tle,[years Nobs mutildeH sigmaH deltaBoth eta]);
disp ' ';
