Functional Lagged Regression with Sparse Noisy Observations

TL;DR

This paper develops a spectral approach for lagged functional regression with sparse, noisy data, enabling estimation of model dynamics and latent functions from limited observations.

Contribution

It introduces a novel spectral estimation method combined with regularization techniques for functional lagged regression under sparse, noisy sampling conditions.

Findings

Effective spectral density estimation from sparse data

Successful recovery of latent functions via prediction

Validated approach through simulations and meteorological data

Abstract

A functional (lagged) time series regression model involves the regression of scalar response time series on a time series of regressors that consists of a sequence of random functions. In practice, the underlying regressor curve time series are not always directly accessible, but are latent processes observed (sampled) only at discrete measurement locations. In this paper, we consider the so-called sparse observation scenario where only a relatively small number of measurement locations have been observed, possibly different for each curve. The measurements can be further contaminated by additive measurement error. A spectral approach to the estimation of the model dynamics is considered. The spectral density of the regressor time series and the cross-spectral density between the regressors and response time series are estimated by kernel smoothing methods from the sparse observations. The impulse response regression coefficients of the lagged regression model are then estimated by means of ridge regression (Tikhonov regularisation) or PCA regression (spectral truncation). The latent functional time series are then recovered by means of prediction, conditioning on all the observed observed data. The performance and implementation of our methods are illustrated by means of a simulation study and the analysis of meteorological data.

Tables

Table 1. Table 1: The mean square prediction error (MSE) and R 2 superscript 𝑅 2 R^{2} coefficients of determination of each of the model with either truncation or Tikhonov regularisation based estimation of the filter coefficients. Both the MSE’s and the R 2 superscript 𝑅 2 R^{2} coefficients were determined on the test partition of the response { Z t } t = 823 1095 superscript subscript subscript 𝑍 𝑡 𝑡 823 1095 \{Z_{t}\}_{t=823}^{1095}

Model	Truncation		Tikhonov
	MSE	$R^2$	MSE	$R^2$
(E)	541	0.14	481	0.23
(T)	407	0.35	379	0.39
(E+T)	392	0.37	335	0.46

Equations

$$Z_t=a+\sum _{k\in \mathbb{Z}}{\mathcal{B}}_k{X}_{t-k}+e_t,t\in \mathbb{Z},$$

$$Y_{tj}=X_t(x_{tj})+{ϵ}_{tj},j=1,\dots ,N_t,t=1,...,T.$$

$$\sum _{h\in \mathbb{Z}}\Vert R_h^X{\Vert }_{\infty }=\sum _{h\in \mathbb{Z}}\underset{x,y\in [0,1]}{\sup }|R_h^X(x,y)|<\infty ,\sum _{h\in \mathbb{Z}}\Vert {\mathcal{R}}_h^X{\Vert }_1=\sum _{h\in \mathbb{Z}}\text{mbox}trace\left\{\sqrt{({\mathcal{R}}_h^X{)}^{\ast }{\mathcal{R}}_h^X}\right\}<\infty .$$

$$Z_t=a+\sum _{k\in \mathbb{Z}}{\mathcal{B}}_k{X}_{t-k}+e_t,t\in \mathbb{Z},$$

$$\sum _{k\in \mathbb{Z}}\Vert {\mathcal{B}}_k\Vert <\infty .$$

$$Z_t=\sum _{k\in \mathbb{Z}}{\mathcal{B}}_k{X}_{t-k}+e_t,t\in \mathbb{Z},\left(\mathbb{E}{X}_t=0,\mathbb{E}{Z}_t=0\right).$$

$$R_h^{ZX}(x)=\mathbb{E}\left[Z_h{X}_0(x)\right],x\in [0,1],{\mathcal{R}}_h^{ZX}=\mathbb{E}\left[Z_h⟨X_0,\cdot ⟩\right]:\mathcal{H}\to \mathbb{R}.$$

$$f_{\omega }^X(\cdot ,\cdot )=\frac{1}{2\pi }\sum _{h\in \mathbb{Z}}{R}_h^X(\cdot ,\cdot )e^{-i\omega h},{\mathcal{F}}_{\omega }^X=\frac{1}{2\pi }\sum _{h\in \mathbb{Z}}{\mathcal{R}}_h^X{e}^{-i\omega h},$$

$${\mathcal{F}}_{\omega }^X=\sum _{m=1}^{\infty }{\lambda }_m^{\omega }{\phi }_m^{\omega }\otimes {\phi }_m^{\omega }=\sum _{m=1}^{\infty }{\lambda }_m^{\omega }⟨{\phi }_m^{\omega },\cdot ⟩{\phi }_m^{\omega }$$

$$R_h^X(\cdot ,\cdot )={\int }_{-\pi }^{\pi }{f}_{\omega }^X(\cdot ,\cdot )e^{i\omega h}d\omega ,{\mathcal{R}}_h^X={\int }_{-\pi }^{\pi }{\mathcal{F}}_{\omega }^X{e}^{i\omega h}d\omega ,h\in \mathbb{Z}.$$

$${\mathcal{F}}_{\omega }^{ZX}=\frac{1}{2\pi }\sum _{h\in \mathbb{Z}}{\mathcal{R}}_h^{ZX}{e}^{-i\omega h},\omega \in [-\pi ,\pi ],$$

$$f_{\omega }^{ZX}(\cdot )=\frac{1}{2\pi }\sum _{k\in \mathbb{Z}}{R}_h^{ZX}(\cdot )e^{-i\omega h},\omega \in [-\pi ,\pi ].$$

$${\mathcal{B}}_{\omega }=\sum _{h\in \mathbb{Z}}{\mathcal{B}}_k{e}^{-ih\omega },\omega \in [-\pi ,\pi ].$$

$${\mathcal{F}}_{\omega }^{ZX}={\mathcal{B}}_{\omega }{\mathcal{F}}_{\omega }^X,\omega \in [-\pi ,\pi ],$$

$${\mathcal{B}}_k=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{\mathcal{B}}_{\omega }{e}^{ih\omega }d\omega ,k\in \mathbb{Z}.$$

$$\begin{gathered} \left({\hat{c}}_0^{(1)},{\hat{c}}_1^{(1)},{\hat{c}}_2^{(1)}\right)=\underset{c_0^{(1)},c_1^{(1)},c_2^{(1)}}{\mathrm{argmin}}\sum _{t=1}^T\stackrel{[}{j}\ne k]\sum _{j=1}^{N_t}\sum _{k=1}^{N_t}K\left(\frac{x_{tj}-x}{B_R}\right)K\left(\frac{x_{tk}-y}{B_R}\right)\times \\ \times {\left\{G_{0,t}^X(x_{tj},x_{tk})-c_0^{(1)}-c_1^{(1)}(x-x_{tj})-c_2^{(1)}(y-x_{tk})\right\}}^2 \end{gathered}$$

$$\begin{gathered} \left(c_0^{(2)},c_1^{(2)},c_2^{(2)}\right)=\underset{c_0^{(2)},c_1^{(2)},c_2^{(2)}}{\mathrm{argmin}}\stackrel{[}{j}\ne k]\sum _{j=1}^{N_t}\sum _{k=1}^{N_t}K\left(\frac{x_{tj}-x}{B_R}\right)K\left(\frac{x_{tk}-x}{B_R}\right)\times \\ \times {\left\{G_{0,t}^X(x_{tj},x_{tk})-c_0^{(2)}-c_1^{(2)}\Delta (x_{tj},x_{tk})-c_2^{(2)}\Delta (x_{tj},x_{tk}{)}^2\right\}}^2 \end{gathered}$$

$$\Delta (x_{tj},x_{tk})=\mathrm{sign}\left(x_{tk}-x_{tj}\right)\sqrt{{\left(P(x_{tj},x_{tk})-x_{tj}\right)}^2+{\left(P(x_{tj},x_{tk})-x_{tk}\right)}^2}$$

$$\left(c_0^{(3)},c_1^{(3)}\right)=\underset{c_0^{(3)},c_1^{(3)}}{\mathrm{argmin}}\sum _{t=1}^T\sum _{j=1}^{N_t}K\left(\frac{x_{tj}-x}{B_V}\right){\left\{G_{0,t}^X(x_{tj},x_{tj})-c_0^{(3)}-c_1^{(3)}(x-x_{tj})\right\}}^2.$$

$${\hat{\sigma }}^2={\int }_0^1\left({\hat{V}}^X(x)-{\bar{R}}_0^X(x)\right)dx.$$

$${\hat{f}}_{\omega }^X(x,y)=\frac{L}{2\pi }{c}_0^{(4)}\left(\in \mathbb{C}\right)$$

$$\left(c^{(4)}_{0},\,c^{(4)}_{1},\,c^{(4)}_{2}\right)=\operatorname*{arg\,min}_{\left(c^{(4)}_{0},\,c^{(4)}_{1},\,c^{(4)}_{2}\right)\in\mathbb{C}^{3}}\sum_{h=-L}^{L}\frac{1}{\mathcal{N}_{h}}\sum_{t=\max(1,1-h)}^{\min(T,T-h)}\stackrel{{\scriptstyle[}}{{j}}\neq k\text{ if }h=0]{}{\sum_{j=1}^{N_{t+h}}\sum_{k=1}^{N_{t}}}\big{|}G_{h,t}^{X}(x_{t+h,j},x_{tk})e^{-\operatorname{i}h\omega}-\\ -c^{(4)}_{0}-c^{(4)}_{1}(x_{t+h,j}-x)-c^{(4)}_{2}(x_{tk}-y)\big{|}^{2}W_{h}\frac{1}{B_{R}^{2}}K\left(\frac{x_{t+h,j}-x}{B_{R}}\right)K\left(\frac{x_{tk}-y}{B_{R}}\right)$$

$${\hat{R}}_h^X(\cdot ,\cdot )={\int }_{-\pi }^{\pi }{\hat{f}}_{\omega }^X(\cdot ,\cdot )e^{ih\omega }d\omega ,{\hat{\mathcal{R}}}_h^X={\int }_{-\pi }^{\pi }{\hat{\mathcal{F}}}_{\omega }^X{e}^{ih\omega }d\omega .$$

$${\hat{f}}_{\omega }^{ZX}(x)=\frac{L}{2\pi }{c}_0^{(5)}\left(\in \mathbb{C}\right)$$

$$\begin{gathered} \left(c_0^{(5)},c_1^{(5)}\right)=\underset{\left(c_0^{(5)},c_1^{(5)}\right)\in {\mathbb{C}}^2}{\mathrm{argmin}}\sum _{h=-L}^L\sum _{t=\max (1,1-h)}^{\min (T,T-h)}\sum _{j=1}^{N_t}\frac{W_h}{B_C}K\left(\frac{x_{tj}-x}{B_C}\right)\times \\ \times {\left|G_{h,t}^{ZX}(x_{tj})e^{-ih\omega }-c_0^{(5)}-c_1^{(5)}(x_{tj}-x)\right|}^2 \end{gathered}$$

$${\hat{\mathcal{F}}}_{\omega }^X=\sum _{j=1}^{\infty }{\hat{\lambda }}_j^{\omega }{\hat{\phi }}_j^{\omega }\otimes {\hat{\phi }}_j^{\omega },\omega \in [-\pi ,\pi ],$$

$$\sum _{j=1}^{K_{\omega }}\frac{1}{{\lambda }_j^{\omega }}{\hat{\phi }}_j^{\omega }\otimes {\hat{\phi }}_j^{\omega },\omega \in [-\pi ,\pi ],$$

$${\hat{\mathcal{B}}}_{\omega }^{trunc}=\sum _{j=1}^{K_{\omega }}\frac{1}{{\lambda }_j^{\omega }}⟨{\hat{\phi }}_j^{\omega },\cdot ⟩{\hat{\mathcal{F}}}_{\omega }^{ZX}{\hat{\phi }}_j^{\omega },\omega \in [-\pi ,\pi ].$$

$${\left({\hat{\mathcal{F}}}_{\omega }^X+\rho \mathcal{I}\right)}^{-1}=\sum _{j=1}^{\infty }\frac{1}{{\lambda }_j^{\omega }+{\rho }^{\omega }}{\hat{\phi }}_j^{\omega }\otimes {\hat{\phi }}_j^{\omega },\omega \in [-\pi ,\pi ],$$

Full text

62M10,

62M15, 60G10,

autocovariance operator,

lagged regression,

functional data analysis,

nonparametric regression,

spectral density operator,

and

1. 1 Introduction

2. 2 Model and estimation methodology

3. 2.1 Functional time series regression model

4. 2.2 Spectral analysis of functional lagged regression

5. 2.3 Nonparametric estimation in the lagged regression

6. 2.4 Forecasting the response process

7. 3 Asymptotic results

8. 4 Numerical experiments

9. 4.1 Simulation setting

10. 4.2 Estimation procedure

11. 4.3 Evaluation criteria

12. 4.4 Results of numerical experiments

13. 5 Data analysis

14. 6 Discussion of the considered regularisation techniques

15. A Lagged regression with functional response

16. A.1 Functional response model and used notation

17. A.2 Nonparametric estimation of the cross-spectral density and spectral transfer function

18. A.3 Forecasting the response process

19. B Lagged regression with multiple inputs

20. B.1 Joint regression model and its spectral analysis

21. B.2 Nonparametric estimation

22. B.3 Forecasting the response process

23. C Proofs of formal statements

24. C.1 Proof of Proposition 1

25. C.2 Auxiliary Lemmas

26. C.3 Proof of Proposition 2

27. C.4 Proof of Theorem 1

28. C.5 Proof of Theorem 2

1 Introduction

