Your theoretical transfer function is G(z) = 1/(1 + z^-2 + 0.5*z^-4). The frequency response is the value of this function for z = exp(j*w*Ts). Ts=2000. w is the frequency for which you want to calculate the spectrum (say, linspace(0, pi/Ts, 1000))
To compute the spectra of time series models th_fb etc, use freqresp or bode command. However, note that AR expects to model a stationary time series (or an impulse response in deterministic case). So I would apply AR to diff(diff(Y)) to get the various models; you might need to set the noisevariance to 1 after estimation for comparison purposes.
Result using FFT: |G(w)| = fft(Y)./|fft(x)| This will not not good luck for a variety of reasons, in addition to the non-periodic nature of your input signal.
