Multiple–index Nonstationary Time Series Models:
Robust Estimation Theory and Practice

Robust model building and determination of good models for estimation and prediction is a very important part of much empirical research in economics and finance. Since real world data often display characteristics, such as nonlinearity, nonstationarity and spatiality, an important aspect of model building and determination is how to take such characteristics into account. Meanwhile, the increasing availability of economic and financial data in recent years has been accompanied by increasing interest in both theoretical research and empirical data analysis in diverse fields of sciences as well as more broadly within social and medical sciences.

One challenge is how to select the most relevant variables when there are many variables available from the data. One initial step is to explore a simple method associated with principal component analysis. Different from dealing with high-dimensional data issues in medical sciences for example, we do need to pay attention to many econometric issues, such as multicollinearity when aggregating real datasets for model building. Meanwhile, one additional complexity is that many original data appear to be highly dependent and even nonstationary. Simply using a transformed version of the original data may completely ignore some key features involved in the original data, such as trending behaviours, which can be a key interest in modelling time series data in climatology, energy, epidemiology, health and macroeconomics. More importantly, transformed versions of two different datasets may have similar features, but the original datasets can have very different structures. Another challenge is the development of new models and estimation methods that are not only robust in theory, but also capable of offering user-friendly tools for both the computational and implementational procedures of empirical data analysis in practice.

There are several levels of robustness involved in the discussion of this paper. First, the proposed model building and estimation methods are robust to data with outliers, spatiality and temporality, different types of nonstationarity. Second, the class of models we propose are applicable to several types of data, such as stationary time series data with deterministic trends, highly dependent and nonstationary time series data. Third, the proposed robust estimation procedures are iterative and computationally tractable. Fourth, the proposed models and robust estimation methods as well as the associated computational algorithms are user-friendly and easily implementable for empirical researchers and practitioners. While the relevant literature accounts for one or more of the model building and estimation issues we will address in this paper, the novelty of this paper is addressing them all simultaneously to produce new models and feasible estimation methods with computational tractability as well as empirical relevance and applicability.

We now propose a cointegrating relationship among time series sequence $\{y_{t},x_{t},z_{t},\tau_{t},1\leq t\leq n\}$, where $y_{t}$ is a response variable, and we first specify the respective structures of $x_{t}$ and $z_{t}$ by

where $\tau_{t}=\frac{t}{n}$, as defined later, $w_{t}$ is a linear vector process, and $h(\cdot)$ is a known vector of functions driven by the time trend $\tau_{t}$ and stationary $v_{t}$.

In the formulation of (1.1) $x_{t}$ is an $I(1)$ vector, and it is well known that $x_{t}$ may capture ‘large’ trending behaviours. Meanwhile, $z_{t}$ is a function in ($\tau_{t},v_{t}$) that can capture ‘small’ fluctuations. The general form of $z_{t}$ accommodates possible nonlinearity, as it is quite common in nonstationary time series settings. Such general structure covers some important special cases, such as (1) $z_{t}=h(\tau_{t})+v_{t}$, a stationary vector superimposed with a vector of time trends; (2) $z_{t}=\sigma(\tau_{t})^{\top}v_{t}$, a combination of $v_{t}$ with varying coefficients ($\sigma(\cdot)$ may be a matrix); and (3) $z_{t}=h(\tau_{t})+\sigma(\tau_{t})^{\top}v_{t}$, which accommodates both deterministic and stochastic trending behaviours in the respective mean and volatility components.

We therefore propose a parametric multiple-index model of the form:

