# Adaptive Non-parametric Estimation of Mean and Autocovariance in Regression
with Dependent Errors

Paulo Serra 111 Department of Mathematics and Computer Science, Eindhoven
University of Technology, P.O. Box 513, 5600 MB Eindhoven, the Netherlands;
<EMAIL_ADDRESS>Tatyana Krivobokova 222 Institute for Mathematical
Stochastics, Georg-August-Universität-Göttingen, Germany. Francisco Rosales
333 Department of Finance, Universidad del Pacífico, Peru.

###### Abstract

We develop a fully automatic non-parametric approach to simultaneous
estimation of mean and autocovariance functions in regression with dependent
errors. Our empirical Bayesian approach is adaptive, numerically efficient and
allows for the construction of confidence sets for the regression function.
Consistency of the estimators is shown and small sample performance is
demonstrated in simulations and real data analysis. The method is implemented
in the R package eBsc that accompanies the paper.

Key words and phrases. Demmler-Reinsch basis, empirical Bayes, spectral
density, stationary process.

##  1 Introduction

Consider a fixed design non-parametric regression model

$\displaystyle\begin{aligned}
&&Y_{i}=f(t_{i})+\sigma\epsilon_{i},\;\mathbb{E}(\epsilon_{i})=0,\;\sigma>0,\;\;i=1,\ldots,n\\\
&&\mbox{cov}(\epsilon_{i},\epsilon_{j})=r(i-j)=r_{|i-j|}\in[-1,1],\;\;i,j=1,\ldots,n,\end{aligned}$
(1)

where $r(\cdot)$ denotes the autocorrelation function of the noise process,
and $r_{|i|}$ denotes autocorrelation at lag $i$. Typically, design points
$t_{i}\in\mathbb{R}$ are equidistant and represent time points. Functions $f$
and $r$ are unknown, but stationarity of $\\{\epsilon_{i}\\}_{i=1}^{n}$ and
smoothness of $f$ assumed. Observations $\\{Y_{i}\\}_{i=1}^{n}$ might be
measures of some experimental quantity observed with a time dependent
measurement error. In this case estimation of $f$ is of interest, while $r$ is
considered as a nuisance parameter. Another instance is
$\\{Y_{i}\\}_{i=1}^{n}$ being some time or space indexed stochastic process
with a seasonal or other deterministic effect, described by $f$. In this case
the focus is rather on estimation of the autocorrelation $r$.

Obviously, having a consistent estimator of $f$ would deliver a reasonable
estimator of $r$. Unfortunately, it is impossible to get a consistent
estimator of $f$ from a single replication without accounting for $r$ as shown
in Hart and Wehrley, (1986). Consider a time series of hourly loads (kW) for a
US utility (grey line in Figure 1), described in detail in Section 4. This
process has clear seasonal effects over years and possibly over weeks and
days. These deterministic effects can be modelled by a smooth function $f$. If
no parametric assumptions on $f$ is made and errors are treated as i.i.d.,
then all standard non-parametric estimators of $f$ are heavily affected by the
ignored dependence in the errors. In particular, automatic selectors of the
smoothing parameter (e.g., cross-validation) choose a biased smoothing
parameter leading to a nearly interpolating estimator in case of positively
correlated $\\{\epsilon_{i}\\}_{i=1}^{n}$; see the left plot in Figure 1.

|
---|---
Figure 1: Hourly loads (kW) for a US utility estimated ignoring dependence in
the data (left) and assuming that errors follow an autoregressive process of
order one (right). Grey line shows the data and black lines show estimators.
Estimators are obtained by functions gam (left) and gamm (right) of the R
package mgcv.

This problem has been known for several decades; for an overview see Opsomer
et al., (2001). Since then, many approaches have been suggested for regression
with dependent errors. The basic idea is to take smoothing parameter selection
criteria that are geared towards independent data, and modify them to take
into account the dependence in the errors.

One group of methods consists of making an explicit, parametric assumption on
the correlation structure, e.g., assuming that errors follow an ARMA$(p,q)$
process. Once $r$ is parametrised, the usual smoothing parameter selection
criteria are adjusted to incorporate $r$ and both $f$ and $r$ are estimated
simultaneously. For example, Hart, (1994) introduces a time series cross-
validation for $f$ estimated by kernel estimators and assuming errors to
follow an $AR(p)$ process, see also Altman, (1990) and Hall and Keilegom,
(2003). Kohn et al., (1992) use spline smoothing for estimation of $f$ and
general ARMA$(p,q)$ model for the errors, estimating all parameters from
either generalised cross validation or maximum likelihood. Similar ideas are
employed in Wang, (1998), Durban and Currie, (2003) and Krivobokova and
Kauermann, (2007). Smoothing with low-rank splines and an ARMA$(p,q)$ model
for the residuals is implemented in the gamm function of the R package mgcv.
This function has been used to estimate hourly loads in the right plot of
Figure 1, assuming that the errors follow an $AR(1)$ process. This approach
makes for computationally attractive algorithms, and the impact of the
correlation structure on the estimation of the regression function is clearer.
However, there are several drawbacks. First, the true correlation structure
might be much more complex than any parametric one, which is especially
problematic if the focus is on estimation of $r$. Moreover, if the correlation
structure is indeed misspecified, then estimators of $f$ are typically
strongly affected and perhaps even inconsistent. Maximum likelihood methods do
seem to be more robust to misspecification, see Krivobokova and Kauermann,
(2007). More importantly, in practice there is no information on $r$ is
available and some strategies for the selection of a parametric model for $r$
(its type and order) are needed.

Another group of methods makes no parametric assumption on $r$ and tries to
eliminate the influence of the error dependence on the smoothing parameter in
different ways. For example, Chu and Marron, (1991) and Hall et al., (1995)
study the modified cross-validation criterion obtained by leaving out a whole
block of $2l+1$ observations around each observation. Chiu, (1989) and Hurvich
and Zeger, (1990) study Mallow’s $C_{p}$ and cross validation in the frequency
domain. A different route is taken by Herrmann et al., (1992) and Lee et al.,
(2010): in the smoothing parameter selection criteria they incorporate some
sample estimators of autocovariances, that depend on unknown parameters linked
to assumptions on the error process. Using the method of Herrmann et al.,
(1992) for estimating the hourly loads gives an estimator of $f$ very similar
to the one obtained with function gamm in Figure 1.

All together, both groups of methods are rather focused on obtaining a
reasonable estimator for $f$, while $r$ is treated as a nuisance parameter.
Moreover, all these methods depend on unknown parameters that require some
knowledge of $r$. Hence, there is a need for a fully automatic (independent on
any unknown parameters) method for non-parametric estimation of both $f$ and
$r$ simultaneously.

Note that estimation of the correlation structure of a data vector with a
known mean is already a challenging problem. It is important in multivariate
analysis in problems such as clustering, principal component analysis (PCA),
linear- and quadratic discriminant analysis, and regression analysis. Two
frameworks are usually considered when the mean is known: either there are $n$
independent observations of a $p-$dimensional vector with correlated
components, cf. Bickel and Levina, 2008a ; Bickel and Levina, 2008b ; or there
is one observation of an $n-$dimensional vector sampled from a stationary
process, cf. Xiao and Wu, (2012). The second framework is the relevant one for
our problem. Irrespectively of the framework, natural (moment) estimators for
correlation matrix of the observed vectors are not consistent (in, e.g.,
operator-norm) and regularisation is needed (e.g., banding, tapering,
thresholding) to ensure positive definiteness and consistency of the estimate.
The minimax rates of convergence of these estimators can be shown to depend on
the decay of the autocorrelation function as the lag increases; cf. Yang et
al., (2001), Cai et al., (2010), Purahmadi, (2011), Xiao and Wu, (2012), see
also Fan et al., (2016).

In this paper we develop a likelihood based method that provides an adaptive
estimator of the regression function $f$, as well as estimators of the noise
level $\sigma$ and autocorrelations $r$. Starting from an arbitrary guess for
the autocorrelations we iteratively update all other parameters until
convergence. Furthermore, our empirical Bayesian framework provides a
computationally attractive way of construction confidence sets for the
regression function that take the correlation structure into account.

There are quite a few novel points in our work. Contrary to other approaches
in the literature, our method is completely automatic so that no tuning
parameters need to be set by the user, and it is also fully non-parametric. To
the best of our knowledge, our estimate of the autocorrelations is also novel
– we use spline smoothers to estimate both the regression function and the
spectral density of the noise process adaptively. The autocorrelations are
then reconstructed from the estimate of the spectral density, rather than by
tapering or thresholding some empirical estimate.

This paper is structured as follows. In Section 2 we introduce our estimator,
Section 3 contains simulation results, Section 4 presents two real data
examples, and Section 5 closes the paper with some conclusions. All of the
asymptotics and respective proofs, as well as auxiliary results are gathered
in the Appendix.

##  2 Construction of the estimators

Assume that in model (1) the design points are given by $t_{i}=(i-1)/(n-1)$,
regression function $f$ is a function from a Sobolev space
$\mathcal{W}_{\beta}$, $\beta>1/2$ and the noise terms $\epsilon_{i}$ are
sampled from a stationary, Gaussian noise process with zero mean and variance
$\sigma^{2}>0$. Setting $\mbox{\boldmath{$Y$}}=(Y_{1},\dots,Y_{n})^{T}$,
$\mbox{\boldmath{$f$}}=f(\mbox{\boldmath{$t$}})=\\{f(t_{1}),\dots,f(t_{n})\\}^{T}$,
$\mbox{\boldmath{$\epsilon$}}=(\epsilon_{1},\ldots,\epsilon_{n})^{T}$ we can
write this model as

$\mbox{\boldmath{$Y$}}=\mbox{\boldmath{$f$}}+\sigma\mbox{\boldmath{$\epsilon$}},\;\;\;\mbox{\boldmath{$\epsilon$}}\sim
N(\mbox{\boldmath{$0$}},\mbox{\boldmath{$\mathcal{R}$}}).$

The entries of $\mathcal{R}$ are ${\mathcal{R}}_{i,j}=r_{|i-j|}$, so that
$\mathcal{R}$ is a Toeplitz matrix, such that $r_{|i-j|}=0$ for $|i-j|>m$, for
some fixed $m\in\mathbb{N}$ independent on $n$. We further assume that for
some $0<\delta<1$, $\mbox{\boldmath{$\mathcal{R}$}}\in\mathcal{M}_{n,\delta}$,
where $\mathcal{M}_{n,\delta}$ denotes the space of $n\times n$ matrices whose
eigenvalues all lie on the interval $[\delta,\,1/\delta]$. Moreover, the
spectral density of $\\{\epsilon_{i}\\}_{i=1}^{n}$ is at least Lipschitz
continuous.

We estimate $f$ using a smoothing spline, i.e., we find $\hat{f}$ that solves

$\min_{f\in\mathcal{W}_{q}}\left[\frac{1}{n}\big{\\{}\mbox{\boldmath{$Y$}}-f(\mbox{\boldmath{$t$}})\big{\\}}^{T}\mbox{\boldmath{$R$}}^{-1}\big{\\{}\mbox{\boldmath{$Y$}}-f(\mbox{\boldmath{$t$}})\big{\\}}+\lambda\int_{0}^{1}\big{\\{}f^{(q)}(t)\big{\\}}^{2}\,dt\right],$
(2)

for some $q\in\mathbb{N}$, $\lambda>0$ and some correlation matrix
$\mbox{\boldmath{$R$}}\in\mathcal{M}_{n,\delta}$. Note that $\beta$ and
$\mathcal{R}$ are unknown in practice and are replaced by some “working”
values in (2). Subsequently, $\beta$, $\lambda$ and $\mathcal{R}$ are
estimated from the data using the empirical Bayes approach. It is well-known
that for given $R$, $q$ and $\lambda$, the resulting estimator $\hat{f}$ is a
natural spline of degree $2q-1$ with knots at $t$ and can be written as
$\hat{f}(t)=\mbox{\boldmath{$S$}}(t)\mbox{\boldmath{$Y$}}$, where $S$ is a
$n\times n$ smoother matrix. To represent $S$ we choose the so-called Demmler-
Reinsch basis of the natural spline space of degree $2q-1$, which is defined
in Section B.1 in the Appendix,

$\displaystyle\mbox{\boldmath{$S$}}=\mbox{\boldmath{$S$}}_{\lambda,q,\mbox{\boldmath{$R$}}}$
$\displaystyle=\mbox{\boldmath{$\Phi$}}_{q}\left\\{\mbox{\boldmath{$\Phi$}}_{q}^{T}\mbox{\boldmath{$R$}}^{-1}\mbox{\boldmath{$\Phi$}}_{q}+\lambda\mbox{diag}(n{\mbox{\boldmath{$\eta$}}}_{q})\right\\}^{-1}\mbox{\boldmath{$\Phi$}}_{q}^{T}\mbox{\boldmath{$R$}}^{-1},$
(3)

where $\mbox{\boldmath{$\Phi$}}_{q}$ is the $n\times n$ basis matrix and
$\mbox{\boldmath{$\eta$}}_{q}\in\mathbb{R}^{n}$ is a vector of eigenvalues.
This representation makes the dependence of $S$ on the parameters $\lambda$,
$R$ and $q$ more explicit. Subsequently, to keep the notation simple we omit
the dependence on these parameters, except if these are set to a particular
value. Note that estimation of $\beta$ from the data makes our estimator
adaptive to the unknown smoothness of the signal; see Serra and Krivobokova,
(2017).

###  2.1 Bayesian interpretation

The estimator $\hat{f}$ has a Bayesian interpretation which provides us with a
convenient way of estimating all of the unknown parameters by employing the
empirical Bayes approach.

We start by endowing $f$ with a partly informative prior – given
$(\mbox{\boldmath{$t$}},\lambda,q,\sigma^{2},\mbox{\boldmath{$R$}})$, the
prior on $f$ admits a density proportional to

$\displaystyle\left|\frac{\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$S$}}^{-1}-\mbox{\boldmath{$I$}}_{n})}{2\pi\sigma^{2}}\right|^{1/2}_{+}\exp\left\\{-\frac{\mbox{\boldmath{$f$}}^{T}\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$S$}}^{-1}-\mbox{\boldmath{$I$}}_{n})\mbox{\boldmath{$f$}}}{2\sigma^{2}}\right\\},$
(4)

where $|\cdot|_{+}$ denotes the product of non-zero eigenvalues ($S$ has
exactly $q$ eigenvalues equal to $1$). This prior is often used for Bayesian
smoothing splines and corresponds to a non-informative part on the null-space
of
$\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$S$}}^{-1}-\mbox{\boldmath{$I$}}_{n})$
and a proper Gaussian prior on the remaining space. Note that this prior
distribution happens to be independent of $R$. This follows from the identity
$\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$S$}}^{-1}-\mbox{\boldmath{$I$}}_{n})=\mbox{\boldmath{$S$}}_{I}^{-1}-\mbox{\boldmath{$I$}}_{n}$,
where $\mbox{\boldmath{$S$}}_{I}$ denotes the smoother matrix with
$\mbox{\boldmath{$R$}}=\mbox{\boldmath{$I$}}_{n}$; cf. Appendix B.2 for the
derivation of the identity. It is well-known that the corresponding posterior
distribution for
$\mbox{\boldmath{$f$}}|(\mbox{\boldmath{$t$}},\lambda,\sigma^{2},\mbox{\boldmath{$R$}})$
has a mean equal to the smoothing spline estimator
$\hat{\mbox{\boldmath{$f$}}\,\,}=\mbox{\boldmath{$S$}}\mbox{\boldmath{$Y$}}$
and variance $\sigma^{2}\mbox{\boldmath{$S$}}\mbox{\boldmath{$R$}}$; cf.
Speckman and Sun, (2003). The variance $\sigma^{2}$ given
$(\mbox{\boldmath{$t$}},\lambda,q,\mbox{\boldmath{$R$}})$ is endowed with an
inverse-gamma prior $\mbox{IG}(a,b)$, $a,b>0$.

The resulting prior on
$(\mbox{\boldmath{$f$}},\sigma^{2})|(\lambda,q,\mbox{\boldmath{$R$}})$ is
conjugate for model (1) in the sense that the posterior distribution on
$(\mbox{\boldmath{$f$}},\sigma^{2})|(\lambda,q,\mbox{\boldmath{$R$}})$ is a
known distribution. Namely, the marginal posterior for $\sigma^{2}$ given
$(\lambda,q,\mbox{\boldmath{$R$}})$ is an inverse gamma distribution with
shape parameter $(n-q+2a)/2$ and scale parameter
$\\{\mbox{\boldmath{$Y$}}^{T}\mbox{\boldmath{$R$}}^{-1}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$Y$}}+2b\\}/2$.
As for $\mbox{\boldmath{$f$}}|(\lambda,q,\mbox{\boldmath{$R$}})$, its
posterior is a multivariate t-distribution with $n+1$ degrees of freedom, mean
$\hat{\mbox{\boldmath{$f$}}\,\,}=\mbox{\boldmath{$S$}}\mbox{\boldmath{$Y$}}$,
and scale $\hat{\sigma}^{2}\mbox{\boldmath{$S$}}\mbox{\boldmath{$R$}}$.

It remains to estimate unknown $\lambda$, $q$ and $R$. To do so, we employ the
empirical Bayes approach and estimate these parameter from the marginal
distribution of $Y$, given
$(\mbox{\boldmath{$t$}},\lambda,q,\mbox{\boldmath{$R$}})$, which is a
multivariate $t$-distribution whose density is

$\frac{\Gamma\\{a+(n-q)/2\\}\left|\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})(2a-q)/(2b)\right|_{+}^{1/2}}{\\{\pi(2a-q)\\}^{n/2}\Gamma(a-q/2)}\left\\{1+\frac{\mbox{\boldmath{$Y$}}^{T}\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$Y$}}}{2b}\right\\}^{-(n+2a-q)/2};$

see Kotz and Nadarajah, (2004). It remains to set the parameters $a$ and $b$.
From our analysis we concluded that asymptotically the choice of parameters
$a$ and $b$ is irrelevant (as long as $a$ and $b$ are $o(n)$ and lead to
proper a prior and marginal distributions). Hereon after, to simplify the log-
likelihood, we set $b=1/2$ and $a=(q+1)/2$ and obtain

$\displaystyle\ell_{n}(\lambda,q,\mbox{\boldmath{$R$}})=-\frac{n+1}{2}\log\left\\{\mbox{\boldmath{$Y$}}^{T}\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$Y$}}+1\right\\}+\frac{1}{2}\log\left|\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})\right|_{+},$
(5)

up to an additive constant that is independent of the parameters of interest.
With this choice of $a$ and $b$ the posterior mean for $\sigma^{2}$ becomes

$\hat{\sigma}^{2}=\frac{\mbox{\boldmath{$Y$}}^{T}\mbox{\boldmath{$R$}}^{-1}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$Y$}}+1}{n+1}.$
(6)

###  2.2 Estimating equations and algorithm

By differentiating the log-likelihood (5) (see identities for the derivatives
in Appendix B.2) with respect to the smoothing parameter $\lambda$, values of
the spectral density $\rho_{i}=\rho(\pi t_{i})$, $i=1,\ldots,n$ that
correspond to $R$ and penalty order $q$, we obtain the following estimating
equations for these parameters.

$\displaystyle T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})$
$\displaystyle=\frac{1}{n_{\lambda,q}}\left[\mbox{\boldmath{$Y$}}^{T}\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$S$}}\mbox{\boldmath{$Y$}}-\hat{\sigma}^{2}\left\\{\mbox{tr}\left(\mbox{\boldmath{$S$}}\right)-q\right\\}\right],$
$\displaystyle T_{q}(\lambda,q,\mbox{\boldmath{$R$}})$
$\displaystyle=\frac{1}{n_{\lambda,q}^{\prime}}\left[\mbox{\boldmath{$Y$}}^{T}\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$\Phi$}}_{q}\mbox{diag}\left[n\mbox{\boldmath{$\eta$}}_{q}\circ\log(n\mbox{\boldmath{$\eta$}}_{q})\right]\mbox{\boldmath{$\Phi$}}_{q}^{T}\mbox{\boldmath{$Y$}}-\hat{\sigma}^{2}\mbox{tr}\left\\{\mbox{\boldmath{$n$}}\mbox{\boldmath{$\eta$}}_{q}\circ\log(n\mbox{\boldmath{$\eta$}}_{q})\right\\}\right],$
$\displaystyle T_{\rho_{i}}(\lambda,q,\mbox{\boldmath{$R$}})$
$\displaystyle=\mbox{tr}\left[\mbox{\boldmath{$R$}}^{-1}\frac{\partial\mbox{\boldmath{$R$}}}{\partial\rho_{i}}\mbox{\boldmath{$R$}}^{-1}\left\\{(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$YY$}}^{T}(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})^{T}-\hat{\sigma}^{2}(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$R$}}\right\\}\right],\quad
i=1,\dots,n,$

