Programming two-dimensional normal CDF in R

Question

I have a question about coding a function that contains a two-dimensional normal CDF in R. The function that I am trying to encode requires one two-dimensional normal CDF, which must be calculated differently depending on the observation. In particular, depending on the value of a certain variable, the correlation should “switch” between positive and negative, but there should be no difference in the call.

This style of the function was encoded in LIMDEP, and I am trying to replicate it, but could not get it to work in R. The LIMDEP command to calculate a two-dimensional normal CDF is "BVN (x1, x2), r)", in which two variables are explicitly required used to calculate (x1, x2) and correlation (r). LIMDEP uses a 15-square Gauss-Laguerre quadrature to compute a two-dimensional normal CDF.

In R, it seems that two packets compute a multidimensional normal CDF. I tried the mnormt package (although there is the mvtnorm package, but I don’t see a significant difference), which uses the Genz method, which seems similar, but more general than the Gauss-Laguerre 15 quadrature method used in LIMDEP (link to documents under? Pmnorm) .

Every time I tried to use the mnormt package, the pmnorm () command requires a form: pmnorm (data, mean, varcov), which I could not code to switch the correlation.

Any ideas how to make this work?

Here is an example of some trivial code explaining what I'm saying about what I would like to do (except for an internal function without a 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 

Thank you for your help!

Answer 1

It seems like all you need is 1) to make your script function so that it can be applied to arbitrary x, y and q and 2) to get rid of the for loop. If so ?function and ?apply should provide you with 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 what you wrote above. Maybe you want the function to take the matrix r, as in limdep? That would be completely wrong.