for $t=1,2,\cdots,n,$ where $\theta_{0}\in\mathbb{R}^{d_{1}}$ and $\beta_{0}\in\mathbb{R}^{d_{2}}$ are unknown parameters, both $g_{1}(\cdot)$ and $g_{2}(\cdot)$ are parametrically known functions, and $e_{t}$ is an error term. Different from existing settings, model (1.2) takes into account not only the ‘large’ stochastic trending component driven by $x_{t}^{\top}\theta_{0}$, but also the ‘small’ slowly varying component $z_{t}^{\top}\beta_{0}$ mainly represented by the deterministic time trend $\tau_{t}$. In doing so, one may consider $g_{2}(z_{t}^{\top}\beta_{0})$ as a ‘conditional mean’ decomposition of the ‘original’ error term: $\xi_{t}=g_{2}(z_{t}^{\top}\beta_{0})+e_{t}$ to eliminate any potential endogeneity involved in $e_{t}$. As a consequence, one may assume that $e_{t}$ is uncorrelated with $(x_{t}^{\top}\theta_{0},z_{t}^{\top}\beta_{0})$, although $\xi_{t}$ may be correlated with $x_{t}^{\top}\theta_{0}$. Our focus of this paper is that $x_{t}^{\top}\theta_{0}$ is nonstationary (rather than cointegrated), and the primary goal is to estimate these unknown parameters by M-estimation to take the advantage of its robustness for modelling nonlinear functional forms of nonstationary time series data.

To the best of our knowledge, model (1.2) has not been proposed and studied in the relevant literature. For the univariate regressor setting, Chang et al., (2001) investigate a parametric additive model to accommodate a time trend, and both univariate stationary and unit root regressors, and the authors estimate the unknown parameters by a nonlinear least-square method. Chang and Park, (2003) study a parametric model that specifies the regression function as a mixture of linear function and nonlinear function of unit root process that appears as a single-index structure. Meanwhile, Dong et al., (2016) consider the estimation for a partially linear single-index model of I(1) process vector. In addition, there are a couple of other models that may be considered as nonparametric counterparts of model (1.2). Dong and Linton, (2018) and Dong et al., (2021) extend Chang et al., (2001), respectively, to additive nonparametric and fully nonparametric models involving time trend as well as stationary and univariate nonstationary regressors.

With respect to M-estimation, there is an extensive literature. Starting from the seminal work of Huber, (1964), M-estimation methods gradually extend to have general loss functions that may circumvent some drawbacks caused by outliers and so on. The resulting estimators yielded by such approaches are coined as $M$- and $Z$-estimators (see, for example, Chapter 3 of van der Vaart, (1996)).

Using M-estimation, the unknown parameters $(\theta^{\prime}_{0},\beta^{\prime}_{0})^{\prime}$ are estimated by minimizing

where $g(u,v)=g_{1}(u)+g_{2}(v)$ and $\rho(\cdot)$ is some loss function. Possible choices of $\rho(\cdot)$ include, but not limits to,

1. (1)

   Squared error loss: $\rho(u)=u^{2}$.

2. (2)

   $L_{p}$ norm: $\rho(u)=|u|^{p}$ with $p>1$.

3. (3)

   Huber’s loss (Huber, , 1964): For fixed $c>0$,

   $\displaystyle\rho_{c}(u)=\begin{cases}\frac{1}{2}u^{2},&|u|\leq c,\\ c|u|-\frac{1}{2}c^{2},&|u|>c.\end{cases}$

4. (4)

   Least absolute deviation (LAD): $\rho(u)=|u|$.

5. (5)

   Quantile loss (Koenker and Bassett, , 1978): $\rho_{\tau}(u)=u(\tau-I(u<0))$ for $\tau\in(0,1)$.

6. (6)

   Expectile loss (Newey and Powell, , 1986): $\rho_{\tau}(u)=|\tau-I(u<0)|u^{2}$ for $\tau\in(0,1)$.

The loss function $\rho(\cdot)$ here may not be differentiable at every point on $\mathbb{R}$. Often in the literature the subgradient needs to be introduced to tackle these potential erratic points (Chen et al., , 2010). Alternatively, some papers consider the population criterion function like ${\mathbb{E}}[\ell(Z,\theta)]$ instead of $\ell(Z,\theta)$ itself when $Z$ is stationary variable. The rationale why the population criterion is workable is that although $\ell(Z,\theta)$ may not be differentiable at some points, ${\mathbb{E}}[\ell(Z,\theta)]$ may be smooth at these points (see, for example, Chen et al., (2014, p. 640)).

However, in this paper the model has nonstationary regressors that make the population criterion approach fail to work; also the so-called “subgradient approach” depends on particular forms of the loss functions under consideration. We therefore propose a generalized function approach.

The generalized function approach has the following advantages: (1) Generalized functions include ordinary functions that are locally integrable. (2) More importantly, in generalized function sense, generalized functions always have derivatives of all orders. This makes it possible to use Taylor expansion of generalized functions for these loss functions although they may not be necessarily smooth at some points (see Stanković, , 1996).

Phillips, (1991, 1995) are among the first in econometrics to develop a generalized function approach for studying asymptotic properties for LAD and M-estimators in the linear model setting. In particular, the two papers show Taylor expansion for generalized functions, establish (weak and strong) law of large numbers and central limit theorem in generalized function context. In addition, Zinde-Walsh, (2008) uses a generalized function approach in kernel estimation to estimate a density that may not be differentiable at some points. Notice also that the idea of generalized function (concretely, delta-convergent sequence) approaches already have been used in earlier studies in statistics, such as Walter and Blum, (1979) and Susarla and Walter, (1981), where the authors use the property of delta sequence converging to Dirac delta function, $\delta(x)$, to estimate density function. Moreover, the rationale behind the conventional kernel estimation method is also the delta-convergent sequence $K_{h}(x)\equiv h^{-1}K(x/h)\to\delta(x)$ as $h\to 0$ in distribution sense where $K(\cdot)$ is a density (the theorem in Kanwal, , 1983, p. 62). Thus, if $g(\cdot)$ is continuous at $x_{0}$, $\int K_{h}(x-x_{0})g(x)dx\to g(x_{0})$ when $h\to 0$, while kernel estimation in econometrics and statistics is merely a sample version of this integral.

In summary, this paper makes the following contributions to the relevant literature. (i) Model (1.2) accommodates three different types of time series vectors in a multiple-index form; (ii) M-estimation is developed for the unknown parameters that deals with a number of loss functions as special cases, such as OLS, LAD, quantile and expectile estimators; (iii) Non-differentiable loss functions are dealt with by a generalized function approach that provides a unified framework in M-estimation; (iv) The empirical findings show that the proposed model outperform its natural competitors in terms of understanding and evaluating stock return predictability from the points of view of both pseudo out-of-sample ${\rm R}^{2}$ and distributional quantiles; and (v) The finite-sample simulation results show that the proposed model and estimation theory works well numerically.

The rest of the paper is organized as follows. Section 2 introduces the estimation method and then the necessary assumptions before an asymptotic theory is established in Section 3. An empirical analysis about stock return predictability is given in Section 4. Simulation results are then given in Section 5. Section 6 concludes the main sections of this paper. Useful definitions and properties about generalized functions are given in Appendix A, and Appendices B and C include all the necessary mathematical technicalities.

Let $\rho(\cdot)$ be a convex nonnegative function defined on $\mathbb{R}$ that plays a role of loss function. Given the objective function $L_{n}(\theta,\beta)$ in equation (1.3), we define the M-estimators $(\widehat{\theta},\widehat{\beta})$ of $(\theta_{0},\beta_{0})$ that satisfy

where $\epsilon_{n}\to 0$ as $n\to\infty$.

Our theory relies on the assumption that $\rho(\cdot)$ is twice differentiable. When $\rho^{\prime}(\cdot)$ and/or $\rho^{\prime\prime}(\cdot)$ may not exist at some points, we shall invoke its derivatives in a generalized function sense, as discussed in Appendix A below. Indeed, our theory depends heavily on the theory of generalized functions.

We now give some examples of the derivatives of generalized functions we discuss in the sequel.

(a) Let $\rho(u)=|u|$ (LAD loss). Again, by the definition of the derivative of generalized functions, $\rho^{\prime}(u)=-H(-u)+H(u)$, and thus $\rho^{\prime\prime}(u)=\delta(-x)+\delta(x)=2\delta(x)$, where we use $\delta(-u)=\delta(u)$ given in Kanwal, (1983, p. 54).

(b) Let $\rho(u)=u(\tau-I(u<0))$ for some $0<\tau<1$ (Quantile loss). Similarly, we have $\rho^{\prime}(u)=\tau-I(u<0)=\tau H(u)+(\tau-1)H(-u)$. Hence, $\rho^{\prime\prime}(u)=\tau\delta(u)-(\tau-1)\delta(-u)=\delta(u)$.

