On a new robust method of inference for general time series models

On a new robust method of inference for general time series models

\def\spacingset#1{ {#1}} \spacingset{1}

\if00 { \affil[1]{\it Department of Statistics and Data Science, Tsinghua University, Beijing, China} \affil[2]{\it Faculty of Business and Economics, The University of Hong Kong, Hong Kong} \affil[3]{\it Department of Statistics, London School of Economics, London, U.K.} } \fi

\if10 { \spacingset{2}

} \fi

\spacingset{1.3}

{\it Keywords:} Conditional heteroscedasticity model; Heavy-tailed innovation; Identifiability condition; Quasi-maximum likelihood estimation; Robust estimation.

\spacingset{1.7}

Introduction

Time series models play a crucial role across diverse fields, such as data science, economics, finance, engineering, environmental science, epidemiology, healthcare, hydrology, and sociology, among others. Researchers and practitioners often need to model time series data to reveal temporal dependent structures, as well as to understand dynamic mechanisms and patterns in data, to predict future trends, and to make decisions. See, for example, box2015time, tsay2010analysis, WSR2024, etc. In most practical scenarios, the modeling is based on the first two conditional moments. Specifically, a general parametric time series model bardet2009asymptotic, ling2010general takes typically the form

where ${\bf Y}_{t-1}=(y_{t-1}, y_{t-2}, \dots)^{{ \mathrm{\scriptscriptstyle T} }}\in\mathbb{R}^{\infty}$, or ${\bf Y}_{t-1}=(y_{t-1},\dots,y_{t-m})^{{ \mathrm{\scriptscriptstyle T} }}\in\mathbb{R}^m$, $\boldsymbol \theta\in\mathbb{R}^d$ is the unknown parameter vector of interest, $g({\bf Y},\boldsymbol \theta):\mathbb{R}^\infty\ ({\rm or\ }\mathbb{R}^{m})\times\mathbb{R}^d\to\mathbb{R}$ and $\sigma({\bf Y},\boldsymbol \theta):\mathbb{R}^\infty\ ({\rm or\ }\mathbb{R}^{m})\times\mathbb{R}^d\to\mathbb{R}^+$ are measurable functions, the innovation sequence $\{\eta_t\}$ is independent and identically distributed (i.i.d.) with zero mean, and $\eta_t$ is independent of $\{y_s: s<t\}$.

Model ((ref)) is quite general, containing many classical time series models as special cases. Specifically, when $\sigma({\bf Y}_{t-1},\boldsymbol \theta)$ is a fixed positive constant, i.e., $\sigma({\bf Y}_{t-1},\boldsymbol \theta)\equiv \sigma>0$, it includes the best known linear time series model, namely the autoregressive moving average (ARMA) model BrockwellDavis, box2015time. When $g({\bf Y}_{t-1},\boldsymbol \theta)\equiv 0$, model (ref) reduces to a general conditionally heteroscedastic time series model

which includes the celebrated autoregressive conditional heteroscedasticity (ARCH) model Engle1982 and the generalized ARCH (GARCH) model Bollerslev1986. A fairly comprehensive review of GARCH models is available in francq2019garch.

In the literature, the Gaussian quasi-maximum likelihood estimation (GQMLE) is widely used for statistical inference of parametric time series models: for example, bardet2009asymptotic and ling2010general for model ((ref)), and RZ2006 for model ((ref)). Apart from these, there are numerous studies on the GQMLE of many specific parametric time series models, including linear and nonlinear ones. For instance, the GQMLE (asymptotically equivalent to the least squares estimator and the Whittle estimator) of ARMA models has been frequently studied; see, e.g., BrockwellDavis, Straumannbook, yao2006gaussian, boubacar2018diagnostic, and wilms2023sparse, among others. For GARCH models, their GQMLE was studied by Lumsdaine1996, HallYao, JensenRahbek, francq2004maximum, FZ2012, straumann2006quasi, jiang2021adaptive, yu2024matrix, etc. francq2019garch provided a quite comprehensive review of the GQMLE of various conditionally heteroscedastic time series models. For the GQMLE of ARMA-GARCH models, see ling1997fractionally, ling2003asymptotic, francq2004maximum, Ling2007 and the references therein.

Although robust in respect of distributional misspecification of the innovation to some extent, the GQMLE often results in a loss of efficiency francq2011two, fan2014quasi, particularly for heavy-tailed time series models.\footnote{Heavy-tailed time series data are ubiquitous in the real world. The modeling of heavy-tailed time series data is a long-lasting topic that has been attracting much attention, resulting in a growing number of monographs; see, e.g., IIWbook, BDMbook, PengQi, KulikSoulier, Nolanbook, NWZbook, and Leipus, among others.} Substantial documented evidence has demonstrated heavy-tailness of the innovation in financial returns, thus ruling out Gaussianity BW1992, Bai2003, li2023maximum. To improve the efficiency of the GQMLE, many remedial methods were proposed for some specific parametric time series models, including the least absolute deviation estimation li2008least, zhu2015lade,zhu2019statistical, zhang2022lade, the conditional quantile estimation li2015quantile,zheng2018hybrid,zhu2018linear, wang2022hybrid,zhu2023quantile, and general $M$-type estimation mukherjee2020bootstrapping, etc. A selective overview of contemporary robust statistics is available in loh2025theoretical. However, these robust estimation methods are often limited in the sense of Huber due to the non-smoothness of the objective function and the unboundedness of the influence function wooldridge2020consistency, and therefore still incur potential loss of efficiency and poor performance to a certain extent.

To overcome the shortcomings, wooldridge2020consistency proposed an ingenious logistic QMLE (LQMLE) similar in principle to the GQMLE. Unfortunately, he focused exclusively on the linear regression settings, omitting completely the time series settings. Perhaps the omission is due to the recognition of a serious obstacle to do with the identifiability of the parameters, that is relevant for time series analysis, but not for linear regression. In this article, after overcoming the obstacle, we extend the LQMLE to cover the general time series model ((ref)) and develop a comprehensive procedure based on LQMLE. As far as we know, this article is the first one that fully explores the LQMLE in the context of time series analysis.

We shall assume a different moment condition on the innovation rather than the usual condition $\mathbb{E}\eta_t^2=1$. BerkesHorvath made a similar point when studying the non-Gaussian QMLE of GARCH models. First, we define

