% APPENDIX.M
%
% Produces the appendix 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_appendix/, which is created if it does not exist.
% Main text figures are produced by main.m and written to figures/.
%
%   Section  Output file                                Paper
%   -------  ----------------------------------------   ----------------
%      1     IndustryContribution.png                   Figure A.1
%      2     scatterplot.png                            Figure A.2
%      3     outputresponsedecomp_ucmarkup.eps/.png     Figure A.3
%      4     RBT_supplysideeffectdecomp.png/.eps        Figure A.4
%      5     supplysideeffectdecomp_apmarkup.eps        robustness, not
%            supplysideeffectdecomp_pfmarkup.eps        reported in the
%            inflationbiasdecomp_apmarkup.eps           current draft
%            inflationbiasdecomp_pfmarkup.eps
%
% Data inputs, all under data_clean/:
%   FPA/FPA_sum.mat          FPA                       66 x 1
%   IO_table/IO_sum.mat      IO_matrix_sum, indname    79 x 94 x 25, 79 x 1
%   markup/markup_sum.mat    mu_sum                    66 x 19 x 3
%
% Functions called, all under functions/:
%   supplyside.m, optimalpolicy.m, outputresponse.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'];

%% Load Data

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

%% Add Functions Path

addpath('functions');

figdir   = [cd, '/figures_appendix'];
mkp_name = {'ucmarkup', 'apmarkup', 'pfmarkup'};

if ~exist(figdir, 'dir')
    mkdir(figdir);
end

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



%% Section 1: Figure A.1 - Industry Contribution to the Supply-Side Effect
% Decomposes the supply-side effect of Proposition 1 in the calibration year
% into the contribution of each industry, and splits every contribution into
% the two channels of the proposition: substitution between sectoral
% products (SSP_sec) and substitution between labor and intermediate inputs
% (SLI_sec). Households appear as the first bar: labor does not enter
% consumption (b_L = 0), so their entire contribution SC is substitution
% between sectoral products and the second channel is zero for them.

year           = 2006;
temp           = IO(:,:,year-1996);
temp           = temp';
Omega          = bsxfun(@rdivide, temp, grosscost(:,year-1996));
alpha          = labor(:,year-1996);
beta           = Final(1:66,year-1996);
beta(beta<0)   = 0;
beta           = beta/sum(beta);

mu             = mu_sum(:, year - 1996, 1);

[T_ind, ~, ~, SC_ind, ~, ~, SSP_sec, SLI_sec] = supplyside(mu, theta, sigma_c, sigma, Omega, alpha, beta);

% Column 1: substitution between sectoral products
% Column 2: substitution between labor and intermediate inputs
bar_dataset        = [SC_ind, 0; SSP_sec, SLI_sec];
percentages        = (sum(bar_dataset, 2) / T_ind) * 100;
secname            = indname(1:66,:);
secname            = [{'Households'}; secname];
percentages        = flipud(percentages);
bar_dataset        = flipud(bar_dataset);
secname            = flipud(secname);

colors  = [102, 204, 0; 128, 128, 255; 255, 128, 128] / 255;

figure('Position', [100, 100, 700, 1500]);
set(gcf, 'DefaultAxesFontName', 'Times New Roman');
b = barh(bar_dataset, 'stacked', 'EdgeColor', 'k', 'LineWidth', 1.5);
b(1).FaceColor = colors(2,:);
b(2).FaceColor = colors(3,:);
set(gca, 'YTickLabel', secname, 'YTick', 1:numel(secname), 'FontSize', 8);
l = legend('Substitution b/w Sectoral Products', ...
    'Substitution b/w Labor & Intermediates', ...
    'Location', 'southoutside', 'Orientation', 'horizontal');
l.FontSize = 8;

% Annotate bars with each industry's share of the total supply-side effect,
% offset to the right of the bar tip so the label clears the bar edge. A
% negative segment is drawn to the left of zero, so the right-hand tip of a
% stacked bar is the sum of its positive segments, not the row total.
bar_right = sum(max(bar_dataset, 0), 2);
offset    = 0.015 * max(bar_right);