where $\circ$ denotes the Hadamard product and $\hat{\sigma}^{2}$ as given in
(6). The scaling factors in the first two equations are given by
$n_{\lambda,q}=\lambda^{-1/(2q)}+n\lambda^{1/(2q)}$ and
$n_{\lambda,q}^{\prime}=n_{\lambda,q}(\log\lambda)^{2}$, $\lambda>0$,
$q\in\mathbb{N}$. These equations might be difficult to evaluate and in
practice we use approximate expressions, which are given in Appendix B.5.

These estimating equations need to be solved simultaneously. In practice we
proceed as follows. For each fixed $q$ we start with an initial guess
$\hat{\mbox{\boldmath{$R$}}\,}^{(0)}$ (typically just an identity matrix),
obtain a preliminary estimate $\hat{\lambda}^{(0)}$ and iterate to get
$\hat{\lambda}_{q}$ and $\hat{\mbox{\boldmath{$R$}}\,}_{q}$. Finally, $q$ is
chosen to solve
$T_{q}(\hat{\lambda}_{q},q,\hat{\,\mbox{\boldmath{$R$}}_{q}})=0$.

Since $\rho_{i}=\rho(\pi t_{i})$ are values of a smooth function $\rho$ at
given points, estimation of $\rho$ should be carried out over a space of
smooth functions. In the Bayesian framework this can be accomplished by
introducing a suitable prior on $R$, which then acts as a penalty term in the
posterior. For simplicity, we perform a two-step procedure instead. First,
$\tilde{\rho}_{i}$ are obtained as solutions of the corresponding estimating
equations. Second, $\tilde{\rho}_{i}$ are smoothed, i.e.,
$\hat{\mbox{\boldmath{$\rho$}}\,}={\mbox{\boldmath{$S$}}}_{\xi,p,\mbox{\boldmath{$I$}}}\tilde{\mbox{\boldmath{$\rho$}}\,}$,
where ${\mbox{\boldmath{$S$}}}_{\xi,p,{\mbox{\boldmath{$I$}}}}$ is a smoother
matrix (3) with parameters $\xi$, $p$ and $\mbox{\boldmath{$I$}}_{n}$.

All together, for each fixed $q$, at each step $\hat{\lambda}^{(j)}$ and
$\tilde{\rho}_{i}^{(j)}$ are obtained as solutions of the corresponding
estimating equations and $\tilde{\rho}_{i}^{(j)}$ are smoothed to get
$\hat{\rho}_{i}^{(j)}$. The algorithm is iterated until convergence of
$\hat{\lambda}$ and $\hat{\rho}$, which are then used to get $\hat{q}$.
Finally, $\hat{\mbox{\boldmath{$R$}}\,}$ is recovered from
$\hat{\mbox{\boldmath{$\rho$}}\,}$ by the discrete Fourier transform; for the
details see Appendix B.4. The summary of the estimation procedure is given in
Algorithm 1.

for  _$q$ in $Q_{n}$_ do

set $j=1$ and $\hat{\mbox{\boldmath{$\rho$}}\,}^{(0)}$

while  _stopping criterium not met_ do

set $\hat{\lambda}^{(j)}$ to a solution of
$T_{\lambda}(\lambda,q,\hat{\mbox{\boldmath{$\rho$}}}^{(j-1)})=0$ ;

set $\tilde{\mbox{\boldmath{$\rho$}}\,}^{(j)}$ to a solution of
$T_{\mbox{\boldmath{$\rho$}}}(\lambda^{(j)},q,\hat{\mbox{\boldmath{$\rho$}}}^{(j-1)})=0$;

compute $\hat{\mbox{\boldmath{$\rho$}}\,}^{(j)}$ by smoothing
$\tilde{\mbox{\boldmath{$\rho$}}\,}^{(j)}$;

set $j=j+1$;

end while

set $\hat{\lambda}_{q}=\hat{\lambda}^{(j)}$ and
$\hat{\mbox{\boldmath{$\rho$}}\,}_{q}=\hat{\mbox{\boldmath{$\rho$}}}^{(j)}$;

end for

set $\hat{q}$ to a solution of
$T_{q}(\hat{\lambda}_{q},q,\hat{\mbox{\boldmath{$\rho$}}\,}_{q})=0$ over
$Q_{n}$;

set $\hat{\lambda}$ to $\hat{\lambda}_{\hat{q}}$;

set
$\hat{\mbox{\boldmath{$\rho$}}\,}=\hat{\mbox{\boldmath{$\rho$}}\,}_{\hat{q}}$;

set $(\hat{\mbox{\boldmath{$R$}}\,})_{i,j}=\hat{r}_{|i-j|}$, with
$\hat{r}_{k}=n^{-1}\sum_{l=1}^{n}\cos\left(k\pi\\{l-1\\}/\\{n-1\\}\right)\hat{\rho}_{l}$;

Algorithm 1 Recursive estimation procedure.

Here $Q_{n}$ denotes collections of values for $q$. The stopping rule is
standard: after each iteration we compare the change in the value of the
estimate of $\lambda$ and the norm of the change in the estimate of $\rho$; if
these fall bellow a threshold, then we stop iterating. The consistency of
estimators $\hat{\mbox{\boldmath{$R$}}}$, $\hat{\lambda}$ and $\hat{q}$ is
addressed in Appendix A.

###  2.3 Confidence sets

Once estimates $\hat{\lambda}$, $\hat{\mbox{\boldmath{$R$}}\,}$, $\hat{q}$,
and $\hat{\sigma}^{2}$ for respectively $\lambda$, $\mathcal{R}$, $q$, and
$\sigma^{2}$ are available, these can be plugged into the marginal posterior
for $f$ to obtain the so called empirical, marginal posterior for $f$:

$\hat{\Pi}^{\mbox{\boldmath{$f$}}}\\!\big{(}\,\cdot\mid\mbox{\boldmath{$Y$}}\big{)}=\Pi^{\mbox{\boldmath{$f$}}}\\!\big{(}\,\cdot\mid\mbox{\boldmath{$Y$}},\hat{\sigma}^{2},\hat{\lambda},\hat{q},\hat{\mbox{\boldmath{$R$}}\,}\big{)}=\Pi^{\mbox{\boldmath{$f$}}}\\!\big{(}\,\cdot\mid\mbox{\boldmath{$Y$}},\sigma^{2},\lambda,q,\mbox{\boldmath{$R$}}\big{)}\big{|}_{\big{(}\lambda,q,\mbox{\boldmath{$R$}},\sigma^{2}\big{)}=\big{(}\hat{\lambda},\hat{q},\hat{\mbox{\boldmath{$R$}}\,},\hat{\sigma}^{2}\big{)}}.$
(7)

Given $Y$, this is just a
$t_{n+1}\big{(}\hat{\mbox{\boldmath{$f$}}\,\,},\hat{\sigma}^{2}\hat{\mbox{\boldmath{$S$}}\,}\hat{\mbox{\boldmath{$R$}}\,}\big{)}$
distribution, which is centred at the spline estimate
$\hat{\mbox{\boldmath{$f$}}\,}$, and whose covariance matrix depends on the
(random) smoother $\hat{\mbox{\boldmath{$S$}}\,}$ – which is just the smoother
$S$ with $(\hat{\lambda},\hat{q},\hat{\mbox{\boldmath{$R$}}\,})$ plugged in
for $(\lambda,q,\mbox{\boldmath{$R$}})$ – and the estimates $\hat{\sigma}^{2}$
and $\hat{\mbox{\boldmath{$R$}}\,}$. From this we can easily construct a
credible set for the regression function.

If $f$ is distributed according to
$\Pi^{\mbox{\boldmath{$f$}}}\\!(\,\cdot\mid\mbox{\boldmath{$Y$}},\lambda,q,\mbox{\boldmath{$R$}},\sigma^{2})$,
then $\mbox{\boldmath{$f$}}-\hat{\mbox{\boldmath{$f$}}\,\,}$ is distributed
according to
$t_{n+1}(\mbox{\boldmath{$0$}},\sigma^{2}\mbox{\boldmath{$S$}}\mbox{\boldmath{$R$}})$.
This means that
$\|\mbox{\boldmath{$f$}}-\hat{\mbox{\boldmath{$f$}}\,\,}\|_{2}^{2}$ has the
same distribution as the random variable
$\sigma^{2}\mbox{\boldmath{$Z$}}^{T}\mbox{\boldmath{$S$}}\mbox{\boldmath{$R$}}\mbox{\boldmath{$Z$}}/U$,
where $\mbox{\boldmath{$Z$}}\sim
N(\mbox{\boldmath{$0$}},{\mbox{\boldmath{$I$}}_{n}})$, and
$U\sim\mathcal{X}^{2}_{n+1}$, with $Z$ and $U$ independent. From this it
follows that there exits a (known) sequence of quantiles
$s_{n}(\lambda,q,\mbox{\boldmath{$R$}})$ such that for every
$n,\lambda,q,\mbox{\boldmath{$R$}}$, and each $\alpha\in(0,1)$,

$\Pi\Big{\\{}\mbox{\boldmath{$f$}}:\|\mbox{\boldmath{$f$}}-\hat{\mbox{\boldmath{$f$}}\,\,}\|_{2}^{2}\leq\sigma^{2}\,s_{n}(\lambda,q,\mbox{\boldmath{$R$}})\mid\mbox{\boldmath{$Y$}},\lambda,q,\mbox{\boldmath{$R$}},\sigma^{2}\Big{\\}}=1-\alpha.$

We can then define the sets

$\hat{C}_{n}(L)=\Big{\\{}\mbox{\boldmath{$f$}}:\|\mbox{\boldmath{$f$}}-\hat{\mbox{\boldmath{$f$}}\,\,}\|_{2}^{2}\leq
L\,\hat{\sigma}^{2}\,s_{n}(\hat{\lambda},\hat{q},\hat{\mbox{\boldmath{$R$}}\,})\Big{\\}},\qquad
L\geq 1,$ (8)

which (simply by definition of the quantile sequence $s_{n}$) satisfy

$\hat{\Pi}\big{\\{}\hat{C}_{n}(L)\mid\mbox{\boldmath{$Y$}}\big{\\}}\geq
1-\alpha,\qquad L\geq 1,$

and are therefore a credible set – a small, high probability region of the
empirical posterior. These sets can be easily simulated by sampling a large
number of draws from the posterior and keeping the $(1-\alpha)$-fraction of
them that are the closest (in $\ell_{2}$-norm) to
$\hat{\mbox{\boldmath{$f$}}\,\,}$. These will give a good visual
representation of the uncertainty in the estimate of the regression function.
One can of course construct confidence sets for $f$ (or functionals of $f$) in
terms of other norms (such as $\ell_{\infty}$\- norm) by simply adjusting the
quantile sequence $s_{n}$. However, the coverage properties of such sets is
still not so well understood outside the setting of regression with Gaussian
i.i.d. noise, and of the Gaussian white noise model.

The construction above can be found in Szabó et al., (2015), and is also used
in Serra and Krivobokova, (2017), where the behaviour of this set as a
confidence set is studied, under a frequentist assumption on the distribution
of the data. In the latter paper, for the case where
$\mbox{\boldmath{$\mathcal{R}$}}={\mbox{\boldmath{$I$}}_{n}}$, it is shown
(cf. Serra and Krivobokova,, 2017, Theorem 3) that this set has two important
properties if $L$ is taken appropriately large. Firstly, the set contains the
true underlying regression function with probability converging to 1,
uniformly over a large subset of functions; secondly, with probability
converging to 1, the (random) radius of $\hat{C}_{n}(L)$ is of the order of
the minimax risk corresponding to the smoothness class to which the regression
function belongs. In other words, uniformly over a large set of regression
function, this credible set covers the true regression function (honest
coverage), and has a size (radius) that adapts to the smoothness of the
underlying regression function (adaptive coverage).

The asymptotic results from Appendix A show that the presence of correlation
has a relatively simple scaling effect on the smoothing parameter. So, in
principle, the theoretical results of Serra and Krivobokova, (2017) can be
extended to cover the correlated noise case as well, by simply picking larger
values of the multiplier $L$, but this problem will be studied in more
generality elsewhere. In this paper, we focus instead on the implementation,
and numerical aspects of the procedure.

##  3 Numerical simulations

In this section we investigate the small sample performance of our estimation
procedure and compare it to two alternatives. Among approaches that make a
parametric assumption on $\mathcal{R}$ we consider the well-established method
based on splines that is implemented in the statistical software R in the
function gamm of package mgcv, see Wood, (2017). Among approaches that make no
parametric assumption on $\mathcal{R}$ we consider the method by Herrmann et
al., (1992). This kernel based method uses sample autocorrelation estimators
to improve bandwidth selection and is developed under assumptions of
$m$-dependence in the residuals. However, the authors state that the method
still works if the residuals satisfy “some mixing conditions”.

Our simulation set up is as follows. We consider three regression functions

$\displaystyle f_{1}(x)$ $\displaystyle=$
$\displaystyle\sum_{i=3}^{n}\psi_{3,i}(x)\\{\pi(i-1)\\}^{-3}\cos(2i)$
$\displaystyle f_{2}(x)$ $\displaystyle=$
$\displaystyle\sum_{i=4}^{n}\psi_{5,i}(x)\\{\pi(i-1)\\}^{-5}\cos(2i)$
$\displaystyle f_{3}(x)$ $\displaystyle=$ $\displaystyle 2\sin(4\pi
x),\phantom{\sum_{i=4}^{n}}$

where $\psi_{j,i}$, $j\in\\{3,5\\}$ is the $i$-th Demmler-Renisch basis of
${\cal{W}}_{j}$ given explicitly in Appendix B.1, $n=500$ and $x$ is a
sequence of $n$ equally spaced data on $[0,1]$. All three functions are
subsequently scaled to have standard deviation $1$. The standard deviation of
the residuals is taken to be $\sigma=0.33$ to imply a medium signal-to-noise
ratio of $3$. All reported results are based on the Monte Carlo sample
$M=500$. We consider $9$ types of the residual processes: independent and
identically distributed, an $AR(1)$ process with the parameters
$\phi\in\\{-0.9,-0.5,0.5,0.9\\}$, a $MA(1)$ process with the parameters
$\theta\in\\{-0.5,0.5\\}$, an ARMA$(2,2)$ process with $\phi=(0.7,-0.4)$ and
$\theta=(-0.2,0.2)$ and a zero-mean Gaussian process with the correlation
matrix with $(i,j)$ entries $\cos(6.5j)\exp(-|i-j|/20)$.

|
---|---
|
Figure 2: Regression functions (black lines) $f_{1}$ (top left), $f_{2}$ (top
right) and $f_{3}$ (bottom left) with the data (grey lines). The data are
simulated adding a zero-mean Gaussian process noise with the correlation
matrix (first row) shown in the bottom right plot.

The top row and bottom left plot in Figure 2 show all three regression
functions (black lines). The data (grey lines) are simulated using a zero-mean
Gaussian noise process with the correlation matrix
${\mathcal{R}}_{i,j}=\cos(6.5\cdot j)\exp(-|i-j|/20)$. Elements
${\mathcal{R}}_{1,j}$, $j=1,\ldots,200$ are shown in the bottom right plot of
Figure 2. All three functions are to be estimated with the method of Herrmann
et al., (1992) (further denoted by HER), with function gamm (denoted by GAM)
and with our approach (denoted by BAS).

|
---|---
|
Figure 3: Estimating equation $T_{q}$ in the top left and the corresponding
estimator of $f_{1}$ (black line) together with the data (black points) and
confidence intervals (grey area) in the top right. First $50$ elements of the
first row of the correlation matrix in grey with its estimator in black
(bottom left) and the true spectral density in grey with its estimator in
black (bottom right).

Before we summarise the simulation results, we demonstrate how our method
works in practice. The data are simulated as described above with the
regression function $f_{1}\in{\cal W}_{3}$ and ARMA$(2,2)$ dependence in the
residuals. After Algorithm 1 has converged, one checks the estimate for $q$.
The top left plot of Figure 3 shows estimating equation
$T_{q}(\hat{\lambda},q,\hat{\mbox{\boldmath{$R$}}\,})$, which is zero at $q$
which is very close to $3$, as it should for $f_{1}\in{\cal{W}}_{3}$. Next,
the corresponding estimate for $f$ is obtained; it is shown in the top right
plot of Figure 3, together with the data (black dots) and confidence bands
(grey area) constructed as described in Section 2.3. The estimate for
${\mathcal{R}}_{1,j}$, $j=1,\ldots,50$ is shown in the bottom left plot of the
same figure in black, which nearly coincides with the truth shown in grey.
This estimate was reconstructed by the discrete Fourier transform of the
estimate of the spectral density, shown in black in the bottom right plot of
Figure 3. The true spectral density is shown in the same plot in grey. Hence,
our method allows for simultaneous, fully automatic and non-parametric
estimation of $q$, $\lambda$ and ${\mathcal{R}}$. However, in the simulation
below we fixed $q=2$ to ensure comparability of all three estimators, which
all should have the same rate of convergence $n^{-{2/5}}$.

For the method of Herrmann et al., (1992) the parameter $m$ is set according
to method (i) described in that paper (p. 787). Namely, we choose $m$ to be
the largest integer such that $\hat{h}_{m}\geq 6/5\;\hat{h}_{m-1}$ and $m\leq
0.2\sqrt{n}\approx 4.5$ for $n=500$, where $\hat{h}_{m}$ is a selected
bandwidth with the parameter $m$. This parameter $m$ is linked to the
assumption of $m$-dependence in the residuals. In our experiments we noticed
that the influence of $m$ on the estimator of $f$ is not very pronounced, but
it does affect the estimator of the autocovariance quite strongly. There is no
simple data-driven approach to choose $m$ such that both mean and
autocovariance estimators are optimal in some sense. Hence, we expect that HER
performs comparably to BAS in estimation of $f$, but should be inferior in
estimation of $R$. In our implementation we used function glkerns of package
lokern by Eva Herrmann for the estimation of $f$ with the second order kernel
and for the estimation of $f^{{}^{\prime\prime}}$ with the fourth order
kernel.

In the function gamm we used low-rank splines with number of knots $n/4=125$,
B-spline basis of degree $3$, penalisation order $q=2$ and specified the
correlation structure according to the true dependence structure in the
residuals for the first $8$ residual processes. Since the parametric
estimators of $\mathcal{R}$ have faster convergence rates, we expect that GAM
will outperform our non-parametric estimator for $\mathcal{R}$ with BAS. Note
that in practice no information on the true correlation structure is given and
some model validation and selection procedures need to be employed. Since the
$9$th process is a non-parametric one, we set the correlation structure in the
call of gamm to be an $AR(1)$ process, mimicking the real situation in
practice where no information on the residual process is available.

The results are summarised in Table 1. For all dependence structures and all
functions we calculate

$\displaystyle A(\hat{f}_{j})$ $\displaystyle=$
$\displaystyle\frac{1}{Mn}\sum_{k=1}^{n}\sum_{i=1}^{M}\\{f_{j}(x_{k})-\hat{f}_{j,i}(x_{k})\\}^{2}$
$\displaystyle A(\hat{R}_{j})$ $\displaystyle=$
$\displaystyle\frac{1}{Mn}\sum_{k=1}^{n}\sum_{i=1}^{M}\\{{\mathcal{R}}_{j}(k)-\hat{R}_{j,i}(k)\\}^{2},\;\;j=1,2,3,$

Here $\hat{f}_{j,i}$ denotes an estimator of $f_{j}$ in $i$th Monte Carlo run
with the residuals following one of the nine processes. Similarly,
${\mathcal{R}}_{j}(k)$ denotes the $k$-th entry of the first row of one of
nine true residual correlation matrices added to the $j$th regression function
$f_{j}$ and $\hat{R}_{j,i}(k)$ is its estimator in the $i$th Monte Carlo run.

