function [g,dg] = mlgmmgx(Theta);
%MLGMMGX  A user-supplied function for MLGMM estimation.
%         [g,dg]=mlgmmgx(theta) is called by MLGMM and MLGMMJ. 
%         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, 1-15-93
%         Revised, ERM, 2-2-93

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

global xgmm pshocks
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

alpha0 = 500;
alpha1 = 0.5;
beta   = 0.909;
delta  = 0.1;
gamma0 = 0.4;
gamma1 = 0.7;
f1     = 0.001;
f3     = 0.001;
f5     = 0.001;
f7     = 0.001;
gg     = 0.85;
rho1   = 0.0;
rho2   = 0.0;
rho3   = 0.6;
sigma1 = 1;
sigma2 = 1;
sigma3 = 15;
muh    = 50.0;
mum    = 100.0;
x0     = zeros(8,1);
x0(2:5)= [17.56*ones(3,1);1];

no     = 4;
grate  = 1;
varm   = ones(no,1);
Gam    = zeros(no);

beta   = 0.96;
delta  = 0.0;
sig    = 0;
f1     = 1.0247e-03;
f3     = 1.9193e-03;
f5     = 3.8389e-05;
f7     = 5.7707e-03;


if no>3; varm(4)= 1e-7; end

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


%---------------------------------------------------------------------
% 3. Assign elements of Theta to underlying parameters.  Use the vector
%    ind to indicate which parameters are to be estimated and which
%    are to be fixed in estimation.  For example, if the vector of
%    parameters is [beta,gamma,delta] and gamma is to be fixed during
%    estimation then set:
%       beta  = Theta(1);
%       delta = Theta(2);
%       ind   = [1,0,1];
%

%    For the current example, the ordering of the parameter vector is:
%       alpha0,alpha1,beta,delta,gamma0,gamma1,f1,f3,f5,f7,
%       g,rho1,rho2,rho3,sigma1,sigma2,sigma3,muh,mum,
%       grate,sig,varm
%         1      2     3     4     5      6    7  8  9  10
%       11 12   13   14    15     16     17   18  19
%        20    21  25
%

ind    = [1,1,0,0,1,1, 0,0,0,0, 1, 1,0,1, 1,0,1, 1,1, 0, 0, 1,1,1,0];

alpha0 = Theta(1);
alpha1 = Theta(2);
gamma0 = Theta(3);
gamma1 = Theta(4);
gg     = Theta(5);
rho1   = Theta(6);
rho3   = Theta(7);
sigma1 = Theta(8);
sigma3 = Theta(9);
muh    = Theta(10);
mum    = Theta(11);
varm(1)= Theta(12);
varm(2)= Theta(13);
varm(3)= Theta(14);

f2     = 0.5*gamma0*gg/f1;
f4     = 0.5*gamma1*gg/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,gg;1,0,0;0,1,0];
gamma  = [0 0 0; f5*(1-delta) 0 gg*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);
no     = length(varm);

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


%---------------------------------------------------------------------
% 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 and sigin.

if ind(3);  betain=sum(ind(1:3));   else; betain=0;  end;
if ind(20); gratein=sum(ind(1:20)); else; gratein=0; end;
if ind(21); sigin=sum(ind(1:21));   else; sigin=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,
%           gg,rho1,rho2,rho3,sigma1,sigma2,sigma3,muh,mum,
%           grate,sig,varm
Ub = [10000;100;.999;1;10;10;1.;1.;1.;1.;1;.999;.999;.999;1000;1000;
      1000;5000;5000;1.1;100;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;-100;0;0;0;0];
Ub = Ub(ind);
Lb = Lb(ind);
if any(Theta>Ub) | any(Theta<Lb);
  g = 1e+20*ones(size(Theta)); dg=1e+20*ones(size(Theta),sum(ind));
  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')+sig*c*c';   % make sure that b0 is symmetric and p.s.d.
P = q;
PT= (v+sig*P*c*c')\P;
k0= a'*PT*b/(b'*PT*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;
  PT = (v+sig*P*c*c')\P;
  k1 = a'*PT*b/(b'*PT*b+R);
  dd = norm(k1-k0);
  it = it+1;
  a0 = a1;
  k0 = k1;
end
F  = k1'+R\W';

%=====================
%% 7b. Direct Iteration
%=====================
%
%v = eye(nw);
%P = -eye(nx);
%PT= P-sig*P*c/(v+sig*c'*P*c)*c'*P;
%f0= (R+beta*B'*PT*B)\(beta*B'*PT*A+W');
%d=1;
%tol=1e-8;       % reset tolerance if necessary
%while (d > tol)
%  P  = beta*A'*PT*A + Q -(beta*A'*PT*B+W)*f0;
%  PT = P-sig*P*c/(v+sig*c'*P*c)*c'*P;
%  f1 = (R+beta*B'*PT*B)\(beta*B'*PT*A+W');
%  d  = norm(f1-f0);
%  f0 = f1;
%end;
%F  = f1;

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

Ao = A-B*F;

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

%data xgmm(t)=[total stock(t),breeding stock(t),slaughter rate(t),  
%              beef price(t),asset prices(t)]';

T   = length(xgmm);
tem = grate.^([1:T]');  
y   = [xgmm(:,1)./tem,xgmm(:,2)./tem,xgmm(:,3)./tem];
if no>3;
  y = [y,xgmm(:,4)+65-xgmm(1,4)];
end;
%---------------------------------------------------------------------

T   = length(y);
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);
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);
sd  = [zeros(nd,nh+nk),ud];
pay = [ ];
D   = [[0,1-delta+gg,gg,gg,0,0,0,0]-F;
       [0,1-delta,0,gg,0,0,0,0]-F;
       F];
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);

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,
%     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);
end;
%-----------------------------------------------------------------------


%-------------------------------------------------------------------
% 11. Specify grada2,gradc2,graddh,graddk,gradg,gradl,gradpc,gradpg,
%             gradpi,gradph,gradth,gradtk,gradub,gradud,gradx0,
%             gradG,gradO,gradc  
%     which are gradients for
%             a22,c2,deltah,deltak,gamma,lambda,phic,phig,phii,pih,
%             thetah,thetak,ub,ud,x0,Gam,Omega,c,
%     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)=gg;
                            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/gg;
                            gradud(3,18)=f2;
                            gradud(8,5)= f2/gamma0;
                            gradud(8,7)=-f2/f1;
                            gradud(8,11)=f2/gg;
                            gradud(4,6)= f4*muh/gamma1;
                            gradud(4,8)=-f4*muh/f3;
                            gradud(4,11)=f4*muh/gg;
                            gradud(4,18)=f4;
                            gradud(9,6)= f4/gamma1;
                            gradud(9,8)=-f4/f3;
                            gradud(9,11)=f4/gg;
                            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,22)=1; gradO(no+2,23)=1;
                            gradO(2*no+3,24)=1;
                            if no>3; gradO(3*no+4,25)=1; end;
