function [H,dH]=gmmh(theta);
%GMMH   A user-supplied function for GMM estimation.
%       [H,dH]=gmmh(theta) is called by gmm.m and gmmj.m.
%       H is nobs x m 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']
%

%       Code for H, dH for TAX1 model
%
l=7; lags=3; lead=4; tm4=3; tm3=4; tm2=5; tm1=6; t=7; t1=8; t2=9; t3=10; t4=11;
c=xgmmc; 
g=xgmmg; 
k=xgmmk; 
n=xgmmn; 
tk=xgmmtk; 
tn=xgmmtn; 
y=xgmmy; 


% parameters estimated
gamma=theta(1); 
thet=theta(2);
phi1=theta(3);
phi2=theta(4); 
phi3=theta(5);
beta=theta(6);
a0=theta(7);

% parameters not estimated:
delta=0.023; 
itc=0; 
phi4=1-phi1-phi2-phi3;
unorm=1304.5; 


H=zeros(nobs,mgmm);
e=ones(nobs,1);
%x(:,4)=zeros(nobs,1);
%x(:,10)=zeros(nobs,1);
%X=[y(t)./c(t),n(t),n(tm1),taun(t),c(t4)./c(t3),c(t3)./c(t2),c(t2)./c(t1), ..
%   c(t1)./c(t),k(t3)./y(t4),tauk(t4),g(t1)./y(t1),g(t)./y(t),g(tm1)./  ..
%   y(tm1),n(t1),k(t)./y(t1),k(tm1)./y(t),k(tm2)./y(t1),y(t1)./y(t), ..
%   g(t1)./g(t),k(t)./k(tm1)];
%xgmm=[X(5:N,:),X(4:N-1,:),X(3:N-2,:),X(2:N-3,:),X(1:N-4,:)];


H(:,1)=a0*gamma*e./(unorm-a0*n(:,t)-(1-a0)*n(:,tm1))-thet*y(:,t).*(e-tn(:,t)) ..
       ./(c(:,t).*n(:,t))+beta*(1-a0)*gamma*e./(unorm-a0*n(:,t1)-(1-a0)*n(:,t));
H(:,2)=beta^4*e./(c(:,t4)./c(:,t)).*(tk(:,t4)*delta+ ..
       y(:,t4)./k(:,t3)*(1-thet).*(e-tk(:,t4))+(1-itc)*(1-delta)*phi1);
H(:,2)=H(:,2)+(1-itc)*(beta^3*((1-delta)*phi2-phi1)*e./(c(:,t3)./c(:,t)) ..
       +beta^2*((1-delta)*phi3-phi2)*e./(c(:,t2)./c(:,t))+  ..
       beta*((1-delta)*phi4-phi3)*e./(c(:,t1)./c(:,t))-phi4);

H(:,3)=H(:,2).*c(:,t)./c(:,tm1);
H(:,4)=H(:,2).*c(:,tm1)./c(:,tm2);
H(:,5)=H(:,2).*c(:,tm2)./c(:,tm3);
H(:,6)=H(:,2).*k(:,tm1)./k(:,tm2);
H(:,7)=H(:,2).*k(:,tm2)./k(:,tm3);
H(:,8)=H(:,2).*k(:,tm3)./k(:,tm4);
H(:,9)=H(:,2).*y(:,tm1)./y(:,tm2);
H(:,10)=H(:,2).*y(:,tm2)./y(:,tm3);
H(:,11)=H(:,2).*y(:,tm3)./y(:,tm4);
H(:,12)=H(:,1).*c(:,t)./c(:,tm1);
H(:,13)=H(:,1).*c(:,tm1)./c(:,tm2);
H(:,14)=H(:,1).*c(:,tm2)./c(:,tm3);
H(:,15)=H(:,1).*k(:,tm1)./k(:,tm2);
H(:,16)=H(:,1).*k(:,tm2)./k(:,tm3);
H(:,17)=H(:,1).*k(:,tm3)./k(:,tm4);
H(:,18)=H(:,1).*y(:,tm1)./y(:,tm2);
H(:,19)=H(:,1).*y(:,tm2)./y(:,tm3);
H(:,20)=H(:,1).*y(:,tm3)./y(:,tm4);
if nargout>1;
  dH=zeros(nobs,mgmm*l);
  dev=diag(.00001*max(abs(theta),3.24e-8*ones(theta)));
% dev=diag(3.24e-8*ones(theta)));
  for i=1:l;
    dH(:,[i:l:mgmm*l])=(gmmh(theta+dev(:,i))-gmmh(theta-dev(:,i)))/(2*dev(i,i));
  end;
end;

