function [U,UC,UM,UL,UCC,UMM,ULL,UCM,UCL,ULM] = util(C,M,L,param);
%UTIL    computes utility and its derivatives.
%
%        Utype=1:
%          U(C,M,L) = [CES(C,M) * (1-L)^psi]^(1-sigma)/(1-sigma)
%
%        Utype=2:
%          U(C,M,L) = CES(C,M)^(1-sigma)/(1-sigma) +psi *(1-L)^(1-xi)/(1-xi)
%
%        Utype=3:
%          U(C,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,M,L) = [C^(1-kappa)/(1-kappa) +psi * (1-L)^(1-xi)/(1-xi)
%                                    +omega * M^nu/nu]^(1-sigma)/(1-sigma)
%
%        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);
beta     = param(11);
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;
  UC       = sigm1*omega *U.*C.^(nu-1)./x;
  UM       = sigm1*omeg1 *U.*M.^(nu-1)./x;
  UL       = -sigm1*psi *U./(1-L);
  UCC      = UC.*UC./U* (1-nu/sigm1) + (nu-1)*UC./C;
  UMM      = UM.*UM./U* (1-nu/sigm1) + (nu-1)*UM./M;
  UCM      = UC.*UM./U* (1-nu/sigm1);
  UCL      = UC.*UL./U;
  ULM      = UL.*UM./U;
  ULL      = UL.*UL./U + UL./(1-L);

elseif Utype == 2;

  U1       = x.^(sigm1/nu)/sigm1;
  U        = U1+psi* (1-L).^(1-xi)/(1-xi);
  UC       = sigm1*omega *U1.*C.^(nu-1)./x;
  UM       = sigm1*omeg1 *U1.*M.^(nu-1)./x;
  UL       = -psi*(1-L).^(-xi);
  UCC      = UC.*UC./U1* (1-nu/sigm1) + (nu-1)*UC./C;
  UMM      = UM.*UM./U1* (1-nu/sigm1) + (nu-1)*UM./M;
  UCM      = UC.*UM./U1* (1-nu/sigm1);
  UCL      = zeros(length(C),1);
  ULM      = UCL;
  ULL      = xi*UL./(1-L);

elseif Utype == 3;

  o        = ones(length(C),1);
  Psi      = C.^(1-kappa)/(1-kappa)-psi*L.^(1+xi)/(1+xi)+omega*M.^nu/nu+100;
  U        = Psi.^sigm1/sigm1;
  UC       = Psi.^(-sigma) .* C.^(-kappa);
  UM       = Psi.^(-sigma) .* M.^(nu-1)*omega;
  UL       = Psi.^(-sigma) .* L.^(xi) *(-psi);
  UCC      = -UC.*(sigma*(C.^(-kappa))./Psi+kappa*o./C);
  UMM      = -UM.*(sigma*(M.^(nu-1))./Psi*omega+(1-nu)*o./M);
  ULL      =  UL.*(sigma*(L.^(xi))./Psi*psi+xi*o./L);
  UCL      =  UC.*(L.^(xi))./Psi*psi*sigma;
  ULM      =  UM.*(L.^(xi))./Psi*psi*sigma;
  UCM      = -UC.*(M.^(nu-1))./Psi*omega*sigma;

elseif Utype == 4;

  o        = ones(length(C),1);
  Phi      = C.^(1-kappa)/(1-kappa)+psi* (1-L).^(1-xi)/(1-xi)+ ...
             omega*M.^nu/nu+1000;
  U        = Phi.^sigm1/sigm1;
  UC       = Phi.^(-sigma) .* C.^(-kappa);
  UM       = Phi.^(-sigma) .* M.^(nu-1)*omega;
  UL       = Phi.^(-sigma) .* (1-L).^(-xi) *(-psi);
  UCC      = -UC.*(sigma*(C.^(-kappa))./Phi+kappa*o./C);
  UMM      = -UM.*(sigma*(M.^(nu-1))./Phi*omega+(1-nu)*o./M);
  ULL      =  UL.*(sigma*((1-L).^(-xi))./Phi*psi+xi*o./(1-L));
  UCL      =  UC.*((1-L).^(-xi))./Phi*psi*sigma;
  ULM      =  UM.*((1-L).^(-xi))./Phi*psi*sigma;
  UCM      = -UC.*(M.^(nu-1))./Phi*omega*sigma;

else

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

end;
