function [L,Lt,dL,dLt]=lqse(tstart,tend,param,ind);
%LQSE_1B   Standard Errors for linear quadratic cattle cycle example, v. 1b.

%          Ellen R. McGrattan,  8-24-93

global ZVAR
no     = 4;
x0     = zeros(8,1);
x0(2:5)= [17.56*ones(3,1);1];

alpha0 = param(1);
alpha1 = param(2);
beta   = param(3);
delta  = param(4);
gamma0 = param(5);
gamma1 = param(6);
f1     = param(7);
f3     = param(8);
f5     = param(9);
f7     = param(10);
g      = param(11);
rho1   = param(12);
rho2   = param(13);
rho3   = param(14);
sigma1 = param(15);
sigma2 = param(16);
sigma3 = param(17);
muh    = param(18);
mum    = param(19);
grate  = param(20);
varm(1)= param(21);
varm(2)= param(22);
varm(3)= param(23);
varm(4)= param(24);

f2     = 0.5*gamma0*g/f1;   
f4     = 0.5*gamma1*g/f3;  
f6     = 0.5/f5;
f8     = 0.5/f7;
a22    = [1 0 0 0; 0 rho1 0 0; 0 0 rho2 0; 0 0 0 rho3];
c2     = [0 0 0;sigma1 0 0; 0 sigma2 0; 0 0 sigma3];
deltah = 0.0;
deltak = [(1-delta),0,g;1,0,0;0,1,0];
gamma  = [0 0 0; f5*(1-delta) 0 g*f5; f1 0 0; 0 f3 0; 0 0 0];
lambda = 0.0;
phic   = [1 f5 0 0 -f7]';
phig   = [0 0 0 0; 1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1];
phii   = [1 0 0 0 0]';
pih    = sqrt(1/alpha1);
thetah = 0.0;
thetak = [1 0 0]';
ub     = [pih*alpha0 0 0 pih];
ud     = [0 0 0 0; [f6;f2;f4]*[muh,1,0,0]; f8*[mum,0,1,0]];
Omega  = diag(varm);
Gam    = zeros(no);

Theta  = param(ind);
np     = sum(ind);
fin    = find(ind);
%---------------------------------------------------------------------


%---------------------------------------------------------------------
% 4. The position of the discount factor (beta) and the growth rate 
%    (grate) are set so that their derivatives can be computed directly.
%    Set betain = i   if the ith element of Theta corresponds to beta 
%               = 0   if beta is not to be estimated               
%    and do the same for gratein.

if ind(3);  betain=sum(ind(1:3));   else; betain=0;  end;
if ind(20); gratein=sum(ind(1:20)); else; gratein=0; end;
%---------------------------------------------------------------------



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

%---------------------------------------------------------------------

% ordering: alpha0,alpha1,beta,delta,gamma0,gamma1,f1,f3,f5,f7,
%           g,rho1,rho2,rho3,sigma1,sigma2,sigma3,muh,mum,
%           grate,varm
Ub = [10000;100;.999;1;10;10;1.;1.;1.;1.;1;.999;.999;.999;1000;1000;
      1000;5000;5000;1.1;10000;10000;10000;10000];
Lb = [0;.00000001;.001;0;0;0;0.0000001;0.0000001;0.0000001;0.0000001;0;
     -.999;-.999;-.999;0;0;0;0;0;.999;0;0;0;0];
Ub = Ub(ind);
Lb = Lb(ind);
if any(Theta>Ub) | any(Theta<Lb);
  L  = 1e+20; dL =1e+20*ones(size(Theta));
  Lt = 1e+20*ones(tend-tstart+1); dLt=1e+20*ones(length(Theta),tend-tstart+1);
  return
end; 

%---------------------------------------------------------------------
% 6. Specify element in x corresponding to x(t)=1, all t: ct.  
              
ct=5;
%---------------------------------------------------------------------

[nb,nh] = size(lambda);
[nd,nc] = size(phic);
[nz,nw] = size(c2);
[tem,ng]= size(phig);
[nk,ni] = size(thetak);
nx      = nh+nk+nz;
n       = ni+nx;

phiin   = inv([phic phig]);
phiinc  = [eye(nc) zeros(nc,ng)]*phiin;
phiing  = [zeros(ng,nc) eye(ng)]*phiin;

