API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Apr 7th, 2020
340
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
R 13.15 KB | None | 0 0
raw download clone embed print report
  1. # get_os <- function(){
  2. #   sysinf <- Sys.info()
  3. #   if (!is.null(sysinf)){
  4. #     os <- sysinf['sysname']
  5. #     if (os == 'Darwin') os <- "osx"
  6. #   } else { ## mystery machine
  7. #     os <- .Platform$OS.type
  8. #     if (grepl("^darwin", R.version$os)) os <- "osx"
  9. #     if (grepl("linux-gnu", R.version$os)) os <- "linux"
  10. #   }
  11. #   as.character(tolower(os))
  12. # }
  13. # sourceName <- paste0("pacman_0.5.0", ifelse(get_os()=="windows",".zip",".tgz"))
  14. # if("pacman" %in% rownames(installed.packages()) == F) install.packages(sourceName, repos = NULL, type="win.binary")
  15. if("pacman" %in% rownames(installed.packages()) == F) install.packages("pacman")
  16. pacman::p_load(fUnitRoots, urca, vars, aod, zoo, tseries, igraph, visNetwork, abind, Hmisc)
  17. TIME <- Sys.time()
  18. Sig_95 <- 0.05 # p-value for stationarity testing with ADF and KPSS
  19. Sig_90 <- 0.1 # p-value for Wald test hypotheses of G-causality
  20. maxLag <- 12 # max lag number if VAR lag selection
  21. nAhead <- 24 # horizon for IRF and FEVD plots
  22. totOrd <- 3 # maximum order of integration to try for stationarity
  23. maxChr <- 25
  24. G_test <- c('chi2','Ftest') # Wald statistic for causality: 'chi2' or 'Ftest'
  25. names(G_test) <- G_test
  26. useLinReg <- FALSE # Should Wald statistic be estimated from linear regression?
  27. savePlots <- FALSE
  28.  
  29. #Load data
  30. setwd(dirname(sys.frame(1)$ofile)) # this line doesn't work on MAC OS
  31. csvFileName <- "LT_2003M07_2013M02"
  32. #csvFileName <- "coffee_data_from_1976"
  33. myData <- read.csv('sutvarkytas_ir_isglaudytas.csv', header=T)
  34. fixData <- function(data_arr, coef){
  35.   index = 1
  36.   res_arr = array()
  37.  
  38.   for (j in 1:6) {
  39.     for (i in 1:12) {
  40.       if(i == 12 && index != 6 * 12) {
  41.         gr = data_arr[index]
  42.         saus = data_arr[index + 1]
  43.         res = (gr + saus) / 2
  44.       } else if(i == 1 && index != 1) {
  45.         res = res_arr[index - 1]
  46.       } else {
  47.         res = data_arr[index]
  48.       }
  49.       res_arr[index] = res
  50.       index = index + 1
  51.     }
  52.   }
  53.   index = index - 1
  54.   res_arr[index] = res_arr[index] * coef
  55.   return(res_arr)
  56. }
  57.  
  58. myData$Y1_GPM = fixData(myData$Y1_GPM, 0.83)
  59. myData$Y2_GPM = fixData(myData$Y2_GPM, 0.83)
  60. cat(myDataN <- names(myData))
  61. myDataX <- as.Date(paste(sub("M","-",myData$Date),"-01",sep=""))
  62. #myDataX <- as.Date(as.yearqtr(sub("Q","-",myData$Date)))
  63. #myDataX <- as.Date(sub("\\.","-",sub("\\.","-",myData$Date)))
  64. #myDataX <- as.Date(paste(sub(".","-",myData$Date),sep=""))
  65. if (savePlots) dir.create(file.path(csvFileName),showWarnings=F)
  66.  
  67. testingPairs <- t(combn(myDataN[-1],2))
  68. #mainCols <- c('ESI','NetMig') # leave only variables of interest
  69. #testingPairs <- testingPairs[testingPairs[,1] %in% mainCols | testingPairs[,2] %in% mainCols,]
  70. v <- myDataN[-1]
  71. vMat <- matrix(0, nrow = length(v), ncol = length(v), dimnames = list(FROM=v, TO=v))
  72. crossImpact <- lapply(G_test, function(x) vMat)
  73. sumOfCoefs <- crossImpact
  74.  
  75. # create R variables from 'csv' columns and make pairs to test
  76. for (i in 1:length(myDataN)) {
  77.   eval(parse(text=paste(myDataN[i],' <- myData$',myDataN[i],sep='')))
  78. }
  79.  
  80. for (nr in 1:nrow(testingPairs)) {
  81.   myPair <- testingPairs[nr,]
  82.   eval(parse(text=paste('data1 <- ',myPair[1],sep='')))
  83.   eval(parse(text=paste('data2 <- ',myPair[2],sep='')))
  84.   varPair <- cbind(data1,data2)
  85.   colnames(varPair) <- myPair
  86.  
  87.   # visualize time series
  88.   ylim1 <- c(min(data1)*(1-sign(min(data1))*0.9),max(data1)*(1+sign(max(data1))/10))
  89.   plot(x=myDataX,y=data1,ylim=ylim1,col='red',type='l',lwd=1,xaxt='n',axes=F,xlab='Time (quarterly periodicity)',ylab='',main='')
  90.   axis(1,pretty(myDataX),labels=format(pretty(myDataX),'%Y'))
  91.   axis(2, pretty(ylim1), col='red')
  92.   par(new=T)
  93.   ylim2 <- c(min(data2)*(1-sign(min(data2))*0.9),max(data2)*(1+sign(max(data2))/10))
  94.   plot(x=myDataX,y=data2,ylim=ylim2,col='blue',type='l',lwd=1,xaxt='n',axes=F,ylab='',xlab='')
  95.   axis(4, pretty(ylim2), col='blue')
  96.   abline(h='', v=myDataX[which(as.integer(format(myDataX,'%m'))==1)], col="gray", lty=3);
  97.   box()
  98.   legend('bottom', myPair,col=c('red','blue'),lwd=c(1,1),horiz=(nchar(myPair[1])<maxChr && nchar(myPair[2])<maxChr))
  99.  
  100.   options(warn=-1) # suppress 'p-value greater than printed p-value' warning
  101.  
  102.   # 1. Test for integration. ADF / KPSS tests for unit-root / stationarity
  103.  
  104.   p_ADF <- matrix(rep(0,2*(totOrd+1)),2,(totOrd+1))
  105.   p_KPSS <- matrix(rep(0,2*(totOrd+1)),2,(totOrd+1))
  106.   p_ADF[1,1] <- adf.test(data1)$p.value
  107.   p_ADF[2,1] <- adf.test(data2)$p.value
  108.   p_KPSS[1,1] <- kpss.test(data1)$p.value
  109.   p_KPSS[2,1] <- kpss.test(data2)$p.value
  110.  
  111.   for (i in 1:totOrd) { #Test for unit roots / stationarity
  112.     eval(parse(text=paste('d_data1 <- diff(data1,',i,')',sep='')))
  113.     eval(parse(text=paste('d_data2 <- diff(data2,',i,')',sep='')))
  114.     p_ADF[1,i+1] <- adf.test(d_data1)$p.value
  115.     p_ADF[2,i+1] <- adf.test(d_data2)$p.value
  116.     p_KPSS[1,i+1] <- kpss.test(d_data1)$p.value
  117.     p_KPSS[2,i+1] <- kpss.test(d_data2)$p.value
  118.   }
  119.  
  120.   ordADF <- NULL
  121.   ordADF <- c(which(p_ADF[1,]<=Sig_95)[1],which(p_ADF[2,]<=Sig_95)[1])
  122.   ordADF[is.na(ordADF)] <- 1
  123.   ordKPSS <- NULL
  124.   ordKPSS <- c(which(p_KPSS[1,]>=Sig_95)[1],which(p_KPSS[2,]>=Sig_95)[1])
  125.   ordKPSS[is.na(ordKPSS)] <- 1
  126.  
  127.   options(warn=1) # restore warning setting back to default  
  128.  
  129.   # 2. Determine max order of integration (m).If none of the series
  130.   # in integrated, the usual Granger-causality test can be done. (m=0)
  131.  
  132.   maxOrd <- max(c(ordADF,ordKPSS)) - 1
  133.  
  134.   # 3. Set up VAR models to select lag order (DO NOT difference the data for VAR).
  135.  
  136.   bestLag <- VARselect(varPair,lag=maxLag,type="both")
  137.  
  138.   varResid <- NULL
  139.  
  140.   for (i in 1:length(bestLag$selection)) { #VAR Model
  141.     varModel <- VAR(varPair,p=bestLag$selection[i],type="both")
  142.     # Portmanteau- and Breusch-Godfrey test for autocorrelation
  143.     result <- try(serial.test(varModel))
  144.     if (class(result) == "try-error") {
  145.       varResid[i] <- 0
  146.     } else {
  147.       varResid[i] <- serial.test(varModel)$serial$p.value
  148.     }
  149.   }
  150.  
  151.   # 4. From best models (found by VARselect) choose the one,
  152.   # which has least autocorrelated residuals (most correct).
  153.  
  154.   optLag <- bestLag$selection[which(varResid == max(varResid))[1]]
  155.  
  156.   # 5. Augmented VAR model with additional lag VAR(p+m) is set up.
  157.   varModel <- VAR(varPair,p=optLag+maxOrd,type="both")
  158.  
  159.   if (useLinReg) {
  160.    
  161.     # VAR model is seperately set up as a linear model
  162.     # which turns the Wald testing more comprihensive
  163.    
  164.     #lag variables
  165.     zoo1 <- NULL; zoo1 <- zoo(varPair[,1])
  166.     zoo2 <- NULL; zoo2 <- zoo(varPair[,2])
  167.     idx <- 2:(optLag+maxOrd+1)
  168.     zoo1.l<-lag(zoo1,-(0:(optLag+maxOrd)),na.pad=T)
  169.     zoo2.l<-lag(zoo2,-(0:(optLag+maxOrd)),na.pad=T)
  170.    
  171.     lm1<-lm(zoo1~zoo2.l[,idx]+zoo1.l[,idx]+index(zoo1))
  172.     lm1$fitted.values <- na.omit(lm1$fitted.values)
  173.     lm2<-lm(zoo2~zoo1.l[,idx]+zoo2.l[,idx]+index(zoo2))
  174.     lm2$fitted.values <- na.omit(lm2$fitted.values)
  175.    
  176.   }
  177.  
  178.   # 6. Wald-test for the first p variables only with p degrees of freedom.
  179.   myTitle <- ''; arTitle <- paste('\n\nAR(',optLag,'+',maxOrd,')',sep='')
  180.  
  181.   wt1 <- NULL; wt2 <- NULL
  182.  
  183.   # Wald.test.1 (H0: 2 does not Granger-cause 1)
  184.   # Wald.test.2 (H0: 1 does not Granger-cause 2)
  185.   # Sum of coefficients to determine response direction
  186.   if (useLinReg) {
  187.     wt1 <- wald.test(b=coef(lm1), Sigma=vcov(lm1), Terms=c(2:(optLag+1)), df=(length(myDataX)-length(coef(lm1))-optLag+maxOrd*2)[[1]] )
  188.     wt2 <- wald.test(b=coef(lm2), Sigma=vcov(lm2), Terms=c(2:(optLag+1)), df=(length(myDataX)-length(coef(lm2))-optLag+maxOrd*2)[[1]] )
  189.     sum1 <- sum(coef(lm1)[c(2:(optLag+1+maxOrd))])
  190.     sum2 <- sum(coef(lm2)[c(2:(optLag+1+maxOrd))])
  191.   } else {
  192.     wt1 <- wald.test(b=coef(varModel$varresult[[1]]), Sigma=vcov(varModel$varresult[[1]]), Terms=seq(2,2*optLag,by=2),df=(length(myDataX)-length(coef(varModel$varresult[[1]]))-optLag+maxOrd*2)[[1]])          
  193.     wt2 <- wald.test(b=coef(varModel$varresult[[2]]), Sigma=vcov(varModel$varresult[[2]]), Terms=seq(1,2*optLag,by=2),df=(length(myDataX)-length(coef(varModel$varresult[[2]]))-optLag+maxOrd*2)[[1]])  
  194.     sum1 <- sum(coef(varModel$varresult[[1]])[seq(2,2*optLag+2*maxOrd,by=2)])
  195.     sum2 <- sum(coef(varModel$varresult[[2]])[seq(1,2*optLag+2*maxOrd,by=2)])    
  196.   }
  197.  
  198.   myTitle <- paste0(myTitle,sprintf('\n(coefs): %s {%5.3f} & %s {%5.3f}',myPair[1],sum2,myPair[2],sum1))
  199.  
  200.   p_wt1 <- NULL; p_wt2 <- NULL
  201.   for (t in 1:length(G_test)) {
  202.     eval(parse(text=paste('p_wt1[t] <- wt1$result$',G_test[t],'[[\'P\']]',sep='')))
  203.     eval(parse(text=paste('p_wt2[t] <- wt2$result$',G_test[t],'[[\'P\']]',sep='')))
  204.     G <- ' \227 '
  205.     if ((p_wt1[t]<Sig_90 && p_wt2[t]<Sig_90)) {
  206.       crossImpact[[G_test[t]]][myPair[1],myPair[2]] <- 1 - p_wt2[t]
  207.       crossImpact[[G_test[t]]][myPair[2],myPair[1]] <- 1 - p_wt1[t]
  208.       sumOfCoefs[[G_test[t]]][myPair[1],myPair[2]] <- sum2
  209.       sumOfCoefs[[G_test[t]]][myPair[2],myPair[1]] <- sum1
  210.       G <- ' \U2194 '
  211.       #G <- ' <-> '
  212.     } else if (p_wt1[t]<Sig_90) {
  213.       crossImpact[[G_test[t]]][myPair[2],myPair[1]] <- 1 - p_wt1[t]
  214.       sumOfCoefs[[G_test[t]]][myPair[2],myPair[1]] <- sum1
  215.       G <- ' \U2190 '
  216.       #G <- ' <- '
  217.     } else if (p_wt2[t]<Sig_90) {
  218.       crossImpact[[G_test[t]]][myPair[1],myPair[2]] <- 1 - p_wt2[t]
  219.       sumOfCoefs[[G_test[t]]][myPair[1],myPair[2]] <- sum2
  220.       G <- ' \U2192 '
  221.       #G <- ' -> '
  222.     }
  223.     myTitle <- paste0(myTitle,sprintf('\n(%s): %s [%5.3f]%s%s [%5.3f]',G_test[t],myPair[1],p_wt2[t],G,myPair[2],p_wt1[t]))
  224.     #myTitle <- paste0(myTitle,sprintf('\n(%s): ',G_test[t]),myPair[1],' [',sprintf('%5.3f',p_wt2[t]),']',G,myPair[2],' [',sprintf('%5.3f',p_wt1[t]),']')
  225.   }
  226.   cat(paste(arTitle,myTitle,sep=''))
  227.   title(myTitle,cex.main=0.8)
  228.  
  229.   # 7. IRF and FEVD plots
  230.   if (!useLinReg) {
  231.     if (any(p_wt1<Sig_90) && any(p_wt2<Sig_90)) {
  232.       if (savePlots) win.metafile(paste0(csvFileName,'/',myPair[2],'_response_to_',myPair[1],'.wmf'),width=8,height=6)
  233.       plot(irf(varModel, impulse = myPair[1], response = myPair[2], ci=1-Sig_90, runs=1000, n.ahead=nAhead))
  234.       if (savePlots) dev.off()
  235.       if (savePlots) win.metafile(paste0(csvFileName,'/',myPair[1],'_response_to_',myPair[2],'.wmf'),width=8,height=6)
  236.       plot(irf(varModel, impulse = myPair[2], response = myPair[1], ci=1-Sig_90, runs=1000, n.ahead=nAhead))
  237.       if (savePlots) dev.off()
  238.       if (savePlots) win.metafile(paste0(csvFileName,'/',myPair[1],'_',myPair[2],'_combined_FEVD.wmf'),width=6,height=8)
  239.       plot(fevd(varModel,n.ahead=nAhead))
  240.       if (savePlots) dev.off()
  241.     } else if (any(p_wt1<Sig_90)) {
  242.       if (savePlots) win.metafile(paste0(csvFileName,'/',myPair[1],'_response_to_',myPair[2],'.wmf'),width=8,height=6)
  243.       plot(irf(varModel, impulse = myPair[2], response = myPair[1], ci=1-Sig_90, runs=1000, n.ahead=nAhead))
  244.       if (savePlots) dev.off()
  245.       if (savePlots) win.metafile(paste0(csvFileName,'/',myPair[1],'_',myPair[2],'_combined_FEVD.wmf'),width=6,height=8)
  246.       plot(fevd(varModel,n.ahead=nAhead))
  247.       if (savePlots) dev.off()
  248.     } else if (any(p_wt2<Sig_90)) {
  249.       if (savePlots) win.metafile(paste0(csvFileName,'/',myPair[2],'_response_to_',myPair[1],'.wmf'),width=8,height=6)
  250.       plot(irf(varModel, impulse = myPair[1], response = myPair[2], ci=1-Sig_90, runs=1000, n.ahead=nAhead))
  251.       if (savePlots) dev.off()
  252.       if (savePlots) win.metafile(paste0(csvFileName,'/',myPair[1],'_',myPair[2],'_combined_FEVD.wmf'),width=6,height=8)
  253.       plot(fevd(varModel,n.ahead=nAhead))
  254.       if (savePlots) dev.off()
  255.     }
  256.   }
  257. }
  258.  
  259. # 8. construct and display graph from the cross-impact matrix
  260. for (t in 1:length(G_test)) {
  261.   iGraph <- graph_from_adjacency_matrix(crossImpact[[G_test[t]]],weighted=TRUE)
  262.   visGraph <- toVisNetworkData(iGraph)
  263.   iiGraph <- graph_from_adjacency_matrix(sumOfCoefs[[G_test[t]]],weighted=TRUE)
  264.   vvisGraph <- toVisNetworkData(iiGraph)
  265.   visGraph$nodes$shadow <- TRUE
  266.   visGraph$edges$color <- "red"
  267.   visGraph$edges$color[vvisGraph$edges$weight > 0] <- "green"
  268.   visGraph$edges$width <- 2 * visGraph$edges$weight ^ 9
  269.   visGraph$edges$title <- paste0('p-value = ',format(1 - visGraph$edges$weight, digits = 1))
  270.   #mutate(visGraph$edges, width = weight/5 + 1)
  271.   #plot.igraph(iGraph,main=G_test[t])
  272.   visnet <- visNetwork(nodes = visGraph$nodes, edges = visGraph$edges, main=paste0("G-causality by ",G_test[t])) %>%
  273.     #visIgraphLayout(layout = "layout_nicely") %>%
  274.     visLayout(improvedLayout=T) %>% visNodes(size=20) %>% visEdges(arrows="middle", physics=T) %>%
  275.     visOptions(highlightNearest=list(enabled=T, labelOnly=F, algorithm="hierarchical"), nodesIdSelection=list(enabled=T, values=v))
  276.   print(visnet)
  277. }
  278.  
  279. # 9. construct a scatterplot with active and passive sums from cross-impact matrix
  280. crossImpactArray <- abind(crossImpact,along=3)
  281. crossImpactBest <- apply(crossImpactArray,c(1,2),max)
  282. correlationResult <- rcorr(as.matrix(myData[,-1]))
  283. pooledCrossImpact <- (correlationResult$P < Sig_95) + (crossImpactBest > (1-Sig_90)) + (crossImpactBest > (1-Sig_95))
  284. XY <- data.frame(Asum=rowSums(pooledCrossImpact,na.rm=T),Psum=colSums(pooledCrossImpact,na.rm=T))
  285. plot(Asum ~ Psum, XY, ylim=c(0,max(XY$Asum)+1), xlim=c(0,max(XY$Psum)+1), main=csvFileName)
  286. abline(v = mean(XY$Psum), h = mean(XY$Asum), col="red", lwd=1, lty=2)
  287. text(Asum ~ Psum, XY, labels = row.names(XY), pos = 4)
  288.  
  289. cat('\n\n')
  290. print(Sys.time() - TIME)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!