#-----------------------------#
#     Ueberpruefung von       #
#       Theorem 12.1          #
#-----------------------------#
        

x = seq(from=-5,to=5,by=1)
x = seq(from=-5,to=5,by=2)
x
n = length(x)

x0 = rep(1,n)
x0

X = cbind(x,x^2)
X

X0 = cbind( x0 , X )
X0

XTXinv = solve( t(X0)%*%X0 )
XTXinv 

XTXinv00 = XTXinv[1,1] 
XTXinv11 = XTXinv[2,2]
XTXinv22 = XTXinv[3,3]
sd0 = sqrt(XTXinv00)
sd1 = sqrt(XTXinv11)
sd2 = sqrt(XTXinv22)

N = 10000

T0 = rep(0,N)
T1 = rep(0,N)
T2 = rep(0,N)
beta0 = -6
beta1 = -1
beta2 = 0.5
sigma = 2

for(k in 1:N)
{
  eps = rnorm(n,mean=0,sd=sigma)
  y = beta0*x0 + beta1*x + beta2*x^2 + eps
  
  res = lm( y ~ X )

  hatbeta = res$coef
  hats = sqrt( 1/(n-3) * sum(res$res^2) )     # res$res: das erste res ist unsere Variable res, das zweite res ist 
                                              # Abkuerzung fuer residuals, tut die lm()-Funktion immer mitberechnen
  T0[k] = (hatbeta[1]-beta0)/(hats*sd0)
  T1[k] = (hatbeta[2]-beta1)/(hats*sd1)
  T2[k] = (hatbeta[3]-beta2)/(hats*sd2)
}



hist(T0,prob=TRUE,breaks=50)
curve(dt(x,df=n-3),col="red",add=TRUE)
curve(dnorm(x,mean=0,sd=1),col="blue",add=TRUE)

hist(T1,prob=TRUE,breaks=50)
curve(dt(x,df=n-3),col="red",add=TRUE)
curve(dnorm(x,mean=0,sd=1),col="blue",add=TRUE)

hist(T2,prob=TRUE,breaks=50)
curve(dt(x,df=n-3),col="red",add=TRUE)
curve(dnorm(x,mean=0,sd=1),col="blue",add=TRUE)


# fuer kleine n:
hist(T0,prob=TRUE,xlim=c(-10,10),ylim=c(0,0.4),breaks=100)
curve(dt(x,df=n-3),col="red",add=TRUE)
curve(dnorm(x,mean=0,sd=1),col="blue",add=TRUE)