A (lagged) time series regression model is perhaps the most basic –and certainly one of the most and longest studied (Kolmogoroff, 1941; Wiener, 1950)– forms of coupled analysis of two time series. In a general context, given two discrete-time stationary time series $\{X_{t}\}_{t\in\mathbb{Z}}$ and $\{Z_{t}\}_{t\in\mathbb{Z}}$, the input (or regressor) and output (or response), valued in some vector spaces $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$, such a model postulates that

[TABLE]

for some constant $a\in\mathcal{H}_{2}$, a sequence of random disturbances $\{e_{t}\}_{t\in\mathbb{Z}}$ valued in $\mathcal{H}_{2}$ and a sequence of linear mappings $\mathcal{B}_{k}:\mathcal{H}_{1}\rightarrow\mathcal{H}_{2},\,k\in\mathbb{Z}$. This linear coupling would be the typical dependence model, for instance, if $\{(X_{t},Z_{t})\}_{t\in\mathbb{Z}}$ were a jointly Gaussian stationary process in $\mathcal{H}_{1}\times\mathcal{H}_{2}$, and is also known as a time-invariant linearly filtered time series model. The estimation problem is then to estimate the unknown transformations $\{\mathcal{B}_{k}\}_{k\in\mathbb{Z}}$ given the realisation of a finite stretch of the joint series $\{(Z_{t},X_{t})\}_{t\in\mathbb{Z}}$ (a problem also known as system identification, particularly in signal processing).

This problem is very well understood and has been extensively studied in the classical context where the spaces $\mathcal{H}_{j},j=1,2,$ coincides with the Euclidean spaces of potentially different dimensions (Brillinger, 1981; Priestley, 1981; Shumway and Stoffer, 2000). Nevertheless, generalising these results to the case where either space, and particularly the regressor space $\mathcal{H}_{1}$, may be an infinite dimensional vector space is far from straightforward, and has been comparatively less studied. The difficulty in this case is that one needs to manipulate operations involving the inverses of compact or even trace-class operators, which fail to exist boundedly on the entire codomain. Consequently, analysis of such models requires drawing on tools from functional analysis, and developing novel methodology that incorporates suitable regularising schemes. Such a setting has only recently been considered for functional time series regression models, for example by Hörmann, Kidziński and Kokoszka (2015), who treat the problem of estimation of the filter coefficients by means of spectral truncation regularisation (PCA regression), and Pham and Panaretos (2018), who deduce convergence rates for the estimated coefficients when using Tikhonov regularisation (ridge regression).

While this may seem to be an artificial abstraction at first sight, such infinite-dimensional time series are becoming increasingly prominent for applications. Indeed the abstraction of an infinite dimensional (Hilbert) space, captures the scenario where the value of a series at each time is a square-integrable function, e.g. a curve. Examples include time series of DNA minicircles evolving in solution (seen as a time series of closed curves in 3D indexed by discrete time (Tavakoli and Panaretos, 2016)) or the data constructed by dividing a continuously observed scalar time series into segments of an obvious periodicity, usually days. Examples of the latter form are particularly prominent in environmental applications, for example the analysis of particulate matter atmospheric pollution (Hörmann and Kokoszka, 2010; Hörmann, Kidziński and Hallin, 2015; Hörmann, Kokoszka and Nisol, 2018; Aue, Norinho and Hörmann, 2015), traffic data modelling (Klepsch, Klüppelberg and Wei, 2017), or financial applications of intra-day trading (Müller, Sen and Stadtmüller, 2011; Kokoszka et al., 2017). Another promising financial application of functional time series emerges in the yield curve modelling (Hays et al., 2012; Kowal, Matteson and Ruppert, 2017a; Sen and Klüppelberg, 2019).

We focus here primarily on the case $\mathcal{H}_{2}=\mathbb{R}$, i.e. when the response process is scalar. This is indeed the case that is most often studied in the literature, due in large part due to the many examples it covers, but also because it is the simplest version of the problem that captures the essence of higher generality: the difficulty in estimating the filter coefficients lies in the ill-possessedness of the spectral density operator inversion that requires regularisation. For a scenario involving $\mathcal{H}_{2}$ being infinite dimensional, i.e. with a functional response, see Hörmann, Kidziński and Kokoszka (2015) or Appendix A.

In practice, one can never observe the regressor time series $\{X_{t}\}_{t\in\mathbb{Z}}$ in its “Platonic” continuum form. For example, if $\mathcal{H}_{1}$ is a space of functions on $[0,1]$ and each $X_{t}:[0,1]\rightarrow\mathbb{R}$ is a curve, then one might be able to observe evaluations of $X_{t}(\cdot)$ at various locations of its domain. In some cases, the sampled locations are sufficiently dense, and the measurement instrumentation sufficiently precise to be free of any additional noise contamination, so that one can disregard the effects of sampling. This is essentially the approach taken in Hörmann, Kidziński and Kokoszka (2015) and Pham and Panaretos (2018) where the regressors are treated as being fully observed as elements of $\mathcal{H}_{1}=L^{2}[0,1]$. However, it may well happen that $\{X_{t}\}_{t\in\mathbb{Z}}$ is only measured at few and randomly positioned locations in its domain, indeed varying with time, and that the measurements are themselves contaminated by noise. That is, instead of observing $X_{t}(u)$ for all $u$ in $[0,1]$, we instead observe

[TABLE]

for a sequence of point processes $(x_{t1},...,x_{tN_{t}}),t=1,\dots,T$, independent in time, and a white noise measurement error sequence. This observation scheme is illustrated on Figure 1 and is exhibited, for example, when dealing with fair weather atmospheric electricity data (Rubín and Panaretos, 2020; Tammet, 2009).

Such a setting escapes both the methods and the theory developed at the level of continuum, and instead requires methodology that accounts for the observation scheme, as well as theory that incorporates the latent/emission processes explicitly. The purpose of this paper is precisely to construct such a methodology and develop the associated asymptotic theory.

Our approach consists in estimating the complete space-time covariance structure of the data by estimating the spectral density operator of the regressor time series, and the cross-spectral density operator between the regressor and response time series. Our estimators are based on the kernel smoothing methods as discussed for example in (Yao, Müller and Wang, 2005a, b; Hall, Müller and Wang, 2006; Li and Hsing, 2010) and extended to the time series case in Rubín and Panaretos (2020). Once these are estimated, one obtains estimating equations whose solution yields the estimators of the filter coefficients. The solution of the estimating equations, however, comprises an ill-posed inverse problem, hence regularisation is required. We offer two regularisation techniques, spectral truncation regularisation (Hörmann, Kidziński and Kokoszka, 2015) and Tikhonov regularisation (Pham and Panaretos, 2018). The forecasting of the response process is then implemented by first predicting the latent functional regressor data (using their estimated spectral characteristics) and then plugging-in these predictions into the estimated lagged regression model.

Sparsely observed functional time series have only recently received attention (Kowal, Matteson and Ruppert, 2017a, b; Rubín and Panaretos, 2020; Sen and Klüppelberg, 2019), and our results appear to be the first in the context of the lagged regression model where the regressor process is functional. A related problem of dynamic function-on-scalar regression was studied by Kowal (2018) by means of Bayesian factor models.

Our presentation is organised as follows. Section 2 defines the framework of sparsely observed functional time series as well as the functional lagged regression model and its analysis in the spectral domain. The section explains the estimation methodology for the model components, the spectral and the cross-spectral densities, and the regularised estimator of the filter coefficients of the lagged regression. Furthermore, the forecasting algorithm is introduced. Section 3 presents the asymptotic results of the proposed method: consistency and convergence rates are established. Section 4 probes the finite sample properties of the proposed methodology on a simulation study. As a proof-of-concept illustration, we included an example featuring meteorological data in Section 5 of the paper. The properties of the two considered regularisation methods are summarised in Section 6. We present the extension of the proposed methodology for lagged regression models with a functional response in Appendix A and for lagged regression models with multiple regressors time series in Appendix B. The proofs of the formal statements presented in Appendix C.

2 Model and estimation methodology

2.1 Functional time series regression model

The classical theory of functional data analysis views its data, random curves, as random elements in a suitable Hilbert space. Throughout our presentation we assume this space to be the Hilbert space of real square integrable function, $\mathcal{H}_{1}\equiv\mathcal{H}=L^{2}([0,1])$, equipped with the inner product $\langle\cdot,\cdot\rangle$ and the norm $\|\cdot\|$. The functional time series is then a sequence of random elements in such space indexed by the integer variable $t\in\mathbb{Z}$, interpreted as times (for example days), and denoted as $\{X_{t}\}\equiv\{X_{t}(\cdot)\}_{t\in\mathbb{Z}}$. The argument of the individual functions, for the most part of this paper denoted as $x\in[0,1]$, can be interpreted as either the spatial location or intra-day time (for the sequence of intra-day data order day-by-day into a sequence). The basic assumption we impose is:

1. (A1)

The functional time series $\{X_{t}\}$ is assumed to be square-integrable, i.e. $\mathbb{E}\|X_{t}\|^{2}<\infty,\,t\in\mathbb{Z}$, and second-order stationary in the variable $t$. Furthermore, the sample paths (trajectories) of $\{X_{t}\}$ are smooth.

The functional time series $\{X_{t}\}$ is assumed to be only latent and therefore not accessible directly. Denote $N_{t}$ the number of measurements available for the random curve at time $t=1,\dots,T$ where $T\in\mathbb{N}$ is the time of the latest curve with at least one measurement. We denote $Y_{tj}$ the $j$-th measurement on the $t$-th curve at spatial position $x_{tj}\in[0,1]$ where $j=1,\dots,N_{t}$ and $t=1,\dots,T$. The measurements are assumed to be contaminated by an additive measurement errors $\epsilon_{tj},\,j=1,\dots,N_{t},\,t=1,\dots,T$ which constitute an independent identically distributed zero-mean array of scalar random variables with variance $\sigma^{2}>0$. The number of measurement locations $\{N_{t}\}$, the measurement errors $\{\epsilon_{tj}\}$, and the underlying functional time series $\{X_{t}\}$ are assumed to be independent. The above outlined sparse observation regime leads to the model equation

[TABLE]

The assumption (A1) allows us to define the mean function as $\mu(x)=\mathbb{E}X_{0}(x),\,x\in[0,1],$ and the lag-$h$ autocovariance kernel $R_{h}^{X}(x,y)=\mathbb{E}\left[\left(X_{h}(x)-\mu(x)\right)\left(X_{0}(y)-\mu(y)\right)\right],\,x,y\in[0,1],$ for $h\in\mathbb{Z}$. The integral operator $\mathscr{R}_{h}^{X}$ associated with this kernel, called the lag-$h$ autocovariance operator, is defined by the right integration $\left(\mathscr{R}_{h}^{X}f\right)(x)=\int_{0}^{1}R_{h}^{X}(x,y)f(y)\operatorname{d}\!{y},\,f\in L^{2}([0,1])$. The autocovariance kernels and operators are assumed to be summable in the supremum norm (denoted $\|\cdot\|_{\infty}$) and the nuclear norm (denoted $\|\cdot\|_{1}$) respectively,

1. (A2)

[TABLE]

Throughout the paper we assume that the response time series $\{Z_{t}\}$ is scalar. We make this assumption to facilitate the presentation while keeping in mind the fact that the scalar case already demonstrates the difficulty of the ill-posed estimation problem and the requirement of regularisation. The extension to a functional response time series is presented in Appendix B.

The functional lagged regression model with a scalar response time series $\{Z_{t}\}$ is given by the equation

[TABLE]

where $a\in\mathbb{R}$ is a constant, called the intercept; $\mathcal{B}_{k}:\mathcal{H}\to\mathbb{R},\,k\in\mathbb{Z},$ are fixed continuous linear functionals called the filter coefficients, with Riesz-representers $\{b_{k}\}_{k\in\mathbb{Z}}$ (i.e. $\mathcal{B}_{k}f=\langle f,b_{k}\rangle$ for $f\in\mathcal{H}$); and $\{e_{t}\}_{t\in\mathbb{Z}}$ is a sequence of independent identically distributed zero-mean real random variables with variance $\tau^{2}>0$. The time series $\{X_{t}\}$, the sampling counts $\{N_{t}\}$, the observation locations $\{x_{tj}\}$, and the errors $\{\epsilon_{tj}\}$ and $\{e_{t}\}$ are assumed to be independent. The filter coefficients are assumed to be summable in the norm

1. (A3)

[TABLE]

where given any continuous linear functional $\mathcal{G}:\mathcal{H}\to\mathbb{R}$ the notation $\|\mathcal{G}\|$ will be used interchangeably with $\|g\|$ for $g$ the Riesz-representer of $\mathcal{G}$, i.e. $\mathcal{G}f=\langle g,f\rangle$.

We make here a few simplifying assumptions concerning the first order structure of the data, to keep focus on the more challenging second order structure. Firstly, we assume that the regressor time series $\{X_{t}\}$ is centred, i.e. $\mu=0$. If that was not the case, the mean function $\mu(\cdot)$ could be estimated by the local-linear kernel smoother (Yao, Müller and Wang, 2005a; Rubín and Panaretos, 2020) and all the sparse observations centred. Secondly, we assume that the intercept is null, i.e. $a=0$. Otherwise it could be estimated using the relation $a=\mathbb{E}\left[Z_{0}\right]-\sum_{k\in\mathbb{Z}}\mathcal{B}_{k}\mu$ where $\mathbb{E}\left[Z_{0}\right]$ is estimated by the sample mean and the coefficients $\mathcal{B}_{k}$ using the methods of this paper. As a consequence of these assumptions, the response time series is also centred, i.e. $\mathbb{E}\left[Z_{0}\right]=0$, leading to the model

