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.
65,915 characters · 11 sections · 77 citation commands
Subgeometrically ergodic autoregressions with autoregressive conditional heteroskedasticity
Let $X_{t}$ ($t=0,1,2,\ldots$) be a Markov chain on the state space $\mathsf{X}$ and initialized from an $X_{0}$ following some initial distribution. If the $n$-step probability measures $P^{n}(x\,;\,\cdot)=\Pr(X_{n}\in\cdot\mid X_{0}=x)$ converge in total variation norm $\lVert\,\cdot\,\rVert_{TV}$ to the stationary probability measure $\pi$ at rate $r^{n}$ (for some $r>1$), that is,
the Markov chain is said to be geometrically ergodic. When the convergence in ((ref)) takes place at a suitably defined rate $r(n)$ slower than geometric, that is,
the Markov chain is called subgeometrically ergodic. Examples of common rates (where $c$ denotes a positive constant) include geometric (or exponential) when $r(n)=e^{cn}=r^{n}$ ($r>1$), subexponential when $r(n)=e^{cn^{\gamma}}$ ($0<\gamma<1$), polynomial when $r(n)=(1+n)^{c}$, and logarithmic when $r(n)=(1+\ln(n))^{c}$. The authoritative and classic reference to Markov chain theory is the monograph of meyn2009markov, while an up-to-date treatment of subgeometric ergodicity can be found in Chapters 16 and 17 of \citet*{douc2018markov}.
To give some background, the notion of subgeometric ergodicity was introduced in the Markov chain literature in the 1980s when nummelin1983rate and tweedie1983criteria obtained the first subgeometric ergodicity results for general state space Markov chains. Subsequent work by tuominen1994subgeometric, fort2000vsubgeometric, jarner2002polynomial, fort2003polynomial, and \citet*{douc2004practical} lead to a formulation of a so-called drift condition to ensure subgeometric ergodicity, paralleling the use of a Foster-Lyapunov drift condition to establish geometric ergodicity (see, e.g., meyn2009markov). Various topics in probability theory and statistics have also been considered under subgeometric assumptions; for instance, \citet*{douc2008bounds} considered the central limit theorem and Berry-Esseen bounds, atchade2010limit the convergence of Markov chain Monte Carlo algorithms, \citet*{merlevede2011bernstein} a Bernstein-type inequality, and meitz2021JAP the rate of $\beta$-mixing. In this paper we are interested in autoregressive time series models. Results regarding the subgeometric ergodicity of first-order autoregressions were obtained by tuominen1994subgeometric, veretennikov2000polynomial, fort2003polynomial, douc2004practical, \nocite{klokov2004sub,klokov2005subexponential}Klokov and Veretennikov (2004, 2005), and klokov2007lower, among others, whereas results for more general higher-order autoregressions were obtained by meitz2022subgear.
In this paper we consider subgeometric (specifically, polynomial) ergodicity of autoregressive models with autoregressive conditional heteroskedasticity (ARCH; engle1982autoregressive). The previous works on subgeometrically ergodic autoregressions listed above only considered the case of independent and identically distributed (IID) errors, and allowing for conditionally heteroskedastic errors considerably widens the scope of potential applications. This is particularly important in applications using economic and financial time series data. In the subgeometrically ergodic AR\textendash ARCH models we consider, the conditional mean is similar to the (homoskedastic) AR models already considered in meitz2022subgear. The precise model formulation will be given and motivated further in Section 2, but we already note that the models we consider accommodate for behavior similar to a unit root process for large values of the observed series but almost no restrictions are placed on their dynamics for moderate values of the observed series. The conditional variance is allowed to follow a rather general nonlinear ARCH process. In our main result, we show that the considered AR\textendash ARCH processes are, under appropriate conditions, subgeometrically ergodic at a polynomial rate;\textcolor{red}{ }the convergence rate of $\beta$-mixing coefficients and finiteness of certain moments are also obtained (for details, see Section 3.2).
The inclusion of ARCH (instead of IID) errors considerably complicates the proofs of (sub)geometric ergodicity of nonlinear autoregressions. Papers considering subgeometric ergodicity of homoskedastic autoregressions were already listed above. Geometric ergodicity of nonlinear autoregressive models with ARCH (or generalized ARCH) errors has previously been considered by numerous authors; see, e.g., cline2004stability, \nocite{meitz2008stability,meitz2010SPL} Meitz and Saikkonen (2008, 2010), and the many references therein. Compared to these two strands of previous literature, the combination of the subgeometrically ergodic type of nonlinear dynamics in the conditional mean with ARCH errors leads to additional complications in the proofs. To appropriately separate these two sources of dynamics we make use of a (relatively unknown) extension of Bernoulli\textquoteright s inequality due to FeffermanShapiro1972 (combined with Young's inequality), and to control terms arising due to conditional heteroskedasticity we devise a special matrix norm that is of a more complicated type than the norms typically used when analysing the stability of nonlinear time series models.
The rest of the paper is organized as follows. Section 2 introduces the nonlinear AR\textendash ARCH model considered and states the assumptions we employ. Results on subgeometric ergodicity are given in Section 3. In Section 4 we consider an empirical application of our model to a daily time series of an energy sector volatility index. Section 5 concludes. All proofs are collected in an Appendix.
We consider the univariate process $y_{t}$ ($t=1,2,\ldots$) generated by
where $p\geq1$ is the autoregressive order, $u_{t}=y_{t}-\pi_{1}y_{t-1}-\cdots-\pi_{p-1}y_{t-p+1}$, $g$ is a real-valued function, $\varepsilon_{t}$ is an IID error term, and $\sigma_{t}=\sigma(\boldsymbol{y}_{t-1})$ is a positive volatility term that depends on $p+q$ lagged values of $y_{t}$, $\boldsymbol{y}_{t-1}=(y_{t-1},\ldots,y_{t-p-q})$, where $q\geq1$ is an ARCH order. For now, one concrete example of the volatility term is a linear ARCH process, where $\sigma_{t}$ satisfies
and $e_{t}=y_{t}-\pi_{1}y_{t-1}-\cdots-\pi_{p-1}y_{t-p+1}-g(u_{t-1})$, $\omega>0$, and $\alpha_{i}\geq0$ ($i=1,\ldots,q$); a more general formulation for the conditional variance will be considered below. Note that a compact expression for $e_{t}$ is $e_{t}=u_{t}-g(u_{t-1})$ so that equation ((ref)) can be expressed as $u_{t}=g(u_{t-1})+\sigma_{t}\varepsilon_{t}$. If $\pi_{1}=\cdots=\pi_{p-1}=0$ in equation ((ref)), we have $u_{t}=y_{t}$ so that the autoregressive order $p$ reduces to one and equation ((ref)) reduces to $y_{t}=g(y_{t-1})+\sigma_{t}\varepsilon_{t}$.
Our first assumption contains basic requirements for the error term $\varepsilon_{t}$ and makes clear that the squared volatility, $\sigma_{t}^{2}$, is the conditional variance of $y_{t}$ (when appropriate moments exist).
Later on we introduce an assumption on the conditional variance $\sigma_{t}^{2}$ which further restricts the moments of $\varepsilon_{t}$.
To further describe the conditional mean of the autoregressions we consider, we next specify the conditions needed for the function $g$ in equation ((ref)). The following assumption is a simplification of Assumption 1 in meitz2022subgear (the somewhat more general formulation used therein is briefly discussed at the end of this subsection).
Assumption 2(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 2(ii) is needed to prove the subgeometric ergodicity of the process $y_{t}$, as already done by fort2003polynomial and douc2004practical in the first-order case $p=1$ and by meitz2022subgear for higher-order autoregressions.
We next provide some intuition and motivation for our model in ((ref)). To clarify the role of inequality ((ref)) restricting the function $g(\cdot)$, suppose Assumptions 1 and 2(i) hold but instead of Assumption 2(ii) suppose the function $g(\cdot)$ were linear with $g(u)=\pi_{0}u$ and $\pi_{0}\in[-1,1]$. Using the lag operator $L$, equation ((ref)) could then be written as
that is, as the familiar linear AR($p$) model (with autoregressive heteroskedasticity). Given Assumptions 1 and 2(i), the case $\pi_{0}\in(-1,1)$ corresponds to geometric ergodicity of $y_{t}$ and the cases $\pi_{0}=\pm1$ to non-ergodicity. Nonlinear functions $g(\cdot)$ satisfying Assumption 2(ii) provide a middle ground between these extreme cases of geometric ergodicity and non-ergodicity. For instance, if $g(u)=(1-r\left|u\right|^{-\rho})u$ for $|u|>r^{1/\rho}$ and $g(u)=0$ otherwise ($r>0$, $0<\rho<2$), then for any fixed $\pi_{0}\in(-1,1)$ and for all $u$ sufficiently large in absolute value (i.e., for the values of $u$ that are crucial for determining ergodicity), \[ |\pi_{0}||u|<|g(u)|<|u|. \] The subgeometrically ergodic autoregressions we consider thus provide one possibility for modeling small departures from unit root autoregressions. Assumption 2(ii) implies that for large values of $|u_{t-1}|$, the conditional mean of model ((ref)) is close to that of an integrated process (of order one). On the other hand, as inequality ((ref)) restricts the function $g(\cdot)$ only for large values of its argument, no restrictions (apart from the boundedness condition in ((ref))) are imposed when the argument takes values inside some bounded set of values. Thus the autoregressions we consider may exhibit rather arbitrary (stationary, unit root, explosive, nonlinear, etc.) behavior for moderate values of the observed series.
The autoregressions we consider are to some extent related to existing models that have autoregressive roots near unity. To illustrate, when $g(u)$ is as in the previous paragraph and we further set $p=r=\rho=1$, the model in ((ref)) simplifies to \[ y_{t}=\Bigl(1-\frac{1}{|y_{t-1}|}\Bigr)y_{t-1}+e_{t}\quad\text{when }|y_{t-1}|>1\quad\text{and}\quad y_{t}=e_{t}\quad\text{otherwise} \] where $e_{t}=\sigma_{t}\varepsilon_{t}$. In comparison, a prototypical local-to-unity autoregression could be expressed as \[ y_{t}=\Bigl(1-\frac{1}{T}\Bigr)y_{t-1}+e_{t},\quad t=1,\ldots,T,\quad\text{where \ensuremath{T} denotes the sample size}. \] Both of the above formulations involve an autoregressive coefficient near unity, the former when the observed process takes on large (absolute) values and the latter when the sample size is large. However, the fact that the sample size is an essential part of local-to-unity autoregressions makes them quite different from the autoregressions we consider \textemdash in particular, the autoregressions we consider are ergodic. For more details on local-to-unity autoregressions and other related models, we refer the reader to the recent contributions of lieberman2020hybrid and phillips2023estimation and the references therein.
Homoskedastic subgeometrically ergodic autoregressions satisfying (a somewhat more general version of) Assumption 2 were already considered by meitz2022subgear. As many time series in economics, finance, and other fields exhibit conditional heteroskedasticity, in this paper we consider an extension to ARCH errors. In the homoskedastic case considered in meitz2022subgear, the term $g(u_{t-1})$ in ((ref)) was replaced with the more general formulation $u_{t-1}+\tilde{g}(y_{t-1},\ldots,y_{t-p})$ (with $\tilde{g}$ a real-valued function) to allow for more general dependence on the past through the variables $y_{t-1},\ldots,y_{t-p}$ (and not only through the linear combination $u_{t-1}=y_{t-1}-\pi_{1}y_{t-2}-\cdots-\pi_{p-1}y_{t-p}$). The present simpler formulation worked well in the empirical application of Section 4 and in some other examples we tried out, and leads to more transparent assumptions and streamlined proofs.
To establish ergodicity, we need the companion form of the $(p+q)$-dimensional process $\boldsymbol{y}_{t}=(\boldsymbol{y}_{1,t},\boldsymbol{y}_{2,t})$ with a $p$-dimensional $\boldsymbol{y}_{1,t}=(y_{t},\ldots,y_{t-p+1})$ and a $q$-dimensional $\boldsymbol{y}_{2,t}=(y_{t-p},\ldots,y_{t-p-q+1})$. First we formulate the $p$-dimensional companion form related to equation ((ref)), which reads as \[ \left[
\right]=
\left[
\right]+g(u_{t-1})\left[
\right]+\sigma_{t}\varepsilon_{t}\left[
\right] \] or, denoting the matrix in this equation with $\boldsymbol{\Phi}$ and setting $\boldsymbol{\iota}_{p}=(1,0,\ldots,0)$ ($p\times1$), as
(when $p=1$, $\boldsymbol{\Phi}=0$ and $u_{t-1}=y_{t-1}$). As $\sigma_{t}=\sigma(\boldsymbol{y}_{t-1})$ depends on the whole $(p+q)$-dimensional vector $\boldsymbol{y}_{t-1}$, we have to expand ((ref)) to the $(p+q)$-dimensional companion form
where $I_{q}$ is the $(q\times q)$ identity matrix and $\boldsymbol{0}_{*\times*}$ denotes a matrix of zeros with the indicated dimensions (and $\boldsymbol{\iota}_{p+q}$ is defined in the obvious way). This shows that $\boldsymbol{y}_{t}$ is a Markov chain on $\mathbb{R}^{p+q}$.
In order to establish ergodicity we further transform the $p$-dimensional companion form ((ref)) in a way already used in meitz2022subgear. 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 equation ((ref)) can be transformed into
where $\mathbf{A}\boldsymbol{y}_{1,t}=(u_{t},y_{t-1},\ldots,y_{t-p+1})$. Now, for any $p$-dimensional vector $\boldsymbol{x}_{1}$, form the partition $\boldsymbol{x}_{1}=(x_{1,1},\ldots,x_{1,p})=(x_{1,1},\boldsymbol{x}_{1,2})$ and define
(when $p=1$, $\boldsymbol{x}_{1,2}$ and $\boldsymbol{z}_{2}(\boldsymbol{x}_{1})$ are dropped). Using this notation equation ((ref)) can be expressed as $\boldsymbol{z}(\boldsymbol{y}_{1,t})=\mathbf{\Pi}\boldsymbol{z}(\boldsymbol{y}_{1,t-1})+g(z_{1}(\boldsymbol{y}_{1,t-1}))\boldsymbol{\iota}_{p}+\sigma(\boldsymbol{y}_{t-1})\varepsilon_{t}\boldsymbol{\iota}_{p}$, that is, as
Here the first equation is in a form where the autoregressive order is one and the volatility term is a function of the $(p+q)$-dimensional vector $\boldsymbol{y}_{t-1}=(\boldsymbol{y}_{1,t-1},\boldsymbol{y}_{2,t-1})$ whereas the second equation involves the $p$-dimensional vector $\boldsymbol{y}_{1,t-1}$ only.
By Assumption 2(i), the roots of the polynomial $\varpi(z)$ lie outside the unit circle, so that the eigenvalues of the matrix $\boldsymbol{\Pi}_{1}$ in the second equation in ((ref)) are smaller than one in absolute value. As is well known, this implies the existence of a matrix norm of $\boldsymbol{\Pi}_{1}$ that is also smaller than one. Specifically, for any vector norm $\Vert\cdot\Vert$, denote by ${\lvert\kern-0.25ex\lvert\kern-0.25ex\lvert \,\cdot\, \rvert\kern-0.25ex\rvert\kern-0.25ex\rvert}$ the corresponding induced matrix norm (HornJohnson2013); that is, for any conformable square matrix $A$, set \[ {\lvert\kern-0.25ex\lvert\kern-0.25ex\lvert A \rvert\kern-0.25ex\rvert\kern-0.25ex\rvert}=\max_{\Vert\boldsymbol{x}\Vert=1}\Vert A\boldsymbol{x}\Vert. \] Then we obtain the following result (HornJohnson2013).
The existence of an induced matrix norm with the property in the above lemma is essential in our proofs. (When $p=1$, Assumption 2(i) and Lemma 1 are redundant.) The norms $\Vert\cdot\Vert_{*}$ and ${\lvert\kern-0.25ex\lvert\kern-0.25ex\lvert \,\cdot\, \rvert\kern-0.25ex\rvert\kern-0.25ex\rvert}_{*}$ are defined on $\mathbb{R}^{p-1}$ and $\mathbb{R}^{(p-1)\times(p-1)}$, respectively, and they have been commonly used in time series models. In the next subsection we introduce norms which are of a different type.
The root condition of Assumption 2(i) and inequality ((ref)) of Assumption 2(ii) are of major importance for establishing the stability of our model. However, as these conditions only concern the conditional mean, we need additional assumptions restricting the conditional variance $\sigma_{t}^{2}$. As an extension of the basic ARCH model ((ref)) we consider a nonlinear formulation of the conditional variance defined as
where $\zeta_{i,t-1}=\zeta_{i}(\boldsymbol{y}_{t-1})$ is a function of $\boldsymbol{y}_{t-1}$ ($i=0,\ldots,q$) and otherwise the notation is as in equation ((ref)) (including the conditions $\omega>0$ and $\alpha_{1},\ldots,\alpha_{q}\geq0$). When the functions $\zeta_{i,t-1}$ are the same for all $i=0,\ldots,q$ we remove the index $i$ and use the notations $\zeta_{t-1}$ and $\zeta(\cdot)$. This is the case in our empirical example where $\zeta_{t-1}=\zeta(y_{t-1})=1/(1+e^{-\gamma(y_{t-1}-a)})$ is a logistic function depending only on $y_{t-1}$. For possible alternatives we consider a more general formulation and introduce the following assumption.
The above assumption includes the case $\zeta_{i}\equiv1$ for all $i$, which corresponds to the linear ARCH model ((ref)). It covers also the above-mentioned logistic function.
Consider the $q$-dimensional process $\boldsymbol{\xi}_{t}=(e_{t}^{2},e_{t-1}^{2},\ldots,e_{t-q+1}^{2})$ ($t\geq1$) with initial values $\boldsymbol{\xi}_{0}=(e_{0}^{2},\ldots,e_{-q+1}^{2})$ where $e_{0}^{2},\ldots,e_{-q+1}^{2}$ are functions of $\boldsymbol{y}_{0}$. Inspired by cline2004stability we now introduce the following equation which is a straightforward implication of equation ((ref)) and the fact $\sigma_{t}\varepsilon_{t}=e_{t}$: \[ {3pt}\left[
\right]=
\left[
\right]+\left[
\right] \] ($t=1,2,\ldots$) or, more briefly,
as $\boldsymbol{\xi}_{t}$ is a function of $\boldsymbol{y}_{t}$, we occasionally write $\boldsymbol{\xi}_{t}=\boldsymbol{\xi}(\boldsymbol{y}_{t})$. For later purposes we also note that due to the identities $\sigma_{t}\varepsilon_{t}=e_{t}$ and $e_{t}^{2}=\boldsymbol{\iota}'_{q}\boldsymbol{\xi}(\boldsymbol{y}_{t})$ we have
When there is need to make the dependence of $\Lambda_{\zeta,t}$ on $\boldsymbol{y}_{t-1}$ explicit we use the notation $\Lambda_{\zeta,t}(\boldsymbol{y}_{t-1})$ and replace the (random) argument $\boldsymbol{y}_{t-1}$ by a fixed counterpart when needed. Specifically, $\Lambda_{\zeta,t}(\boldsymbol{x})$ means that the functions $\zeta_{i,t-1}=\zeta_{i,t-1}(\boldsymbol{y}_{t-1})$ used in $\Lambda_{\zeta,t}(\boldsymbol{y}_{t-1})$ are replaced by $\zeta_{i,t}(\boldsymbol{x})$ for all $i=1,\ldots,q$, and the notations $\sigma^{2}(\boldsymbol{x})$ and $\boldsymbol{\xi}(\boldsymbol{x})$ are used similarly.
We also define the matrices
Note that $\Lambda_{t}$ is obtained from the matrix $\Lambda_{\zeta,t}$ by choosing $\zeta_{i,t}=1$ for all $i=1,\ldots,q$. Similarly, we denote $\boldsymbol{\omega}_{t}=(\omega_{t},0,\ldots,0)'$ with $\omega_{t}=\omega\varepsilon_{t}^{2}$ and $\boldsymbol{\omega}=E[\boldsymbol{\omega}_{t}]=(\omega,0,\ldots,0)'$.
In our proofs we need to appropriately control the size of the random matrix $\Lambda_{t}$, and not just the size of the non-random matrix $\Lambda=E[\Lambda_{t}]$. This is the reason why we next consider vector and matrix norms more complicated than those in Lemma 1. To this end, we first recall the definition of an $L$$^{p}$-norm (for convenience, in this subsection only, we use the notation $p$ in $L^{p}$-norms; elsewhere in the paper $p$ stands for the autoregressive order in model ((ref))). If $\Vert\cdot\Vert$ is any vector norm on $\mathbb{R}^{q}$ and $\boldsymbol{v}$ is a $q$-dimensional random vector, equation \[ \Vert\boldsymbol{v}\Vert_{L^{p}}=(E[\Vert\boldsymbol{v}\Vert^{p}])^{1/p}\qquad(1\leq p<\infty) \] defines an $L^{p}$-norm on the set of (equivalence classes of almost surely equal) $q\times1$ random vectors that are $p$-integrable (see, e.g., Dudley2004). It may be worth noting that for nonrandom vectors there is no difference between the norms $\Vert\cdot\Vert$ and $\Vert\cdot\Vert_{L^{p}}$ but for random vectors the outcome of $\Vert\cdot\Vert$ is random and that of $\Vert\cdot\Vert_{L^{p}}$ is nonrandom. This $L^{p}$-norm can be used to induce a norm for random matrices; for the conventional non-random matrix case and for the terminology used below, see HornJohnson2013. Specifically, to define a generalized (non-submultiplicative) matrix norm ${\lvert\kern-0.25ex\lvert\kern-0.25ex\lvert \,\cdot\, \rvert\kern-0.25ex\rvert\kern-0.25ex\rvert}_{L^{p}}$, for any $q\times q$ random matrix $A$ set
where the latter equality holds as $\boldsymbol{x}$ is nonrandom. This defines a generalized matrix norm\footnote{Axioms (1), (1a), (2), and (3) of a generalized matrix norm (see HornJohnson2013) can be checked similarly as in the proof of Theorem 5.6.2(c) of the same reference (replacing the norms $\Vert\cdot\Vert$ and ${\lvert\kern-0.25ex\lvert\kern-0.25ex\lvert \,\cdot\, \rvert\kern-0.25ex\rvert\kern-0.25ex\rvert}$ therein with $\Vert\cdot\Vert_{L^{p}}$ and ${\lvert\kern-0.25ex\lvert\kern-0.25ex\lvert \,\cdot\, \rvert\kern-0.25ex\rvert\kern-0.25ex\rvert}_{L^{p}}$, replacing appropriate statements therein with their almost sure counterparts, and using Minkowski's inequality as an additional justification for axiom (3)).} on the set of (equivalence classes of almost surely equal) $q\times q$ random matrices with $p$-integrable entries; moreover, the norms $\Vert\cdot\Vert$, $\Vert\cdot\Vert_{L^{p}}$, and ${\lvert\kern-0.25ex\lvert\kern-0.25ex\lvert \,\cdot\, \rvert\kern-0.25ex\rvert\kern-0.25ex\rvert}_{L^{p}}$ are related by the inequality\footnote{Inequality ((ref)) can be verified analogously to Theorem 5.6.2(b) of HornJohnson2013.}
We next state a high-level condition that assumes the existence of a vector norm on $\mathbb{R}^{q}$ with particular additional properties. (Primitive conditions ensuring this high-level assumption will be given momentarily.) One of these properties is monotonicity in the sense of Definition 5.4.18 of HornJohnson2013: a vector norm $\Vert\cdot\Vert$ is monotone if $\boldsymbol{x},\boldsymbol{y}\in\mathbb{R}^{q}$ satisfying $|x_{i}|\leq|y_{i}|$ for $i=1,\ldots,q$ always implies that $\Vert\boldsymbol{x}\Vert\leq\Vert\boldsymbol{y}\Vert$. For clarity, we use the notation $\Vert\cdot\Vert_{\bullet}$ for the specific vector norm in the assumption below; similarly, we denote the related $L^{p}$-norm by $\Vert\cdot\Vert_{\bullet L^{p}}$ and the generalized matrix norm by ${\lvert\kern-0.25ex\lvert\kern-0.25ex\lvert \,\cdot\, \rvert\kern-0.25ex\rvert\kern-0.25ex\rvert}_{\bullet L^{p}}$. We also introduce two constants, $s_{0}\geq1$ and $b\geq1$, such that \[ b=1\text{ when }s_{0}=1\quad\text{and}\quad b>(2s_{0}-\rho)/[s_{0}(2-\rho)]>1\text{ when }s_{0}>1 \] (recall from Assumption 2 that $\rho\in(0,2)$ so that $2s_{0}>\rho$). These constants are used in the next section where we establish our ergodicity result and there the size of $s_{0}$ will have an effect on the rate of convergence obtained and the order of moments that are finite. The rather complex conditions required from the constant $b$ are due to the connection between the conditional mean and ARCH errors (this connection disappears when $s_{0}=1$ as it also does in subgeometric homoskedastic autoregressions).
This assumption tacitly requires that $E[|\varepsilon_{t}|^{2bs_{0}}]$ is finite, thereby strengthening Assumption 1 (when $s_{0}>1$). Assumption 4 is formulated in a way that is convenient in our proofs but is not very transparent. The following lemma gives primitive conditions ensuring that Assumption 4 holds (for a proof, see the appendix).
To illustrate, consider the case $q=s_{0}=b=1$ so that the condition in Lemma 2 reduces to the requirement $\alpha_{1}<1$. In geometrically ergodic AR models with linear ARCH(1) errors, $\alpha_{1}<1$ is the usual requirement for covariance stationarity while geometric ergodicity can hold under even weaker conditions, such as $E[\ln(\alpha_{1}\varepsilon_{t}^{2})]<0$ (see, e.g., meitz2010SPL for further details). In the present setting the situation is different: as will be seen in Section 3.2, condition $\alpha_{1}<1$ does not guarantee a finite variance.
We now consider the stability of the model introduced in the previous section. We begin with a brief account of some necessary Markov chain concepts (for more comprehensive discussions, see meyn2009markov and douc2018markov and also meitz2022subgear). 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. For an arbitrary fixed measurable function $f:\mathsf{X}\rightarrow[1,\infty)$ 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))). When 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, \footnotetext[\numexpr\thefootnote+1\relax]{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$.}
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$. In the conventional geometrically ergodic case, $r(n)=r^{n}$ for some $r>1$. To establish ($f,r$)-ergodicity, we use a so-called drift condition defined as follows (here $\boldsymbol{1}_{S}(x)$ denotes the indicator function taking value one when $x$ belongs to the set $S$ and zero elsewhere).
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 $\tilde{b}$ such that
The idea is to verify this condition with suitable functions $V$ and $\phi$, which together with some additional conditions ensures the ($f,r$)-ergodicity of the process $X_{t}$; for more details, see meitz2022subgear.
Now consider the stability of the Markov chain $\boldsymbol{y}_{t}$ on $\mathbb{R}^{p+q}$ given in ((ref)). To define the function $V$ in ((ref)), we use the functions $z_{1}(\cdot)$, $\boldsymbol{z}_{2}(\cdot)$, and $\boldsymbol{\xi}(\cdot)$ in ((ref))\textendash ((ref)) and ((ref)) and the norms $\lVert\cdot\rVert_{*}$ and $\lVert\cdot\rVert_{\bullet}$ in Lemma 1 and Assumption 4. Set $\boldsymbol{x}=(x_{1},\ldots,x_{p+q})\in\mathbb{R}^{p+q}$ and decompose $\boldsymbol{x}$ to its $p$- and $q$-dimensional components as $\boldsymbol{x}=(\boldsymbol{x}_{1},\boldsymbol{x}_{2})$. We define the function $V$ as
where $s_{0}$ and $b$ are defined above Assumption 4, $s_{1}$ and $s_{2}$ are positive constants to be specified later (with $s_{1}$ small and $s_{2}$ large), and $\alpha=1-\rho/2s_{0}$ (recall from Assumption 2 that $\rho\in(0,2)$ so that $\alpha\in(0,1)$). It may be clarifying to note that when $p=1$, model ((ref)) reduces to $y_{t}=g(y_{t-1})+\sigma_{t}\varepsilon_{t}$; then we can set $s_{1}=0$ and drop $\boldsymbol{z}_{2}(\boldsymbol{x}_{1})$ so that the function $V$ in ((ref)) becomes $V(\boldsymbol{x})=1+|x_{1}|^{2s_{0}}+s_{2}\lVert\boldsymbol{\xi}(\boldsymbol{x})\rVert_{\bullet}^{bs_{0}}$.
To verify Condition D, we need to consider the conditional expectation
bound the conditional expectations on the right hand side of ((ref)), and express these bounds in a way which conforms to inequality ((ref)) with the function $\phi$ satisfying the conditions required in Condition D. These considerations, combined with the checking of some additional technical conditions, lead to the following theorem (the proof can be found in the Appendix).
Theorem 1 provides the first subgeometric ergodicity results for autoregressions with autoregressive conditional heteroskedasticity. In this theorem, the convergence rate $r(n)$ shows the speed at which the $n$-step transition probability measures of the process $\boldsymbol{y}_{t}$ converge to the stationary probability measure. Due to the polynomial convergence rate we therefore call the process $\boldsymbol{y}_{t}$ polynomially ergodic. Note also that the choice of $\delta$ in Theorem 1 allows for a trade-off between the rate of convergence and the size of the $f$-norm.
\paragraph{Geometric ergodicity.}
In previous literature, geometric ergodicity of nonlinear autoregressions with ARCH errors has been considered using a variety of different assumptions for the allowed nonlinear dynamics and for the required moment conditions for the innovations; see, e.g., cline2004stability, meitz2010SPL, and the many references therein.
\paragraph{Homoskedastic case.}
Theorem 1 remains valid also in the homoskedastic case (obtained by setting $\alpha_{1}=\cdots=\alpha_{q}=0$). Previous polynomial ergodicity results for homoskedastic autogressions were obtained by fort2003polynomial and meitz2022subgear, and the above Theorem 1 provides partial improvements over these earlier results in certain cases. Assumptions and notation are slightly different in all the papers, but (in the notation of the present paper) Theorem 1 improves earlier results when $1\leq\rho<2$ and $1<s_{0}<2$.
\paragraph{Proof strategy.}
The proof of Theorem 1 is also somewhat different from the previous polynomial ergodicity results in fort2003polynomial and meitz2022subgear. A rather obvious difference is that these earlier results deal with homoskedastic autoregressions whereas our model contains a nonlinear ARCH term, the size of which is controlled with the special matrix norm defined in Assumption 4. Regarding the conditional expectation, the mentioned earlier results rely on Lemma 3 in fort2003polynomial while our proof of Theorem 1 avoids the use of this lemma, and instead makes use of a (relatively unknown) extension of Bernoulli\textquoteright s inequality due to FeffermanShapiro1972 (combined with Young's inequality).
\paragraph{Mixing and moment results. }
As already indicated in the Introduction, the polynomial ergodicity result of Theorem 1 also implies that the process $\boldsymbol{y}_{t}$ is $\beta$-mixing (and hence $\alpha$-mixing). Moreover, the convergence rate of the $\beta$-mixing coefficients $\beta(n)$ is given by the fastest convergence rate, that is, $\lim_{n\rightarrow\infty}n^{2s_{0}/\rho-1}\beta(n)=0$. For further details and justifications of these mixing results, see meitz2021JAP and meitz2022subgear.
Another consequence of Theorem 1 is that the stationary distribution of $\boldsymbol{y}_{t}$ has finite moments up to order $2s_{0}-\rho$ (for a proof, see the Appendix). Note that depending on the values of $s_{0}\geq1$ and $\rho\in(0,2)$, the order of these finite moments may be very small; in particular, when $s_{0}=1$ we do not obtain a finite variance.
\paragraph{Subexponential ergodicity.}
Theorem 1 concerns only polynomial ergodicity of subgeometric AR\textendash ARCH models, and does not consider subexponential ergodicity (where the rate $r(n)$ in ((ref)) equals, say, $e^{cn^{\gamma}}$ with $c>0$ and $0<\gamma<1$). The reason for this is that the properties of ARCH-type models do not seem compatible with the moment requirements needed for subexponential ergodicity. To elaborate on this, first note that the previous results of douc2004practical and meitz2022subgear on subexponential ergodicity of homoskedastic nonlinear autoregressions (i) require the IID error term to possess moments of all orders and (ii) imply that the observed process $y_{t}$ also has finite moments of all orders. (To provide some further details, (ii) is given as Corollary to Theorem 2 in meitz2022subgear. As for (i), see Assumptions 3.3 and 2(a) of douc2004practical and meitz2022subgear, respectively. These assumptions require the IID error terms to be sub-Weibull random variables, which in turn entails they possess moments of all orders; see vladimirova2020sub or wong2020lasso.)
The abovementioned moment requirements are in stark contrast to ARCH-type models. For instance, in the simplest ARCH(1) model ($e_{t}=\sigma_{t}\varepsilon_{t}$, $\sigma_{t}^{2}=\omega+\alpha_{1}e_{t-1}^{2}$, and $\varepsilon_{t}$ IID N(0,1)), the finiteness of moments of order $2r$ for $e_{t}$ ($E[|e_{t}|^{2r}]<\infty$) is known to require the condition $\alpha_{1}^{r}E[|\varepsilon_{t}|^{2r}]<1$ (see, e.g., ling2002necessary and ling1999probabilistic). For integer values of $r$, this condition is equivalent with $\alpha_{1}<[(2r-1)!!]^{-1/r}=[1\cdot3\cdot\ldots\cdot(2r-1)]^{-1/r}$ and consequently all moments of the ARCH process $e_{t}$ cannot be finite unless $\alpha_{1}=0$. The situation is similar also in more complicated (G)ARCH and AR\textendash (G)ARCH models (see, e.g., meitz2008ergodicity and meitz2008stability, respectively). This suggests that ARCH-type heteroskedastic errors may not be compatible with the moment requirements needed for subexponential ergodicity.
\paragraph{Potential extensions.}
Extending our results to allow for GARCH (and not only ARCH) errors would be interesting. However, previous literature suggests that studying the stability of nonlinear AR\textendash GARCH models can be challenging. Geometric ergodicity of nonlinear AR\textendash GARCH models has previously been studied by \citet*{liu1997threshold}, ling1999probabilistic, cline2007stability, and meitz2008stability; of these articles, the former two are confined to threshold AR\textendash GARCH models, whereas the latter two consider more general nonlinear autoregressions. In the present setting, the autoregressive part of the model we consider is rather general (the restrictions imposed on function $g(\cdot)$ in Assumption 2(ii) are quite mild, essentially restricting $g(\cdot)$ only for large values of its argument) and techniques used for threshold models can not be applied. Using an approach similar to Cline's (2007) appears challenging as the assumptions he employs are quite general and appear difficult to verify (in fact, a threshold AR\textendash GARCH model is the only example that is explicitly treated in his article). On the other hand, meitz2008stability require certain structure and smoothness of the conditional mean (see Assumption 2 of their paper) and it is not clear how to apply these results in the current setting. As the extension to GARCH errors appears challenging, we leave it for future research.
Another useful extension would be to consider the subgeometric ergodicity of multivariate autoregressions with autoregressive conditional heteroskedasticity. fort2003polynomial and douc2004practical already studied multivariate first-order autoregressions with IID errors and obtained results for polynomial and subexponential ergodicity, respectively. In principle, generalizing these results to the higher-order case with multivariate ARCH errors should be possible but it is not immediate how to formulate a general model that would be both theoretically manageable as well as useful in practical applications. We hope to return to this issue in subsequent work.
The conditional mean of the model we have so far discussed is very general, and we next consider some concrete illustrating examples. The following two special cases were introduced in meitz2022subgear in the case of a homoskedastic error term. We first consider a model with a time-varying intercept term based on a logistic function and specified as
where $L(u;\gamma,a)=1/(1+e^{-\gamma(u-a)})$ is the logistic function and the parameters $\gamma$, $a_{1}$, $a_{2}$ are assumed to satisfy $\gamma>0$ and $a_{1}\leq a_{2}$, and $\nu_{1}$, $\nu_{2}$ are assumed to satisfy $\nu_{1}<0<\nu_{2}$. Moreover, $\Delta$ signifies the difference operator (so that $\Delta y_{t-1}=y_{t-1}-y_{t-2}$) and the remaining notation is as in model ((ref)). Arguments similar to those in meitz2022subgear can now be used to prove the following result (for details, see the Appendix).
The convergence rate presented in Proposition 1 also shows the rate of $\beta$-mixing coefficients.
As another special case, we consider a model with a time-varying slope term defined as
where $S(u_{t-1})$ is either $S_{1}(u_{t-1})=1-r_{0}/h(u_{t-1})$ or $S_{2}(u_{t-1})=\exp\{-r_{0}/h(u_{t-1})\}$ (with $r_{0}>0$) and the function $h:\mathbb{R}\rightarrow(0,\infty)$ as defined in Proposition 2 of meitz2022subgear. In addition to a general formulation of the function $h$ that proposition provides six special cases of which two are $h(u)=1+|u-a|^{\rho}$ and $h(u)=(1+(u-a)^{2})^{\rho/2}$ (where $a\in\mathbb{R}$ and $\rho\in(0,2)$; see Assumption 2). Regarding the remaining notation, it is as in model ((ref)).
The following result can be established by using arguments similar to those in the proof of Proposition 2 in meitz2022subgear (for details, see the Appendix).
The rate of $\beta$-mixing coefficients coincides with the rate given in the proposition. As the function $h$ depends on the parameter $\rho\in(0,2)$, the convergence rate in Proposition 2 differs from that obtained in Proposition 1 except in the case $\rho=1$.
Although theoretical work on subgeometric ergodicity has been ongoing for four decades, practical illustrations of (homoskedastic) subgeometrically ergodic autoregressions have been scarce; we are not aware of any previous empirical applications of subgeometrically ergodic autoregressions using real data. A small illustration of simulated data from one subgeometrically ergodic autoregression is given in fort2003polynomial. meitz2022subgear provide examples of some concrete subgeometrically ergodic autoregressive time series models and illustrations of a few simulated data series from them. These simulation exercises suggest that subgeometrically ergodic autoregressions could be useful when the observed time series bears some resemblance to unit root type behavior and the autocorrelation function indicates very strong persistence, but when the time series nevertheless exhibits eventual mean-reverting behavior.\textcolor{red}{ }The discussion in Section 2.1 around equation ((ref)) had a similar message, suggesting these models could be seen as a middle ground between the extreme cases of geometric ergodicity and non-ergodicity.
We next illustrate the use of subgeometrically ergodic AR\textendash ARCH models in a small empirical example. Our aim is simply to provide a proof of concept for the applicability of subgeometrically ergodic AR\textendash ARCH models, illustrating that the model used fits the data well. Further work is certainly needed to judge the usefulness of these models in practical applications but we leave such more comprehensive empirical applications for future research.
The data we employ consists of daily observations on the Chicago Board Options Exchange energy sector volatility index (\url{fred.stlouisfed.org/series/VXXLECLS}) over the period 16 March 2011 through 31 December 2021 (a total of 2719 observations). This data series reflects energy sector risk and is displayed in the top left graph of Figure 1 (the solid graph; the dashed horizontal line shows the estimate $\hat{a}=25.366$, see ((ref)) and ((ref)) below). The time series plot shows signs of strong persistence, which is also reflected in the autocorrelation function of the data shown in the top right graph of Figure 1.
We model this data series using the parametric specification in ((ref)). As for the error distribution, after some experimentation a skew version of the $t$-distribution due to jones2003skew was found to provide a good fit (in contrast, estimation with normal errors lead to a distinct discrepancy between the residual distribution and the Gaussian one). The density function of this distribution is \[ f(x;c,d)=C_{c,d}^{-1}\biggl\{1+\frac{x}{(c+d+x^{\text{2}})^{1/2}}\biggr\}^{c+1/2}\biggl\{1-\frac{x}{(c+d+x^{\text{2}})^{1/2}}\biggr\}^{d+1/2}, \] where $c$ and $d$ are positive parameters and $C_{c,d}=2^{c+d-1}B(c,d)(c+d)^{\text{1/2}}$ (with $B(\cdot,\cdot)$ denoting the beta function); the case $c=d$ results in a symmetric $t$-distribution with $2c$ degrees of freedom and the cases $c<d$ and $c>d$ imply skewness to the left and right, respectively. In our application, we use this distribution centralized to have mean zero and standardized to have unit variance (i.e., in the density function $x$ is replaced by $sx+m$ and $C_{c,d}^{-1}$ by $sC_{c,d}^{-1}$ where $m$ and $s^{2}$ denote the mean and variance, see jones2003skew; this requires that $c>1$ and $d>1$, for a moment of order $k$ is finite when $c>k/2$ and $d>k/2$).
We estimate the model parameters using the method of maximum likelihood and employ optimization routines in R. (We simply assume that standard properties of maximum likelihood estimators hold and calculate standard errors based on the standard formulas.) Trying out different model orders lead to model ((ref)) with order $p=1$ and with a nonlinear ARCH term of order $q=3$ (with these choices the residual diagnostics shown in Figure 2 in Appendix B indicated a very good fit). Specifically, the considered model is
where the errors $\varepsilon_{t}$ are IID$(0,1)$ and follow the above described (centralized and standardized) skew $t$-distribution, $L(y;\gamma,a)=1/(1+e^{-\gamma(y-a)})$ is the logistic function, the parameters $\nu$ and $\gamma$ are positive, and $a\in\mathbb{R}$. (We also tried a model where the logistic functions in the conditional expectation and in the ARCH term were different but this extension had only a minor effect on the results.) ML estimation (with the constraints $\nu,\gamma,\omega>0$, $\alpha_{1},\alpha_{2},\alpha_{3}\geq0$, and $\alpha_{1}+\alpha_{2}+\alpha_{3}<1$ that ensure polynomial ergodicity by Proposition 1) leads to the following results:
where the numbers in parenthesis are standard errors; estimates for the parameters in the error distribution are $\hat{c}=3.551\:(0.422)$ and $\hat{d}=2.138\:(0.197)$.
To illustrate the conditional mean of the estimated model, consider the function $I(x)=-\nu L(x;\gamma,a)+\nu(1-L(x;\gamma,a))$ and the corresponding time-varying intercept term $I(y_{t-1})$ based on the above parameter estimates. These are shown in the middle row of Figure 1. In the left panel, the two horizontal dashed lines show the minimum and maximum $I(x)$ attains, while the three vertical dashed lines indicate the minimum of the observed data series $y_{t}$ ($11.71$), the estimate $\hat{a}=25.366$, and the maximum of $y_{t}$ ($130.61$). On the right, the three horizontal dashed lines show the minimum and maximum $I(y_{t-1})$ attains ($-0.187$ and $0.154$) and the origin. Intuitively, when $y_{t-1}$ is close to $\hat{a}$ the time-varying intercept term $I(y_{t-1})$ is close to zero and the conditional mean of ((ref)) corresponds to unit root type behavior (without drift); when $y_{t-1}$ takes values clearly below/above $\hat{a}$ the intercept $I(y_{t-1})$ is positive/negative and behavior akin to a unit root process with increasing/decreasing drift occurs.
The left panel in the bottom row of Figure 1 displays the estimated volatility series $\hat{\sigma}_{t}$. The variation of the volatility over time is strong, and the large spikes in the volatility series coincide with the large values in the observed series. The logistic formulation of the conditional variance in ((ref)) makes it possible for large observations to amplify volatility more than a standard linear ARCH model would allow for.
The right panel in the bottom row of Figure 1 shows the residual series $\hat{\varepsilon}_{t}$. Four additional graphs analyzing the residuals are available in Figure 2 in Appendix B: autocorrelation functions of the residuals and of the squared residuals, together with a histogram and a Q-Q plot. The autocorrelation functions reveal that the very strong persistence present in the original series has been quite well captured by the estimated subgeometrically ergodic AR\textendash ARCH model (only three of the shown 100 autocorrelation coefficients are barely outside the displayed critical values). The histogram and the Q-Q plot indicate that the employed skew version of the $t$-distribution fits well as only a few outlying observations deviate from the estimated density function and the 45 degree line.
Note also that the estimated AR\textendash ARCH model satisfies the requirements of a stationary and\textcolor{red}{ }polynomially ergodic process with finite absolute moments.\footnote{That is, the parameter estimates in ((ref)) correspond to a process satisfying the requirements of Proposition 1 with $s_{0}=1$. (Note that these requirements are not satisfied with $s_{0}=1.5$ which would correspond to finite second moments of $y_{t}$.)} It may be interesting to note that estimation attempts using standard linear ARMA(1,1)\textendash GARCH(1,1) models (with skew $t$ errors) lead to estimated autoregressive coefficients in excess of $0.999$, reflecting the very persistent nature of the data series apparent from the time series and autocorrelation plots in the top row of Figure 1.
The primary purpose of this small empirical example was to demonstrate what kind of time series could be modeled with subgeometrically ergodic AR\textendash ARCH models. It is worth pointing out that such models may work well even in cases where the graphs of the employed time series and related autocorrelation functions look very different from those displayed in Figure 1.
In this paper, we examined the subgeometric ergodicity of nonlinear autoregressive models with autoregressive conditional heteroskedasticity. We provided conditions that ensured polynomial ergodicity of the considered AR\textendash ARCH models. Our results generalized existing results that assumed the error terms to be IID. The use of subgeometrically ergodic AR\textendash ARCH models was illustrated in an empirical example using energy sector volatility index data.
Several future research topics could be entertained. In this paper we have only considered ARCH-type conditional heteroskedasticity, and extending the results to the generalized ARCH (GARCH) case would be of interest. Subgeometric ergodicity of multivariate autoregressions with autoregressive conditional heteroskedasticity is another interesting topic left for future work. On the empirical side, further applied work is certainly needed to judge the usefulness of subgeometrically ergodic autoregressions (with or without ARCH) in practical applications. For instance, providing more concrete advice on when to use subgeometrically (rather than geometrically) ergodic autoregressions would be useful for practitioners. Another question future applications should address is whether subgeometrically ergodic autoregressions can outperform relevant competing models in out-of-sample forecasting exercises.