$\bf{f_{1}}$ |
---|---
| $A(\hat{f_{1}})$ | $A(\hat{R}_{1})$ | $\hat{q}=3$
Correlation | BAS | GAM | HER | BAS | GAM | HER |
$i.i.d.$ | 3.187 | 3.182 | 3.840 | 0.000 | 0.000 | 0.000 | 0.870
$AR1(-0.9)$ | 1.170 | 1.130 | 3.889 | 0.874 | 0.177 | 5.885 | 0.890
$AR1(-0.5)$ | 1.732 | 1.674 | 2.924 | 0.033 | 0.010 | 0.081 | 0.940
$AR1(0.5)$ | 10.737 | 10.262 | 9.786 | 0.063 | 0.012 | 0.037 | 0.658
$AR1(0.9)$ | 210.067 | 150.041 | 328.900 | 3.709 | 0.426 | 5.088 | 0.010
$MA1(-0.5)$ | 1.464 | 1.126 | 2.057 | 0.015 | 0.001 | 0.064 | 0.572
$MA1(0.5)$ | 8.836 | 6.251 | 6.059 | 0.040 | 0.001 | 0.037 | 0.638
$AR2MA2$ | 6.130 | 5.897 | 5.598 | 0.038 | 0.032 | 0.047 | 0.802
$GP$ | 2.192 | 5.656 | 16.396 | 1.420 | 15.714 | 5.565 | 0.892
$\bf{f_{2}}$ |
| $A(\hat{f_{2}})$ | $A(\hat{R}_{2})$ | $\hat{q}=5$
Correlation | BAS | GAM | HER | BAS | GAM | HER |
$i.i.d.$ | 4.559 | 4.506 | 6.143 | 0.000 | 0.000 | 0.000 | 0.332
$AR1(-0.9)$ | 2.185 | 1.585 | 2.380 | 0.973 | 0.183 | 5.877 | 0.456
$AR1(-0.5)$ | 2.428 | 2.353 | 1.661 | 0.033 | 0.009 | 0.077 | 0.422
$AR1(0.5)$ | 15.007 | 14.083 | 9.033 | 0.075 | 0.013 | 0.032 | 0.122
$AR1(0.9)$ | 267.605 | 191.323 | 338.288 | 4.260 | 0.572 | 5.067 | 0.000
$MA1(-0.5)$ | 1.882 | 1.695 | 1.252 | 0.015 | 0.004 | 0.065 | 0.166
$MA1(0.5)$ | 12.844 | 8.849 | 5.356 | 0.041 | 0.001 | 0.033 | 0.228
$AR2MA2$ | 8.477 | 8.138 | 4.682 | 0.043 | 0.034 | 0.048 | 0.226
$GP$ | 3.907 | 7.028 | 14.673 | 1.548 | 17.964 | 5.663 | 0.330
$\bf{f_{3}}$ |
| $A(\hat{f_{3}})$ | $A(\hat{R}_{3})$ | $\hat{q}=6$
Correlation | BAS | GAM | HER | BAS | GAM | HER |
$i.i.d.$ | 3.270 | 3.264 | 4.058 | 0.000 | 0.000 | 0.000 | 0.422
$AR1(-0.9)$ | 1.093 | 1.186 | 2.735 | 0.902 | 0.189 | 5.874 | 0.800
$AR1(-0.5)$ | 1.832 | 1.753 | 1.981 | 0.032 | 0.009 | 0.072 | 0.500
$AR1(0.5)$ | 11.051 | 10.330 | 9.217 | 0.080 | 0.013 | 0.033 | 0.024
$AR1(0.9)$ | 215.062 | 224.555 | 325.371 | 3.746 | 0.803 | 5.064 | 0.000
$MA1(-0.5)$ | 1.502 | 1.142 | 1.487 | 0.016 | 0.001 | 0.061 | 0.234
$MA1(0.5)$ | 8.353 | 6.688 | 5.821 | 0.029 | 0.001 | 0.036 | 0.118
$AR2MA2$ | 6.535 | 6.224 | 5.360 | 0.036 | 0.034 | 0.046 | 0.108
$GP$ | 2.495 | 6.008 | 15.125 | 1.461 | 16.095 | 5.662 | 0.312
Table 1: Simulation results. Values of $A(\hat{f_{j}})$ and $A(\hat{R}_{j})$,
$j=1,2,3$ are multiplied by $10^{3}$.

Consider first a parametric dependence structures of the residuals (first
eight types). Non-parametric methods HER and BAS perform similarly for the
estimation of the regression function, with HER performing better for
$MA$-processes and BAS being better for $AR$-processes. As expected,
performance of HER for the autocovariance estimation is uniformly worse
(except for an $AR1(0.5)$ process). The GAM approach performs very much like
BAS in estimating the regression function, but, as expected, outperforms in
covariance matrix estimation. Of course, this comes at the cost of using the
true specification for the dependence process, which is unavailable in
practice. In case of a non-parametric dependence in the residuals (the ninth,
Gaussian process), BAS clearly outperforms both GAM and HER. Overall, BAS
shows very good small sample properties, while being fully automatic and non-
parametric.
The last column of Table 1 shows the proportion of correctly estimated $q$ out
of $M$. In general, $q$ is estimated more reliably for smaller $q$s (for the
analytic function $f_{3}$ we set $q=6$ if
$T_{q}(\hat{\lambda},q,\hat{\mbox{\boldmath{$R$}}\,})$ remains negative as it
should). Also, high autoregressive dependence in the data and high smoothness
of the signal make it more difficult to identify $q$. Note also that $n=500$
is rather moderate and performance of
$T_{q}(\hat{\lambda},q,\hat{\mbox{\boldmath{$R$}}\,})$ becomes better with the
growing sample size.

##  4 Examples

###  4.1 Exchange rates of Russian Ruble

The first example is on the official exchange rates of Russian Ruble to Euro
and US Dollars. The data are freely available from the Central bank of Russian
Federation under https://www.cbr.ru/eng/currency_base/dynamics/. We obtained
the data on prices of Euro and US Dollar in Rubles for the time period between
03 June 2008 and 01 June 2018 and calculated the dual currency basket weighing
US Dollar prices by $0.55$ and Euro prices by $0.45$, i.e.,
$B=0.55\;USD+0.45\;EUR$.

Figure 4: Estimators of the mean (top), first row of the correlation matrix
(first 50 values, bottom left) and the estimating equation for the smoothness
degree $q$ (bottom right) of the dual currency basket. Data are shown in grey,
BAS estimators are the black solid lines, HER estimators are doted lines and
GAM estimators are dashed lines.

The top plot in Figure 4 shows the data in grey. The signal-to-noise ratio of
the data is very high and the number of observations is large: $n=2477$.
Fitting the data with our algorithm delivers $\hat{q}=2$ (bottom right plot of
Figure 4). The corresponding mean estimator is shown as a black bold line on
the top plot of Figure 4 and the estimator of the first row of the correlation
matrix (first $50$ values) is shown on the bottom left plot of Figure 4. Since
$q$ is estimated to $2$ (bottom right plot of Figure 4) one can expect that
both competing methods HER and GAM should deliver similar results. Method HER
indeed delivers the mean estimator visually indistinguishable from the one
obtained with BAS (with $m$ set as described in the previous section). The
dotted line on the bottom left plot shows the corresponding sample
autocorrelation. We observed that the shape of the estimated autocorrelation
appears is more influenced by $m$ than the mean estimator. Recall that the
choice of $m$ is not data driven. Method GAM has been applied assuming an
$AR(1)$ process for the residuals. The estimate of the mean seems to be way
too smooth with the most variation moved to the residuals; the corresponding
estimated autocorrelation decays extremely slowly. The model validation via
residual analysis suggests that the estimators HER and BAS are more
appropriate.

###  4.2 Hourly loads at a US facility

In the second example we consider data from the load forecasting track of the
Global Energy Forecasting Competition 2012
(http://www.drhongtao.com/gefcom/2012). These are data on hourly loads of a US
facility at $20$ zones from the 1st hour of January 1st, 2004 to the 6th hour
of June 30th, 2008. The goal of the competition was to consider each time
series of the $20$ zones, as well as their mean, in order to make a one week
out-of-sample forecast, as well as backcast certain values set to be missed
within the observational period. Here we are not interested in forecasting the
data, but rather understanding their structure. We consider the mean over all
$20$ zones over the whole time period – all together $n=1650$ observations.
Missing values were imputed using R package Hmisc; omitting these missing
values lead to the same estimators and same conclusions.

Figure 5: Estimators of the mean (top), first row of the correlation matrix
(first 50 values, bottom left) and the estimating equation for the smoothness
degree $q$ (bottom right) of the hourly loads of a US facility. Data are shown
in grey, BAS estimators are the black solid lines, GAM estimators are dashed
lines.

Estimating with BAS delivered a positive and increasing function in the
estimating equation for $q$. This suggests that the mean function has either
$1$ continuous derivative or is even less smooth. In any case, this suggests a
fit with $q=1$ to be more appropriate than a fit with $q=2$. This agrees with
the nature of the data, where many peak loads might happen; in the literature
similar data are treated e.g., by a mixture of smoothing splines and wavelets,
where the later pick up the peaks, see Amato et al., (2017). The estimator of
the mean with $q=1$ is shown in the top plot of Figure 5 as black line. The
corresponding autocorrelations are shown as solid black line in the bottom
left top. Since the HER method is defined only for even order kernels, we can
not obtain the fit that would be comparable with the fit by BAS with $q=1$.
Setting $q=1$ in GAM with an $AR(2)$ process for the residuals results in an
estimator of the mean that is very close to the one obtained with BAS, having
only slightly less pronounced peaks (not visible in the plot). The
corresponding autocorrelation estimator is shown as a dashed line in the
bottom left plot of Figure 5.

Estimation using BAS with $q=2$ and by HER method gives estimates of the mean
that are very close to the one obtained by GAM, setting the correlation
structure to an $AR(1)$ process; this is shown in Figure 1, right plot. The
corresponding autocorrelation estimators are also reasonably close. The
residual analysis suggests, however, that the fit with $q=1$ is more
appropriate.

##  5  Conclusions

Correlation is ubiquitous in applications, particularly when data collected
sequentially over time, and ignoring this can have severe consequences for
inference. Covariance is typically taken into account by either making
parametric assumptions on the dependence structure, or by relying on
introducing hyperparameters that then have be set heuristically. We propose a
fully automatic, non-parametric, adaptive method to estimate both mean
function and autocovariance function, and further supply confidence sets for
the mean function that quantify the uncertainty of the estimator. The approach
is implemented in the R package eBsc (available from the authors), and
delivers results in a quick and numerically stable way.

Our method is iterative. We start with an arbitrary guess of the correlation
structure (e.g., independent) and recursively update the our data-driven
choice for the smoothing parameter and correlation structure. We can show
that, under appropriate, mild assumptions, the initial choice of the smoothing
parameter is consistent for an oracle, which is of the correct magnitude, but
not optimal. After further iterations, the choice of the smoothing parameter
is consistent for the oracle that has access to the true covariance structure
of the data. Also our estimate of the covariance structure is consistent. The
order of the splines to be used in the estimator is also estimated from the
data.

Our numerical simulations suggest that when the covariance structure is
parametric our (non-parametric) method is typically outperformed at estimating
the autocovariance function by parametric methods, but still delivers
comparable results for the mean. Under the same situation, compared to other
non-parametric methods we seem to perform better at estimating the covariance
structure. When the covariance structure is non-parametric we strongly
outperform competing methods.

Acknowledgements.

T.K. would like to thank the Isaac Newton Institute for Mathematical Sciences,
Cambridge, for support and hospitality during the programme “Statistical
Scalability”, where work on this paper was undertaken and supported by EPSRC
grant numbers EP/K032208/1 and EP/R014604/1. P.S. carried out part of this
research while he was a postdoctoral researcher at the Institute of
Mathematical Stochastics, Georg-August-Universität Göttingen and would like to
acknowledge together with T.K. the support of the German Research Foundation
(Deutsche Forschungsgemeinschaft) as part of the Institutional Strategy of the
University of Göttingen. The authors are greatly indebted to Dr. Christoph
Lehrenfeld for helpful discussions.

## Appendix

##  Appendix A Asymptotics

In this section we describe the asymptotic behaviour of the estimators from
Section 2. The proofs for the results below are in Appendix C. Henceforth let
${\mathbb{P}}$, (resp. ${\mathbb{E}}$, ${\mathbb{V}}$) represent probability
(resp. expectation, variance) with respect to
$N(\mbox{\boldmath{$f$}},\sigma^{2}\mbox{\boldmath{$\mathcal{R}$}})$, where
$f\in L_{2}$, $\sigma^{2}>0$, $\mathcal{R}$ are the true values of the
parameters of interest which determine the distribution of the data. The
matrix $\mathcal{R}$ denotes the true underlying correlation matrix of the
data, while $R$ denotes the “working” correlation matrix used in the smoother
$S$. We assume thought that
$\mbox{\boldmath{$\mathcal{R}$}}\in\mathcal{M}_{n,\delta}$, so that the
eigenvalues of $\mathcal{R}$ are bounded away from zero and from infinity. The
numerical procedure from Section 2.2 ensures that this is also the case for
our estimates of the spectral density of the noise.

Our first result describes the uniform convergence (over $\lambda$ and $R$) of
the criterium $T_{\lambda}$ to it its expectation.

######  Proposition 1 (Uniform convergence of $T_{\lambda}$)

Assume that $\Lambda_{n}$ (a collection such that $\lambda\in\Lambda_{n}$
satisfies $0<\lambda<1$, $\lambda\to 0$, and $n\lambda\to\infty$ as
$n\to\infty$), is such that
$|\Lambda_{n}|\exp\\{-\inf_{\lambda\in\Lambda_{n}}\lambda^{-1/(2q)}\\}=o(1)$.
Then,

$\sup_{\mbox{\boldmath{$R$}}\in\mathcal{M}_{n,\delta}}\mathbb{P}\bigg{\\{}\sup_{\lambda\in\Lambda_{n}}\big{|}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})-\mathbb{E}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})\big{|}>\epsilon\bigg{\\}}\to
0,\quad\epsilon>0,\quad n\to\infty.$ (9)

A consequence of this result is that, for any working correlation
$\mbox{\boldmath{$R$}}\in\mathcal{M}_{n,\delta}$, the solution to the
estimating equation $T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})=0$ converges
to the solution of $\mathbb{E}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})=0$.
In particular we have the following proposition.

######  Proposition 2 (Consistency of the preliminary estimate of $\lambda$)

Let $f\in\mathcal{W}_{\beta}(M)$, $\beta>1/2$, and assume that
$\|f^{(\beta)}\|^{2}>0$. Assume that $\sigma^{2}>0$. Assume that the first row
of $\mbox{\boldmath{$\mathcal{R}$}}\in\mathcal{M}_{n,\delta}$ is absolutely
summable, and denote by $\varrho$ the underlying spectral density. Let
$\Lambda_{n}$ be as in Proposition 1. Assume $q>1/2$ is such that
$\kappa_{q}({\mbox{\boldmath{$I$}}_{n}})$, as defined in (24), is stricly
positive. (Note that this constant depends on $\varrho$.) Denote by
$\lambda_{q,{\mbox{\boldmath{$I$}}}}$ the solution to
$\mathbb{E}T_{\lambda}(\lambda,q,{\mbox{\boldmath{$I$}}_{n}})=0$, $\lambda>0$,
and by $\hat{\lambda}_{q,{\mbox{\boldmath{$I$}}}}$ any solution to
$T_{\lambda}(\lambda,q,{\mbox{\boldmath{$I$}}_{n}})=o(1)$,
$\lambda\in\Lambda_{n}$.

Under these assumptions, $\hat{\lambda}_{q,{\mbox{\boldmath{$I$}}}}$ is
consistent for the oracle $\lambda_{q,{\mbox{\boldmath{$I$}}}}$, in that
$\hat{\lambda}_{q,{\mbox{\boldmath{$I$}}}}/\lambda_{q,{\mbox{\boldmath{$I$}}}}\to
1$, in $\mathbb{P}$-probability, as $n\to\infty$.

As for the oracle, if $q\leq\max\\{\beta>1/2:f\in\mathcal{W}_{\beta}(M)\\}$,
then

$\lambda_{q,{\mbox{\boldmath{$I$}}}}=\bigg{[}\frac{n\|f^{(q)}\|^{2}}{\sigma^{2}\,\kappa_{q}({\mbox{\boldmath{$I$}}_{n}})}\\{1+o(1)\\}\bigg{]}^{-\frac{2q}{2q+1}},$
(10)

and, if $f\in\mathcal{W}_{\beta}$, $\beta>1/2$, and $q>\beta$, then

$\lambda_{q,\mbox{\boldmath{$I$}}}\geq\bigg{[}\frac{n\|f^{(\beta)}\|^{2}}{\sigma^{2}\,\kappa_{q}({\mbox{\boldmath{$I$}}_{n}})}\\{1+o(1)\\}\bigg{]}^{-\frac{2q}{2\beta+1}}.$
(11)

If the noise is white, then $\varrho=1$ in which case the constant in the
denominator becomes $\kappa_{q}(0,2)$ and one just recuperates the oracles for
the i.i.d. noise setting; cf. Serra and Krivobokova, (2017).

The constant $\kappa_{q}(\mbox{\boldmath{$I$}}_{n})$ will not match the
constant $\kappa_{q}(\mbox{\boldmath{$\mathcal{R}$}})$ from the oracle that we
would get if $\mathcal{R}$ were known (this is the oracle from Theorem 1).
Comparing the two constants one can characterise which types of noise
processes lead to under-smoothing and which lead to over-smoothing if the
correlation is ignored. Note also that the assumption that
$\kappa_{q}(\mbox{\boldmath{$I$}}_{n})$ has to be positive means that if the
correlation structure is ignored, then the type of correlation in the noise
introduces constraints on the values that $q$ can take. To the best of our
knowledge this has never been reported in the literature and justifies the
need to pick $q$ in a data driven way if the correlation structure is unknown.
(The constant $\kappa_{q}(\mbox{\boldmath{$\mathcal{R}$}})$ is always
positive.)

Although the working correlation is not the true correlation of the data, the
oracle is of the correct order and will lead to a risk for the resulting
spline estimator of the correct (minimax) order (compare the oracle from
Proposition 2 and the one from Theorem 1). Finite dimensional performance,
though, leaves much to be desired, so the next result describes the behaviour
of the estimator $\hat{\mbox{\boldmath{$R$}}\,}$ when $\lambda$ is set to
$\hat{\lambda}_{{\mbox{\boldmath{$I$}}}}$ and $R$ is set to
${\mbox{\boldmath{$I$}}_{n}}$.

######  Proposition 3 (Consistency of the preliminary estimator of
$\mathcal{R}$)

Let the coordinates of $\tilde{\mbox{\boldmath{$\rho$}}\,}$ solve
$T_{\rho_{i}}(\hat{\lambda}_{q,\mbox{\boldmath{$I$}}},q,1)=0$, and let
$\bar{\mbox{\boldmath{$S$}}\,}=\mbox{\boldmath{$S$}}_{\xi,p,\mbox{\boldmath{$I$}}}$.
Let $\hat{\mbox{\boldmath{$\rho$}}\,}$ be
$\bar{\mbox{\boldmath{$S$}}\,}\tilde{\mbox{\boldmath{$\rho$}}\,}$ with entries
first truncated to $[\delta,1/\delta]$ and then scaled to add up to $n$.
Further, assume that the entries of $\mathcal{R}$ satisfy
$\mathcal{R}_{1,i}=O(|i|^{-2\alpha})$, $\alpha>1/2$, so that in particular
$\sum_{i=1}^{n}|\mathcal{R}_{1,i}|$ is finite. Then,

$\mathbb{E}\|\hat{\sigma}^{2}\hat{\mbox{\boldmath{$\rho$}}\,}-\sigma^{2}\mbox{\boldmath{$\varrho$}}\|_{\infty}=o(1).$

We do now establish how subsequent iterates of the estimate of $\rho$ behave,
but replacing the constant function $1$ with a consistent estimate of the
spectral sensity should further improve the estimator.

This result ensures, in particular, that for any sequences
$\xi,\lambda\in\Lambda_{n}$, the mean squared error of the estimate of the
spectral density of the noise process is consistent in operator norm. The
smoothing parameter $\xi$ is picked using the estimating equation
$T_{\xi}(\xi,p,1)$.

The following theorem specifies the behaviour of the estimator for $\lambda$
when $R$ is set to any estimator $\hat{\mbox{\boldmath{$R$}}\,}$ of
$\mathcal{R}$ that is consistent (e.g., the estimator from Proposition 3).

######  Theorem 1 (Consistency of estimates of $\lambda$)

Let $f\in\mathcal{W}_{\beta}(M)$, $\beta>1/2$, and assume that
$\|f^{(\beta)}\|^{2}>0$. Assume that $\sigma^{2}>0$. Assume that the first row
of $\mathcal{R}$ is absolutely summable, and denote by $\varrho$ the
underlying spectral density. Let $\hat{\mbox{\boldmath{$R$}}\,}$ be any a.s.
positive definite estimator for $\mathcal{R}$ that is consistent in operator
norm.