[TABLE]

Finally, we define the lag-$h$ cross-covariance kernel and the lag-$h$ cross-covariance operator between the scalar time series $\{Z_{t}\}$ and the functional time series $\{X_{t}\}$, for $h\in\mathbb{Z}$, as

[TABLE]

2.2 Spectral analysis of functional lagged regression

Under the assumptions (A1) and (A2), the sums

[TABLE]

converge in the supremum and the trace-class norm respectively, define the spectral density kernel and the spectral density operator at frequency $\omega\in[-\pi,\pi]$ respectively (Panaretos and Tavakoli, 2013a). The spectral density operator $\mathscr{F}^{X}_{\omega}$ is a non-negative self-adjoint trace-class operator for all $\omega\in[-\pi,\pi]$. Hence it admits the spectral representation

[TABLE]

where $\otimes$ is the tensor product in $\mathcal{H}$, defined by the inner product formula on the right-hand side of (2.4). The elements of the sequence $\lambda_{1}^{\omega}\geq\lambda_{2}^{\omega}\geq\dots\geq 0$ are called the harmonic eigenvalues and the corresponding functions $\{\varphi_{m}^{\omega}\}_{m=1}^{\infty}$ are called the harmonic eigenfunctions.

Furthermore, the lagged autocovariance kernels and operators can be recovered by the inversion formula (Panaretos and Tavakoli, 2013a) that holds in the supremum and the nuclear norm, respectively:

[TABLE]

Hörmann, Kidziński and Kokoszka (2015) defined the cross-spectral density operator between $\{Z_{t}\}_{t\in\mathbb{Z}}$ and $\{X_{t}\}_{t\in\mathbb{Z}}$ by the formula

[TABLE]

and showed that this sum converges in the Hilbert-Schmidt norm, provided the autocovariance function $\mathscr{R}_{h}^{X}$ is summable in the Hilbert-Schmidt norm and the sequence of filter coefficient operators are also summable in the operator norm. Since our assumptions (A1) — (A3) are stronger and we consider only scalar response, (2.6) is well defined in the setting of our paper.

Likewise, one defines the cross-spectral density kernel between $\{Z_{t}\}_{t\in\mathbb{Z}}$ and $\{X_{t}\}_{t\in\mathbb{Z}}$ by

[TABLE]

The sum on the right-hand side of (2.7) converges in the supremum norm.

Further, Hörmann, Kidziński and Kokoszka (2015) introduced the frequency response operator

[TABLE]

and obtained the relation between the spectral density operators, cross-spectral density operators, and the frequency response operators

[TABLE]

which provides the basis for the estimation of the filter coefficients introduced in the next section. In our case of scalar response, $\mathscr{B}_{\omega}$ is in fact a functional, i.e. $\mathscr{B}_{\omega}:\mathcal{H}\to\mathbb{R}$. The filter coefficients $\{\mathcal{B}_{k}\}_{k\in\mathbb{Z}}$ can be recovered by the formula

[TABLE]

2.3 Nonparametric estimation in the lagged regression

The nonparametric estimation in functional models given sparse noisy measurements has been frequently based on local-polynomial kernel smoothers (Yao, Müller and Wang, 2005a; Hall, Müller and Wang, 2006; Li and Hsing, 2010; Rubín and Panaretos, 2020). Throughout the paper we use the Epanechnikov kernel $K(v)=\frac{3}{4}(1-v^{2})$ for $v\in[-1,1]$, and [math] otherwise, because it is known to minimise the mean square error in the classical local-polynomial smoothing problem (Fan and Gijbels, 1996, Section 3.2.6). Moreover, the Epanechnikov kernel’s bounded support speeds up the numerical calculations. Nevertheless, the choice of the kernel function is not crucial for the statistical performance of the smoothing estimators (Fan and Gijbels, 1996, Section 3.2.6) and any other commonly used kernel function could be adopted.

We start by defining the “raw” covariances of $\{X_{t}\}$ as $G_{h,t}^{X}\left(x_{t+h,j},x_{tk}\right)=Y_{t+h,j}Y_{tk}$ for $h=-T,\dots,T,\,j=1,\dots,N_{t+h},\,k=1,\dots,N_{t}$. The lag-[math] covariance kernel $R_{0}^{X}(\cdot,\cdot)$ is estimated by the local-linear surface smoother, setting $\hat{R}^{X}_{0}(x,y)=\hat{c}_{0}^{(1)}$ where

[TABLE]

for $x,y\in[0,1]$ and where $B_{R}>0$ is the bandwidth parameter for the surface smoother.

In the next step we aim to estimate the measurement error variance $\sigma^{2}$ by the approach suggested by (Yao et al., 2003; Yao, Müller and Wang, 2005a) for which we need the two following ingredients: the estimator of the diagonal of the lag-[math] covariance kernel of $\{X_{t}\}$ with and without the measurement noise contamination. Firstly, we estimate the diagonal of $R^{X}_{0}$ without the measurement noise contamination by the local-quadratic smoother along the direction perpendicular to the diagonal. For $x\in[0,1]$ we set $\bar{R}_{0}^{X}(x)=\bar{c}_{0}$ where

[TABLE]

where $\Delta(x_{tj},x_{tk})$ is the distance of the point $(x_{tj},x_{tk})$ from the diagonal equipped with the positive sign if the point $(x_{tj},x_{tk})$ is above the diagonal, and with the negative sign if below. Formally

[TABLE]

where $\operatorname{sign}(\cdot)\in\{-1,0,1\}$ is the sign function and $P(x_{tj},x_{tk})$ is the first coordinate of the point $(x_{tj},x_{tk})$ projected onto the diagonal of $[0,1]^{2}$.

Secondly, we estimate the function $x\mapsto R^{X}_{0}(x,x)+\sigma^{2},\,x\in[0,1]$, i.e. the noise contaminated diagonal of the lag-[math] covariance operator. For $x\in[0,1]$ and a bandwidth parameter $B_{V}>0$, we use the local-linear line smoother and set $\hat{V}^{X}(x)=c^{(3)}_{0}$ where

[TABLE]

Having the estimates $\bar{R}^{X}_{0}(\cdot)$ and $\hat{V}^{X}(\cdot)$, the measurement error variance $\sigma^{2}$ is estimated by integrating the difference

[TABLE]

In case the right-hand side of (2.10) is negative, it is recommended to replace it with a small positive number (Yao, Müller and Wang, 2005a).

The spectral density kernels $\{f^{X}_{\omega}\}_{\omega\in[-\pi,\pi]}$ are estimated by Bartlett’s approach (Hörmann, Kidziński and Hallin, 2015; Rubín and Panaretos, 2020), weighing down higher order lags using Barlett’s (triangular) weights defined as $W_{h}=(1-|h|/L)$ for $|h|<L$ and [math] otherwise. The parameter $L\in\mathbb{N}$ is called Bartlett’s span parameter and controls the amount of regularisation involved in the spectral density estimation. For a fixed $\omega\in[-\pi,\pi]$, the spectral density kernel at frequency $\omega$ and $(x,y)\in[0,1]^{2}$ is estimated as

[TABLE]

where $c^{(4)}_{0}\in\mathbb{C}$ is obtained by minimising the following weighted sum of squares

[TABLE]

where $\mathcal{N}_{h}=(T-|h|)(\bar{N})^{2}$ for $h\neq 0$, $\mathcal{N}_{0}=T(\overline{N^{2}}-\bar{N})$, and where $\bar{N}=(1/T)\sum_{t=1}^{T}N_{t}$ and $\overline{N^{2}}=(1/T)\sum_{t=1}^{T}N_{t}^{2}$.

The spectral density estimate $\hat{f}^{X}_{\omega}(\cdot,\cdot)$ and its operator counterpart $\hat{\mathscr{F}}^{X}_{\omega}$ enable us to estimate all autocovariance kernels and operators, thus to estimate the complete second-order structure of $\{X_{t}\}$, by the inversion formula

[TABLE]

In the following paragraphs we extend these kernel smoothing techniques to estimate the cross-spectral density between the response time series $\{Z\}_{t\in\mathbb{Z}}$ and the regressor time series $\{X_{t}\}_{t\in\mathbb{Z}}$.

Define the “raw” lag-$h$ cross covariances $G^{ZX}_{h,t}(x_{tj})=Z_{t+h}Y_{tk}$ where $h=-T+1,\dots,T-1$, $t=\max(1,T-h),\dots,\min(T,T-h)$, and $j=1,\dots,N_{t}$. Using the local linear kernel smoothing techniques, we estimate the cross-spectral density at frequency $\omega\in[-\pi,\pi]$ and at $x\in[0,1]$ as

[TABLE]

where $c^{(5)}_{0}$ is realised as the minimiser of the following weighted sum of squares

[TABLE]

where $B_{C}>0$ is a bandwidth parameter.

The solutions to all above least square optimisation problems can be found explicitly by using a standard argument in local-polynomial regression (Fan and Gijbels, 1996, Section 3.1) or (Rubín and Panaretos, 2020, Section B.2). Moreover, the solutions to the spectral density estimator (2.12) and (2.15) depend on a handful of terms independent of the frequency $\omega\in[-\pi,\pi]$, that can be precalculated, and multiplication by complex exponentials. This allows a computationally feasible evaluation even on a fine grid of frequencies.

Once the estimates of the spectral density kernels $\{\mathscr{F}^{X}_{\omega}\}_{\omega\in[-\pi,\pi]}$ and the cross-spectral density $\{\mathscr{F}^{ZX}_{\omega}\}_{\omega\in[-\pi,\pi]}$ have been constructed, we focus attention on the estimation of the frequency response operators $\{\mathscr{B}_{\omega}\}_{\omega\in[-\pi,\pi]}$. Heuristically, from relation (2.8), we would like to write $\mathscr{B}_{\omega}=\mathscr{F}^{ZX}_{\omega}\left(\mathscr{F}^{X}_{\omega}\right)^{-1},\>\omega\in[-\pi,\pi]$. This formula is indeed only heuristic because the operator $\mathscr{F}^{X}_{\omega}$, being trace class, is not boundedly invertible. The same issue is present also for its empirical counterpart $\hat{\mathscr{F}}^{ZX}_{\omega}\left(\hat{\mathscr{F}}^{X}_{\omega}\right)^{-1}$. Therefore, to achieve consistent estimation, a regularisation of the inverse $\left(\hat{\mathscr{F}}^{X}_{\omega}\right)^{-1}$ is required.

Being a self-adjoint trace class operator, $\hat{\mathscr{F}}^{X}_{\omega}$ admits the spectral representation

[TABLE]

which can be viewed as the empirical version of (2.4). The difficulty in inverting $\hat{\mathscr{F}}^{X}_{\omega}$ can be seen from the fact that $\sum_{j}\hat{\lambda}_{j}^{\omega}=\mathrm{trace}\{\hat{\mathscr{F}}^{X}_{\omega}\}<\infty$, implying that $\hat{\lambda}_{j}^{\omega}$ decays at least as fast as $j^{-(1+\delta)}$, $\delta>0$. It is the small values of $\lambda_{j}^{\omega}$ that cause problems and there are two classical strategies to overcome the issue: spectral truncation and the Tikhonov regularisation.

Spectral truncation. The inverse $\left(\hat{\mathscr{F}}^{X}_{\omega}\right)^{-1}$ is replaced by

[TABLE]

where $K_{\omega}\in\mathbb{N},\,\omega\in[-\pi,\pi],$ is the spectral truncation parameter that needs to grow to infinity sufficiently slowly to allow for the consistency. It may or may not depend on the frequency $\omega\in[-\pi,\pi]$.

The estimator of the spectral transfer function becomes

[TABLE]

We opt to implement the spectral truncation by relying on eigenvalue thresholding approach (Hörmann, Kidziński and Kokoszka, 2015) where we implement the eigenvalue threshold selection by cross-validation (more in Section 4). 2. 2.

Tikhonov regularisation. Here, the inverse of $\hat{\mathscr{F}}^{X}_{\omega}$ is replaced by

[TABLE]

where $\mathcal{I}$ is the identity operator on $\mathcal{H}$ and the Tikhonov regularisation parameter $\rho>0$ tends to zero as $T\to\infty$ slowly enough to allow for consistency. Even though the parameter $\rho$ may, in general, depend on $\omega$ we carry out further analysis with the frequency independent parameter. We do so because in the implementation (more in Section 4) we select the tuning parameter $\rho$ using the cross-validation where it is feasible to optimize over a single (frequency independent) tuning parameter.

The estimator of the spectral transfer function becomes

[TABLE]

This form of regularisation adopted and studied by Pham and Panaretos (2018).

Once the estimators of the spectral transfer operator $\hat{\mathscr{B}}_{\omega}$ have been constructed by either of the above regularisation techniques, the filter coefficients are estimated by

[TABLE]

2.4 Forecasting the response process

In Section 2 we assumed the data to be available up to time $T$ for both the regressor time series $\{X_{t}\}$ as well as the response time series $\{Z_{t}\}$. It may very well happen, though, that the measurement of the response variable is terminated at an earlier time $S$ where $1<S<T$ but the measurements of the regressor time series $\{X_{t}\}$ are available until time $T$. In this case it is of interest to forecast the unobserved values of $Z_{S+1},\dots,Z_{T}$. It turns out that the forecast of $Z_{S+1},\dots,Z_{T}$ is surprisingly straightforward, once the model dynamics have been estimated. We present our forecasting method firstly by assuming the model dynamics are known, and then we plug-in its estimates.

Since the filter coefficients satisfy the assumption (A3), their norms converge to zero as $|k|\to\infty$. Choose a constant $M\in\mathbb{N}$ such that $\|\mathcal{B}_{k}\|$ is negligible for $|k|>M$. Denote $\mathbb{X}=[X_{-M+1},\dots,X_{T+M}]\in\mathcal{H}^{T+2M}$ the random element representing “stacked” curves of the latent regressors time series.

The vector $\mathbb{Y}=(Y_{11},\dots,T_{1N_{1}},\dots,Y_{T,1},\dots,Y_{T,N_{T}})\in\mathbb{R}^{\mathscr{N}_{1}^{T}}$ consists of all observed data of the regressor time series where $\mathscr{N}_{1}^{T}=\sum_{t=1}^{T}N_{t}$ is the total number of observations up to time $T$. The measurement errors $\{\epsilon_{tj}\}$ are stacked into a vector denoted $\mathcal{E}\in\mathbb{R}^{\mathscr{N}_{1}^{T}}$. Further define the evaluation operators $H_{t}:\mathcal{H}\to\mathbb{R}^{N_{t}},\xi\mapsto(\xi(x_{t1}),\dots,\xi(x_{tN_{t}}))$ for each $t=1,\dots,T$ and the “joint” evaluation operator $\mathbb{H}:\mathcal{H}^{T+2M}\to\mathbb{R}^{\mathscr{N}_{1}^{T}},[\xi_{-M+1},\dots,\xi_{T+M}]\mapsto(H_{1}\xi_{1},\dots,H_{T}\xi_{T})$. Hence the observation scheme (2.1) becomes $\mathbb{Y}=\mathbb{H}\mathbb{X}+\mathcal{E}$.

By the symbol $\Pi(\cdot|\mathbb{Y})$ we denote the best linear unbiased predictor of the term represented by the dot, given the data $\mathbb{Y}$. A key observation that simplifies the forecasting is the fact that we may predict the latent functional data first, and then plug them into the filter coefficients. The following proposition summarises this assertion formally.

Proposition 1**.**

The best linear unbiased predictor of $Z_{s}$ given $\mathbb{Y}$, denoted as $\Pi(Z_{s}|\mathbb{Y})$, is equivalent to constructing the best linear unbiased predictors of $X_{t}$ given $\mathbb{Y}$, denoted as $\Pi(X_{t}|\mathbb{Y})$, for all $t\in\mathbb{Z}$ and then applying the filter coefficients $\{\mathcal{B}_{k}\}_{k\in\mathbb{Z}}$ to these predictions:

[TABLE]

In the following paragraph we explain how to construct the predictors $\Pi\left(X_{t}|\mathbb{Y}\right),\,t=-M+1,\dots,T+M$, following Rubín and Panaretos (2020). First, we comment on the distributional properties of the above defined terms. The random element $\mathbb{X}$ inherits the second order structure of $\{X_{t}\}$, thus $\mathbb{E}\mathbb{X}=0$ and

[TABLE]

understood as a linear operator acting on $\mathcal{H}^{T+2M}$ and where $*$ denotes the adjoint operator. Moreover, due to independence of the measurement errors, $\textrm{Var}(\mathcal{E})=\sigma^{2}I_{\mathscr{N}_{1}^{T}}$ where $I_{\mathscr{N}_{1}^{T}}$ is the identity matrix of size $\mathscr{N}_{1}^{T}\times\mathscr{N}_{1}^{T}$.

The best linear unbiased predictor of $\mathbb{X}$ given $\mathbb{Y}$ is then given by

[TABLE]

where the term $\mathbb{H}\mathbb{S}\mathbb{H}^{*}+\sigma^{2}I_{\mathcal{N}_{1}^{T}}$ is a matrix of size $\mathscr{N}_{1}^{T}\times\mathscr{N}_{1}^{T}$ and its inverse is well conditioned since $\sigma^{2}>0$. The best linear predictor of each functional datum $X_{t},t=-M+1,\dots,T+M$ given the observed data $\mathbb{Y}$, denoted as $\Pi(X_{t}|\mathbb{Y})$, is then given by the projection $\Pi(X_{t}|\mathbb{Y})=P_{t}\Pi(\mathbb{X}_{T}|\mathbb{Y})$ where $P_{t}:\mathcal{H}^{T+2M}\to\mathcal{H},[\xi_{-M+1},\dots,\xi_{T+M}]\mapsto\xi_{t}$ is the projection operator for $t=-M+1,\dots,T+M$.

The predictor (2.21) requires the knowledge of the unknown dynamics of the regressor time series through the autocovariance operators (2.20) as well as the measurement error variance $\sigma^{2}$. Instead of these parameters we plug-in the estimated counterparts (2.13) and (2.10) and denote these estimated predictors as $\hat{\Pi}(\cdot|\mathbb{Y})$.

In summary, the forecasting algorithm consists of the following steps:

From the measurements $\mathbb{Y}$ realised on the regressor time series $\{X_{t}\}$ estimate the spectral density $\{\hat{\mathscr{F}}^{X}_{\omega}\}_{\omega\in[-\pi,\pi]}$ and the measurement error variance $\hat{\sigma}^{2}$. Using the formula (2.13), integrate the estimated spectral density to obtain the complete space time covariance $\{\hat{\mathscr{R}}^{X}_{h}\}_{h\in\mathbb{Z}}$ of the regressor time series $\{X_{t}\}$. 2. 2.

From the measurements $\mathbb{Y}$ and the observed response times series $Z_{1},\dots,Z_{S}$, estimate the cross-spectral density $\{\hat{\mathscr{F}}^{ZX}_{\omega}\}_{\omega\in[-\pi,\pi]}$. Using either truncation regularisation (2.16) or Tikhonov regularisation (2.17), estimate the spectral transfer function $\{\mathscr{B}_{\omega}\}_{\omega\in[-\pi,\pi]}$. By means of formula (2.18) or (2.19), depending on the regularisation scheme, integrate the spectral transfer function to obtain the filter coefficients $\{\hat{\mathcal{B}}_{k}^{trunc}\}_{k\in\mathbb{Z}}$ or $\{\hat{\mathcal{B}}_{k}^{Tikh}\}_{k\in\mathbb{Z}}$. 3. 3.

Choose $M$ such that the estimated filter coefficients $\hat{\mathcal{B}}_{k}^{trunc}$ (or $\hat{\mathcal{B}}_{k}^{Tikh}$) are negligible for $|k|>M$. Using the methodology explained at the beginning of this section, construct the prediction of the latent functional data $\hat{\Pi}(X_{-M+1}|\mathbb{Y}),\dots,\hat{\Pi}(X_{T+M}|\mathbb{Y})$. 4. 4.

For each $s=S+1,\dots,T$, construct the forecast $\hat{\Pi}(Z_{s}|\mathbb{Y})=\sum_{k=-M}^{M}\hat{\mathcal{B}}_{k}^{trunc}\hat{\Pi}({X}_{s-k}|\mathbb{Y})$, in the case of spectral truncation or $\hat{\Pi}(Z_{s}|\mathbb{Y})=\sum_{k=-M}^{M}\hat{\mathcal{B}}_{k}^{Tikh}\hat{\Pi}({X}_{s-k}|\mathbb{Y})$), in the case of Tikhonov regularisation.

3 Asymptotic results

In this section we turn to establishing consistency of the proposed estimators of the spectral density (2.11) and the cross-spectral density (2.14), and the regularised estimators of the filter coefficients (2.18) and (2.19). A key requirement for consistency is to control the temporal dependence of $\{X_{t}\}$. Here we employ cumulant mixing conditions which have been successfully used in functional time series before (Panaretos and Tavakoli, 2013a; Rubín and Panaretos, 2020). Alternatives to the cumulant assumptions are $L^{p}$-$m$-approximability (Hörmann and Kokoszka, 2010; Hörmann, Kidziński and Hallin, 2015; Hörmann, Kidziński and Kokoszka, 2015) and strong mixing conditions (Rubín and Panaretos, 2020).

Bellow we list the assumptions that we will made use of to study the asymptotics of the spectral density estimator $\{\hat{f}_{\omega}^{X}\}_{\omega\in[-\pi,\pi]}$ and the cross-spectral density estimator $\{\hat{f}_{\omega}^{ZX}\}_{\omega\in[-\pi,\pi]}$.

1. (B1)

The number of measurements $N_{t}$ at time $t$ are independent random variables, identically distributed with a random variable $N$ satisfying $N\geq 0$, $\mathbb{E}\left[N\right]<\infty$ and $\operatorname{\mathbb{P}}(N>1)>0$. 2. (B2)

The measurement locations $x_{tj},j=1,\dots,N_{t},t=1,\dots,T$ are independent random variables generated from the density $g(\cdot)$ and are independent of the number of measurements $\{N_{t}\}_{t=1,\dots,T}$. The density $g(\cdot)$ is assumed to be twice continuously differentiable and strictly positive on $[0,1]$. 3. (B3)

The autocovariance kernels, $R_{h}(\cdot,\cdot)$, are twice continuously differentiable on $[0,1]^{2}$ for each $h\in\mathbb{Z}$. Moreover, the terms

[TABLE]

are uniformly bounded in $h$ for all combinations of $\alpha_{1},\alpha_{2}\in\mathbb{N}_{0}$ satisfying $\alpha_{1}+\alpha_{2}=2$. 4. (B4)

The 4-th order cumulant kernel of $\{X_{t}\}$ is summable in the supremum norm

[TABLE]

where

[TABLE]

for $t_{1},t_{2},t_{3},t_{4}\in\mathbb{N}$ and $x_{t_{1}},x_{t_{2}},x_{t_{3}},x_{t_{4}}\in[0,1]$. 5. (B5)

The sequence $\{R^{X}_{h}\}_{h\in\mathbb{Z}}$ satisfies the weak dependence condition

[TABLE]

Finally, we state assumptions on the bandwidth parameters $B_{R},B_{V},B_{C}$ and the Bartlett span $L$:

1. (B6)

[TABLE] 2. (B7)

[TABLE] 3. (B8)

[TABLE] 4. (B9)

[TABLE]

Under the assumptions (A1) — (A3), (B1) — (B6), (B9) the spectral density estimator $\{\hat{f}_{\omega}^{X}\}_{\omega\in[-\pi,\pi]}$ is consistent uniformly in the supremum norm (Rubín and Panaretos, 2020). Moreover an upper bound for the convergence rate was established. Assuming further (B7) ensures the consistency of the measurement error variance estimator $\hat{\sigma}^{2}$ given by formula (2.10).

Likewise, the estimator of the cross-spectral density $\{\hat{f}_{\omega}^{ZX}\}_{\omega\in[-\pi,\pi]}$ is consistent by the following proposition.

Proposition 2**.**

Under the conditions (A1) — (A3), (B1) — (B3), (B8), (B9), the cross-spectral density is estimated consistently:

[TABLE]

Assuming further the condition (B5), we obtain the convergence rate:

[TABLE]

The estimators of the spectral density and the cross-spectral density are essential building blocks for the estimation of the filter coefficients. In the following paragraphs we list the conditions for the consistency of the filter coefficients obtained via Tikhonov regularisation (2.19) and via truncation regularisation (2.18).

First, the following condition is required for the regression model (2.2) to be identifiable, regardless of the regularisation method used.

1. (C)

For all $\omega\in[-\pi,\pi]$ the operators $\mathscr{F}^{X}_{\omega}:\mathcal{H}\to\mathcal{H}$ satisfy $\ker\left(\mathscr{F}^{X}_{\omega}\right)=0$.

To ensure the consistency of the filter coefficients estimator by the Tikhonov method we only need to guarantee that the regularisation parameter vanishes slowly.

1. (D)

The Tikhonov regularisation parameter satisfies

[TABLE]

The following theorem establishes the consistency of the Tikhonov filter coefficient estimators.

Theorem 1**.**

Under the conditions (A1) — (A3), (B1) — (B6), (B8), (B9), (C), (D), the filter coefficient estimators (2.19) constructed by means of Tikhonov regularisation are consistent in the sense that:

[TABLE]

We now turn to the truncation estimator (2.18) of the filter coefficients, whose consistency requires more technical assumptions. We use the result by Hörmann, Kidziński and Kokoszka (2015, Theorem 1) relies on having consistent estimators of the spectral density and cross-spectral density operators with a known rate of convergence, on a condition on the eigenvalue spacing, and on an assumption that the spectral truncation parameter $K_{\omega}$ grows sufficiently slowly. In what follows, we review their conditions and adapt them to the setting when the spectral density kernels and the cross-spectral density are estimated by the kernel smoothing methods from sparse noisy observations.

Recall the eigendecomposition of the spectral frequency operator (2.4) and that its harmonic eigenvalues and harmonic eigenfunction are denoted $\{\lambda^{\omega}_{k}\}_{k=1}^{\infty}$ and $\{\varphi^{\omega}_{k}\}_{k=1}^{\infty}$ respectively. Define

[TABLE]

The following condition guarantees that the eigenspaces belonging to each of the eigenvalues $\{\lambda_{m}^{\omega}\}_{m=1}^{\infty}$ are one-dimensional, hence the eigenfunctions $\{\varphi_{m}^{\omega}\}_{m=1}^{\infty}$ can be identified (up to multiplication by a complex number with modulus 1).

1. (E1)

For all $k\geq 1$ we assume $\inf_{\omega\in[-\pi,\pi]}\Lambda_{k}^{\omega}>0$.

Furthermore we need to assume that the truncation parameter $K_{\omega}$ needs to grow sufficiently slowly.

1. (E2)

[TABLE]

where

[TABLE]

and where we further define

[TABLE]

and $\{\hat{\Lambda}_{m}^{\omega}\}$ are the empirical counterparts of $\{\Lambda_{m}^{\omega}\}$ where the estimates $\{\hat{\lambda}_{j}^{\omega}\}$ are plugged-in.

