The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
71,804 characters
Simultaneous Inference of a Partially Linear Model in Time Series
\maketitle
\begin{abstract}
We introduce a new methodology to conduct simultaneous inference of the nonparametric component in partially linear time series regression models where the nonparametric part is a multivariate unknown function. In particular, we construct a simultaneous confidence region (SCR) for the multivariate function by extending the high-dimensional Gaussian approximation to dependent processes with continuous index sets. Our results allow for a more general dependence structure compared to previous works and are widely applicable to a variety of linear and nonlinear autoregressive processes. We demonstrate the validity of our proposed methodology by examining the finite-sample performance in the simulation study. Finally, an application in time series, the forward premium regression, is presented, where we construct the SCR for the foreign exchange risk premium from the exchange rate and macroeconomic data.
\vspace{1cm}
\noindent\textit{Keywords:} time series, simultaneous inference, simultaneous confidence region, partially linear model, Gaussian approximation, forward premium regression
\vspace{1cm}
\end{abstract}
\setlength{\parindent}{15pt}
\newtheorem{definition}{Definition}
\newtheorem{assumption}{Assumption}
\newtheorem{theorem}{Theorem}
\newtheorem{proposition}{Proposition}
\newtheorem{corollary}{Corollary}
\newtheorem{lemma}{Lemma}
\newtheorem{remark}{Remark}
\newtheorem{example}{Example}
\newpage
\section{Introduction}\label{sec_intro}
Partially linear models are of interest in many practical problems. For example, in econometrics, \textcite{engle_semiparametric_1986} modeled the electricity sales as the combination of a smooth function of temperature and a linear function of price and income; in materials science, \textcite{green_semi-parametric_1985} used a semi-parametric generalized linear model to analyze the bioassay data for the study of flame retardants; in biology, \textcite{liang_empirical_2009} applied generalized partially linear models to investigate the relationship between viral load and CD4$^+$ cell counts to understand AIDS pathogenesis. See other applications in \textcite{hardle_partially_2000}.
In this paper, we consider a partially linear time series regression model
\begin{equation}
\label{eq_model}
Y_i=Z_i^{\top}\boldsymbol \beta+\mu(X_i)+\sigma(X_i)\epsilon_i,\quad i=1,\ldots,n,
\end{equation}
where $(Z_i,X_i,Y_i)$ are observed stationary processes with $Z_i\in\mathbb R^{l}$, $X_i\in\mathbb R^d$ and $Y_i\in\mathbb R$, for $l,d\ge1$. Here $\boldsymbol \beta\in\mathbb R^{l}$ is a fixed vector of unknown parameters and $\mu(\cdot)$ [resp. $\sigma^2(\cdot)$] is an unknown smooth regression function (resp. conditional variance or volatility function) from $\mathbb R^d$ to $\mathbb R$. In addition, $\epsilon_i\in\mathbb R$ is an unobserved random error with mean zero, independent of the covariates $X_i$ and $Z_i$. In this work, we'll primarily concentrate on the conditional volatility $\sigma(X_i)$ for clarity's sake although our findings can be extended to $\sigma(X_i,Z_i)$. Compared to completely parametric or nonparametric specifications, a partially linear model in (\ref{eq_model}) enjoys a flexible semi-parametric structure. The parametric components can provide easier interpretations of each variable to better characterize the underlying data-generating mechanism, while the additional nonparametric part allows a data-driven approximation with no specific structures imposed on the true regression function, which can avoid inconsistent estimators and faulty inferences due to model mis-specification in purely parametric models.
In the past decades, much attention has been directed to estimating and testing partially linear models. See, for instance, \textcite{engle_semiparametric_1986,rice_convergence_1986,Robinson,speckman_kernel_1988,schick_root-n_1996} on the $\sqrt{n}$-consistent estimators of $\boldsymbol \beta$; \textcite{gao_convergence_1995,fan_profile_2005,xie_scad-penalized_2009} on the inferences of $\boldsymbol \beta$. It is crucial to include the parametric component in equation (\ref{eq_model}) because the parametric component is derived out of relevant theories and it is practically useful to identify the parametric component. One notable example is the Phillips curve with a time-varying natural unemployment rate (\cite{KHK:2014}), where the parameter $\boldsymbol \beta$ captures the effect of the unemployment rate on the price inflation. Moreover, $\boldsymbol \beta$ can be employed to measure the impact of various demographic variables, such as age, gender, family size, and the residency type, on the household gasoline consumption, as demonstrated in the U.S. case study by \textcite{kim_simultaneous_2021}. In fact, this paper significantly extends the scope of both \textcite{KHK:2014} and \textcite{kim_simultaneous_2021} by introducing a multivariate, time-dependent, and stochastic nonlinear component into (\ref{eq_model}). Furthermore, the parameter $\boldsymbol \beta$ in equation (\ref{eq_model}) can represent the factor that determines the efficiency of the foreign currency market, as discussed in Section \ref{sec_app}. Considering these diverse roles played by the parametric component, it is essential and useful to include the linear parametric part for the estimation and the statistical inference of model (\ref{eq_model}), instead of relying solely on the nonparametric part. It is worth noting that the purely nonparametric model corresponds to a special case of model (\ref{eq_model}) where $\boldsymbol \beta$ equals zero.
The estimation of the nonparametric part has also been intensively studied, including the methods based on kernel, local linear and spline smoothers (\cite{hamilton_local_1997,yu_penalized_2002,fan_kernel-based_2003,aneiros-perez_local_2008}). Several attempts have been made to the inference of $\mu$ in partially linear models, such as
the consistency and asymptotic normality for the estimator of $\mu$ by \textcite{liang_asymptotic_1997}, the point-wise confidence intervals of $\mu$ based on empirical likelihood by \textcite{liang_empirical_2009}, and simultaneous confidence bands of multivariate function $\mu$ by \textcite{KHK:2014,kim_simultaneous_2021}. However, all the aforementioned literature focused on the {\it independent} or non-stochastic observations. No previous work investigated the simultaneous inference of $\mu$ in a partially linear time series model with {\it dependence} as (\ref{eq_model}) that is commonly encountered in real data (\cite{hardle_partially_2000}). Further, most of the studies on partially linear models assume the error terms in model (\ref{eq_model}) to be homoskedastic, where the error term $\sigma(X_i)\epsilon_i$ is simply reduced to $\epsilon_i$ and the conditional variance is constant over time. This can be a shortcoming since it would rule out most macroeconomic and financial time series, at least for the application to asset pricing (\cite{nelson_conditional_1991}), where returns may be uncorrelated but feature stochastic volatility. The current paper aims to fill in these gaps by providing theory for the simultaneous inference of the mean trend $\mu$ in (\ref{eq_model}) under a general dependency structure, which also allows for conditional heteroskedasticity.
Specifically, we allow both the errors $\epsilon_i$ and the covariates $X_i$ in (\ref{eq_model}) to be dependent over $i$. Let $X_i=(X_{i1},X_{i2},\ldots,X_{id})^{\top}$ be a stationary process of the form
\begin{equation}
\label{eq_X_structure}
X_i = H(\ldots,v_{i-1},v_i),
\end{equation}
where $v_i$ are independent and identically distributed (i.i.d.) random vectors in $\mathbb R^{d'}$ for
$d'\ge1$ and $H=(H_1,H_2,\ldots,H_d)^{\top}$ is a measurable function such that $X_i$ is a well-defined. The nonlinear Wold representation (\ref{eq_X_structure}) allows a very general class of stationary processes, including linear processes such as vector autoregressive models (VAR) and autoregressive moving average (ARMA) models and nonlinear transforms such as bilinear models, Volterra processes, Markov chain models, threshold/exponential autoregressive models (TAR/EAR) and (generalized) autoregressive conditionally heteroscedastic (ARCH/GARCH) type models, etc. Within this framework, $v_i$ can be viewed as independent inputs of a physical system, and all the dependencies among the outputs $X_i$ result from the underlying data-generating mechanism $H(\cdot)$.
For the identification of model (\ref{eq_model}), we assume that the error $\epsilon_i$ is independent of both covariates $X_i$ and $Z_i$. In particular, we shall proceed with the conditional expectation $\mathbb E(Y_i\mid X_i,Z_i)=Z_i^{\top}\boldsymbol \beta + \mu(X_i)$. Following the fixed design case in \textcite{hardle_partially_2000}, we let $Z_i$ be a function of $X_i$ plus another noise term (cf. Assumption \ref{asm_iden}). This means $\mathbb E(Y_i\mid X_i,Z_i)$ cannot be reduced to $\mathbb E(Y_i\mid X_i)$ and indicates that it is nontrivial and also challenging to extend the inference of $\mu(\cdot)$ in a purely nonparamteric model to that under a partially linear setting. If $Z_i=0$, then model (\ref{eq_model}) is simply a nonparametric regression process as a special case, that is,
\begin{equation}
\label{eq_model_example}
Y_i= \mu(X_i)+\sigma(X_i)\epsilon_i,\quad i=1,\ldots,n.
\end{equation}
When $X_i=Y_{i-1}$ and $\epsilon_i$ are i.i.d. random noises, model (\ref{eq_model_example}) incorporates many interesting linear and nonlinear autoregressive processes (AR), such as AR processes if $\mu(x)=ax$ for some real parameter $a$, and autoregressive conditional heteroscedastic (ARCH) processes if $\mu(x)=0$ and $\sigma^2(x)=\alpha_0+\alpha_1x^2$ for some non-negative real parameters $\alpha_0,\,\alpha_1\in\mathbb R$.
Many contributions have been made to
constructing the SCR of $\mu(\cdot)$ in completely nonparametric models. For example, concerning independent data, \textcite{johnston_probabilities_1982} was among the first to investigate the inferences of univariate mean regression functions; \textcite{hardle_asymptotic_1989} derived simultaneous confidence bands for one-dimensional kernel M-estimators; \textcite{HS:2010,GH:2012} constructed uniform confidence bands for conditional quantile and expectile functions, respectively. With regard to dependent cases, see inference of trends in a fixed design with $X_i=i/n$ by \textcite{WZ:2007}; confidence bands for the mean function in functional time series with physical dependence by \textcite{CS2015}; nonlinear regression model with nonstationary regressors by \textcite{Li2017} and a time-varying nonlinear regression model by \textcite{ZW2015}. In particular, \textcite{ZW:2008,LW:2010} proposed inferences of the univariate mean and volatility functions in a similar time series regression model in (\ref{eq_model_example}), where they assumed that the error terms $\epsilon_i$ are i.i.d.. Our work can be viewed as a generalization of their results by extending the dependence structure of $\epsilon_i$ and by including an additional parametric part to accommodate a broader class of data-generating mechanisms.
That is, our work is distinct from \textcite{ZW:2008,LW:2010} in that our model framework is semi-parametric with the multivariate covariate $X_i$ and time-dependent $\epsilon_i$, while the framework in \textcite{ZW:2008,LW:2010} is purely nonparametric with an univariate $X_i$ and an i.i.d. noise $\epsilon_i$.
\textbf{Contributions:} Here we summarize our three main contributions to the literature: Firstly, we extend the dependence structure of the error terms to a more general case by allowing $\epsilon_i$ to be {\it dependent over} $i$, while also accounting for the dependence among the covariates $X_i$ and the conditional heteroscedasticity. Secondly, different from the relevant studies relying on the Gumbel convergence to achieve the asymptotics of the statistics (\cite{zhao_kernel_2006,LW:2010}), we provide a new testing methodology based on the multiplier bootstrap enlightened by \textcite{chernozhukov_central_2017}. This allows one to {\it avoid the notoriously slow convergence issue} associated with the Gumbel distribution.
Thirdly, there is no previous work performing simultaneous inference of the {\it multivariate} $\mu(\cdot)$ in (\ref{eq_model}), allowing the multivariate covariate $X_i\in\mathbb{R}^d$ with $d\geq 2$ under some general {\it time dependence} setting. This paper could be a complement to the non-parametric model validation problem with a general dependence structure and conditional heteroscedasticity, applicable in various multivariate scenarios.
\textbf{Notation:} For a vector $v=(v_1,...,v_d)\in\mathbb R^d$ and $q>0$, we denote $|v|_q=(\sum_{i=1}^d|v_i|^q)^{1/q}$ and $|v|_{\infty}=\max_{1\le i\le d}|v_i|$. For $s>0$ and a random vector $X$, we say $X\in\mathcal L^s$ if $\lVert X\rVert_s=[\mathbb E(|X|_2^s)]^{1/s}<\infty$. For two positive number sequences $(a_n)$ and $(b_n)$, we say $a_n=O(b_n)$ or $a_n\lesssim b_n$ (resp. $a_n\asymp b_n$) if there exists $C>0$ such that $a_n/b_n\le C$ (resp. $1/C\le a_n/b_n\le C$) for all large $n$, and say $a_n=o(b_n)$ if $a_n/b_n\rightarrow0$ as $n\rightarrow\infty$. We set $(X_n)$ and $(Y_n)$ to be two sequences of random variables. Write $X_n=O_{\mathbb P}(Y_n)$ if for $\forall \epsilon>0$, there exists $C>0$ such that $\mathbb P(|X_n/Y_n|\le C)>1-\epsilon$ for all large $n$, and say $X_n=o_{\mathbb P}(Y_n)$ if $X_n/Y_n\rightarrow 0$ in probability as $n\rightarrow\infty$. We denote the centered random variable $X$ by $\mathbb E_0(X)$, that is, $\mathbb E_0(X)=X-\mathbb E(X)$.
\textbf{Roadmap:} The rest of the paper is structured as follows. Section \ref{sec_SCR} introduces the overall methodology to perform simultaneous inference of $\mu(\cdot)$ in (\ref{eq_model}) with $d\geq 1$. The asymptotic properties of the proposed statistics and estimators as well as the implementation are provided in Sections \ref{sec_asym} and \ref{sec_est}. Section \ref{sec_simul} is devoted to a simulation study to evaluate the performance of our methods and Section \ref{sec_app} offers an empirical application, the forward premium anomaly, to demonstrate the validity of the proposed methodology in practice. Section \ref{conclusion} concludes the paper and discusses potential extensions for future research. The technical proofs are deferred to the Supplementary Materials.
\section{Simultaneous Confidence Region (SCR)}\label{sec_SCR}
In this section, we first introduce the definition of the simultaneous confidence region (SCR). Then, we shall follow with the estimator of the nonparametric trend function $\mu(\cdot)$ in (\ref{eq_model}). Further, we illustrate our new methodology on constructing the SCR of $\mu(\cdot)$ based on this estimated $\mu(\cdot)$. The theoretical intuition of the proposed simultaneous inference is also provided.
To conduct simultaneous inference of the trend $\mu(\cdot)$ in model (\ref{eq_model}), we shall construct the nonparametric simultaneous confidence region (SCR) for $\mu(\cdot)$.
In particular, we consider deriving asymptotic SCR for $\mu(\cdot)$ over the region $\mathcal T_d=[T_{11},T_{12}]\times [T_{21},T_{22}]\times \cdots \times [T_{d1},T_{d2}] \subset \mathbb R^d$ with confidence level $100(1-\alpha)\%$, $\alpha\in(0,1)$. To this end, we shall find two functions $l_n(\cdot)$ and $r_n(\cdot)$ based on the observations $(Z_i,X_i,Y_i)$, $1\le i\le n$, such that
\begin{equation}
\label{eq_UCBdef}
\lim_{n\rightarrow\infty}\mathbb P\Big(l_n(x)\le \mu(x)\le r_n(x), \text{ for all } x\in\mathcal T_d\Big)=1-\alpha.
\end{equation}
Given the SCR for $\mu(\cdot)$, we can verify whether $\mu(\cdot)$ is of some certain parametric form by testing the null hypothesis
\begin{equation}
\label{eq_testnull}
\mathcal H_0:\,\mu(\cdot)=\mu_{\theta}(\cdot),
\end{equation}
against the alternative $\mathcal H_{\mathcal A}:\, \mu(\cdot)\neq\mu_{\theta}(\cdot)$, where $\theta\in\Theta$ for some parametric space $\Theta$ and $\mu_{\theta}(\cdot)$ is a multivariate parametric function. Specifically, one can test (\ref{eq_testnull}) by checking whether the condition $l_n(x)\le\mu_\theta(x)\le r_n(x)$ holds for all $x\in\mathcal T_d$. If this condition does not hold for some $x\in\mathcal T_d$,
then we reject the null hypothesis at level $\alpha$.
The SCR-based inference is more preferred to other standard inferential procedures utilizing mean-integrated-squared-error (MISE) type statistics, for example, since it is more effective in suggesting the right function form of $\mu(\cdot)$ in (\ref{eq_model}). When the null hypothesis in (\ref{eq_testnull}) is rejected via an MISE-type test statistic, it would be rather difficult to figure out the reason for rejection, which, however, can be easily dealt with under our approach by locating graphically where the SCR is violated by $\mu_{\theta}(\cdot)$ under the null hypothesis.
Next, we provide an estimator for $\mu(\cdot)$ given the observed sample $(Z_i,X_i,Y_i)$, $1\le i\le n$. Let $x=(x_1,\ldots,x_d)^{\top}\in\mathbb R^d$ and set $K(\cdot)\ge0$ to be some kernel function with support $[-1,1]^d$. We consider the following optimization problem:
\begin{equation}
\label{eq_mu_hat_star_goal}
\hat\mu^*(x)=\text{argmin}_{\theta}n^{-1}\sum_{i=1}^nK_h(x-X_i)\big(Y_i-Z^{\top}_{i}\hat\boldsymbol \beta-\theta\big)^2,
\end{equation}
where $K_h(\cdot)=K(\cdot/h)/h^d$, $h$ is a bandwidth parameter with $h\rightarrow0$ and $h^dn\rightarrow\infty$, and $\hat\boldsymbol \beta$ is a consistent estimator of the unknown parameters $\boldsymbol \beta$ in (\ref{eq_model}). We shall defer the details of $\hat\boldsymbol \beta$ to Section \ref{sec_est}. Here in (\ref{eq_mu_hat_star_goal}), we adopt the local constant estimator for the simplicity of notation. One can achieve similar results by applying the local linear estimator introduced in \textcite{fan_local_1996}. By solving (\ref{eq_mu_hat_star_goal}), we can obtain the Nadaraya-Watson estimator for $\mu(\cdot)$ which has the expression
\begin{equation}
\label{eq_mu_hat_star}
\hat\mu^*(x) =\sum_{i=1}^nw_h(x,X_i)\big(Y_i-Z^{\top}_{i}\hat\boldsymbol \beta\big),
\end{equation}
where the weight function $w_h(x,X_i)$ is defined as
\begin{equation}
\label{eq_weight}
w_h(x,X_i) = \frac{K_h(x-X_i)}{\sum_{i=1}^nK_h(x-X_i)}.
\end{equation}
Moreover, we denote the consistent estimator of the volatility function $\sigma(\cdot)$ in (\ref{eq_model}) by $\hat\sigma(\cdot)$, and we shall provide the detailed definition and consistency results of $\hat\sigma(\cdot)$ in Section \ref{sec_est}. Given $\hat\boldsymbol \beta$ and $\hat\sigma(\cdot)$, we consider the statistic $\sup_{x\in\mathcal T_d}\big|\hat\mu^*(x)- \mu(x)\big|/\hat\sigma(x)$ to construct the SCR of $\mu(x)$. We shall note that, due to the smoothness of $\mu(\cdot)$ and the consistency of $\hat\boldsymbol \beta$, this statistic can be approximately written into the supremum of a sum of dependent random fields conditioned on the covariates $X_i$, that is
\begin{equation}
\label{eq_GA_form}
\sup_{x\in\mathcal T_d}\big|\hat\mu^*(x)- \mu(x)\big|/\hat\sigma(x) \approx \sup_{x\in\mathcal T_d}\Big|\sum_{i=1}^nw_h(x,X_i)\sigma(X_i)\epsilon_i\Big|/\hat\sigma(x).
\end{equation}
It is non-trivial to investigate the asymptotic properties of (\ref{eq_GA_form}) when $X_i$ and $\epsilon_i$ are dependent over $i$. \textcite{zhao_kernel_2006,LW:2010} have dealt with the case where the covariates $X_i$ are dependent while the errors $\epsilon_i$ are independent by establishing the Gumbel convergence. To address the more general dependency structure in our study, we propose an extension of the high-dimensional Gaussian approximation theorem introduced by \textcite{chernozhukov_central_2017} to dependent processes with continuous index sets. By this generalized high-dimensional Gaussian approximation, we shall expect the limit distribution of our proposed statistic to be approximated by the one of the maximum of a centered Gaussian random vector $\hat\mathcal Z=(\hat \mathcal Z_1,\ldots,\hat\mathcal Z_n)^{\top}\in\mathbb R^n$, that is
\begin{equation}
\label{eq_GA_intuitive}
\mathbb P\big(\sup_{x\in\mathcal T_d}\sqrt{h^dn}\big|\hat\mu^*(x)- \mu(x)\big|/\hat\sigma(x)<u\big) \approx \mathbb P\big(\max_{1\le j\le n}|\hat\mathcal Z_j|<u\big).
\end{equation}
We defer the detailed definition of the covariance matrix for $\hat\mathcal Z$ to (\ref{eq_Q_hat}). Intuitively, the result in (\ref{eq_GA_intuitive}) would enable us to find the critical value of our proposed statistic, and consequently facilitates the construction of the simultaneous confidence region for $\mu(\cdot)$. Specifically, we can approximate $l_n(x)$ and $r_n(x)$ in (\ref{eq_UCBdef}) by the estimators
\begin{equation}
\label{eq_ln_rn}
\hat l_n(x) = \hat\mu^*(x) - \hat q_{\alpha}\hat\sigma(x), \quad \hat r_n(x) = \hat\mu^*(x) + \hat q_{\alpha}\hat\sigma(x),
\end{equation}
respectively, where $\hat q_{\alpha}$ is the $(1-\alpha)$-th empirical quantile of $\max_{1\le j\le n}|\hat\mathcal Z_j|/\sqrt{h^dn}$ given the significance level $\alpha\in(0,1)$, and it can be evaluated by the multiplier bootstrap (\cite{chernozhukov_central_2017}). Based on the SCR in (\ref{eq_ln_rn}), one can test whether the trend function $\mu(\cdot)$ is of any particular parametric form $\mu_\theta(\cdot)$, such as quadratic or cubic patterns, by evaluating whether $l_n(x)\le\mu_\theta(x)\le r_n(x)$ is satisfied for all $x\in\mathcal T_d$. If the SCR fails to {\it entirely} contain this parametric form, then we reject the null hypothesis (\ref{eq_testnull}) at level $\alpha$. We shall provide the detailed steps for implementing the SCR construction at the end of Section \ref{sec_est} after we introduce our main theorems.
\section{Asymptotic Properties}\label{sec_asym}
This section is devoted to our main results on the asymptotic properties for the statistic $\sup_{x\in\mathcal T_d}\big|\hat\mu^*(x)- \mu(x)\big|/\hat\sigma(x)$, which provides the theoretical foundation for the construction of SCR. In Section \ref{subsec_GA}, we shall first establish the asymptotic distribution of the proposed statistic under the oracle setting, that is, assuming that the unknown parameters $\boldsymbol \beta$ and $\sigma(\cdot)$ in the statistic are the true ones. The case with $\boldsymbol \beta$ and $\sigma(\cdot)$ replaced by their consistent estimators $\hat\boldsymbol \beta$ and $\hat\sigma(\cdot)$, respectively, are dealt with in Section \ref{sec_est}, where we provide the consistency results for $\hat\boldsymbol \beta$ and $\hat\sigma(\cdot)$ and show the similar asymptotic distribution of the statistic.
\subsection{Technical Assumptions}\label{subsec_asm}
We shall start with some regularity conditions which will be useful to establish our main theorems. First, we impose the smoothness condition on the trend function $\mu(\cdot)$ and assume that the volatility function $\sigma(\cdot)$ also varies smoothly and is bounded on the support $\mathcal T_d$.
\begin{assumption}[Trend and volatility]\ \\
\label{asm_sigma}
(i) (Smoothness)
Assume that the trend function $\mu(\cdot)$ and volatility function $\sigma(\cdot)$ defined in (\ref{eq_model}) are both Lipschitz continuous on $\mathcal T_d$. \\
(ii) (Bounds) Assume that for some constants $c_{\sigma},c_{\sigma}'>0$, $c_{\sigma}\le \inf_{x\in\mathcal T_d}\sigma(x)\le \sup_{x\in\mathcal T_d}\sigma(x)\le c_{\sigma}'$.
\end{assumption}
\begin{assumption}[Kernel]\ \\
\label{asm_kernel}
(i) The kernel function $K(\cdot)$ in (\ref{eq_mu_hat_star_goal}) is defined on $\mathbb I=[-1,1]^d$ and is continuously differentiable up to order two. \\
(ii) Assume that $\sup_{x\in \mathbb I}|K(x)|<\infty$ and $\int_{\mathbb I} K(x)dx =1$. Also assume that $K(x)$ has first-order derivative with $ \sup_{x \in \mathbb I}\max_{1\le i\le d}|\partial K(x)/\partial x_i|<\infty$. \\
(iii) Assume that the bandwidth parameter $h\rightarrow0$ and $h^dn\rightarrow\infty$.
\end{assumption}
\begin{remark}[Choice of the bandwidth parameter $h$]
\label{remark_bandwidth}
The bandwidth parameter $h$ can be selected by some data-driven method such as the generalized cross-validation (GCV) introduced by \textcite{CG:1977}. Specifically, the bandwidth $h$ is firstly chosen by GCV and then adjusted downward for undersmoothing, which means $h$ is chosen to converge to 0 faster than the optimal rate yielded by GCV. We defer the detailed discussion of bandwidth selection to Remark \ref{rmk_thm1} after we introduce some additional conditions on $h$ in Theorem \ref{thm1_GA}. Furthermore, we show the results of our simulation study in Section \ref{sec_simul} where we benchmark the effects of multiple bandwidths $h$ on the coverage probabilities.
\end{remark}
Next, we shall specify the dependence structure of the error process $\{\epsilon_i\}_{1\le i\le n}$ and the covariates $X_1,\ldots,X_n$ in model (\ref{eq_model}). In particular, throughout this paper, we assume that $\{\epsilon_i\}_{1\le i\le n}$ is MA($\infty$), which can be formulated as follows:
\begin{equation}
\label{eq_epsilon}
\epsilon_i = \sum_{k=0}^{\infty}a_k\eta_{i-k},
\end{equation}
where $\eta_i\in\mathbb R$ are i.i.d. random variables with mean zero and unit variance, independent of $(X_i,Z_i)$ in (\ref{eq_model}). The coefficients $a_k$, $k\ge0$, take values in $\mathbb R$ such that $\epsilon_i$ is a proper random variable. We assume that the innovations $\eta_i$ have finite $q$-th moment, for some $q\ge 4$. We define the absolute sum of the coefficients as $S=\big|\sum_{k\ge0}a_k\big|$. Then, the long-run variance of $\{\epsilon_i\}_{1\le i\le n}$ is
\begin{equation}
\label{eq_longrun}
S^2=\Big(\sum_{k\ge0}a_k\Big)^2=\sum_{k=-\infty}^{\infty}\gamma(k),
\end{equation}
where $\gamma(k)=\mathbb E(\epsilon_i\epsilon_{i+k})$ is the autocovariance of innovations with lag $k\in\mathbb Z$.
\begin{assumption}[Finite moment]
\label{asm_moment}
Assume that innovations $\{\eta_k\}$ defined in (\ref{eq_epsilon}) are i.i.d. with $\lVert \eta_1\rVert_q<\infty$, for some constant $q\ge4$.
\end{assumption}
\begin{assumption}[Dependence of $\epsilon$]\label{asm_dep_epsilon}
Assume that for any integer $l\ge0$,
$$\sum_{k\ge l}|a_k|/S=O\big\{(1\vee l)^{-\zeta}\big\},$$
where the constant $\zeta>1$ and $S$ is the long-run standard deviation of $\{\epsilon_i\}_{1\le i\le n}$ defined in (\ref{eq_longrun}).
\end{assumption}
Assumptions \ref{asm_moment} and \ref{asm_dep_epsilon} post conditions on the moment and dependency structure of errors $\{\epsilon_i\}_{1\le i\le n}$. Specifically, the moment condition in Assumptions \ref{asm_moment} depends on $q$ which characterizes the heavy-tailedness of the noise, and a larger $q$ means a thinner tail. Assumption \ref{asm_dep_epsilon} requires that the dependency strength of $\{\epsilon_i\}_i$ decays at a polynomial rate, which also ensures that the long-run variance of $\{\epsilon_i\}_i$ is finite.
For the covariates $X_1,\ldots,X_n$, we denote the density function of $X_{i}=(X_{i1},\ldots, X_{id})^{\top}$ by $g(x_1,\ldots, x_d)$. For $s=(s_1,s_2,\ldots,s_d)^{\top} \in \mathbb{R}^d$, we define
\begin{equation}
\label{eq_density_x}
g(s\mid \mathcal{F}_{i-1})=\frac{\partial^d\mathbb P(X_{i}\leq s\mid \mathcal{F}_{i-1})}{\partial s_1\ldots \partial s_d},
\end{equation}
where the filtration $\mathcal{F}_i = (\ldots ,v_{i-1},v_{i})$, $i\in \mathbb Z$.
We let $v_i'$ be an i.i.d. copy of $v_i$ and $\mathcal{F}_{i,\{k\}}$ be $\mathcal{F}_{i}$ with $v_k$ therein replaced by $v_k'$. We define $X_{ij,\{i-k\}}$ as a coupled version of $X_{ij}$ with the form
$$X_{ij,\{i-k\}}=H_j(\ldots ,v_{i-k-1},v'_{i-k},v_{i-k+1},\ldots ,v_i).$$
In addition, we denote $X_{i,\{i-k\}}=(X_{i1,\{i-k\}},\ldots ,X_{id,\{i-k\}})^\top$. We impose assumptions on the moments and dependency structures of covariates $X_1,\ldots,X_n$ as follows.
\begin{assumption}[Covariates $X_i$] \ \\
\label{asm_dep}
(i) (Finite moment). Assume that the conditional density function $g(x\mid\mathcal F_{i-1})$ defined in (\ref{eq_density_x}) has finite $s$-th moment, i.e. $\|g(x\mid\mathcal F_{i-1})\|_s<\infty$, for some constant $s>2$. \\
(ii) (Dependence strength). Following \textcite{Wu:2005}, for all $k\ge1$, we define the physical dependence measure
$$\theta_{k,s} =\sup_{x \in \mathcal T_d} \big\lVert g(x\mid \mathcal F_{i-1})-g(x\mid \mathcal F_{i-1,\{i-k\}})\big\rVert_s.$$
We assume that for some constant $\xi > 0$ and positive integer $m$,
$$\sum_{k\geq m} \theta_{k,s} =O(m^{-\xi}).$$
\end{assumption}
Roughly speaking, the functional dependence measure $\theta_{k,s}$ defined in Assumption \ref{asm_dep} quantifies the dependence of $X_i$ on $v_{i-k}$ by measuring the distance between $g(x\mid\mathcal F_{i-1})$ and its coupled version $g(x\mid\mathcal F_{i-1,\{i-k\}})$. The physical dependence measure accounts for any measurable function of cumulative i.i.d. noises, which makes it applicable to a wide range of linear and nonlinear time series processes. To elaborate such conditional dependency structure, we shall consider a simple moving average (MA) example as follows.
\begin{example}
Suppose that the covariates $(X_1,X_2,\ldots,X_n)$ is a linear process, which takes the form
$$X_{i} = \sum_{l\geq 0}A_{l}v_{i-l},$$
where $v_i \in \mathbb{R}^d$ are i.i.d random vectors with continuous density function $f(\cdot): \mathbb R^d\mapsto\mathbb R$, and $A_l\in\mathbb R^{d\times d}$ are coefficient matrices, $l\ge0$. Without loss of generality, we assume that $A_0$ is an identity matrix, and then, $\theta_{k,s}$ in this case can be written into
$$\theta_{k,s} = \sup_{x \in \mathcal T_d} \Big\| f\Big(x-\sum_{l\ge 1}A_l v_{i-l}\Big)-f\Big(x-\sum_{l\ge 1}A_l v_{i-l}+A_k(v_{i-k}-v'_{i-k})\Big)\Big\|_s,$$
which, by the continuity of the density function $f(\cdot)$, can be bounded as follows
\begin{align*}
\theta_{k,s}\lesssim \max_{1\le j\le d}|A_{k,j,\cdot}|_2\|v_{i-k}-v'_{i-k}\|_s,
\end{align*}
where $A_{k,j,\cdot}$ is the $j$-th row of matrix $A_k$. In Assumption \ref{asm_dep}(ii), this condition actually indicates an algebraic decay rate of the temporal dependence of $X_i$. Namely, for any $m\geq 1,$ there exists some constant $\xi^*>0,$ such that
\begin{equation*}
\max_{1\le j\le d}\sum_{k\geq m} |A_{k,j,\cdot}|_2 \lesssim m^{-\xi^*}.
\end{equation*}
\end{example}
\begin{assumption}[Bounds and smoothness]
\label{asm_smooth}
For some constants $c_g,c_g'>0$, assume that
$$c_g\leq \inf_{x \in \mathcal T_d} g(x)\leq \sup_{x \in \mathcal T_d} g(x) \leq c_g',$$
and
$$\sup_{x \in \mathcal T_d}\max_{1\leq j\leq d} \big|\partial g(x)/\partial x_j\big| <\infty.$$
\end{assumption}
Assumption \ref{asm_smooth} imposes conditions on the boundary and smoothness of the density function $g(x)$ on the support $\mathcal T_d$. In the literature, conditions similar to Assumption \ref{asm_smooth} have been commonly used; see, for instance, \textcite{aneiros-perez_local_2008,kim_simultaneous_2021}.
\subsection{Gaussian Approximation}\label{subsec_GA}
This subsection is devoted to our main results, the asymptotic distribution of the proposed statistic derived by Gaussian approximation. To explicitly illustrate our theory on the simultaneous inference of trend $\mu(\cdot)$, we shall first assume that $\boldsymbol \beta$ and $\sigma(\cdot)$ in model (\ref{eq_model}) are both known. In particular, we define
\begin{equation}
\label{eq_mu_hat_sol}
\hat\mu(x) =\sum_{i=1}^nw_h(x,X_i)\big(Y_i-Z^{\top}_{i}\boldsymbol \beta\big),
\end{equation}
and consider the statistic $\sup_{x\in\mathcal T_d}\big|\hat\mu(x)- \mu(x)\big|/\sigma(x)$. By Slutsky's theorem, the asymptotic distribution of the statistic still holds when replacing the true $\boldsymbol \beta$ and $\sigma(\cdot)$ therein by their consistent estimators $\hat\boldsymbol \beta$ and $\hat\sigma(\cdot)$, respectively. Hence, in this subsection, we shall first proceed under the oracle setting, and we refer to the implementation with theoretical guarantees to Section \ref{sec_est}.
As mentioned in Section \ref{sec_SCR}, we aim to provide the limiting distribution for the statistic $\sup_{x\in\mathcal T_d}\big|\hat\mu(x)- \mu(x)\big|/\sigma(x)$ by applying the high-dimensional Gaussian approximation. Specifically, we shall generate a centered Gaussian random field $\{\mathcal Z_t\}_{t\in\mathcal T_d}$ and utilize the distribution of $\sup_{t\in\mathcal T_d}|\mathcal Z_t|$ for the approximation. We denote the covariance matrix of $\sqrt{h^dn}\big|\hat\mu(x)- \mu(x)\big|/\sigma(x)$ conditioned on the covariates $X_i$, $1\le i\le n$, by $Q=(Q_{t,s})_{t,s\in\mathcal T_d}$, where $Q_{t,s}$ takes the form
\begin{equation}
\label{eq_cov_Z}
Q_{t,s}= h^dn\sum_{k=-\infty}^{\infty}\sum_{i=1\vee (1-k)}^{n\wedge (n-k)}c_{t,s,i,k}w_h(x_t,X_i)w_h(x_s,X_{i+k})\gamma(k),
\end{equation}
where $c_{t,s,i,k} = \sigma(x_t)^{-1}\sigma(x_s)^{-1}\sigma(X_i)\sigma(X_{i+k})$. Further, in practical usage, to evaluate the covariance matrix $Q$, we would need to estimate the autocovariance $\gamma(k)$, for each $k\ge0$, while this is unrealistic when $k$ goes to infinity. Therefore, to establish a consistent estimator for $Q$, we define a truncated version of $Q$ as $Q^{(L)}=(Q_{t,s}^{(L)})_{t,s\in\mathcal T_d}$, for some large positive integer $L$, where
$Q_{t,s}^{(L)}$ is defined as
\begin{equation}
\label{eq_cov_Z_truncate}
Q^{(L)}_{t,s}= h^dn\sum_{k=1-L}^{L-1}\sum_{i=1\vee (1-k)}^{n\wedge (n-k)}c_{t,s,i,k}w_h(x_t,X_i)w_h(x_s,X_{i+k})\gamma(k).
\end{equation}
When applied to real data, the autocovariance $\gamma(k)$ can be estimated by the notable existing methods and we defer the details of this estimation to Section \ref{sec_est}. Now we let $\{\mathcal Z_t\}_{t\in\mathcal T_d}$ be a Gaussian random field with mean zero and conditional covariance matrix $Q^{(L)}$ given the covariates $X_1,X_2,\ldots,X_n$. The first main theorem is stated as follows.
\begin{theorem}[Gaussian approximation]
\label{thm1_GA}
Suppose that Assumptions \ref{asm_sigma}-\ref{asm_smooth} hold. Then, for $$\Delta_1= (h^dn)^{-q/2}n\log^q(n)+1/n, \quad \Delta_2=(h^dn)^{-1/6} \log^{7/6} (n) + (n^{2/q}/(h^dn))^{1/3} \log(n),$$
$$\Delta_3=h\log(n)+h^{2+d/2}n^{1/2
}\sqrt{\log(n)}+n^{-\zeta}\log(n), \quad \Delta_4=L^{-\zeta/3}\log^{2/3}(n),$$
we have
$$\sup_{u\in\mathbb R}\Big|\mathbb P\big(\sup_{x\in\mathcal T_d}\sqrt{h^dn}\big|\hat\mu(x)- \mu(x)\big|/\sigma(x)<u\big) - \mathbb P\big(\sup_{t\in\mathcal T_d}|\mathcal Z_t|<u\big) \Big| \lesssim \Delta_1+\Delta_2+\Delta_3 +\Delta_4,$$
where the constants in $\lesssim$ are independent of $n$ and $h$.
If in addition,
\begin{align}
\label{thm1_o1}
&(h^d)^{2-q}n^{2-q}\log^{3q}(n)\rightarrow0, \quad h\log(n)\rightarrow0, \quad h^{4+d}n\log(n) \rightarrow0,\nonumber \\
& \qquad n^{-\zeta}\log(n) \rightarrow 0 \quad \text{and} \quad L^{-\zeta}\log^2(n) \rightarrow0,
\end{align}
then we have
$$\sup_{u\in\mathbb R}\Big|\mathbb P\big(\sup_{x\in\mathcal T_d} \sqrt{h^dn}\big|\hat\mu(x)- \mu(x)\big|/\sigma(x)<u\big) - \mathbb P\big(\sup_{t\in\mathcal T_d}|\mathcal Z_t|<u\big) \Big| \rightarrow 0.$$
\end{theorem}
\begin{remark}[Comments on the convergence rate]
\label{rmk_thm1}
In Theorem \ref{thm1_GA}, condition (\ref{thm1_o1}) imposes conditions on the strength of dependency and the bandwidth parameter $h$. The first part $(h^d)^{2-q}n^{2-q}\log^{3q}(n)\rightarrow0$ requires that $h$ should not be too small, while the second and third parts, $h\log(n)\rightarrow0$ and $h^{4+d}n\log(n)$, aim to keep $h$ from being too large. Therefore, one can apply cross-validation to choose the bandwidth $h$, such as the order of $n^{-1/3}$ or $n^{-1/5}$. The last two parts $n^{-\zeta}\log(n) \rightarrow 0$ and $L^{-\zeta}\log^2(n) \rightarrow0$ suggest that the dependency strength of $\epsilon_i$ cannot be too strong, and the integer $L$ in the conditional covariance matrix $Q^{(L)}$ should be large to approximate the true covariance matrix of the proposed statistic, respectively. For a general moving average model MA$(\infty)$ defined in (\ref{eq_epsilon}), we could choose $L$ to be some large positive integer such as $L=\sqrt{n}$. For some particular MA$(p)$ model, one could simply let $L=p$. We refer to Proposition \ref{prop_longrun} and the corresponding discussion for more details.
\end{remark}
As shown by Theorem \ref{thm1_GA}, our approach differs from the one based on the Gumbel convergence in extreme value theory adopted by \textcite{zhao_kernel_2006,LW:2010}. We extend the high-dimensional Gaussian approximation in \textcite{chernozhukov_central_2017} to dependent processes and provide a Gaussian limit distribution for $\sup_{x\in\mathcal T_d}|\hat\mu(x)-\mu(x)|/\sigma(x)$. Our method does not need the density estimate of covariate $X_i$ to build the test statistic, and thus offers easier implementation (see details in Section \ref{sec_est}) than the competing ones.
\section{Implementation of SCR}\label{sec_est}
In the previous sections, we assumed that $\boldsymbol \beta$, $\sigma(\cdot)$ in model (\ref{eq_model}) and the conditional long-run covariance matrix of $\sup_{t\in\mathcal T_d}\big|\hat\mu(x)- \mu(x)\big|/\sigma(x)$ are all known, which, however, is not realistic in practice. Hence, we shall introduce the estimators $\hat\boldsymbol \beta$, $\hat\sigma(\cdot)$ and $\hat Q^{(L)}$ in this subsection, followed by the consistency results and the limit distribution of the statistic built on these estimates. The detailed instructions on how to construct SCR in practice are also provided in Section \ref{subsec_scr_construct}.
\subsection{Consistent Estimators}
The literature on estimating the parametric components in model (\ref{eq_model}) has a long history. \textcite{Robinson} was among the first contributing to this problem, where he derived a least square estimator of $\boldsymbol \beta$ based on a Nadaraya-Waston kernel estimator of $\mu$ and provided the $\sqrt{n}$-consistency result of the estimate. In this study, we shall consider the same estimator for $\boldsymbol \beta$, and we will show that the $\sqrt{n}$-consistent rate still holds under our setting. Define the Robinson's estimator $\hat\boldsymbol \beta$ as follows,
\begin{equation}
\label{eq_betahat}
\hat\boldsymbol \beta=\Big(\sum_{i=1}^n(Z_i-\tilde Z_i)(Z_i-\tilde Z_i)^{\top}\mathbf{1}_{X_i\in\mathcal T_d}\Big)^{-1}\Big(\sum_{i=1}^n(Y_i-\tilde Y_i)(Z_i-\tilde Z_i)^{\top}\mathbf{1}_{X_i\in\mathcal T_d}\Big),
\end{equation}
where $\tilde Y_i$ and $\tilde Z_i$ are the kernel estimators of $Y_i$ and $Z_i$, respectively, which are
\begin{equation}
\label{eq_tilde_YZ}
\tilde Y_i=\sum_{t=1}^nw_h(X_i,X_t)Y_t, \quad \text{and} \quad \tilde Z_i=\sum_{t=1}^nw_h(X_i,X_t)Z_t.
\end{equation}
To establish the consistency result of the Robinson's estimator $\hat\boldsymbol \beta$, we shall introduce some regularity conditions on the covariates $Z_i$ of the parametric part in (\ref{eq_model}). Specifically, following \textcite{hardle_asymptotic_1989}, we consider $Z_i$ as random design points and related to $X_i$ in the following way.
\begin{assumption}[Random design]\ \\
\label{asm_iden}
(i) Assume that there exists some bounded function $h(\cdot)$: $\mathbb{R}^d \rightarrow \mathbb{R}^l$, such that
$$Z_{i} = h(X_i)+u_i,$$
where $h(\cdot)$ is Lipschitz continuous on $\mathcal T_d$ and $u_i$ are centered random vectors in $\mathbb{R}^l$ independent of $X_i$. The minimum eigenvalue of the covariance matrix $$\Sigma_u = \mathbb E(u_iu_i^{\top})$$
is lower bounded, that is $\lambda_{\text{min}}(\Sigma_u) \ge c$, for some constant $c>0$. \\
(ii) Assume that for each $1\le i\le n$ and $1\le j\le l$, $0<\|u_{ij}\|_p<\infty$, for some constant $p\ge4$, where $u_{ij}$ is the $j$-th element of $u_i$. \\
(iii) We suppose that for some measurable function $f$,
$$u_i=f(\ldots,e_{i-1},e_i), \quad \text{and} \quad u_{i,\{i-k\}}=f(\ldots,e_{i-k-1},e_{i-k}',e_{i-k+1},\ldots,e_i),$$
where $e_i$ are i.i.d. random vectors in $\mathbb R^l$ independent of $X_i$ and $\epsilon_i$, and $e_i'$ is an i.i.d. copy of $e_i$. We assume that for $p\ge 4$,
$$\sum_{k\geq 0} \Big\| \max_{1\le j\le l}\big|u_{ij}- u_{ij,\{i-k\}}\big|\Big\|_p <\infty,$$
where $u_{ij,\{k\}}$ is the $j$-th element of $u_{i,\{k\}}$.
\end{assumption}
Assumption \ref{asm_iden} can be considered as a generalization of the conditions with weaker requirements compared to the ones introduced in \textcite{hardle_asymptotic_1989}, \textcite{gao_convergence_1995} and \textcite{Sun}, where they assumed $u_i$ to be i.i.d. random vectors, while in this paper, we allow $u_i$ to be weakly dependent over $i$. In fact, Assumption \ref{asm_iden}(iii) implies the short-range dependence (SRD) of $u_i$. The following proposition asserts that, under such dependence conditions, we can consistently estimate $\boldsymbol \beta$ in the time series regression model (\ref{eq_model}) by the Robinson's estimator at the rate of $\sqrt{n}$.
\begin{proposition}[Consistency of $\hat\boldsymbol \beta$]
\label{prop1}
Under Assumptions \ref{asm_sigma}-\ref{asm_iden}, if $h^4n\rightarrow0$, we have
$$|\hat\boldsymbol \beta-\boldsymbol \beta|_{\infty}=O_{\mathbb P}\big\{1/\sqrt{n}\big\}.$$
\end{proposition}
Concerning the conditional volatility function $\sigma(X_i)$ in (\ref{eq_model}), for each $1\le i\le n$, we define the local constant estimator
\begin{equation}
\label{eq_sigmax_est}
\hat\sigma^2(X_i)=\sum_{t=1}^nw_h(X_i,X_t)(Y_t-Z_t^{\top}\hat\boldsymbol \beta - \hat\mu^*(X_i))^2.
\end{equation}
Note that the bandwidth parameter $h$ in (\ref{eq_sigmax_est}) can be selected differently from the one used in the estimator for $\mu(\cdot)$. Similar methods have been applied in the existing studies; see, for example,
\textcite{fan_local_1996,zhao_kernel_2006}. In this paper, we adopt the same bandwidth parameters in both $\hat\mu^*(\cdot)$ and $\hat\sigma(\cdot)$ for brevity. The consistency result of $\hat\sigma(\cdot)$ is stated as follows.
\begin{proposition}[Consistency of $\hat\sigma$]
\label{prop2}
Assume that the conditions in Theorem \ref{thm1_GA} hold. Then,
we have
$$\sup_{x\in\mathcal T_d}|\hat\sigma^2(x)-\sigma^2(x)|=O_{\mathbb P}\Big\{h+\frac{1}{n}+\sqrt{\frac{\log(n)}{h^dn}}\Big\}.$$
\end{proposition}
Recall that in Theorem \ref{thm1_GA}, $\mathcal Z_t$, $t\in\mathcal T_d$, is a Gaussian random field with mean zero and conditional covariance matrix $Q^{(L)}$ as defined in (\ref{eq_cov_Z}). In practice, one shall estimate $Q^{(L)}$ by $\hat Q^{(L)}=(\hat Q_{j,j'}^{(L)})_{1\le j,j'\le n}$, that is
\begin{equation}
\label{eq_Q_hat}
\hat Q_{j,j'}^{(L)}=h^dn\sum_{k=1-L}^{L-1}\sum_{i=1\vee(1-k)}^{n\wedge (n-k)}\hat c_{j,j',i,k}w_h(X_j,X_i)w_h(X_{j'},X_{i+k})\hat\gamma(k),
\end{equation}
where $\hat c_{j,j',i,k} = \hat\sigma^{-1}(X_j)\hat\sigma^{-1}(X_{j'})\hat\sigma(X_i)\hat\sigma(X_{i+k})$, and $\hat\gamma(k)$ is a consistent estimate of $\gamma(k)$ which takes the form
\begin{equation*}
\hat\gamma(k)=\sum_{i=1}^{n-k}\hat\epsilon_i\hat\epsilon_{i+k}/n,
\end{equation*}
with $\hat\epsilon_i=(Y_i-Z_i^{\top}\hat\boldsymbol \beta-\hat\mu^*(X_i))/\hat\sigma(X_i)$. The precision of the estimated autocovariance $\hat\gamma(k)$ has been well investigated in the literature; see, for example, \textcite{shumway}. Also, note that we adopt the estimator $\hat Q_{j,j'}^{(L)}$ which includes the term $\hat c_{j,j',i,k}$ as the estimator of $c_{j,j',i,k}$ appeared in the true $Q_{t,s}^{(L)}$ defined in (\ref{eq_cov_Z}). Since $\hat\sigma(\cdot)$ is a consistent estimator of $\sigma(\cdot)$ by Proposition \ref{prop2} and $\sigma(\cdot)$ is bounded on the support $\mathcal T_d$ according to Assumption \ref{asm_sigma}, it can be shown that $\hat Q_{j,j'}^{(L)}$ is a consistent estimator of $Q_{j,j'}^{(L)}$. We state this result and provide the exact consistency rate as follows.
\begin{proposition}[Precision of long-run covariance matrix]
\label{prop_longrun}
Suppose that Assumptions \ref{asm_sigma}-\ref{asm_iden} hold. Then, for some large positive integer $L\rightarrow\infty$ satisfying $hL^2/n\rightarrow0$, $L^3/n^2\rightarrow0$ and $L^4\log(n)/(h^dn^3)\rightarrow0$, we have
$$\max_{1\le j,j'\le n}\big|\hat Q_{j,j'}^{(L)} - Q_{j,j'}\big| = O_{\mathbb P}\Big\{\frac{L^2}{n}\Big(h+\sqrt{\log(n)/h^dn} + 1/\sqrt{n}\Big) + L/n + L^{-\zeta}\Big\},$$
where the constant $\zeta>1$ is defined in Assumption \ref{asm_dep_epsilon}.
\end{proposition}
In real-data applications, one can simply take $L=\sqrt{n}$ to obtain a consistent estimator $\hat Q^{(L)} =(\hat Q_{j_1,j_2}^{(L)})_{1\le j_1,j_2\le n}$ to estimate the covariance matrix of the statistic $\sqrt{h^dn}|\hat\mu^*(x)-\mu(x)|/\hat\sigma(x)$ conditioned on the covariates $X_i$, $1\le i\le n$. In next section, we shall provide the Gaussian approximation for the test statistic $\sup_{x\in\mathcal T_d}|\hat\mu^*(x)-\mu(x)|/\hat\sigma(x)$ with the estimated covariance matrix applied.
\subsection{Gaussian Approximation for Implementation}
Provided the consistent estimated long-run covariance matrix $\hat Q^{(L)}$ in Proposition \ref{prop_longrun}, we now introduce a sequence of Gaussian random variables $\hat\mathcal Z_j$, $1\le j\le n$, with mean zero and the conditional covariance matrix $\hat Q^{(L)}$, i.e., $\hat Q_{j,j'}^{(L)} = \mathbb E(\hat\mathcal Z_j\hat\mathcal Z_{j'})$, for each $1\le j,j'\le n$. We shall show that the Gaussian approximation result in Theorem \ref{thm1_GA} still holds when we replace $Q_{t,s}^{(L)}$ therein by its estimator $\hat Q_{j,j'}^{(L)}$. In addition, Propositions \ref{prop1} and \ref{prop2} ensure that even when $\boldsymbol \beta$ and $\sigma(x)$ in model (\ref{eq_model}) are both unknown, one shall still achieve the similar asymptotic distribution in Theorem \ref{thm1_GA} for the test statistic $\sup_{x\in\mathcal T_d}\big|\hat\mu^*(x)- \mu(x)\big|/\hat\sigma(x)$ by employing the proposed estimators $\hat\boldsymbol \beta$ and $\hat\sigma(x)$. As a result, combining Propositions \ref{prop1}--\ref{prop_longrun} facilitates the construction of SCR for practical use. For the theoretical guarantee, we show the Gaussian approximation result with $\hat\boldsymbol \beta$, $\hat\sigma(x)$ and $\hat Q^{(L)}$ plugged in as follows.
\begin{theorem}[Gaussian approximation for implementation]
\label{thm2_GA}
Under the conditions in Theorem \ref{thm1_GA} and Propositions \ref{prop1}--\ref{prop_longrun}, for the same $\Delta_1$, $\Delta_2$ and $\Delta_3$ defined in Theorem \ref{thm1_GA} and
$$\Delta_4'=\Big(hL^2/n+L^2\sqrt{\log(n)/(h^dn^3)} + L/n + L^{-\zeta}\Big)^{1/3}\log^{2/3}(n),$$
we have
$$\sup_{u\in\mathbb R}\Big|\mathbb P\big(\sup_{x\in\mathcal T_d}\sqrt{h^dn}\big|\hat\mu^*(x)- \mu(x)\big|/\hat\sigma(x)<u\big) - \mathbb P\big(\max_{1\le j\le n}|\hat\mathcal Z_j|<u\big) \Big| \lesssim \Delta_1+\Delta_2+\Delta_3+\Delta_4',$$
where the constants in $\lesssim$ are independent of $n$ and $h$.
If in addition, the conditions in expression (\ref{thm1_o1}) hold, and
\begin{align}
\label{thm2_o1}
hL^2n^{-1}\log^2(n)\rightarrow0, \quad L^4h^{-2d}n^{-6}\log^5(n)\rightarrow0, \quad Ln^{-1}\log^2(n)\rightarrow0,
\end{align}
then we have
$$\sup_{u\in\mathbb R}\Big|\mathbb P\big(\sup_{x\in\mathcal T_d}\sqrt{h^dn}\big|\hat\mu^*(x)- \mu(x)\big|/\hat\sigma(x)<u\big) - \mathbb P\big(\max_{1\le j\le n}|\hat\mathcal Z_j|<u\big) \Big| \rightarrow 0.$$
\end{theorem}
\begin{remark}[Comparison of convergence rates in Theorems \ref{thm1_GA} and \ref{thm2_GA}]
Note that the convergence rate in Theorem \ref{thm2_GA} differs from the one in Theorem \ref{thm1_GA} only in the last part $\Delta_4'$ (appeared as $\Delta_4$ in Theorem \ref{thm1_GA}). This difference results from the replacement of the covariance matrix $Q^{(L)}$ by its consistent estimator $\hat Q^{(L)}$, which leads to an additional approximation error due to the Gaussian comparison as shown in Lemma \ref{lemma_comparison} in the Supplementary Material. The other approximation errors introduced by the estimators $\hat\boldsymbol \beta$ and $\hat\sigma(x)$ can be shown to be negligible compared to $\Delta_1$--$\Delta_3$ and $\Delta_4'$. In particular, when we let $L=\sqrt{n}$, the convergence rates in Theorems \ref{thm1_GA} and \ref{thm2_GA}
are of the same magnitude up to a logarithmic factor.
We refer to a rigorous proof of these convergence rates to the Supplementary Material. Based on Theorem \ref{thm2_GA}, we can construct the $(1-\alpha)$-th percentile SCR defined in (\ref{eq_ln_rn}) by using the $(1-\alpha)$-th quantile of $\max_{1\le j\le n}|\hat\mathcal Z_j|/\sqrt{h^dn}$.
\end{remark}
\begin{remark}[The necessity of partially linear structure]
Generally, the linear part in (\ref{eq_model}) is derived out of its underlying theories, instead of being added to the model in an ad-hoc fashion; see, for example, the demand equation in \textcite{kim_simultaneous_2021} and the forward premium regression in (\ref{fama}). Secondly, the purely nonparametric approach
suffers from the curse of dimensionality with a large number of model predictors, while our partially linear model is obviously more efficient in dealing with the issue (\cite{wang_conditional_2016}) than the purely nonparametric counterpart. Thirdly, it is nontrivial to extend the asymptotic theory of the purely nonparametric model to the partially linear case, because $Z_i$ in (\ref{eq_model}) is a function of $X_i$ plus a random noise $u_i$ (i.e. Assumption \ref{asm_iden}). As a consequence, establishing the consistency results for $\hat\boldsymbol \beta$ and $\hat\sigma$ is different than in the purely nonparametric case.
\end{remark}
\subsection{SCR Construction}\label{subsec_scr_construct}
Recall the upper and lower bounds of the simultaneous confidence region, i.e., $\hat l_n(x)$ and $\hat r_n(x)$ defined in (\ref{eq_ln_rn}). For the convenience of implementation, we provide the detailed steps of our proposed method for the construction of SCR as follows. Suppose that we have observed a sample $(Z_i,X_i,Y_i)$, $1\le i\le n$. Given the significance level $\alpha\in(0,1)$, we aim to construct the $(1-\alpha)$-th percentile SCR for $\mu(\cdot)$ in model (\ref{eq_model}).
\begin{itemize}
\item \textbf{Step 1.} Calculate the weight function $w_h(x,X_t)$ at each covariate $X_i$:
$$w_h(X_i,X_t) = \frac{K_h(X_i-X_t)}{\sum_{t=1}^nK_h(X_i-X_t)},$$
where $K_h(x-X_t)=K\big((x-X_t)/h\big)/h$. One can choose the kernel function $K(\cdot)$ to be any commonly used kernels, such as the Gaussian kernel $K(u)=(2\pi)^{-1/2}e^{-u^2/2}$ and the Epanechnikov kernel $K(u)=3/4(1-u^2)\mathbf{1}_{|u|\le1}$.
\item \textbf{Step 2.} Estimate the coefficients $\boldsymbol \beta$ of the linear part:
$$\hat\boldsymbol \beta= \Big(\sum_{i=1}^n(Z_i-\tilde Z_i)^2\Big)^{-1}\Big(\sum_{i=1}^n(Y_i-\tilde Y_i)(Z_i-\tilde Z_i)\Big),$$
where $\tilde Y_i=\sum_{t=1}^nw_h(X_i,X_{t})Y_{t}$, and $\tilde Z_i=\sum_{t=1}^nw_h(X_i,X_{t})Z_{t}$.
\item \textbf{Step 3.} Calculate the point estimator of the nonparametric function $\mu(\cdot)$ at each covariate $X_i$:
$$\hat\mu^*(X_i) = \sum_{t=1}^nw_h(X_i,X_t)(Y_t-Z_t^{\top}\hat\boldsymbol \beta).$$
\item \textbf{Step 4.} Estimate the conditional variance/volatility function $\sigma(\cdot)$ at each covariate $X_i$:
$$\hat\sigma^2(X_i)=\sum_{t=1}^nw_h(X_i,X_t)(Y_t-Z_t^{\top}\hat\boldsymbol \beta - \hat\mu^*(X_i))^2.$$
\item \textbf{Step 5.} Let $N=1000$ (or some other large positive number). Apply the multiplier bootstrap to generate a sequence of centered Gaussian random variables $\hat\mathcal Z_j$, $1\le j\le N$, with covariance matrix $\hat Q^{(L)}=(\hat Q_{j,j'}^{(L)})_{1\le j,j'\le N}$, where
$$\hat Q_{j,j'}^{(L)}=h^dn\sum_{j-j'=-L}^L\sum_{i=1}^nc_{j,j',i}\hat\gamma(|j-j'|),$$
with $c_{j,j',i}=w_h(X_j,X_i)w_h(X_{j'},X_{i+j-j'})$, $\hat\gamma(k)=\sum_{i=1}^{n-k}\hat\epsilon_i\hat\epsilon_{i+k}/n$, and $\hat\epsilon_i=(Y_i-Z_i^{\top}\hat\boldsymbol \beta-\hat\mu^*(X_i))/\hat\sigma(X_i)$.
Regarding the choice of $L$, see Proposition \ref{prop_longrun}, and we can choose $L$ to be some large positive integer such as $L=\sqrt{n}$. In some special case such as MA(2), one can simply choose $L$ to be the maximum lag, i.e., $L=2$.
\item \textbf{Step 6.} Construct the $(1-\alpha)$-th percentile SCR for $\mu(\cdot)$:
\begin{equation}
\label{eq_SCR_final}
\hat\mu^*(X_i)-\hat q_{\alpha}\hat\sigma(X_i) \le \mu(X_i)\le \hat\mu^*(X_i) + \hat q_{\alpha}\hat\sigma(X_i), \quad 1\le i\le n,
\end{equation}
where $\hat q_{\alpha}$ is the $(1-\alpha)$-th empirical quantile of $\max_{1\leq j\leq n}|\hat\mathcal Z_j|/\sqrt{h^dn}$.
\end{itemize}
\begin{remark}[Dynamic width of SCR]
Our proposed SCR in (\ref{eq_SCR_final}) takes into account the conditional heteroskedasticity in model (\ref{eq_model}). The conditional heteroskedasticity allows the conditional variance of the error to be time-varying, and also accounts for the dependence between the covariate $X_i$ and the model error $\sigma(X_i)\epsilon_i$ (\cite{zhao_inference_2013}). Additionally, the conditional heteroskedasticity allows the width of the confidence band to vary over the covariate term, and thus, the relative width can indicate the different amounts of information used for different regions of the band. For example, the lack of information at each tail of the curve, due to fewer observations, is implied by the relatively wide confidence intervals around the tails. This dynamic-width SCR is adopted by \textcite{kim_simultaneous_2021} for the i.i.d. data. We extend it to the time series setting, so that one can accommodates macroeconomic and/or financial time series models such as the forward premium regression; see Section \ref{sec_app}.
\end{remark}
\section{Simulation Studies}\label{sec_simul}
We present the simulation study with four different time series models in this section to illustrate the performance of our proposed SCR. Consider the following data-generating process:
\begin{eqnarray}\label{smodel}
y_i=\beta z_i+\mu(x_i)+\sigma(x_i)\epsilon_i,~~~~~~~~i=1,\cdots,n
\end{eqnarray}
where $\beta=0.5$, $\mu(x)=0.3+0.4\,x$ and $\sigma(x)=\big(0.1+0.1\,x^2\big)^{1/2}$. We shall compute the {\it coverage probability} of the SCR for $\mu(\cdot)$ in (\ref{smodel}) to evaluate our method. We set the covariate terms $x_i$ and $z_i$ to be
\begin{eqnarray*}
x_i&=&\sum_{k=0}^{\infty}0.1^k\,\delta_{i-k}\\
z_i&=&0.2+0.4 x_i+u_i,~~~~~~~~i=1,\cdots,n,
\end{eqnarray*}
where $\delta_i$ and $u_i$ are both i.i.d. random variables following the standard normal distribution, and $\delta_i$ are independent of $u_i$. Note that this setting is in accordance with Assumption \ref{asm_iden}. For the specifications of $\epsilon_i$ in (\ref{smodel}), we consider four commonly used models: (i) i.i.d. standard normal, (ii) an autoregressive (AR) process, (iii) a moving-average (MA) process and (iv) an autoregressive moving-average (ARMA) process, which take the forms of:
\begin{eqnarray*}
&&(i)\,\,\mbox{Std. Normal}:\,\epsilon_i=\eta_i,\\
&&(ii)\,\,\mbox{AR(1)}:\,\epsilon_i=0.1\,\epsilon_{i-1}+\eta_i,\\
&&(iii)\,\,\mbox{MA(1)}:\,\epsilon_i=\eta_i+0.2\eta_{i-1},\\
&&(iv)\,\,\mbox{ARMA(1,1)}:\,\epsilon_i=0.1\,\epsilon_{i-1}+\eta_i+0.2\,\eta_{i-1},~~~~~~~~i=1,\cdots,n
\end{eqnarray*}
where $\eta_i$ are i.i.d. random variables following the standard normal distribution. For a range of bandwidths (see Table \ref{tbl_cov_prob}), we follow Steps 1--6, the procedures described in the previous section, to construct the 95\% SCR of $\mu(\cdot)$ in (\ref{smodel}). The process is repeated $M$ times for each chosen bandwidth parameter. Given these $M$ different SCRs, we count how many of those SCRs contain the true $\mu(\cdot)$ in them, which gives us the desired coverage probabilities. Here we let $M=1,000$ (i.e. number of iterations) and $n=200$ (i.e. sample size).
The simulation results under (i)--(iv) are summarized by Table \ref{tbl_cov_prob}.
\begin{table}[tbp]
\centering
\caption[] {Coverage probability of 95\% SCR for $\mu(\cdot)$}
\label{tbl_cov_prob}
\begin{tabular}{cccccccccc}
\hline\hline
bandwidth && \multicolumn{1}{c}{(i) Std. Normal} && \multicolumn{1}{c}{(ii) AR(1)} && \multicolumn{1}{c}{(iii) MA(1)} && \multicolumn{1}{c}{(iv) ARMA(1,1)} & \\ \hline
0.30 && 0.874 && 0.858 && 0.860 && 0.852 & \\
0.32 && 0.898 && 0.890 && 0.884 && 0.888 & \\
0.34 && 0.918 && 0.898 && 0.920 && 0.902 & \\
0.36 && 0.938 && 0.914 && 0.932 && 0.918 & \\
0.38 && 0.948 && 0.924 && 0.934 && 0.936 & \\
0.40 && 0.950 && 0.934 && 0.942 && 0.948 & \\
0.42 && 0.958 && 0.950 && 0.956 && 0.954 & \\
0.44 && 0.963 && 0.951 && 0.958 && 0.958 & \\
0.46 && 0.970 && 0.960 && 0.960 && 0.960 & \\
0.48 && 0.970 && 0.958 && 0.964 && 0.964 & \\
0.50 && 0.972 && 0.964 && 0.964 && 0.964 & \\ \hline
\end{tabular}
\end{table}
From Table \ref{tbl_cov_prob}, we see that the coverage probabilities under various bandwidths are reasonably close to $0.95$, which is the nominal coverage rate at a 5 percent level. For some bandwidths, we have the coverage probabilities that are almost identical to the nominal level $95\%$. However, as the bandwidth gets extreme, the coverage probability deviates from the nominal $0.95$, as we can observe from Table \ref{tbl_cov_prob}.
\section{Application: Forward Premium Regression}\label{sec_app}
One of the useful applications for (\ref{eq_model}) in time series analysis is the {\it forward premium anomaly}. Consider the celebrated forward premium regression (\cite{Fm:1984}):
\begin{equation}\label{fama}
s_{t+1}-s_t=\mu +\beta\cdot(f_{1,t}-s_{t})+u_{t+1},
\end{equation}
where $s_{t}$ is log of monthly {\it spot} exchange rate at time $t$ and $f_{1,t}$ is log of monthly {\it forward} exchange rate with one-month maturity at time $t$, respectively. Here $u_t$ is a model error. A vast amount of literature illustrates that the constant term $\mu$ in (\ref{fama}), the foreign exchange risk premium, is related to {\it fundamentals} (\cite{DH:1985,HS:1986,Hod:1989,Mark:1995,BK:2006,Alv:2009,LRV:2014,BK:2015,KKK:2022}), such as money growth, output growth, interest rates, conditional variance of money growth and so on. Given this, one can possibly model (\ref{fama}) in the following flexible manner:
\begin{equation}\label{fama2}
s_{t+1}-s_t=\mu(x_{1,t},\,x_{2,t}) + \beta\cdot(f_{1,t}-s_{t}) + \sigma(x_{1,t},\,x_{2,t})\,\epsilon_{t+1},
\end{equation}
where
\begin{equation*}
x_{1,t}=m_t-m^*_t,\,\,\,\,\,\,\,x_{2,t}=y_t-y^*_t.
\end{equation*}
Here, $m_t$ and $y_t$ are log of domestic money stock and log of domestic production, respectively, and $m^*_t$ and $y^*_t$ are their foreign counterparts. Note here that $\mu$ in (\ref{fama}) is replaced by $\mu(x_{1,t},\,x_{2,t})$, and that the conditional heteroscedasticity is introduced in (\ref{fama2}). Interestingly, Assumption \ref{asm_iden} for (\ref{fama2}) implies, quite intuitively, that the forward premium $f_{1,t}-s_t$ does depend on the fundamentals $x_{1,t}$ and $x_{2,t}$.
Although the theory of Uncovered Interest Parity (UIP) in international finance implies that $\mu(\cdot)=0$, numerous empirical studies actually show $\mu(\cdot)\not=0$. We refer to \textcite{FT:1990} and \textcite{Eg:1996} for an excellent review of the literature on this. Recently, \textcite{BDKK:2023} revisited the issue and tested the UIP hypothesis to show the effect of modeling lagged spot returns and lagged forward premiums in (\ref{fama}). Interestingly, (\ref{fama2}) is a specification of our partially linear model framework in (\ref{eq_model}). Hence, the methodology developed in this study can be readily applied to decide whether or not the UIP condition holds for (\ref{fama2}).
The exchange rate data used in this study are monthly spots and 30-day forward exchange rate data for Australian Dollar (AUD), British Pound (GBP) and U.S. Dollar (USD), where USD is the numeraire currency. That is, we employ the AUD/USD and GBP/USD currency pairs in this empirical study. For the money and the output variables, we employ monthly M1 and monthly industrial production index home (i.e. U.S.) and abroad (i.e. Australia and U.K.). Note that the usual GDP series cannot be used here because the sampling frequency is monthly in this application. End-of-month observations from December 1988 through December 2016, including the 2008 financial crisis, are used in this study.
The estimation and inference results, including the simultaneous confidence region, are reported by Figures \ref{fig_aud_mean_vol}--\ref{fig_gbp_y}. First of all, the semi-parametric Robinson estimate of the UIP coefficient $\beta$ in (\ref{fama2}) and its standard error are $\hat\beta_R=0.362$ and $s.e.(\hat\beta_R)=0.753$ for AUD/USD, and $\hat\beta_R=0.802$ and $s.e.(\hat\beta_R)=0.776$ for GBP/USD. In contrast, the traditional OLS estimates of $\beta$ in (\ref{fama}) are $\hat\beta_{OLS}=-0.029$ and $s.e.(\hat\beta_{OLS})=0.718$ for AUD/USD, and $\hat\beta_{OLS}=0.689$ and $s.e.(\hat\beta_{OLS})=0.641$ for GBP/USD. Interestingly, our semi-parametric estimates of the UIP coefficient for both AUD/USD and GBP/USD are closer to one, the value implied by UIP, with lower $t$-statistics than their OLS counterparts.
Figure \ref{fig_aud_mean_vol} shows the estimates of $\hat{\mu}(\cdot,\cdot)$ and $\sigma(\cdot,\cdot)$ (i.e. the nonlinear surfaces) for AUD/USD, while Figure \ref{fig_gbp_mean_vol} shows the mean and volatility estimates for GBP/USD. As we can see from Figures \ref{fig_aud_mean_vol} and \ref{fig_gbp_mean_vol}, the estimates of both the mean and the volatility functions for each currency pair are {\it highly nonlinear}, which implies that (\ref{fama}) may not be an appropriate parametric form for the underlying processes.
Figure \ref{fig_aud_m} represents the two-dimensional relationship between $\mu(\cdot,\cdot)$ and the relative output growth for AUD/USD when the relative money growth is fixed at a certain percentile. Similarly, Figure \ref{fig_aud_y} represents the relationship between $\mu(\cdot,\cdot)$ and the relative money growth for AUD/USD when the relative output growth is fixed at a certain percentile. Each panel in Figures \ref{fig_aud_m} and \ref{fig_aud_y} also includes the 95\% SCR of the corresponding function estimate. The horizontal line represents the null hypothesis of $\mu(\cdot,\cdot)=0$.
In order to not reject the null hypothesis of zero risk premium at a 95\% confidence level, the horizontal line fixed at zero must be {\it entirely} contained by $all$ SCRs in Figures \ref{fig_aud_m} and \ref{fig_aud_y}. Clearly, Panels (a)--(c) in Figure \ref{fig_aud_m} and Panels (a)--(c) in Figure \ref{fig_aud_y} fail to contain the horizontal line entirely. That is, all of the panels in Figures \ref{fig_aud_m} and \ref{fig_aud_y} fail to contain the horizontal line entirely. Hence, the null hypothesis gets rejected at a 5\% level, which is the key implication of Figures \ref{fig_aud_m} and \ref{fig_aud_y}. Even a more general null hypothesis of {\it constant} risk premium is also rejected clearly at a 5 percent level because the SCRs in Figure \ref{fig_aud_y} cannot entirely contain any horizontal line in them.
The main reason for the rejection is because of the tendency that the local-linear estimate of $\mu(\cdot,\cdot)$ changes $nonlinearly$ in both the output growth and the money growth. Similarly, Figures \ref{fig_gbp_m} and \ref{fig_gbp_y} show the results corresponding to the GBP/USD pair. Again, the null hypothesis gets rejected at a 5\% level because the 95\% SCRs in all panels of Figures \ref{fig_gbp_m} and \ref{fig_gbp_y}, except for Panel (b) of Figure \ref{fig_gbp_m}, fail to entirely contain the horizontal line fixed at zero.
Interestingly, these findings highlight the relative advantage of our SCR-based inference because, with the usual $p$-value of a MISE-type test statistic (\cite{HM:1993}), it is rather difficult to figure out why the null gets rejected. Under our approach, however, we can easily locate exactly where the SCR fails to contain the null, so that we can propose an appropriate function for the underlying process more easily.
Intuitively, the trend estimates in this study make sense because either a relative increase in the domestic money supply (i.e. an increase in $x_{1,t}$) or a relative decrease in the domestic production (i.e. a decrease in $x_{2,t}$) will cause the foreign currency to appreciate, yielding a decrease in the risk premium for holding the appreciating currency. Hence the risk premium, $\mu(\cdot,\cdot)$, should increase with a relative increase in the domestic production or a relative decrease in the domestic money supply. This intuition appears to be generally in line with the results in Figures \ref{fig_aud_mean_vol}--\ref{fig_gbp_y}.
\section{Concluding Remarks}
\label{conclusion}
This paper proposed a new methodology to conduct simultaneous inference of trend in a semi-parametric partially linear time series model and the number of covariate terms in the nonlinear part is allowed to be two or higher. We generalized the high-dimensional Gaussian approximation theory (\cite{chernozhukov_central_2017}) to construct the simultaneous confidence region (SCR) for the multivariate unknown trend in time series. This work can be viewed as an extension of the two-dimensional uniform confidence band (\cite{johnston_probabilities_1982,KHK:2014}) and a generalization to the {\it time series} setting compared to \textcite{kim_simultaneous_2021}. The relating asymptotic properties of the introduced methodology are investigated and the finite-sample properties are studied through a simulation experiment. The developed methodology is applied to the forward premium regression (\cite{Fm:1984}). The empirical analysis based on simultaneous inference confirms that the zero-risk-premium hypothesis for the AUD/USD and GBP/USD currency pairs is rejected at a 5$\%$ level, mainly due to the underlying nonlinear nature, as shown by Figures \ref{fig_aud_mean_vol}--\ref{fig_gbp_y}.
The current study can be extended to the case where all the model parameters are functions of random processes. For instance, the UIP coefficient in (\ref{fama2}) can be modeled as a function of fundamentals as well. Similarly, both the pricing error and the beta coefficient from a standard factor pricing model in finance can be modeled as functions of relevant state variables. The methodology developed in this study can be easily extended to handle such a framework. Additionally, one can conduct simultaneous inference of the unknown function in (\ref{eq_model}) when the covariate terms are non-stationary processes. This would significantly generalize the results in this study. Further insight can be gained by extending the current work in these and other directions, and the authors are currently working on these issues.
\begin{figure}[tbp]
\centering
\begin{tabular}{ccc}
\setlength{\itemsep}{-1.0cm} \psfig{file = figs/3D_mean_aud20230117.pdf, width = 8.0cm,
angle=0} & \psfig{file = figs/3D_vol_aud20230117.pdf, width = 8.0cm, angle=0} & \\
(a) Conditional mean & (b) Conditional volatility & \\
& &
\end{tabular}
{\small {} }
\caption{Conditional mean and volatility in the money growth and the output growth for Australian Dollar (AUD)--US Dollar(USD); The bandwidth is obtained through under-smoothing of the GCV-chosen one.}
\label{fig_aud_mean_vol}
\end{figure}
\begin{figure}[tbp]
\centering
\begin{tabular}{ccc}
\setlength{\itemsep}{-1.0cm} \psfig{file = figs/20230117aud_y_rp_q25m.pdf, width = 6.0cm,
angle=0} & \psfig{file = figs/20230117aud_y_rp_q50m.pdf, width = 6.0cm, angle=0} & \\
(a) 25th-percentile of money growth & (b) 50th-percentile of money growth & \\
\psfig{file = figs/20230117aud_y_rp_q75m.pdf, width = 6.0cm, angle=0} & \\
(c) 75th-percentile of money growth & \\
& &
\end{tabular}
{\small {} }
\caption{The $solid$ curve is the local-linear regression estimate of $\mu(\cdot)$ for AUD/USD. The $dashed$ band is the 95\% SCR of $\mu(\cdot)$ and the dotted horizontal line is $H_0:\,\mu(\cdot)=0$. The bandwidth is obtained through under-smoothing of the GCV-chosen one.}
\label{fig_aud_m}
\end{figure}
\begin{figure}[tbp]
\centering
\begin{tabular}{ccc}
\setlength{\itemsep}{-1.0cm} \psfig{file = figs/20230117aud_m_rp_q25y.pdf, width = 6.0cm,
angle=0} & \psfig{file = figs/20230117aud_m_rp_q50y.pdf, width = 6.0cm, angle=0} & \\
(a) 25th-percentile of output growth & (b) 50th-percentile of output growth & \\
\psfig{file = figs/20230117aud_m_rp_q75y.pdf, width = 6.0cm, angle=0} & \\
(c) 75th-percentile of output growth & \\
& &
\end{tabular}
{\small {} }
\caption{The $solid$ curve is the local-linear regression estimate of $\mu(\cdot)$ for AUD/USD. The $dashed$ band is the 95\% SCR of $\mu(\cdot)$ and the dotted horizontal line is $H_0:\,\mu(\cdot)=0$. The bandwidth is obtained through under-smoothing of the GCV-chosen one.}
\label{fig_aud_y}
\end{figure}
\begin{figure}[tbp]
\centering
\begin{tabular}{ccc}
\setlength{\itemsep}{-1.0cm} \psfig{file = figs/3D_mean_gbp20230117.pdf, width = 8.0cm,
angle=0} & \psfig{file = figs/3D_vol_gbp20230117.pdf, width = 8.0cm, angle=0} & \\
(a) Conditional mean & (b) Conditional volatility & \\
& &
\end{tabular}
{\small {} }
\caption{Conditional mean and volatility in the money growth and the output growth for British Pound (GBP)--US Dollar (USD); The bandwidth is obtained through under-smoothing of the GCV-chosen one.}
\label{fig_gbp_mean_vol}
\end{figure}
\begin{figure}[tbp]
\centering
\begin{tabular}{ccc}
\setlength{\itemsep}{-1.0cm} \psfig{file = figs/20230117gbp_y_rp_q25m.pdf, width = 6.0cm,
angle=0} & \psfig{file = figs/20230117gbp_y_rp_q50m.pdf, width = 6.0cm, angle=0} & \\
(a) 25th-percentile of money growth & (b) 50th-percentile of money growth & \\
\psfig{file = figs/20230117gbp_y_rp_q75m.pdf, width = 6.0cm, angle=0} & \\
(c) 75th-percentile of money growth & \\
& &
\end{tabular}
{\small {} }
\caption{The $solid$ curve is the local-linear regression estimate of $\mu(\cdot)$ for GBP/USD. The $dashed$ band is the 95\% SCR of $\mu(\cdot)$ and the dotted horizontal line is $H_0:\,\mu(\cdot)=0$. The bandwidth is obtained through under-smoothing of the GCV-chosen one.}
\label{fig_gbp_m}
\end{figure}
\begin{figure}[tbp]
\centering
\begin{tabular}{ccc}
\setlength{\itemsep}{-1.0cm} \psfig{file = figs/20230117gbp_m_rp_q25y.pdf, width = 6.0cm,
angle=0} & \psfig{file = figs/20230117gbp_m_rp_q50y.pdf, width = 6.0cm, angle=0} & \\
(a) 25th-percentile of output growth & (b) 50th-percentile of output growth & \\
\psfig{file = figs/20230117gbp_m_rp_q75y.pdf, width = 6.0cm, angle=0} & \\
(c) 75th-percentile of output growth & \\
& &
\end{tabular}
{\small {} }
\caption{The $solid$ curve is the local-linear regression estimate of $\mu(\cdot)$ for GBP/USD. The $dashed$ band is the 95\% SCR of $\mu(\cdot)$ and the dotted horizontal line is $H_0:\,\mu(\cdot)=0$. The bandwidth is obtained through under-smoothing of the GCV-chosen one.}
\label{fig_gbp_y}
\end{figure}
\clearpage
\printbibliography
\clearpage