Fitting a Copula from Scratch

Question

I am trying to learn about how to work with Copulas. I find that I am often getting lost in the notations and distributions, and wanted to try and solidify my understanding.

As it stands, here is my current view on the general Copula fitting procedure.

Part A) General Copula Fitting Process: Suppose we have random variables $x_1, x_2, ...x_n$ . Each of these random variables has some probability distribution $f(x_1; \theta_1), f(x_2; \theta_2), ... f(x_n; \theta_n)$. We want to create a joint probability distribution of these random variables with a certain correlation structure (e.g. assign a correlation of $\rho_{i,j}$ between every pair of $x_i$ and $x_j$) - however this might not be directly possible. Thus, we achieve this using a Copula. Using a Copula, we can simulate random numbers from this correlated joint distribution which was not possible earlier.

As I understand, this is done with the following steps:

1. Define the Marginal Distributions: For each random variable $X_i$, we have a distribution $f(x_i; \theta_i)$, where $\theta_i$ are the parameters of the distribution.

2. Generate Uniform Random Variables: Simulate random numbers from $f(x_1; \theta_1), f(x_2; \theta_2), ... f(x_n; \theta_n)$ . Then, plug these simulated numbers into the corresponding CDF's $F(x_1; \theta_1), F(x_2; \theta_2), ... F(x_n; \theta_n)$ . The outputs of these will all have Uniform Probability Distributions $u_1, u_2, ... u_n$ (since all CDF's are in effect Uniformly Distributed - this is a classic statistical property).

3. Choose a Copula: Select a copula function $C()$. This function will contain the desired correlation structure $\rho_{ij}$ between each pair of variables. Common choices include Gaussian, Clayton, Gumbel, etc. (I personally don't know the advantages/disadvantages of choosing one of these Copula functions over the others)

4. Apply the Copula: Apply the copula function to the uniform random variables. This transforms the uniform random variables into a new set of uniform random variables with the desired correlation structure. This is represented as $V_i = C(U_1, U_2, ..., U_n)$.

5. Transform Back to Original Scale: Finally, transform these new uniform random variables back to the original scale of your data. We do this by applying the inverse of the CDF of each marginal distribution to the corresponding uniform random variable. Mathematically, this is represented as $X_i = F^{-1}(V_i; \theta_i)$, where $F^{-1}(.; \theta_i)$ is the inverse CDF of the $i$-th marginal distribution.

6. Repeat Step 2 - Step 5 as many times as needed depending on how many samples you want

Part B) Example of Fitting a Copula: To test my understanding, I want to simulate from 3 binomial random variables that are correlated amongst themselves. Here is my attempt to do this using a Gaussian Copula.

1) Define all marginal distributions and correlations: $$X_1 \sim Binom(n, p_1), X_2 \sim Binom(n, p_2), X_3 \sim Binom(n, p_3)$$

$$\rho = \begin{bmatrix} 1 & \rho_{12} & \rho_{13} \\ \rho_{12} & 1 & \rho_{23} \\ \rho_{13} & \rho_{23} & 1 \end{bmatrix}$$

where $\rho_{ij}$ is the correlation between $X_i$ and $X_j$.

2) Set up a mechanism to evaluate the CDF (Cumulative Density Function) of the Binomial Distribution: (This was not initially obvious to me) :

$$ \text{PMF: } P(X=k) = \binom{n}{k} p^k (1-p)^{n-k} $$ $$ \text{CDF: } F(k; n, p) = \sum_{i=0}^{k} \binom{n}{i} p^i (1-p)^{n-i} $$

Thus, we generate 3 random numbers from the 3 Binomial Distributions [this can easily be done in many statistical computing software, e.g. rbinom() in R], and then plug these 3 random into the 3 CDFs. The 3 new numbers that will come out of the CDF will be from a Uniform Distribution.

Note: I have heard that the CDF of the Binomial Distribution is often handled using the Incomplete Beta Function.

