function [gmm se asyvar] = expiv1(Y,X,Z,gmm0)
%% Calculates GMM1 point estimator and standard error 
%% for a two-way exponential regression model on a (balanced) n-by-m panel data set
%% using instrumental variables 
% INPUT: 
% Y = n-by-m array of outcomes; 
% X = d-dimensional cell of n-by-m matrices of regressors ; 
% Z = q-dimensional cell of n-by-m matrices of instruments;  Z must include exogenous X
% starting value for optimization (gmm0).
% OUTPUT: 
% point estimate (gmm); 
% standard error (se);
% asymptotic covariance matrix (asyvar)

% NOTES: 
% This algorithm is optimized for speed in both optimization and calculation of the standard error

%% Demean regressors to deal with possible non-negativity of covariates
[d] = size(gmm0,1); for k=1:d, X{k} = X{k}-mean(mean(X{k})); end
[q] = size(   Z,2);

% Minimization of GMM problem by Newton's method, starting at gmm0
% one-step estimator
[gmm1, condition, numiter, S] = Newton(@QuadraticForm,gmm0,Y,X,Z,eye(q)); 
% Estimation of variance of moment conditions
[n m] = size(Y); nn = nchoosek(n,2); mm = nchoosek(m,2); rho = nn*mm; J = (S/rho);
[V] = Vmatrix(gmm1,Y,X,Z); 
% two-step estimator
[gmm2, condition, numiter, S] = Newton(@QuadraticForm,gmm1,Y,X,Z,inv(V)); 
% Re-estimation of variance of moment conditions
[n m] = size(Y); nn = nchoosek(n,2); mm = nchoosek(m,2); rho = nn*mm; J = (S/rho);
[V] = Vmatrix(gmm2,Y,X,Z); 
% Construction of standard errors
gmm = gmm2; 
Upsilon = inv(J'*inv(V)*J); asyvar = Upsilon/(n*m); se = sqrt(abs(diag(Upsilon))/(n*m));

%% Evaluation of GMM problem at fixed parameter value psi
function [criterion score Hessian H] = QuadraticForm(psi,Y,X,Z,invV)
% dimensions
[n, m] = size(Y); % sample size
[d, ~] = size(X); % number of regressors
[q, ~] = size(Z); % number of moment conditions

% variable definitions
index = zeros(n,m); for k=1:d, index = index+X{k}*psi(k); end % linear index 
phi   = exp(index); % exponential transform  
error =     Y./phi; % disturbance
d_error = cell(d,1); for k=1:d, d_error{k} = error.*X{k}; end % derivative of disturbance
% averages
error_i = sum(error,2);  
error_j = sum(error,1); 
d_error_i =  cell(d,1); for k=1:d, d_error_i{k} = sum(d_error{k},2); end
d_error_j =  cell(d,1); for k=1:d, d_error_j{k} = sum(d_error{k},1); end 
m_error = sum(sum(error)); for k=1:d, m_derror{k} = sum(sum(d_error{k})); end
% score vector
S = zeros(q,1); for k=1:q, S(k) = sum(sum(error.*Z{k}))*sum(sum(error)) - sum(sum((error_i*error_j).*Z{k})); end
% Jacobian matrix
H = zeros(q,d); 
for k=1:q,
    for j=1:d,
        H(k,j) = sum(sum(Z{k}.*error.*(X{j}*m_error+m_derror{j}) - Z{k}.*(error_i*d_error_j{j}+d_error_i{j}*error_j)));
    end
end
H = -H;
%if q>d, V = Vmatrix(psi,Y,X,Z); else V=eye(q); end; invV = inv(V);
% objective function
criterion = -  S'*invV*S;
score     = -2*H'*invV*S;
Hessian   = -2*H'*invV*H;

%% Estimation of the asymptotic variance of the moment conditions
function [mVar] = Vmatrix(psi,Y,X,Z)
% dimensions
[n, m] = size(Y); % sample size
[d, ~] = size(X); % number of regressors
[q, ~] = size(Z); % number of moment conditions
% variable definitions
index = zeros(n,m); for k=1:d, index = index+X{k}*psi(k); end % linear index 
phi   = exp(index); % exponential transform  
error =     Y./phi; % disturbance
xerror = cell(q,1); for k=1:q, xerror{k} = error.*Z{k}; end 

uZu = cell(q,1); for k=1:q, uZu{k} = error*Z{k}'*error; end

u   = sum(sum(error)); xu = cell(q,1); for k=1:q, xu{k} = sum(sum(xerror{k})); end
u_i = sum(error,2); xu_i =  cell(q,1); for k=1:q, xu_i{k} = sum(xerror{k},2) ; end 
u_j = sum(error,1); xu_j =  cell(q,1); for k=1:q, xu_j{k} = sum(xerror{k},1) ; end 
xuu_j = cell(q,1); for k=1:q, xuu_j{k} = sum(Z{k}.*(u_i*ones(1,m)),1); end
xuu_i = cell(q,1); for k=1:q, xuu_i{k} = sum(Z{k}.*(ones(n,1)*u_j),2); end

xi = cell(q,1); mVar = zeros(q,q);
for k=1:q,
    xi{k} = error.*(Z{k}*u+xu{k})-(Z{k}.*(u_i*u_j)+uZu{k}) + (xu_i{k}*u_j+u_i*xu_j{k}) - error.*(xuu_i{k}*ones(1,m)+ones(n,1)*xuu_j{k});
    xi{k} = 4*xi{k}/((n-1)*(m-1));
end
for k=1:q,
    for j=1:q,
        mVar(k,j) = mean(mean(xi{k}.*xi{j})); 
    end
end

%% Newton algorithm used for optimization
function [x condition it J]=Newton(FUN,x,varargin) % varargout
% maximises FUN, starting at x by Newton-Raphson method
tol=1e-5; maxit=100; smalleststep=.5^20;
it=1; condition=1; improvement=1; k=length(x);
[f g H J] =feval(FUN,x,varargin{:}); %varargout
while it<=maxit && condition==1 && improvement==1;
    [s1 s2]=size(H); if s1==s2 && s2>1 d=-inv(H)*g; else d=-g./H; end      
    step=1; improvement=0;
    while step>=smalleststep && improvement==0;
        [ff gg HH JJ] =feval(FUN,x+step*d,varargin{:}); %varargout
        bounded = sum(sum(isnan(HH)))==0 & sum(sum(isinf(HH)))==0;
        if (ff-f)/abs(f)>=-1e-5 & bounded==1;
            improvement=1; condition=sqrt(step*step*(d'*d))>tol & (ff-f)>tol;
            x=x+step*d; f=ff; g=gg; H=HH; J = JJ;
        else
            step=step/2;
        end
    end
    it=it+1;
end
it=it-1;