% Figure 2: Equilibrium and Planner Allocations
%
%   solves model for calibration corresponding to Mu = 1.15 
%   
%   compares equilibrium q(z),l(z) to planner's q*(z),l*(z)
%

clear;
close all;
clc;

set(groot, 'DefaultAxesLineWidth', 1.5);
set(groot, 'DefaultLineLineWidth', 3);
set(groot, 'DefaultAxesTickLabelInterpreter','latex'); 
set(groot, 'DefaultLegendInterpreter','latex');
set(groot, 'DefaultAxesFontSize', 22);

optset('bisect','tol', 1e-18); 

mshare       = 0.45;                      % share of intermediate inputs in total sales

Markuptarget = 1.15;                      % our target for aggregate markup

p.xi         =  6.8409724072258;          % Pareto tail productivity
p.sigma      = 10.8610835311753;          % demand elasticity: match level of markups 

% Assigned Parameters

p.epsi     =  0.162*p.sigma; 
p.phi      =  0.4;                        % doesn't matter since reset below to match material share  

p.beta     = 0.96;                        % period is 1 year
p.varphi   = 0.04;                        % exit rate 
p.delta    = 0.06;                        % depreciation rate

p.r        = 1/p.beta - 1; 
p.nu       = 1;                           % inverse labor supply elasticity
p.alpha    = 1/3;                         % capital share in value added
p.theta    = 0.5;                         % elasticity of substitution between VA and intermediates

p.xie      = 0;                           % entry subsidy
p.xis      = 0;                           % sales subsidy
p.taus     = 0; 

% Initial Aggregates (Normalization)

Y          = 1; 
N          = 1; 
R          = p.r + p.delta; 

printr     = 0; 

if printr == 1
    
    options.Display = 'iter';

else
    
    options.Display = 'none';

end


% Quadrature
 
nz         = 5000;                         % Gauss-Legendre quadrature for exponential r.v.

zmax       = -1/p.xi*log(1e-22);

[z, w]     = qnwlege(nz, 0, zmax); 
w          = p.xi.*w.*exp(-p.xi.*z);
w          = w./sum(w); 

z          = exp(z); 

p.w        = w; 
p.z        = z; 

%%%%% 1. Solve D and firm's choices for a given N %%%%%% 

%%%% Start solving

D              = fsolve('findequilibrium', 1, options, w, z, p, 'market', 'old');

[~, ~, ~, Omega, ~,~, ~, ~, ~, ~, ~, ~, Mu]  = findequilibrium(D, w, z, p, 'market', 'old');

p.phi          = 1 - mshare*Omega^(1 - p.theta)*Mu;                                       % choose to match materials share in sales = mshare 

[~, q, mu, Omega, Z, W, Lp, K, B, U, Uq, C, Mu]  = findequilibrium(D, w, z, p, 'market', 'old');

p.F            = Y/N/W/(1/p.beta - 1 + p.varphi)*(1 - 1/Mu);     % uses normalization N/Y = 1 to back out F

L              = Lp + p.F*p.varphi*N;

p.psi          = W/(C*L^p.nu); 

