I found out a possible reason for the DC gain by myself:
The system identification toolbox gives an estimation of the noise variance. So although I generate my data using randn with variance 1, the toolbox assumes (estimates) another noise variance, which I should use to scale my estimated transfer function to get the right DC gain.
So if in the above code, the armax model estimation part, I use
idout_armax = armax(td_data, [2, 0, 1, 1]);
scale = sqrt(idout_armax.NoiseVariance);
idsys_armax = tf(idout_armax.c, idout_armax.a, Ts)*scale;
figure, bode(sysdt, idsys_armax)
the DC gain should (almost) match sometimes.
I hope this is the correct reason.