% This file plots cross-subsidization in interest rates and generates Figure 2

%% Import all debts and associated prices in benchmark ergodic distribution

jj = 0;  % initialize education group counter

% initialize arrays
qqdbyI_behavioral = cell(1,3);
qqdbyI_rational = qqdbyI_behavioral;
Price = struct([]);
AmountBorrowed = Price;

for group = {'hs','col'}
    jj = jj+1;

    % find DebtDump files that contain all credit taken in ergodic distribution
    for i = 1:length(experimentFolders)
        if experimentFolders(i).isdir && contains(experimentFolders(i).name,'dump') && contains(experimentFolders(i).name,'benchmark') && contains(experimentFolders(i).name,group)
            dumpDir = experimentFolders(i).name;
            break
        end
    end

    type = {'Goofy','Realist'};   % behavioral and rational types
    delimiter = ' ';
    startRow = 3;

    % Format for each line of text:
    %   column1: double (%f)
    %	column2: double (%f)
    %   column3: double (%f)
    %   column4: double (%f)
    %   column5: double (%f)
    formatSpec = '%3f%3f%8f%8f%f%[^\n\r]';


    % Load both DebtDumps
    for i = 1:2

        % Open the text file.
        fileID = fopen(fullfile(runDir,dumpDir,['/DebtDump',type{i},'.txt']),'r');

        % Read columns of data according to the format.
        dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'TextType', 'string', 'EmptyValue', NaN, 'HeaderLines' ,startRow-1, 'ReturnOnError', false, 'EndOfLine', '\r\n');

        % Close the text file.
        fclose(fileID);

        % Allocate imported array to column variable names for both types
        FairPrice.(type{i}) = dataArray{:, 5};
        AmountBorrowed(jj).(type{i}) = dataArray{:, 3};
        Price(jj).(type{i}) = dataArray{:, 4};

        %drop goofies that could not have borrowed without crosssubsidization
        dropindex = (FairPrice.(type{i})==0);

        FairPrice.(type{i})(dropindex) = [];
        AmountBorrowed(jj).(type{i})(dropindex) = [];
        Price(jj).(type{i})(dropindex) = [];

    end




    %Load avg. income
    temp = load(fullfile(runDir,dumpDir,'NewOrderPlus.txt'));
    income_avgannual = temp(25)./3;

    % Clear temporary variables
    clearvars filenamedump filename temp delimiter startRow formatSpec fileID dataArray ans;


    %% Calculate subsidization: (q-qfair)*d

    qqd_behavioral = (Price(jj).Goofy - FairPrice.Goofy) .* AmountBorrowed(jj).Goofy;
    qqd_rational = (Price(jj).Realist - FairPrice.Realist) .* AmountBorrowed(jj).Realist;


    % Normalize with average annual income
    qqdbyI_behavioral{jj} = qqd_behavioral./income_avgannual.*100;
    qqdbyI_rational{jj} = qqd_rational./income_avgannual.*100;

end

% calculate joint x-subs
qqdbyI_rational{3} = [qqdbyI_rational{1} ; qqdbyI_rational{2}(1:round(.38/.62*length(qqdbyI_rational{2})))];
qqdbyI_behavioral{3} = [qqdbyI_behavioral{1} ; qqdbyI_behavioral{2}(1:round(.38/.62*length(qqdbyI_behavioral{2})))];

% define bins
bins = linspace(-.102,.102,51);

% winsorize
qqdbyI_realist_wins = qqdbyI_rational{3};
qqdbyI_realist_wins(qqdbyI_realist_wins<min(bins)) = min(bins);
qqdbyI_realist_wins(qqdbyI_realist_wins>=0) = -1e-99;

qqdbyI_goofy_wins = qqdbyI_behavioral{3};
qqdbyI_goofy_wins(qqdbyI_goofy_wins<=0) = 1e-99;
qqdbyI_goofy_wins(qqdbyI_goofy_wins>max(bins)) = max(bins);


% calculate CV(r)

% size sample to represent true population sizes
cut_sample = shareCollege/(1-shareCollege);  %how much college sample to keep

end_col_goofy = round(cut_sample.*length(AmountBorrowed(2).Goofy));
end_col_real = round(cut_sample.*length(AmountBorrowed(2).Realist));

% collect price and debts: NoCol Behav; NoCol Rational; Col Behav; Col Rational
debts = [AmountBorrowed(1).Goofy;AmountBorrowed(1).Realist;AmountBorrowed(2).Goofy(1:end_col_goofy);AmountBorrowed(2).Realist(1:end_col_real)];   
weights = debts ./ sum(debts);    % weighting with amount borrowed

r = (1./[Price(1).Goofy;Price(1).Realist;Price(2).Goofy(1:end_col_goofy);Price(2).Realist(1:end_col_real)] - 1);    % convert q-notation to interest rate

% calculate CV(r)
avg_r= sum(r.*weights);

var_r= sum(weights.*(r-avg_r).^2);

moments.CVr = sqrt(var_r)./avg_r;   

clear AmountBorrowed Price FairPrice qqd_behavioral qqd_rational qqdbyI_rational qqdbyI_behavioral debts weights r

%% Plot Figure 2: Distribution of Cross-Subsidization
Figure2 = figure('Visible', 'off');
histogram(qqdbyI_realist_wins,bins,'EdgeColor','none','Normalization','probability','FaceColor',coloruse(1,:))
yax = ylim;
hold on
histogram(qqdbyI_goofy_wins,bins,'EdgeColor','none','Normalization','probability','FaceColor',coloruse(2,:))
xlabel('\fontsize{14}(q-q_{fair})*d (% of mean income)')
legend('\fontsize{14}Rational','\fontsize{14}Behavioral','Location','NW')
set(gco,'LineWidth',3);
legend boxoff
set(gca,'FontSize',14,'LineWidth',2);
box('off')
ylim(yax)
xt = xticks;
xl = xticklabels;
xl(1) = strcat({'\leq '} , xl(1));
xl(end) = strcat({'\geq '} , xl(end));
xticklabels(xl)


saveas(Figure2,fullfile(outFolder,'Figure2.png'))


%% Append Table 2 with CV calculation
fileID = fopen([outFolder,'Table2.tex'],'a');
fprintf(fileID,'%25s & %15s & %6s   & & %18s & %6.2f   & %6.3f    \\\\ \r\n','','','','CV of Borrowing r', 0.54, moments.CVr);
fprintf(fileID,'%1s \r\n','\bottomrule');
fprintf(fileID,'%1s \r\n','\end{tabular}');
fprintf(fileID,'%1s \r\n','\end{table}');
fclose(fileID);

disp('Success!')