3) Define the Copula Function: For this problem, I am using the Gaussian Copula - the Copula is the mechanism which allows to take into consideration the correlation structure. We can define the Gaussian Copula as (I used the same notation as https://en.wikipedia.org/wiki/Copula_(probability_theory)):

$$ c_R^{\text{Gauss}}(u) = \frac{1}{\sqrt{\det{R}}}\exp\left(-\frac{1}{2} \begin{pmatrix}\Phi^{-1}(u_1)\\ \vdots \\ \Phi^{-1}(u_d)\end{pmatrix}^T \cdot \left(R^{-1}-I\right) \cdot \begin{pmatrix}\Phi^{-1}(u_1)\\ \vdots \\ \Phi^{-1}(u_d)\end{pmatrix} \right) $$

• $c_R^{\text{Gauss}}(u)$ is the Gaussian copula function.
• $u$ is a vector of uniform random variables.
• $\Phi^{-1}(u_i)$ is the inverse of the cumulative distribution function (CDF) of the standard normal distribution, applied to each element $u_i$ of the vector $u$.
• $R$ is the correlation matrix.

4) Take the output of the Copula function and apply the Inverse Binomial CDF

Again, figuring out how to do this was not immediately obvious to be - but I think an approximate iterative procedure can be used to do this:

Suppose we have a binomial distribution with parameters $n$ (number of trials) and $p$ (probability of success), and we want to find the smallest integer $k$ such that the cumulative distribution function (CDF) $F(k; n, p)$ is greater than or equal to a specific point $w$. This is essentially finding the inverse of the binomial CDF at $w$.

I think we can use the following process:

• Step 1: Initialize $k$ to 0.
• Step 2: Calculate the binomial CDF $F(k; n, p)$.
• Step 3: Check if $F(k; n, p) \geq w$.
  • If true, then $k$ is the correct
  • If false, increment $k$ by 1 and go back to step 2.

A more concrete example of this:

Suppose $n=10$, $p=0.5$, and $w=0.75$. We want to find the smallest $k$ such that $F(k; 10, 0.5) \geq 0.75$.

• Start with $k=0$.
• Calculate $F(0; 10, 0.5)$. This is the probability of getting 0 successes in 10 trials, which is very small.
• Since $F(0; 10, 0.5) < 0.75$, increment $k$ to 1.
• Calculate $F(1; 10, 0.5)$. This is the probability of getting 1 or fewer successes in 10 trials, which is still less than 0.75.
• Continue this process until you find a $k$ such that $F(k; 10, 0.5) \geq 0.75$.

Thus, we would apply this inverse CDF process to each individual output from the Copula. The results after applying the inverse CDF process would be the final answer - this final answer would effectively be from a Multivariate Binomial Distribution with the desired correlation structure.

Part C) Conclusion:

• While I am still learning about why the Copula are mathematically correct (i.e. why are they able to actually do what we believe they are doing, how can we verify to make sure that they are doing what they are supposed to do correctly) .... is the process I have described correct?

Note: In the case of the Gaussian Copula, I read that there is a much more efficient procedure for doing this. However, I wanted to write the general procedure for fitting Copulas.

Answer 1

The steps to fit copula models:

1. Estimate the margins.
2. Transform the data from original to copula data; you can use the pobs function from the copula package in R.
3. Fit the most appropriate copula function to your data. You can use the copula function. Or you can fit multiple copulas, then select the best-fit model using selection criteria such as AIC.

Each copula can deal with a specific type of dependency structure.

1. Gaussian and t-student --> elliptical copulas
2. Clayton --> lower tail
3. Joe and Gumbel --> upper tails
4. Frank ---> very similar to the Gaussian, but not identical
5. BB7, BB8 ---> able to deal with two different types of dependency simultaneously (you can plot them to see the differences).

There are many other copulas. But be careful; a copula model is straightforward in continuous cases as the type of copula is unique, while it is not in the discrete case, as it is not unique.

I think this article is helpful. Copula