System Identification › Parametric Identification

Parametric methods commit to a function class with finitely many parameters and optimize those parameters against data. The commitment is a bargain: you accept misspecification risk in exchange for statistical efficiency, interpretability, and computational tractability. When the class contains (or nearly contains) the truth, parametric methods extract an extraordinary amount from modest data; when it does not, they converge efficiently to the wrong answer. This section builds the two workhorse approaches: regression for linear systems and the prediction-error method as the general template.

Start with the friendliest case, friendly enough that we can see the entire estimation logic laid bare before complications arrive. Suppose we have full-state data $\{x_t, u_t\}_{t=0}^{T-1}$ and posit the linear model class:

$$x_{t+1} = A x_t + B u_t + \varepsilon_t .$$

where:

• $A \in \mathbb{R}^{n \times n}$: the unknown state-transition matrix, propagating today's state forward under the system's internal dynamics
• $B \in \mathbb{R}^{n \times m}$: the unknown input matrix, converting commanded actions into state change
• $\varepsilon_t \in \mathbb{R}^n$: the process-noise residual at time $t$, the part of $x_{t+1}$ the deterministic part cannot account for

The crucial observation is that for a linear model class, each row of the unknown matrices appears in one state equation independently, so what looks like a dynamical-systems problem decomposes into ordinary multiple regression, one row at a time. Stack the unknowns row-wise into a matrix $\Theta = [A \;\; B]$ of shape $n \times (n+m)$, and define the regressor $\varphi_t = [x_t^\top, u_t^\top]^\top$. Every transition is a linear equation:

$$x_{t+1} = \Theta \varphi_t + \varepsilon_t .$$

Arranging all $T$ transitions into matrices $\Phi \in \mathbb{R}^{(n+m)\times T}$ (regressors as columns) and $X_+ \in \mathbb{R}^{n \times T}$ (successor states as columns), ordinary least squares solves

$$\hat\Theta = X_+ \Phi^\top (\Phi \Phi^\top)^{-1},$$

which is exactly $n$ ordinary linear regressions sharing one regressor matrix, one per state dimension, each using the same $T$ transition observations. Read the middle expression physically: $\Phi \Phi^\top$ is the empirical second moment of the stacked state-input regressors used in this fit. Positive definiteness of this finite-record matrix is exactly the condition for a unique least-squares estimate. Without full rank, the displayed inverse does not exist and some parameter directions remain unidentified, although a pseudoinverse still returns one of many minimizing solutions. Persistent excitation is the related asymptotic richness condition on the input sequence: it is designed to keep the relevant empirical information matrices nonsingular as data accumulates, but it is stronger than full rank in one finite record and is not literally identical to the joint state-input rank condition here.

Under i.i.d. Gaussian innovations, least squares coincides with maximum likelihood and is asymptotically efficient under the usual regularity and excitation conditions. With lagged stochastic states as regressors, finite-sample estimates can retain an $O(1/T)$ dynamic-regression bias even though current innovations are uncorrelated with current regressors. Under non-Gaussian white noise with the same exogeneity conditions, least squares remains consistent and asymptotically unbiased but generally loses Gaussian efficiency. The clean distinction for this chapter is therefore consistency versus persistent misspecification bias, not exact finite-sample unbiasedness.

When only inputs and outputs are available, we posit an autoregressive model, the classical ARX form ("auto-regressive with exogenous input"). For a second-order system with scalar output,

$$y_{t+1} = \alpha_1 y_t + \alpha_2 y_{t-1} + \beta_1 u_t + \beta_2 u_{t-1} + e_t ,$$

where:

• $\alpha_1, \alpha_2$: coefficients on the lagged outputs; together they form the characteristic polynomial $\lambda^2 - \alpha_1 \lambda - \alpha_2$, whose roots are the poles we will shortly compare against the true system eigenvalues
• $\beta_1, \beta_2$: coefficients on the lagged inputs
• $e_t$: the one-step-ahead prediction error. Plain least squares is consistent only when this composite error is orthogonal to the lagged-output and lagged-input regressors; process noise or measurement noise entering those lagged outputs can violate that condition. Whiteness is an additional condition used by standard uncertainty and efficiency calculations.

and again least squares applies with regressor $[y_t, y_{t-1}, u_t, u_{t-1}]$. The philosophical move here is elegant: instead of estimating a hidden state and then the dynamics on it, we promote the lagged output history itself to be the state. The model implicitly treats its own past predictions as state, so ARX is a one-step predictor even when the physical system was written in state-space form. For minimal observable systems these two descriptions contain identical input-output information; they differ only in coordinates and in what happens to noise, which is exactly the caveat coming next. Fitting the ARX denominator coefficients and comparing their roots to the true system eigenvalues reveals something important: the recovery happens in input-output coordinates, a different (but equivalent, for minimal systems) description than the internal state.

