Robust estimation for Threshold Autoregressive Moving-Average models

Robust estimation for Threshold Autoregressive Moving-Average models

\affil[1]{Faculty of Economics and Management, Free University of Bozen-Bolzano, Italy} \affil[2]{Department of Statistical Sciences, University of Bologna, Italy} \affil[3]{Department of Data Science and Analytics, BI Norwegian Business School, Norway}

{\it Keywords:} Threshold Autoregressive Moving-Average models; Non-linear time series; Robust estimation; Outliers; Commodity prices

Introduction

Threshold models are popular tools used to describe complex phenomena in many fields, including economics, finance, ecology, epidemiology Ton90,Cha09, giordani2007unified, Ton11,Han11, Cha17b. Non-linearity is introduced by a thresholding mechanism which implies multiple linear regimes; this enables the description of complex non-linear dynamical features, such as jumps, limit cycles, time irreversibility, while retaining good interpretability. Since their introduction by Ton78, threshold models have been widely studied, especially in their autoregressive specification, the so-called threshold autoregressive (TAR) models. In analogy with autoregressive moving-average (MA) models, threshold autoregressive moving-average (TARMA) models extend TAR models by including moving-average components in each regime Ton17.

Although technically more challenging than TAR models due to their non-Markovian nature, TARMA models provide a powerful yet simple framework for many research problems involving non-linear phenomena; e.g., see Ton17, Gor20 and Gor21. Nonetheless, their theoretical development has halted for many years and only recently Cha19 solved the long-standing open problem regarding the probabilistic structure of the first order TARMA model. TARMA models possess a number of desirable features, including the following: they include moving-average (MA) components within a parametric non-linear setting; they naturally account for measurement errors; they are able to describe a wide range of long-run probabilistic behaviors spanning from transience to ergodicity, and even geometric ergodicity. Also, the threshold framework provides a natural way to describe series that appear to behave like random walks, when this behavior is incompatible with the theory underlying the data generating process. One example in economics is the well-known purchasing power parity puzzle, which has motivated the development of unit-root tests where the alternative hypothesis is a stationary threshold model with a local unit-root regime End98,Bec04,Kap06,Bec08,Cha20b. Tests for TARMA nonlinearity have been developed in Li11b and Gor23. Ang22 extend the results in Gor23 and develop a test for non-linear effects in the conditional mean for series with conditional heteroscedasticity by incorporating a GARCH specification.

Compared to the simpler TAR models, estimation for the TARMA models is challenging due to the lack of a linear parameterization conditional on the threshold parameter. Current inference methods mostly rely on the least squares (LS) approach Li11a, which is knowingly influenced by outliers and heavy tails. Although aberrant observations are ubiquitous and appear in many real applications giordani2007unified, the issue of robust estimation and outlier detection for TARMA models has yet to be addressed from either theoretical or methodological viewpoints. On the other hand, robust estimation for the linear ARMA model has been extensively studied maronna2019robust. For the special case of the TAR model, chan1994robust study the effect of additive outliers in LS estimation and propose a generalized M-estimation to mitigate the severe bias of the estimates. zhang2009note consider a general class of robust estimators for threshold autoregressive models and show consistency under regularity conditions. grossi2015robust compare the relative efficiency of M-estimators for the TAR models to that of the LS estimator, showing that the former perform well when the error follows heavy tailed or non-Gaussian distributions. van1999smooth derive a robust estimation method for the parameters in smooth threshold TAR models, using generalized maximum likelihood estimation. Motivated by such a gap in the literature, we consider robust inference for the TARMA model using an M-estimation approach. Our approach consists in replacing the residual sum of squares criterion of Li11a by a function with bounded derivative. This is a crucial feature which is necessary to gain stability of the estimates in the presence of different types of outliers. The resulting estimator for the autoregressive and moving-average parameters is shown to be strongly consistent and asymptotically normal, under standard regularity conditions. To the best of our knowledge, we provide the first results for robust estimation in parametric non-linear time series models with moving-average components. Similarly to the least squares estimator, the threshold parameter is found to be super-consistent with convergence rate of $n^{-1}$, while the autoregressive and moving-average parameters are root-$n$ consistent and asymptotically normal. The methodology is implemented using the special family of objective functions considered in Fer12 and la2015robust, which include the LS estimator as a special case. While common contamination types are shown to increase, sometimes dramatically, bias and variance of the least squares estimator, our estimator mitigates the effect of observations that are incompatible with the assumed model, thus reducing the overall mean squared error. We showcase the performance of our new methodology by analyzing a set of commodity price time series. Commodities are important in economics and finance due to their ability to anticipate the behavior of other macroeconomic variables; e.g., see Hamilton2011 and RR2013. The predictive relationship between commodities and a number of macroeconomic variables is often non-linear, with asymmetric behavior depending on whether prices increase or decrease, e.g. see KV2011 and KV2013. Although TARMA models appear suitable for such data, heavy tails, outliers and non-Gaussian innovations make the least squares estimate untrustworthy. The robust TARMA estimates generally provide a better fit with smaller standard errors compared to non-robust estimates. Our robust TARMA specification confirms the existence of two dynamical regimes, separated by the threshold invariably located at zero, and corresponding to a slow, persistent growth (upper regime) and fast contractions (lower regime). The superior predictive performance of the robust TARMA approach can result in a key advancement in modelling the commodity market. The remainder of the paper is organized as follows. In Section (ref), we describe the general M-estimation approach. In Section (ref), we study the asymptotic behavior and robustness properties of the new estimator. In Section (ref), we study the finite-sample behavior of the new estimator and compare its robustness to the standard LS approach under common contamination models. In Section (ref), we apply the new method for robust estimation and outlier detection for several commodity price series. Conclusions and possible extensions of this work are presented in Section (ref). Further results from the analysis of commodity time series and technical proofs are reported in the Supplementary Material.

