function constrained_152018(eqm_selection_criterion,xi,dlower,dupper)
%constrained_152018(eqm_selection_criterion,xi,dlower,dupper) 
%uses a policy-function iteration procedure to approximate the equilibrium of  the open economy with a flow 
%collateral constraint presented in 
%``Multiple Equilibria in Open Economies with  Collateral Constraints,''  
% by Stephanie Schmitt-Grohé and Martín Uribe. 
%Inputs
%eqm_selection_criterion to pick an equilibrium selection criterion, set the  string variable  follows:
%'b' favors the least severe of the two constrained outcomes (point B in the graphs shown in the paper cited above). 
%'c' favors the most severe of the two constrained outcomes (point C in the graphs) 
%'ab' and 'ac' favor unconstrained outcomes (point A in the graphs of the cited paper).  
% 'ab' is a criterion under which if at a given current state  the equilibrium conditions cannot be satisfied with a nonbinding collateral constraint but all equilibrium conditions are satisfied at two binding points, then criterion 'ab'  favor the least severe of the two constrained outcomes (point B in the graphs) . 
%Criterion 'ac' is defined in a similar fashion. Criteria 'ab' and 'ac' deliver the same equilibrium in the calibration considered here. 
%xi intratemporal elasticity of substitution of tradables for nontradables.
%dlower lower bound of debt grid.
%dupper upper bound of debt grid.
%tpmname transition probability matrix of exogenous state
%©  Stephanie Schmitt-Grohé and Martín Uribe, January 2018.

format compact

filename = ['constrained_152018' eqm_selection_criterion];

load tpm.mat yTgrid yNgrid pai

ny = numel(yTgrid); %number of grid points for the joint process of tradable output and  nontradable output

%Calibration
rstar = 0.04; %interest rate
betta = 0.91; %subjective discount factor
a = 0.31;%weight on traded consumption in the CES aggregator 
sigg = 2;%intertemporal elasticity of consumption substitution
kapa = 0.32*(1+rstar); %fraction of output that is pledgeable as collateral

%Grids in levels
yTgridlevel = exp(yTgrid(:)); %grid of level of tradable output
yNgridlevel = exp(yNgrid(:)); %grid of level of country interest rate

%natural debt limit
NdL = min(yTgridlevel)*(1+rstar)/rstar; 

%The debt grid: 
nd = 800;  %number of points

if nargin<2 | isempty(xi)
xi = 0.5;  %Elasticity of substitution between traded and nontraded goods
end

if nargin<3|isempty(dlower)
dlower =0.2;  %lower bound of debt grid
end

if nargin<4|isempty(dupper)
dupper = 1.02*(1+rstar); %upper bound of debt grid
end

%debt grid
dgrid = linspace(dlower,dupper,nd)';
dd = dgrid(2) - dgrid(1); %debt grid step

n = ny*nd; %total number of states

