proc iml;

/***** (1) *****/
/* Read data into matrix language */
use hw3;
read all var{crimes} into y;
read all var{total_award} into x;

/* Add column of ones */
x = j(nrow(x),1,1)||(x/10000);

/* Perform OLS */
b_hat = inv(x`*x)*x`*y;

/* Variance-covariance matrix and vector of standard errors */
var_b = ((y-x*b_hat)`*(y-x*b_hat))/(nrow(y)-ncol(x))*inv(x`*x);
se_b = sqrt(vecdiag(var_b));

/* Calculate vector of t-statistics */
t_stat_b = b_hat / se_b;
p_value = 2*(1 - cdf("t",abs(t_stat_b),nrow(y)-ncol(x)));

print b_hat se_b t_stat_b p_value;


/**** (2) ****/
/* Show asymptotic consistency */

/* 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;

xx = j(nrow(x),1,1)||x;
yxx = y||xx;

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); 
finish;

iter1 = 100;
b_ = j(iter1,2,.);
v_ = j(iter1,2,.);

do i=1 to iter1;
	s = samplereplace(yxx,i*iter1);
	y_ = s[,1];
	xx_ = s[,2:3];
	b_[i,] = t(inv(xx_`*xx_)*xx_`*y_);
	v_[i,] = t(vecdiag(((y_-xx_*b_[i,]`)`*(y_-xx_*b_[i,]`))/
			   (nrow(y_)-ncol(xx_))*inv(xx_`*xx_)));
end;
print b_ v_;