Moving Horizon Estimation for ARMAX process with t-Distribution Noise

TL;DR

This paper introduces a robust moving horizon estimator for ARMAX processes with t-distribution noise, improving outlier resistance and computational efficiency over traditional methods.

Contribution

It proposes a novel recursive estimator based on influence functions for ARMAX models with t-distribution noise, enhancing robustness and speed.

Findings

Estimator has smaller variance than least-squares

More robust to outliers than traditional methods

Significantly faster than particle filters

Abstract

In this paper, instead of the usual Gaussian noise assumption, $t$-distribution noise is assumed. A Maximum Likelihood Estimator using the most recent N measurements is proposed for the Auto-Regressive-Moving-Average with eXogenous input (ARMAX) process with this assumption. The proposed estimator is robust to outliers because the `thick tail' of the t-distribution reduces the effect of large errors in the likelihood function. Instead of solving the resulting nonlinear estimator numerically, the Influence Function is used to formulate a computationally efficient recursive solution, which reduces to the traditional Moving Horizon Estimator when the noise is Gaussian. The formula for the variance of the estimate is derived. This formula shows explicitly how the variance of the estimate is affected by the number of measurements and noise variance. The simulation results show that the proposed estimator has smaller variance and is more robust to outliers than the Moving Window Least-Squares Estimator. For the same accuracy, the proposed estimator is an order of magnitude faster than the particle filter.

Tables

Table 1. Table 1: Models and Parameters Used in the Examples

Example	1	2	3
$N$	6,9,12	3
$e_k$	$t_3(0,0.5)$	$t_3(0,1)$
$A\left(q^{-1}\right)$	${\left(1+0.855q^{-1}\right)}^5$	$1-0.9q^{-1}$
$B\left(q^{-1}\right)$	$0.1q^{-1}$	$q^{-1}$
$C\left(q^{-1}\right)$	${\left(1+0.8q^{-1}\right)}^5$	$1-0.85q^{-1}$

Table 2. Table 2: Variance of y ^ k subscript ^ 𝑦 𝑘 \hat{y}_{{}_{k}} in Example 1

$k$			6	9	12	40	60
Proposed estimator Eq. (20)	Theoretical results from (24)	N
		$6$	0.6124	0.7055	0.7675	0.8617	0.8620
		$9$	-	0.6879	0.7499	0.8441	0.8444
		$12$	-	-	0.7237	0.8180	0.8183
	Simulation results	$6$	0.5946	0.6834	0.7416	0.8486	0.8729
		$9$	-	0.6723	0.7543	0.8338	0.8546
		$12$	-	-	0.7065	0.7881	0.8343
MWLSE Eq. (22)	Simulation results	$6$	0.9773	1.0705	1.1320	1.2270	1.2270
		$9$	-	1.0360	1.0980	1.1920	1.1920
		$12$	-	-	1.0460	1.1400	1.1400
ARMAX filter[7]	Simulation results	-	0.5946	0.6723	0.7065	0.4883	0.4889
Kalman filter	Simulation results	-	0.9773	1.0360	1.0460	0.4881	0.4889

Equations

$$A\left(q^{-1}\right)y_k=B\left(q^{-1}\right)u_k+C\left(q^{-1}\right)e_k,$$

$$A\left(q^{-1}\right)$$

$$B\left(q^{-1}\right)$$

$$C\left(q^{-1}\right)$$

$$f(e_k)=\frac{\mathbf{\Gamma }\left(\frac{\nu +1}{2}\right)}{\sqrt{\nu \pi }\sigma \mathbf{\Gamma }\left(\frac{\nu }{2}\right)}{\left(1+\frac{e_k^2}{{\sigma }^2\nu }\right)}^{-\frac{\nu +1}{2}}\triangleq t_{\nu }\left(0,\sigma \right)$$

$$x_{k+1}$$

$$y_k$$

$${\Phi }_a$$

$$\Omega$$

$${\hat{x}}_{{}_{\mathrm{MLE}}}=\arg \underset{{\hat{x}}_1}{\min }J,$$

$$J$$

$$e_k$$

$$\Phi$$

$$s_k$$

$${\sum }_{k=1}^T{\psi }_k(e_k)=0,$$

$$\mathrm{IF}\left(y;G,F\right)={\lim }_{h\to 0}\frac{G\left((1-h)F+h{\Delta }_y\right)-G(F)}{h},$$

$${\sum }_{k=1}^T{\int }_{-\infty }^{\infty }{\psi }_k\left(e_k\right)f\left(e_k\right)d{e}_k=0.$$

$${\sum }_{k=1}^T{\int }_{-\infty }^{\infty }{\psi }_k\left(e_k\right)\left[(1-h)f\left(e_k\right)+h{\Delta }_{e_k}\right]d{e}_k=0.$$

$${\hat{x}}_1(h)={\hat{x}}_1(0)+{\frac{d{\hat{x}}_1(h)}{dh}|}_{h=0}h+\frac{1}{2}{\frac{d^2{\hat{x}}_1(h)}{d{h}^2}|}_{h=0}{h}^2+\cdots$$

$$\mathrm{IF}(e)={\frac{d{\hat{x}}_1(h)}{dh}|}_{h=0}$$

$$\mathrm{IF}(e)$$

$$\text{mbox}\mathrm{IF}(e)$$

$${\hat{x}}_{{}_{\mathrm{MLE}}}\approx {\bar{x}}_1+\mathrm{IF}(e)\triangleq x_{{}_{\mathrm{IF}}}.$$

$$x_{{}_{\mathrm{IF}}}$$

$$\mathrm{IF}(e)$$

$$x_{{}_{\mathrm{IF}}}=\mathrm{IF}(e).$$

$$\mathrm{IF}(e)$$

$${\hat{x}}_{{}_T}$$

$${\hat{y}}_{{}_T}$$

$$\Psi$$

$${\mathcal{W}}_{T,N}$$

$$w_{{}_{T-N+i}}$$

$${\hat{x}}_{{}_{T+1}}$$

$${\hat{y}}_{{}_{T+1}}$$

$$z_k$$

$$w_k$$

Full text

††thanks: The authors acknowledge support by the Singapore National Research Foundation (NRF) under its Campus for Research Excellence And Technological Enterprise (CREATE) programme, specifically the Cambridge Centre for Advanced Research and Education in Singapore (Cambridge CARES, http://www.cares.cam.ac.uk), project C4T. This paper was not presented at any IFAC meeting. Corresponding author Dexiang Zhou. Tel. +65 8376 0221.

, , ,

1 Introduction

Gaussian noise is often assumed in state estimation. However, the Gaussian noise assumption is an approximation to reality. The occurrence of outliers, transient data in steady-state measurements, instrument failure, model nonlinearity, etc. can all induce non-Gaussian data [16].

It is well-known that the $\it t$-distribution has the property of ‘ thick tail’  to better model the occurence of outliers [10]. In addition, as a special case, the $\it t$-distribution reduces to the Gaussian distribution when its degree of freedom tends to infinity. Thus the $\it t$-distribution has the flexibility to characterize noise with Gaussian or non-Gaussian statistical properties. In the literature the heavy tail property of the $\it t$-distribution has been used to improve the performance of estimators in [13, 2, 1, 15]. These papers are based on the Bayesian method and the heavy tail property of $t$-distribution was used to reduce the effect of outliers. Our method differs from these estimators in that the probability density function (pdf) of the $\it t$-distribution is used directly in the likelihood function for a Maximum Likelihood Estimator (MLE).

The sensitivity of an estimator when the underlying noise assumption (such as Gaussian noise) is violated has been extensively studied in the robust statistics literature [9, 6, 8, 12]. In particular, Huber [9] studied the effect of outliers by contaminating the underlying noise distribution with data from an arbitrary unknown distribution. Hampel [6] proposed the Influence Function (IF) approach to describe the effect of an infinitesimal contamination of the underlying noise distribution. If the IF of an estimator is bounded and/or decreasing, or increasing slowly for large magnitude of noise, then the estimator is robust to outliers. It can be shown that the IF of the least-squares estimator increases linearly with the magnitude of noise and is unbounded [6]. This confirms the well-known fact that standard least-squares estimation is not robust to outliers. Other approaches for handling non-Gaussian noise include the particle filter [3] which is based on the Monte Carlo method, but at the expense of a heavier computational load.

IF has been mainly used as an analysis tool. The main contribution of this paper, Theorem 6, is to employ the IF for state estimation. We formulate a computationally efficient recursive algorithm that gives an approximate solution to the moving horizon MLE for Auto-Regressive-Moving-Average with eXogenous inputs (ARMAX) process with $t$-distribution noise.

ARMAX processes are popularly used in the statistical analysis of time series [4]. The Kalman filter for an ARMAX process is well-known, see Example 11.6 in [4]. In [7] a MLE of an ARMAX process with generalized $t$-distribution noise was developed. However, in [7] the IF was used to approximate the MLE by using all the measurements. This makes the estimator insensitive to plant parameter changes [11]. In this paper we derive the moving horizon version of the ARMAX filter. We derive a computationally efficient recursive algorithm in Section 3. We also analyze the statistical properties of the proposed estimator in Section 4. More specifically, the formula for the variance of the estimate is derived.

2 Maximum Likelihood Estimation and its approximation for ARMAX process with $t-$distribution noise

Consider the following single-input single-output ARMAX process:

[TABLE]

where $k$ is the sampling instance, $u_{k}$ and $y_{k}$ are the input and output respectively. The polynomials are assumed to be co-prime and given as

[TABLE]

where $n_{B}\leqslant n$ and $q^{-1}$ is the backward shift operator, i.e., $q^{-1}y_{k}=y_{k-1}$. The zeros of the polynomials $A(q^{-1})$ and $C(q^{-1})$ are inside the unit disc. $e_{k}$ is an independent and identically distributed random variable associated with the zero mean $t$-distribution pdf [10]

[TABLE]

The degree of freedom and scale parameters are given by $\nu$ and $\sigma$ respectively, and ${\bf\Gamma}$ is the gamma function.

For convenience, we will derive the estimator in state-space form. The ARMAX process (1) in minimal state-space form is given as [4]

[TABLE]

where

[TABLE]

2.1 Maximum Likelihood Estimation

Given measurements $y_{1},~{}y_{2},~{}\cdots,~{}y_{{}_{T}}$, we denote the MLE of $x_{1}$ for the ARMAX process (3) with $t$-distribution noise as

[TABLE]

where

[TABLE]

Then the solution to (4) can be found by solving the equation

[TABLE]

where $\psi_{k}(e_{k})\triangleq\frac{d}{d\hat{x}_{1}}\left(-\ln f\left(e_{k}\right)\right)=-\frac{\left(\nu+1\right)\left(H\Phi^{k-1}\right)^{{}^{\prime}}e_{k}}{\sigma^{2}\nu+(e_{k})^{2}}$.

Definition 1**.**

Influence Function [6].

The IF of an estimator $G$ at distribution $F$ is

[TABLE]

where $h\in[0,1]$ and $\Delta_{y}$ denotes the probability measure which puts mass 1 at the point $y$.

The IF describes the effect of an infinitesimal contamination at the point $y$ on the estimate, standardized by the mass of the contamination. In this paper the noise is assumed to be $t$ distributed with pdf $f(e_{k})$. Taking the expectation of (6) gives

[TABLE]

Now we wish to obtain the IF of our estimator. As a first step, we replace $f(e_{k})$ in (7) with $(1-h)f\left(e_{k}\right)+h\Delta_{e_{k}}$

[TABLE]

Equation (8) implicitly defines our estimator as a function of $h$ and we denote it as $\hat{x}_{1}(h)$. In general, $\hat{x}_{1}(h)$ is a nonlinear state estimator. The Taylor series expansion of $\hat{x}_{1}(h)$ at $h=0$ is given by

[TABLE]

Definition 1 suggests that the IF of our estimator can be obtained by differentiating equation (8) wrt. $h$ and setting $h=0$

[TABLE]

where we denote $\hat{x}_{1}(0)$ as $\bar{x}_{1}$.

The derivation of (10) is given in Appendix A.

2.2 IF Approximation

Lemma 2**.**

The IF of MLE for ARMAX process (3) with $t$-distribution noise (2) is

[TABLE]

When the noise $e_{k}$ is Gaussian, (11) becomes

[TABLE]

{pf}

See Appendix B.

Since (8) reduces to (6) when $h=1$, $\hat{x}_{1}(1)$ is the solution that we are looking for, namely the solution to (6). We can get an approximate solution to this by retaining only the first two terms of (9). Using (10),

[TABLE]

Remark 3**.**

When the noise is Gaussian, substituting (12) into (13) gives

[TABLE]

which is the least-squares estimator. It is re-assuring that our estimator reduces to the well-known result when the noise is Gaussian. Note that $x_{{}_{\rm IF}}=\hat{x}_{{}_{\rm MLE}}$ in this case since under the Gaussian noise assumption, the MLE and least-squares estimator coincide.

In the sequel we shall assume that $\bar{x}_{1}=0$, for simplicity. Then, since $\left.e_{k}\right|_{\hat{x}_{1}=\bar{x}_{1}}=y_{k}-Hs_{k}$, (11) becomes

[TABLE]

and (13) becomes

[TABLE]

Remark 4**.**

*We consider only cases when $C\left(q^{-1}\right)$ is asymptotically stable. In this case, $\bar{x}_{1}=0$ is a simplifying assumption since error due to a wrong initial state estimate will reduce to zero after an initial transient. *

3 The Proposed Estimator

Instead of using all the measurements, we adopt the moving horizon principle and use the most recent $N$ measurements $\left\{y_{k}\right\}_{T-N+1}^{T}$.

Therefore, the IF of MHE-TD can be obtained by replacing $k=1$ with $k=T-N+1$ in (15)

[TABLE]

We denote the proposed estimator as Moving Horizon Estimation of ARMAX process with $t$-distribution noise (MHE-TD). By recursively applying the ARMAX model, we can obtain an expression for the estimate of the current state $x_{{}_{T}}$, and this is given in the next theorem.

Theorem 5**.**

IF approximation of Batch MHE-TD**

Given the ARMAX process (1) where $C(q^{-1})$ is asymptotically stable and noise $e_{k}$ with $\it t$-distribution pdf of (2), the IF approximation of MHE-TD is

[TABLE]

where $s_{{}_{T}}$ is given by (5d) and

[TABLE]

{pf}

See Appendix C.

In many cases the observations are obtained sequentially. It is desirable to derive a recursive solution for the IF approximation of MHE-TD as follows.

Theorem 6**.**

IF Approximation of Recursive MHE-TD*    The IF approximation of MHE-TD satisfies the following recursive algorithm*

[TABLE]

where

[TABLE]

{pf}

See Appendix D.

Remark 7**.**

When the noise $e_{k}$ in ARMAX process (1) is Gaussian (i.e., $\nu=\infty$), $w_{k}$ in (21b) reduces to $y_{k}-Hs_{k}$. Then $z_{k}$ reduces to $y_{k}$. Thus the recursive MHE-TD reduces to the following recursive solution of the standard Moving Window Least-Squares Estimator (MWLSE).

[TABLE]

Comparing (20a) and (22a) the recursive MHE-TD has the same structure as the standard MWLSE except for a nonlinear transformation, i.e., (21b). When $\Omega=0,~{}\Gamma=0$, Equation (22a) reduces to (10) of [11].

The MHE uses only a finite number of past measurement samples; the oldest measurement is discarded as a new sample becomes available. The ARMAX filter [7] uses the so-called growing memory strategy where the past measurement is not discarded as new measurement becomes available. Therefore, the ARMAX filter can be obtained if the last term of (20a) is discarded and $L_{N}$ and $P_{N}$ become

[TABLE]

4 Variance of the Proposed Estimator

Let the actual noise have pdf $g(e)$ which is different from $f(e)$ used in the design of MHE-TD. It is well-known that the mean and variance of an estimator can be analysed by its IF [6, 14, 5]. The statistical properties of the IF approximation of MHE-TD are derived by using IF as follows.

Theorem 8**.**

Variance of IF Approximation of MHE-TD *    Consider the estimator (18). Assume that the data are generated from the ARMAX process (1) where the $t$-distribution pdf of $e_{k}$ is $g(e)$. The variance of the estimate is given by*

[TABLE]

where

[TABLE]

With $\hat{y}_{{}_{T}}=H\hat{x}_{{}_{T}}$,

[TABLE]

{pf}

See Appendix E.

The third term of (24) which depends on $T$ is due to the dynamic of the ARAMX model. For a fixed $T$, when $N$ increases $M_{N}$ decreases. Thus the variance decreases with increasing $N$. This is also shown in Table 2.

Corollary 9**.**

The expectation of Moving Window Least-Squares Estimation is given by

[TABLE]

and its variance is

[TABLE]

where

[TABLE]

With $\hat{y}_{{}_{T}}=H\hat{x}_{{}_{T}}$,

[TABLE]

{pf}

By setting $\nu=\infty$ in (23).

The assumption that $g(e)$ is the same at each sampling instant $k$ is commonly made. In the following we extend to the case where $g(e)$ could be different at different $k$ denoted by $g_{k}(e)$. This is useful for the analysis of outliers (see Example 2).

Theorem 10**.**

Property of IF approximation of MHE-TD with outlier*    Consider the estimator (18). Assume that the data are generated from the ARMAX process (1) where $x_{1}=0$ and the noise distribution contains an outlier at $k_{1}$, i.e.,*

[TABLE]

then the expectation of the estimator (18) is

[TABLE]

where

[TABLE]

{pf}

See Appendix F.

Remark 11**.**

The expectation of Moving Window Least-Squares Estimation with outlier $e_{k_{1}}$ can be obtained by setting $\nu=\infty$ in (29)

[TABLE]

where $E\left(s_{{}_{T}}\right)$ is given in (10).

5 Examples

The models and parameters used in the examples are summarized in Table 1. The statistical results in Example 1,2 and 3 are obtained from 1000 Monte Carlo simulation runs.

5.1 Example 1:Estimation Variance

Table 2 shows the variances from 1000 simulation runs. The proposed estimator (Equation (20)), ARMAX filter [7], Kalman filter and MWLSE (Equation (22)) are compared for different horizon length, $N$, and different instance, $k$. An example used to illustrate the calculation can be found in Appendix G.

Firstly, it can be seen that for the proposed estimator the theoretical and simulation variance matched approximately. Secondly, when $N$ increases, the variance decreases and therefore $N$ can be used as a tuning parameter for estimator performance. Thirdly, the variance of the proposed estimator is less than that of the MWLSE. Fourthly, for small $k$ ($=6,9,12$) our proposed estimator gives smaller variances than the Kalman Filter because of more accurate noise modeling. However, for large $k$ ($=40,60$), the estimators with growing memory (ARMAX filter and Kalman filter) give smaller estimate variance than the estimators with fixed memory (the proposed estimator and MWLSE). Finally, since the proposed estimator is the moving horizon version of the ARMAX filter, the results for both are the same when $N=k$. Likewise, the MWLSE is the moving horizon version of the Kalman filter, the results for both are also the same for $N=k$.

Figure 1 shows that the estimation error of the proposed estimator is smaller than that of MWLSE.

5.2 Example 2: Performance with outlier

An outlier of $-10$ is introduced at $k=30$.

Figure 2 shows 1000 simulation runs of the proposed estimator (Equation (20)) and MWLSE (Equation (22)). The yellow curves are the mean values of the 1000 runs. The proposed estimator is less affected by the outlier than the MWLSE. These yellow curves in the top and bottom figures can be calculated from (29) and (31) respectively.

5.3 Example 3: Comparison with particle filter

The proposed estimator is compared with a particle filter (bootstrap implementation)[3]. Taking the numerical solution of MLE (4) as the ground truth, we compared the estimators using the index

[TABLE]

where, for the $j$-th run, $\hat{x}_{k}^{(MLE,j)}$ is the numerical solution of the MLE (4) and $\hat{x}_{k}^{(j)}$ is the estimate of $x_{k}$. The indexes, at $k=50$, are $0.0082$, $0.0142$ and $0.0016$ for the proposed estimator, and particle filters with 100 and 1000 particles respectively. The computational times per run using MATLAB R2015a with i7-5500U processor @2.4GHz, 16 GB RAM, are of the order of $0.1$ millisecond, one millisecond and ten milliseconds respectively. For about the same performances ($0.0082$ and $0.0142$) as given by index $I$, the computational time of the proposed estimator is an order of magnitude faster than the particle filter.

6 Conclusion

The IF is employed to give an approximate solution to the moving horizon MLE for ARMAX process with $t$-distribution noise. The approximate solution can be formulated as a recursive MHE-TD algorithm which makes it suitable for on-line and real-time implementation at high sampling rates. We also used the IF to derive a formula for the estimates. The examples show that the proposed estimator gives an estimate with a smaller variance than the MWLSE. It is also less affected by outliers than the MWLSE. For about the same performance index, the computational time of the proposed estimator is an order of magnitude faster than the particle filter.

Appendix A Appendix: Derivation of (10)

Differentiating Equation (8) wrt. $h$ gives

[TABLE]

Substituting (7) and

[TABLE]

into (32) gives

[TABLE]

Setting $h=0$ in (33) gives

[TABLE]

We denote $\hat{x}_{1}(0)$ as $\bar{x}_{1}$.

Thus we have

[TABLE]

Replacing the integration variable $e_{k}$ by $e$ in (34) gives

[TABLE]

which is (10).

Appendix B Appendix: Proof of Lemma 2

Substituting (2) into (5a)

[TABLE]

Differentiating $-{\rm ln}f\left(e_{k}\right)$ with respect to $\hat{x}_{1}$ gives

[TABLE]

Differentiating $\psi_{k}\left(e\right)$ with respect to $\hat{x}_{1}$ gives

[TABLE]

Substituting (35) and (36) into (10) gives

[TABLE]

which is (11).

When $\nu=\infty$ we have

[TABLE]

Substituting (37a) and (37b) into (11) and setting e_{k}=y_{k}-H\Phi^{k-1}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bar{x}_{1}}-H{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}_{k} give

[TABLE]

which is (12).

Appendix C Appendix: Proof of Theorem 5

Rewriting the state space model (3) of ARMAX as

[TABLE]

where $\Phi=\Phi_{a}-\Omega H$. Then using (16) the current state estimate $x_{{}_{T}}$ is obtained by recursively applying the model (38)

[TABLE]

where $s_{{}_{T}}$ is given by (5d).

Multiplying (17) by $\Phi^{T-1}$ and simplifying give

[TABLE]

Rewriting (40) gives

[TABLE]

where

[TABLE]

Equation (18a) is obtained by substituting (41) into (39).

Appendix D Appendix: Proof of Theorem 6

Define

[TABLE]

Thus the estimator (18a) can be written as

[TABLE]

Notice that Equation (42) gives $\Delta x_{{}_{T}}$ from the minimization of the Moving Window Least-Squares loss function

[TABLE]

The recursive solution of the above Moving Window Least-Squares [11] is

[TABLE]

where

[TABLE]

In the next time instant

[TABLE]

Expressing (5d) in recursive form we have

[TABLE]

and

[TABLE]

Define z_{k}=w_{k}+H{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}_{k}. Substituting w_{{}_{T+1}}=z_{{}_{T+1}}-H{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}_{{}_{T+1}} and w_{{}_{T-N+1}}=z_{{}_{T-N+1}}-H{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}_{{}_{T-N+1}} into (43) gives

[TABLE]

Substituting (44) and (45) into (LABEL:xhatTT1) gives

[TABLE]

which is Equation (20a), and (20b) follows.

Appendix E Appendix: Proof of Theorem 8

By substituting $y_{k}=Hx_{k}+e_{k}$ into (5d) and assuming $x_{1}=0$ for simplicity, $s_{{}_{T}}$ can be rewritten as

[TABLE]

Substituting (47) into (39) gives

[TABLE]

Taking expectation of (48), recognizing that $e_{k}$ is zero mean and $\textstyle{\int_{-\infty}^{\infty}}{\rm{IF}}(e)g(e)de=0$ from (17), we get

[TABLE]

Therefore the variance of $\hat{x}_{{}_{T}}$ is given by

[TABLE]

Substituting (17) into (LABEL:varyhatapp1), noticing that $e_{k}$ is i.i.d and simplifying, give

[TABLE]

where

[TABLE]

which is (23), and (24) follows.

Appendix F Appendix: Proof of Theorem 10

Recall (47) and $x_{1}=0$

[TABLE]

Then

[TABLE]

which is (10). Let $k_{2}=k_{1}+N-1$.

Substituting (50) into (40) gives

[TABLE]

From the estimator (39) we have

[TABLE]

Using (52) and (18b) gives

[TABLE]

Substituting (51) and (28) into the second term of (53) and simplifying give

[TABLE]

With

[TABLE]

equation (LABEL:intIF2) becomes

[TABLE]

Equation (29) is obtained by substituting (LABEL:intIF3) into (53).

Appendix G Appendix: Numerical Illustration Example

For example, for $N=6$ and $k=9$, substituting the parameters for Example 1 into (24) gives

[TABLE]

and

[TABLE]

$\rho_{1}=2.6667,~{}\rho_{2}=1,~{}\rho_{3}=0.7414,~{}\rho_{4}=2.6667$. Substituting the parameters for Example 1 into (LABEL:varyhatgaussian) gives

[TABLE]