## R Ue2.1 okuns law
## Reproduktion der Folien aus Kap 2

#setwd("C:/MH/WU/LV/OEKONOMETRIE_BA/Oe1_WS23/Chp2/EXERCISES/")
setwd("C:/Users/hoersaal/Downloads/")
# Daten einlesen (nicht ganz einfach)
dat <- read.table("okuns_law.csv", sep=";", dec=",", header=TRUE, na.strings = "#NV", fill = TRUE, comment.char="")
## Daten anschauen
dim(dat)     # Es sind Daten: 66 Zeitpunkte, 2 Variable (+ Zeitachse)
             # 0-te Zeile sind die names
head(dat,4)     
tail(dat)

# 1950 - 2015
bipr.l <- ts(dat$bipr, start=1950)  # .l steht für lang
ur.l   <- ts(dat$ur  , start=1950)

op <- par(mfrow=c(2,1))
par(mar=c(2, 4, 4, 2) +0.1) 
plot.ts(bipr.l)
#
par(mar=c(5, 4, 1, 2) +0.1) 
plot.ts(ur.l)
par(op)

# 1976 - 2010 Datenfenster 

bipr <- window(bipr.l, start=1976, end=2010)
ur   <- window(ur.l,   start=1976, end=2010)

op <- par(mfrow=c(2,1))
par(mar=c(2, 4, 4, 2) +0.1) 
plot.ts(bipr)
#
par(mar=c(5, 4, 1, 2) +0.1) 
plot.ts(ur)
par(op)

## PLOT 1 LV  -  bipr und ur
# 2 Reihen in einem plot mit unterschiedlichen Skalen (nicht ganz einfach)
par(mar=c(5, 4, 4, 5) +0.1) # Groesse des Plots, 
   # +0.1 ist der zusaetzliche Platz rechts fuer "ur"
plot.ts(bipr, axes=FALSE)
axis(2)
par(new=TRUE)   # naechster Plot wird daruebergelegt
plot.ts(ur, col="blue", axes=FALSE, ylab="")
axis(4)
mtext("ur",side=4,col="blue",line=2) 
axis(1)
title("bipr und ur")
legend("bottomright", legend=c("bipr","ur"), text.col=c("black","blue")) 

## Es schaut so aus, als gaebe es einen positiven Zusammenhang 
## dem realen BIP und der Arbeitslosenrate !?

## Wirtschaftswachstum in %
rho <- ( log(bipr) - log(lag(bipr,-1)) )*100
# Veraenderung der Arbeitslosenrate
dur <- ur - lag(ur,-1) 

## PLOT 2 LV  - rho und dur    (copy & paste von oben)
par(mar=c(5, 4, 4, 5) +0.1) 
# R: Groesse Raender des Plots. default: mar=c(5, 4, 4, 2) +0.1

plot.ts(rho, axes=FALSE)
axis(2)
par(new=TRUE)   # naechster Plot wird daruebergelegt
plot.ts(dur, col="blue", yaxt="n", ylab="")
axis(4)
mtext("dur",side=4,col="blue",line=2) 
axis(1)
title("rho und dur")
legend("bottomleft", legend=c("rho","dur"), text.col=c("black","blue")) 

## Ein negativer Zusammenhang scheint zu ueberwiegen. 

## Wir eklaeren dur durch rho. Daher rho x dur.
## PLOT 3 LV  -  Streudiagramm rho x dur
plot(as.numeric(rho),as.numeric(dur), type="p", xlab="rho", ylab="dur")
# Für das Streudiagramm muessen wir die Eigenschaft Zeitreihe 
# entfernen.
## Gibt es einen positiven oder negativen Zusammenhang?

## OLS Schaetzung von dur = a + b rho + u
reg <- lm(dur ~ rho + 1)  # +1 zeigt an, dass die Konstante mit dabei ist
summary(reg)              # Zeigt die wichtigsten Statistiken der Schaetzung
## Vergleichen sie R Output mit dem Ergebnis auf Folie K2.14    

## In R wird n, die Anz d Beobachtungen, nicht angegeben, sondern die 
## degrees of freedom, die Freiheitsgrade, df. 
## df + Anz der geschaetzen Paramter = n
## 32 + 2 = 34
## Statt RSS wird der Residual standard error ausgegeben.

reg$coefficients       # Vektor der geschaetzten Koeffizienten von reg
a <- reg$coefficients[1]
b <- reg$coefficients[2]
# das geschätzte Modell ist dur-hat = a + b*rho

## PLOT 4 LV  -  Streudiagramm rho x dur mit Regressionsgeraden
plot(as.numeric(rho),as.numeric(dur), type="p", xlab="rho", ylab="dur")
abline(a,b, col="blue")       # zeichnet die Gerade  y=a+b*x
   
## PLOT 5 LV  - Datenreihe (dur) , Fitted, Residuen
# fitted und residuals beginnen erst 1977, da dur
# durch die Differenzenbildung die 1-te Beobachtung verliert
observed <- ts(dur, start=1977)
fitted   <- ts(reg$fitted, start=1977)
res      <- ts(reg$residuals, start=1977)
n <- length(res)
null   <- ts(rep(0,n), start=1977)
#
source("BasicStatistics_R.txt")
plot_ser_fit_res.ts(dur,reg,"dur",1977)