function [alpha,se1,se2] = tsls(y,X,Z);
% [alpha,se1,se2] = TSLS(y,X,Z)
% y = vector of outcomes, n by 1
% X = matrix of regressors, n by k
% Z = matrix of instruments, n by k2, k2 >= k

k = size(X,2);
n = length(y);

% 2SLS estimate
Xhat = Z*inv(Z'*Z)*Z'*X;
alpha = inv(Xhat'*Xhat)*Xhat'*y;

% standard error calculation under homoskedasticity
e = y - X*alpha;
s2 = e'*e/n;
var1 = s2 * inv(Xhat'*Xhat);
se1 = sqrt(diag(var1));

% part d - variance under heteroskedasticity
Xh_e2 = (e*ones(1,k)).* Xhat;
var2 = inv(X'*Xhat)*(Xh_e2'* Xh_e2)*inv(Xhat'*X);
se2 = sqrt(diag(var2));

% end tsls.m