function [u,uj,uh,uc,du,dj,dh]=hp4au(x,p)
%HP4AU   Utility function for home production economy (version 4a).
%        [u,uj,uh,uc,du,dj,dh]=hp4au(x,p) takes variable vector (x) 
%        and parameter vector (p) as inputs and computes u=u(x,p), 
%        uj=du(x,p)/dx, uh=d2u(x,p)/dx2, uc=grad_x (d2u(x,p)/dx2), 
%        du=du(x,p)/dp, dj=d2u(x,p)/dxdp, dh=grad_p (d2u(x,p)/dx2). 

%        Ellen McGrattan, 5-8-91
%        Revised, 5-22-91, ERM


H=1134;
km = x(1);  hm = x(2);  i  = x(3);  hn = x(4); kn = x(5); 
g  = x(6);  sm = x(7);  sn = x(8);  tk = x(9); tn = x(10);
Km = x(11); Hm = x(12); 
nx = 12;

%parameters
a0=p(1); b0=p(2); a1=p(3); b1=p(4); a2=p(5); b2=p(6); a3=p(7); b3=p(8);   
a4=p(9); b4=p(10); b=p(11); risk=p(12); delta=p(13);
np = 13;
e   = eye(nx);
s   = eye(np);

[Y,dY,d2Y,d3Y,dpY,dxpY,d2pY]        = ces3([sm;Km;Hm],[a4;b4]);
[cn,dcn,d2cn,d3cn,dpcn,dxpcn,d2pcn] = ces3([sn;kn;hn],[a3;b3]);
[ht,dht,d2ht,d3ht,dpht,dxpht,d2pht] = ces([hm;hn],[a0;b0]);
cg  = g*Y;
v   = a4*Km^b4;
w   = (1-a4)*(exp(sm)*Hm)^b4;
z   = v+w;
t1  = v/z;
t2  = w/z;
cm  = (1-tk)*Y*t1*km/Km+(1-tn)*Y*t2*hm/Hm+tk*Y*t1+tn*Y*t2-cg-i+delta*tk*(km-Km);
cp  = (a2*cm^b2 + (1-a2)*cn^b2)^(1/b2);
ct  = (a1*cp^b1 + (1-a1)*cg^b1)^(1/b1);
ht  = (a0*hm^b0 + (1-a0)*hn^b0)^(1/b0);
ht  = H-2*ht; 
cs  = 1000;
hs  = 600;
u   = ( (ct/cs)^b*(ht/hs)^(1-b) )^(1-risk)/(1-risk);
scale = ((1/cs)^b*(1/hs)^(1-b))^(1-risk);

