My specific application is analysis of dynamic textures using linear dynamics systems of the form
$$ I(t) = Cz(t) + w(t) \\ z(t + 1) = Az(t) + Bv(t), $$
where $I(t)$ is the original signal, $z(t)$ is its low-dimensional representation, and $C$ is the orthogonal basis for the subspace. For comparing two systems, the discrete Lyapunov equation can be used, taking the form
$$ \mathcal{A} = \begin{bmatrix} A_1 & \textbf{0} \\ \textbf{0} & A_2 \end{bmatrix} $$ $$ \mathcal{Q} = \big[C_1 C_2\big]^T\big[C_1 C_2\big] $$ $$ \mathcal{A}\mathcal{X}\mathcal{A}^T - \mathcal{X} + \mathcal{Q} = 0. $$
However, this only works for first-order systems; that is, systems which for $z(t + 1)$, only the previous time point $z(t)$ is required. How can I generalize this formulation to work for n-th order systems, e.g. $z(t + 1) \approx A_1z(t) + A_2z(t - 1) + ... + A_dz(t - d + 1)$?
Just to clarify: in the definition of $\mathcal{A}$, the terms $A_1$ and $A_2$ are meant to refer to $A$ from system 1 and system 2, respectively. In the approximation of n-th order $z(t + 1)$, the terms $A_1$ through $A_d$ all belong to one system; thus, comparing multiple systems would somehow involve comparing the $d$ parameters $A$ that each separately defines. How exactly that could be done is the question I'm asking here.
You can augment your states as the following:
$$\begin{bmatrix} z(t+1) \\ z(t) \\ z(t-1) \\ \vdots \\ z(t-d+2) \end{bmatrix} = \begin{bmatrix} A_1 & A_2 & \dots & A_{d-1} & A_d \\ I & 0 & \dots & 0 & 0 \\ 0 & I & \dots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & I & 0 \end{bmatrix} \begin{bmatrix} z(t) \\ z(t-1) \\ z(t-2) \\ \vdots \\ z(t-d+1) \end{bmatrix}$$
Then, you get $\tilde{z}(t+1) = \tilde{A} \tilde{z}(t)$.