where $F(x)=1/\{1+\exp(-x)\}$ is the cumulative distribution function of the standard logistic distribution. Clearly, $h(x)$ is nonnegative and even; see Figure (ref)(a). We claim that all parameters are identifiable if $\psi(\eta_t)=1,$ under which the LQMLE is strongly consistent with some extra mild assumptions. Lemma {\color{blue}S.4} of the supplementary material justifies the use of $\psi(\eta_t)=1$. Such a condition relaxes the restriction of $\mathbb{E}\eta_t^2<\infty$ to $\mathbb{E}|\eta_t|<\infty$ since $2F(x)-1\in[-1,1]$ for all $ x\in\mathbb{R}.$ Further, it is worth noting that the conditional variance of $y_t$ (if exists) is proportional to $\sigma^2({\bf Y}_{t-1},\boldsymbol \theta)$, without assuming that $\mathbb{E}\eta_t^2=1$. This means that $\sigma^2({\bf Y}_{t-1},\boldsymbol \theta)$ can still be interpreted as the volatility of models (ref)-(ref), similar to the interpretation of re-parameterized GARCH models in fan2014quasi. For an intuitive understanding of $\psi(\eta_t)=1$, we give an example as follows.

The main contributions of this article are fourfold.

First, we propose a novel LQMLE for general parametric time series model (ref). Under some mild conditions, it is shown that the LQMLE is strongly consistent and asymptotically normal. Compared to the GQMLE and existing robust estimations in the sense of Huber, our LQMLE enjoys the following advantages: (i) it is robust in respect of distributional misspecification and heavy-tailedness of the innovation; (ii) it is more resilient to outliers; (iii) the log logistic quasi-likelihood function is smooth and strictly concave; and (iv) the influence function is bounded. These appealing properties simplify computation and statistical inference.

Second, we provide a new condition for model parameter identifiability, namely $\psi(\eta_t)=1$, which relaxes the restriction $\mathbb{E}\eta_t^2=1$ that is almost routinely used in the GQMLE. The new condition allows for heavy-tailed distributions such as the $\alpha$-stable distributions with $\alpha\in(1, 2)$ and Student's $t_{\nu}$-distributions with degrees of freedom $\nu\in(1, 2]$. Additionally, for the asymptotic normality of the LQMLE, it suffices to assume a finite second moment of the innovation, equivalent to Assumption (ref). This suggests that the innovation following Student's $t_{\nu}$-distribution with $\nu\in(2, 4]$ is also admissible.

Third, based on the LQMLE, we study the Wald test and the Lagrange multiplier test for the unknown parameters, and develop their limiting distributions. Leveraging the smoothness and strict concavity of the log logistic quasi-likelihood function, we can construct consistent estimators of asymptotic covariance matrices in the above test statistics with ease.

Last, to illustrate the applicability of the proposed methodology, we verify technical assumptions for several nonlinear time series models, namely the DAR, GARCH, ARMA-GARCH, DTARMACH, and EXPAR models. Meanwhile, Monte Carlo simulation studies are conducted to examine the performance of our methodology and an empirical example is analyzed.

The remainder of this article is organized as follows. Section (ref) presents the estimation and inference methodology as well as testing, and establishes related theoretical results. Section (ref) provides concrete applications of our methodology in several specific nonlinear time series models. Section (ref) assesses the finite-sample performance of our methodology by Monte Carlo simulation studies. Section (ref) gives an empirical study on treasury yield curve rates. Section (ref) concludes. All proofs of main theoretical results as well as additional simulation results are relegated to the supplementary material.

{\bf Notations}. For any vector ${\bf a},$ we let $\Vert{\bf a}\Vert=\sqrt{{\bf a}^{{ \mathrm{\scriptscriptstyle T} }}{\bf a}}$. $\boldsymbol 0_d$ is a $d$-dimensional vector of zeros. For any matrix ${\bf A},$ we let $\Vert{\bf A}\Vert=\sqrt{\lambda_{\max}({\bf A}^{{ \mathrm{\scriptscriptstyle T} }}{\bf A})}$, where $\lambda_{\max}(\cdot)$ denotes the largest eigenvalue of a matrix, and let $\mathrm{rank}({\bf A})$ be the rank of ${\bf A}$. $\mathbb{R}^k$ is the $k$-dimensional Euclidean space ($1\leq k\leq \infty$), and we write $\mathbb{R}=\mathbb{R}^1$ and $\mathbb{R}^+=(0, \infty)$. $\mathbb{Z}=\{0,\pm1,\pm2,...\}$ and $\mathbb{C}$ denote the sets of integer and complex numbers, respectively. For $x,y \in {\mathbb R},$ we use $x \wedge y = \min(x,y).$ For two positive sequences $\{a_n\}$ and $\{b_n\}$, we write $a_n\lesssim b_n$ or $a_n=O(b_n)$ or $b_n\gtrsim a_n$ if there exists a positive constant $c$ such that $a_n/b_n \leq c$. We write $a_n \asymp b_n$ if and only if $a_n \lesssim b_n$ and $b_n\lesssim a_n$ hold simultaneously. The symbols `$\to_p$' and `$\to_d$' mean convergence in probability and convergence in distribution, respectively. $\mathcal{F}_{t}$ denotes the sigma-algebra generated by random variables $\{y_j: j\leq t\}$ for each $t$, i.e., $\mathcal{F}_{t}=\sigma(\{y_j: j\leq t\})$. For a given random variable $X$, $X\in \mathcal{F}_{t}$ means that $X$ is measurable with respect to (w.r.t.) $\mathcal{F}_{t}$.

Methodology

Logistic QMLE with Asymptotics

Let $\boldsymbol \theta$ be the parameter and $\boldsymbol \Theta$ be the parameter space. Suppose that the observations $\{y_1,\dots,y_n\}$ are from model (ref) with the true parameter $\boldsymbol \theta_0.$ To handle the initial value problem involved in optimizing the objective function ((ref)) below, we first let $\widetilde{{\bf Y}}_0=(\widetilde{y}_0,\widetilde{y}_{-1},\dots)^{{ \mathrm{\scriptscriptstyle T} }}$ be an initial value and then define $\widetilde{{\bf Y}}_{t-1}=(y_{t-1},\dots,y_1, \widetilde{{\bf Y}}_0^{{ \mathrm{\scriptscriptstyle T} }})^{{ \mathrm{\scriptscriptstyle T} }}$ for $t\geq1$. Denote $\widetilde{g}_t(\boldsymbol \theta)=g(\widetilde{{\bf Y}}_{t-1},\boldsymbol \theta)$ and $\widetilde{\sigma}_t(\boldsymbol \theta)=\sigma(\widetilde{{\bf Y}}_{t-1},\boldsymbol \theta)$ for $t\geq 1$. The (conditional) log logistic quasi-likelihood function is defined as

where $f(x)$ is the density of the standard logistic distribution, i.e.,

The LQMLE of $\boldsymbol \theta_0$ is defined as

To facilitate the study on the asymptotic properties of $\widehat{\boldsymbol \theta}_n$, we define the theoretical log logistic quasi-likelihood function as

where $g_{t}(\boldsymbol \theta)=g({\bf Y}_{t-1},\boldsymbol \theta)$ and $\sigma_{t}(\boldsymbol \theta)=\sigma({\bf Y}_{t-1},\boldsymbol \theta)$ for simplicity.