Methodology

Model setup and notation

Let $\{X_t \}_{t \in \mathds{Z}}$ be the TARMA process defined by the difference equation

where: $p \in \mathds{N}$ and $q\in \mathds{N}$ are, respectively, the autoregressive and moving-average orders; the $\phi$'s and $\theta$'s are the autoregressive and moving-average parameters, respectively; $1\leq d\leq D_0$, where $D_0\in \mathcal{D} \subset \mathds{N}$ is the delay parameter; $r\in \mathcal{R} \subseteq \mathds{R}$ is the threshold parameter; and $\{\varepsilon_{t}\}$ is the innovation process, with $E(\varepsilon_t)=0$ and $E(\varepsilon^2_t)=\sigma^2 <\infty$, which is usually assumed to be Gaussian white noise. Note that the TARMA model reduces to a linear ARMA as $|r|\rightarrow \infty$.

Equation ((ref)) defines two regimes, which will be referred to as lower and upper regimes corresponding to $X_{r-d}\le r$ and $X_{r-d}> r$, respectively. For each regime we have specific parameter vectors defined as $ \boldsymbol{\phi}_1=(\phi_{1,0},\dots,\phi_{1,p})^\intercal$, $\boldsymbol\theta_1=(\theta_{1,1},\ldots,\theta_{1,q})^\intercal$, $\boldsymbol{\phi}_2=(\phi_{2,0},\dots,\phi_{2,p})^\intercal$, $\boldsymbol\theta_2=(\theta_{2,1},\ldots,\theta_{2,q})^\intercal$, while $\boldsymbol{\phi}=(\boldsymbol{\phi}_1^\intercal,\boldsymbol{\phi}_2^\intercal)^\intercal$ and $\boldsymbol{\theta}=(\boldsymbol{\theta}_1^\intercal,\boldsymbol{\theta}_2^\intercal)^\intercal$ are used to denote autoregressive and moving-average parameters. The vector collecting all autoregressive and moving-average parameters is denoted by $\boldsymbol{\lambda}=(\boldsymbol{\phi}^\intercal_1,\boldsymbol{\phi}^\intercal_2, \boldsymbol{\theta}^\intercal_1,\boldsymbol{\theta}^\intercal_2)^\intercal\in\mathcal{Q}\subseteq\mathds{R}^{2(1+p+q)}$, while the overall parameter vector including also the threshold parameter is denoted by $\boldsymbol{\eta}=(\boldsymbol{\lambda}^\intercal,r,d)\in\mathcal{Q}\times\mathcal{R} \times \mathcal{D} :=\mathcal{H}$. We assume the parameter space $\mathcal{H}$ to be compact and equipped with product metric. For simplicity of exposition, in this work we focus on TARMA models with a full model structure containing all lags up to order $p$ and $q$ for both regimes. However, our methodology can be applied without loss of generality to more complex order structures, including missing lags and specific orders for the upper and lower regimes.

