ペアワイズ相関-Rコード

Question

私は、各列に異なる株価（宇宙内）の過去の価格を示す時系列（動物園オブジェクト）の行列をRに入れました。私は20の列（20の異なる株式）を持っています。私はこの最初の行列を毎週のリターンの行列に変換します。 私の目標は、<の平均を計算する手助けとなる迅速なRコードを見つけることです。これは、毎週、株式は株価（1）の週別収益と他の19株の相関を1つずつ計算することを意味します。 19x18x17x ....相関時系列。アイデアは、これらのすべての相関の中央値を毎週抽出することです。 あなたの助けをありがとう

2016-20-10 - フォローアップ ここでは、平均ペアワイズ相関を計算するが、私は中央値を計算するために探しているコードです。

#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))   
} 

Answer 1

以下は、10株のローリングメジアン相関の計算であり、これはN株に容易に換算できます。

重要な機能は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)