Under the above stated assumptions, the filter coefficient estimator (2.18) obtained by means of truncation regularisation is consistent.

Theorem 2**.**

Under the conditions (A1) — (A3), (B1) — (B6), (B8), (B9), (C), (E1), (E2), the filter coefficients estimates (2.18) constructed by the spectral truncation regularisation are consistent:

[TABLE]

4 Numerical experiments

4.1 Simulation setting

In this simulation study we asses the performance of the proposed methodology on the basis of two criteria: the estimation error of the filter coefficients estimator (2.19), and the prediction error of the forecasts of the response process (Section 2.4). We also compare the performance of the two regularisation techniques, and corroborate that neither dominates the other. To illustrate this, we introduce two different filter coefficient function for the lagged regression, see (4.4) and (4.5), in which one technique is expected to perform better than the other, and vice versa.

The MATLAB code and the results of the simulation are openly available on GitHub (Rubín and Panaretos, 2019).

We simulate the functional regressor series $\{X_{t}\}_{t\in\mathbb{Z}}$ as functional linear processes: either a functional autoregressive process or a functional moving average process. We define $\{E_{t}\}_{t\in\mathbb{Z}}$ to be the stochastic innovation term for these linear processes. It is assumed to be a sequence of independent identically distributed zero-mean Gaussian random variables in $\mathcal{H}=L^{2}([0,1])$ with the covariance kernel $K(x,y)$ given by

[TABLE]

We simulate realisations of the functional autoregressive process of order 1 and the functional moving average process of order 4 defined in the following two settings:

• (FAR(1))

The process $\{X_{t}\}_{t\in\mathbb{Z}}$ is a functional autoregressive process of order 1 (Bosq, 2012) defined by the iteration

[TABLE]

The operator $\mathcal{A}$ is assumed to be a Hilbert-Schmidt operator and we define its kernel as $A(x,y)=\kappa\sin(x-y),x,y\in[0,1]$, where $\kappa>0$ is chosen such that $\|\mathcal{A}\|_{L(\mathcal{H})}=0.7$ and $\|\cdot\|_{L(\mathcal{H})}$ is the operator norm in the space of linear operators on $\mathcal{H}$.

• (FMA(4))

The process $\{X_{t}\}_{t\in\mathbb{Z}}$ is considered to be the functional moving average process of order 4 defined by

[TABLE]

The operators $\mathcal{M}_{1},\dots,\mathcal{M}_{4}$ are assumed to be Hilbert-Schmidt and given by their kernels $M_{1}(x,y)=\kappa_{1}\sin(x+y),M_{2}(x,y)=\kappa_{2}\sin(1-x+y),M_{3}(x,y)=\kappa_{3}\sin(1+x-y),M_{4}(x,y)=\kappa_{4}\sin(2-x-y),$ for $x,y\in[0,1]$, respectively. The constants $\kappa_{1}>0,\dots,\kappa_{4}>0$ are chosen so that $\|\mathcal{M}_{1}\|_{L(\mathcal{H})}=0.8,\|\mathcal{M}_{2}\|_{L(\mathcal{H})}=0.6,\|\mathcal{M}_{3}\|_{L(\mathcal{H})}=0.4,\|\mathcal{M}_{4}\|_{L(\mathcal{H})}=0.2$ respectively.

The functional autoregressive process (FAR(1)), defined uniquely by the equation (4.2), and the functional moving average process (FMA(4)) are stationary and Gaussian (Bosq, 2012).

Each of the above defined functional processes is simulated with a varying time length $T\in\{300,600,900,1200\}$. The sparse observations (2.1) are generated by fixing the maximal number of observations per curve $N^{max}\in\{10,20,40,60\}$. For each curve, an integer valued random variable is drawn with uniform distribution on $\{0,\dots,N^{max}\}$ corresponding to the number of spatial locations where the $X_{t}$ is observed (with measurement error, to be defined). The measurement locations $x_{tj}$ are sampled as uniform random variables on $[0,1]$. At each $x_{tj}$ location, the measurement error is added as a realisation of a centred Gaussian random variable with variance $\sigma^{2}>0$. The variance $\sigma^{2}>0$ is chosen so the signal-to-noise ratio is $\operatorname{tr}(\mathscr{R}_{0}^{X})/\sigma^{2}=20$.

For the lagged regression model (2.2) we consider regression models where only certain filter coefficients $b_{k}$ are nonzero functions. In particular, the nonzero filter coefficients are considered to be either of the two following options

• (A)

[TABLE]

• (B)

[TABLE]

The considered kernels and filter functions are visualised on Figure 2.

We consider 3 lagged regression schemata with a varying set of nonzero filter coefficients.

• (reg1)

The filter coefficients $b_{0},b_{1}$ are set to either $\beta_{A}$ or $\beta_{B}$.

• (reg2)

The filter coefficients $b_{0},b_{3}$ are set to either $\beta_{A}$ or $\beta_{B}$.

• (reg3)

The filter coefficients $b_{0},b_{1},b_{2},b_{3},b_{4},b_{5}$ are nonzero but with decaying magnitude. They are set to either

[TABLE]

or

[TABLE]

depending on the chosen shape of the filter coefficients.

The variance of the measurement error $\{e_{t}\}_{t\in\mathbb{Z}}$ is set to be $\tau^{2}=0.001$.

For each combination of the settings, i.e. each of the 2 linear processes of $\{X_{t}\}_{t\in\mathbb{Z}}$, each of 4 length parameters $T\in\{300,600,900,1200\}$, each of 4 sampling density parameters $N^{max}\in\{10,20,40,60\}$, each of 3 regression schemata, and 2 shapes of the filter coefficients, we run 90 independent runs. The simulations have to be obviously performed in a finite dimension. We approximate the infinite dimensional dynamics of the process $\{X_{t}\}_{t\in\mathbb{Z}}$ by the B-Spline basis of dimension 21. This follows the same scheme as Rubín and Panaretos (2020), where further details on the definition of the basis as well as the approximation procedure can be found. Moreover, we consider also the regime of complete functional observations in the setting of Hörmann, Kidziński and Kokoszka (2015) in order to compare how much information is lost due to sparse sampling.

4.2 Estimation procedure

The proposed methodology requires the selection of the tuning parameters. We implemented the choice of the bandwidths $B_{R}$ and $B_{C}$ for the estimation of $\{\mathscr{F}_{\omega}^{X}\}_{\omega\in[-\pi,\pi]}$ and $\sigma^{2}$ by means of K-fold cross validation, as explained in detail in Rubín and Panaretos (2020). The Bartlett span parameter for the estimation of the spectral density is set to $L=\lfloor 2T^{1/3}\rfloor$.

In order to select the regularisation parameters for either of the two proposed regularsation methods we resort to holdout cross-validation. We split the response time series $Z_{1},\dots,Z_{T}$ into the training set $Z_{1},\dots,Z_{S}$ and the test set $Z_{S+1},\dots,Z_{T}$. The split is set to be 80:20 in favour of the training set, i.e. $S=0.8T$. The cross-spectral density $\{\mathscr{F}_{\omega}^{ZX}\}_{\omega\in[-\pi,\pi]}$ is estimated from the data in the training set $Z_{1},\dots,Z_{S}$ and the bandwidth parameter $B_{C}$ is selected by K-fold cross-validation within the training set.

The truncation estimator (2.16) is constructed by means of eigenvalue thresholding (Hörmann, Kidziński and Kokoszka, 2015). Specifically, we set $K_{\omega}(\upsilon)=\operatorname*{arg\,max}_{m\geq 1}\left\{\hat{\lambda}_{m}^{\omega}>\upsilon\right\}$ where $\upsilon>0$ is a parameter to be chosen by holdout cross-validation in the following way. Having estimated $\{\mathscr{F}_{\omega}^{X}\}_{\omega\in[-\pi,\pi]}$ and $\sigma^{2}$ from the sparsely observed regressor time series $\{X_{t}\}_{t=1}^{T}$, and the cross-spectral density $\{\mathscr{F}_{\omega}^{ZX}\}_{\omega\in[-\pi,\pi]}$ from $\{X_{t}\}_{t=1}^{T}$ and the training partition of the response $Z_{1},\dots,Z_{S}$, the spectral transfer function is estimated by the formula (2.16) using a candidate value of $K_{\omega}(\upsilon)$. The forecasts $\hat{Z}_{S+1},\dots,\hat{Z}_{T}$ are produced by the methodology outlined in Section 2.4. Comparing the forecasts with the true values of $Z_{S+1},\dots,Z_{T}$ yields a mean square forecast error on the holdout partition which we minimse with respect to $\upsilon$. The Tikhonov estimator (2.17) involves the selection of the parameter $\rho$. In the same way as for the truncation regularisation, we chose $\rho$ by holdout cross-validation based on the mean square forecast error on $Z_{S+1},\dots,Z_{T}$.

In the case of complete functional observations we again use holdout cross-validation for the selection of the eigenvalue thresholding parameter as well as the Tikhonov parameter.

4.3 Evaluation criteria

We asses the estimation error of the filter coefficients by the mean square error criterion:

[TABLE]

Next we want to asses the forecasting performance of the proposed methodology. Because the entire sample was used for fitting the model dynamics, we simulate an independent copy of the regressor time series, denoted as $\{X^{copy}_{t}\}_{t=1}^{T}$, and the response process $\{Z^{copy}_{t}\}_{t=1}^{T}$. Using the estimates of the model dynamics from the original data, we produce the predictions $\{\hat{Z}^{copy}_{t}\}_{t=1}^{T}$ from the sparsely observed $\{X^{copy}_{t}\}_{t=1}^{T}$ and compare with the true values $\{Z^{copy}_{t}\}_{t=1}^{T}$. The prediction relative mean square error is then defined

[TABLE]

Moreover, we include the prediction error of the oracle estimator that assumes that both the dynamics of the regressor time series $\{X_{t}\}_{t\in\mathbb{Z}}$ and the filter coefficients $\{\mathcal{B}_{k}\}_{k\in\mathbb{Z}}$ are known. The oracle estimator completes the steps 3 and 4 of the algorithm of Section 2.4 where the estimates $\{\hat{R}^{X}_{h}(\cdot,\cdot)\}_{h\in\mathbb{Z}}$ and $\{\hat{\mathcal{B}}^{trunc}_{k}\}_{k\in\mathbb{Z}}$ (or $\{\hat{\mathcal{B}}^{Tikh}_{k}\}_{k\in\mathbb{Z}}$) are replaced by the true values of $\{R^{X}_{h}(\cdot,\cdot)\}_{h\in\mathbb{Z}}$ and $\{\mathcal{B}_{k}\}_{k\in\mathbb{Z}}$.

4.4 Results of numerical experiments

Due to the large number of simulation settings considered, we display the results in an aggregated form. Figures 3 and 5 present the results for the filter coefficient estimation and the prediction performance as a function of the sample size parameters $T$ and $N^{max}$. The results are aggregated over both types of simulated dynamics of $\{X_{t}\}_{t}$, i.e. the functional autoregressive process (FAR(1)) and the functional moving average process (FMA(4)), and over all three regression schemes (reg1), (reg2), and (reg3). Figures 4 and 6, on the other hand, present the results for different the simulated dynamics of the process and the different considered regression schemes separately, aggregated over all time series length parameters $T\in\{300,600,900,1200\}$ and all sparse observation regimes $N^{max}\in\{10,20,40,60\}$.

An inspection of figures 3 and 5 reveals that there is no clear winner between the truncation and the Tikhonov methods. The numerical experiments with the shape (A) defined by (4.4) show that the Tikhonov method dominates the truncation method in all considered settings for the estimation of the filter coefficients. The simulations with the shape (B) defined by (4.5) yield the opposite behaviour: the truncation regularisation prevails. We attribute this dichotomy to the following reasons:

• •

The shape (B) corresponds to the leading eigenfunction of the covariance kernel (4.1). Even though the functional regression is performed in the spectral domain, the spectral transfer function $\mathscr{B}_{\omega}$ is still well aligned with the first eigenfunction of $\mathscr{F}_{\omega}^{X}$ and therefore it is enough to cut off after the first eigenvalue, thus favouring truncation regularisation. See Figure 7 to visualise the alignment in the spectral domain.

• •

The shape (A) corresponds to the fourth eigenfunction of the covariance kernel (4.1). Moreover the fourth eigenvalue is tied with the fifth and the sixth one. This structure is preserved also in the spectral domain, c.f. Figure 7. Since the Tikhonov regularisation does not discard the eigenspace corresponding to any of these eigenvalues, it achieves lower estimation error in this non-aligned case. Moreover, the Tikhonov regularisation enjoys the advantage of being stable to spectral eigenvalue ties (Hall and Horowitz, 2007; Pham and Panaretos, 2018).

• •

Note that the shape (A) is generally more difficult to estimate than shape (B). This is not surprising because the signal in the fourth eigenfunction is much weaker than in the first one.

There is no notable difference in the difficulty of estimation between the two simulated dynamics of the regressor time series $\{X_{t}\}_{t}$, i.e. the functional autoregressive process (FAR(1)) and the functional moving average process (FMA(4)). The considered regression schemes do not reveal any surprises: the longer the lagged dependence is, the more difficult estimation becomes. Therefore the (reg3) scheme produces the largest errors while (reg1) the lowest.

