%plot_policy_dp.m 
%Produces Figure 11 `Debt Policy Functions with Fundamental Uncertainty, \xi=0.5' in 
%``Multiple Equilibria in Open Economies with  Collateral Constraints,''  
% by Stephanie Schmitt-Grohé and Martín Uribe. 
%©  Stephanie Schmitt-Grohé and Martín Uribe, January 2018.

clear all
clf
step = 1; 
rows = 2;
cols = 2;

fname = {'constrained_152018ab',  'constrained_152018b', 'constrained_152018c', 'ramsey_1122018'};
%produced by running simu('constrained_152018ab'), simu('constrained_152018b'), simu('constrained_152018c'), simu(ramsey_1122018')
%the input, constrained_152018ab.mat, is produced by running constrained_152018.m  and by setting
% eqm_selection_criterion = 'ab';
%same for 'ac', 'b', and 'c'. Or ramsey_1122018.m

nf = numel(fname); 

load simu_constrained_152018ab ny nd YT YN yTgridlevel yNgridlevel yT yN a rstar kapa xi dgrid
%produced by running simu('constrained_152018ab'), 

%Find the grid point for the  mean of yT and yN
x = 0*std(YT);
y = mean(YT); 
[zT, zTix] = min(abs(yT(:,1)-(y-x)));
yT_m = yT(zTix); 
x = 0*std(YN);
y = mean(YN); 
[zN, zNix] = min(abs(yN(:,1)-( y-x)));
yN_m = yN(zNix); 
s1 = find(yT(:,1)==yT_m & yN(:,1)==yN_m);

ddpp = zeros(ny,nd,3); 
ccrisis = ddpp; 
sslope = ddpp;
for i = 1:nf
eval(['load simu_' fname{i} '.mat  lad dgrid d dp  crisis '])
ddpp(:,:,i) = dp; 
end

dpA=squeeze(ddpp(s1,:,1));
dpB=squeeze(ddpp(s1,:,2));
dpC=squeeze(ddpp(s1,:,3));
dpR=squeeze(ddpp(s1,:,4));

crit = {'Criterion (A)','Criterion (B)','Criterion (C)', 'Ramsey'};
k=1;
%plot dpA
subplot(rows, cols, k)
%find point of discontinuity of dpA
ixA = find(abs(diff(dpA))==max(abs(diff(dpA))));
%Now plot dpA in two parts
plot(dgrid(1:ixA), dpA(1:ixA), 'g')
hold on 
plot(dgrid(ixA+1:end), dpA(ixA+1:end), 'g')
plot(dgrid,dgrid, ':')
xlim([min(dpA) dgrid(end)])
ylim([min(dpA) dgrid(end)])

title(crit(1))
xlabel('$d_t$', 'Interpreter', 'LaTeX')
ylabel('$d_{t+1}$', 'Interpreter', 'LaTeX')

k=k+step;
%plot dpB
subplot(rows, cols, k)
%find point of discontinuities of dpB
[~,ixB] = sort(abs(diff(dpB(:))));
ixB = ixB(end-1:end); 
%Now plot dpB in three parts
plot(dgrid(1:ixB(1)), dpB(1:ixB(1)), 'b')
hold on 
plot(dgrid(ixB(1)+1:ixB(2)), dpB(ixB(1)+1:ixB(2)), 'b')
plot(dgrid(ixB(2)+1:end), dpB(ixB(2)+1:end), 'b')
plot(dgrid,dgrid, ':')
xlim([min(dpB) dgrid(end)])
ylim([min(dpB) dgrid(end)])

title(crit(2))
xlabel('$d_t$', 'Interpreter', 'LaTeX')
ylabel('$d_{t+1}$', 'Interpreter', 'LaTeX')

k=k+step;
%plot dpC
subplot(rows, cols, k)
%find point of discontinuities of dpB
[~,ix] = sort(abs(diff(dpC(:))));
ix = ix(end); 
%Now plot dpB in three parts
plot(dgrid(1:ix), dpC(1:ix),'c')
hold on 
plot(dgrid(ix+1:end), dpC(ix+1:end),'c')
plot(dgrid,dgrid, ':')
xlim([min(dpC) dgrid(end)])
ylim([min(dpC) dgrid(end)])

title(crit(3))
xlabel('$d_t$', 'Interpreter', 'LaTeX')
ylabel('$d_{t+1}$', 'Interpreter', 'LaTeX')


%plot dpR
%for the Ramsey case find  the largest D such that yT and yN are at the  mean. 
load simu_ramsey_1122018  YT YN D Dp yT yN dgrid
find(YT==yT(zTix)& YN==yN(zNix));
d_high = max(D(ans));
d_low = min(D(ans));

k=k+step;
subplot(rows, cols, k)

%find point of discontinuities of dpB
[~,ix] = sort(abs(diff(dpR(:))));
ix = ix(end); 
%Now plot dpR in  parts
plot(dgrid(1:ix), dpR(1:ix))
hold on 
plot(dgrid(ix+1:end), dpR(ix+1:end))
plot(dgrid,dgrid, ':')
xlim([min(dpR) dgrid(end)])
ylim([min(dpR) dgrid(end)])

title(crit(4))
xlabel('$d_t$', 'Interpreter', 'LaTeX')
ylabel('$d_{t+1}$', 'Interpreter', 'LaTeX')

shg