Advertisement
mjaniec

Testing the significancy of correlations

Aug 31st, 2012
278
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
R 2.92 KB | None | 0 0
  1. # http://reakkt.com
  2.  
  3. library(fBasics) # stableFit
  4.  
  5. k <- 1000
  6.  
  7. n <- 250
  8.  
  9. ###
  10.  
  11. assets <- c("pko","peo")
  12.  
  13. market.data <- list()
  14.  
  15. for (asset in assets)
  16.  
  17.   market.data[[asset]] <- read.csv(paste(asset,"_d.csv",sep=""),header=TRUE)
  18.  
  19. available.data <- min(unlist(lapply(market.data, nrow)))
  20.  
  21. # stableFit:
  22.  
  23. # ceny aktywow
  24. x <- do.call(cbind, lapply(market.data, function(i) tail(i[,"Close"],available.data)))
  25.  
  26. d <- cbind(log(x[-1,1]/x[-available.data,1]),log(x[-1,2]/x[-available.data,2]))
  27.  
  28. sf1 <- stableFit(d[,1])
  29. sf2 <- stableFit(d[,2])
  30.  
  31. ###
  32.  
  33. nv <- c(5,10,25,50,100,250,500,1000)
  34.  
  35. nl <- length(nv)
  36.  
  37. cor.mc.norm             <- matrix(NA,k,nl)
  38. cor.mc.stable           <- matrix(NA,k,nl)
  39. colnames(cor.mc.norm)   <- nv
  40. colnames(cor.mc.stable) <- nv
  41.  
  42. distFun <- function(dist.type,ni) {
  43.  
  44.   if (dist.type=="norm")  x <- rnorm(ni,0,6)
  45.   if (dist.type=="sf1")   x <- rstable(ni,sf1@fit$estimate["alpha"],sf1@fit$estimate["beta"],sf1@fit$estimate["gamma"],sf1@fit$estimate["delta"])
  46.   if (dist.type=="sf2")   x <- rstable(ni,sf2@fit$estimate["alpha"],sf2@fit$estimate["beta"],sf2@fit$estimate["gamma"],sf2@fit$estimate["delta"])
  47.    
  48.   return( x )  
  49.    
  50. }
  51.  
  52. for (ni in nv) {
  53.  
  54.   for (i in 1:k) {
  55.  
  56.     x <- distFun("sf1",ni)
  57.     y <- distFun("sf2",ni)
  58.    
  59.     xn <- distFun("norm",ni)
  60.     yn <- distFun("norm",ni)
  61.    
  62.     cor.mc.stable[i,which(nv==ni)] <- cor(x,y)
  63.     cor.mc.norm[i,which(nv==ni)]   <- cor(xn,yn)
  64.  
  65.   }
  66.  
  67. }
  68.  
  69. boxplot(cor.mc.stable,main="Random correlations depending on series length",ylab="correlation",xlab="period")
  70. lines(apply(cor.mc.norm, 2, max),col="Blue",lty="dotted")
  71. lines(apply(cor.mc.norm, 2, min),col="Blue",lty="dotted")
  72.  
  73. #####
  74.  
  75. ### biezaca korelacja dla okresu n:
  76.  
  77. # ceny aktywow
  78. x <- do.call(cbind, lapply(market.data, function(i) tail(i[,"Close"],available.data)))
  79.  
  80. x <- tail(x,(n+1))
  81.  
  82. # zmiany procentowe
  83. d <- cbind(log(x[-1,1]/x[-n,1]),log(x[-1,2]/x[-n,2]))
  84.  
  85. (cor.current <- cor(d[,1],d[,2]))
  86.  
  87. ###
  88.  
  89. cor.sample <- numeric(k)
  90.  
  91. s1 <- round(runif(k,(n+1),available.data))
  92.  
  93. x <- do.call(cbind, lapply(market.data, function(i) tail(i[,"Close"],available.data)))
  94.  
  95. for (i in 1:k) {
  96.  
  97.   xs <- x[(s1[i]-n):s1[i],]
  98.  
  99.   d <- cbind(log(xs[-1,1]/xs[-n,1]),log(xs[-1,2]/xs[-n,2]))
  100.  
  101.   cor.sample[i] <- cor(d[,1],d[,2])
  102.  
  103. }
  104.  
  105. density.sample <- density(cor.sample)
  106. density.random <- density(cor.mc.stable[,paste(n)])
  107.  
  108. ylim=range(density.sample$y,density.random$y)
  109. xlim=range(density.sample$x,density.random$x)
  110.  
  111. plot(density.sample,"Historical Correlation Density",xlim=xlim,ylim=ylim)
  112. lines(density.random,lty="dotted",col="Blue")
  113. abline(v=cor.current,lty="dashed",col="Red")
  114. abline(v=mean(cor.sample),lty="dotted",col="Orange")
  115. mtext(paste("period=",n,
  116.             "mean=",format(mean(cor.sample),digits=4),
  117.             "current=",format(cor.current,digits=4)),
  118.       line=0.5,cex=0.8)
  119. legend("topleft",legend=c("mean","current"),text.col=c("Orange","Red"))
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement