Problem with the formulation of a gaussian copula likelihood function

I recently got to hear about copulas which to me sounded like a nice tool to model relationships between variables. I decided to try to implement the likelihood function for a bivariate Gaussian copula with normally distributed marginals in R for use in MLE estimation or Bayesian estimation. In this case I would expect that this likelihood function would have the same likelihood as a bivariate normal distribution. This doesn't seem to be the case (as demonstrated below) and my question is:

Where do I go wrong implementing the likelihood function for the Gaussian copula with normal marginals? Am I maybe misunderstanding something fundamental?

Here is the code for my likelihood function in R:

# x1 and x2 is the bivariate data, rho is the correlation
loglike_fun <- function(x1, x2, mu1, mu2, sigma1, sigma2, rho) {
# transforming the data according to a Gaussian copula
qx1 <- qnorm(pnorm(x1, mu1, sigma1, log.p=T), 0, 1, log.p=T)
qx2 <- qnorm(pnorm(x2, mu2, sigma2, log.p=T), 0, 1, log.p=T)
qx <- cbind(qx1, qx2)

loglike <- sum(dmvnorm(qx, c(0,0), matrix(c(1, rho, rho, 1), ncol=2), log=T))
loglike <- loglike + sum(dnorm(x1, mu1, sigma1, log=T))
loglike <- loglike + sum(dnorm(x2, mu2, sigma2, log=T))
return(loglike)
}

Now if I generate bivariate normal randomly distributed variables with mu = 0, sigma = 1 and rho = 0.75 as

x <- rmvnorm(100000, c(mu, mu), matrix(c(sigma, rho, rho, sigma),ncol=2))

I would expect the following to have the maximum likelihood:

loglike_fun(x[,1], x[,2], mu1 = 0, mu2 = 0, sigma1 = 1, sigma2 = 1, rho = 0.75)
## -526812.6

But this is not the case as, for example, the following parameters consistently have a larger likelihood:

loglike_fun(x[,1], x[,2],mu1 = 0, mu2 = 0, sigma1 = 1.7, sigma2 = 1.7, rho = 0.9)
## -484621.1

Any pointers to where I go wrong are much appreciated!

Your qx1 is just the same as (x1-mu1)/sigma1, as you can check, and the same for qx2. The dmvnorm function evaluated at $(x,y)$ calculates $$\frac{1}{2\pi\sqrt{1-\rho^2}\sigma_1\sigma_2}\exp\left( \frac{-1}{2(1-\rho^2)}\left[\frac{(x-\mu_1)^2}{\sigma_1^2} + \frac{(y-\mu_2)^2}{\sigma_2^2} - \frac{2\rho (x-\mu_1)(y-\mu_2)}{\sigma_1\sigma_2}\right]\right)$$

but your function calculates (the log of)

$$\frac{1}{2\pi\sqrt{1-\rho^2}\sigma_1\sigma_2}\exp\left( \frac{-1}{2(1-\rho^2)}\left[\frac{(\frac{(x-\mu_1)}{\sigma_1}-\mu_1)^2}{\sigma_1^2} + \frac{(\frac{(y-\mu_2)}{\sigma_2}-\mu_2)^2}{\sigma_2^2} - \frac{2\rho (\frac{x-\mu_1}{\sigma_1}-\mu_1)(\frac{y-\mu_2}{\sigma_2}-\mu_2)}{\sigma_1\sigma_2}\right]\right)$$

times $$\frac{1}{\sqrt{2\pi}}\exp^{(\frac{x-\mu_1}{\sigma_1})^2} \frac{1}{\sqrt{2\pi}}\exp^{(\frac{y-\mu_2}{\sigma_2})^2}$$

which is not the same. I think your intention was to define

loglike_fun <- function(x1, x2, mu1, mu2, sigma1, sigma2, rho)
sum(dmvnorm(cbind(x1,x2), c(0,0), matrix(c(1, rho, rho, 1), ncol=2), log=T))

and then it does what you want.

x <- rmvnorm(100000, c(mu, mu), matrix(c(sigma, rho, rho, sigma),ncol=2))
> loglike_fun(x[,1], x[,2], mu1 = 0, mu2 = 0, sigma1 = 1, sigma2 = 1, rho = 0.75)
 -242503.6
> loglike_fun(x[,1], x[,2],mu1 = 0, mu2 = 0, sigma1 = 1.7, sigma2 = 1.7, rho = 0.9)
 -271719.2

Basically, you accidentally "standardized" x1 and x2 twice.

• Thanks for your explanation! I know that I can use the dmvnorm function to get the correct likelihood, but now I wanted to try to arrive at the same likelihood function using copulas instead. Would that be possible? – Rasmus Bååth May 31 '13 at 7:42
• I think you want to do the calculation given in part A of the following paper arxiv.org/pdf/1011.4997v2.pdf. The two-variable version of the formula you want appears as the first line in equation (17). I don't think the first term is actually a multivariate normal density, and also the constant is different. In other words, you need to replace the call to dmvnorm by something else. Does this help? – Flounderer May 31 '13 at 20:52
• Also, there are a few errors in the formulas in the above answer, and I notice that it's not relevant anyway, since you are using $\mu=0$ and $\sigma=1$. – Flounderer May 31 '13 at 20:56
• I guess what I want is in the paper you linked to. However, what I really wanted, and what I tried to do, was to understand how to construct a Gaussian copula using likelihood functions such as dnorm and dmvnorm. But maybe this is not possible? The reason why I started with the Gaussian copula with normal margins is that I thought it would be an easy case. If I got that working I could later extend it by using other marginal distributions... – Rasmus Bååth Jun 3 '13 at 7:20