Consider the constant $\kappa_{q}(\mbox{\boldmath{$\mathcal{R}$}})$ defined in
(27). Assume also that $\Lambda_{n}$, as specified before is such that
$|\Lambda_{n}|\exp\\{-\inf_{\lambda\in\Lambda_{n}}\lambda^{-1/(2q)}\\}=o(1)$.
For any $q>1/2$ denote by $\hat{\lambda}_{q}$ the solution to
$T_{\lambda}(\lambda,q,\hat{\mbox{\boldmath{$R$}}\,})=o(1)$,
$\lambda\in\Lambda_{n}$, and by $\lambda_{q,\mbox{\boldmath{$\mathcal{R}$}}}$
the solution to
$\mathbb{E}T_{\lambda}(\lambda,q,\mbox{\boldmath{$\mathcal{R}$}})=0$,
$\lambda>0$.

Under these assumptions, $\hat{\lambda}_{q}$ is consistent for the oracle
$\lambda_{q,\mbox{\boldmath{$\mathcal{R}$}}}$, in that
$\hat{\lambda}_{q}/\lambda_{q,\mbox{\boldmath{$\mathcal{R}$}}}\to 1$, in
$\mathbb{P}$-probability, as $n\to\infty$.

As for the value of the oracle, if
$q\leq\max\\{\beta>1/2:f\in\mathcal{W}_{\beta}(M)\\}$, then

$\lambda_{q,\mbox{\boldmath{$\mathcal{R}$}}}=\bigg{[}\frac{n\|f^{(q)}\|^{2}}{\sigma^{2}\,\kappa_{q}(\mbox{\boldmath{$\mathcal{R}$}})}\\{1+o(1)\\}\bigg{]}^{-\frac{2q}{2q+1}},$
(12)

and, if $f\in\mathcal{W}_{\beta}$, $\beta>1/2$, and $q>\beta$, then

$\lambda_{q,\mbox{\boldmath{$\mathcal{R}$}}}\geq\bigg{[}\frac{n\|f^{(\beta)}\|^{2}}{\sigma^{2}\,\kappa_{q}(\mbox{\boldmath{$\mathcal{R}$}})}\\{1+o(1)\\}\bigg{]}^{-\frac{2q}{2\beta+1}}.$
(13)

Note that these oracles match the oracles for $\lambda$ when $\mathcal{R}$ is
known.

The oracle provided by this result is the same as if $\mathcal{R}$ were known.
This means that the resulting spline estimator, with data driven choice for
$\lambda$, will attain the same rate as if $\mathcal{R}$ were known. In turn,
since the noise is short range dependent, this corresponds to the same rate as
one would obtain in the i.i.d. setting for estimating a $\beta$-smooth signal;
cf. Section 3.3 of Yang et al., (2001).

Our final theoretical result addresses the behaviour of the estimator for $q$.

######  Theorem 2 (Consistency of estimate of $q$)

Assume that the conditions of Theorem 1 are met. Assume in addition that the
set $Q_{n}$ from the definition of $\hat{q}$ is such that

$|Q_{n}|\exp\Big{\\{}-\inf_{q\in Q_{n}}\lambda_{q}^{-1/(2q)}\Big{\\}}=o(1),$

where $\lambda_{q}$ is the oracle from the previous theorem. Let
$\mathcal{M}=\mathcal{M}(L,N,\rho)$, $L>0$, $N\in\mathbb{N}$, $s\geq 2$,
corresponds the set of all square integrable sequences $b_{q,1},b_{q,2},\dots$
such that

$\frac{1}{n}\sum_{i=j}^{n}b_{q,i}^{2}\leq\frac{L}{n}\sum_{i=j}^{s\,j}b_{q,i}^{2},\quad
N\leq j\leq n/s,$ (14)

and define
$\bar{\beta}=\max\big{\\{}\beta>1/2:f\in\mathcal{W}_{\beta}(M)\cap\mathcal{M}\big{\\}}$.

If for some $\beta>1/2$, $f\in\mathcal{W}_{\beta}(M)$, then

$\mathbb{P}\big{\\{}\beta<\hat{q}\leq\log(n)\big{\\}}\to 1,\quad n\to\infty.$

If furthermore $\beta=\bar{\beta}$, and $f\in\mathcal{M}$, then

$\hat{q}\stackrel{{\scriptstyle P}}{{\longrightarrow}}\bar{\beta},\quad
n\to\infty,$

If for some $d\in\mathbb{N}$, $f\in\mathcal{P}_{d}$, the space of polynomials
of degree strictly smaller than $d$, then

$\mathbb{P}\big{\\{}d\leq\hat{q}\leq\log(n)\big{\\}}\to 1,\quad n\to\infty.$

##  Appendix B Auxiliary results

In this appendix we collect a number of results that are used throughout this
paper.

###  B.1 Demmler-Reinsch basis

Let $\\{\psi_{i}\\}_{i=1}^{\infty}$ denote the Demmler-Reinsch basis of
${\cal{W}}_{\beta}(M)$, such that

$\int_{0}^{1}\psi_{\beta,i}(x)\psi_{\beta,j}(x)dx=\delta_{ij}=\nu_{\beta,i}^{-1}\int_{0}^{1}\psi_{\beta,i}^{(\beta)}(x)\psi_{\beta,j}^{(\beta)}(x)dx.$

Rosales Marticorena, (2016) found explicit expressions for $\psi_{\beta,i}$
and $\nu_{\beta,i}$ as a solution to

$\displaystyle(-1)^{q}\psi_{\beta,i}^{(2\beta)}=\nu_{\beta,i}\psi_{\beta,i}$
$\displaystyle\psi_{\beta,i}^{(l)}(0)=\psi_{\beta,i}^{(l)}(1)=0,\;\;l=\beta,\beta+1,\ldots,2\beta-1.$

In particular, $\nu_{\beta,1}=\cdots=\nu_{\beta,\beta}=0$ and

$\displaystyle\nu_{\beta,i}$ $\displaystyle=$
$\displaystyle\left\\{\pi\left(i-\frac{\beta+1}{2}\right)\right\\}^{2\beta},\;\;i=\beta+1,\beta+2,\ldots$
$\displaystyle\psi_{\beta,i}(x)$ $\displaystyle=$
$\displaystyle\sqrt{2}\left[\cos\left\\{\pi\left(i-\frac{\beta+1}{2}\right)x+\pi\,\frac{\beta-1}{4}\right\\}+T_{i}(x)\right]$
(15)

where

$T_{i}(x)=\sum_{a_{j}\in
S(\beta)}r_{j}\left[\exp\left\\{-a_{j}\pi\left(i-\frac{\beta+1}{2}\right)x\right\\}+(-1)^{i+1}\exp\left\\{-a_{j}\pi\left(i-\frac{\beta+1}{2}\right)(1-x)\right\\}\right],$

for
$S(\beta)=\cup_{j}\left\\{(-1)^{j/(2\beta)},\overline{(-1)^{j/(2\beta)}}\right\\}$,
with $0\leq j\leq\beta-2$ taking odd values for $\beta$ odd and even values
for $\beta$ even. Constants $r_{j}$ are known and depend on $\beta$ only. Note
that $T_{i}(x)$ vanish exponentially fast away from the boundaries.
The Demmler-Reinsch basis of the natural spline space of degree $2q-1$ with
knots at observations ${\cal{S}}_{2q-1}(\mbox{\boldmath{$x$}})$ is uniquely
defined via

$\sum_{k=1}^{n}\phi_{q,i}(x_{k})\phi_{q,j}(x_{k})=\delta_{ij}=\eta_{q,i}^{-1}\int_{0}^{1}\phi_{q,i}^{(q)}(x)\phi_{q,j}^{(q)}(x)dx$

and
$\mbox{\boldmath{$\Phi$}}_{q}=\mbox{\boldmath{$\Phi$}}_{q}(\mbox{\boldmath{$x$}})=[\phi_{q,1}(\mbox{\boldmath{$x$}}),\dots,\phi_{q,n}(\mbox{\boldmath{$x$}})]=[\phi_{q,j}(x_{i})]_{i,j=1}^{n}$
is the corresponding basis matrix.
Utreras Diaz, (1980) used the results of Fix, (1972) to show that
$|n\eta_{q,i}-\nu_{q,i}|=O(n^{-2})$. From Fix, (1972) and Fix, (1973) also
follows that $\|\sqrt{n}\phi_{q,i}-\psi_{q,i}\|_{{\cal{W}}_{q}}=O(n^{-1})$,
or, equivalently, that
$\|\phi_{q,i}-\psi_{q,i}/\sqrt{n}\|_{L_{2}}=O(n^{-3/2})$.
To apply results by Fix, (1972) and Fix, (1973) the standard result (see Lemma
3.2 in Utreras Diaz,, 1980) is used.

###### Lemma 1

Let $f\in W_{\beta}(M)$, $\beta\geq 2$, then

$\left|\frac{1}{n}\sum_{i=1}^{n}f(x_{i})-\int_{0}^{1}f(x)dx\right|\leq\frac{c}{n^{2}}\|f\|_{{\cal{W}}_{\beta}},$

where $x_{i}=(i-1)/(n-1)$, $i=1,\ldots,n$.

###  B.2 Matrix identitites

We derive some matrix identities that are needed in our derivations. The
smoother matrix $S$ can also be written as

$\mbox{\boldmath{$S$}}=\mbox{\boldmath{$R$}}^{1/2}\left\\{\mbox{\boldmath{$I$}}_{n}+\lambda\mbox{\boldmath{$R$}}^{1/2}\mbox{\boldmath{$\Phi$}}\mbox{diag}(n\mbox{\boldmath{$\eta$}}_{q})\mbox{\boldmath{$\Phi$}}^{T}\mbox{\boldmath{$R$}}^{1/2}\right\\}^{-1}\mbox{\boldmath{$R$}}^{-1/2}.$

From this expression is it clear that
$\mbox{\boldmath{$R$}}^{-1}\mbox{\boldmath{$S$}}=\mbox{\boldmath{$S$}}^{T}\mbox{\boldmath{$R$}}^{-1}$
so that in particular

$\mbox{\boldmath{$R$}}^{-1}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})=({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})^{T}\mbox{\boldmath{$R$}}^{-1}.$

From the definition of $S$ is is also simple to check the scaling relation

$\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$S$}}^{-1}-{\mbox{\boldmath{$I$}}_{n}})=\lambda{\mbox{\boldmath{$\Phi$}}_{q}}\mbox{diag}(n\mbox{\boldmath{$\eta$}}_{q}){\mbox{\boldmath{$\Phi$}}_{q}}^{T}=\mbox{\boldmath{$S$}}_{\mbox{\boldmath{$I$}}}^{-1}-\mbox{\boldmath{$I$}}_{n}.$

The estimating equation for $\lambda$ is obtained in a strightforward way by
noting that since
$\partial\mbox{\boldmath{$S$}}/\partial\lambda=-({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$S$}}/\lambda$,
then

$\frac{\partial\mbox{\boldmath{$R$}}^{-1}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})}{\partial\lambda}=\frac{1}{\lambda}\mbox{\boldmath{$R$}}^{-1}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$S$}},$

and the estimating equation for $q$ is based on the fact that

$\displaystyle\frac{\partial\mbox{\boldmath{$R$}}^{-1}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})}{\partial
q}$
$\displaystyle=\mbox{\boldmath{$R$}}^{-1}\mbox{\boldmath{$S$}}\frac{\partial\mbox{\boldmath{$S$}}^{-1}}{\partial
q}\mbox{\boldmath{$S$}}=\frac{\lambda}{q}\mbox{\boldmath{$R$}}^{-1}\mbox{\boldmath{$S$}}\mbox{\boldmath{$R$}}\mbox{\boldmath{$\Phi$}}_{q}\mbox{diag}\big{\\{}n\mbox{\boldmath{$\eta$}}_{q}\circ\log(n\mbox{\boldmath{$\eta$}}_{q})\big{\\}}\mbox{\boldmath{$\Phi$}}_{q}^{T}\mbox{\boldmath{$S$}}$
$\displaystyle=\frac{\lambda}{q}\mbox{\boldmath{$S$}}^{T}\mbox{\boldmath{$\Phi$}}_{q}\mbox{diag}\big{\\{}n\mbox{\boldmath{$\eta$}}_{q}\circ\log(n\mbox{\boldmath{$\eta$}}_{q})\big{\\}}\mbox{\boldmath{$\Phi$}}_{q}^{T}\mbox{\boldmath{$S$}}$
$\displaystyle=\frac{1}{q}\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$I$}}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$\Phi$}}_{q}\mbox{diag}\big{\\{}\log(n\mbox{\boldmath{$\eta$}}_{q})\big{\\}}\mbox{\boldmath{$\Phi$}}_{q}^{T}\mbox{\boldmath{$S$}},$

where $\circ$ denotes the Hadamard product. For the derivation of the
estimating equations for the $\rho_{i}$ we use the fact that

$\frac{\partial\mbox{\boldmath{$S$}}}{\partial\rho_{i}}=-\mbox{\boldmath{$S$}}\frac{\partial\mbox{\boldmath{$S$}}^{-1}}{\partial\rho_{i}}\mbox{\boldmath{$S$}}=-\mbox{\boldmath{$S$}}\frac{\partial\mbox{\boldmath{$R$}}}{\partial\rho_{i}}(\mbox{\boldmath{$S$}}_{{\mbox{\boldmath{$I$}}_{n}}}^{-1}-{\mbox{\boldmath{$I$}}_{n}})\mbox{\boldmath{$S$}}=-\mbox{\boldmath{$S$}}\frac{\partial\mbox{\boldmath{$R$}}}{\partial\rho_{i}}\mbox{\boldmath{$R$}}^{-1}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}}),$

such that, combining the above, we see that by using the chain rule

$\displaystyle\frac{\partial\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})}{\partial\rho_{i}}=-\mbox{\boldmath{$R$}}^{-1}\frac{\partial\mbox{\boldmath{$R$}}}{\partial\rho_{i}}\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})+\mbox{\boldmath{$R$}}^{-1}\frac{\partial\mbox{\boldmath{$S$}}}{\partial\rho_{i}}$
$\displaystyle\quad=-\mbox{\boldmath{$R$}}^{-1}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})\frac{\partial\mbox{\boldmath{$R$}}}{\partial\rho_{i}}\mbox{\boldmath{$R$}}^{-1}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})=-({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})^{T}\mbox{\boldmath{$R$}}^{-1}\frac{\partial\mbox{\boldmath{$R$}}}{\partial\rho_{i}}\mbox{\boldmath{$R$}}^{-1}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}}).$

###  B.3 Traces and quadratic forms

We compute some traces and quadratic forms involving the smoother matrix $S$.
Let $\rho$ be the spectral density corresponding to $R$.

Define the matrix

$\mbox{\boldmath{$S$}}^{*}=\mbox{\boldmath{$S$}}^{*}_{\lambda,q,\mbox{\boldmath{$R$}}}={\mbox{\boldmath{$\Phi$}}_{q}}\left\\{{\mbox{\boldmath{$I$}}_{n}}+\lambda\,\mbox{diag}(n\mbox{\boldmath{$\eta$}}_{q}\circ\mbox{\boldmath{$\rho$}})\right\\}^{-1}{\mbox{\boldmath{$\Phi$}}_{q}}^{T},$
(16)

where $\mbox{\boldmath{$\rho$}}=(\rho_{1},\dots,\rho_{n})$, for
$\rho_{i}=\rho(\pi t_{i})$, $i=1,\ldots,n$, and $\circ$ represents the
Hadamard product. From the definition of $S$ and $\mbox{\boldmath{$S$}}^{*}$,

$\displaystyle\mbox{\boldmath{$S$}}$
$\displaystyle=\left[{\mbox{\boldmath{$I$}}_{n}}+\lambda\mbox{\boldmath{$R$}}{\mbox{\boldmath{$\Phi$}}_{q}}\mbox{diag}(\mbox{\boldmath{$\rho$}})^{-1}\mbox{diag}(n\mbox{\boldmath{$\eta$}}_{q}\mbox{\boldmath{$\rho$}}){\mbox{\boldmath{$\Phi$}}_{q}}^{T}\right]^{-1}=\left[{\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$R$}}{\mbox{\boldmath{$\Phi$}}_{q}}\mbox{diag}(\mbox{\boldmath{$\rho$}})^{-1}{\mbox{\boldmath{$\Phi$}}_{q}}^{T}(\mbox{\boldmath{$S$}}^{*-1}-{\mbox{\boldmath{$I$}}_{n}})\right]^{-1}$
$\displaystyle=\mbox{\boldmath{$S$}}^{*}\left[{\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$S$}}^{*}-{\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$R$}}{\mbox{\boldmath{$\Phi$}}_{q}}\mbox{diag}(\mbox{\boldmath{$\rho$}})^{-1}{\mbox{\boldmath{$\Phi$}}_{q}}^{T}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}}^{*})\right]^{-1}=\mbox{\boldmath{$S$}}^{*}\left[{\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$\Delta$}}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}}^{*})\right]^{-1},$

where the perturbation matrix $\Delta$ is given by

$\mbox{\boldmath{$\Delta$}}=\mbox{\boldmath{$\Delta$}}(\mbox{\boldmath{$R$}})=\mbox{\boldmath{$R$}}{\mbox{\boldmath{$\Phi$}}_{q}}\mbox{diag}(\mbox{\boldmath{$\rho$}})^{-1}{\mbox{\boldmath{$\Phi$}}_{q}}^{T}-{\mbox{\boldmath{$I$}}_{n}}.$
(17)

With this notation, the relation between $S$ and $\mbox{\boldmath{$S$}}^{*}$
is equivalent to

$\mbox{\boldmath{$S$}}({\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$\Delta$}})=({\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$S\Delta$}})\mbox{\boldmath{$S$}}^{*}.$
(18)

(Note that $S$, $\mbox{\boldmath{$S$}}^{*}$, and
${\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$\Delta$}}$ are positive
definite, so that also
${\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$S$}}\Delta$ is positive
definite.)

Let $M$ be an arbitrary, positive definite $n\times n$ matrix with bounded
singular values. Our goal is to compute traces and quadratic forms involving
matrices of the type
$\mbox{\boldmath{$M$}}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})^{m}\mbox{\boldmath{$S$}}^{l}$,
$m\in\mathbb{N}\cup\\{0\\}$. We first show that if all but a finite number of
singular values of $\Delta$ are $o(1)$, then:

$\mbox{tr}\big{[}\mbox{\boldmath{$M$}}\\{\mbox{\boldmath{$S$}}({\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$\Delta$}})\\}^{m}\big{]}=\mbox{tr}\big{[}\mbox{\boldmath{$S$}}^{r}\mbox{\boldmath{$M$}}\\{\mbox{\boldmath{$S$}}({\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$\Delta$}})\\}^{m-r}\big{]}\\{1+o(1)\\},\quad
r=0,1,\dots,m.$

The equality in the previous display is trivial for $r=0$. Let us then assume
that equality above holds for a certain $r=p\geq 0$ and show that it also
holds for $r=p+1$, so that the previous display follows by induction. By Von
Neumann’s trace inequality,

$\displaystyle\Big{|}\mbox{tr}\big{[}\mbox{\boldmath{$S$}}^{p}\mbox{\boldmath{$M$}}\\{\mbox{\boldmath{$S$}}({\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$\Delta$}})\\}^{m-p}\big{]}-\mbox{tr}\big{[}\mbox{\boldmath{$S$}}^{p+1}\mbox{\boldmath{$M$}}\\{\mbox{\boldmath{$S$}}({\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$\Delta$}})\\}^{m-(p+1)}\big{]}\Big{|}=$
$\displaystyle=$
$\displaystyle\Big{|}\mbox{tr}\big{[}\mbox{\boldmath{$S$}}^{p}\mbox{\boldmath{$M$}}\\{\mbox{\boldmath{$S$}}({\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$\Delta$}})\\}^{m-(p+1)}\mbox{\boldmath{$S$}}\mbox{\boldmath{$\Delta$}}\big{]}\Big{|}\leq\mbox{tr}\big{[}\mbox{\boldmath{$S$}}^{p+1}\mbox{\boldmath{$M$}}\\{\mbox{\boldmath{$S$}}({\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$\Delta$}})\\}^{m-(p+1)}\big{]}o(1)+O(1).$

The singular values of all the matrices above are bounded, so we conclude that
as long as the trace above converges to infinity, then