To obtain the strong consistency of $\widehat{\boldsymbol \theta}_n$, the following assumptions are needed.

The following theorem states the strong consistency of $\widehat{\boldsymbol \theta}_n$.

To obtain the asymptotic distribution of $\widehat{\boldsymbol \theta}_n$, we define {

} where $\dot{{\bf g}}(\boldsymbol \theta)=\partial g_t(\boldsymbol \theta)/\partial\boldsymbol \theta$ and $\dot{\boldsymbol \sigma}_t^2(\boldsymbol \theta)=\partial\sigma_t^2(\boldsymbol \theta)/\partial\boldsymbol \theta$, and define {

} where $\ddot{\boldsymbol \sigma}_t^2(\boldsymbol \theta)=\partial^2\sigma_t^2(\boldsymbol \theta)/\partial\boldsymbol \theta\partial\boldsymbol \theta^{{ \mathrm{\scriptscriptstyle T} }}$ and $\ddot{{\bf g}}_t(\boldsymbol \theta)=\partial^2g_t(\boldsymbol \theta)/\partial\boldsymbol \theta\partial\boldsymbol \theta^{{ \mathrm{\scriptscriptstyle T} }}$. For the reduced model (ref), ${\bf s}_t(\boldsymbol \theta)$ and ${\bf H}_t(\boldsymbol \theta)$ can be simplified as

Let $\widetilde{{\bf s}}_t(\boldsymbol \theta)=\partial\widetilde{\ell}_t(\boldsymbol \theta)/\partial\boldsymbol \theta$ and $\widetilde{{\bf H}}_t(\boldsymbol \theta)=-\partial^2\widetilde{\ell}_t(\boldsymbol \theta)/\partial\boldsymbol \theta\partial\boldsymbol \theta^{{ \mathrm{\scriptscriptstyle T} }}$. To obtain the asymptotic distribution of $\widehat{\boldsymbol \theta}_n$, more assumptions are needed.

The following theorem and corollary state the asymptotic distributions of $\widehat{\boldsymbol \theta}_n$ defined in (ref) for models (ref) and (ref), respectively.

In empirical applications, to make statistical inference on $\boldsymbol \theta_0$, we usually need to estimate the asymptotic covariance matrix ${\bf A}_0^{-1}{\bf B}_0{\bf A}_0^{-1}$. To this end, it suffices to construct weakly consistent estimators of ${\bf A}_0$ and ${\bf B}_0$, respectively. Let

In empirical applications, we are often interested in the statistical significance of some specific component of $\boldsymbol \theta_0$ in modeling, i.e., testing the hypothesis: $H_{0j}:\theta_{0j}=0$ v.s. $H_{1j}:\theta_{0j}\neq 0$ for some $j\in\{1,\dots,d\}.$ Then the Student's $t$-type of test statistic can be constructed in the form $\mathcal{T}_{nj}=\sqrt{n}\big({\bf e}_j^{{ \mathrm{\scriptscriptstyle T} }}\widehat{{\bf A}}_n^{-1}\widehat{{\bf B}}_n\widehat{{\bf A}}_n^{-1}{\bf e}_j\big)^{-1/2}\widehat{\theta}_{nj},$ where the $j$-th entry of ${\bf e}_j$ takes 1 and the rest 0. By Theorems (ref) and (ref), under $H_{0j}$, $\mathcal{T}_{nj}\to_d\mathcal{N}(0,1)$ as $n\to\infty$. Thus, for a given significance level $\alpha\in(0, 1)$, we can reject $H_{0j}$ if $|\mathcal{T}_{nj}|>\Phi^{-1}(1-\alpha/2)$, where $\Phi^{-1}(\cdot)$ is the quantile function of the standard normal distribution.

Testing

Consider the hypothesis

where ${\bf R}$ is a $q\times d$ matrix with full row rank, and ${\bf r}$ is a $q$-dimensional vector.

Based on the LQMLE in Section (ref), we develop the Wald test and the Lagrange multiplier test for the hypothesis (ref). Under $H_0$, by Theorems (ref) and (ref), it follows that

Then, the Wald test statistic is the inner product of the left-hand side of (ref), that is,

For the Lagrange multiplier test, we consider the maximization problem (ref) subject to the constraint ${\bf R}\boldsymbol \theta={\bf r}.$ The constrained LQMLE $\widetilde{\boldsymbol \theta}_n$ of $\boldsymbol \theta_0$ is defined as

The Lagrangian of this constrained problem is defined as $\mathcal{J}_n(\boldsymbol \theta,\boldsymbol \lambda)=\widetilde{\mathcal{L}}_{n}(\boldsymbol \theta)+n\boldsymbol \lambda^{{ \mathrm{\scriptscriptstyle T} }}({\bf R}\boldsymbol \theta-{\bf r}),$ where $\boldsymbol \lambda$ is the vector of Lagrange multipliers. By the Lagrange duality theory, there exists $\widetilde{\boldsymbol \lambda}_n\in\mathbb{R}^q$ such that the solution of (ref) is exactly the solution to the unconstrained problem $\arg\max_{\boldsymbol \theta\in\boldsymbol \Theta}\mathcal{J}_n(\boldsymbol \theta,\widetilde{\boldsymbol \lambda}_n).$ Define $\boldsymbol \Lambda=({\bf R}{\bf A}_0^{-1}{\bf R}^{{ \mathrm{\scriptscriptstyle T} }})^{-1}{\bf R}{\bf A}_0^{-1}{\bf B}_0{\bf A}_0^{-1}{\bf R}^{{ \mathrm{\scriptscriptstyle T} }}({\bf R}{\bf A}_0^{-1}{\bf R}^{{ \mathrm{\scriptscriptstyle T} }})^{-1}.$ It can be shown that $\boldsymbol \Lambda^{-1/2}\sqrt{n}\widetilde{\boldsymbol \lambda}_n\to_d\mathcal{N}(\boldsymbol 0_q,\,{\bf I}_q).$ Once we obtain the constrained LQMLE $\widetilde{\boldsymbol \theta}_n$, a plug-in estimator of $\boldsymbol \Lambda$ follows, namely $\widetilde{\boldsymbol \Lambda}_n=({\bf R}\widetilde{{\bf A}}_n^{-1}{\bf R}^{{ \mathrm{\scriptscriptstyle T} }})^{-1}{\bf R}\widetilde{{\bf A}}_n^{-1}\widetilde{{\bf B}}_n\widetilde{{\bf A}}_n^{-1}{\bf R}^{{ \mathrm{\scriptscriptstyle T} }}({\bf R}\widetilde{{\bf A}}_n^{-1}{\bf R}^{{ \mathrm{\scriptscriptstyle T} }})^{-1},$ where $\widetilde{{\bf A}}_n=n^{-1}\sum_{t=1}^{n}\widetilde{{\bf H}}_t(\widetilde{\boldsymbol \theta}_n)$ and $\widetilde{{\bf B}}_n=n^{-1}\sum_{t=1}^{n}\widetilde{{\bf s}}_t(\widetilde{\boldsymbol \theta}_n)\widetilde{{\bf s}}_t(\widetilde{\boldsymbol \theta}_n)^{{ \mathrm{\scriptscriptstyle T} }}.$ Then, the Lagrange multiplier test statistic can be defined as

