Spatio-Temporal Dual Kriging with Adaptive Coefficient Drift Function

1. Introduction

Spatio-temporal interpolation is a computational technique that estimates the value of an interested variable at a given location and time, provided available information in the sample at other locations and times. It is a combination of spatial interpolation, which deals with the estimation of values in geographic space, and temporal interpolation, which focuses on estimating values over a specific time period. A fundamental assumption that underlies spatio-temporal interpolation is that close observations exhibit greater similarity than distant ones [1]. Spatio-temporal interpolation methods can be categorized into two distinct approaches: the reduction approach and the extension approach [2,3,4]. Within the reduction approach, time is considered an independent dimension distinct from the spatial dimensions. This methodology is a two-step process in which data are first interpolated in space at each time step, and the resulting time series at each spatial location is subsequently interpolated. In contrast, the extension approach is formulated to interpolate simultaneously in both spatial and time dimensions. The reduction methods, such as shape function-based methods, inverse distance weighted (IDW), tension spline, and global polynomial interpolations, offer a simpler and less computationally demanding procedure. However, their suitability may be limited in scenarios where the temporal component of the data plays a significant role [5,6]. On the contrary, the extension methods, such as parallel IDW, random forests, and spatio-temporal kriging methods, are more complex and computationally intensive but possess enhanced capabilities for handling situations where the temporal component is substantial and can potentially yield higher accuracy, particularly in cases where the spatio-temporal relationships are complex [7].

Spatio-temporal kriging methodologies have emerged as powerful and effective geostatistical tools applied in several fields, including environmental sciences [8], pollution demography [9], and meteorology [10,11,12]. These methods expand upon spatial kriging techniques by incorporating both spatial and temporal information and enabling the construction of an optimal linear estimator that minimizes estimation variance under the constraint of unbiasedness. This estimator is a linear combination of spatio-temporal observed values and their corresponding weights. The two most commonly employed kriging techniques in the spatio-temporal domain are spatio-temporal simple kriging (ST-SK) and spatio-temporal ordinary kriging (ST-OK). These methods do not utilize any additional information or auxiliary variables during the estimation process. Their weights are determined by the covariance function or the variogram, which is dependent on the assumed degree of stationarity ascribed to the underlying random function. The ST-SK technique assumes second-order stationarity in the random function, characterized by a constant and known mean and a known covariance function that is dependent only on the separation vector between a pair of measured points in space and time. The ST-OK method, on the other hand, necessitates the assumption of intrinsic stationarity in the random function. This implies that the mean of the random function is constant and the variance of the difference is the same between the random function at two geographic locations lagged by the same distance [13,14].

An alternative spatio-temporal kriging incorporates additional variables, referred to as auxiliary variables, into the kriging model to enhance its predictive abilities, particularly when integrating relevant auxiliary variables that influence the target variable. The spatio-temporal regression kriging (ST-RK) technique involves decomposing the target variable into two distinct components: a deterministic trend component and a stochastic residual component [8]. The trend component encapsulates the systematic variation in the target variable, which can be modeled using auxiliary variables. Conversely, the residual component represents the remaining unexplained variability, which can be modelled using either ST-OK or ST-SK. Similar to the ST-RK method, the target variable in the spatio-temporal kriging with external drift (ST-KED) process can be divided into deterministic drift and stochastic residual terms. However, unlike the ST-RK method, the ST-KED weights are directly derived from the kriging equations, which are formulated under the unbiasedness and minimum error variance conditions [15].

In recent years, these two techniques have been successfully applied in diverse applications. The ST-RK method was used to predict chlorophyll-a [8], nitrogen dioxide (${\mathrm{NO}}_2$) concentrations [9], rainfall [10], and mean daily temperature [11]. Whereas, the ST-KED modeling technique was employed to interpolate soil water content [16], rainfall data [17], daily precipitation, and temperature [18]. Furthermore, extensive research has consistently shown that ST-KED outperforms ST-OK in terms of prediction accuracy and reliability [16,17,18]. Notwithstanding its superior performance, the ST-KED method is impeded by its computational complexity. Therefore, this study aims to tackle this limitation by developing an efficient algorithm for the ST-KED, thereby enhancing its applicability in broader practical scenarios.

In the context of spatial kriging methods, kriging with external drift (KED) is computationally expensive because it requires solving a new kriging system to determine the weights of the linear estimator for each interpolated point. In contrast, dual kriging (DK), a technique derived from KED, provides a unified approach to determining the weights of the estimator for all estimated points. Hence, DK is an advancement of KED that alleviates its computational challenges [19,20]. Although KED and DK produce the same estimated values, there is currently a dearth of research on spatio-temporal interpolation.

Motivated by this observation, we introduce a spatio-temporal dual kriging (ST-DK) algorithm to mitigate the computational burden of the ST-KED method by expanding upon the principles of spatial KED interpolation. The proposed ST-DK algorithm is designed to accommodate both fixed and adaptive coefficients of drift functions.

The remainder of this paper is structured as follows: Section 2 provides background information on ST-RK and ST-KED. Section 3 describes the construction of the proposed ST-DK methodology for fixed and adaptive coefficients of trend functions. In Section 4, a detailed description of the data and the selection criteria for auxiliary variables are presented. Subsequently, an evaluation of the proposed ST-DK method’s performance is conducted using meteorological data, including temperature and air pressure measurements. This evaluation then draws comparisons between the performance of the ST-DK method and the ST-RK method. Conclusions and discussions are drawn in Section 5.

2. Methodological Framework of Spatio-Temporal Kriging Incorporating Auxiliary Variables

Suppose that $\left\{Z\right(s,t)$, $(s,t)\in D\times \mathcal{T}\}$ is a real-value stochastic process or random function on the space-time domain $D\times \mathcal{T}$, where $D\subset {\mathbb{R}}^d$ and $\mathcal{T}\subset \mathbb{R}$ or $\mathbb{Z}$. The random variable $Z(s,t)$ can be disaggregated into two components [21,22]:

$$\begin{array}{c}\hfill Z(s,t)=\mu (s,t)+ϵ(s,t),\end{array}$$

(1)

where $\mu (s,t)$ denotes the deterministic drift function, also known as the trend function, while $ϵ(s,t)$ is a stochastic residual with a mean of zero.

The drift function can typically be represented as a linear combination of the functions evaluated at a spatio-temporal point $(s,t)$, as expressed by the following equation [10,22,23]:

$$\begin{array}{c}\hfill \mu (s,t)=\sum _{l=0}^L{\alpha }_l{f}_l\left(X(s,t)\right),\end{array}$$

(2)

where $X(s,t)={\left[X_1(s,t),...,X_{\eta }(s,t)\right]}^T$ is the vector of $\eta$ auxiliary variables at point $(s,t)$ and ${\alpha }_l$ is a non-zero constant coefficient of a known deterministic function $f_l$, which is represented in terms of auxiliary variables for $l=0,1,...,L$. The function $f_l$ is often referred to as a basis function in certain studies [24,25,26].

