Попарная корреляция-R-код

#r #correlation

#r #корреляция

Вопрос:

у меня есть в R матрица временных рядов (объект zoo), которая показывает в каждом столбце исторические цены для разных акций (в моей вселенной); У меня есть 20 столбцов (20 разных акций). Я преобразую эту первую матрицу в матрицу еженедельного дохода. Моя цель — найти быстрый R-код, который поможет мне вычислить среднее значение < попарной «корреляции с еженедельным возвратом за 52 недели»>. Это означает, что каждую неделю код будет вычислять одну за другой 1y корреляцию недельной доходности акций (1) с другими 19 акциями. 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))         
}
  

Комментарии:

1. Добро пожаловать в SO!, chart.RollingCorrelation из PerformanceAnalytics пакета было бы полезно, что является функцией-оболочкой для ?rollapply

2. Привет, Оссан, большое спасибо за ваше сообщение, я нашел код в R bloggers, который почти соответствует моим потребностям, но я немного пытаюсь внести небольшие изменения,

Ответ №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)
  

введите описание изображения здесь

Комментарии:

1. Дорогой «Osssan», я только что понял, что вы отвечаете на мой «вопрос». Большое спасибо, это идеально! Лучшие

2. Если решение решает вашу проблему, вы можете принять ответ и отметить вопрос как решенный, нажав на стрелку слева