In many real-world applications, the process $X_t$ does not follow exactly the model specified in Equation ((ref)). Although the majority of the observations may be compatible with such model assumptions, real data may diverge substantially from the assumed process due to the presence of heavy-tailed or asymmetric errors and aberrant observations. There are several models that may be used to represent the contamination process in the time series context. One common model is the additive outlier (AO) model, which defines the contaminated process $X_t^\epsilon$ according to $ X_t^\epsilon = X_t + Z^{\epsilon}_{t} W_t$, where $X_t$ is the TARMA process defined in ((ref)); $W_t$ is the contaminating process, independent of $X_t$; and $Z^\epsilon_t$ is a binary process where $P(Z_t^{\epsilon} = 1) = \epsilon$ such that $\epsilon$ is the contamination level. Another common model is the replacement outlier (RO) model, where $X_t^\epsilon = (1-Z^{\epsilon}_{t})X_t + Z^{\epsilon}_{t} W_t$, with $W_t$, $X_t$ and $Z^\epsilon_t$ defined above. Finally, in the innovation outlier (IO) model the outliers affect not only the current observation, but also subsequent observations. IOs are obtained when the innovation $\varepsilon_t$ follows a process different from the assumed nominal model. For example, $\varepsilon_t$ is assumed to follow a Gaussian white noise process while the actual innovation process has the normal mixture distribution $ (1-\epsilon)N(0, \sigma_0^2) + \epsilon N(0, \sigma_1^2) $, with $\sigma_0^2 \ll \sigma_1^2$. Outliers may also differ in their temporal structure. For example, patchy outliers arise from the AO and RO models by letting $Z_t$ be a Markov process remaining in one state for multiple time periods of fixed or random duration.

Robust estimation

For the time series $\{X_1,\dots,X_n\}$, define the residual function $\varepsilon_t(\boldsymbol{\eta}) = X_t-E_{\boldsymbol{\eta}}[X_t|\mathcal{F}_{t-1}]$, $t=1,\dots, n$, where $\mathcal{F}_{t}$ is the sigma algebra generated by $\{X_{t}, X_{t-1},\dots\}$. $E_{\boldsymbol{\eta}}[X_t|\mathcal{F}_{t-1}]$ denotes the expectation of $X_t$ conditional on the process history up to time $t-1$ and is computed with respect to the TARMA model described ((ref)) with parameter $\boldsymbol{\eta}$. From ((ref)), we have

An M-estimate $\hat \boldsymbol{\eta}_n$ of the parameter vector $\boldsymbol{\eta}$ is found by minimizing the objective function

where $\rho:\mathds{R} \mapsto \mathds{R}$ is a loss function often referred to as $\rho$-function in the literature of robust statistics, and $\hat \sigma$ is a robust estimate of scale which is obtained simultaneously with $\boldsymbol{\eta}$ as an $M$-scale estimate. To obtain robustness of $\hat \boldsymbol{\eta}_n$, we require the following standard conditions on $\rho$: (i) $\rho(z)$ is non-decreasing function of $|z|$; (ii) $\rho(0)=0$; (iii) $\rho(z)$ is increasing for $z>0$ such that $\rho(z)<\rho(\infty)$; and (iv) the derivative $\psi(z) = \partial \rho(z)/\partial z$ satisfies $|\psi(z)|<c$ for some finite constant $c>0$ and all $z \in \mathds{R}$. There is a number of functions satisfying the above requirements.

Here we study the $\rho$ function considered in Fer12 by taking $\rho(z) = -(f(z)^\alpha - 1)/\alpha$ for $\alpha>0$, and $\rho(z)= - \log(f(z))$ for $\alpha=0$, where $f$ is the assumed probability density function for the innovations. For the special case of Gaussian innovations, the objective function can be written as

for $\alpha>0$. The limit case $\alpha \rightarrow 0$ corresponds to the maximum likelihood objective with

When $\sigma^2$ is taken as known, minimizing ((ref)) is equivalent to minimizing the residual sum of squares $\sum_{t=1}^{n} \varepsilon_t^2(\boldsymbol{\eta})$. In this respect, the function of Equation ((ref)) represents a robust generalization of the well-established LS estimator for the TARMA model of Li11a. The tuning parameter $\alpha$ controls the trade-off between efficiency and robustness of the underlying estimator; this makes this example particularly useful for analyzing the properties of the estimator for various degrees of robustness. For $\alpha>0$, the derivative $ \psi(z) = \partial \rho(z)/\partial z = f^\alpha(z) \partial \log f(z)/\partial z $ is bounded for common family of density functions and $\psi(z) \rightarrow 0$ as $|z| \rightarrow \infty$, and corresponds to a re-descending M-estimator. On the other hand, for the limit case $\alpha \rightarrow 0$, we have $\rho(z) \rightarrow \log(f(z))$ and $\psi(z) = \partial \log f(z)/\partial z$. This case corresponds to the maximum likelihood estimator and the derivative $\psi(z)$ is typically unbounded, which leads to estimators that are sensitive to the presence of outliers. One practical hurdle in the derivation of $\hat \boldsymbol{\eta}_n$ is the discontinuity of ${\rho}_n(\boldsymbol{\eta})$ in $r$. To cope with this issue, the minimization is carried out in two steps. First, given $r$ and $d$, we take the profile estimator of $\boldsymbol{\lambda}$

