Optimal designs for series estimation in nonparametric regression with correlated data

TL;DR

This paper develops optimal experimental designs for series estimators in nonparametric regression with correlated data, improving estimator performance through explicit solutions and optimal design point selection.

Contribution

It provides explicit solutions for the best linear oracle estimator in continuous time for Markovian error processes and introduces data-driven estimators with optimized design points.

Findings

New series estimator outperforms traditional methods

Optimal design points enhance estimator accuracy

Simulation confirms improved performance

Abstract

In this paper we investigate the problem of designing experiments for series estimators in nonparametric regression models with correlated observations. We use projection based estimators to derive an explicit solution of the best linear oracle estimator in the continuous time model for all Markovian-type error processes. These solutions are then used to construct estimators, which can be calculated from the available data along with their corresponding optimal design points. Our results are illustrated by means of a simulation study, which demonstrates that the new series estimator has a better performance than the commonly used techniques based on the optimal linear unbiased estimators. Moreover, we show that the performance of the estimators proposed in this paper can be further improved by choosing the design points appropriately.

Tables

Table 1. Table 1: Optimal time points for series estimation minimizing the criterion ( 4.22 ). The covariance kernel is given by exp ⁡ ( − L ​ | s − t | ) 𝐿 𝑠 𝑡 \exp(-L|s-t|) with L = 1 𝐿 1 L=1 (first row) and L = 5 𝐿 5 L=5 (second row).

$L$	$n=4$	$n=7$
$1$	$0.00,\mathrm{\hspace{0.17em}0.25},\mathrm{\hspace{0.17em}0.52},\mathrm{\hspace{0.17em}1.00}$	$0.00,\mathrm{\hspace{0.17em}0.12},\mathrm{\hspace{0.17em}0.27},\mathrm{\hspace{0.17em}0.45},\mathrm{\hspace{0.17em}0.57},\mathrm{\hspace{0.17em}0.77},\mathrm{\hspace{0.17em}1.00}$
$5$	$0.00,\mathrm{\hspace{0.17em}0.25},\mathrm{\hspace{0.17em}0.51},\mathrm{\hspace{0.17em}1.00}$	$0.00,\mathrm{\hspace{0.17em}0.12},\mathrm{\hspace{0.17em}0.27},\mathrm{\hspace{0.17em}0.45},\mathrm{\hspace{0.17em}0.57},\mathrm{\hspace{0.17em}0.76},\mathrm{\hspace{0.17em}1.00}$

Table 2. Table 2: Simulated mean integrated squared error of the estimators f ^ ( J ) , n superscript ^ 𝑓 𝐽 𝑛 \hat{f}^{(J),n} and f ˇ ( J ) , n superscript ˇ 𝑓 𝐽 𝑛 \check{f}^{(J),n} defined in ( 5.1 ) and ( 5.2 ) for different regression functions. The covariance kernel is given by exp ⁡ ( − | s − t | ) 𝑠 𝑡 \exp(-|s-t|) . Third column: optimal design; Fourth column: comparative design in ( 5.3 ). Left part: n = 4 𝑛 4 n=4 observations; right part: n = 7 𝑛 7 n=7 observations.

		$n=4$		$n=7$
		design		design
$f$	estimator	optimal	(5.3)	optimal	(5.4)
(5.5)	${\hat{f}}^{\left(J\right),n}$	1.72	2.06	1.58	1.59
	${\check{f}}^{\left(J\right),n}$	1.89	2.22	1.76	1.77
(5.6)	${\hat{f}}^{\left(J\right),n}$	1.67	2.04	1.54	1.56
	${\check{f}}^{\left(J\right),n}$	1.89	2.21	1.76	1.79

Table 3. Table 3: Simulated mean integrated squared error of the estimators f ^ ( J ) , n superscript ^ 𝑓 𝐽 𝑛 \hat{f}^{(J),n} and f ˇ ( J ) , n superscript ˇ 𝑓 𝐽 𝑛 \check{f}^{(J),n} defined in ( 5.1 ) and ( 5.2 ) for different regression functions. The covariance kernel is given by exp ⁡ ( − 5 ​ | s − t | ) 5 𝑠 𝑡 \exp(-5|s-t|) . Third column: optimal design; Fourth column: comparative design in ( 5.3 ). Left part: n = 4 𝑛 4 n=4 observations; right part: n = 7 𝑛 7 n=7 observations.

		$n=4$		$n=7$
		design		design
$f$	estimator	optimal	(5.3)	optimal	(5.4)
(5.5)	${\hat{f}}^{\left(J\right),n}$	0.65	2.13	0.47	0.51
	${\check{f}}^{\left(J\right),n}$	0.77	2.30	0.58	0.62
(5.6)	${\hat{f}}^{\left(J\right),n}$	0.64	2.09	0.43	0.43
	${\check{f}}^{\left(J\right),n}$	0.81	2.30	0.59	0.59

Table 4. Table 4: Simulated mean integrated squared error of the estimators f ^ ( J ) , n superscript ^ 𝑓 𝐽 𝑛 \hat{f}^{(J),n} and f ˇ ( J ) , n superscript ˇ 𝑓 𝐽 𝑛 \check{f}^{(J),n} defined in ( 5.1 ) and ( 5.2 ) for different regression functions. The error process is a Brownian motion. Third column: optimal design; Fourth column: comparative design in ( 5.3 ). Left part: n = 4 𝑛 4 n=4 observations; right part: n = 7 𝑛 7 n=7 observations.

		$n=4$		$n=7$
		design		design
$f$	estimator	optimal	(5.3)	optimal	(5.4)
(5.5)	${\hat{f}}^{\left(J\right),n}$	0.16	0.41	0.13	0.14
	${\check{f}}^{\left(J\right),n}$	0.15	0.43	0.12	0.12
(5.6)	${\hat{f}}^{\left(J\right),n}$	0.13	0.45	0.11	0.11
	${\check{f}}^{\left(J\right),n}$	0.15	0.48	0.12	0.13

Equations

$$Y_i=f(X_i)+{\epsilon }_i,i=1,\dots ,n,$$

$$Y_{t_i}=f(t_i)+{\epsilon }_{t_i},i=1,\dots ,n,$$

$$L^{2}([0,1])=\Big{\{}g:[0,1]\rightarrow\mathbb{R}:\int_{0}^{1}g^{2}(t)dt<\infty\Big{\}},$$

$$f(t)=\sum _{j\in \mathbb{N}}{\theta }_j{\phi }_j(t),$$

$${\theta }_j=⟨f,{\phi }_j⟩={\int }_0^{1}f(t){\phi }_j(t)dtj\in \mathbb{N}.$$

