4
$\begingroup$

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!

$\endgroup$
2
+50
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Rasmus Bååth May 31 '13 at 7:42
  • $\begingroup$ 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? $\endgroup$ – Flounderer May 31 '13 at 20:52
  • $\begingroup$ 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$. $\endgroup$ – Flounderer May 31 '13 at 20:56
  • $\begingroup$ 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... $\endgroup$ – Rasmus Bååth Jun 3 '13 at 7:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.