How to validate an estimated model in case of output-only data (in frequency domain)? Before moving to the actual question, I would like to emphasize on the following points (maybe they are obvious to some of you, but I still would like to list them, since they make the difference): 


I don't have any reliable information on the input, thus I have only
the measurements of the output $y(t)$.
I have only 1 set of recorded data, thus no repetitive measurements are in this case available.


Objective: Identifying a reliable transfer function model for this data in the frequency domain (e.g. rational fraction polynomial with numerator $B$ and denominator $A$)

Question: After obtaining a transfer function model, how does one validate this model and thus can be assured that it does not contain any systematic/modelling errors? Please, keep in mind the above 2 points!  

To my limited experience, calculating the uncertainty of the transfer function residuals will be not possible, since the missing information on the input (and thus also not having the variance of the input). A second approach that I thought of is to do a 'simple' whiteness test of the residuals (sum of squares due to errors (SSE) or the R-square both mentioned in the Matlab manual), but is this approach even significant? I mean, sure... One obtains a result, but it's complex valued due to the data ($Y(k)$) which is also complex valued. Consequently, there is always an imaginary part which makes it fail the criterion (i.e. model is ok if SSE is close to 0). Any thoughts/solutions on this? Or am I mixing things up?

Any help, advise/guidance (tips) or references that I can consult in order to solve this issue is very much appreciated! Because I'm really lost and totally confused.
I don't know of a standard way to do this, but I think taking a leaf out of the training of artificial intelligence algorithms might help:  


divide your output data into two disjoint sets, 
perform whatever identification procedure you want to obtain the ARMA parameters, and 
then compare the two transfer functions using whatever distance metric makes sense for your application.


In the figure below, I've decided that my distance metric is that the poles and zeros of each ARMA model are "not too far" from each other (or the ARMA poles and zeros found using all the data).


The $\color{black}{\bf black} $ markers are the pole ($\times$) and zero ($\circ$) positions for the ARMA model calculated from first half of the data ($G_1$).
The $\color{red}{\bf red} $ markers are the pole ($\times$) and zero ($\circ$) positions for the ARMA model calculated from second half of the data ($G_2$).
The $\color{grey}{\bf grey} $ markers are the pole ($\times$) and zero ($\circ$) positions for the ARMA model calculated from all of the data ($G_{\rm tot}$).




Another way to do it might be to look at the transfer function
$$
E = G_1 - G_2
$$
That is plotted below using the freqz command in R.





R Code Below

#30332
library(signal)
library(tseries)

# Create a time series that is not just white noise
T &lt;- 1000

do_plot &lt;- 1

#for (k in seq(1,100))
{
  x &lt;- rnorm(T, 0, 1)
  bf &lt;- butter(4,0.41)
  y &lt;- filter(bf$b, bf$a, x)

  # Split data in two and compare ARIMA models
  L2 &lt;- length(y)/2
  y1 &lt;- y[100:L2] # Start away from 1 to avoid startup transients
  y2 &lt;- y[(L2+1):length(y)]

  orders &lt;- c(3,0,3)

  # First half
  tryCatch(
  {arima1 &lt;- arima(y1, order = orders)},
  warning = function(w) {  },
  error = function(e) {  },
  finally = {});

  all_1 &lt;- coefficients(arima1)
  num_1 &lt;- all_1[4:6]
  den_1 &lt;- all_1[1:3]

  # Second half
  tryCatch(
  {arima2 &lt;- arima(y2, order = orders)},
  warning = function(w) {  },
  error = function(e) {  },
  finally = {});

  all_2 &lt;- coefficients(arima2)
  num_2 &lt;- all_2[4:6]
  den_2 &lt;- all_2[1:3]

  # All data
  arima_tot &lt;- arima(y, order = orders)
  all_tot &lt;- coefficients(arima_tot)
  num_tot &lt;- all_tot[4:6]
  den_tot &lt;- all_tot[1:3]

  # Plot  boundaries
  xs &lt;- c(-1,1)
  ys &lt;- xs

  if (do_plot == 1)
  {
    par(mfrow=c(1,1), pty='s')
    plot(roots(den_1), xlim=xs, ylim=xs, lwd = 6, pch=4, xlab="Re(z)", ylab="Im(z)")
    do_plot &lt;- 0
  }
  else
  {
    points(roots(den_1), lwd = 6, pch=4)
  }
  points(roots(den_2), lwd = 4, col="red", pch=4)
  points(roots(den_tot), lwd = 2, col="grey", pch=4)

  points(roots(num_1), lwd = 6, pch=1)
  points(roots(num_2), lwd = 4, col="red", pch=1)
  points(roots(num_tot), lwd = 2, col="grey", pch=1)
}


par(pty="m")
subtract_num &lt;- conv(num_1, den_2) - conv(num_2, den_1)
subtract_den &lt;- conv(den_1, den_2)

freqz(subtract_num, subtract_den)