%CCRC_SIMU.M
%This program produces a simulation of the equilibrium dynamics
function [dm, A, Z, R, LA, H, Y, C, BETTA, u] = ccrc_simu(zL, zH, DL, DH, aL, aH, n, RSTAR, SIG, OMEGA, BETTA1, ALFA, PAIZ, KAPA, Kss,QQss);

%The following command sets a different state of the seed for random numbers each time this program is used
rand('state',sum(100*clock)); 

%Length of simulation in periods
T = 5000;

%Number of initial periods to be discarded in computing the denHaan-Marcet (RES, 1994) statistic
Tcut = 1000;

%Number of usable periods for computing the denHaan-Marcet statistic
Tleft = T-Tcut;

%Initial value of debt
A0 = 5.10;%initial debt position 

b=0;h=0;aa=0;%initial conditions for construction of den Haan and Marcet (RES, 1994) accuracy test

%initial state of technology
x = rand<0.5; 
sL = x; %low state
sH = 1-x; %high state

for t=1:T

%Debt
A(t,1) = A0;


%Productivity shock
x = rand<PAIZ(1,1);

sL = sL * x + (1-sL) * (1-x);
sH = sH * x + (1-sH) * (1-x);

Z(t,1) = sL * zL + sH * zH;

%coefficient for parameterizing the variable r_t*c_t
a = sL * aL + sH * aH;
 
%Interest rate under the conjectue of debt limit not binding
R(t,1) = RSTAR;

%Transformed state
X(t,1:n) = chebyreg2(A(t,1),DL,DH,n);

%consumption
C(t,1) = exp(X(t,:)*a)  / R(t,1);

%Hous
H(t,1) = ( (1-ALFA) * exp(Z(t)) * Kss^(ALFA) ) .^(1/(OMEGA+ALFA-1));

%Output
Y(t,1) = exp(Z(t)) * Kss^ALFA * H(t).^(1-ALFA);

%Next period's debt
Ap(t,1) = R(t) .* (A(t) + C(t) - Y(t));

%Impose limits on debt
Ap(t,1) = max(Ap(t),DL);
Ap(t,1) = min(Ap(t),  KAPA*QQss*Kss);

%If debt limit is binding
if Ap(t)==KAPA*QQss*Kss

%adjust interest rate
R(t,1) = (Ap(t)-exp(X(t,:)*a)) / (A(t)-Y(t));
end

%adjust consumption
C(t) = Ap(t,1) / R(t) - A(t)  + Y(t);

%Marginal utility of consumption
LA(t,1) = (C(t) - H(t)^OMEGA/OMEGA)^(-SIG);

%Discount factor
BETTA(t,1) = (1 + C(t) - H(t)^OMEGA/OMEGA)^(-BETTA1);

A0 = Ap(t);

%Prepare for computing denHaan-Marcet (RES, 1994) statistic
if t>Tcut;
tt = t-Tcut;
u(tt,1) = LA(t-1)-R(t-1).* BETTA(t-1).* LA(t);
h = [A(t-5:t-1)' 1 Z(t-5:t-1)'];
b = b + u(tt)*h/Tleft  ;
aa = aa + u(tt)^2 * h'*h/Tleft;
end %if t>Tcut
end %for t=1:T

%denHaan-Marcet (RES, 1994) statistic
dm = Tleft * b*inv(aa)*b';
dm = Tleft * b*inv(aa)*b';