Solved – Problem with the formulation of a gaussian copula likelihood function

copulalikelihoodmultivariate analysisnormal distributionr

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!

Best Answer

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)
[1] -242503.6
> loglike_fun(x[,1], x[,2],mu1 = 0, mu2 = 0, sigma1 = 1.7, sigma2 = 1.7, rho = 0.9)
[1] -271719.2

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

Related Question