%TVCCI.M
%This program computes the equilibrium (parameter vectors aL, aH, a1L, and a1H) in a model with an individual time-varying collateral constraint. Collateral is given by  kappa q_tk_{t+1}, where kappa is a parameter, k_t is a fixed factor of production (land) and q_t is the market price of this factor.  The model economy is developed in ``Individual Versus Aggregate Collateral Constraints and the Overborrowing Syndrome,'' by Martin Uribe. 
%(c) Martin Uribe, May 1, 2006.  

clc
clear all
format compact

%Deep parameters and steady-state values
[RSTAR, SIG, OMEGA, BETTA1, ALFA, MU, z, PAIZ, KAPA, Kss,  QQss, Css, Hss, LAss, Ass, Yss, BETTAss] = tvcci_ss;

%Realizations of productivity shock
zL = z(1);
zH = z(2);

DH = 7.25;%upper bound on debt

DL = -6.5; %lower bound on debt

m = 1000; %number of grid points

n = 30; %number of regressors

a = chebygrid(DL,DH,m); %debt grid in levels 

X = chebyreg(m,n); %regressors
D = inv(X'*X);

%Regression Parameters Low state
aL = zeros(n,1);
aL(1) =  log(LAss/RSTAR);

a1L = zeros(n,1);
a1L(1) =  log(1);

%Regression Parameters High state
aH = aL;
a1H = a1L;

%Stopping criterion
dist=1;

%smoothing parameter
lambda = 0.01;

while dist>1e-12

%marginal utility of consumption-to-interest rate ratio
laL_o_rL = exp(X*aL) ;
laH_o_rH = exp(X*aH);

%price of land
qL = exp(X*a1L) ;
qH = exp(X*a1H);

%Obtain remaining endogenous variables
[rL, rH, laL, laH, hL, hH, yL, yH, cL, cH, apL, apH, bettaL, bettaH] = tvcci_variables(laL_o_rL, laH_o_rH, qL, qH, a, DL, DH);

%Lagrange multiplier on collateral constraint
xiL = 1/RSTAR - 1./rL;
xiH = 1/RSTAR - 1./rH;


%Compute the RHS of the consumption Euler equation
WL = chebyreg2(apL,DL,DH,n);
WH = chebyreg2(apH,DL,DH,n);

laLL_o_rLL = exp(WL*aL);
laLH_o_rLH = exp(WL*aH); 
laHL_o_rHL = exp(WH*aL);
laHH_o_rHH = exp(WH*aH); 

qLL = exp(WL*a1L);
qLH = exp(WL*a1H); 
qHL = exp(WH*a1L);
qHH = exp(WH*a1H); 

[rLL, rLH, laLL, laLH, hLL, hLH, yLL, yLH, cLL, cLH, apLL, apLH, bettaLL, bettaLH] = tvcci_variables(laLL_o_rLL, laLH_o_rLH, qLL, qLH, apL, DL, DH);

[rHL, rHH, laHL, laHH, hHL, hHH, yHL, yHH, cHL, cHH, apHL, apHH, bettaHL, bettaHH] = tvcci_variables(laHL_o_rHL, laHH_o_rHH, qHL, qHH, apH, DL, DH);


ElapL = [laLL laLH] * PAIZ(1,:)';
ElapH = [laHL laHH] * PAIZ(2,:)';

rhsL = log (bettaL .* ElapL);
rhsH = log (bettaH .* ElapH);

%Compute the RHS of Tobin's Q Euler equation
xLL = laLL .* (qLL + exp(zL) * ALFA * (hLL/Kss).^(1-ALFA) );

xLH = laLH .* (qLH + exp(zH) * ALFA * (hLH/Kss).^(1-ALFA) );

xHL = laHL .* (qHL + exp(zL) * ALFA * (hHL/Kss).^(1-ALFA) );

xHH = laHH .* (qHH + exp(zH) * ALFA * (hHH/Kss).^(1-ALFA) );

ExpL = [xLL xLH] * PAIZ(1,:)';
ExpH = [xHL xHH] * PAIZ(2,:)';

rhs1L = log(bettaL  ./ (1 - KAPA * xiL) .* ExpL ./ laL);

rhs1H = log (bettaH ./ (1 - KAPA * xiH) .* ExpH ./ laH);

%Store old parameters
aLold = aL;
aHold = aH;
aold = [aLold; aHold];

a1Lold = a1L;
a1Hold = a1H;
a1old = [a1Lold; a1Hold];

Aold = [aold; a1old];

%Compute new parameters by OLS
aLnew = D*X'*rhsL;
aHnew = D*X'*rhsH; 
anew = [aLnew;aHnew];

a1Lnew = D*X'*rhs1L;
a1Hnew = D*X'*rhs1H; 
a1new = [a1Lnew;a1Hnew];

Anew = [anew; a1new];
dist_old = dist;
dist = max(abs(Anew-Aold));

%Update smoothing parameter
if dist<dist_old
lambda = 1 + 0.99 * (lambda-1);
else
lambda = 0.01;
end

%Compute parameters
aL = aLnew*lambda+aLold*(1-lambda);
aH = aHnew*lambda+aHold*(1-lambda);
a1L = a1Lnew*lambda+a1Lold*(1-lambda);
a1H = a1Hnew*lambda+a1Hold*(1-lambda);

dist 
lambda


end