function  ramsey_1122018(xi,dlower,dupper)
% ramsey_1122018(xi,dlower,dupper) 
% uses a modified value function  iteration procedure to approximate the Ramsey equilibrium 
% of a two-sector open economy in which borrowing is limited by a fraction of output as presented in 
% ``Multiple Equilibria in Open Economies   with Collateral Constraints,''  
% by Stephanie Schmitt-Grohé and Martín Uribe.
%Inputs:
%xi intratemporal elasticity of substitution of tradables for nontradables.
%dlower lower bound of debt grid.
%dupper upper bound of debt grid.
%©  Stephanie Schmitt-Grohé and Martín Uribe, May 2016.

format compact

filename = 'ramsey_1122018';

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<1 | isempty(xi)
xi = 0.5;  %Elasticity of substitution between traded and nontraded goods
end

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

if nargin<3|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 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));
clear cTtry
ctry = (ans).^(1/(1-1/xi));
clear ans

%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)
utry = (ctry.^(1-sigg)-1)/(1-sigg);
clear ctry


%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)]);
bindCix = find(aux1.*aux2==-1 & slacktry<0); %Points of type C in the graphs of the paper cited in the preamble of this program
clear aux2
aux3 = sign([ slacktry(:,2:end) slacktry(:,end)]);
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  aux3
slacktry(bindCix) = 0;
clear bindCix
slacktry(bindBix) = 0;
clear bindBix 


utry(slacktry<0) = NaN;

%initializations
v = zeros(ny,nd); %marginal utilliy of tradable consumption
vnew = zeros(ny,nd);
dpix  = zeros(ny,nd); %debt policy function in index form (entries take values between 1 and nd)
dpix(:) = dix; %set initial guess for debt policy funciton (index form) to current  state 
dpixnew  = zeros(ny,nd); %iupdate  guess 
distTol = 1e-8; %convergence criterion
dist = inf;  %distance criterion

while dist>distTol 

aux = zeros(ny,nd);
for i=1:ny
vaux = v;
pai1 = pai(i,:);
find(pai1==0);
vaux(ans,:) = 0;
aux(i,1:nd) = pai1*vaux;
end

utry +  betta*repmat(aux,nd,1); 
[vnew(:), dpixnew(:)] = max(ans,[],2);
dpixnew(isnan(vnew)) = NaN; 

find(isnan(vnew(:)));
aux1 = vnew; aux1(ans) = 0;
find(isnan(v(:)));
aux2 = v; aux2(ans) = 0;
dist_v = max(abs(aux1(:)-aux2(:)));

find(isnan(dpixnew(:)));
aux1 = dpixnew; aux1(ans) = 0;
find(isnan(dpix(:)));
aux2 = dpix; aux2(ans) = 0;
dist_dpix = max(abs(aux1(:)-aux2(:)));

dist = dist_dpix + dist_v;

xR = numel(find(dpix~=dpixnew)); %number of states that have not yet converged

v = vnew;
dpix = dpixnew;


end %while dist>distTol

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

dp = nan(ny,nd);
dp(pickerok) = dgrid(dpix(pickerok));


%equilibrium slackness
slack = nan(ny,nd);
slack(pickerok) = slacktry(picker(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 = kapa * 1/xi ./(1+rstar) .* (p./cT).* yN; %slope of the RHS of the collateral constraint

la = a*c.^(-sigg) .* (cT./c).^(-1/xi);


%Capital control tax rate
Elap = nan(ny,nd); 
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
repmat(aux,nd,1); 
Elap(pickerok) = ans(picker(pickerok));

tau=nan(ny,nd);
tau(pickerok) = 1- betta*(1+rstar).*Elap(pickerok)./la(pickerok);

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


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