Given the assumption of second-order stationarity for the residual, the spatio-temporal covariance function $\mathbf{C}({\mathbf{h}}_s,h_t)$ exists. This covariance function quantifies the dependence between any pair of residuals separated by a spatial lag ${\mathbf{h}}_s$ and temporal lag $h_t$, defined as follows [27]:

$$\begin{array}{c}\hfill \mathbf{C}({\mathbf{h}}_s,h_t)=\mathrm{Cov}\left[ϵ(s+{\mathbf{h}}_s,t+h_t),ϵ(s,t)\right],\end{array}$$

(3)

where $(s+{\mathbf{h}}_s,t+h_t)$, $(s,t)\in D\times \mathcal{T}$. Furthermore, the spatio-temporal variogram of the residual is given by:

$$\begin{array}{cc}\hfill \gamma ({\mathbf{h}}_s,h_t)& =\frac{1}{2}\mathrm{Var}\left[ϵ(s+{\mathbf{h}}_s,t+h_t)-ϵ(s,t)\right],\hfill \\ \hfill & =\frac{1}{2}\mathbb{E}{\left[ϵ(s+{\mathbf{h}}_s,t+h_t)-ϵ(s,t)\right]}^2,\hfill \end{array}$$

(4)

which results in

$$\begin{array}{cc}\hfill \mathbf{C}({\mathbf{h}}_s,h_t)& =\mathbf{C}(\mathbf{0},0)-\gamma ({\mathbf{h}}_s,h_t),\hfill \end{array}$$

(5)

where $\mathbf{C}(\mathbf{0},0)$ is the global sill value of spatio-temporal variogram.

The classical empirical spatio-temporal variogram estimator [28], which estimates the theoretical spatio-temporal variogram function (presented in Equation (4)), can be formally defined as:

$$\begin{array}{c}\hfill {\gamma }^{*}({\mathbf{h}}_s,h_t)=\frac{1}{2\left|N_s\left({\mathbf{h}}_s\right)\right|\left|N_t\left(h_t\right)\right|}\sum _{s_i,s_u\in N_s\left({\mathbf{h}}_s\right)}\sum _{t_j,t_v\in N_t\left(h_t\right)}{\left(ϵ(s_i,t_j)-ϵ(s_u,t_v)\right)}^2,\end{array}$$

(6)

where $N_s\left({\mathbf{h}}_s\right)$ represents the set of all spatial location pairs separated by the spatial lag vector ${\mathbf{h}}_s$ and $N_t\left(h_t\right)$ denotes the set of all time point pairs separated by the temporal lag $h_t$. The sizes of these sets, $\left|N_s\left({\mathbf{h}}_s\right)\right|$ and $\left|N_t\left(h_t\right)\right|$, indicate the number of distinct pairs in each set.

Several spatio-temporal covariance models have been proposed in the literature, including the sum model [29], the metric model [30], and the product model [27]. The present study focuses on the product-sum model, which is expressed as follows [31]:

$$\begin{array}{c}\hfill \mathbf{C}({\mathbf{h}}_s,h_t)=k_1{\mathbf{C}}_s\left({\mathbf{h}}_s\right){\mathbf{C}}_t\left(h_t\right)+k_2{\mathbf{C}}_s\left({\mathbf{h}}_s\right)+k_3{\mathbf{C}}_t\left(h_t\right),\end{array}$$

(7)

where ${\mathbf{C}}_s$ and ${\mathbf{C}}_t$ are spatial and temporal covariance functions, respectively. The spatial and temporal covariance functions are defined as

$$\begin{array}{cc}\hfill {\mathbf{C}}_s\left({\mathbf{h}}_s\right)& =\mathbf{C}({\mathbf{h}}_s,0)={\mathbf{C}}_s\left(\mathbf{0}\right)-{\gamma }_s\left({\mathbf{h}}_s\right),\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill {\mathbf{C}}_t\left(h_t\right)& =\mathbf{C}(\mathbf{0},h_t)={\mathbf{C}}_t\left(0\right)-{\gamma }_t\left(h_t\right),\hfill \end{array}$$

(9)

where ${\mathbf{C}}_s\left(\mathbf{0}\right)$ and ${\mathbf{C}}_t\left(\mathbf{0}\right)$ are the values of sill of the spatial variogram, ${\gamma }_s$, and temporal variogram, ${\gamma }_t$, respectively. This study uses the exponential model [32], a parametric variogram, to represent both the spatial and temporal variograms, which are fit to the corresponding empirical variograms by weighted least squares [33].

Based on the relationship between the covariance function and variogram established in Equation (5), (7) - (9), the spatio-temporal variogram can be expressed as:

$$\begin{array}{c}\hfill \gamma ({\mathbf{h}}_s,h_t)=[k_1{\mathbf{C}}_t\left(0\right)+k_2]{\gamma }_s\left({\mathbf{h}}_s\right)+[k_1{\mathbf{C}}_s\left(\mathbf{0}\right)+k_3]{\gamma }_t\left(h_t\right)-k_1{\gamma }_s\left({\mathbf{h}}_s\right){\gamma }_t\left(h_t\right),\end{array}$$

(10)

where

$$\begin{array}{cc}\hfill k_1& =\frac{{\mathbf{C}}_s\left(\mathbf{0}\right)+{\mathbf{C}}_t\left(0\right)-\mathbf{C}(\mathbf{0},0)}{{\mathbf{C}}_s\left(\mathbf{0}\right){\mathbf{C}}_t\left(0\right)},\hfill \\ \hfill k_2& =\frac{\mathbf{C}(\mathbf{0},0)-{\mathbf{C}}_t\left(0\right)}{{\mathbf{C}}_s\left(\mathbf{0}\right)},\hfill \\ \hfill k_3& =\frac{\mathbf{C}(\mathbf{0},0)-{\mathbf{C}}_s\left(\mathbf{0}\right)}{{\mathbf{C}}_t\left(0\right)}.\hfill \end{array}$$

The parameters of the spatio-temporal variogram model can be determined by fitting the empirical variogram provided in Equation (7). Subsequently, the spatio-temporal covariance function of residuals is derived by leveraging the relationship between the variogram and covariance function as provided in Equation (5).

The objective of the spatio-temporal kriging with drift function is to estimate the value of the target variable at an unobserved spatio-temporal point $(s_0,t_0)$, represented as $Z(s_0,t_0)$, where $s_0\in D$ and $t_0\in \mathcal{T}$. The estimation is achieved using the principles of unbiasedness and minimum error variance. Two classical spatio-temporal kriging methods, spatio-temporal regression kriging (ST-RK) and spatio-temporal kriging with external drift (ST-KED), are covered in detail in the following.

2.1. Spatio-Temporal Regression Kriging