The limiting distributions of the Wald test statistic and the Lagrange multiplier test statistic follow from the asymptotics of $\widehat{\boldsymbol \theta}_n$ and $\widetilde{\boldsymbol \theta}_n$ and the continuous mapping theorem.

Applications

In this section, we apply the proposed method and theory to five important time series models in nonlinear time series analysis. We focus on checking relevant technical assumptions used to ensure the asymptotics of the LQMLE.

DAR Model

A time series $\{y_t\}$ is said to follow a double autoregressive (DAR) model of order $(p, q)$ if it satisfies the stochastic recurrence equation:

where $\phi_i\in\mathbb{R}$, $\alpha_0>0$, $\alpha_j\ge0$ for $1\leq j\leq q$, $\phi_p\alpha_q\neq0$, $\{\eta_t\}$ is a sequence of i.i.d. random variables with $\eta_t$ independent of $\mathcal{F}_{t-1}$. For model (ref) with symmetric innovations, guegan1994probabilistic proved it is strictly stationary and ergodic if the top Lyapunov exponent $\gamma<0$. Particularly, for $p=q=1$, it is shown that $\gamma=\mathbb{E}\log|\phi_1+\eta_t\sqrt{\alpha_1}|$. By Theorem 2.1 in cline2004stability, a sufficient condition for strict stationarity and ergodicity of $\{y_t\}$ is that $\sum_{i=1}^{p\vee q}(|\phi_i|^2+\alpha_i\mathbb{E}|\eta_t|)<1$ with $\phi_i=0$ for $i>p$ and $\alpha_i=0$ for $i>q,$ and $\eta_t$ having a continuous and positive density over $\mathbb{R}$.

Let $d=p+q+2$, $m=\max\{p, q\}$, $\boldsymbol \theta=(\boldsymbol \phi^{{ \mathrm{\scriptscriptstyle T} }},\boldsymbol \alpha^{{ \mathrm{\scriptscriptstyle T} }})^{{ \mathrm{\scriptscriptstyle T} }}$ with $\boldsymbol \phi=(\phi_0,\dots,\phi_p)^{{ \mathrm{\scriptscriptstyle T} }}$ and $\boldsymbol \alpha=(\alpha_0,\dots,\alpha_q)^{{ \mathrm{\scriptscriptstyle T} }}$, ${\bf Y}_t=(y_t,\dots,y_{t-m+1})^{{ \mathrm{\scriptscriptstyle T} }},g_t(\boldsymbol \theta)=\phi_0+\sum_{i=1}^{p}\phi_iy_{t-i}$, and $\sigma_t(\boldsymbol \theta)=(\alpha_0+\sum_{j=1}^{q}\alpha_jy_{t-j}^2)^{1/2}.$ Assumption (ref) is guaranteed by $\gamma<0$. Assumption (ref) is satisfied due to the compactness of $\boldsymbol \Theta$. Assumptions (ref), (ref), (ref), (ref), and (ref) hold automatically. Assumption (ref) is satisfied by assuming that $\mathbb{E}|y_t|^{\iota}<\infty$ for some $\iota>0$. In addition, note that $\dot{{\bf g}}_t(\boldsymbol \theta)=(1,y_{t-1},\dots,y_{t-p},\boldsymbol 0_{q+1}^{{ \mathrm{\scriptscriptstyle T} }})^{{ \mathrm{\scriptscriptstyle T} }}$, and $\dot{\boldsymbol \sigma}_t^2(\boldsymbol \theta)=(\boldsymbol 0_{p+1}^{{ \mathrm{\scriptscriptstyle T} }},1,y_{t-1}^2,\dots,y_{t-q}^2)$. Then, Assumptions (ref) and (ref) are met by the compactness of $\boldsymbol \Theta$ and $\mathbb{E} y_t^2<\infty$, which is implied by Assumption (ref) and is not necessary if $q\ge p.$ Note that the dimension of the initial value is finite for model (ref). Thus, together with Remark (ref), the proposed LQMLE for the DAR model (ref) is strongly consistent and asymptotically normal.

GARCH Model

Consider the classical GARCH$(p, q)$ model

where $\alpha_0>0$, $\alpha_i\ge0$, $\beta_j\ge0$ for $i,j\geq 1$, $\alpha_p\beta_q>0$, $\{\eta_t\}$ is a sequence of i.i.d. random variables and $\eta_t$ is independent of $\mathcal{F}_{t-1}$. By bougerol1992stationarity, $\{y_t\}$ is strictly stationary and ergodic if and only if the top Lyapunov exponent $\gamma<0.$ Particularly, when $p=q=1$, it is known that $\gamma=\mathbb{E}\log(\beta_1+\alpha_1\eta_t^2)$.

