## R Ue5.3 dats_01_22 Prognose des privaten Konsums 
## Bezeichnungen siehe dat_01_22.docx 
## cpn, cpr ... Privater Konsum nominell, real (Mio Euro)
## yhn, yhr ... Haushaltseinkommen nom, real (Mio Euro)
## wln      ... Geldvermoegen nom (Mio Euro ?)
## pcd      ... Konsumdeflator 2015=100

#setwd("C:/MH/WU/LV/OEKONOMETRIE_BA/Oe1_WS23/Chp5/EXERCISES/")
setwd("C:/Users/hoersaal/Downloads/")

## BasicStatistics_R.txt neu herunterladen
source("BasicStatistics_R.txt")   
# Daten einlesen (nicht ganz einfach)
dat <- read.table("dats_01_22.csv", sep=";", dec=",", header=TRUE, na.strings = "#NV", fill = TRUE, comment.char="")

## Daten anschauen
#head(dat)
#tail(dat)
#dim(dat)   # Zeitreihen n=47 k=9 (Variablen)  von 1976 bis 2021
names(dat)  # Namen der Variablen

## Datenaufbereitung, 
cpr <- ts(dat$cpr, start=1976, end=2021)
cpn <- ts(dat$cpn, start=1976, end=2021)
yhr <- ts(dat$yhr, start=1976, end=2021)
yhn <- ts(dat$yhn, start=1976, end=2021)
pcd <- ts(dat$pcd, start=1976, end=2021)
n_cpr <- length(cpr)
const <- ts(rep(1,n_cpr), start=1976, end=2021)
year <- ts(c(1976:(1976 + n_cpr - 1)), start=1976, end=2021)

## Erzeugen der gelaggten Variablen cpr1 = cpt_(t-1) und wlr1 = wlr_(t-1)
cpr1_h  <- c(NA, cpr[1:(n_cpr-1)])     # cpr_(t-1): muessen den Lag selbst erzeugen
cpr1    <- ts(cpr1_h, start=1976, end=2021)
wlr     <- dat$wln/dat$pcd
wlr1    <- ts(c(rep(NA,5), wlr[1:(n_cpr-1)]), start=1976, end=2021)    
ix      <- c(1:n_cpr)

# gemeinsames neues data.frame df
df <- as.data.frame(cbind(ix,year,yhr,const,cpr,cpr1,wlr1))
head(df)
tail(df)

## plot von cpr und yhr
mincpr <- min(cpr, na.rm=TRUE); maxyhr <- max(yhr, na.rm=TRUE)
plot.ts(yhr, ylab="yhr, cpr", ylim=c(mincpr,maxyhr) )
lines(cpr, col=2)
legend("topleft", legend = c("yhr","cpr"), text.col = c("black","red") )

## HIER
## ENTWEDER Modell 1 ODER Modell 2 hier auswaehlen
##-------------------------------------------------

## Modell 1: cpr = a + b yhr + c cpr(-1) + u 
## Schaetzung für 1977(76) bis 2010
mod <- lm(cpr ~ yhr + cpr1 + 1, data=df[df$year %in% c(1977:2010),])  
summary(mod) 
#
## Abfolge von 1-Schritt-Prognosen, ex-post, die zukuenftigen x-Wrte sind bekannt (!)
## wobei stets das Modell für 1977 -- 2010 verwendet wird
coef <- mod$coefficients
fcp  <- df$ix[df$year %in% c(2011:2020)]                  # fcp ... forecast period
cpr_fc <- coef[1] + coef[2]*df$yhr[fcp] + coef[3]*df$cpr1[fcp]        
lab_main <- "Modell 1"
##-------------------------------------------------

## Modell 2: cpr = a + b yhr + c wln(-1)/pcd(-1) + u 
mod <- lm(cpr ~ yhr + wlr1 + 1, data=df[df$year %in% c(1981:2010),])  
summary(mod)     
#
## Abfolge von 1-Schritt-Prognosen, ex-post, die zukuenftigen x-Wrte sind bekannt (!)
## wobei stets das Modell für 1977 -- 2010 verwendet wird
coef <- mod$coefficients
fcp  <- df$ix[df$year %in% c(2011:2020)]                  # fcp ... forecast period
cpr_fc <- coef[1] + coef[2]*df$yhr[fcp] + coef[3]*df$wlr1[fcp]        
lab_main <- "Modell 2"    
##-------------------------------------------------

