%% Informal Borrowing Completes Formal Borrowing Frontier
% *<https://fanwangecon.github.io/ThaiForInfLuuRobFan/ Go back to LuuRobFan 
% Project Main Page>*
% 
% In <https://fanwangecon.github.io/ThaiForInfLuuRobFan/menu/choice_set.html 
% 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.
%% Formal Borrowing Complemented by Informal Borrowing
% Staying inside the same environment as the one from <https://fanwangecon.github.io/ThaiForInfLuuRobFan/menu/choice_set.html 
% 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 <https://fanwangecon.github.io/ThaiForInfLuuRobFan/menu/choice_set.html 
% 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 <https://fanwangecon.github.io/ThaiForInfLuuRobFan/menu/choice_set.html 
% 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,:))
%% 
%