Let $d=p+q+1$, $m=\infty$, $\boldsymbol \theta=(\boldsymbol \alpha^{{ \mathrm{\scriptscriptstyle T} }},\boldsymbol \beta^{{ \mathrm{\scriptscriptstyle T} }})^{{ \mathrm{\scriptscriptstyle T} }}$ with $\boldsymbol \alpha=(\alpha_0,\dots,\alpha_p)^{{ \mathrm{\scriptscriptstyle T} }}$ and $\boldsymbol \beta=(\beta_1,\dots,\beta_q)^{{ \mathrm{\scriptscriptstyle T} }}$, ${\bf Y}_t=(y_{t},y_{t-1},\dots)^{{ \mathrm{\scriptscriptstyle T} }}\in\mathbb{R}^{\infty}$, $g_t(\boldsymbol \theta)\equiv0,$ and $\sigma_t(\boldsymbol \theta)=\big\{\alpha_0+\sum_{i=1}^{p}\alpha_iy_{t-i}^2+\sum_{j=1}^{q}\beta_j\sigma_{t-j}^2(\boldsymbol \theta)\big\}^{1/2}.$ On assuming that $\alpha_i$'s and $\beta_j$'s are bounded away from zero and $\gamma<0$, Assumptions (ref)--(ref) and (ref)--(ref) can be verified by arguments similar to Section (ref). Since the dimension of the initial value $\widetilde{{\bf Y}}_0$ in optimizing the objective function (ref) for model (ref) is infinite, we focus on checking the initial conditions. For simplicity, we here consider a GARCH(1,1) model as an example, noting that the discussion on higher-order cases is similar. We know that $\beta_1\in(0,1)$ since $\gamma<0.$ Let $\sigma_t^2(\boldsymbol \theta)=\alpha_0+\alpha_1y_{t-1}^2+\beta_1\sigma_{t-1}^2(\boldsymbol \theta)$ with $\boldsymbol \theta=(\alpha_0,\alpha_1,\beta_1)^{{ \mathrm{\scriptscriptstyle T} }}$. By iterations, we have $\sigma_t^2(\boldsymbol \theta)=\alpha_0\sum_{k=1}^{\infty}\beta_1^{k-1}+\alpha_1\sum_{h=1}^{\infty}\beta_1^{h-1}y_{t-h}^2$. Without loss of generality, suppose that the initial value is $\widetilde{{\bf Y}}_0=\bf0$, which is most commonly used in optimizing the objective function. Then $\widetilde{\sigma}_t^2(\boldsymbol \theta)=\alpha_0\sum_{k=1}^{\infty}\beta_1^{k-1}+\alpha_1\sum_{h=1}^{t-1}\beta_1^{h-1}y_{t-h}^2$. Thus $$ \sigma_t^2(\boldsymbol \theta)-\widetilde{\sigma}_t^2(\boldsymbol \theta)=\alpha_1\sum_{h=t}^{\infty}\beta_1^{h-1}y_{t-h}^2=\alpha_1\sum_{h=0}^{\infty}\beta_1^{t+h-1}y_{-h}^2=\beta_1^t\left(\alpha_1\sum_{h=0}^{\infty}\beta_1^{h-1}y_{-h}^2\right). $$ Let $C_1(\boldsymbol \theta)=\alpha_1\sum_{h=0}^{\infty}\beta_1^{h-1}y_{-h}^2$ and $C_1=\sup_{\boldsymbol \theta\in\boldsymbol \Theta}|C_1(\boldsymbol \theta)|=\bar{\alpha}_1\sum_{h=0}^{\infty}\bar{\beta}_1^{h-1}y_{-h}^2$, where $\bar{\alpha}_1=\sup\{\alpha_1|\boldsymbol \theta\in\boldsymbol \Theta\}<\infty$ and $\bar{\beta}_1=\sup\{\beta_1|\boldsymbol \theta\in\boldsymbol \Theta\}<1$ due to the compactness of $\boldsymbol \Theta$ and $\gamma<0$. Clearly, $C_1\in \mathcal{F}_{0}$ is a nonnegative random variable independent of $\boldsymbol \theta$. Thus, Assumption (ref) is satisfied by letting $\rho_1=\bar{\beta}_1$. Then, by simple algebraic calculations, we have $$ \frac{\partial\sigma_t^2(\boldsymbol \theta)}{\partial\boldsymbol \theta}-\frac{\partial\widetilde{\sigma}_t^2(\boldsymbol \theta)}{\partial\boldsymbol \theta}=\Big(0,\,\sum_{h=0}^{\infty}\beta_1^{t+h-1}y_{-h}^2,\,\sum_{h=0}^{\infty}\alpha_1(t+h-1)\beta_1^{t+h-2}y_{-h}^2\Big)^{{ \mathrm{\scriptscriptstyle T} }}, $$ and thus $\big\Vert\partial\sigma_t^2(\boldsymbol \theta)/\partial\boldsymbol \theta-\partial\widetilde{\sigma}_t^2(\boldsymbol \theta)/\partial\boldsymbol \theta\big\Vert\le \beta_1^t\big[\sum_{h=0}^{\infty}\{\beta_1+(t+h-1)\alpha_1\}\beta_1^{h-2}y_{-h}^2\big]$. Similarly, it can be shown that $\big\Vert\partial^2\sigma_t^2(\boldsymbol \theta)/\partial\boldsymbol \theta\partial\boldsymbol \theta^{{ \mathrm{\scriptscriptstyle T} }}-\partial^2\widetilde{\sigma}_t^2(\boldsymbol \theta)/\partial\boldsymbol \theta\partial\boldsymbol \theta^{{ \mathrm{\scriptscriptstyle T} }}\big\Vert\le \beta_1^t\big[\sum_{h=0}^{\infty}\{2(t+h-1)+(t+h-1)(t+h-2)\alpha_1/\beta_1\}\beta_1^{h-2}y_{-h}^2\big]$. Note that there exists a constant $\rho_2$ with $0<\rho_2<1$ such that $t^2\bar{\beta}_1^t<\rho_2^t$ for $t$ large enough. Then Assumption (ref) is satisfied by adopting the sufficient condition discussed in Section {\color{blue}S.3.1} of the supplementary material. Thus, the proposed LQMLE for the GARCH model (ref) is strongly consistent and asymptotically normal.

ARMA-GARCH Model

Consider an ARMA$(p,q)$-GARCH$(r,s)$ model

where $\phi_i\in\mathbb{R}$, $\varphi_j\in\mathbb{R}$, $\phi_p\varphi_q\neq0$, $\alpha_0>0$, $\alpha_k\ge0$, $\beta_l\ge0$, $\alpha_r\beta_s>0$, $\{\eta_t\}$ is i.i.d. and $\eta_t$ is independent of $\mathcal{F}_{t-1}$. The strict stationarity and ergodicity of the ARMA-GARCH model were studied in ling1997fractionally, ling2003asymptotic, etc.

For illustration, model (ref) can be written in the form of model (ref). To this end, we let $\phi(z)=1-\phi_1z-\dots-\phi_pz^p$ and $\varphi(z)=1+\varphi_1z+\dots+\varphi_qz^q$ be the $p$-th and $q$-th degree characteristic polynomials for any $z\in\mathbb{C}$, respectively. Suppose that (i) $\phi(z)\neq 0$, $\varphi(z)\neq 0$ for all $|z|\leq 1$, and $\phi(z)$ and $\varphi(z)$ have no common factors; and (ii) $\sum_{k=1}^{r}\alpha_k+\sum_{l=1}^{s}\beta_l<1$. Then, under (ii), $\{\sigma_t\}$ is stationary (strictly and weakly) and ergodic. Further, $\{y_t\}$ is (causal) stationary, ergodic, and invertible under (i). Thus, by iterations, we can see that $ \varepsilon_t=y_t-\phi_0-\sum_{i=1}^{p}\phi_iy_{t-i}-\sum_{j=1}^{q}\varphi_j\varepsilon_{t-j}\in\mathcal{F}_{t}$. Then the conditional mean function $\phi_0+\sum_{i=1}^{p}\phi_iy_{t-i}+\sum_{j=1}^{q}\varphi_j\varepsilon_{t-j}\in\mathcal{F}_{t-1}$, i.e., there exists a measurable function $g_t(\boldsymbol \theta)\in\mathcal{F}_{t-1}$ such that $g_t(\boldsymbol \theta)=\phi_0+\sum_{i=1}^{p}\phi_iy_{t-i}+\sum_{j=1}^{q}\varphi_j\varepsilon_{t-j}$. Similarly, $\sigma_t^2\in\mathcal{F}_{t-1}$ by noting that $\sigma_t^2$ is a measurable function of $\{\varepsilon^2_{t-j}: j\geq 1\}$. Due to $\sigma_t>0$ a.s., we have $\sigma_t\in\mathcal{F}_{t-1}$. Thus, there exists $\sigma_t(\boldsymbol \theta)\in\mathcal{F}_{t-1}$ such that $\sigma_t(\boldsymbol \theta)=\big(\alpha_0+\sum_{k=1}^{r}\alpha_k\varepsilon_{t-k}^2+\sum_{l=1}^{s}\beta_l\sigma_{t-l}^2\big)^{1/2}$. Therefore, model (ref) has the form of model (ref).

