/* ********************************************************* */
/* COMPARISON OF SIMULTANEOUS EQUATION ESTIMATION TECHNIQUES */
/* ********************************************************* */
cls;
/* ********** */
/* Parameters */
/* ********** */

nobs = 228;
rounds = 50;
k = 2;
b0 = 2;
b1 = 0.8;
sigma = 1;

/* *********** */
/* Create Data */
/* *********** */

g = rndu(nobs,1);
i = rndu(nobs,1);
nx = rndu(nobs,1);
constant = ones(nobs,1);
b = zeros(k,rounds);
tols = zeros(k,rounds);
betahat = zeros(k,rounds);
t2sls = zeros(k,rounds);
 
/* ***************** */
/* Beginning of Loop */
/* ***************** */

for ii(1,rounds,1);
  c = (1/(1-b1))*(b0 + b1*(i+nx+g) + sigma*rndn(nobs,1));
  y = c + i + nx + g;
  xmat = constant~y;

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

  b[.,ii] = inv(xmat'*xmat)*(xmat'*c);
  resids = c - xmat*b[.,ii];
  s2 = (resids'*resids)/(nobs - k);
  varb = s2*inv(xmat'*xmat);
  tols[.,ii] = b[.,ii] ./ diag(sqrt(varb));

  /* *************** */
  /* 2SLS Estimation */
  /* *************** */

  wmat = constant~i~nx~g;
  yhat = wmat*inv(wmat'*wmat)*(wmat'*y);
  zmat = constant~yhat;
  betahat[.,ii] = inv(zmat'*zmat)*(zmat'*C);
  resids = c - xmat*betahat[.,ii];
  sigma2hat = (1/nobs)*(resids'*resids);
  varbetahat = sigma2hat*inv(zmat'*zmat);
  t2sls[.,ii] = betahat[.,ii] ./ diag(sqrt(varbetahat));
endfor;

/* **************************** */
/* Calculate and Print Averages */
/* **************************** */

avgb = meanc(b');
avgtols = meanc(tols');
avgbetahat = meanc(betahat');
avgt2sls = meanc(t2sls');

print "OLS Estimates                    tstats";
print avgb~avgtols;
print;

print "2SLS Estimates                    tstats";
print avgbetahat~avgt2sls;
print;