Simulate a Gaussian Copula with t margins

Question

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:

require(mvtnorm)
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)
plot(Z_gauss[,1],Z_gauss[,2])

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

Answer 1

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.

Answer 2

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)

Answer 3

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

THEORY

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!

PROCEDURE

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
           else:
                 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!