## 3 Review of Gaussian Process Regression Models
### GP-Rlaxx Models  
Ref. [12] proposed an alternative to the GP-NARX model. Their GP-RLARX model assumes a _latent_ autoregressive structure for the lagged inputs, leading to the description below:  
\[y_{t} =x_{t}+\varepsilon_{t}^{(y)} \tag{19a}\] \[x_{t} =f(x_{t-1},\ldots,x_{t-L_{x}},c_{t-1},\ldots,c_{t-L_{x}})+ \varepsilon_{t}^{(x)}\] (19b) \[f(\cdot) =\mathcal{GP}(\mu(\cdot),k(\cdot\,\cdot))\] (19c) \[\varepsilon_{t}^{(x)} \sim\mathcal{N}(\varepsilon_{t}^{(x)}|0,\sigma_{x}^{2})\] (19d) \[\varepsilon_{t}^{(y)} \sim\mathcal{N}\left(\varepsilon_{t}^{(y)}|0,\tau_{t}\right)\] (19e) \[\tau_{t} \sim\mathcal{IG}(\tau_{t}|\alpha,\beta). \tag{19f}\]  
where \(c_{t-1},\ldots,c_{t-L_{x}}\) are lagged exogenous inputs with maximum lag \(L_{x}\), and \(x_{t-1},\ldots,x_{t-L_{x}}\) are the lagged latent states with maximal lag \(L_{x}\). This framework is reminiscent of a state-space model in which (19a) denotes the observation equation at time \(t\) and (19b) is the corresponding state equation, where \(x_{t}\) is an autoregressive function of the preceding \(L_{x}\) lags of the latent state. To facilitate inference in the GP-RLARX model, Ref. [12] used a sparse variational approximation similar to that described in Section 3.1, where \(\mathbf{u}\in\mathbb{R}^{m}\) are inducing points generated by evaluating the GP prior over _pseudo-inputs_\(\mathbf{Z}=\{\mathbf{z}_{1},\ldots,\mathbf{z}_{m}\},\ \mathbf{z}_{i}\in\mathbb{R}^{L_{x}+L_{x}},i=1, \ldots,m\). It follows that \(p(\mathbf{u})=\mathcal{N}(\mathbf{u}|\mathbf{0},\mathbf{K}_{\mathbf{zz}})\), where \(\mathbf{K}_{\mathbf{zz}}\) denotes the kernel covariance matrix evaluated over the _pseudo-inputs_\(\mathbf{Z}\). Then, the GP-RLARX hierarchical model takes the form  
\[p(\mathbf{u}) =\mathcal{N}(\mathbf{u}|\mathbf{0},\mathbf{K}_{\mathbf{zz}}) \tag{20a}\] \[p(f_{t}|\mathbf{u},\mathbf{x}) =\mathcal{N}(f_{t}|a_{\mathbf{x}})_{t},[\mathbf{\Sigma}_{\mathbf{ xx}}]_{tt})\] (20b) \[p(x_{t}) =\mathcal{N}(x_{t}|\mu_{t},\lambda_{t}),\ \forall\ t\in\{1, \ldots,L_{x}\}\ \text{(initial state)}\] (20c) \[p(x_{t}|f_{t}) =\mathcal{N}\left(x_{t}|f_{t},\sigma_{x}^{2}\right),\ \forall\ t\in\{L_{x}+1, \ldots,T\}\] (20d) \[x_{t} =f(x_{t-1},\ldots,x_{t-L_{x}},c_{t-1},\ldots,c_{t-L_{x}})+ \varepsilon_{t}^{(x)}\]0 (20e) \[x_{t} =f(x_{t-1},\ldots,x_{t-L_{x}},c_{t-1},\ldots,c_{t-L_{x}})+ \varepsilon_{t}^{(x)}\]1  
where \(\mathbf{a}_{\mathbf{x}}=\mathbf{K}_{\mathbf{zz}}\mathbf{K}_{\mathbf{zz}}^{-1}\), \(\mathbf{\Sigma}_{\mathbf{xx}}=\mathbf{K}_{\mathbf{xx}}-\mathbf{K}_{\mathbf{ xz}}\mathbf{K}_{\mathbf{zz}}^{-1}\mathbf{K}_{\mathbf{zz}}\), and \(f_{t}=f(x_{t-1},\ldots,x_{t-L_{x}},w_{t-1},\ldots,\allowbreak w_{t-L_{w}})\), with \(f(\cdot)\sim\mathcal{GP}(0,k(\cdot,\cdot))\). For brevity, we denote \(\tilde{\mathbf{x}}_{t}=(x_{t-1},\ldots,x_{t-L_{x}},w_{t-1},\ldots,\allowbreak w _{t-L_{w}})^{{}^{\prime}}\) in the remainder of the paper. The joint distribution is succinctly expressed as  
\[x_{t} =f(x_{t-1},\ldots,x_{t-L_{x}},c_{t-1},\ldots,c_{t-L_{x}})+ \varepsilon_{t}^{(x)}\]2  
Ref. [12] used a variational inference approach [20] to estimate the latent variables, adopting the variational approximation  
\[x_{t} =f(x_{t-1},\ldots,x_{t-L_{x}},c_{t-1},\ldots,c_{t-L_{x}})+ \varepsilon_{t}^{(x)}\]3where \(q(x_{t})=\mathcal{N}(x_{t}|\mu_{t}^{(x)},\lambda_{t}^{(x)})\), \(q(\tau_{t})=\mathcal{IG}(\tau_{t}|a_{t},b_{t})\), \(p(f_{t}|\mathbf{u},\tilde{\mathbf{x}}_{t})=\mathcal{N}(\mathbf{f}|[\mathbf{a}_ {\mathbf{x}}]_{t},[\mathbf{\Sigma}_{\mathbf{xx}}]_{tt})\), and \(q(\mathbf{u})=\mathcal{N}(\mathbf{u}|\mathbf{m}_{\mathbf{z}},\mathbf{\Sigma}_{ \mathbf{zz}})\). In this framework, \(\mu_{t}^{(x)},\lambda_{t}^{(x)},\mathbf{m}_{\mathbf{z}},\mathbf{\Sigma}_{\mathbf{zz }},a_{t},b_{t}\) are variational parameters that are optimized according to a variational inference strategy, similar to that found in [26]. We refer readers to Section 5.1 of [12] for more details, including the exact expression of the ELBO.