The prediction error of the response process $\delta^{pred}$, which is presented in Figures 5 and 6, follows the same conclusions as the estimation of the filter coefficients. The shape (A) is more challenging to predict and the Tikhonov regularisation is seen to be preferable not only for estimation, but for prediction too. The shape (B) is easier to predict using the truncation regularisation. The predictions by either of the two techniques, feature twice to thrice greater prediction error $\delta^{pred}$ than the oracle estimator, i.e. the prediction assuming the model of the data to be known and the uncertainty coming only from sparse noisy sampling regime.

5 Data analysis

In this section we illustrate the proposed methodology on measurements recorded (Tammet, 2009) at the scientific observatory located at mount Wank located in southern Germany. We remark that this analysis should be seen primarily as an example of the type of data that fall in our framework, rather than a complete data analysis, since there are presumably further important covariates that otherwise should be included. The considered period is January 1, 1977 – December 31, 1979 consisting of $T=1095$ days. In particular, we analyse the interdependence of three time series:

• •

Atmospheric electricity. The ionisation processes in the atmosphere causes the air to be conductive and the conductivity can be measured in terms of electric potential difference per distance, expressed in volts per meter ($V/m$). The atmospheric electricity is an important indicator for climate research (Tammet, 2009) and air pollution (Israelsson and Tammet, 2001). However, the atmospheric electricity can be reliably measured only under fair-weather conditions, otherwise the atmospheric ionisation processes are changed and a different quantity is recorded. The standard meteorological methodologies (Xu et al., 2013; Israelsson and Tammet, 2001) suggest to discard the observations under unfair conditions and analyse only those observations recorded under fair weather. Given these guidelines, we take into account only those hourly observations of atmospheric electricity where the wind speed was below $20\,km/h$ and the atmospheric electricity $E$ itself satisfies $0<E<250\,V/m$. The meteorological community also advocate discarding the data based on cloud coverage (Xu et al., 2013; Israelsson and Tammet, 2001), but, unfortunately, the mount Wank dataset (Tammet, 2009) does not contain cloud coverage information.

Based on the above criteria, we eventually retain an unevenly sampled scalar time series which we consequently decatenate into individual days. This technique is useful (Aue, Norinho and Hörmann, 2015; Hörmann, Kidziński and Hallin, 2015; Hörmann and Kokoszka, 2010; Hörmann, Kokoszka and Nisol, 2018) in separating the intra-day variability and the temporal dependence across days. Thus we arrive at a sparsely observed functional time series, where the latent functional data are interpreted as the “atmospheric electricity had the weather been fair”.

In what follows, we denote the fair weather electricity time series as $\{X^{(E)}_{t}\}_{t=1}^{T}$. The time series features total of 18326 measurements or 16.7 measurement per day on average.

• •

Temperature. The temperature was recorded hourly at Wank over the considered period. First, we remove the yearly periodicity in the data, then divide the time domain into individual days, and finally convert the hourly observations into functional data using B-splines. The produced fully observed functional time series is denoted as $\{X^{(\tau)}_{t}\}_{t=1}^{T}$. The time series includes 21 missing days which we treat as missing completely at random.

• •

Recorded visibility. The reported visibility was recorded hourly at a range of locations. We define the response time series $\{Z_{t}\}_{t=1}^{T}$ as the average visibility on the given day. The time series includes 42 missing values which we treat as missing completely at random.

Since the goal of our analysis is to illustrate the lagged regression methodology and compare the Tikhonov and the truncation regularisation, we split the response time series into two parts, the training component $Z_{1},\dots,Z_{822}$ and the test component $Z_{823},\dots,Z_{1095}$ consisting of roughly $75\%$ and $25\%$ of the observations respectively. We fit three models on the entire time span of $\{X^{(E)}_{t}\}_{t=1}^{T}$ and/or $\{X^{(\tau)}_{t}\}_{t=1}^{T}$ and the training set $\{Z_{t}\}_{t=1}^{822}$:

• •

Atmospheric electricity model (E). In this model we use the sparsely observed functional time series $\{X^{(E)}_{t}\}_{t=1}^{T}$ as the regressor time series for the response $\{Z_{t}\}_{t=1}^{T}$ exploiting the methodology outlined in this paper.

• •

Temperature model (T). This model handles the fully observed functional time series $\{X^{(\tau)}_{t}\}_{t=1}^{T}$ as the regressor for the response $\{Z_{t}\}_{t=1}^{T}$ using the methodology developed by Hörmann, Kidziński and Kokoszka (2015).

• •

Joint model (E+T). Finally this model includes the information from both the sparsely observed functional time series $\{X^{(E)}_{t}\}_{t=1}^{T}$ and the fully observed functional time series $\{X^{(\tau)}_{t}\}_{t=1}^{T}$ to predict the response process $\{Z_{t}\}_{t=1}^{T}$ outlined by the equation

[TABLE]

which extends the model (2.2).

The estimation and the prediction in a model that includes multiple functional time series, some under sparse and others under full observation, is explained in detail in Appendix B. Figure 9 displays a schematic visualisation of the prediction for a specific day.

We estimated the filter coefficients in all three models. Figure 8 displays the estimated filter coefficients in the joint model (E+T). The filter coefficients estimates in the marginal models (E) and (T) (not presented here) are very similar to the corresponding filters estimated in (E+T) thanks to the fact that the time series $\{X^{(E)}_{t}\}_{t=1}^{T}$ and $\{X^{(\tau)}_{t}\}_{t=1}^{T}$ are essentially uncorrelated.

A first look at Figure 8 reveals that truncation-based estimates (depicted via their Riesz-representers) feature more spikes. This, combined with a worse prediction performance commented upon later, suggests that the Tikhonov regularisation provides a better fit. Therefore we comment only on the interpretation of the filter coefficients estimated by the Tikhonov regularisation. The filter coefficient $\mathcal{B}_{0}^{(E)}$ is negative, especially in the morning hours, suggesting an obvious interpretation: high atmospheric electricity (which is linked to pollution (Israelsson and Tammet, 2001)) implies a reduction of visibility. The filter coefficients $\mathcal{B}_{k}^{(\tau)}$ corresponding to the temperature time series reveal an opposite effect. The filters $\mathcal{B}_{0}^{(\tau)}$ and $\mathcal{B}_{-1}^{(\tau)}$ are positive, therefore high temperatures on the same day and the next day predict a higher visibility today. We recall that our model is not causal, therefore the filter coefficient $\mathcal{B}_{-1}^{(\tau)}$ indeed predicts an effect backwards in time.

Table 1 presents the prediction performance of the considered models. We calculate the mean square error (MSE) on the test partition $\{Z_{t}\}_{t=823}^{1095}$ and calculate the $R^{2}$ coefficient of determination. The table reveals that the Tikhonov regularisation delivers a better prediction performance for all considered models.

6 Discussion of the considered regularisation techniques

This article presents the methodology for the functional lagged regression problem where the regressor time series is observed sparsely and with noise contamination. We have shown how to estimate the (cross)-spectral density using surface smoothers. The estimation of the spectral transfer function and consequently the filter coefficients using the estimated (cross)-spectral density are ill-posed problems and therefore require regularisation. We considered two regularisation strategies, namely spectral truncation and Tikhonov regularisation, and compared them on a simulation study and the analysis of a data set. In the following we summarise some observations on the differences, strengths, and weaknesses of the two approaches.

The simulation study presented in Section 4 illustrates that neither of the two regularisation method can dominate the other. In one of the considered scenarios the spectral transfer function $\mathscr{B}_{\omega}$ is well-aligned with the leading eigenfunction of the spectral density operator $\mathscr{F}_{\omega}^{X}$ thus being estimated better by means of truncation. In the other considered setting, $\mathscr{B}_{\omega}$ is explained by the fourth, the fifth and the sixth leading eigenvalue of $\mathscr{F}_{\omega}^{X}$, which are moreover nearly tied, resulting in a more challenging estimation task where the stability of Tikhonov regularisation to ties leads to better results than truncation.

The data analysis illustration given in Section 5 analysed the dependence of visibility on atmospheric electricity and temperature. The comparative analysis of the the two regularisation methods revealed that the estimates obtained by the Tikhonov regularisation feature better predictive performance. Moreover we found that the filter coefficients estimated by the Tikhonov regularisation were easier to interpret. At least in this application setting, the spectral transfer function $\mathscr{B}_{\omega}$ does not seem to be well-aligned with the spectral density operator $\mathscr{F}_{\omega}^{X}$.

We conclude that both of the regularisation techniques should belong to the statistician’s repertoire as neither can dominate the other. However, if we were to choose only one to broadly recommend, this would be the Tikhonov approach, as it seems more robust to “spectral misalignment” and eigenvalue ties, and its theoretical treatment requires fewer assumptions.

Appendix A Lagged regression with functional response

A.1 Functional response model and used notation

In the main body of the article we commented on the fact that the main challenge of the methodology lies with the ill-posed inversion of the spectral density operator of the regressor time series $\{\mathscr{F}^{X}_{\omega}\}_{\omega\in[-\pi,\pi]}$, and less so on whether the response is scalar or functional time series. Therefore we opted to present the model with a scalar response in the main body for of the article easier reading and postpone the treatment of the functional response to this appendix.

Now we consider the case where both the regressor time series $\{X_{t}\}\equiv\{X_{t}(\cdot)\}_{t\in\mathbb{Z}}$ and the response time series $\{Z_{t}\}\equiv\{Z_{t}(\cdot)\}_{t\in\mathbb{Z}}$ are zero-mean functional time series. For simplicity of presentation we assume that they take values in the same Hilbert space $\mathcal{H}=L^{2}([0,1])$. The lagged regression model (without intercept) becomes

[TABLE]

where the filter coefficients $\mathcal{B}_{k}$ are now assumed to be Hilbert-Schmidt operators from $\mathcal{H}$ to $\mathcal{H}$, and $\{e_{t}\}_{t\in\mathbb{Z}}$ is a sequence of independent identically distributed zero-mean random elements in $\mathcal{H}$ that is further independent of $\{X_{t}\}_{t\in\mathbb{Z}}$. We denote by $B_{k}(\cdot,\cdot)$ the integral kernel corresponding to the operator $\mathcal{B}_{k}$ for $k\in\mathbb{Z}$. The assumption (A3) is replaced by

1. S.(A3)’

[TABLE]

where $\|\cdot\|_{2}$ is the Hilbert-Schmidt norm.

Both the regressor time series $\{X_{t}\}$ and the response time series $\{Z_{t}\}$ are assumed to be observed sparsely

[TABLE]

where $x_{tj}$ and $z_{tj}$ are the measurement locations, $N_{t}^{X}$ and $N_{t}^{Z}$ are the number of measurement locations at time $t$, and the $\epsilon_{tj}^{X}$ and $\epsilon_{tj}^{Z}$ are additive measurement error contaminants. The $\{\epsilon_{tj}^{X}\}_{tj}$ and $\{\epsilon_{tj}^{Z}\}_{tj}$ are assumed to be two ensembles of independent identically distributed zero-mean real random variables with variances $\sigma^{2}_{X}>0$ and $\sigma^{2}_{Z}>0$ respectively. Moreover, the two ensembles are independent of each other and the underlying pair of functional time series $\{X_{t}\}$ and $\{Z_{t}\}$.

A.2 Nonparametric estimation of the cross-spectral density and spectral transfer function

In this section we describe how to estimate the cross-spectral density operator $\{\mathscr{F}_{\omega}^{ZX}\}_{\omega\in[-\pi,\pi]}$. For $h\in\mathbb{Z},|h|<T,t=\max(1,1+h),\dots,\min(T,T-h),\quad j=1,\dots,N^{Z}_{t+h},\quad k=1,\dots,N^{X}_{t},$ define the raw covariances as

[TABLE]

The cross-spectral density is estimated by smoothing the raw covariances using a local-linear surface smoother. For fixed frequency $\omega\in[-\pi,\pi]$, Bartlett span parameter $L\in\mathbb{N}$, and smoothing bandwidth $B_{R}>0$ we set

[TABLE]

where

[TABLE]

Once the cross-spectral density $\{\mathscr{F}_{\omega}^{ZX}\}_{\omega\in[-\pi,\pi]}$ and the spectral density $\{\mathscr{F}_{\omega}^{ZX}\}_{\omega\in[-\pi,\pi]}$ have been estimated we turn our attention to the estimation of the spectral transfer function.

The spectral function is estimated by the formulae (2.16) and (2.17) where $\hat{\mathscr{F}}_{\omega}^{ZX}\left(\hat{\varphi}_{j}^{\omega}\right)$ is now the application of the operator to the function $\hat{\varphi}_{j}^{\omega}$, yielding a function as output. Thus $\hat{\mathscr{B}}_{\omega}^{trunc}$ and $\hat{\mathscr{B}}_{\omega}^{Tikh}$ are operators. The filter coefficients, now operators, are again recovered by the formulae (2.18) and (2.19).

A.3 Forecasting the response process

We extend the setting outlined in Section 2.4 in the main body of the article. We assume that the response functional time series $\{Z_{t}\}$ is sparsely observed only up to point $S$ for $S<T$ and we wish to forecast the values $Z_{S+1},\dots,Z_{T}$.

The forecasting algorithm follows in fact the same steps as presented in Section 2.4. The only difference is that the formulae in step 4 for the forecast construction $\hat{\Pi}(Z_{s}|\mathbb{Y})=\sum_{k=-M}^{M}\hat{\mathcal{B}}_{k}^{trunc}\hat{\Pi}({X}_{s-k}|\mathbb{Y})$ (or alternatively $\hat{\Pi}(Z_{s}|\mathbb{Y})=\sum_{k=-M}^{M}\hat{\mathcal{B}}_{k}^{Tikh}\hat{\Pi}({X}_{s-k}|\mathbb{Y})$) output elements of the function space $\mathcal{H}$ because $\hat{\mathcal{B}}_{k}^{trunc}$ and $\hat{\mathcal{B}}_{k}^{Tikh}$ are now operators.

