function [grid_rh,change_rh]=income_quantiles_rh_v10(paths,ss0, ell_z,Z, distribution, Parameter, NumPoints, smooth_method)

% - pareto parameters for intiial distribution ----------------------------
invzetaH=distribution.invzetaH;
invzetaP=distribution.invzetaP;
invzetaN=distribution.invzetaN;

% - paths------------------------------------------------------------------
wzt=paths.wz;
rt=paths.r;
time=paths.t;
N=paths.N;
kt=paths.k;

% initial steady state
rss0 =ss0.r;
wzss0=ss0.wz;

%% Compute Delta_z as pdv of wzt discounted by rt minus w/r---------------- 

% - extend the time paths--------------------------------------------------
wzt_4ever = [wzt;repmat(wzt(end,:),2*N,1)];
rt_4ever  = [rt;repmat(rt(end),2*N,1)];
time_4ever= [time(1:end-1);linspace(time(end-1),10*time(end-1),2*N+1)'];

% - compute discount rates-------------------------------------------------
trapezoids=(time_4ever(2:end)-time_4ever(1:end-1)).*(rt_4ever(2:end)+rt_4ever(1:end-1))/2;
discount=exp(-[0; cumsum(trapezoids)]);

% - numerical integration--------------------------------------------------
discounted_wages=repmat(discount,1,Z).*wzt_4ever;
trapezoids=repmat((time_4ever(2:end)-time_4ever(1:end-1)),1,Z).*(discounted_wages(2:end,:)+discounted_wages(1:end-1,:))/2;
hz0=sum(trapezoids,1);

% - obtain revaluation effects for human wealth----------------------------
Delta_z=hz0-wzss0./rss0;

%% Compute M---common growth of effective wealth---------------------------
X_initial_RA=kt(1)+ell_z*hz0';
X_final_RA=kt(end)+ell_z*wzt(N,:)'/rt(end);
M=X_final_RA/X_initial_RA;

%% Compute initial dist of income as y=rx----------------------------------

% - compute counter CDF at a grid of points -------------------------------
ygrid=[wzss0 M*(wzss0+rss0*Delta_z) linspace(0,1.99*max(wzss0), 5*NumPoints) exp(linspace(log(2*max(wzss0)),8*log(max(wzss0)), 2*NumPoints))];
ygrid=sort(ygrid);
GridSize=numel(ygrid);
ymat=repmat(ygrid,Z,1);

% - households ------------------------------------------------------------
prob_y_households=100*(ymat./wzss0').^(-1/invzetaH); % using single-sided Pareto tail
prob_y_households=ell_z'.*min(prob_y_households,100*ones(Z,GridSize));
prob_y_households=sum(prob_y_households,1);
prob_y_households=(1-Parameter.chi)*(100-prob_y_households);

% - investors -------------------------------------------------------------
prob_y_investors=100*invzetaP/(invzetaP-invzetaN)*min((ymat./wzss0').^(-1/invzetaP),ones(Z,GridSize))+ ...                 
                 +100*(invzetaN/(invzetaP-invzetaN))*(ymat<repmat(wzss0',1,GridSize)).*((ymat./wzss0').^(-1/invzetaN)-1); 
prob_y_investors=ell_z'.*min(prob_y_investors,100*ones(Z,GridSize));
prob_y_investors=sum(prob_y_investors,1);
prob_y_investors=Parameter.chi*(100-prob_y_investors);

prob_y_ss0=prob_y_households+prob_y_investors;

%% Compute final dist of income as y'=r*(Mx+Delta_z)-----------------------

% - households at new st st------------------------------------------------
Delta_mat=repmat(Delta_z',1,GridSize);

prob_y_households_ssf=100*((ymat/M-rss0*Delta_mat)./wzss0').^(-1/invzetaH); % using single-sided Pareto tail
prob_y_households_ssf=ell_z'.*min(prob_y_households_ssf,100*ones(Z,GridSize));
prob_y_households_ssf=sum(prob_y_households_ssf,1);
prob_y_households_ssf=(1-Parameter.chi)*(100-prob_y_households_ssf);

% - investors at new st st-------------------------------------------------
prob_y_investors_ssf=100*invzetaP/(invzetaP-invzetaN)*min(((ymat/M-rss0*Delta_mat)./wzss0').^(-1/invzetaP),ones(Z,GridSize))+ ...                 
                     +100*(invzetaN/(invzetaP-invzetaN))*((ymat/M-rss0*Delta_mat)<repmat(wzss0',1,GridSize)).*(((ymat/M-rss0*Delta_mat)./wzss0').^(-1/invzetaN)-1); 
prob_y_investors_ssf=ell_z'.*min(prob_y_investors_ssf,100*ones(Z,GridSize));
prob_y_investors_ssf=sum(prob_y_investors_ssf,1);
prob_y_investors_ssf=Parameter.chi*(100-prob_y_investors_ssf);

prob_y_ssf=prob_y_households_ssf+prob_y_investors_ssf;

% - remove 1s -------------------------------------------------------------
index_last=find(prob_y_ss0>=99.999,1,'first');
index_init=find(prob_y_ss0>0.1,1,'first');
prob_y_ss0=prob_y_ss0(index_init:index_last);
log_y_grid_ss0=log(ygrid((index_init:index_last)));

index_last=find(prob_y_ssf>=99.999,1,'first');
index_init=find(prob_y_ssf>0.1,1,'first');
prob_y_ssf=prob_y_ssf(index_init:index_last);
log_y_grid_ssf=log(ygrid((index_init:index_last)));

grid_quantiles=[prob_y_ss0 prob_y_ssf];
grid_quantiles=sort(grid_quantiles);

%Interpolate findings: 
log_y_grid_ss0 = interp1(prob_y_ss0, log_y_grid_ss0 ,grid_quantiles,smooth_method);
log_y_grid_ssf = interp1(prob_y_ssf, log_y_grid_ssf ,grid_quantiles,smooth_method);

change_rh=100*(exp(log_y_grid_ssf)-exp(log_y_grid_ss0))./exp(log_y_grid_ss0);
grid_rh=grid_quantiles;

%Compute change over coarser grid for the one-panel figure
if Parameter.onepanel==1
   int_grid = [1 5 10:10:90 95 97.5 99 99.5 99.9 99.99 99.99999];
   change_rh=interp1(grid_rh, change_rh , int_grid, smooth_method); 
   change_rh= (change_rh(1:end-1)+ change_rh(2:end))/2;
   grid_rh=int_grid(2:end);
end    

end