and define ${\rho}^*_n(r)={\rho}_n(\hat{\boldsymbol{\lambda}}_n(r),r, d)$. Second, since ${\rho}^*_n(r, d)$ can only take a finite number of values, it can be minimized by searching over some grid $\widetilde{\mathcal{R}} \times \mathcal{D}$, i.e., $$(\hat{r}_n, \hat{d}_n)= \underset{(r,d) \in \widetilde{\mathcal{R}} \times \mathcal{D} }{\text{argmin}} \ {\rho}^*_n(r,d), $$ where $\widetilde{\mathcal{R}}$ may be data-dependent. The final estimator is obtained by the plug-in method as $$\hat{\boldsymbol{\eta}}_n=(\hat{\boldsymbol{\lambda}}^\intercal_n(\hat{r}_n, \hat{d}_n),\hat{r}_n, \hat{d}_n)^\intercal :=(\hat{\boldsymbol{\lambda}}^\intercal_n,\hat{r}_n, \hat{d}_n)^\intercal.$$ Solving the minimization problem in ((ref)) is equivalent to finding the zeros of the weighted least squares estimating equations

where weights take the form $w(\varepsilon_t(\boldsymbol{\eta})) := \psi(|\varepsilon_t(\boldsymbol{\eta})|^{1/2})$. To ensure robustness, such weights must be relatively small when the residual is incompatible with the assumed distribution for the innovations, such as the Gaussian distribution. Solving directly ((ref)) in $\boldsymbol{\lambda}$ may be computationally difficult due to the presence of multiple local minima. This is typically the case for re-descending estimators for which the derivative $\psi(u) = \rho'(u)$ is not monotone. To solve the above computational issues, we propose an iteratively re-weighted least squares (IRLS) approach to compute the estimates. The IRLS algorithm alternates two steps until convergence: (i) computing the weights $\tilde w_t = w(\varepsilon_t(\tilde\boldsymbol{\eta}))$ using the current parameter value, say $\tilde \boldsymbol{\eta}$, and (ii) updating the parameters by solving $ \sum_{t=1}^n \tilde w_t \partial \varepsilon^2_t(\boldsymbol{\eta})/\partial \boldsymbol{\lambda} = \boldsymbol{0} $, which is equivalent to minimizing the weighted residual sum of squares $ \sum_{t=1}^n \tilde w_t \varepsilon^2_t(\boldsymbol{\eta}) $. Note that, for fixed $r$, the parameter update from Step (ii) is just a weighted least squares problem which can be solved efficiently using existing algorithms for TARMA estimation. The above IRLS approach is fast in execution, typically requiring only a few iterations to converge. In all our numerical applications we use the following approach to obtain the initial estimate for the IRLS algorithm. We begin by trimming a percentage of the data corresponding to the most extreme observations; here we choose 10%. Then we run the LS estimator on the trimmed sample. This allows us to obtain a fairly robust initial estimate not affecting the convergence properties of the algorithm. Standard errors for $\hat \boldsymbol{\lambda}_n $ are computed using the asymptotic distribution of the estimator derived in Section (ref). Particularly, $\sqrt{n}\hat \boldsymbol{\lambda}_n$ converges in distribution to a multivariate normal distribution with zero mean and covariance matrix $\boldsymbol{H}(\boldsymbol{\eta})^{-1} \boldsymbol{J}(\boldsymbol{\eta}) \boldsymbol{H}(\boldsymbol{\eta})^{-1}$, where $\boldsymbol{H}(\boldsymbol{\eta})$ and $\boldsymbol{J}(\boldsymbol{\eta})$ are, respectively, the sensitivity and variability matrices whose expression is given in Theorem (ref). The asymptotic variance can be estimated consistently using the sandwich estimator $\hat \boldsymbol{H}(\hat \boldsymbol{\eta}_n)^{-1} \hat \boldsymbol{J}(\hat \boldsymbol{\eta}_n) \hat \boldsymbol{H}(\hat \boldsymbol{\eta}_n)^{-1}$ where

are estimates of the sensitivity and variability matrices $\boldsymbol{H}(\boldsymbol{\eta})$ and $\boldsymbol{J}(\boldsymbol{\eta})$.

Large sample properties

In this section, we study the behavior for the estimator $\hat \boldsymbol{\eta}_n$ as $n$ diverges. We use $\boldsymbol{\eta}_0$ to denote the minimizer of the population objective

