function [xi,xi_earn] = consumption_equivalents(bet,r,invher,ylook,...
    ye,yedtb,pens,vt,Qylife,sig)

global dy TL T TR a sur dye


azero=find(a==0,1,'first');

% Consumption equivalents

% Discount factors

cumsur=cumprod(sur);
dfact=[1; cumsur(1:end-1)].*(bet.^(0:T-1)');
disc=(1/(1+r)).^(0:T-1)';

% Average income at each age

avru_bench=zeros(T,1);

for t=1:TL
   avru_bench(t)=sum(kron(invher,yedtb(:,t)).*kron(ylook(:,t),ones(dye,1)).*...
       kron(ones(dy,1),ye(:,t)));
end

% Plus pensions

for t=TL+1:T
   avru_bench(t)=sum(kron(invher,yedtb(:,TL)).*pens);
end

% Present value of lifetime resources

pvalue_bench=sum(disc.*avru_bench);
evalue_bench=sum(kron(invher,yedtb(:,1)).*reshape(vt(azero,:,:,1,1),dye*dy,1));


% Sum of dfacts

denu=sum((disc.^(1-1./sig))./(dfact.^(-1./sig)));

% For the NL process...

consim2=((disc.^(-1/sig)).*pvalue_bench./denu)./(dfact.^(-1/sig));
actus2=consim2.*disc;
go2=complete_markets(actus2,dfact,sig,T,r);

xi=1-(evalue_bench/(-go2)).^(1/(1-sig));


% - - Next, by starting earnings gridpoint

evalue_bench_earn=permute(yedtb(:,1)'*permute(vt(azero,:,:,1,1,1),[2 3 1 4 5 6]),[2 1]);


% Present value of earnings for each starting gridpoint
pvalue_earn=zeros(dy,1);
for i=1:dy
    pvalue_earn(i)=yedtb(:,1)'*ye(:,1).*ylook(i,1);
    Qyacc=Qylife(i,:,1);
    for t=2:TL
    pvalue_earn(i)=disc(t).*(yedtb(:,t)'*ye(:,t)).*(Qyacc*ylook(:,t))+pvalue_earn(i);
    if t<TL
    Qyacc=Qyacc*Qylife(:,:,t);
    end
    end
    TR=T-TL;
    pvalue_earn(i)=pvalue_earn(i)+sum(disc(TL+1:T).*kron(Qyacc',yedtb(:,TL))'*pens);
end


consim2_earn=((disc.^(-1/sig)).*pvalue_earn'./denu)./(dfact.^(-1/sig));
go2_earn=zeros(dy,1);
for i=1:dy
actus2_earn=consim2_earn(:,i).*disc;
go2_earn(i)=complete_markets(actus2_earn,dfact,sig,T,r);
end

xi_earn=1-(evalue_bench_earn./(-go2_earn)).^(1/(1-sig));

end