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
$$X_{Y;\rho} = \rho\, \operatorname{SD}(Y^\perp)Y + \sqrt{1-\rho^2}\,\operatorname{SD}(Y)Y^\perp.$$
("$\operatorname{SD}$" stands for any calculation proportional to a standard deviation.)
Here is working R
code. If you don't supply $X$, the code will draw its values from the multivariate standard Normal distribution.
complement <- function(y, rho, x) {
if (missing(x)) x <- rnorm(length(y)) # Optional: supply a default if `x` is not given
y.perp <- residuals(lm(x ~ y))
rho * sd(y.perp) * y + y.perp * sd(y) * sqrt(1 - rho^2)
}
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.

There's a remarkable similarity among the plots, isn't there :-).
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.)
y <- rnorm(50, sd=10)
x <- 1:50 # Optional
rho <- seq(0, 1, length.out=6) * rep(c(-1,1), 3)
X <- data.frame(z=as.vector(sapply(rho, function(rho) complement(y, rho, x))),
rho=ordered(rep(signif(rho, 2), each=length(y))),
y=rep(y, length(rho)))
library(ggplot2)
ggplot(X, aes(y,z, group=rho)) +
geom_smooth(method="lm", color="Black") +
geom_rug(sides="b") +
geom_point(aes(fill=rho), alpha=1/2, shape=21) +
facet_wrap(~ rho, scales="free")
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.)
Here's a sketch of the algorithm in R
, where the $Y_i$ are given as columns of a matrix y
:
y <- scale(y) # Makes computations simpler
e <- residuals(lm(x ~ y)) # Take out the columns of matrix `y`
y.dual <- with(svd(y), (n-1)*u %*% diag(ifelse(d > 0, 1/d, 0)) %*% t(v))
sigma2 <- c((1 - rho %*% cov(y.dual) %*% rho) / var(e))
return(y.dual %*% rho + sqrt(sigma2)*e)
The following is a more complete implementation for those who would like to experiment.
complement <- function(y, rho, x) {
#
# 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) # 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 > 0, 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(cbind(z, y))[1,-1],
'Target correlations' = rho)
#
# Display them.
#
colnames(y) <- paste0("y.", 1:d)
colnames(z) <- "z"
pairs(cbind(z, y))