Estimate transfer function from Bode curve I have measured the magnitude response Y_mag(f) and phase response Y_phase(f) of an unknown physical system.

Is it possible to estimate the s-domain transfer function without a priori information about the order of the polynomials in the numerator or denominator?

I am more or less looking for an out-of-the-box approach that estimates the transfer function in a sufficient way. Performance and computing time are not an issue.

Python and MATLAB (Signal Processing Toolbox, Control System Toolbox) are available.
So this page suggests that you might be able to use this package in Python.

The main issue is that that package seems to require a time series for estimation, rather than the frequency domain data you have.  To get around that, you may be able to use your frequency domain data to generate a time series.

Suppose your frequency domain data is 
$$Y(\omega = 2\pi f) = Y_{\tt mag}(f)\exp(\imath Y_{\tt phase}(f)), \ \ f&gt;0$$ 
and the complex conjugate for $f&lt;0$.

Interpolate your data:

$$
y(t) = {\tt IFFT}(Y(\omega))\\
y_{\tt padded}(t) = \left \{ 
\begin{array}{ll}
y(t) &amp; \mbox{for } 0 \le t \le T_y\\
0 &amp; \mbox{for } T_y \le t \le 2T_y \mbox{  (say)}
\end{array}
\right .\\
Y_{\tt padded}(\omega) = {\tt FFT}(y_{\tt padded}(t))
$$
then just use this to filter a padded white noise sequence of length $T_y$ (padded to $2T_y$).