function [g,dg] = gmmh(Theta);
%GMMH  A user-supplied function for GMM2 estimation.
%      [g,dg]=gmmh(theta) is called by GMM2 and GMMJ2. 
%      g is the m x 1 (or mgmm x 1) vector of orthogonality
%      conditions.  dg is m x p where p is the number of parameters.
%

%      Ellen McGrattan, 11-13-92

%      Code for Hansen/Sargent LQ economies
%
%      Model: max  E[sum beta^t x[t]'Q x[t]+u[t]'R u[t]+2x[t]'W u[t]) | x[0]]
%            {u[t]} t=0:inf 
%
%                  subject to  x[t+1]=A*x[t]+B*u[t]+C*e[t+1],  e ~ N(0,I)
%                              x[0]=x0 given
%
%      where  [Q  W]  = -.5*(G'*G+H'*H) 
%             [W' R]
%
%             G = ug*[phic phig]^{-1}[0, gamma, ud, -phii] 
%
%             H = [lambda, phih*uc*[phic phig]^{-1}[gamma, ud, -phii] ]
%                 - [0, 0, ub, 0]
%
%             A = [deltah thetah*uc*[phic phig]^{-1}*[gamma, ud]  ]
%                 [  0           deltak                       0   ]
%                 [  0             0                         a22  ]
%
%             B = [ -thetah*uc*[phic phig]^{-1} ]
%                 [           thetak            ]  
%                 [             0               ]
%
%             C = [ 0  ]
%                 [ 0  ]
%                 [ c2 ]    
%
%             x(t) = [h(t-1)',k(t-1)',z(t)']'  
%             u(t) = i(t) 
%
%             h(t) = household capital stocks  (nh x 1)
%             k(t) = market capital stocks  (nk x 1) 
%             z(t) = taste (b(t)) and technology shocks (d(t))  (nz x 1)
%             i(t) = investment goods   (ni x 1)
%
%      Optimal Controls:
%             u[t] = -F x[t],   F = inv(R+beta*B'*P*B)*(beta*B'*P*A+W')
%                               P = Q+beta*A'*P*A-(beta*A'*P*B+W)*F 
%
%      Transition function for States:
%             x[t+1] = (A-B*F) x[t] + C*e[t+1]
%                    =  Ao x[t] + C*e[t+1]          E e[t]e[t]' = Sigma
%        
%      Observations:   
%
%             z1[t]  =  D*x[t] + v[t]
%             v[t+1] =  Gam*v[t] + nu[t+1],         E nu[t]nu[t]' = Omega   
%         
%             (z1[t], t=1,...,T  in xgmm)
%
%      Score Vector:
%
%             sc[t]    = d l[t]/ d theta, where l[t] is contribution
%                                         at t to the likelihood function, L
%                      
%             L        = T*( ln(|Sigmaz| + trace(inv(Sigmaz)*Sum) ) 
%
%             Sigmaz   = E[z1[t]-z1hat[t]][z1[t]-z1hat[t]]'
%             z1hat[t] = E[z1[t]|z1[t-1],z1[t-2],...,z1[1],x[1]]
%             Sum      = sum(innov(t)*innov(t)')/T
%                        t=1:T
%
%      Reference:   Hansen, Lars P. and Thomas J. Sargent, "Recursive 
%                   Linear Models of Dynamic Economies," manuscript.

Theta=Theta(:);

%---------------------------------------------------------------------
% 1. Provide default values of the parameters underlying the following
%    matrices:     deltak (nk x nk),thetak (nk x ni), phic (nd x nc),
%                  phig (nd x ng), phii (nd x ni), gamma (nd x nk),
%                  deltah (nh x nh), lambda (nb x nh), thetah (nh x nc),
%                  pih (nb x nc), a22 (nz x nz), c2 (nz x nw), x0 (nx x 1)
%    and values for the selector matrices:    ub (nb x nz), ud (nd x nz)
%    and a growth rate: grate
%    and a coefficient matrix for measurement error process: Gam
%    and a covariance matrix for measurement error process:  Omega

% default values for parameters
beta   = 1/1.05;
deltah = .9;
deltak = .95;
gam    = .12;
lambda = -.1;
phi1   = -1;
pih    = 1;
thetah = (1-deltah);
thetak = 1;
rho1   = .8;
rho2   = .5;
sigma1 = 1;
sigma2 = 0;
sigma3 = 1;
ub1    = 30;
ud11   = 5;
ud12   = 1;
varm   = .001*ones(4,1);
x0     = [10.3951 77.0725 1 0 0]';
grate  = 1;
Gam    = zeros(4);
Omega  = zeros(4);

%    Note: make sure the user provided derivatives are consistent with
%    the defined matrices here.  
%---------------------------------------------------------------------


%---------------------------------------------------------------------
% 3. Assign elements of Theta to underlying parameters 

% ordering: beta,deltah,deltak,gam,lambda,phi1,pih,thetah,thetak,rho1,
%           rho2,sigma1,sigma2,sigma3,ub1,ud11,ud12,grate,varm
%            1     2      3     4    5     6    7    8      9     10
%            11    12     13     14   15   16   17   18    22   

%ind=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1]; 
ind=[1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
beta   = Theta(1);
deltah = Theta(2);
deltak = Theta(3);
gam    = Theta(4);
lambda = Theta(5);
phi1   = Theta(6);
pih    = Theta(7);
thetah = Theta(8);
%thetak= Theta(9);
%rho1  = Theta(10);
%rho2  = Theta(11);
%sigma1= Theta(12);
%sigma2= Theta(13);
%sigma3= Theta(14);
%ub1   = Theta(15);
%ud11  = Theta(16);
%ud12  = Theta(17);
grate  = 1; 
%varm  = Theta(18:21);

gamma  = [gam 0]';
phic   = [1 0]';
phig   = [0 1]';
phii   = [1 -phi1]';
ud     = [ud11 ud12 0; 0 0 0];
a22    = [1 0 0; 0 rho1 0; 0 0 rho2];
c2     = [0 0; 
          sigma1 0;
          sigma2 sigma3];
ub     = [ub1 0 0];
Omega  = diag(varm);
no     = length(varm);

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


%---------------------------------------------------------------------
% 4. 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 any other special parameters.

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

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

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

% ordering: beta,deltah,deltak,gam,lambda,phi1,pih,thetah,thetak,rho1,
%           rho2,sigma1,sigma2,sigma3,ub1,ud11,ud12,grate,varm
Ub = [.9999;.999;.999;.2;.3;2;2;.5;2;.9999;.9999;3;3;3;100;20;5;1.1; ..
      1000*ones(4,1)];
Ub = Ub(ind);
Lb = [-.8;-.8;0;.01;-.5;-3;-2;.01;.01;-.8;-.8;0;0;0;0;0;0;.99; ..
      zeros(4,1)];
Lb = Lb(ind);
if any(Theta>Ub) | any(Theta<Lb);
  g  = -99;
  dg = -99;
  return
end;


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

[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');
%g = q;
%k0= a'*g*b/(b'*g*b+R);
%dd=1;
%it=1;
%maxit=1000;
%tol=1e-10;                         % reset tol if necessary
%while (dd>tol & it<=maxit);
%  v1=(v+b0*g);
%  vv=v1\a0;
%  a1=a0*vv;
%  b0=b0+a0*(v1\(b0*a0'));
%  g=g+a0'*g*vv;
%  k1=a'*g*b/(b'*g*b+R);
%  dd=max(max(abs( (k1-k0)./(1+abs(k0)) )));
%  it=it+1;
%  a0=a1;
%  k0=k1;
%end
%F = k1'+R\W';
%P = g;

%=====================
%% 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
chi = v(1:nx,index(1:nx));
lam = v((nx+1):(2*nx),index(1:nx));
P = real(lam/chi);
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,I];
%[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
%chi = v(1:nx,index(1:nx));
%lam = v((nx+1):(2*nx),index(1:nx));
%P = real(lam/chi);
%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. (Use global statement instead of loading.)

%data xgmm(t)=[k(t-1),i(t),c(t),l(t),  asset prices(t)]';

T   = length(xgmm);
tem = grate.^([1:T]');    % assume that k,i,c on balanced growth path

y=[xgmm(:,1)./tem,xgmm(:,2)./tem,xgmm(:,3)./tem,xgmm(:,4)];
dy=y(2:T,:)-y(1:T-1,:)*Gam;
T=T-1;
%---------------------------------------------------------------------


%---------------------------------------------------------------------
% 9. Specify D.   

skl = [zeros(nk,nh),eye(nk),zeros(nk,nz)];
sc  = phiinc*(gamma*skl+[zeros(nd,nh+nk),ud]+phii*F);
sg  = phiing*(gamma*skl+[zeros(nd,nh+nk),ud]+phii*F);
sd  = [zeros(nd,nh+nk),ud];
pay = sd(1,:);
D   = [0,1,0,0,0;-F;sc;sg];
%---------------------------------------------------------------------

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);

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));

%-----------------------------------------------------------------------
% 10. Set jacobian resulting from observations being a function of Theta

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


%-------------------------------------------------------------------
% 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

p      = length(ind);
grada2 = zeros(nz*nz,p);    grada2(5,10)=1; grada2(9,11)=1;  
gradc2 = zeros(nz*nw,p);    gradc2(2,12)=1; gradc2(3,13)=1; gradc2(6,14)=1;  
graddh = zeros(nh*nh,p);    graddh(1,2)=1;
graddk = zeros(nk*nk,p);    graddk(1,3)=1;
gradg  = zeros(nd*nk,p);    gradg(1,4)=1;
gradl  = zeros(nb*nh,p);    gradl(1,5)=1;
gradpc = zeros(nd*nc,p);
gradpg = zeros(nd*ng,p);
gradpi = zeros(nd*ni,p);    gradpi(2,6)=-1;
gradph = zeros(nb*nc,p);    gradph(1,7)=1;
gradth = zeros(nh*nc,p);    gradth(1,8)=1;
gradtk = zeros(nk*ni,p);    gradtk(1,9)=1;
gradub = zeros(nb*nz,p);    gradub(1,15)=1;
gradud = zeros(nd*nz,p);    gradud(1,16)=1; gradud(3,17)=1;
gradx0 = zeros(nx,p);
gradG  = zeros(16,p);
gradO  = zeros(16,p);       gradO(1,19)=1; gradO(6,20)=1; gradO(11,21)=1;
                            gradO(16,22)=1;
%-------------------------------------------------------------------

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;

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(nc),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 gradD. 

tem1=-kron(D(3:4,:)',phiin)*gradp+kron([0;1;0;0;0],phiin)*gradg+ ..
      kron([zeros(3,2),eye(3)]',eye(nd))*gradud+ ..  
      kron(F',phiin)*gradpi+kron(eye(nx),phiin*phii)*gradF;
gradDp=[zeros(5,p);-Kix*gradF;Kdx*tem1];
%-------------------------------------------------------------------

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)=edlyap(AA(nct,nct)',M2(nct,nct));
M=M1+K'*M3*K;
D1=2*( (-sum3-S*SD'*sum4)*(eye(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(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(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(Ao);
tem2=zeros(D');
tem3=zeros(Sigma);
tem4=zeros(Omega);
tem5=zeros(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 Theta.

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; 
%-------------------------------------------------------------------

g = dL;

% tack on moment conditions for prices:
% for now, assume q(x) = x  (see letter to Tom)

[nas,tem] = size(pay);
if nas > 0;
  nt        = length(pshocks);  % pshocks is global 
  e1        = eye(nc);
  e1        = e1(1,:);
  mc        = 2*beta*thetah'*[eye(nh),zeros(nh,nk)]*P(1:nh+nk,:)*Ao+ ...
               pih'*([zeros(nb,nh+nk),ub]-lambda* ...
               [eye(nh),zeros(nh,nk+nz)]-pih*sc);
  xs        = ltitr(Ao,c,pshocks,x0);
  v1        = edlyap(sqrt(beta)*Ao,Sigma)*beta/(1-beta);
  nn        = floor(nt/T);
  tem1      = zeros(T,nn);
  tem2      = zeros(T*nas,nn);
  for i=1:nas;
    v2   = doublej(sqrt(beta)*Ao',pay(i,:)'*e1*mc);
    for j=1:nt;
      apr(j,i) = xs(j,:)*v2*xs(j,:)';
    end;
    %apr(:,i)  = diag(xs*v2*xs');    %  <--  use this if nt not too large
    apr(:,i)   = (apr(:,i) + ones(nt,1)*e1*mc*v1*pay(i,:)')./(xs*mc'*e1');
    tem1(:)    = apr(1:T*nn,i);
    tem2((i-1)*T+1:i*T,:) = tem1;
  end;

  % data: xgmm(:,5:4+nas) 

  tem2    = sum(tem2')'/nn;
  tem1    = zeros(T,nas);
  tem1(:) = tem2;
  tem     = xgmm(2:T+1,5:4+nas)-tem1;
  g       = [g;sum(tem)'/T];
end;

if nargout > 1;
  del    = diag(max(abs(Theta)*1e-5,1e-8));
  for i=1:p;
    dg(:,i) = .5*(gmmh2(Theta+del(:,i))-gmmh2(Theta-del(:,i)))/del(i,i);
  end;
  dg = dg';
end;
