% Figure 4: Optimal Entry Subsidy, Mu = 1.15
%
%   calculates welfare for each entry subsidy 
% 
%   chi_e in (0.1:0.01:0.7);
%
%   this grid is just for illustrative purposes
%   optimal entry calculations reported in tables etc find optimal chi_e 
%   by golden search

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

subsidies = (0.1:0.01:0.7); % grid of entry subsidies
 
welfare   = zeros(size(subsidies));

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

% Calibrated Parameters

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

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


%start_planner
%start_subsidy

for tt=1:numel(subsidies)
    
 p.xie        = subsidies(tt);        
 
 Xnew         = fsolve('findequilibrium', [D*1.1; Y*1.5; N*1.2], 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');

% Now start doing transition dynamics

% See if we can use dynare: parameterize Z(N) with chebyshev polynomials

global cfun fspace; 

npts      = 11; 

optset('chebnode', 'nodetype', 2);              % Lobatto nodes (include endpoints)

fspace    = fundefn('cheb', npts, 1, 2); 

Ngrid     = funnode(fspace); 


Zgrid     = zeros(npts, 1); 
Mugrid    = zeros(npts, 1); 

options.Display = 'none';

for i = 1 : npts

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

    [~, ~, ~, ~, Zgrid(i), ~, ~, ~, ~, ~, ~, ~, Mugrid(i)]  = findequilibrium(D, w, z, p, 'market', 'old', Ngrid(i));
    
end


Ggrid  = (1 + p.xie)*(1 + p.xis - 1./Mugrid);               % flow gains from varieties

cfun   = funfitxy(fspace, Ngrid, [Zgrid, Ggrid, Mugrid]);

xpar   = [p.beta; p.varphi; p.delta; p.psi; p.F; p.nu; p.alpha; p.theta; p.phi];                              save xpar xpar; 

x1_ss  = [K;    N];                                                                                           save x1_ss x1_ss; 
x2_ss  = [Ynew; Cnew; Bnew; Lpnew; Knew; Nnew; Wnew; Rnew; Qnew; Znew; Gnew; Znew/Munew; Munew];              save x2_ss x2_ss; 


dynare transition noclearall;

% Cleanup

filename = 'transition';

eval(['delete ' filename '*.m']);
eval(['delete ' filename '*.mat']);
eval(['delete ' filename '*.swp']);
eval(['delete ' filename '*.log']);
eval(['delete ' filename '*.eps']);
eval(['delete ' filename '*.pdf']);
eval(['delete ' filename '*.fig']);
eval(['delete ' filename '*.jnl']);


% Populate first period values from original steady state

Yt(1)  = Y; 
Ct(1)  = C; 
Bt(1)  = B; 
Lpt(1) = Lp; 
Kt(1)  = K; 
Nt(1)  = N; 
Wt(1)  = W; 
Rt(1)  = R; 
Zt(1)  = Z; 
Mut(1) = Mu; 

Lt     = [L; Lpt(2:end) + p.F*(Nt(2:end) - (1 - p.varphi)*Nt(1 : end-1))]; 
Xt     = [p.delta*K; Kt(2:end) - (1 - p.delta)*Kt(1 : end - 1)]; 


T    = length(Ct) - 1; 

Wold = 1/(1 - p.beta)*(log(C)  - p.psi/(1 + p.nu)*L^(1 + p.nu)); 

Wnew = (p.beta.^(0 : 1 : T - 1))*(log(Ct(2 : end)) - p.psi/(1 + p.nu)*Lt(2:end).^(1 + p.nu)) + p.beta^T/(1 - p.beta)*(log(Ct(end)') - p.psi/(1 + p.nu)*Lt(end)'.^(1 + p.nu));

dW   = (exp((1 - p.beta)*(Wnew - Wold)) - 1)*100;

fprintf('\n')
fprintf('Subsidy                             = %9.5f \n',  p.xie);
fprintf('\n')
fprintf('Welfare Gains                       = %9.5f \n',  dW);
fprintf('\n')



welfare(tt) = dW;




end

figure(4)
hold on
plot(subsidies,welfare)
plot(subsidies,zeros(size(subsidies)),'k--')
hold off
grid on 
box on
xlim([0    ,0.75]) 
ylim([-0.15,0.60])
xlabel('entry subsidy, $\chi_e$','Interpreter','latex')
ylabel('welfare, CEV, \%','Interpreter','latex')

print -dpng figures/Figure4
 