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.
77,157 characters · 11 sections · 77 citation commands
Subgeometrically ergodic autoregressions
Markov chain theory and the notion of geometric ergodicity have become standard tools in econometrics and statistics when analyzing the stationarity and ergodicity of nonlinear autoregressions or other nonlinear time series models. A detailed discussion of the relevant Markov chain theory will be given in Section (ref). For now, consider a Markov chain $X_{t}$ ($t=0,1,2,\ldots$) on the state space $\mathsf{X}$ and initialized from $X_{0}$ following some initial distribution\textcolor{red}{ }(that is not necessarily the stationary distribution). Geometric ergodicity of $X_{t}$ entails that the $n$-step probability measures $P^{n}(x\,;\,\cdot)=\Pr\left(X_{n}\in\cdot\,|\,X_{0}=x\right)$ converge in total variation norm $\left\Vert \cdot\right\Vert _{TV}$ to the stationary probability measure $\pi$ at rate $r^{n}$ (for some $r>1$), that is,
(the definition of $\left\Vert \cdot\right\Vert _{TV}$ and a formulation of ((ref)) using a more general norm are given in Section (ref)). A common and convenient way to establish geometric ergodicity involves the verification of a so-called drift condition. Useful implications obtained with this approach include the existence of a stationary probability distribution $\pi$ of $X_{t}$ as well as the geometric $\beta$-mixing of $X_{t}$. (For a definition of $\beta$-mixing, see doukhan1994mixing or bradley2007introduction.)\textcolor{red}{ }The authoritative and classic reference to Markov chain theory is the monograph of meyn1993markov,meyn2009markov. Recent papers establishing geometric ergodicity of different nonlinear time series models include francq2006mixing, ling2007double, meitz2008ergodicity, and \citet*{fokianos2009poisson}, among others.
In this paper we consider autoregressions that may exhibit rather arbitrary (stationary, unit root, explosive, nonlinear, etc.) behavior for moderate values of the observed series and that behave similarly to a unit root process for large values of the observed series. What this exactly means will be clarified shortly, but first we would like to emphasize that the autoregressions we consider will not necessarily be geometrically ergodic. Under appropriate conditions they will, nevertheless, satisfy a weaker form of so-called subgeometric ergodicity. A Markov chain is said to be subgeometrically ergodic when the convergence in ((ref)) takes place at a rate $r(n)$ slower than geometric, that is,
In the geometric case $r(n)=r^{n}$ with $r>1$ or, equivalently, $r(n)=e^{cn}$ with $c>0$. Examples of rates slower than geometric include subexponential rates (say, $r(n)=e^{cn^{\gamma}}$ with $c>0$ and $\gamma\in(0,1)$) and polynomial rates (say, $r(n)=(1+n)^{\beta}$ with $\beta>0$). For an up-to-date treatment of subgeometric ergodicity we refer to Chapters 16 and 17 of \citet*{douc2018markov} (further references will be given below).
As will be discussed in Section 2, subgeometric ergodicity can conveniently be established by verifying a suitably formulated drift condition and useful implications analogous to those in the case of geometric ergodicity again follow. In particular, the existence of a stationary probability distribution $\pi$ of $X_{t}$ as well as the finiteness of certain moments are obtained. Moreover, in a companion paper meitz2019subgemix we show that subgeometric ergodicity implies $\beta$-mixing with subgeometrically decaying mixing coefficients. Subgeometric ergodicity therefore allows one to use limit theorems developed for $\beta$-mixing processes.
The main aims of this paper are to establish subgeometric ergodicity of certain higher-order nonlinear autoregressions and to illustrate the potential of the concept of subgeometric ergodicity for nonlinear time series models. To facilitate discussion, first consider a simple special case at an informal level. Specifically, consider the univariate first-order nonlinear autoregressive model
where the error term $\varepsilon_{t}$ is a sequence of independent and identically distributed (IID) zero-mean random variables and $g$ is a real-valued function. For now, assume that $g$ is such that
and that $g(x)$ is bounded for $\left|x\right|\leq M_{0}$. A concrete example where ((ref)) can be easily verified is
The model defined in ((ref)) can be thought of as first-order autoregression with a time-varying autoregressive coefficient. For large values of $\left|y_{t-1}\right|$ the autoregressive coefficient takes values that are close to one and the generation mechanism of $y_{t}$ is close to a random walk while for small values of $\left|y_{t-1}\right|$ the autoregressive coefficient is close to $1-r_{0}$ and, for $r_{0}$ not very close to zero, $y_{t}$ is generated from a less persistent stationary autoregressive process. Overall, the generation mechanism of $y_{t}$ fluctuates between these two borderline cases. Simulated examples in Section 5 demonstrate that processes of the type described in ((ref)) can exhibit behavior close to a random walk for rather long times before returning to a less persistent regime.
When $\rho\geq1$ the model defined by equation ((ref)) can be viewed as a special case of the model
where the function $\tilde{g}$ is bounded (but not constant).\footnote{This model belongs to a class of models referred to as “random-walk-type Markov chains” by jarner2003necessary; see particularly equation (3) of their paper.} In Section (ref), a higher-order version of equation ((ref)) (without assuming boundedness) is used as a starting point of the formulation of our general model. Our main results in Section (ref) show that, depending on the assumptions made, either geometric, subexponential, or polynomial ergodicity is obtained.
The preceding discussion illustrates what kind of behavior the autoregressions we consider may exhibit for large values of the observed series. However, it should be emphasized that inequality ((ref)) restricts the regression function $g$ only for large values of its argument. As long as the assumed boundedness condition imposed on the function $g$ is satisfied, no restrictions are required when the process evolves in the vicinity of the origin. Allowing unit root type behavior for large (absolute) values of the observed series is the main feature which distinguishes the models we consider from most previous nonlinear autoregressions where stationary behavior is related to large (absolute) values of the process.\footnote{For instance, lu1998geometric, gourieroux2006stochastic, and bec2008acr, among others, establish geometric ergodicity (and thus the existence of a stationary distribution) for autoregressions whose behavior approaches stationarity when the process moves away from the origin while in the vicinity of the origin its behavior can be rather arbitrary.}
Previously results on subgeometric ergodicity of nonlinear autoregressions have been obtained in the probability literature by tuominen1994subgeometric, veretennikov2000polynomial, fort2003polynomial, douc2004practical, \nocite{klokov2004sub,klokov2005subexponential}Klokov and Veretennikov (2004, 2005), and klokov2007lower, among others (further discussion on these and some related papers will be provided in Section 4). To our knowledge, all of the previous results concern only first-order models. We contribute to this literature by obtaining results for more general higher-order autoregressions. This is achieved using techniques similar to those in the aforementioned papers, especially in fort2003polynomial and douc2004practical. Depending on the assumptions imposed on the moments of the error term, the resulting rate of ergodicity is either geometric or subexponential or polynomial.
The rest of the paper is organized as follows. Section (ref) contains basic concepts of Markov chains and summarizes existing results on subgeometric ergodicity. Section (ref) introduces the nonlinear autoregressive model considered and states the assumptions used to obtain the results of the paper. The main results on subexponential and polynomial ergodicity are given in Section (ref). In Section 5 we provide examples of our general model. Section 6 concludes. All proofs are collected in an Appendix and a Supplementary Appendix.
Finally, a few notational conventions are given. The minimum (maximum) of the real numbers $x$ and $y$ is denoted by $x\wedge y$ ($x\lor y$), and $L$ and $\Delta$ signify the lag operator and the difference operator, respectively (so that $\Delta x{}_{t}=(1-L)x_{t}=x_{t}-x_{t-1}$). The notation $\boldsymbol{1}_{S}(x)$ is used for the indicator function which takes the value one when $x$ belongs to the set $S$ and zero elsewhere, and $\left|\cdot\right|$ is used for both an absolute value and Euclidean norm. Furthermore, $\text{\textbf{0}}{}_{k}$ denotes a $k\times1$ vector of zeros and $\boldsymbol{\iota}_{k}=(1,0,\ldots,0)$ ($k\times1$).
In this section we discuss basic concepts of Markov chains needed to obtain our results. More comprehensive discussions can be found in meyn2009markov and douc2018markov. Let $X_{t}$ ($t=0,1,2,\ldots$) be a Markov chain on a general measurable state space $(\mathsf{X},\mathcal{B}(\mathsf{X}))$ (with $\mathcal{B}(\mathsf{X})$ the Borel $\sigma$-algebra) and let $P^{n}(x\,;\,A)=\Pr(X_{n}\in A\mid X_{0}=x)$ signify its $n$-step transition probability measure. As in fort2003polynomial and douc2004practical our goal is to establish the convergence \textemdash in a suitably defined norm and at rate $r(n)$ \textemdash of the $n$-step probability measures $P^{n}(x\,;\,\cdot)$ to the stationary distribution $\pi$. To this end, let $f:\mathsf{X}\rightarrow[1,\infty)$ be an arbitrary fixed measurable function and, for any signed measure $\mu$, define the $f$-norm $\left\Vert \mu\right\Vert _{f}$ as
where $\mu(f_{0})=\int_{x\in\mathsf{X}}f_{0}(x)\mu(dx)$ (and the supremum in ((ref)) runs over all measurable functions $f_{0}:\mathsf{X}\to\mathbb{R}$ such that $\left|f_{0}(x)\right|\leq f(x)$ for all $x\in\mathsf{X}$). When $f\equiv1$, the $f$-norm $\left\Vert \mu\right\Vert _{f}$ reduces to the total variation norm $\left\Vert \mu\right\Vert _{TV}=\sup_{f_{0}:\left|f_{0}\right|\leq1}\left|\mu(f_{0})\right|$ used in ((ref)) and ((ref)).
We aim to establish that the $n$-step probability measures $P^{n}(x\,;\,\cdot)$ converge in $f$-norm and at rate $r(n)$ to the stationary probability measure $\pi$ satisfying $\pi(f)<\infty$, that is, that
\textcolor{red}{\footnotetext{That is, the convergence in ((ref)) is required to hold for all $x\in\mathsf{X}$ except for those $x$ in a set that has probability zero with respect to the stationary measure $\pi$.}}If ((ref)) holds we say that the Markov chain $X_{t}$ is ($f,r$)-ergodic; this implicitly entails the existence of $\pi$ as well as certain moments as $\pi(f)<\infty$. (For instance, if $\mathsf{X}=\mathbb{R}$ and $f(x)=1+x^{2}$, then ($f,r$)-ergodicity implies that the stationary distribution of $X_{t}$ has finite second moments.)
In the probability literature the preceding definition of ($f,r$)-ergodicity is standard. However, an equivalent and more transparent formulation is obtained by replacing equation ((ref)) with
\[ \lim_{n\to\infty}r(n)\sup_{f_{0}:\left|f_{0}\right|\leq f}\left|E[f_{0}(X_{n})\mid X_{0}=x]-\pi(f_{0})\right|=0\qquad\text{for }\pi\text{-almost all }x\in\mathsf{X} \] (see tuominen1994subgeometric). For instance, if $f(x)=1+\left|x\right|$ the above equation shows that, for almost any initial value $x,$ the conditional expectation $E[X_{n}\mid X_{0}=x]$ converges to $\int_{x\in\mathsf{X}}x\thinspace\pi(dx)$, the expectation of the stationary distribution of $X_{t}$, and the rate of the convergence is given by $r(n)$.
Most of the recent ergodicity results obtained for nonlinear autoregressions have established geometric ergodicity so that the rate of convergence in ((ref)) is given by $r(n)=r^{n}$, $r>1$. The subgeometric rate functions we consider are defined as follows (cf., e.g., nummelin1983rate and douc2004practical). Let $\Lambda_{0}$ be the set of positive nondecreasing functions $r_{0}\,:\,\mathbb{N}\rightarrow[1,\infty)$ such that $\ln[r_{0}(n)]/n$ decreases to zero as $n\rightarrow\infty$. The class of subgeometric rate functions, denoted by $\Lambda$, consists of positive functions $r\,:\,\mathbb{N}\rightarrow(0,\infty)$ for which there exists some $r_{0}\in\Lambda_{0}$ such that \[ 0<\liminf_{n\rightarrow\infty}\frac{r(n)}{r_{0}(n)}\leq\limsup_{n\rightarrow\infty}\frac{r(n)}{r_{0}(n)}<\infty. \] Typical examples are obtained of rate functions $r$ for which these inequalities hold with (for notational convenience, we set $\ln(0)=0$) \[ r_{0}(n)=(1+\ln(n))^{\alpha}\,\cdot\,(1+n)^{\beta}\,\cdot\,e^{cn^{\gamma}},\qquad\alpha,\beta,c\geq0,\,\gamma\in(0,1). \] The rate function $r_{0}(n)$ is called subexponential when $c>0$, polynomial when $c=0$ and $\beta>0$, and logarithmic when $\beta=c=0$ and $\alpha>0$. douc2004practical consider subexponential convergence rates whereas fort2003polynomial consider polynomial convergence rates in model ((ref)) (see also the related references mentioned in these papers).
The proofs of our results make use of the following condition adapted from douc2018markov.\footnote{A somewhat more general version which allows $V$ to be extended-real-valued (i.e., $V\,:\,\mathsf{X}\rightarrow[1,\infty]$) is given in douc2004practical.}
Condition D. There exist a measurable function $V\,:\,\mathsf{X}\rightarrow[1,\infty)$, a concave increasing continuously differentiable function $\phi\,:\,[1,\infty)\rightarrow(0,\infty)$, a measurable set $C$, and a finite constant $b$ such that
Conditions of this kind are known as drift conditions; when $\phi(v)=\lambda v$ for some $\lambda>0$ the so-called Foster-Lyapunov drift condition used to establish geometric ergodicity is obtained. For ease of discussion and reference, the following theorem summarizes geometric, subexponential, and polynomial ergodicity results that can be obtained using Condition D. (For the definitions of irreducibility, aperiodicity, and petite sets appearing in the theorem we refer the reader to meyn2009markov.)
In the geometric case the result of Theorem 1 is given in meyn2009markov, and in the subexponential and polynomial cases the result can be obtained from douc2004practical; some further details are provided in the proof of Theorem 1 in the Supplementary Appendix. Note that in the subexponential case choosing $v_{0}$ sufficiently large ensures the concavity of $\phi$ required in Condition D and also that (in the subexponential case) results with a faster rate of convergence and/or larger $f$-norm could be obtained at the expense of more complex notation; see douc2004practical.
An essential feature of the subgeometric ergodicity results in Theorem 1 is that there is a trade-off between the rate of convergence and the size of the $f$-norm; in Theorem 1 the choice of $\delta$ reflects this. If a fast rate of convergence is desired one has to accept a small $f$-norm (recall from ((ref)) that the size of the $f$-norm is directly proportional to the order of finite moments the stationary distribution is guaranteed to have). For instance, in the polynomial case choosing $\delta=1/(1-\alpha)$ gives the fastest rate of convergence and with this choice the $f$-norm reduces to the total variation norm (so that $f\equiv1$); the extreme case $\alpha=0$ results in $r(n)\equiv1$ and standard ergodicity. In the subexponential case values of $\delta$ that are close to zero (one) correspond to small (large) $f$-norms.
It is also worth noting that Condition D is only sufficient, not necessary, for ($f,r$)-ergodicity. It is therefore possible that with another drift condition (not necessarily a special case of Condition D) a better rate function could be obtained, but presumably at the cost of a smaller norm. Being able to obtain necessary conditions for particular subgeometric ergodicity rates would be of interest but we will not pursue this issue. Necessary conditions for geometric and polynomial ergodicity in the context of random-walk-type Markov chains (see ((ref))) are given in jarner2003necessary (for an application of this result to a threshold autoregressive model, see meitz2019subgemix).
As already indicated in the Introduction, the ergodicity results of Theorem 1 imply results on $\beta$-mixing or, more specifically, on convergence rates of $\beta$-mixing coefficients $\beta(n)$ ($n=1,2,\ldots$) (for a definition of $\beta(n)$ and properties of $\beta$-mixing, see doukhan1994mixing, bradley2007introduction, or meitz2019subgemix). To illustrate this point, let $\mu$ signify the distribution of $X_{0}$, the initial value of the Markov chain $X_{t}$, and assume that$\int_{x\in\mathsf{X}}V(x)\mu(dx)<\infty$ (with $V$ as in Theorem 1). Then, using Theorems 1 and 2 of meitz2019subgemix the three cases in Theorem 1 imply the following convergence rates for $\beta$-mixing coefficients (here $c$ and $\alpha$ are as in Theorem 1):\footnote{See also the discussion following Theorem 2 of meitz2019subgemix, and note that their Theorem 2(e) is also used to obtain the subexponential rate shown here.}
Thus, the convergence rates of the $\beta$-mixing coefficients are qualitatively similar to the fastest convergence rates of ergodicity obtained in Theorem 1 (as indicated above, a slight improvement can be achieved in the subexponential case). These results, combined with the fact that the $(f,r)$-ergodicity given in Theorem 1 implies finiteness of moments, make possible to use limit theorems developed for $\beta$-mixing processes (and also for $\alpha$-mixing processes because $\beta$-mixing is known to imply $\alpha$-mixing).
We now introduce a higher-order generalization of the model discussed in the Introduction. Suppose the process $y_{t}$ ($t=1,2,\ldots$) is generated by
where $\tilde{g}$ is a real-valued function, the error term $\varepsilon_{t}$ is a sequence of IID random variables, and exactly one of the roots of the polynomial $\varphi(z)=1-\varphi_{1}z-\cdots-\varphi_{p}z^{p}$ is equal to unity and (when $p\geq2$) all others lie outside the unit circle. Thus, the regression function of the model has a linear part and a nonlinear part, and without the nonlinear part we have a standard linear $p$th order autoregression with a single unit root (cf. model ((ref))).
To express ((ref)) in a different way, set $\pi_{j}=-\sum_{i=j+1}^{p}\varphi_{i}$ ($j=1,\ldots,p-1$; when $p=1$, set $\pi_{1}=\cdots=\pi_{p-1}=0$) so that we can express the polynomial $\varphi(z)$ as \[ \varphi(z)=(1-z)(1-\pi_{1}z-\cdots-\pi_{p-1}z^{p-1}), \] where the roots of the polynomial $\varpi(z)=1-\pi_{1}z-\cdots-\pi_{p-1}z^{p-1}$ lie outside the unit circle. This shows that we can write equation ((ref)) alternatively as
Denoting $u_{t}=y_{t}-\pi_{1}y_{t-1}-\cdots-\pi_{p-1}y_{t-p+1}$ equation ((ref)) can be written as
or as $u_{t}=u_{t-1}+\tilde{g}(y_{t-1},\ldots,y_{t-p})+\varepsilon_{t}$; when $p=1$ we obtain $y_{t}=y_{t-1}+\tilde{g}(y_{t-1})+\varepsilon_{t}$ as in ((ref)). The formulation in ((ref)) is convenient in our theoretical developments and will therefore be used instead of ((ref)). One reason for this convenience is that in cases where the function $\tilde{g}$ depends on $y_{t-1},\ldots,y_{t-p}$ only through the linear combination $u_{t-1}=y_{t-1}-\pi_{1}y_{t-2}-\cdots-\pi_{p-1}y_{t-p}$ we can write equation ((ref)) (with a slight abuse of notation) in a more compact way as $u_{t}=u_{t-1}+\tilde{g}(u_{t-1})+\varepsilon_{t}$. Then the process $u_{t}$ can be treated as the first-order model ((ref)) and, as will be discussed shortly, with a suitable assumption, we can make use of results in fort2003polynomial and douc2004practical (this turns out to be the case even when $\tilde{g}$ is not a function of the process $u_{t-1}$ only).
Next we introduce the assumptions needed to prove our results. Our first assumption restricts the dynamics in equation ((ref)).
Assumption (ref)(i) corresponds to the conventional stationarity condition of a linear autoregression in that it requires the roots of the polynomial $\varpi(z)$ to lie outside the unit circle. In the first-order case $p=1$, this condition becomes redundant because then $\pi_{1}=\cdots=\pi_{p-1}=0$.
Assumption (ref)(ii) requires the function $\tilde{g}$ to be bounded on compact subsets and links it to another function $g$. Condition (ii.a) controls the difference between the functions $u+\tilde{g}(\boldsymbol{x})$ and $g(u)$ or, in model ((ref)), the difference between the processes $u_{t-1}+\tilde{g}(y_{t-1},\ldots,y_{t-p})$ and $g(u_{t-1})$. In the special case where the function $\tilde{g}$ depends on $u$ only, condition (ii.a) becomes obvious because then one can choose $u+\tilde{g}(\boldsymbol{x})=g(u)$ and $\epsilon(\boldsymbol{x})=0$, and it suffices to check condition (ii.b) only. In this case we can use results in fort2003polynomial and douc2004practical directly in our proofs. However, we can do the same, albeit in a more complicated way, also when the function $\tilde{g}$ depends on the whole $p$-dimensional vector $\boldsymbol{x}$, but then the difference between the functions $u+\tilde{g}(\boldsymbol{x})$ and $g(u)$ may not increase “too fast” when $\left|\boldsymbol{x}\right|$ gets large. What is “too fast” is controlled by the function $\epsilon(\boldsymbol{x})$, and when $d\geq1$ the difference between $u+\tilde{g}(\boldsymbol{x})$ and $g(u)$ becomes negligible when $\left|\boldsymbol{x}\right|$ increases.
Condition (ii.a) implies that $\left|u+\tilde{g}(\boldsymbol{x})\right|\leq\left|g(u)\right|+\left|\epsilon(\boldsymbol{x})\boldsymbol{x}\right|$, which combined with condition (ii.b) yields
This fact is used in our proofs. Note also that condition (ii.a) is implied by the equality $u+\tilde{g}(\boldsymbol{x})=g(u)+\tilde{\epsilon}(\boldsymbol{x})\theta'\boldsymbol{x}$ where $\theta$ is a $p$-dimensional parameter vector and, if $\tilde{\epsilon}(\boldsymbol{x})=o(\left|\boldsymbol{x}\right|^{-d})$ is assumed, condition (ii.a) holds with $\epsilon(\boldsymbol{x})=\left|\theta\right|\tilde{\epsilon}(\boldsymbol{x})$. This approach for checking condition (ii.a) is illustrated in Section 5.
Condition (ii.b) is similar to its first-order counterpart ((ref)) to which it reduces when $p=1$. Note that apart from the boundedness condition, no restrictions are placed on $g(u)$ for moderate values of $u$. In the higher-order case this assumption concerns the filtered process $u_{t}=\varpi(L)y_{t}$. In the first-order case we also have $\boldsymbol{x}=u$ and the easiest way to verify Assumption 1(ii) may then be to define the function $g$ as $g(x)=x+\tilde{g}(x)$ (and $\epsilon(\boldsymbol{x})=0$), and verify condition (ii.b) directly. When $p\geq2,$ the fact that the domain of the function $\tilde{g}$ is larger than that of $g$ complicates the situation in that then no simple connection between inequalities ((ref)) and ((ref)) can generally be found. An example of this case is provided in Section 5.
Our second assumption gives conditions required of the error term in equation ((ref)).
Assumption 2(a) corresponds to Assumption 3.3 of douc2004practical, whereas Assumption 2(b) is a combination of the conditions imposed in (NSS 1), (NSS 4), and Lemma 3 of fort2003polynomial. The boundedness condition imposed in Assumption 2 on the density of the error term is stronger than would be needed but is used for simplicity (see Assumption 3.3 of douc2004practical for a more general alternative).
Note that finiteness of the first expectation in Assumption (ref)(a) holds with $\kappa_{0}=1$ if the distribution of $\varepsilon_{1}$ has a moment generating function in some interval of the origin. Although many widely used distributions satisfy this condition some heavy tailed distributions are ruled out (this applies to distributions whose densities cannot be bounded by a term of the form $c_{1}e^{-c_{2}\left|x\right|}$ with $c_{1}$ and $c_{2}$ positive constants). An example is Student's t-distribution irrespective of the value of the degrees of freedom parameter. The condition in Assumption (ref)(b) is used to address this issue. In this condition the case $0<\rho<s_{0}<1$ is rather extreme in that not even the expectation $E[\varepsilon_{1}]$ is assumed to exist.
We now present our ergodicity results which we base on model ((ref)). In Section 4.1 the rate of ergodicity established is subexponential whereas a slower polynomial rate of ergodicity is obtained in Section 4.2. The difference between these two cases stems from the assumed moment conditions: in Section 4.1 the condition in Assumption 2(a) is assumed whereas in Section 4.2 the weaker condition in Assumption 2(b) is employed. First we have to present the companion form of model ((ref)) which applies to both of these cases and will be needed in the proofs of our theorems.
To simplify notation, denote $\boldsymbol{y}_{t}=(y_{t},\ldots,y_{t-p+1})$ and define the function $\overline{g}\,:\,\mathbb{R}^{p}\rightarrow\mathbb{R}$ as \[ \overline{g}(\boldsymbol{x})=x_{1}-\pi_{1}x_{2}-\cdots-\pi_{p-1}x_{p}+\tilde{g}(\boldsymbol{x})=u+\tilde{g}(\boldsymbol{x}) \] so that $\overline{g}(\boldsymbol{y}_{t-1})=u_{t-1}+\tilde{g}(\boldsymbol{y}_{t-1})$. It is readily seen that the companion form related to equation ((ref)) reads as \[ \left[
\right]=
\left[
\right]+\overline{g}(\boldsymbol{y}_{t-1})\left[
\right]+\varepsilon_{t}\left[
\right] \] or, with obvious matrix notation,
(when $p=1$, $\boldsymbol{\Phi}=0$). Thus, Assumption (ref) ensures that $\boldsymbol{y}_{t}$ is a Markov chain on $\mathbb{R}^{p}$. For later purposes it is convenient to transform the companion form ((ref)). To this end, we define the matrices
where $\mathbf{A}$ is nonsingular and $\boldsymbol{\Pi}_{1}$ is the $(p-1)\times(p-1)$ dimensional lower right hand corner of $\mathbf{\Pi}$ (when $p=1$, $\mathbf{A}=1$ and $\mathbf{\Pi}=0$). With these definitions ((ref)) can be transformed into
where $\mathbf{A}\boldsymbol{y}_{t}=(u_{t},y_{t-1},\ldots,y_{t-p+1})$. Now, for any $p$-dimensional vector $\boldsymbol{x}$, form the partition $\boldsymbol{x}=(x_{1},\ldots,x_{p})=(x_{1},\boldsymbol{x}_{2})$ and define $\boldsymbol{z}(\boldsymbol{x})=(z_{1}(\boldsymbol{x}),\boldsymbol{z}_{2}(\boldsymbol{x}))=\boldsymbol{A}\boldsymbol{x}$, where (due to ((ref))) $z_{1}(\boldsymbol{x})=x_{1}-\pi_{1}x_{2}-\cdots-\pi_{p-1}x_{p}$ and $\boldsymbol{z}_{2}(\boldsymbol{x})=\boldsymbol{x}_{2}$ (when $p=1$, $\boldsymbol{x}_{2}$ and $\boldsymbol{z}_{2}(\boldsymbol{x})$ are dropped). With this notation equation ((ref)) can be expressed as
The first equation in ((ref)) is now in a form that can be analyzed by using the results in fort2003polynomial and douc2004practical. As for the second equation, by Assumption (ref)(i) the roots of the polynomial $\varpi(z)$ lie outside the unit circle, so that the eigenvalues of the matrix $\boldsymbol{\Pi}_{1}$ are smaller than one in absolute value. As is well known, this implies the existence of a matrix norm $\left\Vert \cdot\right\Vert _{*}$ induced by a vector norm, also denoted by $\left\Vert \cdot\right\Vert _{*}$, such that $\left\Vert \boldsymbol{\Pi}_{1}\right\Vert _{*}\leq\eta$ for some $\eta<1$ (see, e.g., Definition 5.6.1 and Lemma 5.6.10 in horn2013matrix). These facts will be useful in our proofs.
Our results make use of Condition D which requires choosing the function $V$. To this end, let $b_{1}$, $b_{2}$, and $b_{3}$ be positive constants whose values (to be specified later) depend on the constants $\beta_{0}$, $\kappa_{0}$, and $\rho$ introduced in Assumptions 1 and 2; for $b_{3}$ we already mention that it will always satisfy $b_{3}\in(0,1]$. When $p\geq2$, we define the function $V$ as
and when $p=1$, we define $V(x)=\exp\{b_{1}\left|x\right|^{b_{3}}\}$.
Now we can state the following theorem which makes use of the stronger moment requirement in Assumption 2(a). (The proof is given in the Supplementary Appendix.)
In Theorem 2, the case $\rho=\kappa_{0}$ represents a qualitative change in the ergodic behavior of the considered Markov chain: For $\rho>\kappa_{0}$ a slower subexponential convergence rate is obtained and $\rho=\kappa_{0}$ is the borderline case where a change to the faster geometric rate occurs. For $\rho<\kappa_{0}$, geometric ergodicity could also be established but we omit this case for brevity (in the first-order case this result is an immediate consequence of Theorem 3.3 of douc2004practical). Note also that, by the definition of the constant $b_{3}$, the rate of ergodicity in the subexponential case decreases as the value of $\rho$ increases.
Previously, douc2004practical obtained the results of Theorem 2 in the first-order case; our primary purpose here is to provide higher-order analogs of their results. klokov2004sub and klokov2007lower have also studied the first-order model ((ref)) satisfying inequality ((ref)) with $1<\rho<2$ but otherwise their assumptions are rather different from ours. They obtain subexponential bounds for ergodicity in total variation norm (i.e., ($1,r$)-ergodicity) and for $\beta$-mixing coefficients but they do not discuss general ($f,r$)-ergodicity.
As discussed in Section 2, we can also establish $\beta$-mixing and, in contrast to klokov2004sub and klokov2007lower, we can permit all initial values with distribution $\mu$ such that $\int_{x\in\mathbb{R}^{p}}V(x)\mu(dx)<\infty$ (and $V$ as in Theorem 2). Specifically, the discussion at the end of Section 2 and Theorem 2 imply $\beta$-mixing with the following rates: In case (i) the rate is subexponential, i.e., $\lim_{n\rightarrow\infty}e^{\tilde{k}n^{b_{3}/\rho}}\beta(n)=0$ with any $\tilde{k}\in(0,\left\{ c\rho/2b_{3}\right\} ^{b_{3}/\rho})$, and in case (ii) the rate is geometric, i.e., $\lim_{n\rightarrow\infty}\tilde{r}^{n}\beta(n)=0$ for some $\tilde{r}>1$ (or, equivalently, $\lim_{n\rightarrow\infty}e^{\tilde{c}n}\beta(n)=0$ for some $\tilde{c}>0$).
The following corollary is an immediate consequence of the discussion after equality ((ref)) (for a formal result, see Theorem 14.0.1 in meyn2009markov).
As we remarked after Theorem 1, in the subexponential case of Theorem 2 results with a faster rate of convergence and/or larger $f$-norm could be obtained at the expense of more complex notation. This means that in the above corollary finiteness of slightly larger moments could be obtained and the subexponential $\beta$-mixing rate discussed above could similarly be slightly improved.
Next we consider ergodicity results relying only on the weaker moment requirement in Assumption 2(b). This will below lead to a slower polynomial rate of ergodicity. The key result used to relax the moment requirement is Lemma 3 of fort2003polynomial (which the authors use in conjunction with their analog of Condition D; we depart from their approach and use Condition D which corresponds to an analogous condition described in Section 1.2 in fort2003polynomial).
The function $V$ employed is now different from the subexponential case. When $p\geq2$, we define the function $V$ as
where $s_{0}$ is as in Assumption 2(b), $\alpha=1-\rho/s_{0}$ with $\rho$ as in Assumption 1(ii.b), and $s_{1}$ is a positive constant (to be specified later); when $p=1$, we define $V(x)=1+\left|x\right|^{s_{0}}$.
The following theorem presents the ergodicity result obtained when using the weaker moment condition in Assumption 2(b). (The proof is given in the Supplementary Appendix.)
Options (i)\textendash (iii) in Theorem 3 represent the combinations of the values of $\rho$ and $s_{0}$ for which the result can be obtained by relying on the corresponding cases (i)\textendash (iii) in Lemma 3 of fort2003polynomial. Unlike in Theorem 2, the case $\rho=2$ is allowed, but then an additional and rather intricate moment condition is required. A further departure from Theorem 2 is that the same polynomial rate of ergodicity is obtained in all cases. However, similarly to the subexponential case in Theorem 2, the rate of ergodicity decreases as the value of $\rho$ increases. Also, from the discussion at the end of Section 2 we can conclude that the rate of $\beta$-mixing implied by Theorem 3 is polynomial and, specifically, $\lim_{n\rightarrow\infty}n^{s_{0}/\rho-1}\beta(n)=0$.
The first-order case of Theorem 3 was obtained by fort2003polynomial (with slightly different assumptions). Polynomial ergodicity results for first-order autoregressions similar to that in ((ref)) have previously appeared also in tuominen1994subgeometric (in the case $0<\rho<1$), tanikawa2001markov (in the case $\rho=1$), and \textcolor{black}{veretennikov2000polynomial and klokov2007lower (in the case $\rho=2$; these authors also obtain polynomial bounds for }$\beta$-mixing coefficients\textcolor{black}{ }but do not consider general ($f,r$)-ergodicity\textcolor{black}{).}
The following corollary on the moments of the stationary distribution is proved in the Supplementary Appendix. (In contrast to the subexponential case, using the ($f,r$)-ergodicity result of Theorem 3 would here yield a weaker moment result; hence, some extra steps are needed.)
In this section we discuss special cases of the general model introduced in Section (ref). Using the three equivalent formulations ((ref))\textendash ((ref)) in Section 3, the model considered can be written as \[ y_{t}-\varphi_{1}y_{t-1}-\cdots-\varphi_{p}y_{t-p}=\Delta y_{t}-\pi_{1}\Delta y_{t-1}-\cdots-\pi_{p-1}\Delta y_{t-p+1}=u_{t}-u_{t-1}=\tilde{g}(y_{t-1},\ldots,y_{t-p})+\varepsilon_{t}, \] where the polynomial $\varphi(z)$ has precisely one unit root, can be decomposed as $\varphi(z)=(1-z)\varpi(z)$, and $u_{t}=\varpi(L)y_{t}$. Further formal assumptions will be stated in Propositions 1 and 2 below.
First-order subgeometrically ergodic autoregressions were already exemplified, albeit at a rather general level, in fort2003polynomial and douc2004practical. In meitz2019subgemix we study rates of subgeometric ergodicity and $\beta$-mixing in a first-order multi-regime self-exciting threshold autoregressive (SETAR) model; the proof of Theorem 3 in that paper illustrates how Theorems 2 and 3 of the present paper can be applied in a first-order case. In what follows we focus on examples of higher-order subgeometrically ergodic autoregressive models.
We consider three main examples. We first give a heuristic overview of them and then present the formalities. The first example we consider can be expressed as
where $I(u_{t-1})$ is not constant and will be interpreted as a time-varying drift or intercept term (the case where $I(u_{t-1})$ is a constant is not of interest here because then model ((ref)) reduces to a nonstationary unit root model (with or without a drift)). We consider the case where $I(u_{t-1})$ takes values in a bounded interval and fluctuates suitably between increasing and decreasing drifts. This ensures that model ((ref)) will be (subgeometrically) ergodic and stationary under appropriate conditions (see Proposition 1 below) even though its dynamics involve a unit root component and a (time-varying) drift.
The second example we consider is
where $\nu\in\mathbb{R}$ is an intercept term and $S(u_{t-1})$ will be interpreted as a time-varying slope coefficient. In the cases $S(u_{t-1})\equiv0$ and $S(u_{t-1})\equiv1$ model ((ref)) reduces to the (linear) stationary model $\varpi(L)y_{t}=u_{t}=\varepsilon_{t}$ and to the nonstationary unit root model $\varphi(L)y_{t}=\varepsilon_{t}$, respectively. We consider the case of $S(u_{t-1})$ being time-varying, taking values in some interval $[s,1)$ ($s<1$), and attaining values arbitrarily close to 1 for values of $u_{t-1}$ large in absolute value. Under appropriate conditions (see Proposition 2 below), model ((ref)) will be (subgeometrically) ergodic and stationary although exhibiting features similar to a unit root process.
Our third example generalizes the second one and illustrates that one can allow nonlinear dependence not only on $u_{t-1}$ but on the entire $\boldsymbol{y}_{t-1}$. Specifically, we consider model
where we have omitted the intercept for simplicity, $S(u_{t-1})$ again represents a time-varying slope coefficient, and the term $F(\boldsymbol{y}_{t-1})$ captures the nonlinear dependence on the entire $\boldsymbol{y}_{t-1}$.
We consider example ((ref)) with the time-varying intercept term $I(u_{t-1})$ specified as in logistic smooth transition autoregressive (LSTAR) models (see, e.g., vandijk2002smooth). Specifically, we choose
with $L(u;b,a)=1/(1+e^{-b(u-a)})$ denoting the logistic function. The parameters $b,a_{1},a_{2}$ are assumed to satisfy $b>0$ and $a_{1}\leq a_{2}$ as usual, and $\nu_{1},\nu_{2}$ are assumed to satisfy $\nu_{1}<0<\nu_{2}$ to obtain ergodicity below. The time-varying intercept term $I(u_{t-1})$ now takes values in the interval $(\nu_{1},\nu_{2})$. Note that for large values of $b$ the logistic function $L(u;b,a)$ is close to the indicator function and then this model provides a close approximation to the above-mentioned threshold autoregressive model where $L(u_{t-1};b,a_{i})$ is replaced with an indicator function.
The following proposition shows the ergodicity of this model.
The proof of Proposition 1 is a straightforward application of Theorems 2 and 3 in the case $\rho=1$ (see the Appendix). Depending on moment assumptions the rate of ergodicity is geometric, subexponential, or polynomial, and similar rates also apply to $\beta$-mixing coefficients (see the discussions following Theorems 2 and 3 as well as the Corollaries to these theorems for existence of finite moments of the stationary distribution).
We now provide intuitive and graphical illustrations of the behavior processes covered by Proposition 1 can exhibit. An informal description captures the main features concisely: For `values of $u_{t-1}$ in the extreme left tail', the intercept term $I(u_{t-1})$ is close to $\nu_{2}>0$ resulting in increasing drift towards `central values of $u_{t-1}$'; conversely, for `values of $u_{t-1}$ in the extreme right tail', $I(u_{t-1})$ is close to $\nu_{1}<0$ and decreasing drift towards `central values of $u_{t-1}$' takes place; finally, for such `central values of $u_{t-1}$', unit root behavior without drift occurs. This informal description is illustrated in the top row of Figure (ref) using three different cases of the function $I(\cdot)$ in ((ref)); the precise parameter values used can be found in the caption of Figure (ref).
The graph on the left illustrates the simplest possible case. When $u$ takes negative values the function $I(u)$ takes positive values; for $u<-2$ the values of $I(u)$ are very close to its supremum, in this example $0.08$. As $u_{t}=I(u_{t-1})+u_{t-1}+\varepsilon_{t}$, the behavior of the process $u_{t}$ is for $u_{t-1}<0$ (and particularly for $u_{t-1}<-2$) close to that of a unit root process with an increasing drift. A similar behavior occurs when $u_{t-1}$ takes positive values (and particularly when $u_{t-1}>2$) but then the drift is decreasing. Thus, large absolute values of $u_{t-1}$ induce a drift towards the origin which provides intuition why the process $u_{t}$ can be ergodic and stationary even though its behavior shows resemblance to a unit root process with a (negligible) drift near the origin where $I(u)\approx0$. The other two graphs in the top row of Figure (ref) illustrate similar but somewhat more involved possibilities for the function $I(u)$. The middle graph illustrates that the `unit root regime' can take place over a wider range of values of $u$ and these values need not be centered at the origin; in this graph $I(u)\approx0$ for $u$ roughly between $1.5$ and $3.5$. The graph on the right illustrates a possibility in which also the `middle regime' corresponds to a unit root with positive drift.
The second row of Figure (ref) presents three examples of simulated time series of $y_{t}$. These are generated from second-order versions of model ((ref)) using the three time-varying intercept terms $I(u_{t-1})$ depicted in the first row. In all cases the autoregressive coefficient $\pi_{1}$ is equal to 0.75 and the error terms $\varepsilon_{t}$ are generated from rescaled Student's $t$-distribution with five degrees of freedom (and different variances, see the caption for details). The third and fourth rows of Figure (ref) show the time series graphs of the time-varying intercept terms $I(u_{t})$ and the autocorrelation functions of $y_{t}$, respectively, in these three examples. The time series graphs of $y_{t}$ in the second row bear a resemblance to those of unit root processes but nevertheless exhibit mean-reverting behavior. However, the mean reversion takes place slower than it would for a geometrically ergodic process. The autocorrelation functions also show very strong persistence. These features are related to the fact that the considered models are only polynomially ergodic with rate $r(n)=n^{4-\delta}$ for some $\delta>0$ (this follows from Proposition 1 as the error terms used only have moments of order smaller than five).
We notice from Figure (ref) that the periods when $u_{t}$ takes large absolute values can be rather long and they can contain both increasing and decreasing periods for $y_{t}$. An example is the time series in the first column around $t\approx800$ corresponding to the largest peak of the time series of $y_{t}$; similar features occur also in the time series in the second and third columns, often around peaks or troughs of the series. We also note from the second and third columns of Figure (ref) that the time series of $I(u_{t})$ stays close to 0 or 0.03, respectively, for some time (corresponding to behavior of $u_{t}$ close to that of a unit root process without a drift or with a drift). An example is the rather long decreasing period in the time series of $y_{t}$ in the third column, roughly between $t\approx200$ and $t\approx500$, that is to large extent due to unit root type behavior.
Next we consider example ((ref)) with the time-varying slope term $S(u_{t-1})$ being either
where $r_{0}>0$ and the positive-valued function $h$ is such that $h(u)$ is large whenever $u$ is large in absolute value (formal requirements for $h$ are given in Proposition 2 below). Then the time-varying slope term $S(u_{t-1})$ takes values in some interval $[s,1)$ ($s<1$) and attains values arbitrarily close to 1 for values of $u_{t-1}$ large in absolute value. The shapes of $S_{1}(u)$ and $S_{2}(u)$ as functions of $u$ resemble the `inverted bell curve form' commonly employed in exponential STAR (ESTAR) models (see, e.g., vandijk2002smooth). The main difference between the two functions in ((ref)) is that $S_{2}(u)$ takes values in the unit interval $(0,1)$ whereas $S_{1}(u)$ can also take negative values.
Before providing concrete examples we state the following proposition which imposes conditions on the function $h$ above to ensure that the results of Theorems 2 and 3 hold for model ((ref)) with $S(\cdot)$ as in ((ref)). The proof is straightforward and available in the Appendix.
The obtained rate of ergodicity is again either geometric, subexponential, or polynomial depending on the moment assumptions made (the precise convergence rates can be obtained from Theorems 2 and 3; for existence of moments, see the corollaries to Theorems 2 and 3). The last part of the proposition lists several potential concrete choices for the function $h$, of which cases (iii) and (vi) may be convenient if differentiability of the function $h$ is desired.
To illustrate the type of behavior processes covered by Proposition 2 may exhibit, consider as a simple example model ((ref)) with $\nu=0$, the slope term $S_{1}(u_{t-1})$ in ((ref)), and $h(u)=1+\lvert u\rvert^{\rho}$ (cf. ((ref))). The considered model then becomes
The shape of the slope coefficient $S_{1}(u_{t-1})$ as a function of $u_{t-1}$ is now similar to that of an inverted bell curve increasing monotonically to unity as $\lvert u_{t-1}\rvert$ increases, and $S_{1}(u_{t-1})$ takes values within the interval $[1-r_{0},1)$. Note that depending on the value of $r_{0}$, the slope coefficient takes values potentially only near the unity (say, in the interval $[0.95,1)$) or even in rather extreme ranges (say, in the interval $[-100,1)$) so that very different behaviors can be accommodated.
To explain why such processes can still be ergodic and stationary, note that for large values of $\lvert u_{t-1}\rvert$ the slope $S_{1}(u_{t-1})$ always takes values within $(-1,1)$ (or a much smaller subset of it near unity). This prevents the process $u_{t}$ from exploding and ensures mean-reverting behavior. Note that for larger values of $\rho$ (within the permitted range $0<\rho\leq2$), the slope $S_{1}(u_{t-1})$ approaches unity faster as $\lvert u_{t-1}\rvert$ increases, intuitively corresponding to more `wandering' behavior of the observed series. This is reflected in Proposition 2 as slower rates of ergodicity being related to larger values of the parameter $\rho$ (see also Theorems 2 and 3).
Figure (ref) illustrates the preceding discussion. The top row depicts examples of the function $S(u)$ with some particular choices of the function $h$ (see the caption of Figure 2 for the details). In the figures on the left and in the middle $S(u)=S_{1}(u)$ whereas in the figure on the right $S(u)=S_{2}(u)$. In each example the shape of the function $S(u)$ is similar to that of an inverted bell curve increasing monotonically to unity as $\lvert u\rvert$ increases.
Rows 2\textendash 4 of Figure (ref) present simulated time series of $y_{t}$, time series graphs of the time-varying slope coefficients $S(u_{t-1})$, and the autocorrelation functions of $y_{t}$, respectively, in these three examples. In the first and second columns the time series graphs of $y_{t}$ and the related autocorrelation functions indicate unit root type behavior and strong persistence; in these cases also the slope coefficient $S(u_{t-1})$ takes values mostly larger than $0.90$ and often very close to one. Compared to these to cases, the time series graph of $y_{t}$ in the third column appears less \textquoteleft wandering\textquoteright and this feature is also reflected in the related autocorrelation function which decays faster and in a slope coefficient $S(u_{t-1})$ taking values further away from unity. Despite the three examples exhibiting somewhat different behaviours, all of them exhibit mean-reverting behavior and are subexponentially ergodic (this is due to Proposition 2 because $\rho>1$ and the error terms used are normally distributed and hence satisfy Assumption 2(a) with $\kappa_{0}=1$). The mean reversion again takes place slower than would be the case for geometrically ergodic processes.
Finally, we briefly consider example ((ref)) and for simplicity set $p=2$. The time-varying slope term $S(u_{t-1})$ can be either one of the two options in ((ref)). As for $F(\bm{y}_{t-1})$, we set $F(\bm{y}_{t-1})=\exp\{-\gamma\lvert\boldsymbol{y}_{t-1}\rvert^{2}\}(\theta_{1}y_{t-1}+\theta_{2}y_{t-2})$, where $\boldsymbol{y}_{t-1}=(y_{t-1},y_{t-2})$, $\gamma>0$, and $\theta=(\theta_{1},\theta_{2})$ can take any values in $\mathbb{R}^{2}$. That is, the considered model reads as
On the right hand side of ((ref)) the term $\exp\{-\gamma\lvert\boldsymbol{y}_{t-1}\rvert^{2}\}$ has a bell-shaped form (as a function of $\lvert\boldsymbol{y}_{t-1}\rvert$) with maximum at the origin while for choices such as $h(u_{t-1})=1+\lvert u\rvert^{\rho}$ the shape of $S(u_{t-1})$ is that of an inverted bell curve. Thus, given the shape of the terms $S(u_{t-1})$ and $\exp\{-\gamma\lvert\boldsymbol{y}_{t-1}\rvert^{2}\}$, model ((ref)) can be viewed as a certain type of three-regime ESTAR model (see, e.g., vandijk2002smooth).
It is straightforward to check that model ((ref)) satisfies Assumption 1(ii) for both of the two options of $S$ and with $h$ as in Proposition 2 (details are available in the Appendix; note that in this example the function $\tilde{g}$ in equation ((ref)) does not depend on $u_{t-1}$ only). It may be worth noting that Assumption 1(ii.a) does not hold if we replace the norm $\lvert\boldsymbol{y}_{t-1}\rvert$ in ((ref)) with a linear function of $\boldsymbol{y}_{t-1}$ such as $u_{t-1}$. Another point worth noting is that Assumption 1(ii) does not rule out the possibility of setting $\pi_{1}=0$ (and $u_{t-1}=y_{t-1}$) in ((ref)), and similarly for its higher-order counterparts where some or even all the coefficients $\pi_{1},\ldots,\pi_{p-1}$ may be equal to zero.
Allowing the parameters $\theta_{1}$ and $\theta_{2}$ in model ((ref)) to be totally unrestricted highlights the fact that the autoregressions we consider may exhibit rather arbitrary behavior for moderate values of the observed series. As indicated in the Introduction, geometrically ergodic nonlinear autoregressions with features of this kind have previously been considered by lu1998geometric, gourieroux2006stochastic, bec2008acr, and others. However, in most of these previous models stationarity is approached the further away the process moves from the origin whereas in our model unit root type behavior prevails for large absolute values of the process.
In this paper we examined the subgeometric ergodicity of certain higher-order nonlinear autoregressive models. Generalizing existing first-order results, we provided conditions that ensure subexponential and polynomial ergodicity of the considered autoregressions. These results were established by utilizing suitably formulated drift conditions. Relying on results in a companion paper meitz2019subgemix, useful conclusions on the convergence rates of $\beta$-mixing coefficients were also obtained.
After obtaining theoretical results for rather general models we considered concrete examples and illustrated them with simulation. However, further work is needed to judge the usefulness of these models in practical applications. Several extensions could also be envisioned. For instance, subgeometric ergodicity of multivariate higher-order autoregressions or of models with conditional heteroskedasticity are interesting topics left for future work.