The ST-RK estimator briefly described in the previous section consists of drift function and residual components, in which the implementation of each term is obtained independently. The drift function is conventionally modeled using multiple linear regression (MLR), while the residuals are interpolated using spatio-temporal simple kriging (ST-SK). Given a set of $n\times m$ observed values of the target variables and auxiliary variables at distinct spatio-temporal points $(s_1,t_1),...,(s_n,t_m)$, the ST-RK estimator at the point $(s_0,t_0)$ denoted as $Z_{\mathrm{ST}-\mathrm{RK}}(s_0,t_0)$, can be written as [23]:

$$\begin{array}{cc}\hfill Z_{\mathrm{ST}-\mathrm{RK}}(s_0,t_0)& =\mu (s_0,t_0)+{ϵ}_{\mathrm{ST}-\mathrm{SK}}(s_0,t_0),\hfill \\ \hfill & =\sum _{l=0}^L{\alpha }_l{f}_l\left(X(s_0,t_0)\right)+\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}^{{ϵ}_0}ϵ(s_i,t_j),\hfill \\ \hfill & ={\alpha }^T{\mathit{F}}_0+{ϵ}^T{\mathsf{\Lambda }}^{{ϵ}_0},\hfill \end{array}$$

(11)

where $\alpha ={\left[{\alpha }_0,...,{\alpha }_L\right]}^T$ is a vector of coefficients,

${\mathit{F}}_0={\left[f_0\left(X(s_0,t_0)\right),...,f_L\left(X(s_0,t_0)\right)\right]}^T$ is a vector of the basis functions,

$ϵ={\left[ϵ(s_1,t_1),...,ϵ(s_n,t_m)\right]}^T$ is a residual vector, and

${\mathsf{\Lambda }}^{{ϵ}_0}={\left[{\lambda }_{11}^{{ϵ}_0},...,{\lambda }_{nm}^{{ϵ}_0}\right]}^T$ is the ST-SK weight vector.

The optimal weight vector for the ST-SK method is determined by minimizing the prediction error variance under the constraint of unbiasedness. This minimization problem is solved utilizing the Lagrange multiplier method, resulting in a system of linear equations known as the ST-SK system [34]:

$$\begin{array}{c}\hfill \sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}^{{ϵ}_0}\mathbf{C}\left(s_u-s_i,t_v-t_j\right)=\mathbf{C}\left(s_u-s_0,t_v-t_0\right),foru=1,...,n\mathrm{and}v=1,...,m.\end{array}$$

(12)

The ST-SK equations can be expressed in the following matrix form:

$$\begin{array}{c}\hfill \mathsf{\Omega }{\mathsf{\Lambda }}^{{ϵ}_0}={\mathsf{\Omega }}_0,\end{array}$$

(13)

where $\mathsf{\Omega }=\left[\begin{array}{ccc}\mathbf{C}\left(s_1-s_1,t_1-t_1\right)& \cdots & \mathbf{C}\left(s_1-s_n,t_1-t_m\right)\\ ⋮& \ddots & ⋮\\ \mathbf{C}\left(s_n-s_1,t_m-t_1\right)& \cdots & \mathbf{C}\left(s_n-s_n,t_m-t_m\right)\end{array}\right]$ and

   ${\mathsf{\Omega }}_0=\left[\begin{array}{c}\mathbf{C}\left(s_1-s_0,t_1-t_0\right)\\ ⋮\\ \mathbf{C}\left(s_n-s_0,t_m-t_0\right)\end{array}\right]$.

2.2. Spatio-Temporal Kriging with External Drift

Let $Z_{\mathrm{ST}-\mathrm{KED}}(s_0,t_0)$ be the estimated value of the target variable at the point $(s_0,t_0)$. The ST-KED estimator is expressed as a linear combination of input data, $Z(s_i,t_j)$ for $i=1,...,n$ and $j=1,...,m$ [34]:

$$\begin{array}{cc}\hfill Z_{\mathrm{ST}-\mathrm{KED}}(s_0,t_0)& =\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}Z(s_i,t_j),\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill & ={\mathit{Z}}^T\mathsf{\Lambda },\hfill \end{array}$$

(15)

where $\mathit{Z}={\left[Z(s_1,t_1),...,Z(s_n,t_m)\right]}^T$ is the vector of observed data for target variable and

$\mathsf{\Lambda }={\left[{\lambda }_{11},...,{\lambda }_{nm}\right]}^T$ is the ST-KED weight vector.

The optimal weight vector for the ST-KED technique can be found using the Lagrange multiplier method to minimize the prediction error variance while satisfying the unbiasedness constraint. This approach is based on the ST-KED estimator and the form of the target variable Z, as represented by Equation (1). The ST-KED system is written as [15]:

$$\begin{array}{c}\hfill \left[\begin{array}{ccc}\mathsf{\Omega }& |& \mathit{F}\\ --& +& --\\ {\mathit{F}}^T& |& {\mathbf{0}}_{L+1}\end{array}\right]\left[\begin{array}{c}\mathsf{\Lambda }\\ --\\ \mathit{N}\end{array}\right]=\left[\begin{array}{c}{\mathsf{\Omega }}_0\\ --\\ {\mathit{F}}_0\end{array}\right],\end{array}$$

(16)

where ${\mathbf{0}}_{L+1}$ is the $(L+1)\times (L+1)$ zero matrix,

$\mathit{N}={\left[{\eta }_0,...,{\eta }_L\right]}^T$ represents the vector of Lagrange multipliers,

$\mathit{F}=\left[\begin{array}{c}{\mathit{F}}_1\\ ⋮\\ {\mathit{F}}_m\end{array}\right]$, and

${\mathit{F}}_j=\left[\begin{array}{ccc}{f}_0\left(X(s_1,t_j)\right)& \cdots & f_L\left(X(s_1,t_j)\right)\\ ⋮& \ddots & ⋮\\ f_0\left(X(s_n,t_j)\right)& \cdots & f_L\left(X(s_n,t_j)\right)\end{array}\right]$, for $j=1,...,m$.

Since the ST-KED weights need to be recalculated for each unobserved data point, this leads to a significant increase in computational resources. Employing a single weight matrix to estimate the target variable across all ungauged locations can effectively address this challenge. The spatio-temporal dual kriging (ST-DK) method is therefore proposed to overcome this limitation.

3. Methodology of Spatio-Temporal Dual Kriging

The formulation of the ST-DK interpolation is inspired by the principles of the DK [35]. The ST-DK algorithm can be described as a refinement of the ST-KED method, in which the ST-DK weight coefficient matrix is performed only once for every unobserved site. In this section, the derivation of both fixed and adaptive coefficients of the ST-DK estimation is presented. To simplify subsequent references, we will use the abbreviations FST-DK and AST-DK to designate the ST-DK method with a fixed coefficient drift function and the ST-DK method with an adaptive coefficient drift function, respectively.

3.1. Spatio-Temporal Dual Kriging with Fixed Coefficient Drift Function

