################################################################################
# Univariate Zeitreihenanalyse
################################################################################

library(urca)
library(tseries)
library(vars)
library(MASS) # for dataset

data(EuStockMarkets)


dax <- EuStockMarkets[, "DAX"]
smi <- EuStockMarkets[, "SMI"]
cac <- EuStockMarkets[, "CAC"]
ftse <- EuStockMarkets[, "FTSE"]

# 0. Visualisierung 

plot(EuStockMarkets)



# 1. Einheitswurzeltests
 
# a) ADF Test

adf.y <-  ur.df(y,type="drift",selectlags="AIC")
summary(adf.y)    # Nicht ablehnen

adf.dy <-  ur.df(dy,selectlags="AIC")
summary(adf.dy)   # Ablehnen

# Ergebnis:
# Log GDP ist I(1) (enth?lt eine Einheitswurzel),
# W'rate des GDP ist I(0)

# b) KPSS Test

kpss.test(y) # Ablehnen: Nichtstation?r
kpss.test(dy) # Auch ablehnen: Nichtstation?r 

# c) weitere:
# ur.pp (Phillips & Perron Unit Root Test)
# ur.sp (Schmidt & Phillips Unit Root Test)
# ur.ers (Elliott, Rothenberg & Stock Unit Root Test)
# ur.za (Zivot & Andrews Unit Root Test)


# 2. Autokorrelationen und partielle ACF
acf(dy,lag.max=60)
pacf(dy,lag.max=60)

# 3. Sch?tzen AR Modell
(ar.4 <- ar(dy,order=4,aic=FALSE,method="ols"))

# Elemente des Sch?tzobjekts
summary(ar.4)

# AR-Parameter und Standardfehler
cbind(ar.4$ar,
ar.4$asy.se.coef$ar)

# Residuen Plotten und ACF anschauen
u.ar4 <- na.omit(ar.4$resid) 
plot.ts(u.ar4)


# Lagl?ngenwahl mit AIC (sehr viele Parameter zu sch?tzen!)
(ar.aic <- ar(dy,order.max=40,aic=TRUE,method="ols"))
 
 
# 4. Sch?tzen eines ARMA Prozesses 

(arma.33 <- arima0(dy,order=c(3,0,3),method="ML" )) 
# AR(3)-Term ist insignifikant, AIC wird hier kleiner:
(arma.23 <- arima0(dy,order=c(2,0,3),method="ML" ))  

# Residuen generieren
u.23 <- arma.23$resid


# 5. Autokorrelation in Residuen?

acf(u.23,lag.max=60)
acf(u.ar4,lag.max=60)

Box.test(u.23,lag=20,type="Ljung-Box",fitdf=5)
Box.test(u.23,lag=20,type="Box-Pierce",fitdf=5)

Box.test(u.ar4,lag=20,type="Ljung-Box",fitdf=4)
Box.test(u.ar4,lag=20,type="Box-Pierce",fitdf=4)    
# Keine Signifikante Autokorrelation

# 6. ARCH Effekte?

acf(u.23^2)
acf(u.ar4^2)

order <- 3
U_sq <- as.data.frame(embed(u.23^2,order+1) )
names(U_sq) <- c("u.sq",paste("u.sq.l",1:order,sep=""))

arch.eq <- lm(u.sq~.,data=U_sq)
summary(arch.eq)

Arch.LM <- summary(arch.eq)$r.squared*nrow(U_sq)
Arch.LM
pchisq(Arch.LM,df=order,lower.tail=FALSE)  
# Ablehnen auf 5%-Level: Signifikante bedingte Heteroskedastie
 
# 7. Normalit?t?

hist(u.23,probability=TRUE,breaks=25)
x <-  seq(min(u.23),max(u.23),0.001)
lines(x,dnorm(x,sd=sd(u.23),mean=mean(u.23)),col="darkred" ) 

jarque.bera.test(u.23) # Wird klar abgelehnt: Keine Normalit?t der Fehler

# 8. Prognose

n.ahead <- 24
pr.ar <- predict(ar.4,n.ahead=n.ahead)
pr.arma <- predict(arma.23,n.ahead=n.ahead)

plot(c(dy,pr.ar$pred),type="n",ylab="GDP Growth, Prediction")
lines(dy)
lines((length(dy)+1):(length(dy)+n.ahead),pr.ar$pred,col=2)
lines((length(dy)+1):(length(dy)+n.ahead),pr.ar$pred-2*pr.ar$se,col=2,lty=2)
lines((length(dy)+1):(length(dy)+n.ahead),pr.ar$pred+2*pr.ar$se,col=2,lty=2)


plot(c(dy,pr.ar$pred),type="n",ylab="GDP Growth, Prediction")
lines(dy)
lines((length(dy)+1):(length(dy)+n.ahead),pr.arma$pred,col=3)
lines((length(dy)+1):(length(dy)+n.ahead),pr.arma$pred-2*pr.ar$se,col=3,lty=2)
lines((length(dy)+1):(length(dy)+n.ahead),pr.arma$pred+2*pr.ar$se,col=3,lty=2)

 

################################################################################
# Multivariate Zeitreihen (VARS)
################################################################################

