Simultaneous Inference for Time Series Functional Linear Regression

In order to facilitate the formulation of roughness penalization estimation, we first prepare some notations. Throughout this paper, we assume that $X_{ij}(t)$ for $j=1,...,p,~{}i=1,...,n$ is continuous on $[0,1]$ a.s. and hence admits the following expansion $X_{ij}(t)=\sum_{k=1}^{\infty}\widetilde{x}_{ij,k}\alpha_{k}(t),$ where $\{\alpha_{k}(t)\}_{k=1}^{\infty}$ is a set of pre-selected orthonormal basis functions of $H$. From Theorem 1 of [49], $X_{ij}(t)$ has the standard Karhunen-Loève type expansion as follows

$$X_{ij}(t)=\sum_{k=1}^{\infty}f_{jk}x_{ij,k}\alpha_{k}(t),$$

(2)

where $f_{jk}={\rm Std}(\widetilde{x}_{ij,k})$ with ${\rm Std}$ denoting standard deviation and $x_{ij,k}=\widetilde{x}_{ij,k}/f_{jk}$ if $f_{jk}\neq 0$. Set $x_{ij,k}=0$ if $f_{jk}=0$. Notice that $f_{jk}$ captures the decay speed of $\widetilde{x}_{ij,k}$ as $k$ increases and the random coefficient $\{x_{ij,k}\}_{k=1}^{\infty}$ remains at the same magnitude with variance $1$ as $k$ increases if $f_{jk}\neq 0$. Similarly, write $\beta_{j}(t)=\sum_{k=1}^{\infty}\beta_{jk}\alpha_{k}(t).$ The following assumption restricts the decay speed of the basis expansion coefficients of $\{\bm{X}_{i}(t)\}$ and $\{\bm{\beta}(t)\}$.

Let $d\geq 0$ be an integer. It is well-known that for a general $\mathcal{C}^{d}([0,1])$ function the fastest decay rate for its $k$-th basis expansion coefficient is $O(k^{-d-1})$ for a wide class of basis functions ([11]). For instance, the Fourier basis (for periodic functions), the weighted Chebyshev polynomials ([53]) and the orthogonal wavelets with degree $m\geq d$ ([39]) admit the latter decay rate under some extra mild assumptions on the behavior of the function’s $d$-th derivative. Hence Assumption 2.1 essentially requires that the basis expansion coefficients of $\beta_{j}(t)$ and $X_{ij}(t)$ decay at the fastest rate. On the other hand, we remark that the basis expansion coefficients may decay at slower speeds for some orthonormal bases. An example is the Legendre polynomials where the coefficients decay at an $O(k^{-d-1/2})$ speed ([55]). For basis functions whose coefficients decay at slower rates, following the proofs of this paper it is obvious to see that the multiplier bootstrap method for the JSCB construction is still asymptotically valid under the corresponding restrictions on the tuning parameters. However, in this case the estimates of $\bm{\beta}(t)$ will converge at a slower speed and the bootstrap approximation may be less accurate. For the sake of brevity we shall stick to the fastest decay Assumption 2.1 for our theoretical investigations throughout this paper.

The roughness penalization approach to the FLM (1) solves the following penalized least squares problem

$$\widetilde{\bm{\beta}}(t)=\mathop{\arg\min}_{\bm{\beta}(t)}\left\{\frac{1}{n}\sum_{i=1}^{n}\left[Y_{i}-\sum_{j=1}^{p}\int_{0}^{1}\beta_{j}(t)X_{ij}(t)\mathrm{d}t\right]^{2}+\lambda\sum_{j=1}^{p}\int_{0}^{1}[\beta^{\prime\prime}_{j}(t)]^{2}\mathrm{d}t\right\},$$

(3)

see [50] and references therein. On the other hand, in functional data analysis, researchers usually work on the optimization problem on a finite-dimensional subspace ([34]), which makes the procedure easily implementable. Following the method proposed by [45], we truncate the coefficient functions to finite (but diverging) dimensional spans of a priori set of basis functions $\{\alpha_{k}(t)\}_{k=1}^{\infty}$ while involving a smoothness penalty. Specifically, assume the truncation number of $\beta_{j}(t)$ equals $c_{j}$ ($c_{j}\rightarrow\infty$). Then the right hand side of (3) can be approximated by

$$\frac{1}{n}\sum_{i=1}^{n}\left[Y_{i}-\sum_{j=1}^{p}\int_{0}^{1}\beta_{j,c_{j}}(t)X_{ij}(t)\mathrm{d}t\right]^{2}+\lambda\sum_{j=1}^{p}\int_{0}^{1}[\beta^{\prime\prime}_{j,c_{j}}(t)]^{2}\mathrm{d}t,$$

where $\beta_{j,c_{j}}(t)=\sum_{k=1}^{c_{j}}\beta_{jk}\alpha_{k}(t)$. Simplifying the above expression, the estimation can be achieved by minimizing the following penalized least squares criterion function

		$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left[Y_{i}-\sum_{j=1}^{p}\sum_{k=1}^{c_{j}}\beta_{jk}\widetilde{x}_{ij,k}\right]^{2}+\lambda\sum_{j=1}^{p}\sum_{k,l=1}^{c_{j}}\beta_{jk}\beta_{jl}\widetilde{R}_{j}(k,l)$
	$\displaystyle=$	$\displaystyle\frac{1}{n}[\bm{Y}-\bm{X}_{c}\bm{\theta}_{c}]^{\top}[\bm{Y}-\bm{X}_{c}\bm{\theta}_{c}]+\bm{\theta}_{c}^{\top}\bm{R}(\lambda)\bm{\theta}_{c},$		(4)

where $\lambda$ ($\lambda\rightarrow 0$) is a common smoothing parameter that measures the rate of exchange between fit to the data and smoothness of the estimator, as measured by the residual sum of squares in the first term, and variability of the functions $\beta_{j,c_{j}}(t)$ in the second term. Here we impose a roughness penalty associated with the same smoothing parameter $\lambda$ for every $j$, which is commonly used in practice, see [18].

In (4), $c=c(n)=\sum_{j=1}^{p}c_{j}$ and $\bm{\theta}_{c}$ is a $c$-dimensional block vector where $\bm{\theta}_{j}^{c}=(\theta_{j1},...,\theta_{jc_{j}})^{\top}$ is the $j$-th block with $\theta_{jk}=f_{jk}\beta_{jk}$ for $k=1,...,c_{j}$. Furthermore, $\bm{X}_{c}$ is the $n\times c$ block design matrix, $\bm{R}(\lambda)$ is a $c\times c$ block diagonal matrix with $\lambda\bm{R}_{j}$ as its diagonals and $\bm{R}_{j}$ is a $c_{j}\times c_{j}$ matrix with elements $R_{j}(k,l)=\int_{0}^{1}\alpha_{k}^{\prime\prime}(t)\alpha_{l}^{\prime\prime}(t)\mathrm{d}t/(f_{jk}f_{jl})$. Then the penalized least squares estimator turns out to be

$$\widetilde{\bm{\theta}}_{c}=\left[\frac{\bm{X}_{c}^{\top}\bm{X}_{c}}{n}+\bm{R}(\lambda)\right]^{-1}\frac{\bm{X}_{c}^{\top}\bm{Y}}{n}.$$

Consequently, it is easy to find the matrix representation of the roughness penalized estimator $\widetilde{\bm{\beta}}(t)=\bm{C}_{f}(t)\widetilde{\bm{\theta}}_{c},$ where $\bm{C}_{f}(t)=(\bm{C}_{1}^{\top}(t),...,\bm{C}_{p}^{\top}(t))$ is a $p\times c$ matrix and $\bm{C}_{j}(t)$ is a $c$-dimensional block vector with $(\alpha_{1}(t)/f_{j1},...,\alpha_{c_{j}}(t)/f_{jc_{j}})^{\top}$ as its $j$-th block and other elements being 0.

JSCB construction of $\bm{\beta}(t)$ boils down to evaluating the distributional behavior of the weighted maximum deviation

$$\Xi_{n,\bm{g}_{n}}:=\sqrt{n}\sup_{t\in[0,1]}|\widetilde{\bm{\beta}}(t)-\bm{\beta}(t)|_{\bm{g}_{n}(t)},$$

where $|\bm{V}(t)|_{\bm{g}(t)}=\max_{1\leq j\leq p}|V_{j}(t)/g_{j}(t)|$ for any function $\bm{V}(t)=(V_{1}(t),\\ ...,V_{p}(t))^{\top}$ and weight function $\bm{g}(t)=(g_{1}(t),...,g_{p}(t))^{\top}$. The weight function $\bm{g}_{n}(t)$ is assumed to belong to a class $\mathcal{G}$ with

	$\displaystyle\mathcal{G}$	$\displaystyle=\{\bm{f}(t):[0,1]\to\mathbb{R}^{p},~{}\bm{f}=(f_{1},...,f_{p})^{\top}~{}\text{is a continuous}~{}p\text{-dimensional}$
		$\displaystyle\text{~{}vector~{}function~{}satisfying}\inf_{t\in[0,1]}\min_{1\leq j\leq p}f_{j}(t)\geq\kappa\text{~{}for some constant}~{}\kappa>0\}.$

For a given $\alpha\in(0,1)$, denote by $\xi_{n,\bm{g}_{n}}(\alpha)$ the $(1-\alpha)$-th quantile of $\Xi_{n,\bm{g}_{n}}$. Then a JSCB with coverage probability $1-\alpha$ can be constructed as $\widetilde{\beta}_{j}(t)\pm\xi_{n,\bm{g}_{n}}(\alpha)g_{nj}(t)/\sqrt{n}$, $t\in[0,1]$, $j=1,2,...,p$.

Observe that the width of the JSCB is proportional to the weight function $\bm{g}_{n}(t)$. In practice one could simply choose some fixed weight functions such as $g_{nj}(t)\equiv 1$ which yields equal JSCB width at each $t$ and $j$. Alternatively, when the sample size is sufficiently large and temporal dependence is weak or moderately strong, we recommend choosing $g_{nj}(t)={\rm Std}(\widetilde{\beta}_{j}(t))/\int_{0}^{1}{\rm Std}(\widetilde{\beta}_{j}(s))\mathrm{d}s,~{}j=1,...,p.$ There are two advantages for this data-driven choice of weights. First, the resulting width of the JSCB reflects the standard deviation of $\widetilde{\bm{\beta}}(t)$ which gives direct visual information on the estimation uncertainty at every $t\in[0,1]$ and $j=1,...,p$. Second, this choice of weight function yields much smaller average width of the JSCB compared to some fixed choices such as $\bm{g}_{n}(t)\equiv 1$; see our simulations in Section 4 for a finite-sample illustration. On the other hand, ${\rm Std}(\widetilde{\beta}_{j}(t))$ has to be estimated in practice. Later in this article, we shall discuss its estimation and also asymptotic robustness of our multiplier bootstrap methodology when the weight function is inconsistently estimated.

To motivate the multiplier bootstrap, denote $\widetilde{\bm{\Sigma}}_{c}(\lambda)=\bm{X}_{c}^{\top}\bm{X}_{c}/n+\bm{R}(\lambda)$. Note that $\widetilde{\bm{\Sigma}}_{c}(\lambda)\approx\frac{\mathbb{E}(\bm{X}_{c}^{\top}\bm{X}_{c})}{n}+\bm{R}(\lambda):=\bm{\Sigma}_{c}(\lambda)$ under mild conditions. For simplicity, we assume that each $\beta_{j}(t)$ has the same degree of smoothness, here we let each $c_{j}$ have the same rate of divergence. By elementary calculations and basis expansions, we have

$$\widetilde{\bm{\beta}}(t)-\bm{\beta}(t)=\bm{C}_{f}(t)\widetilde{\bm{\Sigma}}_{c}^{-1}(\lambda)\frac{\bm{X}_{c}^{\top}\widetilde{\bm{\epsilon}}}{n}-\bm{C}_{f}(t)\widetilde{\bm{\Sigma}}_{c}^{-1}(\lambda)\bm{R}(\lambda)\bm{\theta}_{c}+\mathcal{O}(c^{-d_{1}}),$$

(5)

