The task is the following:

Given is $Z_1,...Z_{50}$ different hypothetical assets.

Each $Z_k \sim t_3$ with standard deviation $\sigma=0.01$ and $\tau(Z_i,Z_k)=0.4$ for $j\neq k$.

I want to simulate from the distribution $(Z_1,...,Z_{50})$ assuming that it has a Gaussian Copula with $t_3$ margins.

For what I understand, from the kendals tau $\tau$ I can estimate the correlation parameter $\rho=sin(\pi\tau/2)=0.5878$. From this I can determina the covariance matrix when I do want to do the simulation. Here is the code:

N <- 1000
tau <- 0.4
correlation <- sin(pi*tau/2)
std <- 0.01

S <- matrix((std^2)*correlation,50,50
for (i in 1:50) {

        S[i,i] <- std^2


gauss <- rmvt(N, sigma = S, df = 3)
U_norm <- pnorm(gauss)
Z_gauss <- qt(U_norm,df=3)

However, this method is not correct. I cannot see any plot that resembles a gaussian copula. Any thoughts, corrections?

  • $\begingroup$ Could you elaborate on what you mean by "cannot see any plot that resembles a gaussian copula"? I am puzzled by this because (1) there are no plotting commands or images in your post and (2) this copula exists in 50 dimensions. I am also curious why in your penultimate line pnorm(gauss) you do not refer to S. Incidentally, typos and missing commands (such as one that defines rmvt) in your code will make it impossible to execute. $\endgroup$
    – whuber
    Dec 3, 2015 at 15:36
  • $\begingroup$ rmvt is for generating from multidimensional t distribution. I added the library in the top. Also, added plotting code. $\endgroup$
    – Elekko
    Dec 3, 2015 at 15:48
  • $\begingroup$ Great, thank you. May I direct your attention to the one part of my comment you have not yet addressed? It may hold the key to the problem: why are you applying pnorm, which is the distribution function of a standard Normal variate, when gauss contains variables whose marginal distributions are given by the nonstandard covariance matrix S? $\endgroup$
    – whuber
    Dec 3, 2015 at 17:44
  • $\begingroup$ Not sure... how should it in that case be done? I use pnorm for since I want to create a Gaussian Copula, however, i just checked that the histogram of U_norm is NOT uniformly distributed. Any help what to do? $\endgroup$
    – Elekko
    Dec 3, 2015 at 17:50

3 Answers 3


There's an R package called "copula" that will let you do exactly this.

The process goes:

  1. Specify a copula

  2. Specify the population distribution, including whatever marginals you want. From the documentation: "A user-defined distribution, for example, fancy, can be used as margin provided that dfancy, pfancy, and qfancy are available."

  3. Generate samples from that multivariate distribution.

For you, you would specify a Gaussian copula in step 1 and then say that you want t-distributed marginals in step 2.

# Step 1
my_copula <- normalCopula(0.8)

# Step 2
my_population <- mvdc(my_copula, c("t","t"),list(t=3,t=3))

# Step 3
my_sample <- rMvdc(1000,my_population)

Caveat: I don't have access to this package right now, so I can't swear that this will compile, though it gives the gist of what to do.


I believe the approach you're attempting to implement would proceed as follows:

  1. Generate sample from a multivariate normal with the desired correlation matrix

  2. transform each margin from normal to uniform (in R, pnorm is suitable), to obtain a sample from the required copula

  3. transform each margin from uniform to the desired $t$ distribution (in R, qt will do that)

  • $\begingroup$ +1. Watch out, though: qt will produce infinite values especially with small degrees of freedom. The extreme t tails will make it practically impossible to gauge whether the resulting multivariate distribution follows a Gaussian copula. One should study the distribution of the uniform versions (in 2) rather than the distribution of the t variates (in 3). $\endgroup$
    – whuber
    Dec 4, 2015 at 15:06
  • $\begingroup$ @whuber the question asks for t-margins, so I don't see how to avoid attempting some transform from normal to $t_3$; of course one might try to implement a more direct approach there. However, I just generated 10 million values from a uniform and transformed via qt without getting a single non-finite value returned. How often are you seeing it? $\endgroup$
    – Glen_b
    Dec 4, 2015 at 20:39
  • $\begingroup$ It looks like making use of lower.tail=FALSE in qt with a check for being in the upper half will make it an event that would happen less than once in a googol times for each tail. When I get back I'll look at whether to lay that out in detail. $\endgroup$
    – Glen_b
    Dec 4, 2015 at 20:50
  • $\begingroup$ With df=3, I was seeing it relatively frequently. I didn't check, but it was on the order of a percent or so. I was using a minor variation of the OP's code, corrected to use the proper SD and to generate 6-D vectors rather than 50-D vectors (which are much harder to examine in detail!) $\endgroup$
    – whuber
    Dec 5, 2015 at 0:21
  • $\begingroup$ @whuber But nothing so complex is required to check the effect of qt -- each marginal will be uniform and we're only transforming the margins, so we only have to check the frequency with which qt(runif(...),3) has a problem. $\endgroup$
    – Glen_b
    Dec 5, 2015 at 6:27

Here I will try to answer this. I am using Python code for ease.


First, lets explicitly define your correlation matrix $$R^{50\times50} = \begin{bmatrix} .01&.4 &\cdots&.4 \\ .4 & \ddots& &.4 \\ \vdots&&\ddots &\vdots\\ .4 &\cdots &\cdots &.01 \end{bmatrix}$$ Then, since we are using a Gaussian copula to model the joint distrubtion of $z_1,...,z_{50}$, we have: $C_{\textit{Gauss}}(u_1,...,u_{50}) = Pr[U_1 \leq u_1, ..., U_{50}\leq u_{50}] = \Phi_R(x_1, ..., x_{50})$, where $\Phi_R$ is the multivariate normal CDF with correlation matrix $R$ and mean $\mathbf{0}$, $x_i$ is a Gaussian standard normal random variable, and $\Phi(x_i) = u_i$ is the standard normal CDF. Often this copula is written as

$$C_{\textit{Gauss}} (u_1,...,u_{50}) = \Phi_R(\Phi^{-1}(u_1),...,\Phi^{-1}(u_{50})),$$ where $\Phi^{-1}$ is the inverse standard normal CDF, thus $\Phi^{-1}(u_i)$ gives us a standard normal random variable $x_i$!

Lastly, we note that $u_i = F(Z_i)$, where $Z_i \sim t_3$ is your t-distributed random variable and $F$ is the CDF of your t-distrubted random variable. Thus, with a Gaussian copula, we first calculate the CDF of your random variable ($Z_i$ in this case), which gives $u_i$. We then take the inverse standard normal CDF of $u_i$, which gives us a standard normal random variable $\Phi^{-1}(u_i) = x_i \sim \mathcal{N}(0,1)$. We do this for each $Z_i$. Then after getting all the standard normal variables $x_1,...x_{50}$, we look at the CDF of their joint distribution. Because each $x_i$ is a standard normal random variables, the joint distribution of all $x_i$'s is a multivariate normal distribution! Thus, the CDF of this joint distribution is $\Phi_R$, the reason we can use a correlation matrix as our covariance matrix is because each $x_i$ has mean 0, and variance 1, so the covariance matrix is already normalized to 1!


With all of this in mind. We can now write the procedure for sampling from this copula.

  1. Sample $u_1, ..., u_{50}$ from your Gaussian copula. $$u_1, ..., u_{50} \sim C_{\text{Gauss}}(u_1,...,u_{50}).$$ Here are the steps to do this as your copula is Gaussian

    a. Since $C_{\text{Gauss}}(u_1,...,u_{50}) = \Phi_R(x_1, ..., x_{50})$, as described above, first sample $x_1,...,x_{50}$ from a joint normal multivariate distribution with mean $\mathbf{0}$ and correlation $R$ $$ x_1,...,x_{50} \sim \mathrm{Multi}(0,R).$$ There are many Python and R packages that can do this quickly and easily.

    b. Once you have $x_1,...,x_{50}$, take the standard normal CDF of each $x_i$ $$ \Phi(x_1),...,\Phi(x_{50}) = u_1,...,u_{50}.$$ This will give you your uniform marginals $u_1,...,u_{50}$

  2. Now that you have $u_1,...,u_{50}$, note for each of your t-distributed $Z_i$, that $F(Z_i) = u_i$, where $F$ is the CDF of the t-distribution. Thus to get $Z_1, ...,Z_{50}$ simply take the inverse of the CDF for each $u_i$. $$F^{-1}(u_1), ... , F^{-1}(u_{50}) = Z_1, ..., Z_{50}$$

Implementation Here is the implementation of the above steps in Python

'''Simulate from a Gaussian Copula with t-margins
variance = .01
variance amongst covariates (correlation)  = .4
import numpy as np
import scipy.stats as ss

n = 1 # The number of samples
p = 50 # The number of covariates
# the mean of the multivariate Gaussian CDF is zero
mu = np.zeros(p)
The correlation matrix R with R_ii = .01, R_ij = .4
correlation_matrix_r = np.zeros((p,p))
for i in range(n):
     for j in range(n):
           if (i != j):
                 correlation_matrix_r[i,j] = .4
                 correlation_matrix_r[i,j] = .01

# Step 1: sample u_1, ..., u_50 from copula

# a. Get standard normal random variables x_1,...,x_50 ~ N(0,1)
x_is = np.random.multivariate_normal(mu, correlation_matrix_r, size = (n,p))

# b. Get the uniform marginals u_1,..,u_50 from the standard normal random variables
# Do this by taking the standard normal CDF for each x_i
# ss.norm(x_is) will take the the standard normal CDF of each element in x_is
uniform_marginals_u = ss.norm(x_is)

# Step 2: The invese t-distrubtion CDF of each uniform marginal to get your samples
# F^{-1}(u_1),...,F^{-1}(u_50) = Z_1, ...,Z_50

# get t-distrubtion object from scipy, with specified parameters
t_3 = ss.t(df = 3)
t_dist_samples_z = t_3.ppf(uniform_marginals_u) #ppf is the inverse cdf 

This code will give you $n$ samples of $Z_1,..Z_{50}$. $$\begin{bmatrix} \{(Z_1,...Z_{50})_1 \\ \vdots \\ (Z_1,...Z_{50})_n\end{bmatrix}$$

Hope this helped!


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.