$${\hat{f}}^{(J)}(t)=\sum _{j=1}^J{\hat{\theta }}_j{\phi }_j(t),$$

$$\sum _{i=2}^n(t_i-t_{i-1}){\phi }_j(t_{i-1})Y_{t_{i-1}}.$$

$$\mathbb{E}\Big{[}\int_{0}^{1}\big{(}\hat{f}^{(J)}(t)-f(t)\big{)}^{2}\Big{]}=\sum_{j=1}^{J}\mathbb{E}\big{[}(\hat{\theta}_{j}-\theta_{j})^{2}\big{]}+\sum_{j=J+1}^{\infty}\theta_{j}^{2},$$

$$\displaystyle\sum_{j=1}^{J}\mathbb{E}\big{[}(\hat{\theta}_{j}-\theta_{j})^{2}\big{]}$$

$$Y_t=f(t)+{\epsilon }_t,t\in [0,1],$$

$$d{Y}_t=\dot{f}(t)dt+d{\epsilon }_t,t\in [0,1],$$

$$\mathbb{E}[{\epsilon }_s{\epsilon }_t]=K(s,t)=u(s)v(t)\text{mbox}fors\le t,$$

$$q(t)=\frac{u(t)}{v(t)},$$

$${\hat{\theta }}_j={\int }_0^1{Y}_t{\xi }_j(dt),j\in \mathbb{N},$$

$$\sum_{j=1}^{\infty}\big{\{}(\mathbb{E}[\hat{\theta}_{j}])^{2}+\mbox{Var}(\hat{\theta}_{j})\big{\}}=\sum_{j=1}^{\infty}\Big{(}\int_{0}^{1}f(t)d\xi_{j}(t)\Big{)}^{2}+\sum_{j=1}^{\infty}\int^{1}_{0}\int^{1}_{0}K(s,t)d\xi_{j}(s)d\xi_{j}(t)<\infty.$$

$$\hat{f}(t)=\sum _{j=1}^{\infty }{\hat{\theta }}_j{\phi }_j(t).$$

$$\displaystyle\lim_{J\rightarrow\infty}\mathbb{E}\Big{[}\int_{0}^{1}\big{(}\hat{f}^{(J)}(t)-f(t)\big{)}^{2}dt\Big{]}=\lim_{J\rightarrow\infty}\sum_{j=1}^{J}\mathbb{E}[(\hat{\theta}_{j}-\theta_{j})^{2}]=\sum_{j=1}^{\infty}\mathbb{E}[(\hat{\theta}_{j}-\theta_{j})^{2}]<\infty,$$

$$\text{mbox}MISE(\hat{f})$$

$$\begin{split}\Psi_{j}(\xi_{j}):&=\mathbb{E}\Big{[}\Big{(}\int_{0}^{1}Y_{t}\xi_{j}(dt)-\theta_{j}\Big{)}^{2}\Big{]}\\ &=\int_{0}^{1}\int_{0}^{1}\big{[}f(s)f(t)+K(s,t)\big{]}\xi_{j}(ds)\xi_{j}(dt)-2\theta_{j}\int_{0}^{1}f(s)\xi_{j}(ds)+\theta^{2}_{j}\\ &=\int_{0}^{1}\int_{0}^{1}K(s,t)\xi_{j}(ds)\xi_{j}(dt)+\Big{(}\int_{0}^{1}f(s)\xi_{j}(ds)-\theta_{j}\Big{)}^{2}\,.\end{split}$$

$${\tilde{\theta }}_j={\int }_0^1{Y}_t{\phi }_j(t)dt(j\in \mathbb{N}),$$

$$\sum _{j=1}^{\infty }{\int }_0^1{\int }_0^{1}K(s,t)d{\xi }_j(s)d{\xi }_j(t)<\infty .$$

$$\sum _{j=1}^{\infty }{\int }_0^1{\int }_0^{1}K(s,t){\phi }_j(s){\phi }_j(t)dsdt=\sum _{j=1}^{\infty }{\lambda }_i{\int }_0^1{\phi }_j(t){\phi }_j(t)dt=\sum _{j=1}^{\infty }{\lambda }_j<\infty .$$

$${\xi }_j^{\ast }(dt)=\frac{{\theta }_j}{1+c}\left(P_0{\delta }_0(dt)+P_1{\delta }_1(dt)+p(t)dt\right),$$

$$c$$

$$P_0$$

$$P_1$$

$$p(t)$$

$$c$$

$$P_0$$

$$P_0=\frac{1}{f(0)}.$$

$$f^{\ast }(t)=\sum _{j=1}^{\infty }{\hat{\theta }}_j^{\ast }{\phi }_j(t),$$

$${\hat{\theta }}_j^{\ast }={\int }_0^1{Y}_t{\xi }_j^{\ast }(dt),j\in \mathbb{N},$$

Full text

Keywords: Optimal design, nonparametric regression, integrated mean squared error, optimal estimator

AMS Subject classification: 62K05

1 Introduction

Nonparametric regression is a common tool of statistical inference with numerous applications [see the monographs of Fan and Gijbels, (1996), Efromovich, (1999), Fan and Yao, (2003), Tsybakov, (2009) among many others]. The basic model is formulated in the form

[TABLE]

where one usually distinguishes between random and fixed predictors $X_{i}$. In the latter case a natural question is how to choose $X_{1},\ldots,X_{n}$ to obtain the most precise estimates of the regression function $f$ and several authors have worked on this problem. For example, Müller, (1984), Biedermann and Dette, (2001) and Zhao and Yao, (2012) derived optimal designs with respect to different criteria for kernel estimates, while Dette and Wiens, 2008a and Dette and Wiens, 2008b considered the design problem for series estimation in terms of spherical harmonics and Zernike polynomials, respectively. We also refer to the work of Efromovich, (2008), who proposed a sequential allocation scheme in a nonparametric model of the form (1.1) with random predictors and heteroscedastic errors. A common feature of the literature in this field is the fact that all authors investigate the design problem in a model (1.1) with independent errors. However, there are many situations, where this assumption is not satisfied, in particular, when the explanatory variable represents time.

