Guest User

Untitled

a guest
Jul 18th, 2018
73
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 4.20 KB | None | 0 0
  1. ---
  2. title: "CCM without Tidyverse"
  3. output: html_document
  4. ---
  5.  
  6.  
  7. ```{r setup, include=FALSE}
  8. knitr::opts_chunk$set(echo = TRUE, cache = TRUE)
  9.  
  10. indent1 = ' '
  11. indent2 = paste(rep(indent1, 2), collapse='')
  12. indent3 = paste(rep(indent1, 3), collapse='')
  13.  
  14. doeval = TRUE
  15. doecho = FALSE
  16.  
  17. library(knitr) # for kable
  18. #library(tidyverse) # column as nested data frame
  19. library(ggplot2)
  20. library(rEDM)
  21. library(TransferEntropy)
  22. library(forecast)
  23. library(lmtest)
  24. library(vars)
  25. knitr::opts_chunk$set(echo = TRUE)
  26. ```
  27.  
  28. ```{r read data, message=FALSE}
  29. # read the data from the file
  30. TeslaSearch <- read.csv("Tesla Search.txt",header = FALSE,col.names=c("Tesla Search")) #%>% as.tbl
  31. TeslaSale <- read.csv("Tesla Sale.txt",header = FALSE, col.names=c("Tesla Sale")) #%>% as.tbl
  32.  
  33. AudiSale <- read.csv("AudiSale.txt",header = FALSE,col.names=c("Audi Sale")) #%>% as.tbl
  34. AudiSearch <- read.csv("AudiSearch.txt",header = FALSE,col.names = c("Audi Search")) #%>% as.tbl
  35. ```
  36.  
  37.  
  38.  
  39.  
  40. ```{r remove na, message=FALSE}
  41. TeslaSale <- na.omit(TeslaSale)
  42. TeslaSearch <- na.omit(TeslaSearch)
  43.  
  44. AudiSale <- na.omit(AudiSale)
  45. AudiSearch <- na.omit(AudiSearch)
  46.  
  47.  
  48.  
  49. # CONVERT THE DATA TO STATIONARY TIME SERIES
  50. N1 <- nrow(AudiSale)
  51. N2 <- nrow(AudiSearch)
  52.  
  53. Audi.ts1 <- cbind(log(AudiSearch[2:N1,] /AudiSearch[1:N1-1,]),log(AudiSale[2:N1,]/ AudiSale[1:N1-1,])) #%>% as.tbl
  54. Tesla.ts1 <- cbind(log(TeslaSearch[2:N2,] /TeslaSearch[1:N2-1,]),log(TeslaSale[2:N2,]/ TeslaSale[1:N2-1,])) #%>% as.tbl
  55. ```
  56.  
  57.  
  58. ```{r added if not writing in `tidyverse`, message=FALSE echo=FALSE}
  59. # **must** be added if not writing in `tidyverse`
  60. colnames(Audi.ts1) <- c("Audi.Search","Audi.Sale")
  61. colnames(Tesla.ts1) <- c("Tesla.Search","Tesla.Sale")
  62. ```
  63.  
  64.  
  65. ```{r granger test, message=FALSE}
  66.  
  67. audi.var1 <- VAR(Audi.ts1, p =1 , type = "const")
  68. causality(audi.var1, cause = "Audi.Search")
  69.  
  70.  
  71. tesla.var1 <- VAR(Tesla.ts1, p =1 , type = "const")
  72. causality(tesla.var1, cause = "Tesla.Search")
  73.  
  74. ```
  75.  
  76.  
  77. 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.
  78.  
  79. ```
  80. ## Compute the TE from AudiSearch to AudiSale
  81. computeTE(AudiSearch,AudiSale,3,1,"MI_diff") ## should be circa 0.16
  82. ```
  83.  
  84.  
  85. ```{r}
  86. block_Tesla <- cbind(TeslaSearch,TeslaSale)
  87.  
  88. block_Audi <- cbind(AudiSearch,AudiSale)
  89.  
  90. vars <- colnames(block_Audi)
  91.  
  92.  
  93. # generate all combinations of lib_column, target_column, tp
  94. params <- expand.grid(lib_column = vars,
  95. target_column = vars,
  96. tp = 1:12)
  97.  
  98. # throw out cases where lib == target
  99. params <- params[params$lib_column != params$target_column, ]
  100.  
  101. # E = 3 is optimal or very close to optimal for both vars
  102. # In other circumstances, we should use the best univariate E for each lib_column
  103. E <- 3
  104.  
  105.  
  106. audi.output <- do.call(rbind, lapply(seq_len(NROW(params)), function(i) {
  107. ccm(block_Audi, E = 3,
  108. lib_sizes = NROW(block_Audi), random_libs = FALSE,
  109. lib_column = params$lib_column[i],
  110. target_column = params$target_column[i],
  111. tp = params$tp[i], silent = TRUE)
  112. }
  113. )
  114. )
  115. kable(audi.output)
  116. ```
  117.  
  118.  
  119. ```
  120. audi.output$direction <- paste(output$lib_column, "xmap to", audi.output$target_column)
  121.  
  122. # visualize
  123. ggplot(output, aes(x = tp, y = rho, color = direction)) +
  124. geom_line() + theme_bw()
  125. ```
  126.  
  127. # Exaple for CCM test
  128.  
  129. Following is the standard example from the r-blogger.
  130.  
  131. ```{r example, message=FALSE}
  132. data(paramecium_didinium)
  133.  
  134.  
  135. vars <- names(paramecium_didinium)[2:3] # c("paramecium", "didinium")
  136.  
  137. # generate all combinations of lib_column, target_column, tp
  138. params <- expand.grid(lib_column = vars,
  139. target_column = vars,
  140. tp = -10:10)
  141.  
  142. # throw out cases where lib == target
  143. params <- params[params$lib_column != params$target_column, ]
  144.  
  145. # E = 3 is optimal or very close to optimal for both vars
  146. # In other circumstances, we should use the best univariate E for each lib_column
  147. E <- 3
  148.  
  149. # Perform cross mapping runs
  150.  
  151. output <- do.call(rbind, lapply(seq_len(NROW(params)), function(i) {
  152. ccm(paramecium_didinium, E = 3,
  153. lib_sizes = NROW(paramecium_didinium), random_libs = FALSE,
  154. lib_column = params$lib_column[i],
  155. target_column = params$target_column[i],
  156. tp = params$tp[i], silent = TRUE)
  157. }))
  158.  
  159. kable(output)
  160.  
  161. ```
Add Comment
Please, Sign In to add comment