Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
88,617 characters · 11 sections · 100 citation commands
Robust M--Estimation for Additive Single--Index Cointegrating Time Series Models
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 therefore propose a parametric additive single-index model of the form:
for $1\leq t\leq n$, where $\gamma_{ij}^{0}\in \mathbb{R}^{1}$ and $\theta_{ij}^{0}\in \mathbb{R}^{d_i}$, for $1\leq j\leq p_i$ and $1\leq i\leq 2$, are all unknown parameters, $g_{ij}(\cdot)$, for $1\leq j\leq p_i$ and $1\leq i\leq 2$, are all known functions, $p_1$ and $p_2$ are known positive integers, $x_t$ is a $d_1$--dimensional vector of nonstationary time series, $z_t$ is a $d_2$--dimensional vector of trending stationary time series, and $e_t$ is the error term.
Before we discuss about how to estimate model ((ref)) and then summarize the main contributions of this paper, we would like to first point out some novel features of model ((ref)) and explain why we start with model ((ref)) in this paper.
(i) Model ((ref)) offers a simple way of dimension reduction by involving a partial sum of a number of linear combinations of the form: $\left(x_t^{\top} \theta_{11}^{0}, \cdots, x_t^{\top} \theta_{1, p_1}^{0}; z_t^{\top} \theta_{21}^{0}, \cdots, z_t^{\top} \theta_{2, p_2}^{0}\right)$, which naturally use a principal component analysis approach to obtaining the first $p_1$ leading components from $x_t$, and also the first $p_2$ leading components from $z_t$. Even though in this paper we focus on the case where $d_1, p_1, d_2$ and $p_2$ do not diverge with $n$, the actual values of $d_1$ and $d_2$ can be very large in both theory and practice, and much larger than $p_1$ and $p_2$, respectively.
(ii) Model ((ref)) may be considered as a parametric neural network (NN) model with one hidden layer. In this case, $\{g_{ij}(\cdot): 1\leq j\leq p_i, 1\leq i\leq 2\}$ may be generated by a sequence of activation functions of the form: $g_{ij}(\cdot)= \sigma_i(\cdot + c_{ij})$, in which each $\sigma_i(\cdot)$ may be chosen as a commonly used activation function, and each $c_{ij}$ is the location parameter. The relevant NN literature includes cg1989, cs1998, cw1999, crs2001, chen2007, kk2017, bk2019, and sh2020.
(iii) Note that model ((ref)) covers a wide class of linear combinations of single--index parameters and basis functions commonly used in the NN literature for approximating unknown nonparametric functions by such linear combinations, as discussed in the NN references cited above. Note also that the choice of $\{g_{ij}(\cdot): 1\leq j\leq p_i, 1\leq i\leq 2\}$ allows for many commonly used polynomials, such as Hermite polynomials used in the empirical analysis in Section 5.2 below.
(iv) The functional forms in model ((ref)) considerably extend those studied in the relevant literature by dgp2015, dgd2016, degui2016, degui2021 and others. As discussed below, moreover, we are able to consistently estimate both the single--index parameters, $\{\theta_{ij}^{0}\}$, and their coefficients, $\{\gamma_{ij}^{0}\}$, under the data structure outlined in model ((ref)) below.
(v) The parametric form: $m(x, z; \theta_0, \gamma_0) = \sum_{j=1}^{p_1}\gamma_{1j}^{0}\,g_{1j}(x^\top\theta_{1j}^{0})+ \sum_{j=1}^{p_2}\gamma_{2j}^{0}\,g_{2j}(z^\top\theta_{2j}^{0})$ represents the true conditional mean when the null hypothesis $H_0: \, E\left[y_t|(x_t=x, z_t=z)\right] = m(x, z; \theta_0, \gamma_0)$ is true, as studied in dgdy2017 for a special case of model ((ref)). We therefore think that it is both important and necessary to investigate model ((ref)) before we may study estimation and inferential problems for various semi--parametric forms of model ((ref)).
We now come back to depict time series sequence $\{x_t, z_t, \tau_t, 1\leq t \leq n\}$ in model (ref). We specify the respective structures of $x_t$ and $z_t$ by
where $\tau_t=\frac{t}{n}$, as defined later, vector $w_t$ is a linear 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 (ref), $x_t$ is a vector of linearly integrated processes of order one. It is 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.
Different from existing settings, model ((ref)) takes into account not only the `large' stochastic trending component driven by $x_t^{\top} \theta_{1j}^{0}$, but also the `small' slowly varying component $z_t^{\top} \theta_{2j}^{0}$ mainly represented by the deterministic time trend $\tau_t$.
Note that the products of $\gamma_{ij}$ and $g_{ij}(\cdot)$ may give raise to an identification issue. For example, if $g_{11}(u)=u$, $\gamma_{11}g_{11}(x_{t}^\top\theta_{11}^0)=\gamma_{11}x_{t}^\top\theta_{11}^0$. For the sake of identification, we shall assume that $\|\theta_{ij}^0\|=1$ for all $1\leq j \leq p_i$ and $1\leq i\leq 2$, and possibly the first elements of $\theta_{ij}^0$ are positive that is not necessary in the case that $g_{ij}(\cdot)$ are nonlinear.
We shall adopt robust M--estimation (robust estimation hereafter for simplicity) method to estimate all unknown parameters. With regard to robust estimation, there is an extensive literature. Starting from the seminal work of huber1964 with Huber's loss, robust estimation methods gradually extend to some general loss functions extensively studied in the literature, such as bai1992, fanjq1994, fan2018 and degui2021.
By robust estimation, the unknown parameters $\theta, \gamma$ are estimated by minimizing
where $\rho(\cdot)$ is some loss function, $\theta$ and $\gamma$ stand for generic vectors in a compact sets including the true vectors $\gamma_i^{0}$ and $\theta_{ij}^{0}$, $i=1,2$, $j=1,\cdots,p_i$, as interior points.
Though there are several loss functions in the robust estimation literature, the most important ones are: (1) Least absolute deviation (LAD), $\rho(u)=|u|$; (2) Quantile loss koenker1978, $\rho_{\tau}(u)=u(\tau-I(u<0))$ for $\tau\in (0,1)$; and (3) Huber's loss huber1964, i.e. for fixed $c>0$,
The first two loss functions are not smooth while the third does not have second derivative at $\pm c$. This feature makes aymptotic analysis of the estimators much difficult than usual where the loss is sufficiently smooth.
We shall adopt generalized function approach to dealing with the nonsmooth loss function $\rho(\cdot)$. Generalized functions started with the use in physics, and the most famous one is Dirac delta function that is actually not a function in ordinary sense. Mathematicians make the concept of generalized functions rigorous by defining them as functionals acting on tempered test function space $S$ (another test space is $D$ but $D\subset S$), introduced in the supplementary material of the paper. Because of the property of functions in $S$, any local integrable function with no faster than a polynomial increase at infinity is a generalized function, and all loss functions in the literature belong to the category. See, for example, gelfand1964 and kanwal1983. In particular, LAD and quantile losses have Dirac delta function $\delta(\cdot)$ as their second-order generalized derivatives.
As is well known, generalized functions can not be regarded as ordinary functions because they might not be well-defined at each point. Nonetheless, their operation is valid under integration with certain conditions. For example, for any random variable $e$ with density function $f_e(x)$,
provided $f_e(x)$ is continuous at $x=0$. Thus, in our proposed mechanism we use generalized functions only in expectations, whereas if we have to consider the derivatives of loss functions, we use their regular sequences defined later.
The focus of this paper is on non-differentiable convex loss functions $\rho(\cdot)$. It is well known that the subgradient $\psi(\cdot)$ of convext $\rho(\cdot)$ always exists, see rockafellar1970; subgradient will coincide with derivative wherever it is differentiable. The step stone in dealing with non-differentiable convex loss $\rho(\cdot)$ is using a so-called regular sequence $\{\rho_m(\cdot), m=1,2, \cdots\}$ that approximates to $\rho(\cdot)$, as $m\to\infty$; and $\rho_m(\cdot)$ are infinitely differentiable, $\rho_m^{(k)}(u)\to \rho^{(k)}(u)$ for any $k$ in generalized function sense. More importantly, we establish the rates of $\rho_m(u)-\rho(u)$, $\rho_m'(e_t)-\psi(e_t)$ in probability and ${\mathbb E}[\rho_m''(e_t)-\rho''(e_t)]$ in Theorem 2.1 below that are crucial for our theory, where $\rho''(\cdot)$ is understood in generalized function sense such that ${\mathbb E}[\rho''(e_t)]$ exists. As far as we are aware, this is the debut of these rates in the literature.
Closely related to our generalized function approach, phillips1991, phillips1995 use generalized function approach to studying asymptotic properties for LAD and M-estimators in linear models. From then on, to the best of our knowledge, there has been no attempts for using generalized functions to tackle robust estimation. The Achilles heel to the approach is that the approximation rate of regular sequence to nonsmooth loss functions is not established in phillips1991, phillips1995, and this approximation rate nevertheless is crucial in the theoretical side that, as what aforementioned, we shall establish in Section two below. Notice also that the regular sequence (also delta-convergent sequence) approach already has been used in earlier studies in statistics to estimate density function, such as walter1979 and walter1981. Thus, we believe that the rationale of the estimation method we propose considerably enriches the relevant literature.
Our methodology is to approximate $L_n$ by a quadratic form $Q_{n}$ of parameters whose minimizer has an explicit expression, and we show that the difference between the minimizers of $L_n$ and $Q_{n}$ is negligible, so that the limit of $L_{n}$'s minimizer is the same as that of $Q_{n}$'s. This mechanism is much simpler and more generally applicable than those developed in the existing literature, such as knight1989, phillips1991, phillips1995, pollard1991, bai1992, fanjq1994, vaart1996, gao2009, fan2018 and degui2021. In addition, the estimator's explicit expression derived in our paper is comparable with Bahadur representation in the relevant literature; under the same setting, our remainder term is $O_P(n^{-1/2+\lambda})$ for any $\lambda\in (0,1/2)$, whereas bahadur1966 has the remainder $O_{a.s.}(n^{-1/4}\log(n))$.
To illustrate our methodology of the proposed generalized function approach, we use a linear model as an exemplar in Section 2 to show the essence of the approach; in Section 3 we present all assumptions. In Section 4 we establish the corresponding asymptotic distributions of the proposed estimators of all the parameters in model (ref) in two scenarios where $g_{1j}(\cdot)$ are $H$-regular and $I$-regular, respectively. Section 5.1 gives Monte Carlo simulation results to illustrate the performance of our estimators in finite sample situations; an empirical study about stock returns is shown in Section 5.2. Technical lemmas are given in Appendix A, whose proofs however are shown in supplementary materials; all main results are proven in Appendix B.
Let us finally introduce the following notation. $\|\cdot\|$ is either Euclidean norm for vectors or element-wise norm for matrices; $\equiv$ means equal by definition; $\int$ means an integration on $\mathbb{R}$; $C$, $C_1$ and so on represent absolute constants that may be different at each appearance; $\dot{g}, \ddot{g}$ signify the first and second derivatives of $g$, respectively.
This section will show the rationale of robust estimation based on generalized function approach, temporarily letting model (ref) alone. To begin with, we study the convergence rate of regular sequence of loss functions where regular sequence is a tool by which generalized functions are defined. See Chapter 5 of milne1980.
The loss functions $\rho(\cdot)$ we use satisfy the following assumption. Recall that $\psi(\cdot)$ is the subgradient of $\rho(\cdot)$. Though the value of $\psi(\cdot)$ at some points may not be unique, we allow $\psi(\cdot)$ to take any of them while the analysis in the sequel remains unchanged. The focus then rests on the second derivative of $\rho(\cdot)$ that may not exist in ordinary sense.
Assumption 2.1.
By convexity, $\rho(\cdot)$ is a continuous function. See, Corollary 10.1.1 of rockafellar1970. The other conditions in (a) (local integrability and no increasing faster than a polynomial) make operation of generalized function in the sequel valid, and these are certainly fulfilled by all loss functions encountered in the literature.
The Lipschitz condition in (b) of the assumption is the center that depicts a technical requirement in the following analysis, and it is guaranteed by the boundedness of subgradient. Notice also that it is satisfied by several most popular loss functions. For example, when $\rho(u)=|u|$, we have $C=1$; when $\rho(u)$ is the check function with parameter $\tau\in (0,1)$, we have $C=\max(\tau, 1-\tau)$; when $\rho(u)$ is Huber's loss with parameter $c$, we have $C=c$.
We next give crucial results about regular sequences of loss functions that is used as a bridge between nonsmooth loss and its approximation counterpart in generalized function context.
The proof is given in Appendix B. Note that each $\phi_m(x)$ is the density of variables $N(0, (2m)^{-1})$ which form a delta-convergent sequence, i.e. $\phi_m(x)\to \delta(x)$ as $m\to \infty$ in generalized function sense; see kanwal1983 and stein2003. Note also that in the last assertion, because of convexity of $\rho(u)$ and $\rho_m''(u)\to \rho''(u)$ in generalized function sense, ${\mathbb E}[\rho_m''(e)]>0$ for large $m$.
{\bf Remark 2.1}. \ Theorem (ref) is of independent interest. The importance of Theorem (ref) is the convergence rate $\sup_{u\in \mathbb{R}}|\rho_m(u)-\rho(u)|\le Cm^{-1/2}$. To the best of our knowledge, this is first shown in the literature although regular sequences are discussed in several papers, such as phillips1991, phillips1995. It is due to this rate that our analysis below is on a solid and rigorous ground. In Figure (ref) above three losses and their regular sequences are plotted to visualize the approximations established in Theorem (ref). Note from the proof that the Lipschitz condition in Assumption 2.1 can be relaxed as $|\rho(x)-\rho(y)|\le C|x-y|^\alpha$ for $\alpha\in (0,1]$. All the results below under the relaxation hold with a change on the choice of $m$. \qed
To show the essence of our generalized function approach, in the sequel we shall illustrate through a linear parametric regression. To do so, let us consider a very simple regression model of the form:
where unknown parameter $\theta_0\in \Theta$, a compact subset of $\mathbb{R}^d$. Notice that the notation in this section has a different meaning from other sections.
The estimator of $\theta_0$ in the regression (ref) is defined by
where $\rho(\cdot)$ satisfies Assumption 2.1.
To analyse $\widehat{\theta}$, instead of considering minimization over $\theta\in\Theta$, we shall follow the literature to consider vectors $\beta$ in the tangent cone $T_{\Theta}(\theta_0)$ of $\Theta$ at $\theta_0$, that is, $\beta=\sqrt{n}(\theta-\theta_0)$ where $\theta\in\Theta$. See among others, bickel1974, badu1989, davis1992, phillips1995 and gao2009. In particular, charles1994 argues that $\beta=\sqrt{n}(\theta-\theta_0)\in T_{\Theta}(\theta_0)$, and shows the convergence of $\widehat\theta$ through that of $\widehat\beta$.
Therefore, we focus on the minimization of $\Pi_n(\theta)-\Pi_n(\theta_0)$ instead of $\Pi_n(\theta)$ in (ref), and write
Hereby, the objective function is reparametrized in $\beta\in \mathbb{R}^d$, and when $\widehat\beta$ minimizes $\tilde{\Pi}_n(\beta)$, $\widehat\theta$ minimizes $\Pi_n(\theta)$ where $\widehat\beta=\sqrt{n}(\widehat\theta -\theta_0)$.
Notice that, to illustrate our method manifestly, we impose very strong conditions on the regressors and error term. Such strong conditions are used only in this section. It is clear that
and it is the unique minimizer of $Q_n(\beta)$.
{\bf Remark 2.2}. \ Note that pollard1991 establishes `Convexity Lemma' and points that if convex sequence $\lambda_n(\theta)$ converges to $\lambda(\theta)$ in probability for each $\theta\in \Theta$ where $\Theta$ is convex and open subset of $\mathbb{R}^d$, then $\sup_{\theta\in K}|\lambda_n(\theta)- \lambda(\theta)|\to 0$ in probability for any compact $K\subset \Theta$. This is much weaker than our results (ref) and (ref). In addition, fan2003 use quadratic approximation for smooth objective function and derive limit theory from this quadratic function. By contrast, the generalized function approach we propose mainly deals with nonsmooth loss functions although it can also be used for smooth losses. \qed
{\bf Remark 2.3}. \ Interestingly, equation (ref) gives $\widehat{\beta}=\widehat\beta_Q+o_P(n^{-1/2+\lambda})$ for any $\lambda\in (0, 1/2)$ that is better than Bahadur representation in some sense. See, for example, bahadur1966 where the reminder term is of order $O_{a.s.}(n^{-1/4}\log(n))$. \qed
{\bf Remark 2.4}. \ As argued in pollard1991, quadratic approximation to objective function of M estimation avoids technical difficulty such as stochastic equicontinuity in an asymptotic proof. Indeed, once we establish (ref), the limit of $\widehat{\beta}$ is the same as $\widehat\beta_Q$ that is quite easily derived from its explicit expression and some regular conditions such as ${\mathbb E}[\psi(e_t)]=0$. \qed
We shall use this mechanism for model (ref) in the following section.
Now we turn back to model (ref), and we shall use the generalized function approach for the establishment of the asymptotic limits of the proposed robust estimators. We give the following assumptions.
{\bf Assumption 3.1}
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 phillips1999,phillips2001 and dgd2016. It is straightforward to calculate that $d_t^2\equiv\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_DB(r)$ as $n\to\infty$ where $B(r)$ is a standard Brownian motion on $[0,1]$ with covariance $AA^\top$.
In the case where some of $x_t^{\top} \theta_{1j}^{0}$ or even all of them reduce to stationary processes, the estimation method of $\theta_{1j}^{0}$, for $1\leq j \leq p_1$, and the other parameters remains the same. The corresponding asymptotic properties, however, become different. In such cases, the corresponding assumptions on $g_{1j}(\cdot)$ are roughly the same as those on $g_{2j}(\cdot)$. Since theoretical developments for such cases require substantially different technologies, we wish to leave this for future research.
Assumption 3.2
Basically, $\{v_t\}$ is a vector sequence of strictly stationary time series 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}_{2j}^2(\theta_{2j}^{0\top} z_t)z_t z_t^\top$ in probability by virtue of Davydov's inequality, as shown in the supplementary material of the paper. Since this relies on a known function $h(\cdot, \cdot)$, so we impose in 3.2(i) the continuity on $h(\cdot, v)$ to make sure that $h(r,v)$ is bounded in $r\in [0,1]$ for each $v$.
Recall that $\psi(\cdot)$ is the subgradient of $\rho(\cdot)$ that satisfies Assumption 2.1.
Assumption 3.3
which is a vector of standard Brownian motion processes of $(d_1+1)$ dimension, where $U_\rho(r)$ and $W(r)$ have variance and covariance $a_1$ and $I_{d_1}$, respectively.
In Assumption 3.3(i), $\psi(\cdot)$ and $\rho''(\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 3.3(i) can be taken as $\mathcal{F}_{t}=\sigma(e_1,\cdots, e_t; x_s, z_s, s\le t+1)$ that is less restrictive than mutual independence between $e_t$ and $(x_t, z_t)$. The martingale difference sequence condition renders $\mathbb{E}[\psi(e_t)|\mathcal{F}_{t-1}]=0$ that is easily fulfilled for the loss functions of LAD and Huber's when the conditional distribution of $e_t$ given $\mathcal{F}_{t-1}$ is symmetric, because these functions $\psi(\cdot)$ are odd. Also, this condition is equivalent to $F_e(0|\mathcal{F}_{t-1}) =\tau$ for quantile loss (see, for example, Assumption 2 in degui2019). Some authors obtain this condition by adjusting the intercept that we also have to deal with when it is violated, like hexuming2000 and xiao2009a.
When $\rho''(\cdot)$ in 3.3(i) is genuinely a generalized function, 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 points, $\mathbb{E}[\psi(e_t)|\mathcal{F}_{t-1}]$ and $\mathbb{E}[\rho''(e_t)|\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 C in the supplementary material); this is due to the use of generalized functions that extend the ordinary derivatives. The condition D.2 in hexuming2000 defines $\mathbb{E}[\rho''(e_t)]$ in the same fashion as for non-differentiable loss functions.
In addition, the positiveness $\mathbb{E}[\rho''(e_t)|\mathcal{F}_{t-1}] =a_2>0$ is due to the convexity of the loss functions under consideration. In particular, $a_2=f(0)>0$ for both check loss and LAD loss with $f(\cdot)$ being the density of $e_t$, while for Huber's loss, $a_2= P(|e_t|\le c)>0$ in exogenous situation.
Condition 3.3(ii) is used as a technical requirement in related calculations. See Lemma (ref). This holds in particular when the regressors are independent of the error terms.
Note that 3.3(iii) is quite commonly encountered in the cointegrating regression literature when there is no additional stationary regressor involved (see phillips2001, xiao2009a and qiying2018). Nevertheless, here the subscript $\rho$ indicates that the joint convergence relies on the loss function. For example, equation (3) of phillips1995 shows the functional invariance principle for the derivative of LAD loss, i.e. sign$(\cdot)$. For simplicity of notation we suppress the subscript, viz., denote $U_\rho(r)$ by $U(r)$, if there is no confusion raised.
Note that 3.3(ii), in particular, ensures that as $n\rightarrow \infty$
where $B(r)$ is a vector of Brownian motion with covariance matrix $AA^\top$.
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 phillips1999,phillips2001. However, in our paper we only consider a simpler version of $H$-regular class, that is, power functions, because our model incorporates an additive form of $H$-regular functions.
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 phillips1999, phillips2001.
Assumption 3.4 \ \ Let all $g_{ij}(\cdot)$, for $1\leq j \leq p_i$ and $1\leq i\leq 2$, be continuously differentiable up to the second order. Suppose further that
Assumption 3.4 requires all $g_{ij}(\cdot)$, for $1\leq j \leq p_i$ and $1\leq i\leq 2$ to be continuously differentiable up to the second order, while considering two classes of functions for all $g_{1j}(u)$, $j=1,\ldots,p_1,$ that is, $H$-regular and $I$-regular. This is because, when $\theta'x_t$ is a unit root process, the behaviour of $g(\theta'x_t)$ depends heavily on whether $g(\cdot)$ is $H$-regular or $I$-regular (see, for example, dgd2016). More importantly, in $I$-regular case we require that all index vectors be the same $\theta_{1j}^0=\theta_1^0$, $j=1,\ldots, p_1$. That is, all integrable functions share the common stochastic trend $x_t'\theta_1^0$, since otherwise there would be a technical bottle neck. Under this setting, the first part of the regression function can be rewritten as
where $g_1(x^\top \theta_{1}^0; \gamma_{1}^0)$ is known up to $\theta_{1}^0$ and $\gamma_{1}^0$.
The asymptotic behaviours of the proposed estimators of the parameters associated with the nonstationary part in our regression equation depend heavily on the functional functions of $g_{1j}(\cdot)$, $j=1,\cdots, p_1$. We therefore establish our asymptotic theory according to $H$-regular and $I$-regular functions, respectively.
We firstly consider the case of $H$-regular functions, that is, Assumption 3.4(i) holds. For better exposition, denote $x_{nt}=n^{-1/2}x_t$ and two vectors associated with the regression functions and regressors by
In addition, we define a symmetric matrix $B$ and a vector $B_1$ given by
where the positive constants $a_1, a_2$ are given by Assumption 3.3 and processes $B(r)$ and $U(r)$ are defined by (ref). Also, let
where $\widehat\theta_{1j}$, $\widehat\theta_{2j}$, $\widehat\gamma_{1j}$ and $\widehat\gamma_{2j}$ are the minimizers of the objective function $L_n(\theta, \gamma)$ defined by (ref).
The proof of the theorem is given in Appendix B. Note that in the proof we make use of the regular sequence of $\rho(\cdot)$ defined in Theorem 2.1, and then similar to Theorem 2.2, the objective function can be approximated by a quadratic function. Hence, the consistency and convergence rate of all estimators are implied by $\widehat{\Lambda}=O_P(1)$.
As can be seen from the asymptotics in (ref), all convergence rates of $\widehat\theta_{1j}$ and $\widehat\gamma_{1j}$ are affected by the $H$-regular property of $g_{1j}(\cdot)$, $j=1,\cdots, p_1$. Precisely, $\widehat\theta_{1j}$ has convergence rate $[\dot{\nu}_j(\sqrt{n})n]^{-1}$ where $\nu_j(\cdot)$ is the asymptotic order of $g_{1j}(\cdot)$. For example, if $g_{11}(u)=u$, $\nu_1(\lambda)=\lambda$, then $\widehat\theta_{11}$ has rate $n^{-1}$, so-called super rate in nonstationary context; if $g_{12}(u)=u^2$, $\nu_1(\lambda)=\lambda^2$, then $\widehat\theta_{12}$ has rate $n^{-3/2}$.
Notice that, given the identification condition $\|\theta_{1j}^0\|=1$, we are able to identify $\gamma_{1j}$. Moreover, the rate of $\widehat\gamma_{1j}$ is affected by $g_{1j}(\cdot)$ too, as it converges at rate of $[\nu_j(\sqrt{n})\sqrt{n}]^{-1}$. Since $x_t^\top \theta_{1j}^0$ is still a unit root process, $g_{1j}(x_t^\top \theta_{1j}^0)$ has asymptotic order $\nu_j(\sqrt{n})$ that results in a rapid rate for $\widehat\gamma_{1j}$ by a factor $[\nu_j(\sqrt{n})]^{-1}$, comparing with stationary case.
By contrast, all the estimators $\widehat\theta_{2j}$ and $\widehat\gamma_{2j}$, $j=1,\cdots, p_2$, have usual square-root-$n$ convergence rate, although these parameters are involved in the regression accommodating both stationary and nonstationary variables. This is simply because the regression function takes additive form so that the stationary part and the nonstationary part are separated, and thus the parameters in the stationary part retain the usual square-root-$n$. As can be seen below, however, the limits of all estimators in both parts include all ingredients in the regression. This was also found in DL2018.
In order to obtain the asymptotic properties of $\widehat{\Lambda}_1$ and $\widehat{\Lambda}_2$ separately, we partition the matrix $B=(B_{ij})_{2\times 2}$ and vector $B_1=(b_{11}^\top, b_{12}^\top)^\top$ conformably,
Then, (ref) implies that
To understand this, take $p_1=p_2=1$ as example. Consider
where $g_1(\cdot)$ is $H$-regular with asymptotic order $\nu_1(\cdot)$, $\mathbb{E}z_t=0$ that implies $B_{12}=B_{21}^\top =0$. Theorem (ref), in particular (ref), gives
where all notation $B_{11}, B_{22}, b_{11}, b_{12}$ can be explicitly given from the above discussion.
It is noteworthy that the impact of the loss function to estimation is through the constants $a_1$ and $a_2$ that are easily separated from the matrices and vectors, so the limits of $\widehat{\Lambda}_1$ and $\widehat{\Lambda}_2$ have factors $1/a_2$ and $\sqrt{a_1}/a_2$, respectively.
In order to make a comparison with the relevant literature, here we simply suppose the model is exogenous. (1) When $\rho(u)=|u|$, $\rho'(u)=-H(-u)+H(u)$ with $H(\cdot)$ being Heaviside function and $\rho''(u)=2\delta(u)$. We have $a_1=\mathbb{E}[\rho'(e_t)]^2=1$ and $a_2=\mathbb{E}[\rho''(e_t)] =2f_e(0)>0$ where $f_e(\cdot)$ is the density of $e_t$. Moreover, if $g_1(u)=u$, $\gamma_1^0=1$ needed not to estimate, and $g_2(v)\equiv 0$, the model reduces to a usual linear cointegrating model. Our result implies that $$n(\widehat \theta_1-\theta_1^0)\to_D [2f_e(0)]^{-1}\left[\int_0^1B(r)B(r)^\top dr\right]^{-1} \int_0^1 B(r)dU(r),$$ which is the same as phillips1995 without correlation between $x_t$ and $e_t$; if $g_2(v)=v$, $\gamma_2^0=1$ needed not to estimate, 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 \theta_2-\theta_2^0)\to_D[2f_e(0)]^{-1} [\mathbb{E}(v_1v_1^\top)]^{-1/2} N(0,I_{d_2})$ as shown in pollard1991.
(2) When $\rho(u)$ is the quantile loss, $\rho'(u)=\tau H(u)+(\tau-1)H(-u)$ and $\rho''(u)=\delta(u)$. We have $a_1=\mathbb{E}[\rho'(e_t)]^2 =\tau(1-\tau)$ and $a_2=\mathbb{E}[\rho''(e_t)]=f_e(0)>0$ where $f_e(\cdot)$ is the density of $e_t$. Further, if $g_1(u)=u$, $\gamma_1^0=1$ needed not to estimate and $g_2(v)\equiv 0$, the model reduces to a linear cointegrating model. Our result implies that $n(\widehat \theta_1-\theta_1^0)\to_D [f_e(0)]^{-1} [\int_0^1B(r)B(r)^\top dr]^{-1} \int_0^1 B(r)dU(r)$ that is the same result as in xiao2009a without correlation between $x_t$ and $e_t$; if $g_2(v)=v$, $\gamma_2^0=1$ need not to be estimated and $g_1(u)\equiv 0$, $h(r,v)\equiv v$, the model reduces to a linear model with one stationary regressor. Our result implies that $\sqrt{n}(\widehat \theta_2- \theta_2^0) \to_D\sqrt{\tau(1-\tau)}[f_e(0)]^{-1} [\mathbb{E}(v_1v_1^\top)]^{-1/2} N(0,I_{d_2})$, which is the result of quantile regression for linear model. See koenker1978.
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
where $\widehat{e}_t=y_t-\sum_{j=1}^{p_1}\widehat\gamma_{1j}g_{1j} (\widehat{\theta}_{1j}^{\,\top} x_t) -\sum_{j=1}^{p_2}\widehat\gamma_{2j}g_{2j}(\widehat{\theta}_{2j}^{\,\top} z_t)$ for $1\leq t\leq n$, and $m=[n^{2+\varepsilon}]$ for some $\varepsilon>0$. Although both $a_1$ and $a_2$ are defined in terms of conditional expectation, they are the corresponding expectations as well, so we simply define their estimators as the corresponding sample averages. We then have the following corollary.
In the proof of $\widehat{a}_1\to_Pa_1$ we use the ergodicity of $\{\psi(e_t)\}$ that is implied from its martingale difference sequence property; similarly in the proof of $\widehat{a}_2\to_Pa_2$, we need the ergodicity of $\{\rho_m''(e_t)\}$. Certainly, the ergodicity of the two process is ensured by that of $\{e_t\}$, but we do not impose this condition in order to keep the condition as weak as possible.
In this subsection, we shall consider the case where all $g_{1j}(\cdot)$ are $I$-regular satisfying Assumption 3.4(ii). This assumption will change the asymptotic theory considerably. Note that in 3.4(ii), $\theta_{1j}^0= \theta_1^0, j=1, \cdots, p_1$, so all $I$-regular functions have the same stochastic trend $x_t^\top\theta_{1}^0$ as their argument. Moreover, the limit theory depends on the local time defined for scalar Brownian motion $W(r)$ as
which measures the sojourning time of $W(r)$ around $s$ during the time period $(0,p)$. In our context, the scalar Brownian motion is $B_{1}(r)\equiv (\theta_{1}^0)^\top B(r)$ where $B(r)$ is defined by (ref), and we denote its local time process by $L_1(p,s)$.
Another feature of the limit theory about this kind of regression functions is that we derive the asymptotics in a new coordinate systems first, then recover the limits in the original coordinate. This is the same as phillips2000, and dgd2016. To do so, using $\theta_{1}^0$ that is a unit vector, we construct an orthogonal matrix $P=(\theta_{1}^0, P_{2})_{d_1\times d_1}$ to rotate all the vectors of interest, including $\theta_{1}^0$, $x_t$ and $B(r)$. Hence, $P^\top \theta_{1}^0 =(1,0,\cdots,0)^\top$; other relevant rotated vectors are $P^\top x_t$ that have the first element $x_{1,t}\equiv (\theta_{1}^0)^\top x_t$ and sub-vector $x_{2,t}\equiv P_{2}^\top x_t$; due to continuous mapping theorem, after normalized by $\sqrt{n}$ they converge jointly as follows:
To proceed, denote $D_n=\text{diag}(\sqrt[4]{n}, \sqrt[4]{n}^3I_{d_1-1})$, $p_n=\sqrt{n}$ and $q_n=\sqrt[4]{n}$, and define
Here, $\theta_1, \theta_{2j}, \gamma_{1j}, \gamma_{2j}$ are generic parameters, from which $\alpha_1, \beta_{1j}, \alpha_{2j}, \beta_{2j}$ are centered and rescaled, and in particular $\alpha_1$ is the centered vector $\theta_{1}-\theta_{1}^0$ rotated by $P$ and rescaled by $D_n$.
The role that the rotation takes is for asymptotic analysis only, while in Monte Carlo experiments and empirical study the rotation is not necessary; more details can be seen in the relevant literature such as phillips2000, and dgd2016.
To make the notation compact, we also define $\Lambda\equiv(\Lambda_1^\top, \Lambda_2^\top)^\top$ where
We shall give the asymptotics of $\widehat\Lambda=(\widehat\Lambda_1^\top, \widehat\Lambda_2^\top)^\top$ where
that are defined in terms of $\widehat\theta_1, \widehat\theta_{2j}, \widehat\gamma_{1j}, \widehat\gamma_{2j}$ through the relationship in (ref).
To state the asymptotic properties of these estimators, we need to define several much more complicated vectors and matrices. First, let $x_{nt}\equiv D_n^{-1}P^\top x_t$, $\dot{g}(u)\equiv\sum_{j=1}^{p_1} \gamma_{1j}^0 \dot{g}_{1j}(u)$ and $Z_t\equiv(Z_1(x_{nt})^\top, Z_2(z_t)^\top)^\top$, where
where $p_n$ and $q_n$ are the same as in (ref).
Second, to show the limits of $\sum_{t=1}^nZ_tZ_t^\top$ and $\sum_{t=1}^n\psi(e_t)Z_t$ we need to define a matrix. Let $\mathcal{R}$ be a square matrix of order $(d_1+p_1)+p_2(d_2+1)$, standing for the limit of $\sum_{t=1}^nZ_tZ_t^\top$, that we divide into blocks $\mathcal{R}=(\mathcal{R}_{ij})_{2\times 2}$, conformably with the blocks in $Z_tZ_t^\top$, i.e. $Z_1Z_1^\top$, $Z_1Z_2^\top$, $Z_2Z_1^\top$ and $Z_2Z_2^\top$. Thus, $\mathcal{R}_{11}=(R_{ij}^{11})$ is a symmetric matrix of order $d_1+p_1$ that has elements,
where $R_{2:d_1, 2:d_1}^{11}$ denotes all elements $R_{i, j}^{11}$ with $i,j =2, \cdots, d_1$; other symbols, such as $R_{2:d_1, d_1+1}^{11}, \cdots$ and $R_{2:d_1, d_1+p_1}^{11}$, are explained similarly.
Note that $\mathcal{R}_{22}=\Sigma$, the same as given in Theorem 4.1, while $\mathcal{R}_{12}=\mathcal{R}_{12}^\top=0$ a zero matrix of order $(d_1+p_1)\times (d_2p_2+p_2)$. Thus, $\mathcal{R}$ is a diagonal block matrix.
It can be seen that there may be many zeros even in the block $\mathcal{R}_{11}$ if $g_{11}(u), \cdots, g_{1p_1}(u)$ are taken from an orthogonal sequence such that $\int g_{1j}(u)g_{1k}(u)du=0$ for $j\ne k$. This is most likely to be true because the orthogonality among $g_{11}(u), \cdots, g_{1p_1}(u)$ is the primary motivation for their choice.
The proof of the theorem is given in Appendix B. Similar to the comment right below Theorem 4.1, the consistency and convergence order of all estimators are implied by $\widehat{\Lambda}=O_P(1)$.
Due to diagonal block structure of $\mathcal{R}$, the assertion (ref) implies
as $n\to\infty$. It follows from (ref) that all estimators $\widehat\theta_{2j}$ and $\widehat\gamma_{2j}$, $j=1,\cdots, p_2$, have conventional $\sqrt{n}$-convergent rate, whereas the rates of $\widehat\theta_{1}$ and $\widehat\gamma_{1j}$, $j=1,\cdots, p_1$, affected by the unit root property of $x_t$ and $I$-regularity of functions $f_{1j}$, are distorted, similar to the results in dgd2016. This feature coincides with DL2018 since these models have additive form where these nonstationary and stationary variables are separated in different components.
In particular, $\widehat\gamma_{1j}$ have a very slow rate $n^{-1/4}$, because they are the coefficients of $g_{1j}(x_t^\top \theta_1^0)$, but $g_{1j}(x_t^\top \theta_1^0)$ attenuates to zero when $t$ gets large. This rate is the same as classical papers such as phillips1999, phillips2001. However, the rate of $\widehat\theta_{1}$ is blurred due to the involvement of $D_n$ and $P$. To find out its exact rate, denote by $r_{11}$ the left-top $d_1\times d_1$ submatrix of $\mathcal{R}^{-1}_{11}$. Thus,
The following corollary recovers the asymptotic distribution of $\widehat\theta$ from the above limit of its rotation.
Though indicated by (ref) that the coordinates of $\widehat\theta_1$ on $\theta_1^0$ and $P_2$ have rates $n^{-1/4}$ and $n^{-3/4}$, respectively, the estimator $\widehat\theta_1$ eventually has slower rate $n^{-1/4}$. This is the same as dgd2016 and one may find more detailed explanation therein. In addition, if we normalize $\widehat\theta_1$ to be $\widehat\theta_{1,unit}=\widehat\theta_1/\|\widehat\theta_1\|$, one may find $\widehat\theta_{1,unit}$ converges to $\theta_1^0$ with a quicker rate $n^{-3/4}$. Since this is exactly the same as dgd2016, we omit the derivation and refer the readers to the reference paper.
Similar to Corollary (ref), we may define $\widehat{a}_1$ and $\widehat{a}_2$ and then establish their consistency. This however is omitted due to similarity.
This subsection presents simulation experiments to evaluate the finite sample performances of the robust M estimators in the multiple index cointegration model. For space consideration, we only entertain two examples. The first example considers homogeneous cointegrating functions, while the second example studies the case with an integrable cointegrating function.
The time series $\{x_t, z_t\}$ used in the following two examples are generated as follows. $x_t=\rho_{1} x_{t-1}+\sigma_1 w_{t}$, $z_t=h(t/n)+v_t$, $h(\tau)=(\tau, \tau)^\top$, $v_t=\rho_{2} v_{t-1}+\sigma_2 \epsilon_{t}$, both $x_t$ and $v_t$ are bivariate autoregressive vector processes with $x_0=0$, $v_0=0$, $\rho_{1}=I_2$, $\rho_{2}=0.5\times I_2$, $\sigma_1={\rm diag}\{0.2,0.5\}$, $\sigma_2=I_2$, and $(w _{t}^\top, \epsilon_{t}^\top)$ is a vector of independent and normally distributed random variables with zero mean and identity covariance matrix.
{ \bf Example 5.1}. We first consider the following model:
where ${\bf\theta}^\top_1= (\theta_{11},\theta_{12})=(1,1)/\sqrt{2}$, ${\bf\theta}^\top_2=(\theta_{21},\theta_{22})=(1, 1)/\sqrt{2})$, ${\bf\theta}^\top_3=(\theta_{31},\theta_{32})=(1,1)/\sqrt{2}$, $\gamma_1=\gamma_2=2$, $\gamma_3=1$. The regression residual is set as $e_t=0.5\cdot u_t$ with $u_t$ independently generated according to four distributions: (D1) standard normal, $N(0,1)$; (D2) mixed normal, $0.9\cdot N(0,1)+0.1\cdot N(0,4)$; (D3) $t(2)$, $t$ distribution with 2 degrees of freedom; (D4) $t(1)$, standard Cauchy distribution.
{ \bf Example 5.2}. We next consider the following model:
where $\phi(\cdot)$ is the standard normal density function, ${\bf\theta}^\top_1=(\theta_{11},\theta_{12})=(1,1)/\sqrt{2}$, ${\bf\theta}^\top_2=(\theta_{21},\theta_{22})=(1,1)/\sqrt{2}$, $\gamma_1=2$, $\gamma_2=1$. $e_t$ is generated the same as in Example 5.1.
For demonstration, three loss functions are entertained in both examples: (L1) Huber's loss $\rho_c (e)=\frac{1}{2}e^2\cdot 1\{|e|\leq c\}+(\delta (|e|-\frac{1}{2} c^2) \cdot 1\{|e|> c\}$ with $c=1.25$; (L2) absolute errors loss $\rho(e)=|e|$, and (L3) quantile loss $\rho_\tau(e)=e\cdot (\tau-I\{e<0\}$ with $\tau=0.3$. The simulations are conducted for $n=50, 100, 200$ with 5,000 replications. For the quantile regression, the regression residual is normalized such that its $\tau$-th quantile is set as 0. To measure the estimation accuracy of the nonlinear least-squares estimates for the index parameters, we compute the bias, estimated standard deviations and mean squared errors for each element of the indices. For space limitation, we only report the mean squared errors in Tables (ref)-(ref) for Examples 5.1-5.2, respectively. Each table contains the estimation results under the three loss functions L1-L3, and four types of error distributions D1-D4.
The main findings from Table (ref) are summarized as follows. First, the mean squared errors are decreasing as the sample size $n$ increases, which indicate the consistency of the robust M estimators. Second, the index parameter estimates associated with the nonstationary variables enjoy super-rate of convergence, as shown in all the cases, in comparison with the standard $\sqrt{n}$ rate of the estimators associated with the stationary variable. In particular, the estimators of $\gamma_1$, $\theta_{11}$ and $\theta_{12}$ are close to $n^2$-consistent, while those of $\gamma_2$, $\theta_{21}$ and $\theta_{22}$ are close to $n^3$-consistent. These results are consistent with our theoretical developments. Finally, the three estimates considered are quite robust in terms of involving the different error distributions. It is found that the Huber's estimator is the most efficient, except for the case with the Cauchy ($t(1)$) errors. The least absolute error estimate becomes the most efficient one when the errors are Cauchy, which is consistent with the theory.
Turning to Table (ref) where the cointegrating function is integrable, there are a few new findings which are summarized as follows. First, the estimator of $\gamma_1$ is consistent but converges at a relatively slow rate. This is primarily due to the integrable nature of the cointegrating function. Second, the normalized estimators of $\theta_{11}$ and $\theta_{12}$ are converging at a relatively faster rate than that of $\gamma_1$ does. This finding corroborates the comment below Corollary 4.2 and that discovered in dgd2016.
{
}
To demonstrate the practical relevance of our proposed model over some natural competitors, we investigate its applicability in stock return prediction. It is common to use a linear predictive mean regression in the literature, which has led to considerable disagreements in the empirical findings as to whether stock returns are predictable or not (CT2008; WG2008). Recently, KAP2015a use a nonparametric mean regression, Lee2016, and FL2019 consider a quantile linear regression, and TLW2021 adopt a nonparametric quantile framework to re-investigate this important issue. We shall demonstrate the nonlinearity in return predictability through the use of the proposed multiple-index model with both stationary and nonstationary predictors.
The data sets to be used for return prediction are obtained from WG2008. They have been extensively used in the predictive literature, including the recent balanced predictive mean regression model by Ren et al. (2019), the linear quantile predictive regression by Lee2016, and FL2019, the linear prediction with cointegrated variables by kasy2020, the LASSO predictive regression of LSG2021, and the nonparametric quantile predictive regression of TLW2021, among others. The monthly data spans from January 1927 to December 2005. The dependent variable, excess stock return, is defined as the difference between the S&P 500 index return, including dividends and the one month Treasury bill rate. The nonstationary (persistent) predictors include dividend-price ($dp$), dividend-payout ratio ($de$), long term yield ($lty$), book to market ($bm$) ratios, T-bill rate ($tbl$), default yield spread ($dfy$), net equity expansion ($ntis$), earnings-price ($ep$), term spread ($tms$), while the stationary variables include default return spread ($dfr$), long term rate of return ($ltr$), stock variance ($svar$) and inflation ($infl$). The first-order serial correlations of these variables and their time series plot are available from, for example, LSG2021. For detailed description on each series and the data construction, please refer to WG2008.
To measure the performance of the predictive regression models under discussion, we use the pseudo out-of-sample $R^2$, defined as $$PR^2=\frac{\sum_{t=T_1+1}^{T_2} \rho(\widehat e_t)}{\sum_{t=T_1+1}^{T_2} \rho(\overline{e}_t)},$$ where $\widehat{e}_t$ and $\overline{e}_t$ are the respective out-of-sample prediction errors at time $t$ for a given model and that for the base model with only a constant predictor, for the forecasting sample spanning from $T_1$ to $T_2$. Four loss functions are entertained for $\rho(\cdot)$: (SE) squared errors loss $\rho(e)=e^2$; (AE) absolute errors loss $\rho(e)=|e|$, (HL) Huber's loss $\rho_c (e)=\frac{1}{2}e^2\cdot 1\{|e|\leq c\}+c\cdot (e-\frac{1}{2} c) \cdot 1\{|e|> c\}$ with $c=1.25$; and (QL) quantile loss $\rho_\tau (e)=e\cdot (\tau-I(e<0))$. Positive $PR^2$ indicates the better performance of the given model over the base model, and the larger the value is, the better the performance is.
For space limitation, we follow KAP2015a, Lee2016, and TLW2021 to consider two subsamples: (i) the period spanning from January 1927 to December 2005, where significant predictability has been discovered, and (ii) the tranquil period starting from January 1952 to December 2005, where mixed evidences are found for predictability. We consider a rolling-window scheme for the out-of-sample return prediction and let the rolling window in-sample size (RWS) be 120 (10 years), 240 (20 years), and 360 (30 years), respectively. For example, the forecast starts from January 1927 and ends at December 2005 for the first subsample when RWS equals 120.
In contrast to TLW2021, our proposed model allows for the use of multivariate nonstationary predictors together with those stationary ones to capture nonlinearity in the return predictability. There are numerous choices of combinations of the nonstationary and stationary predictors. For demonstration purposes, we use the stationary predictor $z=\{inf, svar\}$, and choose the nonstationary predictors as either Case A: $x=\{bm,lty\}$ or Case B: $x=\{dp,tms\}$. For the nonstationary index, we consider either a linear transformation $g_1(u)=\gamma u$ or an integrable transformation $g_1(u)=\ell(u)=\gamma_1 \cdot exp(-u^2)+\gamma_2 \cdot u \cdot exp(-u^2)$. Similarly, we use either $g_2(v)=\delta v$ or $g_2(v)=\ell(v)=\delta_1 \cdot exp(-v^2)+\delta_2 \cdot v \cdot exp(-v^2)$ for the stationary index. The choice for the nonlinear transformation is made based on the (normalized) Hermite polynomials in the sieve literature where Hermite polynomials form a basis in the related function space. This amounts to four possible combinations. The extensive analysis with a comprehensive set of predictors noted above and all possible specifications of the nonlinear link functions remains an interesting but a challenging work, which is left as future study.
{
}
{
}
Table (ref) collects the pseudo out-of-sample $R^2$ for excess return prediction with the above two sets of selected predictors and specifications of the link functions, respectively, and for the loss functions SE, AE and HL. For each loss function and each RWS, the combination of the link functions that produces the largest $PR^2$ is bolded. There are several interesting findings. First, the proposed models demonstrate significant return predictability compared to the baseline model with only an intercept. For Case A with $x=\{bm,lty\}$, the linear predictive regression outperform other specifications with nonlinearities for all three loss functions for the first sample period January 1927 to December 2005 (Panel A). When it comes to the second sample period in Panel B, nonlinear predictability is discovered, even though the predictability is seen to decline, as that discovered by CY2006. For Case B with $x=\{dp,tms\}$, we find no predictability from the linear predictive model for both sample periods. The pseudo out-of-sample $R^2$ can reach 10.8%, 4.35%, 11.2% for the first sample period under the squared error loss, the least absolute error loss, and the Huber loss, respectively. The nonlinear predictability is also seen to decline in the “tranquil” period, but remains significant with the pseudo out-of-sample $R^2$ being as large as 10.2%, when both indices enter the model in the nonlinear manner. This finding is in contrast to the weak nonlinear predictability discovered by KAP2015a with only a single nonstationary predictor. To summarize, significant return predictability has been revealed in the nonlinear predictive regression with both stationary and nonstationary indices we introduced.
As the recent literature has witnessed a growing interest in investigating the return prediction in quantiles (Lee, 2016; Fan and Lee, 2019; Tu, et al., 2021), we next study the performance of our proposed model in predicting the return quantiles. We consider the quantile level $\tau=0.05,0.1,0.2,\ldots,0.9,0.95$. The analysis will only focus on the nonlinear effect of the nonstationary index as discovered from Tables (ref). Therefore, we only present the results for $x=\{dp,tms\}$ and $z=\{inf, svar\}$ in Table (ref) at all the quantile levels specified above, for the two subsamples mentioned earlier. The specifications for transformation on the stationary and nonstationary components are the same as above. From Table (ref), it is observed that linear quantile predictability seems to only exist for the lower quantile levels, such as those at $\tau=0.05,0.1$ and $0.2$. However, nonlinear quantile predictability prevails at all quantile levels. Even at the lower quantiles, the proposed nonlinear index model outperforms the linear model except for a very few cases. In addition, it seems that the nonlinear predictability becomes stronger as the quantile level moves to the two ends, especially for the first subsample.
The paper has studied a class of additive single index models with diverse regressors including I(1) process, time trend and $\alpha$-mixing stationary process. To deliver the asymptotic theory of robust estimation, a generalized function approach has been proposed that has independent interest, as shown in Section two with a simple linear model. This approach associated with Theorem 2.1 has answered a call from the literature to establish the convergence rate of regular sequences to nonsmooth loss functions. Thus, we therefore believe that Theorem 2.1 fills a gap in the literature.
The generalized function approach is then applied in general robust estimation for a class of additive single--index cointegrating time series models where the objective function is constructed with possible nonsmooth loss, such as LAD, quantile loss and Huber's loss; the corresponding asymptotic theory has been established, respectively, according to H-regular and I-regular classes of the regression functional forms involving I(1) processes. Monte Carlo simulations are conducted to verify the performance of estimators proposed in finite sample situations; and an empirical study on stock returns has also been implemented to exhibit both the relevance and applicability of the proposed model and estimation method developed in the paper.
Dong would like to thank the financial support from National Natural Science Foundation of China (Grant 72073143) and the Fundamental Research Funds for the Central Universities, Zhongnan University of Economics and Law (2722022EG001); Gao acknowledges financial support from the Australian Research Council Discovery Grants Program under Grant Number: DP200102769; Peng acknowledges the Australian Research Council Discovery Grants Program for its financial support under Grant Number DP210100476; and Tu (Corresponding Author) would like to thank support from National Natural Science Foundation of China (Grant 72073002, 12026607, 92046021), the Center for Statistical Science at Peking University, and Key Laboratory of Mathematical Economics and Quantitative Finance (Peking University), Ministry of Education.
{
}