
clc;
clear;

%% Set Paths

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 Data

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

%% Add Functions Path

addpath('functions');

%% Transform Monthly Data to Quarterly

theta = 1 - exp(-3*FPA);

%% Baseline calibration

% Number of industries
N              = 66;

% Preference parameters
gamma          = 1;
varphi         = 2;

% Within-sector elasticities
varepsilon     = 6 * ones(N,1);

% Cross-sector elastiticities
sigma        = 0.5 * ones(N,1);
sigma_c      = 1;

% Process 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,:) = [];

% Set years range
yr             = 1997:2015;
start_date     = yr(1);


% Preallocate Matrices for Efficiency
T           = zeros(length(yr),2); % Total supply-side effect
SSP         = zeros(length(yr),2); % Substitution between sectoral products
SLI         = zeros(length(yr),2); % Substitution between labor & intermediates
pi_star_aw  = zeros(length(yr),2); % Inflation bias due to aggregate wedge
pi_star_ss  = zeros(length(yr),2); % Inflation bias due to supply-side effect


%% Compute Supply-Side Effects and Inflation Bias

for t = 1:length(yr)

    % Time-varing markups
    mu           = mu_sum(:, t, 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)

    % Baseline IO calibration (Time-varing IO)
    temp         = IO(:,:,t);
    temp         = temp';
    Omega        = bsxfun(@rdivide, temp, grosscost(:,t));
    alpha        = labor(:,t);
    beta         = Final(1:66,t);
    beta(beta<0) = 0;
    beta         = beta/sum(beta);
    [T(t,1), SSP(t,1), SLI(t,1)] = supplyside(mu, theta, sigma_c, sigma, Omega, alpha, beta);

    % Compute Supply-Side Effect
    [T(t,1), SSP(t,1), SLI(t,1)] = supplyside(mu, theta, sigma_c, sigma, Omega, alpha, beta);

    % Compute Inflation Bias
    [phi_star,~,pi_star1,pi_star2] = optimalpolicy(mu, theta, sigma_c, sigma, varepsilon, Omega, alpha, beta, gamma, varphi);
    pi_star_aw(t,1) = pi_star1 / sum(phi_star);
    pi_star_ss(t,1) = pi_star2 / sum(phi_star);


    % No IO case
    Omega        = zeros(N,N);
    alpha        = ones(N,1);
    [T(t,2), SSP(t,2), SLI(t,2)] = supplyside(mu, theta, sigma_c, sigma, Omega, alpha, beta);
    [phi_star,~,pi_star1,pi_star2] = optimalpolicy(mu, theta, sigma_c, sigma, varepsilon, Omega, alpha, beta, gamma, varphi);
    pi_star_aw(t,2) = pi_star1 / sum(phi_star);
    pi_star_ss(t,2) = pi_star2 / sum(phi_star);

end

%% Convert Inflation Rates to Percentage

pi_star_aw = pi_star_aw * 100;
pi_star_ss = pi_star_ss * 100;

%% Figure 1: Plot Supply-Side Effect Decomposition

figure('Position', [100, 100, 560, 180]) % [left, bottom, width, height] in pixels

% Set the default font for the figure to Arial
set(gcf, 'DefaultAxesFontName', 'Times New Roman');

% Define a purple color using RGB values
purple = [0.5, 0, 0.5];  % This is purple