GDP <- read.table("GDP.csv", sep = ",", header = TRUE)
names(GDP)[2] <- "gdp"
CPI <- read.table("CPIAUCSL.csv", sep = ",", header = TRUE)
names(CPI)[2] <- "cpi"

GDPCPI <-  merge(GDP,CPI,by.y="DATE",by.x="DATE")
GDPCPI$DATE <- as.Date(GDPCPI$DATE)
GDPCPI$gdp.real <- GDPCPI$gdp/GDPCPI$cpi
GDPCPI$dp <- c(NA,diff(log(GDPCPI$cpi)))
GDPCPI$dy <- c(NA,diff(log(GDPCPI$gdp.real)))
GDPCPI <- subset(GDPCPI,select=c("dy","dp"),subset=DATE>=as.Date("1954-01-01"))

plot.ts(GDPCPI)
attach(GDPCPI)

summary(ur.df(dp,type="drift",selectlags="AIC"))
summary(ur.df(dy,type="drift",selectlags="AIC"))

varaic <- VAR(cbind(dy,dp),lag.max=12,ic="AIC")
summary(varaic)

serial.test(varaic,type="BG")
arch.test(varaic)       # Abgelehnt
normality.test(varaic)  # Abgelehnt


Acoef(varaic)
predict(varaic,n.ahead=12)
ImpResp1 <- irf(varaic,n.ahead=45,ortho=FALSE,cumulative=TRUE)
plot(varaic)











################################################################################
# Kointegrationsanalyse
################################################################################

# 0.1 Daten Einlesen: CPI und 3-Monats T-Bill
rsetwd("c:/dokumente und einstellungen/roland/eigene dateien/lehre/r kurs/datens?tze")
CPI <- read.table("CPIAUCSL.csv", sep = ",", header = TRUE)
TBILL <-  read.table("TB3MS.csv", sep = ",", header = TRUE)

# 0.2 Generiere Jahresinflationsrate; Quartalsdaten
US_Fisher <- data.frame( date=as.Date(CPI$DATE), tbill=TBILL$VALUE, lcpi=log(CPI$VALUE) )
US_Fisher$infl <- c(rep(NA,12),diff(US_Fisher$lcpi,12))*100
US_Fisher <- subset(US_Fisher,
    subset=c(TRUE,FALSE,FALSE) & date>=as.Date("1954-01-01"))

plot.ts(US_Fisher[,-1])
attach(US_Fisher)

# 1. Einheitswurzeltests
      
adf.tbill <-  ur.df(tbill,type="drift",selectlags="AIC")
summary(adf.tbill)  # Nicht ablehnen
adf.dtbill <-  ur.df(diff(tbill),selectlags="AIC")
summary(adf.dtbill) # Ablehnen

adf.infl <-  ur.df(infl,type="drift",selectlags="AIC")
summary(adf.infl)  # Nicht ablehnen
adf.dinfl <-  ur.df(diff(infl),selectlags="AIC")
summary(adf.dinfl) # Ablehnen

# Also: Beide Reihen I(1)

# 2.1 OLS Regression (Statische Regression)

fisher.ols <- lm(tbill~infl)
summary(fisher.ols)


# 2.2 Test auf Kointegration (Residuenbasiert)

ec.ols <- fisher.ols$resid
plot.ts(ec.ols)
summary(ur.df(ec.ols))  # Achtung: Kritische Werte sind falsch!!!
 

PO.test <- ca.po(cbind(tbill,infl),demean="const",type="Pz")
summary(PO.test) 

# Ergebnis:  Kointegration (auf 5% Signifikanzniveau)

# 3. Johansen Analyse


# 3.1. VAR Sch?tzung und Wahl der Lag-Ordnung

fisher.var <- VAR(cbind(tbill,infl),ic="SC",lag.max=12)
summary(fisher.var)
p <- fisher.var$p

# Problematisch:  
serial.test(fisher.var,type="BG", lags.bg=p+6)
arch.test(fisher.var, lags.multi=p+6)

# 3.2. Test auf den Kointegrationsrang

fisher.MaxEigen <- ca.jo(cbind(tbill,infl),K=p,ecdet="const",spec="transitory",
  type="eigen")
summary(fisher.MaxEigen)


fisher.Trace <- ca.jo(cbind(tbill,infl),K=p,ecdet="const",spec="transitory",
  type="trace")
summary(fisher.Trace)

# Ergebnis: 
# H_0: Keine Kointegrationsbeziehung abgelehnt nur bei Trace Test
# und nur auf 10% Niveau
# H_0: Genau eine KI-Beziehung zumeist nicht abgelehnt
# W?hle EINE Kointegrationsbeziehung

# 3.3.  Sch?tzung des VECM (r=1: eine CI Beziehung)

fisher.VECM <- cajorls(fisher.Trace, r=1)
fisher.VECM$beta
summary(fisher.VECM$rlm)

# Koeffizienten in Levels
vec2var(fisher.Trace,r=1)


### Ergebnis: Existenz von Kointegration ist nicht eindeutig
# Koeffizienten (Kointegrationsvektor) unterscheidet sich deutlich zwischen
# OLS
c(1,-fisher.ols$coef)
# Johansen
fisher.VECM$beta[c(1,3,2)]

#