where $\widetilde{\bm{\epsilon}}=(\epsilon_{1}+\mathcal{O}_{\mathbb{P}}(c^{-(d_{1}+d_{2}+1)}),...,\epsilon_{n}+\mathcal{O}_{\mathbb{P}}(c^{-(d_{1}+d_{2}+1)}))^{\top}$. Hence if $c$ is sufficiently large and $\lambda$ is relatively small such that $\widetilde{\bm{\beta}}(t)$ is under-smoothed, i.e., the standard deviation of the estimation (captured by the first term on the right hand side of (5)) dominates the bias asymptotically, Eq.(5) reveals that the maximum deviation of $\sqrt{n}\widetilde{\bm{\beta}}(t)$ on $[0,1]$ is determined by the uniform probabilistic behavior of $\bm{Q}_{n}^{z}(t,\lambda):=\bm{C}_{f}(t)\bm{\Sigma}_{c}^{-1}(\lambda)\bm{Z}_{n}^{c}$, where $\bm{Z}_{n}^{c}:=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\bm{z}_{ci}$ and $\bm{z}_{ci}=\bm{x}_{ci}\epsilon_{i}$ with $\bm{x}_{ci}$ being the $i$-th column of $\bm{X}_{c}^{\top}$.

There are two major difficulties in the investigation of $\bm{Q}_{n}^{z}(t,\lambda)$ uniformly in $t$. Firstly, $\{\bm{z}_{ci}\}_{i\in\mathbb{Z}}$ is typically a moderately high-dimensional time series whose dimensionality $c$ diverges slowly with $n$ and $\bm{Q}_{n}^{z}(t,\lambda)$ is not a tight sequence of stochastic processes on $[0,1]$. Consequently, deriving the explicit limiting distribution of the maximum deviation of $\bm{Q}_{n}^{z}(t,\lambda)$ is a difficult task. Second, the convergence rate of $\sup_{t\in[0,1]}|\bm{Q}_{n}^{z}(t,\lambda)|_{\bm{g}_{n}(t)}$ depends on many nuisance parameters such as the smoothness of $\bm{X}_{i}(t)$ and $\bm{\beta}(t)$, and the diverging rate of the truncation parameters $c_{j}$, which are difficult to estimate in practice. To circumvent the aforementioned difficulties, one possibility is to utilize certain bootstrap methods to avoid deriving and estimating the limiting distributions and nuisance parameters explicitly. In this article, we resort to the multiplier/wild/weighted bootstrap ([57]) to mimic the probabilistic behavior of the process $\bm{Q}_{n}^{z}(t,\lambda)$ uniformly over $t$. In the literature, the multiplier bootstrap has been used for high dimensional inference on hyper-rectangles and certain classes of simple convex sets that can be well approximated by (possibly higher dimensional) hyper-rectangles after linear transformations; see for instance [12] and [14] for independent data and [64] for functional time series. The inference of $\sup_{t\in[0,1]}|\bm{Q}_{n}^{z}(t,\lambda)|_{\bm{g}_{n}(t)}$ uniformly over all quantiles and weight functions in $\mathcal{G}$ can be transformed into investigating the probabilistic behavior of $\bm{Z}_{n}^{c}$ over a large class of moderately high-dimensional convex sets. However, these convex sets have complex geometric structures for which results that are based on approximations on hyper-rectangles and their linear transformations are not directly applicable. As a result, in this article we shall extend the uniform Gaussian approximation and comparison results over all high-dimensional convex sets for sums of independent and $m$-dependent data established in, for instance, [5], [20] and [19] to sums of stationary and short memory time series in order to validate the multiplier bootstrap. These results may be of wider applicability in other moderately high-dimensional time series problems.

To be more specific, we will consider the bootstrapped sum given a block size $m$: $\bm{U}_{n}^{boots}=\frac{1}{\sqrt{n-m+1}}\sum_{j=1}^{n-m+1}\left(\frac{1}{\sqrt{m}}\sum_{i=j}^{j+m-1}\bm{z}_{ci}\right)u_{j},$ where $\{u_{j}\}_{j=1}^{n-m+1}$ is a sequence of i.i.d. standard normal random variables which is independent of $\bm{Z}_{1}^{n}:=\{\bm{z}_{c1},...,\bm{z}_{cn}\}$. Define $\bm{Q}_{n}^{boots}(t,\lambda)=\bm{C}_{f}(t)\widetilde{\bm{\Sigma}}_{c}^{-1}(\lambda)\bm{U}_{n}^{boots}$, then we will show that, conditional on the data, $\bm{U}_{n}^{boots}$ approximates $\bm{Z}_{n}^{c}$ in distribution in large samples with high probability. Naturally, we can use the conditional distribution of $\sup_{t\in[0,1]}|\bm{Q}_{n}^{boots}(t,\lambda)|_{\bm{g}_{n}(t)}$ to approximate the law of $\sup_{t\in[0,1]}|\bm{Q}_{n}^{z}(t,\lambda)|_{\bm{g}_{n}(t)}$ uniformly over all quantiles and weight functions in ${\cal G}$.

