% Code to calclaute markups costs
clear all
close all
clc

addpath('\Model\Tables\CompEcon2017\CEtools')

% Save outputs
outputs = zeros(3,1,2);
outputs_distr = zeros(5,12,2);
welfare = zeros(7,12,2);

% Calibrated Parameters
pre_param = xlsread('output\out_full_pre.xlsx');
post_param = xlsread('output\out_full_post.xlsx');
pre_data=xlsread('\Inputs for model.xlsx','pre_main','d2:i13');
post_data=xlsread('\Inputs for model.xlsx','post_main','d2:i13');

%Loop over periods
for pre = 1:2

if pre == 1 
input = pre_data;
param = pre_param;
else
input = post_data;
param=post_param;
end

%loop over specifciation - as in calibration
for option = 1:2
% Output Selection
if input(option,1)>0

% Select report - T3 is misallcation cost main metric. T4 is broader markup cost
report       = 'T4';          % 'T1', 'T3', 'T4' 
markup       = string(option);        % '1.05', '1.15', '1.25', '1.35'
%slect comaprison point for braoder cost - efficient is total cost.
experiment   = 'efficient';   % 'efficient', 'uniform', 'sizedep', 'entry'

% Options

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

options.Display = 'none';

Markuptarget = input(option,3);                      % our target for aggregate markup

p.xi  = param(1,option);
p.sigma = param(2,option);

% Assigned Parameters


mshare     = input(option,6);                        % share of intermediate inputs in total sales

p.epsi     =  input(option,1)*p.sigma;

p.phi      =  mshare;                        % 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    = input(option,5)/(input(option,5)+input(option,4));                         % 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; 

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

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

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

% compute value added misallocation

GDPp  = p.phi^(1/(p.theta - 1))*Zp/(1 - (1 - p.phi)*Zp^(p.theta - 1))^(1/(p.theta - 1))*K^p.alpha*Lp^(1 - p.alpha);
Yp    = p.phi^(1/(p.theta - 1))*Zp/(1 - (1 - p.phi)*Zp^(p.theta - 1))^(p.theta/(p.theta - 1))*K^p.alpha*Lp^(1 - p.alpha);
Bp    = (1 - p.phi)*Zp.^(p.theta - 1)*Yp; 

GDP   = p.phi^(1/(p.theta - 1))*(1 - (1 - p.phi)*Z^(p.theta - 1)*Mu^(-p.theta))*Z/(1 - (1 - p.phi)*Z^(p.theta - 1)*Mu^(1 - p.theta))^(p.theta/(p.theta - 1))*K^p.alpha*Lp^(1 - p.alpha);
GDP1  = p.phi^(1/(p.theta - 1))*(1 - (1 - p.phi)*Z^(p.theta - 1)              )*Z/(1 - (1 - p.phi)*Z^(p.theta - 1)                 )^(p.theta/(p.theta - 1))*K^p.alpha*Lp^(1 - p.alpha);

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

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

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

profits            = (P*Uq.*q - Omega*(q./z));

transfer           = (U - Uq.*q)*P;

net_profits        = profits+transfer;

% Output 

