# R 3.4.1
# Ex4_R.txt
## Dynamic Dependency: US stock to other stock indices
## wrt returns or volatility

rm(list=ls())
library("tseries")  
library("zoo")
library("vars")     # help(package="urca")

# OK ret((FTSE, DJA/DJI)), ret(DAX, DJA)
# NOT ret(FTSE, NDX/GSPC), ret(BSESN, DJA)

## read data
start0 <- "1998-10-01" ;   end0   <- "2010-07-30"      # choose period  
require(quantmod) # Finanzdaten download

spc <- new.env()
dax <- new.env()

setDefaults(getSymbols,src="yahoo")
getSymbols("^GSPC", env = spc, from = as.Date(start0), to = as.Date(end0))
getSymbols("^GDAXI", env = dax, from = as.Date(start0), to = as.Date(end0))
# read in series are of xts type,  plot is different for xts and ts

x.c <-as.zoo( spc$GSPC[,4]);       y.c <- as.zoo(dax$GDAXI[,4]);                                        
###
nx <- length(x.c); ny <- length(y.c);                       

## plot data
xy.c <- na.omit(merge(x.c, y.c,all=FALSE))     
x.c <- xy.c[,1];  y.c <- xy.c[,2]                       
plot( xy.c, main = "^GSPC ,  ^GDAXI")    # plot series
plot(x.c, y.c, main = "^GSPC x  ^GDAXI")   # plot scatter
plot(x.c, y.c, main = "^GSPC x  ^GDAXI", type="l")   # plot scatter

# not used
#x <-        log(x.c) ;  
#y <- zoo( log(y.c),  order.by= index(x.c));  # fewer trading days in x.c
#xy <- cbind(x,y);

x <- log(x.c) ;   y <- log(y.c)
xy <- merge.zoo(x,y, all=FALSE)    
plot(xy, main = "log(^GSPC), log(^GDAXI)")   # plot series
plot(xy[,1], xy[,2], main = "^GSPC x  ^GDAXI")   # plot scatter
plot(xy[,1], xy[,2], main = "^GSPC x  ^GDAXI", type="l")   # plot scatter

## omit period shortly after 9/11
#start0_1 = "1998-10-01"; end0_1 = "2001-09-10"; 
#start0_2 = "2001-10-01"; end0_2 = "2010-07-30";    

## choose your subperiod
start0_1 = "1998-10-01"; end0_1 = "2001-09-10";  
## ci  e.g.: 00-01, 09-10, 01-02, 02-03, 03-04, 04-05, 05-06,
start0_2 = "2005-10-01"; end0_2 = "2006-07-30";     

 i.x   <- index(xy)
 i.x1 <- i.x[ as.Date(start0_1) <= i.x & i.x <= as.Date(end0_1) ]
 i.x2 <- i.x[ as.Date(start0_2) <= i.x & i.x <= as.Date(end0_2) ]

x1  <- zoo(x[i.x1])
## alternatively: 
## x1  <- zoo(xts(window(x, start = start0_1, end = end0_1)))
## colnames(x1) <- colnames(x)

y1  <- zoo(y[i.x1])
xy1 <- merge.zoo(x1,y1, all=FALSE)   
mainlab1 <- paste("log(^GSPC), log(^GDAXI):", start0_1,"to",end0_1)
plot(xy1, main = mainlab1 )   # plot series
plot(x1,y1, main = mainlab1 )   # plot scatter
plot(x1,y1, main = mainlab1, type="l" )   # plot scatter

x2  <- zoo(x[i.x2])
y2  <- zoo(y[i.x2])
xy2 <- merge.zoo(x2,y2, all=FALSE)   
mainlab2 <- paste("log(^GSPC), log(^GDAXI):", start0_2,"to",end0_2)
plot(xy2, main = mainlab2 )   # plot series
plot(x2,y2, main = mainlab2 )   # plot scatter
plot(x2,y2, main = mainlab2, type="l" )   # plot scatter

nxy1 <- length(xy1[,1]);   nxy2 <- length(xy2[,1]);   # no of observations 

##Dickey-Fuller Unit root test: level and 1st difference
# tau* test for  H0:                    (1-rho) = 0  in various models
# phi1  test for H0: const=0 & (1-rho) = 0  in model w const/drift
# phi3  test for H0: no trend & (1-rho) = 0  in trend model
## type = "none" or "drift" or "trend"
## selectlags= "Fixed" or "AIC", ... 
## relevant statistic: tau

x.df <- ur.df(x, type="drift", lags=5, selectlags="Fixed")     
summary(x.df)
y.df <- ur.df(y, type="drift", lags=5, selectlags="BIC")
summary(y.df)

dx.df <- ur.df(diff(x), type="drift", lags=5, selectlags="Fixed")     
summary(dx.df)
dy.df <- ur.df(diff(y), type="drift", lags=5, selectlags="BIC")
summary(dy.df)

# 1
x1.df <- ur.df(x1, type="drift", lags=0, selectlags="Fixed")    
summary(x1.df)
y1.df <- ur.df(y1, type="drift", lags=0, selectlags="Fixed")
summary(y1.df)

dx1.df <- ur.df(diff(x1), type="drift", lags=5, selectlags="BIC")    
summary(dx1.df)
dy1.df <- ur.df(diff(y1), type="drift", lags=5, selectlags="BIC")
summary(dy1.df)

# 2
x2.df <- ur.df(x2, type="drift", lags=0, selectlags="BIC")    
summary(x2.df)
y2.df <- ur.df(y2, type="drift", lags=0, selectlags="BIC")
summary(y2.df)

dx2.df <- ur.df(diff(x2), type="drift", lags=5, selectlags="BIC")    
summary(dx2.df)
dy2.df <- ur.df(diff(y2), type="drift", lags=5, selectlags="BIC")
summary(dy2.df)


##Cointegration  - test
# spec="transitory"  is our specification: Delta x_t = ... Pi x_(t-1) ...
# const in Delta x always included
# help("ca.jo-class")   reference: ci_test@Z0
## type = c("eigen", "trace"),  ... tests
## ecdet = c("none", "const", "trend")   ... deterministic part in EC part
## (K-1) no of lags of Delta x included
# r in output is our m

xy.cit <- ca.jo(xy, ecdet = "none", type="eigen", K=4, spec="transitory" )
summary(xy.cit)

xy1.cit <- ca.jo(xy1, ecdet = "none", type="eigen", K=2, spec="transitory" )
summary(xy1.cit)

xy2.cit <- ca.jo(xy2, ecdet = "none", type="eigen", K=2, spec="transitory" )
summary(xy2.cit)
plot(xy2.cit)
plotres(xy2.cit)
 
####

## CI VAR
xy2.civar <- cajorls(xy2.cit, r = 1)    # r (in our notation m) ... chosen coint rank
xy2.civar                               # details: xy.civar$beta 
summary(xy2.civar$rlm)                  # no t-values for beta, AIC
#pix <- xy2.cit@ZK %*% t(xy.cit@PI)     # Pi * x(t-1)
#varlevel <- vec2var(xy2.cit, r=1)      # coeff matrices of the lagged endogenous vars

## EXERCISE: 
# (1) Compare series and scatter plots for different periods.
# (2) Test for unit roots for different periods. 
# (3) Test for cointegration for different periods.
# (4) Estimate the CIVAR and comment the result.
# (5) Write down at least one estimated CIVAR model in conventional notation.
#
# Do this analysis for 2 indices, and 3 indices of your choice.
# For the latter you have to modify the R script.