(c) For Huber’s loss function $\rho(u)$, it has the following derivative,

Moreover, $\rho^{\prime\prime}(u)=I(|u|\leq c)$ in generalized function sense.

(d) For expectile loss function $\rho(u)=|\tau-I(u<0)|u^{2}$, it is a smooth function in ordinary sense with $\rho^{\prime}(u)=2|\tau-I(u<0)|u$, whereas $\rho^{\prime}(u)$ is not smooth at $u=0$. By definition, $\rho^{\prime\prime}(u)=2|\tau-I(u<0)|$.

(e) For $L_{p}$ norm loss $\rho(u)=|u|^{p}$ with $p>1$, we have $\rho^{\prime}(u)=p|u|^{p-1}\text{sgn}(u)$. However, if $1<p<2$, as an ordinary function $\rho^{\prime\prime}(u)$ does not exist at $u=0$. In generalized function sense, $\rho^{\prime\prime}(u)=p(p-1)|u|^{p-2}$, while the meaning of generalized function $|u|^{p-2}$ is discussed in Gel’fand and Shilov, (1964) and Kanwal, (1983) for example.

Before stating assumptions needed for our theoretical development, we introduce some notations. In what follows, $\|\cdot\|$ signifies Euclidean norm for vector or element-wise norm for matrix; $\int$ means an integral (maybe multiple integral) over entire real set $\mathbb{R}$ ( $\mathbb{R}^{p}$).

The structure of the unit root process $x_{t}$ is commonly imposed in the literature since $w_{t}$ has a considerably general form that covers many special cases. See Park and Phillips, (1999, 2001) and Dong et al., (2016). It is straightforward to calculate that $d_{t}^{2}:=\mathbb{E}(x_{t}x_{t}^{\top})=AA^{\top}t(1+o(1))$, and it is well known that, for $r\in[0,1]$, $n^{-1/2}x_{[nr]}\to_{D}B(r)$ as $n\to\infty$ where $B(r)$ is a Brownian motion on $[0,1]$ with covariance $AA^{\top}$.

1. (B.1)

   Let $h(r,v)$ be a $d_{2}$-dimensional vector of continuous functions in $r$ on $[0,1]$.

2. (B.2)

   Suppose that $v_{t}=q(\eta_{t},\cdots,\eta_{t-d_{0}+1};\tilde{v}_{t})$ for $d_{0}\geq 1$, where $\{\tilde{v}_{t},1\leq t\leq n\}$ is also a vector sequence of strictly stationary and $\alpha$-mixing time series, and independent of $\eta_{j}$’s introduced in Assumption A, and the vector function $q(\cdots)$ is measurable such that

   $$\sup_{r\in[0,1]}{\mathbb{E}}\|\dot{g}_{2}^{2}(\beta_{0}^{\top}h(r,v_{1}))\;h(r,v_{1})h(r,v_{1})^{\top}\|^{4}<\infty.$$

   Moreover, the $\alpha$-mixing coefficient $\alpha(k)$ of $\{\tilde{v}_{t},1\leq t\leq n\}$ satisfies $\sum_{k=1}^{\infty}\alpha^{1/2}(k)<\infty$.

Basically, $\{v_{t}\}$ is a vector sequence of strictly stationary variables and $v_{t}$ can be correlated with $x_{t}$ through a finite number of innovations of $\eta$’s. The existence of the fourth moment and the condition on the mixing coefficients guarantee the convergence of the average of $\dot{g}_{2}^{2}(\beta_{0}^{\top}z_{t})z_{t}z_{t}^{\top}$ in probability by virtue of Davydov’s inequality, as shown in Appendix B. All these are quit common but in this paper relies on a known function $h(\cdot)$, so we impose in B.1 the continuity on $h(\cdot,v)$ to make sure that $h(r,v)$ is bounded in $r\in[0,1]$ for each $v$. Let $\xi_{1t}=\rho^{\prime}(e_{t})$, $\xi_{2t}=\rho^{\prime}(e_{t})z_{t}$ and $\xi_{3t}=\rho^{\prime\prime}(e_{t})$.

