Numerical instabilities basically occur whenever the probability mass of a conditional distribution is largely concentrated on a very tiny intervall, for example in the lower of upper corner. This means, that points outside of the interval basically have no likelihood of occurrence, so that the conditional cdf will have horizontal vicinities. Inverting these horizontal vicinities would lead to a violation of the rule of functions that each $x$ value must be mapped to a unique $y$ value. In these cases the inverse function simply picks one of all possible underlying $x$ values, and all others can never be achieved through the inverse function. For simulation purposes this should not make a difference, since all of these points should never appear anyways.

Let's look at a couple of examples that threw some errors during unit testing.

BB7

The following BB7 copula exhibits a horizontal $v$-function in a local area for given $u_{1}=0.1$.

This way, different $u_{2}$ values very close to 1 get mapped to the same value:

As we can see, the problem occurs since values in upper quantiles have a likelihood of almost 0:

This of course causes problems with inversion, as the inverse $v$-function in principle would have multiple function values for a single $u_{2}$ value. In this case, the function gets assigned a single value for input 1.0:

Remark: if two $x$-values are too close to each other, Winston does show only one of them. Hence, simply attaching 1.0 to testGrid did not show the value for 1.0, so that I needed to replace the last grid point.

Hence, although at least all values above 0.99 get mapped to 1 in the $v$-function, inverting this operation will always lead to a value of 1. Hence, $v^{-1}(v(x))\neq x$ for these values, but will always return value 1:

Hence, a value $x$ with basically no probability of occurrence gets mapped to another value that should never occur.

Through this numerical instability, some values will never occur during simulation. In this case, values above 0.99 are impossible to simulate.

In general, however, it should not matter much, as these values would have a too small likelihood anyways.

BB1

We can get the same problem for the BB1 copula if we plug in values of highly improbable regions.

The conditional distribution gives rise to a highly degenerated cdf:

This way, $v^{-1}(v(x))$ will again map highly unlikely $x$ values to somewhere else. In this case, for $u_{1}=2.0\times 10^{-5}$ a value of $u_{2}=0.00098$ is already too unlikely and will be mapped all the way to a value of 1:

Only incredibly small numbers get a $v$-function value unequal to 1:

Looking a $v^{-1}$ we can see that almost all values will get mapped very close to 0, while only a very tine interval next to 1 will get mapped to different values:

Although we can not see this in the picture, the function exhibits discontinuity, as there are values close to 1 that can never occur.

Still, simulation does not lead to unreliable results. In principle, some values in upper quantiles now have a probability of occurrence of 0. However, we should not observe them anyways.

BB6

For the BB6 copula, the unit tests brought up a much worse numerical instability. Here, the $h$-function does never reach a value of 1:

In this case, plugging values larger than this threshold into the inverse function will throw an error in the underlying C++ root searching algorithm and leads Julia to crash.

This most likely is an artifact of shifting values away from the unit range bounds. The conditional distribution seems to have a lot of probability in the interval very close to 1 which can not be reached if $h(0.9999999999999|0.999)$ in fact is evaluated at $h(0.99999|0.999)$.

Even more suspicious becomes the behaviour for conditional distributions in more extreme quantiles:

Here the $h$-function exhibits non-monotonic behaviour! One reason might be that the pdf function exhibits infinite values: