/* PRETEST ESTIMATORS */

/* Create Artificial Data */

nobs = 50;
loops = 2000;
x0 = ones(nobs,1);
x1 = seqa(1,0.05,nobs);
x2 = rndn(nobs,1);
xmat = x0~x1~x2;
invxx = inv(xmat'*xmat);
b = {0.5, 1, 2};
k = rows(b);

y = zeros(nobs,loops);
for ii(1,loops,1);
  	y[.,ii] = xmat*b + 10*rndn(nobs,1);
endfor;

/* Unrestricted Estimator */

urb = zeros(k,loops);
err = zeros(nobs,loops);
s2 = zeros(loops,1);
for ii(1,loops,1);
	urb[.,ii] = invxx*(xmat'*y[.,ii]); 
	err[.,ii] = y[.,ii] - xmat*urb[.,ii];
	s2[ii,1] = (err[.,ii]'*err[.,ii])/(nobs - k);
endfor;

/* Restricted Estimator */

rb = zeros(k-1,loops);
rxmat = xmat[.,1:2];
for ii(1,loops,1);
	rb[.,ii] = inv(rxmat'*rxmat)*(rxmat'*y[.,ii]);
endfor;

/* Pretest Estimator */

tcrit = 2.00;
seurb = zeros(loops,1);
tstat = zeros(loops,1);
ptb = zeros(loops,1);
for ii(1,loops,1);
	seurb[ii,1] = s2[ii,1]*invxx[3,3];
	tstat[ii,1] = urb[k,ii]/sqrt(seurb[ii,1]);
	if tstat[ii,1] < tcrit;
		ptb[ii,1] = rb[k-1,ii];
	else;
		ptb[ii,1] = urb[k-1,ii]; 
	endif;
endfor;

/* Mean Squared Error Criterion */

urmse = zeros(loops,1);
rmse = zeros(loops,1);
ptmse = zeros(loops,1);
for ii(1,loops,1);
	urmse[ii,1] = (urb[2,ii] - b[2,1])^2;
	rmse[ii,1] = (rb[2,ii] - b[2,1])^2;
	ptmse[ii,1] = (ptb[ii,1] - b[2,1])^2;
endfor;

avgurmse = meanc(urmse);
avgrmse = meanc(rmse);
avgptmse =  meanc(ptmse);

print "Unrestricted Average MSE = " avgurmse;
print;
print "Restricted Average MSE = " avgrmse;
print;
print "Pretest Average MSE = " avgptmse;