if nargout>1;
  e1    = e([7,11,12],:);
  e2    = s([9,10],:);
  dY    = e1'*dY;
  d2Y   = e1'*d2Y*e1;
  d3Y   = kron(e1',e1')*d3Y*e1;
  dpY   = e2'*dpY;
  dxpY  = e1'*dxpY*e2;
  d2pY  = kron(e1',e1')*d2pY*e2;
  e1    = e([8,5,4],:);
  e2    = s([7,8],:);
  dcn   = e1'*dcn;
  d2cn  = e1'*d2cn*e1;
  d3cn  = kron(e1',e1')*d3cn*e1;
  dpcn  = e2'*dpcn;
  dxpcn = e1'*dxpcn*e2;
  d2pcn = kron(e1',e1')*d2pcn*e2;
  e1    = e([2,4],:);
  e2    = s([1,2],:);
  dht   = -2*e1'*dht;
  d2ht  = -2*e1'*d2ht*e1;
  d3ht  = -2*kron(e1',e1')*d3ht*e1;
  dpht  = -2*e2'*dpht;
  dxpht = -2*e1'*dxpht*e2;
  d2pht = -2*kron(e1',e1')*d2pht*e2;
  % 
  %  Construct cm and its derivatives 
  %
  dcm    = b4*v*w*[-1;1/Km;-1/Hm]/(z*z);
  d2cm   = b4*v*w*[b4*(w-v)  ,b4*(v-w)/Km,b4*(w-v)/Hm; 
                 b4*(v-w)/Km,(b4*(w-v)-z)/Km^2,b4*(v-w)/(Km*Hm);
                 b4*(w-v)/Hm,b4*(v-w)/(Km*Hm),(b4*(w-v)+z)/Hm^2]/z^3;
  d3cm       = zeros(9,3);
  d3cm(1,1)  = b4*b4*b4*v*w*(4*v*w-w*w-v*v)/z^4;
  d3cm(2,1)  =-d3cm(1,1)/Km;
  d3cm(3,1)  = d3cm(1,1)/Hm;
  d3cm(1,:)  = d3cm(1:3,1)';
  d3cm(2,2)  = b4*b4*v*w*(4*b4*v*w-(b4-1)*w*w-(b4+1)*v*v)/(Km*Km*z^4);
  d3cm(2,3)  =-b4*b4*b4*v*w*(4*v*w-w*w-v*v)/(Km*Hm*z^4);
  d3cm(3,3)  = b4*b4*v*w*(4*b4*v*w-(b4+1)*w*w-(b4-1)*v*v)/(Hm*Hm*z^4);
  d3cm(3,2)  = d3cm(2,3);
  d3cm(5,1)  = (d3cm(1,1) - b4*b4*v*w*(v-w)/z^3)/Km^2;
  d3cm(5,2)  = ( (b4+2)*(b4+1)*v*v-4*(b4*b4-1)*v*w+(b4-2)*(b4-1)*w*w )/ ...
                 (Km^3*z^4)*b4*v*w;
  d3cm(6,1)  = d3cm(3,2);
  d3cm(6,2)  = ( -(b4+1)*v*v+4*b4*v*w-(b4-1)*w*w )/(Hm*Km^2*z^4)*b4*b4*v*w;
  d3cm(9,1)  = (d3cm(1,1) + b4*b4*v*w*(v-w)/z^3)/Hm^2;
  d3cm(9,2)  = ( (b4-1)*v*v-4*b4*v*w+(b4+1)*w*w )/(Hm^2*Km*z^4)*b4*b4*v*w;
  d3cm(5:6,3)= d3cm([6,9],2);
  d3cm(9,3)  =-( (b4-2)*(b4-1)*v*v-4*(b4*b4-1)*v*w+(b4+2)*(b4+1)*w*w )/ ...
                 (Hm^3*z^4)*b4*v*w;
  d3cm(4,:)  = d3cm(2,:);
  d3cm(7,:)  = d3cm(3,:);
  d3cm(8,:)  = d3cm(6,:);
  dpcm       = v*w*[1/(a4*(1-a4));log(Km/(exp(sm)*Hm))]/(z*z);
  t          = [(w-v)/(a4*(1-a4)),(log(Km/(exp(sm)*Hm))*b4*(w-v)+z)/b4]/z;
  dxpcm      = dcm*t;
  dt         = -2*b4*v*w*[-1;1/Km;-1/Hm]*[1/(a4*(1-a4)),log(Km/(exp(sm)*Hm))- ...
                .5*(w-v)*z/(b4*v*w)]/(z*z);
  d2pcm      = vec(d2cm)*t+kron(eye(3),dcm)*dt;
  e1         = e([7,11,12],:);
  e2         = s([9,10],:);
  dcm        = e1'*dcm;
  d2cm       = e1'*d2cm*e1;
  d3cm       = kron(e1',e1')*d3cm*e1;
  dpcm       = e2'*dpcm;
  dxpcm      = e1'*dxpcm*e2;
  d2pcm      = kron(e1',e1')*d2pcm*e2;
  [t,dt,d2t,d3t,dpt,dxpt,d2pt]=dmult(t1,dcm,d2cm,d3cm,dpcm, ...
                         dxpcm,d2pcm,Y,dY,d2Y,d3Y,dpY,dxpY,d2pY);  %t=t1*Y
  e1    = e(:,1);
  e2    = e(:,11);
  dcm   = e1/Km-km*e2/(Km*Km);
  d2cm  = -(e1*e2'+e2*e1')/(Km*Km)+2*km*e2*e2'/Km^3;
  d3cm  = 2*vec(e1*e2'+e2*e1')/Km^3*e2'+2*vec(e2*e2')*(e1'-3*km*e2'/Km)/Km^3;
  dpcm  = zeros(np,1);
  dxpcm = zeros(nx,np);
  d2pcm = zeros(nx*nx,np);
  [cm,dcm,d2cm,d3cm,dpcm,dxpcm,d2pcm]=dmult(t,dt,d2t,d3t,dpt,dxpt,d2pt, ...
                         km/Km,dcm,d2cm,d3cm,dpcm,dxpcm,d2pcm); %cm=t1*Y*km/Km
  e1    = e(:,9);
  dcg   = (1-tk)*dcm-cm*e1+tk*dt+t1*Y*e1;
  d2cg  = (1-tk)*d2cm-dcm*e1'-e1*dcm'+tk*d2t +dt*e1'+e1*dt';
  d3cg  = (1-tk)*d3cm-vec(d2cm)*e1'-(kron(e1,e)+kron(e,e1))*d2cm+ ...
          tk*d3t +vec(d2t)*e1'+(kron(e1,e)+kron(e,e1))*d2t;
  dpcg  = (1-tk)*dpcm+tk*dpt;
  dxpcg = (1-tk)*dxpcm-e1*dpcm'+tk*dxpt+e1*dpt';
  d2pcg = (1-tk)*d2pcm-(kron(e1,e)+kron(e,e1))*dxpcm+ ...
          tk*d2pt+(kron(e1,e)+kron(e,e1))*dxpt;   %cm=(1-tk)*t1*Y*km/Km+tk*Y*t1
  %
  %  Add the second term to cm: 
  % 
  dcm    = b4*v*w*[1;-1/Km;1/Hm]/(z*z);
  d2cm   =-b4*v*w*[b4*(w-v)  ,b4*(v-w)/Km,b4*(w-v)/Hm; 
                 b4*(v-w)/Km,(b4*(w-v)-z)/Km^2,b4*(v-w)/(Km*Hm);
                 b4*(w-v)/Hm,b4*(v-w)/(Km*Hm),(b4*(w-v)+z)/Hm^2]/z^3;
  d3cm       = zeros(9,3);
  d3cm(1,1)  = b4*b4*b4*v*w*(4*v*w-w*w-v*v)/z^4;
  d3cm(2,1)  =-d3cm(1,1)/Km;
  d3cm(3,1)  = d3cm(1,1)/Hm;
  d3cm(1,:)  = d3cm(1:3,1)';
  d3cm(2,2)  = b4*b4*v*w*(4*b4*v*w-(b4-1)*w*w-(b4+1)*v*v)/(Km*Km*z^4);
  d3cm(2,3)  =-b4*b4*b4*v*w*(4*v*w-w*w-v*v)/(Km*Hm*z^4);
  d3cm(3,3)  = b4*b4*v*w*(4*b4*v*w-(b4+1)*w*w-(b4-1)*v*v)/(Hm*Hm*z^4);
  d3cm(3,2)  = d3cm(2,3);
  d3cm(5,1)  = (d3cm(1,1) - b4*b4*v*w*(v-w)/z^3)/Km^2;
  d3cm(5,2)  = ( (b4+2)*(b4+1)*v*v-4*(b4*b4-1)*v*w+(b4-2)*(b4-1)*w*w )/ ...
                 (Km^3*z^4)*b4*v*w;
  d3cm(6,1)  = d3cm(3,2);
  d3cm(6,2)  = ( -(b4+1)*v*v+4*b4*v*w-(b4-1)*w*w )/(Hm*Km^2*z^4)*b4*b4*v*w;
  d3cm(9,1)  = (d3cm(1,1) + b4*b4*v*w*(v-w)/z^3)/Hm^2;
  d3cm(9,2)  = ( (b4-1)*v*v-4*b4*v*w+(b4+1)*w*w )/(Hm^2*Km*z^4)*b4*b4*v*w;
  d3cm(5:6,3)= d3cm([6,9],2);
  d3cm(9,3)  =-( (b4-2)*(b4-1)*v*v-4*(b4*b4-1)*v*w+(b4+2)*(b4+1)*w*w )/ ...
                 (Hm^3*z^4)*b4*v*w;
  d3cm(4,:)  = d3cm(2,:);
  d3cm(7,:)  = d3cm(3,:);
  d3cm(8,:)  = d3cm(6,:);
  d3cm       = -d3cm;
  dpcm       =-v*w*[1/(a4*(1-a4));log(Km/(exp(sm)*Hm))]/(z*z);
  t          = [(w-v)/(a4*(1-a4)),(log(Km/(exp(sm)*Hm))*b4*(w-v)+z)/b4]/z;
  dxpcm      = dcm*t;
  dt         = -2*b4*v*w*[-1;1/Km;-1/Hm]*[1/(a4*(1-a4)),log(Km/(exp(sm)*Hm))- ...
                .5*(w-v)*z/(b4*v*w)]/(z*z);
  d2pcm      = vec(d2cm)*t+kron(eye(3),dcm)*dt;
  e1         = e([7,11,12],:);
  e2         = s([9,10],:);
  dcm        = e1'*dcm;
  d2cm       = e1'*d2cm*e1;
  d3cm       = kron(e1',e1')*d3cm*e1;
  dpcm       = e2'*dpcm;
  dxpcm      = e1'*dxpcm*e2;
  d2pcm      = kron(e1',e1')*d2pcm*e2;
  [t,dt,d2t,d3t,dpt,dxpt,d2pt]=dmult(t2,dcm,d2cm,d3cm,dpcm, ...
                         dxpcm,d2pcm,Y,dY,d2Y,d3Y,dpY,dxpY,d2pY);  %t=t2*Y
  e1    = e(:,2);
  e2    = e(:,12);
  dcm   = e1/Hm-hm*e2/(Hm*Hm);
  d2cm  = -(e1*e2'+e2*e1')/(Hm*Hm)+2*hm*e2*e2'/Hm^3;
  d3cm  = 2*vec(e1*e2'+e2*e1')/Hm^3*e2'+2*vec(e2*e2')*(e1'-3*hm*e2'/Hm)/Hm^3;
  dpcm  = zeros(np,1);
  dxpcm = zeros(nx,np);
  d2pcm = zeros(nx*nx,np);
  [cm,dcm,d2cm,d3cm,dpcm,dxpcm,d2pcm]=dmult(t,dt,d2t,d3t,dpt,dxpt,d2pt, ...
                         hm/Hm,dcm,d2cm,d3cm,dpcm,dxpcm,d2pcm); %cm=t2*Y*hm/Hm
  e1    = e(:,10);
  dcg   = dcg+(1-tn)*dcm-cm*e1+tn*dt+t2*Y*e1;
  d2cg  = d2cg+(1-tn)*d2cm-dcm*e1'-e1*dcm'+tn*d2t +dt*e1'+e1*dt';
  d3cg  = d3cg+(1-tn)*d3cm-vec(d2cm)*e1'-(kron(e1,e)+kron(e,e1))*d2cm+ ...
          tn*d3t +vec(d2t)*e1'+(kron(e1,e)+kron(e,e1))*d2t;
  dpcg  = dpcg+(1-tn)*dpcm+tn*dpt;
  dxpcg = dxpcg+(1-tn)*dxpcm-e1*dpcm'+tn*dxpt+e1*dpt';
  d2pcg = d2pcg+(1-tn)*d2pcm-(kron(e1,e)+kron(e,e1))*dxpcm+ ...
          tn*d2pt+(kron(e1,e)+kron(e,e1))*dxpt;  
                     %cm=(1-tk)*t1*Y*km/Km+tk*Y*t1+(1-tn)*t2*Y*hm/Hm+tn*Y*t2
  dcm   = dcg;
  d2cm  = d2cg;
  d3cm  = d3cg;
  dpcm  = dpcg;
  dxpcm = dxpcg;
  d2pcm = d2pcg;
  e1    = e(:,6);
  dcg   = g*dY+Y*e1;
  d2cg  = g*d2Y+dY*e1'+e1*dY';
  d3cg  = g*d3Y+vec(d2Y)*e1'+(kron(e1,e)+kron(e,e1))*d2Y;
  dpcg  = g*dpY;
  dxpcg = g*dxpY+e1*dpY';
  d2pcg = g*d2pY+(kron(e1,e)+kron(e,e1))*dxpY;
  cm    = (1-tk)*Y*t1*km/Km+(1-tn)*Y*t2*hm/Hm+tk*Y*t1+tn*Y*t2-cg-i+ ...
          delta*tk*(km-Km);
  dcm   = dcm  - dcg - e(:,3)+delta*tk*(e(:,1)-e(:,11))+delta*(km-Km)*e(:,9);
  d2cm  = d2cm+delta*(e(:,9)*(e(1,:)-e(11,:))+(e(:,1)-e(:,11))*e(9,:))- d2cg;
  d3cm  = d3cm - d3cg;
  dpcm  = dpcm - dpcg  + tk*(km-Km)*s(:,13);
  dxpcm = dxpcm- dxpcg + (tk*(e(:,1)-e(:,11))+(km-Km)*e(:,9))*s(13,:);
  d2pcm = d2pcm+vec(e(:,9)*(e(1,:)-e(11,:))+(e(:,1)-e(:,11))*e(9,:))* ...
          s(13,:) - d2pcg;
  [cp,dcp,d2cp,d3cp,dpcp,dxpcp,d2pcp]=cesvec(cm,dcm,d2cm,d3cm,dpcm, ....
             dxpcm,d2pcm,cn,dcn,d2cn,d3cn,dpcn,dxpcn,d2pcn,a2,b2,5,6);
  [ct,dct,d2ct,d3ct,dpct,dxpct,d2pct]=cesvec(cp,dcp,d2cp,d3cp,dpcp, ....
             dxpcp,d2pcp,cg,dcg,d2cg,d3cg,dpcg,dxpcg,d2pcg,a1,b1,3,4);
  br   = b*(1-risk);
  b1r  = (1-b)*(1-risk);
  b1r1 = b1r-1;
  br1  = br-1;
  b1r2 = b1r-2;
  br2  = br-2;
  uj   = (br*dct/ct+b1r*dht/ht)*u;
  uh   = br *u*(d2ct+br1*dct*dct'/ct)/ct+ ....
         b1r*u*(d2ht+b1r1*dht*dht'/ht)/ht+ ....
         br*b1r*u*(dct*dht'+dht*dct')/(ct*ht);
  uc   = br*br1*u/(ct*ct)*  ...
         ( vec(dct*dct')*(br2*dct'/ct+b1r*dht'/ht)+ ...
           (kron(dct,eye(nx))+kron(eye(nx),dct))*d2ct  ); 
  uc   = uc+br*b1r*u/(ct*ht)* ...
         ( vec(dct*dht'+dht*dct')*(br1*dct'/ct+b1r1*dht'/ht)+ ... 
           (kron(dct,eye(nx))+kron(eye(nx),dct))*d2ht+ ...
           (kron(dht,eye(nx))+kron(eye(nx),dht))*d2ct );
  uc   = uc+b1r*b1r1*u/(ht*ht)* ...
         ( vec(dht*dht')*(br*dct'/ct+b1r2*dht'/ht)+ ...
           (kron(dht,eye(nx))+kron(eye(nx),dht))*d2ht  );
  uc   = uc+br*u/ct* ....
         ( vec(d2ct)*( br1*dct'/ct+b1r*dht'/ht) + d3ct );
  uc   = uc+b1r*u/ht* ...
         ( vec(d2ht)*( br*dct'/ct+b1r1*dht'/ht) + d3ht);

  du       = (br*dpct/ct+b1r*dpht/ht)*u;
  du(11,1) = u*(1-risk)*log(ct*hs/ht/cs);
  du(12,1) = u/(1-risk)-u*log( (ct/cs)^b*(ht/hs)^(1-b) );

  dj   = br* u*(dxpct+br1*dct*dpct'/ct)/ct+ ....
         b1r*u*(dxpht+b1r1*dht*dpht'/ht)/ht+ ....
         br*b1r*u*(dct*dpht'+dht*dpct')/(ct*ht);
  dj   = dj + (1-risk)*(uj*log(ct*hs/ht/cs)+u*dct/ct-u*dht/ht)* ...
               [zeros(1,10),1,0,0];
  dj   = dj + ( uj/(1-risk)-uj*log((ct/cs)^b*(ht/hs)^(1-b))- ...
                u*(b*dct/ct+(1-b)*dht/ht) )*[zeros(1,11),1,0];
  dh = br*br1*u/(ct*ct) *  ...
       ( vec(dct*dct')*(br2*dpct'/ct+b1r*dpht'/ht)+ ...
         (kron(dct,eye(nx))+kron(eye(nx),dct))*dxpct );
  dh = dh+br*b1r*u/(ct*ht)*  ...
       ( vec(dct*dht'+dht*dct')*(br1*dpct'/ct+b1r1*dpht'/ht)+ ...
         (kron(dct,eye(nx))+kron(eye(nx),dct))*dxpht+ ...
         (kron(dht,eye(nx))+kron(eye(nx),dht))*dxpct );
  dh = dh+b1r*b1r1*u/(ht*ht)* ....
       ( vec(dht*dht')*(br*dpct'/ct+b1r2*dpht'/ht)+ ...
         (kron(dht,eye(nx))+kron(eye(nx),dht))*dxpht );
  dh = dh+br*u/ct* ...
       ( vec(d2ct)*( br1*dpct'/ct+b1r*dpht'/ht)+d2pct);
  dh = dh+b1r*u/ht* ...
       ( vec(d2ht)*( br*dpct'/ct+b1r1*dpht'/ht)+d2pht);
  dh = dh + (1-risk)*vec(uh*log(ct*hs/ht/cs)+uj*dct'/ct-uj*dht'/ht+ ...
                         dct*uj'/ct-dht*uj'/ht+u*d2ct/ct-u*d2ht/ht- ...
                         u*dct*dct'/ct/ct+u*dht*dht'/ht/ht )* ...
                         [zeros(1,10),1,0,0];
  dh = dh + vec( uh/(1-risk)-uh*log((ct/cs)^b*(ht/hs)^(1-b))- ...
                 uj*(b*dct'/ct+(1-b)*dht'/ht)- ...
                 (b*dct/ct+(1-b)*dht/ht)*uj'- ...
                 u*(b*d2ct/ct-b*dct*dct'/ct/ct+ ...
                    (1-b)*d2ht/ht-(1-b)*dht*dht'/ht/ht) )* ...
                 [zeros(1,11),1,0];
end;