The reason for this gap in the existing literature is that the design problem for models with correlated errors (even parametric models) is substantially harder compared to the uncorrelated case. In contrast to the latter case, where a very well developed and powerful methodology for the construction of optimal designs has been established [see, for example, the monograph of Pukelsheim, (2006)], optimal designs for models with correlated observations are only available in rare circumstances considering parametric models [see, for example, Pázman and Müller, (2001), Näther and Simák, (2003), Müller and Stehlík, (2004), Dette et al., (2009), Zhigljavsky et al., (2010); Pázman, (2010), Harman and Stulajter, (2010), Amo-Salas et al., (2012), Stehlik et al., (2015), Rodríguez-Díaz, (2017) among others]. Some general results on optimal designs for linear models with correlated observations can be found in the seminal work of Sacks and Ylvisaker, (1966, 1968), while more recently in a series of papers Dette et al., (2013, 2016, 2017) provided a general approach for the problem of designing experiments in linear models with correlated observations by considering the problem of optimal (unbiased linear) estimation and optimal design simultaneously. Usually, authors use asymptotic arguments to embed the discrete (non-convex) optimization problem in a continuous (or approximate) one. However, unlike the uncorrelated case, in the context of correlated observations this approach does not simplify the problem substantially and due to the lack of convexity the resulting approximate optimal design problems for regression models with correlated observations are still extremely difficult to solve.

In this paper we consider optimal design theory for series estimation in the nonparametric regression model (1.1) with correlated data. The basic notation and the general design problem are introduced in Section 2. In order to address the particular difficulties in design problems for series estimation from correlated data, in Section 3 we consider a continuous time version of the discrete model. We first determine optimal oracle estimators for the coefficients in a Fourier expansion of the regression function $f$. These are shrinkage estimators and not unbiased.

Section 4 is devoted to the implementation of the results from Section 3 for the construction of an efficient estimator with a corresponding optimal design. In particular, we determine an optimal approximation of the Fourier coefficients in the continuous model (which requires the full trajectory of the process) by an estimator which can be calculated from the available data $\{Y_{t_{1}},\ldots,Y_{t_{n}}\}$ and determine the designs points $t_{1},\ldots,t_{n}$ such that the approximation has minimal mean squared error with respect to the solution in the continuous time model. The resulting estimator is a two stage estimator shrinking the best linear unbiased estimator when the design points are chosen in an optimal way. The superiority of our approach is demonstrated in Section 5 by means of a small simulation study, while all technical details are given in Section 6.

2 Optimal designs for series estimation

Throughout this paper we consider the nonparametric regression model with a fixed design, that is,

[TABLE]

where $f:[0,1]\to\mathbb{R}$ is the regression function, $0\leq t_{1}<t_{2}<\ldots<t_{n}\leq 1$ are $n$ distinct time points in the interval $[0,1]$, $\mathbb{E}[\varepsilon(t_{j})]=0$ and $K(t_{i},t_{j})=\mathbb{E}[\varepsilon_{t_{i}}\varepsilon_{t_{j}}]$ denotes the covariance between observations at the points $t_{i}$ and $t_{j}$ ($i,j=1,\dots,n$). Let

[TABLE]

denote the space of square integrable (real valued) functions with inner product $\langle g_{1},g_{2}\rangle=\int_{0}^{1}g_{1}(t)g_{2}(t)dt$ and norm \|g\|_{2}=\big{(}\int_{0}^{1}g^{2}(t)dt\big{)}^{1/2}. Let $\{\varphi_{j}(\cdot):j\in\mathbb{N}\}$ be an orthonormal basis, then any function $f\in L^{2}([0,1])$ admits a series expansion of the form

[TABLE]

in $L^{2}{([0,1])}$ with Fourier coefficients

[TABLE]

Moreover, the coefficients are squared summable, that is, $\sum_{j\in\mathbb{N}}\theta_{j}^{2}<\infty$. In order to estimate the unknown function $f$ we now follow the idea of projection estimators [see Tsybakov, (2009), pp.47] and estimate the truncated series ${f}^{(J)}(t)=\sum_{j=1}^{J}\theta_{j}\varphi_{j}(t)$ by

[TABLE]

where $\hat{\theta}_{j}$ is an appropriate estimator for the Fourier coefficient $\theta_{j}$ $(j=1,\ldots,J$). For example, if $\max^{n}_{i=2}(t_{i}-t_{i-1})\to 0$, as $n\to\infty$, an asymptotically unbiased estimator of $\theta_{j}$ is given by

[TABLE]

More general estimators will be specified later on. At this point it is only important to note that the performance of any reasonable estimator will depend on the design points $t_{1},\ldots,t_{n}$. We are interested in choosing these design points such that the mean integrated squared error

[TABLE]

is minimal. We also note that any solution of this discrete optimization problem depends on the unknown regression function $f$, the truncation point $J$ used in (2.4) and on the covariance kernel $K$, which is assumed to be known throughout this paper. On the other hand, the term $\sum_{j=J+1}^{\infty}\theta_{j}^{2}$ does not depend on the design points which can therefore, be determined by minimizing \sum_{j=1}^{J}\mathbb{E}\big{[}(\hat{\theta}_{j}-\theta_{j})^{2}\big{]} with respect to the choice of $t_{1},\ldots,t_{n}$. For example, if $\hat{\theta}_{j}=\sum_{i=1}^{\ell_{j}}\alpha_{ji}Y_{t_{i}}$ is a linear estimator of $\theta_{j}$ ($j=1,\ldots,J$) we have that

[TABLE]

which has to be minimized with respect to the choice of the time points $t_{1},\ldots,t_{n}$.

3 Optimal estimation in the continuous time model

The discrete optimization problem (2.6) stated in the previous section is extremely difficult to be solved. In this section in order to derive efficient designs, we investigate a simpler problem and consider the continuous time nonparametric regression model of the form

[TABLE]

where $f$ is an unknown square integrable function and the error process $\varepsilon=\{\varepsilon_{t}:t\in[0,1]\}$ is a centered Gaussian process with covariance kernel $K(s,t)=\mathbb{E}[\varepsilon_{s}\varepsilon_{t}]$. As we assume that the full trajectory of the process is available, there is in fact no optimal design problem but only the issue of optimal estimation of the regression function $f$. The optimal design question will appear later, when we return to the discrete model (2.1). The main result of this section provides an oracle solution of the optimal estimation problem. In particular, the optimal estimator depends on the unknown function $f$ in model (3.1) and is therefore, not implementable (even if the full trajectory of the process $\{Y_{t}:t\in[0,1]\}$ is available). However, our solution serves as benchmark and actually provides a clear hint how good estimators and corresponding optimal designs can be constructed. This will be formulated precisely in Section 4.

Model (3.1) is often written in terms of a stochastic differential equation (provided that the regression function $f$ is differentiable with derivative $\dot{f}$), that is

[TABLE]

