%% This code generates Figure 3


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


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


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

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

% across-sector elasticities

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


start_date     = 1997;

pi_star_aw     = zeros(19,6);
pi_star_ss     = zeros(19,6);


for year = 1997:2015


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


    % baseline calibration

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


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

    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(year-1996,1)    = pi_star1/sum(phi_star);
    pi_star_ss(year-1996,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(year-1996,2)    = pi_star1/sum(phi_star);
    pi_star_ss(year-1996,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(year-1996,3)    = pi_star1/sum(phi_star);
    pi_star_ss(year-1996,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(year-1996,4)    = pi_star1/sum(phi_star);
    pi_star_ss(year-1996,4)    = pi_star2/sum(phi_star);
end

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

figure;

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

yr = 1997:2015;

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

subplot(2,2,1)
hold on;
plot(yr,pi_star_aw(:,1),'--r', 'LineWidth', 1);
plot(yr,pi_star_ss(:,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(:,2),'--r', 'LineWidth', 1);
plot(yr,pi_star_ss(:,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(:,3),'--r', 'LineWidth', 1);
plot(yr,pi_star_ss(:,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(:,4),'--r', 'LineWidth', 1);
plot(yr,pi_star_ss(:,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;


figdir = [cd, '/figures'];
print(fullfile(figdir, 'RBT_inflationbiasdecomp_preferences'), '-depsc');



