% hw10_4.m
% suggested solutions to hw10, part 4

% load data and calculate sample size
load hw10.mat
n=length(Hosp);

% part a: means of treated and untreated groups
m1=mean(Hosp(logical(Flushot)))
m0=mean(Hosp(logical(1-Flushot)))

% part b: add Age as a regressor, calculate robust SE's
X = [ones(n,1) Flushot Age];
b = inv(X'*X)*X'*Hosp;

e = Hosp - X*b;
s2 = e'*e/n;

var_ols_standard = s2 * inv(X'*X);
Xtilde = (e * ones(1,3) ).*X;
Delta = (1/n)* Xtilde'*Xtilde;
var_ols_robust = n*inv(X'*X)*Delta*inv(X'*X);

se_ols_standard = sqrt(diag(var_ols_standard));
se_ols_robust = sqrt(diag(var_ols_robust));

%report results for part b

[b se_ols_standard se_ols_robust]

% part d: 2SLS analysis
Z = [ones(n,1) Letter Age];

% first stage: Regress Flushot on Z
Fitted = Z*inv(Z'*Z)*Z'*Flushot;
R2 = sum( (Fitted - mean(Flushot)).^2 )/ sum( (Flushot - mean(Flushot)).^2 );
Fstat = ((n-3)/2)*(R2/(1-R2))

% call tsls.m to do TSLS and calculate SEs

[b_2sls, se_2sls_1, se_2sls_2] = tsls(Hosp,X,Z);

[b_2sls, se_2sls_1, se_2sls_2]

% end hw10_4.m