#----------------#
#   Aufgabe 1    #
#----------------#


x = (-5):5
x
n = 11
alpha = 0.9
xq = qt(0.95,df=9)

beta0 = 3
beta1 = -0.5
sigma = 2
eps = rnorm(n,mean=0,sd=sigma)
y = beta0 + beta1*x + eps
plot(x,y)    


N = 10000
hatbeta0   = rep(0,N)
hatbeta1   = rep(0,N)
hatbeta2   = rep(0,N)
deltabeta0 = rep(0,N)
deltabeta1 = rep(0,N)
deltabeta2 = rep(0,N)

X0 = cbind(rep(1,n),x)        # Matrix der Regressoren mit x0, brauchen 
                              # wir zur Berechnung der deltabeta's

X0TX0inv = solve(t(X0)%*%X0)  # die Matrix (X0^T*X0)^{-1}, mit Vektor x0

sqrt00 = sqrt( X0TX0inv[1,1] )
sqrt11 = sqrt( X0TX0inv[2,2] )


for(k in 1:N)
{
  eps = rnorm(n,mean=0,sd=sigma)
  y = beta0 + beta1*x + eps

  # die lineare Regression:
  res = lm(y ~ x)       

  betas = res$coeff
  hatbeta0[k] = betas[1]
  hatbeta1[k] = betas[2]

  # res$residuals = gegebene Daten y - Regression-Fit
  #               = y - hat(y),   wir brauchen sigma-Hut:
  hatsigma = sqrt( 1/(n-2) * sum(res$residuals^2) )    

  deltabeta0[k] = xq * hatsigma * sqrt00
  deltabeta1[k] = xq * hatsigma * sqrt11
}


# obere und untere Vertrauensintervallgrenzen:
beta0up   = hatbeta0 + deltabeta0
beta0down = hatbeta0 - deltabeta0
beta1up   = hatbeta1 + deltabeta1
beta1down = hatbeta1 - deltabeta1


# Schauen wir uns die Intervallgrenzen mal an:
plot( rep(beta0,N), col="red", ylim=c(0,6) )         # wirkliches beta
points(beta0up)                  
points(beta0down,col="blue")         # ok, sieht nicht unplausibel aus..


# jetzt etwas quantitativer:
test0 = ifelse( beta0<beta0up & beta0>beta0down , 1 , 0 )
sum(test0)                     # ok, nahe bei N*alpha = 9000, all good..


# dasselbe fuer beta1: 
plot( rep(beta1,N), col="red", ylim=c(-2,1) )        # wirkliches beta
points(beta1up)                  
points(beta1down,col="blue")             # sieht nicht unplausibel aus
test1 = ifelse( beta1<beta1up & beta1>beta1down , 1 , 0 )
sum(test1)                                   # nahe bei N*alpha = 9000