I would like to generate correlated data draws across two unrelated distributions (or across two empirical datasets that are from unknown distributions). The correct way to do this is to use a copula. I'm not a copula expert.

Mechanically, my steps are as follows:

  1. Do multivariate normal draws with a fixed correlated structure.
  2. Map the resulting draws to their quantiles.
  3. Use these quantiles on the CDF (or eCDF) of the desired target distribution
  4. Take a rank correlation of the resulting data to verify the correlation structure is preserved. (Simulate this many times to ensure that finite sample bias is zero)

I wrote a bunch of custom code to do this for my use case, only to discover that I had a slight (3-5%) attenuation bias. Suppose I was drawing for 2 variables with correlation rho=0.6, I'd get data correlated at 0.58. Suppose rho=0.2, I'd get 0.19. Suppose rho=-0.2, I'd get -0.19. The magnitude of the bias is multiplicative and setting rho=1 or rho=0 worked the way you'd expect (1/0 observed correlation)

Figuring there was something wonky with my code, I tried using the copula package in R to do it parametrically:

# Load copula library

# Normal copula, rho=0.6
copula = normalCopula(0.6, dim = 2, dispstr="ex")

# Output distributions: uniform, normal, parameterized
decl = mvdc(copula = copula,
     margins = c("unif", "norm"),
     paramMargins = list(list(min=1, max=1000),
                         list(mean=75, sd=8)))

# Run 1,000 simulations of 1,000 observations
results = replicate(1000, 
          rMvdc(1000, decl))

# Get the average correlation coefficient across all 1,000 sims
summary(t(apply(results, 3, cor, method="spearman")))

The result is 0.58 correlation.

Before writing this question, I did a search on whether this is a known property of copula. This article suggests it is the case -- under "Using Rank Correlation Coefficients":

The correlation parameter, $\rho$, of the underlying bivariate normal determines the dependence between X1 and X2 in this construction. However, the linear correlation of X1 and X2 is not $\rho$. For example, in the original lognormal case, a closed form for that correlation is: $cor(X1, X2) = \frac{e^{\rho\sigma^2} - 1}{e^{\sigma^2} - 1}$ which is strictly less than $\rho$, unless $\rho$ is exactly 1. In more general cases such as the Gamma/t construction, the linear correlation between X1 and X2 is difficult or impossible to express in terms of $\rho$, but simulations show that the same effect happens.

Well, sugar. I'm in the "more general case" and I would like a closed form for the correlation so I can develop a correction to fix this. My end goal is writing a method in a package for simulating particular data in R (think wakefield but a little different). It's not the end of the world if this is a known property that I cannot work around, but I would at least like to understand a little bit better what's going on to be able to document the behaviour.

Can anyone help intuit what's going on here or point me to a discussion of this phenomenon or, even better, let me know what kind of correction I can apply to resolve this given arbitrary marginal distributions and/or empirical distributions?

Thanks in advance.


Your Answer

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

Browse other questions tagged or ask your own question.