hold on;
for i = 1:length(bar_right)
    text(bar_right(i) + offset, i, sprintf('%.1f%%', percentages(i)), ...
        'HorizontalAlignment', 'left', 'VerticalAlignment', 'middle', ...
        'FontSize', 7)
end
hold off;

% Leave room on the right so the labels are not clipped
xl = xlim;
xlim([xl(1), max(bar_right) + 10*offset]);

exportgraphics(gcf, fullfile(figdir, 'IndustryContribution.png'),'Resolution',600);



%% Section 2: Figure A.2 - Upstream Markup vs. Wage Pass-Through
% Scatter of the wage pass-through to sectoral prices (Psi_f_hat) against
% the inverse upstream markup (Psi_f) in the calibration year, one panel per
% markup measure. The weighted covariance between these two objects is the
% sufficient statistic for the supply-side effect in Proposition 1.

figure('Position', [100, 100, 900, 350]);
set(gcf, 'DefaultAxesFontName', 'Times New Roman');

for k = 1:3

    year         = 2006;

    temp         = IO(:,:,year-1996);
    temp         = temp';

    Omega        = bsxfun(@rdivide, temp, grosscost(:,year-1996));
    alpha        = labor(:,year-1996);
    beta         = Final(1:66,year-1996);
    beta(beta<0) = 0;
    beta         = beta/sum(beta);

    mu           = mu_sum(:, year-1996, k);

    % Revenue-based IO variables: Psi_f is the inverse upstream markup
    Omega_p        = diag(mu)^(-1) * Omega;
    Omega_f        = diag(mu)^(-1) * alpha;
    Psi_p          = (eye(N)-Omega_p)^(-1);
    Psi_f          = Psi_p * Omega_f;

    % Sticky-price IO variables: Psi_f_hat is the wage pass-through
    Omega_p_hat    = diag(theta) * Omega;
    Omega_f_hat    = diag(theta) * alpha;
    Psi_p_hat      = (eye(N)-Omega_p_hat)^(-1);
    Psi_f_hat      = Psi_p_hat * Omega_f_hat;

    x = Psi_f;
    y = Psi_f_hat;

    % Broad industry groups used only for the marker colors
    cate        = zeros(66,1);
    cate(1:2)   = 1;  % Agriculture
    cate(3:26)  = 2;  % Manufacturing
    cate(27:66) = 3;  % Services

    colors  = [102, 204, 0; 128, 128, 255; 255, 128, 128] / 255;
    markers = {'^', 's', 'o'};

    subplot(1,3,k);
    hold on;

    for i = 1:3
        scatter(x(cate == i), y(cate == i), 20, colors(i,:), 'filled', 'Marker', markers{i});
    end

    % Fitted line
    p     = polyfit(x, y, 1);
    x_fit = linspace(min(x), max(x), 100);
    y_fit = polyval(p, x_fit);
    plot(x_fit, y_fit, 'k-', 'LineWidth', 1);

    xlabel('Upstream markup $\Psi_{iL}$', 'Interpreter', 'latex');
    ylabel('Wage pass-through $\hat{\Psi}_{iL}$', 'Interpreter', 'latex');
    if k == 1
        title('UC Markups', 'Interpreter', 'latex');
    elseif k == 2
        title('AP Markups', 'Interpreter', 'latex');
    else
        title('PF Markups', 'Interpreter', 'latex');
    end

    box on;
    grid on;
    axis tight;

    hold off;

end

legend('Agriculture', 'Manufacturing', 'Services', 'Location', 'best');

exportgraphics(gcf, fullfile(figdir, 'scatterplot.png'));



%% Section 3: Figure A.3 - Output Response Decomposition
% Splits the elasticity of real output with respect to the money supply into
% the supply-side effect (the TFP response) and the demand-side effect (the
% labor response), following Proposition 2. Panel (a) uses the observed IO
% table; panel (b) shuts down input-output linkages.

T_or          = zeros(length(yr),2);
supplyside_or = zeros(length(yr),2);
demandside_or = zeros(length(yr),2);

