function [theta,Jstat,df1,thetac,Rstat,df2,se]=mlgmm(theta0,n,wgt,nlag, ..
                                                     restrict,induc)
%MLGMM Main file for ML/GMM estimation. 
%      [theta,Jstat,df1,thetac,Rstat,df2,se]=gmm(theta0,n,wgt,restrict,induc)
%      finds the GMM estimates for theta where:
%
%             E [ h(x(t+n),theta) | I(t) ] = 0,  t=1,...,T,  h: m x 1
%
%      and x(t) is a vector of variables in the information set, I(t).
%      Estimates are chosen to minimize:    J = g(theta)' W g(theta) where
%       g(theta) =  1/T sum_{t=1:T} kron( h(x(t+n),theta), z(t) )
%          W     =  R(0)+ sum_{j=1:n-1} (R(j) + R(j)')       if wgt=0 and W>0
%                =  R(0)+ sum_{j=1:n-1} (1-j/n) (R(j)+R(j)') if wgt=1 or W<0 
%                =  W described in Eichenbaum, et.al. (1986) if wgt=2
%         R(j)   =  1/T sum_{t=1+j:T}  kron(h(x(t+n),theta),z(t))*  ..
%                                               kron(h(x(t+n),theta),z(t))
%         z(t)   = a q x 1 vector of instrumental variables in I(t)
%      Also returned are Hansen's J statistic, Jstat, (distributed as chi-
%      square with df1 degrees of freedom) testing overidentifying restric-
%      tions, the statustic T*(Jc-J), Rstat, (distributed as chi-square with 
%      df2 degrees of freedom) testing parameter restrictions across two 
%      runs with the same W, and the standard errors of theta, se. The vector, 
%      induc, gives the indices of the theta vector that are not constrained 
%      in computing Rstat. Use MLGMMHC for moments with restrictions imposed.
%      MLGMM calls MLGMMJ, GMM calls GMMJ. See GMMJ,MLGMMJ for differences.
%

%      Ellen R. McGrattan, 12-07-92
%
%      References:
%      ----------
%        [1]  Eichenbaum, Martin S. and Lars P. Hansen, ``Estimating
%             Models with Intertemporal Substitution Using Aggregate
%             Time Series Data,'' NBER Working Paper 2181, 1987. 
%        [2]  Eichenbaum, Martin S. and Lars P. Hansen, ``Estimating
%             Models with Intertemporal Substitution Using Aggregate
%             Time Series Data,'' JOURNAL OF BUSINESS AND ECONOMIC
%             STATISTICS, 8:53-70, 1990. 
%        [3]  Eichenbaum, Martin S., Lars P. Hansen, and Kenneth Singleton,
%             ``A Time Series Analysis of Representative Agent Models of 
%             Consumption and Leisure Choice Under Uncertainty,'' NBER 
%             Working Paper No. 1981, 1986.
%        [4]  Hansen, L.P. and K.J. Singleton, "Generalized Instrumental
%             Variables Estimation of Nonlinear Rational Expectations
%             Models," ECONOMETRICA, 50:1269-1286, 1982.
%        [5]  Hansen, Lars P., "Large-sample Moments Properties of 
%             Generalized Method of Moments Estimators,'' ECONOMETRICA, 
%             50:1029-1054, 1982.
%        [6]  Newey, W.K. and K.D. West, "A Simple, Positive Semi-
%             Definite, Heteroskedasticity and Autocorrelation Consistent 
%             Covariance Matrix," ECONOMETRICA, 55:703-708, 1987.
%
format short e
if nargin<5; restrict=0; induc=0; end;
if nargin<3; wgt=1; end;
if nargin<2; n=1; end;
if n==1; wgt=0; end;
[nobs,tem]=size(xgmm);
if nargin<4; nlag=min(10,nobs-10); end;
[tem,qgmm]=size(zgmm);
lgmm=length(theta0);
if nobs~=tem; error('Incompatible number of observations in x and z'); end;
disp(' ')
disp(' ')
disp('STARTING PARAMETER VALUES')
disp('-------------------------')
disp(' ')
disp([theta0])
disp(' ')
disp(' ')
disp('OPTIONS SPECIFIED')
disp('-----------------')
disp(' ')
fprintf('wgt               = %g \n',wgt)                      
fprintf('nlag              = %g \n',nlag)                      
fprintf('n                 = %g \n',n)                      
fprintf('restrict          = %g \n',restrict)                      
fprintf('induc             = \n')                      
disp([induc])

H=mlgmmhx(theta0);
[tem,mgmm]=size(H);
if nobs~=tem; error('Incompatible number of observations in x and H'); end;
r=mgmm*qgmm;
if r<lgmm; 
  error(['The number of orthogonality conditions must exceed the number', ..
         ' of parameters being estimated'])
end;

Wgmm=eye(r);
[theta,J,grad,code,status]=uncmin(theta0,'mlgmmjx');

disp(' ')
disp(' ')
disp('RESULTS FROM OPTIMIZATION ROUTINE -- STAGE 1   ')
disp('--------------------------------------------   ')
disp(' ')
fprintf('J                 = %12.4e \n',J)                      
disp('theta, dJ/dtheta  = ')
disp([theta,grad])
disp(' ')
if code==1; 
  disp('An optimal solution was found.')
elseif code==2; 
  disp('A lower point could not be found.  Check tolerances in UNCMIN.')
elseif code==3;
  disp('The iteration limit in UNCMIN was reached.')
elseif code==4;
  disp('Too many large steps taken in UNCMIN.  Function may be unbounded.')
end;
disp('Current convergence status:')
if  ~sum(status) ; disp('  Nothing converged.           '); end;
if status(1) == 1; disp('  norm(g)               < epsa '); end;
if status(2) == 1; disp('  norm(g)/(1+f(x))      < gtol '); end;
if status(3) == 1; disp('  norm(dx)/(1+norm(x))  < xtol '); end;
if status(4) == 1; disp('  df/(1+f(x))           < ftol '); end;

H=mlgmmhx(theta);
F=[ ];
for i=1:mgmm;
  F=[F;(ones(qgmm,1)*H(:,i)').*zgmm'];
end;
R=zeros(2*n,r*r);
Wgmm=zeros(r);
for j=1:n;
  tem1=F(:,j:nobs)*F(:,1:nobs-j+1)';
  tem2=tem1';   
  R(2*j-1,:)=tem2(:)';
  R(2*j,:)=tem1(:)';
end;
R(1,:)=zeros(1,r*r);
if ~wgt | n==1;
  tem1=sum(R)';
  Wgmm(:)=tem1;
  Wgmm=inv(Wgmm/nobs);
  if min(eig(Wgmm))<=0; 
    disp('WARNING: wgt set to 1 because W not positive definite')
    wgt=1; 
  end;
end;
if wgt==2;
  X=[ ];
  for i=nlag:-1:1;
    X=[X,F(:,i:nobs-nlag+i-1)'];
  end;
  [beta,sb,resid]=ols(X,F(:,nlag+1:nobs)');
  X=[ ];
  for i=n-1:-1:1;
    X=[X,resid(i:nobs-nlag-n+i,:)];
  end;
  [beta,sb,resid]=ols(X,F(:,nlag+n:nobs)');
  beta=[eye(r);beta];
  for i=1:r;
    BB(:,i)=sum(beta([i:r:n*r],:))';
  end;
%  br=BB'*resid;
  br= BB'*resid';
  Wgmm=br*br'/(nobs-nlag-n+1);
  Wgmm=inv(Wgmm);
  if min(eig(Wgmm))<=0; 
    disp('WARNING: wgt set to 1 because W not positive definite')
    wgt=1; 
  end;
end;  
if wgt==1;
  tem1=([0;ones(2*n-1,1)-sort(abs([-n+1:n-1]'/n))])*ones(1,r*r);
  R=tem1.*R;
  tem2=sum(R)';
  Wgmm(:)=tem2;
  Wgmm=inv(Wgmm/nobs);
end;

[theta,J,grad,code,status]=uncmin(theta,'mlgmmjx');   

disp(' ')
disp(' ')
disp('RESULTS FROM OPTIMIZATION ROUTINE -- STAGE 2   ')
disp('--------------------------------------------   ')
disp(' ')
fprintf('J                 = %12.4e \n',J)                      
disp('theta, dJ/dtheta  = ')
disp([theta,grad])
disp(' ')
if code==1; 
  disp('An optimal solution was found.')
elseif code==2; 
  disp('A lower point could not be found.  Check tolerances in UNCMIN.')
elseif code==3;
  disp('The iteration limit in UNCMIN was reached.')
elseif code==4;
  disp('Too many large steps taken in UNCMIN.  Function may be unbounded.')
end;
disp('Current convergence status:')
if  ~sum(status) ; disp('  Nothing converged.           '); end;
if status(1) == 1; disp('  norm(g)               < epsa '); end;
if status(2) == 1; disp('  norm(g)/(1+f(x))      < gtol '); end;
if status(3) == 1; disp('  norm(dx)/(1+norm(x))  < xtol '); end;
if status(4) == 1; disp('  df/(1+f(x))           < ftol '); end;

Jstat=nobs*J;
df1=r-lgmm;

if restrict;
  [thetac,Jc,grad,code,status]=uncmin(theta(induc),'mlgmmjc');
  disp(' ')
  disp(' ')
  disp('RESULTS FROM OPTIMIZATION ROUTINE -- WITH CONSTRAINED OBJECTIVE   ')
  disp('---------------------------------------------------------------   ')
  disp(' ')
  fprintf('Jc                = %12.4e \n',Jc)                      
  disp('thetac, dJc/dthetac = ')
  disp([thetac,grad])
  disp(' ')
  if code==1; 
    disp('An optimal solution was found.')
  elseif code==2; 
    disp('A lower point could not be found.  Check tolerances in UNCMIN.')
  elseif code==3;
    disp('The iteration limit in UNCMIN was reached.')
  elseif code==4;
    disp('Too many large steps taken in UNCMIN.  Function may be unbounded.')
  end;
  disp('Current convergence status:')
  if  ~sum(status) ; disp('  Nothing converged.           '); end;
  if status(1) == 1; disp('  norm(g)               < epsa '); end;
  if status(2) == 1; disp('  norm(g)/(1+f(x))      < gtol '); end;
  if status(3) == 1; disp('  norm(dx)/(1+norm(x))  < xtol '); end;
  if status(4) == 1; disp('  df/(1+f(x))           < ftol '); end;

  Rstat=(Jc-J)*nobs;
  df2=lgmm-sum(induc>0);
else;
  thetac=0;
  Rstat=0;
  df2=0;
end;

if nargout>6;
  [H,dH]=mlgmmhx(theta);
  F=[ ];
  for i=1:mgmm;
    F=[F;(ones(qgmm,1)*H(:,i)').*zgmm'];
  end;
  R=zeros(2*n,r*r);
  Wgmm=zeros(r);
  for j=1:n;
    tem1=F(:,j:nobs)*F(:,1:nobs-j+1)';
    tem2=tem1';   
    R(2*j-1,:)=tem2(:)';
    R(2*j,:)=tem1(:)';
  end;
  R(1,:)=zeros(1,r*r);
  if ~wgt;
    tem1=sum(R)';
    Wgmm(:)=tem1;
    Wgmm=inv(Wgmm/nobs);
    if min(eig(Wgmm))<0; wgt=1; end;
  end;
  if wgt;
    tem1=([0;ones(2*n-1,1)-sort(abs([-n+1:n-1]'/n))])*ones(1,r*r);
    R=tem1.*R;
    tem2=sum(R)';
    Wgmm(:)=tem2;
    Wgmm=inv(Wgmm/nobs);
  end;

  D=zeros(r,lgmm);
  for i=1:nobs;
    tem1=zeros(lgmm,mgmm);
    tem1(:)=dH(i,:);
    D=D+kron(tem1',zgmm(i,:)')/nobs;
  end;
  if rank(D)<lgmm; disp('WARNING: D does not have full rank.'); end;
  ivcov=D'*Wgmm*D;                         
  if rank(ivcov)<lgmm; 
    disp(['WARNING: D''*W*D is not nonsingular.  ', ..
          'Standard errors can not be computed.  ']) 
    se=0;
  else;
    se=sqrt(diag(inv(ivcov)/nobs));       % vcov=inv(ivcov)/nobs; 
  end;
  Dgmm=D;
end;
