%solveEGVFI.m
% Use a value-function-iteration procedure to obtain the policy function of  the
% Eaton and Gersovitz default model. 
% output:
% dpc is the policy function for debt under continuation. It is an ny-by-nd matrix, where ny is the number of output grid points and nd is the number of debt grid points
% dpix is the policy function for debt under continuation (dpc) but expressed in terms of the index in the debt grid. That is, dpc=d(dpix). 
% vc is the value function under continuation, an ny-by-nd matrix.
% vb value function under bad standing, an ny-by-nd matrix.
% vg value function under good financial standing, an ny-by-nd matrix.
% q is the price of debt, an ny-by-nd matrix.  Note, the rows of q indicate y_t and the columns d_{t+1} (not d_t). 
% tauc capital control tax under continuation, an ny-by-nd matrix. This fiscal instrument is the one that decentralizes the Eaton-Gersovitz model.
% lac, lag, marginal utility of consumption under continuation and good standing, respectively, both  ny-by-nd matrices. 
% Elagp  the expected value of  next period's marginal utility of consumption,  la_{t+1}, conditional on  continuation in t. 
% This object is useful to compute the capital control tax rate, tauc, in period t.

clear all

% outputloss = 'quadratic';  %the options are 'constant' or 'quadratic'

% load specification of the function
load loss_spec

%% Parameters & Grids
%Exogenous Process for Endowment: 
load grid.mat ygrid pai     % ygrid=endowment grid; associated pai=transition probability matrix
ny = numel(ygrid);          % number of grid points for log of ouput
y = exp(ygrid(:));          % level of output

%Calibrated parameters
rstar = 0.01;     % quarterly risk-free interest rate (Chatterjee and Eyigungor, AER 2012). 
theta = 0.0385;   % probability of reentry (USG  and Chatterjee and Eyigungor, AER 2012). 
sig   = 2;        % intertemporal elasticity of consumption substitution
beta  = 0.85;     % discount factor, from Na, Schmitt-Grohe, Uribe, and Yue (2014)

%Output loss function 
if strcmp(outputloss,'constant')
    load cond_av
    a0 =    0;
    a1 = loss; % set it to match conditional average
    a2 = 0;
    clear loss
end
if strcmp(outputloss,'quadratic')
    a0 = 0;
    a1 = -0.35; 
    a2 = ((1-a1)/2)*(1/max(y));
end

yl = max(0,a0+ a1*y + a2*y.^2); % output loss
ya =  y - yl;                   % output in autarky

%debt grid
dupper  = 2.5;
dlower  = 0;
nd      = 300; % # of grid points for debt
d       = linspace(dlower,dupper,nd);
d       = d';

%force the element closest to zero to be exactly zero
[~,nd0] = min(abs(d)); 
d(nd0)  = 0; 

n = ny*nd; % total number of states 

% matrix for  indices of output as a function of the current state (size ny-by-nd)
yix = repmat((1:ny)',1,nd);

% matrix for  indices of debt as a function of the current state (size ny-by-nd)
dix = repmat(1:nd,ny,1);

% Consider a generic n-by-nd matrix Xtry. Each row of Xtry indicates one of
% the n possible states in the current period  and columns correspond to all
% (nd) possible values for debt assumed in the current peirod and due in 
% the next period under continuation. 
% The variable X could be d_t, y_t, d_{t+1}, etc. 
dtry = repmat(d',[ny 1  nd]);
dtry = reshape(dtry,n,nd);

dptryix = repmat(1:nd,n,1);
dptry = d(dptryix);

ytryix = repmat(yix,nd,1);
ytry = y(ytryix);

% AUTARKY
% Consumption under autarky
ca = repmat(ya, [1 nd]);            % consumption of under bad standing
ua = ( ca.^(1-sig)-1)  / (1-sig);   % period utility under autarky

%Initialize the Value functions
vc      = zeros(ny,nd);         % continue repaying
vcnew   = vc;
vg      = zeros(ny,nd);         % good standing
vb      = zeros(ny,nd);         % bad standing
vr      = zeros(ny,nd);         % reentry
dpix    = zeros(ny,nd);         % debt policy function (expressed in indices)  
dpixnew = zeros(ny,nd);
f       = zeros(ny,nd);         % default indicator 1 if default, 0 otherwise; f maps (y_t,d_t)->{1,0}
q       = ones(ny,nd)/(1+rstar);% q is price of debt; it is a function of  (y_t, d_{t+1}) 

dist = 1;
while dist>1e-8

qtry = repmat(q,[nd 1]);
ctry =  dptry .* qtry - dtry + ytry;
utry = (ctry.^(1-sig) -1)  / (1-sig);
utry(ctry<=0) = -inf;

evgptry = pai *  vg;
Evgptry = repmat(evgptry,nd,1);

[vcnew(:), dpixnew(:)] = max(utry+beta*Evgptry,[],2); % maximum value for each row

vbnew = ua+beta*pai*(theta*vr + (1-theta) * vb);

f = vc<vb; 

qnew = (1- pai*f)/(1+rstar);

dist = max(abs(qnew(:)-q(:))) + max(abs(vcnew(:)-vc(:))) + max(abs(vbnew(:)-vb(:))) + max(abs(dpixnew(:)-dpix(:))); 
distance = ['dist = ', num2str(dist)];
disp(distance);

q = qnew;
vc = vcnew;
vb = vbnew;
vg = max(vc,vb);
vr = repmat(vg(:,nd0),1,nd);
dpix = dpixnew;

end %while dist>...


% Policy functions under continuation;

% debt choice under continuation
dpc = d(dpix);

% consumption of tradables  under continuation
I   = sub2ind([ny nd],yix,dpix);
cc  = q(I).*dpc+y(yix)-d(dix); 

% consumption  under good standing
cg = ca.*f + cc.*(1-f); 

% marginal utility of consumption under continuation
lac = cc.^(-sig);

% marginal utility of consumption  in good standing
lag = cg.^(-sig);

% marginal utility of consumption in autarky
laa = ca.^(-sig); 

% marginal utility of consumption under reentry;
lar = repmat(lag(:,nd0),1,nd); 

% Elagp=E_tla_{t+1} if continuation in t
pai*lag;
Elagp =  ans(I);

% capital controls under continuation
beta* Elagp ./ lac ./q(I);
tauc = 1-ans;

% price of debt under autarky
beta*theta*pai*lar + beta*(1-theta)*pai*laa;
qa = ans./laa;

clear ans *try *new *tryix I

eval(['save eg_' outputloss])