function [ylife,ylook,ageyn2,earnginiage,ageyn_dye]= stats_earnings_process_dye(ylife,Qylife,invherhip,normtype,ye,yedtb)

global dy dye dbet dhip TL gpop eff sur pfrac


% Much of this is not relevant for the loop - the normalization is
% important though and happens here (need ageyn for it)

% form the transition matrix with the blocks for beta on the main diag.:
%                                    age 1        age 2    .... age TL
%                              prod1....prod dy  .....         ........
% age 1;HIP1 productivity level 1     
%        productivity level 2
%        ....................
%        productivity level dy
%       HIP2...
%       dhip
% age 2; productivity level 1
%        productivity level 2
%        ....................
%        productivity level dy
% ...........................
% age TL productivity level 1
%        ....................
%        productivity level dy

EM=zeros(dbet*TL*dhip*dy,dbet*TL*dhip*dy);
clear j k
for j=1:dbet
    for k=1:TL-1
        for j2=1:dhip
        EM((j-1)*TL*dhip*dy+(k-1)*dhip*dy+(j2-1)*dy+1:(j-1)*TL*dhip*dy+...
            (k-1)*dhip*dy+j2*dy,(j-1)*TL*dhip*dy+k*dhip*dy+(j2-1)*dy+1:...
            (j-1)*TL*dhip*dy+k*dhip*dy+j2*dy)=sur(k).*Qylife(:,:,k);
        end
    end
end % for j=1:dbet; 

% form the vector for the newborns
eq=zeros(dbet*TL*dhip*dy,1);
tempvec=kron((0:(dbet-1))'*TL*dhip*dy,ones(dhip*dy,1))+kron(ones(dbet,1),...
    (0:(dhip*dy-1))')+1;
eq(tempvec)=invherhip;         

if isempty(gpop)
    gpop=1;
end

% now compute the invariant distribution
epsi=1;
eminv=eq;
while epsi>1e-08
   eminv1=eminv;
   eminv=EM'*eminv/gpop+eq;
   epsi=max(abs(eminv1-eminv));
end % for while epsi>1e-08

% renormalize the invariant distribution
eminv=eminv./sum(eminv);


eminvt=reshape(eminv,size(eminv,1)/dbet,dbet); % put betas into columns
if dbet>1 % if there are heterogenous betas, aggregate people over beta
   eminvt=(sum(eminvt,2));
end

% form a matrix with age on the rows and productivity on the columns
% containing the number of people in each cell in order to compute
% the age-efficiency profile IMPLIED by our income process
ageyn2=(reshape(eminvt,dhip*dy,size(eminvt,1)/(dhip*dy)));
% transform the matrix in terms of fractions of people of a given age
ageyn=(ageyn2./(ones(dhip*dy,1)*sum(ageyn2)))';

ageyn_dye=kron(ageyn,ones(1,dye)).*kron(ones(1,dy),yedtb');


%% Adjustment for age-efficiency profile
% Now we do it for all of the options!

% readjust the income process so that the implied age-efficiency profile
% is the one given in the eflife matrix. Be careful: we need to use ageyn2
% for this. First make every age mean 1:

ylook=reshape(ylife,dy,TL);
ylook_l=kron(ylook,ones(dye,1)).*kron(ones(dy,1),ye);
meanage=sum(ylook_l.*ageyn_dye',1);
%menage_simple=sum(ylook.*ageyn',1);

%Perform normalization only in persistent component, but taking
%transitory component into account
ylook=ylook./repmat(meanage,dy,1);
%ylook_noage=ylook;
ylife=reshape(ylook,dy*TL,1);
ylife=kron(eff,ones(dhip*dy,1)).*ylife; %Now input eff
%Now normalize mean to 1 if desired
if normtype==1
ylife=ylife./(sum(pfrac(1:TL)));
end
ylook=reshape(ylife,dy,TL);

% to compute the adjusted gini index, now construct y*eff for each possible
% y and eff and sort them in ascending order, keeping track of the 
% corresponding number of people in that cell. Note y includes both
% hip+stochastic income process.

aearn=ylife;
%Find total earnings for
    %each point - note that dy is the most internal variable and dhip is
    %the previous one, then we change by age. Betas are already summed
    %inside here.

[aearn3,i]=sort(aearn); %we need aearn later
aearn2=[0; aearn3];
adist2=eminvt(i);
adist=[0; adist2];

cumadist=cumsum(adist); % this is the cumulative number of agents
cumaearn=cumsum(aearn2.*adist)./sum(aearn2.*adist); % this is cumulative adjusted productivity

%plot(cumadist,cumaearn);
%axis([0 1 0 1]);
%title('effective productivity Lorenz Curve');
%xlabel('cumulative fraction of agents');
%ylabel('cumulative fraction of effective productivity');
% print -append -dps unobr.ps

% compute the gini index
dgini=size(cumaearn,1);
aarea=0;
for i=1:dgini-1
   aarea=aarea+(cumaearn(i+1)+cumaearn(i))*(cumadist(i+1)-cumadist(i));
end
agini=1-aarea;
fprintf('Agini is  %.4f\n',agini)

%------- now compute the Gini index per age group ----------%
ginidist=zeros(dy+1,TL);
giniy=zeros(dy+1,TL);
earnginiage=zeros(TL,1);

for i=1:TL
ginidist(:,i)=[0; eminvt((i-1)*dhip*dy+1:i*dy*dhip)];
giniy(:,i)=[0; ylook(:,i)];
dyy=dy+1;
%Now note that these matrices are not necessarily in order - change it:
[giniy(:,i),KK]=sort(giniy(:,i));
ginidist(:,i)=ginidist(KK,i);

% lorenz curve plot (cumulative productivity plotted against cumulative
% number of agents; both normalized to one)
cumdist1=cumsum(ginidist(:,i))/sum(ginidist(:,i)); % this is the cumulative number of agents
cumearn1=cumsum(giniy(:,i).*ginidist(:,i))/sum(giniy(:,i).*ginidist(:,i)); % this is cumulative productivity

% plot(cumdist1,cumearn1);
% axis([0 1 0 1]);
% title('productivity Lorenz Curve for new agents');
% xlabel('cumulative fraction of agents');
% ylabel('cumulative fraction of productivity');
% print -append -dps unobr.ps

% plot(cumdist2,cumearn2);
% axis([0 1 0 1]);
% title('productivity Lorenz Curve for TL years old agents');
% xlabel('cumulative fraction of agents');
% ylabel('cumulative fraction of productivity');
% print -append -dps unobr.ps

% compute the gini index for new agents
area=0;
for j=1:dyy-1
   area=area+(cumearn1(j+1)+cumearn1(j))*(cumdist1(j+1)-cumdist1(j));
end
earnginiage(i)=1-area;

end
fprintf('Gini for new agents is %.4f \n',earnginiage(1));
fprintf('Gini for TL-old agents is %.4f \n',earnginiage(TL));
end