%% This code generates Table 1

% For user-cost markups use mu_sum(~, ~, 1)
% For account profit markups use mu_sum(~, ~, 2)
% For production function markups use mu_sum(~, ~, 3)

clc;
clear;

Path_FPA          = [cd, '/data_clean', '/FPA', '/FPA_sum.mat'];
Path_IO           = [cd, '/data_clean', '/IO_table', '/IO_sum.mat'];
Path_markup       = [cd, '/data_clean', '/markup', '/markup_sum.mat'];
Path_productivity = [cd, '/data_clean', '/productivity', '/MFP.mat'];


load(Path_FPA);
load(Path_IO);
load(Path_markup);
load(Path_productivity);

addpath('functions');


global N mu theta sigma_c sigma varepsilon omega_tilde alpha_tilde beta_tilde gamma varphi

%% Baseline calibration

% preference parameters
gamma          = 1;
varphi         = 2;

% number of industries

N              = 66;

% transforming the data from monthly to quarterly and calculating the
% degree of price flexibility of each industry
theta          = 1 - exp(-3*FPA);

% Calibrate the IO Table

IO             = IO_matrix_sum;
Final          = IO(:,93,:);
IO(isnan(IO))  = 0;
IO(IO < 0)     = 0;

% Labor as single input
grosscost      = squeeze(sum(IO(1:66,1:66,:)))+squeeze(IO(75,1:66,:));
labor          = squeeze(IO(75,1:66,:))./grosscost;
IO(67:end,:,:) = [];
IO(:,67:end,:) = [];


Omega_sum      = zeros(66,66);
beta_sum       = zeros(66,1);

for year = 1997:2015

    temp         = IO(:,:,year-1996);
    temp         = temp';
    Omega_sum    = Omega_sum + bsxfun(@rdivide, temp, grosscost(:,year-1996));
    beta_tilde   = Final(1:66,year-1996);
    beta_tilde(beta_tilde<0) = 0;
    beta_tilde   = beta_tilde/sum(beta_tilde);
    beta_sum     = beta_sum + beta_tilde;
end

omega_tilde  = Omega_sum/19;
alpha_tilde  = mean(labor(:,1:19),2);
beta_tilde   = beta_sum/19;
mu           = mean(mu_sum(:,:,1),2);

varepsilon   = 6;
sigma        = 0.5 * ones(N,1);
sigma_c      = 1;

% %% Calculating the weights corresponding to various policies
% 
% [phi_star,pi_star,~,~,phi_og,~,~,~,phi_dc, phi_dc_eff]  = optimalpolicy(mu, theta, sigma_c, sigma, varepsilon, omega_tilde, alpha_tilde, beta_tilde, gamma, varphi);
% 
% phi_dc         = phi_dc./sum(phi_dc);
% phi_dc_eff     = phi_dc_eff./sum(phi_dc_eff);
% phi_star       = phi_star./sum(phi_star);
% phi_og         = phi_og./sum(phi_og);
% 
% phi_dc_sec     = sector_aggregate_flexwage(phi_dc);
% phi_dc_eff_sec = sector_aggregate_flexwage(phi_dc_eff);
% phi_star_sec   = sector_aggregate_flexwage(phi_star);
% phi_og_sec     = sector_aggregate_flexwage(phi_og);

% bar_dataset       = [phi_dc_sec, phi_dc_eff_sec,phi_og_sec, phi_star_sec];
% secname = {   
%     'AGRICULTURE, FORESTRY, FISHING, AND HUNTING';
%     'MINING';
%     'UTILITIES';
%     'CONSTRUCTION';
%     'DURABLE GOODS';
%     'NONDURABLE GOODS';
%     'WHOLESALE TRADE';
%     'RETAIL TRADE';
%     'TRANSPORTATION AND WAREHOUSING, EXCLUDING POSTAL SERVICE';
%     'INFORMATION';
%     'FINANCE AND INSURANCE';
%     'REAL ESTATE AND RENTAL AND LEASING';
%     'PROFESSIONAL AND TECHNICAL SERVICES';
%     'MANAGEMENT OF COMPANIES AND ENTREPRISES';
%     'ADMINISTRATIVE AND WASTE SERVICES';
%     'EDUCATIONAL SERVICES';
%     'HEALTH CARE AND SOCIAL ASSISTANCE';
%     'ARTS, ENTERTAINMENT, AND RECREATION';
%     'ACCOMMODATION AND FOOD SERVICES';
%     'OTHER SERVICES, EXCEPT GOVERNMENT'};
% bar_dataset = flipud(bar_dataset);
% secname     = flipud(secname);
% 
% % figure;
% % bar(bar_dataset);
% % set(gca, 'XTickLabel', secname, 'XTick', 1:numel(secname));
% % legend('Optimal Policy', 'CPI Targetting');
% % title('Optimal Policy and CPI Targetting');
% % grid on;
% 
% figure;
% barh(bar_dataset);
% set(gca, 'YTickLabel', secname, 'YTick', 1:numel(secname));
% legend('DC index in Distorted Economy', 'DC Index in Efficient Economy', 'Output-Gap Targeting', 'Optimal Policy', 'Location', 'northeast');
% title('Inflation Index');
% 
% 
% % Rotate x-axis labels for clarity
% xtickangle(30);