In this subsection, we will first discuss the issue of tuning parameter selection for the roughness penalization regression. We have three parameters to choose, that is the auxiliary truncation parameters $c_{j}$, the smoothing parameter $\lambda$ and the window size $m$. We recommend choosing $c_{j}=2d_{j}$ where $d_{j}$ can be selected via the cumulative percentage of total variance (CPV) criterion. Specifically, we choose $d_{j}$ such that the quantity $\sum_{k=1}^{d_{j}}\rho_{k}/\sum_{k=1}^{\infty}\rho_{k}$ exceeds a pre-determined high percentage value (e.g., 85% used in the simulations), where $\{\rho_{k}\}_{k=1}^{\infty}$ are the eigenvalues of ${\rm Cov}(X_{ij}(s),X_{ij}(t))$. The rationale is that with the aid of the roughness penalization, $c_{j}$ can be chosen at a relatively large value to reduce the sieve approximation bias without blowing up the variance of the estimation.

In addition, the generalized cross validation (GCV) method can be used to choose $\lambda$, see for examples [8, 47]. To be more specific, the GCV criterion for the smoothing parameter is defined as ${\rm GCV}(\lambda)=\frac{1}{n}\sum_{i=1}^{n}\frac{(Y_{i}-\widehat{Y}_{i})^{2}}{(1-{\rm Trace}(\bm{H})/n)^{2}},$ where $\bm{H}=\bm{X}_{c}[\bm{X}_{c}^{\top}\bm{X}_{c}/n+\bm{R}(\lambda)]^{-1}\bm{X}_{c}^{\top}/n$ and $\widehat{Y}_{i}$ is the $i$-th element of the vector $\widehat{\bm{Y}}=\bm{H}\bm{Y}$. Thus one can select $\lambda$ over a range by minimizing the above function.

For the window size $m$, we suggest using the minimum volatility (MV) method, which was proposed by [44]. Denote the estimated conditional covariance matrix $\widehat{\bm{\Xi}}^{c}=\widehat{\bm{\Xi}}^{c}(m)=\frac{1}{(n-m+1)m}\sum_{j=1}^{n-m+1}\left(\sum_{i=j}^{j+m-1}\widehat{\bm{z}}_{ci}\right)\left(\sum_{i=j}^{j+m-1}\widehat{\bm{z}}_{ci}^{\top}\right).$ The rationale behind the MV method is that the estimator $\widehat{\bm{\Xi}}^{c}(m)$ becomes stable as a function of $m$ when $m$ is in an appropriate range. Let the grid of candidate window sizes be $\{m_{1},...,m_{M_{1}}\}$. The MV criterion selects window size $m_{j_{0}}$ such that it minimizes the function $L(m_{j}):={\rm SE}\left(\left\{\widehat{\bm{\Xi}}^{c}(m_{j+k})\right\}_{k=-2}^{2}\right)$, where ${\rm SE}$ denotes the standard error

${\rm SE}\left(\left\{\widehat{\bm{\Xi}}^{c}(m_{j+k})\right\}_{k=-2}^{2}\right)=\left[\frac{1}{4}\sum_{k=-2}^{2}\left|\widehat{\bm{\Xi}}^{c}(m_{j+k})-\bar{\bm{\Xi}}^{c}(m_{j})\right|_{F}^{2}\right]^{1/2},$
with $\bar{\bm{\Xi}}^{c}(m_{j})=\sum_{k=-2}^{2}\widehat{\bm{\Xi}}^{c}(m_{j+k})/5$.

Next, we will describe the detailed steps of the multiplier bootstrap procedure for JSCB construction when the weight function $g_{nj}(t)={\rm Std}(Q_{nj}^{z}(t,\lambda))/\int_{0}^{1}{\rm Std}(Q_{nj}^{z}(s,\lambda))\,ds$, $j=1,...,p$. Note that ${\rm Std}(\sqrt{n}{\widetilde{\beta}_{j}(t)})/{\rm Std}(Q_{nj}^{z}(t,\lambda))\rightarrow 1$ in probability under some mild conditions.

1. (a)

   Select the window size $m$, such that $m\to\infty,~{}m=o(n)$.

