% hw73.m
% Q3, part b
load hw6_nohead.dat
y = hw6_nohead(:,10);
n = length(y);
treat = hw6_nohead(:,11);
x3 = hw6_nohead(:,3) > 12;

X = [ones(n,1) treat x3];
betahat = inv(X'*X)*X'*y;
e = y-X*betahat;
s2 = e'*e/(n-3);
Vhat = s2*inv(X'*X);
SE = sqrt(diag(Vhat));

% report results:
betahat
SE
disp('Confidence Intervals')
[ (betahat - 1.96*SE), (betahat + 1.96*SE)]

% test beta_2 = beta_3 = 0
R = [0 1 0
     0 0 1];
r = zeros(2,1);
W = (R*betahat - r)'*inv(s2 * R * inv(X'*X)* R')*(R*betahat - r)
% end hw73.m