If $\varepsilon=\{\varepsilon_{t}:t\in[0,1]\}$ is a Brownian motion, the model (3.2) is called Gaussian white noise model and has found much attention in the statistical literature [see, for example, Ibragimov and Hasminskiĭ, (1981) or Tsybakov, (2009) among many others]. In particular, the model is asymptotically equivalent to the nonparametric regression model $Z_{i}=\dot{f}(i/n)+\eta_{i}$ ($i=1,\ldots,n$), where $\eta_{1},\ldots,\eta_{n}$ are independent standard normally distributed random variables [see Brown and Low, (1996)]. Note that the focus in the aforementioned publications is on the optimal estimation of the function $\dot{f}$, whereas in this section we are interested in the estimation the function $f$ in model (3.1). Nevertheless, under additional assumptions we can investigate the properties of the derivative of the oracle estimator developed in what follows and a brief discussion of these relations is given in Example 3.1.

Another important difference between model (3.1) and the Gaussian white noise model commonly discussed in the literature of mathematical statistics lies in the fact that we consider a general error process $\{\varepsilon_{t}:t\in[0,1]\}$. In particular, we concentrate on Markovian Gaussian error processes with a covariance kernel of the form

[TABLE]

where $u(\cdot)$ and $v(\cdot)$ are some (known) functions defined on the interval $[0,1]$, such that $v(t)\neq 0$ for $t\in[0,1]$. Kernels of this form generalize the Brownian motion, which is obtained for $u(t)=t$, $v(t)=1$, and are called triangular kernels in the literature. The property (3.3) essentially characterizes a Gaussian process to be Markovian [see Doob, (1949) or Mehr and McFadden, (1965) for more details]. We assume that the process $\{\varepsilon_{t}:t\in[0,1]\}$ is non-degenerate on the open interval $(0,1)$, which implies that the function

[TABLE]

is positive on the interval $(0,1)$ and strictly increasing and continuous on $[0,1]$.

Regarding the estimation of the unknown function $f$, we propose to estimate the coefficients $\theta_{j}$ in the projection estimator (2.2) using statistics of the form [see Grenander, (1950)]

[TABLE]

where $\xi_{j}$ is a signed measure on the interval $[0,1]$ such that

[TABLE]

Obviously, this condition implies for the sequence of estimators $(\hat{\theta}_{j})_{j\in\mathbb{N}}$ that $\sum_{j=1}^{\infty}\mathbb{E}[\hat{\theta}^{2}_{j}]<\infty,$ and thus we can define the random variable

[TABLE]

In particular, if $\hat{f}^{(J)}(t)=\sum_{j=1}^{J}\hat{\theta}_{j}\varphi_{j}(t)$ is the truncated series from (3.7), we have that

[TABLE]

and the mean integrated squared error of the estimator $\hat{f}$ in (3.7) is given by

[TABLE]

We conclude that the optimal linear oracle estimator $\hat{f}$ of the function $f$ minimizing (3.8) can be determined minimizing the individual mean squared errors $\mathbb{E}[(\hat{\theta}_{j}-\theta_{j})^{2}]$ separately. Due to the definition of linear estimators in (3.5), this problem corresponds to the determination of a signed measure $\xi^{*}_{j}$ on the interval $[0,1]$, which minimizes the functional

[TABLE]

Remark 3.1

• (1)

Note that - in contrast to most of the literature - we do not assume that $\hat{\theta}_{j}$ is an unbiased estimator of the Fourier coefficient $\theta_{j}$ $(j\in\mathbb{N})$. A prominent unbiased estimator for $\theta_{j}$ is given by

[TABLE]

and for general unbiased estimates of the form (3.5) the condition (3.6) reduces to

[TABLE]

Moreover, if the kernel $K$ is continuous on $[0,1]\times[0,1]$ and if $\varphi_{1},\varphi_{2},\ldots$ are the eigenfunctions of the integral operator associated with the covariance kernel $K$ with corresponding eigenvalues $\lambda_{1},\lambda_{2},\ldots$, then condition (3.6) further reduces to

[TABLE]

• (2)

Under the additional assumption that the estimator (3.5) is unbiased for $\theta_{j}$, the second term in (3.9) vanishes and the resulting optimization problem corresponds to the problem of finding the best linear estimator in the location scale model $Y_{t}=\theta_{j}+\varepsilon_{t}$, which has been first studied in a seminal paper of Grenander, (1950). This author showed that under the additional constraint $\int_{0}^{1}d\xi_{j}(dt)=1$ the optimal solution $\xi_{j}^{*}$ minimizing $\int_{0}^{1}\int_{0}^{1}K(s,t)\xi_{j}(ds)\xi_{j}(dt)$ can be characterized by the property that the function $t\to\int_{0}^{1}K(s,t)\xi_{j}^{*}(ds)$ is constant on the interval $[0,1]$.

The following theorem provides a complete solution of the optimization problem (3.9) and is proven in the appendix. For a precise statement of the result we denote by $\delta_{x}$ the Dirac measure at the point $x$ and distinguish the following cases for the triangular kernel (3.3).

• (A)

$u(0)\neq 0$.

• (B)

$u(0)=0,f(0)=0$.

• (C)

$u(0)=0,f(0)\neq 0$.

Theorem 3.1

Consider the functional $\Psi_{j}$ in (3.9) with a twice differentiable regression function $f$ and a triangular covariance kernel of the form (3.3), where the functions $u$ and $v$ are also twice differentiable. For any $j\in\mathbb{N}$ the signed measure $\xi^{*}_{j}(dt)$ minimizing the functional $\Psi_{j}$ in the class of all signed measures on the interval $[0,1]$ is given by

[TABLE]

where $\theta_{j}$ is the $j$-th Fourier coefficient in the Fourier expansion (2.2). The values for $c$, $P_{0}$, $P_{1}$ and the function $p(\cdot)$ do not depend on the index $j$ and take different values corresponding to the properties of the functions $u(\cdot)$ and $f(\cdot)$. In particular, we have the following cases

1. (A)

If $u(0)\neq 0$, the quantities $c$, $P_{0}$, $P_{1}$ and $p$ are given by

[TABLE]

where the function $q$ is defined in (3.4). 2. (B)

If $u(0)=0$ and $f(0)=0$, the quantities $c$ and $P_{0}$ are given by

[TABLE]

and $P_{1}$ and $p$ are given by (3.15) and (3.16), respectively. 3. (C)

If $u(0)=0$ and $f(0)\neq 0$, the quantities $c$, $p(t)$ and $P_{1}$ are equal to zero, whereas $P_{0}$ is given by

[TABLE]

Corollary 3.1

Consider the regression model (3.1) with a twice differentiable regression function $f$ and a non-degenerate centered Gaussian error process $\{\varepsilon_{t}:t\in[0,1]\}$ with a triangular covariance kernel of the form (3.3), where the functions $u$ and $v$ are twice differentiable. The best linear oracle estimator minimizing the mean integrated squared error in (3.8) in the class of all linear estimators of the form (3.7) satisfying (3.6) is defined by