%% Drawing log-productivity shocks with covariance matrix Sigmaz

% number of iterations
iteration = 100;

% initialize the random number generator to make the results reproducable
rng('default');
rng(1);

% drawing productivity shocks with variance-covariance matrix Sigmaz
dlogz = mvnrnd(zeros(N,1),Sigmaz,iteration);
dlogz = dlogz';


%% Monte Carlo simulation

[phi_star,pi_star,~,~,phi_og,~,~,~,phi_dc, phi_dc_eff]  = optimalpolicy(mu, theta, sigma_c, sigma, varepsilon, omega_tilde, alpha_tilde, beta_tilde, gamma, varphi);

phi_cpi      = beta_tilde';

iteration = size(dlogz,2);


w_cpi        = zeros(iteration,7);
w_og         = zeros(iteration,7);
w_star       = zeros(iteration,7);
w_dc         = zeros(iteration,7);

for i=1:iteration
    y_tilde_cpi          = solveoutputgap(phi_cpi, 0, dlogz(:,i));
    y_tilde_og           = solveoutputgap(phi_og, 0, dlogz(:,i));
    y_tilde_star         = solveoutputgap(phi_star, pi_star, dlogz(:,i));


    [w_cpi(i,1),w_cpi(i,2),w_cpi(i,3),w_cpi(i,4),w_cpi(i,5),w_cpi(i,6),w_cpi(i,7)]...
        = welfarefunction (y_tilde_cpi, dlogz(:,i));

    [w_og(i,1),w_og(i,2),w_og(i,3),w_og(i,4),w_og(i,5),w_og(i,6),w_og(i,7)]...
        = welfarefunction (y_tilde_og, dlogz(:,i));

    [w_star(i,1),w_star(i,2),w_star(i,3),w_star(i,4),w_star(i,5),w_star(i,6),w_star(i,7)]...
        = welfarefunction (y_tilde_star, dlogz(:,i));

end

%% cosine_similarity to optimal policy

CS_og             = dot(phi_star, phi_og) / (norm(phi_star) * norm(phi_og));
CS_cpi            = dot(phi_star, phi_cpi) / (norm(phi_star) * norm(phi_cpi));


%% Calculating the welfare gain for each policy

EWG_cpi  = mean(w_cpi);
EWG_og   = mean(w_og);
EWG_star = mean(w_star);

%% Expected welfare gain measured as a fraction of steady-state consumption

EWG_GDP_cpi   = expectedwelfaregain(EWG_cpi,gamma);
EWG_GDP_og    = expectedwelfaregain(EWG_og,gamma);
EWG_GDP_star  = expectedwelfaregain(EWG_star,gamma);

input.data    = [EWG_GDP_star; EWG_GDP_og; EWG_GDP_cpi];
input.data    = input.data';
input.data    = [input.data; 1, CS_og, CS_cpi];



%% Generates a LaTeX table

% Optional fields (if not set default values will be used):

% Set the position of the table in the LaTex document using h, t, p, b, H or !
input.tablePositioning = 'h';

% Set column labels (use empty string for no label):
input.tableColLabels = {'Optimal Policy', 'Output-Gap Targeting', 'CPI Targeting'};
% Set row labels (use empty string for no label):
input.tableRowLabels = {'Welfare Gain (Percent GDP)','\hline \quad First order welfare gain\\ \quad \quad Aggregate Wedge',...
    '\quad \quad Supply-side effect','\hline \quad Second order welfare gain\\ \quad \quad Output gap volatility',...
    '\quad \quad Within-sector price dispersion','\quad \quad Cross-sector price dispersion',...
    '\quad \quad Variation in allocative efficiency','\hline Cosine similarity to optimal policy'};

% Switch transposing/pivoting your table:
input.transposeTable = 0;
input.dataFormat = {'%.3f',3};
input.dataNanString = '-';
input.tableColumnAlignment = 'c';
input.tableBorders = 1;
input.booktabs = 1;

% LaTex table:
input.tableCaption = 'Welfare Gain under Various Policies';
input.tableLabel = 'welfare-gain-decomp';
input.makeCompleteLatexDocument = 1;
latex = latexTable(input);





















