成对相关R代码 [英] Pairwise correlation-R code

问题描述

我在R中有一个时间序列矩阵（动物园对象），该矩阵在每列中显示了不同资产（在我的宇宙范围内）的历史价格；我有20列（20种不同的股票）。我将第一个矩阵转换为每周收益矩阵。
我的目标是找到快速的R代码，以帮助我计算<成对每周滚动52周回报相关性。这意味着该代码每周将以一个股票（1）的每周收益与其他19个股票的1y关联性逐一计算。 19x18x17x ....相关时间序列。想法是每周提取所有这些相关性的中值。
感谢您的帮助

i've got in R a matrix of time series ( zoo object) which shows in each column historical prices for different equity ( within my universe); I have 20 columns ( 20 different stocks ). I transform this first matrix into a matrix of weekly return. My goal is to find a quick R code that help me to calcultate the mean of < the pairwise "rolling 52 weeks weekly return correlation">. It imply that every week, the code will compute one by one the 1y correlation of the stocks(1)'weekly return with the other 19 stocks . 19x18x17x.... Correlation time series. The idea is to extract every week the median of these all correlation. Thanks for your help

2016-20-10-跟进
这是下面的代码，用于计算平均成对相关性，但是我

2016-20-10 - follow up here is the code below that calculate the average pairwise correlation, but i am looking for to calculate the median.

#EMU_NAV_D : Matrix of daily prices
#EMU_ret_d : : Matrix of daily prices
#index : to count iteration 
index = which((count(t(EMU_NAV_D)) > 49 ))
index = index[-c(1:252)]   
# average correlation among SX5E components ( Eurostoxx 50 => 50  securities)
#avg.cor => initialisation matrice finalle des average correl
avg.cor = NA * EMU_NAV[,1]
for(i in index) {
    hist = EMU_ret_d[ (i- 252 +1):i, ]
    hist = hist[ , count(hist)==252, drop=F]
    nleft = ncol(hist)
    correlation = cor(hist, use='complete.obs',method='pearson')
    avg.cor[i,] = (sum(correlation) - nleft) / (nleft*(nleft-1))         
}

推荐答案

以下是对10只股票的滚动中值相关性的计算，可以很容易地将其缩放为N只股票

The following is calculation of rolling median correlation for 10 stocks, which could easily be scaled to N stocks.

重要函数是 getSymbols ， Return.calculate 和 rollapply

加载数据：

library(xts)       # for time series manipulation
library(zoo)       # for na.locf,na.approx
library(quantmod)  #for downloading closing prices from yahoo
library(PerformanceAnalytics) # for performance calculations and plotting
library(reshape2)  #for melt function



symbolList = c("AAPL", "GOOG", "MSFT", "CSCO", "AMZN", "YHOO", "EBAY","NVDA", "BIDU", "FB")


#Create new environment to save price data from yahoo

dataStocks <- new.env()
getSymbols(symbolList, src = 'yahoo', from = '2005-01-01', env = dataStocks, auto.assign = TRUE)

价格汇总：

#function to extract closing prices from OHLC object ( Open,High,Low,Close)

getClosePx = function(envObj=env,symbolName="AAPL") {

closePx = Cl(envObj[[symbolName]])
colnames(closePx) = symbolName

cat("Symbol is",symbolName,"\n")

return(closePx)
}

closePxList = lapply(symbolList,function(x) getClosePx(envObj=dataStocks,symbolName=x) )


#Create merged time series of all prices

dailyPx = Reduce(function(x,y) merge.xts(x,y),closePxList)



head(dailyPx)
#            AAPL     GOOG  MSFT  CSCO  AMZN  YHOO   EBAY     NVDA BIDU FB
#2005-01-03 63.29 202.7104 26.74 19.32 44.52 38.18 114.11 23.58000   NA NA
#2005-01-04 63.94 194.5003 26.84 18.56 42.14 36.58 111.31 22.47000   NA NA
#2005-01-05 64.50 193.5103 26.78 18.57 41.77 36.13 110.90 22.68000   NA NA
#2005-01-06 64.55 188.5503 26.75 18.85 41.05 35.43 106.18 22.46001   NA NA
#2005-01-07 69.25 193.8503 26.67 18.72 42.32 35.96 106.58 22.02999   NA NA
#2005-01-10 68.96 195.0603 26.80 18.72 41.84 36.32 107.31 22.08000   NA NA
tail(dailyPx)
#             AAPL   GOOG  MSFT  CSCO   AMZN  YHOO  EBAY  NVDA   BIDU     FB
#2016-10-12 117.34 786.14 57.11 30.34 834.09 42.36 31.50 66.43 175.41 129.05
#2016-10-13 116.98 778.19 56.92 30.17 829.28 41.62 31.51 65.35 174.61 127.82
#2016-10-14 117.63 778.53 57.42 30.18 822.96 41.44 31.89 65.99 175.51 127.88
#2016-10-17 117.55 779.96 57.22 30.22 812.95 41.79 31.81 65.61 175.15 127.54
#2016-10-18 117.47 795.26 57.66 30.44 817.65 41.68 31.64 66.61 175.65 128.57
#2016-10-19 117.12 801.50 57.53 30.35 817.69 42.73 32.52 66.47 176.20 130.11


#Omit Rows with missing values
dailyPx = dailyPx[complete.cases(dailyPx),]

#If missing values are to retained, you can use na.locf to carry forward previous non-missing value
#dailyPx = na.locf(dailyPx)

回报计算：

#Use na.approx function from zoo package to convert daily time series to weekly frequency
weeklyPx = na.approx(dailyPx, xout = seq(start(dailyPx), end(dailyPx), by = "week"))

#Calculate weekly returns
weeklyRet = na.omit(Return.calculate(weeklyPx,method="discrete"))

滚动相关：

#lookbackPeriod = 52 weeks
lookbackPeriod = 52


#Function to calculation median correlation from correlation matrix
#we take lower triangle of correlation matrix which contains all required
#pair-wise correlations and calculate median 

medianCor = function(returnMatrix=NULL) {

corMatrix = cor(returnMatrix)

medianCorValue = median(corMatrix[lower.tri(corMatrix)],na.rm=TRUE)

return(medianCorValue)

}


#Use rollapply to calculate median correlation on moving window of 52 weeks

rollMedianCorr = rollapply(weeklyRet,width= lookbackPeriod ,FUN=function(x) medianCor(x) ,by=1,fill=NA,align='right',by.column=FALSE)
colnames(rollMedianCorr) = "1Year_Rolling_Median_Correlation"

图：

#Plot series: base graphics
chart.TimeSeries(rollMedianCorr)

#Plot series: ggplot

DF = data.frame(date=index(rollMedianCorr),coredata(rollMedianCorr))
DF_melt = melt(DF,id.vars= "date")

customTitle = "1Year_Rolling_Median_Correlation"

gg<-ggplot(DF_melt,aes(x=date,y=value,color=variable))+geom_line() + 
   ggtitle(customTitle) + theme(plot.title = element_text(lineheight=.8, face="bold"),legend.position = "none")

print(gg)

这篇关于成对相关R代码的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！