[TABLE]

where the coefficients $\hat{\theta}_{j}^{*}$ are given by

[TABLE]

and the signed measure $\xi^{*}_{j}(dt)$ is defined in Theorem 3.1. Moreover, the corresponding mean integrated squared error is given by

[TABLE]

where $c$ is defined in (3.13).

Note that Theorem 3.1 is a theoretical result as it requires knowledge of the unknown regression function $f$. Nevertheless, we will use it extensively in the following section to construct good estimators and corresponding optimal designs for series estimation in model (2.1).

Remark 3.2

• (1)

In model (2.1) with covariance kernel (3.3) and $u(0)=0$, the observation $Y_{0}$ at $t=0$ does not contain any error. Therefore, the value of $f(0)$ is known so that it can be checked whether case (B) or (C) of Theorem 3.1 holds.

• (2)

The estimator given in Theorem 3.1 depends on the orthonormal system of the series expansion via the parameter $\theta_{j}$.

• (3)

Using integration by parts the resulting estimator $\hat{\theta}^{*}_{j}$ in Theorem 3.1 can be represented as stochastic integral. For example, in case (A) (where $u(0)\neq 0$) the estimator can be represented as

[TABLE]

where the constant $c$ is defined in (3.13). Similarly in case (B) (where $u(0)=0$ and $f(0)=0$), the estimator can be represented by

[TABLE]

Finally, in case (C) (where $u(0)=0$ and $f(0)\neq 0$), the estimator directly reduces to

[TABLE]

In the latter case the estimator in (3.23) is not random, but fixed to the true - but unknown - parameter $\theta_{j}$.

Example 3.1

A very popular orthonormal basis of $L^{2}([0,1])$ is given by the trigonometric functions

[TABLE]

Under the assumptions of Theorem 3.1 we assume that $f$ and its derivative $\dot{f}$ can be represented as a trigonometric series, that is,

[TABLE]

Note that (under suitable assumptions) the Fourier coefficients in (3.25) and (3.26) are related by the equations

[TABLE]

If the error process $\{\varepsilon_{t}:t\in[0,1]\}$ in model (2.1) is given by a Brownian motion, we have $u(t)=t$, $v(t)=1$ in the definition of the triangular kernel (3.3) and thus $q(t)=t$. A straightforward application of Corollary 3.1 (case (B)) yields for the optimal oracle estimator of the function $f$

[TABLE]

where the estimated Fourier coefficients are given by

[TABLE]

(note that $f(0)=f(1)=0$). We thus also obtain an estimator of the function $\dot{f}$ in model (3.2) by taking the derivative of $f^{*}$ given n (3.28), that is,

[TABLE]

Using the relation (3.27), the estimator in (3.29) can be rewritten as

[TABLE]

and the mean integrated squared error of the estimator $\dot{f^{*}}$ in (3.30) is given by

[TABLE]

where we have used the representation c=\int_{0}^{1}\big{(}\dot{f}(t)\big{)}^{2}dt=\sum_{j=1}^{\infty}\bar{\theta}^{2}_{j}=\sum_{j=2}^{\infty}\bar{\theta}^{2}_{j} in the last equality.

It might be of interest to compare this estimator with the linear oracle estimator

[TABLE]

proposed in Tsybakov, (2009)[p. 67], where

[TABLE]

is used as the estimator of the Fourier coefficient $\bar{\theta}_{j}$ ($j=1,2,\ldots).$ This estimator is a shrinkage version of the unbiased estimator in (3.10) and the mean integrated squared error of $\dot{\tilde{f}}$ is given by

[TABLE]

Comparing (3.1) and (3.33), we observe that the oracle estimator $\dot{f^{*}}$, which is constructed by an application of Corollary 3.1, has a smaller mean integrated squared error than the estimator $\dot{\tilde{f}}$ defined in (3.32). **

4 Efficient series estimation from correlated data

In this section we apply the results from the continuous time model to construct optimal designs for series estimation of the function $f$ in model (2.1). In this transition from the continuous to the discrete model we are faced with several challenges. First, the signed measure defining the optimal oracle estimator $\hat{\theta}_{j}^{*}$ depends on the unknown function $f$ through its Fourier coefficients and through the constant $c$, and the function $f$ also appears in the stochastic integrals in (3.21) and (3.22). Secondly, we need to address the problem that even with preliminary knowledge of the function $f$, the stochastic integrals can not be computed since as the continuous time process $\{Y_{t}:t\in[0,1]\}$ is not observable. In order to overcome these difficulties and construct an implementable estimator, which does not require preliminary knowledge of $f$, we proceed to several steps, which are explained in detail below. Roughly speaking, these steps consist of a two stage estimation procedure, a truncation and an appropriate approximation of the stochastic integrals by sums, which can be calculated from the available data. In the latter step of this procedure we also determine the optimal design points.

Throughout this section we will restrict ourselves to the cases (A) and (B) of Theorem 3.1. For the case (C) we simply propose to replace the parameter value (3.23) by the best linear unbiased estimator derived in Dette et al., (2017).

4.1 Truncation in the continuous time model

In model (2.1) with $n$ observations, only a finite number, say $J$, of Fourier coefficients in the series expansion (2.2) can be estimated. For this reason, we consider for fixed $J\in\mathbb{N}$ the best $L^{2}$-approximation

[TABLE]

of the function $f$ by functions from the $\mbox{span}\{\varphi_{1},\ldots,\varphi_{J}\}$, space where the vectors $\theta^{(J)}$ and $\Phi^{(J)}$ are defined by $\theta^{(J)}=(\theta_{1},\ldots,\theta_{J})^{T}$ and $\Phi^{(J)}(t)=(\varphi_{1}(t),\ldots,\varphi_{J}(t))^{T}$, respectively. We now replace the function $f$ by the function $f^{(J)}$ in the estimators $\hat{\theta}_{1}^{*},\ldots,\hat{\theta}_{J}^{*}$ defined in (3.21) and (3.22) for cases (A) and (B) respectively. In case (A) this gives the vector

[TABLE]

where

[TABLE]

and the $J\times J$ matrix $C^{(J)}$ is defined by

[TABLE]

Similarly, in case (B) we obtain

[TABLE]

where

[TABLE]

and the $J\times J$ matrix $M^{(J)}$ is given by

[TABLE]

The resulting estimators (4.2) and (4.5) still depend on the first $J$ unknown Fourier coefficients $\theta_{1},\ldots,\theta_{J}$ and also depend on the full trajectory of the process $\{Y_{t}:t\in[0,1]\}$. This dependence will be removed in the following sections.

