Can it be so hard to estimate a simple transfer function?

I have bode plot data from a circuit simulator (Cadence SpectreRF). To first order, the system should be an integrator but it has much high frequency stuff going on. According to the simulator's .pz analysis (estimation), it has 17 poles and 14 zeros.

I use 100 points per decade in my data (which is a lot).

Now no matter which number of poles/zeros I take, this ftest function just gives garbage:

gfr = idfrd(vod.V, 2*pi*vod.freq, 0);
save gfr;
figure, bode(gfr, tfest(gfr,3,2));

(blue: original data, red: estimate; see code).

Am I doing something wrong? Can it be so hard to estimate a simple transfer function?

For reference: Here is the gfr object data with the data: https://www.dropbox.com/s/zn3pkz9ckihjdbl/gfr.mat?dl=0