%% This code generates Figure 2, Figure A.6 and Figure A.7

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

dlogz_quarterly = log_MFP_linear_detrended(38:113,:);
dlogz_quarterly_reshaped = reshape(dlogz_quarterly, 4, []);
meanValues = mean(dlogz_quarterly_reshaped, 1);
dlogz = reshape(meanValues, 19, 66);
dlogz = dlogz';

% transforming the data from monthly to quarterly and calculating the
% degree of price flexibility of each industry

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;

% 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,:) = [];
grosscost = grosscost(:,1:19);


% Initialize a year for IO Matrix


yr = 1997:2015;
pi_star_aw = zeros(19,2);
pi_star_ss = zeros(19,2);

for k = 1

    for t = 1: length(yr)



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


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



        % Benchmark
        mu           = mu_sum(:, t, k);

        [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

        Omega        = zeros(N,N);
        alpha        = ones(N,1);

        [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


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


    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


    figdir = [cd, '/figures'];

    if k == 1
        print(fullfile(figdir, 'inflationbiasdecomp_ucmarkup'), '-depsc');
    elseif k == 2
        print(fullfile(figdir, 'inflationbiasdecomp_apmarkup'), '-depsc');
    else
        print(fullfile(figdir, 'inflationbiasdecomp_pfmarkup'), '-depsc');
    end


end