All this is easier with a theory of causality. For example, let's use here the Structural Causal Models (which includes the Potential Outcomes) approach.
A Structural Causal Model (SCM) is triplet $M = \langle V, U, F\rangle$ where $U$ is a set of exogeneous variables, $V$ a set of endogenous variables and $F$ is a set of structural equations that determines the values of each endogenous variable. The structural equations are in the sense of assignments not equalities. For example, consider the simple structural equation $Y = X^2$. This is meant to be read $Y\leftarrow X^2$, in the sense that if I experimentally set $X=2$ then this causally determines the value of $Y = 4$ but experimentally setting $Y = 4$ does nothing to $X$. The asymmetry is important/fundamental in causality: rain causes the floor to be wet, but making the floor wet does not cause rain.
So our causal model can be thought as functional relationships among variables and we are considering these relationships as autonomous. You can think of it as an idealized representation of the real world, where the variables $V$, the endogenous variables, are what we choose to model, and the variables $U$ are the aspects we chose to ignore. Since we chose not to model the $U$, what we usually do is to represent our ignorance about $U$ with a probability distribution $P(U)$ over the domain of $U$, giving us a probabilistic SCM which is pair $\langle M, P(U) \rangle$. Notice this means that causal relationships are ultimately functional relationships, therefore causal relationships may or may not translate to specific probabilistic dependencies. Finally, every causal model can be associated with a directed (acyclic) graph $G(M)$.
Hence, one way to simulate from a probabilistic causal model is by specifying: (i) the endogenous variables $V$ you are going to model; (ii) the exogenous variables $U$ which are usually the "disturbances", along with their joint probability distribution; and, (iii) the (causal) structural relationships among the variables. It might be easier to start this process qualitatively by first drawing the causal DAG with the main features that you want to illustrate and then add the details of the simulation (functional forms) later.
To see how this can be easily done in practice, let's simulate a simple causal model that illustrates simpson's paradox (for more see Pearl). Suppose our model $M$ is given by the following causal DAG, where the variables in parenthesis are "unobserved" and each variable has an associated exogenous disturbance $U$ which is omitted for convenience:
More specifically we will assume the following structural equations $F$:
$$
\begin{aligned}
W_1 &= U_{W_1}\\
W_2 &= U_{W_2}\\
Z &= W_1 + W_2 + U_{Z}\\
X &= W_1 + U_{x}\\
Y &= X+ 10W_2+ U_{y}\\
\end{aligned}
$$
Finally, assume all disturbances in $U$ are independent standard normal random variables. Now it's easy to simulate from our causal model. In R for instance:
rm(list = ls())
set.seed(1)
n <- 1e5
w1 <- rnorm(n)
w2 <- rnorm(n)
z <- w1 + w2 + rnorm(n)
x <- w1 + rnorm(n)
y <- x + 10*w2 + rnorm(n)
This example is interesting because if you run the regression $Y \sim X$ you get $1$:
lm(y ~ x)
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept) x
0.01036 1.00081
But if you further "control" for $Z$, which is a pre-treatment variable correlated with both $X$ and $Y$ --- and some people still erroneously would say it's a confounder --- you will actually see a sign reversal of the estimate and get $-1$:
lm(y ~ x + z)
Call:
lm(formula = y ~ x + z)
Coefficients:
(Intercept) x z
0.00845 -1.01127 4.00041
In this example, since we simulated the data, we know the true causal effect is $1$ which is captured by the first regression. But you can only know that if you know the true causal structure. There's nothing in the data itself that tells you which one is the correct answer. Hence if you simulate this and give to a researcher only the variables $x$, $y$ and $z$ he can't tell the right answer just from looking at the correlations. If you want further play/simulate causal models with multi-stage Simpson's paradox reversals, you can check it here.
Simulating correlations/dependencies
To simulate correlations/dependencies you can take a similar approach. You can simply create a causal model that gives you the correlations/dependencies you want (adding latent variables if needed), simulate as above and the resulting data will have the desired correlations/dependencies. To make it easier, you can start by drawing the causal DAG (bayesian network) and read from the graph if the desired conditional dependencies/independencies are implied by your model. After that you might think of specific functional forms to get other quantitative aspects that you want. Notice that several models with different causal interpretations can give you the same correlations.
