Let a bivariate random vector $(X_1,X_2)\sim N(\mathbf{\mu},\Sigma)$ with correlation $\rho$. Let us apply probability integral transform (PIT) on each of the marginals to make them uniformly distributed and call the resulting random vector $(U_1, U_2)$.
Question: Given $\rho$, how can we obtain the correlation between the PITs of the marginals, $\text{Corr}(U_1, U_2)$?
I have run a simulation and noticed $\text{Corr}(U_1,U_2)$ tends to be very close to (just slightly below) $\rho$. I thought maybe this has to do with symmetry of both the normal and the uniform marginals, so I tried transforming $(U_1, U_2)$ to skew-Normal, producing $(Y_1, Y_2)$, and to exponential, producing $(Z_1, Z_2)$. The correlation of the skew-Normal marginals, $\text{Corr}(Y_1,Y_2)$, was still very close to $\rho$. The correlation of the exponential marginals, $\text{Corr}(Y_1,Y_2)$, was a noticeably further away from $\rho$. This is perhaps unsurprising given the not-so-strong asymmetry of the skew-Normal distribution and the strong asymmetry of the exponential distribution; see the scatterplots below. Correlation between the same marginal transformed into different distributions is in line with this thinking.
library(MASS)
library(fGarch)
# Generate bivariate Normal
n=1e4
set.seed(1); x=mvrnorm(n=n, mu=c(0,0), Sigma=matrix(c(1,-0.8,-0.8,1),nrow=2), tol=1e-6, empirical=TRUE)
cor(x)
# Normal to uniform via PIT
u=cbind( pnorm(x[,1]), pnorm(x[,2]))
cor(u)
# Normal to uniform to skew-Normal
y=cbind( qsnorm(u[,1], mean=0, sd=1, xi=1.5), qsnorm(u[,2], mean=0, sd=1, xi=1.5) )
cor(y)
# Normal to uniform to exponential
z=cbind( qexp(u[,1], rate=1), qexp(u[,2], rate=1) )
cor(z)
# Correlations between Normal and uniform, skew normal, exponential for the first marginal
cor(x[,1],u[,1])
cor(x[,1],y[,1])
cor(x[,1],z[,1])
# Scatterplots
par(mfrow=c(2,2),mar=c(2,2,2,1))
plot(x=x[,1],y=x[,2],main="Bivariate normal")
plot(x=u[,1],y=u[,2],main="Normal copula + uniform marginals")
plot(x=y[,1],y=y[,2],main="Normal copula + skew-Normal marginals")
plot(x=z[,1],y=z[,2],main="Normal copula + exponential marginals")
par(mfrow=c(1,1))
