%% Load markups
clear;
clc;

% Need to run IO_import.m first to import IO data
[ParentPath, ~, ~] = fileparts(pwd);
Path = fullfile(ParentPath, 'IO_table', 'IO_sum.mat');

load(Path);

mu_types = {'mu_cfg1', 'mu_li', 'markup_cobb_douglas'};

% Define the path to the raw data directory (without changing any logic)
rawDataDir = fullfile(cd, '..', '..', 'data_raw', 'markup'); % Access myfile/data_raw/markup

% Scales compustat firms to take up the whole industry, or only a part of
% the industry with a blob of unobserved firms left over. For micro impact of macro shocks
% we want there to be a blob. For markup reduction, we don't want there to
% be a blob.

blob = 0;
start_date = 1997;

%% Calibrate the IO Table
IO = IO_matrix_sum;
IO(isnan(IO)) = 0;
IO(IO < 0) = 0;

grossout = squeeze(sum(IO(1:66,1:66,:))) + squeeze(IO(77,1:66,:)) + squeeze(IO(79,1:66,:));

for i = 1:3
    mu_type = mu_types{i};

    for year = 1997:2015
        error = 0;

        for T = 1:66
            try
                % Construct absolute file path for CSV
                csvFilePath = fullfile(rawDataDir, sprintf('Firm_level_mu_%d_%d_%s.csv', T, year, mu_type));

                % Use csvread with the absolute file path
                eval(strcat('data_matrix_', num2str(T), ' = csvread(''', csvFilePath, ''',1,0);'));               
                eval(strcat('temp = data_matrix_', num2str(T), ';'));
                
                if blob == 1
                    eval(strcat('temp = [temp temp(:,3)./grossout(',num2str(T),',year-start_date+1)];'));
                    temp = [temp(:,1) temp(:,3) temp(:,4) temp(:,5)];
                    if sum(temp(:,4)) >1
                        temp(:,4) = temp(:,4)/sum(temp(:,4));
                    end
                else
                    eval(strcat('temp = [temp temp(:,3)./sum(temp(:,3))];'));
                    temp = [temp(:,1) temp(:,3) temp(:,4) temp(:,5)];
                end

                average_mu(T) = 1./sum((1./temp(:,3)).*temp(:,2)/sum(temp(:,2)));
                eval(strcat('data_matrix_', num2str(T),'=temp'));
            catch
                error = error +1;
            end

        end
        %replace missings with the industry-weighted average markup
        average_mu(average_mu==0) = (average_mu(average_mu~=0)*grossout(average_mu~=0,year-start_date+1))/sum(grossout(average_mu~=0,year-start_date+1));
        
        mu_temp(:,year - 1996) = average_mu';

    end
    mu_sum(:,:,i)= mu_temp;
end

% Save processed markup data in the current directory (data_clean/markup)
save('markup_sum.mat', 'mu_sum');