0
$\begingroup$

I am trying to implement the (average) total correlation of assets, which is discussed here and here in R. Specifically, I am looking at $\rho_{av(1)}$ and $\rho_{av(2)}$:

$$ \rho_{av(1)} = \frac{2 \sum_{i=1}^N \sum_{j>i}^N w_i w_j \rho_{i, j}}{1 - \sum_{i=1}^N w_i^2} $$

$$ \rho_{av(2)} = \frac{\sigma^2 - \sum_{i=1}^N w_i^2\sigma_i^2}{2 \sum_{i=1}^N \sum_{j>i}^N w_i w_j \sigma_i \sigma_j} $$

So far, I have the following code:

set.seed(123)
a = rnorm(100, 0, 0.02)
b = rnorm(100, 0, 0.03)
c = rnorm(100, 0, 0.04)
dt = cbind(a, b, c)
wg = runif(ncol(dt))
wg = wg/sum(wg)
wg = round(wg, digits = 3)

options(scipen = 999)
cv = cov(dt)
VARs = diag(cv)
cr = cor(dt)
length(wg)

sum_1 = 0
for (i in 1 : length(wg) ) {
  for (j in 2 : length(wg) ) {
    sum_1 = sum_1 + 2 * wg[i] * wg[j] * cr[i, j]
  }}
rho_1 = sum_1/(1 - sum(wg^2))

sum_2 = 0
for (i in 1 : length(wg) ) {
  for (j in 2 : length(wg) ) {
    sum_2 = sum_2 + 2 * wg[i] * wg[j] * sqrt(VARs[i]) * sqrt(VARs[j])
  }}

PV = t(wg) %*% cv %*% wg
ss = sum(wg^2 * VARs)
rho_2 = (PV - ss)/sum_2

rhos = c(rho_1, rho_2)
rhos

And the results are:

[1]  1.388473356 -0.004104948

I expected these two to be close to each other. I think there might be an error in my code. I would appreciate if somebody could verify the code.

Thanks.

$\endgroup$
  • $\begingroup$ the inner of each of your for loops should be for(j in i+1 : length(wg)), you are double counting the total variance, and since one is a numerator, an the other a denominator, it's going to push the two father away from each other. I would expect that to only account for a factor of 4 difference though, which is much less than what you have. Additionally, the first number is clearly wrong -> average correlation > 1? i don't think so. $\endgroup$ – will Jun 19 '17 at 10:06
  • $\begingroup$ in that case my correlation matrix ref gets out of bounds ... Exactly, it should be between [-1; 1]. $\endgroup$ – AK88 Jun 19 '17 at 10:08
  • $\begingroup$ try for(j in min(i+1, length(wg)) : length(wg))... $\endgroup$ – will Jun 19 '17 at 10:09
  • $\begingroup$ got closer: 0.10906278; -0.01044804. I think there is still something wrong ... $\endgroup$ – AK88 Jun 19 '17 at 10:12
0
$\begingroup$

Here is the corrected original code (thanks @will):

 sum_1 = 0
  for (i in 1 : (length(wg) - 1) ) {
    for (j in  min(i+1, length(wg)) : length(wg) ) {
      sum_1 = sum_1 + 2 * wg[i] * wg[j] * cr[i, j]
    }}
  rho_11 = sum_1/(1 - sum(wg^2))

  sum_2 = 0
  for (i in 1 : (length(wg) - 1) ) {
    for (j in  min(i+1, length(wg)) : length(wg) ) {
      sum_2 = sum_2 + 2 * wg[i] * wg[j] * sqrt(VARs[i]) * sqrt(VARs[j])
    }}

  PV = t(wg) %*% cv %*% wg
  ss = sum(wg^2 * VARs)
  rho_22 = (PV - ss)/sum_2

A vectorized version:

  # double summation in numerator (including multiplier 2)
  diag(cr) <- 0
  double_sum.1 <- c(crossprod(wg, cr %*% wg))

  # single summation in denominator
  single_sum.1 <- c(crossprod(wg))
  rho_1 = double_sum.1/(1 - single_sum.1)


  # single summation in numerator
  single_sum.2 <- c(crossprod(wg * sqrt(diag(cv))))
  # double summation in denominator
  double_sum.2 <- sum(tcrossprod(wg * sqrt(diag(cv)))[upper.tri(diag(length(wg)))])
  PV = t(wg) %*% cv %*% wg
  rho_2 = (PV - single_sum.2)/ (2 * double_sum.2)

And the comparison of the two methods (initially they were way off): enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.