function [b,s,V,s0,V0,J,pv,iter] = iterated_gmm_nl_cen(b0)

% Calculates iterated GMM estimator for non-linear model
% Calculates covariance matrix and standard errors robust to misspecification

% Inputs:
%	y	nx1 vector of observations
%	x	nxk matrix of regressors
%	z	nxl matrix of instruments, l>=k (includes exogenous components of x)

% Outputs:
%	b	kx1 vector of coefficient estimates
%	s	kx1 vector of robust asymptotic standard errors
%	V	kxk robust asymptotic covariance matrix estimate
%	s0	kx1 vector of classic asymptotic standard errors
%	V0	kxk classic asymptotic covariance matrix estimate
%	J	J-statistic for overidentifying restrictions
%	pv	Asymptotic chi-square p-value for J-statistic

% Output variables = NaN if the iteration does not converge
% 
% Required functions (save in the same folder)
% moment.m    m(theta), the moment function
% dmoment.m   Q(theta), the Jacobian of the moment
% d2moment.m  R(theta), the second derivative of the moment
% dwc.m       S(theta), the derivative of the vec(W(theta))
% crit.m      J(theta,w) evaluated at the parameter and the weight matrix 

tolerance = 1e-5;
maxit = 1e+3;

m0 = moment(b0);

n = length(m0(:,1));
k = length(b0);
l = length(m0(1,:));

w = eye(l);

critw = @(b)crit(b,w);

options = optimoptions(@fminunc,'Algorithm','trust-region','GradObj','on','Hessian','on','Display','off');
b1 = fminunc(critw,b0,options);

for iter = 1:maxit
    m1 = mean(moment(b1))';
    w = ((moment(b1)-repmat(m1',n,1))'*(moment(b1)-repmat(m1',n,1)))/n;

    critw = @(b)crit(b,w);
    
    b = fminunc(critw,b1,options);
    if (norm(b-b1) < tolerance)
      break
   end   
   b1 = b;
end

if iter == maxit
    
    b = NaN;
    s = NaN;
    V = NaN;
    s0 = NaN;
    V0 = NaN;
    J = NaN;
    pv = NaN;
    
    return
end
    
m = moment(b);
Q = dmoment(b);
R = d2moment(b);

mu = mean(m)';
Qn = mean(Q,3);
Rn = mean(R,3);
Sn = dwc(b);

H = Qn'/w*Qn + kron(mu'/w,eye(k))*Rn - kron(mu'/w,Qn'/w)*Sn;

Omega = zeros(k,k);
for i = 1:n
    h = Qn'/w*m(i,:)'+Q(:,:,i)'/w*mu-Qn'/w*(m(i,:)'-mu)*(m(i,:)'-mu)'/w*mu;
    Omega = Omega + h*h';
end

Omega = Omega/n;

V = H\Omega/H';
s = sqrt(diag(V/n));

V0 = inv(Qn'/w*Qn);
s0 = sqrt(diag(V0/n));

if l>k
  J = (mu'/w*mu)*n;
  pv = chi2cdf(J,l-k,'upper');
else
  J = 0;
  pv = 1;
end