What is the advantage of modelling dependency using Copulas? I have been trying to understand why modeling dependency using Copulas is widely used - specifically, what are the advantages of using Copulas? 
Here is my understanding: 


*

*The marginals of the data may not be normally distributed, hence, we cannot write it in the form of a multi-variate normal distribution. Therefore, it may be difficult to write its density in a closed form. However, if we map our data using a percentile-to-percentile basis into a multivariate Gaussian distribution, we can represent this new Gaussian distribution in a closed-form. But, still what is the advantage of representing so? 


Can you please give a complete example of why modeling becomes easier using copulas? 
 A: Before I write the answer to my question, I would like to post some excellent answers already given for a similar question
First timers, usually, do not appreciate the modeling flexibility associated with a Gaussian copula vis-a-vis a Gaussian distribution. What I mean is, say, we have 10000 data points and we want to fit a distribution to these available 10k datapoints. One of the reasons why multi-variate Gaussian distributions are widely used is because of its simple closed form expression. We can easily use Maximum likelihood estimation to fit a mutlivariate Gaussian distribution to the available 10k datapoints. However, in the real world, the data is rarely multi-variate Gaussian. In such real world cases, it does not make sense to try and fit a Multivariate Gaussian distribution to the data.  Here is where Copulas offer the flexibility. 
Some theory on Copulas: 
The CDF of a bi-variate Gaussian Copula distribution is given by $C(u,v) = \Phi_{\rho}(\phi^{-1}(u),\phi^{-1}(v))$ where the $U$ and $V$ are uniformly distributed random variables, $\Phi_{\rho}$ - CDF of bivariate Gaussian with a correlation $\rho$ and $\phi$ is the CDF of univariate Gaussian.
It is straightforward to see that any continuous cummulative distribution $U = F(X)$ produces a uniformly distributed random variable $U$. The magic lies in the $\phi^{-1}(u),\phi^{-1}(v)$ parts - Copula formulation preserves the marginal distribution of the $u = F(x)$ and $v=G(y)$ 
How?  


*

*We know $\Phi_{\rho}(p,q)$ is bi-variate Gaussian distribution in $P-Q$ space. Since the joint distribution is Gaussian, the marginal distributions have to be Gaussian. Essentially, the marginal distribution of $Q$ is given by $\phi(q)$. 

*Now, if we do a non-linear transformation of the $Q$-axis in the $P-Q$ space, the distribution may no longer be Gaussian. Specifically, if we do a non-linear transformation $q = \phi^{-1}(v)$, the marginal CDF along the new transformed axis becomes $\phi(\phi^{-1}(v))$, which essentially is $v$ - a uniform random variable. Since the CDF of any continuous distribution is $v$, why not have that $ v = G(y)$. Essentially, we moved from a $P-Q$ space to a $P-Y$ space using the non-linear transform $q= \phi^{-1}(G(y))$.  


Let me demonstrate with an example.   
library(mvtnorm)
library(MASS)
library(ggplot2)
Sigma <- matrix(c(1,0,0,1),2,2)
Sigma
sim1 = mvrnorm(n = 10000, c(0,0), Sigma)
ggplot()+ geom_point(aes(sim1[,1],sim1[,2])) + labs(x = "X-Axis",y = "Y-Axis")

A simple bi-variate normal distribution. 
The same Bi-variate distribution with Y-axis scaled by a factor of 2. Note that the overall distribution is still Gaussian.

#Y-axis is dilated - meaning mulitplied by 2
ggplot()+ geom_point(aes(sim1[,1],2*sim1[,2])) + labs(x = "X-Axis",y = "Y-Axis")

Only Y-axis is given a squared transformation. Clearly this is not bi-variate Gaussian.

ggplot()+ geom_point(aes(sim1[,1],sim1[,2]^2)) + labs(x = "X-Axis",y = "Y-Axis")
# It is no longer bi-variate normal. 

Both X-axis Y-axis is given a squared transformation

ggplot()+ geom_point(aes(sim1[,1]^2,sim1[,2]^2)) + labs(x = "X-Axis",y = "Y-Axis")
# It is no longer bi-variate normal.

Copulas give the flexibility to capture these additional cases where the marginals are not necessarily normal. Using Copulas we can get a closed form expression for the distributions of all these cases. 
For e.g. For the last figure, using Copula theory, we can write the density function \begin{align}  c(u,v) &= \frac{f_\rho(\phi^{-1}(u),\phi^{-1}(v))}{f(\phi^{-1}(u)) f(\phi^{-1}(u))} & \text{Where } f_\rho & \text{ is the bivariate Gaussian distribution}  \\
&= \frac{f_\rho(\phi^{-1}(F(x)),\phi^{-1}(G(y)))}{f(\phi^{-1}(F(x))) f(\phi^{-1}(G(y)))} & \text{and } f & \text{ is the univariate Gaussian distribution} 
\end{align}
In this case, $F,G$ are the CDF of Chi-square (kind of Squared normal) distribution. 
