Estimate transfer function from Bode curve
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$).