function [U,U1,U2,UM,UL] = util(C,Clag,M,L,shock,param);
%UTIL    computes utility and its derivatives.
%
%        Utype=1:
%          U(C,Clag,M,L) = [CES(C,M) * (1-L)^psi]^(1-sigma)/(1-sigma)
%
%        Utype=2:
%          U(C,Clag,M,L) = CES(C,M)^(1-sigma)/(1-sigma)+psi*(1-L)^(1-xi)/(1-xi)
%
%        Utype=3:
%          U(C,Clag,M,L) = [C^(1-kappa)/(1-kappa) -psi * L^(1+xi)/(1+xi)
%                                    +omega * M^nu/nu]^(1-sigma)/(1-sigma)
%        Utype=4:
%          U(C,Clag,M,L) = [C^(1-kappa)/(1-kappa) +psi * (1-L)^(1-xi)/(1-xi)
%                                    +omega * M^nu/nu]^(1-sigma)/(1-sigma)
%        Utype=5:
%          U(C,Clag,M,L) = shock*log(C-kappa*Clag)-psi * L^(1+xi)/(1+xi)
%                                    + omega * M^nu/nu
%
%        where 
%          CES(C,M) = [omega*C^nu+(1-omega)*M^nu]^(1/nu)
%          nu       = (eta-1)/eta
%

%        Ellen McGrattan, 2-1-00
%        Revised, ERM, 4-27-00
%

Utype    = param(4);
eta      = param(5);
kappa    = param(6);
omega    = param(7);
psi      = param(8);
sigma    = param(9);
xi       = param(10);
nu       = (eta-1)/eta;
omeg1    = 1-omega;
sigm1    = 1-sigma;

x        = omega * C.^nu + omeg1 * M.^nu;

if Utype == 1;
  U        = (x.^(1/nu) .* (1-L).^psi).^sigm1/sigm1;
  U1       = sigm1*omega *U.*C.^(nu-1)./x;
  U2       = 0*U1;
  UM       = sigm1*omeg1 *U.*M.^(nu-1)./x;
  UL       = -sigm1*psi *U./(1-L);

elseif Utype == 2;

  U0       = x.^(sigm1/nu)/sigm1;
  U        = U0+psi* (1-L).^(1-xi)/(1-xi);
  U1       = sigm1*omega *U1.*C.^(nu-1)./x;
  U2       = 0*U1;
  UM       = sigm1*omeg1 *U1.*M.^(nu-1)./x;
  UL       = -psi*(1-L).^(-xi);

elseif Utype == 3;

  Psi      = C.^(1-kappa)/(1-kappa)-psi*L.^(1+xi)/(1+xi)+omega*M.^nu/nu+100;
  U        = Psi.^sigm1/sigm1;
  U1       = Psi.^(-sigma) .* C.^(-kappa);
  U2       = 0*U1;
  UM       = Psi.^(-sigma) .* M.^(nu-1)*omega;
  UL       = Psi.^(-sigma) .* L.^(xi) *(-psi);

elseif Utype == 4;

  Phi      = C.^(1-kappa)/(1-kappa)+psi* (1-L).^(1-xi)/(1-xi)+ ...
             omega*M.^nu/nu+1000;
  U        = Phi.^sigm1/sigm1;
  U1       = Phi.^(-sigma) .* C.^(-kappa);
  U2       = 0*U1;
  UM       = Phi.^(-sigma) .* M.^(nu-1)*omega;
  UL       = Phi.^(-sigma) .* (1-L).^(-xi) *(-psi);

elseif Utype == 5;

  U        = shock.*log(C-kappa*Clag)-psi*L.^(1+xi)/(1+xi)+omega*M.^nu/nu+100;
  U1       = shock./(C-kappa*Clag);
  U2       = -kappa*shock./(C-kappa*Clag);
  UM       = omega*M.^(nu-1);
  UL       = -psi*L.^xi;

else

  error(' WARNING: Utype is not in {1,2,3,4,5}.')

end;
