function muhat = imbens_newey(Y,X,Z,L);
% muhat = imbens_newey(Y,X,Z)
% assumes Y,X,Z are all n x 1 vectors 
% L - vector of values for x to evaluate mu(x) at
% uses just a quadratic function of z for eta-hat estimation
% uses just a linear function of x, eta-hat, and interaction for 2nd stage 

% estimate eta

n = length(Y);
Q = [ones(n,1) Z Z.^2];

etahat=NaN*ones(n,1);
QQinv = inv(Q'*Q);

for i=1:n,
    etahat(i,1) = [1 Z(i) Z(i)^2]*QQinv*Q'*(X <= X(i));
end;

etahat = max(etahat,zeros(n,1));
etahat = min(etahat,ones(n,1));


% estimate beta([1,x,eta,x*eta])

P = [ones(n,1) X etahat X.*etahat];
PPinv = inv(P'*P);
gammahat = PPinv*P'*Y;

etagrid=0.01:.01:.99;
K = length(etagrid);

J=length(L);
muhat = NaN*ones(J,1);
for j=1:J,
	x_j = L(j);
	P_j = [ones(K,1) x_j*ones(K,1) etagrid' x_j*etagrid'];
	muhat(j,1) = mean(P_j*gammahat);
end;