Prediction-Error Models, Order Selection & Residual Validation

A black-box model contains two stories: how the measured input drives the output, and how unmeasured disturbances make the output unpredictable. ARX, ARMAX, output-error and Box–Jenkins structures differ mainly in which poles and zeros those two stories are forced to share. Choosing among them is therefore a claim about disturbance physics, not a choice of file format.

Plant, noise and prediction #

Write a discrete-time linear model as

$$y_k=G(q,\theta)u_k+H(q,\theta)e_k,$$

where $q^{-1}u_k=u_{k-1}$, $G$ is the deterministic input/output dynamics, $H$ is a stable monic disturbance filter, and $e_k$ is the innovation: the part of $y_k$ that cannot be predicted from past measured data under the model.

The one-step prediction error is

$$\varepsilon_k(\theta)=H^{-1}(q,\theta)\left[y_k-G(q,\theta)u_k\right],$$

and the one-step predictor is

$$\hat y_{k|k-1}=y_k-\varepsilon_k =\left(1-H^{-1}\right)y_k+H^{-1}G u_k.$$

The predictor uses past measured output whenever $H\ne1$. Changing the noise model therefore changes the criterion used to fit $G$; it does not merely add a decorative disturbance output after plant estimation.

Prediction-error estimation chooses

$$\hat\theta_N=\arg\min_\theta V_N(\theta),\qquad V_N=\frac1N\sum_{k=1}^{N}\varepsilon_k^T\Lambda^{-1}\varepsilon_k,$$

with output weighting or innovation covariance $\Lambda$. Under a correctly specified innovation distribution and likelihood loss, this is maximum likelihood. Otherwise it is still a defined prediction-error criterion, but not automatically a probabilistic optimum.

One polynomial family, four common structures #

A general SISO polynomial form is

$$A(q)y_k=\frac{B(q)}{F(q)}u_{k-n_k}+\frac{C(q)}{D(q)}e_k,$$

with delay $n_k$. Different choices set some polynomials to one:

The polynomial orders $n_a,n_b,n_c,n_d,n_f$ count free lag coefficients; $n_k$ is a pure input delay in samples. Their numerical values are not comparable across structures without the equations: a “second-order ARX” and a “second-order OE” make different noise assumptions.

ARX #

Expanding

$$A(q)y_k=B(q)u_{k-n_k}+e_k$$

expresses $y_k$ as a linear combination of past $y$ and $u$. The parameter estimate is ordinary or regularized least squares. This makes ARX fast for scanning delays and orders, and useful when its disturbance structure is plausible or signal-to-noise ratio is high.

Its weakness is the shared $A$: if plant and disturbance have different poles, ARX may increase order to approximate both. Because past noisy outputs are regressors, arbitrary output noise can also bias the plant estimate when the assumed $1/A$ noise path is wrong. A low loss does not repair the structural coupling.

ARMAX #

The moving-average polynomial $C$ lets the innovation shape output disturbance beyond $1/A$. This is useful when disturbance enters through dynamics that plausibly share the plant denominator, such as a load entering near the plant input. Estimation is no longer a simple linear regression because past innovations depend on the unknown parameters; initialization and nonlinear optimization matter.

Output-error #

OE separates plant poles into $F$ and assumes $H=1$. Its prediction error equals the output simulation error for the deterministic plant model. When input is independent of output noise and the plant class is appropriate, OE can target $G$ without forcing plant poles to explain colored output disturbance.

But it does not make colored disturbance disappear. Residuals can remain strongly autocorrelated, one-step prediction can be poor, and uncertainty estimates relying on white innovations become misleading. In closed loop, ignoring the disturbance model is especially risky because $u$ and disturbance are correlated through feedback.

Box–Jenkins #

BJ gives $B/F$ and $C/D$ independent poles and zeros. It can represent a plant and colored measurement/load disturbance with genuinely different dynamics. The price is a larger, more correlated parameter vector and a nonconvex numerical problem. A flexible BJ model needs more informative data, careful initialization and stronger validation; otherwise it can trade plant dynamics against noise dynamics without a unique physical interpretation.

