## R Ue6.5 dats_01_22 Modell-Spezifikation Konsumgleichung 
## 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/Chp6/EXERCISES/")
setwd("C:/Users/hoersaal/Downloads/")

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

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)

## cpr1 ist das gelaggte cpr, cpr(-1) bzw L(cpr_t) = cpr_(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 ist das gelaggte wlr, wlr(-1) bzw L(wlr_t)
wlr1    <- ts(c(rep(NA,5), wlr[1:(n_cpr-1)]), start=1976, end=2021)  
yhr1_h  <- c(NA, yhr[1:(n_cpr-1)])     # yhr_(t-1): muessen den Lag selbst erzeugen
yhr1    <- ts(yhr1_h, start=1976, end=2021)
  
ix      <- c(1:n_cpr)

# df enthaelt die Variablen ix (Index) ,year, yhr, yhr (yhr(t-1)), 
#    const,cpr,cpr1 (cpr(-1)=L(cpr_t)) ,wlr1 (wlr(-1)=L(wlr_t))
df <- as.data.frame(cbind(ix,year,yhr,yhr1,const,cpr,cpr1,wlr1))  # ~~ in K5Ue3
head(df)
tail(df)

## plot von cpr und yhr                     # wie in K5Ue3 zum Anschauen
#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") )

## Es sollen 5 verschiedene Modelle durchprobiert und verglichen werden
## Schaetzung hier jeweils für 1985 bis 2021 
## Sie koennen auch eine andere Periode wählen. 
## Sie muss aber für alle Modelle dieselbe sein.
n_mod <- 5      
mod_ICs <- rep(NA,n_mod*2);  dim(mod_ICs) <- c(n_mod,2)

## Modell 1: cpr = a + u     (Modell nur mit Konstanter)
mod_1 <- lm(cpr ~ 1, data=df[df$year %in% c(1985:2021),])  
summary(mod_1)
mod_ICs[1,1] <- AIC(mod_1);   mod_ICs[1,2] <- BIC(mod_1);  
cat("Modell 1: ","AIC= ", AIC(mod_1), "BIC= ", BIC(mod_1), "\n")

## ... Probieren sie 3 weitere Modellspezifikationen

## Modell 2: cpr = a + ... + u                                   ## HIER IHRE MODELLSPEZIFIKATION
mod_2 <- lm(cpr ~ ... + 1, data=df[df$year %in% c(1985:2021),])  ## HIER IHRE MODELLSPEZIFIKATION
summary(mod_2)
mod_ICs[2,1] <- AIC(mod_2);   mod_ICs[2,2] <- BIC(mod_2);  
cat("Modell 2: ","AIC= ", AIC(mod_2), "BIC= ", BIC(mod_2), "\n")

## Modell 3: cpr = a + ... + u                                   ## HIER IHRE MODELLSPEZIFIKATION
mod_3 <- lm(cpr ~ ... + 1, data=df[df$year %in% c(1985:2021),])  ## HIER IHRE MODELLSPEZIFIKATION
summary(mod_3)
mod_ICs[3,1] <- AIC(mod_3);   mod_ICs[3,2] <- BIC(mod_3);  
cat("Modell 3: ","AIC= ", AIC(mod_3), "BIC= ", BIC(mod_3), "\n")

## Modell 4: cpr = a + ... + u                                   ## HIER IHRE MODELLSPEZIFIKATION
mod_4 <- lm(cpr ~ ... + 1, data=df[df$year %in% c(1985:2021),])  ## HIER IHRE MODELLSPEZIFIKATION
summary(mod_4)
mod_ICs[4,1] <- AIC(mod_4);   mod_ICs[4,2] <- BIC(mod_4);  
cat("Modell 4: ","AIC= ", AIC(mod_4), "BIC= ", BIC(mod_4), "\n")

##

## Modell 5: cpr = a + b yhr + c wlr1 + d cpr1 + e Trend + f yhr1 + u 
##                                                   (Modell mit allen Variablen) 
mod_5 <- lm(cpr ~ yhr + wlr1 + cpr1 + year + yhr1 + 1, data=df[df$year %in% c(1985:2021),])  
summary(mod_5) 
mod_ICs[5,1] <- AIC(mod_5);   mod_ICs[5,2] <- BIC(mod_5);  
cat("Modell 5: ","AIC= ", AIC(mod_5), "BIC= ", BIC(mod_5), "\n")


## Welche Modellspezifikation scheint die geeignetste zu sein?
cat("         AIC      BIC","\n")
mod_ICs