proc iml;

/* Set up the simulation */
x = j(10000,1,.);
call randgen(x,"Normal",2,2);

e = j(10000,1,.);
call randgen(e,"Normal",0,1);

y = 2 + 0.5*x + e;

/* Determine the population regression line */
xx = j(nrow(x),1,1)||x;
b = inv(xx`*xx)*xx`*y;
print b;

yxx = y||xx;

/* Look at the bias properties */
/* Specify module that will sample with replacement */
start samplereplace(A, nSamples);
   results = j(nSamples,1);            /** allocate result matrix **/
   call randgen(results, "Uniform");    /** fill with random U(0,1) **/
   if ncol(A)>nrow(A) then n=ncol(A);
   else n=nrow(A);
   results = ceil(n*results);           /** convert to integers 1,2,...n **/
   ss = A[results[1,],];
   do i=2 to nrow(results);
	 ss = ss//A[results[i,],];
   end;
   return(ss); /** reshape and return from A **/
finish;

iter = 1000;
b_ = j(iter,2,.);

/* Perform simulation; use many samples of 200 observations and 
	keep the parameter values */
do i=1 to iter;
	s = samplereplace(yxx,100);
	y_ = s[,1];
	xx_ = s[,2:3];
	b_[i,] = t(inv(xx_`*xx_)*xx_`*y_);
end;

/* Check whether you have bias */
b_mean = b_[:,];
print b_mean;