Condition (C.1) is mild that ensures the eligibility of $\rho(\cdot)$ as a generalized function, apart from being a loss function. In Assumption (C.2), $\rho^{\prime}(\cdot)$ and $\rho^{\prime\prime}(\cdot)$ are the first two derivatives of $\rho(\cdot)$ possibly in the generalized function sense, while they are the same as the derivatives in ordinary sense whenever exist. The filtration in (C.2) can be taken as $\mathcal{F}_{t}=\sigma(e_{1},\cdots,e_{t};x_{s},z_{s},s\leq t+1)$ that is general than series independence between $e_{t}$ and $(x_{t},z_{t})$. When $\rho^{\prime}(\cdot)$ exists in ordinary sense in (C.2), the martingale difference sequence condition renders $\mathbb{E}[\xi_{1t}|\mathcal{F}_{t-1}]=0$ that is easily fulfilled for the loss functions of OLS, LAD, $p$-norm ($p>1$) and Huber’s when the conditional distribution of $e_{t}$ given $\mathcal{F}_{t-1}$ is symmetric, because these functions have odd derivatives. Also, this is the same as $F_{e}(0|\mathcal{F}_{t-1})=\tau$ for quantile loss (see, for example, Assumption 2 in Chen et al., , 2019). Some authors obtain this condition by adjusting the intercept that we also have to deal with when it is violated, like He and Shao, (2000, p. 124) and Xiao, (2009, p. 251). In addition, the positiveness $\mathbb{E}[\xi_{3t}|\mathcal{F}_{t-1}]=a_{2}>0$ is due to the convexity of the loss function under consideration. In particular, $a_{2}=f_{e}(0)>0$ for both check loss and LAD loss, while for Huber’s loss, $a_{2}=P(|e_{t}|\leq c)>0$ in exogenous situation.

When $\rho^{\prime}(\cdot)$ and $\rho^{\prime\prime}(\cdot)$ in (C.2) are genuinely generalized functions, the conditional expectations can be defined as

under some conditions on the conditional density $f_{t,e}(u)$ of $e_{t}$ given $\mathcal{F}_{t-1}$. Though $\rho(\cdot)$ may not be differentiable at some isolated points, $\mathbb{E}[\xi_{1t}|\mathcal{F}_{t-1}]$ and $\mathbb{E}[\xi_{3t}|\mathcal{F}_{t-1}]$ are well-defined given that the conditional density is smooth and is a member in $S$ (the rapid decay test function space, see Appendix A); this is due to the use of generalized functions that extend the ordinary derivatives. The condition D.2 in He and Shao, (2000, p. 125) defines $\mathbb{E}[\xi_{3t}]$ in the same fashion as for non-differentiable loss functions.

Note that (C.3) is quite commonly encountered in the cointegrating regression literature when there is no additional stationary regressor involved (see Park and Phillips, , 2001, Xiao, , 2009 and Wang et al., , 2018). Nevertheless, here the subscript $\rho$ indicates that the joint convergence relies on the loss function. See, for example, equation (3) of Phillips, (1995) shows the functional invariance principle for the derivative of LAD loss, i.e. sign$(\cdot)$.

For simplicity of notation we would mostly suppress these subscript, viz., denote $U_{\rho}(r)$ by $U(r)$, and $V_{\rho}(r)$ by $V(r)$, if there is no confusion raised.

Note that (C.3), in particular, ensures that as $n\rightarrow\infty$

where $B(r)$ has a different covariance function from that of $W(r)$.

Because of the involvement of the unit root process $x_{t}$, we need a similar classification on the function classes, such as $H$-regular and $I$-regular in Park and Phillips, (1999, 2001). However, in our paper we only consider a simpler version of $H$-regular class though a general form can be adopted straightforwardly.

The class of $H$-regular functions mainly contains power functions while the class of $I$-regular functions has all integrable functions on the entire real line as its member, such as probability density functions. More detailed discussion on these definitions can be found in Park and Phillips, (1999, 2001).

It is possible to relax the restriction on the function forms but they are the benchmark. For example, mixtures of $H$-regular and $I$-regular in terms of argument $u$ are feasible for the development in what follows, while we do not pursue it since this would make notation much complicated.

In the process of establishing our asymptotic theory, we shall consider in the proofs that each of the generalized derivatives of $\rho(\cdot)$ is a limit of regular sequence under the setting of generalized functions (see, for example, Phillips, (1991, 1995)).