Notice what changed, because it quietly invalidates the clean theory. With full-state data, the regression noise is the current process disturbance and the current regressor is built from variables available before that disturbance. With output data, eliminating the hidden state produces a composite equation error containing earlier process disturbances; sensor noise can also enter the lagged outputs. Either source can therefore appear in both today's regressor and today's equation error, biasing plain ARX least squares. In the drone example below, process noise alone is enough: the eliminated equation contains $e_t=w_{t,1}-0.75w_{t-1,1}+0.35w_{t-1,2}$, while $y_t$ contains $w_{t-1,1}$. Such correlation can persist with more data, so the estimator may converge to a shifted value. Correcting that bias leads to the general template.

The prediction-error method (PEM) generalizes linear regression into a principle for essentially any model class. Rather than regressing outputs on raw lagged quantities and hoping the implicit noise bookkeeping works out, PEM asks each candidate model to do its own optimal prediction, noise structure included. Define a one-step-ahead predictor associated with candidate model $\theta$: $\hat y_{t+1|t}(\theta)$, the conditional expectation of the next output given everything observed through time $t$, assuming the world behaves like $\theta$. The definition matters: this conditional expectation uses only data available through time $t$, so when past measurements are noisy, the optimal predictor weighs recent observations less than it would trust clean states. For a linear state-space model with Gaussian noises, computing this predictor is exactly the Kalman filter recursion from our previous chapter, run with $\theta$'s matrices. Here the two chapters meet operationally: the filter, previously the consumer of a known model, becomes an internal subroutine of the estimator searching for the model. PEM then optimizes:

$$\hat\theta = \arg\min_\theta \frac{1}{T}\sum_{t=0}^{T-1} \ell\!\big(y_{t+1} - \hat y_{t+1|t}(\theta)\big),$$

typically with quadratic loss $\ell(e) = \|e\|^2$, whereupon PEM reduces to negative log-likelihood maximization for Gaussian output noise. Three properties justify why PEM became the classical centerpiece:

• Statistical grounding. With structural containment, identifiability, informative data, and standard stability and regularity assumptions, a correctly specified maximum-likelihood PEM is consistent and asymptotically efficient. Under misspecification it instead converges to the best approximation selected by its loss, which need not be the physical system.
• Noise-model honesty. Unlike bare ARX regression, PEM lets you declare whether disturbances corrupt the state dynamics (equation error), only the sensor (output error), or something in between. A wrong disturbance model can shift the pseudo-true parameters and forfeits the clean maximum-likelihood efficiency guarantee, even when some parameter estimates remain consistent under weaker conditions. When a network minimizes squared one-step error against noisy observations, it adopts a Gaussian one-step-prediction scoring rule on noisy targets; that loss alone does not specify a classical output-error model structure.
• Model-quality framing. PEM scores exactly the quantity a filter consumes: predictive likelihood of new data. Models good for filtering are judged on their own terms. This alignment between the training objective and the downstream consumer is rarer than it sounds; later we will see planning consumes models quite differently, and the mismatch will matter.

Its costs: PEM optimization is nonlinear in $\theta$ (each loss evaluation runs a filter pass), initialization-sensitive, and locally convergent at best. Each gradient or search step requires filtering the entire record, and the loss surface over multi-dimensional parameter spaces can harbor local minima corresponding to qualitatively wrong dynamics that nonetheless predict the data mediocre-but-decently. Multi-start strategies, careful initialization from cheaper methods, and structure-exploiting implementations mitigate but do not eliminate this friction. This computational friction motivated the family we visit next.

Consider a crude lateral dynamics model of a quadrotor hovering, linearized about hover and discretized at 50 Hz:

$$x_{t+1} = \begin{bmatrix} 0.90 & 0.35 \\ -0.15 & 0.75 \end{bmatrix} x_t + \begin{bmatrix} 0 \\ 0.6 \end{bmatrix} u_t + w_t, \qquad y_t = x_{t,1},$$

where $x_1$ is the lateral velocity proxy, $x_2$ is its rate, $u$ is the differential motor command, and $y_t = x_{t,1}$ records that only the first component of the state reaches the sensor. The complex-conjugate eigenvalues of this $A$ have modulus about $0.853$, describing a moderately damped side-to-side oscillation: errors rotate as they decay rather than disappearing in one step. We simulate it under a spectrally rich shaped-noise input, then fit (a) the full-state linear regression and (b) a second-order ARX model on the scalar output, checking that both recover the dynamics, one up to nothing (states observed, exact coordinates) and one up to input-output equivalence. We then evaluate the fitted full-state model recursively on overlapping starts from the held-out tail, under the recorded exogenous action sequence. Comparing its horizon-30 error with the true-model floor isolates the excess caused by estimation imperfection; the code reports that gap explicitly rather than asking you to infer it from the plot.

Let's code this. Everything numerical in this chapter uses NumPy first, PyTorch only in the neural section.