$\mbox{tr}\big{[}\mbox{\boldmath{$S$}}^{p}\mbox{\boldmath{$M$}}\\{\mbox{\boldmath{$S$}}({\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$\Delta$}})\\}^{m-p}\big{]}=\mbox{tr}\big{[}\mbox{\boldmath{$S$}}^{(p+1)}\mbox{\boldmath{$M$}}\\{\mbox{\boldmath{$S$}}({\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$\Delta$}})\\}^{m-(p+1)}\big{]}\\{1+o(1)\\},$

so that
$\mbox{tr}\big{[}\mbox{\boldmath{$M$}}\\{\mbox{\boldmath{$S$}}({\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$\Delta$}})\\}^{m}\big{]}=\mbox{tr}\big{[}\mbox{\boldmath{$S$}}^{p+1}\mbox{\boldmath{$M$}}\\{\mbox{\boldmath{$S$}}({\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$\Delta$}})\\}^{m-(p+1)}\big{]}\\{1+o(1)\\}$.

By taking $r=m$ above, and using the identity (18),

$\mbox{tr}\big{[}\\{({\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$S$}}\mbox{\boldmath{$\Delta$}})\mbox{\boldmath{$S$}}^{*}\\}^{m}\mbox{\boldmath{$M$}}\big{]}=\mbox{tr}\big{[}\mbox{\boldmath{$M$}}\\{\mbox{\boldmath{$S$}}({\mbox{\boldmath{$I$}}_{n}}+\mbox{\boldmath{$\Delta$}})\\}^{m}\big{]}=\mbox{tr}\big{[}\mbox{\boldmath{$M$}}\mbox{\boldmath{$S$}}^{m}\big{]}\\{1+o(1)\\}.$

Using exactly the same argument with $\Delta$ replaced with $S$$\Delta$, and
using the fact that the singular values of $S$ belong to $(0,1]$, we conclude
that the left-hand-side of the previous display is equal to
$\mbox{tr}\big{\\{}\mbox{\boldmath{$M$}}(\mbox{\boldmath{$S$}}^{*})^{m}\big{\\}}\\{1+o(1)\\}$,
so that

$\mbox{tr}\big{(}\mbox{\boldmath{$M$}}\mbox{\boldmath{$S$}}^{m}\big{)}=\mbox{tr}\big{\\{}\mbox{\boldmath{$M$}}(\mbox{\boldmath{$S$}}^{*})^{m}\big{\\}}\\{1+o(1)\\}.$

(Note that the $o(1)$ term can be taken uniform in $m$.)

Using the previous display and the linearity of the trace,

$\displaystyle\mbox{tr}\\{\mbox{\boldmath{$M$}}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})^{m}\mbox{\boldmath{$S$}}^{l}\\}=\sum_{k=0}^{m}(-1)^{k}{m\choose
k}\mbox{tr}\big{(}\mbox{\boldmath{$M$}}\mbox{\boldmath{$S$}}^{k+l}\big{)}=$
$\displaystyle=\sum_{k=0}^{m}(-1)^{k}{m\choose
k}\mbox{tr}\big{\\{}\mbox{\boldmath{$M$}}(\mbox{\boldmath{$S$}}^{*})^{k+l}\big{\\}}\\{1+o(1)\\}=\mbox{tr}\\{\mbox{\boldmath{$M$}}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}}^{*})^{m}(\mbox{\boldmath{$S$}}^{*})^{l}\\}\\{1+o(1)\\}.$

Let the superscript ∘ represent the Hadamard power. Define the matrices

$\mbox{\boldmath{$\Delta$}}_{t,s}=\mbox{\boldmath{$\Delta$}}_{t,s}(\mbox{\boldmath{$\mathcal{R}$}},\mbox{\boldmath{$R$}})=\mbox{\boldmath{$\mathcal{R}$}}^{t}\mbox{\boldmath{$R$}}^{s}{\mbox{\boldmath{$\Phi$}}_{q}}\mbox{diag}(\mbox{\boldmath{$\varrho$}}^{\circ
t}\circ\mbox{\boldmath{$\rho$}}^{\circ
s})^{-1}{\mbox{\boldmath{$\Phi$}}_{q}}^{T}-{\mbox{\boldmath{$I$}}_{n}},$ (19)

where $\mathcal{R}$ is a Toeplitz correlation matrix with underlying spectral
density $\varrho$, and $\varrho$ is defined in analogy to $\rho$. Note that in
particular
$\mbox{\boldmath{$\Delta$}}(\mbox{\boldmath{$R$}})=\mbox{\boldmath{$\Delta$}}_{0,1}({\mbox{\boldmath{$I$}}_{n}},\mbox{\boldmath{$R$}})$.
We can also write
$\mbox{\boldmath{$\mathcal{R}$}}^{t}\mbox{\boldmath{$R$}}^{s}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})^{m}\mbox{\boldmath{$S$}}^{l}=(\mbox{\boldmath{$\Delta$}}_{t,s}+{\mbox{\boldmath{$I$}}_{n}}){\mbox{\boldmath{$\Phi$}}_{q}}\mbox{diag}(\mbox{\boldmath{$\varrho$}}^{\circ
t}\circ\mbox{\boldmath{$\rho$}}^{\circ
s}){\mbox{\boldmath{$\Phi$}}_{q}}^{T}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})^{m}\mbox{\boldmath{$S$}}^{l}$,
so that if we assume that all but a finite number of singular values of
$\mbox{\boldmath{$\Delta$}}_{t,s}$ are $o(1)$, we can use the argument from
before to get rid of the $\mbox{\boldmath{$\Delta$}}_{t,s}$ perturbation;
setting
$\mbox{\boldmath{$M$}}={\mbox{\boldmath{$\Phi$}}_{q}}\mbox{diag}(\mbox{\boldmath{$\varrho$}}^{\circ
t}\circ\mbox{\boldmath{$\rho$}}^{\circ s}){\mbox{\boldmath{$\Phi$}}_{q}}^{T}$
we get

$\mbox{tr}\\{\mbox{\boldmath{$\mathcal{R}$}}^{t}\mbox{\boldmath{$R$}}^{s}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})^{m}\mbox{\boldmath{$S$}}^{l}\\}=\sum_{i=q+1}^{n}\frac{\varrho_{i}^{t}\rho_{i}^{m+s}(\lambda
n\eta_{q,i})^{m}}{\big{(}1+\lambda
n\eta_{q,i}\rho_{i}\big{)}^{m+l}}\\{1+o(1)\\}.$

In the same way, if we now set
$\mbox{\boldmath{$M$}}={\mbox{\boldmath{$\Phi$}}_{q}}\mbox{diag}(\mbox{\boldmath{$\varrho$}}^{\circ
t}\circ\mbox{\boldmath{$\rho$}}^{\circ
s}){\mbox{\boldmath{$\Phi$}}_{q}}^{T}\mbox{\boldmath{$ff$}}^{T}$, we get

$\mbox{\boldmath{$f$}}^{T}\mbox{\boldmath{$\mathcal{R}$}}^{t}\mbox{\boldmath{$R$}}^{s}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})^{m}\mbox{\boldmath{$S$}}^{l}\mbox{\boldmath{$f$}}=\sum_{i=q+1}^{n}\frac{\varrho_{i}^{t}\rho_{i}^{m+s}(\lambda
n\eta_{q,i})^{m}B_{q,i}^{2}}{\big{(}1+\lambda
n\eta_{q,i}\rho_{i}\big{)}^{m+l}}\\{1+o(1)\\},$

where $\mbox{\boldmath{$B$}}_{q}=(B_{q,1},\dots,B_{q,n})^{T}$ denotes
${\mbox{\boldmath{$\Phi$}}_{q}}^{T}\mbox{\boldmath{$f$}}$.

We proceed like in Lemma 1 of Serra and Krivobokova, (2017). To deal with the
trace, let

$g_{q,n}(y,\lambda)=\rho\left(\frac{\\{y/\lambda\\}^{1/\\{2q\\}}+(q+1)\pi/2}{\pi(n-1)}\right),\quad
h_{q,n}(y,\lambda)=\varrho\left(\frac{\\{y/\lambda\\}^{1/\\{2q\\}}+(q+1)\pi/2}{\pi(n-1)}\right),$

and define for $m\in\mathbb{N}\cup\\{0\\}$, $l\in\mathbb{N}$,
$t,s\in\mathbb{Z}$,

$\kappa_{q}(m,l,t,s,\varrho,\rho)=\frac{1}{2\pi
q}\lim_{n\to\infty}\int_{\lambda\pi^{2q}}^{\lambda\\{\pi(n-q)\\}^{2q}}\frac{y^{\frac{1}{2q}+m-1}g_{q,n}(y,\lambda)^{m+s}h_{q,n}(y,\lambda)^{t}}{\big{\\{}1+y\,g_{q,n}(y,\lambda)\big{\\}}^{m+l}}\,dy.$
(20)

Note that for any $\lambda>0$, within the integration range the argument of
$\rho$ and $\varrho$ above lies in $(0,1/2]$. Note also that if the spectral
densities satisfy $\delta\leq\varrho(t),\rho(t)\leq 1/\delta$, for some
$0<\delta<1$, $t\in[0,1]$, then

$\delta^{2m+l+t+s}\kappa_{q}(m,l)\leq\kappa_{q}(m,l,t,s,\varrho,\rho)\leq\delta^{-(2m+l+t+s)}\kappa_{q}(m,l),$

where $\kappa_{q}(m,l)$ are the constants from Lemma 1 of Serra and
Krivobokova, (2017), so that indeed the limit above is well defined.

With this notation we have

$\mbox{tr}\\{\mbox{\boldmath{$\mathcal{R}$}}^{t}\mbox{\boldmath{$R$}}^{s}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})^{m}\mbox{\boldmath{$S$}}^{l}\\}=\lambda^{-1/(2q)}\,\kappa_{q}(m,l,t,s,\varrho,\rho)\\{1+o(1)\\}.$
(21)

We deal now with the quadratic forms. From Serra and Krivobokova, (2017) it is
known that for $m\in\mathbb{N}\cup\\{0\\}$, $l\in\mathbb{N}$,

$\displaystyle\sum_{i=q+1}^{n}\frac{\lambda
n\eta_{q,i}B_{q,i}^{2}}{\big{(}1+\lambda n\eta_{q,i}\big{)}^{m}}$
$\displaystyle=n\lambda\|f^{(q)}\|^{2}\\{1+o(1)\\},$ $\displaystyle
f\in\mathcal{W}_{q},$ $\displaystyle\sum_{i=q+1}^{n}\frac{\lambda
n\eta_{q,i}B_{q,i}^{2}}{\big{(}1+\lambda n\eta_{q,i}\big{)}^{l}}$
$\displaystyle\leq n\lambda^{\beta/q}\|f^{(\beta)}\|^{2}\\{1+o(1)\\},$
$\displaystyle f\in\mathcal{W}_{\beta},$ $\displaystyle\quad q>\beta.$

In the case where $\delta\leq\rho_{i}\leq 1/\delta$, for some $0<\delta<1$,
exactly the same argument can be used to show that for
$m\in\mathbb{N}\cup\\{0\\}$, $l\in\mathbb{N}$,

$\displaystyle\sum_{i=q+1}^{n}\frac{\lambda
n\eta_{q,i}B_{q,i}^{2}}{\big{(}1+\lambda n\eta_{q,i}\rho_{i}\big{)}^{m}}$
$\displaystyle=n\lambda\|f^{(q)}\|^{2}\\{1+o(1)\\},$ $\displaystyle
f\in\mathcal{W}_{q},$ $\displaystyle\sum_{i=q+1}^{n}\frac{\lambda
n\eta_{q,i}B_{q,i}^{2}}{\big{(}1+\lambda n\eta_{q,i}\rho_{i}\big{)}^{l}}$
$\displaystyle\leq n\lambda^{\beta/q}\|f^{(\beta)}\|^{2}\\{1+o(1)\\},$
$\displaystyle f\in\mathcal{W}_{\beta},$ $\displaystyle\quad q>\beta.$

Combining the previous display with the representation for the quadratic form
from before, we conclude in particular that for $m\in\mathbb{N}\cup\\{0\\}$,

$\displaystyle\mbox{\boldmath{$f$}}^{T}\mbox{\boldmath{$R$}}^{-1}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$S$}}^{m}\mbox{\boldmath{$f$}}$
$\displaystyle=n\lambda\|f^{(q)}\|^{2}\\{1+o(1)\\},$ $\displaystyle
f\in\mathcal{W}_{q},$ (22)
$\displaystyle\mbox{\boldmath{$f$}}^{T}\mbox{\boldmath{$R$}}^{-1}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$S$}}^{m}\mbox{\boldmath{$f$}}$
$\displaystyle\leq n\lambda^{\beta/q}\|f^{(\beta)}\|^{2}\\{1+o(1)\\},$
$\displaystyle f\in\mathcal{W}_{\beta},$ $\displaystyle\quad q>\beta,$

which covers all the quadratic forms that we need to control in this paper.

###  B.4 $R$ and its spectral density

Let $R$ be a Toeplitz matrix in that
$(\mbox{\boldmath{$R$}})_{i,j}=R_{i,j}=r_{i-j}$ for some sequence
$\\{r_{i}\\}_{i\in\mathbb{Z}}$. Let us further assume that $R_{i,j}=0$, if
$|i-j|>m$, $m\in\mathbb{N}$. Denote the Fourier spectrum of $R$ as

$\rho(x)=\sum_{k=-m}^{m}r_{k}\exp(\mbox{\boldmath{$i$}}kx),\quad\text{so
that}\quad
r_{k}=\frac{1}{2\pi}\int_{-\pi}^{\pi}\rho(x)\exp(-\mbox{\boldmath{$i$}}kx)\,dx,$

where $i$ denotes the imaginary unit. If $R$ is the correlation matrix of a
stationary process, then additionally we have that $r_{i-j}=r_{j-i}$ so that
we can write

$\rho(x)=1+2\sum_{k=1}^{m}r_{k}\cos(kx)\quad\text{and}\quad
r_{k}=\frac{1}{2\pi}\int_{-\pi}^{\pi}\cos(kx)\rho(x)\,dx=\int_{0}^{1}\cos(k\pi
x)\rho(\pi x)\,dx.$

In practice we approximate $r_{k}$ via

$r_{k}=\frac{1}{n}\sum_{i=1}^{n}\cos\left(k\pi\frac{i-1}{n-1}\right)\rho_{i}+O(n^{-1}),\quad\mbox{where}\;\;\rho_{i}=\rho\left(\pi\frac{i-1}{n-1}\right).$

For more details see Gray, (1972).

The following lemmas are crucial for our approach.

###### Lemma 2

Let $R$ be a $n\times n$ covariance matrix of a stationary process with a
Lipschitz continuous spectral density $\rho$. If $R_{i,j}=0$, for $|i-j|>m$,
$m\in\mathbb{N}$, then

$\\{{\mbox{\boldmath{$\Phi$}}}_{q}^{T}\mbox{\boldmath{$R$}}{\mbox{\boldmath{$\Phi$}}}_{q}\\}_{i,j}={\rho}_{j}\delta_{i,j}+\mathbb{I}\\{|i-j|\mbox{
is even}\\}O(n^{-1}),\;\;i,j=q+1,\ldots,n.$

Proof
Noting that $x_{i}=(i-1)/(n-1)$ and denoting $x_{j,q}=\\{j-(q+1)/2)\\}/(n-1)$
we get for $i,j=q+1,\ldots,n$

$\displaystyle\begin{aligned}
\\{{\mbox{\boldmath{$\Phi$}}}_{q}^{T}\mbox{\boldmath{$R$}}{\mbox{\boldmath{$\Phi$}}}_{q}\\}_{i,j}&=\sum_{l=1}^{n}\phi_{q,i}(x_{l})\sum_{k=1}^{n}\sqrt{\frac{2}{n}}\cos\\{\pi(k-1)x_{j,q}+\pi(q-1)/4\\}r_{|l-k|}\\\
&+\sum_{l=1}^{n}\phi_{q,i}(x_{l})\sum_{k=1}^{n}\left[\phi_{q,j}(x_{k})-\sqrt{\frac{2}{n}}\cos\\{\pi(k-1)x_{j,q}+\pi(q-1)/4\\}\right]r_{|l-k|}\\\
&=\rho(\pi
x_{j,q})\delta_{i,j}-\\{1+(-1)^{|i-j|}\\}\sqrt{\frac{2}{n}}\sum_{l=1}^{m}\phi_{q,i}(x_{l})\sum_{k=1}^{m-l+1}\cos(\pi
kx_{j,q})r_{k+l-1}\\\ &-\rho(\pi
x_{j,q})\sum_{l=1}^{n}\phi_{q,i}(x_{l})\left[\phi_{q,j}(x_{l})-\sqrt{\frac{2}{n}}\cos\\{\pi(l-1)x_{j,q}+\pi(q-1)/4\\}\right]\\\
&+\sum_{l=1}^{n}\phi_{q,i}(x_{l})\sum_{k=1}^{n}\left[\phi_{q,j}(x_{k})-\sqrt{\frac{2}{n}}\cos\\{\pi(k-1)x_{j,q}+\pi(q-1)/4\\}\right]r_{|l-k|}.\end{aligned}$
(23)

We used that $R_{i,j}=r_{|i-j|}=0$ for $|i-j|>m$, for some fixed $m$
independent on $n$ and that cosine and sine are even and odd functions,
respectively. Three last terms in (23) vanish for $|i-j|$ odd; for the last
two terms this follows since $\phi_{q,j}(x_{k})$ has the same number of sign
changes as $\cos\\{\pi(k-1)x_{j,q}+\pi(q-1)/4\\}$ (see section B.1) and
$\phi_{q,i}(x_{k})$ is an even function for $i$ odd and odd function for $i$
even. If $|i-j|$ is even, then the second term is of order $O(n^{-1})$ since
$\phi_{q,i}(x)=O(n^{-1/2})$ and $m$ is fixed, small and independent on $n$. To
see that the last two terms are of order $O(n^{-1})$ for $|i-j|$ even, use the
Cauchy-Schwarz inequality, Lemma 1 and the result from Section B.1 that
$\|\phi_{q,i}-\psi_{q,i}/\sqrt{n}\|_{L_{2}}=O(n^{-3/2})$, where $\psi_{q,i}$
is the Demmler-Reinsch basis of the Sobolev space ${\cal{W}}_{\beta}(M)$ given
in (15), with the tail parts $T_{i}(x)$ vanishing exponentially away from the
boundaries. Finally, since $\rho$ is Lipschitz continuous, it follows
$\rho(\pi x_{j,q})=\rho(\pi x_{j})+O(n^{-1})=\rho_{j}+O(n^{-1})$.

###### Lemma 3

Under assumptions of Lemma 2, as $n$ goes to infinity, the perturbation
matrices $\Delta$ from (17) has a bounded number of eigenvalues of order
$O(1)$, and the remaining eigenvalues are $o(1)$.

Proof
Let
${\mbox{\boldmath{$\Phi$}}}_{q}^{T}\mbox{\boldmath{$R$}}{\mbox{\boldmath{$\Phi$}}}_{q}=\mbox{\boldmath{$U$}}\mbox{diag}(\mbox{\boldmath{$\zeta$}})\mbox{\boldmath{$U$}}^{t}$.
On the one hand, $\zeta_{i}=\rho_{i}+O(n^{-1})$, on the other hand, from Lemma
2 we know that
$\\{{\mbox{\boldmath{$\Phi$}}}_{q}^{T}\mbox{\boldmath{$R$}}{\mbox{\boldmath{$\Phi$}}}_{q}\\}_{i,j}={\rho}_{j}\delta_{i,j}+O(n^{-1})$,
$i=q+1,\ldots,n$, so that

$\\{{\mbox{\boldmath{$\Phi$}}}_{q}^{T}\mbox{\boldmath{$R$}}{\mbox{\boldmath{$\Phi$}}}_{q}\\}_{i,i}=\\{\mbox{\boldmath{$U$}}\mbox{diag}(\mbox{\boldmath{$\zeta$}})\mbox{\boldmath{$U$}}^{t}\\}_{i,i}=U_{i,i}^{2}\zeta_{i}+\sum_{j\neq
i}U_{i,j}^{2}\zeta_{j}=\rho_{i}+O(n^{-1})=\zeta_{i}+O(n^{-1}).$

Hence,

$(U_{i,i}^{2}-1)\zeta_{i}+\sum_{j\neq
i}U_{i,j}^{2}\zeta_{j}=O(n^{-1}),\;i=q+1,\ldots,n.$

