Simulating highly correlated data for 25 variables via multivariate normal distribution

212 views Asked by At

I try to simulate highly correlated (average absolute correlation of ~0.7) variables for 25 variables using 4,000 observations.

However, I can't get Sigma positive definite, and thus I need to use a trick to get sigma PD. And hereby I lose most of the initial correlation, I go from ~0.7 to ~0.25.

The reason might be that I need positive ánd negative correlated variables. I am able to simulate highly correlated variables for 25 variables if all correlations are only positive, however I need 50/50 positive/negative on average.

I make a distinction in the code for informative (relevant) & uninformative (noisy) variables.

Is there a solution to this issue?

Thank you for your time.

library("matrixcalc")
library("MASS")
library("Matrix")

p_inf <- 4
p_unf <- 21
N <- 4000

# Number of unique correlations in correlation matrix
n_inf_sample = (p_inf^2-p_inf)/2
n_unf_unf_and_inf_unf = p_inf*p_unf+(p_unf^2-p_unf)/2

Corr_inf_sample <- sample(runif(n_inf_sample,-0.3,0.3))

Corr_unf_sample <- sample(c(runif(n_unf_unf_and_inf_unf, -0.99,0.99),
                            runif(n_unf_unf_and_inf_unf,0.7,0.9),
                            runif(n_unf_unf_and_inf_unf,-0.9,-0.7)
                            ), n_unf_unf_and_inf_unf)

Corr_matrix <- matrix(ncol = p_inf+p_unf, nrow = p_inf+p_unf,0)
Corr_matrix[upper.tri(Corr_matrix, diag= FALSE)] <- c(Corr_inf_sample,Corr_unf_sample)
Corr_matrix[lower.tri(Corr_matrix, diag= FALSE)] <- t(Corr_matrix)[lower.tri(t(Corr_matrix))]
diag(Corr_matrix) <- 1

#  Mu and Sigma 
mu<- c(runif(p_inf+p_unf, -3,3)) 
sigma<- Corr_matrix

if (is.positive.definite(as.matrix(sigma))) {
  df <- as.data.frame(mvrnorm(n=N, mu=mu, Sigma=sigma))
  print("Sigma is Positive Definite")
}
if (!is.positive.definite(as.matrix(sigma))) {
  sigma2 <- nearPD(sigma)
  df<-as.data.frame(mvrnorm(n=N, mu=mu, Sigma=sigma2$mat)) 
  print("Sigma is not Positive Definite, we create PD matrix")
}

print(c("Mean absolute accorelation sigma: ", round(mean(abs(sigma)),3)))
print(c("Mean absolute correlation variables: " , round(mean(abs(cor(df))),3)))
0

There are 0 answers