For simplicity, we here consider an ARMA(1,1)-GARCH(1,1) model as a benchmark. Specifically, $y_t=\phi_0+\phi_1y_{t-1}+\varphi_1\varepsilon_{t-1}+\varepsilon_t$ with $\varepsilon_t=\sigma_t\eta_t$ and $\sigma_t^2=\alpha_0+\alpha_1\varepsilon_{t-1}^2+\beta_1\sigma_{t-1}^2.$ Let $\boldsymbol \theta=(\phi_0,\phi_1,\varphi_1,\alpha_0,\alpha_1,\beta_1)^{{ \mathrm{\scriptscriptstyle T} }}.$ Suppose that (i) $^{\prime}$\, $|\phi_1|<1$, $|\varphi_1|<1$, and $\phi_1+\varphi_1\neq0$; (ii)$^{\prime}$\, $\alpha_1+\beta_1<1$. By Lemma 2.2 in francq2019garch and the fact $\mathbb{E}|\eta_t|<\infty$ implied by Assumption (ref), there exists some constant $\iota>0$ such that $\mathbb{E}\sup_{\boldsymbol \theta\in\boldsymbol \Theta}\sigma_t^{\iota}(\boldsymbol \theta)<\infty.$ By iterations, it follows that

Let $g_t(\boldsymbol \theta)=\phi_0\sum_{k=1}^{\infty}(-\varphi_1)^{k-1}+(\phi_1+\varphi_1)\sum_{h=1}^{\infty}(-\varphi_1)^{h-1}y_{t-h}\in\mathcal{F}_{t-1}$ and $\varepsilon_t(\boldsymbol \theta)=y_t-g_t(\boldsymbol \theta)\in\mathcal{F}_{t}$. Similar to Section (ref), we can get that

Next, we focus on verifying Assumptions (ref), (ref), (ref), and (ref), while the others are straightforward to check. By tedious algebraic calculations, we have { $$

$$} Then, Assumption~\ref{ass.dot_g} is satisfied by noting that the conditions (i)$^{\prime}$--(ii)$^{\prime}$, the compactness of $\boldsymbol \Theta$, and the fact $\mathbb{E}|y_t|<\infty$ implied by Assumption~\ref{ass.eta}. Similarly, Assumption~\ref{ass.Ht} is also satisfied by noting that $\mathbb{E} y_t^2<\infty$ implied by Assumption~\ref{ass.eta_2} and the conditions (i)$^{\prime}$--(ii)$^{\prime}$.

Without loss of generality, suppose that the initial value is $\widetilde{{\bf Y}}_0=\boldsymbol 0.$ Then let $\widetilde{g}_t(\boldsymbol \theta)=g(\widetilde{{\bf Y}}_{t-1},\boldsymbol \theta),$ $\widetilde{\sigma}_t(\boldsymbol \theta)=\sigma(\widetilde{{\bf Y}}_{t-1},\boldsymbol \theta)$ and $\widetilde{\varepsilon}_t(\boldsymbol \theta)=\widetilde{y}_t-\widetilde{g}_t(\boldsymbol \theta)$. It can be shown that { $$

$$} which implies Assumption~\ref{ass.init_sigma} by similar arguments in Section~\ref{subsec.GARCH}. By tedious algebraic calculations, it follows that {\small $$

$$} Similarly, we can also obtain that the expressions of $\partial^2 g_t(\boldsymbol \theta)/\partial\boldsymbol \theta\partial\boldsymbol \theta^{{ \mathrm{\scriptscriptstyle T} }}-\partial^2 \widetilde{g}_t(\boldsymbol \theta)/\partial\boldsymbol \theta\partial\boldsymbol \theta^{{ \mathrm{\scriptscriptstyle T} }}$ and $\partial^2 \sigma_t^2(\boldsymbol \theta)/\partial\boldsymbol \theta\partial\boldsymbol \theta^{{ \mathrm{\scriptscriptstyle T} }}-\partial^2 \widetilde{\sigma}_t^2(\boldsymbol \theta)/\partial\boldsymbol \theta\partial\boldsymbol \theta^{{ \mathrm{\scriptscriptstyle T} }}$. Then Assumption (ref) is satisfied by the sufficient condition discussed in Section {\color{blue}S.3.1} of the supplementary material. Thus, the proposed LQMLE for the ARMA-GARCH model (ref) is strongly consistent and asymptotically normal.

DTARMACH Model

Consider a double threshold ARMA conditional heteroscedasticity (DTARMACH) model, proposed by ling1999probabilistic and defined as

where $\phi_i^{(u)}\in\mathbb{R}$, $\varphi_j^{(u)}\in\mathbb{R}$, $\alpha_0^{(v)}>0$, $\alpha_k^{(v)}\ge0$, $\beta_l^{(v)}\ge0$ for $k,l>0$ and $u=1,\dots,U$, $v=1,\dots,V,$ the positive integers $b,h$ are the delay parameters, the threshold parameters satisfy $-\infty=a_0<a_1<\cdots<a_U=\infty$, and $-\infty=c_0<c_1<\cdots<c_V=\infty,\{\eta_t\}$ is a sequence of i.i.d. random variables and $\eta_t$ is independent of $\mathcal{F}_{t-1}$. The DTARMACH model includes many well-known models such as GARCH, TARMA brockwell1992existence, and DTARCH li1996double. The strict stationarity and ergodicity of the DTARMACH model were studied in ling1999probabilistic and a sufficient condition is

