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

% - pareto parameters for initial 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;

% - Volatility of investment parameter ------------------------------------

kappanu_ss0 = distribution_ss0.kappanu;
kappanu_ssf = distribution_ssf.kappanu;
alpha_capitalrisk = Parameter.alpha_capitalrisk;

% - Plotting option -------------------------------------------------------

plot_intermediate_graphs = Parameter.plot_intermediate_graphs;

% - 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;


%% INITIALIZE GRIDS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% - compute counter CDF at a grid of points -------------------------------
ygrid=[ -linspace(0.001,max(wzssf),100) 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(unique(ygrid));

GridSize=numel(ygrid);
ymat=repmat(ygrid,Z,1);

% Format grid for wage
wzss0mat = repmat(wzss0,length(ygrid),1);
wzssfmat = repmat(wzssf,length(ygrid),1);

% Compute average reinjection point
reinjection0 = wzss0(1:10:end)/rHss0;
reinjectionf = wzssf(1:10:end)/rHss0;


%% COMPUTE INITIAL DIST OF INCOME USING y=x0(r+sigma_y W), W~N(0,1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% - households ------------------------------------------------------------
prob_y_households=100*((ymat./wzss0mat').^(-1/invzetaH_ss0).*(ymat>=0)+(ymat<0)); % using single-sided Pareto tail (from lemma A7)
prob_y_households=ell_z*min(prob_y_households,100*ones(Z,GridSize));
prob_y_households=(1-Parameter.chi)*(100-prob_y_households);

% - investors -------------------------------------------------------------

ygrid1 = linspace(-3,5,1000);
ygrid1 = sort(unique(ygrid1));

dy = diff(ygrid1);
dy = [dy,dy(end)];

integrationgrid = [linspace(1e-6,100,1000) linspace(100,5000,1000)];
integrationgrid = sort(unique(integrationgrid));
d = diff(integrationgrid);
d = [d,d(end)];

f_x0 = 1/(invzetaP_ss0-invzetaN_ss0)*integrationgrid.^(-1/invzetaP_ss0-1).*(integrationgrid>1)+ ...                 
                 +(1/(invzetaP_ss0-invzetaN_ss0))*(integrationgrid<1).*(integrationgrid>0).*(integrationgrid.^(-1/invzetaN_ss0-1)); %lemma A7

clear y1_cdf

% Compute CCDF using lemma A8

for i = 1:length(ygrid1) 
    y = ygrid1(i);
    if alpha_capitalrisk~=0
        fnorm = normcdf(y./integrationgrid,rIss0,alpha_capitalrisk*kappanu_ss0);
    else
        fnorm = heaviside(y./integrationgrid-rIss0);
    end    
    % Compute convolution
    integrand = f_x0.*(1-fnorm);
    y1_ccdf(i) = sum(integrand.*d);
end

% Compute PDF
y1_dens = [0,diff((1-y1_ccdf)./dy)];

prob_y_investors = 100*interp1(ygrid1,y1_ccdf,rHss0*ymat./wzss0mat');
index_low = isnan(prob_y_investors);
index_high = (rHss0*ymat./wzss0mat'>ygrid1(end));
prob_y_investors(index_low) = 100; prob_y_investors(index_high) = 0;
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./wzss0mat').^(-1/invzetaH_ss0-1)./wzss0mat'.*(ymat>=0);

pi_y = interp1(ygrid1,y1_dens,rHss0*ymat./wzss0mat');
index = isnan(pi_y);
pi_y(index) = 0; % the nans come from interpolating outside of the bounds of ymat(i,:)

delta_ymat = [0,diff(ygrid)];
integration_scaling = pi_y*delta_ymat';
integration_scaling = repmat(integration_scaling,1,length(ygrid));
pi_y = pi_y./integration_scaling; %Rescale so that the pdf integrates to 1

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 USING y=x0(r+sigma_y W), W~N(0,1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% - households ------------------------------------------------------------
prob_y_households_ssf=100*((ymat./wzssfmat').^(-1/invzetaH_ssf).*(ymat>=0)+(ymat<0)); % using single-sided Pareto tail (from lemma A7)
prob_y_households_ssf=ell_z*min(prob_y_households_ssf,100*ones(Z,GridSize));
prob_y_households_ssf=(1-Parameter.chi)*(100-prob_y_households_ssf);

% - investors -------------------------------------------------------------

f_x0 = 1/(invzetaP_ssf-invzetaN_ssf)*integrationgrid.^(-1/invzetaP_ssf-1).*(integrationgrid>1)+ ...                 
                 +(1/(invzetaP_ssf-invzetaN_ssf))*(integrationgrid<1).*(integrationgrid>0).*(integrationgrid.^(-1/invzetaN_ssf-1)); %lemma A7

clear y1_cdf

% Compute CCDF using lemma A8

for i = 1:length(ygrid1)   
    y = ygrid1(i);
    if alpha_capitalrisk~=0
        fnorm = normcdf(y./integrationgrid,rIssf,alpha_capitalrisk*kappanu_ssf);
    else
        fnorm = heaviside(y./integrationgrid-rIssf);
    end
    % Compute convolution
    integrand = f_x0.*(1-fnorm);
    y1_ccdf(i) = sum(integrand.*d);
end

% Compute PDF
y1_dens = [diff((1-y1_ccdf)./dy),0];

prob_y_investors_ssf = 100*interp1(ygrid1,y1_ccdf,rHssf*ymat./wzssfmat');
index_low = isnan(prob_y_investors_ssf);
index_high = (rHssf*ymat./wzssfmat'>ygrid1(end));
prob_y_investors_ssf(index_low) = 100; prob_y_investors_ssf(index_high) = 0;
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./wzssfmat').^(-1/invzetaH_ssf-1)./wzssfmat';

pi_y = interp1(ygrid1,y1_dens,rHssf*ymat./wzssfmat');
index = isnan(pi_y);
pi_y(index) = 0; % the nans come from interpolating outside of the bounds of ymat(i,:)

delta_ymat = [0,diff(ygrid)];
integration_scaling = pi_y*delta_ymat';
integration_scaling = repmat(integration_scaling,1,length(ygrid));
pi_y = pi_y./integration_scaling; %Rescale so that the pdf integrates to 1

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;



%% COMPUTE CHANGE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% - remove duplicates in the CDF (required for interpolation) -------------
[interp_prob_y_ss0,index] = unique(prob_y_ss0);
interp_ygrid_ss0 = ygrid(index);
interp_wgrid_ss0 = wgrid_ss0(index);

[interp_prob_y_ssf,index] = unique(prob_y_ssf);
interp_ygrid_ssf = ygrid(index);
interp_wgrid_ssf = wgrid_ssf(index);

% Get reinjection points on the quantile grids

reinjection0_quantile = interp1(interp_ygrid_ss0,interp_prob_y_ss0,reinjection0);
reinjectionf_quantile = interp1(interp_ygrid_ssf,interp_prob_y_ssf,reinjectionf);


%Squeeze the data using arcshinh for better interpolation

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

asinh_ygrid_ss0 = asinh(interp_ygrid_ss0);
asinh_ygrid_ssf = asinh(interp_ygrid_ssf);
ashinh_wgrid_ss0 = asinh(interp_wgrid_ss0);
ashinh_wgrid_ssf = asinh(interp_wgrid_ssf);

% Interpolate findings to get quantile functions
asinh_y_grid_ss0 = interp1(interp_prob_y_ss0, asinh_ygrid_ss0 ,grid_quantiles,smooth_method);
asinh_y_grid_ssf = interp1(interp_prob_y_ssf, asinh_ygrid_ssf ,grid_quantiles,smooth_method);
asinh_w_grid_ss0 = interp1(interp_prob_y_ss0, ashinh_wgrid_ss0 ,grid_quantiles,smooth_method);
asinh_w_grid_ssf = interp1(interp_prob_y_ssf, ashinh_wgrid_ssf ,grid_quantiles,smooth_method);

%Smooth using integration on a less fine grid

%threshold = 0.01*interp1(interp_prob_y_ss0,interp_ygrid_ss0,50);
min_grid = ceil(interp1(interp_ygrid_ss0,interp_prob_y_ss0,1)); %drop bellow 1% of the median income


if Parameter.onepanel==1
    int_grid = [5 10:10:90 95 97.5 99 99.5 99.9 99.99 99.99999];
else
    int_grid = [linspace(min_grid,3,5) linspace(3.1,70,40) linspace(70,99,15) linspace(99.01,99.99,99)];
end

int_grid = unique(sort(int_grid));
int_grid_diff = diff(int_grid);
int_grid_diff = [int_grid_diff(1),int_grid_diff];

%int_grid=[linspace(1,99,50) linspace(99.01,99.99,99)];
%{
asinh_y_grid_ss0_smth = interp1(grid_quantiles,asinh_y_grid_ss0,int_grid,smooth_method);
asinh_y_grid_ssf_smth = interp1(grid_quantiles,asinh_y_grid_ssf,int_grid,smooth_method);
asinh_w_grid_ss0_smth = interp1(grid_quantiles,asinh_w_grid_ss0,int_grid,smooth_method);
asinh_w_grid_ssf_smth = interp1(grid_quantiles,asinh_w_grid_ssf,int_grid,smooth_method);
%}



%Smooth by integrating

f_grid=sinh(asinh_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)];
asinh_y_grid_ss0_smth=asinh(final_integral./int_grid_diff);

f_grid=sinh(asinh_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)];
asinh_y_grid_ssf_smth=asinh(final_integral./int_grid_diff);


f_grid= min(max(sinh(asinh_w_grid_ss0),-1e6),1e6);%censor unreasonably high values
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)];
asinh_w_grid_ss0_smth=asinh(final_integral./int_grid_diff);


f_grid=min(max(sinh(asinh_w_grid_ssf),-1e6),1e6); %censor unreasonably high values
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)];
asinh_w_grid_ssf_smth=asinh(final_integral./int_grid_diff);

%%
% SAVE DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%change = (sinh(asinh_y_grid_ssf)-sinh(asinh_y_grid_ss0))./((sinh(asinh_y_grid_ss0)+sinh(asinh_y_grid_ssf))/2)*100;
%change_labor = (sinh(asinh_w_grid_ssf)-sinh(asinh_w_grid_ss0))./((sinh(asinh_w_grid_ss0)+sinh(asinh_w_grid_ssf))/2)*100;
%change_capital = change-change_labor;

change = (sinh(asinh_y_grid_ssf_smth)-sinh(asinh_y_grid_ss0_smth))./sinh(asinh_y_grid_ss0_smth)*100;
change_labor = (sinh(asinh_w_grid_ssf_smth)-sinh(asinh_w_grid_ss0_smth))./sinh(asinh_y_grid_ss0_smth)*100;
change_capital = change-change_labor;

change_nosmooth = (sinh(asinh_y_grid_ssf)-sinh(asinh_y_grid_ss0))./sinh(asinh_y_grid_ss0)*100;
change_labor_nosmooth = (sinh(asinh_w_grid_ssf)-sinh(asinh_w_grid_ss0))./sinh(asinh_y_grid_ss0)*100;
change_capital_nosmooth = change_nosmooth-change_labor_nosmooth;

grid=int_grid;


%% PLOT INTERMEDIATE GRAPHS

if plot_intermediate_graphs ==1


%%
close
figure
hold on
plot(grid,sinh(asinh_y_grid_ss0_smth),'LineWidth',2)
plot(grid,sinh(asinh_y_grid_ssf_smth),'LineWidth',2)
legend({'1980','2012'},'Interpreter','Latex','FontSize',13)
xlabel('Quantiles','Interpreter','Latex','FontSize',18)
ylabel('Income','Interpreter','Latex','FontSize',18)
set(gca, 'YScale', 'log')
export_fig '../output_nov/income_quantiles'.pdf
%%
close
figure
hold on
plot(grid,sinh(asinh_y_grid_ss0_smth),'LineWidth',2)
plot(grid,sinh(asinh_y_grid_ssf_smth),'LineWidth',2)
xlim([99 100])
legend({'1980','2012'},'Interpreter','Latex','FontSize',13)
xlabel('Quantiles','Interpreter','Latex','FontSize',18)
ylabel('Income','Interpreter','Latex','FontSize',18)
set(gca, 'YScale', 'log')
%export_fig '../output_nov/income_quantiles_top1p'.pdf 

%%
close
figure
hold on
plot(grid,sinh(asinh_w_grid_ss0_smth),'LineWidth',2)
plot(grid,sinh(asinh_w_grid_ssf_smth),'LineWidth',2)
legend({'1980','2012'},'Interpreter','Latex','FontSize',13)
xlabel('Quantiles','Interpreter','Latex','FontSize',18)
ylabel('Wage','Interpreter','Latex','FontSize',18)
export_fig '../output_nov/wage_quantiles'.pdf

close
figure
hold on
plot(grid,sinh(asinh_w_grid_ss0_smth),'LineWidth',2)
plot(grid,sinh(asinh_w_grid_ssf_smth),'LineWidth',2)
xlim([99 100])
legend({'1980','2012'},'Interpreter','Latex','FontSize',13)
xlabel('Quantiles','Interpreter','Latex','FontSize',18)
ylabel('Wage','Interpreter','Latex','FontSize',18)
%export_fig '../output_nov/wage_quantiles_top1'.pdf

%%
close
figure
hold on
plot(grid_quantiles,change_labor_nosmooth,'LineWidth',2)
plot(grid_quantiles,change_capital_nosmooth,'LineWidth',2)
plot(grid_quantiles,change_nosmooth,'LineWidth',3,'Color','black')
legend({'labor','capital','total'},'Interpreter','Latex','FontSize',13)
ylim([-20 80])
xlabel('Quantiles','Interpreter','Latex','FontSize',18)
ylabel('Income','Interpreter','Latex','FontSize',18)
%export_fig '../output_nov/fine_change_quantiles'.pdf 

%%

close
figure
hold on
plot(grid_quantiles,change_labor_nosmooth,'LineWidth',2)
xlim([85 100])
ylim([-40 80])
xlabel('Quantiles','Interpreter','Latex','FontSize',18)
ylabel('Wage','Interpreter','Latex','FontSize',18)
%export_fig '../output_nov/wage_upper_tail'.pdf 

%%

close
figure
hold on
plot(grid,change,'LineWidth',2)
scatter(reinjection0_quantile,interp1(grid,change,reinjection0_quantile));
scatter(reinjectionf_quantile,interp1(grid,change,reinjectionf_quantile),'*');
ylim([-15 100])
legend({'Total growth','Reinjection 1980','Reinjection 2012'},'Interpreter','Latex','FontSize',13)
xlabel('Quantiles','Interpreter','Latex','FontSize',18)
ylabel('Income','Interpreter','Latex','FontSize',18)
export_fig '../output_nov/income_quantiles_reinjection'.pdf


end

end

