################################################################################ # 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)] #