for k = 1   % set to 1:3 to also produce the AP and PF markup versions

    for t = 1:length(yr)

        mu           = mu_sum(:, t, k);

        % Benchmark, time-varying 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);
        [supplyside_or(t,1), demandside_or(t,1)] = outputresponse(mu, theta, sigma_c, sigma, Omega, alpha, beta, gamma, varphi);

        % No IO
        Omega        = zeros(N,N);
        alpha        = ones(N,1);
        [supplyside_or(t,2), demandside_or(t,2)] = outputresponse(mu, theta, sigma_c, sigma, Omega, alpha, beta, gamma, varphi);

    end

    T_or = supplyside_or + demandside_or;


    figure('Position', [100, 100, 560, 180])
    set(gcf, 'DefaultAxesFontName', 'Times New Roman');

    purple = [0.5, 0, 0.5];

    subplot(1,2,1);
    hold on;
    plot(yr, supplyside_or(:,1), '-b', 'LineWidth', 1);
    plot(yr, demandside_or(:,1), '--', 'Color', purple, 'LineWidth', 1);
    plot(yr, T_or(:,1), ':o', 'Color', 'r', 'LineWidth', 1, 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r', 'MarkerSize', 3);
    box on;
    xlim([1997, 2015]);
    title("(a) Benchmark", 'FontSize', 10, 'FontName', 'Times New Roman','FontWeight', 'normal');
    xlabel('Year', 'FontSize', 10);
    ylabel('Output Response', 'FontSize', 10, 'Interpreter','latex');

    subplot(1,2,2);
    hold on;
    plot(yr, supplyside_or(:,2), '-b', 'LineWidth', 1);
    plot(yr, demandside_or(:,2), '--', 'Color', purple, 'LineWidth', 1);
    plot(yr, T_or(:,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('Output Response', 'FontSize', 10, 'Interpreter','latex');
    l = legend('Supply-Side Effect', 'Demand-Side Effect', 'Total', 'FontSize', 8);
    set(l, 'Location', 'best');
    legend('boxon');

    % Put both panels on a common vertical scale
    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, ['outputresponsedecomp_' mkp_name{k}]), '-depsc');
    saveas(gcf, fullfile(figdir, ['outputresponsedecomp_' mkp_name{k} '.png']));

end



%% Section 4: Figure A.4 - Varying Homogeneous Markups and Constant IO
% Two robustness variants of the supply-side effect decomposition. Panel (a)
% replaces the observed markup vector with the Domar-weighted harmonic mean
% markup, so markups are common across industries but still vary over time.
% Panel (b) holds the IO table, labor shares and final expenditure shares
% fixed at their 1997 values.

T_rbt          = zeros(length(yr),2);
SC_rbt         = zeros(length(yr),2); % substitution between sectoral products
SP_rbt         = zeros(length(yr),2); % substitution between labor & intermediates

% Constant (1997) IO calibration
year           = 1997;
temp           = IO(:,:,year-1996);
temp           = temp';
Omega_cons     = bsxfun(@rdivide, temp, grosscost(:,year-1996));
beta_cons      = Final(1:66,year-1996);
beta_cons(beta_cons<0) = 0;
beta_cons      = beta_cons/sum(beta_cons);
alpha_cons     = labor(:,year-1996);

for t = 1:length(yr)

    mu           = mu_sum(:, t, 1);

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

    % Domar-weighted harmonic mean markup, imposed on every industry
    Domar  = (beta' * (eye(N)-diag(mu)^(-1) * Omega)^(-1))';
    Domar  = Domar./(sum(Domar));
    mu_hm  = (Domar' * mu.^(-1))^(-1);
    mu_hm  = mu_hm * ones(N,1);

    [T_rbt(t,1), SC_rbt(t,1), SP_rbt(t,1)] = supplyside(mu_hm, theta, sigma_c, sigma, Omega, alpha, beta);

    % Constant IO
    [T_rbt(t,2), SC_rbt(t,2), SP_rbt(t,2)] = supplyside(mu, theta, sigma_c, sigma, Omega_cons, alpha_cons, beta_cons);

end


figure('Position', [100, 100, 560, 180])
set(gcf, 'DefaultAxesFontName', 'Times New Roman');

purple = [0.5, 0, 0.5];

subplot(1,2,1);
hold on;
plot(yr, SC_rbt(:,1), '-b', 'LineWidth', 1);
plot(yr, SP_rbt(:,1), '--', 'Color', purple, 'LineWidth', 1);
plot(yr, T_rbt(:,1), ':o', 'Color', 'r', 'LineWidth', 1, 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r', 'MarkerSize', 3);
box on;
xlim([1997, 2015]);
title("(a) Varing Homog. Markups ", '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', 'northeast');
legend('boxon');

subplot(1,2,2);
hold on;
plot(yr, SC_rbt(:,2), '-b', 'LineWidth', 1);
plot(yr, SP_rbt(:,2), '--', 'Color', purple, 'LineWidth', 1);
plot(yr, T_rbt(:,2), ':o', 'Color', 'r', 'LineWidth', 1, 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r', 'MarkerSize', 3);
box on;
xlim([1997, 2015]);
title("(b) Constant IO", 'FontSize', 10, 'FontName', 'Times New Roman','FontWeight', 'normal');
xlabel('Year', 'FontSize', 10);
ylabel('$d \log TFP/d \log w$', 'FontSize', 10, 'Interpreter','latex');

% Put both panels on a common vertical scale
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

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



%% Section 5: Alternative Markup Measures (AP and PF)
% Repeats the supply-side effect decomposition and the inflation bias
% decomposition of the main text under the accounting profit and production
% function markups. Not reported in the current draft.

for k = 2:3

    T_k          = zeros(length(yr),2); % Total supply-side effect
    SSP_k        = zeros(length(yr),2); % Substitution between sectoral products
    SLI_k        = zeros(length(yr),2); % Substitution between labor & intermediates
    pi_aw_k      = zeros(length(yr),2); % Inflation bias due to aggregate wedge
    pi_ss_k      = zeros(length(yr),2); % Inflation bias due to supply-side effect

    for t = 1:length(yr)

        mu           = mu_sum(:, t, k);

        % Benchmark, time-varying 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_k(t,1), SSP_k(t,1), SLI_k(t,1)] = 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_aw_k(t,1) = pi_star1 / sum(phi_star);
        pi_ss_k(t,1) = pi_star2 / sum(phi_star);

        % No IO
        Omega        = zeros(N,N);
        alpha        = ones(N,1);

        [T_k(t,2), SSP_k(t,2), SLI_k(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_aw_k(t,2) = pi_star1 / sum(phi_star);
        pi_ss_k(t,2) = pi_star2 / sum(phi_star);

    end

    pi_aw_k = pi_aw_k * 100;
    pi_ss_k = pi_ss_k * 100;


    % Supply-side effect decomposition

    figure('Position', [100, 100, 560, 180])
    set(gcf, 'DefaultAxesFontName', 'Times New Roman');

    purple = [0.5, 0, 0.5];

    subplot(1,2,1);
    hold on;
    plot(yr, SSP_k(:,1), '-b', 'LineWidth', 1);
    plot(yr, SLI_k(:,1), '--', 'Color', purple, 'LineWidth', 1);
    plot(yr, T_k(:,1), ':o', 'Color', 'r', 'LineWidth', 1, 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r', 'MarkerSize', 3);
    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_k(:,2), '-b', 'LineWidth', 1);
    plot(yr, SLI_k(:,2), '--', 'Color', purple, 'LineWidth', 1);
    plot(yr, T_k(:,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');
    legend('boxon');

    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, ['supplysideeffectdecomp_' mkp_name{k}]), '-depsc');


    % Inflation bias decomposition

    figure('Position', [100, 100, 560, 180])
    set(gcf, 'DefaultAxesFontName', 'Times New Roman');

    subplot(1,2,1)
    hold on;
    plot(yr, pi_aw_k(:,1), '--o', 'LineWidth', 1, 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r', 'MarkerSize', 3, 'Color', 'r');
    plot(yr, pi_ss_k(:,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_aw_k(:,2), '--o', 'LineWidth', 1, 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r', 'MarkerSize', 3, 'Color', 'r');
    plot(yr, pi_ss_k(:,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;

    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_' mkp_name{k}]), '-depsc');

end
