rm(list=ls());
set.seed(12345789);
nrep=1000;
nobs=2000;
nvar=3;
X = matrix(rnorm(nobs*nvar,mean=10,sd=2),ncol=nvar);
b= c(1,2,3);
Bhat=matrix(0,ncol=nvar,nrow=nrep);
that=matrix(0,ncol=nvar,nrow=nrep);
for (i in 1:nrep)
{ e=0.1*rnorm(nobs);
  y = X %*% b + e ;
  Bhat[i,]= solve(t(X) %*% X ) %*% t(X) %*% y 
  sig2= t(y-X %*% Bhat[i,]) %*% (y-X %*% Bhat[i,])/nrow(X);
  that[i,1]=1/sqrt(sig2)*(Bhat[i,1]-b[1])/sqrt(solve(t(X) %*% X)[1,1]);
  that[i,2]=1/sqrt(sig2)*(Bhat[i,2]-b[2])/sqrt(solve(t(X) %*% X)[2,2]);
  that[i,3]=1/sqrt(sig2)*(Bhat[i,3]-b[3])/sqrt(solve(t(X) %*% X)[3,3]);
  }
  
  hist((Bhat[,1]-b[1]),prob=TRUE);
  hist((that[,1]),prob=TRUE);
   curve(dnorm(x),-3,3,add=T,col="red");
   
 hist((Bhat[,2]-b[2]),prob=TRUE);
 hist((Bhat[,3]-b[3]),prob=TRUE);