Ledger Nano X - The secure hardware wallet
SHARE
TWEET

Untitled

a guest Apr 7th, 2020 265 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  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)
RAW Paste Data
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
Top