## (a) "In-sample Prognose"             
# Plot v observed fitted, residuals
plot_ser_fit_res.ts(df$cpr,mod,"cpr", 1977, lab_main) 


## (b) ex post Prognose fuer 2011 bis 2020   
## Um welche Prognosen handelt es sich? 1-Schritt, 2-Schritt, ... Prognosen

## plot von cpr und cpr_fc  1977 - 2020
##
cpr_fc_h <- ts( c(rep(NA,2010-1976),cpr_fc), start=1977, end=2020)  
cpr_h    <- ts( df$cpr[df$year %in% c(1977:2020) ], start=1977, end=2020)  
year_h   <- c(1977:2020) ;   ix_h <- c(1: length(year_h))
ix_2010  <- ix[year_h == 2010];  cpr_fc_h[ix_2010] <- cpr_h[ix_2010]      

#cbind(year_h, cpr_h, cpr_fc_h)   # check

mincpr <- min(cpr_h, na.rm=TRUE);  maxfc <- max(cpr_fc_h, na.rm=TRUE)
plot.ts(cpr_h, ylab="cpr", ylim=c(mincpr,maxfc), main=lab_main)
lines(cpr_fc_h, col=2)
legend("topleft", legend = c("cpr","cpr_fc"), text.col = c("black","red") )

## plot von cpr und cpr_fc  nur fuer 2000 - 2020
##
ix_0020 <- ix[year_h %in% c(2000:2020)]
cpr_fc_h1 <- ts(cpr_fc_h[ix_0020], start=2000, end=2020)  
cpr_h1    <- ts(cpr_h[ix_0020], start=2000, end=2020)              

mincpr <- min(cpr_h1, na.rm=TRUE); maxfc <- max(cpr_fc_h1, na.rm=TRUE)
plot.ts(cpr_h1, ylab="cpr", ylim=c(mincpr,maxfc) , main=lab_main)
lines(cpr_fc_h1, col=2)
legend("topleft", legend = c("cpr","cpr_fc"), text.col = c("black","red") )


## Prognoseintervall
## Berechnung der Varianz des Prognosefehlers Kap 5 Folie 16 nur 
## fuer das bivariate Modell, 
## fuer multiple Regression siehe zB Verbeek p.52
##
#  The prediction error has the variance for a single fc-step
#  V[y0_hat - y0] = sigma^2[ 1 + x0'(X'X)^(-1)x0 ]
sigma <- summary(mod)$sigma
X     <- as.matrix( df[df$year %in% c(1977:2010),c("const","cpr","cpr1")] )
XtXinv  <-  solve( t(X) %*% X )          # berechnet die Inverse von (X'X)

se_fc <- rep(NA,10);   
for (i in 2011:2020) {  j <- i - 2010
   x0  <- as.numeric( df[df$year %in% c(i),c("const","cpr","cpr1")] )
   se_fc[j] <- sigma * ( 1 + t(x0) %*% XtXinv %*% x0)**(0.5) 
                     }
PI_links  <- cpr_fc - 1.96*se_fc
PI_rechts <- cpr_fc + 1.96*se_fc
cbind(PI_links,PI_rechts)


## (c) plot von cpr und cpr_fc  mit Prognosintervall fuer 2000 - 2020
##
PI_links_h1  <- ts(PI_links, start=2011, end=2020)
PI_rechts_h1 <- ts(PI_rechts, start=2011, end=2020)            

mincpr <- min(cpr_h1, na.rm=TRUE); maxfc <- max(PI_rechts_h1, na.rm=TRUE)
plot.ts(cpr_h1, ylab="cpr", ylim=c(mincpr,maxfc), main= lab_main )
lines(cpr_fc_h1, col=2)
lines(PI_links_h1, col="blue", type="l")
lines(PI_rechts_h1, col="blue", type="l")
legend("topleft", legend = c("cpr","cpr_fc"), text.col = c("black","red") )


## (d) RMSE, MAPE
diff <- cpr_fc - df$cpr[df$year %in% c(2011:2020)] 
MSE <- sum(diff^2) /10
RMSE <- MSE^0.5
#
MAPE <- 100*sum( abs(diff/df$cpr[df$year %in% c(2011:2020)]) )/10
#
cat(RMSE,MAPE, "\n")    
## Modell 1: 4969.03  1.805646
## Modell 2: 5169.426 1.710254


## (e) Wie wuerde eine ex ante Progose durchzufuehren sein?