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.
66,634 characters · 7 sections · 72 citation commands
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 (ref) $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 ((ref)) 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 ((ref)) has not been proposed and studied in the relevant literature. For the univariate regressor setting, chang2001 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. chang2003 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, dgd2016 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 ((ref)). DL2018 and dlp2021 extend chang2001, 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 huber1964, 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 vaart1996).
Using M-estimation, the unknown parameters $(\theta'_0, \beta'_0)'$ 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,
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 (chenjia2010). 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, chenxh2014).
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 stankovic1996).
phillips1991, phillips1995 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, walsh2008 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 walter1979 and walter1981, 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 kanwal1983). 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 (ref) 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 (ref), 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.
{\bf Examples}.\ \ The function given by
is called Heaviside function. Its value at zero usually is defined as $1/2$, while sometimes it is taken as $c$ for $0<c<1$, written as $H_c(x)$. Simply by the definition of the derivative of generalized functions, $H'(x)= \delta(x)$ the so-called Dirac delta function.
(a) Let $\rho(u)=|u|$ (LAD loss). Again, by the definition of the derivative of generalized functions, $\rho'(u)= -H(-u)+H(u)$, and thus $\rho''(u)=\delta(-x) +\delta(x)= 2\delta(x)$, where we use $\delta(-u)=\delta(u)$ given in kanwal1983.
(b) Let $\rho(u)=u(\tau-I(u<0))$ for some $0<\tau <1$ (Quantile loss). Similarly, we have $\rho'(u)=\tau-I(u<0)=\tau H(u)+(\tau-1) H(-u)$. Hence, $\rho''(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''(u)=I(|u|\le 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'(u)=2|\tau-I(u<0)|u$, whereas $\rho'(u)$ is not smooth at $u=0$. By definition, $\rho''(u)=2|\tau-I(u<0)|$.
(e) For $L_p$ norm loss $\rho(u)=|u|^p$ with $p>1$, we have $\rho'(u)=p|u|^{p-1}\text{sgn}(u)$. However, if $1<p<2$, as an ordinary function $\rho''(u)$ does not exist at $u=0$. In generalized function sense, $\rho''(u)=p(p-1)|u|^{p-2}$, while the meaning of generalized function $|u|^{p-2}$ is discussed in gelfand1964 and kanwal1983 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$).
{\bf Assumption A}\ \
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:=\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 Brownian motion on $[0,1]$ with covariance $AA^\top$.
Assumption B.
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'(e_t)$, $\xi_{2t} = \rho'(e_t)z_t$ and $\xi_{3t} = \rho''(e_t)$.
Assumption C. \
which is a Brownian motion of $(d+1)$ dimension, where $U_\rho(r), V_\rho(r)$ and $W(r)$ have variance and covariance, $a_1$, $a_1 \int {\mathbb E}(h(r, v_1)h(r, v_1)^\top)dr$ and $I_{d_1}$, respectively.
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'(\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 (C.2) can be taken as $\mathcal{F}_{t}=\sigma(e_1,\cdots, e_t; x_s, z_s, s\le t+1)$ that is general than series independence between $e_t$ and $(x_t,z_t)$. When $\rho'(\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 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. 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|\le c)>0$ in exogenous situation.
When $\rho'(\cdot)$ and $\rho''(\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 hexuming2000 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 phillips2001, xiao2009a and qiying2018). Nevertheless, here the subscript $\rho$ indicates that the joint convergence relies on the loss function. See, 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 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 phillips1999,phillips2001. 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 phillips1999, phillips2001.
{\bf Assumption D}\; \; Let both $g_1(u)$ and $g_2(v)$ be continuously differentiable up to the second order. Suppose further that
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, phillips1991, phillips1995).
{\bf 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 phillips2001. 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'(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$ 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 [2f_e(0)]^{-1}\left[\int_0^1B(r)B(r)^\top dr\right]^{-1} \int_0^1 B(r)dU(r),$$ that is the result of phillips1995 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_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$ 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^1B(r)B(r)^\top dr]^{-1} \int_0^1 B(r)dU(r)$ that is the result of xiao2009a 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_1v_1^\top)]^{-1/2} N(0,I_{d_2})$, which is the result of quantile regression for linear model. See koenker1978. $\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'(\widehat{e}_t)]^2$ and $\widehat{a}_2=\frac{1}{n}\sum_{t=1}^n \rho''(\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 phillips1995).
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 (ref).
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 (ref) 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 (ref) 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 pollard1991 and koenker1978). 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, phillips2001 and dgd2016.
The convergence of $\eqref{int1}$ is similar to Theorem 3.1 in dgd2016 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 phillips2001 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 phillips2000.
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
The following corollary recovers the asymptotic distribution of $\widehat\theta$ from the limit of its rotation.
Though indicated by (ref) that the coordinates of $\widehat\theta$ on $\theta_0$ and $P_1$ have rates $\sqrt[4]{n}$ and $\sqrt[4]{n}^3$, respectively, the estimator $\widehat\theta$ eventually has slower rate $\sqrt[4]{n}$. This is the same as dgd2016 and one may find more detailed explanation therein. Interestingly notice that since the model in dgd2016 is nonparametric, $\theta_0$ is identified up to its direction, so the condition $\|\theta_0\|=1$ is imposed. Thus, the authors normalize the estimator $\widehat\theta$ to be a unit vector and find that the normalized estimator has a much faster rate than $\sqrt[4]{n}$. By contrast, no identification condition is needed in this paper as the model is parametrically nonlinear. Consequently we would not be able to achieve any rates faster than what we have obtained. However, if one could know the length of $\theta_0$, say $\|\theta_0\|=q$, the estimator $q\;\widehat\theta/\|\widehat\theta\|$ would converge to $\theta_0$ with a rate of an order $\sqrt[4]{n}^3$, basically because the normalization stretches $\widehat\theta$ directly on the unit circle with radius $q$ (see, Theorem 3.2 in dgd2016, for example). Moreover, one may have the estimates of $a_1$ and $a_2$ similarly to Corollary (ref) and that of the local time in the limit similarly to dgd2016. The details are omitted due to the similarity.
Before we conclude this section, we briefly discuss another important issue related to model (ref). If the model has conditional heteroscedasticity, such as $e_t=\sigma(\vartheta_0^\top x_t, \pi_0^\top z_t) \varepsilon_t$, where $\varepsilon_t$ is independent of $(z_t,x_t)$, and if $\sigma(u,v)=\sigma_1(u)\sigma_2(v)$ is known and separable, we will have
where $\zeta_t$ is the centralized version of $\log(\varepsilon^2_t)$ with $\mu_0=\mathbb{E}[\log(\varepsilon^2_t)]$.
In this case, the proposed estimation method is readily applicable for us to derive the corresponding M-estimators of $\vartheta_0$ and $\pi_0$ after we replace $e_t$ by $\widehat{e}_t=y_t-g_1(\widehat{\theta}^{\;\top} x_t)-g_2(\widehat{\beta}^{\;\top} z_t)$. Under some same assumptions on $\log(\sigma_1^2(u))$ and $\log(\sigma_2^2(u))$ as in Assumption D, we may establish asymptotic properties corresponding to the above theorems and corollaries for the estimators of $\vartheta_0$ and $\pi_0$.
To demonstrate the practical relevance and superiority of our proposed model and estimation method over some natural competitors, we investigate its applicability in stock return predictability. It is common to use a linear predictive mean regression in the literature, which has led to considerable disagreements in the empirical finding 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)$: (L1) squared errors loss $\rho(e)=e^2$; (L2) absolute errors loss $\rho(e)=|e|$, (L3) Huber's loss $\rho_\delta (e)=\frac{1}{2}e^2\cdot 1\{|e|\leq \delta\}+\delta\cdot (e-\frac{1}{2} \delta) \cdot 1\{|e|> \delta\}$ with $\delta=1.25$; and (L4) 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 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 $x=\{bm,lty\}$ or $x=\{dp,tms\}$. When $x=\{bm,lty\}$, we consider an exponential link function $g_1(u)=u\cdot e^{-0.4u^2}$ and a linear link $g_2(v)=v$ outside of the nonstationary (resp. stationary) and stationary (resp. nonstationary) index variables, respectively. This amounts to four possible combinations. When $x=\{dp,tms\}$, we alternatively consider the normal cumulative distribution function $\Phi$ to capture possible nonlinearity. The logistic distribution function has also been experimented in this case and makes little difference to subsequent main findings, the results of which are therefore not reported below. 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.
Tables (ref) and (ref) collect 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, for the loss functions specified in (L1)-(L3). For each loss function and each RWS, the combination of the link function that produces the largest $PR^2$ is bolded. There are several interesting findings. First, noticeable nonlinearity exists in return prediction, for both choices of the nonstationary predictors, and for both subsamples considered. It is especially found that for the squared error loss and the Huber's loss, the linear specification is always worse than the nonlinear alternatives. When $x=\{bm,lty\}$, nonlinearity has been discovered for both the stationary index and the nonstationary component. However, when $x=\{dp,tms\}$, only the nonstationary index presents nonlinearity in return prediction. Second, nonlinear predictability is found in forecasting return when linear predictability is not. The linear predictive regression (nearly) outperforms the constant model without any predictor when $x=\{bm,lty\}$, but the former underperforms the latter when $x=\{dp,tms\}$. Although the linear model does not reveal predictability in this case, it is uniformly observed that the nonlinear link function specification improves over the constant model for both choices of subsamples. Third, the conventional “tranquil” period 1952-2005 reveals strong nonlinear predictability, which is quite comparable to the 1927-2005 subsample. This finding is in contrast to the declining linear predictability discovered by CY2006, or the weak nonlinear predictability with only a single nonstationary predictor by KAP2015a. 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) and (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 nonlinear specification for the nonstationary component includes not only the normal distribution function, but also the popularly used logistic distribution function (Logit), as a robustness check. 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. To compare the normal with logistic distributions, the latter seems to enjoy better performance for both sample periods. In addition, it seems that the nonlinear predictability becomes stronger as the quantile level moves to the two ends, especially for the first subsample.
{
}
We consider the following data generating process:
where $x_t=\rho_{n1} x_{t-1}+\sigma_1 w_{t}$, $z_t=h(t/n)+v_t$, $h(\tau)=\tau$, $v_t=\rho_{n,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_{n1}=I_2$, $\rho_{n2}=0.5\times I_2$, $\sigma_1=diag\{0.2,0.5\}$, $\sigma_2=I_2$, and $(w _{t}^\top,\epsilon_{t}^\top)$ is a series of independent four dimensional normal random vector with zero mean and identity covariance matrix, $\theta=(1,0)^\top$, $\beta=(2,1)^\top$. The regression residual $e_t=0.1\cdot u_t$ with $u_t$ independently generated according to four distributions: (D1) standard normal; (D2) mixed normal $0.9\cdot N(0,1)+0.1\cdot N(0,4)$; (D3) $t$ distribution with 2 degrees of freedom; (D4) standard Cauchy distribution.
For the bivariate nonlinear function $g(\cdot,\cdot)$, we consider the following designs:
where $\phi(\cdot)$ is the standard normal density function, $\Phi(\cdot)$ is the standard normal distribution function. Three loss functions are entertained: (L1) squared errors loss $\rho(e)=e^2$; (L2) absolute errors loss $\rho(e)=|e|$, and (L3) Huber's loss $\rho_\delta (e)=\frac{1}{2}e^2\cdot 1\{|e|\leq \delta\}+\delta\cdot (e-\frac{1}{2} \delta) \cdot 1\{|e|> \delta\}$ with $\delta=1.25$. The simulations are conducted for $n=100, 200, 400$ with 5,000 replications.
To measure the estimation accuracy of the nonlinear least-square 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 M1-M4, 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 are summarized as follows. First, the mean squared errors are decreasing as the sample size $n$ increases, except when the errors are generated according to the Cauchy distribution and when the least squares loss is implemented. Under Cauchy errors, the least-square estimator is known to be inconsistent. However, the median estimator and the Huber's estimator remain consistent. As a result, we shall exclude this scenario when we further comment on the estimator performance below. Second, the index parameter estimate of the nonstationary variables enjoys super-rate of convergence, as shown in all the cases, as compared to the standard $\sqrt{n}$ rate for the stationary index estimator. Finally, the three estimates considered are quite competitive in terms of the different error distributions. It is found that, as expected, the least squares estimator is the most efficient estimator when the errors are drawn from normal distributions. The Huber's estimator becomes the most efficient, in general, for the mixed normal errors and $t(2)$ errors. The least absolute error estimate is the most efficient one when the errors are Cauchy.
{\samepage {
}}
In order to cater for the practical usefulness this paper proposes a class of multiple index time series parametric models that accommodate time trend as well as both stationary and nonstationary vectors. An $M$-estimation approach is adopted where the loss functions can be, but not limited to, six popular ones, such as the squared loss, LAD, Huber's loss, quantile loss, $L_p$ and expectile loss. Meanwhile, two categories of link functions are investigated due to the different behaviour of functions of nonstationary vector variables, that is, $I$-regular and $H$-regular classes in the related literature. Accordingly, our asymptotic theory dwells on two categories of estimators and the rates of convergence are discussed under different classes of loss functions. Moreover, these models are used in the analysis of predictability where we find that our models are competitive comparing with the literature in terms of predictability. Finally, we conduct Monte Carlo simulations that reveal the satisfactory performance of our estimators proposed in finite sample situations.
Dong would like to thank the financial support from National Natural Science Foundation of China (Grant 72073143); Gao acknowledges financial support from the Australian Research Council Discovery Grants Program under Grant Number: DP170104421; Peng acknowledges the Australian Research Council Discovery Grants Program for its financial support under Grant Number DP210100476; and Tu 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.
{
}
{