function y = sscumcb(x,param)

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

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

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*k^theta;
c0     = c1*(1-hrs)^(hwgt/cwgt1);
if itype==1;
  pn   = pa*emp+pb*emp*emp+pc;
  dpn  = pa+2*pb*emp;
end;
if itype==2;
  d    = (emp-pa)*pb;
  pn   = d*d*d+abs(d*d*d);
  dpn  = 3.d0*d*d*(1+sign(d))*pb;
end;
if itype==3;
  pn   = (exp(pa*emp)-1)*pb+pc;
  dpn  = pa*pb*exp(pa*emp);
end;
y(1,1) =  gamm1/gamma*emp^theta*c1-zkk*(1-hrs);
y(2,1) =  c1^cwgt*(1-hrs)^hwgt/ome1-c0^cwgt/ome1- ...
          gamma*c1^cwgt1*(1-hrs)^hwgt*(c1-c0-zkk*thet1*emp^(-theta)*hrs-dpn);
y(3,1) =  emp*c1+(1-emp)*c0+delta*k+bg-zkk*hrs*emp^thet1+pn;
y(4,1) =  beta*(1-delta+theta*zkk/k*emp^thet1*hrs)-1;