4.2 Discrete approximation of stochastic integrals

In concrete applications the integrals in (4.2) and (4.5) cannot be evaluated and have to be approximated from the given data. For this purpose we assume that $n$ observations $Y_{t_{1}},\ldots,Y_{t_{n}}$ from model (2.1) at $n$ distinct time points $0=t_{1}<t_{2}<\ldots<t_{n-1}<t_{n}=1$ are available and we consider the estimators

[TABLE]

as approximations of the quantities in (4.2) and (4.5), respectively. Note that $\hat{\theta}^{(J),*}$ depends on the full trajectory $\{Y_{t}:t\in[0,1]\}$, while $\hat{\theta}^{(J),n}$ is an approximation based on the sample $\{Y_{t_{i}}:i=1,\ldots,n\}$. In (4.8) and (4.9) $\mu_{2},\dots,\mu_{n}$ denote $J$-dimensional weights which depend on the time points $0=t_{1}<t_{2}<\ldots<t_{n-1}<t_{n}=1$ and will be chosen in an optimal way. In particular we propose to determine the weights $\mu_{2},\ldots,\mu_{n}$ such that the expected $L^{2}$-distance

[TABLE]

between $\hat{\theta}^{(J),*}$ and its discrete analogue $\hat{\theta}^{(J),n}$ is minimized, where $\|\cdot\|$ denotes the Euclidean norm in $\mathbb{R}^{J}$.

The following result provides an alternative expression of the expectation of this distance.

Proposition 4.1

Assume that the conditions of Theorem 3.1 are satisfied. The Euclidean distance between the estimators $\hat{\theta}^{(J),*}$ and $\hat{\theta}^{(J),n}$ can be represented as

[TABLE]

where the quantities $V$ and $B$ are defined by

[TABLE]

and the constant $k^{(J)}$ is given by

[TABLE]

Note that the expected $L^{2}$-distance in (4.11) only differs in the multiplicative factor $k^{(J)}$ for the different cases (A) and (B) and this factor does not depend on the vector-weights $\mu_{2},\ldots,\mu_{n}$. Therefore optimal weights minimizing the expected $L^{2}$-distance can be determined without distinguishing between the two cases (A) and (B).

The function $B$ in the criterion (4.11) still depends on the unknown regression function $f$ which we replace again by its truncation $f^{(J)}$ defined in (4.1). The resulting criterion is given by

[TABLE]

where

[TABLE]

We now determine the optimal weights such that the term $B^{(J)}(\mu_{2},\ldots,\mu_{n})$ in (4.13) vanishes for all potential Fourier coefficients $\theta_{1},\ldots,\theta_{J}$ in the function $f^{(J)}$. Therefore, the optimal weights are obtained by minimizing $\Phi$ in (4.13) under the constraint

[TABLE]

In this situation the criterion (4.13) reduces to the minimization of

[TABLE]

with respect to the weights $\mu_{2},\ldots,\mu_{n}$ (depending on the time points $0=t_{1}<t_{2},\ldots,t_{n-1}<t_{n}=1$). In order to simplify this optimization we introduce the following notation

[TABLE]

which however does not reflect the dependence on the time points. Using the notation in (4.17), the approximation of the expected $L^{2}$-distance in (4.16) can be rewritten in terms of the quantities $\gamma_{2},\ldots,\gamma_{n}$ as

[TABLE]

and the constraint (4.15) is given by

[TABLE]

where $M^{(J)}$ is the matrix defined in (4.7) (for both cases (A) and (B)). Note that both the function $\Psi$ and the constraint in (4.19) do not involve the function $f$ and only include assumptions concerning the first $J$ basis functions $\varphi_{1},\ldots,\varphi_{J}$ used in the approximation $f^{(J)}$.

The resulting optimization problem (4.18) with constraint (4.19) has the same structure as an optimization problem considered in Dette et al., (2017) and from the results in this paper we obtain the solution

[TABLE]

where the matrix $B$ is given by

[TABLE]

and $M^{(J)}$ and $\beta_{i}$ are defined in (4.7) and in (4.17), respectively. If the matrix $B$ is singular, we replace the inverse $B^{-1}$ in (4.20) by a generalized inverse $B^{-}$. Using the relation between $\gamma_{i}$ and $\mu_{i}$ in (4.17), we obtain the optimal weights

[TABLE]

Note that these weighs still depend on the design points $t_{2},\ldots,t_{n-1}$ which will be determined next.

4.3 Optimal designs for series estimation

Using the optimal $\gamma^{*}_{2},\ldots,\gamma_{n}^{*}$ given in (4.20) in the expression for the function $\Psi$ defined in (4.18), we obtain an appropriate optimal design criterion for the choice of the time points $0=t_{1}<t_{2}<\ldots t_{n-1}<t_{n}=1$. More precisely, for the optimal weights, the function $\Psi$ depends only on the design points and can be represented as the function

[TABLE]

where the matrices $B$ and $M^{(J)}$ are defined in (4.21) and (4.7), respectively and depend on $0=t_{1}<t_{2},\ldots,t_{n-1}<t_{n}=1$. The optimal design is now determined by minimizing the function $\tilde{\Psi}$, which is different from the criterion considered in Dette et al., (2017) for unbiased linear estimation in the linear regression model

[TABLE]

The optimal time points only depend on the first $J$ basis functions which are used for the estimator of the regression function $f$ and have to be determined numerically in all cases of practical interest. We will present some examples in Section 5.

4.4 The final estimate

With the optimal weights $\mu_{2}^{*},\ldots,\mu_{n}^{*}$ determined in Section 4.2 and the optimal time points $t_{2}^{*},\ldots t_{n-1}^{*}$ determined in Section 4.3, the estimators in (4.8) and (4.9) corresponding to the cases (A) and (B) are given by

[TABLE]

and

[TABLE]

respectively, where

[TABLE]

For their application we still require knowledge of the vector of Fourier coefficients $\theta^{(J)}$ and the constants $c^{(J)}$ and $m^{(J)}$ defined in (4.3) and (4.6) (note that these quantities also depend on $\theta^{(J)}$). For this purpose we propose to use the linear unbiased estimate derived by Dette et al., (2017) for the linear model (4.23). This estimate is defined as

[TABLE]

and the quantity $c^{(J)}$ in (4.3) is estimated by $\check{c}^{(J),n}=(\check{\theta}^{(J),n})^{T}C^{(J)}\check{\theta}^{(J),n}$. A straightforward calculation shows that the resulting estimator for the case (A) is given by

[TABLE]

which is a shrinkage version of the estimator $\check{\theta}^{(J),n}$ in (4.26).

