% -------------------------------------------------------------------------
% Project: Uneven Growth
%
% -------------------------------------------------------------------------

close all
clear all

% change location to the folder where this m file is saved
mfile_name          = mfilename('fullpath');
[pathstr,name,ext]  = fileparts(mfile_name);
cd(pathstr);

% add all routines in the subfolders
cd ..
addpath 'output'
addpath 'code'

% latex interpreter
set(groot,'defaulttextinterpreter','latex');  
set(groot, 'defaultAxesTickLabelInterpreter','latex');  
set(groot, 'defaultLegendInterpreter','latex');

% - solver default options ------------------------------------------------
options=optimset('disp','off','LargeScale','off', ...
                 'TolFun',.00000001,'MaxFunEvals',10^7,'MaxIter',10^7);

%% - Choice of Parameters -------------------------------------------------

binding=0; %is the borrowing constraint binding (binding when risk aversion low). set to 0 or 1. Baseline is 0, with target of kappa=1.3.

if (binding==1)
    nu = 0.06;              % Capital income risk
    gamma = 0.1;            % Risk aversion of investors
    Parameter.chi = 0.045;  % Share of Investors (target 1/zeta=0.54)
    Parameter.lambda    = 2/3;   % Lambda=1-theta (Credit constraint)
else
    nu = 0.077;              % Capital income risk
    gamma = 2;              % Risk aversion of investors
    Parameter.chi = 0.066;  % Share of Investors (target 1/zeta=0.54)
    Parameter.lambda    = 1/2;   % Lambda=1-theta (Credit constraint)
end

%Parameter.lambda    = 2/3;              % Lambda=1-theta (Credit constraint)
Parameter.gamnu     = gamma * nu^2;     % Gamma x nu^2
Parameter.rho       = 0.000;            % Rho (Time preference)
Parameter.sigma     = 2.000;            % Sigma (Inverse of Intertemporal Elasticity of Substitution)
Parameter.g         = 0.015;            % Economy wide growth rate along the BGP
Parameter.p         = 0.045;            % Reset probability 
Parameter.govbonds  = 0.000;            % Total Government bonds (zero in notes)
Parameter.delta     = 0.050;            % Depreciation rate

Parameter.alpha_capitalrisk = 1;        % Auxiliary control, 0 => shuts down capital risk
Parameter.plot_intermediate_graphs =0;  % Plot the auxiliary graphs
Parameter.onepanel = 1;                 % set to 1 to plot figure on one panel only


% - feed in the shock to the alphas and A----------------------------------

    data_sim = csvread('Data/data_calibration_matlab_pi30.csv',1,0);

    % alpha's
    cal_alpha_0  = data_sim(:,3)';
    cal_alpha_f  = data_sim(:,5)'; 
        
    % eta's
    cal_eta_0  = data_sim(:,2)'/sum(data_sim(:,2)); %note: we are keeping etas fixed at 1980 values throughout! 
    cal_eta_f  = data_sim(:,2)'/sum(data_sim(:,2)); 
    
    % psi_z's
    cal_psi_L     = data_sim(:,8)';
    
    % A term
    cal_A_0 = data_sim(1,7);       
    cal_A_f = data_sim(1,7); 
    
    % Skill groups
    cal_Z           = numel(data_sim(:,1));   % number of skill groups
    cal_ell_z       = ones(1,cal_Z)/cal_Z;            % equal sized skill groups 
    types           = linspace(1,cal_Z,cal_Z);
    
    % Target for r_k
    rk_target=data_sim(1,6)-Parameter.delta;

%% - Routine to compute equilibrium wages and returns in extended model----

% - short routine to calibrate rho----------------------------------------- 
Parameter.alpha=cal_alpha_0*cal_eta_0';
Parameter.rho=calibration_routine_v10(rk_target, Parameter);
assert(Parameter.rho+(Parameter.sigma-1)*Parameter.g>0); % Check usual condition for finite utility/wealth being finite

% - compute initial steady state ------------------------------------------
Parameter.alpha=cal_alpha_0*cal_eta_0';
alpha_ss0     = Parameter.alpha;



