Let's say I want to produce a series of random numbers that (i) follow a Laplace (or Cauchy) distribution and (ii) are correlated with a correlation rho. In the case of normal samples, if $r_{1}\sim U(0,1)$ and $r_{2}\sim U(0,1)$, then $ z=\rho r_{1}+\sqrt{1-\rho^{2}}r_{2}$ are correlated with correlation $\rho$.
But how does one go about doing the same for other distributions?
A general approach to generating correlated random variables with given marginals is to use copulas. The basic idea is to define the cdf of the joint distribution as $$F(x_1,\ldots,x_n) = C(F_1(x_1), \ldots, F_n(x_n))$$ for a suitable function $C$ which should be a joint cdf with uniform marginals on the unit cube. There are all kind of copulas, the linked Wikipedia article shows some common ones.
While Pearson correlation is not well defined in this case, because it depends on the marginal distribution (and might not even exist), the rank correlations of the generated variables are defined by the copula function. Specifically for a bivariate copula, Spearman's rho $r_S=12\int\int [C(u,v)-uv]du dv$, and Kendall's tau $\tau=4\int\int C(u,v)dC(u,v)-1$. Since most of the common copula families have a free parameter, you can adjust it to get the desired rank correlation.
Once you selected a copula, you can generate your sample using the following steps:
Step 1 is the most difficult, but it is known for the common copulas. For example, the gaussian copula would require generation from the multivariate normal distribution. Marius Hofert, Martin Maechler (2011) Nested Archimedean Copulas Meet R: The nacopula Package. Journal of Statistical Software, 39(9) describes the logic of generation from some copulas which is also implemented in an R package.
