function [WG,Ujf1,Ujc1] = gains(curpolt,newpolt);
% GAINS  compute gains by cohort for the current policy (curpolt)
%        and the new policy (newpolt) for the transition.
%
%        Ellen McGrattan, 9-23-13


beta          = curpolt(4);
gam           = curpolt(9);
nj            = curpolt(20);
ne            = curpolt(21);
T             = curpolt(22);
nh            = ne*nj;
ps            = zeros(T,nh);
ps(:)         = curpolt(28+nh:27+nh+T*nh);
ntx           = curpolt(24);
ntr           = curpolt(25);
i1            = 27+nh+3*T*nh+2*nh+32*T+3*T*(ntx+ntr);
cc            = zeros(nh,T);
lc            = zeros(nh,T);
cc(:)         = curpolt(i1+1:i1+nh*T);
lc(:)         = curpolt(i1+1+nh*T:i1+2*nh*T);
cc            = cc';
lc            = lc';


for j=2:nj
  Discl(j,1)  = sum((beta.^[j-1:nj-1]').* ...
                cumprod([diag(ps(1:nj-j+1,j-1:nj-1))]));
end;
for t=1:T-nj+1;
  Disc(t,1)   = sum((beta.^[0:nj-1]').* ...
                cumprod([1;diag(ps(t:t+nj-1,1:nj-1))]));
end;
Disc(T-nj+2:T)= ones(nj-1,1)*Disc(T-nj+1);
Disc          = [Discl(101:-1:2);Disc];

for i=1:ne
  s           = int2str(i);
  c           = cc(:,(i-1)*nj+1:i*nj);
  l           = lc(:,(i-1)*nj+1:i*nj);
  for j=2:nj;
    Cj          = diag(c(1:nj-j+1,j:nj));        
    Lj          = diag(l(1:nj-j+1,j:nj));        
    Ujlagc(j,1) = sum((beta.^[j-1:nj-1]').* ...
                  cumprod([diag(ps(1:nj-j+1,j-1:nj-1))]).* ...
                  (log(Cj)+gam*log(1-Lj)));
  end;
  for t=1:T-nj+1;
    Cj          = diag(c(t:t+nj-1,:));        
    Lj          = diag(l(t:t+nj-1,:));        
    Ujc(t,1)    = sum((beta.^[0:nj-1]').* ...
                  cumprod([1;diag(ps(t:t+nj-1,1:nj-1))]).* ...
                  (log(Cj)+gam*log(1-Lj)));
  end;
  Ujc(T-nj+2:T)  = ones(nj-1,1)*Ujc(T-nj+1);
  eval(['Ujc',s,'=Ujc;']);
  eval(['Uc',s,'= [Ujlagc(101:-1:2);Ujc];'])
end;


cf            = zeros(nh,T);
lf            = zeros(nh,T);
ntx           = newpolt(24);
ntr           = newpolt(25);
i1            = 27+nh+3*T*nh+2*nh+32*T+3*T*(ntx+ntr);
cf(:)         = newpolt(i1+1:i1+nh*T);
lf(:)         = newpolt(i1+1+nh*T:i1+2*nh*T);
cf            = cf';
lf            = lf';

for i=1:ne
  s           = int2str(i);
  c           = cf(:,(i-1)*nj+1:i*nj);
  l           = lf(:,(i-1)*nj+1:i*nj);
  for j=2:nj;
    Cj          = diag(c(1:nj-j+1,j:nj));        
    Lj          = diag(l(1:nj-j+1,j:nj));        
    Ujlagf(j,1) = sum((beta.^[j-1:nj-1]').* ...
                  cumprod([diag(ps(1:nj-j+1,j-1:nj-1))]).* ...
                  (log(Cj)+gam*log(1-Lj)));
  end;
  for t=1:T-nj+1;
    Cj          = diag(c(t:t+nj-1,:));        
    Lj          = diag(l(t:t+nj-1,:));        
    Ujf(t,1)    = sum((beta.^[0:nj-1]').* ...
                  cumprod([1;diag(ps(t:t+nj-1,1:nj-1))]).* ...
                  (log(Cj)+gam*log(1-Lj)));
  end;
  Ujf(T-nj+2:T)  = ones(nj-1,1)*Ujf(T-nj+1);
  eval(['Ujf',s,'=Ujf;']);
  eval(['Uf',s,'= [Ujlagf(101:-1:2);Ujf];'])
end;

WG = [];
for i=1:ne;
  s           = int2str(i);
  eval(['WG',s,'= (exp((Uf',s,'-Uc',s,')./Disc)-1)*100;'])
  eval(['WG=[WG,WG',s,'];'])
end;