In the rest of this section, we assume that $\boldsymbol{\eta}_0$ exists and is unique. Note that here, differently from previous works on TARMA estimation, the true process generating the data does not necessarily coincide with the nominal TARMA model described in Section (ref) and the expectation in ((ref)) may be taken with respect to a process outside the TARMA model family. In this case, the population parameter $\boldsymbol{\eta}_0$ should be regarded as the optimal process in terms of minimizing the density divergence implied by $\rho$ between the parametric TARMA model and the actual process underlying the data.

For the results presented in the remainder of this section, we require the following regularity conditions:

Assumptions (A1), (A2) and (A4) are standard requirements in the threshold framework. Regarding Assumption (A1), more details on the conditions ensuring stationarity and ergodicity of TARMA models are given by Lin99 and Cha19, while invertibility is studied in Cha10. A discussion on the invertibility of threshold moving-average models can also be found in Lin05 and Lin07. Assumption (A3) is a basic requirement for robustness. For instance, the re-descending estimator in Fer12 satisfies these properties for common families of distributions for the innovation process. Another possible choice for $\rho$ leading to similar robustness properties is Tukey's bisquare function maronna2019robust. Assumption (A4) is the same condition considered in Li11a in order to ensure threshold identification. Assumption (A5) is stronger than Assumption (A3), and is needed to guarantee a regular behavior of the expansion leading to asymptotic normality for the ARMA parameter estimator in the two regimes. The loss function $\rho$ should satisfy at least the Fisher consistency property. Namely, when the data are generated by a TARMA process with parameter $\boldsymbol{\eta}_0$, then $\boldsymbol{\eta}_0$ should be also the minimizer of the population objective $\rho^\dag(\boldsymbol{\eta}) = E_{\boldsymbol{\eta}_0}[\rho(\varepsilon_t(\boldsymbol{\eta}))]$, where expectation is taken with respect to the TARMA process with parameter $\boldsymbol{\eta}_0$. Following steps analogous to Lemma 1 in Fer12, one can show that the re-descending estimator minimizing ((ref)) is Fisher consistent for the parameter $\boldsymbol{\eta}$ for any $\alpha > 0$. The special case $\alpha=0$ corresponds to the maximum likelihood estimator, which is clearly Fisher consistent, but does not satisfy Assumption (A3) and leads to estimates that are influenced by outliers. The next theorem shows the strong consistency of the estimator $\hat{\boldsymbol{\eta}}_n$.

In the following, we derive the convergence rates of $\hat{r}_n$ and $\hat{\boldsymbol{\lambda}}_n$ and prove the uniform asymptotic normality of $\hat{\boldsymbol{\lambda}}_n$. To this end, let $\partial\rho(\boldsymbol{\eta}_0)/\partial\boldsymbol{\lambda}$ and $\partial^2\rho(\boldsymbol{\eta}_0)/\partial\boldsymbol{\lambda}\partial\boldsymbol{\lambda}^\intercal$ be the first and the second derivative of the function $\rho(\boldsymbol{\eta})$ with respect to $\boldsymbol{\lambda}$ evaluated at the parameter vector $\boldsymbol{\eta}_0$. Moreover, define

The matrices $\boldsymbol{H}(\boldsymbol{\eta})$ and $\boldsymbol{J}(\boldsymbol{\eta})$, evaluated at $\boldsymbol{\eta}_0$ form the asymptotic variance $\boldsymbol{H}(\boldsymbol{\eta}_0)^{-1}\boldsymbol{J}(\boldsymbol{\eta}_0)\boldsymbol{H}(\boldsymbol{\eta}_0)^{-1}$ for $\hat{\boldsymbol{\lambda}}_n$. Hence, we require the following assumptions.

The proofs of Theorems (ref) and (ref) follow an approach similar to Kou03 and Li11a with some notable differences. While Kou03 also focus on general M-estimators, their results rely heavily on the simpler structure of the autoregressive process, while here we also take into account the moving-average component. Li11a consider both autoregressive and moving-average components, but their proofs are only valid for the specific case of the least squares objective function, which is much simpler to handle than generic M-estimating functions. Finally, note that the estimator of the threshold $\hat r_n$ is super-consistent. In practical terms, this means that the threshold can be taken as given, provided the sample size is adequate. For this reason we have omitted the derivation of the robust asymptotic distribution for $\hat r_n$.

Monte Carlo study

In this section, we perform a Monte Carlo study to assess the performance of our robust estimator. We consider the four parameter settings shown in Table (ref) for the following TARMA$(1,1)$ process:

where $\varepsilon_t \sim N(0,1)$ and $d=1$. The choice of parameters reflects different long-run probabilistic behaviors of the TARMA process. In particular, Cases 1 and 3 correspond to ergodic processes, whereas Cases 2 and 4 correspond to geometrically ergodic processes. Also, Case 3 has unit roots in both regimes but is globally stationary; this is a challenging case laying on the boundary of the ergodicity region; see Cha19 for more details.

