function [wgini,ginage,dist,cons_aworkers,cons_aold,lcons_aworkers,lcons_aold,...
    ltransiti,lincwork,lvarlabincwork,lthirdincwork,lfourthincwork,...
    lvarcons,percentable,...
    cons_perctiles,lvarlabincgroup,lvarconsgroup,negwealth,varage,savbyage]=...
    compute_cw_statistics(ylook,u,invher,copt,colr,sopt,invm,invmt,dimTL,dimTR)

global da dy dbet dhip TL T TR a agetime ageoldyoung pfrac

ndim1=length(dimTL);
ndim2=length(dimTR);
        
s1 = prod(dimTL);
s2 = prod(dimTR);

if ndim1==5
   ye=u.ye;
   ylife=u.ylife;
else
    ylife=u;
end


%% 1. Statistics on wealth inequality

% We want to obtain a matrix with assets on rows and the rest of dimensions
% on the columns

invdist=reshape(invm,da,dbet*(s1+s2)/da);
dist=sum(invdist,2);

% Now aggregate over ages

agew=zeros(da,T);

% construct the part for the retired agents
agew(:,TL+1:T)=permute(sum(reshape(invmt(s1+1:s1+s2),da,s2/(da*TR),TR),2),[1 3 2]);

% construct the part for the working agents
temp1=reshape(invmt(1:s1),da,s1/(da*TL),TL);
temp1=permute(temp1,[1 3 2]);
temp1=reshape(temp1,da*TL,s1/(da*TL));
temp1=sum(temp1,2)';
agew(:,1:TL)=reshape(temp1,da,TL);

azero=find(a==0,1,'first');
dist(azero)=sum(agew(azero,2:T)); %eliminate the people just entering the labour
                          %market from the fraction with zero wealth bec.
                          %their wealth is zero by construction


%%%%%%%%%%%%% lorenz curve plot (cumulative wealth plotted against
% cumulative number of agents; both normalized to one)
cumdist=cumsum(dist); % this is the cumulative number of agents
cumwealth=cumsum(a.*dist)/sum(a.*dist); % this is cumulative wealth
% plot(cumdist,cumwealth);
% axis([0 1 0 1]);
% title('Wealth Lorenz Curve');
% xlabel('cumulative fraction of agents');
% ylabel('cumulative fraction of wealth');
% print -append -dps hugg4y.ps


%%%%%%%%%%%%%% gini index
% the area under the lorenz curve is composed by trapeziums:
%          |c
% b |      |
%   |      |
% a |      |d
% whose area is (cd+ba)*ad/2.
% the formula for the gini index given by the ratio of the
% area between the 45 degree line and the the lorenz curve divided by
% the total area under the 45 degree line: 
% (1/2 - \sum_{i=1}^{da-1}(w(i+1)+w(i))*(cumdist(i+1)-cumdist(i))/2)/(1/2)=
% 1-\sum_{i=1}^{da-1}(w(i+1)+w(i))*(cumdist(i+1)-cumdist(i))
area=0;
for i=1:da-1
   area=area+(cumwealth(i+1)+cumwealth(i))*(cumdist(i+1)-cumdist(i));
end
wgini=1-area;

%%%%%%%%%%%%%% top percentiles for wealth for the population
topwel=[cumsum(flipud(a.*dist/sum(a.*dist))) cumsum(flipud(dist))];

%%%%%%%%%%%%%% bottom percentiles for wealth for the population
botwel=[cumsum(a.*dist/sum(a.*dist)) cumsum(dist)];

%%%%%%%%%%%%%% wealth profiles within age groups
%%%implied gini for newborns


%%%%%%%%%%%%%% wealth gini within age groups 
cumaw=cumsum(agew(:,2:T))./(ones(da,1)*sum(agew(:,2:T)));
cumaw0=cumsum(agew(:,1))./(ones(da,1)*sum(agew(:,1)));
% this is the cumulative number of agents for each age  normalized to 1
cumawea=cumsum(agew(:,2:T).*(a*ones(1,T-1)))./...
        (ones(da,1)*sum(agew(:,2:T).*(a*ones(1,T-1))));
    cumawea0=cumsum(agew(:,1).*(a))./...
        (ones(da,1)*sum(agew(:,1).*(a)));
    
    piece1=cumsum(agew(:,2:T).*(a*ones(1,T-1)));
    
% this is cumulative wealth
%%%%%%%%%%%%% lorenz curve plot within age groups
%for j=1:T-1,
%   plot(cumaw(:,j),cumawea(:,j));
%   axis([0 1 0 1]);
%   title(['Wealth Lorenz Curve for age ' num2str(j+1)]);
%   xlabel('cumulative fraction of agents');
%   ylabel('cumulative fraction of wealth');
%   print -append -dps hugg4y.ps
% end
%%%%%%%%%%%%%% gini index within age groups
area=zeros(1,T-1);
for i=1:da-1
   area(1,:)=area(1,:)+...
             (cumawea(i+1,:)+cumawea(i,:)).*(cumaw(i+1,:)-cumaw(i,:));