subplot(1,2,1);
hold on;
plot(yr, SSP(:,1), '-b', 'LineWidth', 1); % Blue solid line
plot(yr, SLI(:,1), '--', 'Color', purple, 'LineWidth', 1); % Purple dash line
plot(yr, T(:,1), ':o', 'Color', 'r', 'LineWidth', 1, 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r', 'MarkerSize', 3); % Red dot line
box on;
xlim([1997, 2015]);
title("(a) Benchmark", 'FontSize', 10, 'FontName', 'Times New Roman','FontWeight', 'normal');
xlabel('Year', 'FontSize', 10);
ylabel('$d \log TFP/d \log w$', 'FontSize', 10, 'Interpreter','latex');


subplot(1,2,2);
hold on;
plot(yr, SSP(:,2), '-b', 'LineWidth', 1); % Blue solid line
plot(yr, SLI(:,2), '--', 'Color', purple, 'LineWidth', 1); % Purple dash line
plot(yr, T(:,2), ':o', 'Color', 'r', 'LineWidth', 1, 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r', 'MarkerSize', 3);
box on;
xlim([1997, 2015]);
title("(b) No IO", 'FontSize', 10, 'FontName', 'Times New Roman','FontWeight', 'normal');
xlabel('Year', 'FontSize', 10);
ylabel('$d \log TFP/d \log w$', 'FontSize', 10, 'Interpreter','latex');
l = legend('Substitution b/w Sectoral Products', 'Substitution b/w Labor & Intermediates', 'Total', 'FontSize', 8);
set(l, 'Location', 'best'); % Place the legend in the best position
legend('boxon');

% Initialize min and max y-limits
minY = Inf;
maxY = -Inf;

for t = 1:2
    subplot(1, 2, t);
    currentYLim = ylim;
    minY = min(minY, currentYLim(1));
    maxY = max(maxY, currentYLim(2));
end

for t = 1:2
    subplot(1, 2, t);
    ylim([minY, maxY]);
end

figdir = [cd, '/figures'];


print(fullfile(figdir, 'supplysideeffectdecomp_ucmarkup'), '-depsc');
saveas(gcf, fullfile(figdir, 'supplysideeffectdecomp.png'));



%% Figure 2: Plot Inflation Bias Decomposition

figure('Position', [100, 100, 560, 180]) % [left, bottom, width, height] in pixels

% Set the default font for the figure to Arial
set(gcf, 'DefaultAxesFontName', 'Times New Roman');


subplot(1,2,1)
hold on;
plot(yr, pi_star_aw(:,1), '--o', 'LineWidth', 1, 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r', 'MarkerSize', 3, 'Color', 'r');
plot(yr,pi_star_ss(:,1),'-b', 'LineWidth', 1);
box on;
xlim([1997, 2015]);
title("(a) Benchmark", 'FontSize', 10, 'FontName', 'Times New Roman','FontWeight', 'normal');
xlabel('Year', 'FontSize', 10);
ylabel('Inflation Rate (%)', 'FontSize', 10);

subplot(1,2,2)
hold on;
plot(yr, pi_star_aw(:,2), '--o', 'LineWidth', 1, 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r', 'MarkerSize', 3, 'Color', 'r');
plot(yr,pi_star_ss(:,2),'-b', 'LineWidth', 1);
box on;
xlim([1997, 2015]);
title("(b) No IO ", 'FontSize', 10, 'FontName', 'Times New Roman','FontWeight', 'normal');
xlabel('Year', 'FontSize', 10);
ylabel('Inflation Rate (%)', 'FontSize', 10);

l = legend('Aggregate Wedge','Supply-Side Effect','Orientation','vertical','Location','best');
l.FontSize = 8;

% Initialize min and max y-limits
minY = Inf;
maxY = -Inf;

for i = 1:2
    subplot(1, 2, i);
    currentYLim = ylim;
    minY = min(minY, currentYLim(1));
    maxY = max(maxY, currentYLim(2));
end

for i = 1:2
    subplot(1, 2, i);
    ylim([minY, maxY]);
end

print(fullfile(figdir, 'inflationbiasdecomp_ucmarkup'), '-depsc');



%% Robustness Analysis

pi_star_aw_rbt     = zeros(19,6);
pi_star_ss_rbt     = zeros(19,6);

for t = 1:length(yr)

    % Time-varing markups
    mu           = mu_sum(:, t, 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)

    % Baseline IO calibration (Time-varing IO)
    temp         = IO(:,:,t);
    temp         = temp';
    Omega        = bsxfun(@rdivide, temp, grosscost(:,t));
    alpha        = labor(:,t);
    beta         = Final(1:66,t);
    beta(beta<0) = 0;
    beta         = beta/sum(beta);

    % Preference parameters
    gamma  = 1;
    varphi = 2;

    % High Wealth Effect
    gamma  = 2;

    [phi_star,~,pi_star1,pi_star2]    = optimalpolicy(mu, theta, sigma_c, sigma, varepsilon, Omega, alpha, beta, gamma, varphi);

    pi_star_aw_rbt(t,1)    = pi_star1/sum(phi_star);
    pi_star_ss_rbt(t,1)    = pi_star2/sum(phi_star);

    % Low Wealth Effect
    gamma = 0.5;
    [phi_star,~,pi_star1,pi_star2]    = optimalpolicy(mu, theta, sigma_c, sigma, varepsilon, Omega, alpha, beta, gamma, varphi);
    pi_star_aw_rbt(t,2)    = pi_star1/sum(phi_star);
    pi_star_ss_rbt(t,2)    = pi_star2/sum(phi_star);

    gamma  = 1;
    % High Frisch elasticity
    varphi = 0.5;
    [phi_star,~,pi_star1,pi_star2]    = optimalpolicy(mu, theta, sigma_c, sigma, varepsilon, Omega, alpha, beta, gamma, varphi);
    pi_star_aw_rbt(t,3)    = pi_star1/sum(phi_star);
    pi_star_ss_rbt(t,3)    = pi_star2/sum(phi_star);


    % Low Frisch elasticity
    varphi = 5;
    [phi_star,~,pi_star1,pi_star2]    = optimalpolicy(mu, theta, sigma_c, sigma, varepsilon, Omega, alpha, beta, gamma, varphi);
    pi_star_aw_rbt(t,4)    = pi_star1/sum(phi_star);
    pi_star_ss_rbt(t,4)    = pi_star2/sum(phi_star);
end

pi_star_aw_rbt = pi_star_aw_rbt * 100;
pi_star_ss_rbt = pi_star_ss_rbt * 100;

%% Figure 3

figure;

% Set the default font for the figure to Arial
set(gcf, 'DefaultAxesFontName', 'Times New Roman');

yr = 1997:2015;

% Define a deep red color using RGB values
deepRed = [139, 0, 0] / 255;

subplot(2,2,1)
hold on;
plot(yr,pi_star_aw_rbt(:,1),'--r', 'LineWidth', 1);
plot(yr,pi_star_ss_rbt(:,1),'-b', 'LineWidth', 1);
box on;
xlim([1997, 2015]);
title("(a) High Wealth Effect ($\gamma = 2$)", 'FontSize', 10, 'FontName', 'Times New Roman','FontWeight', 'normal', 'interpreter', 'latex');
xlabel('Year', 'FontSize', 10);
ylabel('Inflation Rate (%)', 'FontSize', 10);

subplot(2,2,2)
hold on;
plot(yr,pi_star_aw_rbt(:,2),'--r', 'LineWidth', 1);
plot(yr,pi_star_ss_rbt(:,2),'-b', 'LineWidth', 1);
box on;
xlim([1997, 2015]);
ylim([-0.1, 0.8]);
title("(b) Low Wealth Effect ($\gamma = 0.5$)", 'FontSize', 10, 'FontName', 'Times New Roman','FontWeight', 'normal', 'interpreter', 'latex');
xlabel('Year', 'FontSize', 10);
ylabel('Inflation Rate (%)', 'FontSize', 10);

subplot(2,2,3)
hold on;
plot(yr,pi_star_aw_rbt(:,3),'--r', 'LineWidth', 1);
plot(yr,pi_star_ss_rbt(:,3),'-b', 'LineWidth', 1);
box on;
xlim([1997, 2015]);
title("(c) High Frisch elasticity ($1/\varphi = 2$)", 'FontSize', 10, 'FontName', 'Times New Roman','FontWeight', 'normal', 'interpreter', 'latex');
xlabel('Year', 'FontSize', 10);
ylabel('Inflation Rate (%)', 'FontSize', 10);

subplot(2,2,4)
hold on;
plot(yr,pi_star_aw_rbt(:,4),'--r', 'LineWidth', 1);
plot(yr,pi_star_ss_rbt(:,4),'-b', 'LineWidth', 1);
box on;
xlim([1997, 2015]);
title("(d) Low Frisch elasticity ($1/\varphi = 0.2$)", 'FontSize', 10, 'FontName', 'Times New Roman','FontWeight', 'normal', 'interpreter', 'latex');
xlabel('Year', 'FontSize', 10);
ylabel('Inflation Rate (%)', 'FontSize', 10);


% Initialize min and max y-limits
minY = Inf;
maxY = -Inf;

for i = 1:4
    subplot(2, 2, i);
    currentYLim = ylim;
    minY = min(minY, currentYLim(1));
    maxY = max(maxY, currentYLim(2));
end

for i = 1:4
    subplot(2, 2, i);
    ylim([minY, maxY]);
end

l = legend('Aggregate Wedge','Supply-Side Effect','Orientation','vertical','Location','best');
legend('boxon');
l.FontSize = 8;


print(fullfile(figdir, 'RBT_inflationbiasdecomp_preferences'), '-depsc');


%% Table 1: Welfare Gains under Various Policies

% Preference parameters
gamma          = 1;
varphi         = 2;

% Computes sample averages of input-output (IO) coefficients, final expenditure shares, and markups over the period 1997–2015. 

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

% Drawing log-productivity shocks with covariance matrix Sigmaz

% number of iterations
iteration = 10000;

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


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

% Generates a LaTeX table

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


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