The data generated from the clean TARMA process are contaminated using a fraction $\epsilon$ of outliers. In particular, we consider both additive outliers (AOs) and innovation outliers (IOs) corresponding to the two Monte Carlo experiments described below.

To assess the performance of our robust methodology, we compute Monte Carlo estimates of the bias, $ \Vert E(\hat \boldsymbol{\eta}_n) - \boldsymbol{\eta}_0 \Vert^2_2, $ and variance, $ \Vert \text{var}(\hat \boldsymbol{\eta}_n)\Vert^2_2$ of our estimator, where $\boldsymbol{\eta}_0$ represents the parameter vector for the clean TARMA process and $\Vert \cdot \Vert_2$ is the Euclidean norm. Estimates are based on 1000 Monte Carlo replications with sample size $n=100, 200$. In practice, in each contamination setting, we add 10% of equally spaced outliers of size $k=10$ with random sign depending on $\xi_t$. Note that, differently from AOs, IOs are much harder to treat since they enter the state equation and interact non-trivially with the non-linearity of the TARMA process. This can exert a long-term influence upon the series and even produce a qualitative change in the dynamics. The above experiments aim to mimic real scenarios encountered in economics and finance where contamination may occur in both tails but are prevalent in one. Figure (ref) shows bias and variance for the four TARMA specifications under AO contamination for values of the robustness parameter $\alpha = 0,0.3,0.6,0.9,1.2,1.5$, where $\alpha=0$ corresponds to the special case of the non-robust LS estimator. In all the settings, the bias decreases significantly when $\alpha$ moves away from zero and stabilizes for values of $\alpha$ larger than $1$. Cases 2 and 4 show the sharpest decrease, which may be an effect due to the geometric ergodicity. Interestingly, the variance also stabilizes starting for a value of $\alpha$ larger $1$; however, differently from the bias, here Cases 1 and 3 show the sharpest decrease. While both bias and variance of the LS estimator are considerably affected by outliers, the robust estimator with $\alpha \geq 1$ is generally successful in mitigating their influence.

Figure (ref) shows bias and variance for the four TARMA specifications under IO contamination. The findings are consistent with those reported in Figure (ref). In all the scenarios, a value of $\alpha>0$ suffices to improve both bias and variance. The improvement is dramatic in most cases for sufficiently large $\alpha$. Note that the reduction in bias and variance is less marked for Case 3, which sits at the boundary of the parametric region of ergodicity. Even if the process is globally stationary, its regimes are both $I(1)$ so that the outlier effect decays very slowly.

Asymptotic bias under contamination

The asymptotic bias under contamination is a common measure of robustness for the time series framework. Other measures include the influence curve, introduced by hampel1974influence in the i.i.d. framework, which measures the influence of infinitesimal outlier contamination on the parameter estimates. Also, martin1986influence consider a generalization of influence functionals in time-series models based on a replacement outlier model. Here we focus on the asymptotic bias since it does not assume infinitesimal contaminations and provides a realistic representations of the behavior of the estimator in practical situations. Let $\hat \boldsymbol{\eta}_\infty(F)$ be the almost sure limit of the estimator $\hat \boldsymbol{\eta}_n = (\hat \boldsymbol{\lambda}_n^\intercal, \hat r_n , \hat d_n)$ applied to a process with distribution $F$. The asymptotic squared bias for $\hat \boldsymbol{\eta}_\infty$ applied to the contaminated process $\{X_t^{\epsilon, k}\}$ is given by

where $F(X_t^{\epsilon, k})$ is the distribution of the contaminated process $\{X_t^{\epsilon, k}\}$ and $\left\Vert \cdot \right\Vert_2$ is the Euclidean norm. In Figures (ref) and (ref) we show the behavior of the asymptotic bias against outlier size $k$, for contamination levels $\epsilon=0.05, 0.1, 0.15, 0.2$ and different values for the robustness parameter $\alpha$. The asymptotic values are computed using series of size $n=20000$ and the plots summarize the four cases through the median. The behavior for both additive and innovation outliers is similar. The bias of the non-robust estimator ($\alpha=0$) diverges quickly as $k$ increases. On the other hand, for contamination levels up to 10%, small values of $\alpha$ are enough to achieve robustness. As the contamination level increases, larger values of $\alpha$ are needed to stabilize the asymptotic bias. A value of $\alpha$ close to one achieves a remarkable robustness even when 20% of the data are contaminated and in case of large outliers ($\epsilon=0.2$, lower right panels).

