function Lt = gmmh(Theta);
%GMMH  A user-supplied function for GMM estimation.
%      [H,dH]=gmmh(theta) is called by GMM and GMMJ. H is a
%      nobs x m (or mgmm) matrix with rows containing the m x 1 
%      vectors h(x(t+n),theta), where it is assumed that
%
%          E[ h(x(t+n),theta) | I(t) ] = 0,  t=1,...,nobs.
%
%      dH is a nobs x (m*l) matrix containing the derivatives
%      of H with respect to the elements of theta, i.e.
%      dH(t,:)= [dh1(x(t+n),theta)/dtheta', ..., dhm(x(t+n),theta)/dtheta']
%

%      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    = .1;
lambda = 0;
phi1   = .00001;
pih    = 1;
thetah = (1-deltah);
thetak = 1;
rho1   = .8;
rho2   = .5;
sigma1 = 1;
sigma2 = 0;
sigma3 = 1;
ub1    = 30;
ud11   = 5;
ud12   = 1;
varm   = 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.  
%---------------------------------------------------------------------


%---------------------------------------------------------------------
% 2. 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]; 
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;
%---------------------------------------------------------------------


%---------------------------------------------------------------------
% 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));
lambda = v((nx+1):(2*nx),index(1:nx));
P = real(lambda/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));
%lambda = v((nx+1):(2*nx),index(1:nx));
%P = real(lambda/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)]';

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.   

D = [0,1,0,0,0;
     -F;
      phiin*(gamma*[0,1,0,0,0]+ud*[zeros(3,2),eye(3)] ..
       +phii*F)];
%---------------------------------------------------------------------

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';
Lt(1,1)    = innov(1,:)*Sigmazi*innov(1,:)';
for i=2:T;
  xt(i,:)    = xt(i-1,:)*Ao'+innov(i-1,:)*K';
  innov(i,:) = dy(i,:)-xt(i,:)*Dbar';
  Lt(i,1)    = innov(i,:)*Sigmazi*innov(i,:)';  
end;
sum1=innov(1:T,:)'*innov(1:T,:)/T;

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

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

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