Given the ST-KED system represented by Equation (16), where the coefficient matrix is assumed to be invertible, the solution can be expressed in the following form [35,36]:

$$\begin{array}{c}\hfill \left[\begin{array}{c}\mathsf{\Lambda }\\ --\\ \mathit{N}\end{array}\right]=\left[\begin{array}{ccc}\mathit{U}& |& \mathit{V}\\ --& +& --\\ {\mathit{V}}^T& |& \mathit{W}\end{array}\right]\left[\begin{array}{c}{\mathsf{\Omega }}_0\\ --\\ {\mathit{F}}_0\end{array}\right],\end{array}$$

(17)

where $\mathit{U}$, $\mathit{V}$, and $\mathit{W}$ are matrices with dimensions of $nm\times nm$, $nm\times (L+1)$, and $(L+1)\times (L+1)$, respectively. Substituting the weight solution of Equation (17) for the weight vector $\mathsf{\Lambda }$ in Equation (15) yields the following expression for the estimated value:

$$\begin{array}{c}\hfill Z_{\mathrm{ST}-\mathrm{KED}}(s_0,t_0)={\mathit{Z}}^T{\mathit{VF}}_0+{\mathit{Z}}^T\mathit{U}{\mathsf{\Omega }}_0.\end{array}$$

(18)

To establish a formal representation of the dual formulation of ST-KED, two matrices are introduced

$$\begin{array}{c}\hfill {\mathit{a}}^T=\left[a_0,...,a_L\right]={\mathit{Z}}^T\mathit{V}\hspace{1em}\hspace{1em}\mathrm{and}\hspace{1em}\hspace{1em}{\mathit{b}}^T=\left[b_{11},...,b_{nm}\right]={\mathit{Z}}^T\mathit{U}.\end{array}$$

(19)

Consequently, Equation (18) undergoes a transformation, resulting in a novel form known as the FST-DK estimator, as depicted below:

$$\begin{array}{cc}\hfill Z_{\mathrm{FST}-\mathrm{DK}}(s_0,t_0)& =\sum _{l=0}^L{a}_l{f}_l\left(X(s_0,t_0)\right)+\sum _{j=1}^m\sum _{i=1}^n{b}_{ij}\mathbf{C}\left(s_i-s_0,t_j-t_0\right),\hfill \\ \hfill & ={\mathit{a}}^T{\mathit{F}}_0+{\mathit{b}}^T{\mathsf{\Omega }}_0.\hfill \end{array}$$

(20)

As a result, the vectors $\mathit{a}$ and $\mathit{b}$ become the weight vectors of the FST-DK estimator, exhibiting invariance with respect to the target point $(s_0,t_0)$. According to Equation (19), the weight vectors $\mathit{a}$ and $\mathit{b}$ can be represented in matrix form as:

$$\begin{array}{c}\hfill \left[\begin{array}{c}\mathit{b}\\ --\\ \mathit{a}\end{array}\right]=\left[\begin{array}{ccc}\mathit{U}& |& \mathit{A}\\ --& +& --\\ {\mathit{V}}^T& |& \mathit{B}\end{array}\right]\left[\begin{array}{c}\mathit{Z}\\ --\\ {\underline{\mathbf{0}}}_{L+1}\end{array}\right],\end{array}$$

(21)

where $\mathit{A}$, $\mathit{B}$ are arbitrary matrices with dimensions of $nm\times (L+1)$ and $(L+1)\times (L+1)$, respectively, and ${\underline{\mathbf{0}}}_{L+1}$ represents the zero vector of size $L+1$.

Given the substitutions $\mathit{A}=\mathit{V}$ and $\mathit{B}=\mathit{W}$, the vectors of coefficients $\mathit{a}$ and $\mathit{b}$ are solutions of the FST-DK system,

$$\begin{array}{c}\hfill \left[\begin{array}{ccc}\mathsf{\Omega }& |& \mathit{F}\\ --& +& --\\ {\mathit{F}}^T& |& {\mathbf{0}}_{L+1}\end{array}\right]\left[\begin{array}{c}\mathit{b}\\ --\\ \mathit{a}\end{array}\right]=\left[\begin{array}{c}\mathit{Z}\\ --\\ {\underline{\mathbf{0}}}_{L+1}\end{array}\right].\end{array}$$

(22)

The estimated value of the target variable at the unobserved point $(s_0,t_0)$ is therefore obtained by substituting the vectors $\mathit{a}$ and $\mathit{b}$ into Equation (20).

3.2. Spatio-Temporal Dual Kriging with Adaptive Coefficient Drift Function

In the previous case, the coefficient ${\alpha }_l$, where $l=0,1,...,L$, of the drift function $\mu (s,t)$ in Equation (2) remain constant across all time steps. This property of the coefficients makes it simple to use in practice, as it does not require us to keep track of its value at each time step. Although this can be a significant advantage, especially when dealing with large datasets, but it is a trade-off between the estimation accuracy and a reduction in computational cost. To address this limitation and incorporate a more flexible drift function, this study proposes a novel drift function approach that incorporates dynamically adjusted coefficients. The proposed drift function, denoted by ${\mu }_A$, is hereby introduced as follows:

$$\begin{array}{c}\hfill {\mu }_A(s,t)=\sum _{l=0}^L\sum _{\tau =t_1}^{t_m}{\alpha }_{l\tau }{f}_{l\tau }\left(X(s,t)\right),\end{array}$$

(23)

where $f_{l\tau }$ is l-th basis function of auxiliary variables at time step $\tau$ with ${\alpha }_{l\tau }\in \mathbb{R}-\left\{0\right\}$ being the constant coefficient. Furthermore, $t_1$ and $t_m$ represent the initial and final time steps, respectively.

For $l=0,1,...,L$, function $f_{l\tau }$ is represented in term of auxiliary variables and is defined by the equation:

