function y = sscu(x,param)

hrs    = x(1);
emp    = x(2);
c1     = x(3);
k      = x(4);

beta   = param(1);
delta  = param(2);
gamma  = param(3);
lam    = param(4);
omega  = param(5);
pa     = param(6);
pb     = param(7);
phi    = param(8);
psi    = param(9);
rho    = param(10);
sig    = param(11);
theta  = param(12);
ig     = param(13);
bg     = param(14);
army   = param(15);

kg     = ig/delta;
z      = 1; 
ome1   = 1-omega;
thet1  = 1-theta;
gamm1  = 1-gamma;
hwgt   = (1-gamma)*(1-omega);
cwgt   = gamma*(1-omega);
cwgt1  = cwgt-1;
hwgt1  = hwgt-1;

zkk    =  z*lam*(phi*k^rho+(1-phi)*kg^rho)^(theta/rho);
c0     =  c1*(1-hrs)^(hwgt/cwgt1);
%dpn   =  pa+2*pb*emp;
d      =  (emp-pa)*pb;
dpn    =  3.d0*d*d*(1+sign(d))*pb;
y(1,1) =  gamm1/gamma*emp^theta*c1-zkk*(1-hrs);
y(2,1) =  (1-hrs)^hwgt/ome1-(c0/c1)^cwgt/ome1-dpn/(c1^cwgt)- ...
          gamma*(1-hrs)^hwgt/c1*(c1-c0-zkk*thet1*emp^(-theta)*hrs);
y(3,1) =  emp*c1+(1-emp)*c0+delta*k+ig+bg-zkk*hrs*emp^thet1;
y(4,1) =  beta*(1-delta+theta*phi*k^(rho-1)*zkk/(phi*k^rho+ ...
          (1-phi)*kg^rho)*emp^thet1*hrs)-1;