% - Initialize Guess ------------------------------------------------------
% - Note: equations_v10 accepts log(k), b, log(rb-g) ----------------------
Solution0    = zeros(3,1);
Solution0(1) = log((Parameter.alpha/(1-Parameter.alpha)) ...
              /(Parameter.delta+Parameter.rho+Parameter.sigma*Parameter.g)); 
             % log of normalized capital (k in notes)
Solution0(2) = Parameter.govbonds;                                          
             % normalized bond holdings of households (b in notes)
Solution0(3) = log(Parameter.rho+(Parameter.sigma-1)*Parameter.g);            
             % log of return to bonds (rb in notes) - g

% - load solver -----------------------------------------------------------
Solution = fsolve(@(Solution)equations_v10(Solution,Parameter),Solution0,options);
Check=equations_v10(Solution,Parameter);
assert(max(abs(Check))<0.0001);


% - extract steady state variables-----------------------------------------

    k_norm_ss0       = exp(Solution(1));
    b_h_ss0          = Solution(2);
    b_e_ss0          = Parameter.govbonds - b_h_ss0;
    r_b_ss0          = Parameter.g+exp(Solution(3));
    r_k_ss0          = (Parameter.alpha/(1-Parameter.alpha))*(r_b_ss0-Parameter.g)/k_norm_ss0 - Parameter.delta;
    r_dif_ss0        = r_k_ss0 - r_b_ss0;
    kappas_ss0       = min(1/Parameter.lambda,(r_dif_ss0)/Parameter.gamnu);
    r_bar_ss0        = kappas_ss0 * r_k_ss0 + (1-kappas_ss0)*r_b_ss0;
    r_w_ss0           =r_bar_ss0+ 0.5 * (Parameter.sigma-1) * Parameter.gamnu * kappas_ss0^2;
    
% - wages in initial steady state -----------------------------------------
    
    Lwtd=prod((cal_psi_L.*cal_ell_z./(cal_eta_0.*(1-cal_alpha_0))).^(cal_eta_0.*(1-cal_alpha_0)));
    k_hour_ss0=Parameter.alpha*(cal_A_0*Lwtd/(r_k_ss0+Parameter.delta))^(1/(1-Parameter.alpha));
    y_hour_ss0=cal_A_0*((k_hour_ss0/Parameter.alpha)^Parameter.alpha)*Lwtd;
    w_hour_ss0=cal_eta_0.*(1-cal_alpha_0).*y_hour_ss0./cal_ell_z;
    w_ave_hour_ss0=mean(w_hour_ss0);
   
% - wealth distribution in initial steady state----------------------------
    
    muI_ss0=(r_w_ss0-Parameter.rho)/Parameter.sigma-Parameter.g-0.5*nu^2*kappas_ss0^2;
    muH_ss0=(r_b_ss0-Parameter.rho)/Parameter.sigma-Parameter.g;
    
    distribution_ss0.invzetaH=muH_ss0/Parameter.p;
    distribution_ss0.invzetaP=(muI_ss0+sqrt(muI_ss0^2+2*nu^2*kappas_ss0^2*Parameter.p))/(2*Parameter.p);
    distribution_ss0.invzetaN=(muI_ss0-sqrt(muI_ss0^2+2*nu^2*kappas_ss0^2*Parameter.p))/(2*Parameter.p);
    distribution_ss0.wz=w_hour_ss0;
    distribution_ss0.rI=r_bar_ss0;
    distribution_ss0.rH=r_b_ss0;
    distribution_ss0.kappanu=kappas_ss0*nu;
    
   
% - compute final steady state --------------------------------------------
Parameter.alpha=cal_alpha_f*cal_eta_f';
alpha_ssf      = Parameter.alpha;

% - Initialize Guess ------------------------------------------------------
Solution0    = Solution;
% - load solver -----------------------------------------------------------
Solution = fsolve(@(Solution)equations_v10(Solution,Parameter),Solution0,options);
Check=equations_v10(Solution,Parameter);
assert(max(abs(Check))<0.0001);

