I am looking into conducting some simulations. For those simulations I need to generate observations of variables $X_1$ and $X_2$ (suppose that both variables have mean zero) such that $X_1$ and $X_2$ have non-zero correlation but at the same time the second conditional moment $$V[X_1 |X_2=x_2]$$ does not depend on $x_2$.
What would be easy ways to do that?
EDIT. I am now starting to think that something like is not even possible.
I'm answering to this question: how to generate variate $Y$ or to correct existing variate $Y$ so that it shows correlations precisely (or near so) as the specified correlations, with one or more existing variables $X$, and at the same time $Y$ is homoscedastic with those $X$ as much as possible.
Homo/heteroscedasticity can be defined concretely in various ways. The OP defines homoscedasticity - as I interpreted their question and their comment to me - that it is the (near) zero correlation of the squared residuals of the linear dependency of $Y$ on $X$ with $X$. This is one of common technical understandings of homoscedasticity.
I'm suggesting an iterative approach which is the direct extension of the one given by me as a way to generate/correct a variate for it to have the specified correlations with the given one(s). The current task is also the task to bring a variate to some wanted correlations with existing variates, but now the wish to attain or to secure homoscedasticity is added. For example, variable $Y$ already exists or is an randomly generated "ingot", and it correlates with the $X$(s) not as precisely as we want; we'll transform it to reach the wanted correlations and simultaneously to cure an observed heteroscedasticity with them (or to insure against losing homoscedasticity).
Maybe a more elegant/short and more exact noniterative algorithm could be proposed by another respondent, yet I'm giving what I see currently as "not very bad" solution.
So, do everything in stages as described in the above linked answer, but on iterations in point 8, use, in place of this:
$$E_i[\text{corrected}]=E_i-\frac{\sum_{j=1}^m \frac{C_j X_{ij}^3}{\sum_{i=1}^n X_{ij}^2}} {\sum_{j=1}^m X_{ij}^2}$$
this modified formula:
$$E_i[\text{corrected}]=E_i-\frac{\sum_{j=1}^m \frac{C_j X_{ij}^3}{\sum_{i=1}^n X_{ij}^2} + \text{sgn}_i \sum_{i=1}^m \frac{K_j W_{ij}^3}{\sum_{i=1}^n W_{ij}^2} } {\sum_{j=1}^m X_{ij}^2 + \sum_{i=1}^m W_{ij}^2}$$
where $\text{sgn}_i$ is the sign of the currently observed $E_i$ value ($i$ is data point index); $W_j$ is the centered $X_j$ variable; $K_j$ is the sum cross-product of $\text{abs}(E)$ with $W_j$, observed at the beginning of an iteration.
You can see that it is the formula like the previous, only with term (more burden) added, to drift towards homoscedasticity. Everything else of the algorithm remains the same. The procedure cures / insures against fan-like heteroscedasticity, like that:
If more than one $X$ variables, homoscedasticity is pursued in relations of $Y$ with all of them, but on a simultaneous, regression $Y \sim X$s basis, not pairwise "marginal" basis.
There were n=400 data with variables $X_1$, $X_2$ (r=0.086) and to-be transformed variable $Y$ (the "ingot"). $Y$ is pronouncedly heteroscedastic with $X_1$ (notice on the scatterplot below). We require $Y$ to correlate with the two $X$s as $0.3$ and $0.2$, respectively, while getting rid of the heteroscedasticity as much as possible. Here's what we had and what we got after 20 iterations (I omit showing $X_2$ on the plot because heteroscedastisity were not an issue with it):
Result: correlations $0.3$ and $0.2$ were (almost) reached while heteroscedasticity of $Y$ with $X_1$ fell from $0.245$ to $0.060$.
It is possible to correct specifically symmetric heteroscedasticities such as
where data cloud is either bow tie or rhombic about the mean of $X$ variable(s). In order to transform such shape to more homoscedastic, ellipsoidal one, re-define $W_j$ in the formula to be the centered $\text{abs(centered}X_j)$ variable. This is more radical training of $Y$ for homoscedasticity.
Example (the $X$ variables were used from the previous example, $Y$ was simulated different; notice the bow tie shape heteroscedasticity):
The heteroscedasticity decreased from $0.435$ to $0.111$. More than 20 iterations might be able to yield more success. The described iterative approach (starting from the one presented in the linked answer) are more successful if n is large and if there are many unique values in the data (i.e. input variables are continuous).