Appendix B Lagged regression with multiple inputs

B.1 Joint regression model and its spectral analysis

It may very well happen that we want to analyse the relationship between multiple time series, including functional time series observed either sparsely or fully, multivariate time series, and scalar time series. To keep the notation simple we focus our presentation on the setting of a scalar response and two functional time series regressors, one of which is observed sparsely and the other fully. This setting corresponds to the model (E+T) of the data analysis in Section 5. The generalisation to a higher number of functional time series, observed sparsely or fully, is clear, treating multivariate time series as time series in a finite-dimensional Hilbert space; and, the further incorporation of multivariate or scalar time series is uncomplicated because the regularisation of the spectral density inversion is no longer required.

Consider the model with two functional time series regressors $\{X_{t}^{(1)}\}\equiv\{X_{t}^{(1)}(\cdot)\}_{t\in\mathbb{Z}}$ and $\{X_{t}^{(2)}\}\equiv\{X_{t}^{(2)}(\cdot)\}_{t\in\mathbb{Z}}$. To simplify the notation we assume the both functional time series to take values in the same Hilbert space $\mathcal{H}$. The lagged regression with a scalar response $\{Z_{t}\}\equiv\{Z_{t}\}_{t\in\mathbb{Z}}$ becomes:

[TABLE]

where $a\in\mathbb{R}$ is the intercept, $\{\mathcal{B}_{k}^{(1)}\}_{k\in\mathbb{Z}}$ and $\{\mathcal{B}_{k}^{(1)}\}_{k\in\mathbb{Z}}$ are two sequences of linear mappings from $\mathcal{H}$ to $\mathbb{R}$, and $\{e_{t}\}$ is a sequence of zero mean independent identically distributed real random variables. Denote $\{b_{k}^{(1)}\}_{k\in\mathbb{Z}}$ and $\{b_{k}^{(1)}\}_{k\in\mathbb{Z}}$ the filter functions corresponding to the functionals $\{\mathcal{B}_{k}^{(1)}\}$ and $\{\mathcal{B}_{k}^{(1)}\}$ by the Riesz representation theorem. We assume that the first functional time series, $\{X_{t-k}^{(1)}\}$, is sparsely observed, that is we have access to only the observations generated by

[TABLE]

where $N_{t}$ is a number of observation locations $\{x_{tj}\}$ at time $t=1,\dots,T$. The second functional time series is fully observed, therefore we have access to $X_{1}^{(2)},\dots,X_{T}^{(2)}$. The mean functions of $\{X_{t}^{(1)}\}$ and $\{X_{t}^{(2)}\}$ are denoted as $\mu^{(1)}(\cdot)$ and $\mu^{(2)}(\cdot)$ respectively.

Denote $R_{h}^{\ast}$ and $\mathscr{R}_{h}^{\ast}$ the (auto)covariance kernel and (auto)covariance operator where $\ast$ is substituted by a single functional time series or a pair theoreof, as the case may be. Likewise we denote the (cross)-spectral density kernel and the (cross)-spectral density operator as $f_{\omega}^{\ast}$ and $\mathscr{F}_{\omega}^{\ast}$.

Assuming

[TABLE]

and $\sum_{k\in\mathbb{Z}}\|b_{k}^{(1)}\|_{\infty}<\infty$, $\sum_{k\in\mathbb{Z}}\|b_{k}^{(2)}\|_{\infty}<\infty$, $\sum_{k\in\mathbb{Z}}\|\mathcal{B}_{k}^{(1)}\|_{\mathcal{H}}<\infty$, $\sum_{k\in\mathbb{Z}}\|\mathcal{B}_{k}^{(2)}\|_{\mathcal{H}}<\infty$ implies that the (cross)-spectral density kernels and operators $\{f_{\omega}^{X^{(1)}}\}_{\omega}$, $\{f_{\omega}^{X^{(2)}}\}_{\omega}$, $\{f_{\omega}^{X^{(1)}X^{(2)}}\}_{\omega}$, $\{f_{\omega}^{ZX^{(1)}}\}_{\omega}$, $\{f_{\omega}^{ZX^{(2)}}\}_{\omega}$, $\{\mathscr{F}_{\omega}^{X^{(1)}}\}_{\omega}$, $\{\mathscr{F}_{\omega}^{X^{(2)}}\}_{\omega}$, $\{\mathscr{F}_{\omega}^{X^{(1)}X^{(2)}}\}_{\omega}$, $\{\mathscr{F}_{\omega}^{ZX^{(1)}}\}_{\omega}$, $\{\mathscr{F}_{\omega}^{ZX^{(2)}}\}_{\omega}$ are well defined.

The joint frequency response operator

[TABLE]

satisfies the relation

[TABLE]

The filter coefficients can be recovered by the following formula, for $k\in\mathbb{Z}$,

[TABLE]

B.2 Nonparametric estimation

The mean function of the sparsely observed $\{X_{t}^{(1)}\}$ can be estimated by the local linear smoother (Rubín and Panaretos, 2020). The mean function of the fully observed $\{X_{t}^{(2)}\}$ is estimated by the classical empirical pointwise mean. We denote these estimates $\hat{\mu}^{(1)}(\cdot)$ and $\hat{\mu}^{(2)}(\cdot)$ respectively.

The estimation of the spectral density $\{\mathscr{F}_{\omega}^{X^{(1)}}\}_{\omega\in[-\pi,\pi]}$ and the cross-spectral density $\{\mathscr{F}_{\omega}^{ZX^{(1)}}\}_{\omega\in[-\pi,\pi]}$ from the sparse observations is explained in the main body of the article. For the fully observed functional case we use functional version of Bartlett’s estimator (Hörmann, Kidziński and Hallin, 2015; Hörmann, Kidziński and Kokoszka, 2015)

[TABLE]

where $\hat{\mathscr{R}}_{h}^{X^{(2)}}$ and $\hat{\mathscr{R}}_{h}^{ZX^{(2)}}$ are the classical empirical lag-$h$ autocovariance operators of $\{X^{(2)}_{t}\}$ and the lag-$h$ covariance operator of $\{Z_{t}\}$ and $\{X^{(2)}_{t}\}$ respectively.

It remains to estimate the cross-spectral density $\{\mathscr{F}_{\omega}^{X^{(1)}X^{(2)}}\}_{\omega\in[-\pi,\pi]}$. Aiming to estimate $f_{\omega}^{X^{(1)}X^{(2)}}(x,y)$ for $x,y\in[0,1],\omega\in[-\pi,\pi]$, we define the raw covariances

[TABLE]

and for fixed $(x,y)\in[0,1]^{2}$ we estimate the cross-spectral density as

[TABLE]

where $\hat{d}_{0}\in\mathbb{C}$ is the minimiser of the following weighted sum of squares

[TABLE]

Once we have estimated all the above (cross-)spectral densities, we wish to recover the filter coefficients from the equation (B.1). The inversion of the joint spectral density of $\{X^{(1)}_{t}\}$ and $\{X^{(1)}_{t}\}$ is ill-conditioned due to the reasons explained in Section 2.3 of the main body of the paper where we overcame this issue by two regularisation techniques: spectral truncation and Tikhonov regularisation. Here we do the same while allowing each of the regressor time series to have a different degree of regularisation. This is important because, generally speaking, the estimation from sparse data will require more regularisation.

Denote the spectral decompositions of the (cross)-spectral density operators by

[TABLE]

where $\{\hat{\lambda}_{i}^{\omega}\}_{i=1}^{\infty}$ and $\{\hat{\eta}_{i}^{\omega}\}_{i}^{\infty}$ are the harmonic eigenvalues of $\mathscr{F}_{\omega}^{X^{(1)}}$ and $\mathscr{F}_{\omega}^{X^{(2)}}$ respectively, for given $\omega\in[-\pi,\pi]$. The sequences of functions $\{\hat{\varphi}^{\omega}\}_{i}^{\infty}$ and $\{\hat{\psi}^{\omega}\}_{i}^{\infty}$ are orthonormal bases of $\mathcal{H}$ and are called the harmonic eigenfunctions of $\mathscr{F}_{\omega}^{X^{(1)}}$ and $\mathscr{F}_{\omega}^{X^{(2)}}$ respectively. The complex numbers $\{\hat{\gamma}_{ij}^{\omega}\}_{i,j=1}^{\infty}$ are the basis coefficients of $\mathscr{\hat{F}}_{\omega}^{X^{(1)}X^{(2)}}$ with respect to $\{\hat{\varphi}_{i}^{\omega}\otimes\hat{\psi}_{j}^{\omega}\}_{i,j=1}^{\infty}$, a basis of the space of Hilbert-Schmidt operators on $\mathcal{H}$.

Truncation regularisation. Select two (possibly different) truncation parameters $K_{\omega}^{(1)}\in\mathbb{N}$ and $K_{\omega}^{(2)}\in\mathbb{N}$ that may depend on $\omega\in[-\pi,\pi]$. The idea of the truncation regularisation is to replace the empirical joint spectral density operator on the right hand side of the empirical version of (B.1) by

[TABLE]

Alternatively, the operator in the form (B.2) can be expressed with respect to the reduced basis of $\mathcal{H}\times\mathcal{H}$ composed of $\{[\hat{\varphi}_{1}^{\omega},\textbf{0}],\dots,[\hat{\varphi}_{K_{\omega}^{(1)}}^{\omega},\textbf{0}],[\textbf{0},\hat{\psi}_{1}^{\omega}],\dots,[\textbf{0},\hat{\psi}_{K_{\omega}^{(2)}}^{\omega}]\}$ where 0 is the zero element of $\mathcal{H}$. Denote $M\in\mathbb{C}^{({K_{\omega}^{(1)}+K_{\omega}^{(2)}})\times({K_{\omega}^{(1)}+K_{\omega}^{(2)}})}$ the complex matrix given by the inverse of the following matrix, assuming it is invertible,

[TABLE]

Denote the elements of $M$ according to the blocks using the following scheme

[TABLE]

Then, the frequency response operators can be recovered by

[TABLE] 2. 2.

Tikhonov regularisation. We consider again two possibly different regularisation parameters $\rho^{(1)}>0$ and $\rho^{(2)}>0$.

The empirical joint spectral density operator on the right hand side of the empirical version of (B.1) is replaced by

[TABLE]

and the joint frequency response operator is estimated as

[TABLE]

Once the estimates of the frequency response operators have been constructed, the filter coefficients are estimated by integrating back into the temporal domain

[TABLE]

B.3 Forecasting the response process

In this section we extend the forecasting algorithm introduced in Section 2.4 of the main body of the article. As in that Section, we assume here that the response time series is observed only at times $1,\dots,S$ where $S<T$ and we wish to predict $Z_{S+1},\dots,Z_{T}$.

From the measurements $Y_{tj},t=1,\dots,T,j=1,\dots,N_{t}$, of sparsely observed $\{X_{t}^{(1)}\}$, the fully observed functional observations $X_{1}^{(1)},\dots,X_{T}^{(1)}$, and the response time series $Z_{1},\dots,Z_{S}$, construct the estimates of the model: the mean functions $\hat{\mu}^{(1)},\hat{\mu}^{(2)}$, the spectral densities $\{\hat{\mathscr{F}}^{X^{(1)}}_{\omega}\}_{\omega\in[-\pi,\pi]},\{\hat{\mathscr{F}}^{X^{(2)}}_{\omega}\}_{\omega\in[-\pi,\pi]}$ and the cross-spectral densities $\{\hat{\mathscr{F}}^{X^{(1)}X^{(2)}}_{\omega}\}_{\omega\in[-\pi,\pi]},\{\hat{\mathscr{F}}^{ZX^{(1)}}_{\omega}\}_{\omega\in[-\pi,\pi]},\{\hat{\mathscr{F}}^{ZX^{(2)}}_{\omega}\}_{\omega\in[-\pi,\pi]}.$ Using either the truncation regularisation or the Tikhonov regularisation, estimate the spectral transfer functions $\{\hat{\mathscr{B}}_{\omega}^{(1)}\}_{\omega\in[-\pi,\pi]}$ and $\{\hat{\mathscr{B}}_{\omega}^{(2)}\}_{\omega\in[-\pi,\pi]}$ and the filter coefficients $\{\hat{\mathcal{B}}_{k}^{(1)}\}_{k\in\mathbb{Z}}$ and $\{\hat{\mathcal{B}}_{k}^{(2)}\}_{k\in\mathbb{Z}}$. 2. 2.

Denoting the ensemble of the measurements $Y_{tj},t=1,\dots,T,j=1,\dots,N_{t}$ as $\mathbb{Y}^{(1)}$, predict the latent values of $X_{-M+1}^{(1)},\dots,X_{T+M}^{(1)}$ using the method recalled in Section 2.4 of the main body of the article. Denote these predictions as $\hat{\Pi}(X_{-M+1}^{(1)}|\mathbb{Y}^{(1)}),\dots,\hat{\Pi}(X_{T+M}^{(1)}|\mathbb{Y}^{(1)})$. The constant $M$ is determined in such a way that the filter coefficients $\hat{\mathcal{B}}_{k}^{(1)}$ and $\hat{\mathcal{B}}_{k}^{(2)}$ are negligible for $|k|>M$. 3. 3.

Pad the fully observed functional time series $\{X_{t}^{(2)}\}$ by the mean function $\hat{\mu}^{(2)}$, i.e. set $X_{t}=\hat{\mu}^{(2)}$ for $t<1$ or $t>T$. This approach was suggested in functional lagged prediction for fully functional regressor (Hörmann, Kidziński and Kokoszka, 2015) where the zero-mean functional regressor time series is padded by zeros. 4. 4.

