# R: 3.4.1
# Ex1_6_Granger_R.txt    Granger causality: US stock to other stock indices wrt returns.
# For volatility use e.g. rx <- volx <- log((x.h - x.l)/x.c) instead of return 

#setwd("C:/LVs/QF_FinancialEconometrics/R_Test") 

library("tseries")
library("zoo")
library("dynlm")

require(quantmod) # Finanzdaten download
#####################################
LjungBox <- function(y,maxlag) { lb <- c(1:(3*maxlag)); dim(lb) <- c(maxlag,3); 
    for (ii in 1:maxlag) { 
       bb <- Box.test(y, lag =  ii, type = "Ljung-Box", fitdf = 0)
       for (jj in 1:3) { lb[ii,jj] <- as.numeric(bb[jj]) }
                          }
    print(c("Ljung-Box statistics:"))
    print(c("Chi^2", "lag", "p-value"))
    lb 
                                }
#####################################
spc <- new.env()
dax <- new.env()
nasdaq <- new.env()

setDefaults(getSymbols,src="yahoo")
getSymbols("^GSPC", env = spc, from = "1995-01-01", to = "2013-07-18", return.class = "xts")
# to = "2006-12-31")
getSymbols("^GDAXI", env = dax, from = "1995-01-01", to = "2013-07-18", return.class = "xts")
getSymbols("^NDX", env = nasdaq, from = "1995-01-01", to = "2013-07-18", return.class = "xts")
# read in series are of xts type,  plot is different for xts and ts

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

### Plotting univariate data
plot( merge(as.zoo(x.c),as.zoo(y.c)) , main = "^GSPC , ^GDAXI")  # plot price series, one scale, indiv NAs     

rx  <- diff(log(x.c))[2:nx];         
ry  <- diff(log(y.c))[2:ny];         
plot(merge(as.zoo(rx),as.zoo(ry)), main = "^GSPC , ^GDAXI")  # plot returns  

# simultaneous analysis
xy  <- merge(x.c,y.c, all=FALSE)        # one scale no indiv NAs
nxy <- length(xy[,1])
rxy <- diff(log(xy))[2:nxy] 

### Testing for white noise
rx_v <- as.vector(rxy[,1]); ry_v <- as.vector(rxy[,2]);    
acf(rx_v);                                               # acf of returns
acf(ry_v, na.action=na.pass); 

LjungBox(rx_v, 10)                                       # Ljung-Box tests
LjungBox(ry_v, 10)  

### Granger causality , choose type of causality 
### choose maximal lags (preset as 3, 3)
rx_z <- as.zoo(rx);    ry_z <- as.zoo(ry);  # dynlm needs zoo data

## model for Granger causality 
#form0 <- " ry_z ~    lag(ry_z,-1) + lag(ry_z,-2) + lag(ry_z,-3) "
#formb <- "         + lag(rx_z,-1) + lag(rx_z,-2) + lag(rx_z,-3) "
#form1 <- paste(form0,formb,sep="")
#df.2 <- 3

## model for instantaneous Granger causality 
form0 <- " ry_z ~       lag(ry_z,-1) + lag(ry_z,-2) + lag(ry_z,-3) "
formb <- "     + rx_z + lag(rx_z,-1) + lag(rx_z,-2) + lag(rx_z,-3)  "
form1 <- paste(form0,formb,sep="")
df.2 <- 3 + 1

## a model for feedback causality 
#form0 <- " rx_z ~     lag(rx_z,-1) + lag(rx_z,-2) + lag(rx_z,-3) "
#formb <- "          + lag(ry_z,-1) + lag(ry_z,-2) + lag(ry_z,-3)  "
#form1 <- paste(form0,formb,sep="")
#df.2 <- 3

## define the subperiods
## September 11, 2001
start0_1 = "2001-06-10"; end0_1 = "2001-09-10"; 
start0_2 = "2001-09-13"; end0_2 = "2001-12-31";
## Lehman Brothers, Chapter 11 bankruptcy protection Sep 15,2008
#start0_1 = "2006-01-01"; end0_1 = "2008-07-01"; 
#start0_2 = "2008-10-01"; end0_2 = "2009-12-31";

## models for the whole period: without(0)/with causality(1)
m0 <- dynlm(formula = as.formula(form0) , start= start0_1, end=end0_2)
m1 <- dynlm(formula = as.formula(form1) , start= start0_1, end=end0_2)
c_meas <- log( summary(m0)$sigma^2 / summary(m1)$sigma^2 )     # causality measure of strength         
p_lrtest <- 1 - pchisq( length( summary(m1)$residuals ) * c_meas, df.2 ) #LR-test
p_lrtest
summary(m0);
summary(m1); 

## models for the 1st subperiod: without(0)/with causality(1)
m0_1 <- dynlm(formula = as.formula(form0) , start=start0_1, end=end0_1)
m1_1 <- dynlm(formula = as.formula(form1) , start=start0_1, end=end0_1)
c_meas_1 <- log(summary(m0_1)$sigma^2 / summary(m1_1)$sigma^2 )  # causality measure of strength    
p_lrtest_1 <- 1 - pchisq( length( summary(m1_1)$residuals ) * c_meas_1, df.2) #LR-test
p_lrtest_1
summary(m0_1); 
summary(m1_1); 

## models for the 2nd subperiod: without(0)/with causality(1)
m0_2 <- dynlm(formula = as.formula(form0) , start=start0_2, end=end0_2)
m1_2 <- dynlm(formula = as.formula(form1) , start=start0_2, end=end0_2)
c_meas_2 <- log(summary(m0_2)$sigma^2 / summary(m1_2)$sigma^2 )  # causality measure of strength    
p_lrtest_2 <- 1 - pchisq( length( summary(m1_2)$residuals ) * c_meas_2, df.2) #LR-test
p_lrtest_2
summary(m0_2); 
summary(m1_2); 

print( ' Causality measures: M, M_1, M_2')
print( c(c_meas, c_meas_1, c_meas_2) )


## Comment:
# (A) dynlm does not use the data for the subset(start, end) only, if lags are 
# introduced, but also uses the observations before start for the lagged variables. 
# (B) if class(x) of a data-vector x is 
# "zoo"          : lag(x,-1) means lag 1 of x
# "xts" and "zoo": lag(x, 1) means lag 1 of x !


### EXERCISE:
#   Comment on every step. Interprete the results.
#   (1) Select the data and period for your investigation.
#   (2) Plot the data (prices, returns) and test for white noise.
#   (3) Check for Granger causality 
#   (3.a) Repeat it for instantanenous Granger causality 
#   (3.b) Consider the possibility of feedback causality 
#   (4) Test for Granger causality. Interprete the causality measures.
#   (5!) Write down the models and interprete the parameters.
#   (6) Select an appropriate maximal lag for the AR model for the y variable,
#       and the maximal lag for the augmented AR term for the x variable.