# Question on GHK algorithm

After implementing the GHK algorithm, I noticed an unexpected behaviour that lead me to question the validity of the algorithm, but since this is widely used I am sure I must be missing something.

Specifically, even choosing a bivariate normal with zero-mean and symmetric bound (e.g. (-Inf,0) and (0, Inf) respectively, for the two variables), I would expect the two variables to have the same but symmetrical distribution, but this is not the case, since the variances of the two variables are different.

I initially thought this was a coding mistake but using the packages tmvnsim (GHK-based) and tmvtnorm (rejection-based) reveals the same issue.

library(tmvnsim); library(tmvtnorm)

Sigma <- matrix(c(1, .8, .8, 1), 2, 2)

x <- rtmvnorm(10000, rep(0, 2),
sigma = Sigma,
lower = c(-Inf,0),
upper = c(0, Inf),
algorithm = "rejection")

cov(x)


[,1] [,2]

[1,] 0.10147623 0.01712324

[2,] 0.01712324 0.10071812

x <- tmvnsim(10000, 2, lower= c(-Inf,0), upper= c(0, Inf),
imod=rep(FALSE, 2), means=rep(0, 2), sigma= Sigma)

cov(x\$samp)


[,1] [,2]

[1,] 0.3645155 0.05037320

[2,] 0.0503732 0.08500028

The variances of the two components are different, with the GHK one being clearly wrong as the two variables are symmetrical and they should have the same variance.

Following the steps of the GHK algorithm, the first component is always simulated as a univariate truncated normal distribution, but this also does not seem correct as in general marginal of a multivariate truncated normal distribution are not univariate truncated normal distribution, which further makes me wonder what I am missing about the algorithm.

Does anyone have any intuition for why the GHK does not give symmetric draws even if the bounds are symmetric (and the covariance matrix has 1 on the diagonal)?