%Let yT, yN, and d  be 3 ny-by-nd matrices containing values  of the current state
yT = repmat(yTgridlevel,1,nd);  
yN = repmat(yNgridlevel,1,nd);  
d = repmat(dgrid',ny,1);
dix = repmat(1:nd,ny,1); %ny-by-nd matrix of indices of the current debt state

%n-by-nd matrix of tradable consumption as a function of the current state (n in total) and  possible choices of next-period debt (nd in total)
%bsxfun(@rdivide,dgrid',1+r(:));
cTtry = bsxfun(@plus,dgrid'/(1+rstar),yT(:)-d(:));

if min(max(cTtry,[],2))<=0
error('Natural debt limit violated')
end

cTtry(cTtry<=0) = NaN; %rule out negative consumption

%n-by-nd matrix of relative price of nontradables as a function of the current state (n in total) and possible values for next-period debt (nd in total)
bsxfun(@rdivide, cTtry, yN(:)); 
ptry  = (1-a)/a*ans.^(1/xi);

%n-by-nd matrix of the value of collateral  as a function of the current state (n in total) and possible values for next-period debt (nd in total)
bsxfun(@times, ptry, yN(:)); 
clear ptry
collateraltry = kapa*bsxfun(@plus,yT(:),ans);

%n-by-nd matrix of the slackness of the collateral constraint  as a function of the current state (n in total) and possible values for next-period debt (nd in total). 
slacktry = bsxfun(@minus,collateraltry,dgrid');
clear collateraltry


%Gridize  values of  slacktry matrix for   which the collateral constraint is binding  by setting them equal to 0
%Also Identify the indices  at which slacktry is nil (binding collateral constraint). 
aux1 = sign(slacktry);
aux2 = sign([slacktry(:,1) slacktry(:,1:end-1)]);
aux3 = sign([ slacktry(:,2:end) slacktry(:,end)]);
bindCix = find(aux1.*aux2==-1 & slacktry<0); %Points of type C in the graphs of the paper cited in the preamble of this program
bindBix = find(aux1.*aux3==-1 & slacktry<0); %Points of type B in the graphs of the paper cited in the preamble of this program
clear aux1 aux2 aux3
slacktry(bindCix) = 0;
slacktry(bindBix) = 0;

%n-by-nd matrix of consumption of the composite good  as a function of the current state (n in total) and possible values for next-period debt (nd in total)
bsxfun(@plus, a * cTtry.^(1-1/xi),(1-a) * yN(:).^(1-1/xi));
ctry = (ans).^(1/(1-1/xi));

%n-by-nd matrix of marginal utility of tradable consumption  as a function of the current state (n in total) and possible values for next-period debt (nd in total)
latry = a*ctry.^(-sigg) .* (cTtry./ctry).^(-1/xi);
clear ctry cTtry 

%n-by-nd matrix of lambda_t/(1+r_t)  as a function of the current state (n in total) and possible values for next-period debt (nd in total)
%LHSEulertry = bsxfun(@rdivide,latry,1+r(:));
%LHSEulertry = latry/(1+rstar);

%initializations
cT0 = min(yT(:))- rstar / (1+rstar)*d;
c0 = (a * cT0.^(1-1/xi) + (1-a) * yN.^(1-1/xi)).^(1/(1-1/xi));
la0 = a*c0.^(-sigg) .* (cT0./c0).^(-1/xi);
la = la0; %marginal utilliy of tradable consumption
laA = zeros(ny,nd); %marginal utility of tradable consumption when unconstrained

dist = inf;

while  dist>0

%This loop computes  pai*la but allowing for the product pai(i,j)*la(j,k)=0 if pai(i,j)=0 and la(j,k)=NaN
aux = zeros(ny,nd);
for i=1:ny
laaux = la;
pai1 = pai(i,:);
find(pai1==0);
laaux(ans,:) = 0;
aux(i,1:nd) = pai1*laaux;
end

%n-by-nd matrix containing values of  mu_t*la_t as a function of the current state (n in total) and possible values of next-period debt (nd in total)


MUtry = latry/(1+rstar) -  betta*repmat(aux,nd,1);

%Gridize iMUtry=0, by identifying sign change along rows and setting them to 0
aux1 = sign(MUtry);
aux2 = sign([MUtry(:,1) MUtry(:,1:end-1)]);
%MUtry(aux1~=aux2 & ~isnan(aux1.*aux2)) = 0;
%alternative, make the last positive MU already 0
MU1ix = find(aux1.*aux2==-1 & MUtry<0);
clear aux1 aux2 
MUtry(MU1ix) = 0;
clear MU1ix

%Find policy function assuming a nonbinding  collateral constraint
[auxA,dpixA] = min(abs(MUtry),[],2);
pickerA = sub2ind([n nd],(1:n)',dpixA(:));  %In general pickerx (in this case x=u) picks the policy choice out of any n-by-nd matrix  containing for each state (n in total) possible policy choices (nd in total)
laA(:) = latry(pickerA); %marginal utility of tradabels when the collateral constraint is  nonbinding 
find(slacktry(pickerA)<0 | isnan(MUtry(pickerA))  );
laA(ans) = NaN; %rule out violations of the collateral constraint and  negative consumption 
pickerAok = find(~isnan(laA));
muA=MUtry(pickerA);
muA(isnan(laA))=NaN; 

%marginal utility of tradables under the assumption that 'c' points are chosen. 
laC = NaN(ny,nd);
find(MUtry>=0 & ~isnan(MUtry));
pickerCok = intersect(ans,bindCix);
[pickerCiok,~] = ind2sub([n nd],pickerCok);
laC(pickerCiok) = latry(pickerCok);

%marginal utility of tradables under the assumption that 'b' points are chosen. 
laB = NaN(ny,nd);
find(MUtry>=0 & ~isnan(MUtry));
pickerBok = intersect(ans,bindBix);
[pickerBiok,ans] = ind2sub([n nd],pickerBok);
laB(pickerBiok) = latry(pickerBok);

laold = la;

%Construct la under the different equilibrium selection criteria

if strcmp(eqm_selection_criterion,'ab')
la = laC;
la(pickerBiok) = laB(pickerBiok);
la(pickerAok) = laA(pickerAok);
end

if strcmp(eqm_selection_criterion,'ac')
la = laB;
la(pickerCiok) = laC(pickerCiok);
la(pickerAok) = laA(pickerAok);
end

if strcmp(eqm_selection_criterion,'b')
la = laA;
la(pickerCiok) = laC(pickerCiok);
la(pickerBiok) = laB(pickerBiok);
end

if strcmp(eqm_selection_criterion,'c')
la = laA;
la(pickerBiok) = laB(pickerBiok);
la(pickerCiok) = laC(pickerCiok);
end

%Compute the distance between the old and new values of la to determine convergence. This is longer than just computing max(abs(la(:)-laold(:))) because NaN enries must be accounted for

find(isnan(la));
aux1 = la; aux1(ans)=0;
find(isnan(laold));
aux2 = laold; aux2(ans)=0;
dist = max(abs(aux1(:)-aux2(:)));
end  %while  dist>0



clear pickerA   pickerAok   pickerBiok   pickerBok   pickerCiok     pickerCok 


%debt policy function in index form (entries take values between 1 and nd)
dpix  = zeros(ny,nd); 
abs(bsxfun(@minus,latry,la(:)));
clear latry
[~,dpix(:)] = min(ans,[],2);
dpix(isnan(la)) = NaN;

picker = sub2ind([n nd],(1:n)',dpix(:));
pickerok = find(~isnan(picker));

%policy function for the slackness of the collateral constraint
slack = nan(ny,nd);
slack(pickerok) = slacktry(picker(pickerok));
clear slacktry

%equilibrium mu
mu = nan(ny,nd);
mu(pickerok) = MUtry(picker(pickerok));
clear MUtry

%Debt policy function
dp = nan(ny,nd);
dp(pickerok) = dgrid(dpix(pickerok));

%Policy functions for cT, c,  and p
cT = yT + dp/(1+rstar) -d;
c =  a * cT.^(1-1/xi)+ (1-a)*yN.^(1-1/xi);
c = c.^(1/(1-1/xi));
p = (1-a)/a *(cT./yN).^(1/xi);

%crisis in the present period
crisis = slack;
crisis(slack==0) = 1;
crisis(slack~=0) = 0;
crisis(isnan(slack)) = nan;

%slope of the RHS of the collateral constraint
slope = kapa /xi ./(1+rstar) .* (p./cT) .* yN;

if xi==0.83
filename = [filename '_xi83'];
end

%eval(['save ' filename '.mat filename nd ny n betta sigg   pai xi a dgrid yTgrid yNgrid yTgridlevel yNgridlevel  laA  la dpix dp v crisis yT yN cT c kapa  picker  dupper dlower dd rstar  distTol d p  slack mu  bindCix bindBix   slope  '])
eval(['save ' filename '.mat filename nd ny n betta sigg   pai xi a dgrid yTgrid yNgrid yTgridlevel yNgridlevel  laA  la dpix dp  crisis yT yN cT c kapa  picker  dupper dlower dd rstar   d p  slack mu  bindCix bindBix   slope  '])


