function [L,dL]=mle(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=3.08; 
h1=.432; 
h2=.0639; 
theta=.365;
aa=3.62; 
delta=.0241; 
lambda=1.005;    
rho0=5.20; 
rho1=.7; 
rho2=.277;
sigma=1.64;
beta=.995; 
varm=[1000;.001;.001;1000;1000];
coefm=zeros(5,1);

% parameter vector
param=[d h1 h2 theta aa lambda rho0 rho1 rho2 delta sigma beta]';
param=[param;varm;coefm];
%---------------------------------------------------------------------

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

%---------------------------------------------------------------------
% 2. Scale Theta.  I suggest having L and Theta of order unity. 
Tsize=[ ..
     3.717135508567333e+00
     1.816439845465921e-01
     3.592631231510675e-01
     3.809023987557150e-01
     5.087117940771078e+00
     1.005250731575861e-01
     9.701359011931011e-01
     2.942454377737908e-02
     2.354265203786171e-02
     2.557150144332708e+00
     9.925760857991500e-01
     8.029968809863021e+01
     7.274133929082318e-05
     1.506317232755622e-04
     1.549631527499441e+03
     3.318470548594628e+03
     1.419938856306720e-06
     8.501089002209466e-06
     1.391246931936938e-05
     1.198196067843829e-05
     1.407381554935088e-05];

Lsize=1100;

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



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

%parameters to be estimated
ind=[1:5,7:22];
d=Theta(1); h1=Theta(2); h2=Theta(3); theta=Theta(4); aa=Theta(5);
rho0=Theta(6);     rho1=Theta(7);    rho2=Theta(8); 
delta=Theta(9);    sigma=Theta(10);   beta=Theta(11);   
varm=Theta(12:16);
coefm=Theta(17:21);
Omega=diag(varm);
G=diag(coefm);

param(ind)=Theta;
param(6)=lambda;
np=length(ind);
%---------------------------------------------------------------------



%---------------------------------------------------------------------
% 4. Set betain = i   if the ith element of Theta corresponds to beta 
%               = 0   if beta is not to be estimated               
betain=11;
%---------------------------------------------------------------------

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=[-10;1e-8;1e-8;.2;.01; .99;.0001; .6;-.5;.02;.0075; .96;    1;1e-8;
    1e-8;    1;    1;-.99;-.99;-.99;-.99;-.99];
ub=[ 10;  .8;  .5;.8; 10;1.0065;   10;1.1; .3;.03;   10;.999;10000;  .1;
      .1;10000;10000; .99; .99; .99; .99; .99];
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);'])
  du=[du(:);zeros(np-5,1)];
  db=[db,zeros(6,np-5)];
  dq=[dq,zeros(36,np-5)];
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;

