How can I generate (two) gaussian variables with and without heteroskedasticity simultaneously I'm trying to generate two gaussian variables: one that has heteroskedasticity using the rnorm function like this answer. I have trying to use rmvnorm but it generates the same sigma simultaneously.
 A: I think a simple example would be to simulate using a model like this (assuming the $i^{th}$ value of variable $Y$ is a quadratic heteroscedastic function—such as that pictured in the link you provided, copied here—of the $i^{th}$ value of variable $X$):

$Y_i = \beta_{0i} + \beta_{Xi} X_{i} + \varepsilon_{i}$, where
$X \sim \mathcal{N}\left(\mu_{X}, \sigma^{2}_{X}\right)$, 
$\beta_{0i} = \beta_{0} + \varepsilon_{0}$,
$\beta_{Xi} = \beta_{X} + \varepsilon_{X}$, and
$\varepsilon_{i} = \left[\begin{array}{c}
\varepsilon_{0}\\
\varepsilon_{X}\\
\end{array}\right] \sim \Omega \left(0,\begin{array}{cc}
\sigma^{2}_{\varepsilon_{0}} &\\
\sigma_{\varepsilon_{0} \varepsilon_{X}} & \sigma^{2}_{\varepsilon_{X}}
\end{array}\right)$.
This will produce values of $Y$ with a mean of $\beta_{0} + \beta_{X}X$, and with a variance of $\sigma^{2}_{\varepsilon0} + 2\sigma_{\varepsilon_{0} \varepsilon_{X}}X + \sigma^{2}_{\varepsilon_{X}}X$ at a given $X$.
To simulate draws of $Y$:


*

*Simulate draws of $X_{i}$, given $\mu_{X}$ and $\sigma^{2}_{X}$.

*Simulate of $\varepsilon_{i}$, given $\sigma^{2}_{\varepsilon0}$, $\sigma_{\varepsilon_{0} \varepsilon_{X}}$, and $\sigma^{2}_{\varepsilon_{X}}$, and using $X_{i}$.

*Finally, simulate draws of $Y_{i}$ using $\beta_{0}$, $\beta_{X}$, $X_{i}$, and $\varepsilon_{i}$
