%% Informal Borrowing Completes Formal Borrowing Frontier
% Go back to LuuRobFan Project Main Page
% In Borrowing Choice Set Discretization in a Deterministic Two Period Model, we discussed the effects of discretization on optimal choices and wealth/welfare. If informal borrowing options exist, they will help complete the choice set.
% Table of Contents
% Formal Borrowing Complemented by Informal Borrowing
% Choice Frontier and Informal Borrowing
% Optimal Formal and Informal Choices
% Expenditure Minimization to Calculate Wealth Loss due to Discretization
% Graphical Results
% Grid Constrained Borrowing Choice Set at Four Wealth Levels
% Summary Tabulate Graph Results
% Shifting Wealth and Wealth Loss Due to Discretization
% Fraction of Total Loan that is Informal
% Tabulate All Results
%
%% Formal Borrowing Complemented by Informal Borrowing
% Staying inside the same environment as the one from Borrowing Choice Set Discretization in a Deterministic Two Period Model, we have:
% Exogenous endowments in both periods: z_1, z_2
% Discrete Formal Borrowing Grid with equi-distance gaps: \gamma
% Discrete Formal borrowing choices: b^F_i, where i is the grid index
% Continuous Informal borrowing choices: b^I
% Formal and informal borrowing have different interest rates
% If r^I < r^F, will only borrow informally.
% If r^I \ge r^F, could borrow only formally or jointly formally and informally
%
% Optimization with Discrete Formal and Continuous Informal Choices
% Ignore savings for now, assume z_2 is sufficiently larger than z_1 that households want to borrow if facing the unconstrained problem.
% The new maximization problem can be written as follows, with discrete formal borrowing and continuous informal borrowing choices:
% \max_{i\in\left\{ 0,..., n\left(z_2, \gamma\right) \right\} } \left\{ \max_{b^I \in \left[ 0, B^{I}_i\left(b^F_i\right) \right)} \left\{ \log \left(z_1 - b^F_i - b^I\right) + \beta \log(z_2 + b^F_i \left(1+r^F\right) + b^I \left(1+r^I\right)) \right\} \right\}
% where: B^{I}_i\left(b^F_i\right) = \frac{z_2 + b_i^F\left(1+r^F\right)}{1+r^I}
% and b^F_i is from sequence \left\{ b^F_i  \right\}_{i=0}^{n\left(z_2, \gamma\right)}
% where b^F_i = - \gamma \cdot i, and n\left(z_2, \gamma\right) = \lfloor_{\text{floor}}\frac{z_2}{ \left(1+r^F\right) \gamma}\rfloor
%
%% Choice Frontier and Informal Borrowing
% The choice set above ignores that some of the informal choice points would be within the choice frontier. To only consider the informal borrowing choices that are on the frontier:
% \widehat{B^{I}_i}\left(b^F_i\right) = -\gamma\left(\frac{1+r^F}{1+r^I}\right)
% In solving the model, we will use \widehat{B^{I}_i}\left(b^F_i\right), rathern than B^{I}_i\left(b^F_i\right).
% Note that with just discrete formal borrowing and continuous informal borrowing, the choice frontier is not "complete" yet, as can be seen below. To close the choice frontier, savings is also required.
%
%% Optimal Formal and Informal Choices
% We can find the optimal choices using the codes below. For the graphs and simulations below, parameters are kept consistent with Borrowing Choice Set Discretization in a Deterministic Two Period Model.

% Parameters
z1 = 0.2;
z2 = 2;
beta = 1;

% Informal and Formal Intrest Rate
rI = 0.50; % 50 percent interest rate
rF = 0;

% There is a grid of b choices
bF_grid_gap = 0.6;
bI_eps = 0.01; % less than
bI_grid_n = 50; % this should be very dense to very close to exact/continuous choice

% Same Utility Function as Before, but now b=bINF + bFOR
f_utility = @(b) log(z1-b) + beta*log(z2+b*(1+r));

% This is the max natural borrowing constraint
b_max_borrow = z2/(1+rF); 
% The is the Formal Borrowing Choice Grid
bF_grid = (-1)*(0:bF_grid_gap:b_max_borrow);
bF_grid_n = length(bF_grid);

% Show Choice Grids
disp(table(bF_grid));

