I am trying to use the GHK simulator to estimate the probabilities $F(\mathbf{x} > k\mathbf{a})$ that the values of a high dimensional ($n>1000$), correlated random vector $\mathbf{x}$ will exceed some threshold vector $\mathbf{a}$, linearly scaled by $k$ with $k \rightarrow \infty$. The GHK approach seems simple enough (e.g. here), however the results of my MATLAB implementation are pretty disappointing.
For my scenario I know that for large values of $k$, $F$ will be fairly small (e.g. $<10^{-6}$). (I can verify this using a brute force Monte-Carlo approach.) However, after some good initial estimates for small $k$, the probabilities returned by GHK simulator begin fluctuating around some value (e.g. $ \sim 10^{-3}$) and eventually drop to zero (as the calculated CDFs become ones).
My questions are:
Is the GHK simulator fundamentally unsuitable for use with high dimensional data, or is there some workaround which avoids numerical problems?
How many sample vectors should I be generating? Currently I use 1000 samples for each value of $k$. Using a larger number improves the results slightly, however then, in terms of computation time, the advantage of the GHK approach compared to a brute-force MC approach is negated.
Is there a more suitable approach for this type of problem?