For the case (B) similar arguments show that the estimator in (4.25) can also be rewritten in terms of the linear unbiased estimate $\check{\theta}^{(J),n}$, that is,

[TABLE]

where $\check{m}^{(J),n}=(\check{\theta}^{(J),n})^{T}M^{(J)}\check{\theta}^{(J),n}$. Here the structure of the estimator $\check{\theta}^{(J),n}$ depends on the structure of the basis functions contained in the vector $\Phi^{(J)}$ [see Section 5 in Dette et al., (2017) for more details].

5 Numerical results

In this section, we illustrate the properties of the estimator and the corresponding optimal design derived in Section 4 by means of a small simulation study. We consider a Gaussian process assuming both an exponential kernel and a Brownian motion as the error process in model (2.1). In both cases, we present the numerically calculated optimal time points with respect to the criterion defined in (4.22) and the corresponding simulated integrated mean squared errors for the estimator

[TABLE]

proposed in this paper and the estimator

[TABLE]

which is based on the best linear unbiased estimates in the tuncated Fourier expansion.

Throughout this section, we will use the trigonometric series defined in (3.24) as orthonormal basis of $L^{2}([0,1])$. We further assume that the unknown function $f$ is symmetric on the interval $[0,1]$ such that it is sufficient to use only the cosine functions in the series expansions of $f$. Consequently, the orthonormal system is given by

[TABLE]

In Section 5.1 we consider the exponential kernel, whereas in Section 5.2 we concentrate on the Brownian motion.

5.1 The Exponential kernel

We assume that the error process $\{\varepsilon_{t}:t\in[0,1]\}$ is a centered Gaussian process with an exponential kernel of the form $K(s,t)=\exp(-L|s-t|)$, where $L\in\mathbb{R}^{+}$ is a given constant. This can be represented in the triangular form (3.3) with $u(t)=\exp(Lt)$ and $v(t)=\exp(-Lt)$ and the function $q$ is obtained as $q(t)=u(t)/v(t)=\exp(2Lt)$. Therefore, we have $u(0)\neq 0$ (which corresponds to case (A)) and the preliminary estimator $\check{\theta}^{(J),n}$ in (4.26) is given by

[TABLE]

where the matrices $M^{(J)}$, $C^{(J)}$ and $B$ become

[TABLE]

The estimator $\hat{\theta}^{(J),n}$ proposed in this paper is given in equation (4.27) and the corresponding estimators of the function $f$ are defined in (5.1) and (5.2).

We first consider the exponential covariance kernel with $L=1$ and assume that three basis functions $\varphi_{1}(t)=1,\ \varphi_{2}(t)=\cos(2\pi t),\ \varphi_{3}(t)=\cos(4\pi t)$ are used in the series estimator, where $n=4$ and $n=7$ observations at different time points $0=t_{1}<t_{2}<\ldots<t_{n-1}<t_{n}=1$ can be taken. Note that one needs at least $n=4$ observations at different time points to guarantee that the matrix $B$ in the preliminary estimator $\check{\theta}^{(J),n}$ is non-singular. The optimal points are determined minimizing the criterion (4.22) by particle swarm optimization [see Clerc, (2006) for details] and the results are presented in the first row of Table 1. The second row shows the results for $L=5$ and interestingly the optimal points do not change substantially for different values of the constant $L$. Also note that all designs are nearly equidistant.

We now evaluate the performance of the different estimators and the optimal time points by means of a simulation study. For the sake of comparison we also consider non-optimized time points for the simulation, which are given by

[TABLE]

for the case $n=4$ and $n=7$, respectively.

In the simulation study we generate data according to model (2.1) with two regression functions

[TABLE]

(note that both proposed functions are symmetric with $f(0)=f(1)=0$). For each model the mean integrated squared error of the estimators $\hat{f}^{(J),n}$ and $\check{f}^{(J),n}$ defined in (5.1) and (5.2) respectively is determined. More precisely, if $S$ denotes the number of simulation runs and $\bar{f}_{\ell}$ is the estimator based on the $\ell$-th run (either $\hat{f}^{(J),n}$ and $\check{f}^{(J),n}$), the simulated mean integrated squared error, $\text{MISE}_{n}$, is given by

[TABLE]

where $f$, the “true” regression function under consideration, is either given by (5.5) or by (5.6). All results are based on $S=1000$ simulation runs.

For the case of the sample size $n=4$, the resulting mean integrated squared error of the different estimators (and corresponding optimal time points) is shown in the left part of Table 2. For instance, the mean integrated squared error of the estimator $\hat{f}^{(J),n}$ (based on the on the optimal design) is $1.72$, if the true function is given by (5.5), whereas it is $2.06$ if the observations are taken according to the non-optimized design (5.3). Thus, the optimal design yields a reduction by $17\%$ in the mean integrated squared error. The optimal design also yields a reduction of $15\%$ of the mean squared error of the preliminary estimator $\check{f}^{(J),n}$ (although it is not constructed for this purpose). We also observe that the new estimator $\hat{f}^{(J),n}$ clearly outperforms the estimator $\check{f}^{(J),n}$ in all cases under consideration (reduction of the mean squared error between $9\%$ and $12\%$).

For the case of the sample size $n=7$, the corresponding results are presented in the right part of Table 2 and we observe a similar behavior. The new estimator $\hat{f}^{(J),n}$ clearly outperforms $\check{f}^{(J),n}$ regardless of the design and model under consideration. On the other hand the improvement by the choice of the design is less visible compared to the case where the sample size is $n=4$. This means that the influence of design on the performance of the estimators decreases with increasing sample size. The reason for this observation lies in the fact that in the models under consideration the discrete model (2.1) already provides a good approximation of the continuous model (3.1) for the sample size $n=7$. As in this model the full trajectory is available the impact of the design is negligible for sample sizes larger than $10$. As a consequence, a larger sample size would not decrease the integrated mean squared error substantially either. A similar effect was also observed by Dette et al., (2017) in the linear regression model with correlated observations.

Next we consider a situation where the correlation between the different observations is smaller and so we use the constant $L=5$ for the exponential kernel. The time points minimizing the criterion $\tilde{\Psi}$ in (4.22) are depicted in the second row of Table 1 for $n=4$ and $n=7$. The simulated mean integrated squared error of the estimators $\hat{f}^{(J),n}$ and $\check{f}^{(J),n}$ defined in (5.1) and (5.2) are displayed in Table 3 for the cases of sample size $n=4$ and $n=7$. When the sample size is $n=4$, we observe that the optimal design yields a substantial reduction in the mean squared errors of both estimators (between $65\%$ and $70\%$). Compared to the case $L=1$ (see Table 2) the reduction is larger. When the sample size is $n=7$ the mean integrated squared error of the estimators based on the optimal time points are slightly smaller compared to the non-optimized time points. We observe again that the influence of the position of the time points, and thus of the design, decreases if the sample size $n$ increases (see Table 3). A comparison of the two estimators (5.1) and (5.2) shows again that the new estimator $\hat{f}^{(J),n}$ outperforms the estimator $\check{f}^{(J),n}$ in all cases under consideration (reduction of the mean squared error between $16\%$ and $27\%$).