Application: Analysis of commodity prices

Commodities are raw materials or primary agricultural products used as inputs in the production of other goods and are commonly traded in the cash market or as derivatives. Commodity prices are extremely important in individual, country-level economies: since they respond quickly to economic shocks, such as increase in demand, they are often used to predict the behavior of other economic variables. We consider $336$ monthly observations for the price of five commonly traded energy or precious commodities. The energy commodities are the WTI crude oil price index, the US natural gas spot price at the Henry Hub in Louisiana, the average of the Australian coal price at Port Thermal in Newcastle and the South African coal price at Richards Bay. The precious commodities are the gold and silver prices traded in London, afternoon fixing. All the series are sampled in the period February 1994 -- December 2021, and are obtained from the World Bank website \url{https://www.worldbank.org/en/research/commodity-markets}. For each commodity, we model log returns of their prices, that is $x_{t,i} = \nabla \log(y_{t,i}) = \log(y_{t,i}/y_{t-1,i})$, where $y_{t,i}$ denotes the price of commodity $i$ ($i=1,\dots,5$) at time $t$. The time plot reported in Figure (ref) of the Supplementary Material highlights that the series have different volatility, which is lower for gold and coal while is more pronounced for oil and natural gas. We estimate TARMA$(1,1)$ models for the five commodities using our robust estimation method and the least squares approach, the latter corresponding to the special case $\alpha = 0$. We also include the linear ARMA$(1,1)$ model, estimated through full maximum likelihood (ML). In preliminary analyses not reported here, we found estimates for the threshold parameter $r$ consistently close to zero for most values of $\alpha$ ranging from 0 to 1, which confirms the general asymmetric behavior in growth and contraction periods, see DL2013. Motivated by these findings, we set $r=0$ to obtain our final TARMA estimates. Moreover, based on macroeconomic theory, we set $d=1$. The overall model accuracy is assessed through the mean absolute percentage error (MAPE) (see Section (ref) if the Supplementary Material), using 12 out-of-sample observations from January to December 2021 as the test set, while the remaining 324 observations are used as the training set. Table (ref) shows parameter estimates for the five series with standard errors in parentheses below the estimates. The third column shows the values of the tuning parameter $\alpha$, computed by minimizing the MAPE over a grid of equally spaced values in the interval $(0,1)$. In all the series, we note that the autoregressive and moving-average parameters change, sometimes dramatically, when using our robust method compared to the LS approach. Moreover, the standard errors from the TARMA models based on the LS method are generally larger than the robust standard errors. Thus, using the non-robust method can hinder the discovery of separate regimes and make it impossible to assess the actual significance of many parameters. For instance, for the silver series, the LS method does not show significantly different estimates in the two regimes and its MAPE is even larger than that of the linear ARMA. On the other hand, the robust TARMA reveals the existence of two dynamical regimes with clearly different autoregressive and moving-average behaviors. The prediction error of the resulting model is 30% smaller than the LS fit. All the estimated intercepts for the robust TARMA are close to zero and this suggests that the transition between the two regimes is not discontinuous (see also Figure (ref), last row). Moreover, in absolute value, the parameters for the lower regime are almost always smaller than those of the upper regime, which highlights the asymmetric behavior of the commodity series characterized by periods of persistent growth and sharp contraction. In particular, the difference in the moving-average parameters denotes the different reaction to shocks in the two regimes: with the exception of coal, in the upper regime the shocks exert a stronger and more persistent influence.

As already mentioned, the higher estimation accuracy of the robust TARMA is also beneficial for prediction since this model always outperforms the least square TARMA. Gains in terms of MAPE are sizeable for precious commodities (up to 67% for gold and 33% for silver); they are moderate for coal (4%) and oil (6%), and small for natural gas (1%). The linear ARMA model is the least accurate, except for gold, where many parameter estimates are not significant. In order to detect the most influential outliers, for each series we compute the robust weights $$ \hat w_t = \dfrac{\exp\{ - \hat \alpha \times \hat \varepsilon^2_{t}/ (2\hat \sigma^2)\}}{\sum_{s=1}^n\exp\{ - \hat \alpha \times \hat \varepsilon^2_{s}/ (2\hat \sigma^2)\}}, \ \ t=1,\dots, 324, $$ where $\hat \varepsilon_t$ is the residual at time $t$ from the robust fit, $\hat \sigma^2$ is the estimated error variance and $\hat \alpha$ is the data-driven tuning parameter obtained by minimizing the MAPE. Smaller weights correspond to observations that are further from the assumed clean model, i.e. the TARMA model with Gaussian errors.

Figure (ref) (top) shows the histograms of the residuals from the robust TARMA for oil, coal and gold prices. All the histograms appear to be different from the nominal standard Gaussian density (superimposed in red) due to heavy tails or asymmetry. The residuals translate into robust weights mostly concentrated on larger values above $0.30$, although a number of observations receives smaller weights closer to zero (see Figure (ref), second row), indicating the presence of strong outliers. Figure (ref) (third row) highlights with circles the most influential outliers corresponding to the smallest robust weights (5% of the sample, or 15 values) in the time plot of the original log-return series. Many of these extreme outliers are evident and correspond to the shocks which occurred during the financial crisis in 2008--2009 and at the beginning of the COVID pandemic, although some of the abrupt changes appear to be compatible with the assumed TARMA model. The last row of Figure (ref) shows the same outliers in the state space (lag plot of $x_{t}$ versus $x_{t-1}$), where we have also added estimated piecewise linear autoregression lines. Note that the 15 most extreme observations are, in fact, outliers in the state space while this is not so evident from the time plot (Figure (ref), third row). Moreover, the placement of such observations marked as outliers appears to be linked to the commodity type. For precious commodities (gold and silver) the outliers tend to fall in the upper regime, while for oil they are found in the lower regime. Finally, for gas and coal, there is roughly the same proportion of outliers in both regimes. Figure (ref) in the Supplementary Material reports the plots for the two remaining commodities (oil and silver).

Conclusions

TARMA models have attracted considerable interest due to their ability to parsimoniously describe complex dynamical features such as jumps, asymmetric limit cycles, time irreversibility, and chaos. They are unique in that they provide a natural interpretation for phenomena that change qualitatively across regimes and react differently to shocks. Nonetheless, estimation for TARMA model is currently limited to the least squares method, which is known to be severely influenced by the presence of outliers. In this paper we provide the first theoretical framework for robust M-estimation for TARMA models and also study its practical relevance. Theorems (ref) and (ref) extend the results of Li11a and establish an asymptotic theory for a wide class of estimators found as the solution of M-estimating equations with bounded derivatives. We establish the superconsistency for $\hat r_n$ and defer to future research the derivation of the limit distribution of the threshold estimator, which is a challenging task. Our results can be used to derive other robust inference and model selection tools for TARMA processes. For example, following ronchetti1997robustness and muller2009robust, a robust model-selection criterion for TARMA models may be formulated as $2 {\rho}_n(\hat \boldsymbol{\eta}_n) + 2 \text{trace}(\hat \boldsymbol{H}^{-1} \hat \boldsymbol{J})$, where $\hat \boldsymbol{H}$ and $\hat \boldsymbol{J}$ are the plug-in estimators based on the sensitivity and variability estimators defined in ((ref)). These can also be used to derive Wald and score statistics to test hypotheses on the parameters. Focusing on the re-descending estimator of Fer12, we study the robustness properties of the proposed M-estimator in a range of scenarios involving both additive and innovation outliers. The results from our Monte Carlo experiments show that moving away from the LS estimator even by a small amount already achieves robustness both in terms of bias and variance. Overall, our estimator reduces considerably the asymptotic bias also in the presence of severe contaminations and high fractions of outliers, where the least squares estimator fails. The findings suggest that robust M-estimation should be generally preferred to the least squares method, even when the actual data deviate only slightly from the nominal TARMA model. The analysis of the five time series of commodity prices shows that the robust TARMA estimates present smaller standard errors and lead to superior forecasting accuracy compared to the least squares fit. This enables us to detect regime changes with confidence and support the hypothesis of a two-regime, asymmetric nonlinearity around zero, characterised by slow expansions and fast contractions. Although a thorough analysis of the price dynamics for different commodities is beyond the scope of the present work, the robust TARMA framework could be used as the foundation for future modelling approaches, possibly leading to important advancements in the field. An interesting direction for future investigations could be the study of the performance in the presence of specific contamination processes. For example IOs are generally more challenging to handle and would require the development of some ad-hoc estimating function. One possible approach is to introduce a robust filtering step within the residual function, as in the bounded innovation propagation ARMA (BIP-ARMA) of muler2009robust. For the time being, we note that there is a fundamental difference in the way non-linear processes react to perturbations compared to linear processes. In general, the presence of dynamic noise can alter qualitatively and non trivially the nature of the process, see e.g., Cha01 for a discussion. For instance, for linear processes the response function to noise is flat, whereas non-linear processes can act both as noise amplifiers and noise suppressors, producing a plethora of characteristic phenomena, such as resonances, or the state-dependence predictability, which is well known in the forecasting literature, see e.g., Fan05, Ch. 10.

\setcounter{table}{0} \setcounter{section}{0} \setcounter{page}{1}