$$\begin{array}{c}\hfill f_{l\tau }\left(X(s,t)\right)=\left\{\begin{array}{cc}{f}_l\left(X(s,t)\right),& \tau =t,\\ 0,& \tau \ne t.\end{array}\right.\end{array}$$

(24)

This expression of $f_{l\tau }$ indicates that the drift function ${\mu }_A$ can change over time step $\tau$.

Incorporating these adaptive coefficients into the drift function of the ST-KED, the value at the target point $(s_0,t_0)$ is estimated using the ST-KED estimator, as presented in Equation (14). This estimator is denoted as ${\widehat{Z}}_{\mathrm{ST}-\mathrm{KED}}(s_0,t_0)$. Moreover, the matrix equation that arises from employing the Lagrange multiplier method to minimize the variance of prediction error under the constraint of unbiasedness is expressed as follows (Refer to Appendix A for a detailed explanation of the derivation process):

$$\begin{array}{c}\hfill \left[\begin{array}{ccc}\mathsf{\Omega }& |& \mathcal{F}\\ --& +& --\\ {\mathcal{F}}^T& |& {\mathbf{0}}_{(L+1)m}\end{array}\right]\left[\begin{array}{c}\mathsf{\Lambda }\\ --\\ \mathcal{N}\end{array}\right]=\left[\begin{array}{c}{\mathsf{\Omega }}_0\\ --\\ {\mathcal{F}}_0\end{array}\right],\end{array}$$

(25)

where $\mathcal{N}={\left[{\eta }_{0t_1},...,{\eta }_{L{t}_m}\right]}^T$ is a vector of the Lagrange multiplier,

$\mathcal{F}=\left[\begin{array}{cccc}{\mathit{F}}_1& {\mathbf{0}}_{n\times (L+1)}& \cdots & {\mathbf{0}}_{n\times (L+1)}\\ {\mathbf{0}}_{n\times (L+1)}& {\mathit{F}}_2& \cdots & {\mathbf{0}}_{n\times (L+1)}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{0}}_{n\times (L+1)}& {\mathbf{0}}_{n\times (L+1)}& \cdots & {\mathit{F}}_m\end{array}\right]$,

${\mathbf{0}}_{n\times (L+1)}$ represents the $n\times (L+1)$ zero matrix, and

${\mathcal{F}}_0=\left[\begin{array}{c}{f}_{0t_1}\left(X\left(s_0,t_0\right)\right)\\ ⋮\\ f_{L{t}_m}\left(X\left(s_0,t_0\right)\right)\end{array}\right]$, see detail in Appendix A.

By applying the same procedure as in SubSection 2.2 and SubSection 3.1, the AST-DK estimator takes the following form:

$$\begin{array}{cc}\hfill Z_{\mathrm{AST}-\mathrm{DK}}(s_0,t_0)& =\sum _{l=0}^L\sum _{\tau =t_1}^{t_m}{a}_{l\tau }{f}_{l\tau }\left(X(s_0,t_0)\right)+\sum _{j=1}^m\sum _{i=1}^n{b}_{ij}\mathbf{C}\left(s_i-s_0,t_j-t_0\right),\hfill \\ \hfill & ={\widehat{\mathit{a}}}^T{\mathcal{F}}_0+{\mathit{b}}^{\mathbf{T}}{\mathsf{\Omega }}_0,\hfill \end{array}$$

(26)

where $\widehat{\mathit{a}}$ and $\mathit{b}$ are the vector of coefficients such that ${\widehat{\mathit{a}}}^T=\left[a_{0t_1},...,a_{L{t}_m}\right]$ and ${\mathit{b}}^T=\left[b_{11},...,b_{nm}\right]$.

Analogously to the approach outlined in SubSection 3.1, we can derive the AST-DK system for determining the vectors $\widehat{\mathit{a}}$ and $\mathit{b}$, which is given by

$$\begin{array}{c}\hfill \left[\begin{array}{ccc}\mathsf{\Omega }& |& \mathcal{F}\\ --& +& --\\ {\mathcal{F}}^T& |& {\mathbf{0}}_{(L+1)m}\end{array}\right]\left[\begin{array}{c}\mathit{b}\\ --\\ \widehat{\mathit{a}}\end{array}\right]=\left[\begin{array}{c}\mathit{Z}\\ --\\ {\underline{\mathbf{0}}}_{(L+1)m}\end{array}\right].\end{array}$$

(27)

To evaluate the effectiveness of the proposed AST-DK techniques, these methods are applied in the subsequent section to interpolate meteorological data, encompassing temperature and air pressure. Additionally, the performance of these methods will be compared to that of the ST-RK with fixed coefficient drift function (denoted as FST-RK) and the ST-RK with adaptive coefficient drift function method (denoted as AST-RK).

4. Application of the Proposed ST-DK Techniques for Temperature and Air Pressure Interpolations in Thailand

4.1. Study Area and Data Description

Our domain of study covers the entire country of Thailand, which lies between $5^{\circ }{37}^{\prime }$ N and ${20}^{\circ }{27}^{\prime }$ N latitude and ${97}^{\circ }{22}^{\prime }$ E and ${105}^{\circ }{37}^{\prime }$ E longitude. The total area of the country is 513,115 square kilometres, with a coastline of 3,219 kilometres [37,38]. The present study leverages observational data from the Hydro-Informatics Institute (HII) [39] including temperature (in °C), air pressure (in hPa), relative humidity (in %RH), and the Digital Elevation Model (DEM) (in metre). A comprehensive network of 226 weather stations encloses the study area, as depicted in Figure 1. The performances of the kriging schemes are evaluated based on monthly mean data over the period of January 2018 to December 2018. Prior to the interpolation process, data cleaning and preprocessing procedures are conducted to ensure data quality.

4.2. Selection of Auxiliary Variables

The estimation accuracy of both the ST-RK and ST-DK methods mainly depends on the careful selection of auxiliary variables for which a strong correlation with target variables is found [40]. This stage involves an interaction assessment between external factors and target variables using the Pearson correlation coefficient [41]. If the absolute value of the correlation coefficient exceeds a predefined threshold of 0.5 (as established in [42]), the candidate variable is considered statistically significant and deemed suitable to include in the basis functions of the drift terms. In this work, air pressure, relative humidity, and DEM are considered potential auxiliary variables for the temperature interpolation, whereas temperature, relative humidity, and DEM are those of the air pressure estimation.

4.3. Accuracy Assessment

A comparison of the performances between the proposed ST-DK and the traditional ST-RK methods is demonstrated in the case of the fixed and adaptive coefficient drift functions. The accuracy of all algorithms is evaluated using a $10-$folds cross validation technique. The process is done by dividing the entire dataset into ten equal-sized folds. Each fold is then used once as the test set, while the remaining 9 folds serve as the training set [43]. After each of the 10 folds was once the left-out fold, the overall accuracy was computed by averaging prediction error values obtained from test folds. The root mean square error (RMSE) [44] and mean absolute percentage error (MAPE) [45] are the performance evaluation statistics, which can be written as follows:

$$\mathrm{RMSE}=\sqrt{\frac{1}{N}\sum _{n=1}^N{\left(Z_n-Z_n^{*}\right)}^2}$$

and

$$\mathrm{MAPE}=\frac{1}{N}\sum _{n=1}^N\frac{|Z_n-Z_n^{*}|}{Z_n}\times 100,$$

where N is the number of the observed data in each fold and $Z_n$ and $Z_n^{*}$ are the n-th actual and estimated values, respectively.

4.4. Results

4.4.1. Comparison of ST-RK and ST-DK with Fixed Coefficient Drift Function

Table 1 presents the correlation coefficients between the target and potential auxiliary variables. It is found that air pressure is correlated with temperature more significantly than others, with the correlation coefficients being 0.5775. On the other hand, a significant negative correlation is also observed between air pressure and DEM. Consequently, air pressure and temperature are selected as each other’s auxiliary variables, with DEM being another auxiliary variable for air pressure estimation.