P              = (N*w'*(Uq.*q))^(-1);       % Demand index so we can compute prices


% Individual Choices

py             = Uq*P;                      % producer price
y              = q*Y; 
l              = py.*y/Y.*Mu./mu.*Lp; 
b              = py.*y/Y.*Mu./mu.*B; 
k              = py.*y/Y.*Mu./mu.*K; 
mc             = Omega./z; 


% Report Some Aggregate Statistics
    
% fprintf('\n');
% fprintf('\n')
% fprintf('Entry cost, kappa                       = %9.3f \n',  p.F);
% fprintf('Disutility from work, psi               = %9.3f \n',  p.psi);
% fprintf('\n');
% fprintf('\n');
% fprintf('Measure firms, N                        = %9.3f \n',  N);
% fprintf('Capital-Output rato, K/Y                = %9.3f \n',  K/Y);
% fprintf('N-Output rato, N/Y                      = %9.3f \n',  N/Y);
% fprintf('Wage rate                               = %9.3f \n',  W);
% fprintf('Output, Y                               = %9.3f \n',  Y);
% fprintf('Consumption, C                          = %9.3f \n',  C);
% fprintf('Materials, B                            = %9.3f \n',  B);
% fprintf('GDP, Y - B                              = %9.3f \n',  Y - B);
% fprintf('Investment, X                           = %9.3f \n',  p.delta*K);
% fprintf('Employment, L                           = %9.3f \n',  L);
% fprintf('Employment (production), Lp             = %9.3f \n',  Lp);
% fprintf('Labor share (production)                = %9.3f \n',  W*Lp/Y);
% fprintf('Agg Labor share                         = %9.3f \n',  W*L/Y);
% fprintf('Variable Cost share                     = %9.3f \n',  (W*Lp + B)/Y);
% fprintf('\n');
% fprintf('Investment to GDP                       = %9.4f \n',  p.delta*K/(Y - B));
% fprintf('Materials Share in Sales                = %9.4f \n',  B/Y);
% fprintf('Profits to GDP                          = %9.4f \n',  (Y - W*Lp - R*K - B)/(Y - B));
% 
% fprintf('\n');

% weight = w.*l/sum(w.*l); 

% fprintf('Aggregate Markup (cost-weighted)        = %9.3f \n',  Mu);
% fprintf('\n');
% 
% fprintf('10 p.c Markup (cost-weighted)           = %9.3f \n',  wprctile(mu, 10, weight));
% fprintf('25 p.c Markup (cost-weighted)           = %9.3f \n',  wprctile(mu, 25, weight));
% fprintf('50 p.c Markup (cost-weighted)           = %9.3f \n',  wprctile(mu, 50, weight));
% fprintf('75 p.c Markup (cost-weighted)           = %9.3f \n',  wprctile(mu, 75, weight));
% fprintf('90 p.c Markup (cost-weighted)           = %9.3f \n',  wprctile(mu, 90, weight));
% fprintf('95 p.c Markup (cost-weighted)           = %9.3f \n',  wprctile(mu, 95, weight));
% fprintf('99 p.c Markup (cost-weighted)           = %9.3f \n',  wprctile(mu, 99, weight));
% 
% fprintf('\n');
% fprintf('\n');
% 
% fprintf('Various Checks of the Code\n');
% fprintf('\n');
% fprintf('Error in Kimball aggregator             = %9.2e \n',  N*w'*U - 1);
% fprintf('Zero Profits Final Goods Firm           = %9.2e \n',  N*w'*(py.*y) - Y);
% fprintf('Error in Consumption-Leisure Choice     = %9.2e \n',  p.psi - W/(C*L^p.nu));
% fprintf('Error in Capital Resource constraint    = %9.2e \n',  N*w'*k - K);
% fprintf('Error in Labor Resource constraint      = %9.2e \n',  N*w'*l - Lp);
% fprintf('Error in Materials  Res constraint      = %9.2e \n',  N*w'*b - B);
% fprintf('1 - Markup x Marginal Cost Aggreg       = %9.2e \n',  1 - Mu*Omega/Z);
% fprintf('Error in Agg Prod Function              = %9.2e \n',  Y - Z*(p.phi^(1/p.theta)*(K^p.alpha*Lp^(1 - p.alpha))^(1 - 1/p.theta) + (1 - p.phi)^(1/p.theta)*B^(1 - 1/p.theta))^(p.theta/(p.theta - 1)));
% 
% fprintf('Error in Markups                        = %9.2e \n',  norm(py(q > 0)./mc(q > 0) - mu(q > 0)));
% 
% fprintf('\n');
% fprintf('\n');


Dp       = fsolve('findequilibrium', D, options, w, z, p, 'planner', 'old');

[~, qp, ~, ~, Zp, ~,~, ~, ~, Up]  = findequilibrium(Dp, w, z, p, 'planner', 'old');

lp = qp./z;

lbar = Lp/(N*w'*lp); % planner's allocation of l(z) = market Lp at given N

lp   = lp*lbar;       

% fprintf('\n');
% fprintf('Error Kimball aggregator planner             = %9.2e \n',  N*w'*Up - 1);
% fprintf('\n');
% fprintf('Losses from Misallocation: Z, p.p.           = %9.2f \n',  (Zp/Z-1)*100);
% fprintf('\n');
% fprintf('\n');


figure(2)
subplot(1,2,1)
title('relative size')
loglog(z,[q,qp])
xlabel('$\log z$','interpreter','latex')
ylabel('$\log q(z)$','interpreter','latex') 
xlim([10^0 ,10^2])
ylim([10^-2,10^6])
legend('equilibrium','planner','Location','SouthEast')

subplot(1,2,2)
title('employment')
loglog(z,[l,lp])
xlabel('$\log z$','interpreter','latex')
ylabel('$\log l(z)$','interpreter','latex')
xlim([10^0 ,10^2])
ylim([10^-2,10^6])
legend('equilibrium','planner','Location','SouthEast')

print -dpng figures/Figure2




