The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
77,137 characters
Subgeometrically ergodic autoregressions
\title{\vspace{20pt}
Subgeometrically ergodic autoregressions\thanks{The authors thank the Academy of Finland for financial support, and
Co-Editor Donald Andrews and two anonymous referees for useful comments
and suggestions. Contact addresses: Mika Meitz, Department of Economics,
University of Helsinki, P. O. Box 17, FI\textendash 00014 University
of Helsinki, Finland; e-mail: [email removed]. Pentti Saikkonen,
Department of Mathematics and Statistics, University of Helsinki,
P. O. Box 68, FI\textendash 00014 University of Helsinki, Finland;
e-mail: [email removed].}\vspace{20pt}
}
\author{Mika Meitz\\\small{University of Helsinki} \and Pentti Saikkonen\\\small{University of Helsinki}\vspace{20pt}
}
\date{First version April 2019, revised February 2020}
\maketitle
\begin{abstract}
\noindent In this paper we discuss how the notion of subgeometric
ergodicity in Markov chain theory can be exploited to study stationarity
and ergodicity of nonlinear time series models. Subgeometric ergodicity
means that the transition probability measures converge to the stationary
measure at a rate slower than geometric. Specifically, we consider
suitably defined higher-order nonlinear autoregressions that behave
similarly to a unit root process for large values of the observed
series but we place almost no restrictions on their dynamics for moderate
values of the observed series. Results on the subgeometric ergodicity
of nonlinear autoregressions have previously appeared only in the
first-order case. We provide an extension to the higher-order case
and show that the autoregressions we consider are, under appropriate
conditions, subgeometrically ergodic. As useful implications we also
obtain stationarity and $\beta$-mixing with subgeometrically decaying
mixing coefficients.
\bigskip{}
\bigskip{}
\bigskip{}
\noindent\textbf{JEL classification:} C22
\bigskip{}
\noindent \textbf{Keywords:} Nonlinear autoregressive model, subgeometric
ergodicity, Markov chain, $\beta$-mixing.
\end{abstract}
\vfill{}
\pagebreak{}
\section{Introduction}
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{sec:markov}. 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,
\begin{equation}
\lim_{n\to\infty}r^{n}\lVert P^{n}(x\,;\,\cdot)-\pi\rVert_{TV}=0\label{eq:Geom-erg}
\end{equation}
(the definition of $\left\Vert \cdot\right\Vert _{TV}$ and a formulation
of (\ref{eq:Geom-erg}) using a more general norm are given in Section
\ref{sec:markov}). 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 \citet[Sec 1.1]{doukhan1994mixing} or \citet[Ch 3]{bradley2007introduction}.)\textcolor{red}{{}
}The authoritative and classic reference to Markov chain theory is
the monograph of \citet{meyn1993markov,meyn2009markov}. Recent papers
establishing geometric ergodicity of different nonlinear time series
models include \citet{francq2006mixing}, \citet{ling2007double},
\citet{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{eq:Geom-erg}) takes place at a rate $r(n)$ slower than
geometric, that is,
\begin{equation}
\lim_{n\to\infty}r(n)\lVert P^{n}(x\,;\,\cdot)-\pi\rVert_{TV}=0.\label{eq:SubGeom-erg}
\end{equation}
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 \citet{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
\begin{equation}
y_{t}=g(y_{t-1})+\varepsilon_{t},\quad t=1,2,\ldots,\label{NLAR(1)}
\end{equation}
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
\begin{equation}
\left|g(x)\right|\leq\left(1-r\left|x\right|^{-\rho}\right)\left|x\right|\quad\textrm{for }\left|x\right|\geq M_{0}\qquad\qquad[r>0,\,M_{0}>r^{1/\rho},\,0<\rho\leq2],\label{Ineq g(x)_p=00003D1}
\end{equation}
and that $g(x)$ is bounded for $\left|x\right|\leq M_{0}$. A concrete
example where (\ref{Ineq g(x)_p=00003D1}) can be easily verified
is
\begin{equation}
y_{t}=\left(1-\frac{r_{0}}{1+\left|y_{t-1}\right|^{\rho}}\right)y_{t-1}+\varepsilon_{t}\qquad\qquad[r_{0}>0,\,0<\rho\leq2].\label{Example_1}
\end{equation}
The model defined in (\ref{Example_1}) 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{Example_1}) 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{Example_1})
can be viewed as a special case of the model
\begin{equation}
y_{t}=y_{t-1}+\tilde{g}(y_{t-1})+\varepsilon_{t},\label{RW-type example_2}
\end{equation}
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 \citet{jarner2003necessary}; see particularly
equation (3) of their paper.} In Section \ref{sec:model}, a higher-order version of equation (\ref{RW-type example_2})
(without assuming boundedness) is used as a starting point of the
formulation of our general model. Our main results in Section \ref{sec:results}
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{Ineq g(x)_p=00003D1})
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, \citet{lu1998geometric}, \citet{gourieroux2006stochastic},
and \citet{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 \citet{tuominen1994subgeometric},
\citet{veretennikov2000polynomial}, \citet{fort2003polynomial},
\citet{douc2004practical}, \nocite{klokov2004sub,klokov2005subexponential}Klokov
and Veretennikov (2004, 2005), and \citet{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 \citet{fort2003polynomial} and \citet{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{sec:markov}
contains basic concepts of Markov chains and summarizes existing results
on subgeometric ergodicity. Section \ref{sec:model} 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{sec:results}.
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$).
\section{Markov chains and subgeometric ergodicity\label{sec:markov} }
In this section we discuss basic concepts of Markov chains needed
to obtain our results. More comprehensive discussions can be found
in \citet{meyn2009markov} and \citet{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 \citet{fort2003polynomial}
and \citet{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
\begin{equation}
\left\Vert \mu\right\Vert _{f}=\sup_{f_{0}:\left|f_{0}\right|\leq f}\left|\mu(f_{0})\right|,\label{eq:f-norm}
\end{equation}
where $\mu(f_{0})=\int_{x\in\mathsf{X}}f_{0}(x)\mu(dx)$ (and the
supremum in (\ref{eq:f-norm}) 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{eq:Geom-erg}) and (\ref{eq:SubGeom-erg}).
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
\begin{equation}
\lim_{n\to\infty}r(n)\lVert P^{n}(x\,;\,\cdot)-\pi\rVert_{f}=0\qquad\text{for }\pi\text{-almost all }x\in\mathsf{X}.~\footnotemark\label{f-ergodicity}
\end{equation}
\textcolor{red}{\footnotetext{That is, the convergence in (\ref{f-ergodicity}) 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{f-ergodicity}) 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{f-ergodicity})~with\hspace*{-10pt}
\[
\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 \citet[p.�776]{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{f-ergodicity}) is given by $r(n)=r^{n}$, $r>1$. The subgeometric
rate functions we consider are defined as follows (cf., e.g., \citet{nummelin1983rate}
and \citet{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$. \citet[Sec 3.3]{douc2004practical} consider subexponential
convergence rates whereas \citet[Sec 2.2]{fort2003polynomial} consider
polynomial convergence rates in model (\ref{NLAR(1)}) (see also the
related references mentioned in these papers).
The proofs of our results make use of the following condition adapted
from \citet[Defn 16.1.7]{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 \citet{douc2004practical}.}
\bigskip{}
\vspace*{\fill}
\pagebreak{}
\noindent \textbf{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
\begin{equation}
E\left[V(X_{1})\,\left|\,X_{0}=x\right.\right]\leq V(x)-\phi\left(V(x)\right)+b\boldsymbol{1}_{C}(x),\qquad x\in\mathsf{X}.\label{Drift condition}
\end{equation}
\bigskip{}
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 \citet{meyn2009markov}.)
\begin{thm}[(\citet{meyn2009markov}, \citet{douc2004practical})]
Suppose $X_{t}$ is a $\psi$-irreducible and aperiodic Markov chain
on $(\mathsf{X},\mathcal{B}(\mathsf{X}))$ and that Condition D holds
with a petite set $C$ such that $\sup_{x\in C}V(x)<\infty$ and the
function $\phi$ being either\vspace*{-2pt}
\begin{lyxlist}{(ii;�subexponential�case) }
\item [{(i;~geometric~case)}] $\phi(v)=\lambda v$ for some $\lambda>0$,\vspace*{-6pt}
\item [{(ii;~subexponential~case)}] $\phi(v)=c(v+v_{0})/[\ln(v+v_{0})]^{\alpha}$
for some $c,\alpha,v_{0}>0$, or\vspace*{-6pt}
\item [{(iii;~polynomial~case)}] $\phi(v)=cv^{\alpha}$ for some $\alpha\in[0,1)$
and $c\in(0,1]$.\vspace*{-2pt}
\end{lyxlist}
Then $X_{t}$ is $(f,r)$-ergodic with either\vspace*{-2pt}
\begin{lyxlist}{(iii)}
\item [{(i)}] $f=V$ and $r(n)=r^{n}$ for some $r>1$ (or, equivalently,
$r(n)=(e^{c})^{n}$ for some $c>0$),\vspace*{-6pt}
\item [{(ii)}] $f=V^{\delta}$ and $r(n)=(e^{d})^{n^{1/(1+\alpha)}}$ for
any $\delta\in(0,\negthinspace1)$ and any $d\in(0,(1\negthinspace-\negthinspace\delta)\negthinspace\left\{ c(1\negthinspace+\negthinspace\alpha)\right\} ^{1/(1+\alpha)})$,
or\hspace*{-4pt}\vspace*{-6pt}
\item [{(iii)}] $f=V^{1-\delta(1-\alpha)}$ and $r(n)=n^{\delta-1}$ for
any $\delta\in[1,1/(1-\alpha)]$.\bigskip{}
\end{lyxlist}
\end{thm}
In the geometric case the result of Theorem 1 is given in \citet{meyn2009markov},
and in the subexponential and polynomial cases the result can be obtained
from \citet{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 \citet[Sec 2.3]{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{f-ergodicity}) 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{RW-type example_2}))
are given in \citet{jarner2003necessary} (for an application of this
result to a threshold autoregressive model, see \citet{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 \citet[Sec 1.1]{doukhan1994mixing}, \citet[Ch 3]{bradley2007introduction},
or \citet{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 \citet{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 \citet{meitz2019subgemix},
and note that their Theorem 2(e) is also used to obtain the subexponential
rate shown here.}
\begin{lyxlist}{XXXXXXXXXXXXXXX}
\item [{(i)~geometric~case:}] $\lim_{n\rightarrow\infty}\tilde{r}^{n}\beta(n)=0$
for some $\tilde{r}>1$;
\item [{(ii)~subexponential~case:}] $\lim_{n\rightarrow\infty}(e^{\tilde{d}})^{n^{1/(1+\alpha)}}\beta(n)=0$
for any $\tilde{d}\in(0,\left\{ c(1\negthinspace+\negthinspace\alpha)/2\right\} ^{1/(1+\alpha)})$;
\item [{(iii)~polynomial~case:}] $\lim_{n\rightarrow\infty}n^{\alpha/(1-\alpha)}\beta(n)=0$.
\end{lyxlist}
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).
\section{Model and assumptions\label{sec:model}}
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
\begin{equation}
y_{t}=\varphi_{1}y_{t-1}+\cdots+\varphi_{p}y_{t-p}+\tilde{g}(y_{t-1},\ldots,y_{t-p})+\varepsilon_{t},\label{NLAR(p)_phi}
\end{equation}
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{RW-type example_2})).
To express (\ref{NLAR(p)_phi}) 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{NLAR(p)_phi}) alternatively as
\begin{equation}
y_{t}=y_{t-1}+\pi_{1}\Delta y_{t-1}+\cdots+\pi_{p-1}\Delta y_{t-p+1}+\tilde{g}(y_{t-1},\ldots,y_{t-p})+\varepsilon_{t}.\label{DNLAR(p)}
\end{equation}
Denoting $u_{t}=y_{t}-\pi_{1}y_{t-1}-\cdots-\pi_{p-1}y_{t-p+1}$ equation
(\ref{DNLAR(p)}) can be written as
\begin{equation}
y_{t}=\pi_{1}y_{t-1}+\cdots+\pi_{p-1}y_{t-p+1}+u_{t-1}+\tilde{g}(y_{t-1},\ldots,y_{t-p})+\varepsilon_{t}\label{NLAR(p)_pi}
\end{equation}
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{RW-type example_2}). The formulation in (\ref{NLAR(p)_pi})
is convenient in our theoretical developments and will therefore be
used instead of (\ref{DNLAR(p)}). 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{NLAR(p)_pi}) (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{RW-type example_2}) and, as will be discussed shortly, with
a suitable assumption, we can make use of results in \citet[Sec 2.2]{fort2003polynomial}
and \citet[Sec 3.3]{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{NLAR(p)_pi}).
\begin{assumption}
\noindent \label{assu:dynamics}Suppose the polynomial $\varpi(z)=1-\pi_{1}z-\cdots-\pi_{p-1}z^{p-1}$
and the function $\tilde{g}\,:\,\mathbb{R}^{p}\rightarrow\mathbb{R}$
in (\ref{NLAR(p)_pi}) satisfy the following conditions:\vspace*{-2pt}
\begin{lyxlist}{(ii)}
\item [{(i)}] \noindent The roots of $\varpi(z)$ lie outside the unit
circle.\vspace*{-2pt}
\item [{(ii)}] The function $\tilde{g}$ is measurable, bounded on compact
subsets of $\mathbb{R}^{p}$, and there exists a measurable function
$g\,:\,\mathbb{R}\rightarrow\mathbb{R}$ with the property $\left|g(x)\right|\rightarrow\infty$
as $\left|x\right|\rightarrow\infty$ such that the following two
conditions hold.\vspace*{-2pt}
\begin{lyxlist}{(ii.b)}
\item [{(ii.a)}] \noindent With $\boldsymbol{x}=(x_{1},\ldots,x_{p})$
and $u=x_{1}-\pi_{1}x_{2}-\cdots-\pi_{p-1}x_{p}$, the function $\tilde{g}$
satisfies
\begin{equation}
\left|u+\tilde{g}(\boldsymbol{x})-g(u)\right|\leq\left|\epsilon(\boldsymbol{x})\boldsymbol{x}\right|,\label{Inequality Ass 2a}
\end{equation}
where $\epsilon(\boldsymbol{x})$ is a real-valued function such that
$\left|\epsilon(\boldsymbol{x})\right|=o(\left|\boldsymbol{x}\right|^{-d})$
as $\left|\boldsymbol{x}\right|\rightarrow\infty$ for some $d>0$.\vspace*{-2pt}
\item [{(ii.b)}] There exist positive constants $r$, $M_{0}$, $K_{0}$,
and $0<\rho\leq2$ such that for all $u\in\mathbb{R}$
\begin{equation}
\left|g(u)\right|\leq\begin{cases}
(1-r\left|u\right|^{-\rho})\left|u\right| & \textrm{for }\left|u\right|\geq M_{0},\\
K_{0} & \textrm{for }\left|u\right|\leq M_{0}.
\end{cases}\label{Inequality_Ass 2}
\end{equation}
\end{lyxlist}
\end{lyxlist}
\end{assumption}
Assumption \ref{assu:dynamics}(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{assu:dynamics}(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{NLAR(p)_pi}),
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 \citet[Sec 2.2]{fort2003polynomial} and \citet[Sec 3.3]{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
\begin{equation}
\left|u+\tilde{g}(\boldsymbol{x})\right|\leq\left(1-r\left|u\right|^{-\rho}\right)\left|u\right|+o(\left|\boldsymbol{x}\right|^{-d})\left|\boldsymbol{x}\right|\quad\textrm{for }\left|u\right|\geq M_{0}.\label{Implication of Ass 2}
\end{equation}
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{Ineq g(x)_p=00003D1})
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{Inequality Ass 2a})
and (\ref{Inequality_Ass 2}) 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{NLAR(p)_pi}).
\begin{assumption}
\label{assu:errors}$\{\varepsilon_{t},\,t=1,2,\ldots\}$ is a sequence
of IID random variables that is independent of $(y_{0},\ldots,y_{-p+1})$
(with $p$ as in Assumption 1), the distribution of $\varepsilon_{1}$
has a (Lebesgue) density that is bounded away from zero on compact
subsets of $\mathbb{R}$, and either
\begin{lyxlist}{00}
\item [{(a)}] $E\bigl[e^{\beta_{0}\left|\varepsilon_{1}\right|^{\kappa_{0}}}\bigr]<\infty$
for some $\ensuremath{\beta_{0}>0}$ and $\ensuremath{\kappa_{0}\in(0,1]}$,
and $E[\varepsilon_{1}]=0$; or
\item [{(b)}] $E\left[\left|\varepsilon_{1}\right|^{s_{0}}\right]<\infty$
for some $s_{0}>\rho$ (with $\rho$ as in Assumption 1), and $E[\varepsilon_{1}]=0$
holds if $\rho\geq1$.
\end{lyxlist}
\end{assumption}
Assumption 2(a) corresponds to Assumption 3.3 of \citet[Sec 3.3]{douc2004practical},
whereas Assumption 2(b) is a combination of the conditions imposed
in (NSS 1), (NSS 4), and Lemma 3 of \citet[Sec 2.2]{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 \citet[Sec 3.3]{douc2004practical} for a more
general alternative).
Note that finiteness of the first expectation in Assumption \ref{assu:errors}(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{assu:errors}(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.
\section{Results\label{sec:results}}
We now present our ergodicity results which we base on model (\ref{NLAR(p)_pi}).
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{NLAR(p)_pi}) 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{NLAR(p)_pi})
reads as
\[
\left[\begin{array}{c}
y_{t}\\
y_{t-1}\\
\vdots\\
\vdots\\
y_{t-p+1}
\end{array}\right]=\begin{bmatrix}\pi_{1} & \pi_{2} & \cdots & \pi_{p-1} & 0\\
1 & 0 & \cdots & 0 & 0\\
0 & \ddots & \ddots & \vdots & \vdots\\
\vdots & \ddots & \ddots & 0 & 0\\
0 & \cdots & 0 & 1 & 0
\end{bmatrix}\left[\begin{array}{c}
y_{t-1}\\
y_{t-2}\\
\vdots\\
\vdots\\
y_{t-p}
\end{array}\right]+\overline{g}(\boldsymbol{y}_{t-1})\left[\begin{array}{c}
1\\
0\\
\vdots\\
\vdots\\
0
\end{array}\right]+\varepsilon_{t}\left[\begin{array}{c}
1\\
0\\
\vdots\\
\vdots\\
0
\end{array}\right]
\]
or, with obvious matrix notation,
\begin{equation}
\boldsymbol{y}_{t}=\boldsymbol{\Phi}\boldsymbol{y}_{t-1}+\overline{g}(\boldsymbol{y}_{t-1})\boldsymbol{\iota}_{p}+\varepsilon_{t}\boldsymbol{\iota}_{p}\label{Companion form}
\end{equation}
(when $p=1$, $\boldsymbol{\Phi}=0$). Thus, Assumption \ref{assu:errors}
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{Companion form}). To this end, we define the matrices
\begin{equation}
\negthickspace\negthickspace\negthickspace\mathbf{A}=\setlength{\arraycolsep}{2pt}
\global\long\begin{bmatrix}1 & -\pi_{1} & -\pi_{2} & \cdots & -\pi_{p-1}\\
0 & 1 & 0 & \cdots & 0\\
\vdots & \:\:\:\:\ddots & \ddots & \ddots & \vdots\\
0 & \:\:\: & \ddots & \ddots & 0\\
0 & \cdots & \cdots & 0 & 1
\end{bmatrix}\quad\text{and}\quad\mathbf{\Pi}=\mathbf{A}\boldsymbol{\Phi}\mathbf{A}^{-1}=\begin{bmatrix}0 & 0 & 0 & \cdots & 0\\
1 & \pi_{1} & \pi_{2} & \cdots & \pi_{p-1}\\
0 & 1 & 0 & \cdots & 0\\
\vdots & \ddots & \ddots & \ddots & \vdots\\
0 & \cdots & 0 & 1 & 0
\end{bmatrix}=\begin{bmatrix}0 & \text{\textbf{0}}'_{p-1}\\
\boldsymbol{\iota}_{p-1} & \boldsymbol{\Pi}_{1}
\end{bmatrix},\negthickspace\negthickspace\negthickspace\negthickspace\label{Matrix Pi}
\end{equation}
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{Companion form}) can be transformed into
\begin{equation}
\mathbf{A}\boldsymbol{y}_{t}=\mathbf{\Pi}\mathbf{A}\boldsymbol{y}_{t-1}+\overline{g}(\boldsymbol{y}_{t-1})\boldsymbol{\iota}_{p}+\varepsilon_{t}\boldsymbol{\iota}_{p},\label{Companion form_A}
\end{equation}
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{Matrix Pi})) $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{Companion form_A})
can be expressed as
\begin{equation}
\begin{bmatrix}z_{1}(\boldsymbol{y}_{t})\\
\boldsymbol{z}_{2}(\boldsymbol{y}_{t})
\end{bmatrix}=\setlength{\arraycolsep}{2pt}\begin{bmatrix}0 & \text{\textbf{0}}'_{p-1}\\
\boldsymbol{\iota}_{p-1} & \boldsymbol{\Pi}_{1}
\end{bmatrix}\begin{bmatrix}z_{1}(\boldsymbol{y}_{t-1})\\
\boldsymbol{z}_{2}(\boldsymbol{y}_{t-1})
\end{bmatrix}+\overline{g}(\boldsymbol{y}_{t-1})\boldsymbol{\iota}_{p}+\varepsilon_{t}\boldsymbol{\iota}_{p}=\begin{bmatrix}\overline{g}(\boldsymbol{y}_{t-1})+\varepsilon_{t}\\
\boldsymbol{\Pi}_{1}\boldsymbol{z}_{2}(\boldsymbol{y}_{t-1})+z_{1}(\boldsymbol{y}_{t-1})\boldsymbol{\iota}_{p-1}
\end{bmatrix}.\label{Companion form_A2}
\end{equation}
The first equation in (\ref{Companion form_A2}) is now in a form
that can be analyzed by using the results in \citet{fort2003polynomial}
and \citet{douc2004practical}. As for the second equation, by Assumption
\ref{assu:dynamics}(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
\citet{horn2013matrix}). These facts will be useful in our proofs.
\subsection{Subexponential case}
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
\begin{equation}
V(\boldsymbol{x})=\tfrac{1}{2}\exp\bigl\{ b_{1}\left|z_{1}(\boldsymbol{x})\right|^{b_{3}}\bigr\}+\tfrac{1}{2}\exp\bigl\{ b_{2}\left\Vert \boldsymbol{z}_{2}(\boldsymbol{x})\right\Vert _{*}^{b_{3}}\bigr\}\label{Def. V_1}
\end{equation}
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.)
\begin{thm}
\noindent Suppose $p\geq2$ and consider the Markov chain $\boldsymbol{y}_{t}$
defined in equation (\ref{Companion form}). Let Assumptions 1 and
2(a) hold, suppose that in Assumption 1 the constants $\rho$ and
$d$ satisfy $0<\rho<2$ and $d=\rho/b_{3}$, and let $V(\boldsymbol{x})$
be as in (\ref{Def. V_1}).
\begin{lyxlist}{000}
\item [{(i)}] \noindent If $\rho>\kappa_{0}$, then $\boldsymbol{y}_{t}$
is ($f,r$)-ergodic with the subexponential convergence rate $r(n)=e^{kn^{b_{3}/\rho}}$
and the function $f$ given by $f(\boldsymbol{x})=V(\boldsymbol{x})^{\delta}$;
this result holds for any choice of $\delta\in(0,1)$, for any $k$
such that $0<k<(1-\delta)\left\{ c\rho/b_{3}\right\} ^{b_{3}/\rho}$,
and for some (small) $b_{1},\,b_{2}\in(0,\beta_{0})$, $b_{3}=\kappa_{0}\wedge(2-\rho)\in(0,1)$,
and some (small) $c>0$.
\item [{(ii)}] If $\rho=\kappa_{0}$, then $\boldsymbol{y}_{t}$ is geometrically
ergodic with the convergence rate $r(n)=e^{cn}$ and the function
$f$ given by $f(\boldsymbol{x})=V(\boldsymbol{x})$; this result
holds for some (small) $b_{1},\,b_{2}\in(0,\beta_{0})$, $b_{3}=\kappa_{0}\in(0,1]$,
and some $c>0$.
\end{lyxlist}
When $p=1$, consider the Markov chain $y_{t}$ defined by $y_{t}=y_{t-1}+\tilde{g}(y_{t-1})+\varepsilon_{t}$.
The above results hold for $y_{t}$ with the function $V$ defined
as $V(x)=\exp\{b_{1}\left|x\right|^{b_{3}}\}$ (and the constant $d$
becomes redundant).
\end{thm}
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 \citet[Sec 3.3]{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, \citet[Sec 3.3]{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. \citet[2005]{klokov2004sub}
and \citet{klokov2007lower} have also studied the first-order model
(\ref{NLAR(1)}) satisfying inequality (\ref{Ineq g(x)_p=00003D1})
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 \citet[2005]{klokov2004sub} and \citet{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{f-ergodicity}) (for a formal result, see Theorem
14.0.1 in \citet{meyn2009markov}).
\begin{cor*}[\textbf{to Theorem 2}]
Let $\pi$ signify the stationary distribution of $\boldsymbol{y}_{t}$
in Theorem 2 and let the function $f$ be as in cases (i) and (ii)
of Theorem 2. Then $\pi(f)=\int_{\boldsymbol{x}\in\mathbb{R}^{p}}f(\boldsymbol{x})\pi(d\boldsymbol{x})<\infty$;
in particular, $\pi(\left|\boldsymbol{x}\right|^{s})<\infty$ for
all $s>0$ so that the stationary distribution has finite moments
of all orders.
\end{cor*}
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.
\subsection{Polynomial case}
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 \citet[Sec 2.2]{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 \citet{fort2003polynomial}).
The function $V$ employed is now different from the subexponential
case. When $p\geq2$, we define the function $V$ as
\begin{equation}
V(\boldsymbol{x})=1+\left|z_{1}(\boldsymbol{x})\right|^{s_{0}}+s_{1}\left\Vert \boldsymbol{z}_{2}(\boldsymbol{x})\right\Vert _{*}^{\alpha s_{0}},\label{Def. V_2}
\end{equation}
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.)
\begin{thm}
\noindent Suppose $p\geq2$ and consider the Markov chain $\boldsymbol{y}_{t}$
defined in equation (\ref{Companion form}). Let Assumptions 1 and
2(b) hold, suppose that in Assumption 1 the constants $\rho$ and
$d$ satisfy $0<\rho\leq2$ and $d=\rho/s_{0}$ when $s_{0}<1$ and
$d=\rho$ when $s_{0}\geq1$, and let $V(\boldsymbol{x})$ be as in
(\ref{Def. V_2}). Assume further that either \vspace*{-4pt}
\begin{lyxlist}{000000}
\item [{~~~(i)}] \noindent $0<\rho<1$ and $s_{0}>\rho$,\vspace*{-6pt}
\item [{~~~(ii)}] \noindent $1\leq\rho<2$ and either $s_{0}=2$ or
$s_{0}\geq4$, or\vspace*{-6pt}
\item [{~~~(iii)}] \noindent $\rho=2$ and $s_{0}\geq4$ with $s_{0}r-\frac{1}{2}s_{0}(s_{0}-1)E[\varepsilon_{1}^{2}]>0$.\vspace*{-4pt}
\end{lyxlist}
\noindent Then $\boldsymbol{y}_{t}$ is ($f,r$)-ergodic with the
polynomial convergence rate $r(n)=n^{\delta-1}$ and the function
$f$ given by $f(\boldsymbol{x})=V(\boldsymbol{x})^{1-\delta\rho/s_{0}}$;
this result holds for any choice of $\delta\in[1,s_{0}/\rho]$ and
for some (small) $s_{1}>0$.
\medskip{}
\noindent When $p=1$, consider the Markov chain $y_{t}$ defined
by $y_{t}=y_{t-1}+\tilde{g}(y_{t-1})+\varepsilon_{t}$. The above
results hold for $y_{t}$ with the functions $V$ and $f$ defined
as $V(x)=1+\left|x\right|^{s_{0}}$ and $f(x)=1+\left|x\right|^{s_{0}-\delta\rho}$
(and the constant $d$ becomes redundant).
\end{thm}
\noindent 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 \citet[Sec 2.2]{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 \citet[Sec 2.2]{fort2003polynomial}
(with slightly different assumptions). Polynomial ergodicity results
for first-order autoregressions similar to that in (\ref{NLAR(1)})
have previously appeared also in \citet[Sec 5.2]{tuominen1994subgeometric}
(in the case $0<\rho<1$), \citet{tanikawa2001markov} (in the case
$\rho=1$), and \textcolor{black}{\citet{veretennikov2000polynomial}
and \citet{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.)
\begin{cor*}[\textbf{to Theorem 3}]
Let $\pi$ signify the stationary distribution of $\boldsymbol{y}_{t}$
in Theorem 3. Then $\pi(f)=\int_{\boldsymbol{x}\in\mathbb{R}^{p}}f(\boldsymbol{x})\pi(d\boldsymbol{x})<\infty$
with $f(\boldsymbol{x})=\left|\boldsymbol{x}\right|^{s_{0}-\rho}$
so that the stationary distribution has finite moments up to order
$s_{0}-\rho$.
\end{cor*}
\section{Illustrative examples}
In this section we discuss special cases of the general model introduced
in Section \ref{sec:model}. Using the three equivalent formulations
(\ref{NLAR(p)_phi})\textendash (\ref{NLAR(p)_pi}) 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 \citet[Sec 2.2]{fort2003polynomial}
and \citet[Secs 3.3, 3.4]{douc2004practical}. In \citet[Sec 5]{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
\begin{equation}
u_{t}=u_{t-1}+I(u_{t-1})+\varepsilon_{t},\label{eq:Ex1}
\end{equation}
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{eq:Ex1}) 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{eq:Ex1}) 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
\begin{equation}
u_{t}-\nu=S(u_{t-1})(u_{t-1}-\nu)+\varepsilon_{t},\label{eq:Ex2}
\end{equation}
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{eq:Ex2}) 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{eq:Ex2}) 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
\begin{equation}
u_{t}=S(u_{t-1})u_{t-1}+F(\boldsymbol{y}_{t-1})+\varepsilon_{t},\label{eq:Ex3}
\end{equation}
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}$.
\subsection{Example with time-varying intercept term of LSTAR type}
We consider example (\ref{eq:Ex1}) with the time-varying intercept
term $I(u_{t-1})$ specified as in logistic smooth transition autoregressive
(LSTAR) models (see, e.g., \citet{vandijk2002smooth}). Specifically,
we choose
\begin{equation}
I(u_{t-1})=\nu_{1}L(u_{t-1};b,a_{1})+\nu_{2}(1-L(u_{t-1};b,a_{2}))\label{Function I(u)}
\end{equation}
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.
\begin{prop}
Consider the process $y_{t}$ defined by $u_{t}=u_{t-1}+I(u_{t-1})+\varepsilon_{t}$
as in (\ref{eq:Ex1}) (with $u_{t}=\varpi(L)y_{t}$ and the roots
of $\varpi(z)$ outside the unit circle), and with $I(u_{t-1})$ as
in (\ref{Function I(u)}) (with $\nu_{1}<0<\nu_{2}$). Assume further
that either\vspace*{-4pt}
\begin{lyxlist}{00}
\item [{\small\ \ (1)}] Assumption 2(a) is satisfied with $\kappa_{0}\in(0,1)$,\vspace*{-6pt}
\item [{\small\ \ (2)}] Assumption 2(a) is satisfied with $\kappa_{0}=1$,
or\vspace*{-6pt}
\item [{\small\ \ (3)}] Assumption 2(b) is satisfied with either $s_{0}=2$
or $s_{0}\geq4$.\vspace*{-4pt}
\end{lyxlist}
Then, under condition {\small (1)/(2)/(3)}, the process $\boldsymbol{y}_{t}=(y_{t},\ldots,y_{t-p+1})$
is either\vspace*{-4pt}
\begin{lyxlist}{00}
\item [{\small\ \ (1)}] subexponentially ergodic with convergence rate
$r(n)=(e^{k})^{n^{\kappa_{0}}}$ (for some $k>0$),\vspace*{-6pt}
\item [{\small\ \ (2)}] geometrically ergodic with convergence rate $r(n)=(e^{c})^{n}$
(for some $c>0$), or\vspace*{-6pt}
\item [{\small\ \ (3)}] polynomially ergodic with convergence rate $r(n)=n^{s_{0}-1}$.
\end{lyxlist}
\end{prop}
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{FigEx1} using three different cases of the function
$I(\cdot)$ in (\ref{Function I(u)}); the precise parameter values
used can be found in the caption of Figure~\ref{FigEx1}.
\begin{figure}[p]
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig2a}
\end{minipage}\hfill{}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig2b}
\end{minipage}\hfill{}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig2c}
\end{minipage}\vspace*{-25pt}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig31a}
\end{minipage}\hfill{}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig32a}
\end{minipage}\hfill{}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig33a}
\end{minipage}\vspace*{-25pt}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig31b}
\end{minipage}\hfill{}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig32b}
\end{minipage}\hfill{}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig33b}
\end{minipage}\vspace*{-25pt}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig31c}
\end{minipage}\hfill{}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig32c}
\end{minipage}\hfill{}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig33c}
\end{minipage}\vspace*{-15pt}
\caption{\label{FigEx1} \uline{Top row.} Graphs of the function $I(u)=\nu_{1}L(u;b,a_{1})+\nu_{2}(1-L(u;b,a_{2}))$
in (\ref{Function I(u)}). Left: $b=2,$ $\nu_{1}=-0.08$, $\nu_{2}=0.08$,
and $a_{1}=a_{2}=0$. Middle: $b=4$, $\nu_{1}=-0.08$, $\nu_{2}=0.08$,
$a_{1}=1$, $a_{2}=4$. Right: $b=5$, $\nu_{1}=-0.05$, $\nu_{2}=0.08$,
$a_{1}=1$, and $a_{2}=4$ (the horizontal dotted line at $0.03$).
\uline{Second row.} Simulated time series of $y_{t}$ corresponding
to the time-varying intercept functions $I(u_{t})$ in the top row.
The 1000 observations are generated from model (\ref{eq:Ex1}) with
$p=2$, $\pi_{1}=0.75$, and $\varepsilon_{t}$ distributed as rescaled
Student's $t$ distribution with five degrees of freedom (and $E[\varepsilon_{1}^{2}]$
equal to $0.300$ on the left, $0.250$ in the middle, and $0.125$
on the right). \uline{Third row.} The corresponding time series
graphs of $I(u_{t})$. The dotted lines show the minimum and maximum,
and the value zero in the first and second columns and $0.03$ in
the third column. \uline{Bottom row.} The corresponding autocorrelation
functions of the simulated series. The first three autocorrelations
are 0.999, 0.996, 0.992 (left); 0.997, 0.991, 0.984 (middle); and
0.996, 0.987, 0.973 (right). }
\end{figure}
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{FigEx1} 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{FigEx1} presents three examples of
simulated time series of $y_{t}$. These are generated from second-order
versions of model (\ref{eq:Ex1}) 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{FigEx1} 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{FigEx1} 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{FigEx1} 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.
\subsection{Example with time-varying slope term of ESTAR type}
Next we consider example (\ref{eq:Ex2}) with the time-varying slope
term $S(u_{t-1})$ being either
\begin{equation}
S_{1}(u_{t-1})=1-\frac{r_{0}}{h(u_{t-1})}\qquad\text{or}\qquad S_{2}(u_{t-1})=\exp\Bigl\{-\frac{r_{0}}{h(u_{t-1})}\Bigr\},\label{Functions S}
\end{equation}
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.,
\citet{vandijk2002smooth}). The main difference between the two functions
in (\ref{Functions S}) 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{eq:Ex2}) with
$S(\cdot)$ as in (\ref{Functions S}). The proof is straightforward
and available in the Appendix.
\begin{prop}
Consider the process $y_{t}$ defined by $u_{t}-\nu=S(u_{t-1})(u_{t-1}-\nu)+\varepsilon_{t}$
as in (\ref{eq:Ex2}) (with $u_{t}=\varpi(L)y_{t}$ and the roots
of $\varpi(z)$ outside the unit circle) with $S(u_{t-1})$ being
either $S_{1}(u_{t-1})$ or $S_{2}(u_{t-1})$ in (\ref{Functions S})
(with $r_{0}>0$) and with the function $h$ therein satisfying\vspace*{-4pt}
\begin{lyxlist}{000}
\item [{(h)}] $h:\mathbb{R}\rightarrow(0,\infty)$ is measurable, bounded
on compact sets, satisfies $h(u)\rightarrow\infty$ as $\left|u\right|\rightarrow\infty$,
and is such that $c_{1}h(u)\leq\left|u\right|^{\rho}$ and $\left|u\right|^{\rho+c_{2}}\leq c_{3}h^{2}(u)$
for $\left|u\right|\geq M_{0}$, some $c_{1},c_{2},c_{3},M_{0}>0$,
and $0<\rho\leq2$.\vspace*{-4pt}
\end{lyxlist}
Assume further that either Assumption\vspace*{-4pt}
\begin{lyxlist}{00}
\item [{\small\ \ (1)}] 2(a) is satisfied with $\kappa_{0}<\rho$,\vspace*{-6pt}
\item [{\small\ \ (2)}] 2(a) is satisfied with $\kappa_{0}=\rho$, or\vspace*{-6pt}
\item [{\small\ \ (3)}] 2(b) is satisfied with an $s_{0}$ such that one
of the conditions (i)\textendash (iii) of Theorem 3 holds.\vspace*{-4pt}
\end{lyxlist}
Then, under condition {\small (1)/(2)/(3)}, the process $\boldsymbol{y}_{t}=(y_{t},\ldots,y_{t-p+1})$
is either\vspace*{-4pt}
\begin{lyxlist}{00}
\item [{\small\ \ (1)}] subexponentially ergodic,\vspace*{-6pt}
\item [{\small\ \ (2)}] geometrically ergodic, or\vspace*{-6pt}
\item [{\small\ \ (3)}] polynomially ergodic.\vspace*{-4pt}
\end{lyxlist}
Moreover, condition (h) above is satisfied for i) $h(u)=1+\lvert u-a\rvert^{\rho}$,
ii) $h(u)=(1+\lvert u-a\rvert)^{\rho}$, iii) $h(u)=(1+(u-a)^{2})^{\rho/2}$,
iv) $h(u)=1+\lvert u-a_{1}\rvert^{\rho_{1}}+\lvert u-a_{2}\rvert^{\rho_{2}}$,
v) $h(u)=1+(1+\lvert u-a_{1}\rvert)^{\rho_{1}}+(1+\lvert u-a_{2}\rvert)^{\rho_{2}}$,
or vi) $h(u)=1+(1+(u-a_{1})^{2})^{\rho_{1}/2}+(1+(u-a_{2})^{2})^{\rho_{2}/2}$
($\rho,\rho_{1},\rho_{2}\in(0,2]$; $a,a_{1},a_{2}\in\mathbb{R}$).
\end{prop}
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{eq:Ex2})
with $\nu=0$, the slope term $S_{1}(u_{t-1})$ in (\ref{Functions S}),
and $h(u)=1+\lvert u\rvert^{\rho}$ (cf.~(\ref{Example_1})). The
considered model then becomes
\begin{equation}
u_{t}=\biggl(1-\frac{r_{0}}{1+\lvert u_{t-1}\rvert^{\rho}}\biggr)u_{t-1}+\varepsilon_{t}.\label{Specific example_1}
\end{equation}
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{FigEx2} 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.
\begin{figure}[p]
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig4a}
\end{minipage}\hfill{}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig4b}
\end{minipage}\hfill{}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[clip,width=1.1\textwidth]{Fig4c}
\end{minipage}\vspace*{-25pt}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig51a}
\end{minipage}\hfill{}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig52a}
\end{minipage}\hfill{}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[bb=0bp 0bp 396bp 287bp,width=1.1\textwidth]{Fig53a}
\end{minipage}\vspace*{-25pt}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig51b}
\end{minipage}\hfill{}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig52b}
\end{minipage}\hfill{}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig53b}
\end{minipage}\vspace*{-25pt}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig51c}
\end{minipage}\hfill{}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig52c}
\end{minipage}\hfill{}
\begin{minipage}[t][1\totalheight][c]{0.33\textwidth}
\includegraphics[width=1.1\textwidth]{Fig53c}
\end{minipage}\vspace*{-15pt}
\caption{\label{FigEx2} \uline{Top row.} Graphs of $S(u)$ in (\ref{Functions S}).
Left: $S_{1}(u)$ with $r_{0}=0.15$ and $h(u)=1+\left|u\right|^{1.5}$.
Middle: $S_{1}(u)$ with $r_{0}=0.5$ and $h(u)=1+\lvert u+4\rvert^{1.25}+\lvert u+8\rvert^{1.25}$.
Right: $S_{2}(u)$\textbf{ }with $r_{0}=1.5$ and $h(u)=1+\left|u+2\right|^{1.5}$.
\uline{Second row.} Simulated time series of $y_{t}$ corresponding
to the time-varying slope functions $S(u_{t})$ in the top row. The
1000 observations are generated from model (\ref{eq:Ex2}) with $p=2$,
$\pi_{1}=0.75$, $\nu=2$, and Gaussian $\varepsilon_{t}$ ($\varepsilon_{t}\sim N(0,0.25)$
on the left; $\varepsilon_{t}\sim N(0,0.75)$ in the middle; $\varepsilon_{t}\sim N(0,1.5)$
on the right). \uline{Third row.} The corresponding time series
graphs of $S_{1}(u_{t-1})=1-r_{0}/h(u_{t-1})$ on the left and in
the middle, and $S_{2}(u_{t-1})=\exp\{-r_{0}/h(u_{t-1})\}$ on the
right. The dotted lines show the minimum and maximum. \uline{Bottom
row.} The corresponding autocorrelation functions of the simulated
series. The first three autocorrelations are 0.997, 0.990, 0.981 (left);
0.996, 0.989, 0.980 (middle); and 0.992, 0.974, 0.950 (right). }
\end{figure}
Rows 2\textendash 4 of Figure \ref{FigEx2} 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.
\subsection{Example of a more general formulation}
Finally, we briefly consider example (\ref{eq:Ex3}) 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{Functions S}). 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
\begin{equation}
y_{t}=\pi_{1}y_{t-1}+S(u_{t-1})(y_{t-1}-\pi_{1}y_{t-2})+\exp\{-\gamma\lvert\boldsymbol{y}_{t-1}\rvert^{2}\}(\theta_{1}y_{t-1}+\theta_{2}y_{t-2})+\varepsilon_{t}.\label{General example_2}
\end{equation}
On the right hand side of (\ref{General example_2}) 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{General example_2}) can be viewed as a certain type of
three-regime ESTAR model (see, e.g., \citet{vandijk2002smooth}).
It is straightforward to check that model (\ref{General example_2})
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{NLAR(p)_phi})
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{General example_2}) 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{General example_2}), 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{General example_2})
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 \citet{lu1998geometric}, \citet{gourieroux2006stochastic},
\citet{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.
\section{Conclusions}
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 \citet{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.