% hw8_3.m: Solution to HW8, Question 3
% calls ols.m and ols_het.m

% load data
load ruud/wage_nohead.dat
y=log(wage_nohead(:,1));
fe = wage_nohead(:,2);
nw = wage_nohead(:,3);
un = wage_nohead(:,4);
ed = wage_nohead(:,5);
age = wage_nohead(:,7);
n = length(y);
X = [ones(n,1) fe nw un ed age];

% part a: ols, conventional and Eicker-White SEs

[b, se, s2] = ols(y,X);
[b_r, se_r] = ols_het(y,X);

[b, se, se_r]

% part b: test heteroskedasticity

fe_index = logical(fe);
y_f = y(fe_index);
X_f = X(fe_index,[1 3 4 5 6]);
[b_f, se_f, s2_f] = ols(y_f, X_f);

y_m = y(~fe_index);
X_m = X(~fe_index, [1 3 4 5 6]);
[b_m, se_m, s2_m] = ols(y_m, X_m);

Fstat = s2_m/s2_f

n_f=sum(fe);
n_m = n-n_f;
dof_num = n_m-6
dof_den = n_f-6

% test should reject if Fstat < 0.86 or Fstat > 1.17

% part c: FGLS
var_vector = s2_f*fe + s2_m*(1-fe);
Omega = diag(var_vector);
Omega_inv = inv(Omega);
V_gls = inv(X'*Omega_inv*X);
b_gls = V_gls * X'*Omega_inv*y;
se_gls = sqrt(diag(V_gls));
[b_gls, se_gls]