function [grid,change, change_labor, change_capital]=income_quantiles_v05(ell_z,Z, distribution_ss0, distribution_ssf, Parameter, NumPoints, smooth_method)

% - pareto parameters for intiial distribution ----------------------------
invzetaH_ssf=distribution_ssf.invzetaH;
invzetaP_ssf=distribution_ssf.invzetaP;
invzetaN_ssf=distribution_ssf.invzetaN;

invzetaH_ss0=distribution_ss0.invzetaH;
invzetaP_ss0=distribution_ss0.invzetaP;
invzetaN_ss0=distribution_ss0.invzetaN;

% - returns ---------------------------------------------------------------
rIss0=distribution_ss0.rI;
rIssf=distribution_ssf.rI;

rHss0=distribution_ss0.rH;
rHssf=distribution_ssf.rH;

% - wages -----------------------------------------------------------------
wzss0=distribution_ss0.wz;
wzssf=distribution_ssf.wz;

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

% - compute counter CDF at a grid of points -------------------------------
ygrid=[wzss0 wzss0*(rIss0/rHss0) wzssf wzssf*(rIssf/rHssf) linspace(0.0001,1.99*max(wzssf), 5*NumPoints) exp(linspace(log(2*max(wzssf)),8*log(max(wzssf)), 2*NumPoints))];
ygrid=sort(ygrid);
GridSize=numel(ygrid);
ymat=repmat(ygrid,Z,1);

