/* SUR ESTIMATION -- AN APPLICATION TO WAGES AND BENEFITS */
cls;
/* Load Data */
load path = c:\gauss8.0\classes\econ5360\data\;
load data[617,12] = fringe.txt;

/* Define Variables */
wage = data[2:617,1];
exper = data[2:617,2];
expersq = exper^2;
married = data[2:617,3];
tenure = data[2:617,4];
tenuresq = tenure^2;
educ = data[2:617,5];
nrthe = data[2:617,6];
nrthc = data[2:617,7];
south = data[2:617,8];
male = data[2:617,9];
white = data[2:617,10];
union1 = data[2:617,11];
hrben = data[2:617,12];
T = 616;
n = 2;

/* ************** */
/* OLS Estimation */
/* ************** */

y1 = wage;
y2 = hrben;
x1 = ones(T,1)~educ~exper~expersq~tenure~tenuresq~
     union1~south~nrthe~nrthc~married~male~white;
x2 = ones(T,1)~educ~exper~expersq~tenure~tenuresq~
     union1~south~nrthe~nrthc~married~male~white;
k = cols(x1);

b1 = inv(x1'*x1)*(x1'*y1);
e1 = y1 - x1*b1;
varb1 = ((e1'*e1)/(T-k))*inv(x1'*x1);
t1 = b1 ./ diag(sqrt(varb1));

b2 = inv(x2'*x2)*(x2'*y2);
e2 = y2 - x2*b2;
varb2 = ((e2'*e2)/(T-k))*inv(x2'*x2);
t2 = b2 ./ diag(sqrt(varb2));

print "OLS Estimates (Wage Eq)                    tstats";
print b1~t1;
print;
print "OLS Estimates (Benefit Eq)                 tstats";
print b2~t2;
print;

/* *************** */
/* FGLS Estimation */
/* *************** */

sigma = zeros(2,2);
s11 = e1'*e1/(T-k);
s12 = e1'*e2/(T-k);
s22 = e2'*e2/(T-k);
sigma[1,1] = s11;
sigma[1,2] = s12;
sigma[2,1] = s12;
sigma[2,2] = s22;

id = eye(T);
invomega = inv(sigma).*.id;
x = (x1|zeros(T,k))~(zeros(T,k)|x2);
y = y1|y2;
bfgls = inv(x'*invomega*x)*(x'*invomega*y);
varb = inv(x'*invomega*x);
tfgls = bfgls ./ diag(sqrt(varb));

print "FGLS Estimates (Wage Eq)                    tstats";
print bfgls[1:k]~tfgls[1:k];
print;
print "FGLS Estimates (Benefit Eq)                 tstats";
print bfgls[k+1:2*k]~tfgls[k+1:2*k];
print;

/* ******************************** */
/* Contemporaneous Correlation Test */
/* ******************************** */

r12 = s12/sqrt(s11*s22);
lam = T*(r12^2);

print "Test Statistic (df=1) for No Contemporaneous Correlation =" lam;
print;

/* ********************************** */ 
/* Test of Common Coefficient Vectors */
/* ********************************** */ 

R = eye(k)~(-eye(k));
q = zeros(k,1);
W = (R*bfgls-q)'*inv(R*varb*R')*(R*bfgls-q);

print "Test Statistic (df=k) for Common Coefficient Vector =" W;