How to estimate a model with a positive impulse response, i.e., with only real poles?

I would like to estimate a system model from input and output data that has a purely positive impulse response. If I understand correctly, this is equivalent to all the poles of the transfer function being real.

One possible solution to the problem is given by the "process estimation" methods, which by default enforce real poles. For instance,

load dryer2;
dry = iddata(y2,u2,0.08);
opt = procestOptions('Focus', 'simulation');
sys = procest(dry, 'P3DZ', opt);
figure;
impulse(sys);

does exactly what I want (see below).

However, the procest method seems to only work for a very limited class of models (max. 1 zero, max. 3 poles). Is there any way to achieve the same thing for more general models?