Programming a Bivariate Normal CDF in R

Question

I have a question regarding coding a function that contains a bivariate normal CDF in R. The function I am trying to code requires one bivariate normal CDF, that should be calculated differently depending on the observation. Specifically, depending upon the value of a certain variable, the correlation should "switch" between being positive and negative, but there should be no difference in the call.

This style of function has been coded in LIMDEP and I am trying to replicate it, but have not been able to get it to work in R. The command in LIMDEP to calculate a bivariate normal CDF is "BVN(x1, x2, r)", which explicitly requires the two variables used for calculation (x1, x2) and the correlation (r). LIMDEP uses the Gauss-Laguerre 15 point quadrature to calculate the bivariate normal CDF.

In R, it appears that two packages calculate the multivariate normal CDF. I have been trying the mnormt package (though there is the mvtnorm package as well--I do not see a major difference) that uses the Genz method, which seems to be similar, but more general than the Gauss-Laguerre 15 quadrature method used in LIMDEP (referencing the papers under ?pmnorm).

Everytime that I have tried to use the mnormt package, the command pmnorm(), requires the form: pmnorm(data, mean, varcov), which I have not been able to code for correlation switching.

Any ideas how to get this to work??

Here is an example of some trivial code to explain what I am talking about what I would like to do (except inside of a function without the for loop):

 library(mnormt)
 A <- c(0,1, 1, 1, 0, 1, 0, 1, 0, 1)
 q <- 2*A-1 
 set.seed(1234)
 x <- rnorm(10)
 y <- rnorm(10, 2, 2)

 #Need to return a value of the CDF for each row of data:
 cdf.results <- 0
 for(i in 1:length(A)){
 vc.mat <- matrix(c(1, q[i]*.7, q[i]*.7, 1.3), 2, 2)
 cdf.results[i] <- pmnorm(cbind(x[i], y[i]), c(0, 0), vc.mat) 
 }
 cdf.results

Thanks for your help!

Answer 1

It sounds like all you need is 1) to make your script into a function so it can apply to arbitrary x,y and q and 2) to get rid of the for loop. If that is the case ?function and ?apply should give you what you need.

BVN=function(x,y,q) {

  cdf.results=apply(cbind(x,y,q),1,FUN=function(X) 
       {
      x=X[1]
      y=X[2]
      q=X[3]
          vc.mat <- matrix(c(1, q*.7, q*.7, 1.3), 2, 2)
          pmnorm(cbind(x, y), c(0, 0), vc.mat)     
                    # I think you may want c(0,2) but not sure
        }
                   )
 cdf.results
}



BVN(x,y,q)

Here x,y and q are how you wrote above. Maybe you want the function to take matrix r, like in limdep? It wouldn't be much different.