How to validate an estimated model in case of output-only data (in frequency domain)?
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)