The drift function used in the temperature estimation process takes the following form:

$$\begin{array}{c}\hfill \mu (s,t)={\alpha }_0+{\alpha }_1{X}_1(s,t),\end{array}$$

(28)

where $X_1(s,t)$ is the air pressure at point $(s,t)$. By comparing the given form with Equation (2), we observe that $f_0\left(X_1(s,t)\right)=1$ and $f_1\left(X_1(s,t)\right)=X_1(s,t)$ for all point $(s,t)$. The drift function for air pressure estimation is represented by

$$\begin{array}{c}\hfill \mu (s,t)={\alpha }_0+{\alpha }_1{X}_2(s,t)+{\alpha }_2{X}_3(s,t),\end{array}$$

(29)

where $X_2(s,t)$ and $X_3(s,t)$ denote temperature and DEM at point $(s,t)$, respectively. Suppose that $X(s,t)={\left[X_2(s,t),X_3(s,t)\right]}^T$ and through a direct comparison with Equation (2), it follows that $f_0\left(X(s,t)\right)=1$, $f_1\left(X(s,t)\right)=X_2(s,t)$, and $f_2\left(X(s,t)\right)=X_3(s,t)$. All unknown parameters of the drift functions are estimated through ordinary least squares (OLS) [46].

The RMSE and MAPE values attained from the FST-RK and FST-DK methods are summarized in Table 2. Regarding accuracy measures, the FST-RK exhibited superior performance in temperature estimation compared to the FST-DK. The FST-RK method reduces the RMSE and MAPE by 2.582% and 2.025%, respectively. On the contrary, the FST-DK excelled in air pressure estimation, achieving lower RMSE (15.6068%) and MAPE (4.4717%) compared to the benchmark scheme.

4.4.2. Comparison of the ST-RK and ST-DK with Adaptive Coefficient Drift Function

Since trends are assumed to be changing across time in the present case, the procedure to compute overall correlation coefficients therefore differs from the fixed coefficients case. Here, monthly correlation coefficients between variables are averaged to obtain an overall correlation. A highly positive correlated relationship is marked between temperature and air pressure, with a correlation coefficient exceeding 0.8. The DEM variable is also included in the drift function for both temperature and air pressure estimations as the absolute values of correlation coefficients are greater than 0.5.

Table 3. Monthly means of Pearson correlation coefficients between candidate auxiliary and target variables.

Variable	Temperature	Air Pressure
Temperature	1.0000	0.8021
Air Pressure	0.8021	1.0000
Relative Humidity	-0.0710	0.3432
DEM	-0.5907	-0.7608

Table 3. Monthly means of Pearson correlation coefficients between candidate auxiliary and target variables.

Variable	Temperature	Air Pressure
Temperature	1.0000	0.8021
Air Pressure	0.8021	1.0000
Relative Humidity	-0.0710	0.3432
DEM	-0.5907	-0.7608

Therefore, we define the drift function based on the selected auxiliary variables for temperature estimation. In this work, we define two such functions below, utilizing the variable definitions stated in the previous subsection. The first function is formulated as:

$$\begin{array}{c}\hfill {\mu }_A(s,t)=\sum _{\tau =1}^{12}{\alpha }_{0\tau }{f}_{0\tau }\left(X_1(s,t)\right)+\sum _{\tau =1}^{12}{\alpha }_{1\tau }{f}_{1\tau }\left(X_1(s,t)\right),\end{array}$$

(30)

