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

% Calibrated Parameters


%Markuptarget = 1.05;                      % our target for aggregate markup

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


%Markuptarget = 1.10;                      % our target for aggregate markup

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


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 


%Markuptarget = 1.20;                      % our target for aggregate markup

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


%Markuptarget = 1.25;                      % our target for aggregate markup

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


%Markuptarget = 1.30;                      % our target for aggregate markup

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


%Markuptarget = 1.35;                      % our target for aggregate markup

%p.xi         = 2.88859963714985;          % Pareto tail productivity
%p.sigma      = 5.66479780771211;          % 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     = 1; 

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

else
    
    options.Display = 'none';

end


% Quadrature
 
nz         = 200;                         % 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; 


if printr
    
loglog(z, [1 - cumsum(w), z.^(-p.xi)]);
ylim([1e-8, 1])

fprintf('Mean z (True, Quadrature)     = %9.4f %9.4f\n',  [p.xi/(p.xi - 1),               w'*z]);
fprintf('Var z (True, Quadrature)      = %9.4f %9.4f\n',  [p.xi/(p.xi - 1)^2/(p.xi - 2),  w'*(z - w'*z).^2]);
fprintf('E z^9 (True, Quadrature)      = %9.4f %9.4f\n',  [(p.xi/(p.xi - 5))^(1/5),       (w'*z.^5)^(1/5)]);

end


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

if printr
    
% 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');


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');

end

moment_data    = [Markuptarget; 0.45; 0.229; 0.573; 0.739; 0.162]; 

    
moment_model       = zeros(size(moment_data)); 

moment_model(1)    = Mu; 
moment_model(2)    = B/Y;


data               = sortrows([w, py.*y], 2); 

data               = data(data(:, 2) > 1e-16, :); 


data(:,1)          = data(:,1)/sum(data(:,1)); 

cumF               = cumsum(data(:,1)); 



top1               = cumF >= 0.99; 
top5               = cumF >= 0.95; 
top10              = cumF >= 0.90; 

moment_model(3)    = data(top1,  1)'* data(top1,  2)/(data(:,1)'*data(:,2));
moment_model(4)    = data(top5,  1)'* data(top5,  2)/(data(:,1)'*data(:,2));
moment_model(5)    = data(top10, 1)'* data(top10, 2)/(data(:,1)'*data(:,2));


keep               = py.*y > 1e-8; 

bhat               = lscov([ones(size(y(keep))), log(py(keep).*y(keep))], 1./mu(keep) + log(1 - 1./mu(keep)));

moment_model(6)    = bhat(2);


if printr
    

fprintf('\n');

fprintf('Moments used in Calibration: Left = Model, Right = Data  \n');   
fprintf('\n');
fprintf('\n');
fprintf('Aggregate Markup                           = %9.3f %9.3f\n',  [moment_model(1),    moment_data(1)]);
fprintf('Share Intermediates in Sales               = %9.3f %9.3f\n',  [moment_model(2),    moment_data(2)]);
fprintf('\n')
fprintf('Sales share largest 1   percent firms      = %9.3f %9.3f\n',  [moment_model(3),    moment_data(3)]);
fprintf('Sales share largest 5   percent firms      = %9.3f %9.3f\n',  [moment_model(4),    moment_data(4)]);
fprintf('Sales share largest 10  percent firms      = %9.3f %9.3f\n',  [moment_model(5),    moment_data(5)]);
fprintf('\n');
fprintf('Regression coefficient                     = %9.3f %9.3f\n',  [moment_model(6),    moment_data(6)]);
fprintf('\n');
end

weights        = zeros(numel(moment_model), 1); 

weights([1; 2; 4])  = 1; 


weights        = weights/sum(weights);


err_mom        = (moment_model - moment_data)./(1 + moment_data);
err_mom        = (weights'*err_mom.^2).^(1/2);

if printr
        
fprintf('\n');
fprintf('\n');   
fprintf('Error in Moments                           = %9.6f \n', err_mom);
fprintf('\n');
fprintf('\n');   

end


%start_planner
%start_subsidy

 p.xie        = 3;        % approximately triples N if Mu=1.15
       
 Xnew         = fsolve('findequilibrium', [D; Y; N], options, w, z, p, 'market', 'new');

 [~, qnew, munew, ~, Znew, Wnew, Lpnew, Knew, Bnew, Unew, Uqnew, Cnew, Munew] = findequilibrium(Xnew, w, z, p, 'market', 'new');

 
 % New Steady State Allocations

Dnew          = Xnew(1); 
Ynew          = Xnew(2); 
Nnew          = Xnew(3); 


Gnew          = (1 + p.xie)*(1 + p.xis - 1/Munew);
Qnew          = 1/(1 - p.beta*(1 - p.varphi))*Gnew*Ynew/Nnew; 

Lnew          = Lpnew + p.F*p.varphi*Nnew; 
Rnew          = 1/p.beta - 1 + p.delta; 


ynew          = qnew*Ynew; 

Pnew          = (Nnew*w'*(Uqnew.*qnew))^(-1);       % Demand index so we can compute prices


% Individual Choices

pynew          = Uqnew*Pnew;                      % producer price
y1             = qnew*Ynew; 
lnew           = pynew.*ynew/Ynew.*Munew./munew.*Lpnew; 






fprintf('Various Checks of the Code\n');
fprintf('\n');
fprintf('Error in Kimball aggregator             = %9.2e \n',  Nnew*w'*Unew - 1);
fprintf('Error in Consumption-Leisure Choice     = %9.2e \n',  p.psi - Wnew/(Cnew*Lnew^p.nu));
fprintf('Error in Agg Prod Function              = %9.2e \n',  Ynew - Znew*(p.phi^(1/p.theta)*(Knew^p.alpha*Lpnew^(1 - p.alpha))^(1 - 1/p.theta) + (1 - p.phi)^(1/p.theta)*Bnew^(1 - 1/p.theta))^(p.theta/(p.theta - 1)));

fprintf('\n');
fprintf('\n');

fprintf('\n');
fprintf('\n');
fprintf('Compare to original allocations: 1. Old, 2. New 3. Change \n');
fprintf('\n');
fprintf('\n');
fprintf('Measure firms, N             = %9.3f %9.3f %9.3f\n',  [N,       Nnew,        Nnew/N   - 1]);
fprintf('Output, Y                    = %9.3f %9.3f %9.3f\n',  [Y,       Ynew,        Ynew/Y   - 1]);
fprintf('Consumption, C               = %9.3f %9.3f %9.3f\n',  [C,       Cnew,        Cnew/C   - 1]);
fprintf('Capital, K                   = %9.3f %9.3f %9.3f\n',  [K,       Knew,        Knew/K   - 1]);
fprintf('Employment, L                = %9.3f %9.3f %9.3f\n',  [L,       Lnew,        Lnew/L   - 1]);
fprintf('Employment (production), Lp  = %9.3f %9.3f %9.3f\n',  [Lp,      Lpnew,       Lpnew/Lp - 1]);
fprintf('TFP, Z                       = %9.3f %9.3f %9.3f\n',  [Z,       Znew,        Znew/Z   - 1]);
fprintf('True Markup, Mu              = %9.3f %9.3f %9.3f\n',  [Mu,      Munew,       Munew/Mu   - 1]);
    
fprintf('\n');
fprintf('\n');



figure(1)

subplot(1,2,1)
semilogx(z,mu,z,munew)
axis([1 100 1 1.6])
xlabel('productivity, $z$')
ylabel('markups $\mu(z)$')
legend('$N=1$','$N=3$','Location','SouthEast')
subplot(1,2,2)
semilogx(z,lnew./l)
axis([1 100 0 1])
xlabel('productivity, $z$')
ylabel('ratio $l_{new}(z)/l(z)$')


print -dpng triple_entry_figure
matlab2tikz('triple_entry_figure.tikz')


 