% Informal Borrowing max for each grid point derived before 
% bI between [0, -1*GAP*((1+rF)/(1+rI))
bI_min = 0;
bI_max = (-1)*bF_grid_gap*((1+rF)/(1+rI));
bI_grid = linspace(bI_min, bI_max, bI_grid_n);
bI_grid_gap = abs(bI_grid(2)-bI_grid(1));
bI_grid(bI_grid_n) = bI_grid(bI_grid_n); % less than

% mesh formal and informal borrowing grids
[bI_mesh, bF_mesh] = meshgrid(bI_grid, bF_grid);

% for the largest bF, BI to b_grid_gap would exceed Natural Borrowing
bFmax_bI_too_much_idx = (bF_mesh(bF_grid_n)*(1+rF) + bI_grid*(1+rI) <= -z2);
bFmax_bI_max = (-1)*(z2 + bF_mesh(bF_grid_n)*(1+rF))/(1+rI) + bI_grid_gap*bI_eps;
bI_mesh(bF_grid_n, bFmax_bI_too_much_idx) = bFmax_bI_max;

% BFI aded up together
bFI_principle = bF_mesh + bI_mesh;
bFI_interests = bF_mesh*rF + bI_mesh*rI;

% C1 and C2 given just Formal Choices
c1_bF = z1 - bF_grid;
c2_bF = z2 + bF_grid.*(1+rF);
% C1 and C2 given B for + Inf choices
c1_bFI = z1 - bFI_principle;
c2_bFI = z2 + bFI_principle + bFI_interests;

% Utility Along Choice Grid
Utility_Grid = log(c1_bFI) + beta*log(c2_bFI);

% Optimal Choice and Value at Choices
[val_og, max_lin_idx] = max(Utility_Grid(:));
[max_bF_idx, max_bI_idx] = ind2sub(size(Utility_Grid),max_lin_idx);
c1_bFI_opti = c1_bFI(max_bF_idx, max_bI_idx);
c2_bFI_opti = c2_bFI(max_bF_idx, max_bI_idx);
bI_bFI_opti = bI_mesh(max_bF_idx, max_bI_idx);
bF_bFI_opti = bF_mesh(max_bF_idx, max_bI_idx);

% Optimal Choices Display
opti_choices = table(bF_bFI_opti, bI_bFI_opti, c1_bFI_opti, c2_bFI_opti, val_og);
disp(opti_choices);

%
%% Expenditure Minimization to Calculate Wealth Loss due to Discretization
% Similar to before, we solve the expenditure minimization problem below: If only the formal option is available, but it is not discrete, what would be the wealth required to achieve the same level of utility as achieved when discrete formal and continuous informal choices are available?

% Solve a Expenditure Minimization Problem and Calculate Wealth Loss (with rF)
syms c1 c2 lambda
lagrangian = (c2 + (1+rF)*c1 - lambda*(log(c1) + beta*log(c2) - val_og));
solu_min = solve(diff(lagrangian, c1)==0,...
                 diff(lagrangian, c2)==0,...
                 diff(lagrangian, lambda)==0,...
                 c1, c2, lambda, 'Real', true);
solu_min_c1 = double(solu_min.c1);
solu_min_c2 = double(solu_min.c2);
wealth_grid = solu_min_c1*(1+rF) + solu_min_c2;
wealth_loss = 1- (wealth_grid/(z1*(1+rF) + z2));
wealth_loss_percent = round(wealth_loss*10000)/100;

% Show Expenditure Minimization Table
exp_min_table = table(solu_min_c1, solu_min_c2, wealth_loss, wealth_loss_percent);
disp(exp_min_table);

%
%% Graphical Results
% Show results graphically. The Graph below will show these:
% Endowment Point: where z_1 and z_2 are
% Unconstrained (complete budget) Budget: Budget given interest rate, frontier when any borrowing and savings are allowed
% Unconstrained (complete budget) Optimal Choice: optimal unconstrained choice
% B Formal Grid Choice Set: Borrow Grid Constrained Choice Set
% B Informal Cts Choice Set: Informal Choices combined with each formal grid point
% For+Inf  Expenditure Minimizing Budget: Given For+Inf optimal value, what is the budget from the expenditure minimization problem
% For+Inf Optimal Value: Given For+Inf  optimal choice, what is the utility.
% For+Inf Optimal Choice: Optimal For+Inf choices

% For Graphing
syms c1
% The grid constrained value's unconstrained optimal budget line
f_budget_grid_constrained = solu_min_c1*(1+rF) + solu_min_c2 - c1*(1+rF);
% Indifference for grid constrained value
f_indiff_grid_constrained = exp((val_og-log(c1))/(beta));

% Graphing
figure();
hold on;

% Endowment Point
scatter(z1, z2, 350, 'k', 's');
plot(linspace(0,z1,10),ones(10,1) * z2, 'k--', 'HandleVisibility','off');
plot(ones(10,1) * z1, linspace(0,z2,10), 'k--', 'HandleVisibility','off');

% Unconstrained Budget Line
syms c1
f_budget = z1*(1+rF) + z2 - c1*(1+rF);
fplot(f_budget, [0, (z1 + z2/(1+rF))*1.25], 'b--');

% Unconstrained Optimal Point
syms b
b_opti = double(solve(diff(log(z1-b) + beta*log(z2+b*(1+rF)), b)==0,...
                      b, 'Real', true));
c1_opti = z1 - b_opti;
c2_opti = z2 + b_opti*(1+rF);
scatter(c1_opti, c2_opti, 200, 'k', 'd');

% Constrained Lines
scatter(c1_bF(:), c2_bF(:), 100, 'k', 'filled', 'c');
scatter(c1_bFI(:), c2_bFI(:), 10, 'k', 'x');
fplot(f_budget_grid_constrained, [0, (z1 + z2/(1+rF))*1.25], 'b-');
fplot(f_indiff_grid_constrained, [0, (z1 + z2/(1+rF))*1.25], 'r-');

% Constrained Optimal Choice
scatter(c1_bFI_opti, c2_bFI_opti, 250, 'r', 'd');

% Labeling
ylim([0, (z1 + z2/(1+rF))*1.25])
title({[num2str(wealth_loss_percent) '% Wealth Loss Due to Discretization'],...
       ['B Grid Gap = ' num2str(bF_grid_gap) ', beta = ' num2str(beta)...
        ', z1 = ' num2str(z1) ', z2 = ' num2str(z2)...
        ', rF = ' num2str(rF) ', rI = ' num2str(rI)],...
       ['Unconstrained: C1*= ' num2str(c1_opti) ', C2*=' num2str(c2_opti)...
        ', B*=' num2str(b_opti)],...
       ['Constrained: C1*= ' num2str(c1_bFI_opti) ', C2*=' num2str(c2_bFI_opti),...
        ', BF*=' num2str(bF_bFI_opti) ', BI*=' num2str(bI_bFI_opti)]});
    
xlabel('consumption today');
ylabel('consumption tomorrow');
legend({'Endowment',...
        'Complete Budget',...
        'Complete Optimal Choice',...
        'B For. Grid Choice Set',...
        'B Inf. Cts Choice Set',...
        'For+Inf Budget with Wealth Loss',...
        'For+Inf Optimal Value',...
        'For+Inf Optimal Choice'...
        })
grid on;

%
%% Grid Constrained Borrowing Choice Set at Four Wealth Levels
% In these two period models, wealth is z_1\cdot (1+r) + z_2. Similar to before, I will now shift wealth and study in several examples of how households' optimal formal and informal choice combinations change. For different households, the relative ratio of z_1 to z_2 are fixed, so their optimal consumption allocation shares should be the same. However, that is not the case here given formal borrowing discretization.
% In most models of formal and informal borrowing, there is a sorting map in which households with varying wealth and productivity sort into formal and informal borrowing choices. We will generate this here as well.
% Following Borrowing Choice Set Discretization in a Deterministic Two Period Model:
% We set \beta=0.95, and r^F = 0.05
% Borrow grid gap is still: 0.6
% Shift Overall Wealth level: z_1 \cdot (1+r) + z_2, from equal to the borrow grid gap to six times the borrow grid gap.
% In the Four Figures Below, when r_^I = 0.375:
% \text{wealth}=0.6: informal borrowing only
% \text{wealth}=1.6: formal borrowing only
% \text{wealth}=2.6: joint formal borrowing and informal borrowing
% \text{wealth}=3.6: formal borrowing only again
% I first show four graphs, than a table summarizing results from the four graphs.

% Wealth Shares and Discount
z1_share = 0.15;
z2_share = 1 - z1_share;
beta = 0.95;

% There is a grid of b choices
bF_grid_gap = 0.6;
bI_eps = 0.01; % less than

% Informal and Formal Intrest Rate
rF = 0.05; % 10 percent interest rate
rI_grid_n = 3;
rI_grid = linspace(rF + 0.05, 0.65, rI_grid_n); % 5 levels of interest rates
rI_grid_graph = [2];

% Wealth Grid
wealth_grid_min = bF_grid_gap;
wealth_grid_max = bF_grid_gap*6;
wealth_grid_n = 500;
wealth_grid = linspace(wealth_grid_min, wealth_grid_max, wealth_grid_n);
wealth_grid_graph = [1 round(wealth_grid_n/3) round((wealth_grid_n)*2/3) wealth_grid_n];

% Store Wealth Loss Along Grid
wealth_loss_percent_mat = zeros(rI_grid_n, wealth_grid_n);
% Store Wealth Loss Along Grid
opti_borrow_mat = zeros(rI_grid_n, wealth_grid_n);
% Informal Fraction of Total Borrowing
opti_inf_frac_mat = zeros(rI_grid_n, wealth_grid_n);

% Main Graphs Table Results Collection
table_results_summary = [];
table_row_names_summary = {};
graph_counter_summary = 0;

% All Results Store
table_results_all = [];
table_rows_all = {};
graph_counter_all = 0;

% Graphing
for rI_grid_i=1:1:length(rI_grid)
    rI = rI_grid(rI_grid_i);
    for wealth_i=1:1:length(wealth_grid)
        
        % Numbers
        z1 = z1_share * wealth_grid(wealth_i);
        z2 = z2_share * wealth_grid(wealth_i);           
        
        % Same Utility Function as Before, but now b=bINF + bFOR
        f_utility = @(b) log(z1-b) + beta*log(z2+b*(1+r));
        
        % This is the max natural borrowing constraint
        b_max_borrow = z2/(1+rF); 
        % The is the Formal Borrowing Choice Grid
        bF_grid = (-1)*(0:bF_grid_gap:b_max_borrow);
        bF_grid_n = length(bF_grid);
        
        % Informal Borrowing max for each grid point derived before 
        % bI between [0, -1*GAP*((1+rF)/(1+rI))
        bI_min = 0;
        bI_max = (-1)*bF_grid_gap*((1+rF)/(1+rI));
        bI_grid_n = bF_grid_gap/0.01;
        bI_grid = linspace(bI_min, bI_max, bI_grid_n);
        bI_grid_gap = abs(bI_grid(2)-bI_grid(1));
        bI_grid(bI_grid_n) = bI_grid(bI_grid_n); % less than
        
        % mesh formal and informal borrowing grids
        [bI_mesh, bF_mesh] = meshgrid(bI_grid, bF_grid);
        
        % for the largest bF, BI to b_grid_gap would exceed Natural Borrowing
        bFmax_bI_too_much_idx = (bF_mesh(bF_grid_n)*(1+rF) + bI_grid*(1+rI) <= -z2);
        bFmax_bI_max = (-1)*(z2 + bF_mesh(bF_grid_n)*(1+rF))/(1+rI) + bI_grid_gap*bI_eps;
        bI_mesh(bF_grid_n, bFmax_bI_too_much_idx) = bFmax_bI_max;
        
        % BFI aded up together
        bFI_principle = bF_mesh + bI_mesh;
        bFI_interests = bF_mesh*rF + bI_mesh*rI;
        
        % C1 and C2 given just Formal Choices
        c1_bF = z1 - bF_grid;
        c2_bF = z2 + bF_grid.*(1+rF);
        % C1 and C2 given B for + Inf choices
        c1_bFI = z1 - bFI_principle;
        c2_bFI = z2 + bFI_principle + bFI_interests;
        
        % Utility Along Choice Grid
        Utility_Grid = log(c1_bFI) + beta*log(c2_bFI);
        
        % Optimal Choice and Value at Choices
        [val_og, max_lin_idx] = max(Utility_Grid(:));
        [max_bF_idx, max_bI_idx] = ind2sub(size(Utility_Grid),max_lin_idx);
        c1_bFI_opti = c1_bFI(max_bF_idx, max_bI_idx);
        c2_bFI_opti = c2_bFI(max_bF_idx, max_bI_idx);
        bI_bFI_opti = bI_mesh(max_bF_idx, max_bI_idx);
        bF_bFI_opti = bF_mesh(max_bF_idx, max_bI_idx);
        
        % Solve a Expenditure Minimization Problem and Calculate Wealth Loss (with rF)
        syms c1 c2 lambda
        lagrangian = (c2 + (1+rF)*c1 - lambda*(log(c1) + beta*log(c2) - val_og));
        solu_min = solve(diff(lagrangian, c1)==0,...
                         diff(lagrangian, c2)==0,...
                         diff(lagrangian, lambda)==0,...
                         c1, c2, lambda, 'Real', true);
        solu_min_c1 = double(solu_min.c1);
        solu_min_c2 = double(solu_min.c2);
        wealth_lvl = solu_min_c1*(1+rF) + solu_min_c2;
        wealth_loss = 1- (wealth_lvl/(z1*(1+rF) + z2));
        wealth_loss_percent = round(wealth_loss*10000)/100;
        
        % Store Results
        wealth_loss_percent_mat(rI_grid_i, wealth_i) = wealth_loss_percent;
        opti_borrow_mat(rI_grid_i, wealth_i) = bI_bFI_opti + bF_bFI_opti;
        opti_inf_frac_mat(rI_grid_i, wealth_i) = bI_bFI_opti/(bI_bFI_opti + bF_bFI_opti);
        
        
        % Counting results
        graph_counter_all = graph_counter_all + 1;            
        
        % Collect Results in Tables for these cases
        bF_bFI_o = str2num(sprintf('%4.3f', (bF_bFI_opti)));
        bI_bFI_o = str2num(sprintf('%4.3f', (bI_bFI_opti)));
        c1_bFI_o = str2num(sprintf('%4.3f', (c1_bFI_opti)));
        c2_bFI_o = str2num(sprintf('%4.3f', (c2_bFI_opti)));
        val_o = str2num(sprintf('%4.3f', (val_og)));
        inc_loss_perc_o = str2num(sprintf('%4.3f', (wealth_loss_percent)));
        opti_choices_table = table(bF_bFI_o, bI_bFI_o,...
                                   c1_bFI_o, c2_bFI_o,...
                                   val_o, inc_loss_perc_o);
        table_results_all = [table_results_all; opti_choices_table];
        table_row_names_all{graph_counter_all} = ...
            ['wealth:' sprintf('%4.3f', wealth_grid(wealth_i)) ' rI:' num2str(rI)];
        
        % Graph Budget and Indiff for a Subset
        if (ismember(wealth_i, wealth_grid_graph) && ismember(rI_grid_i, rI_grid_graph))
            
            % Counting graphs
            graph_counter_summary = graph_counter_summary + 1;            
            % Collect Results in Tables for these cases
            table_results_summary = [table_results_summary; opti_choices_table];
            table_row_names_summary{graph_counter_summary} =...
                ['wealth = ' sprintf('%4.3f', wealth_grid(wealth_i))];
            
            % Graphing
            figure();
            hold on;
            
            % Endowment Point
            scatter(z1, z2, 350, 'k', 's');
            plot(linspace(0,z1,10),ones(10,1) * z2, 'k--', 'HandleVisibility','off');
            plot(ones(10,1) * z1, linspace(0,z2,10), 'k--', 'HandleVisibility','off');
            
            % Unconstrained Budget Line
            syms c1
            f_budget = z1*(1+rF) + z2 - c1*(1+rF);
            fplot(f_budget, [0, (z1 + z2/(1+rF))*1.25], 'b--');
            
            % Unconstrained Optimal Point
            syms b
            b_opti = double(solve(diff(log(z1-b) + beta*log(z2+b*(1+rF)), b)==0,...
                                  b, 'Real', true));
            c1_opti = z1 - b_opti;
            c2_opti = z2 + b_opti*(1+rF);
            scatter(c1_opti, c2_opti, 200, 'k', 'd');
            
            % Constrained Lines
            % For Graphing
            syms c1
            % The grid constrained value's unconstrained optimal budget line
            f_budget_grid_constrained = solu_min_c1*(1+rF) + solu_min_c2 - c1*(1+rF);
            % Indifference for grid constrained value
            f_indiff_grid_constrained = exp((val_og-log(c1))/(beta));
            
            scatter(c1_bF(:), c2_bF(:), 100, 'k', 'filled', 'c');
            scatter(c1_bFI(:), c2_bFI(:), 10, 'k', 'x');
            fplot(f_budget_grid_constrained, [0, (z1 + z2/(1+rF))*1.25], 'b-');
            fplot(f_indiff_grid_constrained, [0, (z1 + z2/(1+rF))*1.25], 'r-');
            
            % Constrained Optimal Choice
            scatter(c1_bFI_opti, c2_bFI_opti, 250, 'r', 'd');
            
            % Labeling
            ylim([0, (z1 + z2/(1+rF))*1.25])
            title({[num2str(wealth_loss_percent) '% Wealth Loss Due to Discretization'],...
                   ['W = ' num2str(wealth_grid(wealth_i))...
                    ', z1 = ' num2str(z1),...
                    ', z2 = ' num2str(z2) ],...
                   ['BF Grid Gap = ' num2str(bF_grid_gap),...
                    ', beta = ' num2str(beta) ...
                    ', rF = ' num2str(rF) ', rI = ' num2str(rI)],...
                   ['Unconstrained: C1*= ' num2str(c1_opti)...
                    ', C2*=' num2str(c2_opti)...
                    ', B*=' num2str(b_opti)],...
                   ['Constrained: C1*= ' num2str(c1_bFI_o)...
                    ', C2*=' num2str(c2_bFI_o)...
                    ', BF*=' num2str(bF_bFI_o)...
                    ', BI*=' num2str(bI_bFI_o)]});
            xlabel('c t');
            ylabel('c t+1');            
            legend({'Endowment',...
                    'Complete Budget',...
                    'Complete Optimal Choice',...
                    'B For. Grid Choice Set',...
                    'B Inf. Cts Choice Set',...
                    'For+Inf Budget with Wealth Loss',...
                    'For+Inf Optimal Value',...
                    'For+Inf Optimal Choice'...
                    })            
            grid on;
        end
    end
end

%
%% Summary Tabulate Graph Results
% The Table below corresponds to the four main graphs above:

table_results_summary.Properties.RowNames = table_row_names_summary;
disp(table_results_summary)

%
%% Shifting Wealth and Wealth Loss Due to Discretization
% Plot also based on stored results, given the parameters here, variations in wealth loss along the total 2 period wealth grid.
% If Informal interest rate is high, informal options do not improve things
% If r^I = r^F, then households can borrow any amount, there are no loss to discretization.
% In the results below, we can see that, along the x-axis:
% For the lowest wealth levels, households borrow only informally
% As wealth increases, households borrow only formally
% Then households borrow jointly formally and informally
% Then just formally again
% When the informal interest rate is high, only small segments of wealth benefit from informal borrowing.

figure();
hold on;
    
% Left Axis Wealth Loss
yyaxis left

% find y max 
inc_loss_max = max(max(wealth_loss_percent_mat));
opti_borr_min = min(min(opti_borrow_mat));

% Plot Each Interest Rate Group One by One
r_legend = {};
legend_counter = 0;
for rI_grid_i=length(rI_grid):-1:1
    rI = rI_grid(rI_grid_i);
    legend_counter = legend_counter + 1;
    r_legend{legend_counter} = ['Informal r = ' num2str(rI)];    
    wealth_loss_percent_grid = wealth_loss_percent_mat(rI_grid_i, :);
    
    if (rI_grid_i == 1)
        plot(wealth_grid, wealth_loss_percent_grid, ':' , 'LineWidth', 3);
    elseif (rI_grid_i == 2)
        plot(wealth_grid, wealth_loss_percent_grid, '--', 'LineWidth', 2);
    elseif (rI_grid_i == 3)
        plot(wealth_grid, wealth_loss_percent_grid, '-', 'LineWidth', 1);
    end
    
%     plot(wealth_grid, ones(size(wealth_grid)) * 0, 'k--');
    ylim([-(inc_loss_max)/8 inc_loss_max*9/8]);
    xlim([wealth_grid_min wealth_grid_max]);
    
    title({['% Wealth Loss to Borrowing Grid Discretization'],...
           ['B Grid Gap = ' num2str(bF_grid_gap) ...
            ', beta = ' num2str(beta) ', rF = ' num2str(rF) ],...
           ['z1 wealth share = ' num2str(z1_share) ...
            ', z2 wealth share = ' num2str(z2_share)]})
    xlabel('wealth = z1*(1+r) + z2 ');
    ylabel({['Wealth Loss Percentage'],...
            ['Computed From Expenditure Minimization'],...
            ['Given Value at Optimal Constrained-For+Inf Choice'],...
            ['Compare to Value at Optimal Un-Constrained Choice']});
    grid on;
    grid minor;
end

% Right Axis Optimal Choices
yyaxis right
% Plot Each Interest Rate Group One by One
for rI_grid_i=1:1:length(rI_grid)
    rI = rI_grid(rI_grid_i);
    opti_borrow_grid = opti_borrow_mat(rI_grid_i, :);
    
    if (rI_grid_i == 1)
        plot(wealth_grid, (-1)*opti_borrow_grid, ':', 'LineWidth', 3);
    elseif (rI_grid_i == 2)
        plot(wealth_grid, (-1)*opti_borrow_grid, '--', 'LineWidth', 2);
    elseif (rI_grid_i == 3)
        plot(wealth_grid, (-1)*opti_borrow_grid, '-',  'LineWidth', 1);
    end    
    
    ylim([-(-1*opti_borr_min)/8 (-1*(opti_borr_min+opti_borr_min/8))]);
    xlim([wealth_grid_min wealth_grid_max]);
    ylabel({['Optimal Discrete Grid Borrowing Choices * (-1)'],...
            ['The sum of Formal and Informal Choices']});            
end
legend(r_legend, 'Location','north');

%
%% Fraction of Total Loan that is Informal
% Based on the same data as above, we can distinguish formal and informal borrowing by calculating the fraction of total borrowing that is informal borrowing.
% Lowest wealth households only borrow informally
% Then there are subsets of households whose unconstrained optimal choices are close to the formal grid points, hence they do not borrow informally
% For households whose unconstrained optimal choices are further away from the formal grid points, they borrow informally as well
% When the informal borrowing interest rate increases, informal borrowing share decreases.

figure();
hold on;
   
% Plot Each Interest Rate Group One by One
r_legend = {};
for rI_grid_i=1:1:length(rI_grid)
    rI = rI_grid(rI_grid_i);
    r_legend{rI_grid_i} = ['Informal r = ' num2str(rI)];
    opti_inf_frac = opti_inf_frac_mat(rI_grid_i, :);
    
    if (rI_grid_i == 1)
        plot(wealth_grid, opti_inf_frac, ':', 'LineWidth', 3);
    elseif (rI_grid_i == 2)
        plot(wealth_grid, opti_inf_frac, '--', 'LineWidth', 2);
    elseif (rI_grid_i == 3)
        plot(wealth_grid, opti_inf_frac, '-',  'LineWidth', 1);
    end    
    
    ylim([0 1]);
    xlim([wealth_grid_min wealth_grid_max]);
    
    title({['Informal Fraction of Total Borrowing'],...
           ['B Grid Gap = ' num2str(bF_grid_gap) ...
            ', beta = ' num2str(beta) ', rF = ' num2str(rF) ],...
           ['z1 wealth share = ' num2str(z1_share) ...
            ', z2 wealth share = ' num2str(z2_share)]})
    xlabel('wealth = z1*(1+r) + z2 ');    
    ylabel({['Informal Borrowing Fraction of Total Borrowing']});
        
    legend(r_legend, 'Location','north');
    grid on;
    grid minor;    
end

%
%% Tabulate All Results
% The Tables below tabulate results, showing 100 rows from the graph above:

table_results_all.Properties.RowNames = table_row_names_all;
[rows, cols] = size(table_results_all);
% Only show 100 rows
rows_display = floor(linspace(1, rows, 100));
disp(table_results_all(rows_display,:))

%