The DTARMACH model is a piecewise ARMA-GARCH model and can be regarded as an extension of model (ref). To incorporate model (ref) into our framework, we assume that the threshold parameters are all known, the sufficient condition ((ref)) holds, and $\sum_{j=1}^{q}\max_{u}|\varphi_j^{(u)}|<1$. Then $\{y_t\}$ is invertible and $\varepsilon_t$ can be written as the infinite summation of $\{y_s:s\le t\}$. For simplicity, let $p=q=r=s=1$, $d=3U+3V,$ and $\boldsymbol \theta=\big(\phi_0^{(1)},\phi_1^{(1)},\varphi^{(1)},\alpha_0^{(1)},\alpha_1^{(1)},\beta^{(1)},\dots,\phi_0^{(U)},\phi_1^{(U)},\varphi^{(U)},\alpha_0^{(V)},\alpha_1^{(V)},\beta^{(V)}\big)^{{ \mathrm{\scriptscriptstyle T} }}\in\mathbb{R}^{d}.$ Denote $I^{(u)}(y_{t-b})=I(a_{u-1}<y_{t-b}\le a_u)$ and $I^{(v)}(y_{t-h})=I(c_{v-1}<y_{t-h}\le c_v)$. By iterations, similar to (ref), it follows that {

} where the convention $\prod_{l=1}^{0}\cdot\equiv1$ is adopted. Further, let {

} and $\varepsilon_t(\boldsymbol \theta)=y_t-g_t(\boldsymbol \theta)\in \mathcal{F}_{t}$. Similarly, we can get {

} Thus, model (ref) has the form of model (ref). By similar arguments to Section (ref), we can verify Assumptions (ref), (ref), (ref), and (ref) hold, while the others are straightforward.

EXPAR Model

Consider an exponential autoregressive (EXPAR) model of order $p$:

where $\phi_i\in\mathbb{R}$, $\varphi_i\in\mathbb{R}$ for $1\leq i\leq p$, $|\phi_p|+|\varphi_p|>0$, $\delta>0$, $\{\eta_t\}$ is a sequence of i.i.d. random variables and $\eta_t$ is independent of $\mathcal{F}_{t-1}$. The EXPAR model Jones1978,haggan1981modelling is a nonlinear time series model partly motivated by amplitude-frequency dependency. By Theorem 2.1 in chen2018generalized, a sufficient condition for strict stationarity and ergodicity of $\{y_t\}$ in (ref) follows that all the roots of the characteristic equation: $z^p-(|\phi_1|+|\varphi_1|)z^{p-1}-\dots-(|\phi_p|+|\varphi_p|)=0$, are inside the unit circle, and $\eta_t$ has a continuous and positive density over $\mathbb{R}$.

Let $d=2p+1$, $m=p$, $\boldsymbol \theta=(\boldsymbol \phi^{{ \mathrm{\scriptscriptstyle T} }},\boldsymbol \varphi^{{ \mathrm{\scriptscriptstyle T} }},\delta)^{{ \mathrm{\scriptscriptstyle T} }}$ with $\boldsymbol \phi=(\phi_1,\dots,\phi_p)^{{ \mathrm{\scriptscriptstyle T} }}$, $\boldsymbol \varphi=(\varphi_1,\dots,\varphi_p)^{{ \mathrm{\scriptscriptstyle T} }}$, ${\bf Y}_t=(y_t,\dots,y_{t-p+1})^{{ \mathrm{\scriptscriptstyle T} }}$, $g({\bf Y}_{t-1},\boldsymbol \theta)=\sum_{i=1}^{p}\{\phi_i+\varphi_i\exp(-\delta y_{t-1}^2)\}y_{t-i}$, and $\sigma({\bf Y}_{t-1},\boldsymbol \theta)=1.$ Assumptions (ref)--(ref), (ref), (ref), and (ref) are satisfied automatically, and Assumption (ref) is satisfied by assuming that $\mathbb{E}|y_t|^{\iota}<\infty$ for some $\iota>0$. Then note that $$ \dot{{\bf g}}_t(\boldsymbol \theta)=\big({\bf Y}_{t-1}^{{ \mathrm{\scriptscriptstyle T} }}, ~\exp(-\delta y_{t-1}^2){\bf Y}_{t-1}^{{ \mathrm{\scriptscriptstyle T} }},~ -y_{t-1}^2\exp(-\delta y_{t-1}^2)\boldsymbol \varphi^{{ \mathrm{\scriptscriptstyle T} }}{\bf Y}_{t-1}\big)^{{ \mathrm{\scriptscriptstyle T} }}, $$ and $\ddot{{\bf g}}_t(\boldsymbol \theta)$ can be accordingly obtained. To guarantee Assumptions (ref) and (ref), it only requires the conditions $\mathbb{E}\Vert\dot{{\bf g}}_t(\boldsymbol \theta_0)\Vert<\infty$ and $\mathbb{E}\Vert\ddot{{\bf g}}_t(\boldsymbol \theta_0)\Vert^2<\infty$, which can be verified by using $\mathbb{E} |y_t|<\infty$ and $\mathbb{E} y_t^2<\infty$, which are implied by Assumption (ref) and (ref), respectively. Thus, the LQMLE of the EXPAR model is strongly consistent and asymptotically normal.

Simulation Studies

Performance of the LQMLE

To assess the finite-sample performance of the LQMLE, we consider an ARMA(1,1)-GARCH(1,1) model and a DAR(1,1) model for illustrations in the following examples.

In each simulation scenario, we use the length of observations $n=100$, 200, and 400 with $1000$ replications. The absolute estimation bias and standard deviation results for each simulation experiment of models (ref) and (ref) are reported in Tables (ref) and (ref), respectively. From the tables, we can find that (i) with a larger $n$, both the absolute estimation bias and standard deviation of each parameter show decreasing trends; (ii) even though it is motivated by the log-likelihood function of a standard logistic distribution when the model innovations come from other distributions, our LQMLE is still consistent, which implies that the LQMLE is robust to distributional misspecifications.

Figures (ref)--(ref) plot the histograms of $\sqrt{n}(\widehat{\boldsymbol \theta}_n-\boldsymbol \theta_0)$ for models (ref) and (ref) with $n=400$, the innovations $\{\eta_t\}$ being generated from i.i.d. $\mathcal{N}(0,1.75^2)$ and $1.25 t_3$, and the true parameter $(\phi_1,\varphi_1,\alpha_0,\alpha_1,\beta_1)=(0.3,0.2,0.2,0.1,0.3)$ and $(\phi_0,\phi_1,\alpha_0,\alpha_1)=(1.0,0.5,0.3,0.5)$, respectively. The asymptotic standard deviations are calculated from the asymptotic covariance matrix in Theorems (ref) through Monte Carlo simulation. It can be seen that the empirical densities of each estimator closely align with the asymptotic normal distributions, despite the absence of finite fourth moments of the innovations, which are generated from a Student's $t_{3}$-distribution. These findings strongly support Theorems (ref) and demonstrate the robustness of the LQMLE compared to the GQMLE.

Finally, we compare the finite-sample performance of our LQMLE and the GQMLE. As discussed in Remark (ref), the identifiability conditions for the LQMLE and GQMLE are different so that it is infeasible to fairly compare the finite-sample performance of the LQMLE of parameters in volatility function directly. Thus, we first focus on the estimators of parameters in mean functions, i.e., $(\phi_1,\varphi_1)$ in ARMA-GARCH model (ref), and $(\phi_0,\phi_1)$ in DAR model (ref). In each simulation scenario, we use $n=100,200,400$ with $1000$ replications. The absolute estimation bias and standard deviation results for two estimators are reported in Table (ref).

