use mroz;
read all into y var{lwage};
read all into x var{one educ exper exper2};
read all into x1 var{educ};

/* instrument matrix includes all external and internal exogenous variables */
read all into z var {one exper exper2 mothereduc fathereduc}; * iv;
read all into z1 var {mothereduc fathereduc}; * external iv;
read all into z2 var {one exper exper2}; * included iv;
rname = {one educ exper exper2};

/* Hausman Test for Endogeneity Module */
*=======================================;
/* first run 2sls module */

start hausman;
b = inv(x`*x)*x`*y;
sig2ols = ssq(y-x*b)/(n-k);  * sigma-squared hat;
covb = sig2ols*inv(x`*x);    * OLS covariance;
h = (biv-b)`*inv(covbiv-covb)*(biv-b); * Constrast test;
df_sas = round(trace(ginv(covbiv-covb)*(covbiv-covb)));
p_sas = 1-probchi(h,df_sas);
p_correct = 1-probchi(h,k1);
print / "Hausman contrast test, df and pvalues", h df_sas p_sas p_correct;


* Durbin-Wu-Hausman chi-square test;
*----------------------------------;
hols = (biv-b)`*ginv(inv(xhat`*xhat)-inv(x`*x))*(biv-b)/sig2ols;
pols = 1-probchi(hols,k1);
print ,, "Hausman test and p-value using OLS variance", hols pols;
v1 = x1-z*inv(z`*z)*z`*z1;  * first stage residual;
xa = v1||x;
ba = inv(xa`*xa)*xa`*y;
sse_u = ssq(y-xa*ba);
sig2a = sse_u/(n-ncol(xa));  * estimated error var;
sse_r = ssq(y-x*b);
F = (sse_r - sse_u)/(k1*sig2a);
p = 1-probf(F,k1,n-ncol(xa));
print ,, "Regression based Hausman Test, df and p-value", F k1 p;
finish;
run hausman;