switch report

    case 'T1'

        symbol      = categorical({sprintf('%c ', 946);sprintf('%c ', 948);sprintf('%c ', 966);sprintf('%c ', 945);sprintf('%c ', 957);sprintf('%c ', 952)});
        description = categorical({'discount factor';'depreciation rate';'exit rate';'elasticity of value-added to capital';'elasticity of labor supply';'elasticity of substitution between value-added and materials'});
        values      = [p.beta; p.delta; p.varphi; p.alpha; p.nu; p.theta];
        
        fprintf('<strong> Table 1: Parameterization, Benchmark Model </strong> \n'); 

        fprintf('\n');
        
        fprintf([' Aggregate Markup at ', num2str(round(Mu,2)), '\n'])
        
        fprintf('\n');
        
        fprintf(' A. Assigned Parameter \n');  
        
        fprintf('\n');
        
        disp(table(symbol, description, round(values, 2), 'VariableNames', {'Symbol','Variable Description','Value'}))
        
        fprintf('\n');

        fprintf(' B. Calibrated Parameter \n');
        
        fprintf('\n');
        
        symbol      = categorical({'M';'-';'-';'b';sprintf('%c ', 958);sprintf('%c ', 963);[sprintf('%c', 949),'/',sprintf('%c ', 963)];sprintf('%c ', 981)});
        description = categorical({'aggregate markup';'top 5% sales share';'materials share';'regression coefficient';'Pareto tail';'demand elasticity';'super-elasticity';'weight on valued-added'});
        values      = [Mu; moment_model(4); B; moment_model(6); p.xi; p.sigma; moment_model(6); p.phi];
        
        disp(table(symbol, description, round(values, 2), 'VariableNames', {'Symbol','Variable Description','Value'}))
        
        fprintf('\n');

    case 'T3'

        fprintf('<strong> Table 3: Markup Disperion and Productivity Losses, Benchmark Model </strong> \n'); 

        fprintf('\n');
        
        fprintf([' Aggregate Markup at ', num2str(round(Mu,2)), '\n'])

        fprintf('\n');
        
        description = categorical({'25th percentile';'50th percentile';'75th percentile';'90th percentile';'99th percentile'});
        values      = [wprctile(mu, 25, weight); wprctile(mu, 50, weight); wprctile(mu, 75, weight); wprctile(mu, 90, weight); wprctile(mu, 99, weight)];
        values_distr = values;
        fprintf(' Cost-weighted distribution of markups: \n');  
        
        fprintf('\n');
        
        disp(table(description, round(values, 2), 'VariableNames', {'Variable Description','Value'}))
        
        fprintf('\n');
        
        description = categorical({'gross output ';'value added';'value added, M = 1'});
        values      = [-(Z/Zp - 1)*100; -(GDP/GDPp  - 1)*100; -(GDP1/GDPp - 1)*100];
        
        fprintf(' Aggregate productivity loss in percent: \n');  
        
        fprintf('\n');
        
        disp(table(description, round(values, 2), 'VariableNames', {'Variable Description','Value'}))
        
        fprintf('\n');

outputs(:,option,pre) = values;
outputs_distr(:, option,pre) = values_distr;
    case 'T4'

        switch experiment

            case {'uniform', 'entry'}

                start_subsidy

            case {'efficient', 'sizedep'}

                start_planner

        end

    cost_vars = [round((Ynew/Y   - 1)*100, 1), round((Cnew/C   - 1)*100, 1), round((Lnew/L   - 1)*100, 1), round((Nnew/N   - 1)*100, 1), round((Knew/K   - 1)*100, 1), round((Znew/Z   - 1)*100, 1), round(dWs, 2)];
    welfare(:,option,pre) = cost_vars;

end

end
end
end

%% Save outputs - elect based on which table and experiment you are running

% Misallcoation cost
%writematrix(outputs(:,:,1), '\cost_full_pre.xlsx')
%writematrix(outputs(:,:,2), 'G\cost_full_post.xlsx')

% Markup distribution
%writematrix(outputs_distr(:,:,1), 'G:\Research\JonathanH\EMX 2022 Replication Package\distr_full_pre.xlsx')
%writematrix(outputs_distr(:,:,2), 'G:\Research\JonathanH\EMX 2022 Replication Package\distr_full_post.xlsx')

% Total cost
%writematrix(welfare(:,:,1), 'G:\Research\JonathanH\EMX 2022 Replication Package\welfare_pre.xlsx')
%writematrix(welfare(:,:,2), 'G:\Research\JonathanH\EMX 2022 Replication Package\welfare_post.xlsx')

% Cost with size depdnet subsidu
%writematrix(welfare(:,:,1), 'G:\Research\JonathanH\EMX 2022 Replication Package\size_welfare_pre.xlsx')
%writematrix(welfare(:,:,2), 'G:\Research\JonathanH\EMX 2022 Replication Package\size_welfare_post.xlsx')

% Cost with uniform subsidy
%writematrix(welfare(:,:,1), 'G:\Research\JonathanH\EMX 2022 Replication Package\uni_welfare_pre.xlsx')
%writematrix(welfare(:,:,2), 'G:\Research\JonathanH\EMX 2022 Replication Package\uni_welfare_post.xlsx')

