My goal is to generate two variable which are correlated and one of them is heteroscedastic with regards to an grouping variable.

To create two variables with a desired correlation the common way to go is to use the cholesky decomposition. I additionally used the residuals of the orthogonal projection and standardized the variables up front to make sure the correlation stays at the desired value.

To introduce heteroscedasticity to one of the variables i tried the following:

  1. I did the cholesky decomposition for $u=(u_1,u_2), \ q[i] = c(1,1.01,1.02,1.03,...)$

    $Var(u) = \begin{pmatrix} 1 & \rho\cdot q[i] \\ \rho\cdot q[i] & q[i]^2 \end{pmatrix}$

    The reason why I choose this covariance matrix is that while the variance of $u_2$ increases with $q[i]$ the correlation stays at the desired value $\rho$.

    But this does obviously not introduce heteroscedasticity to $u_2$....

  2. A second attempt was to decomposed both variables $u_1,u_2$ with sample size n into two parts each of length $\frac{1}{2}\cdot n$: $\ u_1 = (u_{1.1},u_{1.2}), u_2 = (u_{2.1},u_{2.2})$. Now i did the cholesky decomposition for

    $Var(u) = \begin{pmatrix} 1 & 0 & \rho & \rho\sqrt{q[i]} \\ 0 & 1 & \rho & \rho\sqrt{q[i]} \\ \rho & \rho & 1 & 0 \\ \rho\sqrt{q[i]} & \rho\sqrt{q[i]} & 0 & q[i] \end{pmatrix}$

    After I generated the four variables I merged $u_{1.1},u_{1.2}$ back to one $n\times 1$ vector. Analogous for $u_2$. Testing the outcome yields that the heteroscedasticity is now where it should be as the variance of $u_{2.2}$ increases with $q[i]$. But the correlation vanishes as $q[i]$ increases because the variance of $u_2$ increases more rapidly then the covariance between $u_1$ and $u_2$

  3. A third attempt which i end up doing was to just multiply $u_2$ with a dummy-variable

    $z_2 =\begin{cases} 1\\ q[i] \end{cases}$

    This approach does work for small values of $q[i]$ (as my graphic shows) but as soon as $q[i]$ is large enough it dominates the nominator and denominator of the correlation coefficient such that the conditional variance (depending on the value of $z_2$) does no longer increase in $q[i]$.

Here is my code for case (3):

q     <- seq(1.01,10,0.1)
n         <- 100
rho   <- 0.5
sd_u1     <- numeric(0)
sd_u2.1   <- numeric(0)
sd_u2.2   <- numeric(0)
cor_u1_u2 <- numeric(0)

