4
$\begingroup$

How can I write the following function for the normal copula in R?
$$ C_\theta(u, v)=\Phi_\theta\left(\Phi^{-1}(u), \Phi^{-1}(v)\right), $$ where $\Phi$ is the $N(0,1)$ cdf, $\Phi^{-1}$ is the functional inverse of $\Phi$ and $\Phi_\theta$ is the bivariate standard normal CDF with correlation $\theta$.

I have tried writing it like this but it is not good because eventually when I try to find the log likelihood it does not return any value.

pcopula.fam1 <- function(u, v, a) {
  return(exp(-(qnorm(u)^2-2*a*qnorm(u)*qnorm(y)+qnorm(y)^2)/(2*(1-a^2)))/(2*pi))
}

The Frank copula
This is how I was able to write my function for the Frank copula, but I am stuck with the normal copula described above.

pcopula.fam1<-function(u,v,a){
  - (1/a) * log(1 + ((exp(-a*u)-1) * (exp(-a*v)-1))/(exp(-a)-1))
} 
New contributor
Aria is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
$\endgroup$
1
  • 2
    $\begingroup$ I voted to reopen this question because it implicitly asks for a computable formula for the bivariate Normal CDF. $\endgroup$
    – whuber
    2 days ago

2 Answers 2

5
$\begingroup$

By definition, the bivariate CDF for variables $(X,Y)$ is

$$\Phi_\theta(x,y) = \Pr(X\le x,\ Y\le y).$$

It is the area of an infinite rectangle (with upper right vertex at $(x,y)$) weighted by the probability density. Fubini's Theorem states it can be computed as an iterated integral. For a bivariate Normal distribution, one integral will result in an error function; the second integral usually is best computed numerically. (The special case where one of the upper limits is infinite, solved at How can I calculate $\int^{\infty}_{-\infty}\Phi\left(\frac{w-a}{b}\right)\phi(w)\,\mathrm dw$, hints at the difficulty.)

Several well-established R packages offer this capability. The following solution copula.normal uses pmvnorm in the mvtnorm package to compute $\Phi_\theta.$ The inverses $\Phi^{-1}$ of the standard univariate Normal are offered in any statistical computing platform: in R they are called qnorm.

copula.normal <- Vectorize(function(u, v, rho = 0) {
  require(mvtnorm)
  pmvnorm(upper = qnorm(c(u, v)), corr = matrix(c(1, rho, rho, 1), 2))
})

You can see this is a direct translation of the mathematical formula into the corresponding function calls.

To illustrate its use, let's plot some of these CDFs.

u <- seq(0, 1, length.out = 21)
pars <- par()
par(mfrow = c(2, 2), mai = c(0.5,0.4,1,0.2))
for (rho in c(-0.8, 0, 0.5, 0.9)) 
  image(outer(u, u, copula.normal, rho = rho), 
        main = bquote(rho==.(rho)), cex.main = 2)
u <- seq(0, 1, length.out = 201);
pars <- par(no.readonly = TRUE)
par(mfrow = c(2, 2), mai = c(0.5,0.4,1,0.2));
for (rho in c(-0.8, 0, 0.5, 0.9)) 
  image(outer(u, u, copula.normal, rho = rho), 
        main = bquote(rho==.(rho)), cex.main = 2)
par(pars)

enter image description here

This looks about right, because as $\rho\to-1$ (perfect negative correlation) the plot should look like the lower Fréchet–Hoeffding bound $W$, which it approximates for $\rho=-0.8$ in the upper left, and as $\rho\to1$ (perfect positive correlation) the plot should look like the upper Fréchet–Hoeffding bound $M$], which it approximates for $\rho=0.9$ in the lower right, passing through the independence copula $C(u,v)=uv$ for $\rho=0$ (upper right) whose level curves are perfect hyperbolae.

Because this example computes $4\times 201\times 201 = 161,604$ of these integrals, it takes a few seconds. In practice, if you need many such calculations, you can speed it up by precomputing values on a coarse grid and interpolating.

$\endgroup$
0
$\begingroup$

While I am not sure if what I have done is correct, I was able to write a function for $$ C_\theta(u, v)=\Phi_\theta\left(\Phi^{-1}(u), \Phi^{-1}(v)\right), $$

pcopula.fam1 <- function(u, v, a) {
pbinorm(qnorm(u), qnorm(v), mean1 = 0, mean2 = 0, var1 = 1, var2 = 1, cov12=a)
}

This function works in R and has given me a value for my log-likelihood. Thank you for your help.

New contributor
Aria is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
$\endgroup$
2
  • $\begingroup$ pbinorm is not a standard R function. Where does it come from? Regardless, assuming it performs the same calculation as I describe in my answer, your solution is identical to the one I posted. $\endgroup$
    – whuber
    2 days ago
  • $\begingroup$ @whuber Hello! the pbinorm is the cumulative distribution function for the bivariate normal distribution distribution and it is from the R package VGAM. Thank you for your help. Your code also gives me more or less the same answer. $\endgroup$
    – Aria
    yesterday

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.