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

X = [ones(n,1) treat x3 x4];
betahat = inv(X'*X)*X'*y;
e = y-X*betahat;
s2 = e'*e/(n-4);
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 0 1 0
     0 0 0 1];
r = zeros(2,1);
W = (R*betahat - r)'*inv(s2 * R * inv(X'*X)* R')*(R*betahat - r)
% end hw74.m