Note that both $\zeta_{i}$ and $\rho_{i}$ are positive and $U$ is an
orthogonal matrix. Therefore, for each $j$, it holds that
$U_{i,i}^{2}=1+O(n^{-1})$ and there is a bounded number of
$U_{i,j}=O(n^{-1/2})$, while the remaining $n\\{1+o(1)\\}$ elements
$U_{i,j}=O(n^{-1})$. Hence, $U$ is asymptotically a diagonal matrix.
We need to find the bound on the eigenvalues of $\Delta$, which are the same
as the eigenvalues of

$\mbox{diag}\left(\frac{1}{\sqrt{\mbox{\boldmath{$\rho$}}}}\right){\mbox{\boldmath{$\Phi$}}}_{q}^{T}\mbox{\boldmath{$R$}}{\mbox{\boldmath{$\Phi$}}}_{q}\mbox{diag}\left(\frac{1}{\sqrt{\mbox{\boldmath{$\rho$}}}}\right)-\mbox{\boldmath{$I$}}_{n}.$

Now,

${\mbox{\boldmath{$\Phi$}}}_{q}^{T}\mbox{\boldmath{$R$}}{\mbox{\boldmath{$\Phi$}}}_{q}=\mbox{\boldmath{$U$}}\mbox{diag}(\mbox{\boldmath{$\rho$}})\mbox{\boldmath{$U$}}^{t}+\mbox{\boldmath{$U$}}\mbox{diag}(\mbox{\boldmath{$\zeta$}}-\mbox{\boldmath{$\rho$}})\mbox{\boldmath{$U$}}^{t}.$

Since the eigenvalues of the second matrix are $O(n^{-1})$, except for the
first $q$ eigenvalues, it is enough to consider

$\mbox{\boldmath{$V$}}=\mbox{diag}\left(\frac{1}{\sqrt{\mbox{\boldmath{$\rho$}}}}\right)\mbox{\boldmath{$U$}}\mbox{diag}(\mbox{\boldmath{$\rho$}})\mbox{\boldmath{$U$}}^{t}\mbox{diag}\left(\frac{1}{\sqrt{\mbox{\boldmath{$\rho$}}}}\right).$

The diagonal elements $V_{i,i}$, $i=q+1,\ldots,n$ are given by

$\sum_{j:|i-j|\leq
k}U_{i,j}^{2}\frac{\rho_{i}}{\rho_{j}}+\sum_{j:|i-j|>k}U_{i,j}^{2}\frac{\rho_{i}}{\rho_{j}}=\sum_{j:|i-j|\leq
k}U_{i,j}^{2}\\{1+O(n^{-1})\\}+O(n^{-1})=\sum_{j=1}^{n}U_{i,j}^{2}\\{1+O(n^{-1})\\},$

for some fixed bounded $k$. We used that $\rho_{i}/\rho_{j}=1+O(n^{-1})$ for
$|i-j|\leq k$ due to Lipschitz continuity of the spectral density $\rho$,
while for $|i-j|>k$ the ratio $\rho_{i}/\rho_{j}$ is bounded by the assumption
on $R$. In the same way, the off-diagonal elements $V_{i,j}$ are given by

$\sum_{l=1}^{n}U_{i,l}U_{j,l}\frac{\rho_{l}}{\sqrt{\rho_{i}\rho_{j}}}=\sum_{l=1}^{n}U_{i,l}U_{j,l}\\{1+O(n^{-1})\\}.$

Hence, we can write
$V_{i,j}=\\{\mbox{\boldmath{$U$}}\mbox{diag}\\{1+O(n^{-1})\\}\mbox{\boldmath{$U$}}^{t}\\}_{i,j}=\delta_{i,j}+O(n^{-1})$,
$i,j=q+1,\ldots,n$ and the $i$-th eigenvalue of $V$ is $1+O(n^{-1})$,
$i=q+1,\ldots,n$. This completes the proof.
A similar result holds for the matrices $\mbox{\boldmath{$\Delta$}}_{t,s}$
from (19) but we omit the details.

###  B.5 Derivation of the estimating equations

Exact estimating equations given in Section 2.2 follow from the expression for
the log-likelihood in (5) and the identities for the derivatives from Appendix
B.2. These exact expressions for the estimating equations become difficult to
evaluate in a numerically stable way if $n$ and/or $q$ are large. We suggest
to use the approximate expressions that can be evaluated in a numerically
stable way and provide excellent approximations for the exact expressions even
for small sample size.

The approximate estimating equations for
$T_{\lambda}(\lambda,q,\mbox{\boldmath{$\rho$}})$,
$T_{q}(\lambda,q,\mbox{\boldmath{$\rho$}})$, and
$T_{\rho_{i}}(\lambda,q,\mbox{\boldmath{$\rho$}})$ are obtained by using the
approximation (16) of the smoother $S$.

$\displaystyle T_{\lambda}(\lambda,q,\mbox{\boldmath{$\rho$}})$
$\displaystyle=\frac{1}{n_{\lambda,q}}\left\\{\sum_{i=q+1}^{n}\frac{B_{i}^{2}\lambda
n\eta_{q,i}}{\big{(}1+\lambda
n\eta_{q,i}\rho_{i}\big{)}^{2}}-\hat{\sigma}^{2}\sum_{i=q+1}^{n}\frac{1}{1+\lambda
n\eta_{q,i}\rho_{i}}\right\\}\\{1+o(1)\\},$ $\displaystyle
T_{q}(\lambda,q,\mbox{\boldmath{$\rho$}})$
$\displaystyle=\frac{1}{n_{\lambda,q}^{\prime}}\left\\{\sum_{i=q+1}^{n}\frac{B_{i}^{2}\lambda
n\eta_{q,i}\log(n\eta_{q,i})}{\big{(}1+\lambda
n\eta_{q,i}\rho_{i}\big{)}^{2}}-\hat{\sigma}^{2}\sum_{i=q+1}^{n}\frac{\log\big{(}n\eta_{q,i}\big{)}}{1+\lambda
n\eta_{q,i}\rho_{i}}\right\\}\\{1+o(1)\\},$ $\displaystyle
T_{\rho_{i}}(\lambda,q,\mbox{\boldmath{$\rho$}})$
$\displaystyle=\left(\frac{B_{i}^{2}\lambda n\eta_{q,i}\rho_{i}}{1+\lambda
n\eta_{q,i}\rho_{i}}-\rho_{i}\right)\\{1+o(1)\\},\quad i=1,\dots,n,$

where

$\hat{\sigma}^{2}=\frac{1}{n+1}\left(\sum_{i=q+1}^{n}\frac{B_{i}^{2}\lambda
n\eta_{q,i}}{1+\lambda n\eta_{q,i}\rho_{i}}+1\right)$

and the logarithm is taken entry-wise, and with the convention that
$0\cdot\log(0)=0$. The theoretical performance of the solutions to these
equations can be found in Appendix C.

##  Appendix C Proofs

Here we collect the proofs for Proposition 1, Proposition 2, Proposition 3,
Theorem 1, and Theorem 2 in Sections C.1 through C.5, respectively.

###  C.1 Uniform convergence of $T_{\lambda}$ to its expectation

Note that
$\mbox{\boldmath{$W$}}=\mbox{\boldmath{$\mathcal{R}$}}^{-1/2}\mbox{\boldmath{$Y$}}\sim
N(\mbox{\boldmath{$\mathcal{R}$}}^{-1/2}\mbox{\boldmath{$f$}},\sigma^{2}\mbox{\boldmath{$I$}}_{n})$
so that if $L$ is any symmetric, positive semi-definite matrix, then
$\mathbb{E}\big{(}\mbox{\boldmath{$W$}}^{T}\mbox{\boldmath{$L$}}\mbox{\boldmath{$W$}}\big{)}=\mbox{\boldmath{$f$}}^{T}\mbox{\boldmath{$\mathcal{R}$}}^{-1/2}\mbox{\boldmath{$L$}}\mbox{\boldmath{$\mathcal{R}$}}^{-1/2}\mbox{\boldmath{$f$}}+\sigma^{2}\mbox{tr}(\mbox{\boldmath{$L$}})$.
We write

$T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})=\frac{1}{n_{\lambda,q}}\mbox{\boldmath{$Z$}}^{T}\mbox{\boldmath{$L$}}\mbox{\boldmath{$Z$}}-\frac{1}{n_{\lambda,q}}\left\\{\mbox{\boldmath{$Z$}}^{T}\mbox{\boldmath{$L$}}^{\prime}\mbox{\boldmath{$Z$}}+1\right\\}\frac{\mbox{tr}(\mbox{\boldmath{$S$}})-q}{n+1},$

where
$\mbox{\boldmath{$Z$}}=\mbox{\boldmath{$\Phi$}}_{q}^{T}\mbox{\boldmath{$\mathcal{R}$}}^{-1/2}\mbox{\boldmath{$Y$}}\sim
N(\mbox{\boldmath{$z$}},\sigma^{2}\mbox{\boldmath{$I$}}_{n})$, for
$\mbox{\boldmath{$z$}}={\mbox{\boldmath{$\Phi$}}_{q}}^{T}\mbox{\boldmath{$\mathcal{R}$}}^{-1/2}\mbox{\boldmath{$f$}}$,
and the symmetric, positive definite matrices $L$ and
$\mbox{\boldmath{$L$}}^{\prime}$ are

$\mbox{\boldmath{$L$}}=\mbox{\boldmath{$\Phi$}}_{q}^{T}\mbox{\boldmath{$\mathcal{R}$}}^{1/2}\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$S$}}\mbox{\boldmath{$\mathcal{R}$}}^{1/2}\mbox{\boldmath{$\Phi$}}_{q},\quad\quad\mbox{\boldmath{$L$}}^{\prime}=\mbox{\boldmath{$\Phi$}}_{q}^{T}\mbox{\boldmath{$\mathcal{R}$}}^{1/2}\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$\mathcal{R}$}}^{1/2}\mbox{\boldmath{$\Phi$}}_{q}.$

Denote by $\Lambda_{n}$ a finite collection of values for $\lambda$ such that
$\lambda\to 0$, and $n\lambda\to\infty$. For any $R$, we can upper bound
$\sup_{\lambda\in\Lambda_{n}}|T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})-{\mathbb{E}}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})|$
as the sum of two suprema, one involving $L$, and one involving
$\mbox{\boldmath{$L$}}^{\prime}$.

Consider the term involving $L$ (the term involving
$\mbox{\boldmath{$L$}}^{\prime}$ is controlled in exactly the same way, and is
in fact smaller). By Lemma 3 all but a finite number of singular values of
$\mbox{\boldmath{$\Delta$}}_{1/2,-1}$ and of
$\mbox{\boldmath{$\Delta$}}_{1/2,0}$ are $o(1)$, the rest of which are $O(1)$
so then it suffices to control

$\sup_{\lambda\in\Lambda_{n}}\left|\frac{1}{n_{\lambda,q}}\mbox{\boldmath{$Z$}}^{T}\mbox{diag}\left[\left\\{\frac{\lambda\,n\eta_{q,i}\varrho_{i}}{\big{(}1+\lambda\,n\eta_{q,i}\rho_{i}\big{)}^{2}}\right\\}_{i}\right]\mbox{\boldmath{$Z$}}-\mathbb{E}\frac{1}{n_{\lambda,q}}\mbox{\boldmath{$Z$}}^{T}\mbox{diag}\left[\left\\{\frac{\lambda\,n\eta_{q,i}\varrho_{i}}{\big{(}1+\lambda\,n\eta_{q,i}\rho_{i}\big{)}^{2}}\right\\}_{i}\right]\mbox{\boldmath{$Z$}}\right|.$

To control this we use the following inequalities for quadratic form in
Gaussian random variables; cf. Enikeeva et al., (2018). The random variables
$Z_{i}^{2}/\sigma^{2}$, $i=q+1,\dots,n$, are independent, and distributed like
$\mathcal{X}_{1}(z_{i}^{2})$, a non-central $\mathcal{X}^{2}$-distribution
with non-centrality parameter $z_{i}^{2}$, and 1 degree of freedom. If we
abbreviate $M=\sum_{i=q+1}^{n}c_{\lambda,q,i}Z_{i}^{2}/\sigma^{2}$, for
$c_{\lambda,q,i}\geq 0$, and if
$c_{\lambda,q}=\max_{i=q+1,\dots,n}c_{\lambda,q,i}$, then for any $x>0$,

$\displaystyle\mathbb{P}\left\\{M-\mathbb{E}M\leq-\sqrt{2\mathbb{V}M\,x}\right\\}$
$\displaystyle\leq\exp(-x),$
$\displaystyle\mathbb{P}\left\\{M-\mathbb{E}M>\sqrt{2\mathbb{V}M\,x}+2c_{\lambda,q}\,x\right\\}$
$\displaystyle\leq\exp(-x),$

where $\mathbb{E}M=\sum_{i=q+1}^{n}c_{\lambda,q,i}(1+z_{i}^{2})$, and
$\mathbb{V}M=2\sum_{i=q+1}^{n}c_{\lambda,q,i}^{2}(1+2z_{i}^{2})$. Note that
using the bounds on the eigenvalues of autocorrelation matrices,

$2\mbox{tr}\left[\mbox{\boldmath{$\mathcal{R}$}}^{2}\mbox{\boldmath{$R$}}^{-2}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}}^{*})^{2}\\{\mbox{\boldmath{$S$}}^{*}\\}^{2}\right]\leq\mathbb{V}M\leq
2\mbox{tr}\left[\mbox{\boldmath{$\mathcal{R}$}}^{2}\mbox{\boldmath{$R$}}^{-2}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}}^{*})^{2}\\{\mbox{\boldmath{$S$}}^{*}\\}^{2}\right]+4\delta^{-5}\,n\lambda^{1\wedge(\beta/q)},$

such that if we assume that $\beta>1/2$,

$O\big{\\{}\lambda^{-1/(2q)}\big{\\}}\leq\mathbb{V}M\leq
O\big{\\{}\lambda^{-1/(2q)}+n\lambda^{1\wedge(\beta/q)}\big{\\}}\leq
O\big{\\{}\lambda^{-1/(2q)}+n\lambda^{1/(2q)}\big{\\}}.$

Note that this holds, in fact, for the variance of all of the quadratic forms
involved in the definition of the estimating equation for $\lambda$, hence the
scaling used there.

Combining the probability bounds with a union bound we get that for any
$\epsilon>0$,

$\mathbb{P}\bigg{(}\sup_{\lambda\in\Lambda_{n}}\frac{|M-\mathbb{E}M|}{n_{\lambda,q}}>\epsilon\bigg{)}\leq
2|\Lambda_{n}|\sup_{\lambda\in\Lambda_{n}}\exp(-x),$

as long as
$\epsilon\,n_{\lambda,q}\geq\sqrt{2\mathbb{V}M\,x}+2c_{\lambda,q}\,x$. The
largest $x$ such that this inequality holds is

$\frac{\mathbb{V}M}{2c_{\lambda,q}}\left(\sqrt{\frac{4\epsilon\,n_{\lambda,q}c_{\lambda,q}}{\mathbb{V}M}+1}-1\right)^{2}\geq
O\big{\\{}\lambda^{-1/(2q)}\big{\\}},$

where the lower bound (which comes from using the bounds on the variance of
$M$, noting that $c_{\lambda,q}=1+o(1)$, and by definition of $n_{\lambda,q}$)
holds for any $\epsilon>0$.

The conclusion is that as long as
$|\Lambda_{n}|\exp\\{-\inf_{\lambda\in\Lambda_{n}}\lambda^{-1/(2q)}\\}=o(1)$,
i.e., if the grid $\Lambda_{n}$ does not have too many elements, and if the
largest element in $\Lambda_{n}$ is not too large, then the estimating
equation for $\lambda$ converges uniformly over $\Lambda_{n}$ to its
expectation. Note that this holds true over all $R$ that satisfy our
eigenvalue bounds, so we conclude that

$\sup_{\mbox{\boldmath{$R$}}\in\mathcal{M}_{n,\delta}}\mathbb{P}\bigg{\\{}\sup_{\lambda\in\Lambda_{n}}\big{|}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})-\mathbb{E}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})\big{|}>\epsilon\bigg{\\}}\to
0,\quad\epsilon>0,\quad n\to\infty.$ (9)

###  C.2 Consistency of the preliminary estimate of $\lambda$

Convergence of the estimator to the respective oracle follows from Proposition
1 and Theorem 5.9 of van der Vaart, (1998). It remains to compute the oracle.

Firstly, the expression from the expectation of
$T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})$ for arbitrary $R$ is given in
Appendix C.3. The identity matrix $\mbox{\boldmath{$I$}}_{n}$ can be quite off
from the true correlation matrix, so the trace that features in the
expectation of $T_{\lambda}$ is not necessarily positive. By (21), this trace
is $\kappa_{q}({\mbox{\boldmath{$I$}}_{n}})\,\lambda^{-1/(2q)}\\{1+o(1)\\}$,
where

$\kappa_{q}({\mbox{\boldmath{$I$}}_{n}})=\kappa_{q}(0,2,1,0,\varrho,1)+\kappa_{q}(0,1,0,0,\varrho,1)-\kappa_{q}(0,1,1,0,\varrho,1).$
(24)

To ensure that a solution to the estimating equation
${\mathbb{E}}T_{\lambda}(\lambda,q,\mbox{\boldmath{$I$}}_{n})=0$ exists, we
assume that $\kappa_{q}({\mbox{\boldmath{$I$}}_{n}})$ is positive. Again using
(21) and (22), and solving for $\lambda$, we conclude that if
$\kappa_{q}({\mbox{\boldmath{$I$}}_{n}})>0$, then the solution to the
estimating equation when $f\in\mathcal{W}_{q}$,
$\mbox{\boldmath{$R$}}=\mbox{\boldmath{$I$}}_{n}$, is

$\lambda_{q,\mbox{\boldmath{$I$}}}=\bigg{[}\frac{n\|f^{(q)}\|^{2}}{\sigma^{2}\,\kappa_{q}({\mbox{\boldmath{$I$}}_{n}})}\\{1+o(1)\\}\bigg{]}^{-\frac{2q}{2q+1}},$
(10)

and if $f\in\mathcal{W}_{\beta}$, $q>\beta$, then

$\lambda_{q,\mbox{\boldmath{$I$}}}\geq\bigg{[}\frac{n\|f^{(\beta)}\|^{2}}{\sigma^{2}\,\kappa_{q}({\mbox{\boldmath{$I$}}_{n}})}\\{1+o(1)\\}\bigg{]}^{-\frac{2q}{2\beta+1}}.$
(11)

###  C.3 Consistency of the preliminary estimator of $\mathcal{R}$

The result follows by just combining known results. The estimators for the
eigenvalues $\sigma^{2}\mbox{\boldmath{$\varrho$}}$ of the covariance matrix
$\sigma^{2}\mbox{\boldmath{$\mathcal{R}$}}$ are

$Z_{i}=\hat{\sigma}^{2}\tilde{\rho}_{i}=\frac{\lambda
n\eta_{q,i}B_{q,i}^{2}}{1+\lambda n\eta_{q,i}},\qquad i=q+1,\dots,n,$

where $\mbox{\boldmath{$B$}}\sim
N\big{(}\mbox{\boldmath{$\Phi$}}_{q}\mbox{\boldmath{$f$}},\sigma^{2}\mbox{\boldmath{$\Phi$}}_{q}\mbox{\boldmath{$\mathcal{R}$}}\mbox{\boldmath{$\Phi$}}_{q}^{T}\big{)}$.
Based on Lemma 2, as $n$ grows, this is distributed like
$N\big{(}\mbox{\boldmath{$b$}}_{q},\sigma^{2}\mbox{diag}\\{\mbox{\boldmath{$\varrho$}}\\}\big{)}$,
where
$\mbox{\boldmath{$b$}}_{q}=\mathbb{E}\mbox{\boldmath{$B$}}_{q}=\mbox{\boldmath{$\Phi$}}_{q}\mbox{\boldmath{$f$}}$.
Note then that

$\mathbb{E}\big{(}Z_{i}-\sigma^{2}\varrho_{i}\big{)}=\frac{b_{i}^{2}\,\lambda
n\eta_{q,i}}{1+\lambda n\eta_{q,i}}-\frac{\sigma^{2}\varrho_{i}}{1+\lambda
n\eta_{q,i}},\quad\hbox{and}\quad\mathbb{V}Z_{i}=\left(\frac{\lambda
n\eta_{q,i}}{1+\lambda
n\eta_{q,i}}\right)^{2}\left(6\sigma^{2}b_{i}^{2}\varrho_{i}+3\sigma^{4}\varrho_{i}^{2}\right).$