5.2 Brownian motion

We now consider the case where the error process in (2.1) is given by a Brownian motion, that is $K(s,t)=s\wedge t\,$, which can be represented by $K(s,t)=s,\quad s\leq t.$ Therefore, the functions $u$ and $v$ in (3.3) are given by $u(t)=t$ and $v(t)=1$, respectively, and the function $q$ is obtained as $q(t)=u(t)/v(t)=t$. This situation corresponds to case (B), where $u(0)=0$ and $f(0)=0$. The estimator $\hat{\theta}^{(J),n}$ is given by

[TABLE]

where the matrices $M^{(J)}$, $B$ and the constant $m^{(J)}$ are of the form

[TABLE]

Note that both the first row and the first column of the matrices $M^{(J)}$ and $B$ are zero (since $\varphi_{1}(t)=1$), such that both matrices are singular. Consequently, as proposed in Section 4, we use the generalized inverse

[TABLE]

of $B$, where the matrix $\tilde{B}$ is given by

[TABLE]

Here the vector $\bf{0}_{(J-1)}$ is of dimension $(J-1)$ with zero entries. where the matrix $\mathbf{I}_{(J-1)\times(J-1)}$ is the $(J-1)$ dimensional identity matrix. The estimator $\check{\theta}^{(J),n}$ is obtained from Section 5.2 in Dette et al., (2017)

[TABLE]

where the matrix $\underline{C}^{(J)}$ is of the form

[TABLE]

We now analyze the behavior of the resulting estimators of the function $f$ if the first three basis functions are used for the series estimator and $n=4$ or $n=7$ observations at different time points $0=t_{1}<t_{2}<\ldots<t_{n-1}<t_{n}=1$ are available. The optimal time points minimizing the criterion (4.22) derived in Section 4 are given by

[TABLE]

for sample sizes $n=4$ and $n=7$ respectively. Note that the optimal time points (5.8) and (5.9) differ from the optimal time points for the case of the exponential kernel displayed in Table 1. This indicates that the position of the optimal time points depends on the structure of the covariance kernel.

The resulting mean integrated squared errors of the estimators $\hat{f}^{(J),n}$ and $\check{f}^{(J),n}$ are displayed in Table 4, where we again consider the comparative set of time points depicted in (5.3) and (5.4). We obtain similar results as in Section 5.1. More specifically, for the case of sample size $n=4$, we observe that the optimal design yields a substantial reduction in the mean squared errors of both estimators (see the left part of Table 4). When the sample size is $n=7$, the difference between the optimal time points and the design (5.4) is less visible.

A comparison of the two estimators shows a different behavior as in Section 5.1, that is, unlike the case of an exponential Kernel, when the error process is a Brownian motion, both estimators perform well and they have similar (small) mean integrated squared errors (see Table 4).

Acknowledgements This work has been supported in part by the Collaborative Research Center “Statistical modeling of nonlinear dynamic processes” (SFB 823, Teilprojekt C2) of the German Research Foundation (DFG) and by a grant from the National Institute of General Medical Sciences of the National Institutes of Health under Award Number R01GM107639. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

6 Technical details

Proof of Theorem 3.1 We restrict ourselves to the proof of the result in case (A), the other cases can be proved in a similar way. Note that the function $\Psi_{j}$ is convex on the space of all signed measures and therefore, a signed measure $\xi_{j}^{*}$ minimizes $\Psi_{j}$ if and only if the directional derivative from $\xi^{*}_{j}$ in any direction is nonnegative, that is

[TABLE]

for all signed measures $\eta$ on the interval $[0,1]$. A straightforward calculation gives

[TABLE]

Consequently, the signed measure $\xi^{*}_{j}$ minimizes $\Psi_{j}$ if and only the inequality

[TABLE]

is satisfied for all signed measures $\eta$ on the interval $[0,1]$.

In order to check (5.2) for the signed measure $\xi^{*}_{j}$ we calculate each term in (5.2) separately, where we use the following representation of the quantities in (3.13)-(3.16)

[TABLE]

To simplify (5.2) we note that integration by parts yields

[TABLE]

Similarly, we obtain

[TABLE]

where we have used again integration by parts for the third equality. Consequently, we get

[TABLE]

[TABLE]

and thus the left hand side of (5.2) reduces to

[TABLE]

for an arbitrary signed measure $\eta$. This proves that (5.2) holds and the signed measure $\xi^{*}_{j}$ defined in Theorem 3.1 minimizes the function $\Psi_{j}$.

Proof of Proposition 4.1 For the term on the left hand side of equation (4.11) we obtain

[TABLE]

For the determination of the expected value inside the trace

[TABLE]

we use a transformation of the Gaussian process $\{Y_{t}:t\in[0,1]\}$ to a Brownian motion, as it was introduced by Doob, (1949). This result shows that the error process $\{\varepsilon_{t}:t\in[0,1]\}$ with covariance kernel (3.3) can be represented by

[TABLE]

where $W=\{W(s):s\in[q(0),q(1)]\}$ is a Brownian motion on the interval $[q(0),q(1)]$. We use this relationship to represent the process $\{Y_{t}:t\in[0,1]\}$ as

[TABLE]

Dividing by $v(t)$ and using the transformation $s=q(t)$, we get the transformed model

[TABLE]

where

[TABLE]

Consequently, we obtain for arbitrary $0\leq t_{i-1}<t_{i}\leq 1$

[TABLE]

where the function $\tilde{\Phi}^{(J)}(s)$ is given by $\tilde{\Phi}^{(J)}(s)=\frac{\Phi^{(J)}(q^{-1}(s))}{v(q^{-1}(s))}$ and we set $t=q^{-1}(s)$. This gives for the transformed derivatives

[TABLE]

We now introduce the notation

[TABLE]

As $W$ is a Brownian motion, the random variables $X_{2},\ldots,X_{n}$ are independent and the expected value in (5.3) can be rewritten as

[TABLE]

Obviously

[TABLE]

and Itô’s isometry gives

[TABLE]

Inserting these representations in (5.4) results in (4.11), which proves Proposition 4.1.