function [L,dL,Xt,innov,Y,DY,Cbar,C,Ao]=mlehs4(Theta,dummy);
%MLE   Maximum likelihood estimation for nondistorted economies 
%      [L,dL]=mle(Theta) has user-defined sections and thus must
%      be edited.  In mle.m, the class of economies is described
%      and instructions for the user are provided.  See also
%      LQMLE.M (linear-quadratic economies), HSMLE.M (Hansen-
%      Sargent economies), DISTMLE.M (distorted economies).

%      Ellen R. McGrattan,  3-14-90
%      Revised, ERM, 3-14-90

%      Model:      max  E[   sum beta^t r(X[t],u[t]) | X[0] ] 
%                 {u[t]}   t=0:inf 
%
%                  subject to  X[t+1]=A*X[t]+B*u[t]+e[t+1]
%                              X[0] given
%
%      Approximation:  
%                  r(X[t],u[t]) = X[t]'Q X[t]+u[t]'R u[t]+2X[t]'W u[t]  
%                                                 + higher order terms 
%           
%      Equilibrium in Approximated Economy:
%                  u[t] = -F X[t] 
%
%      Transition function for States in Approximated Economy:
%                  X[t+1] = (A-B*F) X[t] + e[t+1]
%                         =  Ao X[t] +e[t+1]
%        
%      Observations:  z(1),z(2),...z(T)
%
%                  Y[t] = C X[t] + v[t]   
%                  Y[t] = g(z[t],Theta) 
%                  v[t] = G v[t-1] + w[t]    
%
%                  E[ e[t] e[s]'] = Sigma, t=s
%                                 = 0      t~=s
%                  E[ w[t] w[s]'] = Omega, t=s
%                                 = 0      t~=s
%
%      Gaussian Likelihood Function:
%                  L= T*( ln(|Sigmaz| + trace(inv(Sigmaz)*Sum) )
%
%                  Sigmaz = E[      ]
%                  Sum    = sum(innov(t)*innov(t)')/T
%                          t=1:T
%
%
%      Miscellaneous: 
%                  x[t]  = set of variables in return function
%                        = D*X[t] + E*u[t], D,E are selector matrices
%                  param = vector of parameters needed to evaluate
%                          the return function 
%                  rfn   = a string giving name of m-file with 
%                          inputs x and param and outputs: 
%                          b   = dr(x,param)/x 
%                          q   = d^2r(x,param)/dx^2
%                          cub = d^3r(x,param)/dx^3
%                          du  = dr(x,param)/dparam
%                          db  = db/dparam 
%                          dq  = dq/dparam
%              
%                  
%      User specifies:    a. Theta
%                         b. 

%      References:
%      [1] Wilson, David A., Alok Kumar, "Derivative Computations
%          for the Log Likelihood Function," IEEE Transactions on
%          Automatic Control, Vol. AC-27, No. 1, pp. 230-232, 1982. 
%      [2] Zadrozny, Peter, "Analytic Derivatives for Estimation of
%          Discrete-Time, Linear-Quadratic, Dynamic, Optimization
%          Models," Econometrica, 56:467-472, 1988. 

Theta=Theta(:);

%---------------------------------------------------------------------
% 1. Provide default values of the parameters.                

d=1e-12; 
h1=.432; 
h2=.0639; 
theta= 3.653089857473236e-01;
aa=3.62; 
delta=.0241; 
lambda=1.004245039075813254;
rho0= 3.537426122486124e+00;
rho1= 1.287634205016144e+00;
rho2=-3.024683392896476e-01;
rho0=     1.076832677675266e+01;
rho1=     6.596511919100013e-01;
rho2=     2.955508450620621e-01;
sigma=1.64;
beta=.995; 
varm=[1000;.001;.001;1000;1000];
coefm=[ 9.949999999999999e-01
        5.859207040009263e-01
        9.495843032926892e-01
        8.696211375605091e-01
        9.505603152629579e-01];
coeft=[ -1.160012222059382e+01 
         0
         0
         1.082789939864638e+00 
        -1.427844644127792e+00]; 
coeft=zeros(5,1);
coefm=zeros(5,1);
varm=[100;1e-8;1e-8;100;100];
%varm=[1;1e-8;1e-8;1;1];
%---------------------------------------------------------------------

Tsize=ones(Theta);    % default values for optional Tsize,Lsize
Lsize=1;

%---------------------------------------------------------------------
% 2. Scale Theta.  I suggest having L and Theta of order unity. 
Tsize=[..
     5.0
     0.45
     0.13
     2.0
    1.475002908821986e+01     
    6.895881052987415e-01     
    2.459122286891504e-01
     0.023
     4.0
     0.99
     1000.0
     1e-5 
     1e-5
     1000.0
     1000.0];


Lsize=1e+3;

Theta=diag(Tsize)*Theta;
%---------------------------------------------------------------------



%---------------------------------------------------------------------
% 3. Assign elements of Theta to underlying parameters (including the
%    parameter vector param).                                   

%parameters to be estimated

% ordering: d,h1,h2,theta,aa,lambda,rho0,rho1,rho2,delta,sigma,beta,
%           varm,coefm,coeft

ind=[1,1,1,0,1, 0,1,1,1,1,1,1, 1,1,1,1,1, 0,0,0,0,0, 0,0,0,0,0];
fin=find(ind);
d=Theta(1); 
dd=d;
h1=Theta(2); 
h2=Theta(3); 
%theta=Theta(4); 
aa=Theta(4);
%lambda=Theta(5);
rho0=Theta(5);     
rho1=Theta(6);    
rho2=Theta(7); 
delta=Theta(8);    
sigma=Theta(9);   
beta=Theta(10);   
varm=Theta(11:15);
%coefm=Theta(13:17);
Omega=diag(varm);
G=diag(coefm);
%coeft([1,4:5])=Theta(21:23);
trend=0;

% parameter vector
param=[ d,h1,h2,theta,aa,lambda,rho0,rho1,rho2,delta,sigma,beta,varm', ..
        coefm',coeft']';
param(fin)=Theta;
np=sum(ind);
if ind(6); lambdain=sum(ind(1:6)); else; lambdain=0; end;
%---------------------------------------------------------------------



%---------------------------------------------------------------------
% 4. Set betain = i   if the ith element of Theta corresponds to beta 
%               = 0   if beta is not to be estimated               
if ind(12); betain=sum(ind(1:12)); else betain=0; end;
%---------------------------------------------------------------------

gL(1)=h1+h2-1;
gL(2)=rho1+rho2-1;

%---------------------------------------------------------------------
% 5. Put bounds on Theta if there are values of the parameters for 
%    which L cannot be computed.
lb=[0;1e-12;1e-12;.2;.01;.99;.0001;-1;-1;0.0;.0075; .8; 1e-10;1e-10;
    1e-10; 1e-10; 1e-10;-.99;-.99;-.99;-.99;-.99;-100*ones(5,1)];
ub=[ 100;  1.0;  1.0;.65; 10;1.1;   1000;10;10;1.0;   10;1.0;50000;   .1;
      .1;10000;50000; (1-1e-12)*ones(5,1); 100*ones(5,1)];
if any(lb>param) | any(ub<param) | param(2)+param(3)>1 | param(8)+param(9)>1;
  disp('WARNING: in restricted zone with:')
  disp(Theta)
  L=1e+10; gL=[0,0]; dL=1e+10*ones(Theta);
  return
end;
%---------------------------------------------------------------------



%---------------------------------------------------------------------
% 6. Specify A,B,D,E,Sigma,ct.                                       
                                                                     
A=[1     0    0   0      0;
   rho0 rho1 rho2 0      0;
   0     1    0   0      0;
   0     0    0 (1-delta)/lambda  0;
   0     0    0   0      0];
B=[0 0 0;0 0 0;0 0 0;0 0 1/lambda;1 0 0];
D=[0 1 0 0 0;
   0 0 0 1 0;
   0 0 0 0 1;
   0 0 0 0 0;
   0 0 0 0 0;
   0 0 0 0 0];
E=[0 0 0;
   0 0 0;
   0 0 0;
   1 0 0;
   0 1 0;
   0 0 1];
Sigma=zeros(5);
Sigma(2,2)=sigma^2;
ct=1;
%---------------------------------------------------------------------



%---------------------------------------------------------------------
% 7. Specify steady states of X, u, x.                      

a1=-aa*log(1-h1);
a2=-aa*log(1-h1-h2);
lbd=(lambda-beta*(1-delta))/(beta*theta);
ti=1/theta;
tt=1-ti;
x(1,1)=rho0/(1-rho1-rho2);
x(2,1)=((h1^ti*a1^tt+h2^ti*(a2-a1)^tt)/lbd)^(-1/tt)*(1-theta)*x(1)^ ..
       (1/(1-theta))/(lbd-lambda+1-delta);
x(3,1)=(a1*(lbd-lambda+1-delta)*x(2)^(1-theta)/(x(1)*h1*(1-theta)))^(-ti);
x(4,1)=x(3);
x(5,1)=(h2*a1/(h1*(a2-a1)))^ti*x(3);
x(6,1)=(lambda-1+delta)*x(2);
X=[1 x(1) x(1) x(2) x(3)]';
u=[x(4) x(5) x(6)]';
%---------------------------------------------------------------------



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

ufn='hsu4';
%---------------------------------------------------------------------




% Compute Q,R,W and gradQ, gradR, gradW. 
% retrieve b=dr(x,theta)/x 
%          q=d^2r(x,theta)/dx^2
% plus     cub=d^3r(x,theta)/dx^3
%          du=dr(x,theta)/dtheta
%          db=db/dtheta 
%          dq=dq/dtheta   if derivatives required
if nargout > 1;
  eval(['[ux,b,q,cub,du,db,dq]=',ufn,'(x,param);'])
  tem=sum(ind(1:5));
  du=du(:); du=[du(find(ind(1:5)));zeros(np-tem,1)];
  db=[db(:,find(ind(1:5))),zeros(6,np-tem)];
  dq=[dq(:,find(ind(1:5))),zeros(36,np-tem)];
else;
  [ux,b,q]=hsu4(x,param);
end;
b=b(:);

% construct Q,R,W 
ns=length(X);
nc=length(u);
nx=length(x);
Q0=zeros(ns);
Q0(ct,ct)=ux-b'*x+x'*q*x/2;
d=(b'-x'*q)';
Q1=zeros(ns);
Q1(ct,:)=d'*D;
Q1(:,ct)=(d'*D)';
Q1=Q1/2;
Q1(ct,ct)=2*Q1(ct,ct);
Q=Q0+Q1+D'*q*D/2;
R=E'*q*E/2;
Q2=zeros(ns,nc);
Q2(ct,:)=(d'*E)/2;
W=Q2+D'*q*E/2;

% if R is invertible, construct undiscounted, no cross-products
% versions of A,B,Q,R and use double.m; otherwise, use olrp.m 
% to compute decision rules and Riccati matrix.
%if rank(R)>=length(R);
%  An=sqrt(beta)*(A-B*(R\W'));
%  Bn=sqrt(beta)*B;
%  Qn=Q-W*(R\W');
%  [Fn,P]=double(An',Bn',Qn,R);
%  F=Fn'+(R\W');
%else;
  [F,P]=dlqrtes1(beta,A,B,Q,R,W);
%end;

% compute matrix for state transition function 
Ao=A-B*F;


%---------------------------------------------------------------------
% 9. Specify C.

[tem1,tem2]=hsf4(x,param(1:4));
C=[0,0,0,1,0;
   -F(1,:);
   -1043.6*(h1*F(1,:)+h2*F(2,:));
   -F(3,:);
   tem2'*(D-E*F)+[tem1-tem2'*x,zeros(1,ns-1)]];
%---------------------------------------------------------------------



%---------------------------------------------------------------------
% 10. Observations.

%data Y(t)=[k(t),n1(t),n1(t)*ha(t),i(t),c(t)]';
T=length(zxvar);
tem=lambda.^([1:T]');

%   test:
%g4 = 1.005404130506354657.^([1:T]');
%g5 = 1.002824227733914497.^([1:T]');
g4=tem; g5=tem;


Y=[zxvar(:,1)./tem,zxvar(:,2),zxvar(:,3),zxvar(:,5)./g4,zxvar(:,4)./g5];
%k=Y(2:T,1); n1=Y(2:T,2); Th=Y(2:T,3); n1l=Y(1:T-1,2);
%c=Y(2:T,5); i=Y(2:T,4);
%n2=(Th-80*13.045*h1*n1)/(80*13.045*h2);
%tech=(c+i)./( k.^theta.*(h1*( n1-(dd/2)*(n1-n1l).^2 ).^(1-theta) +  ..
%        h2*n2.^(1-theta)) );
%lt=length(tech);
%shock=tech(3:lt)-rho0-rho1*tech(2:lt-1)-rho2*tech(1:lt-2);
%plot(tech)
%[bb,se]=autoreg(tech,2)
%save gunk bb se tech
DY=Y(2:T,:)-Y(1:T-1,:)*diag(coefm);
T=T-1;
%---------------------------------------------------------------------

if trend;
  Aoh=Ao;
  Ao=[Ao,zeros(ns,1);zeros(1,ns),1];      
  Ao(ns+1,ct)=1;
  C=[C,coeft];
  Sigma=[Sigma,zeros(ns,1);zeros(1,ns+1)];
  X=[X;0];
end;

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

X0=X;
t0=3;
%---------------------------------------------------------------------

nobs=T-t0+1;
Cbar=C*Ao-G*C;
Obar=Omega+C*Sigma*C';
[K,S]=kfilter(Ao,Cbar,Sigma,Obar,Sigma*C');
Sigmaz=Obar+Cbar*S*Cbar';
Sigmazi=inv(Sigmaz);

Xt(t0,:)=X0';
innov(t0,:)=DY(t0,:)-X0'*Cbar';
for i=t0+1:T;
  Xt(i,:)=Xt(i-1,:)*Ao'+innov(i-1,:)*K';
  innov(i,:)=DY(i,:)-Xt(i,:)*Cbar';
end;
sum1=innov(t0:T,:)'*innov(t0:T,:)/nobs;

L= nobs*(log(det(Sigmaz))+trace(Sigmazi*sum1));

if trend;
  Ao=Aoh;
end;

%-------------------------------------------------------------------
if lambdain;
  L = L + 1.5*log(lambda)*nobs*(nobs+1);
end;
%-------------------------------------------------------------------
L=L/Lsize;

if nargout > 1;


  %-------------------------------------------------------------------
  % 12. Specify gradA, gradB, gradS, gradO, gradG, gradx, gradX
 
  % user-defined
  gradA=zeros(25,np);
  gradB=zeros(15,np);
  gradS=zeros(25,np);
  gradO=zeros(25,np);
  gradG=zeros(25,np);
  grada=zeros(5,np);
  gradx=zeros(nx,np);
  tem=(h1^ti*a1^tt+h2^ti*(a2-a1)^tt)/(lambda-beta*(1-delta));
  tem1=lbd-lambda+1-delta;  
  tem2 =h1^ti*a1^tt;
  tem3=h2^ti*(a2-a1)^tt;
  tem4=-(1-delta+lbd*theta)/(beta*theta);
  tem5=h2*a1/(h1*(a2-a1));

  % derivatives with respect to d
  %i=1;

  % derivatives with respect to h1 
  i=2;
  gradx(2,i)=theta/(1-theta)*x(2)/(tem2+tem3)*(tem2/(h1*theta)+tt*tem2*aa/ ..
     (a1*(1-h1))+tt*tem3*aa*(1/(1-h1-h2)-1/(1-h1))/(a2-a1));
  gradx(3,i)=-ti*x(3)*(aa/(a1*(1-h1))-1/h1)+tt*x(3)*gradx(2,i)/x(2);
  gradx(5,i)=ti*tem5^(-tt)*x(3)*(aa*h2/(1-h1)-a1*h2/h1-a1*h2*aa* ..
     (1/(1-h1-h2)-1/(1-h1))/(a2-a1))/(h1*(a2-a1))+tem5^ti*gradx(3,i);
  gradx(6,i)=(lambda-1+delta)*gradx(2,i);

  % derivatives with respect to h2 
  i=3;
  gradx(2,i)=theta/(1-theta)*x(2)/(tem2+tem3)*(tem3/(h2*theta)+tt*tem3* ..
     aa/((1-h1-h2)*(a2-a1)) );
  gradx(3,i)=tt*x(3)*gradx(2,i)/x(2);
  gradx(5,i)=ti*tem5^(-tt)*x(3)*a1*(1-h2*aa/((a2-a1)*(1-h1-h2)))/ ..
     (h1*(a2-a1))+tem5^ti*gradx(3,i);
  gradx(6,i)=(lambda-1+delta)*gradx(2,i);

  % derivatives with respect to theta 
  % dx/dtheta, if done numerically:
  % del=eye(length(param))*1e-5;
  % gradx(:,4)=(hsgradx(param+del(:,4))-hsgradx(param-del(:,4)))/(2*del(4,4));
  if ind(4);
  gradx(2,i)=x(2)*log(beta*theta*tem*x(1))/((1-theta)*(1-theta))+ ..
     x(2)*(h1^ti*a1^tt*log(a1/h1)+h2^ti*(a2-a1)^tt*log((a2-a1)/h2))/ ..
     ((lambda-beta*(1-delta))*tem*theta*(1-theta)) +x(2)*lbd/(theta*tem1) ;
  gradx(3,i)=x(3)*( log(a1*x(2)*tem1/((1-theta)*x(1)*h1)) + theta*  ..
     (theta-1)*gradx(2,i)/x(2)+ lbd/tem1-theta/(1-theta) )/ ..
     (theta*theta);
  gradx(5,i)=-x(5)*log(h2*a1/(h1*(a2-a1)))/(theta*theta)+x(5)/x(3)* ..
     gradx(3,i);
  gradx(6,i)=(lambda-1+delta)*gradx(2,i);
  end;

  % derivatives with respect to aa 
  i=4;
  gradx(2,i)=theta/(1-theta)*x(2)*tt/aa;
  gradx(3,i)=-ti*x(3)/aa+tt*x(3)*gradx(2,i)/x(2);
  gradx(5,i)=tem5^ti*gradx(3,i);
  gradx(6,i)=(lambda-1+delta)*gradx(2,i);

  % derivatives with respect to lambda 
  if ind(6);
  gradA(ns*3+4,i)=-(1-delta)/(lambda*lambda);
  gradB(ns*2+4,i)=-1/(lambda*lambda);
  gradx(1,i)= 0; 
  gradx(2,i)=1/tt*x(2)/lbd/(beta*theta)-x(2)*(1/(beta*theta)-1)/ ..
              (lbd-lambda+1-delta);
  gradx(3,i)=-ti*x(3)/(lbd-lambda+1-delta)*(1/(beta*theta)-1)-ti*(1-theta)* ..
             x(3)/x(2)*gradx(2,i);
  gradx(5,i)=x(5)/x(3)*gradx(3,i);
  gradx(6,i)=x(2)+(lambda-1+delta)*gradx(2,i);
  end;

  % derivatives with respect to rho0 
  i=5;
  gradA(2,i)=1;
  gradx(1,i)=x(1)/rho0;
  gradx(2,i)=x(2)/((1-theta)*x(1))*gradx(1,i);
  gradx(3,i)=tt*x(3)*gradx(2,i)/x(2)+ti*x(3)*gradx(1,i)/x(1);
  gradx(5,i)=tem5^ti*gradx(3,i);
  gradx(6,i)=(lambda-1+delta)*gradx(2,i);

  % derivatives with respect to rho1 
  i=6;
  gradA(ns+2,i)=1;
  gradx(1,i)=x(1)/(1-rho1-rho2);
  gradx(2,i)=x(2)/((1-theta)*x(1))*gradx(1,i);
  gradx(3,i)=tt*x(3)*gradx(2,i)/x(2)+ti*x(3)*gradx(1,i)/x(1);
  gradx(5,i)=tem5^ti*gradx(3,i);
  gradx(6,i)=(lambda-1+delta)*gradx(2,i);

  % derivatives with respect to rho2 
  i=7;
  gradA(ns*2+2,i)=1;
  gradx(1,i)=x(1)/(1-rho1-rho2);
  gradx(2,i)=x(2)/((1-theta)*x(1))*gradx(1,i);
  gradx(3,i)=tt*x(3)*gradx(2,i)/x(2)+ti*x(3)*gradx(1,i)/x(1);
  gradx(5,i)=tem5^ti*gradx(3,i);
  gradx(6,i)=(lambda-1+delta)*gradx(2,i);

  % derivatives with respect to delta 
  i=8;
  gradA(ns*3+4,i)=-1/lambda;
  gradx(2,i)=-1/(lbd*(1-theta))*x(2)+tt*x(2)/(lbd-lambda+1-delta);
  gradx(3,i)=ti*tt*x(3)/(lbd-lambda+1-delta)+tt*x(3)*gradx(2,i)/x(2);
  gradx(5,i)=tem5^ti*gradx(3,i);
  gradx(6,i)=(lambda-1+delta)*gradx(2,i)+x(2);

  % derivatives with respect to sigma 
  i=9;
  gradS(ns+2,i)=2*sigma;

  % derivatives with respect to beta 
  i=10;
  gradx(2,i)=1/tt*x(2)*tem4/lbd-x(2)*tem4/(lbd-lambda+1-delta);
  gradx(3,i)=-ti*x(3)/(lbd-lambda+1-delta)*tem4-ti*(1-theta)* ..
    x(3)/x(2)*gradx(2,i);
  gradx(5,i)=x(5)/x(3)*gradx(3,i);
  gradx(6,i)=(lambda-1+delta)*gradx(2,i);

  % derivatives with respect to varm 
  i=11;
  gradO(1,i)=1; 
  gradO(7,i+1)=1; 
  gradO(13,i+2)=1; 
  gradO(19,i+3)=1; 
  gradO(25,i+4)=1; 

  % derivatives with respect to coefm 
  if ind(18);
  gradG(1,i)=1; 
  gradG(7,i+1)=1; 
  gradG(13,i+2)=1; 
  gradG(19,i+3)=1; 
  gradG(25,i+4)=1; 
  end;

  % derivatives with respect to coeft 
  %grada(1,i)=1;
  %grada(2,i+1)=1;
  %grada(3,i+2)=1;
  %grada(4,i+3)=1;
  %grada(5,i+4)=1;

  % derivatives of x(4):
  gradx(4,:)=gradx(3,:);
  
  gradX=zeros(5,np);
  gradX(2:5,:)=gradx([1,1:3],:);
  %-------------------------------------------------------------------


  % construct gradb,gradq
  gradb=q*gradx+db;
  gradq=cub*gradx+dq;

  % construct gradQ, etc.
  tem1=zeros(ns,nc);
  tem1(:)=[1:ns*nc]';
  tem2=tem1'; 
  tem2=tem2(:);
  tem1=tem1(:);
  Knk=zeros(ns,nc);
  for i=1:ns*nc;
    Knk(tem1(i),tem2(i))=1;
  end;
  tem1=zeros(ns);
  tem1(:)=[1:ns*ns]';
  tem2=tem1'; 
  tem2=tem2(:);
  tem1=tem1(:);
  Knn=zeros(ns*ns);
  for i=1:ns*ns;
    Knn(tem1(i),tem2(i))=1;
  end;
  tem1=zeros(nc,ns);
  tem1(:)=[1:ns*nc]';
  tem2=tem1'; 
  tem2=tem2(:);
  tem1=tem1(:);
  Kkn=zeros(nc,ns);
  for i=1:ns*nc;
    Kkn(tem1(i),tem2(i))=1;
  end;
  tem=(ns+1)*(ct-1);
  gradQ=[zeros(tem,np);
         du'-x'*gradb+x'*q*gradx+.5*kron(x',x')*gradq;
         zeros((ns*ns)-tem-1,np)];
  tem= [zeros(ns*(ct-1),np);
        .5*D'*gradb-.5*D'*q*gradx-.5*kron(x',D')*gradq;
        zeros(ns*(ns-ct),np)];
  gradQ=gradQ + tem + Knn*tem;
  gradQ=gradQ + .5*kron(D',D')*gradq;
  gradR=kron(E',E')*gradq/2;
  gradWp=[zeros(nc*(ct-1),np);
          (E'*gradb-E'*q*gradx-kron(x',E')*gradq)/2;
          zeros((ns-ct)*nc,np)] +kron(D',E')*gradq/2; 
  gradW=Kkn*gradWp;

  ZO=inv(R+beta*B'*P*B);
  ZB=beta*kron(Ao',ZO*B')/(eye(ns*ns)-beta*kron(Ao',Ao'));  % rcond errors?
  Z1=beta*kron(eye(ns),ZO*B'*P)+ZB*(beta*kron(eye(ns),Ao'*P)+ ..
      beta*kron(Ao'*P,eye(ns))*Knn); 
  Z2=-beta*kron(F',ZO*B'*P)+beta*kron(Ao'*P,ZO)*Knk-ZB*(beta* ..
      kron(F',Ao'*P)+beta*kron(Ao'*P,F')*Knk);
  Z3=ZB;
  Z4=-kron(F',ZO)+ZB*kron(F',F');
  Z5=kron(eye(ns),ZO)*Knk-ZB*(kron(F',eye(ns))+kron(eye(ns),F')*Knk);
  Z6=zeros(ns*nc,np);
  if betain; 
    Z6(:,betain)=vec(ZO*B'*P*Ao)+ZB*vec(Ao'*P*Ao); 
  end;
  gradF=Z1*gradA+ Z2*gradB+ Z3*gradQ+ Z4*gradR +Z5*gradW+ Z6;
  gradAo=gradA-kron(F',eye(ns))*gradB-kron(eye(ns),B)*gradF;

  if trend;
    tem=zeros(Ao);
    for i=1:np;
      tem(:)=gradAo(:,i);
      tem1=[tem,zeros(ns,1);zeros(1,ns+1)];
      gradAn(:,i) = tem1(:);
      tem(:)=gradS(:,i);
      tem1=[tem,zeros(ns,1);zeros(1,ns+1)];
      gradSn(:,i) = tem1(:);
    end;
    gradAo=gradAn;
    gradS=gradSn;
    Ao=[Ao,zeros(ns,1);zeros(1,ns),1];      
    Ao(ns+1,ct)=1;
  end;

  %-------------------------------------------------------------------
  % gradC
  gradFp=Kkn*gradF;
  gradCp=zeros(25,np);
  gradCp(6:5+ns*nc,:)=-gradFp;
  gradCp(11:15,:)=1043.6*(gradCp(6:10,:)*h1+gradCp(11:15,:)*h2);
  if ind(2);
     tem=sum(ind(1:2));
     gradCp(11:15,tem)=gradCp(11:15,tem)-1043.6*F(1,:)';
  end;
  if ind(3);
     tem=sum(ind(1:3));
     gradCp(11:15,tem)=gradCp(11:15,tem)-1043.6*F(2,:)';
  end;
  tem=sum(ind(1:4));
  [f,fj,fh,f,df,dfj]=hsf4(x,param(1:4));
  dfj=[dfj(:,find(ind(1:4))),zeros(6,np-tem)];
  df=df(:); df=[df(find(ind(1:4)));zeros(np-tem,1)];
  gradCp(21:25,:)=[df'-x'*(fh*gradx+dfj);zeros(4,np)]+ ..
                  (D-E*F)'*(fh*gradx+dfj)-kron(fj'*E,eye(ns))*gradFp;
  gradC=Knn*gradCp;

  if trend;
    gradC=[gradC;grada];
  end;
  %-------------------------------------------------------------------

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

  sum2=Xt(t0:T,:)'*innov(t0:T,:)/nobs;
  sum3=Xt(t0:T-1,:)'*lamb(t0+1:T,:)/nobs;
  sum4=DY(t0:T-1,:)'*lamb(t0+1:T,:)/nobs;
  sum5=innov(t0:T-1,:)'*lamb(t0+1:T,:)/nobs;
  sum6=Y(t0:T,:)'*innov(t0:T,:)/nobs;
  sum7=Y(t0:T-1,:)'*lamb(t0+1:T,:)/nobs;

  %-------------------------------------------------------------------
  if lambdain;
    sum8=-( ([1;0;0;1;1]*[t0+1:T+1]).*Y(t0+1:T+1,:)' )*innov(t0:T,:)/..
          (lambda*nobs);
    sum9=-( ([1;0;0;1;1]*[t0:T]).*Y(t0:T,:)' )*innov(t0:T,:)/ ..
          (lambda*nobs);
    sum10=-( ([1;0;0;1;1]*[t0:T-1]).*Y(t0:T-1,:)' )*lamb(t0+1:T,:)/ ..
          (lambda*nobs);
    sum11=-( ([1;0;0;1;1]*[t0+1:T]).*Y(t0+1:T,:)' )*lamb(t0+1:T,:)/ ..
          (lambda*nobs);
  end;

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


  M1=Sigmazi-Sigmazi*sum1*Sigmazi;
  M2=Cbar'*M1*Cbar-2*SC'*sum5*AA;
  M2=.5*(M2+M2');
  nct=ones(ns,1); nct(ct)=0;
  M3=zeros(ns);
  M3(nct,nct)=edlyap(AA(nct,nct)',M2(nct,nct));
  M=M1+K'*M3*K;
  D1=2*( (-sum3-S*SC'*sum5)*(eye(S)-K*C)+S*Cbar'*M*C- ..
     (sum2+S*AA'*sum5')*Sigmazi*C+S*AA'*M3-S*Ao'*M3*K*C );
  D2=2*( Sigma*C'*M+(Ao*S*Cbar'+Sigma*C')*Sigmazi*sum5*K-Ao*S*AA'*(sum5' ..
     *Sigmazi+M3*K)+S*AA'*sum5'*Sigmazi*G-Sigma*(eye(S)-C'*K')*sum5'*Sigmazi);
  D2=D2+2*(-S*Cbar'*(M1+Sigmazi*sum5*K)*G+Ao*S*Cbar'*M1-Ao*sum2* ..
     Sigmazi+sum2*Sigmazi*G-sum3*K*G+Ao*sum3*K-Sigma*M3*K+S*AA'*M3*K*G ); 
  D3=2*( .5*C'*M*C-C'*Sigmazi*sum5*(eye(S)-K*C)+.5*M3-C'*K'*M3 );
  D4=2*Sigmazi*sum5*K+M;
  D5=2*(-C*S*Cbar'*(M+Sigmazi*sum5*K)-(sum6-C*sum2)*Sigmazi-C*sum3*K+ ..
     C*S*(AA'*sum5'*Sigmazi+Ao'*M3*K)+sum7*K );
  tem1=zeros(Ao);
  tem2=zeros(C);
  tem3=zeros(Sigma);
  tem4=zeros(Omega);
  tem5=zeros(G);
  for j=1:np;
    tem1(:)=gradAo(:,j);
    tem2(:)=gradC(:,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(gradX(:,j)*lamb(t0,:)/nobs);
  end; 
  %-------------------------------------------------------------------
  if lambdain;
    dL(lambdain)=dL(lambdain)-2*trace(K*sum11-K*G*sum10)+ ..
            2*trace(sum8*Sigmazi-G*sum9*Sigmazi) +1.5*(nobs+1)/lambda;
  end; 
  %-------------------------------------------------------------------
  dL=nobs*diag(Tsize)*dL/Lsize;
end;