The $Z_{i}$, as estimates of $\sigma^{2}\varrho_{i}$ are not consistent – they
are biased and their variances do not converge to zero. To obtain consistent
estimates we smooth the $Z_{i}$ using a smoothing spline. Let
$\bar{\mbox{\boldmath{$S$}}\,}=\mbox{\boldmath{$S$}}_{\xi,p,\mbox{\boldmath{$I$}}}$
be the spline smoother, and consider the estimates
$\hat{\mbox{\boldmath{$\rho$}}\,}=\bar{\mbox{\boldmath{$S$}}\,}\tilde{\mbox{\boldmath{$\rho$}}\,}$.

The heteroscedasticity of the noise terms is not an issue. As remarked by
Eggermont and LaRiccia, (2009, p.233), as long as the variances of the noise
terms are bounded uniformly over $n$, the usual smoothing spline estimator
(for homoskedastic noise) attains the optimal rate. The authors remark that
weighted spline estimators based on estimates of the variances of the noise
terms do not lead to large differences in the estimator.

Let us assume that the autocorrelations of the stationary noise process in (1)
satisfy $r_{i}=O(|i|^{-2\alpha})$, $\alpha>1/2$, such that in particular
$\sum_{i=0}^{\infty}|r_{i}|$ is finite. The bias can be shown to vanish by
combining this with (22) and using the consistency of the preliminary
estimator for $\lambda$ based on the identity as working correlation. Then the
risk of the estimates of the eigenvalues of the autocovariance matrix is

$\mathbb{E}\|\hat{\sigma}^{2}\hat{\mbox{\boldmath{$\rho$}}\,}-\sigma^{2}\mbox{\boldmath{$\varrho$}}\|_{\infty}=o(1).$

This follows from Eggermont and LaRiccia, (2006) where the authors determine
uniform rates for spline smoothers.

###  C.4 Consistency of estimates of $\lambda$

By Lemma 3 all but a finite number of singular values of
$\mbox{\boldmath{$\Delta$}}_{1/2,-1}(\mbox{\boldmath{$\mathcal{R}$}},\hat{\mbox{\boldmath{$R$}}\,})$,
$\mbox{\boldmath{$\Delta$}}_{1/2,}(\mbox{\boldmath{$\mathcal{R}$}},\hat{\mbox{\boldmath{$R$}}\,})$,
and $\mbox{\boldmath{$\Delta$}}(\hat{\mbox{\boldmath{$R$}}\,})$ are
$o_{P}(1)$, the rest of which are $O_{P}(1)$. By construction,
$\hat{\mbox{\boldmath{$R$}}\,}\in{\mathcal{M}}_{n,\delta}$, with probability
1, for some $0<\delta<1$. The same can be assumed without loss of generality
for $\mathcal{R}$. To control
$T_{\lambda}(\lambda,q,\hat{\mbox{\boldmath{$R$}}\,})-\mathbb{E}T_{\lambda}(\lambda,q,\mbox{\boldmath{$\mathcal{R}$}})$
uniformly over $\lambda\in\Lambda_{n}$ we write this difference as

$\big{\\{}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})-\mathbb{E}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})\\}|_{\mbox{\boldmath{$R$}}=\hat{\mbox{\boldmath{$R$}}\,}}+\big{\\{}\mathbb{E}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})\\}|_{\mbox{\boldmath{$R$}}=\hat{\mbox{\boldmath{$R$}}\,}}-\mathbb{E}T_{\lambda}(\lambda,q,\mbox{\boldmath{$\mathcal{R}$}}),$
(25)

where in the first three terms $\hat{\mbox{\boldmath{$R$}}\,}$ is plugged in
for $R$. In Appendix C.2 we show that

$\sup_{\mbox{\boldmath{$R$}}\in\mathcal{M}_{n,\delta}}\mathbb{P}\Big{\\{}\sup_{\lambda\in\Lambda_{n}}\big{|}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})-\mathbb{E}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})\big{|}>\epsilon\Big{\\}}\to
0,\quad n\to\infty.$ (9)

Using this fact we can control the first difference in (25) by noting that

$\displaystyle\mathbb{P}\left[\sup_{\lambda\in\Lambda_{n}}\left|\big{\\{}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})-\mathbb{E}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})\\}|_{\mbox{\boldmath{$R$}}=\hat{\mbox{\boldmath{$R$}}\,}}\right|>\epsilon\right]\leq$
$\displaystyle\qquad\leq\mathbb{P}\left[\sup_{\lambda\in\Lambda_{n}}\left|T_{\lambda}(\lambda,q,\hat{\mbox{\boldmath{$R$}}\,})-\mathbb{E}\big{\\{}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})\\}|_{\mbox{\boldmath{$R$}}=\hat{\mbox{\boldmath{$R$}}\,}}\right|>\epsilon\mid\hat{\mbox{\boldmath{$R$}}\,}\in{\mathcal{M}}_{n,\delta}\right]+\mathbb{P}\big{(}\hat{\mbox{\boldmath{$R$}}\,}\not\in{\mathcal{M}}_{n,\delta}\big{)}.$

Both terms on the right-hand-side converge to zero, as $n\to\infty$, so it
remains to bound the second difference in (25).

Let us denote
$\mbox{\boldmath{$\mathcal{S}$}}=\mbox{\boldmath{$\mathcal{S}$}}_{\lambda,q}=\mbox{\boldmath{$S$}}_{\lambda,q,\mbox{\boldmath{$\mathcal{R}$}}}$.
These are the ideal spline smoothers that we would have used if the
correlation matrix $\mathcal{R}$ were known. First we look at the expectation
of the estimates for the variance that feature in the expectation of the
estimating equation for $\lambda$. Denote

$\mbox{\boldmath{$\nabla$}}=\mbox{\boldmath{$\nabla$}}(\mbox{\boldmath{$\mathcal{R}$}},\mbox{\boldmath{$R$}})=\mbox{\boldmath{$\mathcal{R}-R$}},$
(26)

and let ${\mathcal{N}}_{n}={\mathcal{N}}_{n}(\mbox{\boldmath{$\mathcal{R}$}})$
be the collection of all matrices $R$ such that the largest eigenvalue of
$\nabla$ is $o(1)$. Then, for any
$\mbox{\boldmath{$R$}}\in{\mathcal{N}}_{n}\cap{\mathcal{M}}_{n,\delta}$, and
any $\lambda\in\Lambda_{n}$,

$\displaystyle{\mathbb{E}}\big{(}\hat{\sigma}^{2}_{\mbox{\boldmath{$R$}}}\big{)}$
$\displaystyle=\frac{1}{n+1}\Big{[}\mbox{\boldmath{$f$}}^{T}\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$f$}}+1+\sigma^{2}\mbox{tr}\big{\\{}(\mbox{\boldmath{$\nabla$}}\mbox{\boldmath{$R$}}^{-1}+{\mbox{\boldmath{$I$}}_{n}})(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})\big{\\}}\Big{]}$
$\displaystyle=\sigma^{2}\frac{n}{n+1}\Big{[}1+O\big{(}\lambda^{1\wedge(\beta/q)}\big{)}+O\big{\\{}\lambda^{-1/(2q)}/n\big{\\}}+O\\{\mbox{tr}(\mbox{\boldmath{$\nabla$}})/n\\}\Big{]}=\sigma^{2}\\{1+o(1)\\},$

so that
${\mathbb{E}}\big{(}\hat{\sigma}^{2}_{\mbox{\boldmath{$R$}}}\big{)}=\sigma^{2}\\{1+o(1)\\}={\mathbb{E}}\big{(}\hat{\sigma}^{2}_{\mbox{\boldmath{$\mathcal{R}$}}}\big{)}$,
if we assume w.l.g., that
$\mbox{\boldmath{$\mathcal{R}$}}\in\mathcal{M}_{n,\delta}$.

The estimating equation for $\lambda$ can be written in matrix notation as

$T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})=\frac{1}{n_{\lambda,q}}\Big{[}\mbox{\boldmath{$Y$}}^{T}\mbox{\boldmath{$R$}}^{-1}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$S$}}\mbox{\boldmath{$Y$}}-\hat{\sigma}^{2}\big{\\{}\mbox{tr}(\mbox{\boldmath{$S$}})-q\big{\\}}\Big{]},$

so simple computations and use of (21) and (22) show that for each $\lambda$,
$q$, $R$,

${\mathbb{E}}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})=\frac{1}{n_{\lambda,q}}\Big{[}\mbox{\boldmath{$f$}}^{T}\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$S$}}\mbox{\boldmath{$f$}}+\sigma^{2}\mbox{tr}\big{\\{}\mbox{\boldmath{$\mathcal{R}R$}}^{-1}(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$S$}}\big{\\}}-{\mathbb{E}}\hat{\sigma}^{2}\big{\\{}\mbox{tr}(\mbox{\boldmath{$S$}})-q\big{\\}}\Big{]},$

For any $R$, up to lower order terms,

$\displaystyle{\mathbb{E}}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})$
$\displaystyle=\frac{1}{n_{\lambda,q}}\Big{[}\mbox{\boldmath{$f$}}^{T}\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$S$}}\mbox{\boldmath{$f$}}+\sigma^{2}\mbox{tr}\big{\\{}(\mbox{\boldmath{$\nabla$}}\mbox{\boldmath{$R$}}^{-1}+{\mbox{\boldmath{$I$}}_{n}})(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$S$}}\big{\\}}-\sigma^{2}\mbox{tr}(\mbox{\boldmath{$S$}})\Big{]}.$

The difference
${\mathbb{E}}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})-{\mathbb{E}}T_{\lambda}(\lambda,q,\mbox{\boldmath{$\mathcal{R}$}})$
involves the difference of two quadratic forms and differences of traces,
which we treat separately. From (22), for any
$\mbox{\boldmath{$R$}}\in\mathcal{M}_{n,\delta}$, the influence of $R$ on the
particular quadratic form above lies only in the lower order terms so we can
directly conclude that

$\mbox{\boldmath{$f$}}^{T}\mbox{\boldmath{$R$}}^{-1}(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$S$}}\mbox{\boldmath{$f$}}-\mbox{\boldmath{$f$}}^{T}\mbox{\boldmath{$\mathcal{R}$}}^{-1}(\mbox{\boldmath{$I$}}_{n}-\mbox{\boldmath{$\mathcal{S}$}})\mbox{\boldmath{$\mathcal{S}$}}\mbox{\boldmath{$f$}}=o\left(n\lambda^{1\wedge(\beta/q)}\right).$

Now we control the differences of traces. If we assume that $\Delta$ is
$o(1)$, and that
$\mbox{\boldmath{$R$}}\in\mathcal{N}_{n}\cap\mathcal{M}_{n,\delta}$, then this
difference is

$\displaystyle\mbox{tr}\big{\\{}\mbox{\boldmath{$\nabla$}}\mbox{\boldmath{$R$}}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$S$}}\big{\\}}+\mbox{tr}\big{(}\mbox{\boldmath{$S$}}^{2}\big{)}-\mbox{tr}\big{(}\mbox{\boldmath{$\mathcal{S}$}}^{2}\big{)}=$
$\displaystyle\quad=o\left[\mbox{tr}\big{\\{}\mbox{\boldmath{$R$}}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$S$}})\mbox{\boldmath{$S$}}\big{\\}}\right]+\mbox{tr}\big{(}\\{\mbox{\boldmath{$S$}}^{*}\\}^{2}\big{)}\\{1+o(1)\\}-\mbox{tr}\big{(}\\{\mbox{\boldmath{$\mathcal{S}$}}^{*}\\}^{2}\big{)}\\{1+o(1)\\}$
$\displaystyle\quad=\mbox{tr}\big{\\{}(\mbox{\boldmath{$S$}}^{*}+\mbox{\boldmath{$\mathcal{S}$}}^{*})(\mbox{\boldmath{$S$}}^{*}-\mbox{\boldmath{$\mathcal{S}$}}^{*})\big{\\}}+o\\{\lambda^{-1/(2q)}\\},$

where $\mbox{\boldmath{$\mathcal{S}$}}^{*}$ is $\mbox{\boldmath{$S$}}^{*}$
where $R$ was replaced with $\mathcal{R}$, and since the eigenspaces of
$\mbox{\boldmath{$S$}}^{*}$ and $\mbox{\boldmath{$\mathcal{S}$}}^{*}$ are the
same so that these two matrices comute. The eigenvalues of
$\mbox{\boldmath{$S$}}^{*}-\mbox{\boldmath{$\mathcal{S}$}}^{*}$ satisfy

$\displaystyle\left|\frac{1}{1+\lambda n\eta_{q,i}\rho_{i}}-\frac{1}{1+\lambda
n\eta_{q,i}\varrho_{i}}\right|$ $\displaystyle\leq\frac{\lambda
n\eta_{q,i}|\rho_{i}-\varrho_{i}|}{\big{(}1+\lambda
n\eta_{q,i}\rho_{i}\big{)}\big{(}1+\lambda n\eta_{q,i}\varrho_{i}\big{)}}$
$\displaystyle\leq o\left\\{\frac{1}{\delta}\min\Big{(}\frac{1}{1+\lambda
n\eta_{q,i}\rho_{i}},\frac{1}{1+\lambda
n\eta_{q,i}\varrho_{i}}\Big{)}\right\\},$

since the largest eigenvalue of $\nabla$ is $o(1)$.

We conclude that

$\mbox{tr}\big{\\{}(\mbox{\boldmath{$S$}}^{*}+\mbox{\boldmath{$\mathcal{S}$}}^{*})(\mbox{\boldmath{$S$}}^{*}-\mbox{\boldmath{$\mathcal{S}$}}^{*})\big{\\}}=o\left[\mbox{tr}\big{(}\\{\mbox{\boldmath{$\mathcal{S}$}}^{*}\\}^{2}\big{)}\right]+o\left[\mbox{tr}\big{(}\\{\mbox{\boldmath{$S$}}^{*}\\}^{2}\big{)}\right]=o\\{\lambda^{-1/(2q)}\\},$

whence uniformly over $\lambda\in\Lambda_{n}$, and for any
$\mbox{\boldmath{$R$}}\in\mathcal{A}$, which denotes the collection of
matrices such that the largest singular value of $\Delta$ is $o(1)$, and so
that
${\mbox{\boldmath{$R$}}}\in{\mathcal{N}}_{n}\cap{\mathcal{M}}_{n,\delta}$, it
holds that

${\mathbb{E}}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})-{\mathbb{E}}T_{\lambda}(\lambda,q,\mbox{\boldmath{$\mathcal{R}$}})=\frac{o\left(n\lambda^{1\wedge(\beta/q)}\right)+o(\lambda^{-1/(2q)})}{n_{\lambda,q}}=\frac{o(n_{\lambda,q})}{n_{\lambda,q}}=o(1),$

because $\beta>1/2$. Since the probability that
$\hat{\mbox{\boldmath{$R$}}\,}\in\mathcal{A}$ converges to $1$, we conclude
that

$\displaystyle\mathbb{P}\left\\{\sup_{\lambda\in\Lambda_{n}}\big{|}\mathbb{E}\big{\\{}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})\\}|_{\mbox{\boldmath{$R$}}=\hat{\mbox{\boldmath{$R$}}\,}}-{\mathbb{E}}T_{\lambda}(\lambda,q,\mbox{\boldmath{$\mathcal{R}$}})\big{|}>\epsilon\right\\}\leq$
$\displaystyle\quad\leq\mathbb{P}\left\\{\sup_{\lambda\in\Lambda_{n}}\big{|}\mathbb{E}\big{\\{}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})\\}|_{\mbox{\boldmath{$R$}}=\hat{\mbox{\boldmath{$R$}}\,}}-{\mathbb{E}}T_{\lambda}(\lambda,q,\mbox{\boldmath{$\mathcal{R}$}})\big{|}>\epsilon\mid\hat{\mbox{\boldmath{$R$}}\,}\in\mathcal{A}\right\\}+\mathbb{P}(\hat{\mbox{\boldmath{$R$}}\,}\not\in\mathcal{A})=o(1).$

From this we conclude that if
$|\Lambda_{n}|\exp\\{-\inf_{\lambda\in\Lambda_{n}}\lambda^{-1/(2q)}\\}=o(1)$,
then

$\mathbb{P}\left\\{\sup_{\lambda\in\Lambda_{n}}\big{|}\big{\\{}T_{\lambda}(\lambda,q,\hat{\mbox{\boldmath{$R$}}\,})\\}-{\mathbb{E}}T_{\lambda}(\lambda,q,\mbox{\boldmath{$\mathcal{R}$}})\big{|}>\epsilon\right\\}\to
0,\quad n\to\infty.$

It remains to compute the oracle for $\lambda$. We have

${\mathbb{E}}T_{\lambda}(\lambda,q,\mbox{\boldmath{$\mathcal{R}$}})=\frac{1}{n_{\lambda,q}}\Big{\\{}\mbox{\boldmath{$f$}}^{T}\mbox{\boldmath{$\mathcal{R}$}}^{-1}({\mbox{\boldmath{$I$}}_{n}}-\mbox{\boldmath{$\mathcal{S}$}})\mbox{\boldmath{$\mathcal{S}$}}\mbox{\boldmath{$f$}}-\sigma^{2}\mbox{tr}\big{(}\mbox{\boldmath{$\mathcal{S}$}}^{2}\big{)}\Big{\\}}.$

Contrary to the case of the preliminary estimator where
$\mbox{\boldmath{$R$}}={\mbox{\boldmath{$I$}}_{n}}$, there always exists a
solution to
${\mathbb{E}}T_{\lambda}(\lambda,q,\mbox{\boldmath{$\mathcal{R}$}})=0$. Using
(21) and (22), we obtain the oracle
$\lambda_{q,\mbox{\boldmath{$\mathcal{R}$}}}$ which satisfies the following:
when $f\in\mathcal{W}_{q}$, and $q>1/2$, then

$\lambda_{q,\mbox{\boldmath{$\mathcal{R}$}}}=\bigg{[}\frac{n\|f^{(q)}\|^{2}}{\sigma^{2}\kappa_{q}(\mbox{\boldmath{$\mathcal{R}$}})}\\{1+o(1)\\}\bigg{]}^{-\frac{2q}{2q+1}},$
(12)

and if $f\in\mathcal{W}_{\beta}$, $q>\beta$, then

$\lambda_{q,\mbox{\boldmath{$\mathcal{R}$}}}\geq\bigg{[}\frac{n\|f^{(\beta)}\|^{2}}{\sigma^{2}\kappa_{q}(\mbox{\boldmath{$\mathcal{R}$}})}\\{1+o(1)\\}\bigg{]}^{-\frac{2q}{2\beta+1}},$
(13)

where

$\kappa_{q}(\mbox{\boldmath{$\mathcal{R}$}})=\kappa_{q}(0,2,0,0,\varrho,\varrho).$
(27)

By Theorem 5.9 of van der Vaart, (1998), we conclude that since
$T_{\lambda}(\lambda,q,\hat{\mbox{\boldmath{$R$}}\,})$ converges uniformly
over $\lambda$ to
$\mathbb{E}T_{\lambda}(\lambda,q,\mbox{\boldmath{$\mathcal{R}$}})$, then the
solution to $T_{\lambda}(\lambda,q,\hat{\mbox{\boldmath{$R$}}\,})=0$ converges
to the solution of
$\mathbb{E}T_{\lambda}(\lambda,q,\mbox{\boldmath{$\mathcal{R}$}})=0$. This
concludes the proof.

###  C.5 Consistency of the estimate of $q$

Let $\hat{\mbox{\boldmath{$R$}}\,}$ be any operator norm consistent estimator
for $\mathcal{R}$. To avoid excessive use of subscripts in what follows,
denote by $\hat{\lambda}_{q}$ the solution to
$T_{\lambda}(\lambda,q,\hat{\mbox{\boldmath{$R$}}\,})=0$, and by $\lambda_{q}$
the oracle that solves
$\mathbb{E}T_{\lambda}(\lambda,q,\mbox{\boldmath{$\mathcal{R}$}})=0$. Denote
$n_{\lambda,q}^{\prime}=n_{\lambda,q}(\log\lambda)^{2}$.

Using matrix notation and the identitites from before, we can write
$T_{q}(\lambda,q,\mbox{\boldmath{$R$}})$ as

$\frac{1}{n_{\lambda,q}^{\prime}}\bigg{[}\mbox{\boldmath{$Y$}}^{T}\mbox{\boldmath{$S$}}^{T}\mbox{\boldmath{$\Phi$}}_{q}\mbox{diag}\big{\\{}\lambda
n\mbox{\boldmath{$\eta$}}_{q}\log(n\mbox{\boldmath{$\eta$}}_{q})\big{\\}}\mbox{\boldmath{$\Phi$}}_{q}^{T}\mbox{\boldmath{$S$}}\mbox{\boldmath{$Y$}}-\hat{\sigma}^{2}\mbox{tr}\Big{[}\mbox{\boldmath{$S$}}\mbox{\boldmath{$\Phi$}}_{q}\mbox{diag}\big{\\{}\log(n\mbox{\boldmath{$\eta$}}_{q})\big{\\}}\mbox{\boldmath{$\Phi$}}_{q}^{T}\Big{]}\bigg{]}.$

