#
#
#
#
#

x1 <- rnorm(20)
x2 <- rnorm(20)
x3 <- rnorm(20)

# wahres Modell
y <- 1 + 2*x1 - 1*x2 + 1*x3 + rnorm(20)

X <- cbind(rep(1,20),x1,x2,x3)

# X'X
t(X) %*% X

# (X'X)^-1
solve( t(X) %*% X)

# (X'X)^-1 X'
solve(t(X) %*% X) %*% t(X)

# (X'X)^-1 X'y
solve(t(X)%*%X) %*% t(X) %*% y

# compare to 
lm(y ~ x1 + x2 +x3)

lm1 <- lm(y ~ x1 + x2 + x3)
summary(lm1)

# predict
beta <- solve(t(X)%*%X) %*% t(X) %*%y

# fitted values
yhat <- X %*% beta

fitted(lm1)
summary(fitted(lm1))
summary(yhat)

summary(lm1)


####
# X'X nicht invertiertbar?

solve(t(X)%*%X)

X2 <- cbind(X, 2 + 3*x1)
cor(X2)

t(X2)%*%X2

solve(t(X2)%*%X2)
x3 <- 1 + 2*x1
lm(y ~ x1 +x2)
#
x <- runif(10000)
x <- x[order(x)]
y <- 2+ 3*x + rnorm(10000)


plot(x,y, pch = 19)

lm1 <- lm(y ~ x)
abline(lm1, col = 2, lwd = 2)

predict(lm1, newdata = NULL, interval = "predict")

predict(lm1, newdata = NULL, interval = "confidence")

lines(x, predict(lm1, newdata = NULL, interval = "predict")[,2], col = 4)
lines(x, predict(lm1, newdata = NULL, interval = "predict")[,3], col = 4)

lines(x, predict(lm1, newdata = NULL, interval = "confidence")[,2], col = 3)
lines(x, predict(lm1, newdata = NULL, interval = "confidence")[,3], col = 3)



y <- 2 +