function discs=NKStatModel(x)
global sigma psi chi alpha gamma beta delta epsilon

c=real(x(1));
h=real(x(2));
q=real(x(3));
mu=epsilon/(epsilon-1);
k=q-c;
y=h^alpha*k^(1-alpha);
muc=c^(-1/sigma)*(1-chi*(1-1/sigma)*h^(1+1/psi));
muh=(-chi*(1+1/psi)*c^(1-1/sigma)-gamma)*h^(1/psi);
w=alpha*h^(alpha-1)*k^(1-alpha)/mu;
r=(1-alpha)*h^alpha*k^(-alpha)/mu-delta;
discs(1)=muh+w*muc;
discs(2)=beta*(1+r)-1;
discs(3)=q-(1-delta)*k-y;
