%NCC_SIMU.M
function [dm, A, Z, R, LA, H, Y, C, BETTA, u] = ncc_simu(zL, zH, DL, DH, aL, aH, n, RSTAR, SIG, OMEGA, BETTA1, ALFA, PAIZ, Kss);
%This program produces a simulation of the equilibrium dynamics. 
%Inputs are obtained by running NCC.M

rand('state',sum(100*clock)); %sets a different state each time.
T = 5000;
Tcut = 1000;%number of periods thrown out in computing the test

Tleft = T-Tcut;

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
x = rand<0.5; 
sL = x; %low state
sH = 1-x; %high state

for t=1:T

A(t,1) = A0;

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;

a = sL * aL + sH * aH;
 
R(t,1) = RSTAR;

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

LA(t,1) = exp(X(t,:)*a) .* R(t,1);

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

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

C(t,1) = LA(t).^(-1/SIG) + H(t)^OMEGA/OMEGA;

Ap(t,1) = RSTAR .* (A(t) + C(t) - Y(t));

Ap(t,1) = max(Ap(t),DL);
Ap(t,1) = min(Ap(t),DH);

C(t) = Ap(t,1) / RSTAR - A(t)  + Y(t);

LA(t,1) = (C(t) - H(t)^OMEGA/OMEGA)^(-SIG);

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;
u = 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*h/Tleft  ;
aa = aa + u^2 * h'*h/Tleft;
end %if t>Tcut
end %for t=1:T

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