## 3 Materials and Methods
### Methodology
#### 3.1.3 Recursive Estimate  
Given the observed data \(\mathbf{Y}_{t}\) and the latent variable \(\widehat{\mathbf{c}}_{t-1}\) at time \(t\), we estimate \(Z_{t}(x)\) to obtain \(\widehat{Y}_{t}(x)\). In this case, we use BLUE, which is widely used in universal Kriging, with expression by Equation (11):  
\[\widehat{Z}_{t}(x)=\widehat{c}_{t-1}(x)+\mathbf{q}_{t}^{T}(x)\widetilde{\mathbf{Y}}_{t}. \tag{13}\]  
Unbiasedness implies that the estimate must satisfy \(E[\widehat{Z}_{t}(x)]=E[\widehat{c}_{t-1}(x)+\mathbf{q}_{t}^{T}(x)\widetilde{\mathbf{Y }}_{t}]=\mathbf{f}^{T}(x)\mathbf{m}_{t}\), which leads to the constraint  
\[\mathbf{F}\mathbf{W}_{t}^{T}\mathbf{q}_{t}(x)=\mathbf{f}(x). \tag{14}\]  
The best estimate \(\widehat{Z}_{t}(x)\), then minimizes the mean squared error \(E[(Z_{t}(x)-\widehat{Z}_{t}(x))^{2}]\), subject to the constraint  
\[ \begin{array}{c}\sigma_{cz}^{2}=E[(Z_{t}(x)-\widehat{c}_{t-1}(x)-\mathbf{q}_{t} ^{T}(x)\widetilde{\mathbf{Y}}_{t})^{2}]\\ =E[(\varepsilon_{t}(x)+c_{t-1}(x)+\mathbf{f}^{T}(x)\mathbf{m}_{t}-\widehat{c}_{t-1}(x) -...\\ \mathbf{q}_{t}^{T}(x)(\mathbf{W}_{t}\mathbf{e}_{t}+\mathbf{W}_{t}\mathbf{e}_{t-1}-\mathbf{W}_{t} \widehat{c}_{t-1}+\mathbf{W}_{t}\mathbf{F}^{T}\mathbf{m}_{t}+\mathbf{v}_{t})^{2}].\end{array}  \tag{15}\]  
Squaring the argument and applying the expectation yields  
\[ \begin{array}{c}\sigma_{cz}^{2}=\sigma_{cz}^{2}(x)+\sigma_{c}^{2}(x)+\mathbf{q} _{t}^{T}(x)(\mathbf{W}_{t}\mathbf{P}_{t-1}\mathbf{W}_{t}^{T}+\mathbf{W}_{t}\mathbf{K}\mathbf{W}_{t}^{T} +\mathbf{R}_{t})\mathbf{q}_{t}(x)-...\\ \mathbf{2}\mathbf{q}_{t}^{T}(x)\mathbf{W}_{t}(\mathbf{p}_{t-1}(x)+\mathbf{k}(x)).\end{array}  \tag{16}\]  
where  
\[\sigma_{cz}^{2}(x)=E[(\varepsilon_{t}(x))^{2}], \tag{17}\]  
\[\sigma_{cz}^{2}(x)=E[(\varepsilon_{t-1}(x)-\widehat{c}_{t-1}(x))^{2}], \tag{18}\]  
Figure 2: Flexible Universal Kriging state-space model flow chart. \[\mathbf{p}_{t-1}(x)=E[(c_{t-1}(x)-\widehat{c}_{t-1}(x))(c_{t-1}-\widehat{c}_{t-1})], \tag{19}\]  
\[\mathbf{P}_{t-1}=E[(c_{t-1}-\widehat{c}_{t-1})(c_{t-1}-\widehat{c}_{t-1})^{T}], \tag{20}\]  
\[\mathbf{k}(x)=E[\epsilon_{t}(x)\epsilon_{t}], \tag{21}\]  
\[\mathbf{K}=E[\epsilon_{t}(\epsilon_{t})^{T}]. \tag{22}\]  
In this result, we know that \(\mathbf{v}_{t}\) is white noise, such that its expectation with the other values is 0. The terms \(\mathbf{q}_{t}^{T}(x)\mathbf{W}\mathbf{F}^{T}\mathbf{m}_{t}\) and \(\mathbf{f}^{T}(x)\mathbf{m}_{t}\) can be removed from the expression by Equation (14). Thus, the question now is how to compute \(\mathbf{q}_{t}(x)\). The \(\mathbf{q}_{t}(x)\) is derived by minimizing the variance of the prediction error, subject to \(\mathbf{F}\mathbf{W}_{t}^{T}\mathbf{q}_{t}(x)=\mathbf{f}(x)\), such that the prediction is unbiased. The constrained minimizing the variance problem can be solved by the Lagrange multiplier method:  
\[\xi(x)=\sigma_{ez}^{2}+\mathbf{\lambda}^{T}\mathbf{F}\mathbf{W}_{t}^{T}\mathbf{g}_{t}(x), \tag{23}\]  
where the Lagrange multipliers \(\mathbf{\lambda}\) are unknown and must be estimated. Setting the partial derivatives of this criterion with respect to \(\mathbf{g}_{t}(x)\) and \(\mathbf{\lambda}\) equal to zero, we obtain the following equations:  
\[\left\{ \begin{array}{c}2(\mathbf{W}_{t}\mathbf{P}_{t-1}\mathbf{W}_{t}^{T}+\mathbf{W}_{t}\bm {K}\mathbf{W}_{t}^{T}+\mathbf{R}_{t})\mathbf{q}_{t}(x)-2\mathbf{W}_{t}(\mathbf{p}_{t-1}(x)+\mathbf{k} (x))+\mathbf{W}_{t}\mathbf{F}^{T}\mathbf{\lambda}=0\\ \mathbf{F}\mathbf{W}_{t}^{T}\mathbf{q}_{t}(x)=\mathbf{f}(x)\end{array} .\right.