(The justification for the scaling follows below.) If we use the fact that we
can write $\log(n\mbox{\boldmath{$\eta$}}_{q})=\log(1/\lambda)+\log(\lambda
n\mbox{\boldmath{$\eta$}}_{q})$, take expectations and note that since the
eigenvalues of $\mathcal{R}$ are bounded away from zero and from infinity, and
since the largest singular value of
$\mbox{\boldmath{$\Delta$}}(\mbox{\boldmath{$R$}})$ is $o(1)$, then both
$\mbox{tr}\big{[}\mbox{\boldmath{$\mathcal{R}$}}\mbox{\boldmath{$S$}}\mbox{\boldmath{$\Phi$}}_{q}\mbox{diag}\big{\\{}\lambda
n\mbox{\boldmath{$\eta$}}_{q}\log(\lambda
n\mbox{\boldmath{$\eta$}}_{q})\big{\\}}\mbox{\boldmath{$\Phi$}}_{q}^{T}\mbox{\boldmath{$S$}}\big{]}$
and also
$\mbox{tr}\big{[}\mbox{\boldmath{$S$}}\mbox{\boldmath{$\Phi$}}_{q}\mbox{diag}\big{\\{}\lambda
n\mbox{\boldmath{$\eta$}}_{q}\log(\lambda
n\mbox{\boldmath{$\eta$}}_{q})\big{\\}}\mbox{\boldmath{$\Phi$}}_{q}^{T}\big{]}$
are $o\\{\log(1/\lambda)\lambda^{-1/(2q)}\\}$ (cf. Lemma 2 of Serra and
Krivobokova, (2017)); consequently we have

$\mathbb{E}T_{q}(\lambda,q,\mbox{\boldmath{$R$}})=\frac{1}{n_{\lambda,q}^{\prime}}\mbox{\boldmath{$f$}}^{T}\mbox{\boldmath{$S$}}\mbox{\boldmath{$\Phi$}}_{q}\mbox{diag}\big{\\{}\lambda
n\mbox{\boldmath{$\eta$}}_{q}\log(\lambda
n\mbox{\boldmath{$\eta$}}_{q})\big{\\}}\mbox{\boldmath{$\Phi$}}_{q}^{T}\mbox{\boldmath{$S$}}\mbox{\boldmath{$f$}}+\frac{1}{\log(1/\lambda)}\mathbb{E}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}}).$

This last expression is useful because it tells us that if we plug
$\mathcal{R}$ into $R$, and $\lambda_{q}$ into $\lambda$, then the second term
on the right-hand-side of the previous display is zero, and so, since the
largest singular values of
$\mbox{\boldmath{$\Delta$}}(\mbox{\boldmath{$\mathcal{R}$}})$ is $o(1)$,

$\displaystyle
n_{\lambda,q}^{\prime}\cdot\mathbb{E}T_{q}(\lambda_{q},q,\mbox{\boldmath{$\mathcal{R}$}})=\mbox{\boldmath{$f$}}^{T}\mbox{\boldmath{$\mathcal{S}$}}\mbox{\boldmath{$\Phi$}}_{q}\mbox{diag}\big{\\{}\lambda_{q}n\mbox{\boldmath{$\eta$}}_{q}\log(\lambda_{q}n\mbox{\boldmath{$\eta$}}_{q})\big{\\}}\mbox{\boldmath{$\Phi$}}_{q}^{T}\mbox{\boldmath{$\mathcal{S}$}}\mbox{\boldmath{$f$}}=$
$\displaystyle\qquad=\sum_{i=q+1}^{n}\frac{b_{q,i}^{2}\lambda_{q}n\eta_{q,i}\log(\lambda_{q}n\eta_{q,i})}{(1+\lambda_{q}n\eta_{q,i}\varrho_{i})^{2}}\\{1+o(1)\\}=\sum_{i=q+1}^{n}\frac{b_{q,i}^{2}\lambda_{q}n\eta_{q,i}\log(\lambda_{q}n\eta_{q,i})}{(1+\lambda_{q}n\eta_{q,i})^{2}}\\{1+o(1)\\},$

since the eigenvalues of $\mathcal{R}$ are bounded away from zero and from
infinity; the last equality follows from the proof of Lemma 3 of Serra and
Krivobokova, (2017). (It is only the order of $\lambda n\eta_{q,i}$ in the
denominator that is relevant to compute the sum.) This lemma specifies, in
fact, the behaviour of the last sum depending on the smoothness of $f$, which
in turn characterises the oracle for $q$; cf. Theorem 2 of Serra and
Krivobokova, (2017).

It remains to show consistency. This follows the same steps as the the uniform
consistency proof for $T_{\lambda}$ (with slightly different weights in the
quadratic forms) so we only highlight the differences. As before, to show that
$\hat{q}$ is consistent for the respective oracle we control, uniformly over
$q\in Q_{n}$, the difference
$T_{q}\big{(}\hat{\lambda}_{q},q,\hat{\mbox{\boldmath{$R$}}\,}\big{)}-\mathbb{E}T_{q}(\lambda_{q},q,\mbox{\boldmath{$\mathcal{R}$}})$
which we write as

$\\{T_{q}\big{(}\lambda,q,\mbox{\boldmath{$R$}}\big{)}-\mathbb{E}T_{q}(\lambda,q,\mbox{\boldmath{$R$}})\\}|_{(\lambda,\mbox{\boldmath{$R$}})=(\hat{\lambda}_{q},\hat{\mbox{\boldmath{$R$}}\,})}+\\{\mathbb{E}T_{q}(\lambda,q,\mbox{\boldmath{$R$}})\\}|_{(\lambda,\mbox{\boldmath{$R$}})=(\hat{\lambda}_{q},\hat{\mbox{\boldmath{$R$}}\,})}-\mathbb{E}T_{q}(\lambda_{q},q,\mbox{\boldmath{$\mathcal{R}$}}).$

The first difference is controlled in the same way as the respective
difference in $T_{\lambda}$. In turn, this term can be written as the
difference between two quadratic forms and their respective expectations. If
$\mbox{\boldmath{$\Delta$}}_{1/2,-1}$ and $\mbox{\boldmath{$\Delta$}}_{1/2,0}$
are $o(1)$, then the first of these is, up to lower order terms, bounded by

$\frac{1}{n_{\lambda,q}^{\prime}}\left|\mbox{\boldmath{$Z$}}^{T}\mbox{diag}\left[\left\\{\frac{\lambda\,n\eta_{q,i}\varrho_{i}\log(n\eta_{q,i})}{\big{(}1+\lambda\,n\eta_{q,i}\rho_{i}\big{)}^{2}}\right\\}_{i}\right]\mbox{\boldmath{$Z$}}-\mathbb{E}\mbox{\boldmath{$Z$}}^{T}\mbox{diag}\left[\left\\{\frac{\lambda\,n\eta_{q,i}\varrho_{i}\log(n\eta_{q,i})}{\big{(}1+\lambda\,n\eta_{q,i}\rho_{i}\big{)}^{2}}\right\\}_{i}\right]\mbox{\boldmath{$Z$}}\right|,$

where
$\mbox{\boldmath{$Z$}}=\mbox{\boldmath{$\Phi$}}_{q}^{T}\mbox{\boldmath{$\mathcal{R}$}}^{-1/2}\mbox{\boldmath{$Y$}}\sim
N(\mbox{\boldmath{$z$}},\sigma^{2}\mbox{\boldmath{$I$}}_{n})$. We want to
control the supremum over $q\in Q_{n}$ of the previous display. (In analogy to
the $T_{\lambda}$ case, the remaining difference is of the same type and is
controlled in the same way.) If we denote by $M$ the quadratic form in the
previous display, then if we assume that $\beta>1/2$,

$O\big{\\{}\lambda^{-1/(2q)}\big{\\}}\leq\frac{\mathbb{V}M}{(\log\lambda)^{2}}\leq
O\big{\\{}\lambda^{-1/(2q)}+n\lambda^{1\wedge(\beta/q)}\big{\\}}\leq
O\big{\\{}\lambda^{-1/(2q)}+n\lambda^{1/(2q)}\big{\\}},$

where we use Lemmas 2 and 3 of Serra and Krivobokova, (2017). (This justifies
the scaling that we use in the criterium for $q$.) Furthermore, uniformly over
$\lambda$ such that $\lambda_{q}/\lambda\in[1-\zeta,1+\zeta]=L_{\zeta}$, and
$R$ such that
$\mbox{\boldmath{$R$}}\in\mathcal{N}_{n}\cap\mathcal{M}_{n,\delta}$, this
variance satisfies

$O\big{\\{}(\log\lambda_{q})^{2}\lambda_{q}^{-1/(2q)}\big{\\}}\leq\mathbb{V}M\leq
O\big{\\{}(\log\lambda_{q})^{2}\lambda_{q}^{-1/(2q)}+(\log\lambda_{q})^{2}n\lambda_{q}^{1/(2q)}\big{\\}}.$

Proceeding as before we have that uniformly over $\lambda$ such that
$\lambda_{q}/\lambda\in L_{\zeta}$ and
$\mbox{\boldmath{$R$}}\in\mathcal{N}_{n}\cap\mathcal{M}_{n,\delta}$,

$\mathbb{P}\bigg{(}\sup_{q\in
Q_{n}}\frac{|M-\mathbb{E}M|}{n_{\lambda,q}^{\prime}}>\epsilon\bigg{)}\leq
2|Q_{n}|\exp\left[-\inf_{q\in
Q_{n}}(\log\lambda_{q})^{2}\lambda_{q}^{-1/(2q)}\right].$

Let $\mathcal{A}$ denote the event that $\lambda_{q}/\hat{\lambda}_{q}\in
L_{\zeta}$, $0<\zeta<1$, and that
$\hat{\mbox{\boldmath{$R$}}\,}\in\mathcal{N}_{n}\cap\mathcal{M}_{n,\delta}$.
From the consistency of $\hat{\lambda}_{q}$ and
$\hat{\mbox{\boldmath{$R$}}\,}$ we known that the probability of this event
converges to 1, as $n\to\infty.$ The remaining term that comes from the
variance is actually smaller, but can be controlled in the same way. From this
we conclude that

$\displaystyle\mathbb{P}\left[\sup_{q\in
Q_{n}}\left|\big{\\{}T_{q}(\lambda,q,\mbox{\boldmath{$R$}})-\mathbb{E}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})\\}|_{(\lambda,\mbox{\boldmath{$R$}})=(\hat{\lambda}_{q},\hat{\mbox{\boldmath{$R$}}\,})}\right|\right]\leq$
$\displaystyle\qquad\leq\mathbb{P}\left[\sup_{q\in
Q_{n}}\left|\big{\\{}T_{q}(\lambda,q,\mbox{\boldmath{$R$}})-\mathbb{E}T_{\lambda}(\lambda,q,\mbox{\boldmath{$R$}})\\}|_{(\lambda,\mbox{\boldmath{$R$}})=(\hat{\lambda}_{q},\hat{\mbox{\boldmath{$R$}}\,})}\right|\mid\mathcal{A}\right]+\mathbb{P}\big{(}\mathcal{A}^{c}\big{)}=o(1).$

It remains to control
$\\{\mathbb{E}T_{q}(\lambda,q,\mbox{\boldmath{$R$}})\\}|_{(\lambda,\mbox{\boldmath{$R$}})=(\hat{\lambda}_{q},\hat{\mbox{\boldmath{$R$}}\,})}-\mathbb{E}T_{q}(\lambda_{q},q,\mbox{\boldmath{$\mathcal{R}$}})$.
As we saw before, under the event that the eigenvalues of $R$ and of
$\mathcal{R}$ are bounded away from zero and infinity (which occurs with
probability going to 1),

$\displaystyle\mathbb{E}T_{q}(\lambda,q,\mbox{\boldmath{$R$}})-\mathbb{E}T_{q}(\lambda_{q},q,\mbox{\boldmath{$\mathcal{R}$}})=$
$\displaystyle\qquad=\left[\frac{1}{n_{\lambda,q}^{\prime}}\sum_{i=q+1}^{n}\frac{b_{q,i}^{2}\lambda
n\eta_{q,i}\log(\lambda n\eta_{q,i})}{(1+\lambda
n\eta_{q,i})^{2}}-\frac{1}{n_{\lambda_{q},q}^{\prime}}\sum_{i=q+1}^{n}\frac{b_{q,i}^{2}\lambda_{q}n\eta_{q,i}\log(\lambda_{q}n\eta_{q,i})}{(1+\lambda_{q}n\eta_{q,i})^{2}}\right]\\{1+o(1)\\}.$

By the fundamental theorem of calculus, we see that

$\left|\mathbb{E}T_{q}(\lambda,q,\mbox{\boldmath{$R$}})-\mathbb{E}T_{q}(\lambda_{q},q,\mbox{\boldmath{$\mathcal{R}$}})\right|\leq\frac{1}{n_{\lambda,q}^{\prime}}\left|\sum_{i=q+1}^{n}\frac{b_{q,i}^{2}\lambda_{q}n\eta_{q,i}\log(\lambda_{q}n\eta_{q,i})}{(1+\lambda_{q}n\eta_{q,i})^{2}}\right|\cdot\left|1-\frac{\lambda}{\lambda_{q}}\right|.$

Under the event $\mathcal{A}$ this is indeed $o(1)$. With this we conclude the
proof.

## References

  * Altman, (1990) Altman, N. (1990).  Kernel smoothing of data with correlated errors.  J. Am. Statist. Ass., 85:749–759.
  * Amato et al., (2017) Amato, U., Antoniadis, A., De Feis, I., and Goude, Y. (2017).  Estimation and group variable selection for additive partial linear models with wavelets and splines.  South African Statist. J., 51(2):235–272.
  * (3) Bickel, P. and Levina, E. (2008a).  Covariance regularization by thresholding.  Ann. Statist., 36:2577–2604.
  * (4) Bickel, P. and Levina, E. (2008b).  Regularized estimation of large covariance matrices.  Ann. Statist., 36:199–227.
  * Cai et al., (2010) Cai, T., Zhang, C.-H., and Zhou, H. (2010).  Optimal rates of convergence for covariance matrix estimation.  Ann. Statist., 38:2118–2144.
  * Chiu, (1989) Chiu, S.-T. (1989).  Bandwidth selection for kernel estimate with correlated noise.  Statist. Probab. Lett., 8:347–354.
  * Chu and Marron, (1991) Chu, C.-K. and Marron, J. (1991).  Comparison of two bandwidth selectors with dependent errors.  Ann. Statist., 19:1906–1918.
  * Durban and Currie, (2003) Durban, M. and Currie, I. (2003).  A note on P-spline additive models correlated errors.  Comput. Statist., 18:89–121.
  * Eggermont and LaRiccia, (2009) Eggermont, P. and LaRiccia, V. (2009).  Maximum penalized likelihood estimation, volume 2 of Springer series in Statistics.  Springer.
  * Eggermont and LaRiccia, (2006) Eggermont, P. P. B. and LaRiccia, V. N. (2006).  Uniform error bounds for smoothing splines.  Lecture Notes-Monograph Series, 51:220–237.
  * Enikeeva et al., (2018) Enikeeva, F., Munk, A., and Werner, F. (2018).  Bump detection in heterogeneous gaussian regression.  Bernoulli, 24(2):1266–1306.
  * Fan et al., (2016) Fan, J., Liao, Y., and Liu, H. (2016).  An overview of the estimation of large covariance and precision matrices.  Econom. J., 19(1).
  * Fix, (1972) Fix, G. J. (1972).  Effects of quadrature errors in finite element approximation of steady state, eigenvalue and parabolic problems.  In Aziz, A., editor, The Mathematical Foundations of the Finite Element Method with Applications to Partial Differential Equations, pages 525 – 556. Academic Press.
  * Fix, (1973) Fix, G. J. (1973).  Eigenvalue approximation by the finite element method.  Advances in Mathematics, 10(2):300 – 316.
  * Gray, (1972) Gray, R. (1972).  On the asymptotic eigenvalue distribution of toeplitz matrices.  IEEE Trans. Inform. Theory, 18(6):725–730.
  * Hall and Keilegom, (2003) Hall, P. and Keilegom, I. V. (2003).  Using difference-based methods for inference in nonparametric regression with time series errors.  J. R. Statist. Soc. B, 65(2):443–456.
  * Hall et al., (1995) Hall, P., Lahiri, S., and Polzehl, J. (1995).  On bandwidth choice in nonparametric regression with both short- and long-range dependent errors.  Ann. Statist., 23:1921–1936.
  * Hart, (1994) Hart, J. (1994).  Automated kernel smoothing of dependent data by using time series cross-validation.  J. R. Statist. Soc. B, 56:529–542.
  * Hart and Wehrley, (1986) Hart, J. and Wehrley, T. (1986).  Kernel regression estimation using repeated measurements data.  J. Am. Statist. Ass., 81(396):1080–1088.
  * Herrmann et al., (1992) Herrmann, E., Gasser, T., and Kneip, A. (1992).  Choice of bandwidth for kernel regression when residuals are correlated.  Biometrika, 79:783–795.
  * Hurvich and Zeger, (1990) Hurvich, C. and Zeger, S. (1990).  A frequency domain selection criterion for regression with autocorrelated errors.  J. Am. Statist. Ass., 85:705–714.
  * Kohn et al., (1992) Kohn, R., Ansley, C., and Wong, C. (1992).  Nonparametric spline regression with autoregressive moving average errors.  Biometrika, 79:335–346.
  * Kotz and Nadarajah, (2004) Kotz, S. and Nadarajah, S. (2004).  Multivariate t-distributions and their applications.  Cambridge University Press.
  * Krivobokova and Kauermann, (2007) Krivobokova, T. and Kauermann, G. (2007).  A note on penalized spline smoothing with correlated errors.  J. Am. Statist. Ass., 102(480).
  * Lee et al., (2010) Lee, Y. K., Mammen, E., and Park, B. U. (2010).  Bandwidth selection for kernel regression with correlated errors.  Statistics, 44(4):327–340.
  * Opsomer et al., (2001) Opsomer, J., Wang, Y., and Yang, Y. (2001).  Nonparametric regression with correlated errors.  Statist. Sci., 16(134–153).
  * Purahmadi, (2011) Purahmadi, M. (2011).  Covariance estimation: the GLM and regularized perspectives.  Statist. Sci., 26:369–387.
  * Rosales Marticorena, (2016) Rosales Marticorena, L. F. (2016).  Empirical Bayesian Smoothing Splines for Signals with Correlated Errors: Methods and Applications.  PhD thesis, Georg-August-Universität Göttingen, http://hdl.handle.net/11858/00-1735-0000-0028-87F9-6.
  * Serra and Krivobokova, (2017) Serra, P. and Krivobokova, T. (2017).  Adaptive empirical Bayesian smoothing splines.  Bayesian Anal., 12(1):219–238.
  * Speckman and Sun, (2003) Speckman, P. and Sun, D. (2003).  Fully Bayesian spline smoothing and intrinsic autoregressive priors.  Biometrika, 90(2):289–302.
  * Szabó et al., (2015) Szabó, B., van der Vaart, A. W., van Zanten, J., et al. (2015).  Frequentist coverage of adaptive nonparametric bayesian credible sets.  Ann. Statist., 43(4):1391–1428.
  * Utreras Diaz, (1980) Utreras Diaz, F. (1980).  Pur ie choix du paramétre d’ajustement dans le lissage par fonctions spline.  Numer. Math., 34:15–28.
  * van der Vaart, (1998) van der Vaart, A. (1998).  Asymptotic Statistics, volume 3 of Cambridge Series in Statistical and Probabilistic Mathematics.  Cambridge University Press.
  * Wang, (1998) Wang, Y. (1998).  Smoothing spline models with correlated random errors.  J. Am. Statist. Ass., 93:341–348.
  * Wood, (2017) Wood, S. (2017).  Generalized Additive Models: An Introduction with R.  Chapman and Hall/CRC.
  * Xiao and Wu, (2012) Xiao, H. and Wu, W. (2012).  Covariance matrix estimation for stationary time series.  Ann. Statist., 40(1):466–493.
  * Yang et al., (2001) Yang, Y. et al. (2001).  Nonparametric regression with dependent errors.  Bernoulli, 7(4):633–655.