end
area0=0;
for i=1:da-1
    area0=area0+(cumawea0(i+1)+cumawea0(i)).*(cumaw0(i+1)-cumaw0(i));
end
ginage0=1-area0;
ginage=ones(T-1,1)-area';

%print -append -dps hugg4y.ps

%%%%%%%%%%%%% quantiles by age
% define the quantiles for which we want to trace the age behavior
alph=[.1; .3; .5; .8; .9];
qalph=zeros(T-1,length(alph));
for j=1:length(alph)
   for i=1:T-1
      f=cumaw(:,i);
      i1=find(f>alph(j),1)-1;
      if i1==0
         qalph(i,j)=a(i1+1);
      else
         i2=i1+1;
         m=(f(i2)-f(i1))/(a(i2)-a(i1));
         n=f(i1)-m*a(i1);
         qalph(i,j)=(alph(j)-n)/m;
      end % for if
   end
end
meana=a'*(agew./(ones(da,1)*sum(agew)));


% Percentage of people who have negative wealth

negwealth=sum(dist(1:azero))./sum(dist);

disp('People with non-positive wealth')
disp(sum(dist(1:azero))./sum(dist))

%% 2. Finding var of consumption statistics

% There are several important reshapes of invm, income and consumption that
% we need

% First, one element per dimension

% Order: da x everything else x age x dbet

copt=permute(copt,[1:(ndim1-1) ndim1+1 ndim1]);
colr=permute(colr,[1:(ndim2-1) ndim2+1 ndim2]);

if ndim1==4
    inv5d_workers=reshape(invm(1:s1),dimTL(1),dimTL(2),dimTL(3),dimTL(4));
elseif ndim1==5
    inv5d_workers=reshape(invm(1:s1),dimTL(1),dimTL(2),dimTL(3),dimTL(4),dimTL(5));
end

if ndim2==3
    inv5d_old=reshape(invm(s1+1:s1+s2),dimTR(1),dimTR(2),dimTR(3));
elseif ndim2==4
    inv5d_old=reshape(invm(s1+1:s1+s2),dimTR(1),dimTR(2),dimTR(3),dimTR(4));
end

if ndim1==4
   inc5d_workers=permute(repmat(reshape(ylife,dy,dhip,TL),1,1,1,da),[ndim1 1:(ndim1-1)]);
elseif ndim1==5
   inc5d_workers=permute(repmat(reshape(ylife,dy,dhip,TL),1,1,1,dimTL(2),da),[ndim1 ndim1-1 1:(ndim1-2)]).*...
   reshape(kron(ones(1,dy),kron(ye,ones(da,1))),dimTL(1),dimTL(2),dimTL(3),dimTL(4),dimTL(5));
end


% Second, age on columns, everything else on rows
% Age is always the last dimension of all of these matrices

invlong_workers=reshape(permute(inv5d_workers,[ndim1 1:(ndim1-1)]),TL,s1/TL);
conslong_workers=reshape(permute(copt,[ndim1 1:(ndim1-1)]),TL,s1/TL);
inclong_workers=reshape(permute(inc5d_workers,[ndim1 1:(ndim1-1)]),TL,s1/TL);

invlong_old=reshape(permute(inv5d_old,[ndim2 1:(ndim2-1)]),TR,s2/TR);
conslong_old=reshape(permute(colr,[ndim2 1:(ndim2-1)]),TR,s2/TR);

%% 2.1. Variances and other moments of labour income

lvarlabincwork=zeros(TL,1);
lthirdincwork=lvarlabincwork;
lfourthincwork=lvarlabincwork;
lvarconswork=lvarlabincwork;
lincwork=lvarlabincwork;
cons_aworkers=lvarlabincwork;
lcons_aworkers=lvarlabincwork;

% Pure labor income variances from the process itself

purelabinc=zeros(TL,1);

for i=1:TL
    purelabinc(i)=var(log(ylook(:,i)),invher);
end


for k=1:TL
    lincwork(k)=sum(log(inclong_workers(k,:)).*invlong_workers(k,:)./sum(invlong_workers(k,:)));
    lvarlabincwork(k)=var(log(inclong_workers(k,:)),invlong_workers(k,:));
    lthirdincwork(k)=sum(((log(inclong_workers(k,:))-...
        sum((log(inclong_workers(k,:)).*invlong_workers(k,:))/sum(invlong_workers(k,:))...
        ))./sqrt(lvarlabincwork(k))).^3.*...
        invlong_workers(k,:)./sum(invlong_workers(k,:)));
    lfourthincwork(k)=sum(((log(inclong_workers(k,:))-...
        sum((log(inclong_workers(k,:)).*invlong_workers(k,:))/sum(invlong_workers(k,:))...
        ))./sqrt(lvarlabincwork(k))).^4.*...
        invlong_workers(k,:)./sum(invlong_workers(k,:)));
    lvarconswork(k)=var(log(conslong_workers(k,:)),invlong_workers(k,:));
    cons_aworkers(k)=sum(conslong_workers(k,:).*invlong_workers(k,:)/sum(invlong_workers(k,:)));
    lcons_aworkers(k)=sum(log(conslong_workers(k,:)).*invlong_workers(k,:)/sum(invlong_workers(k,:)));
