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.

• Welcome to DSP.SE! This is an interesting problem. A comment: not having the variance of the input will mean you can't know the gain, so you'll have to assume some value for the input. Question: What are you actually going to do with the model once you have it? That may have a bearing on how to approach the problem. – Peter K. Apr 25 '16 at 14:00
• @Peter - Thank you for the warm welcome. Regarding your question: once I've the validated model, then I would like to use the poles (and the corresponding zeros) that stems from that model in a subsequent step. However, it doesn't involve any feedback-loop, so it's all within the context of an 'open loop' system where the outcome (i.e poles of the transfer function) is deduced/taken and substituted in another separate step, which will assume that the information is correct. Therefore, the need of some kind of tangible form of model validation, in order to substantiate this made assumption. – WaffleTeX Apr 25 '16 at 15:29
• Just another question: you mention y(t) and Y(k) are they the same? And your data, y(t) is complex? – Peter K. Apr 26 '16 at 19:34

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 <- 1000

do_plot <- 1

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

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

orders <- c(3,0,3)

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

all_1 <- coefficients(arima1)
num_1 <- all_1[4:6]
den_1 <- all_1[1:3]

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

all_2 <- coefficients(arima2)
num_2 <- all_2[4:6]
den_2 <- all_2[1:3]

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

# Plot  boundaries
xs <- c(-1,1)
ys <- 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 <- 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 <- conv(num_1, den_2) - conv(num_2, den_1)
subtract_den <- conv(den_1, den_2)

freqz(subtract_num, subtract_den)