2. (b)

   Choose the number of basis expansion $c_{j}$ for each $j=1,...,p$ and choose the smoothing parameter $\lambda$.

3. (c)

   Estimate FLM (1) and obtain the residuals $\widehat{\epsilon}_{i}=Y_{i}-\sum_{j=1}^{p}\sum_{k=1}^{c_{j}}\widetilde{\theta}_{jk}x_{ij,k}$.

4. (d)

   Generate $B$ (say 1000) sets of i.i.d. standard normal random variables $\{u^{(r)}_{j}\}_{j=1}^{n-m+1}$, $r=1,2,...,B$. For each $r=1,2,...,B$, calculate $\widehat{\bm{Q}}_{n,r}^{boots}(t,\lambda):=\widehat{\bm{C}}_{f}(t)\widehat{\bm{\Sigma}}_{c}^{-1}(\lambda)\widehat{\bm{U}}_{n,r}^{boots}$, where $\widehat{\bm{U}}_{n,r}^{boots}=\frac{1}{\sqrt{n-m+1}}\sum_{j=1}^{n-m+1}\left(\frac{1}{\sqrt{m}}\sum_{i=j}^{j+m-1}\widehat{\bm{z}}_{ci}\right)u^{(r)}_{j}$ with $\widehat{\bm{z}}_{ci}=\bm{x}_{ci}\widehat{\epsilon}_{i}$, $\widehat{\bm{C}}_{f}(t),\widehat{\bm{\Sigma}}_{c}(\lambda)$ have similar definitions to $\bm{C}_{f}(t)$ and $\widetilde{\bm{\Sigma}}_{c}(\lambda)$ with $f_{jk}$ replaced by its estimate $\widehat{f}_{jk}$.

5. (e)

   Estimate the bootstrap sample standard deviation of $\{\widehat{Q}_{nj,r}^{boots}(t,\lambda)\}_{r=1}^{B}$, $\widehat{{\rm Std}}(\widehat{Q}_{nj}^{boots}(t,\lambda))$, where $\widehat{Q}_{nj,r}^{boots}(t,\lambda)$ denotes the $j$-th component of $\widehat{\bm{Q}}_{n,r}^{boots}(t,\lambda)$, $j=1,...,p$. Calculate an estimator of $g_{nj}(t)$ as $\widehat{g}_{nj}(t):=\widehat{{\rm Std}}(\widehat{Q}_{nj}^{boots}(t,\lambda))/\int_{0}^{1}\widehat{{\rm Std}}(\widehat{Q}_{nj}^{boots}(s,\lambda))\,ds$ and obtain $M_{r}=\sup_{t\in[0,1]}|\widehat{\bm{Q}}_{n,r}^{boots}(t,\lambda)|_{\widehat{\bm{g}}_{n}(t)}$, $r=1,2,...,B$ with $\widehat{\bm{g}}_{n}(t)=(\widehat{g}_{n1}(t),\cdots,\widehat{g}_{np}(t))^{\top}$.

6. (f)

   For a given level $\alpha\in(0,1)$, let the $(1-\alpha)$-th sample quantile of the sequence $\{M_{r}\}_{r=1}^{B}$ be $\hat{q}_{n,1-\alpha}$. Then the JSCB of $\bm{\beta}(t)$ can be constructed as $\widetilde{\beta}_{j}(t)\pm\widehat{g}_{nj}(t)\hat{q}_{n,1-\alpha}/\sqrt{n}$ for $t\in[0,1]$, $j=1,2,...,p$.

In the rare case where $\widehat{g}_{nj}(t_{0})$ is close to 0 at some $t_{0}$, one can lift up $\widehat{g}_{nj}(t_{0})$ to a certain threshold (say, $\max_{t\in[0,1]}\widehat{{\rm Std}}(\widehat{Q}_{nj}^{boots}(t,\lambda))/[100\int_{0}^{1}\widehat{{\rm Std}}(\widehat{Q}_{nj}^{boots}(s,\lambda))\,ds]$) while keeping the weight function continuous such that $\widehat{g}_{nj}(t)\in{\cal G}$. As we will show in Section 4, the above manipulations do not influence the asymptotic validity of the bootstrap.

