---
title: "CCM test"
date: "June 27, 2018"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, cache = TRUE)

indent1 = '    '
indent2 = paste(rep(indent1, 2), collapse='')
indent3 = paste(rep(indent1, 3), collapse='')

doeval = TRUE
doecho = FALSE

library(knitr)       # for kable
library(tidyverse)   # column as nested data  frame
library(ggplot2)
library(rEDM)
library(TransferEntropy)
library(forecast)
library(lmtest)
library(vars)
knitr::opts_chunk$set(echo = TRUE)
```

```{r read data, message=FALSE}
# read the data from the file
TeslaSearch <- read.csv("Tesla Search.txt",header = FALSE,col.names=c("Tesla Search")) %>% as.tbl
TeslaSale <- read.csv("Tesla Sale.txt",header = FALSE, col.names=c("Tesla Sale")) %>% as.tbl

AudiSale <- read.csv("AudiSale.txt",header = FALSE,col.names=c("Audi Sale")) %>% as.tbl
AudiSearch <- read.csv("AudiSearch.txt",header = FALSE,col.names = c("Audi Search")) %>% as.tbl
```




```{r remove na, message=FALSE}
TeslaSale <- na.omit(TeslaSale)
TeslaSearch <- na.omit(TeslaSearch)

AudiSale <- na.omit(AudiSale)
AudiSearch <- na.omit(AudiSearch)



# CONVERT THE DATA TO STATIONARY TIME SERIES
N1 <- nrow(AudiSale)
N2 <- nrow(AudiSearch)

Audi.ts1 <- cbind(log(AudiSearch[2:N1,] /AudiSearch[1:N1-1,]),log(AudiSale[2:N1,]/ AudiSale[1:N1-1,])) %>% as.tbl
Tesla.ts1 <- cbind(log(TeslaSearch[2:N2,] /TeslaSearch[1:N2-1,]),log(TeslaSale[2:N2,]/ TeslaSale[1:N2-1,])) %>% as.tbl
```


```{r granger test, message=FALSE}

audi.var1 <- VAR(Audi.ts1, p =1 , type = "const")
causality(audi.var1, cause = "Audi.Search")


tesla.var1 <- VAR(Tesla.ts1, p =1 , type = "const")
causality(tesla.var1, cause = "Tesla.Search")

```


Under the confidance level $\alpha = 0.05$, Audi.Search do not Granger-cause Audi.Sale. But Tesla.Search do not Granger-cause Tesla.Sale with lag = 1.

```
## Compute the TE from AudiSearch to AudiSale
computeTE(AudiSearch,AudiSale,3,1,"MI_diff") ## should be circa 0.16
```


```{r}
block_Tesla <- cbind(TeslaSearch,TeslaSale)

block_Audi <- cbind(AudiSearch,AudiSale)

vars <- colnames(block_Audi)


# generate all combinations of lib_column, target_column, tp
params <- expand.grid(lib_column = vars,
                      target_column = vars,
                      tp = 1:12)

# throw out cases where lib == target
params <- params[params$lib_column != params$target_column, ]

# E = 3 is optimal or very close to optimal for both vars
# In other circumstances, we should use the best univariate E for each lib_column
E <- 3


audi.output <- do.call(rbind, lapply(seq_len(NROW(params)), function(i) {
    ccm(block_Audi, E = 3,
        lib_sizes = NROW(block_Audi), random_libs = FALSE,
        lib_column = params$lib_column[i],
        target_column = params$target_column[i],
        tp = params$tp[i], silent = TRUE)
  }
  )
  )
kable(audi.output)
```


```
audi.output$direction <- paste(output$lib_column, "xmap to", audi.output$target_column)

# visualize
ggplot(output, aes(x = tp, y = rho, color = direction)) +
    geom_line() + theme_bw()
```

# Exaple for CCM test

Following is the standard example from the r-blogger.

```{r example, message=FALSE}
data(paramecium_didinium)


vars <- names(paramecium_didinium)[2:3] # c("paramecium", "didinium")

# generate all combinations of lib_column, target_column, tp
params <- expand.grid(lib_column = vars,
                      target_column = vars,
                      tp = -10:10)

# throw out cases where lib == target
params <- params[params$lib_column != params$target_column, ]

# E = 3 is optimal or very close to optimal for both vars
# In other circumstances, we should use the best univariate E for each lib_column
E <- 3

# Perform cross mapping runs

output <- do.call(rbind, lapply(seq_len(NROW(params)), function(i) {
    ccm(paramecium_didinium, E = 3,
        lib_sizes = NROW(paramecium_didinium), random_libs = FALSE,
        lib_column = params$lib_column[i],
        target_column = params$target_column[i],
        tp = params$tp[i], silent = TRUE)
}))

kable(output)

```