where $f_{0\tau }\left(X_1(s,t)\right)=\left\{\begin{array}{cc}1,& \tau =t,\\ 0,& \tau \ne t\end{array}\right.$ and

  $f_{1\tau }\left(X_1(s,t)\right)=\left\{\begin{array}{cc}{X}_1(s,t),& \tau =t,\\ 0,& \tau \ne t\end{array}\right.$

for $t=1,2,...,12$.

Furthermore, suppose that $\mathbf{X}(s,t)={\left[X_1(s,t),X_3(s,t)\right]}^T$, a second drift function is expressed as:

$$\begin{array}{c}\hfill {\mu }_A(s,t)=\sum _{\tau =1}^{12}{\alpha }_{0\tau }{f}_{0\tau }\left(\mathbf{X}(s,t)\right)+\sum _{\tau =1}^{12}{\alpha }_{1\tau }{f}_{1\tau }\left(\mathbf{X}(s,t)\right)+\sum _{\tau =1}^{12}{\alpha }_{2\tau }{f}_{2\tau }\left(\mathbf{X}(s,t)\right),\end{array}$$

(31)

where $f_{0\tau }\left(\mathbf{X}(s,t)\right)=\left\{\begin{array}{cc}1,& \tau =t,\\ 0,& \tau \ne t,\end{array}\right.$

  $f_{1\tau }\left(\mathbf{X}(s,t)\right)=\left\{\begin{array}{cc}{X}_1(s,t),& \tau =t,\\ 0,& \tau \ne t,\end{array}\right.$ and

  $f_{2\tau }\left(\mathbf{X}(s,t)\right)=\left\{\begin{array}{cc}{X}_3(s,t),& \tau =t,\\ 0,& \tau \ne t,\end{array}\right.$

for $t=1,2,...,12$.

For air pressure estimation, the drift function is formulated as:

$$\begin{array}{c}\hfill {\mu }_A(s,t)=\sum _{\tau =1}^{12}{\alpha }_{0\tau }{f}_{0\tau }\left(X(s,t)\right)+\sum _{\tau =1}^{12}{\alpha }_{1\tau }{f}_{1\tau }\left(X(s,t)\right)+\sum _{\tau =1}^{12}{\alpha }_{2\tau }{f}_{2\tau }\left(X(s,t)\right),\end{array}$$

(32)

where $f_{0\tau }\left(X(s,t)\right)=\left\{\begin{array}{cc}1,& \tau =t,\\ 0,& \tau \ne t,\end{array}\right.$

  $f_{1\tau }\left(X(s,t)\right)=\left\{\begin{array}{cc}{X}_2(s,t),& \tau =t,\\ 0,& \tau \ne t,\end{array}\right.$ and

  $f_{2\tau }\left(X(s,t)\right)=\left\{\begin{array}{cc}{X}_3(s,t),& \tau =t,\\ 0,& \tau \ne t,\end{array}\right.$

for $t=1,2,...,12$.

At each monthly time step, the unknown coefficients are identified through the application of the OLS method, as in the previous subsection.

Table 4 demonstrates the estimation performances of the AST-RK and AST-DK methods for temperature and air pressure based on the drift functions given in Equations (30)-(32). As evidenced by the RMSE and MAPE metrics, the AST-DK method with one and two auxiliary variables consistently outperform the AST-RK method for estimation of both target variables. In particular, the RMSE values attained from the AST-DK method with two auxiliary variables are 0.8336 for temperature and 10.5835 for air pressure, resulting in 12.5544% and 7.4256% improvement compared to the AST-RK method. Furthermore, the AST-RK and AST-DK methods with one auxiliary variable surpass the accuracy of the fixed-coefficient approaches presented in Table 2.

Besides, the role of the adaptive coefficient drift function in the ST-DK method becomes apparent when considering temperature estimation with a single auxiliary variable. In the case of fixed coefficients, the ST-RK model yields higher accuracy estimates than the ST-DK, but the opposite scenario occurs for the adaptive coefficient schemes.

For illustration purposes, the spatial distribution maps of temperature and air pressure obtained from the AST-RK and AST-DK models are generated. The maps were gridded at a resolution of 0.05 degrees, corresponding to an approximate spatial area of 5.5 $k{m}^2$ per grid cell. Thailand’s weather is divided into mainly three seasons; the summer season (March to May), the rainy season (May to October), and the winter season (November to February) [47,48]. One month of each season (March, July, and November) was performed to represent seasonal meteorological fluctuations across the country. The temperature variations across Thailand are visualized in Figure 2. As can be seen in the figures, stronger temperature gradients are produced by the ST-DK method for all three representative months. A more dispersed distribution of high temperatures can be markedly observed in March from the north to the centre regions, as displayed in Figure 2(a).

A comparison of air pressure distributions obtained from the AST-RK and AST-DK methods is presented in Figure 3. The results of both methods show similar distribution patterns for March and November. A notable difference can be found in July, in which a relatively low air pressure produced by the AST-RK model is identified in the northern area as opposed to the AST-DK approach.

5. Conclusions and Discussions

This work introduces the ST-DK technique, a novel methodological approach designed as a computationally efficient alternative to the ST-KED method for spatio-temporal interpolation. The ST-DK leverages the ability to utilize fixed coefficients within its drift functions while offering increased flexibility to incorporate adaptive coefficients. This capability enables the ST-DK to effectively model complex spatio-temporal dynamics. The results of the case study revealed that the ST-DK and ST-RK methods with adaptive coefficient drift functions demonstrated superior performance in the estimation of both temperature and air pressure compared to their counterparts employing fixed coefficients. Furthermore, the ST-DK method consistently outperformed the ST-RK method for both temperature and air pressure estimations. The advantages of the ST-DK model in estimation over the ST-RK method may be attributed to the fact that correlated residuals are considered in the DK trend process, which is not the case for the RK model. Despite the satisfactory results of this study, higher computational resources for a large number of unobserved locations and estimated time steps can be a challenge. To handle the issue of the growing dimensionality of the ST-DK system, specialized solvers designed for large linear systems are used to exploit the sparsity and symmetry of the ST-DK coefficient matrix [49,50,51]. Another practical way to boost simulation efficiency is by localizing ST-DK computations to a neighborhood of each unobserved site, resulting in a model size reduction. Possible future research endeavors include the identification of an optimal drift function and the development of techniques for determining the optimal coefficients within adaptive-coefficient drift functions tailored for the ST-DK framework.

Appendix A.1. Unbiasedness Condition

The estimator ${\widehat{Z}}_{\mathrm{ST}-\mathrm{KED}}(s_0,t_0)$ is said to be unbiased estimator for $Z\left(s_0,t_0\right)$ if the expected value of estimation error is zero [35]. which can be expressed mathematically as:

$$\begin{array}{c}\hfill \mathbb{E}\left[{\widehat{Z}}_{\mathrm{ST}-\mathrm{KED}}(s_0,t_0)-Z\left(s_0,t_0\right)\right]=0.\end{array}$$

(A1)

Utilizing the expression of Equation (1) and the ST-KED estimator presented in Equation (14) with adaptive coefficient drift function as expressed in Equation (23), this expected value can be expanded as:

$$\begin{array}{cc}\hfill \mathbb{E}\left[{\widehat{Z}}_{\mathrm{ST}-\mathrm{KED}}(s_0,t_0)-Z\left(s_0,t_0\right)\right]& =\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}\left(\mathbb{E}\left[{\mu }_A\left(s_i,t_j\right)\right]+\mathbb{E}\left[ϵ\left(s_i,t_j\right)\right]\right)-\left(\mathbb{E}\left[{\mu }_A\left(s_0,t_0\right)\right]+\mathbb{E}\left[ϵ\left(s_0,t_0\right)\right]\right),\hfill \\ \hfill & =\left(\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}{\mu }_A\left(s_i,t_j\right)\right)-{\mu }_A\left(s_0,t_0\right),\hfill \\ \hfill & =\left(\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}\sum _{l=0}^L\sum _{\tau =t_1}^{t_m}{\alpha }_{l\tau }{f}_{l\tau }\left(X(s_i,t_j)\right)\right)-\left(\sum _{l=0}^L\sum _{\tau =t_1}^{t_m}{\alpha }_{l\tau }{f}_{l\tau }\left(X(s_0,t_0)\right)\right),\hfill \\ \hfill & =\left(\sum _{l=0}^L\sum _{\tau =t_1}^{t_m}{\alpha }_{l\tau }\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}{f}_{l\tau }\left(X(s_i,t_j)\right)\right)-\left(\sum _{l=0}^L\sum _{\tau =t_1}^{t_m}{\alpha }_{l\tau }{f}_{l\tau }\left(X(s_0,t_0)\right)\right),\hfill \\ \hfill & =\sum _{l=0}^L\sum _{\tau =t_1}^{t_m}{\alpha }_{l\tau }\left(\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}{f}_{l\tau }\left(X(s_i,t_j)\right)-f_{l\tau }\left(X(s_0,t_0)\right)\right).\hfill \end{array}$$

(A2)

To guarantee unbiasedness condition, the following condition must be satisfied:

$$\begin{array}{c}\hfill \sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}{f}_{l\tau }\left(X(s_i,t_j)\right)-f_{l\tau }\left(X(s_0,t_0)\right)=0,\end{array}$$

(A3)

for all $l=0,1,...,L$ and $\tau =t_1,t_2,...,t_m$.

This impile that

$$\begin{array}{c}\hfill {\mathcal{F}}^T\mathsf{\Lambda }={\mathcal{F}}_0\end{array}$$

(A4)

or

$$\begin{array}{c}\hfill {\mathsf{\Lambda }}^T\mathcal{F}={\mathcal{F}}_0^T.\end{array}$$

(A5)

Appendix A.2. The Variance of the Estimation Error

The variance of the estimation error, denoted by ${\sigma }_{ϵ}^2$, is given by

$$\begin{array}{c}\hfill {\sigma }_{ϵ}^2=\mathrm{Var}\left[{\widehat{Z}}_{\mathrm{ST}-\mathrm{KED}}(s_0,t_0)-Z\left(s_0,t_0\right)\right]=\mathbb{E}\left[{\left({\widehat{Z}}_{\mathrm{ST}-\mathrm{KED}}(s_0,t_0)-Z\left(s_0,t_0\right)\right)}^2\right].\end{array}$$

(A6)

Applying Equations (1), (14), and (23), this expression can be rewritten as:

$$\begin{array}{cc}\hfill {\sigma }_{ϵ}^2& =\mathbb{E}\left[{\left(\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}\left({\mu }_A\left(s_i,t_j\right)+ϵ\left(s_i,t_j\right)\right)-\left({\mu }_A\left(s_0,t_0\right)+ϵ\left(s_0,t_0\right)\right)\right)}^2\right],\hfill \\ \hfill & =\mathbb{E}\left[{\left(\left(\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}{\mu }_A\left(s_i,t_j\right)-{\mu }_A\left(s_0,t_0\right)\right)+\left(\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}ϵ\left(s_i,t_j\right)-ϵ\left(s_0,t_0\right)\right)\right)}^2\right],\hfill \\ \hfill & =\mathbb{E}\left[{\left(\left(\sum _{l=0}^L\sum _{\tau =t_1}^{t_m}{\alpha }_{l\tau }\left(\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}{f}_{l\tau }\left(s_i,t_j\right)-f_{l\tau }\left(s_0,t_0\right)\right)\right)+\left(\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}ϵ\left(s_i,t_j\right)-ϵ\left(s_0,t_0\right)\right)\right)}^2\right].\hfill \end{array}$$