From the table, we can see that (i) when $\eta_t$ follows the logistic distribution, the LQMLE is indeed the MLE and its outperformance is clear; (ii) when $\eta_t$ follows the normal distribution, the GQMLE reduces to the MLE and outperforms the LQMLE, which is unquestionable; (iii) when $\eta_t$ follows the uniform distribution, the GQMLE outperforms the LQMLE since the uniform distribution is lighted-tailed; (iv) when $\eta_t$ follows the Student's $t_\nu$ and stable distributions, the LQMLE achieves significantly better performance than the GQMLE, which shows that the LQMLE is robust to heavy-tailed innovations. Moreover, the heavier the tail of the innovation, the better the performance of the LQMLE.

To fairly compare the finite-sample performance of our LQMLE and the GQMLE of the parameters in the volatility function, a rescaling technique li2018zd can be used to address the issues caused by different identifiability conditions. The related numerical results are presented in Section {\color{blue}S.4.1} of the supplementary material, where we have similar findings.

Performance of the Test Statistics

To examine the finite-sample performance of two tests proposed in Section (ref), we consider the null hypothesis $H_0:\boldsymbol \theta_0=\boldsymbol \theta_{H_0}=(\phi_{1,H_0},\varphi_{1,H_0},\alpha_{0,H_0},\alpha_{1,H_0},\beta_{1,H_0})^{{ \mathrm{\scriptscriptstyle T} }}=(0.3,0.2,0.2,0.1,0.3)^{{ \mathrm{\scriptscriptstyle T} }}$ and ${\bf R}=(1,1,2,3,1)$ for model (ref) , and $H_0:\boldsymbol \theta_0=\boldsymbol \theta_{H_0}=(\phi_{0,H_0},\phi_{1,H_0},\alpha_{0,H_0},\alpha_{1,H_0})^{{ \mathrm{\scriptscriptstyle T} }}=(1.0,0.5,0.3,0.5)^{{ \mathrm{\scriptscriptstyle T} }}$ and ${\bf R}=(1,1,1,1)$ for model (ref). As discussed in Remark (ref), the asymptotic normality of the estimator requires the innovations with finite variance, thus $\eta_t$ is generated by one of the first four random variables discussed in Example (ref). The nominal significance level is $\alpha=5\%.$ To present the empirical powers, we consider three alternatives $H_1:\boldsymbol \theta_0=1.1\boldsymbol \theta_{H_0}$, $H_1:\boldsymbol \theta_0=1.3\boldsymbol \theta_{H_0}$, and $H_1:\boldsymbol \theta_0=1.5\boldsymbol \theta_{H_0}$ for both two models, respectively. In each simulation scenario, we use $n=100,200,400$ with $1000$ replications. The empirical sizes and powers of the hypothesis tests for models (ref) and (ref) are reported in Tables (ref) and (ref), respectively.

From the tables, we can see that (i) for $\eta_t$ generated by four different random variables, the empirical sizes of the Wald test and the Lagrange multiplier test are close to the nominal significance level $\alpha=5\%$ in most cases when $n$ is large; (ii) with a larger $n$, the empirical powers of the two hypothesis tests show increasing trends. For example, under the ARMA-GARCH model (ref), for $H_1:\boldsymbol \theta_0=1.3\boldsymbol \theta_{H_0}$, $\eta_t$ is generated by $\mathcal{N}(0,1.75^2)$, when $n$ increases from 100 to 400, the empirical powers of $\mathcal{T}_n^{{ \mathrm{\scriptscriptstyle W} }}$ and $\mathcal{T}_n^{{ \mathrm{\scriptscriptstyle LM} }}$ increase quickly from 0.491 to 0.975, and from 0.463 to 0.993, respectively. Similar simulation studies on a GARCH(1,1) model will be presented in Section {\color{blue}S.4.2} of the supplementary material. In summary, our numerical results provide strong support to our theoretical results in Sections (ref) and (ref).

Real Data Analysis

To showcase our method, we consider the U.S. monthly 3-month treasury par yield curve rate series from January 1990 to December 2024 with a total of 420 observations, which are available at \url{https://home.treasury.gov/}. Denote the data by $x_t$ and $y_t=x_t-x_{t-1}$ as the first difference of $x_t$. The Phillips-Perron test phillips1988testing suggests that $x_t$ is not stationary, while $y_t$ can be viewed as stationary. Figure (ref) plots the monthly series $\{x_t\}_{t=0}^{419}$ and $\{y_t\}_{t=1}^{419}$, respectively. li2023maximum also studied this dataset (until September 2021) and used an $\alpha$-stable DAR to fit the series $\{y_t\}.$ Their proposed test indicates that the innovation follows a symmetric $\alpha$-stable distribution with $1<\alpha<2$ which corresponds to the heavy-tailed distribution without finite second-order moment.

Then we fit three different models discussed in Section (ref), i.e., DAR(1,1), GARCH(1,1), and ARMA(1,1)-GARCH(1,1) models, respectively, to the series $\{y_t\}_{t=1}^{419}.$ The fitting results and related statistics of all three models are summarized in Table (ref), which provides the following information. First, ARMA(1,1)-GARCH(1,1) model achieves the largest log-likelihood (226.778) and the smallest AIC ($-441.555$) among these three models, and the parameters $\phi_1,\varphi_1,\alpha_0,\alpha_1$ and $\beta_1$ are significant at the significance level $\alpha=5\%$. Second, the Lyapunov exponents of these three models are $-0.6524$, $-0.7863$, and $-0.6437$, respectively, which are all less than 0 and indicate that the series is strictly stationary and ergodic. Third, the Hill estimators of the residuals $\{\widehat{\eta}_t\}$ in these three models are 0.8163, 0.6716, and 0.7377, respectively, which are all less than 1 and indicate that $\{\widehat{\eta}_t\}$ tends to be heavy-tailed. The histograms of the residuals $\{\widehat{\eta}_t\}$ in these three models are plotted in Figure (ref). Thus, it is more appropriate to model the given using our LQMLE.

Conclusion and Discussion

In this article, we have proposed a novel logistic quasi-maximum likelihood methodology in the context of time series analysis. It enjoys numerous advantages. Specifically, it is robust in respect of distributional misspecification and heavy-tailedness of the innovation, and is more resilient to outliers than the Gaussian quasi-maximum likelihood method and the least squares method. In our asymptotic theory, conditional symmetry of the innovation is assumed, which is mild in practice. When there exist asymmetry phenomena in the data, we can introduce asymmetric structures or threshold effects in time series models to circumvent asymmetry issues.

In this article, we only focus on univariate parametric time series models. There are several topics worthy of further study. For example, with our methodology, we can consider spatial autoregressive models, threshold autoregressive models with unknown threshold parameters, multivariate autoregressive models involving a joint multivariate logistic distribution as in malik1973multivariate. We leave these topics for future research.