State-space innovation models #

A multivariable alternative is

$$x_{k+1}=Ax_k+Bu_k+Ke_k,$$

$$y_k=Cx_k+Du_k+e_k.$$

Here $K$ is the innovation gain that shapes how unpredictable output enters state prediction. It is not the controller feedback gain. This form provides a compact common state for MIMO dynamics and noise and avoids high-order polynomial matrices.

Subspace methods first recover a state sequence or observability subspace from block Hankel matrices, then estimate $A,B,C,D$; PEM can refine the result. Singular-value gaps are useful order clues, not automatic truth: colored noise, feedback, short records and weakly excited modes blur them. Similarity transformations mean the estimated state coordinates are not unique even when the input/output model is.

What the optimizer is really fitting #

When the true system lies outside the candidate set, PEM converges toward the model minimizing an experiment-weighted prediction error. In open loop, the input spectrum emphasizes frequencies with more excitation, and $H^{-1}$ reweights the mismatch. A model can fit the experiment excellently while being poor at an unexcited control crossover.

Consequences:

• record and plot the input spectrum beside every frequency-domain fit claim;
• choose weighting from the intended decision, but report it because it changes the estimand;
• use several initial parameterizations for nonlinear structures;
• constrain stability or parameter ranges only when the use demands it and disclose the constraint;
• inspect whether different structures converge to compatible $G$ within uncertainty.

An optimization message saying “converged” means the numerical stopping conditions were met. It does not prove a global minimum, identifiability, correct order or model adequacy.

Delay and order selection #

Estimate gross delay first from physical timing, cross-correlation, impulse response or phase slope, then scan a plausible discrete delay range. An incorrect $n_k$ forces numerator/denominator coefficients to imitate dead time. Too-small delay leaves input-correlated residuals at missing lags; too-large delay throws away real early response.

For a SISO model with prediction-error variance

$$V=\frac1N\sum_{k=1}^{N}\varepsilon_k^2$$

and $p$ free parameters, common criteria, up to constants shared across candidates, are

$$\mathrm{AIC}=N\log V+2p,$$

$$\mathrm{BIC}=N\log V+p\log N,$$

$$\mathrm{AIC_c}=\mathrm{AIC}+\frac{2p(p+1)}{N-p-1},$$

and Akaike's final prediction error is approximately

$$\mathrm{FPE}=V\frac{1+p/N}{1-p/N}.$$

AIC emphasizes expected predictive performance; BIC applies a stronger large-sample complexity penalty; AICc corrects AIC when $p$ is not negligible relative to $N$. These formulas assume comparable data, likelihood/error definitions and parameter counting. Estimated initial states, multiple outputs, weighting and effective sample loss from lags affect the exact implementation.

Information criteria rank the candidate list. If every candidate omits a delay, nonlinear regime or disturbance channel, the winner is merely the least bad member of a wrong set. Use the smallest order whose independent validation, residuals and decision-relevant frequency response are adequate; do not select order from one scalar alone.

Simulation and prediction answer different questions #

Free-run simulation generates output from inputs and initial conditions without repeatedly correcting from measured output. It tests whether $G$ reproduces accumulated plant dynamics.

One-step prediction uses past measured output through $H$ and can track a noisy record well even when the free-running plant model is poor. Longer-horizon prediction moves between these extremes.

For controller design, validate free-run response and frequency response of $G$ over the intended band. For forecasting or monitoring, validate the actual prediction horizons. Reporting one-step fit as “model accuracy” without the horizon and response mode is incomplete.

A common normalized fit display is

$$100\left(1-\frac{\lVert y-\hat y\rVert_2}{\lVert y-\bar y\rVert_2}\right)\%.$$

It can be negative, depends on the tested record and baseline, and hides phase-local or frequency-local failures. Preserve the raw error units and plots.

Residual tests #

For an adequate one-step predictor, residuals should resemble the innovation class assumed by the model.

Whiteness #

For nonzero lag $\ell$,

$$R_{\varepsilon\varepsilon}[\ell]=E[\varepsilon_k\varepsilon_{k-\ell}]\approx0.$$

