Return to Answer

I will describe the most general possible solution. Solving the problem in this generality allows us to achieve a remarkably compact software implementation: just two short lines of R code suffice.

Pick a vector $X$, of the same length as $Y$, according to any distribution you like. Let $Y^\perp$ be the residuals of the least squares regression of $X$ against $Y$: this extracts the $Y$ component from $X$. By adding back a suitable multiple of $Y$ to $Y^\perp$, we may produce a vector having any desired correlation $\rho$ with $Y$. Up to an arbitrary additive constant and positive multiplicative constant--which you are free to choose in any way--the solution is

Here is working R code. If you don't supply $X$, the code will draw its values from the multivariate standard Normal distribution.

To illustrate, I generated a random $Y$ with $50$ components and produced $X_{Y;\rho}$ having various specified correlations with this $Y$. They were all created with the same starting vector $X=(1,2,\ldots, 50)$. Here are their scatterplots. The "rugplots" at the bottom of each panel show the common $Y$ vector.

If you would like to experiment, here is the code that produced these data and the figure. (I didn't bother to use the freedom to shift and scale the results, which are easy operations.)

BTW, this method readily generalizes to more than one $Y$: if it's mathematically possible, it will find an $X_{Y_1,Y_2,\ldots,Y_k;\rho_1,\rho_2,\ldots,\rho_k}$ having specified correlations with an entire set of $Y_i$. Just use ordinary least squares to take out the effects of all the $Y_i$ from $X$ and form a suitable linear combination of the $Y_i$ and the residuals. (It helps to do this in terms of a dual basis for $Y$, which is obtained by computing a pseudo-inverse. The following code uses the SVD of $Y$ to accomplish that.)

I will describe the most general possible solution. Solving the problem in this generality allows us to achieve a remarkably compact software implementation: just two short lines of R code suffice. At the end is a generalization to multiple $Y$ vectors, with working code.

Pick a vector $X$, of the same length as $Y$, according to any distribution you like. Let $Y^\perp$ be the residuals of the least squares regression of $X$ against $Y$: this removes the $Y$ component from $X,$ producing a vector orthogonal to $Y.$ By adding back a suitable multiple of $Y$ to $Y^\perp$, we may produce a vector having any desired correlation $\rho$ with $Y$ (except $\rho=\pm 1$, but then $\pm Y$ works). Up to an arbitrary additive constant and positive multiplicative constant--which you are free to choose in any way--the solution is

Here is working R code. If you don't supply $X$, the code will draw its values randomly from the multivariate standard Normal distribution.

To illustrate, I generated a vector $Y$ with $50$ components and produced various $Z=X_{Y;\rho}$ having specified correlations $\rho$ with this $Y$. They were all created with the same starting vector $X=(1,2,\ldots, 50)$. Here are their $(Y,Z)$ scatterplots. The "rugplots" at the bottom of each panel show the common $Y$ vector.

If you would like to experiment, modify this code that produced the data and the figure. (I didn't bother to use the freedom to shift and scale the results, which are easy operations.)

BTW, this method readily generalizes to more than one $Y$: if it's mathematically possible, it will find an $X_{Y_1,Y_2,\ldots,Y_k;\rho_1,\rho_2,\ldots,\rho_k}$ having specified correlations with an entire set of $Y_i$. Just use ordinary least squares to take out the effects of all the $Y_i$ from $X$ and form a suitable linear combination of the $Y_i$ and the residuals. (It helps to do this in terms of a dual basis for $Y$, which is obtained by computing a pseudo-inverse. The following code uses the SVD of $Y$ to accomplish that.)

complement <- function(y, rho, x, threshold=1e-12) {
  #
  # Process the arguments.
  #
  if(!is.matrix(y)) y <- matrix(y, ncol=1)
  if (missing(x)) x <- rnorm(n)
  d <- ncol(y)
  n <- nrow(y)
  y <- scale(y, center=FALSE) # Makes computations simpler
  #
  # Remove the effects of `y` on `x`.
  #
  e <- residuals(lm(x ~ y))
  #
  # Calculate the coefficient `sigma` of `e` so that the correlation of
  # `y` with the linear combination y.dual %*% rho + sigma*e is the desired
  # vector.
  #
  y.dual <- with(svd(y), (n-1)*u %*% diag(ifelse(d > threshold, 1/d, 0)) %*% t(v))
  sigma2 <- c((1 - rho %*% cov(y.dual) %*% rho) / var(e))
  #
  # Return this linear combination.
  #
  if (sigma2 >= 0) {
    sigma <- sqrt(sigma2) 
    z <- y.dual %*% rho + sigma*e
  } else {
    warning("Correlations are impossible.")
    z <- rep(0, n)
  }
  return(z)
}
#
# Set up the problem.
#
d <- 3           # Number of given variables
n <- 50          # Dimension of all vectors
x <- 1:n         # Optionally: specify `x` or draw from any distribution
y <- matrix(rnorm(d*n), ncol=d) # Create `d` original variables in any way
rho <- c(0.5, -0.5, 0)          # Specify the correlations
#
# Verify the results.
#
z <- complement(y, rho, x)
cbind('Actual correlations' = cor(y, z),
      'Target correlations' = rho)
#
# Display them.
#
colnames(y) <- paste0("y.", 1:d)
colnames(z) <- "z"
pairs(cbind(z, y))

complement <- function(y, rho, x, threshold=1e-12) {
  #
  # Process the arguments.
  #
  if(!is.matrix(y)) y <- matrix(y, ncol=1)
  d <- ncol(y)
  n <- nrow(y)
  y <- scale(y, center=FALSE) # Makes computations simpler
  if (missing(x)) x <- rnorm(n)
  #
  # Remove the effects of `y` on `x`.
  #
  e <- residuals(lm(x ~ y))
  #
  # Calculate the coefficient `sigma` of `e` so that the correlation of
  # `y` with the linear combination y.dual %*% rho + sigma*e is the desired
  # vector.
  #
  y.dual <- with(svd(y), (n-1)*u %*% diag(ifelse(d > threshold, 1/d, 0)) %*% t(v))
  sigma2 <- c((1 - rho %*% cov(y.dual) %*% rho) / var(e))
  #
  # Return this linear combination.
  #
  if (sigma2 >= 0) {
    sigma <- sqrt(sigma2) 
    z <- y.dual %*% rho + sigma*e
  } else {
    warning("Correlations are impossible.")
    z <- rep(0, n)
  }
  return(z)
}
#
# Set up the problem.
#
d <- 3           # Number of given variables
n <- 50          # Dimension of all vectors
x <- 1:n         # Optionally: specify `x` or draw from any distribution
y <- matrix(rnorm(d*n), ncol=d) # Create `d` original variables in any way
rho <- c(0.5, -0.5, 0)          # Specify the correlations
#
# Verify the results.
#
z <- complement(y, rho, x)
cbind('Actual correlations' = cor(y, z),
      'Target correlations' = rho)
#
# Display them.
#
colnames(y) <- paste0("y.", 1:d)
colnames(z) <- "z"
pairs(cbind(z, y))

The following is a more complete implementation for those who would like to experiment.

Another thread provides a detailed explanation of each line of code. The following is a more complete implementation for those who would like to experiment.