/* illustration of the sampling distributions of regression coefficients, hypothesis testing 
samp_dist.txt  */

@bivariate regression
  y = 3 + .5*x + epsilon
  epsilon i.i.d. normal N(0,4)  @

output file = samp_dist.out reset;

/* set sample size */
n = 10;

/*draw x from uniform distribution, sample space [0,4]*/
x = rndu(n,1)*4;

print "x draws";
print x;

/* set number of sample replications */
r = 100000;

e_mat = 2*rndn(n,r);

/* generate y matrix for these samples */
y = 3 + .5*x + e_mat;

/* generate estimates of regression coefficients for each sample */
xmat = ones(n,1) ~ x;
bhat_mat = invpd(xmat'xmat)*xmat'y;

print "mean of estimates of regression coefficients";
print meanc(bhat_mat');
print "standard deviation (standard errors) of regression coefficients";
print stdc(bhat_mat');

/* compute s2 estimate of variance */
s2_vec = zeros(1,r);
i = 1;
do until i gt r;
  ypred = xmat*bhat_mat[.,i];
  s2_vec[i] = (y[.,i]-ypred)'(y[.,i]-ypred)/(n-2);
  i = i+1;
endo;

print "mean of s2";
print meanc(s2_vec');
print "standard error of s2";
print stdc(s2_vec');

/* standard errors of regression coefficients from 
    s2 invpd(x'x)  */

xxinv = invpd(xmat'xmat);
se_slope_mat = sqrt(s2_vec*xxinv[2,2]);

print "mean of slope coef stand error";
print meanc(se_slope_mat');
print "st dev of slope coef stand error";
print stdc(se_slope_mat');
print "actual standard error of slope coefficient";
print sqrt(4*xxinv[2,2]);

/* graph distributions of bhat (slope) and s2 */

library pgraph;
graphset;
margin(1,1,0,0);
_pdate = "";
_ptitlht = .2;
_paxht = .15;
window(2,1,0);
setwind(1);
title("Distribution of Slope Coefficient\LN = 50");
xlabel("Slope Estimate");
histp(bhat_mat[2,.]',100);
nextwind;
title("Distribution of s[2]\LN = 50");
xlabel("s[2]");
histp(s2_vec',100);
endwind;

end;