function [L,Lt,dL,dLt]=hpse(param,ind);
%HPSE      Standard Errors for Home Production Model.

%          Ellen R. McGrattan,  6-12-91
%          Revised, ERM, 1-15-95

global zxvar
TH         = 1134;
a0         = param(1);
b0         = param(2);
a1         = param(3);
b1         = param(4);
a2         = param(5);
b2         = param(6);
a3         = param(7);
b3         = param(8);
a4         = param(9);
b4         = param(10);
b          = param(11);
risk       = param(12);
delta      = param(13);
beta       = param(14);
grate      = param(15);
c1         = zeros(5);
c2         = zeros(5);
c3         = zeros(5);
c4         = zeros(5);
c0         = param(16:20);
c1(:)      = param(21:45);
c2(:)      = param(46:70);
c3(:)      = param(71:95);
c4(:)      = param(96:120);
varm       = param(121:128);
coefm      = param(129:136);
RR         = diag(varm);
D          = diag(coefm);
if isempty(ind); 
  ind = ones(1,136); 
end;
Theta      = param(find(ind));
np         = sum(ind);
fin        = find(ind);
si         = cumsum(ind);
i          = find(si==1);
si(1:i-1)  = ones(1,i-1);
%---------------------------------------------------------------------

% parameter vector:
param= [a0,b0,a1,b1,a2,b2,a3,b3,a4,b4,b,risk,delta,beta,grate, ...
        c0(:)',c1(:)',c2(:)',c3(:)',c4(:)',varm',coefm']';
%---------------------------------------------------------------------
% 4. Set gratein = i   if the ith element of Theta corresponds to grate
%                = 0   if grate is not to be estimated     
          
if ind(15); gratein = si(15); else; gratein=0; end;
%---------------------------------------------------------------------

X     = zeros(14,1);
X(2)  = 1;
tm1   = c1+c2;
X(3:7)= (eye(5)-tm1)\c0;
s=[sprintf('means of exog states = %g,',X(3)),sprintf('%g,',X(4)), ...
   sprintf('%g,',X(5)),sprintf('%g,',X(6)),sprintf('%g',X(7))];
%disp(s);

%---------------------------------------------------------------------
% 5. Put bounds on Theta if there are values of the parameters for 
%    which L cannot be computed.

imax  = 0.998;
lb    = [ .05; -10; 0; -10;  0; -10; 0; -10; 0; -10;     0;  ...
          -10; .001;     0; 0.999; -100;-100;-100;0;0; -20000*ones(100,1); ...
           zeros(8,1);     -.99*ones(8,1)];
ub    = [  1;    1; 1;   1;  1;   1; 1;   1; 1;   1; 0.999;  ...
           7; .040; 0.999;  1.01; 100000*ones(105,1);  ...
          100000*ones(8,1); 0.8*ones(8,1)];

if any(param(129:136)>1) | any(abs(eig(tm1))>1) | any(X([3,4,6,7])<0.05);
  disp('Parameter vector such that likelihood can''t be computed');
  L  = 1e+15;
  dL = 1e+15*ones(size(Theta));
  return
end;
%---------------------------------------------------------------------


%---------------------------------------------------------------------
% 6. Specify Ay,Az,B,Cy,Sx,Su,Thet,Psi,ct.
                                         
Ay         = [(1-delta)/grate,zeros(1,11);
              0,1,zeros(1,10);
              zeros(5,1),  c0,c1,c2;
              zeros(5,2),  eye(5), zeros(5)];
Az         =  zeros(12,2);
By         = [0 0 1/grate 0;zeros(11,4)];
C          = zeros(12);
C(3:7,3:7) = c4;
Sx         = zeros(12,14);
Su         = zeros(12,4);
Sx(5:12,[1,3:7,13:14]) = eye(8);
Su(1:4,:)  = eye(4);
Su(5,1)    = -1;
ct         = 2;
Thet       = zeros(2,12);
Psi        = [eye(2),zeros(2)];
%---------------------------------------------------------------------


%---------------------------------------------------------------------
% 7. Specify steady states of X, u, x.                      
%    states:    k,1,g,sm,sn,tk,tn,g(-1),sm(-1),sn(-1),tk(-1),tn(-1),Km,Hm
%    controls:  km,hm,i,hn
%    x:         km,hm,i,hn,kn=k-km,g,sm,sn,tk,tn,Km,Hm

X(8:12)  = X(3:7);
khm      = ( (-a4 + ( (grate-beta*(1-delta)-beta*delta*X(6))/ ...
           (beta*a4*(1-X(6))) )^(b4/(1-b4)) )/(exp(X(4))^b4*(1-a4)) )^(-1/b4);
z0       = [300;140;20];
z        = newton('sshp4b',z0,[param(1:15);X(3:7);khm],1e-7,100);
hm       = z(1); 
hn       = z(2); 
khn      = z(3);
u(2)     = hm;
u(1)     = khm*u(2);
X(1)     = u(1)+khn*hn;
u(3)     = (grate-1+delta)*X(1);
u(4)     = hn;
u        = u(:);
X(13:14) = u(1:2);
x        = Sx*X+Su*u;
nx       = length(x);
ny       = 12;
nz       = 2;
nX       = 14;
nu       = 4;
%---------------------------------------------------------------------


%---------------------------------------------------------------------
% 8. Specify ufn (or replace this string in eval statement in (9) ).

ufn='hp7u';
%---------------------------------------------------------------------

% Compute Q,R,W and necessary matrices for gradQ, gradR, gradW. 
% retrieve vx  = dr(x,theta)/x 
%          vxx = d^2r(x,theta)/dx^2
% plus     vxxx= d^3r(x,theta)/dx^3
%          dv  = dr(x,theta)/dtheta
%          dvx = dvx/dtheta 
%          dvxx= dvxx/dtheta   if derivatives required
if nargout >2;
  eval(['[v,vx,vxx,vxxx,dv,dvx,dvxx]=',ufn,'(x,param(1:13));'])
  %-------------------------------------------------------------------
  tem1 = np-sum(ind(1:13));
  dv   = [dv(find(ind(1:13)));zeros(tem1,1)];
  dvx  = [dvx(:,find(ind(1:13))),zeros(nx,tem1)];
  dvxx = [dvxx(:,find(ind(1:13))),zeros(nx*nx,tem1)];
  %-------------------------------------------------------------------
else;
  eval(['[v,vx,vxx]=',ufn,'(x,param(1:13));'])
end;
vx=vx(:);

x        = Sx*X+Su*u;
e        = zeros(nX,1);
e(ct)    = 1;
Q        = e*(v-vx'*x+.5*x'*vxx*x)*e'+.5*e*(vx'-x'*vxx)*Sx+ ...
           .5*Sx'*(vx-vxx*x)*e'+.5*Sx'*vxx*Sx;
R        = .5*Su'*vxx*Su;
W        = .5*e*(vx'-x'*vxx)*Su+.5*Sx'*vxx*Su;
Qy       = Q(1:ny,1:ny);
Qz       = Q(1:ny,ny+1:nX);
Wy       = W(1:ny,:);
Wz       = W(ny+1:nX,:);

% Compute the equilibrium decision function F

Ri       = inv(R);
BRi      = By*inv(R);
sbeta    = sqrt(beta);
ay       = sbeta*(Ay-BRi*Wy');
az       = sbeta*(Az-BRi*Wz');
by       = sbeta*By;
q        = Q - W*Ri*W';
qy       = q(1:ny,1:ny);
qz       = q(1:ny,ny+1:nX);
Iny      = eye(ny);
Inz      = eye(nz);
psi      = inv(Inz+Psi*Ri*Wz');
thet     = psi*(Thet-Psi*Ri*Wy');
psi      = psi*Psi;
ahat     = ay+az*thet;
bhat     = by+az*psi;
qhat     = qy+qz*thet;
atil     = ay-by*Ri*psi'*qz';
brb      = bhat*Ri*by';
[v,d]    = eig([ahat,zeros(ny);-qhat,Iny],[Iny,brb;zeros(ny),atil']);
d        = diag(d);
[d,in]   = sort(abs(d));
if (~( (d(ny) < 1) & (d(ny+1) > 1) ) ) error('Can''t order eigenvalues'); end;
W11      = v(1:ny,in(1:ny));
W21      = v(ny+1:2*ny,in(1:ny));
P        = real(W21/W11);
Rbpi     = inv(R+by'*P*bhat);
P        = qhat+atil'*P*ahat-atil'*P*bhat*Rbpi*by'*P*ahat;
f        = Rbpi*(by'*P*ahat);
F        = (R+Wz'*Psi)\(R*f+Wy'+Wz'*Thet);

% Compute the matrix for state transition function
Ao       = Ay+Az*Thet-Az*Psi*F-By*F;

%---------------------------------------------------------------------
% 9. Specify G.
%    G =[ K=Km+Kn; Hm; I; Y; g*Y; Kn; tk; tn];

yk  = a4*X(13)^b4;
yh  = (1-a4)*(exp(X(4))*X(14))^b4;
y   = (yk+yh)^(1/b4);
G   = [1,zeros(1,11);
       -F(2:3,:);
       [0,-y*yh*X(4)/(yk+yh),0,y*yh/(yk+yh),zeros(1,8)]- ...
          y*yk/(yk+yh)/X(13)*F(1,:)-y*yh/(yk+yh)/X(14)*F(2,:); 
       [0,-y*yh*X(4)*X(3)/(yk+yh)-X(3)*y,y,y*yh*X(3)/(yk+yh),zeros(1,8)]- ...
          y*yk*X(3)/(yk+yh)/X(13)*F(1,:)-y*yh*X(3)/(yk+yh)/X(14)*F(2,:);
       [1,zeros(1,11)]+F(1,:);
       0,0,0,0,0,1,zeros(1,6);
       0,0,0,0,0,0,1,zeros(1,5)];
Gp  = G';
%---------------------------------------------------------------------
 

%---------------------------------------------------------------------
% 10. Observations.
%     data = Y = [k ,hm, i, cm+i+g, g, kn, tk, tn]';
%              = [knnr+kncd+knr 
%                 lhours
%                 if+cd
%                 cn+cs+cd+if+g-csh
%                 g
%                 kncd+knr
%                 taxk
%                 taxn]
%
%                                
no  = 8;
T   = length(zxvar);
tem = grate.^([1:T]');
Y   = [(zxvar(:,1)+zxvar(:,7))./tem,zxvar(:,2),zxvar(:,3)./tem, ...
        zxvar(:,5)./tem,zxvar(:,6)./tem,zxvar(:,7)./tem,zxvar(:,8),zxvar(:,9)];
DY  = Y(2:T,:)-Y(1:T-1,:)*D';
T   = T-1;
%---------------------------------------------------------------------


%---------------------------------------------------------------------
% 11. Specify initial state and initial date, X0, t0, T0.

X0  = X(1:ny);     
t0  = 1;
T0  = T;  % choose >=2
%---------------------------------------------------------------------

nct          = [1:ct-1,ct+1:ny];
Sigma        = C*C';
nobs         = T-t0+1;
Gbar         = G*Ao-D*G;
Rbar         = RR+G*Sigma*G';
[K,S]        = kfilter(Ao,Gbar,Sigma,Rbar,Sigma*G');
S0           = zeros(ny);     % alternative initial variance-covariances
S0           = S;
%S0(nct,nct) = dlyap(Ao(nct,nct),Sigma(nct,nct));

S1           = S0;
Xt           = zeros(T,ny);
innov        = zeros(T,no);
Xt(t0,:)     = X0';
innov(t0,:)  = DY(t0,:)-X0'*Gbar';
Omega        = Rbar+Gbar*S1*Gbar';
Omegai       = inv(Omega);
Lt(t0,1)     = log(det(Omega))+innov(t0,:)*Omegai*innov(t0,:)';
for i=t0+1:T;
  K1         = (Ao*S1*Gbar'+Sigma*G')*Omegai;
  Xt(i,:)    = Xt(i-1,:)*Ao'+innov(i-1,:)*K1';
  innov(i,:) = DY(i,:)-Xt(i,:)*Gbar';
  S1         = Ao*S1*Ao'+Sigma-(Ao*S1*Gbar'+Sigma*G')*Omegai* ...
               (Gbar*S1*Ao'+G*Sigma);
  Omega      = Rbar+Gbar*S1*Gbar';
  Omegai     = inv(Omega);
  Lt(i,1)    = log(det(Omega))+innov(i,:)*Omegai*innov(i,:)';
end;

%-------------------------------------------------------------------
% 12. Add Jacobian to L if gratein = 1

if gratein;
  Lt   = Lt + 5*log(grate)*(nobs+1);
end;
%-------------------------------------------------------------------
L            = sum(Lt);

if nargout >2;
  %-------------------------------------------------------------------
  % 13. Specify gradAy, gradAz, gradBy, gradC, gradRR, gradD, 
  %     gradx, gradX.

  % user-defined
  gradAy  = zeros(ny*ny,np);
  gradAz  = zeros(ny*nz,np);
  gradBy  = zeros(ny*nu,np);
  gradGp  = zeros(ny*no,np);
  gradC   = zeros(ny*ny,np);
  gradD   = zeros(no*no,np);
  gradRR  = zeros(no*no,np);

  % derivatives with respect to delta
  if ind(13);
    i                = si(13);
    gradAy(1,i)      = -1/grate;
  end;

  % derivatives with respect to grate
  if ind(15);
    i                = si(15);
    gradAy(1,i)      = -(1-delta)/(grate*grate);
    gradBy(25,i)     = -1/(grate*grate);
  end;

  % derivatives with respect to c0
  gradAy(15,si(16))  = gradAy(15,si(16))  + ind(16);
  gradAy(16,si(17))  = gradAy(16,si(17))  + ind(17);
  gradAy(17,si(18))  = gradAy(17,si(18))  + ind(18);
  gradAy(18,si(19))  = gradAy(18,si(19))  + ind(19);
  gradAy(19,si(20))  = gradAy(19,si(20))  + ind(20);

  % derivatives with respect to c1
  gradAy(27,si(21))  = gradAy(27,si(21))  + ind(21);
  gradAy(28,si(22))  = gradAy(28,si(22))  + ind(22);
  gradAy(29,si(23))  = gradAy(29,si(23))  + ind(23);
  gradAy(30,si(24))  = gradAy(30,si(24))  + ind(24);
  gradAy(31,si(25))  = gradAy(31,si(25))  + ind(25);
  gradAy(39,si(26))  = gradAy(39,si(26))  + ind(26);
  gradAy(40,si(27))  = gradAy(40,si(27))  + ind(27);
  gradAy(41,si(28))  = gradAy(41,si(28))  + ind(28);
  gradAy(42,si(29))  = gradAy(42,si(29))  + ind(29);
  gradAy(43,si(30))  = gradAy(43,si(30))  + ind(30);
  gradAy(51,si(31))  = gradAy(51,si(31))  + ind(31);
  gradAy(52,si(32))  = gradAy(52,si(32))  + ind(32);
  gradAy(53,si(33))  = gradAy(53,si(33))  + ind(33);
  gradAy(54,si(34))  = gradAy(54,si(34))  + ind(34);
  gradAy(55,si(35))  = gradAy(55,si(35))  + ind(35);
  gradAy(63,si(36))  = gradAy(63,si(36))  + ind(36);
  gradAy(64,si(37))  = gradAy(64,si(37))  + ind(37);
  gradAy(65,si(38))  = gradAy(65,si(38))  + ind(38);
  gradAy(66,si(39))  = gradAy(66,si(39))  + ind(39);
  gradAy(67,si(40))  = gradAy(67,si(40))  + ind(40);
  gradAy(75,si(41))  = gradAy(75,si(41))  + ind(41);
  gradAy(76,si(42))  = gradAy(76,si(42))  + ind(42);
  gradAy(77,si(43))  = gradAy(77,si(43))  + ind(43);
  gradAy(78,si(44))  = gradAy(78,si(44))  + ind(44);
  gradAy(79,si(45))  = gradAy(79,si(45))  + ind(45);

  % derivatives with respect to c2 
  gradAy(87,si(46))  = gradAy(87,si(46))  + ind(46);
  gradAy(88,si(47))  = gradAy(88,si(47))  + ind(47);
  gradAy(89,si(48))  = gradAy(89,si(48))  + ind(48);
  gradAy(90,si(49))  = gradAy(90,si(49))  + ind(49);
  gradAy(91,si(50))  = gradAy(91,si(50))  + ind(50);
  gradAy(99,si(51))  = gradAy(99,si(51))  + ind(51);
  gradAy(100,si(52)) = gradAy(100,si(52)) + ind(52);
  gradAy(101,si(53)) = gradAy(101,si(53)) + ind(53);
  gradAy(102,si(54)) = gradAy(102,si(54)) + ind(54);
  gradAy(103,si(55)) = gradAy(103,si(55)) + ind(55);
  gradAy(111,si(56)) = gradAy(111,si(56)) + ind(56);
  gradAy(112,si(57)) = gradAy(112,si(57)) + ind(57);
  gradAy(113,si(58)) = gradAy(113,si(58)) + ind(58);
  gradAy(114,si(59)) = gradAy(114,si(59)) + ind(59);
  gradAy(115,si(60)) = gradAy(115,si(60)) + ind(60);
  gradAy(123,si(61)) = gradAy(123,si(61)) + ind(61);
  gradAy(124,si(62)) = gradAy(124,si(62)) + ind(62);
  gradAy(125,si(63)) = gradAy(125,si(63)) + ind(63);
  gradAy(126,si(64)) = gradAy(126,si(64)) + ind(64);
  gradAy(127,si(65)) = gradAy(127,si(65)) + ind(65);
  gradAy(135,si(66)) = gradAy(135,si(66)) + ind(66);
  gradAy(136,si(67)) = gradAy(136,si(67)) + ind(67);
  gradAy(137,si(68)) = gradAy(137,si(68)) + ind(68);
  gradAy(138,si(69)) = gradAy(138,si(69)) + ind(69);
  gradAy(139,si(70)) = gradAy(139,si(70)) + ind(70);

  % derivatives with respect to c3 (need to fill in if 3rd lag is added)
 
  % derivatives with respect to c4
  gradC(27,si(96))   = gradC(27,si(96))   + ind(96);
  gradC(28,si(97))   = gradC(28,si(97))   + ind(97);
  gradC(29,si(98))   = gradC(29,si(98))   + ind(98);
  gradC(30,si(99))   = gradC(30,si(99))   + ind(99);
  gradC(31,si(100))  = gradC(31,si(100))  + ind(100);
  gradC(39,si(101))  = gradC(39,si(101))  + ind(101);
  gradC(40,si(102))  = gradC(40,si(102))  + ind(102);
  gradC(41,si(103))  = gradC(41,si(103))  + ind(103);
  gradC(42,si(104))  = gradC(42,si(104))  + ind(104);
  gradC(43,si(105))  = gradC(43,si(105))  + ind(105);
  gradC(51,si(106))  = gradC(51,si(106))  + ind(106);
  gradC(52,si(107))  = gradC(52,si(107))  + ind(107);
  gradC(53,si(108))  = gradC(53,si(108))  + ind(108);
  gradC(54,si(109))  = gradC(54,si(109))  + ind(109);
  gradC(55,si(110))  = gradC(55,si(110))  + ind(110);
  gradC(63,si(111))  = gradC(63,si(111))  + ind(111);
  gradC(64,si(112))  = gradC(64,si(112))  + ind(112);
  gradC(65,si(113))  = gradC(65,si(113))  + ind(113);
  gradC(66,si(114))  = gradC(66,si(114))  + ind(114);
  gradC(67,si(115))  = gradC(67,si(115))  + ind(115);
  gradC(75,si(116))  = gradC(75,si(116))  + ind(116);
  gradC(76,si(117))  = gradC(76,si(117))  + ind(117);
  gradC(77,si(118))  = gradC(77,si(118))  + ind(118);
  gradC(78,si(119))  = gradC(78,si(119))  + ind(119);
  gradC(79,si(120))  = gradC(79,si(120))  + ind(120);

  % derivatives with respect to varm
  gradRR(1,si(121))  = gradRR(1,si(121))  + ind(121);
  gradRR(10,si(122)) = gradRR(10,si(122)) + ind(122);
  gradRR(19,si(123)) = gradRR(19,si(123)) + ind(123);
  gradRR(28,si(124)) = gradRR(28,si(124)) + ind(124);
  gradRR(37,si(125)) = gradRR(37,si(125)) + ind(125);
  gradRR(46,si(126)) = gradRR(46,si(126)) + ind(126);
  gradRR(55,si(127)) = gradRR(55,si(127)) + ind(127);
  gradRR(64,si(128)) = gradRR(64,si(128)) + ind(128);


  % derivatives with respect to coefm
  gradD(1,si(129))   = gradD(1,si(129))   + ind(129);
  gradD(10,si(130))  = gradD(10,si(130))  + ind(130);
  gradD(19,si(131))  = gradD(19,si(131))  + ind(131);
  gradD(28,si(132))  = gradD(28,si(132))  + ind(132);
  gradD(37,si(133))  = gradD(37,si(133))  + ind(133);
  gradD(46,si(134))  = gradD(46,si(134))  + ind(134);
  gradD(55,si(135))  = gradD(55,si(135))  + ind(135);
  gradD(64,si(136))  = gradD(64,si(136))  + ind(136);

  % derivatives for Thet, Psi
  gradTh = zeros(24,np);
  gradPs = zeros(8,np); 

  % gradx, gradX
  tm4 = fin(find(fin<96));
  del = diag(max(abs(param(1:95))*1e-5,1e-8));
  for i=1:length(tm4);
    gradx(:,i) = (grad7(param(1:95)+del(:,tm4(i))) -  ...
                  grad7(param(1:95)-del(:,tm4(i))))/(2*del(tm4(i),tm4(i)));
  end; 
  gradx = [gradx,zeros(12,np-length(tm4))];
  gradX([1,3:12],:) = [gradx(1,:)+gradx(5,:); ...
                       gradx(6:10,:);gradx(6:10,:)];
  gradX([13:14],:)  = [gradx(1:2,:)];
  %-------------------------------------------------------------------

  tem1    = zeros(ny,nu);
  tem1(:) = [1:ny*nu]';
  tem2    = tem1'; 
  tem2    = tem2(:);
  tem1    = tem1(:);
  Syu     = zeros(ny,nu);
  for i=1:ny*nu;
    Syu(tem1(i),tem2(i))=1;
  end;
  tem1    = zeros(nu,ny);
  tem1(:) = [1:ny*nu]';
  tem2    = tem1'; 
  tem2    = tem2(:);
  tem1    = tem1(:);
  Suy     = zeros(nu,ny);
  for i=1:ny*nu;
    Suy(tem1(i),tem2(i))=1;
  end;
  tem1    = zeros(ny,nz);
  tem1(:) = [1:ny*nz]';
  tem2    = tem1'; 
  tem2    = tem2(:);
  tem1    = tem1(:);
  Syz     = zeros(ny,nz);
  for i=1:ny*nz;
    Syz(tem1(i),tem2(i))=1;
  end;
  tem1    = zeros(nz,nu);
  tem1(:) = [1:nz*nu]';
  tem2    = tem1'; 
  tem2    = tem2(:);
  tem1    = tem1(:);
  Szu     = zeros(nz,nu);
  for i=1:nz*nu;
    Szu(tem1(i),tem2(i))=1;
  end;
  tem1    = zeros(ny);
  tem1(:) = [1:ny*ny]';
  tem2    = tem1'; 
  tem2    = tem2(:);
  tem1    = tem1(:);
  Syy     = zeros(ny);
  for i=1:ny*ny;
    Syy(tem1(i),tem2(i))=1;
  end;
  tem1    = zeros(no);
  tem1(:) = [1:no*no]';
  tem2    = tem1'; 
  tem2    = tem2(:);
  tem1    = tem1(:);
  Soo     = zeros(no);
  for i=1:no*no;
    Soo(tem1(i),tem2(i))=1;
  end;
  tem1    = zeros(ny,no);
  tem1(:) = [1:ny*no]';
  tem2    = tem1'; 
  tem2    = tem2(:);
  tem1    = tem1(:);
  Syo     = zeros(ny,no);
  for i=1:ny*no;
    Syo(tem1(i),tem2(i))=1;
  end;
  gradvx   = vxx*gradx+dvx;
  gradvxx  = vxxx*gradx+dvxx;
  e        = e(1:ny);
  sx       = Sx(:,1:ny);
  sz       = Sx(:,ny+1:nX);

  gradQy   = kron(e,e)*dv(:)'-kron(e*x',e)*gradvx-kron(e,e*vx')*gradx+ ...
             .5*kron(e*x',e*x')*gradvxx+ ...
             .5*(kron(e*x'*vxx,e)+kron(e,e*x'*vxx))*gradx+ ...
             .5*(kron(sx',e)+kron(e,sx'))*gradvx- ...
             .5*(kron(sx',e*x')+kron(e*x',sx'))*gradvxx- ...
             .5*(kron(sx'*vxx,e)+kron(e,sx'*vxx))*gradx+ ...
             .5*kron(sx',sx')*gradvxx;
  gradQz   = .5*kron(sz',sx')*gradvxx+.5*kron(sz',e)*gradvx- ...
             .5*kron(sz'*vxx,e)*gradx-.5*kron(sz',e*x')*gradvxx;
  gradR    = .5*kron(Su',Su')*gradvxx;
  gradWy   = .5*kron(Su',sx')*gradvxx+.5*kron(Su',e)*gradvx- ...
             .5*kron(Su',e*x')*gradvxx-.5*kron(Su'*vxx,e)*gradx;
  gradWz   = .5*kron(Su',Sx(:,ny+1:nX)')*gradvxx;
  Inu      = eye(nu);
  PRWi     = -inv(Inz+Psi/R*Wz');
  RWPi     = inv(R+Wz'*Psi);
  ftil     = (R+bhat'*P'*by)\bhat'*P'*atil;
  aotil    = atil-by*ftil;
  kwri     = kron(Wy*Ri,Iny);
  graday   = sbeta*(gradAy-kwri*gradBy+ ...
             kron(Wy*Ri,BRi)*gradR-kron(Iny,BRi)*Syu*gradWy);
% gradaz   = sbeta*(gradAz-kron(Wz*Ri,Iny)*gradBy+ ...
%            kron(Wz*Ri,BRi)*gradR-kron(Inz,BRi)*Szu*gradWz);
  gradaz   = sbeta*(kron(Wz*Ri,Iny)*gradBy+ ...
             kron(Wz*Ri,BRi)*gradR-kron(Inz,BRi)*Szu*gradWz);
  gradby   = sbeta*gradBy;

  graday(:,si(14)) = graday(:,si(14)) + 0.5*ay(:)/beta*ind(14); 
  gradaz(:,si(14)) = gradaz(:,si(14)) + 0.5*az(:)/beta*ind(14); 
  gradby(:,si(14)) = gradby(:,si(14)) + 0.5*by(:)/beta*ind(14); 

% gradth   = kron((Wy+thet'*Wz)*Ri,PRWi)*gradPs-kron((Wy+thet'*Wz)*Ri, ...
%            PRWi*Psi*Ri)*gradR+kron(thet',PRWi*Psi*Ri)*Szu*gradWz- ...
%            kron(Iny,PRWi)*gradTh+kron(Iny,PRWi*Psi*Ri)*Syu*gradWy;
  gradth   =-kron((Wy+thet'*Wz)*Ri,PRWi*Psi*Ri)*gradR+ ...
             kron(thet',PRWi*Psi*Ri)*Szu*gradWz+ ...
             kron(Iny,PRWi*Psi*Ri)*Syu*gradWy;
% gradps   = kron(psi'*Wz*Ri,PRWi)*gradPs-kron(psi'*Wz*Ri,PRWi*Psi*Ri)* ...
%            gradR+kron(psi',PRWi*Psi*Ri)*Szu*gradWz-gradPs;
  gradps   =-kron(psi'*Wz*Ri,PRWi*Psi*Ri)*gradR+ ...
             kron(psi',PRWi*Psi*Ri)*Szu*gradWz;
  gradah   = graday+kron(thet',Iny)*gradaz+kron(Iny,az)*gradth;
  gradbh   = gradby+kron(psi',Iny)*gradaz+kron(Inu,az)*gradps;
  gradqy   = gradQy-kwri*gradWy+kron(Wy*Ri,Wy*Ri)*gradR ...
             -kron(Iny,Wy*Ri)*Syu*gradWy;
  gradqz   = gradQz-kron(Wz*Ri,Iny)*gradWy+kron(Wz*Ri,Wy*Ri)*gradR ...
             -kron(Inz,Wy*Ri)*Szu*gradWz;
  gradqh   = gradqy+kron(thet',Iny)*gradqz+kron(Iny,qz)*gradth;
  gradat   = graday-kron(qz*psi*Ri,Iny)*gradby+kron(qz*psi*Ri,by*Ri)* ...
             gradR-kron(qz,by*Ri)*Szu*gradps-kron(Iny,by*Ri*psi')*Syz*gradqz;
  tem      = sparse(eye(ny*ny)-sbeta*kron(Ao',aotil'));
  gradP    = tem\ ...
             (gradqh+kron(Iny,aotil'*P)*gradah+sbeta*kron(Ao'*P',Iny)* ...
             Syy*gradat-kron(f',aotil'*P)*gradbh-sbeta* ...
             kron(Ao'*P',ftil')*Syu*gradby+kron(f',ftil')*gradR);
  gradf    = -kron(f',Rbpi)*gradR+sbeta*kron(Ao'*P',Rbpi)*Syu*gradby ...
             +sbeta*kron(Ao',Rbpi*by')*gradP+kron(Iny,Rbpi*by'*P)* ...
             gradah-kron(f',Rbpi*by'*P)*gradbh;
% gradF    = kron(f'-F',RWPi)*gradR+kron(Iny,RWPi*R)*gradf+kron(Iny,RWPi)* ...
%            Syu*gradWy+kron((Thet-Psi*F)',RWPi)*Szu*gradWz+ ...
%            kron(Iny,RWPi*Wz')*gradTh-kron(F',RWPi*Wz')*gradPs;
  gradF    = kron(f'-F',RWPi)*gradR+kron(Iny,RWPi*R)*gradf+kron(Iny,RWPi)* ...
             Syu*gradWy+kron((Thet-Psi*F)',RWPi)*Szu*gradWz;
  gradAo   = (gradah-kron(f',Iny)*gradbh-kron(Iny,bhat)*gradf)/sbeta;
  gradAo(:,si(14)) = gradAo(:,si(14))-.5*Ao(:)/beta*ind(14);

  %-------------------------------------------------------------------
  % 14.  Specify gradGp.

  gradyk           = zeros(1,np);
  gradyk(si(9))    = yk/a4*ind(9);
  gradyk(si(10))   = gradyk(si(10)) + yk*log(X(13))*ind(10); 
  gradyk           = gradyk + b4*yk/X(13)*gradX(13,:);

  gradyh           = zeros(1,np);
  gradyh(si(9))    = -yh/(1-a4)*ind(9);
  gradyh(si(10))   = gradyh(si(10)) + yh*(X(4)+log(X(14)))*ind(10);
  gradyh           = gradyh + b4*yh/X(14)*gradX(14,:)+b4*yh*gradX(4,:);

  grady            = zeros(1,np);
  grady(si(10))    = -y*log(yk+yh)/b4/b4*ind(10);
  grady            = grady + y/(b4*(yk+yh))*(gradyk+gradyh);

  gradFp           = Suy*gradF;
  gradGp           = zeros(ny*no,np);
  gradGp(13:36,:)  = -gradFp(13:36,:);
  gradGp(61:72,:)  = gradFp(1:12,:); 

  % derivatives of G(4,:)
  tem1             = grady/y+gradyh/yh-(gradyk+gradyh)/(yk+yh);
  tem2             = tem1+gradyk/yk-gradyh/yh;
  g41              = [0,-y*yh*X(4)/(yk+yh),0,y*yh/(yk+yh),zeros(1,8)];
  g42              = -y*yk/(yk+yh)/X(13);
  g43              = -y*yh/(yk+yh)/X(14);
  gradGp(37:48,:)  = g42*F(1,:)'*(tem2-gradX(13,:)/X(13))+ ...
                     g42*gradFp(1:12,:)+g43*F(2,:)'*(tem1- ...
                     gradX(14,:)/X(14))+g43*gradFp(13:24,:);
  gradGp(38,:)     = gradGp(38,:) + g41(2)*(tem1+gradX(4,:)/X(4));
  gradGp(40,:)     = gradGp(40,:) + g41(4)*tem1;
  
  % derivatives of G(5,:)
  g51              = [0,-y*yh*X(4)*X(3)/(yk+yh)-X(3)*y,y,y*yh*X(3)/(yk+yh), ...
                      zeros(1,8)];
  g52              = -y*yk*X(3)/(yk+yh)/X(13);
  g53              = -y*yh*X(3)/(yk+yh)/X(14);
  gradGp(49:60,:)  = g52*F(1,:)'*(tem2-gradX(13,:)/X(13)+gradX(3,:)/X(3))+ ...
                     g53*F(2,:)'*(tem1-gradX(14,:)/X(14)+gradX(3,:)/X(3))+ ...
                     g52*gradFp(1:12,:)+g53*gradFp(13:24,:);
  gradGp(50,:)     = gradGp(50,:) + (X(3)*y+g51(2))*(tem1+gradX(4,:)/X(4)+ ...
                     gradX(3,:)/X(3))-y*gradX(3,:)-X(3)*grady; 
  gradGp(51,:)     = gradGp(51,:) + grady;
  gradGp(52,:)     = gradGp(52,:) + g51(4)*(tem1+gradX(3,:)/X(3));
  gradG            = Syo*gradGp;
  %-------------------------------------------------------------------
  
  gradX       = gradX(1:ny,:);
  dAo         = zeros(ny);
  dG          = zeros(no,ny);
  dC          = zeros(ny);
  dD          = zeros(no);
  dR          = zeros(no);
  dX          = zeros(ny,1);
  dS          = zeros(ny);
  dK          = zeros(ny,no);
  gradS0      = zeros(ny*ny,np);
  gradX0      = gradX;
  tem         = zeros(ny);
  sum1        = zeros(np,1);
  sum2        = sum1;
  tem0        =-[1;0;1;1;1;1;0;0]/grate;

  % derivative of S0, where  S0 = Ao*S0*Ao'+C*C'
% for j=1:np;
%   dAo(:)     = gradAo(:,j);
%   dC(:)      = gradC(:,j);
%   tem2       = dAo*S0*Ao'+dC*C';
%   tem1       = tem2+tem2';
%   tem(nct,nct)= dlyap(Ao(nct,nct),tem1(nct,nct));
%   gradS0(:,j) = tem(:);
% end

  % derivative of S0, where  S0 = S (from kfilter)
  AA          = Ao-K*Gbar;
  IKG         = Iny-K*G;
  for j=1:np
    dAo(:)       = gradAo(:,j);
    dC(:)        = gradC(:,j);
    dG(:)        = gradG(:,j);
    dD(:)        = gradD(:,j);
    dR(:)        = gradRR(:,j);
    tem2         = IKG*dAo*S*AA'+IKG*dC*C'*IKG'-K*dG*S+K*D*dG*S*AA'+ ...
                   K*dD*G*S*AA'+.5*K*dR*K'; 
    tem1         = tem2+tem2';
    tem(nct,nct) = dlyap(AA(nct,nct),tem1(nct,nct));
    gradS0(:,j)   = tem(:);
  end;


  S1        = S0;
  Omega     = Rbar+Gbar*S1*Gbar';
  Omegai    = inv(Omega);
  OG        = Omegai*Gbar;
  M         = Omegai-Omegai*innov(t0,:)'*innov(t0,:)*Omegai;
  xio       = Xt(t0,:)'*innov(t0,:)*Omegai;
  yio       = Y(t0,:)'*innov(t0,:)*Omegai;
  for j=1:np;
    dX(:)   = gradX0(:,j);
    dS(:)   = gradS0(:,j);
    dAo(:)  = gradAo(:,j);
    dG(:)   = gradG(:,j);
    dC(:)   = gradC(:,j);
    dR(:)   = gradRR(:,j);
    dD(:)   = gradD(:,j);
    ss1     = 2*trace(dAo*S1*Gbar'*M*G)+2*trace(dC*C'*G'*M*G) ...
                 +2*trace(dG*(Ao*S1*Gbar'*M-S1*Gbar'*M*D ...
                 +Sigma*G'*M))-2*trace(dD*G*S1*Gbar'*M) ...
                 +trace(dR*M)+trace(dS*Gbar'*M*Gbar);
    ss2     = -2*trace(dAo*xio*G)-2*trace(dG*(Ao*xio-xio*D)) ...
                 -2*trace(dD*(yio-G*xio)) ...
                 -2*trace(dX*innov(t0,:)*OG); 
    if gratein & j==si(15);
      ss2   = ss2+2*trace((((t0+1)*tem0.*Y(t0+1,:)') ...
                            -D*(t0*tem0.*Y(t0,:)'))*innov(t0,:)*Omegai);  
    end;
    dLt(j,t0) = ss1+ss2;
    sum1(j)   = sum1(j) + ss1;
    sum2(j)   = sum2(j) + ss2;
  end
  for i=t0+1:T;
disp(i)
    K1        = (Ao*S1*Gbar'+Sigma*G')*Omegai;
    AA        = Ao-K1*Gbar;
    IKG       = Iny-K1*G;
    tem1      = Ao*S1*Ao'+Sigma-K1*Gbar*S1*Ao'-K1*G*Sigma;
    gradK     = kron(Omegai*G*C,IKG)*gradC ...
               +kron(Omegai*G,IKG*C)*Syy*gradC ...
               +kron(OG*S1,IKG)*gradAo ...
               +kron(Omegai*G,AA*S1)*Syy*gradAo ...
               +kron(OG,AA)*gradS0 ...
               +kron(OG*S1*G',K1)*gradD ...
               -kron(Omegai,AA*S1*G')*Soo*gradD ...
               -kron(K1',K1)*gradG ...
               +kron(OG*S1,K1*D)*gradG ...
               -kron(Omegai*D,AA*S1)*gradGp ...
               +kron(Omegai,tem1)*gradGp ...
               -kron(Omegai,K1)*gradRR;
    gradX0    = AA*gradX0 +kron(Xt(i-1,:),IKG)*gradAo ...
                +kron(innov(i-1,:),Iny)*gradK ...
                -kron(Xt(i-1,:)*Ao',K1)*Syo*gradGp ...
                +kron(Xt(i-1,:),K1*D)*Syo*gradGp ...
                +kron(Xt(i-1,:)*G',K1)*gradD ...
                -kron(Y(i-1,:),K1)*gradD;
    gradX0(:,si(15)) = gradX0(:,si(15))+ ...
                K1*(tem0.*(i*Y(i,:)'-D*(i-1)*Y(i-1,:)'))*ind(15);
    S0        = S1;
    S1        = Ao*S1*Ao'+Sigma-(Ao*S1*Gbar'+Sigma*G')*Omegai* ...
               (Gbar*S1*Ao'+G*Sigma);
    Omega     = Rbar+Gbar*S1*Gbar';
    Omegai    = inv(Omega);
    OG        = Omegai*Gbar;
    M         = Omegai-Omegai*innov(i,:)'*innov(i,:)*Omegai;
    xio       = Xt(i,:)'*innov(i,:)*Omegai;
    yio       = Y(i,:)'*innov(i,:)*Omegai;
    %tem2      = kron(AA*S0,IKG)*gradAo ...
    %            +kron(C-K1*G*C,IKG)*gradC ...
    %            -kron(AA'*S0*Ao'+Sigma*IKG',K1)*Syo*gradGp ...
    %            +kron(AA*S0,K1*D)*Syo*gradGp ...
    %            +kron(AA*S0*G',K1)*gradD ...
    %            +0.5*kron(K1,K1)*gradRR;
    %gradS0    = kron(AA,AA)*gradS0+tem2+Syy*tem2;
    for j=1:np;
      dX(:)   = gradX0(:,j);
      dS(:)   = gradS0(:,j);
      dAo(:)  = gradAo(:,j);
      dG(:)   = gradG(:,j);
      dC(:)   = gradC(:,j);
      dR(:)   = gradRR(:,j);
      dD(:)   = gradD(:,j);
      tem2    = dAo*S0*Ao'+dC*C'-dC*C'*G'*K1'-K1*G*dC*C'-K1*dG*C*C' ...
                -dAo*S0*Gbar'*K1'-K1*dG*Ao*S0*Ao'-K1*G*dAo*S0*Ao' ...
                +K1*dD*G*S0*Ao'+K1*D*dG*S0*Ao'+.5*K1*dR*K1' ...
                +K1*dG*Ao*S0*Gbar'*K1'+K1*G*dAo*S0*Gbar'*K1' ...
                -K1*dD*G*S0*Gbar'*K1' ...
                -K1*D*dG*S0*Gbar'*K1'+K1*dG*C*C'*G'*K1'+K1*G*dC*C'*G'*K1';
      dS      = AA*dS*AA'+tem2+tem2';
      gradS0(:,j) = dS(:);

      ss1     = 2*trace(dAo*S1*Gbar'*M*G)+2*trace(dC*C'*G'*M*G) ...
                   +2*trace(dG*(Ao*S1*Gbar'*M-S1*Gbar'*M*D ...
                   +Sigma*G'*M))-2*trace(dD*G*S1*Gbar'*M) ...
                   +trace(dR*M)+trace(dS*Gbar'*M*Gbar);
      ss2     = -2*trace(dAo*xio*G)-2*trace(dG*(Ao*xio-xio*D)) ...
                   -2*trace(dD*(yio-G*xio)) ...
                   -2*trace(dX*innov(i,:)*OG); 
      if gratein & j==si(15);
        ss2   = ss2+2*trace((((i+1)*tem0.*Y(i+1,:)') ...
                              -D*(i*tem0.*Y(i,:)'))*innov(i,:)*Omegai);  
      end;
      dLt(j,i) = ss1+ss2;
      sum1(j)  = sum1(j) + ss1;
      sum2(j)  = sum2(j) + ss2;
    end
  end
  dL          = sum1+sum2; 
  if gratein;
    dLt(si(15),:)  = dLt(si(15),:) + 5*(nobs+1)/grate;
    dL(si(15))     = dL(si(15))    + 5*(nobs+1)*nobs/grate;
  end;
end;
