## R Ue5.1 Eigene Daten, Prognose, MSE
##         Alle einzelnen Schritte sind zu beachten!
##
#### Geben sie IHRE Daten ein. ####
x     <- c(1,2,3,4,5)
y     <- c(2,1,4,3,5)

const <- c(1,1,1,1,1)
dat <- as.data.frame(cbind(y,const,x))

## Streudiagramm anschauen
plot(dat$x, dat$y)

## Modell  y_i = const + beta x_i + u_i
# Regressionsgerade y_i-hat = const-hat + b * x_i      y-hat steht fuer y-Dach
reg <- lm(y ~ x + 1, data=dat)
summary(reg)

## (a) Berechnen sie haendisch die in-sample Prognose für i=3,5
y_hat <- rep(NA,5)           # Speicherplatz für 5 in-sample Prognosen anlegen

## y_3-hat = const_hat + b * x_3         
y_hat[3] <- reg$coefficients[1] + reg$coefficients[2]*dat$x[3]

## y_5-hat = const_hat * const + b * x_5
y_hat[5] <- reg$coefficients[1] + reg$coefficients[2]*dat$x[5]

## Vergleichen sie y_3-hat mit y_3 und y_5-hat mit y_5
## Sie sehen die in-sample Prognose ist nichts anderes als 
## die fitted Regressionsgerade: y_i = y-hat_i + residual_i
cbind(y, y_hat, reg$residuals)

## Berechnen sie den MSE nur fuer diese beiden Zeitpunkte
## s. Kap 5 Folie 21. Hier gilt T=2 
## Der mittlere quadratische Fehler:
MSE_3u5 <- ( (y[3]-y_hat[3])^2 + (y[5]-y_hat[5])^2 )/2
as.numeric(MSE_3u5)
## Die Wurzel aus MSE, RMSE, sollte ungefähr dem Standardfehler der Residuen 
## aus der Schätzung entsprechen, wenn n groß ist. Ist das der Fall?
as.numeric(MSE_3u5**0.5)


## (b) MSE fuer alle 5 Beobachtungen
## Die in-sample Prognosen sind die fitted Regressionsgerade
y_hat <- dat$y - reg$residuals
cbind(y_hat, reg$fitted.values)

MSE_all <- t(y_hat - dat$y) %*% (y_hat - dat$y) / 5
as.numeric(MSE_all)
as.numeric(MSE_all**0.5)

## (c) Prognose fuer Beobachtung 6 - out of sample
####  Waehlen sie IHRE Werte fuer x_6 und ein y_6 ####
x_6 <- 7;   y_6 <- 9
y_6_hat <- reg$coefficients[1] + reg$coefficients[2]*x_6

## (d) Vergleichen sie y_6_hat mit y_6
cat(y_6,y_6_hat,"\n")
MSE_6 <- (y_6 - y_6_hat)^2
as.numeric(MSE_6)          
## Vergleichen sie die 3 MSE-Werte MSE_3u5, MSE_all, MSE_6
cat(MSE_3u5, MSE_all, MSE_6,"\n")

## (e) Auffuellen der gespeicherten Daten mit den 
##     Beobachtungen x_6, ..., x_10                
## (e.1) mit aehnlichen x-Werten
x2 <- c(x, 0.5 , 2.5 , 3.5 , 1.5 , 5.5 )  # die neuen Beobachtungen werden an 
                                  # die alten anghaengt

## Prognosen für y fuer i=6,...,10
y_hat2 <- reg$coefficients[1] + reg$coefficients[2]*x2[6:10]

# anschauen
cbind(y_hat2, x2[6:10])

## Progonseintervall 
## Berechnung der Varianz des Prognosefehlers Kap 5 Folie 16
Sxx <- as.numeric( t(x) %*% x )
n   <- 5      # 5 Beobachtungen wurden fuer die Regression verwendet
# sigma der Residuen ist in summary(reg)$sigma abgelegt.
V_pf <- summary(reg)$sigma^2 *( 1 + 1/n + x2[6:10] / Sxx )
se_pf <- V_pf^(0.5)
## Die Laenge des Prognosintervalls haengt vom x-Wert ab!
## Die se's der Prognose fuer i=6,..,10 sind
se_pf

## Prognosintervalle, PIs, fuer i=6,...,10
PIs_links  <- y_hat2 - 1.96*se_pf
PIs_rechts <- y_hat2 + 1.96*se_pf
PIs <- cbind(PIs_links,PIs_rechts)
PIs
## Laenge der PIs
PIs_Laenge <- cbind(PIs_links,PIs_rechts, PIs_rechts - PIs_links)
PIs_Laenge

## (e.2) wie (e.1) oben mit untypischen x-Werten im Prognoseintervall
x2 <- c(x, 5 , 25 , 50 , 100 , 500 )  # die neuen Beobachtungen werden an 
                           # die alten anghaengt...
...