% - extract steady state variables-----------------------------------------

    k_norm_ssf       = exp(Solution(1));
    b_h_ssf          = Solution(2);
    b_e_ssf          = Parameter.govbonds - b_h_ssf;
    r_b_ssf          = Parameter.g+exp(Solution(3));
    r_k_ssf          = (Parameter.alpha/(1-Parameter.alpha))*(r_b_ssf-Parameter.g)/k_norm_ssf - Parameter.delta;
    r_dif_ssf        = r_k_ssf - r_b_ssf;
    kappas_ssf       = min(1/Parameter.lambda,(r_dif_ssf)/Parameter.gamnu);
    r_bar_ssf        = kappas_ssf * r_k_ssf + (1-kappas_ssf)*r_b_ssf;
    r_w_ssf           =r_bar_ssf+ 0.5 * (Parameter.sigma-1) * Parameter.gamnu * kappas_ssf^2;
    
% - wages in final steady state -------------------------------------------
    
    Lwtd=prod((cal_psi_L.*cal_ell_z./(cal_eta_f.*(1-cal_alpha_f))).^(cal_eta_f.*(1-cal_alpha_f)));
    k_hour_ssf=Parameter.alpha*(cal_A_f*Lwtd/(r_k_ssf+Parameter.delta))^(1/(1-Parameter.alpha));
    y_hour_ssf=cal_A_f*((k_hour_ssf/Parameter.alpha)^Parameter.alpha)*Lwtd;
    w_hour_ssf=cal_eta_f.*(1-cal_alpha_f).*y_hour_ssf./cal_ell_z;
    w_ave_hour_ssf=mean(w_hour_ssf);
   
% - wealth distribution in final steady state------------------------------

    muI_ssf=(r_w_ssf-Parameter.rho)/Parameter.sigma-Parameter.g-0.5*nu^2*kappas_ssf^2;
    muH_ssf=(r_b_ssf-Parameter.rho)/Parameter.sigma-Parameter.g;
    
    distribution_ssf.invzetaH=muH_ssf/Parameter.p;
    distribution_ssf.invzetaP=(muI_ssf+sqrt(muI_ssf^2+2*nu^2*kappas_ssf^2*Parameter.p))/(2*Parameter.p);
    distribution_ssf.invzetaN=(muI_ssf-sqrt(muI_ssf^2+2*nu^2*kappas_ssf^2*Parameter.p))/(2*Parameter.p);
    distribution_ssf.wz=w_hour_ssf;
    distribution_ssf.rI=r_bar_ssf;
    distribution_ssf.rH=r_b_ssf;
    distribution_ssf.kappanu=kappas_ssf*nu;
    
    distribution_ss0.invzetaP
    distribution_ssf.invzetaP
    kappas_ss0
    kappas_ssf
    
%% - Transitional dynamics in basic version of the model ------------------    

% - calibration
Calibration.A_0      =cal_A_0;
Calibration.A_f      =cal_A_f;
Calibration.eta_0    =cal_eta_0;
Calibration.eta_f    =cal_eta_f;
Calibration.alpha_0  =cal_alpha_0;
Calibration.alpha_f  =cal_alpha_f;
Calibration.ell_z    =cal_ell_z;
Calibration.psi_L    =cal_psi_L;
Calibration.rk_target=rk_target;

% - our model -------------------------------------------------------------
[paths,ss0, elasticity]=transition_v10(Parameter, Calibration);

% - representative housheold version --------------------------------------
Parameter_rh=Parameter;
Parameter_rh.p=0;
[paths_rh,ss0_rh]=transition_v10(Parameter_rh, Calibration);


%% - Plot of transitional dynamics ----------------------------------------

% - plot the shock --------------------------------------------------------
figure_alpha = figure('units','centimeters','position',[10 10 7*(16/9) 7]);

legendinfo = cell(1,2); 
hold on
box on
plot(types, cal_alpha_f*100, 'r--', 'linewidth',1.2)
legendinfo{1} = 'New st st';
plot(types, cal_alpha_0*100, 'k-', 'linewidth', 1.2)
legendinfo{end} = '1980';
ylim([30,50])
xlabel('Wage percentile, $z$')
ylabel('percent')
title('Calibrated share of automated tasks $\alpha_z(t)$')
legend(legendinfo,'location','northeast')
legend('boxoff')

set(gcf, 'Color', 'w');
export_fig 'output/Figure_paper_shock'.pdf 

% - plot comparisson between rep housheold model and ours -----------------