Remark.  It can be seen from the theorem that the imposition of $H$-regular functions $g_{1}(u)$ and its derivatives delivers a fast rate of convergence for $\widehat{\theta}$ due to the divergence of unit root process $x_{t}$. This is comparable with the result in Park and Phillips, (2001). On the other hand, the estimator $\widehat{\beta}$ possesses a usual square-root-$n$ rate, although the regressor $z_{t}$ has a deterministic trending component. Notice that the positive definiteness of $\Sigma$ is easily satisfied as long as $v_{1}$ is a continuous variable. More importantly, it is readily seen that the different loss functions take a role in the asymptotic limits through $a_{1}$ and $a_{2}$ that may affect the covariance but not the rates of convergence.

In order to make a comparison with the literature, here we simply suppose the model is exogenous. (1) When $\rho(u)=|u|$, $\rho^{\prime}(u)=-H(-u)+H(u)$ with $H(\cdot)$ being Heaviside function and $\rho^{\prime\prime}(u)=2\delta(u)$. We have $a_{1}=\mathbb{E}[\rho^{\prime}(e_{t})]^{2}=1$ and $a_{2}=\mathbb{E}[\rho^{\prime\prime}(e_{t})]=2f_{e}(0)>0$ where $f_{e}(\cdot)$ is the density of $e_{t}$. Moreover, if $g_{1}(u)=u$ and $g_{2}(v)\equiv 0$, the model reduces to a linear cointegrating model. Our result implies that

that is the result of Phillips, (1995) without correlation between $x_{t}$ and $e_{t}$; if $g_{2}(v)=v$ and $g_{1}(u)\equiv 0$, $h(r,v)\equiv v$ the model reduces to a linear model with stationary regressors. Our result also naturally covers that $\sqrt{n}(\widehat{\beta}-\beta_{0})\to_{D}[2f_{e}(0)]^{-1}[\mathbb{E}(v_{1}v_{1}^{\top})]^{-1/2}N(0,I_{d_{2}})$ as shown in Pollard, (1991).

(2) When $\rho(u)$ is the quantile loss, $\rho^{\prime}(u)=\tau H(u)+(\tau-1)H(-u)$ and $\rho^{\prime\prime}(u)=\delta(u)$. We have $a_{1}=\mathbb{E}[\rho^{\prime}(e_{t})]^{2}=\tau(1-\tau)$ and $a_{2}=\mathbb{E}[\rho^{\prime\prime}(e_{t})]=f_{e}(0)>0$ where $f_{e}(\cdot)$ is the density of $e_{t}$. Further, if $g_{1}(u)=u$ and $g_{2}(v)\equiv 0$, the model reduces to a linear cointegrating model. Our result implies that $n(\widehat{\theta}-\theta_{0})\to_{D}[f_{e}(0)]^{-1}[\int_{0}^{1}B(r)B(r)^{\top}dr]^{-1}\int_{0}^{1}B(r)dU(r)$ that is the result of Xiao, (2009) without correlation between $x_{t}$ and $e_{t}$; if $g_{2}(v)=v$ and $g_{1}(u)\equiv 0$, $h(r,v)\equiv v$ the model reduces to a linear model with stationary regressor. Our result implies that $\sqrt{n}(\widehat{\beta}-\beta_{0})\to_{D}\sqrt{\tau(1-\tau)}[f_{e}(0)]^{-1}[\mathbb{E}(v_{1}v_{1}^{\top})]^{-1/2}N(0,I_{d_{2}})$, which is the result of quantile regression for linear model. See Koenker and Bassett, (1978). $\Box$

While $\Sigma$ can easily be consistently estimated by the function $h$ and observations of $v_{t}$, we need to construct consistent estimators for both $a_{1}$ and $a_{2}$ in order to use the asymptotic limits for statistical inferences. These estimators are given by $\widehat{a}_{1}=\frac{1}{n}\sum_{t=1}^{n}[\rho^{\prime}(\widehat{e}_{t})]^{2}$ and $\widehat{a}_{2}=\frac{1}{n}\sum_{t=1}^{n}\rho^{\prime\prime}(\widehat{e}_{t}$), respectively, where $\widehat{e}_{t}=y_{t}-g(\widehat{\theta}^{\,\top}x_{t},\widehat{\beta}^{\,\top}z_{t})$ for $1\leq t\leq n$. We then have the following corollary.