Significant autocorrelation means predictable temporal structure remains: missing poles/zeros, wrong delay, colored disturbance, nonlinearity, drift or preprocessing artifacts. A whiteness test is relative to record length and confidence level; isolated crossings among many lags are expected.

Independence from input #

In an appropriate open-loop experiment,

$$R_{\varepsilon u}[\ell]\approx0$$

at causally relevant lags. Correlation means measured input still explains residual output—evidence of missing dynamics, wrong delay or bias. In closed loop, $u$ is inherently correlated with output disturbance, so the expected correlation structure and valid test depend on the identification method and external excitation. Do not apply the open-loop interpretation mechanically.

Distribution and operating dependence #

Plot residuals against predicted output, input amplitude, state, time and operating conditions. Heteroscedasticity suggests amplitude-dependent noise or nonlinearity; clusters by mode suggest an omitted switched regime; heavy tails may reflect outliers, jumps or saturation. Gaussianity is required for some interval and likelihood interpretations, but not for the basic demand that predictable structure be absent.

Residual whiteness is necessary for a declared innovation model, not sufficient for physical truth. A high-order black box can whiten one experiment while placing fragile poles and zeros outside the band where data carry information.

Worked structure diagnosis #

Suppose the plant and output disturbance are approximately

$$G_0(q)=\frac{0.3q^{-1}}{(1-0.8q^{-1})(1-0.2q^{-1})},$$

$$H_0(q)=\frac{1+0.4q^{-1}}{1-0.6q^{-1}}.$$

The plant poles are $0.8$ and $0.2$; the disturbance pole is $0.6$.

• A low-order ARX must use one $A$ polynomial for both, so it cannot exactly represent all three distinct poles. Increasing $n_a$ may improve prediction while adding a pole that belongs to disturbance rather than plant.
• ARMAX adds $C$ but still shares $A$ poles between $G$ and $H$; it remains structurally mismatched unless a higher shared denominator approximates both paths.
• OE can represent $G_0$ with $F=(1-0.8q^{-1})(1-0.2q^{-1})$, but leaves the colored $H_0$ dynamics in residuals.
• BJ can represent $G_0$ with $B/F$ and $H_0$ with $C/D$ independently.

BJ is the structurally correct candidate here, not automatically the best finite-data estimate. With weak noise and short data, OE may estimate the plant more precisely because it spends fewer parameters. The choice is resolved by independent plant-response validation, residual diagnostics and uncertainty—not by naming the most flexible family.

Parameter and response uncertainty #

An estimated covariance matrix describes local uncertainty around $\hat\theta$ under the chosen structure, data and noise assumptions. Strong parameter correlations and a nearly singular information matrix reveal combinations the experiment cannot separate.

Transform parameter uncertainty into the quantities used for decisions:

• frequency-response bands near controller crossover;
• pole/damping and delay uncertainty;
• step/impulse response bands;
• simulation envelopes over initial conditions;
• uncertainty of steady gain and disturbance paths.

A small parameter standard deviation inside a wrong model class is false precision. Compare models across independent runs and operating points; between-run variation often dominates asymptotic covariance from one stationary record.

Validation protocol #

1. Reproduce the estimation criterion and verify numerical convergence from multiple starts when needed.
2. Compare nonparametric and parametric impulse/frequency responses over the excited band with uncertainty.
3. Simulate on independent runs with the initial-state treatment declared.
4. Evaluate one-step and decision-relevant multi-step prediction separately.
5. Test residual autocorrelation and input correlation with multiple-lag interpretation.
6. Check poles, zeros, delays, units and signs against physical constraints without forcing convenient cancellations.
7. Repeat across amplitudes and operating points to define the LTI validity envelope.
8. For control use, close the loop in simulation with uncertainty, sensor/actuator dynamics, saturation and representative disturbances.
9. Publish rejected structures and why they failed; otherwise a single winning fit hides model-selection uncertainty.

The finished artifact is not a coefficient vector. It is a model structure, estimated parameters, disturbance assumptions, uncertainty, validation evidence, data provenance and a clear statement of where the model is not trusted.