# How can I generate random numbers from any given copula?

Suppose that I have a 2-dim copula function C(x_1,x_2).

How can I generate bivariate numbers from this copula?

For specific types of copulas, I can use 'rCopula' function of 'copula' package in R. But I have no idea what to do if I have an arbitrary copula function.

For a copula that corresponds to a known multivariate distribution, you can simulate from that distribution and then make the margins uniform (e.g. Gaussian copula, t-copula).

More generally if you can work out the conditional (either $C(u|v)$ or $c(u|v)$), you can simulate from a uniform for $V$ and then from the conditional, perhaps via inverse-cdf (if you know $C(u|v)$) or perhaps via say accept-reject (maybe an adaptive accept-reject, some version of ziggurat, etc, if you know $c(u|v)$).

In the case of bivariate Archimedean copulas, following Nelsen (1999) or Embrechts et al., (2001), we have a mechanism for then generating from them as follows. Suppose $(U_1,U_2)$ has a two-dimensional Archimedean copula with generator $\phi$. Then:

1. Simulate two independent $U(0,1)$ random variables, $v_1$ and $v_2$

2. Set $t=K_C^{-1}(v_2)\,$, where $K_C(t)=t-\phi(t)/\phi'(t)$

3. The desired simulated values are $u_1=\phi^{-1}(v_1\,\phi(t))$ and $u_2=\phi^{-1}((1-v_1)\phi(t))$.

There are other methods; for example in some cases it might sometimes be practical to do some version of bivariate accept-reject, say, or via transformation to some convenient bivariate distribution on which accept-reject might be applied.

• A paper which outlines this method for the D-dimensional case using approximations is here: ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=4419639 Title: "Analysis and Generation of Random Vectors using Copulas" Authors: Johann Christoph Strelen and Feras Nassaj Dec 2 '15 at 14:09

"

require(mvtnorm)
S <- matrix(c(1,.8,.8,1),2,2) #Correlation matrix
AB <- rmvnorm(mean=c(0,0),sig=S,n=1000) #Our gaussian variables
U <- pnorm(AB) #Now U is uniform - check using hist(U[,1]) or hist(U[,2])
x <- qgamma(U[,1],2) #x is gamma distributed
y <- qbeta(U[,2],1,2) #y is beta distributed
plot(x,y) #They correlate!


"