end

lvarconsold=zeros(T-TL,1);
cons_aold=lvarconsold;
lcons_aold=lvarconsold;
conslong_old(conslong_old<0 & invlong_old<1e-13)=0.1; % gets multiplied by zero anyway
% but avoid problems of infinity*0
for k=1:(T-TL)
lvarconsold(k)=var(log(conslong_old(k,:)),invlong_old(k,:));
cons_aold(k)=sum(conslong_old(k,:).*invlong_old(k,:))/sum(invlong_old(k,:));
lcons_aold(k)=sum(log(conslong_old(k,:)).*invlong_old(k,:))/sum(invlong_old(k,:));
end
lvarcons=[lvarconswork; lvarconsold];

% Now by age groups!


% Number of groups
% 25-29, 30-34, 35-39, 40-44, 45-50, 51-55, 56-60
% TL=36
stgroup=[1 6 11 16 21 27 32];
engroup=[5 10 15 20 26 31 36];
plgroup=stgroup+2;

ng=7;

lvarlabincgroup=zeros(ng,1);
lvarconsgroup=zeros(ng,1);

for k=1:ng
    lvarlabincgroup(k)=var(reshape(log(inclong_workers(stgroup(k):engroup(k),:)),...
    1,(engroup(k)-stgroup(k)+1)*(s1/TL)),...
        reshape(invlong_workers(stgroup(k):engroup(k),:),...
        1,(engroup(k)-stgroup(k)+1)*(s1/TL)));
    lvarconsgroup(k)=var(reshape(log(conslong_workers(stgroup(k):engroup(k),:)),...
         1,(engroup(k)-stgroup(k)+1)*(s1/TL)),...   
        reshape(invlong_workers(stgroup(k):engroup(k),:),...
        1,(engroup(k)-stgroup(k)+1)*(s1/TL)));
end


transiti=[cons_aworkers(TL);cons_aold(1)];
ltransiti=[lcons_aworkers(TL);lcons_aold(1)];

% 
% % Savings averages
% 
% sav_workers=permute(sopt,[1 2 3 5 4]);
% sav_aworkers=reshape(permute(sav_workers,[4 1 2 3 5]),TL,da*dbet*dy*dhip);
% sav_aworkers=sum(sav_aworkers.*inv_aworkers,2);

    ageoldyoung=[agetime;(61:agetime(1)+T-1)'];

 
 %% 2.2. Percentiles of log consumption in the model
 
 invlong_order=zeros(size(invlong_workers));
 conslong_order=zeros(size(conslong_workers));
 
 conslong_workers(invlong_workers==0 & conslong_workers<=0)=1e-10;
 
for i=1:TL
    [conslong_order(i,:),I]=sort(log(conslong_workers(i,:)),2);
    invlong_order(i,:)=invlong_workers(i,I);
end

invlong_cumorder=cumsum(invlong_order,2);
invlong_cumorder=invlong_cumorder./(repmat(invlong_cumorder(:,end),...
    1,s1/TL));

dcp=3;
cons_perctiles=zeros(dcp,TL);

for i=1:TL
   I=find(invlong_cumorder(i,:)>0.1,1,'first');
   I1=find(invlong_cumorder(i,:)>0.05,1,'first');
   I2=find(invlong_cumorder(i,:)>0.15,1,'first');
   cons_perctiles(1,i)=(conslong_order(i,I1)+conslong_order(i,I2))/2;
   I=find(invlong_cumorder(i,:)>0.5,1,'first');
   I1=find(invlong_cumorder(i,:)>0.45,1,'first');
   I2=find(invlong_cumorder(i,:)>0.55,1,'first');
   cons_perctiles(2,i)=(conslong_order(i,I1)+conslong_order(i,I2))/2;
   I=find(invlong_cumorder(i,:)>0.9,1,'first');
   I1=find(invlong_cumorder(i,:)>0.85,1,'first');
   I2=find(invlong_cumorder(i,:)>0.95,1,'first');
   cons_perctiles(3,i)=(conslong_order(i,I1)+conslong_order(i,I2))/2;
end

%% Savbyage and varage for figures

savbyage = permute(sum(sum(sum(inv5d_workers.*...
    permute(sopt,[1 2 3 4 6 5 7]),1),2),3),[5 1 2 3 4]);
savbyage=savbyage./pfrac(1:36);

varage=zeros(TL,0);

for i=1:TL
    w1=inv5d_workers(:,:,:,:,i);
    v1=sopt(:,:,:,:,i);
    varage(i)=var(v1(:),w1(:));
end


%% Now compute table of topwel botwel for paper.

percentable=zeros(7,2);
percentable(:,1)=[0.01; 0.05; 0.2; 0.4; 0.6; 0.8; 0.9];
A=find(topwel>0,1,'first');
B=find(topwel(:,1)<0.9999,1,'last');
percentable(:,2)=interp1(topwel(A:B,2),topwel(A:B,1),percentable(:,1),'linear','extrap');


end