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.