N=paths.N;
timec = [(-10:0.1:0)';paths.t(2:end)];
Nnew = numel(timec);
final_year=2060;

figure_dynamics = figure('units','centimeters','outerposition',[5 5 2*17*0.8 3*9*0.8]);

subplot(3,2,1)
hold on
box on
plot(timec+1980,100*[repmat(paths.ls(1,1),Nnew-N,1);paths.ls],'k-', 'Linewidth',1.2)
plot(timec+1980,100*[repmat(paths_rh.ls(1,1),Nnew-N,1);paths_rh.ls],'r--', 'Linewidth',1.2)
xlim([timec(1)+1980,final_year])
ylim([55,70])
title('Labor share')
ylabel('percent')

subplot(3,2,2)
hold on
box on
plot(timec+1980,[repmat(1,Nnew-N,1);paths.y/paths.y(1,1)],'k-', 'Linewidth',1.2)
plot(timec+1980,[repmat(1,Nnew-N,1);paths_rh.y/paths_rh.y(1,1)],'r--', 'Linewidth',1.2)
xlim([timec(1)+1980,final_year])
ylim([0.95,1.25])
title('Output per worker')
ylabel('relative to initial st st')
x1 = 2030;
y1 = 1.2;
txt1 = 'Rep. household';
text(x1,y1,txt1,'Color','r');
x2 = 2030;
y2 = 1.09;
txt2 = 'Diss. shocks';
text(x2,y2,txt2,'Color','k');


subplot(3,2,3)
hold on
box on
plot(timec+1980,[repmat(paths.ktoy(1,1),Nnew-N,1);paths.ktoy],'k-', 'Linewidth',1.2)
plot(timec+1980,[repmat(paths_rh.ktoy(1,1),Nnew-N,1);paths_rh.ktoy],'r--', 'Linewidth',1.2)
title('Capital-output ratio')
ylabel('ratio')
xlim([timec(1)+1980,final_year])
ylim([2.8,3.8])

subplot(3,2,4)
hold on
box on
plot(timec+1980,[repmat(100*ss0.i,Nnew-N,1);100*paths.i],'k-', 'Linewidth',1.2)
plot(timec+1980,[repmat(100*ss0_rh.i,Nnew-N,1);100*paths_rh.i],'r--', 'Linewidth',1.2)
xlim([timec(1)+1980,final_year])
ylim([14,20]);
title('Gross investment rate')
ylabel('percent')

subplot(3,2,5)
hold on
box on
plot(timec+1980,[repmat(100*paths.r(1,1),Nnew-N,1);100*paths.r],'k-', 'Linewidth',1.2)
plot(timec+1980,[repmat(100*paths_rh.r(1,1),Nnew-N,1);100*paths_rh.r],'r--', 'Linewidth',1.2)
xlim([timec(1)+1980,final_year])
ylim([6,8])
title('Return to wealth')
ylabel('percent')

subplot(3,2,6)
hold on
box on
plot(timec+1980,[repmat(1,Nnew-N,1);paths.w_ave/paths.w_ave(1,1)],'k-', 'Linewidth',1.2)
plot(timec+1980,[repmat(1,Nnew-N,1);paths_rh.w_ave/paths_rh.w_ave(1,1)],'r--', 'Linewidth',1.2)
xlim([timec(1)+1980,final_year])
title('Average real wage', 'Interpreter','latex')
ylabel('relative to initial st st', 'Interpreter','latex')

set(gcf, 'Color', 'w');

export_fig 'output/Figure_paper_transition'.pdf 

%% - Summary of st-st comparissons ----------------------------------------

% - report sum stats ------------------------------------------------------

display(['Calibrated rho=',num2str(Parameter.rho)]);
display(['Productivity gains automation=',num2str(mean(w_hour_ss0./cal_psi_L / (r_k_ss0+Parameter.delta)))]);
display(['Initial labor share=',num2str(cal_alpha_0*cal_eta_f')]);
display(['Final labor share=',num2str(cal_alpha_f*cal_eta_f')]);
display(['Initial output per hour=',num2str(y_hour_ss0)]);
display(['Initial capital per hour=',num2str(k_hour_ss0)]);
display(['Initial wage per hour=',num2str(w_ave_hour_ss0)]);

% - wages (Figure XX) -----------------------------------------------------

wages_data = csvread("Data/observed_wages.csv",1,0); %observed wage changes
wage_percentile=linspace(1,100,100);                      %wage percentiles

figure_wages = figure('units','centimeters','outerposition',[10 10 16*1.2 9*1.2]);

hold on
box on
yyaxis left
h1=plot(wage_percentile, 100*(log(w_hour_ssf)-log(w_hour_ss0)),'k-', 'linewidth',1.2);
h2=plot(wage_percentile, zeros(1,100), ':','linewidth',1,'color',[0,0,0 ]+0.5 );
set(gca,'ycolor','k');
ylabel('wage changes in model, percent')
ylim([-14,11])

yyaxis right
h3=plot(wage_percentile, 100*wages_data(:,2),'r--', 'linewidth',1.2);
xlabel('wage percentile, z')
ylabel('observed wage changes, percent')
title('Change in wages by percentile of the wage distribution')
legend([h1 h3],{'Change in model (left axis)',  'Observed change (right axis)'}, 'location', 'best')
legend('boxoff')
set(gcf, 'Color', 'w');
set(gca,'ycolor','k');


if (binding==1)
export_fig 'output/Figure_paper_wages_binding'.pdf 
else
export_fig 'output/Figure_paper_wages_notbinding'.pdf 
end

% - income disttribution (Figure XX) --------------------------------------
smooth_method='spline';
npoints=20000;

% - representative household model after wage change ----------------------
[grid_quantiles_rh, grid_change_rh]=income_quantiles_rh_v10(paths_rh,ss0_rh, cal_ell_z,cal_Z,distribution_ss0, Parameter, npoints, smooth_method);

% - model with dissipation shocks -----------------------------------------
[grid_quantiles, grid_change, grid_change_labor, grid_change_capital]=income_quantiles_v10(cal_ell_z,cal_Z,distribution_ss0, distribution_ssf, Parameter, npoints, smooth_method);

% - set colors and properties for the figures -----------------------------
colx=0.6;
labor_color=[0, 0.4470, 0.7410]*colx+(1-colx)*[1,1,1];
growth_color=[0.9500, 0.4250, 0.40]*colx+(1-colx)*[1,1,1];
cut5pct = 1;    %set to 1 cut everything below 5 pct
leftaxis = 1;   %set to 1 to have the same yaxis for the two panels

if cut5pct==1
    startpct= 'cut5pct';
else
    startpct= 'include5pct';    
end

if leftaxis==1
    plotaxis= 'sameaxis';
else
    plotaxis= 'diffaxis';    
end

if Parameter.onepanel==1
    figtype= 'onepanel';
else
    figtype= '';    
end

% - two-panel figure - %
if Parameter.onepanel==0
    
    figure_quantiles = figure('units','centimeters','outerposition',[10 10 16*1.2 12*1.2]);

    title('Change in income by percentile of the income distribution')
    set(gca,'ycolor','w');
    set(gca,'xcolor','w');

    axes('Position', [0.06 0.12 .65 .8]);
    hold on
    box on
    h1= area(grid_quantiles, grid_change_capital + grid_change_labor.*(grid_change_labor>0),'Facecolor',growth_color,'EdgeColor',growth_color);
    area(grid_quantiles, grid_change_capital,'Facecolor',growth_color,'EdgeColor',growth_color); % needed so that the mass of capital and labor add to total on the graph (see when they overlap, eg around 90)
    h2= area(grid_quantiles, grid_change_labor ,'Facecolor',labor_color,'EdgeColor',labor_color);
    h3= plot(grid_quantiles,grid_change,'k-', 'linewidth',1.2);
    h4= plot(grid_quantiles_rh,grid_change_rh,'r--', 'linewidth',1.2);

    if leftaxis==1
        ylim([-15 75])
    else
        ylim([-15 40])
    end

    if cut5pct==1
        xlim([5,100])
    end

    xlabel('income percentile')
    ylabel('percent change')

    legend([h3 h2 h1 h4], ...
        {'Total income growth', 'Part due to wage income',...
        'Part due to capital income' ,'Rep. household model'} ...
        , 'location', 'northwest');  
    legend('boxoff')

    axes('Position', [0.75 0.12 .2 .8],'Box','on');
    hold on
    box on
    h1= area(grid_quantiles, ...
        grid_change_capital+grid_change_labor.*(grid_change_labor>0) , ...
        'Facecolor',growth_color,'EdgeColor',growth_color);
    h2= area(grid_quantiles, grid_change_labor ,'Facecolor',...
        labor_color,'EdgeColor',labor_color);
    h3= plot(grid_quantiles,grid_change,'k-', 'linewidth',1.2);
    h4= plot(grid_quantiles_rh,grid_change_rh,'r--', 'linewidth',1.2);
    ylim([-15,75])


    xlim([99,100])
    xlabel('top tail')
    set(gcf,'color','w');
    
% - single panel figure - %    
else

    figure_quantiles = figure('units','centimeters','outerposition',[10 10 16*1.2 12*1.2]);

    xaxis = 1:length(grid_quantiles);
    
    title('Change in income by percentile of the income distribution')
    set(gca,'ycolor','w');
    set(gca,'xcolor','w');


    axes('Position', [0.06 0.22 .9 .7]);
    hold on
    box on
    h1= area(xaxis, grid_change_capital + grid_change_labor.*(grid_change_labor>0),'Facecolor',growth_color,'EdgeColor',growth_color);
    area(xaxis, grid_change_capital,'Facecolor',growth_color,'EdgeColor',growth_color); % needed so that the mass of capital and labor add to total on the graph (see when they overlap, eg around 90)
    h2= area(xaxis, grid_change_labor ,'Facecolor',labor_color,'EdgeColor',labor_color);
    h3= plot(xaxis,grid_change,'k-', 'linewidth',1.2);
    h4= plot(xaxis,grid_change_rh,'r--', 'linewidth',1.2);

    ylim([-15 60])
    xlim([1.99 16.01])

    set(gca,'XTick',1:length(xaxis),'XTickLabel',{'P0-P5','P5-P10','P10-P20','P20-P30','P30-P40','P40-P50','P50-P60','P60-P70','P70-P80' ...
        ,'P80-P90','P90-P95','P95-P97.5','P97.5-P99','P99-P99.5','P99.5-P99.9','P99.9-P99.99','P99.99-P100'}, ...
        'XTickLabelRotation',45)

    xlabel('income percentile')
    ylabel('percent change')

    legend([h3 h2 h1 h4], ...
        {'Total income growth', 'Part due to wage income',...
        'Part due to capital income' ,'Rep. household model'} ...
        , 'location', 'northwest');  
    legend('boxoff')
    set(gcf,'color','w');
      
end


if (binding==1)
    eval(['export_fig output/Figure_paper_quantiles_binding_' figtype '_' plotaxis '_' startpct '.pdf'])
else
    eval(['export_fig output/Figure_paper_quantiles_nobinding_' figtype '_' plotaxis '_' startpct '.pdf'])
end


%% - Comparative statics following change in alphas -----------------------

N=100; % number of grid points and placeholders 
alphas           = zeros(N,1);
r_k_comp         = zeros(N,1);
r_b_comp         = zeros(N,1);
r_bar_comp       = zeros(N,1);
r_dif_comp       = zeros(N,1);
r_w_comp         = zeros(N,1);
kappas_comp      = zeros(N,1);
mu_adj_e_comp    = zeros(N,1);
mu_adj_e_comp_p  = zeros(N,1);
mu_adj_h_comp    = zeros(N,1);
pareto_full_comp = zeros(N,1);
pareto_mu_comp   = zeros(N,1);
k_norm_comp      = zeros(N,1);
b_h_comp         = zeros(N,1);
b_e_comp         = zeros(N,1);
lambda_comp      = zeros(N,1);
chialpha_comp    = zeros(N,1);
rprem_adj_comp   = zeros(N,1);
adj_return_comp  = zeros(N,1);
k_hour_comp      = zeros(N,1);
y_hour_comp      = zeros(N,1);
k_to_y_comp      = zeros(N,1);
tfp_comp         = zeros(N,1);
w_ave_hour_comp  = zeros(N,1);
w_med_hour_comp  = zeros(N,1);

% - Loop through values and compute equilibrium -
for i = 1:N
    
    % - Set parameter value -
    aalpha=cal_alpha_0+(i-1)*(cal_alpha_f-cal_alpha_0)/(N-1);
    eeta=cal_eta_0;
    A_neutral=cal_A_0;
    Parameter.alpha = aalpha*eeta';
    alphas(i)= aalpha*eeta';
    
    if (i == 1) % Initialize
        Solution0(1) = log((Parameter.alpha/(1-Parameter.alpha)) ...
                        /(Parameter.delta+Parameter.rho+Parameter.sigma*Parameter.g)); % Capital stock (relative to wages)
        Solution0(2) = Parameter.govbonds;                         % Bond holdings of households (relative to wages)
        Solution0(3) = log(Parameter.rho+Parameter.sigma*Parameter.g);    % Return to bonds
    else % Use previous result
        Solution0 = Solution;
    end
    
    % - Solve system -
    Solution = fsolve(@(Solution)equations_v10(Solution,Parameter),Solution0,options);
    
    % - Compute variables of interest -
    lambda_comp(i)       = Parameter.lambda;
    chialpha_comp(i)     = Parameter.chi/Parameter.alpha;
    k_norm_comp(i)       = exp(Solution(1));
    b_h_comp(i)          = Solution(2);
    b_e_comp(i)          = Parameter.govbonds - b_h_comp(i);
    r_b_comp(i)          = exp(Solution(3))+Parameter.g;
    r_k_comp(i)          = (Parameter.alpha/(1-Parameter.alpha))*(r_b_comp(i)-Parameter.g)/k_norm_comp(i) - Parameter.delta;
    r_dif_comp(i)        = r_k_comp(i) - r_b_comp(i);
    kappas_comp(i)       = min(1/Parameter.lambda,(r_dif_comp(i))/Parameter.gamnu);
    r_bar_comp(i)        = kappas_comp(i) * r_k_comp(i) + (1-kappas_comp(i))*r_b_comp(i);
    mu_adj_e_comp(i)     = (1/Parameter.sigma)*(r_bar_comp(i) - Parameter.rho - Parameter.sigma*Parameter.g + 0.5 * (Parameter.sigma-1) * Parameter.gamnu * kappas_comp(i)^2)-0.5*kappas_comp(i)^2*nu^2;
    mu_adj_h_comp(i)     = (1/Parameter.sigma)*(r_b_comp(i) - Parameter.rho - Parameter.sigma*Parameter.g);
    pareto_full_comp(i)  =(1/(2*Parameter.p))*(mu_adj_e_comp(i) + sqrt(mu_adj_e_comp(i)^2 + 2 * Parameter.p*kappas_comp(i)^2*nu^2));
    mu_adj_e_comp_p(i)   = (1/Parameter.sigma)*(r_bar_comp(i) - Parameter.rho - Parameter.sigma*Parameter.g + 0.5 * (Parameter.sigma-1) * Parameter.gamnu * kappas_comp(i)^2)-0.5*kappas_comp(1)^2*nu^2;
    pareto_mu_comp(i)    =(1/(2*Parameter.p))*(mu_adj_e_comp_p(i) + sqrt(mu_adj_e_comp_p(i)^2 + 2 * Parameter.p*kappas_comp(1)^2*nu^2));
    rprem_adj_comp(i)    = r_dif_comp(i)/Parameter.gamnu;
    adj_return_comp(i)   =r_bar_comp(i)-Parameter.rho- Parameter.sigma*Parameter.g + 0.5 * (Parameter.sigma-1) * Parameter.gamnu * kappas_comp(i)^2;

    Lwtd=prod((cal_psi_L.*cal_ell_z./(eeta.*(1-aalpha))).^(eeta.*(1-aalpha)));
    k_hour_comp(i)=(A_neutral*Lwtd*Parameter.alpha^(1-Parameter.alpha)/(r_k_comp(i)+Parameter.delta))^(1/(1-Parameter.alpha));
    y_hour_comp(i)=A_neutral*((k_hour_comp(i)/Parameter.alpha)^Parameter.alpha)*Lwtd;
    w_ave_hour_comp(i)=mean(eeta.*(1-aalpha).*y_hour_comp(i)./cal_ell_z);
    w_med_hour_comp(i)=median(eeta.*(1-aalpha).*y_hour_comp(i)./cal_ell_z);
    k_to_y_comp(i)=k_hour_comp(i)/y_hour_comp(i);
    
    if i>1
    tfp_comp(i)=log(y_hour_comp(i)/y_hour_comp(i-1))-alphas(i-1)*log(k_hour_comp(i)/k_hour_comp(i-1))+tfp_comp(i-1);
    end
end

%% - Summary of comparative statics with respect to alpha -----------------

% - Plot aggregates -------------------------------------------------------
xaxis = alphas;
xlab = 'Capital share, $\alpha$';

% - Figure 1 -
figure('Name','Aggregates');
hold on
subplot(2,3,1)
plot(xaxis,100*r_k_comp,'-r',xaxis,100*r_b_comp,'-k','LineWidth',1.5)
text(xaxis(20),100*r_k_comp(20)*0.99,'Capital','Color','red')
text(xaxis(1)*1.01,100*r_b_comp(1)*1.03,'Bonds')
axis([-inf inf 100*r_b_comp(1)*0.95 inf]);
title('Returns (in p.p)')
xlabel(xlab)

subplot(2,3,2)
plot(xaxis,100*r_dif_comp,'-k','LineWidth',1.5)
title('Difference $r_K-r_B$ (in p.p)')
axis([-Inf Inf -Inf Inf]);
xlabel(xlab)

subplot(2,3,3)
plot(xaxis,w_med_hour_comp,'-r', xaxis,w_ave_hour_comp,'-k','LineWidth',1.5)
text(xaxis(1)*1.01,w_med_hour_comp(1)*0.95,'Median','Color','red')
text(xaxis(1)*1.01,w_med_hour_comp(1)*0.91,'wage','Color','red')
text(xaxis(1)*1.01,w_ave_hour_comp(1)*1.03,'Mean wage')
title('Wage per hour (in dollars)')
yticks([18 20 22 24 26 28])
yticklabels({'\$18','\$20','\$22','\$24','\$26','\$28'})
axis([-Inf Inf 18 28]);
xlabel(xlab)

subplot(2,3,4)
plot(xaxis,y_hour_comp,'-k','LineWidth',1.5)
title('Output per hour (in dollars)')
yticks([38 40 42 44 46])
yticklabels({'\$38','\$40','\$42','\$44','\$46'})
axis([-Inf Inf 38 46]);
xlabel(xlab)

subplot(2,3,5)
plot(xaxis,k_to_y_comp,'-k','LineWidth',1.5)
title('Capital-output ratio')
axis([-Inf Inf -Inf Inf]);
xlabel(xlab)

subplot(2,3,6)
plot(xaxis,100*tfp_comp,'-k','LineWidth',1.5)
title('TFP gains (percent change)')
yticks([0 0.5 1 1.5])
yticklabels({'0','0.5\%','1\%','1.5\%'})
axis([-Inf Inf -Inf Inf]);
xlabel(xlab)

set(gcf,'units', 'centimeters', 'position',[0,0,17,10]);
set(gcf, 'Color', 'w');

if (binding==1)
export_fig 'output/Figure_paper_alpha_aggregates_binding'.pdf 
else
export_fig 'output/Figure_paper_alpha_aggregates_notbinding'.pdf 
end

% - Plot outcomes related to wealth inequality ----------------------------
figure('Name','Wealth inequality');

hold on
subplot(1,2,1)
plot(xaxis,100*adj_return_comp,'-k','LineWidth',1.5)
title('Return gap $r_W-\rho-\sigma g$ (in p.p)')
axis([-Inf Inf -Inf Inf]);
xlabel(xlab)

subplot(1,2,2)
plot(xaxis,pareto_full_comp,'-k',xaxis, pareto_mu_comp,'--r','LineWidth',1.5)
hl=legend('Full effect', 'Drift effect', 'Location', 'Southeast');
title('Tail index, $1/\zeta$')
axis([-Inf Inf -Inf Inf]);
xlabel(xlab)

set(gcf,'units', 'centimeters', 'position',[0,0,17,7]);
set(gcf, 'Color', 'w');
  
if (binding==1)
export_fig 'output/Figure_paper_alpha_inequality_binding'.pdf 
else
export_fig 'output/Figure_paper_alpha_inequality_notbinding'.pdf 
end