% - households ------------------------------------------------------------
prob_y_households=100*(ymat./wzss0').^(-1/invzetaH_ss0); % 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_ss0/(invzetaP_ss0-invzetaN_ss0)*min(((rHss0/rIss0)*ymat./wzss0').^(-1/invzetaP_ss0),ones(Z,GridSize))+ ...                 
                 +100*(invzetaN_ss0/(invzetaP_ss0-invzetaN_ss0))*((rHss0/rIss0)*ymat<repmat(wzss0',1,GridSize)).*(((rHss0/rIss0)*ymat./wzss0').^(-1/invzetaN_ss0)-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;

% - contribution of labor income ( E[labor income| income=y]---------------

ph_y=(ymat>=repmat(wzss0',1,GridSize))*(1/invzetaH_ss0).*(ymat./wzss0').^(-1/invzetaH_ss0-1)./wzss0';
pi_y=((rHss0/rIss0)*ymat<repmat(wzss0',1,GridSize))*(1/(invzetaP_ss0-invzetaN_ss0)).*((rHss0/rIss0)*ymat./wzss0').^(-1/invzetaN_ss0-1)*(rHss0/rIss0)./wzss0'+ ...
     +((rHss0/rIss0)*ymat>=repmat(wzss0',1,GridSize))*(1/(invzetaP_ss0-invzetaN_ss0)).*((rHss0/rIss0)*ymat./wzss0').^(-1/invzetaP_ss0-1)*(rHss0/rIss0)./wzss0';
 
wph_y= (1-Parameter.chi)*repmat(wzss0'.*ell_z', 1, GridSize).*ph_y;
wpi_y= Parameter.chi*repmat(wzss0'.*ell_z', 1, GridSize).*pi_y;
wp_y=sum(wph_y+wpi_y,1);

lph_y= (1-Parameter.chi)*repmat(ell_z', 1, GridSize).*ph_y;
lpi_y= Parameter.chi*repmat(ell_z', 1, GridSize).*pi_y;
lp_y=sum(lph_y+lpi_y,1);
wgrid_ss0=wp_y./lp_y;

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

% - households ------------------------------------------------------------
prob_y_households_ssf=100*(ymat./wzssf').^(-1/invzetaH_ssf); % 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 -------------------------------------------------------------
prob_y_investors_ssf=100*invzetaP_ssf/(invzetaP_ssf-invzetaN_ssf)*min(((rHssf/rIssf)*ymat./wzssf').^(-1/invzetaP_ssf),ones(Z,GridSize))+ ...                 
                     +100*(invzetaN_ssf/(invzetaP_ssf-invzetaN_ssf))*((rHssf/rIssf)*ymat<repmat(wzssf',1,GridSize)).*(((rHssf/rIssf)*ymat./wzssf').^(-1/invzetaN_ssf)-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;

% - contribution of labor income ( E[labor income| income=y]---------------

ph_y=(ymat>=repmat(wzssf',1,GridSize))*(1/invzetaH_ssf).*(ymat./wzssf').^(-1/invzetaH_ssf-1)./wzssf';
pi_y=((rHssf/rIssf)*ymat<repmat(wzssf',1,GridSize))*(1/(invzetaP_ssf-invzetaN_ssf)).*((rHssf/rIssf)*ymat./wzssf').^(-1/invzetaN_ssf-1)*(rHssf/rIssf)./wzssf'+ ...
     +((rHssf/rIssf)*ymat>=repmat(wzssf',1,GridSize))*(1/(invzetaP_ssf-invzetaN_ssf)).*((rHssf/rIssf)*ymat./wzssf').^(-1/invzetaP_ssf-1)*(rHssf/rIssf)./wzssf';
     
wph_y= (1-Parameter.chi)*repmat(wzssf'.*ell_z', 1, GridSize).*ph_y;
wpi_y= Parameter.chi*repmat(wzssf'.*ell_z', 1, GridSize).*pi_y;
wp_y=sum(wph_y+wpi_y,1);

lph_y= (1-Parameter.chi)*repmat(ell_z', 1, GridSize).*ph_y;
lpi_y= Parameter.chi*repmat(ell_z', 1, GridSize).*pi_y;
lp_y=sum(lph_y+lpi_y,1);
wgrid_ssf=wp_y./lp_y;

% - 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)));
log_w_grid_ss0=log(wgrid_ss0((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)));
log_w_grid_ssf=log(wgrid_ssf((index_init:index_last)));

grid_quantiles=[prob_y_ss0 prob_y_ssf linspace(1,99,99) linspace(99.01,99.99,99)];
grid_quantiles=sort(grid_quantiles);

%Smooth data and construct grid

%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);
log_w_grid_ss0 = interp1(prob_y_ss0, log_w_grid_ss0 ,grid_quantiles,smooth_method);
log_w_grid_ssf = interp1(prob_y_ssf, log_w_grid_ssf ,grid_quantiles,smooth_method);

%log_y_grid_ss0 = csaps(prob_y_ss0, log_y_grid_ss0 ,.01, grid_quantiles);
%log_y_grid_ssf = csaps(prob_y_ssf, log_y_grid_ssf ,.01, grid_quantiles);
%log_w_grid_ss0 = csaps(prob_y_ss0, log_w_grid_ss0 ,.01, grid_quantiles);
%log_w_grid_ssf = csaps(prob_y_ssf, log_w_grid_ssf ,.01, grid_quantiles);

% - integrate data over 98 percentiles 1-2, ..., 98-99
% - and hundreds of percentiles 99.00-99.01, 99.01--99.02, etc...

int_grid=[linspace(1,99,50) linspace(99.01,99.99,99)];

f_grid=exp(log_y_grid_ss0);
trapezoids=(grid_quantiles(2:end)-grid_quantiles(1:end-1)).*(f_grid(1:end-1)+f_grid(2:end))/2;
integral=cumsum(trapezoids);
cumpoints=interp1(grid_quantiles(2:end), integral , int_grid, 'nearest');
final_integral=[cumpoints(1) cumpoints(2:end)-cumpoints(1:end-1)];
log_y_grid_ss0=log(final_integral);


f_grid=exp(log_y_grid_ssf);
trapezoids=(grid_quantiles(2:end)-grid_quantiles(1:end-1)).*(f_grid(1:end-1)+f_grid(2:end))/2;
integral=cumsum(trapezoids);
cumpoints=interp1(grid_quantiles(2:end), integral , int_grid, 'nearest');
final_integral=[cumpoints(1) cumpoints(2:end)-cumpoints(1:end-1)];
log_y_grid_ssf=log(final_integral);


f_grid=exp(log_w_grid_ss0);
trapezoids=(grid_quantiles(2:end)-grid_quantiles(1:end-1)).*(f_grid(1:end-1)+f_grid(2:end))/2;
integral=cumsum(trapezoids);
cumpoints=interp1(grid_quantiles(2:end), integral , int_grid, 'nearest');
final_integral=[cumpoints(1) cumpoints(2:end)-cumpoints(1:end-1)];
log_w_grid_ss0=log(final_integral);


f_grid=exp(log_w_grid_ssf);
trapezoids=(grid_quantiles(2:end)-grid_quantiles(1:end-1)).*(f_grid(1:end-1)+f_grid(2:end))/2;
integral=cumsum(trapezoids);
cumpoints=interp1(grid_quantiles(2:end), integral , int_grid, 'nearest');
final_integral=[cumpoints(1) cumpoints(2:end)-cumpoints(1:end-1)];
log_w_grid_ssf=log(final_integral);

% - save data 

change=100*(exp(log_y_grid_ssf)-exp(log_y_grid_ss0))./((exp(log_y_grid_ss0)+exp(log_y_grid_ssf))/2);
change_labor=100*(exp(log_w_grid_ssf)-exp(log_w_grid_ss0))./((exp(log_y_grid_ss0)+exp(log_y_grid_ssf))/2);
change_capital=max(change, change-change_labor);
grid=int_grid;


end