% make sure R is invertible and construct undiscounted, no cross-products
% versions of A,B,Q,R 
if rank(R)>=length(R);
  An=sqrt(beta)*(A-B*(R\W'));
  Bn=sqrt(beta)*B;
  Qn=Q-W*(R\W');
else;
  disp('WARNING: R is not invertible at Theta equal:')
  Theta
  return
end;

% compute Riccati matrix and decision rules 
[Fn,P]=double(An',Bn',Qn,R);
Fn=Fn';
Aon=An-Bn*Fn;

% construct matrices for original problem
F=Fn+(R\W');
Ao=A-B*F;


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

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



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

%data Y(t)=[k(t),n1(t),n2(t),i(t),c(t)]';
T=length(zxvar);
tem=lambda.^([1:T]');
Y=[zxvar(:,1)./tem,zxvar(:,2),zxvar(:,3),zxvar(:,5)./tem,zxvar(:,4)./tem];
DY=Y(2:T,:)-Y(1:T-1,:)*diag(coefm);
T=T-1;
%---------------------------------------------------------------------



%---------------------------------------------------------------------
% 11. Specify initial state, X0.

X0=X;
%---------------------------------------------------------------------



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);
t0=3;
Xt(t0,:)=X0';
innov(t0,:)=DY(t0,:)-X0'*Cbar';
sum1=innov(t0,:)'*innov(t0,:);
sum2=Xt(t0,:)'*innov(t0,:);
sum6=Y(t0,:)'*innov(t0,:);
for i=t0+1:T;
  Xt(i,:)=Xt(i-1,:)*Ao'+innov(i-1,:)*K';
  innov(i,:)=DY(i,:)-Xt(i,:)*Cbar';
  sum1=sum1+innov(i,:)'*innov(i,:);
  sum2=sum2+Xt(i,:)'*innov(i,:);
  sum6=sum6+Y(i,:)'*innov(i,:);
end;
sum1=sum1/(T-t0+1);
sum2=sum2/(T-t0+1);
sum6=sum6/(T-t0+1);


L= (T-t0+1)*(log(det(Sigmaz))+trace(Sigmazi*sum1))+  ..
   1.5*log(lambda)*(T-t0+1)*(T-t0+2);
L=L/Lsize;



if nargout > 1;


  %-------------------------------------------------------------------
  % 12. Specify gradA, gradB, gradS, gradO, gradG, gradx, gradX
 
  % user-defined
  gradA=zeros(25,np);
  gradA(2,6)=1;
  gradA(7,7)=1;
  gradA(12,8)=1;
  %gradA(19,6)=-(1-delta)/(lambda*lambda);
  gradA(19,9)=-1/lambda;
  gradB=zeros(15,np);
  %gradB(14,6)=-1/(lambda*lambda);
  gradS=zeros(25,np);
  gradS(7,10)=2*sigma;
  gradO=zeros(25,np);
  gradO(1,12)=1; 
  gradO(7,13)=1; 
  gradO(13,14)=1; 
  gradO(19,15)=1; 
  gradO(25,16)=1; 
  gradG(1,17)=1; 
  gradG(7,18)=1; 
  gradG(13,19)=1; 
  gradG(19,20)=1; 
  gradG(25,21)=1; 
  gradx=zeros(6,np);

  % dx/dtheta 
  % if done numerically:
   del=eye(length(param))*1e-10;
  % gradx(:,4)=(hsgradx(param+del(:,4))-hsgradx(param-del(:,4)))/(2*del(4,4));
  fp= hsgradx(param+del(:,8))
  fm=hsgradx(param-del(:,8))
  (fp-fm)/(2*del(8,8))
  fp= hsgradx(param+del(:,9))
  fm=hsgradx(param-del(:,9))
  (fp-fm)/(2*del(9,9))

  tem=(h1^ti*a1^tt+h2^ti*(a2-a1)^tt)/(lambda-beta*(1-delta));
  tem1=lbd-lambda+1-delta;  
  
  gradx(2,4)=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,4)=x(3)*( log(a1*x(2)*tem1/((1-theta)*x(1)*h1)) + theta*  ..
     (theta-1)*gradx(2,4)/x(2)+ lbd/tem1-theta/(1-theta) )/ ..
     (theta*theta);
  gradx(5,4)=-x(5)*log(h2*a1/(h1*(a2-a1)))/(theta*theta)+x(5)/x(3)* ..
     gradx(3,4);
  gradx(6,4)=(lambda-1+delta)*gradx(2,4);

  % dx/dlambda:
  % 1: 0, 
  % 2: 1/tt*x(2)/lbd/(beta*theta)-x(2)*(1/(beta*theta)-1)/(lbd-lambda+1-delta);
  % 3: -ti*x(3)/(lbd-lambda+1-delta)*(1/(beta*theta)-1)-ti*(1-theta)* ..
  %    x(3)/x(2)*gradx(2,6);
  % 4: dx(3)/dlambda
  % 5: x(5)/x(3)*gradx(3,--);
  % 6: x(2)+(lambda-1+delta)*gradx(2,--);

  % derivatives of x(1):
  gradx(1,[6:8])=[x(1)/rho0,[1,1]*x(1)/(1-rho1-rho2)];
  
  % derivatives of x(2):
  tem =h1^ti*a1^tt;
  tem1=h2^ti*(a2-a1)^tt;
  tem2=-(1-delta+lbd*theta)/(beta*theta);
  gradx(2,2)=theta/(1-theta)*x(2)/(tem+tem1)*(tem/(h1*theta)+tt*tem*aa/ ..
     (a1*(1-h1))+tt*tem1*aa*(1/(1-h1-h2)-1/(1-h1))/(a2-a1));
  gradx(2,3)=theta/(1-theta)*x(2)/(tem+tem1)*(tem1/(h2*theta)+tt*tem1* ..
     aa/((1-h1-h2)*(a2-a1)) );
  gradx(2,5)=theta/(1-theta)*x(2)*tt/aa;
  gradx(2,[6:8])=x(2)/((1-theta)*x(1))*gradx(1,[6:8]);
  gradx(2,9)=-1/(lbd*(1-theta))*x(2)+tt*x(2)/(lbd-lambda+1-delta);
  gradx(2,11)=1/tt*x(2)*tem2/lbd-x(2)*tem2/(lbd-lambda+1-delta);

  % derivatives of x(3):
 
  gradx(3,2)=-ti*x(3)*(aa/(a1*(1-h1))-1/h1)+tt*x(3)*gradx(2,2)/x(2);
  gradx(3,3)=tt*x(3)*gradx(2,3)/x(2);
  gradx(3,5)=-ti*x(3)/aa+tt*x(3)*gradx(2,5)/x(2);
  gradx(3,[6:8])=tt*x(3)*gradx(2,[6:8])/x(2)+ti*x(3)*gradx(1,[6:8])/x(1);
  gradx(3,9)=ti*tt*x(3)/(lbd-lambda+1-delta)+tt*x(3)*gradx(2,9)/x(2);
  gradx(3,11)=-ti*x(3)/(lbd-lambda+1-delta)*tem2-ti*(1-theta)* ..
    x(3)/x(2)*gradx(2,11);
 
  % derivatives of x(4):
  gradx(4,:)=gradx(3,:);
 
  % derivatives of x(5):
  tem=h2*a1/(h1*(a2-a1));
 
  gradx(5,2)=ti*tem^(-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))+tem^ti*gradx(3,2);
  gradx(5,3)=ti*tem^(-tt)*x(3)*a1*(1-h2*aa/((a2-a1)*(1-h1-h2)))/ ..
     (h1*(a2-a1))+tem^ti*gradx(3,3);
  gradx(5,[5:9])=tem^ti*gradx(3,[5:9]);
  gradx(5,11)=x(5)/x(3)*gradx(3,11);
 
  % derivatives of x(6):
  gradx(6,2)=(lambda-1+delta);
  gradx(6,[2,3,5:9,11])=(lambda-1+delta)*gradx(2,[2,3,5:9,11]);
  gradx(6,9)=gradx(6,9)+x(2);

gradx(:,7:8)

  gradX=zeros(5,np);
  gradX(2:5,:)=gradx([1,1:3],:);

  V=eye(5); V(2,2)=0; V(3,3)=0;
  sum8=(t0+1)*V*Y(t0+1,:)'*innov(t0,:);
  sum9=t0*V*Y(t0,:)'*innov(t0,:);
  for i=t0+1:T;
    sum8=sum8+(i+1)*V*Y(i+1,:)'*innov(i,:);
    sum9=sum9+i*V*Y(i,:)'*innov(i,:);
  end; 
  sum8=-sum8/(lambda*(T-t0+1));
  sum9=-sum9/(lambda*(T-t0+1));
  %-------------------------------------------------------------------





  % 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')*inv(eye(ns*ns)-beta*kron(Ao',Ao'));
  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;

  % user-defined 
  gradFp=Kkn*gradF;
  gradCp=zeros(25,np);
  gradCp(6:5+ns*nc,:)=-gradFp;
  [f,fj,fh,f,df,dfj]=hsf4(x,param(1:4));
  dfj=[dfj,zeros(6,np-4)];
  df=[df(:);zeros(np-4,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;

  % user-defined
  % sum10,sum11 below

  AA=Ao-K*Cbar;
  lamb(T,:)=innov(T,:)*Sigmazi*Cbar;
  sum3=Xt(T-1,:)'*lamb(T,:);
  sum4=DY(T-1,:)'*lamb(T,:);
  sum5=innov(T-1,:)'*lamb(T,:);
  sum7=Y(T-1,:)'*lamb(T,:);
  sum10=(T-1)*V*Y(T-1,:)'*lamb(T,:);
  sum11=T*V*Y(T,:)'*lamb(T,:);
  for i=T-1:-1:t0+1;
    lamb(i,:)=lamb(i+1,:)*AA+innov(i,:)*Sigmazi*Cbar;  
    sum3=sum3+Xt(i-1,:)'*lamb(i,:);
    sum4=sum4+DY(i-1,:)'*lamb(i,:);
    sum5=sum5+innov(i-1,:)'*lamb(i,:);
    sum7=sum7+Y(i-1,:)'*lamb(i,:);
    sum10=sum10+(i-1)*V*Y(i-1,:)'*lamb(i,:);
    sum11=sum11+i*V*Y(i,:)'*lamb(i,:);
  end;
  lamb(t0,:)=lamb(t0+1,:)*AA+innov(t0,:)*Sigmazi*Cbar;
  sum3=sum3/(T-t0+1);
  sum4=sum4/(T-t0+1);
  sum5=sum5/(T-t0+1);
  sum7=sum7/(T-t0+1);
  sum10=-sum10/(lambda*(T-t0+1));
  sum11=-sum11/(lambda*(T-t0+1));

  M1=Sigmazi-Sigmazi*sum1*Sigmazi;
  M2=Cbar'*M1*Cbar-2*Cbar'*Sigmazi*sum5*AA;
  M3=dlyap(AA',.5*(M2+M2'));
  tem1=zeros(Ao);
  tem2=zeros(C);
  tem3=zeros(Omega);
  for j=1:np;
    tem1(:)=gradAo(:,j);
    d1(j,1)=2*trace(tem1*((-sum3-S*Cbar'*Sigmazi*sum5)*(eye(S)-K*C) .. 
      +S*Cbar'*M1*C-sum2*Sigmazi*C-S*AA'*sum5'*Sigmazi*C+S*Cbar'*K'* ..
      M3*K*C+S*AA'*M3-S*Ao'*M3*K*C));
    tem2(:)=gradC(:,j);
    d2(j,1)=2*trace(tem2*(Ao*S*Cbar'*M1-S*Cbar'*M1*G+Sigma*C'*M1- ..
      Ao*sum2*Sigmazi+sum2*Sigmazi*G-sum3*K*G+Ao*sum3*K-Sigma* ..
      (eye(S)-C'*K')*sum5'*Sigmazi));
    d2(j,1)=d2(j,1)+2*trace(tem2*(Sigma*C'*Sigmazi*sum5*K-Ao*S*AA'*sum5' ..
      *Sigmazi+Ao*S*AA'*sum5'*Sigmazi*G+Ao*S*Cbar'*Sigmazi*sum5*K-S* ..
      Cbar'*Sigmazi*sum5*K*G-Sigma*M3*K-Ao*S*AA'*M3*K+S*AA'*M3*K*G)); 
    d2(j,1)=d2(j,1)+2*trace(tem2*(Sigma*C'*K'*M3*K));
    tem1(:)=gradS(:,j);  
    d3(j,1)=2*trace(tem1*(.5*C'*M1*C-C'*Sigmazi*sum5*(eye(S)-K*C)+.5*M3- ..
      C'*K'*M3+.5*C'*K'*M3*K*C));
    tem3(:)=gradO(:,j);
    d4(j,1)=trace(tem3*(2*Sigmazi*sum5*K+M1+K'*M3*K));
    d5(j,1)=-2*trace(gradX(:,j)*lamb(t0,:)/(T-t0+1));
    tem1(:)=gradG(:,j);
    d6(j,1)=2*trace(tem1*(-C*S*Cbar'*M1-(sum6-C*sum2)*Sigmazi-C*sum3* ..
      K+C*S*AA'*sum5'*Sigmazi-C*S*Cbar'*Sigmazi*sum5*K+sum7*K+C*S*Ao'* ..
      M3*K-C*S*Cbar'*K'*M3*K));
  end; 
  dL=d1+d2+d3+d4+d5+d6;
%  dL(6)=dL(6)-2*trace(K*sum11-K*G*sum10)+2*trace(sum8*Sigmazi- ..
%    G*sum9*Sigmazi)  +1.5*(T-t0+2)/lambda;
  dL=diag(Tsize)*(T-t0+1)*dL/Lsize;
end;

