Transfer function determination from input and output data
if your input excitation signal (call it $x[n]$) is arbitrarily long and so also your output $y[n]$, then you should be able to take as long of block (let's call that length $N$) of each and, using the DFT and division compute
$$ H[k] = \frac{Y[k]}{X[k]} $$

in audio, we call that method of system identification the "dual-channel FFT".  the DFT of $h[n]$ is equal to the the DTFT of $h[n]$ (wrapped-around and overlapped and added every $N$ samples) which is the Z-transform of the same wrapped (and aliased) $h[n]$ sampled at $N$ equally-spaced points around the unit circle.
$$ H(e^{j 2 \pi k/N}) = H[k]  $$

so if your 4th-order IIR transfer function is
$$ G(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + b_3 z^{-3} + b_4 z^{-4}}{1 + a_1 z^{-1} + a_2 z^{-2} + a_3 z^{-3} + a_4 z^{-4}} $$

what you want to do is choose $ \{ b_0, b_1, b_2, b_3, b_4, a_1, a_2, a_3, a_4 \} $ in such a way to minimize
$$ \sum_{k=0}^{N-1} \left| G(e^{j 2 \pi k/N}) - H(e^{j 2 \pi k/N}) \right|^2 \ = \  \sum_{k=0}^{N-1} \left| G(e^{j 2 \pi k/N}) - H[k] \right|^2 $$
or, perhaps a weighted version of that, 
$$ \sum_{k=0}^{N-1} W[k] \left| G(e^{j 2 \pi k/N}) - H[k] \right|^2 $$
where $W[k]$ is known.