function [f,df,d2f,d3f,dfp,dfzp,d2fzp]=ces1(x,p);
%CES   computes derivatives for the following function:
%                        b                    b    1/b
%               f = ( a x   + (1-a) (exp(s)*y)   )        
%     
%      [f,df,d2f,d3f,dfp,dfzp,d2fzp] = ces3(z,p) (z=[s,x,y],p=[a,b])
%      computes first (df), second (d2f), and third (d3f) derivatives 
%      of f with respect to z, the derivative of f with respect to p
%      (dfp), the derivative of f with respect to z and p (dfzp), and 
%      the derivative of d2f with respect to p (d2fzp).
%      Column 1 of d3f is equal to vec( d/dx(d2f) ) and column 2
%      of d3f is vec( d/dy(d2f) ).
%      See also CES, CES1, CES2.

%      Ellen McGrattan, 1-9-91
%      Revised, 6-5-91, ERM

s=x(1); y=x(3); x=x(2); a=p(1); b=p(2);

c = a*x^b; 
r = exp(s);
d = (1-a)*r^b*y^b; 
z = c + d;
f = z^(1/b);

if nargout > 1;   
  I     = eye(3);
  t     = [d;c/x;d/y]/z;
  df    = f*t;
  dt    = [b*c*d,-b*c*d/x,b*c*d/y; 
          -b*c*d/x,c/(x*x)*(b*d-z),-b*c*d/(x*y);
           b*c*d/y,-b*c*d/(x*y),d/(y*y)*(b*c-z)]/(z*z);
  d2f   = df*t' + f*dt;
  d2t   = zeros(9,3);
  bcdpz = b*(c-d)+z;
  bcdmz = bcdpz-2*z;
  d2t(1,1) = b*b*c*d*(c-d);
  d2t(1,2) = b*b*c*d*(d-c)/x;
  d2t(1,3) = b*b*c*d*(c-d)/y;
  d2t(2,1) = d2t(1,2);
  d2t(2,2) = b*c*d*bcdpz/(x*x);
  d2t(2,3) = b*b*c*d*(d-c)/(x*y);
  d2t(3,1:2) = d2t(1:2,3)';
  d2t(3,3) = -b*c*d*(z-b*(c-d))/(y*y);
  d2t(4:6,1) = d2t(1:3,2);
  d2t(5:6,1) = d2t(2,2:3)';
  d2t(5,2) = (-bcdpz*(b*d-z)-(b-1)*z*z)*c/x^3;
  d2t(5,3) = bcdpz*b*c*d/(x^2*y);
  d2t(6,2:3)= b*c*d*[bcdpz/x,-bcdmz/y]/(x*y);
  d2t(4,:) = d2t(2,:);
  d2t(7,:) = d2t(3,:);
  d2t(8,:) = d2t(6,:);
  d2t(9,1) = d2t(3,3);
  d2t(9,2) = -bcdmz*b*c*d/(x*y^2);
  d2t(9,3) = (bcdmz*(b*c-z)-(b-1)*z*z)*d/y^3;
  d3f   = kron(t,I)*d2f+kron(I,df)*dt+vec(dt)*df'+f*d2t/z^3;   

  % to calculate dfp(2), use:  d(x(p)^g(p))/dp = x^g*(g'*ln(x)+g*x'/x) 
  t     = [(x^b-(r*y)^b)/z; -log(z)/b+(c*log(x)+d*(s+log(y)))/z]/b; 
  dfp   = f*t;
  dt    = [[-(r*y)^b-b*d*t(1),d*(s+log(y))-d*(c*log(x)+d*(s+log(y)))/z]/z;
           [1/x;-1/y]*[d*x^b+c*(r*y)^b,c*d*log(x/(r*y))]/(z*z)]; 
  d2t   = zeros(9,2);
  d2t(1,1) = b*(d-c)*(d*x^b+c*(r*y)^b)/z^3; 
  d2t(2,1) =-d2t(1,1)/x;
  d2t(3,1) = d2t(1,1)/y;
  d2t(1,2) = c*d*(b*(d-c)*(log(x)-log(r*y))+z)/(z^3);
  d2t(2,2) =-c*d*(b*(d-c)*(log(x)-log(r*y))+z)/(x*z^3);
  d2t(3,2) = c*d*(b*(d-c)*(log(x)-log(r*y))+z)/(y*z^3);
  d2t(5,1) =-dt(2,1)*bcdpz/(x*z); 
  d2t(6,1) = dt(2,1)*b*(c-d)/(y*z); 
  d2t(9,1) = dt(3,1)*bcdmz/(y*z); 
  d2t(5,2) = c*d*(-log(x/(r*y))*bcdpz+z)/(x^2*z^3);
  d2t(6,2) = c*d*( log(x/(r*y))*b*(c-d)-z)/(x*y*z^3);
  d2t(9,2) = c*d*( log(r*y/x)*bcdmz+z)/(y^2*z^3);
  d2t(4,:) = d2t(2,:);
  d2t(7,:) = d2t(3,:);
  d2t(8,:) = d2t(6,:);
  dfzp  = df*t' + f*dt;
  d2fzp = vec(d2f)*t'+(kron(I,df)+kron(df,I))*dt+f*d2t ;
end;