The ergodicity for the two sequences is quite weak and can therefore be fulfilled under several sufficient conditions, such as requiring $e_{t}$ be a martingale difference sequence or a strictly stationary and $\alpha$-mixing sequence. When the derivations of the loss function $\rho(u)$ do not exist in ordinary sense, its regular sequence may be used for $\rho(u)$ (see Appendix A and Phillips, (1995, p 918)).

Now we consider the regression function $g(u,v)=g_{1}(u)+g_{2}(v)$ where both $g_{1}(u)$ and $g_{2}(v)$ are smooth but $g_{1}(u)$ is integrable on the entire real line stipulated by Assumption D(3). This assumption will change the asymptotic theory drastically. The limit theory depends on the local time of some scalar Brownian motion $W(r)$ defined as

which measures the sojourning time of $W(r)$ around $s$ during the time period $(0,p)$. In our context, the scalar Brownian motion may be $B_{1}(r):=\theta_{0}^{\top}B(r)$ where $B(r)$ is defined by (2.2).

Another feature of the limit theory about this regression function is that we derive in a new coordinate system first, then recover the limit in the original coordinate. To do so, using $\theta_{0}$ that is a nonzero but not necessarily a unit vector, we construct an orthogonal matrix $P=(\theta_{0}/\|\theta_{0}\|,P_{1})_{d_{1}\times d_{1}}$ to rotate all the vectors of interest, including $\theta_{0}$, $x_{t}$ and $B(r)$. The relevant rotated vectors are $x_{1t}\equiv\theta_{0}^{\top}x_{t}$ and $x_{2t}\equiv P_{1}^{\top}x_{t}$ that after normalization converge jointly to $B_{1}(r)\equiv\theta_{0}^{\top}B(r)$ and $B_{2}(r)\equiv P_{1}^{\top}B(r)$; and we shall denote by $L_{1}(p,s)$ the local time of $B_{1}(r)$.

Noting that $M$ is a block diagonal matrix, the result in Theorem 3.2 implies

as $n\to\infty$ where $J=\text{diag}(1/\|\theta_{0}\|,I_{d_{1}-1})$.

Notice also that in the condition of D.2, $g_{2}(v)$ is the same as in D.1. However, the convergence of (3.2) has not been bothered by the first additive component $g_{1}(v)$, because the convergence is similar to the existing literature except our setting for the regressor is more complicated (see Pollard, , 1991 and Koenker and Bassett, , 1978). By sharp contrast, the asymptotic limit of $\widehat{\beta}$ in Theorem 2.1 is affected by the limit of the regressor in the first component function. This is due to the drastic different behavior between the H-regular and I-regular functions of unit root processes. See, for example, Park and Phillips, (2001) and Dong et al., (2016).

The convergence of $\eqref{int1}$ is similar to Theorem 3.1 in Dong et al., (2016) where $J=I_{d}$ due to the identification condition $\|\theta_{0}\|=1$ in semiparametric single-index model, $a_{1}=4\sigma_{e}^{2}$ and $a_{2}=2$ because $\rho(u)=u^{2}$ in the paper. It is also similar to Theorem 5.1 in Park and Phillips, (2001) where the regressor is univariate but the parameter is multivariate. In addition, when the loss function reduces to LS loss, the convergence of $\eqref{int1}$ is the same as Theorem 2 in Park and Phillips, (2000).

If one is concerned about the estimate of the matrix $P_{1}$, note that $P_{1}$ is an orthogonal system of $\theta_{0}^{\bot}$ (orthogonal complement space), so that once $\widehat{\theta}$ is available, we may define $\widehat{P}_{1}=\widehat{\theta}^{\bot}$. Normally, one does not need to rotate the coordinate system in practice that is introduced to facilitate our theory only.

Using the block representation of $M_{1}=(m_{ij})_{2\times 2}$, we have $M_{1}^{-1}=(\tau_{ij})_{2\times 2}$ with