## 2 Problem Description and Identification Model  
For narrative convenience, let us introduce some notation. The nomenclature is displayed in Table 1.  
Consider the following state-space system with colored noise,  
\[\mathbf{x}(k+1) =\mathbf{G}\mathbf{x}(k)+\mathbf{h}u(k), \tag{1}\] \[y(k) =f\mathbf{x}(k)+e(k), \tag{2}\] \begin{table}
\begin{tabular}{c l l l} \hline \hline \multicolumn{2}{c}{**Nomenclature**} & & \\ \hline T & matrix/vector transpose & \(\mathbf{A}^{-1}\) & inverse of the matrix \(\mathbf{A}\) \\ E & expectation operator & \(z\) & unit forward shift operator \\ \(\|\mathbf{X}\|\) & 2-norm of a matrix \(\mathbf{X}\) & \(\mathbf{I}_{m}\) & \(m\times m\) identity matrix \\ \(\mathbf{1}_{n}\) & \(n\)-dimensional column vector whose elements are 1 & \\ \(\hat{\mathbf{\theta}}(k)\) & estimate of the parameter vector \(\mathbf{\theta}\) at instant \(k\) & \\ \(\mathbf{\theta}(1:m)\) & a vector consisting of the first entry to the \(m\)-th entry of \(\mathbf{\theta}\) & \\ \hline \hline \end{tabular}
\end{table} 
Table 1: The nomenclature for symbols. where \(\mathbf{x}(k):=[x_{1}(k),x_{2}(k),\cdots,x_{n}(k)]^{\intercal}\in\mathbb{R}^{n}\) is the state vector, \(u(k)\in\mathbb{R}\) is the system input, \(y(k)\in\mathbb{R}\) is the system output, \(e(k)\in\mathbb{R}\) is a random noise and \(\mathbf{G}\in\mathbb{R}^{n\times n}\), \(\mathbf{h}\in\mathbb{R}^{n}\) and \(\mathbf{f}\in\mathbb{R}^{1\times n}\) are the system parameter matrix and vectors, defined as:  
\[\mathbf{G}:= \left[ \begin{array}{ccccc}0&1&0&\cdots&0\\ 0&0&1&\ddots&\vdots\\ \vdots&\vdots&\ddots&\ddots&0\\ 0&0&\cdots&0&1\\ g_{1}&g_{2}&g_{3}&\cdots&g_{n}\end{array} \right]\in\mathbb{R}^{n\times n}, \quad\mathbf{h}:=\left[ \begin{array}{c}h_{1}\\ h_{2}\\ \vdots\\ h_{n-1}\\ h_{n}\end{array} \right]\in\mathbb{R}^{n},\] \[\mathbf{f}:= [1,\ 0,\cdots,\ 0,\ 0]\in\mathbb{R}^{1\times n}.\]  
The external disturbance \(e(k)\) can be fitted by a moving average process, an autoregressive process or an autoregressive moving average process. Without loss of generality, we consider \(e(k)\) as a moving average noise process:  
\[e(k):=(1+e_{1}z^{-1}+e_{2}z^{-2}+\cdots+e_{n_{e}}z^{-n_{e}})v(k), \tag{3}\]  
where \(v(k)\) is the white noise with zero mean and variance \(\delta\). Assume that \(y(k)=0\), \(u(k)=0\) and \(v(k)=0\) for \(k\leqslant 0\). The objective is to estimate the parameters \(g_{i}\), \(h_{i}\) and \(e_{i}\) and the system state \(\mathbf{x}(k)\) from the available input-output data \(\{u(k),y(k):\ k=0,1,2,\cdots\}\). Note that the system in Equations (1) and (2) is an observability canonical form, and the observability matrix \(\mathbf{T}\) is an identity matrix, i.e.,  
\[\mathbf{T}=[\mathbf{f}^{\intercal},(\mathbf{f}\mathbf{G})^{\intercal},\cdots,(\mathbf{f}\mathbf{G}^{ n-1})^{\intercal}]^{\intercal}=\mathbf{I}_{n}. \tag{4}\]  
From Equations (1) and (2), we have:  
\[y(k) = \mathbf{f}\mathbf{x}(k)+e(k), \tag{5}\] \[y(k+1) = \mathbf{f}\mathbf{x}(k+1)+e(k+1)\] (6) \[= \mathbf{f}\mathbf{G}\mathbf{x}(k)+\mathbf{f}hu(k)+e(k+1),\] \[y(k+2) = \mathbf{f}\mathbf{G}\mathbf{x}(k+1)+\mathbf{f}hu(k+1)+e(k+2)\] (7) \[= \mathbf{f}\mathbf{G}^{2}\mathbf{x}(k)+\mathbf{f}\mathbf{G}hu(k)+\mathbf{f}hu(k+1)+e(k+2),\] \[\vdots\] \[y(k+n-1) = \mathbf{f}\mathbf{G}^{n-1}\mathbf{x}(k)+\mathbf{f}\mathbf{G}^{n-2}\mathbf{h}u(k)+\mathbf{f} \mathbf{G}^{n-3}\mathbf{h}u(k+1)+\cdots\] (8) \[e(k):=(1+e_{1}z^{-1}+e_{2}z^{-2}+\cdots+e_{n_{e}}z^{-n_{e}})v(k), \tag{3}\]0 \[e(k):=(1+e_{1}z^{-1}+e_{2}z^{-2}+\cdots+e_{n_{e}}z^{-n_{e}})v(k), \tag{3}\]1 (9) \[e(k):=(1+e_{1}z^{-1}+e_{2}z^{-2}+\cdots+e_{n_{e}}z^{-n_{e}})v(k), \tag{3}\]2  
Define the parameter matrix \(\mathbf{M}\) and the information vectors \(\mathbf{\phi}_{\text{y}}(k)\), \(\mathbf{\phi}_{\text{u}}(k)\) and \(\mathbf{\phi}_{\text{e}}(k)\) as:  
\[\mathbf{M}:= \left[