So there are oscillation phenomena that depend on large signal behavior and nonlinearity (e.g. relaxation oscillators), but the point here is that the current tools fail to predict ordinary, small signal instability, ones that can be easily seen in an s-plane plot if you work out the transfer function by hand. Suppose you have a feedback loop somewhere in your circuit that at a particular frequency has a loop gain of 2, in other words Vout = Vin + 2*Vout. The solver will find that for Vin = 1, Vout = -1, and completely miss the instability. Why is this? When the circuit was linearized at a particular frequency, an assumption was made which is that the circuit is already in steady state (its internal state and excitations set up in such a way that the output is a steady sinusoid). Since this assumption is broken, everything that follows is rubbish. In the s-domain it's as if you excite the circuit with a impulse and then take the dot product of the response with various decaying and growing sinusoids, which will catch the oscillation when the decay factor is close to the oscillation growth factor.
In spice all circuits are built from basic elements (passives and dependent sources) which all have trivial s-domain transfer functions. Even a transmission line's S parameters are easily calculated in the s-domain instead of frequency domain (e.g. instead of S21 being e**(jwL) it would be e**sL). In a linear simulator like RFSim99 it's the same idea, basic elements and transmission lines, except for device models that are specified by S parameter files.
How do you derive s-domain S parameters from frequency domain measurements? First observation is that as long as the individual circuit element (transistor or MMIC) itself isn't oscillating, it can be completely described by its S parameters as long as there are enough frequency points that go down to near DC. A linear two port network is nothing more than a LTI system with 2 inputs and 2 outputs. The S parameters simply describe the transfer function of each combination of (input, output). If you transform the S parameters to the time domain, what you get is the impulse response of the system for each pair of (input, output), or in other words time domain S parameters. Given the impulse response it's easy to just use the Laplace transform to get s-domain S parameters for any arbitrary s. (btw "s" should be thought of as a combination of frequency and decay/growth factor.)
It's still possible with S parameters that don't go down to near DC. If you think about what the Laplace transform for one particular s is, it's the dot product between a decaying complex sinusoid and the impulse response of the system. A dot product can be taken in the frequency domain as well as time. However, the Fourier transform of the decaying/growing sinusoid doesn't exist, so we have to assume the device is causal, and truncate the decaying sinusoid to only the positive t part. For s with growth rather than decay you can assume a finite length of impulse response for the DUT. In the end what we get is a "convolution" between a kernel and the DUT's frequency domain S parameters, in other words by looking at neighboring frequencies you can obtain the s-domain S parameters at a particular s.
Put all this together and it should be possible to develop a replacement for RFSim99 that can plot in the s-domain.