Transfer Function Identification giving negative static value

I'm trying to fit a complex system with the identification toolbox of matlab. The fit is great for the frequency and phase if I go to an order around 20, with the Transfer Function estimator, but the step response goes on a negative value. This is impossible on my physical system.

If I use other techniques, like state space, the fit is not good enough in the low frequencies, even for higher orders, and I cannot use the "focus" function (to implement a custom filter). I forced my identification to be stable, and so, all the poles are inside the unit circle. My system is a non-minimum phase system, and so, some zeros are outside the unit circle.

What is wrong in my procedure ? Is the step response reliable or 20 poles and zeros is too much ? How can a perfect fit of amplitude and phase in the frequency domain gives a wrong step response ?