If one is interested in constructing JSCB for a group of parameter functions, say $\beta_{i_{1}}(t),\cdots,\beta_{i_{k}}(t)$, then one just need to focus on the $i_{1}$th, $i_{2}$th, $\cdots$, $i_{k}$th elements of the bootstrap process $\widehat{\bm{Q}}_{n,r}^{boots}(t,\lambda)$, $[\widehat{Q}_{ni_{1},r}^{boots}(t,\lambda),\widehat{Q}_{ni_{2},r}^{boots}(t,\lambda),\cdots,\widehat{Q}_{ni_{k},r}^{boots}(t,\lambda)]^{\top}$, to conduct simultaneous inference of those parameter functions. The implementation procedure is very similar to the above and we shall omit the details.

Based on the basis expansion (2), we aim at utilizing a general time series model for $\{x_{ij,k}\}_{k=1}^{\infty}$ from a nonlinear system point of view as follows, which will serve as a preliminary for our theoretical investigations.

Note that in the above definition $\delta_{x}(l,q)$ does not depend on $i$. The above formulation of time series $x_{ij,k}$ can be viewed as a physical system where functions $G_{jk}$ are the underlying data generating mechanisms and $\{\eta_{i}\}$ are the shocks or innovations that drive the system. Meanwhile, $\delta_{x}(l,q)$ measures the temporal dependence of $\{\bm{X}_{i}(t)\}$ by quantifying the corresponding changes in the system’s output uniformly across all basis expansion coefficients when the shock of the system $l$ steps ahead is changed to an i.i.d. copy. We refer to [58] for more discussions of the physical dependence measures with examples on how to calculate them for a wide range of linear and nonlinear time series models.

Definition 1 is related to the class of functional time series formulated in [64]. The difference is that in (2) we separate the standard deviation $f_{jk}$ from $\widetilde{x}_{ij,k}$ and the functional time series model in [64] is formulated without this extra step. Standardization of the basis expansion coefficients is needed in the fitting of the FLM (1) to avoid near singularity of the design matrix. Furthermore, Definition 1 is also related to the concept of $m$-approximable functional time series introduced in [26] as both formulations utilize the concepts of Bernoulli shifts and coupling. The difference lies in our adaptation of the basis expansion similar to that in [64] which separates the functional index $t$ and time index $i$ and hence makes it easier technically to investigate the behavior of various estimators of $\bm{\beta}(t)$ uniformly over $t$. Now, we impose an assumption on the speed of decay for the dependence measures $\delta_{x}(l,q)$.

Assumption 3.1 is a mild short-range dependence assumption which asserts that the temporal dependence of the functional time series $\bm{X}_{i}(t)$ decays at a sufficiently fast polynomial rate. For independent functional data, the condition $\delta_{x}(l,q)\leq C_{2}(l+1)^{-\tau}$ is automatically satisfied as there is no temporal dependence ($\delta_{x}(l,q)=0$). In Section C of the online supplemental material, we will provide two examples on how to calculate $\delta_{x}(l,q)$ for a class of functional MA$(\infty)$ and functional AR(1) processes, respectively.

Adapting the nonlinear system point of view ([58]), we model the stationary time series of errors $\{\epsilon_{i}\}_{i=1}^{n}$ as $\epsilon_{i}=G(\mathcal{F}_{i})$ for some measurable function $G$. Observe that both $\bm{X}_{i}(t)$ and $\epsilon_{i}$ are generated by the same set of shocks $\{\eta_{j}\}_{j\in\mathbb{Z}}$ and hence they could be statistically dependent. To establish the main results, we need the following conditions:

The above assumptions are mild and are needed for establishing a Gaussian approximation and comparison theory for roughness-penalized estimators. 3.2 ensures positive definiteness of the design matrix in order to avoid multi-colinearity. 3.3 is a short range dependent condition on $\{\epsilon_{i}\}_{i=1}^{n}$ in accordance with 3.1. Finally, 3.4 puts some moment restrictions on the random variable $\bm{z}_{ci,j}$.

Let $\widetilde{\bm{R}}_{j}$ be a $c_{j}\times c_{j}$ matrix with its $(k,l)$ element $\widetilde{R}_{j}(k,l)$ $=\int_{0}^{1}\alpha_{k}^{\prime\prime}(t)\alpha_{l}^{\prime\prime}(t)\mathrm{d}t$. The following additional assumptions are needed for the validation of the multiplier bootstrap.