
% MAIN.M
%
% Produces the main text figures of "Production Networks and Monetary Policy
% in Distorted Economies". Run from the root of the replication package
% (the directory that contains data_clean/ and functions/). All output is
% written to figures/. The appendix figures are produced by appendix.m and
% written to figures_appendix/; the welfare gain table is produced by
% tabularwelfaregiandecomp.m.
%
%   Output file                                    Paper
%   --------------------------------------------   ----------------
%   supplysideeffectdecomp_ucmarkup.eps            Figure 1
%   supplysideeffectdecomp.png                     Figure 1
%   inflationbiasdecomp_ucmarkup.eps               Figure 2
%   RBT_inflationbiasdecomp_preferences.eps        Figure 3
%   RBT_inflationbiasdecomp_preferences.png        Figure 3
%
% Data inputs, all under data_clean/:
%   FPA/FPA_sum.mat          FPA                       66 x 1
%   IO_table/IO_sum.mat      IO_matrix_sum             79 x 94 x 25
%   markup/markup_sum.mat    mu_sum                    66 x 19 x 3
%   productivity/MFP.mat     Sigmaz                    66 x 66
%                            loaded here for consistency across the driver
%                            scripts; used by tabularwelfaregiandecomp.m
%
% Functions called, all under functions/:
%   supplyside.m, optimalpolicy.m, weightedcov.m
%
% Markup measures are stored in the third dimension of mu_sum:
%   (:,:,1) user cost (UC), (:,:,2) accounting profit (AP),
%   (:,:,3) production function (PF). The paper reports UC markups.

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 each year the supply-side effect of monetary policy is evaluated from
% the sufficient statistic of Proposition 1, split into its two channels:
% substitution between sectoral products (SSP) and substitution between
% labor and intermediate inputs (SLI), with T the total. The inflation bias
% is then read off the right-hand side of the optimal policy condition in
% Proposition 5, whose two terms are the aggregate wedge (pi_star_aw) and
% the supply-side effect (pi_star_ss); both are normalized by the sum of the
% sectoral inflation weights phi_star so that the weights add to one.
% Column 1 of each matrix uses the observed input-output table, column 2
% shuts down input-output linkages (Omega = 0, alpha = 1).

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
% The two channels of Proposition 1 and their total, 1997-2015. Panel (a)
% uses the observed input-output table, panel (b) shuts down input-output
% linkages. Both panels are drawn on a common vertical scale.

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
% The two sources of inflation bias in Proposition 5, in percent, 1997-2015.
% Panel (a) uses the observed input-output table, panel (b) shuts down
% input-output linkages. Both panels are drawn on a common vertical scale.

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
% Recomputes the inflation bias decomposition of Proposition 5 under four
% alternative preference calibrations, holding the observed input-output
% table fixed. The columns of pi_star_aw_rbt and pi_star_ss_rbt are:
%   1  high wealth effect        gamma  = 2
%   2  low wealth effect         gamma  = 0.5
%   3  high Frisch elasticity    varphi = 0.5
%   4  low Frisch elasticity     varphi = 5
% Columns 5 and 6 are preallocated but unused.

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: Inflation Bias under Different Preference Parameters
% One panel per alternative calibration from the robustness block above, in
% percent, 1997-2015. All four panels are drawn on a common vertical scale.

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');
saveas(gcf, fullfile(figdir, 'RBT_inflationbiasdecomp_preferences.png'));




