function [b,s,V,s0,V0,J,pv] = IV_robust_cluster(y,x,z,mem)

% Calculates 2SLS estimator for linear model   y = X*b + e with instrument matrix Z
% Calculates covariance matrix and standard errors robust to
% misspecification and clusters

% Inputs:
%	y	nx1 vector of observations
%	x	nxk matrix of regressors
%	z	nxl matrix of instruments, l>=k (includes exogenous components of x)
%   mem     nx1 vector of cluster membership set mem=(1:n)' for iid

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

n = size(y,1);
k = size(x,2);
l = size(z,2);
G = length(unique(mem));
b =  ((x'*z)/(z'*z)*(z'*x))\((x'*z)/(z'*z)*(z'*y));

idx = repmat(mem,1,G)==kron(ones(n,1),unique(mem)');

e = y - x*b;
ze = z.*(e*ones(1,l));
mu = mean(ze)';
zz = z'*z;
H = (1/n)*(x'*z)/zz*(z'*x);

Psi = -ze/zz*(z'*x)-repmat(((e'*z)/zz*z')',1,k).*x+repmat(((e'*z)/zz*z')',1,k).*(z/zz*z'*x);

if n == G
    Omega = Psi'*Psi/n;
    we = (ze'*ze)/n;
else
    we = zeros(l,l);
    Psig = zeros(G,k);
    for g = 1:G
        zg = z(idx(:,g),:);
        eg = e(idx(:,g));
        zeg = zg'*eg;
        we = we + zeg*zeg';
        Psig(g,:) = sum(Psi(idx(:,g),:));
    end
    we = we/n;
    Omega = Psig'*Psig/n;    
end

V = H\Omega/H';
s = sqrt(diag(V/n));
V0 = H\((x'*z)/zz*we/zz*(z'*x))/H';
s0 = sqrt(diag(V0/n));


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