gradc  = zeros(nx*nw,p);    gradc(6,15)=1; gradc(15,16)=1; gradc(24,17)=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);
gradc  = gradc(:,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(nx,nw);
tem1(:)=[1:nx*nw]';
tem2=tem1';
tem2=tem2(:);
tem1=tem1(:);
Kxw=zeros(nx,nw);
for i=1:nx*nw;
  Kxw(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);

gradcp = Kxw*gradc;
pci  = inv(eye(nw)+sig*c'*P*c);
tem1 = eye(nx*nx)-sig*kron(P*c*pci*c',eye(nx))- ...
       sig*kron(eye(nx),P*c*pci*c')+sig*sig*kron(P*c*pci*c',P*c*pci*c');
tem2 = -sig*kron(P*c*pci,P)+sig*sig*kron(P*c*pci,P*c*pci*c'*P);
tem3 = -sig*kron(P,P*c*pci)+sig*sig*kron(P*c*pci*c'*P,P*c*pci);
ri=inv(R+beta*B'*PT*B);
ai=inv(eye(nx*nx)-beta*kron(Ao',Ao')*tem1);
gradP = ai*(gradQ+(beta*kron(eye(nx),Ao'*PT)+beta*kron(Ao'*PT,eye(nx)) ...
       *Kxx)*gradA-beta*kron(F',Ao'*PT)*gradB-beta*kron(Ao'*PT,F')*gradBp);
gradP = gradP+ai*(kron(F',F')*gradR-(kron(F',eye(nx))+kron(eye(nx),F')* ...
        Kxi)*gradW);
gradP = gradP+ai*beta*kron(Ao',Ao')*(tem2*gradc+tem3*gradcp);
gradPT= tem1*gradP+tem2*gradc+tem3*gradcp;
gradF = beta*kron(eye(nx),ri*B'*PT)*gradA-beta*kron(F',ri*B'*PT)* ...
        gradB+beta*kron(Ao'*PT,ri)*gradBp-kron(F',ri)*gradR+ ...
        kron(eye(nx),ri)*Kxi*gradW+ beta*kron(Ao',ri*B')*gradPT;
if betain;
  gradF(:,betain) = gradF(:,betain)-beta*kron(Ao',ri*B')*gradPT(:,betain) ...
                   +vec(ri*B'*PT*Ao);
  gradP(:,betain) = gradP(:,betain)+ai*vec(Ao'*PT*Ao);
  gradPT(:,betain)= tem1*gradP(:,betain);
  gradF(:,betain) = gradF(:,betain)+beta*kron(Ao',ri*B')*gradPT(:,betain);
end;
if sigin;
  gradPT(:,sigin) = (eye(nx*nx)-beta*kron(Ao',Ao')+beta*sig* ...
      kron(P*c*pci*c'*Ao',Ao')+sig*beta*kron(Ao',P*c*pci*c'*Ao') ...
      -sig*sig*beta*kron(P*c*pci*c'*Ao',P*c*pci*c'*Ao'))\ ...
      vec(-PT*c*pci*c'*P);
  gradF(:,sigin)  = gradF(:,sigin)+beta*kron(Ao',ri*B')*gradPT(:,sigin);
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,5,12],sum(ind(1:11)))=ones(5,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; 
dL = T*dL;

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

g = dL;

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

Aos    = (eye(nx)-sig*c*pci*c'*P)*Ao;
Sigmas = c*(eye(nw)+pci)*c';
[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)]*PT(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        = dlyap(sqrt(beta)*Aos,Sigmas)*beta/(1-beta);
  nn        = floor(nt/T);
  tem1      = zeros(T,nn);
  tem2      = zeros(T*nas,nn);
  for i=1:nas;
    v2   = doublej(sqrt(beta)*Aos',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*(mlgmmgx(Theta+del(:,i))-mlgmmgx(Theta-del(:,i)))/del(i,i);
  end;
  dg = dg';
end;