For each $s=S+1,\dots,T$, construct the forecast

[TABLE]

where $\bar{Z}$ is the sample mean of the response process.

Appendix C Proofs of formal statements

C.1 Proof of Proposition 1

Proof.

The proof follows directly from the formula for the best linear unbiased predictors.

[TABLE]

∎

C.2 Auxiliary Lemmas

Lemma 1**.**

Assuming (A1) — (A3), the time series $\{Z_{t}\}_{t\in\mathbb{Z}}$ defined by (2.2) is stationary. The cross-covariance function between $\{Z_{t}\}_{t\in\mathbb{Z}}$ and $\{X_{t}(\cdot)\}_{t\in\mathbb{Z}}$ satisfies

[TABLE]

where the series converge in the vector norm and the supremum norm respectively. Furthermore

[TABLE]

[TABLE]

Assuming further (B3), the function $R^{ZX}_{h}(\cdot)$ is twice continuously differentiable in $y$ and

[TABLE]

is uniformly bounded in $h\in\mathbb{Z}$.

Proof.

The claims (C.1) and (C.3) were proven by (Hörmann, Kidziński and Kokoszka, 2015, Lemma 2).

For the claim (C.2) first consider a functional $\mathcal{A}:\mathcal{H}\to\mathbb{R}$ represented by $a(\cdot)\in\mathcal{H}$, i.e.

[TABLE]

and a Hilbert-Schmidt operator $\mathcal{T}:\mathcal{H}\to\mathcal{H}$ with integral kernel $T(\cdot,\cdot)$. Then its composition $\mathcal{A}\mathcal{T}$ is a functional represented by a function $y\mapsto\int a(x)T(x,y)\operatorname{d}\!{x}$ because

[TABLE]

for all $f\in\mathcal{H}$.

The convergence of the sum on the right-hand side of (C.2) and the claim (C.4) are derived analogously to (Hörmann, Kidziński and Kokoszka, 2015, Lemma 2) where the Hilbert-Schmidt norm is replaced by the supremum norm.

The second derivative of $R^{ZX}_{h}(y)$ is given by

[TABLE]

The exchange of the integration-summation and the derivative is justified by the fact that the following bound is integrable in $x$ and absolutely summable in $k$

[TABLE]

where $C$ is the uniform bound from the assumption (B3) with $(\alpha_{1},\alpha_{2})=(0,2)$. ∎

C.3 Proof of Proposition 2

The minimiser of the optimisation problem (2.15), and hence the estimator (2.14), can be expressed explicitly:

[TABLE]

where

[TABLE]

for $r\in\{0,1\}$. The above quantities are defined as functions of $x\in[0,1]$ and all of the operations are understood pointwise, including the division operation.

The derivation of the formula (C.6) is a simple modification of the standard argument in the local-polynomial regression literature (Fan and Gijbels, 1996, Section 3.1) or (Rubín and Panaretos, 2020, Section B.2).

Lemma 2**.**

Under the conditions (B1), (B2), and (B8)

[TABLE]

where the constant $U$ is uniform in $r=0,1,2$, $T\in\mathbb{N}$, $|h|<T$, and $B_{C}$, and where

[TABLE]

Furthermore,

[TABLE]

*uniformly in $x\in[0,1]$. *

Proof.

We have the usual bias-variance decomposition

[TABLE]

The bias term, i.e. the second term on the right-hand side of (C.10), can be Taylor expanded to second order and the formulae in display (C.8) follow. Moreover

[TABLE]

uniformly in $x\in[0,1]$, $h$, and $T$.

To bound the variance term, i.e. the first term on the right-hand side of (C.10), we make use of the Fourier transform. Denote the inverse Fourier transform of the function $u\mapsto K(u)u^{r}$ by $\zeta_{r}(t)=\int e^{-\operatorname{i}ut}K(u)u^{r}du$. Hence

[TABLE]

where

[TABLE]

Thanks to the above introduced notation we may separate the stochastic terms and the dependence on $x$:

[TABLE]

We firstly focus on $\mathbb{E}\left|\phi_{r}(v)-\mathbb{E}\phi_{r}(v)\right|$. By Jensen’s inequality,

[TABLE]

hence it suffices to bound the variance. By the independence of $\{N_{t}\}$ and $\{x_{tj}\}$

[TABLE]

Now we treat the integral on the right hand side of (C.13)

[TABLE]

for a constant $c$. Combining (C.13), (C.15), and (C.16) yields

[TABLE]

Combining the bounds (C.11) and (C.17) yields the existence of a constant $U$ as in the claim (C.7).

We now turn our attention to the claim (C.9). Since $L=o(T)$ we may assume $L<T/2$ and obtain $1/(T-|h|)\leq 2/T$ for $|h|<L$. Noting that $L^{-1}\sum_{h=-L}^{L}W_{h}=1$ we calculate

[TABLE]

Taking the supremum norm and the expectation we bound the first term by $\sqrt{2}U/(\sqrt{T}B_{C})$. The second term is bounded by $LM_{[S_{r}]}/T$ which is a faster rate than the one above. Hence the claim (C.9). ∎

Lemma 3**.**

Under the conditions (A1) — (A3), (B1) — (B3), and (B8)

[TABLE]

where the constant $\tilde{U}$ is uniform in $r=0,1$, $T$, $|h|<T$, and

[TABLE]

Proof.

The standard bias-variance decomposition of the left-hand side of (C.18) yields

[TABLE]

The bias therm, i.e. the second term on the right-hand side of (C.20), is treated by first conditioning on $x_{tj}$ and then expanding in a Taylor series of order 2. The bias term is of order $O(B_{C}^{2})$ uniformly in $|h|<T$. The Taylor expansion is justified by (C.5) in Lemma 1. The formulae (C.19) follow.

Now, using similar steps as in (C.12), we obtain:

[TABLE]

where

[TABLE]

Replicating the steps in (C.13) we arrive at the inequality

[TABLE]

The integral on the right-hand side of (C.21) is treated in (C.16). Using (C.14), it remains to bound the variance of $\varphi_{r}(v)$. First, remark that for an arbitrary stationary time-series $\{W_{t}\}$ with a summable autocovariance function $\{R^{W}_{h}\}_{h\in\mathbb{Z}}$, one has the bound:

[TABLE]

Define $W_{t}=\sum_{j=1}^{N_{t}}G_{h,t}^{ZX}(x_{tj})e^{\operatorname{i}vx_{tj}}$. This sequence of complex-valued random variables is clearly a stationary (scalar) time series. By conditioning on $N_{t}$ and $x_{tj}$, and applying the law of total covariance, we can bound the autocovariance function of $\{W_{t}\}$ by $\left|R_{h}^{W}\right|\leq\mathbb{E}\left[N\right]\sup_{x\in[0,1]}\left|R_{h}^{ZX}(x)\right|$ which is summable in $h$ due to the claim (C.4) of Lemma 1. Hence

[TABLE]

for a constant $K$ uniform in $h$ and $T$. Combining (C.20), (C.21), and (C.22) yields the claim (C.18). ∎

Lemma 4**.**

For $r=0,1$:

under the conditions (A1) — (A3), (B1) — (B3), (B8), (B9)

[TABLE] 2. 2.

if (B5) is further assumed, then

[TABLE]

where all convergences are uniform in $\omega\in[-\pi,\pi]$ and $x\in[0,1]$, and where

[TABLE]

Proof.

We start with the observation that the sum $\sum_{h\in\mathbb{Z}}M_{[Q_{r}^{(h)}]}$ converges absolutely in the supremum norm and $\sum_{h\in\mathbb{Z}}M_{[Q_{r}^{(h)}]}e^{-\operatorname{i}h\omega}=M_{[Q_{r}^{\omega}]}$. Moreover,

[TABLE]

The first term on the right hand side is bounded in the supremum by

[TABLE]

The last three terms on the right-hand side of (C.26) converge uniformly to zero due to the summability of $\sum_{h\in\mathbb{Z}}M_{[Q_{r}^{(h)}]}$ in the supremum norm and by Kroncker’s lemma. Hence the claim (C.23).

If (B5) is further assumed then $\sum_{h\in\mathbb{Z}}|h|\left|M_{[Q_{r}^{(h)}]}\right|<\infty$ and therefore the last three terms on the right-hand side of (C.26) are of order $O(1/L)$, hence (C.24). ∎

Proof of Proposition 2.

(C.6) The proof follows directly from formula (C.6) and Lemmas 2, 3, 4:

[TABLE]

uniformly in $\omega\in[-\pi,\pi]$ and $x\in[0,1]$. Hence

[TABLE]

uniformly in $\omega\in[-\pi,\pi]$ and $x\in[0,1]$. If the condition (B5) is further assumed, we replace $o_{p}(1)$ by $O_{p}(L/(\sqrt{T}B_{C}))$. ∎

C.4 Proof of Theorem 1

Thanks to the fact that

[TABLE]

the proof reduces to establishing the convergence rate of the frequency response operator $\mathscr{B}_{\omega}$.

For the Tikhonov regularisation parameter $\rho=\rho(T)>0$ define

[TABLE]

Further, split the desired difference into three terms:

[TABLE]

We bound each of the terms $\mathscr{S}_{1}$, $\mathscr{S}_{2}$, and $\mathscr{S}_{3}$ and show the convergence of the bound to zero uniformly in $\omega\in[-\pi,\pi]$. We start with bounding $\mathscr{S}_{1}$.

[TABLE]

Now, bounding $\mathscr{S}_{2}$:

[TABLE]

The right-hand sides of (C.28) tend to zero uniformly in $\omega\in[-\pi,\pi]$ as $T\to\infty$ thanks to assumption (D) and Proposition 2. The right-hand side of (C.29) tend to zero uniformly in $\omega\in[-\pi,\pi]$ as $T\to\infty$ thanks to assumption (D) and Rubín and Panaretos (2020, Theorem 2).

It remains to handle the deterministic term $\mathscr{S}_{3}$. Since the spectral density operator $\mathscr{F}^{X}_{\omega}$ is self-adjoint and trace-class, it admits the series decomposition

[TABLE]

where $\{\lambda_{j}^{\omega}\}_{j\in\mathbb{N}}$ and $\{\varphi_{j}^{\omega}\}_{j\in\mathbb{N}}$ are the harmonic (a.k.a. dynamic) eigenvalues and the harmonic (dynamic) eigenfunctions at frequency $\omega\in[-\pi,\pi]$ (Panaretos and Tavakoli, 2013b; Hörmann, Kidziński and Hallin, 2015). From the relation (2.8) we obtain the following series expansions

[TABLE]

Combining the above expansions yields

[TABLE]

Since $\sum_{j=1}^{\infty}\left\|\mathscr{B}_{\omega}\varphi_{j}^{\omega}\right\|^{2}=\sum_{j=1}^{\infty}\|\mathscr{B}_{\omega}\|^{2}<\infty$ and $\rho/(\lambda_{j}^{\omega}+\rho)\to 0$ as $\rho\searrow 0$ for each $j\in\mathbb{N}$, the right-hand side of (C.30) tends to zero, completing the proof.

C.5 Proof of Theorem 2

The proof of this theorem is an application of the result by Hörmann, Kidziński and Kokoszka (2015, Theorem 1). Their theorem requires that the spectral density operators $\{\mathscr{F}^{X}_{\omega}\}_{\omega\in[-\pi,\pi]}$ are estimated with a certain rate, say $(\psi^{X}_{T})_{T=1}^{\infty}$, and the cross-density operators $\{\mathscr{F}^{ZX}_{\omega}\}_{\omega\in[-\pi,\pi]}$ with another rate, say $(\psi_{T}^{ZX})_{T=1}^{\infty}$. For $T\in\mathbb{N}$ we put $\psi^{X}_{T}=LT^{-1/2}B_{R}^{-2}$ and $\psi_{T}^{ZX}=LT^{-1/2}B_{C}^{-1}$ by Rubín and Panaretos (2020, Theorem 2) and Proposition 2 respectively. Note that the convergence rates for the spectral density kernels $\{f_{\omega}^{X}\}_{\omega\in[-\pi,\pi]}$ and the cross-spectral density kernels $\{f_{\omega}^{ZX}\}_{\omega\in[-\pi,\pi]}$ are in the supremum norm which is stronger than the operator norm of the assumptions of Hörmann, Kidziński and Kokoszka (2015, Theorem 1).

We can now replicate all the steps of the proof of Hörmann, Kidziński and Kokoszka (2015, Theorem 1) and see that they only require the above stated rates and the assumptions (C), (E1), and (E2).

It is worth noting that the results by Hörmann, Kidziński and Kokoszka (2015, Theorem 1) are derived under the assumption of $L^{p}$-$m$-approximability as opposed to cumulant mixing conditions considered in this paper. Nevertheless, they use the assumption of the $L^{p}$-$m$-approximability only to prove the rate of convergence of the spectral density and cross-spectral density estimators in (Hörmann, Kidziński and Kokoszka, 2015, Lemma 1). Once these rates are established, the proof of (Hörmann, Kidziński and Kokoszka, 2015, Theorem 1) does not use the $L^{p}$-$m$-approximability. Therefore our proof of Theorem 2 is indeed its simple adaptation which takes as inputs the convergence rates $\psi^{X}_{T}=LT^{-1/2}B_{R}^{-2}$ and $\psi_{T}^{ZX}=LT^{-1/2}B_{C}^{-1}$ proved under the cumulant mixing conditions.

Data availability statement

The MATLAB code and the results that support the findings of this study are openly available on GitHub at DOI 10.5281/zenodo.3190740, reference Rubín and Panaretos (2019).