function discs=NKACCModel(x)
global sigma psi chi alpha beta delta gamma T gs StatVars epsilon theta phi lambda kappa rBar yBar IfSave

nCore=8;
nTerm=1;  %  Number of core variables with terminal conditions
discs=zeros(nCore*(T-1)-nTerm,1);
yLoc=1;
pLoc=yLoc+T-1;
hLoc=pLoc+T-1;
rLoc=hLoc+T-1;
wLoc=rLoc+T-2;
bLoc=wLoc+T-1;
qLoc=bLoc+T-1;
cLoc=qLoc+T-1;
kLoc=cLoc+T-1;
ys=[0;x(yLoc:yLoc+T-2)];
ps=[1;x(pLoc:pLoc+T-2)];  %  Initial value of p is 1
hs=[0;x(hLoc:hLoc+T-2)];
rs=[0;0;x(rLoc:rLoc+T-3)]; %  rs(2) is undefined
ws=[0;x(wLoc:wLoc+T-2)];
bs=[0;x(bLoc:bLoc+T-2)];
qs=[1;x(qLoc:qLoc+T-2)]; % Stationary value of q is 1
cs=[0;x(cLoc:cLoc+T-2)];
ks=[StatVars(9);x(kLoc:kLoc+T-3);StatVars(9)]; % k is k* in periods 1 and T

%  Calculate auxiliary variables

ms=zeros(T,1);
for t=2:T-1;
    mucN=cs(t)^(-1/sigma)*(1-chi*(1-1/sigma)*hs(t)^(1+1/psi));
    mucF=cs(t+1)^(-1/sigma)*(1-chi*(1-1/sigma)*hs(t+1)^(1+1/psi));
    ms(t)=beta*mucF/mucN;
end
% ms(T)=beta;  %  Stationary value of discounter
ptildes=zeros(T,1);
Ds=zeros(T,1);
ptildes(1)=1;
Ds(1)=1;
D=1;
for t=2:T
    ptildes(t)=((ps(t)^(-(epsilon-1))-theta*ps(t-1)^(-(epsilon-1)))/(1-theta))^(-1/(epsilon-1));
    D=(1-theta)*ptildes(t)^(-epsilon)+theta*D;
    Ds(t)=D;
end
zs=zeros(T,1);
mult=alpha^(-alpha)*(1-alpha)^(-(1-alpha));
for t=2:T
    zs(t)=mult*ws(t)^alpha*bs(t)^(1-alpha);
end

%  Backward recursions for A and B

yStar=ys(1);
A=yStar/(1-theta*beta);
B=zs(1)*yStar/(1-theta*beta);
t=T;
As=zeros(T,1);
Bs=zeros(T,1);
for tback=1:T
    A=ys(t)*ps(t)^(epsilon-1)+theta*ms(t)*A;
    As(t)=A;
    B=zs(t)*ys(t)*ps(t)^epsilon+theta*ms(t)*B;
    Bs(t)=B;
    t=t-1; 
end

iDisc=1;

%  Discrepancies except Euler equation and rental price of capital

for t=2:T
    % Labor demand
    discs(iDisc)=hs(t)-alpha*zs(t)*ys(t)*ps(t)^epsilon*Ds(t)/ws(t);
    iDisc=iDisc+1;
    %  Tobin's q
    discs(iDisc)=kappa*(ks(t)-ks(t-1))/ks(t-1) -(qs(t)-1);
    iDisc=iDisc+1;
    % Markup
    discs(iDisc)=ptildes(t)-epsilon*Bs(t)/(As(t)*(epsilon-1));
    iDisc=iDisc+1;
    % Taylor
    if t<T
        discs(iDisc)=rs(t+1)-rBar-phi*ps(t)/ps(t-1)-lambda*ys(t)/yBar;
    else
        discs(iDisc)=ps(T)-ps(T-1);
    end
    iDisc=iDisc+1;
    % Labor supply
    muc=cs(t)^(-1/sigma)*(1-chi*(1-1/sigma)*hs(t)^(1+1/psi));
    muh=(-chi*(1+1/psi)*cs(t)^(1-1/sigma)-gamma)*hs(t)^(1/psi);
    discs(iDisc)=muh+ws(t)*muc;
    iDisc=iDisc+1;
    % Capital demand
    discs(iDisc)=ks(t-1)-(1-alpha)*zs(t)*ys(t)*ps(t)^epsilon*Ds(t)/bs(t);
    iDisc=iDisc+1;
    % Law of motion of capital
    discs(iDisc)=ks(t)+kappa*(ks(t)-ks(t-1))^2/(2*ks(t-1))-(1-delta)*ks(t-1)-ys(t)+cs(t)+gs(t);
    iDisc=iDisc+1;
end

%  Rental price of capital

for t=3:T
     discs(iDisc)=bs(t)-(qs(t-1)*(rs(t)+delta)-(qs(t)-qs(t-1)));
     iDisc=iDisc+1;
end

%  Euler 

for t=2:T-1
    discs(iDisc)=(1+rs(t+1))*ms(t)-1;
    iDisc=iDisc+1;
end

if IfSave==0
    return
end

AllVars=[ys,ps,hs,rs,ws,bs,qs,cs,ks,gs,ms,ptildes,Ds,zs,As,Bs];
save NKACCModel AllVars