a1  = thetah*phiinc*gamma;
a12 = [thetah*phiinc*ud; zeros(nk,nz)];
a11 = [deltah a1 ; zeros(nk,nh) deltak];
b11 = -thetah*phiinc*phii;
b1  = [b11; thetak];
c   = [zeros(nk+nh,nw);c2];
G1  = phiing*[zeros(nd,nh) gamma ud];
G2  = -phiing*phii;
H1  = [lambda pih*phiinc*gamma -ub+pih*phiinc*ud];
H2  = -pih*phiinc*phii; 
Q   = -.5*(G1'*G1 + H1'*H1);
R   = -.5*(G2'*G2 + H2'*H2);
W   = -.5*(G1'*G2 + H1'*H2);
A   = [a11 a12; [zeros(nz,nk+nh) a22]];
B   = [b1; zeros(nz,ni)];
a   = sqrt(beta)*(A-B*(R\W'));
b   = sqrt(beta)*B;
q   = Q-W*(R\W');
Sigma = c*c';

%---------------------------------------------------------------------
% 7. Use one of the following methods for computing F,P,Ao=A-B*F:
%    doubling algorithm, Vaughan's method, direct iteration.
%    Comment out or delete sections not to be used (after checking
%    that all algorithms provide same result).

%=======================
% 7a. Doubling Algorithm
%=======================

%v = eye(nx);
%a0= a;
%b0= b*(R\b');
%P = q;
%k0= a'*P*b/(b'*P*b+R);
%dd=1;
%it=1;
%maxit=1000;
%tol=1e-5;                         % reset tol if necessary
%while (dd>tol & it<=maxit);
%  v1 = (v+b0*P);
%  vv = v1\a0;
%  a1 = a0*vv;
%  b0 = b0+a0*(v1\(b0*a0'));
%  P  = P+a0'*P*vv;
%  k1 = a'*P*b/(b'*P*b+R);
%  dd = norm(k1-k0);
%  it = it+1;
%  a0 = a1;
%  k0 = k1;
%end
%F = k1'+R\W';

%=====================
%% 7b. Vaughan's method
%=====================

%[v,d] = eig([a+b/R*b'/a'*q  -b/R*b'/a'; -a'\q  inv(a)']);
%d     = diag(d);
%[d,index] = sort(abs(d));	 % sort on magnitude of eigenvalues
%if (~((d(nx) < 1) & (d(nx+1)>1)))
%  error('Can''t order eigenvalues, (A,B) may be uncontrollable.')
%end
%%select vectors with eigenvalues inside unit circle
%v11 = v(1:nx,index(1:nx));
%v21 = v((nx+1):(2*nx),index(1:nx));
%P   = real(v21/v11);
%F   = (R+b'*P*b)\b'*P*a+R\W';

%==================================================
%% 7c. Vaughan's method when A-B*R\W' not invertible
%==================================================

tem1 = [eye(nx),b/R*b';zeros(nx),a'];
tem2 = [a,zeros(nx);-q,eye(nx)];
[v,d]=eig(tem2,tem1);
d = diag(d);
[d,index] = sort(abs(d));        % sort on magnitude of eigenvalues
if (~((d(nx) < 1) & (d(nx+1) > 1)))
        error('Can''t order eigenvalues, (A,B) may be uncontrollable.')
end
% select vectors with eigenvalues inside unit circle
v11 = v(1:nx,index(1:nx));
v21 = v((nx+1):(2*nx),index(1:nx));
P   = real(v21/v11);
F   = (R+b'*P*b)\b'*P*a+R\W';

%=====================
%% 7d. Direct Iteration 
%=====================
% 
%P = -eye(nx);
%F0= (R+beta*B'*P*B)\(beta*B'*P*A+W');
%d=1;
%tol=1e-8;       % reset tolerance if necessary 
%while (d > tol)
%  P  = beta*A'*P*A + Q -(beta*A'*P*B+W)*F0;
%  F1 = (R+beta*B'*P*B)\(beta*B'*P*A+W');
%  d  = norm(F1-F0);
%  F0 = F1;
%end;
%F  = F1;
%---------------------------------------------------------------------

Ao = A-B*F;


%---------------------------------------------------------------------
% 8. Specify observables (y). (Use global statement on ZVAR instead of 
%    loading the dataset each time.)

%data y(t)=[total stock(t),breeding stock(t),slaughter rate(t)]
T   = length(ZVAR);
tem = grate.^([1:T]');  

y   = [ZVAR(:,1)./tem,ZVAR(:,2)./tem,ZVAR(:,3)./tem];
if no>3; 
  y = [y,ZVAR(:,4)+100-ZVAR(1,4)];
end;
%---------------------------------------------------------------------

T   = length(y);
dy  = y(2:T,:)-y(1:T-1,:)*Gam; 
T   = T-1;


%---------------------------------------------------------------------
% 9. Specify D which relates observables to states (e.g.,y = D x + v).   

D  = [[0,1-delta+g,g,g,0,0,0,0]-F;
      [0,1-delta,0,g,0,0,0,0]-F; 
      F];
skl = [zeros(nk,nh),eye(nk),zeros(nk,nz)];
sc  = phiinc*(gamma*skl+[zeros(nd,nh+nk),ud]+phii*F);
mc  = 2*beta*thetah'*[eye(nh),zeros(nh,nk+nz)]*P*Ao+ ...
      pih'*([zeros(nb,nh+nk),ub]-lambda* ...
      [eye(nh),zeros(nh,nk+nz)]-pih*sc);
if no>3; D= [D;mc]; end;
%---------------------------------------------------------------------

Dp  = D';
Dbar=D*Ao-Gam*D;
Obar=Omega+D*Sigma*D';
[K,S]=kfilter(Ao,Dbar,Sigma,Obar,Sigma*D');
Sigmaz=Obar+Dbar*S*Dbar';
Sigmazi=inv(Sigmaz);

t0= 1;
xt(1,:)=x0';
innov(1,:)=dy(1,:)-x0'*Dbar';
for i=2:T;
  xt(i,:)=xt(i-1,:)*Ao'+innov(i-1,:)*K';
  innov(i,:)=dy(i,:)-xt(i,:)*Dbar';
end;
sum1=innov(1:T,:)'*innov(1:T,:)/T;

L= T*(log(det(Sigmaz))+trace(Sigmazi*sum1));
for t=tstart:tend;
  Lt(t,1)=log(det(Sigmaz))+innov(t,:)*Sigmazi*innov(t,:)';
end;
Lt=Lt(tstart:tend);

%-----------------------------------------------------------------------
% 10. Set jacobian resulting from observations being a function of Theta,
%     e.g., if n elements of y have a geometric trend removed then add
%     n*log(grate)*T*(T+1) to L.

if gratein;
  L  = L + 3*log(grate)*T*(T+1);
  Lt = Lt + 3*log(grate)*(T+1);
end;
%-----------------------------------------------------------------------


if nargout > 2;

  %-------------------------------------------------------------------
  % 11. Specify grada2,gradc2,graddh,graddk,gradg,gradl,gradpc,gradpg,
  %             gradpi,gradph,gradth,gradtk,gradub,gradud,gradx0,
  %             gradG,gradO  
  %     which are gradients for
  %             a22,c2,deltah,deltak,gamma,lambda,phic,phig,phii,pih,
  %             thetah,thetak,ub,ud,x0,Gam,Omega,
  %     e.g., gradA(:,i) = vec(dA/dTheta(i)).
  %    
  p      = length(ind);
  grada2 = zeros(nz*nz,p);    grada2(6,12)=1; grada2(11,13)=1; grada2(16,14)=1; 
  gradc2 = zeros(nz*nw,p);    gradc2(2,15)=1; gradc2(7,16)=1;  gradc2(12,17)=1; 
  graddh = zeros(nh*nh,p);   
  graddk = zeros(nk*nk,p);    graddk(1,4)=-1; graddk(7,11)=1;
  gradg  = zeros(nd*nk,p);    gradg(2,4)=-f5; gradg(2,9)=1-delta; 
                              gradg(3,7)=1;   gradg(9,8)=1;    gradg(12,9)=g;
                              gradg(12,11)=f5;
  gradl  = zeros(nb*nh,p);    
  gradpc = zeros(nd*nc,p);    gradpc(2,9)=1;  gradpc(5,10)=-1;
  gradpg = zeros(nd*ng,p); 
  gradpi = zeros(nd*ni,p); 
  gradph = zeros(nb*nc,p);    gradph(1,2)=-0.5*alpha1^(-1.5);
  gradth = zeros(nh*nc,p);    
  gradtk = zeros(nk*ni,p);    
  gradub = zeros(nb*nz,p);    gradub(1,1)=pih;
                              gradub(1,2)=-0.5*alpha1^(-1.5)*alpha0;
                              gradub(4,2)=-0.5*alpha1^(-1.5);
  gradud = zeros(nd*nz,p);    gradud(2,9)=-f6*muh/f5;          
                              gradud(2,18)=f6;
                              gradud(7,9)=-f6/f5;         
                              gradud(3,5)= f2*muh/gamma0;    
                              gradud(3,7)=-f2*muh/f1;          
                              gradud(3,11)=f2*muh/g;
                              gradud(3,18)=f2;
                              gradud(8,5)= f2/gamma0;         
                              gradud(8,7)=-f2/f1;         
                              gradud(8,11)=f2/g;
                              gradud(4,6)= f4*muh/gamma1;    
                              gradud(4,8)=-f4*muh/f3;          
                              gradud(4,11)=f4*muh/g;
                              gradud(4,18)=f4;
                              gradud(9,6)= f4/gamma1;         
                              gradud(9,8)=-f4/f3;         
                              gradud(9,11)=f4/g;
                              gradud(5,10)=-f8*mum/f7;          
                              gradud(5,19)=f8;
                              gradud(15,10)=-f8/f7;         
  gradx0 = zeros(nx,p);
  gradG  = zeros(no*no,p);
  gradO  = zeros(no*no,p);    gradO(1,21)=1; gradO(no+2,22)=1; 
                              gradO(2*no+3,23)=1; 
                              if no>3; gradO(3*no+4,24)=1; end;
  %-------------------------------------------------------------------

  grada2 = grada2(:,ind);
  gradc2 = gradc2(:,ind);
  graddh = graddh(:,ind);
  graddk = graddk(:,ind);
  gradg  = gradg(:,ind);
  gradl  = gradl(:,ind);
  gradpc = gradpc(:,ind);
  gradpg = gradpg(:,ind);
  gradpi = gradpi(:,ind);
  gradph = gradph(:,ind);
  gradth = gradth(:,ind);
  gradtk = gradtk(:,ind);
  gradub = gradub(:,ind);
  gradud = gradud(:,ind);
  gradx0 = gradx0(:,ind);
  gradG  = gradG(:,ind);
  gradO  = gradO(:,ind);
  gradp  = [gradpc;gradpg];
  p      = sum(ind);

  tem1=zeros(nx,ni);
  tem1(:)=[1:nx*ni]';
  tem2=tem1';
  tem2=tem2(:);
  tem1=tem1(:);
  Kxi=zeros(nx,ni);
  for i=1:nx*ni;
    Kxi(tem1(i),tem2(i))=1;
  end;
  tem1=zeros(ni,nx);
  tem1(:)=[1:nx*ni]';
  tem2=tem1';
  tem2=tem2(:);
  tem1=tem1(:);
  Kix=zeros(ni,nx);
  for i=1:nx*ni;
    Kix(tem1(i),tem2(i))=1;
  end;
  tem1=zeros(nx);
  tem1(:)=[1:nx*nx]';
  tem2=tem1';
  tem2=tem2(:);
  tem1=tem1(:);
  Kxx=zeros(nx*nx);
  for i=1:nx*nx;
    Kxx(tem1(i),tem2(i))=1;
  end;
  tem1=zeros(ni);
  tem1(:)=[1:ni*ni]';
  tem2=tem1';
  tem2=tem2(:);
  tem1=tem1(:);
  Kii=zeros(ni*ni);
  for i=1:ni*ni;
    Kii(tem1(i),tem2(i))=1;
  end;
  tem1=zeros(nh,nd);
  tem1(:)=[1:nh*nd]';
  tem2=tem1';
  tem2=tem2(:);
  tem1=tem1(:);
  Khd=zeros(nh,nd);
  for i=1:nh*nd;
    Khd(tem1(i),tem2(i))=1;
  end;
  tem1=zeros(nz,nw);
  tem1(:)=[1:nz*nw]';
  tem2=tem1';
  tem2=tem2(:);
  tem1=tem1(:);
  Kzx=zeros(nz,nw);
  for i=1:nz*nw;
    Kzw(tem1(i),tem2(i))=1;
  end;
  tem1=zeros(nk,ni);
  tem1(:)=[1:nk*ni]';
  tem2=tem1';
  tem2=tem2(:);
  tem1=tem1(:);
  Kki=zeros(nk,ni);
  for i=1:nk*ni;
    Kki(tem1(i),tem2(i))=1;
  end;
  tem1=zeros(ng,nx);
  tem1(:)=[1:ng*nx]';
  tem2=tem1';
  tem2=tem2(:);
  tem1=tem1(:);
  Kgx=zeros(ng,nx);
  for i=1:ng*nx;
    Kgx(tem1(i),tem2(i))=1;
  end;
  tem1=zeros(nb,nx);
  tem1(:)=[1:nb*nx]';
  tem2=tem1';
  tem2=tem2(:);
  tem1=tem1(:);
  Kbx=zeros(nb,nx);
  for i=1:nb*nx;
    Kbx(tem1(i),tem2(i))=1;
  end;
  tem1=zeros(nd,nx);
  tem1(:)=[1:nd*nx]';
  tem2=tem1';
  tem2=tem2(:);
  tem1=tem1(:);
  Kdx=zeros(nd,nx);
  for i=1:nd*nx;
    Kdx(tem1(i),tem2(i))=1;
  end;
  tem1=zeros(ng,ni);
  tem1(:)=[1:ng*ni]';
  tem2=tem1';
  tem2=tem2(:);
  tem1=tem1(:);
  Kgi=zeros(ng,ni);
  for i=1:ng*ni;
    Kgi(tem1(i),tem2(i))=1;
  end;
  tem1=zeros(nb,ni);
  tem1(:)=[1:nb*ni]';
  tem2=tem1';
  tem2=tem2(:);
  tem1=tem1(:);
  Kbi=zeros(nb,ni);
  for i=1:nb*ni;
    Kbi(tem1(i),tem2(i))=1;
  end;
  tem1=zeros(nc,nh);
  tem1(:)=[1:nc*nh]';
  tem2=tem1';
  tem2=tem2(:);
  tem1=tem1(:);
  Kch=zeros(nc,nh);
  for i=1:nc*nh;
    Kch(tem1(i),tem2(i))=1;
  end;
  tem1=zeros(nc,nb);
  tem1(:)=[1:nc*nb]';
  tem2=tem1';
  tem2=tem2(:);
  tem1=tem1(:);
  Kcb=zeros(nc,nb);
  for i=1:nc*nb;
    Kcb(tem1(i),tem2(i))=1;
  end;

  tem = kron(phiinc',eye(nh))*gradth-kron(phiin',thetah*phiinc)*gradp;
  grada12 = kron(gamma',eye(nh))*tem+kron(eye(nk),thetah*phiinc)*gradg;
  grada13 = kron(ud',eye(nh))*tem+kron(eye(nz),thetah*phiinc)*gradud;
  gradBp  = [-kron(thetah*phiinc,eye(ni))*Kii*gradpi-kron(eye(nh),phii')* ...
               Khd*tem; Kki*gradtk; zeros(nz*ni,p)];
  grads   = (kron(eye(nz),c2)*Kzw+kron(c2,eye(nz)))*gradc2;
         
  tem1 = zeros(nx);
  tem2 = zeros(nx);
  tem3 = zeros(nh);
  tem4 = zeros(nh,nk);
  tem5 = zeros(nh,nz);
  tem6 = zeros(nk);
  tem7 = zeros(nz);
  for i=1:p;
    tem3(:) = graddh(:,i);
    tem4(:) = grada12(:,i);
    tem5(:) = grada13(:,i);
    tem6(:) = graddk(:,i);
    tem7(:) = grada2(:,i);
    tem1(1:nh,1:nh)             = tem3;
    tem1(1:nh,nh+1:nh+nk)       = tem4;
    tem1(1:nh,nh+nk+1:nx)       = tem5;
    tem1(nh+1:nh+nk,nh+1:nh+nk) = tem6;
    tem1(nh+nk+1:nx,nh+nk+1:nx) = tem7;
    tem7(:) = grads(:,i);
    tem2(nh+nk+1:nx,nh+nk+1:nx) = tem7;
    gradA(:,i)  = vec(tem1);
    gradS(:,i)  = vec(tem2);
  end;
  gradB = Kix*gradBp;
  gradG1=-kron( [zeros(nd,nh),gamma,ud]'*phiin',phiing)*gradp+ ...
          kron(eye(nx),phiing)*[zeros(nd*nh,p);gradg;gradud];
  gradG2= kron(phii'*phiin',phiing)*gradp-kron(eye(ni),phiing)* ...
          gradpi;
  gradH1=[gradl;
          kron(gamma'*phiinc',eye(nb))*gradph-kron(gamma'*phiin', ...
             pih*phiinc)*gradp+kron(eye(nk),pih*phiinc)*gradg;
          kron(ud'*phiinc',eye(nb))*gradph-kron(ud'*phiin', ...
             pih*phiinc)*gradp-gradub+kron(eye(nz),pih*phiinc)*gradud];
  gradH2=-kron(phii'*phiinc',eye(nb))*gradph+kron(phii'*phiin', ...
          pih*phiinc)*gradp-kron(eye(ni),pih*phiinc)*gradpi;
  gradQ=-.5*( (kron(G1',eye(nx))*Kgx+kron(eye(nx),G1'))*gradG1+ ...
              (kron(H1',eye(nx))*Kbx+kron(eye(nx),H1'))*gradH1 );
  gradR=-.5*( (kron(G2',eye(ni))*Kgi+kron(eye(ni),G2'))*gradG2+ ...
              (kron(H2',eye(ni))*Kbi+kron(eye(nc),H2'))*gradH2 );
  gradW=-.5*(  kron(G2',eye(nx))*Kgx*gradG1+kron(eye(ni),G1')*gradG2+ ...
               kron(H2',eye(nx))*Kbx*gradH1+kron(eye(ni),H1')*gradH2 );
  ri=inv(R+beta*B'*P*B);
  ai=inv(eye(nx*nx)-beta*kron(Ao',Ao'));
  gradP = ai*(gradQ+(beta*kron(eye(nx),Ao'*P)+beta*kron(Ao'*P,eye(nx)) ...
         *Kxx)*gradA-beta*kron(F',Ao'*P)*gradB-beta*kron(Ao'*P,F')*gradBp);
  gradP = gradP+ai*(kron(F',F')*gradR-(kron(F',eye(nx))+kron(eye(nx),F')* ...
          Kxi)*gradW);
  gradF = beta*kron(eye(nx),ri*B'*P)*gradA-beta*kron(F',ri*B'*P)* ...
          gradB+beta*kron(Ao'*P,ri)*gradBp-kron(F',ri)*gradR+ ...
          kron(eye(nx),ri)*Kxi*gradW+ beta*kron(Ao',ri*B')*gradP;
  if betain;
    gradF(:,betain)=gradF(:,betain)-beta*kron(Ao',ri*B')*gradP(:,betain) ...
                   +vec(ri*B'*P*Ao);
    gradP(:,betain)=gradP(:,betain)+ai*vec(Ao'*P*Ao);
    gradF(:,betain)=gradF(:,betain)+beta*kron(Ao',ri*B')*gradP(:,betain);
  end;
  gradAo=gradA-kron(F',eye(nx))*gradB-kron(eye(nx),B)*gradF;

  %-------------------------------------------------------------------
  % 12. Specify gradDp where Dp = D'.   In doing so, you may find
  %     the following formulas useful:  
  %
  %     vec(A*B*C) = kron(C',A)*vec(B)
  %     gradA(:,i) = vec(dA/dTheta(i)).

  gradDp=zeros(24,p);
  if ind(4); 
    gradDp(2,sum(ind(1:4)))=-1; 
    gradDp(10,sum(ind(1:4)))=-1; 
  end;
  if ind(11);
    gradDp([2,3,4,12],sum(ind(1:11)))=ones(4,1);
  end;
  gradDp=gradDp+kron([-1;-1;1],eye(nx))*Kxi*gradF;
  tem   =phiin*(gamma*skl+[zeros(nd,nh+nk),ud]+phii*F);
  gradsc=-kron(tem',phiinc)*gradp+kron(skl',phiinc)*gradg+ ...
         kron(eye(nx),phiinc)*[zeros(nd*nh+nd*nk,p);gradud]+ ...
         kron(F',phiinc)*gradpi+kron(eye(nx),phiinc*phii)*gradF;
  tem   =[eye(nh),zeros(nh,nk+nz)];
  gradmc=2*beta*(kron(Ao'*P'*tem',eye(nc))*Kch*gradth+ ...
         kron(Ao',thetah'*tem)*gradP+kron(eye(nx),thetah'*tem*P)*gradAo);
  tem   =[zeros(nb,nh+nk),ub]-lambda*[eye(nh),zeros(nh,nk+nz)]-pih*sc;
  gradmc=gradmc+kron(tem',eye(nc))*Kcb*gradph+kron(eye(nx),pih')* ...
         [zeros(nb*nh+nb*nk,p);gradub]- ...
         kron([eye(nh),zeros(nh,nk+nz)]',eye(nb))*gradl;
  gradmc=gradmc-kron(sc',pih')*gradph-kron(eye(nx),pih'*pih)*gradsc;
  if no>3;
    gradDp(25:32,:)=gradmc;
    if betain;
      tem   =[eye(nh),zeros(nh,nk+nz)];
      gradDp(25:32,betain) = gradDp(25:32,betain)+2*vec(thetah'*tem*P*Ao);
    end;
  end;
  
  %-------------------------------------------------------------------

  AA=Ao-K*Dbar;
  SD=Sigmazi*Dbar;
  lamb(T,:) = innov(T,:)*SD;
  for i=T-1:-1:1;
    lamb(i,:) = lamb(i+1,:)*AA+innov(i,:)*SD;  
  end;

  sum2 = xt(1:T,:)'*innov(1:T,:)/T;
  sum3 = xt(1:T-1,:)'*lamb(2:T,:)/T;
  sum4 = innov(1:T-1,:)'*lamb(2:T,:)/T;
  sum5 = y(1:T,:)'*innov(1:T,:)/T;
  sum6 = y(1:T-1,:)'*lamb(2:T,:)/T;

  M1=Sigmazi-Sigmazi*sum1*Sigmazi;
  M2=Dbar'*M1*Dbar-2*SD'*sum4*AA;
  M2=.5*(M2+M2');
  nct=ones(nx,1); nct(ct)=0;
  M3=zeros(nx);
  M3(nct,nct)=dlyap(AA(nct,nct)',M2(nct,nct));
  M=M1+K'*M3*K;
  D1=2*( (-sum3-S*SD'*sum4)*(eye(size(S))-K*D)+S*Dbar'*M*D- ...
     (sum2+S*AA'*sum4')*Sigmazi*D+S*AA'*M3-S*Ao'*M3*K*D );
  D2=2*( Sigma*D'*M+(Ao*S*Dbar'+Sigma*D')*Sigmazi*sum4*K-Ao*S*AA'*(sum4' ...
     *Sigmazi+M3*K)+S*AA'*sum4'*Sigmazi*Gam-Sigma*(eye(size(S))-D'*K')* ...
     sum4'*Sigmazi);
  D2=D2+2*(-S*Dbar'*(M1+Sigmazi*sum4*K)*Gam+Ao*S*Dbar'*M1-Ao*sum2* ...
     Sigmazi+sum2*Sigmazi*Gam-sum3*K*Gam+Ao*sum3*K-Sigma*M3*K+S*AA'*M3*K*Gam); 
  D3=2*( .5*D'*M*D-D'*Sigmazi*sum4*(eye(size(S))-K*D)+.5*M3-D'*K'*M3 );
  D4=2*Sigmazi*sum4*K+M;
  D5=2*(-D*S*Dbar'*(M+Sigmazi*sum4*K)-(sum5-D*sum2)*Sigmazi-D*sum3*K+ ...
     D*S*(AA'*sum4'*Sigmazi+Ao'*M3*K)+sum6*K );
  tem1=zeros(size(Ao));
  tem2=zeros(size(D'));
  tem3=zeros(size(Sigma));
  tem4=zeros(size(Omega));
  tem5=zeros(size(Gam));
  for j=1:p;
    tem1(:) = gradAo(:,j);
    tem2(:) = gradDp(:,j);
    tem3(:) = gradS(:,j);
    tem4(:) = gradO(:,j);
    tem5(:) = gradG(:,j);
    dL(j,1) = trace(tem1*D1+tem3*D3)+trace(tem2'*D2+tem4*D4+tem5*D5) ...
              -2*trace(gradx0(:,j)*lamb(1,:)/T);
  end; 

  %-------------------------------------------------------------------
  % 13. Specify derivatives resulting from observations being functions
  %     of grate.

  if gratein;
    sum7 =-( ([1;1;1;0]*[2:T+1]).*y(2:T+1,:)' )*innov(1:T,:)/(grate*T);
    sum8 =-( ([1;1;1;0]*[1:T]).*y(1:T,:)' )*innov(1:T,:)/(grate*T);
    sum9 =-( ([1;1;1;0]*[1:T-1]).*y(1:T-1,:)' )*lamb(2:T,:)/(grate*T);
    sum10=-( ([1;1;1;0]*[2:T]).*y(2:T,:)' )*lamb(2:T,:)/(grate*T);
    dL(gratein)=dL(gratein)-2*trace(K*sum10-K*Gam*sum9)+ ...
            2*trace(sum7*Sigmazi-Gam*sum8*Sigmazi) +3*(T+1)/grate;
  end; 
  %-------------------------------------------------------------------

  dL=T*dL;
  for t=tstart:tend;
    sum1=innov(t,:)'*innov(t,:);
    sum2=xt(t,:)'*innov(t,:);
    sum6=y(t,:)'*innov(t,:);
    %-----------------------------------------------------------------
    if gratein;
      %disp('WARNINGam: derivatives in dLt for grate have not been checked');
      sum8=-( (t+1)*[1;1;1;0].*y(t+1,:)' )*innov(t,:)/grate;
      sum9=-( t*[1;1;1;0].*y(t,:)' )*innov(t,:)/grate;
    end;
    %-----------------------------------------------------------------
    M1=Sigmazi-Sigmazi*sum1*Sigmazi;
    d1=2*( S*Dbar'*M1*D-sum2*Sigmazi*D );
    d2=2*( Sigma*D'*M1-S*Dbar'*M1*Gam+Ao*S*Dbar'*M1-Ao*sum2* ...
       Sigmazi+sum2*Sigmazi*Gam );
    d3=D'*M1*D;
    d4=M1;
    d5=2*(-D*S*Dbar'*M1-(sum6-D*sum2)*Sigmazi);
 
    tem1=zeros(size(Ao));
    tem2=zeros(size(D'));
    tem3=zeros(size(Sigma));
    tem4=zeros(size(Omega));
    tem5=zeros(size(Gam));
    for j=1:np;
      tem1(:)=gradAo(:,j);
      tem2(:)=gradDp(:,j);
      tem3(:)=gradS(:,j);
      tem4(:)=gradO(:,j);
      tem5(:)=gradG(:,j);
      OM=tem1*S*AA'-Ao*S*tem1'*D'*K'+K*D*tem1*S*Dbar'*K';
      OM=OM-Ao*S*(Ao'*tem2*K'-tem2*Gam'*K')-Sigma*tem2*K'+K*tem2'*Sigma* ...
         D'*K'+K*(tem2'*Ao-Gam*tem2')*S*Dbar'*K';
      OM=OM-tem3*D'*K'+.5*tem3+.5*K*D*tem3*D'*K' +.5*K*tem4*K';
      OM=OM+Ao*S*D'*tem5'*K'-K*tem5*D*S*Dbar'*K';
      dS=zeros(nx);
      OM=OM+OM';
      dS(nct,nct)=dlyap(AA(nct,nct),OM(nct,nct));
      dDbar=tem2'*Ao+D*tem1-tem5*D-Gam*tem2';
      dK=( (eye(size(S))-K*D)*(tem3*D'+Sigma*tem2)-K*tem2'*Sigma*D'+ ...
         tem1*S*Dbar'+AA*dS*Dbar'+AA*S*dDbar'-K*dDbar*S*Dbar' ...
         -K*tem4 )*Sigmazi;
      s(t0,:)=gradx0(:,j)';
      for i=t0+1:t;
        s(i,:)=s(i-1,:)*AA'+xt(i-1,:)*(tem1-K*(tem2'*Ao+D*tem1-tem5*D- ...
               Gam*tem2'))'+innov(i-1,:)*dK' - y(i-1,:)*tem5'*K';
        %---------------------------------------------------------------
        if gratein & j==gratein;
          s(i,:)=s(i,:)+( [1,1,1,0].*(-i*y(i,:)+(i-1)*y(i-1,:)*Gam'))* ...
                 K'/grate;
        end;
        %---------------------------------------------------------------
      end;
      dLt(j,t)=trace(tem1*d1+tem3*d3)+trace(tem2'*d2+tem4*d4+tem5*d5) ...
              -2*trace(s(t,:)'*innov(t,:)*SD) + trace(dS*Dbar'*M1*Dbar);
    end;
    %-------------------------------------------------------------------
    if gratein;
      dLt(gratein,t)=dLt(gratein,t)+2*trace(sum8*Sigmazi-Gam*sum9*Sigmazi) ...
                   +5*(T+1)/grate;
    end;
    %-------------------------------------------------------------------
  end; 
  dLt=dLt(:,tstart:tend);
end;