(A7)

Under the condition of unbiasedness established in Equation (A3), the variance can be expressed as follows:

$$\begin{array}{cc}\hfill {\sigma }_{ϵ}^2=& \mathbb{E}\left[{\left(\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}ϵ\left(s_i,t_j\right)-ϵ\left(s_0,t_0\right)\right)}^2\right],\hfill \\ \hfill =& \mathbb{E}\left[{\left(ϵ\left(s_0,t_0\right)\right)}^2-2\left(\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}ϵ\left(s_i,t_j\right)ϵ\left(s_0,t_0\right)\right)+{\left(\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}ϵ\left(s_i,t_j\right)\right)}^2\right],\hfill \\ \hfill =& \mathbb{E}\left[{\left(ϵ\left(s_0,t_0\right)\right)}^2\right]-\mathbb{E}\left[2\left(\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}ϵ\left(s_i,t_j\right)ϵ\left(s_0,t_0\right)\right)\right]+\mathbb{E}\left[{\left(\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}ϵ\left(s_i,t_j\right)\right)}^2\right],\hfill \\ \hfill =& \mathbb{E}\left[ϵ\left(s_0,t_0\right)ϵ\left(s_0,t_0\right)\right]-2\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}\mathbb{E}\left[ϵ\left(s_i,t_j\right)ϵ\left(s_0,t_0\right)\right]+\sum _{v=1}^m\sum _{u=1}^n{\lambda }_{uv}\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}\mathbb{E}\left[ϵ\left(s_i,t_j\right)ϵ\left(s_u,t_v\right)\right],\hfill \\ \hfill =& \mathrm{Cov}\left[ϵ\left(s_0,t_0\right),ϵ\left(s_0,t_0\right)\right]-2\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}\mathrm{Cov}\left[ϵ\left(s_i,t_j\right),ϵ\left(s_0,t_0\right)\right]+\sum _{v=1}^m\sum _{u=1}^n{\lambda }_{uv}\sum _{j=1}^m\sum _{i=1}^n{\lambda }_{ij}\mathrm{Cov}\left[ϵ\left(s_i,t_j\right),ϵ\left(s_u,t_v\right)\right].\hfill \end{array}$$

(A8)

The variance can be written as the following matrix:

$$\begin{array}{c}\hfill {\sigma }_{ϵ}^2=\mathbf{C}(\mathbf{0},0)-2{\mathsf{\Lambda }}^T{\mathsf{\Omega }}_0+{\mathsf{\Lambda }}^T\mathsf{\Omega }\mathsf{\Lambda }.\end{array}$$

(A9)

Appendix A.3. Minimization Approach for the ST-KED System Construction

The optimal weight vector $\mathsf{\Lambda }$ of the estimator ${\widehat{Z}}_{\mathrm{ST}-\mathrm{KED}}(s_0,t_0)$, as shown in Equation (15), is determined by solving the following minimization problem:

$$\begin{array}{c}\hfill \begin{array}{cc}\mathrm{minimum}\mathrm{of}& \mathbf{C}(\mathbf{0},0)-2{\mathsf{\Lambda }}^T{\mathsf{\Omega }}_0+{\mathsf{\Lambda }}^T\mathsf{\Omega }\mathsf{\Lambda }\\ \mathrm{subject}\mathrm{to}& {\mathsf{\Lambda }}^T\mathcal{F}={\mathcal{F}}_0^T.\end{array}\end{array}$$

(A10)

To solve this problem using the method of Lagrange multipliers [52], we formulate the Lagrangian function $\varphi$ as follows:

$$\begin{array}{c}\hfill \varphi (\mathsf{\Lambda },\mathcal{N}):=\mathbf{C}(\mathbf{0},0)-2{\mathsf{\Lambda }}^T{\mathsf{\Omega }}_0+{\mathsf{\Lambda }}^T\mathsf{\Omega }\mathsf{\Lambda }-2\left({\mathcal{F}}_0^T-{\mathsf{\Lambda }}^T\mathcal{F}\right)\mathcal{N}.\end{array}$$

(A11)

The first order partial derivatives with respect to $\mathsf{\Lambda }$ and $\mathcal{N}$ yield

$$\begin{array}{c}\hfill \frac{\partial \varphi }{\partial \mathsf{\Lambda }}(\mathsf{\Lambda },\mathcal{N})=-2{\mathsf{\Omega }}_0+2\mathsf{\Omega }\mathsf{\Lambda }+2\mathcal{F}\mathcal{N}\end{array}$$

(A12)

and

$$\begin{array}{c}\hfill \frac{\partial \varphi }{\partial \mathcal{N}}(\mathsf{\Lambda },\mathcal{N})=-2\left({\mathcal{F}}_0^T-{\mathsf{\Lambda }}^T\mathcal{F}\right).\end{array}$$

(A13)

Setting both partial derivative matrices to zero matrices, the ST-KED system for adaptive coefficient drift function can be rewritten in the form:

$$\begin{array}{c}\hfill \mathsf{\Omega }\mathsf{\Lambda }+\mathcal{F}\mathcal{N}={\mathsf{\Omega }}_0\end{array}$$

(A14)

and

$$\begin{array}{c}\hfill {\mathsf{\Lambda }}^T\mathcal{F}={\mathcal{F}}_0^T\end{array}$$

(A15)

with the result that

$$\begin{array}{c}\hfill {\mathcal{F}}^T\mathsf{\Lambda }={\mathcal{F}}_0.\end{array}$$

(A16)

These equations are represented in the form of a matrix equation, as seen in Equation (25).