# How to write a function for the normal copula in R?

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.
• I voted to reopen this question because it implicitly asks for a computable formula for the bivariate Normal CDF.
– whuber
2 days ago

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);
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)


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.

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.
• 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.
– whuber
2 days ago
• @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.
– Aria
yesterday