for(i in 1:length(q)){
u1      <- rnorm(n,0,1)
u1      <- ( u1 - mean(u1) )/sd(u1)
u2      <- rnorm(n,0,1)
z2      <- c(rep(1,0.5*n),rep(q[i],0.5*n))
u2      <- u2*z2
u2      <- as.vector( ( diag(n) - u1%*%solve(t(u1)%*%u1)%*%t(u1) ) %*% u2 )
u2      <- ( u2 - mean(u2) )/sd(u2)
z       <- cbind(u1, rho*u1+sqrt(1-rho^2)*u2) 
sd_u1[i]      <- sd(z[,1])
sd_u2.1[i]    <- sd(z[,2][1:(0.5*n)])
sd_u2.2[i]    <- sd(z[,2][(0.5*n+1):n])
cor_u1_u2[i]  <- cor(z[,1],z[,2])
plot(q,sd_u1, type="l")
plot(q,sd_u2.1, type="l", ylim=c(0,2))
plot(q,cor_u1_u2, type="l")

enter image description here

  • 2
    $\begingroup$ It would help to explain your problem in English and mathematical terms, rather than wholly as an R program. There are several reasons for this, including (1) the program might not be doing what you intend, (2) it would help readers avoid having to work through code to understand the question, and (3) in the course of having to explain the circumstances clearly, you might arrive at a solution yourself. $\endgroup$
    – whuber
    Commented Nov 8, 2012 at 14:59
  • $\begingroup$ i edited a short description on top of the main article $\endgroup$
    – Druss2k
    Commented Nov 8, 2012 at 16:34
  • $\begingroup$ Thanks. I find the notation mysterious: what is the relationship between $u$ and the $q[i]$? Obviously there are many $q[i]$, so how is it possible to specify correlations with only a two by two matrix? If you intend this as a sequence of matrices, then you have only specified a small part of the full covariance matrix (blocks along the diagonal). What should the off-diagonal blocks contain? Perhaps you could back up a little and tell us what exactly you mean by a "variable" and "heteroscedastic with respect to a grouping variable." What is your model for heteroscedasticity? $\endgroup$
    – whuber
    Commented Nov 8, 2012 at 17:27
  • $\begingroup$ Sorry about my crappy notation. By a variable I mean a random variable. In this case I got n realisations of this particular random variable. Since I got two of them I design a 2x2 covariance matrix. The relationship between $u$ and $q[i]$ is that I use a sequence $q[i]$ to increase the variance of one part of $u$ such that it couses heteroscedasitcity. Since both variables have zero mean I can describe the correlationmatrix in this way because $E(uu') = I_2$ With L equal to the cholesky decomposition ill get $E(Luu'L') = LI_2 L' = Q$ with Q the desired covariance matrix- $\endgroup$
    – Druss2k
    Commented Nov 8, 2012 at 17:40
  • 1
    $\begingroup$ Let's check that I understand, because many interpretations still seem possible. Do you wish to generate a sequence of vector-valued variables $(U_i,V_i)$ such that (1) Var$(V_i)=q_i$ (with $q_i$ predetermined) and (2) the correlation coefficient of the realized values $(u_i,v_i)$ equals $\rho$? (Please notice the distinction between the random variables and their realizations in this description.) And perhaps you want $V_i/q_i$ and $U_i$ all to be identically distributed? Maybe you only want the expectation of the correlation coefficient of $(u_i,v_i),i=1,\ldots,n$ to equal $\rho$? $\endgroup$
    – whuber
    Commented Nov 8, 2012 at 20:35

1 Answer 1


Ruscio and Kaczetow (2008) proposed a method for inducing an exact correlation into two random variables with arbitrary distributions. Both variables can be sampled from existing empirical distributions, and the marginal distributions will reflect these distributions.

They have an accompanying R code which works (I used it a couple of times). As some minor typing errors in the original code had to be corrected: here's a cleaned version.

Ruscio, J., & Kaczetow, W. (2008). Simulating Multivariate Nonnormal Data Using an Iterative Algorithm. Multivariate Behavioral Research, 43(3), 355–381. doi:10.1080/00273170802285693

  • $\begingroup$ I hope that is an answer to the question (as it is not entirely clear to me as well ...) $\endgroup$
    – Felix S
    Commented Nov 9, 2012 at 13:46
  • $\begingroup$ sry i was at my moms :) yes this was very helpfull though i did it in a different way. as i said i was looking for a way to generate two variables with a certain correlation from which one is heteroscedastic. one covariance-matrix which does the job is the following: $\begin{pmatrix} 1 & 0 & \rho & rho[l]*\sqrt{q[i]} \\ 0 & 1 & \rho & \rho*\sqrt{2}*\sqrt{1+q[i]}-\rho \\ \rho & \rho & 1 & 0 \\ \rho*\sqrt{q[i]} & \rho*\sqrt{2}*\sqrt{1+q[i]}-\rho & 0 & q[i] \end{pmatrix}$ $\endgroup$
    – Druss2k
    Commented Nov 11, 2012 at 14:01
  • $\begingroup$ Before I was using the Cholesky-Decomposition to get the desired weights but I find it much easier to use the Singular Value Decomposition instead. Though this matrix does control the correlation in each group it tends to get larger in the group $u_{1.2},u_{2.2}$ for a certain $\rho$ if $q[i]$ increases. If $\rho$ reaches 0.5 the overall correlation does does decrease too if $q[i]$ increases. I know that my descriptions are rather cryptic and im very sorry about that. If someone is really interestend in my simulation study i would really like to explained it in more detial perhaps via chat :) $\endgroup$
    – Druss2k
    Commented Nov 11, 2012 at 14:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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