I usually use baseband equivalent models and VerilogA. You can simulate the interferer by mixing an independent complex baseband signal up to the offset frequency using a complex multiplication with exp(jwt). That signal adds to the baseband equivalent of the desired input signal. The sum then goes to the baseband equivalent model of the downconverting mixer pair, which amounts to another complex multiplication but this time you multily with exp(j*phi) where phi is the phase noise. You simulate the phase noise by passing white Gaussian noise through a digitial filter, perhaps implemented with a z-domain block. You generate the white Gaussian noise with a random number generator updated at a sufficiently high rate.

But before you dive in, ask yourself how you plan to measure the desensing effects of reciprocal mixing. Do you plan to run a BER simulation? An EVM simulation? A simple rms noise simulation with only the interferer present? The answer may confine you to a particular simulator, perhaps one that does not support VerilogA. Perhaps it has something better.

If you are using Spectre and VerilogA, you could run two identical models in paralle, one corrupted by the interferer and one not so corrupted. You can then compute EVM based on the normalized difference between the two sets of outputs.

I have also written VerilogA models to sample and record the corrupted outputs to a file so that I could then process the file with the same software the lab equipment uses to compute EVM. This is a good approach if you are dealing with OFDM signals because many times you need to equalize the signal before computing EVM. The lab software has the equalizer built in.