The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
75,626 characters
Testing for threshold regulation in presence of measurement error with an application to the PPP hypothesis.
\def\spacingset#1{\renewcommand{\baselinestretch}
{#1}\small\normalsize} \spacingset{1}
\if00
{
\title{\bf Testing for threshold regulation in presence of measurement error with an application to the PPP hypothesis.}
\author{Kung-Sik Chan
\hspace{.2cm}\\
{Department of Statistics and Actuarial Science, The University of Iowa, USA}\\
and \\
Simone Giannerini \\
Department of Statistical Sciences, University of Bologna, Italy\\
and \\
Greta Goracci \\
Department of Statistical Sciences, University of Bologna, Italy\\
Faculty of Economics and Management, Free University of Bozen/Bolzano, Italy\\
and \\
Howell Tong\\
University of Electronic Science and Technology of China, Chengdu, China;\\
Tsinghua University, China;
London School of Economics and Political Science, U.K.
}
\maketitle
} \fi
\if10
{
\bigskip
\bigskip
\bigskip
\begin{center}
{\LARGE\bf Testing for threshold regulation in presence of measurement error with an application to the PPP hypothesis.}
\end{center}
\medskip
} \fi
\vspace*{-0.95cm}
\begin{abstract}
Regulation is an important feature characterising many dynamical phenomena and can be tested within the threshold autoregressive setting, with the null hypothesis being a global non-stationary process. Nonetheless, this setting is debatable since data are often corrupted by measurement errors. Thus, it is more appropriate to consider a threshold autoregressive moving-average model as the general hypothesis. We implement this new setting with the integrated moving-average model of order one as the null hypothesis. We derive a Lagrange multiplier test which has an asymptotically similar null distribution and provide the first rigorous proof of tightness pertaining to testing for threshold nonlinearity against difference stationarity, which is of independent interest. Simulation studies show that the proposed approach enjoys less bias and higher power in detecting threshold regulation than existing tests when there are measurement errors. We apply the new approach to the daily real exchange rates of Eurozone countries. It lends support to the purchasing power parity hypothesis, via a nonlinear mean-reversion mechanism triggered upon crossing a threshold located in the extreme upper tail. Furthermore, we analyse the Eurozone series and propose a threshold autoregressive moving-average specification, which sheds new light on the purchasing power parity debate.
\end{abstract}
\noindent
{\it Keywords:} Lagrange multiplier test; Random walk hypothesis; Threshold autoregressive moving-average model; Purchasing power parity; Volatility.
\vfill
\spacingset{1.1}
\section{Introduction} \label{sec:intro}
Regulation plays a fundamental role in various fields including economics, finance, biological growth and population fluctuations, etc. Growth processes are generally regulation-free until they enter into extreme phases. For instance, there are strong theoretical economic arguments to support the view that the price gap for goods of a particular kind (or a basket of goods) in different countries should rapidly converge to zero. In other words, currencies should have the same long-run purchasing power implying that the real exchange rates would be regulated, and that there exists a threshold that triggers a mean‐reversion \citep{Tay01b}. However, empirical evidence points to a strong persistence of the price gap and the inability of existing unit-root tests to reject the null hypothesis of a random walk. This is also known as the Purchasing Power Parity (PPP) puzzle, which has been a subject of debate ever since its first formulation. See, e.g., \cite{Tay04} for a discussion.
\par
The random walk is a simple model for regulation-free dynamics. On the other hand, regulation from above (below) may be captured by a first-order threshold autoregressive model (TAR) which follows a random walk until the process crosses a certain threshold above (below) which mean-reversion takes place, while the process as a whole is \emph{stationary}. A nonlinear stationary process generally renders nonlinear and state-dependent the impulse response to a random shock, which is consequential and could be leveraged in economic regulation. Thus, an approach to test for dynamic regulation is to adopt the preceding threshold model as the general model and test whether it reduces to a \emph{global} random walk. This approach has received much attention in the literature \citep{End98,Can01,Bec04,Kap06,Bec08,Seo08,Par16,deJ07,Gio17}. However, data are almost always corrupted by measurement error. In this case, the TAR model is not appropriate and the null hypothesis should be an exponential smoothing model instead, i.e., the integrated moving-average IMA(1,1) model. Then, the general hypothesis may be taken as the first-order threshold autoregressive moving-average model, i.e., TARMA(1,1), which is driven by an IMA(1,1) model in one of its two regimes. See Section~\ref{sec:Merr} of the Supplementary Material for further justification. Above all, we cannot over-emphasize the critical importance of the role of the moving average term for practical applications.
\par
Just as ARMA models provide a parsimonious approximation to some long AR models, so may TARMA models well approximate some high-order TAR models parsimoniously \cite{Gor20,Gor21}. Thus, the TARMA model holds substantial promise as a class of nonlinear time series models for exploring nonlinear dynamics in economics and other fields. Yet, the TARMA model has been under-explored, partly because of a lack of progress in obtaining conditions on stationarity and ergodicity. Unlike the AR-ARMA analogy, the incorporation of a moving-average part in a nonlinear framework poses major theoretical challenges and has non-trivial implications on the probabilistic structure of the process. Recent work by \cite{Cha19} provides, for the first time, a breakthrough in deriving a set of necessary and sufficient conditions for the (multi-regime) TARMA(1,1) model to admit an irreducible and invertible state-space representation. Moreover, they derived a set of necessary and sufficient conditions for stationarity and ergodicity of the TARMA(1,1) model.
\par
By leveraging on the recent results of \cite{Cha19}, we deploy a supremum Lagrange Multiplier test (supLM) for threshold regulation, with the TARMA(1,1) model as the general framework. We specify an IMA$(1,1)$ model as the null hypothesis and a TARMA$(1,1)$ with a unit-root regime as the alternative. A difficulty arising from testing for a unit-root against a TARMA model is that the threshold parameters are absent under the null hypothesis.
This non-standard situation, in the nonlinear time series context, is well recognized both in the TAR setting \citep{Cha90,Han96} and in the TARMA setting \cite{Li11,Gor21c}. Fortunately, the supLM framework overcomes this problem. We derive its asymptotic distribution both under the null hypothesis and local alternatives. We prove that the test is consistent and asymptotically similar in that its asymptotic null distribution does not depend on the value of the MA parameter. Moreover, we provide the first rigorous proof of tightness pertaining to testing for threshold nonlinearity against difference stationarity. The tightness result is of independent interest. It constitutes a general theoretical framework for ARIMA versus TARMA testing. We also introduce a wild bootstrap version of the supLM statistic that, for finite samples, possesses good properties and robustness against heteroskedasticity. We perform a large scale simulation study to compare our tests with existing tests, in which the alternative hypothesis is that of a threshold model. In general, the size of the latter tests is severely biased in a number of cases to the extent that their use in practical applications remains questionable unless additional information on the data generating process is available. The surprisingly good size property of our tests may be owing to the versatility of an IMA(1,1) model in approximating general non-seasonal difference stationary processes. In addition, the comparison includes some of the best performing unit-root tests to date, where the alternative hypothesis does not specify explicitly a nonlinear process.
\par
The paper is structured as follows. In Section~\ref{sec:TARMA} we present some fundamentals of the first-order TARMA model and a parametrization that reduces to the IMA(1,1) process under the null hypothesis. In Section~\ref{sec:LM} we present a supremum Lagrange Multiplier test, which we denote by supLM, including the theoretical framework based on Brownian local time. Section~\ref{sec:null} is devoted to the derivation of the asymptotic distribution of the supLM test statistic under the null hypothesis and we show that it is nuisance-parameter-free and depends only on the search range of the threshold. The results concerning the local power of the proposed test are summarized in Section~\ref{sec:LP}. In Section~\ref{sec:sim} we perform a large scale simulation study to show the performance of the supLM test and a wild bootstrap version of it and compare them with numerous existing tests in the recent literature. Section~\ref{sec:real} presents the application of the tests to the daily real exchange rates of the panel of Eurozone countries, plus the TARMA-GARCH modelling for the global Eurozone series. All the proofs are collected in the Supplementary Material, which contains further results from the Monte Carlo study and from the real data application.
\section{Threshold autoregressive moving-average model}\label{sec:TARMA}
Consider the following first-order threshold autoregressive moving-average (TARMA) model:
\begin{equation}
X_t=\begin{cases}
\phi_{1,0}+\phi_{1,1}X_{t-1} +\varepsilon_t- \theta \varepsilon_{t-1}, & \text{if } X_{t-d} \le r \\
\phi_{2,0}+\phi_{2,1}X_{t-1} +\varepsilon_t- \theta \varepsilon_{t-1} & \text{otherwise},
\end{cases}
\label{tarma.r}
\end{equation}
where $\phi_{2,1}$ is fixed at 1 unless stated otherwise,
the innovations $\{\varepsilon_t\}$ are independent and identically distributed random variables with zero mean and variance $\sigma^2$, $\varepsilon_t$ is independent of $X_{t-j}, j\ge 1$, the delay $d$ is a positive integer which, for simplicity, is taken to be 1 henceforth, $r$ is the real-valued threshold parameter, and the $\phi$'s and $\theta$'s are unknown coefficients. The assumption of independence and identical distribution of the innovations will be relaxed later to a martingale difference sequence. The preceding (constrained) TARMA model assumes that the sub-model in the upper regime is a first-order IMA model while the lower regime specifies a general first-order ARMA model. Statistical inference with a TARMA model hinges on whether the model is invertible. We assume $|\theta|<1$ since it is a necessary and sufficient condition for the invertibility of Model~(\ref{tarma.r}) \citep{Cha10}. By assuming that the innovations admit a positive, continuous probability density function with finite absolute first moment, \cite{Cha19} showed that Model~(\ref{tarma.r}) is an ergodic Markov chain if and only if $\phi_{2,0}<0$ and either $(i)$ $\phi_{1,1}<1$, or $(ii)$ $\phi_{1,1}=1$, $\phi_{1,0}>0$; ergodicity then implies that the first-order TARMA model admits a unique stationary distribution. Furthermore, under the stronger condition that the innovations admit a finite absolute $k$th moment for some $k>2$, \cite{Cha19} provides a complete classification of the parametric regions of Model~(\ref{tarma.r}) into sub-regions of ergodicity, null recurrence and transience. In particular, the (constrained) first-order TARMA model defined by Model~(\ref{tarma.r}) is null-recurrent if any of the following holds:
$(iii) \,\, \phi_{1,1}=1, \phi_{2,0}=0, \phi_{1,0}\geq0; \quad
(iv) \,\, \phi_{1,1}=1, \phi_{2,0}<0, \phi_{1,0}=0; \quad
(v) \,\, \phi_{1,1}<1, \phi_{2,0}=0.$
If none of the conditions $(i)$--$(v)$ holds, then the model is transient. Therefore, Model~(\ref{tarma.r}) is a rich model that encompasses both linear and nonlinear processes spanning a wide spectrum of long-run behaviors including ergodicity, null recurrence and transience.
\section{Lagrange multiplier test for threshold regulation}\label{sec:LM}
We first formulate a framework for testing for threshold regulation from below.
Let $\{X_t, t=0,1,\ldots\}$, be a time series and assume that, for $t\ge 1$, $X_t$ satisfies the equation
\begin{equation}
H:\quad X_t=\phi_0+X_{t-1} + \varepsilon_{t}-\theta\varepsilon_{t-1} + (\phi_{1,0}+\phi_{1,1} X_{t-1})\times I(X_{t-1}\le r),
\label{reparm.model}
\end{equation}
which is a re-parameterization of Model~(\ref{tarma.r})
with
$\phi_0=\phi_{2,0}$ and by an abuse of notation, $\phi_{1,0}$ and $\phi_{1,1}$ represent, respectively, the difference of intercept and slope of the lower regime relative to their upper-regime counterparts;
the initial value $X_0$ can be fixed at, say, 0. Our interest is in testing whether $\phi_{1,0}=\phi_{1,1}=0$, in which case the data are generated by
the IMA(1,1) model
\begin{equation}
H_0:\quad X_t=\phi_0+X_{t-1}+\varepsilon_t-\theta\varepsilon_{t-1},
\label{reparm.modelH0}
\end{equation}
where $|\theta|<1$.
If the intercept $\phi_0\not = 0$, then the IMA(1,1) process has a linear trend. If no such linear trend is apparent in the data, it is reasonable to omit the intercept. Henceforth, we assume that $\phi_0=0$ under $H_0$. The case for $\phi_0\not= 0$ will be studied elsewhere. However, the intercept terms on the two regimes of any competing stationary first-order TARMA model will be required to model the mean of the data. Indeed, even for mean-deleted data, the intercept terms of the first-order TARMA model are not necessarily zero. Thus, the intercept terms are essential and retained in the constrained TARMA model under $H$. Note that testing for threshold regulation from above can be conducted by applying the test to $\{-X_t\}$.
\par
Under the null hypothesis, the threshold parameter is absent thereby complicating the test \citep{Cha90, Han96, Li11, Gor21c}. Our approach is to develop a Lagrange multiplier test statistic for $H_0$ initially with the threshold parameter fixed at some $r$. Denote the test statistic as $T_{n}(r)$. Since $r$ is unknown and indeed absent under $H_0$, we shall compute $T_{n}(r)$ for all $r$ over some data-driven interval, say, $[a,b]$ with the end points being some percentiles of the observed data. For instance, $a$ could be the 20-th percentile and $b$ the 80-th percentile. Then the overall test statistic results in $T_n=\sup_{r\in [a,b]} T_{n}(r)$. Besides taking the supremum, other approaches including integration can be employed to derive an overall test statistic.
\par
For fixed $r$, the Lagrange multiplier test is developed based on the Gaussian likelihood conditional on $X_0$:
\begin{equation}
\ell = -\log (2\pi\sigma^2)\times n/2 -\sum_{t=1}^n \varepsilon_t^2/(2\sigma^2),
\label{log-like}
\end{equation}
where, by an abuse of notation,
\begin{equation}
\varepsilon_{t}=X_t-\{\phi_0+X_{t-1} + (\phi_{1,0}+\phi_{1,1} X_{t-1})\times I(X_{t-1}\le r)\}+\theta\varepsilon_{t-1},\quad \forall t\ge 1, \label{residual}
\end{equation}
with the unknown $\varepsilon_0$ set to be zero; $\varepsilon_t$ in the preceding formula is a function of $\phi_0,\phi_{1,0}, \phi_{1,1}, \theta$ and $r$, but the arguments are generally suppressed for simplicity.
Let $\boldsymbol{\psi}=(\phi_0,\theta, \sigma^2,$ $ \phi_{1,0}, \phi_{1,1})^{\intercal}$, with its components denoted by $\psi_j, j=1,2,\ldots,5$, and let it be partitioned into $\boldsymbol{\psi}_1=(\phi_0,\theta, \sigma^2)^{\intercal}$ and $\boldsymbol{\psi}_2=(\phi_{1,0}, \phi_{1,1})^{\intercal}$. The null hypothesis can be succinctly expressed as $H_0: \boldsymbol{\psi}_2=0$.
The score vector is
\begin{align*}
\frac{\partial \ell}{\partial \psi_j}&= -\sum_{t=1}^n \frac{\varepsilon_t}{\sigma^2} \frac{\partial \varepsilon_t}{\partial \psi_j},\quad
1\le j \le 5, j \not = 3, \\
\frac{\partial \ell}{\partial \psi_3}&= \frac{\partial \ell}{\partial \sigma^2} =\sum_{t=1}^n \frac{\varepsilon_t^2-\sigma^2}{2\sigma^4}
\end{align*}
where for $t> 1$,
\begin{align}
\frac{\partial \varepsilon_t}{\partial \phi_0} &=
-1+\theta \frac{\partial \varepsilon_{t-1}}{\partial \phi_0} = -\sum_{j=0}^{t-1}\theta^j,\label{eq:ppsi0}\\
\frac{\partial \varepsilon_t}{\partial \theta} &= \varepsilon_{t-1}+\theta \frac{\partial \varepsilon_{t-1}}{\partial \theta}
=\sum_{j=0}^{t-1}\theta^j\varepsilon_{t-1-j}, \label{eq:ptheta}\\
\frac{\partial \varepsilon_t}{\partial \phi_{1,0}} &= -I(X_{t-1}\le r)+ \theta \frac{\partial \varepsilon_{t-1}}{\partial \phi_{1,0}}
=-\sum_{j=0}^{t-1}\theta^jI\left(X_{t-1-j}\le r\right), \label{eq:pphi10}\\
\frac{\partial \varepsilon_t}{\partial \phi_{1,1}} &= -X_{t-1}I(X_{t-1}\le r)+ \theta \frac{\partial \varepsilon_{t-1}}{\partial \phi_{1,1}}
=-\sum_{j=0}^{t-1}\theta^jX_{t-1-j}I\left(X_{t-1-j} \le r\right), \label{eq:pphi11}
\end{align}
with initial values given by
$\partial \varepsilon_1/\partial \phi_0=-1, \partial \varepsilon_1/\partial \theta=0,
\partial \varepsilon_1/\partial \phi_{1,0}=-I(X_0\le r)$ and $\partial \varepsilon_1/\partial \phi_{1,1}=-X_0I(X_0\le r)$. Below, we sometimes write, as a typical example,
$\partial \varepsilon_t/\partial \phi_{1,1}=-(1-\theta B)^{-1} \left\{ X_{t-1}I(X_{t-1} \le r)\right\}$, where $B$ is the backshift operator that shifts the indices backward by 1 unit. The IMA(1,1) model under the null hypothesis can be estimated by solving the score equation
$
\partial \ell/\partial \boldsymbol{\psi}_1=0,
$
yielding $\hat{\boldsymbol{\psi}}_1=\hat{\boldsymbol{\psi}}_{1,n}=(\hat{\phi}_{0,n}, \hat{\theta}_n, \hat{\sigma}^2_n)^{\intercal}$. Thus, the overall estimator of $\boldsymbol{\psi}$ under $H_0$ is $\hat{\boldsymbol{\psi}}=(\hat{\phi}_{0,n}, \hat{\theta}_n, \hat{\sigma}^2_n, 0,0)^{\intercal}$, with the residuals given by
\begin{equation}
\hat{\varepsilon}_{t}=X_t-X_{t-1}-\hat{\phi}_0 +\hat{\theta}\hat{\varepsilon}_{t-1},\quad \forall t\ge 1, \label{residualunderH0}
\end{equation}
where $\hat{\varepsilon}_0=0$. The observed Fisher information (excluding the threshold parameter) is given by
$$I_n = -\frac{\partial^2 \ell}{\partial \boldsymbol{\psi} \partial \boldsymbol{\psi}^{\intercal}}, $$
whose $(i,j)$-th element with $i,j\not= 3$ is given by
\begin{equation}
\sum_{t=1}^n
\frac{1}{\sigma^2}\frac{\partial \varepsilon_t}{\partial \psi_i}\frac{\partial \varepsilon_t}{\partial \psi_j} + \sum_{t=1}^n
\frac{\varepsilon_t}{\sigma^2} \frac{\partial^2 \varepsilon_t}{\partial \psi_i \partial\psi_j} = (1+o_p(1))\times\sum_{t=1}^n
\frac{1}{\sigma^2}\frac{\partial \varepsilon_t}{\partial \psi_i}\frac{\partial \varepsilon_t}{\partial \psi_j},
\label{info1}
\end{equation}
its $(3,i)$-th element with $i\not = 3$ equal to
\begin{equation}
\sum_{t=1}^n \frac{\varepsilon_t}{\sigma^4}
\frac{\partial \varepsilon_t}{\partial \psi_i}=o_p(n)
\label{info2}
\end{equation}
and the $(3,3)$-th element equal to
$$
\sum_{t=1}^n \left(\frac{1}{2\sigma^4}-\frac{\varepsilon_t^2}{\sigma^6}\right),
$$
where the $o_p(1)$ and $o_p(n)$ terms hold uniformly in $r$, when the expressions are evaluated at the true parameter value under the null hypothesis; hence they are asymptotically negligible (via arguments similar to those in the proof of Theorem~\ref{thm1}), and omitted in all numerical work reported herein. Partition the Fisher information matrix according to $\boldsymbol{\psi}_i, i=1,2$ into
\begin{equation}\label{eq:In}
I_n =
\begin{pmatrix}
I_{1,1, n} & I_{1,2,n}\\
I_{2,1,n} & I_{2,2,n}
\end{pmatrix}.
\end{equation}
Note $\partial \ell/\partial \psi_j, \partial \varepsilon_t/\partial \psi_i, I_n,
$ depend on $\boldsymbol{\psi}$ and $r$ implicitly. Below, we sometimes write, e.g., $\partial \ell/\partial \psi_j(\boldsymbol{\psi};r)$, to highlight the role of the arguments; we further simplify the notation, for example, from $\partial \ell/\partial \psi_j(\boldsymbol{\psi}_0;r)$ to $\partial \ell/\partial \psi_j(r)$, with $\boldsymbol{\psi}_0$ denoting the true value under $H_0$. Moreover, $I_{1,1,n}(\boldsymbol{\psi}_0;r)$ and $\partial \ell/\partial \boldsymbol{\psi}_1(\boldsymbol{\psi}_0;r)$ are further simplified as $I_{1,1,n}$ and $\partial \ell/\partial \boldsymbol{\psi}_1$ as they do not depend on $r$. By an abuse of notation, the true values of the moving-average coefficient and the innovation variance under $H_0$ are simply denoted by $\theta$ and $\sigma^2$; no confusion should arise as the context will make clear whether they stand for the generic parameters or their true values.
\par
The Lagrange multiplier test statistic is an asymptotic approximation of twice the Gaussian likelihood ratio statistic, based on a second-order Taylor expansion. For fixed $r$, it equals
\begin{equation}
T_{n}(r)=\frac{\partial \hat{\ell} }{\partial \boldsymbol{\psi}_2^{\intercal}}(r) \left\{\hat{I}_{2,2,n}(r)-\hat{I}_{2,1,n}(r)\hat{I}_{1,1,n}^{-1}(r) \hat{I}_{1,2,n}(r)\right\}^{-1}\frac{\partial \hat{\ell} }{\partial \boldsymbol{\psi}_2}(r) \label{LM-stat}
\end{equation}
where $\partial \hat{\ell}/\partial \boldsymbol{\psi}_2(r)$ is equal to $\partial \ell/\partial \boldsymbol{\psi}_2$ evaluated at $\boldsymbol{\psi}_1=\hat{\boldsymbol{\psi}}_1$, $\boldsymbol{\psi}_2=0$ and the threshold parameter $r$. Similarly defined are $\hat{I}_{i,j,n}(r), 1\le i,j,\le 2$. Because the threshold $r$ is unknown, the overall supLM statistic is $T_n=\sup_{r\in [a, b]} T_{n}(r)$ with $a$ and $b$, for instance, being some pre-specified percentiles of the observed data. For theoretical analysis, the threshold range is specified as $R_n=(n^{1/2}(1-\theta) \sigma \times r_L, n^{1/2}(1-\theta) \sigma \times r_U)$ where $r_L<r_U$ are two fixed finite numbers. We now justify this choice of the threshold range. First, some heuristics will be employed.
Under the null hypothesis (with $\phi_0=0$),
\begin{equation*}
X_t =
\varepsilon_t+(1-\theta)\sum_{s=1}^{t-1} \varepsilon_s-\theta\varepsilon_0+ X_0.
\end{equation*}
Hence, $\{n^{-1/2}\,X_{[sn]}, 0\le s\le 1\}$, where $X_{[sn]} = \sum_{t=1}^{[sn]} X_t $ and $[sn]$ is the largest integer less than or equal to $sn$, converges in distribution to $\{(1-\theta)\sigma W_s\}$ where $\{W_s\}$ is the standard Brownian motion. It is well known \citep[Theorems 3.1 and 3.2]{Bjo19} that the Brownian local time $\{L_t^x, t\ge 0, -\infty<x<\infty\}$ defined as follows:
$$
L_t^x=|W_t-x|-|x|-\int_0^t \operatorname{sign}(W_s-x)ds,
$$
where $\operatorname{sign}(x)$ denotes the sign of $x$, is essentially the probability density function of the Brownian realization in the sense that for any bounded real-valued Borel function $f$,
$$
\int_0^1 f(W_s)ds= \int_{-\infty}^\infty f(x) L_1^x dx.
$$
Thus, any quantile of $\{X_t, t=0,\ldots,n\}$ is asymptotically equal to $n^{1/2}(1-\theta)\sigma$ times the corresponding quantile of $\{W_s, 0\le s\le 1\}$. Since the Brownian local time process is a random process, so the quantiles are realization specific! This motivates us to set the threshold to be of the form $r_n=(1-\theta)\tau\sigma n^{1/2}$ for some fixed $\tau$, in which case
\begin{equation}
n^{-1/2} \frac{\partial \ell}{\partial \phi_{1,0} }(r_n)=
n^{-1/2}\sum_{t=1}^n \frac{\varepsilon_t}{\sigma^2} \frac{1}{1-\theta B}\left\{ I\left(\frac{X_{t-1}}{n^{1/2}(1-\theta)\sigma}\le \tau\right)\right\}.
\label{score10}
\end{equation}
The right side of (\ref{score10}) is a Riemann-Stieltjes sum over $[0,1]$, with a step integrator jumping at $t/n$ with jump size $ (n\sigma^{2})^{-1/2} \varepsilon_t$ and the integrand is a piecewise constant function which equals $\sum_{j=0}^{t-1}\theta^j I\left(\{n^{1/2}(1-\theta)\sigma\}^{-1}X_{t-1-j} \le \tau \right)$ over the interval $[n^{-1}(t-1), n^{-1}t]$, for $t=1,2,\ldots,n$. The integrator converges weakly to the standard Brownian motion whereas the integrand to $(1-\theta)^{-1} I(W_s\le \tau)$ as $t, n\to \infty$ such that $t/n\to s$ in $[0,1]$. Thus, heuristically, $ n^{-1/2} \partial \ell/\partial \phi_{1,0}(r_n)$ converges in distribution to $(1-\theta)^{-1} \sigma\int_0^1 I(W_s\le \tau )dW_s $
under $H_0$ and as $n\to\infty$, or in symbol,
\begin{equation}
n^{-1/2} \frac{\partial \ell}{\partial \phi_{1,0} }(r_n) \rightsquigarrow \frac{1}{(1-\theta) \sigma}\int_0^1 I(W_s\le \tau )dW_s.
\label{wconv1}
\end{equation}
This asymptotic result and other heuristic results stated below can be essentially justified using Theorem 7.10 in \cite{Kur96}.
Similarly,
\begin{align}
n^{-1} \frac{\partial \ell}{\partial \phi_{1,1} }(r_n)&= n^{-1/2}\sum_{t=1}^n \frac{\varepsilon_t}{\sigma} \frac{1}{1-\theta B}\left[\frac{X_{t-1}}{n^{1/2}\sigma} I\left\{\frac{X_{t-1}}{n^{1/2}(1-\theta)\sigma}\le \tau\right\}\right]\nonumber \\
&\rightsquigarrow \int_0^1 W_sI(W_s\le \tau )dW_s \label{wconv2} \\
n^{-1/2} \frac{\partial \ell}{\partial \phi_{0} }&=
n^{-1/2} \sum_{t=1}^n \frac{\varepsilon_t}{\sigma^2} \frac{1}{1-\theta B} (1) \rightsquigarrow \frac{1}{(1-\theta)\sigma}\int_0^1 dW_s=\frac{W_1}{(1-\theta)\sigma}.
\label{wconv3}
\end{align}
Note the different rates of normalization. Let $K_n$ be the $5\times 5$ diagonal matrix with the last diagonal elements being $n$ and other diagonal elements all being $n^{1/2}$. We can also show that
$
K_n^{-1} I_n(r_n) K_n^{-1}$ converges in probability to a matrix denoted by $\mathcal{I}(\tau)$ which can be blocked as $I_n$ (see Eq.~\ref{eq:In}). In particular, $\mathcal{I}_{1,1}$ is a diagonal matrix comprising $(1-\theta)^{-2}\sigma^{-2}, (1-\theta^2)^{-1}, (4\sigma^4)^{-1}$ as its diagonal elements,
\begin{equation*}
\mathcal{I}_{2,2}(\tau)=
\begin{pmatrix}
\frac{1}{(1-\theta)^2 \sigma^2}\int_0^1 I(W_s\le \tau) ds &
\frac{1}{(1-\theta) \sigma}\int_0^1 W_sI(W_s\le \tau) ds\\
\frac{1}{(1-\theta) \sigma}\int_0^1 W_sI(W_s\le \tau) ds & \int_0^1 W_s^2 I(W_s\le \tau) ds
\end{pmatrix};
\label{I22}
\end{equation*}
\begin{equation*}
\mathcal{I}_{2,1}(\tau)=
\begin{pmatrix}
\frac{1}{(1-\theta)^2 \sigma^2}\int_0^1 I(W_s\le \tau) ds &
0 & 0\\
\frac{1}{(1-\theta) \sigma}\int_0^1 W_sI(W_s\le \tau) ds & 0 & 0
\end{pmatrix}.
\label{I21}
\end{equation*}
Note that $\mathcal{I}_{1,1}$ does not depend on $\tau$.
Thus, $\theta$ and $\sigma^2$ are locally orthogonal to the other parameters around the true parametric value under $H_0$ . Hence, their estimates are expected to be asymptotically independent of the proposed test statistic, as will be shown below to be the case.
\section{The null distribution}\label{sec:null}
We now derive the asymptotic distribution of $T_{n}(r)$ under the null hypothesis of an IMA(1,1) model with zero intercept.
Using second-order Taylor expansion and after some routine algebra, it holds that
\begin{equation}
\frac{\partial \hat{\ell} }{\partial \boldsymbol{\psi}_2} (r_n)
\approx\frac{\partial \ell}{\partial \boldsymbol{\psi}_2}(r_n)-I_{2,1,n}(r_n)I_{1,1,n}^{-1}\frac{\partial \ell}{\partial \boldsymbol{\psi}_1}. \label{eq:scoreapprox}
\end{equation}
More rigorously, letting
$$Q_n=\begin{pmatrix}
n^{-1/2} & 0 \\
0 & n^{-1}
\end{pmatrix}, \quad
P_n=n^{-1/2}\begin{pmatrix}
1 & 0 &0 \\
0 & 1 & 0\\
0 & 0 & 1
\end{pmatrix},
$$
we shall prove below that uniformly for $ r_n=n^{1/2}(1-\theta) \sigma\tau\in R_n=(n^{1/2}(1-\theta) \sigma \times r_L, n^{1/2}(1-\theta) \sigma \times r_U)$, where $r_L< r_U$ are fixed numbers,
\begin{align}
Q_n\frac{\partial \hat{\ell} }{\partial \boldsymbol{\psi}_2}(r_n)
&=Q_n\frac{\partial \ell}{\partial \boldsymbol{\psi}_2}(r_n)
-I_{2,1}(\tau)I_{1,1}^{-1}
P_n\frac{\partial \ell}{\partial \boldsymbol{\psi}_1}+o_P(1) \nonumber \\
&= Q_n\frac{\partial \ell}{\partial \boldsymbol{\psi}_2}(r_n)
-\tilde{I}_{2,1}(\tau)\tilde{I}_{1,1}^{-1}
P_n\frac{\partial \ell}{\partial \phi_0}+o_P(1),
\label{eq:exactscoreapprox}
\end{align}
where,
owing to the form of $I_{2,1}(\tau)$,
$\tilde{I}_{1,1}= (1-\theta)^{-2}\sigma^{-2}$ and
$$
\tilde{I}_{2,1}=
\begin{pmatrix}
\frac{1}{(1-\theta)^2 \sigma^2}\int_0^1 I(W_s\le \tau) ds \\
\frac{1}{(1-\theta) \sigma}\int_0^1 W_sI(W_s\le \tau) ds
\end{pmatrix}.
$$
The intercept $\hat{\phi}_{0,n}$ admits the asymptotic representation under $H_0$ \citep[c.f. Eqn.(8.11.5)]{Bro01}
$$P_n^{-1}(\hat{\phi}_{0,n}-\phi_0)= (\tilde{I}_{1,1})^{-1} P_n \frac{\partial \ell}{\partial \phi_0}+o_P(1).$$
A key step in deriving the limiting null distribution of the proposed test is then to demonstrate that uniformly for $r_n=n^{1/2}(1-\theta) \sigma\tau\in R_n$
\begin{equation}
Q_n\frac{\partial \hat{\ell} }{\partial \boldsymbol{\psi}_2}(r_n)
= Q_n\frac{\partial \ell}{\partial \boldsymbol{\psi}_2}(r_n)-\tilde{I}_{2,1}(\tau)P_n^{-1}(\hat{\phi}_{0,n} -\phi_0)+o_p(1). \label{essence}
\end{equation}
Let
\begin{equation}
H(\tau)=\left(\int_0^1 dW_s, \int_0^1 I(W_s\le \tau) dW_s, \int_0^1 W_sI(W_s\le \tau) dW_s\right)^{\intercal}
\label{eq:Htau}
\end{equation}
and
\begin{equation}
\Lambda(\tau)
=
\begin{pmatrix}
1 & \int_0^1 I(W_s\le \tau)ds & \int_0^1 W_sI(W_s\le \tau)ds \\
\int_0^1 I(W_s\le \tau)ds & \int_0^1 I(W_s\le \tau)ds & \int_0^1 W_sI(W_s\le \tau)ds \\
\int_0^1 W_sI(W_s\le \tau) ds & \int_0^1 W_sI(W_s\le \tau) ds & \int_0^1 W_s^2I(W_s\le \tau) ds
\end{pmatrix}.\label{eq: Lambdatau}
\end{equation}
Let $\Lambda(\tau)$ be partitioned into a $2\times 2$ block matrix with the $(2,2)$-th block being $2\times 2$. Similarly partitioned is $H(\tau)=(H_{1}(\tau), H_{2}(\tau))^{\intercal}$. It follows from Eq.~(\ref{eq:exactscoreapprox}) and Eqs.~(\ref{wconv1})--(\ref{wconv3}) that the asymptotic null distribution of $T_{n}(r_n)$ can be shown to be the same as that of $\left\|(\{\Lambda^{-1}(\tau)\}_{2,2})^{1/2} \left( H_{2}(\tau)-\Lambda_{2,1}(\tau)H_{1}(\tau)\right)\right\|^2$, where $\|\cdot\|^2$ is the squared Euclidean norm of the enclosed vector. It is readily shown that $\{\Lambda^{-1}(\tau)\}_{2,2}=\{\Lambda_{2,2}(\tau)-\Lambda_{2,1}(\tau)\Lambda_{1,2}(\tau)\}^{-1}$. The asymptotic null distribution of $T_n$ is derived in the following theorem.
\begin{theorem}
\label{thm1}
Suppose $H_0$ holds so that $\{X_t, t=0,1,\ldots,\}$ is an IMA(1,1) process satisfying Eq.~(\ref{reparm.modelH0}), with the intercept $\phi_0=0$, $|\theta|<1$ and the innovations are independent and identically distributed with zero mean and finite positive variance.
Let $r_L< r_U$ be two fixed real numbers.
Let $$\mathcal{T}_n(\tau)= n^{-1/2}\sum_{t=2}^{n} \frac{\varepsilon_t}{\sigma}\sum_{j=0}^{t-2}\theta^j I\left\{ r_L<\frac{X_{t-1-j}}{n^{1/2}(1-\theta)\sigma}\leq\tau\right\},$$ for $ r_L\le \tau\le r_U$.
Suppose $(i)$ there exists a constant $C>0$ such that, for any fixed $r_L\leq\tau_1<\tau_2\leq r_U$,
\begin{align}
E\left\{\left|\mathcal{T}_n(\tau_2)-\mathcal{T}_n(\tau_1)\right|^4\right\}\leq C(|\tau_2-\tau_1|^{3/2}+|\tau_2-\tau_1|/n),
\label{tightness-c1}
\end{align}
and $(ii)$ uniformly for $a\le \tau_1<\tau_2 \le b$,
\begin{equation}
|\mathcal{T}_n(\tau_2)-T_n(\tau_1)| \le K\times L(n)(n\log\log n)^{1/2}|\tau_2-\tau_1|+o_p(1)
\label{tightness-c2}
\end{equation}
as $n\to\infty$ where the $o_p(1)$ term holds uniformly, $K$ is a constant that may depend on $\theta$, and $L(\cdot)$ is some slowly varying function, i.e., for any $\lambda>0, L(\lambda x)/L(x)\to 1$ as $x\to\infty$. Then as $n\to\infty$, $T_n=\sup \{T_{n}(r), r \in [n^{1/2}(1-\theta)\sigma r_L, n^{1/2}(1-\theta)\sigma r_U]\}$ converges in distribution to
\begin{equation}
F(W; r_L, r_U)=\sup_{\tau\in [r_L, r_U]} \left\|\left[\{\Lambda^{-1}(\tau)\}_{2,2}\right]^{1/2} \left\{ H_{2}(\tau)-\Lambda_{2,1}(\tau)H_{1}(\tau)\right\}\right\|^2, \label{limiting-dist}
\end{equation}
whose distribution is parameter-free, although it depends on the search range of the threshold.
\end{theorem}
We remark that the assumption of independence and identical distribution of the innovations in the preceding theorem can be relaxed to $\{\varepsilon_t\}$ being a stationary, ergodic, martingale difference sequence with respect to the $\sigma$-algebra $\mathcal{F}_t$ generated by $\varepsilon_{t-s}, s\le 0$; the proof is essentially the same.
\begin{remark}
Conditions (\ref{tightness-c1})--(\ref{tightness-c2}) provide a new set of general sufficient conditions for the tightness of a sequence of stochastic processes; specifically the tightness of $\{T_{n}(n^{1/2}(1-\theta)\tau), r_L\le \tau \le r_U\}$. These sufficient conditions are motivated by the approach taken by \cite{Bil68}, Theorem 22.1, for studying the tightness of empirical processes for stationary, mixing data, and are tailor made for coping with nonstationarity under the null. To the best of our knowledge, this is the first rigorous proof of tightness for testing threshold nonlinearity against difference stationarity and constitutes a general theoretical framework that can be used in different settings.
\end{remark}
\noindent
The preceding theorem assumes deterministic threshold search interval. It can be readily extended to the case that the end points are fixed quantiles of the data, which are realization specific. We omit the proof as it is based on routine analysis that builds on Theorem~\ref{thm1} and the facts that for any fixed $0<p<1$, $(i)$ the $p$-quantile of $\{W_s, 0\le s \le 1\}$ is $O_p(1)$, which follow from \cite{Bjo19}, Proposition 3.2, and the Markov inequality, and $(ii)$ the $p$-quantile of $\{X_t, t=0,\ldots, n\}$ is asymptotically equal to its counterpart of $\{W_s, 0\le s \le 1\}$ times $n^{1/2}(1-\theta) \sigma$.
The following result shows that Theorem~\ref{thm1} holds for normally distributed innovations.
\begin{theorem}~\label{thm2}
Conditions (\ref{tightness-c1}) and (\ref{tightness-c2}) hold if $(i)$ $|\theta|<1$ and $(ii)$ $\{\varepsilon_t\}$ are independent and identically normally distributed with zero mean and finite positive variance.
\end{theorem}
Since the null distribution of $T_n$ is asymptotically similar, its quantiles can be derived numerically. The tabulated quantiles of the null distribution for different threshold ranges can be found in Section~\ref{SM:tab} of the Supplementary Material.
\section{Local Power}\label{sec:LP}
In this section we derive the asymptotic distribution of the supLM statistic under a sequence of local threshold alternatives and prove its consistency in having power approaching 1 with increasing departure in some direction from the null hypothesis. The mathematical framework is as follows. For each positive integer $n$, the system of hypothesis is:
\vskip .1 in
\noindent $H_{0,n}$: $(X_0,\ldots,X_n)$ follow the IMA(1,1) model: $X_t=X_{t-1}+\varepsilon_t-\theta\varepsilon_{t-1}.$\\
\noindent $H_{1,n}$: $(X_0,\ldots,X_n)$ follow the TARMA(1,1) model:
\begin{equation}
X_t=\begin{cases}
{n^{-1/2}h_{1,0}}{}+\left(1+{n^{-1}h_{1,1}}\right)X_{t-1}+\varepsilon_t-\theta\varepsilon_{t-1} & \text{if $\frac{X_{t-1}}{\sigma n^{1/2}(1-\theta)} \leq\tau_0$} \\
n^{-1/2}h_{2,0}+\left(1+{n^{-1}h_{2,1}}\right)X_{t-1}+\varepsilon_t-\theta\varepsilon_{t-1} & \text{if $\frac{X_{t-1}}{\sigma n^{1/2}(1-\theta)}>\tau_0$,}
\end{cases} \label{TARMAalt}
\end{equation}
\noindent
where $\mathbf{h}=(h_{1,0},h_{2,0},h_{1,1},h_{2,1})^{\intercal}$ is a fixed vector with $h_{i,1}\le 0, i=1,2$ and $\tau_0$ is a fixed threshold. Note that if $h_{1,1}<0$ ($h_{2,1}<0$), then the model is locally stable in the lower (upper) regime, for all $n$ sufficiently large. In order to derive the local power, we henceforth impose the following mild regularity conditions:
\begin{itemize}
\item[C1:] The innovations are assumed to be independent and identically distributed, with zero mean, finite positive standard deviation, $\sigma$, and probability density function
$f(\cdot/\sigma)/\sigma$, where $(i)$ $f$ is a bounded function, $\log(f(x))$ is twice differentiable with Lipschitz continuous first and second derivatives over the support of the probability density function, $(ii)$ the moment generating function of the innovations exists and is finite over some open interval around 0, and $(iii)$ $\mathcal{I}_f=-\int ({\ddot{f}f-\dot{f}^2}/{f^2})(x)\times f(x) dx$ is a finite positive number, where the first (second) derivative of $f$ is denoted by $\dot{f}$ ($\ddot{f}$).
\item[C2:] $-{\pi}/{2}<h_{1,1},h_{2,1}\le 0$ and $h_{1,1}+ h_{2,1} <0.$
\end{itemize}
Note that $\mathcal{I}_f$ is the Fisher information for the location model $f(\cdot-\mu)$ where $\mu$ is the location parameter.
Let $P_{0,n}$ and $P_{1,n}$ be the probability measures induced by $(X_0,\ldots,X_n)$ under $H_{0,n}$ and $H_{1,n}$, respectively.
Condition (C1) holds for many commonly used innovation distributions.
Condition (C2) ensures that the local alternative first-order TARMA model is asymptotically locally stable in at least one regime. These two conditions are imposed to ensure that $\{P_{1,n}\}$ is contiguous to $\{P_{0,n}\}$.
Finally, let $\rho$ be the correlation between $\varepsilon_t$ and $(\dot{f}/{f})(\varepsilon_t)$, i.e., $\rho \surd{\mathcal{I}_f} =\int x \dot{f}(x)dx$.
\begin{theorem}\label{thm:3}
Suppose all the conditions stated in Theorem~\ref{thm1} hold. Assume (C1) and (C2) hold. Under $H_{1,n}$ and as $n\to\infty$, $T_n=\sup \{T_{n}(r), r \in [n^{1/2}(1-\theta)\sigma r_L, n^{1/2}(1-\theta)\sigma r_U]\}$, where $r_L, r_U$ are two fixed numbers, converges in distribution to $F(W; r_L, r_U)$ defined in Eq.~(\ref{limiting-dist}) but with $W$ now being a threshold diffusion process satisfying the following stochastic differential equation (SDE):
\begin{equation}
dW_s=dW_s^\dagger+
\begin{cases} \rho\surd{\mathcal{I}_f}\left\{\frac{h_{1,0}}{\sigma(1-\theta)}+h_{1,1}W_s\right\} ds, & \text{if}\,\,\, W_s\leq\tau_0, \\
\rho\surd{\mathcal{I}_f}\left\{\frac{h_{2,0}}{\sigma(1-\theta)}+h_{2,1}W_s\right\} ds, & \text{otherwise},
\end{cases} \label{TDif}
\end{equation}
where $W_0=0$ almost surely and $\{dW_s^\dagger, s\ge 0\}$ is a standard Brownian motion.
\end{theorem}
\noindent
Henceforth in this section, $W$ denotes the threshold diffusion satisfying Eq.~(\ref{TDif}). Note that if $h_{i,0}=h_{i,1}=0, i=1,2$, then we get back the limiting null distribution for $T_n$. Otherwise, $W$ is a threshold diffusion process \citep{Su15}. Thus, the building block $W$ determining the limiting distribution of the supLM statistic changes from a standard Brownian motion under $H_{0,n}$ to a threshold diffusion under $H_{1,n}$, if $\rho\not =0$. Consequently, the proposed test would have power to detect the local threshold alternatives. Since the functional $F(\cdot; r_L, r_U)$ is quite complex, we examine an example to demonstrate the consistency of the proposed test.
\subsection{An asymptotically ergodic, symmetric TARMA(1,1) alternative.}\label{example}
Suppose the parameters in Eq.~(\ref{TARMAalt}) are such that $\tau_0=0$, $\rho\surd{\mathcal{I}_f}h_{1,1}= \rho\surd{\mathcal{I}_f}h_{2,1}=-1/2$ and
$\rho\surd{\mathcal{I}_f}\frac{h_{1,0}}{\sigma (1-\theta)}= 2h>0$ and
$h_{2,0}=-h_{1,0}$.
Consequently, $W$ is driven by the following SDE:
\begin{equation}
dW_s=\{2hI(W_s\le 0)-2hI(W_s>0)-W_s/2\}ds+ dW_s^\dagger,
\label{STDif}
\end{equation}
with initial condition $W_0=0$ and parameter $h>0$. It is an ergodic diffusion whose stationary marginal probability density function is given by
$$
\pi(x)=k^{-1}\exp[ -\{(x-h)^2I(x\le 0)+(x+h)^2I(x>0)\}/2 ], \quad -\infty<x<\infty,
$$
where
$$
k=2\sqrt{2\pi}\Phi(-h),
$$
and
$\Phi(\cdot)$ is the standard normal cumulative distribution function \citep[see][Theorem 1]{Su15}. It
is well known that
$$
\frac{1}{h}\phi(h)\le \Phi(-h) \le \frac{h}{h^2+1} \phi(h),
$$
where $\phi(\cdot)$ is the standard normal probability density function. Thus, the stationary probability density function equals $\pi(x)=(h/2)\{1+o(1)\} \exp(-x^2/2)\exp(-|x|h)$, which is asymptotically the Laplace distribution with scale parameter $h^{-1}$ (mean 0 and variance $2/h^2$), as $h\to\infty$.
As $h\to\infty$, the $W$ process approaches
stationarity at an increasing rate, implying that solving Eq.~(\ref{STDif}) via discretization can be achieved with step size inversely proportional to $h$. ($\mathcal{W}_h$ below stands for a random variable having the Laplace distribution with scale parameter $h^{-1}$.)
Consequently, it follows from ergodicity and the preceding discussion that, letting $\mathcal{D}_h$ be a $3\times 3$ diagonal matrix with $1,1,h$ as its diagonal,
$$
\mathcal{D}_h\Lambda_0 \mathcal{D}_h \to
\begin{pmatrix}
1 &\frac{1}{2} & -\frac{1}{2}\\
\frac{1}{2} & \frac{1}{2} & -\frac{1}{2} \\
-\frac{1}{2} & -\frac{1}{2} & 1
\end{pmatrix},
$$
a non-singular matrix; this is because, as $h\to \infty$,
$(i)$ $\int_0^1 I(W_s\le 0)ds\approx E\{I(\mathcal{W}_h\le 0)\}=1/2$, $(ii)$ $h\int_0^1 W_sI(W_s\le 0)ds\approx h E\{W_hI(\mathcal{W}_h\le 0)\}=-1/2$ and $(iii)$ $h^2 \int_0^1 W_s^2I(W_s\le 0)ds\approx h^2 E\{\mathcal{W}_h^2I(\mathcal{W}_h\le 0)\}=1$.
It can be similarly checked that
\begin{align}
\mathcal{D}_h H(0) &=
\begin{pmatrix}
\int_0^1 dW_s^\dagger +\int_0^1\{ -2h \times \operatorname{sign}(W_s)-W_s/2\}ds\\
\int_0^1 I(W_s\le 0)dW_s^\dagger + \int_0^1 (2h-W_s/2)I(W_s\le 0)ds\\
h\int_0^1 W_s I(W_s\le 0)dW_s^\dagger + h\int_0^1 (2h-W_s/2)W_sI(W_s\le 0)ds
\end{pmatrix}\\
&= (O_p(1), 2h+O_p(1) , -h+O_p(1))^{\intercal}. \nonumber
\end{align}
Assuming that $r_L\le 0\le r_U$, then
$$F(W; r_L,r_U)\ge \left\|\{(\Lambda_0^{-1})_{2,2}\}^{1/2} \left\{ H_{2}(0)-\Lambda_{2,1,0}H_{1}(0)\right\}\right\|^2\to\infty,$$
in probability as $h\to\infty$. Thus the supLM test statistic has power approaching 1 in rejecting the null hypothesis, as $h\to \infty$.
\section{Finite sample performance} \label{sec:sim}
In this section we compare the finite sample performance of the proposed test with a number of existing relevant tests. To better approximate the finite sample distribution of $T_n$, we have simulated the null distributions for the sample sizes in use. Moreover, since we have found that the finite sample distribution of $T_n$ changes appreciably only when $|\theta|$ is close to one, we have adopted the following, conservative, approach: if $|\hat\theta|>0.3$, we use the quantiles of the simulated null with $\theta = \operatorname{sign}(\hat\theta)\cdot0.9$. Furthermore, wild bootstrap obtained by bootstrapping with randomly signed residuals is also added to improve the empirical size of the test. We denote our asymptotic test and its wild bootstrap version by sLM and sLMb, respectively.
\par
The set of competing tests can be divided into those whose alternative is a threshold autoregressive model and those that do not specify explicitly a nonlinear alternative. The former set includes the tests proposed by \cite{Kap06} (KS), \cite{End98} (EG) and \cite{Bec04} (BBC), with their bootstrap variants (if implemented) signified by appending the extension b to their abbreviations. The latter set includes the ADF test of \cite{Dic79}, the class of M tests of \cite{Ng01} ($\bar{\mathrm{M}}^{\mathrm{g}}$), the $\bar{\mathrm{MP}}_{\mathrm{T}}^{\mathrm{GLS}}$ test of \cite{Ng01} (${\mathrm{MP}}_{\mathrm{T}}$) and the GLS detrended version of the ADF test (ADF$^{\mathrm{g}}$), and the test $\mathrm{M}^{\mathrm{GLS}}$ of \cite{Per07} ($\mathrm{M}^{\mathrm{g}}$). Note that we have obtained the results for all the tests proposed in the above references. Of these tests, we report only the best performing ones. Unreported results are available upon request.
\par
The sample sizes considered are 100, 300 and 500. The rejection percentages are derived with a nominal size $\alpha = 5\%$ and based upon 10000 replications. In order to reduce the computational burden, for the bootstrap tests we select 1000 replications and $B=1000$ bootstrap resamples. The threshold search ranges from the 25\% to the 75\% of the sample distribution.
We simulate data from the following first-order TARMA model
\begin{equation}
X_t=
\begin{cases}
\phi_{1,0}+\phi_{1,1}X_{t-1} +\varepsilon_t- \theta \varepsilon_{t-1}, & \text{if } X_{t-1} \le 0, \\
\phi_{2,0}+\phi_{2,1}X_{t-1} +\varepsilon_t- \theta \varepsilon_{t-1}, & \text{otherwise},
\end{cases}
\label{tarma.sim}
\end{equation}
where $(\phi_{1,0}, \phi_{1,1}, \phi_{2,0}, \phi_{2,1})=\tau\times (0, 0.7,-0.02,0.99)+(1-\tau)\times(0,1,0,1)$ with $\tau$ increasing from 0 to 1.5 with increments 0.5. When $\tau=0$, the model is an IMA(1,1) model with zero intercept. When $\tau>0$, the model becomes a stationary first-order TARMA model that is increasingly distant from the IMA(1,1) model with increasing $\tau$. As for the MA parameter we set $\theta=-0.9,-0.5,0,0.5,0.9$. The empirical sizes of the tests are displayed in Table~\ref{tabT:size}. Note that we have partitioned the set of 11 tests according to their nature: the first 9 are asymptotic and the last 2 are bootstrap tests. Clearly, the ADF, the KS, the BBC and the EG tests are severely oversized as $\theta$ approaches unity. Moreover, the wild bootstrap sLMb test is the only test that shows a correct size in all the settings, whereas both the sLM and the M class of tests show some bias, albeit small. Note that, when $\theta=0$ the TARMA model reduces to a TAR model. In this case, the auxiliary model of the KS, BBC, EG tests is correctly specified and their size is correct; however, when $\theta$ becomes positive their size is severely biased and this raises issues concerning their practical utility.
\begin{table}
\spacingset{1}
\small
\centering
\caption{Rejection percentages from the TARMA model of Eq.(\ref{tarma.sim}), with nominal size at $\alpha = 5\%$. Sizes over 15\% are highlighted in bold font.}\label{tabT:size}
\begin{tabular}{rrrrrrrrrrrr}
& \multicolumn{9}{c}{asymptotic} & \multicolumn{2}{c}{bootstrap}\\
\cmidrule(lr){2-10} \cmidrule(lr){11-12}
$\theta$ & sLM & $\bar{\mathrm{M}}^{\mathrm{g}}$ & $\mathrm{M}^{\mathrm{g}}$ & $\mathrm{MP}_\mathrm{T}$ & ADF & ADF$^{\mathrm{g}}$ & KS & BBC & EG & sLMb & KSb \\
\cmidrule(lr){2-10} \cmidrule(lr){11-12}
\multicolumn{10}{l}{$n=100$} & \multicolumn{2}{c}{}\\
\cmidrule(lr){1-1}
-0.9 & 2.2 & 7.7 & 7.0 & 7.1 & 2.5 & 3.8 & 8.1 & 11.2 & 7.1 & 5.1 & 4.9 \\
-0.5 & 1.6 & 6.3 & 6.1 & 5.8 & 4.8 & 5.1 & 7.0 & 6.1 & 5.7 & 5.0 & 5.6 \\
0.0 & 1.6 & 5.1 & 5.1 & 4.6 & 5.3 & 5.6 & 8.1 & 2.7 & 5.0 & 4.5 & 5.3 \\
0.5 & 1.7 & 5.6 & 5.9 & 5.1 & 6.7 & 7.4 &{\bfseries 64.5} & 10.2 &{\bfseries 57.5} & 5.2 &{\bfseries 58.6} \\
0.9 & 11.3 & 6.5 &{\bfseries 17.7} & 6.4 &{\bfseries 77.9} &{\bfseries 17.8} &{\bfseries 100.0} &{\bfseries 92.4} &{\bfseries 100.0} & 5.7 &{\bfseries 99.8} \\
\multicolumn{10}{l}{$n=300$} & \multicolumn{2}{c}{}\\
\cmidrule(lr){1-1}
-0.9 & 5.5 & 6.7 & 6.3 & 6.1 & 3.3 & 4.2 & 6.3 & 14.0 & 6.5 & 5.3 & 3.8 \\
-0.5 & 4.7 & 5.2 & 5.1 & 4.8 & 4.5 & 4.5 & 5.1 & 8.5 & 5.4 & 5.0 & 4.5 \\
0.0 & 2.9 & 4.9 & 4.9 & 4.4 & 5.1 & 4.6 & 6.9 & 3.2 & 4.4 & 5.6 & 4.3 \\
0.5 & 2.3 & 5.5 & 5.4 & 5.1 & 5.4 & 5.8 &{\bfseries 74.5} &{\bfseries 19.0} &{\bfseries 61.1} & 4.9 &{\bfseries 67.7} \\
0.9 & 4.9 & 1.9 & 2.4 & 1.9 &{\bfseries 86.0} &{\bfseries 15.8} &{\bfseries 100.0} &{\bfseries 99.7} &{\bfseries 100.0} & 4.9 &{\bfseries 100.0} \\
\multicolumn{10}{l}{$n=500$} & \multicolumn{2}{c}{}\\
\cmidrule(lr){1-1}
-0.9 & 8.1 & 6.4 & 6.1 & 6.0 & 7.4 & 4.7 & 5.7 &{\bfseries 16.0} & 6.1 & 5.5 & 4.0 \\
-0.5 & 5.3 & 5.5 & 5.3 & 5.0 & 5.1 & 4.8 & 5.2 & 9.2 & 5.4 & 4.7 & 4.2 \\
0.0 & 3.5 & 4.9 & 4.8 & 4.5 & 4.9 & 4.6 & 7.3 & 3.5 & 5.0 & 3.8 & 4.5 \\
0.5 & 2.5 & 5.2 & 5.1 & 4.8 & 5.1 & 5.3 &{\bfseries 78.4} &{\bfseries 23.7} &{\bfseries 62.3} & 4.5 &{\bfseries 71.7} \\
0.9 & 3.3 & 1.3 & 1.4 & 1.4 &{\bfseries 83.2} & 14.5 &{\bfseries 100.0} &{\bfseries 99.9} &{\bfseries 100.0} & 5.4 &{\bfseries 100.0} \\
\cmidrule(lr){2-10} \cmidrule(lr){11-12}
\end{tabular}
\end{table}
\begin{table}
\spacingset{1}
\centering \small
\caption{Size corrected power of the asymptotic and bootstrap tests at nominal size $\alpha=5\%$
}\label{tabT:SCpower}
\begin{tabular}{crrrrrrrrrrrr}
$n=300$ & \multicolumn{9}{c}{asymptotic} & \multicolumn{2}{c}{bootstrap}\\
\cmidrule(lr){1-1} \cmidrule(lr){2-10} \cmidrule(lr){11-12}
$\tau\;;\;\theta$ & sLM & $\bar{\mathrm{M}}^{\mathrm{g}}$ & $\mathrm{M}^{\mathrm{g}}$ & $\mathrm{MP}_\mathrm{T}$ & ADF & ADF$^{\mathrm{g}}$ & KS & BBC & EG & sLMb & KSb \\
\cmidrule(lr){2-10} \cmidrule(lr){11-12}
0.0;-0.9 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 \\
0.5;-0.9 & 25.7 & 17.6 & 17.8 & 18.2 & 10.2 & 19.4 & 5.1 & 16.5 & 1.6 & 23.7 & 8.3 \\
1.0;-0.9 & 52.5 & 26.7 & 26.9 & 27.7 & 15.8 & 30.4 & 17.4 & 31.9 & 3.8 & 54.3 & 27.0 \\
1.5;-0.9 & 77.1 & 33.5 & 34.0 & 35.1 & 22.1 & 38.2 & 36.9 & 50.1 & 8.2 & 75.6 & 45.5 \\
\cmidrule(lr){2-10} \cmidrule(lr){11-12}
0.0;-0.5 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.1 \\
0.5;-0.5 & 21.7 & 22.8 & 22.8 & 22.7 & 11.6 & 22.4 & 11.2 & 15.6 & 3.3 & 23.5 & 9.1 \\
1.0;-0.5 & 48.3 & 34.5 & 34.9 & 34.9 & 18.2 & 35.0 & 32.3 & 31.1 & 8.2 & 47.8 & 29.1 \\
1.5;-0.5 & 72.6 & 45.0 & 45.8 & 46.1 & 25.9 & 45.7 & 55.5 & 50.1 & 17.1 & 74.9 & 53.6 \\
\cmidrule(lr){2-10} \cmidrule(lr){11-12}
0.0;0.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.1 & 5.1 \\
0.5;0.0 & 22.4 & 25.8 & 26.1 & 26.9 & 11.0 & 26.6 & 37.9 & 15.2 & 22.6 & 22.5 & 40.7 \\
1.0;0.0 & 50.5 & 41.3 & 42.0 & 41.7 & 18.0 & 42.0 & 66.7 & 33.6 & 43.5 & 46.9 & 69.7 \\
1.5;0.0 & 75.3 & 54.7 & 55.7 & 55.7 & 27.7 & 55.7 & 84.8 & 55.8 & 65.8 & 73.8 & 84.7 \\
\cmidrule(lr){2-10} \cmidrule(lr){11-12}
0.0;0.5 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.1 & 0.0 \\
0.5;0.5 & 20.6 & 25.0 & 25.1 & 24.5 & 12.9 & 25.6 & 42.9 & 18.8 & 35.3 & 21.8 & 0.0 \\
1.0;0.5 & 50.1 & 39.5 & 40.1 & 39.2 & 26.2 & 40.9 & 70.9 & 45.1 & 65.8 & 49.2 & 0.0 \\
1.5;0.5 & 76.9 & 49.8 & 51.9 & 49.9 & 43.1 & 53.1 & 88.9 & 72.5 & 88.0 & 77.3 & 0.0 \\
\cmidrule(lr){2-10} \cmidrule(lr){11-12}
0.0;0.9 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 5.0 & 4.9 & 0.0 \\
0.5;0.9 & 24.8 & 16.2 & 19.9 & 16.0 & 14.9 & 18.2 & 6.4 & 34.6 & 29.6 & 14.3 & 0.0 \\
1.0;0.9 & 62.8 & 22.9 & 34.5 & 22.5 & 32.8 & 26.5 & 12.4 & 63.6 & 52.4 & 36.1 & 0.0 \\
1.5;0.9 & 86.3 & 25.3 & 44.9 & 25.0 & 47.8 & 29.5 & 23.5 & 77.2 & 65.8 & 61.7 & 0.0 \\
\cmidrule(lr){2-10} \cmidrule(lr){11-12}
\end{tabular}
\end{table}
The size-corrected power of the tests is presented in Table~\ref{tabT:SCpower}. Here, the sample size is 300; see Section~\ref{SM:power} of the Supplementary Material for results for $n=100, 500$. The rows for $\tau=0$ correspond to the size and other rows give size-corrected power. The size correction for bootstrap tests is achieved by calibrating the $p$-values. In some cases, the corrected size deviates from the nominal 5\% due to discretization effects on the empirical distribution of bootstrap $p$-values. Clearly, the supLM tests are almost always more powerful than the other tests, especially as $\tau$ increases. For instance, when $\tau=1.5$ the sLM test has almost double the power of M tests in several instances. As mentioned before, the case $\theta=0$ (central panel) corresponds to a TAR model and this is one of two instances where the KS tests are slightly more powerful than the supLM tests. The power of the bootstrap version of the KS test is zero in three cases, due to its 100\% oversize. See the Supplementary Material for further simulation results.
\subsection{Measurement error and heteroskedasticity}
In this section we assess the effect of measurement error and heteroskedasticity on the behaviour of the tests. We simulate from the following IMA(1,1) model
\begin{equation}\label{eq:ima11}
X_t= X_{t-1} + \theta \varepsilon_{t-1}+ \varepsilon_{t},
\end{equation}
where $\theta=-0.9$ (model M1), -0.5 (model M2), 0.5 (model M3), 0.9 (model M4). We add measurement noise as follows
\begin{equation}\label{eq:mod2}
Y_t = X_t + \eta_t,
\end{equation}
where the measurement error $\eta_t\sim N(0,\sigma^2_\eta)$ is such that the signal to noise ratio SNR~$=\sigma^2_X/\sigma^2_\eta$ is equal to $\{+\infty, 50,10,5\}$. Here, $\sigma^2_X$ is the variance of $X_t$ computed by means of simulation. The case without noise (SNR~$=+\infty$) is taken as the benchmark. The empirical sizes (rejection percentages) for models M1--M4 are presented in Table~\ref{tab:mes300} for $n=300$ and the results for $n=100$ and 500 can be found in Section~\ref{SM:MCmerr} of the Supplementary Material. Clearly, the measurement noise has little effect upon the size of the supLM tests. On the contrary, the size bias of the tests KS, BBC and EG increases appreciably when $\theta$ is positive (Models M3--M4). Worst still the bias does not reduce when the sample size increases.
\begin{table}
\spacingset{1}
\centering
\caption{Empirical size (rejection percentage) at nominal $\alpha=5\%$ and $n=300$ for the IMA(1,1) models M1--M4 with increasing levels of measurement error.}\label{tab:mes300}
\begin{tabular}{rrrrrrrrrrrrr}
&& \multicolumn{9}{c}{asymptotic} & \multicolumn{2}{c}{bootstrap}\\
\cmidrule(lr){3-11}\cmidrule(lr){12-13}
& \textsc{snr} & sLM & $\bar{\mathrm{M}}^{\mathrm{g}}$ & $\mathrm{M}^{\mathrm{g}}$ & $\mathrm{MP}_\mathrm{T}$ & ADF & ADF$^{\mathrm{g}}$ & KS & BBC & EG & sLMb & KSb \\
\cmidrule(lr){3-11}\cmidrule(lr){12-13}
\multirow{4}{5pt}{M1}
& $\infty$ & 4.4 & 7.1 & 6.7 & 6.8 & 3.4 & 4.8 & 6.2 & 12.5 & 6.7 & 5.0 & 5.3 \\
& 50 & 3.6 & 6.4 & 6.4 & 5.4 & 4.8 & 4.9 & 5.8 & 10.4 & 6.5 & 3.8 & 4.7 \\
& 10 & 2.8 & 5.0 & 5.0 & 4.2 & 6.1 & 4.7 & 4.9 & 5.7 & 5.1 & 5.0 & 4.0 \\
& 5 & 5.5 & 5.3 & 5.0 & 4.8 & 5.2 & 4.8 & 4.9 & 3.1 & 3.1 & 5.3 & 3.8 \\
\cmidrule(lr){3-11}\cmidrule(lr){12-13}
\multirow{4}{5pt}{M2}
& $\infty$ & 4.0 & 6.4 & 6.2 & 5.6 & 5.6 & 5.4 & 4.2 & 5.5 & 5.6 & 5.2 & 3.4 \\
& 50 & 4.7 & 6.3 & 5.9 & 5.9 & 5.8 & 5.4 & 4.0 & 4.8 & 5.4 & 6.1 & 3.2 \\
& 10 & 5.9 & 6.3 & 6.1 & 5.1 & 6.6 & 5.3 & 3.6 & 4.1 & 4.6 & 6.4 & 2.3 \\
& 5 & 5.4 & 5.5 & 5.3 & 5.1 & 6.3 & 5.2 & 5.4 & 2.5 & 5.6 & 5.4 & 3.8 \\
\cmidrule(lr){3-11}\cmidrule(lr){12-13}
\multirow{4}{5pt}{M3}
& $\infty$ & 2.8 & 5.8 & 5.8 & 4.4 & 5.6 & 6.3 & 67.8 & 14.1 & 59.9 & 5.2 & 62.0 \\
& 50 & 3.4 & 5.6 & 5.7 & 4.2 & 5.7 & 5.9 & 68.2 & 15.7 & 60.9 & 5.0 & 62.7 \\
& 10 & 2.4 & 6.2 & 6.0 & 5.0 & 5.8 & 7.0 & 74.2 & 19.7 & 67.4 & 4.8 & 66.7 \\
& 5 & 2.5 & 5.6 & 5.5 & 4.2 & 5.2 & 6.7 & 84.3 & 28.4 & 77.6 & 5.3 & 76.7 \\
\cmidrule(lr){3-11}\cmidrule(lr){12-13}
\multirow{4}{5pt}{M4}
& $\infty$ & 6.1 & 1.2 & 2.1 & 0.9 & 86.6 & 15.4 & 100.0 & 98.5 & 100.0 & 4.3 & 99.5 \\
& 50 & 5.8 & 1.2 & 1.8 & 1.1 & 87.8 & 15.6 & 100.0 & 98.8 & 100.0 & 3.6 & 99.6 \\
& 10 & 4.2 & 2.5 & 2.7 & 1.4 & 89.8 & 17.1 & 100.0 & 99.3 & 100.0 & 2.2 & 99.9 \\
& 5 & 6.7 & 3.5 & 4.5 & 2.9 & 94.9 & 19.5 & 100.0 & 99.8 & 100.0 & 3.6 &100.0 \\
\cmidrule(lr){3-11}\cmidrule(lr){12-13}
\end{tabular}
\end{table}
Since daily financial time series are characterized by volatility and this is known to bias the size of unit root tests \citep[see e.g][and references therein]{Bos18}, we simulate from the following models
\begin{center}
\small
\begin{tabular}{ll}
M5. ARIh &
$
\Delta X_t=-0.6\Delta X_{t-1}+\varepsilon_{t} $, where $V(x_t) = (1.5)^2 V(x_s)$, \\
& $t = 1,\dots,[n/2]$ , $s = [n/2]+1,\dots,n$.\\
M6. IMA-GARCH &
$
\begin{array}{ll}
\Delta X_t = -0.6\varepsilon_{t-1} + \sqrt{h_t}\,\varepsilon_t, & \text{where } h_t = 0.05 + 0.30 \varepsilon^2_{t-1} + 0.65 h_{t-1} \\
\end{array}
$
\\
M7. ARI-GARCH &
$
\begin{array}{ll}
\Delta X_t = 0.3 \Delta X_{t-1} + \sqrt{h_t}\,\varepsilon_t, & \text{where } h_t = 0.05 + 0.30 \varepsilon^2_{t-1} + 0.65 h_{t-1} \\
\end{array}
$
\end{tabular}
\end{center}
Model M5 is an integrated linear AR(1) where the variance of the series in the first half of the sample is different from that of the second half. Models M6 and M7 are integrated MA and AR with GARCH innovations. Measurement errors are added as described previously. The results of Table~\ref{tab:mesh300} confirm that the supLM tests are well behaved in the presence of heteroskedasticity and measurement error, with the sLMb wild bootstrap test being more so. The KS, BBC, and EG tests are severely affected by the combined presence of heteroskedasticity and measurement error. Furthermore, the size bias gets worse with increasing sample size (see the Supplementary Material, Section~\ref{SM:MCmerr}). The class of M tests is also robust in this respect but they can display low power in a number of instances, especially when the DGP is nonlinear. See also \cite{Cha20} and Section~\ref{SM:MCmerr} of the Supplementary Material.
\begin{table}
\spacingset{1}
\centering
\caption{Empirical size (rejection percentage) at nominal $\alpha=5\%$ and $n=300$ for the heteroskedastic models M5--M7 with increasing levels of measurement error.}\label{tab:mesh300}
\begin{tabular}{rrrrrrrrrrrrr}
&& \multicolumn{9}{c}{asymptotic} & \multicolumn{2}{c}{bootstrap}\\
\cmidrule(lr){3-11}\cmidrule(lr){12-13}
& \textsc{snr} & sLM & $\bar{\mathrm{M}}^{\mathrm{g}}$ & $\mathrm{M}^{\mathrm{g}}$ & $\mathrm{MP}_\mathrm{T}$ & ADF & ADF$^{\mathrm{g}}$ & KS & BBC & EG & sLMb & KSb \\
\cmidrule(lr){3-11}\cmidrule(lr){12-13}
\multirow{4}{5pt}{M5}
& $\infty$& 2.6 & 4.9 & 4.9 & 3.9 & 7.5 & 5.4 & 65.8 & 4.7 & 62.2 & 3.7 & 57.3 \\
& 50 & 2.7 & 5.9 & 5.9 & 4.7 & 6.9 & 6.4 & 82.1 & 8.0 & 81.3 & 3.5 & 72.0 \\
& 10 & 3.5 & 6.1 & 6.0 & 4.9 & 8.7 & 7.0 & 98.9 & 30.9 & 100.0 & 5.6 & 92.3 \\
& 5 & 4.0 & 5.7 & 5.4 & 4.7 & 19.2 & 9.5 & 100.0 & 67.2 & 100.0 & 7.4 & 98.9 \\
\cmidrule(lr){3-11}\cmidrule(lr){12-13}
\multirow{4}{5pt}{M6}
& $\infty$& 7.7 & 4.8 & 4.6 & 4.1 & 9.8 & 5.8 & 82.0 & 48.4 & 82.3 & 6.2 & 76.4 \\
& 50 & 7.4 & 5.0 & 5.0 & 3.8 & 10.3 & 5.9 & 91.9 & 53.5 & 92.7 & 7.0 & 86.0 \\
& 10 & 5.7 & 4.9 & 4.7 & 3.2 & 15.1 & 7.0 & 99.6 & 70.0 & 100.0 & 6.6 & 95.5 \\
& 5 & 5.2 & 3.5 & 3.7 & 2.4 & 27.6 & 7.9 & 100.0 & 90.0 & 100.0 & 6.9 & 98.4 \\
\cmidrule(lr){3-11}\cmidrule(lr){12-13}
\multirow{4}{5pt}{M7}
& $\infty$&11.1 & 5.2 & 5.0 & 4.6 & 7.6 & 4.5 & 10.6 & 12.5 & 8.4 & 4.2 & 8.6 \\
& 50 & 9.2 & 5.6 & 5.2 & 5.1 & 7.2 & 4.6 & 36.8 & 15.5 & 19.5 & 8.4 & 30.3 \\
& 10 & 7.0 & 6.5 & 6.2 & 5.8 & 10.0 & 6.3 & 96.2 & 40.5 & 100.0 & 9.1 & 89.4 \\
& 5 & 5.1 & 5.4 & 5.5 & 4.4 & 23.2 & 7.5 & 100.0 & 71.1 & 100.0 & 7.4 & 97.3 \\
\cmidrule(lr){3-11}\cmidrule(lr){12-13}
\end{tabular}
\end{table}
\section{A real application: new light on the PPP hypothesis}\label{sec:real}
In this section we apply our supLM tests to the real exchange rates of the panel of 17 European countries of the eurozone, plus the series of the aggregated eurozone. The idea is to contribute to the widely debated issue of the purchasing power parity (PPP) and show that the threshold ARMA model can be a useful tool to this aim. As mentioned in the introduction, macroeconomic theory suggests that price gaps (measured in a common currency) for the same goods in different countries should rapidly disappear. However, existing unit-root tests have generally failed to reject the null hypothesis of a random walk, suggesting the existence of a strong persistence of the gaps. As also pointed out in \cite{Tay01} the divergence between macroeconomic theory and econometric findings can be ascribed to multiple reasons. First, the way economic data are produced or aggregated can result in severely biased inference.
This is also noted in \cite{Pel15} where the authors show that real exchange rates based on the consumer price index do not preserve the possible stationarity properties of the ratios. Also, deviations from the equilibrium can be produced by the presence of non-traded goods in the basket and price frictions can slow down the mean-reversion mechanism. A third possible reason is that the commonly presumed linear price dynamics is a mis-specification. Indeed, the presence of trading costs implies that the mechanisms governing price adjustments are nonlinear and threshold autoregressive models provide a solution to the problem by allowing a ``band of inaction'' random walk regime, where arbitrage does not occur, and other regimes where mean reversion takes place so that the model is globally stationary \citep[see][and references therein for further discussion]{Bec04}. For a review on how TAR models are used to analyse the exchange rates dynamics see also \cite{Han11}. Among other approaches, \cite{Bec08b} and \cite{Gou06} introduced the switching models to incorporate a random threshold that delineates the regimes where arbitrage takes place. A critical investigation on the practical usefulness of combining unit-root tests and other stationarity tests in the PPP debate was put forward by \cite{Can01a}.
\par
Measurement noise is ubiquitous and can be a further hindrance to the empirical verification of a threshold regulated mechanism governing the dynamics of real exchange rates. To the best of our knowledge, our attempt will represent the first time that measurement noise is incorporated in the PPP debate. By leveraging the threshold ARMA framework, we focus on the daily $\log$ real exchange rates from the exclusive adoption of the Euro for the founding countries from \texttt{2002-01-01} to \texttt{2021-08-16} ($n=5120$). The countries are: Estonia (EE), Spain (ES), Luxembourg (LU), Ireland (IE), Italy (IT), Slovenia (SI), Portugal (PT), Netherland (NL), Malta (MT), Belgium (BE), Cyprus (CY), Slovakia (SK), Finland (FI), France (FR), Greece (GR), Austria (AT), Germany (DE). The series denoted with XM indicates the Eurozone as a whole. We also add to the panel the series for Great Britain (GB) and USA (US). Overall, the panel comprises 20 time series whose time plots are presented in Figure~\ref{fig:S1} of the Supplementary Material. The real exchange rates are produced by the Bank of International Settlements (BIS) by taking the geometric weighted average of a basket of bilateral exchange rates (60 economies), adjusted with the corresponding relative consumer prices. Such weights are constructed from manufacturing trade flows so as to encompass both third-market competition and direct bilateral trade through a double-weighting scheme. See \cite{Kla06} and \url{https://www.bis.org/} for more details on the construction of the indexes.
\par
Previous studies on the Eurozone focused on monthly data and failed to reveal convincing evidence for the PPP hypothesis. In particular, \cite{Bos18} suggested the presence of a non-stationary volatility over the period 1973-2015 and proposed an adaptive likelihood ratio test for the null hypothesis of a unit root against a linear heteroskedastic AR model. However, we do not find signs of non-stationary volatility over the period 2002-2021. As we will show, the key factors behind our results are the specification that includes a threshold regulation with measurement error and the long daily series.
\par
\begin{table}
\spacingset{1}
\centering \small
\caption{$p$-values of the supLM tests applied to the panel of 20 series. The $p$-values for the asymptotic test have been obtained from the quantiles of the simulated null distribution. The third row reports the sample percentile corresponding to the threshold.}\label{tab:ppp1}
\begin{tabular}{rrrrrrrrrrr}
& XM & US & EE & ES & LU & IE & IT & SI & PT & NL \\
\cmidrule(lr){2-11}
sLM & 0.001 & 0.272 & 0.381 & 0.076 & 0.001 & 0.001 & 0.052 & 0.001 & 0.056 & 0.001 \\
sLMb & 0.027 & 0.561 & 0.704 & 0.290 & 0.014 & 0.009 & 0.248 & 0.060 & 0.238 & 0.009 \\
\cmidrule(lr){2-11}
\end{tabular}
\begin{tabular}{rrrrrrrrrrr}
& MT & BE & CY & SK & FI & FR & GR & GB & AT & DE \\
\cmidrule(lr){2-11}
sLM & 0.001 & 0.001 & 0.001 & 0.001 & 0.167 & 0.001 & 0.122 & 0.001 & 0.001 & 0.001 \\
sLMb & 0.002 & 0.009 & 0.046 & 0.114 & 0.432 & 0.023 & 0.373 & 0.048 & 0.025 & 0.023 \\
\cmidrule(lr){2-11}
\end{tabular}
\end{table}
Table~\ref{tab:ppp1} reports the results of the application of the supLM tests on the 20 series. The threshold is searched from the 1st to the 99th percentiles of the sample and the $p$-values for the asymptotic test were obtained from the quantiles of the simulated null distribution. The asymptotic $p$-values indicate a clear rejection of the null hypothesis in 13 out of 20 series. The wild-bootstrap $p$-values are more conservative and account for the possible presence of heteroskedasticity. At level 5\%, they reject the null hypothesis for 11 series. Note that none of the M tests rejects the null hypothesis at the 5\% level (see Table~\ref{tab:S0} of the Supplement). In many situations when the supLM tests reject the null hypothesis the threshold is located in the extreme upper tail of the sample distribution; see also Figure~\ref{fig:S2} of the Supplementary Material for details. The first conclusion that can be drawn is that the mean reversion mechanism is present in most of the series but a large sample size is needed in order to detect it. Another important piece of information, presented in Figure~\ref{fig:1}, is brought by the plot of the values of the sLM statistic versus the threshold $r$. Clearly, two clusters of countries are identifiable. The first one includes those of northern Europe (plus Cyprus and Malta) for which the sLM statistic is well above the 99\%th critical value of the null distribution (red dashed line) and with a higher threshold (around 105). The second cluster identifies the Mediterranean countries plus Finland and Estonia and is characterized by lower thresholds and values of the statistic just below the 90\%th percentile of the null distribution (purple dashed line). This result can reflect the different sectorial composition of trade of the two groups: the northern countries export is mainly directed at medium to high-tech industries while the Mediterranean countries (plus Finland and Estonia) are more focused on low-tech export \citep[see e.g.][for further details]{Dim08}. Finally, we have tested for regulation from below by applying the test to $-X_t$. The results are presented in Figure~\ref{fig:S8} of the Supplementary Material and show no apparent threshold effects in the lower tail, the only exception being the Great Britain.
\begin{figure}
\centering
\includegraphics[width=0.6\linewidth]{fig_1}
\caption{Values of the sLM statistic (in log scale) versus threshold for the eurozone countries. The critical values of the null distribution at levels 90\%, 95\% and 99\% are added as purple, green and red dashed lines, respectively.}\label{fig:1}
\end{figure}
\par
We now turn to estimating a TARMA model on the global Eurozone series (XM). In Figure~\ref{fig:2} we present the time series (in natural scale) with the estimated threshold that minimizes the AIC criterion $\hat r=108.74$ indicated by a red dashed line. The threshold grid ranges from the 1st to the 99th percentiles of the data and the estimated threshold falls around the 96th percentile. It lends support to the suggestion that the mean reversion mechanism is triggered upon crossing a threshold located in the extreme upper tail of the sample distribution. For this reason, a large sample size is needed and in our case the upper regime contains 204 observations. The gray shaded area denotes the days associated with the upper regime.
\begin{figure}
\centering
\includegraphics[width=0.7\linewidth]{fig_2}
\caption{Time series of the real daily exchange rates for the Eurozone series (XM, in blue). The estimated threshold $\hat r=108.74$ that minimises the AIC criterion is indicated with a red dashed line and the gray shaded area indicates the days associated with the upper regime. The estimated volatility from Eq.~(\ref{garch.fit}) is also shown as a green line.}\label{fig:2}
\end{figure}
\begin{equation}
X_t=
\begin{cases}
\underset{(0.19)}{0.67} + \underset{(0.04)}{0.86} X_{t-1} - \underset{(0.07)}{0.41} \varepsilon_{t-1} +\varepsilon_t, & \text{if } X_{t-1} > 108.74\\
\underset{(0.0004)}{0.0135} + \underset{(0.0009)}{0.9971} X_{t-1} + \underset{(0.014)}{0.008} \varepsilon_{t-1} +\varepsilon_t, & \text{if } X_{t-1} \leq 108.74
\end{cases}\label{tarma.fit}
\end{equation}
The parameter estimates are presented in Equation~(\ref{tarma.fit}) (with standard errors in parenthesis). The fitted model, which is globally stationary and ergodic, has a near unit-root lower regime, and an upper TARMA(1,1) regime where the slope of the AR part is clearly smaller than 1. Figure~\ref{fig:2} shows that the ``intervention'' regime is visited mostly in 2008-2009 and between August 2020 and January 2021. A further confirmation of the goodness of the fit is shown in Figure~\ref{fig:3}, where we present the lag plot of $X_t$ vs $X_{t-1}$, zoomed in the proximity of the threshold, and with the 45$^\circ$ dashed line in blue and the nonparametric fit (loess) in continuous red line. Clearly, the nonparametric fit deviates from the 45$^\circ$ line starting from values close to the estimated threshold.
\par
We note that the MA parameter $\theta$ has played a crucial role in greatly enhancing the fitting capability of the model while retaining parsimony. Now we examine the TARMA model fitted above. While there is no autocorrelation in its residuals (see Figure~\ref{fig:S3} of the Supplement), the test based upon the entropy measure $S_k$, introduced in \cite{Gia15} rejects both the null hypothesis of independence and that of linearity of the residuals (see Figure~\ref{fig:S4} of the Supplement). Further, an inspection of the correlograms of the squared residuals (Figure~\ref{fig:S5} of the Supplement) and a test for ARCH effects suggest the possible presence of conditional heteroscedasticity. Hence, we improve the TARMA(1,1) model by fitting a GARCH(1,2) model to the residuals $\hat e_t$ and obtain the following:
\begin{align}\label{garch.fit}
\hat e_t &= \sqrt{h_t}z_t \nonumber\\
h_t &= \underset{(0.004)}{0.056} \hat e^2_{t-1} + \underset{(0.01)}{0.49} h_{t-1} + \underset{(0.01)}{0.45} h_{t-2}.
\end{align}
Here, the conditional distribution of the innovations is assumed to follow the Generalized Error distribution with estimated shape $1.33$ (with standard error 0.03), which accounts for the non-Gaussianity, as shown in Figure~\ref{fig:S6} of the Supplement. The entropy-based test for serial independence applied upon the standardized residuals of the fit shows no residual dependence of any kind (see Figure~\ref{fig:S7} of the Supplement), and this is also confirmed by the battery of diagnostic tests implemented in the R package \texttt{rugarch} \cite{Gha20} and applied to the fitted GARCH model. (The results are available upon request.) As a further confirmation of the results, we have computed both the supLM and the M tests on simulated series from the fitted models. In particular, the empirical sizes reported in Table~\ref{tab:Sfit1} of the Supplement have been obtained from an ARIMA-GARCH model with GED innovations whose parameters have been estimated on the Eurozone series. While the asymptotic sLM test shows some size bias, probably due to the presence of volatility and heavy tails, the wild bootstrap sLMb test has always a correct size and behaves like the M tests. Table~\ref{tab:Sfit2} of the Supplement reports the empirical power computed on simulated series from the fitted TARMA-GARCH process with GED innovations. (See Eq.~(\ref{tarma.fit}) and (\ref{garch.fit}).) Here, the asymptotic sLM test is the most powerful whereas the wild bootstrap test sLMb and M tests have similar power, the former being slightly more conservative so that a rejection can be taken as a genuine indication of a stationary DGP.
Finally, the estimated volatility from Eq.~(\ref{garch.fit}) is added as a green line to Figure~\ref{fig:2}. The inspection of the plot reveals that the series is characterized by consistent fluctuations of conditional variance but the threshold crossings are not trivially associated to periods of high volatility.
\begin{figure}
\centering
\includegraphics[width=0.4\linewidth]{fig_3}
\caption{Lag plot of $X_t$ vs $X_{t-1}$ for the daily time series of real exchange rates for the Eurozone (XM). The 45$^\circ$ line and the nonparametric fit (loess) are shown in the blue dashed line and the red continuous line, respectively. The estimated threshold $\hat r$ is indicated by the vertical red dashed line.}\label{fig:3}
\end{figure}
\section{Conclusion}\label{sec:conc}
In this paper, we argue that the ubiquity of measurement error implies that to test for regulation in dynamics, it is more appropriate and perhaps even crucially important to formulate the test within a TARMA specification. We adopt the TARMA(1,1) model as the general hypothesis and the IMA(1,1) model as the null hypothesis. As far as we know, this is the first time that a TARMA specification is used in the present context although it was previously utilized in a very different context, namely for linearity testing under stationarity \citep{Li11}. We derive a Lagrange multiplier test which is asymptotically similar given the threshold search range. Empirical studies confirm that the proposed approach enjoys much higher power in detecting regulation in dynamics than existing tests that do not address measurement errors. The IMA(1,1) model can mimic well other integrated processes with short-memory first differences. In particular, empirical results reported in \cite{Cha20} and in the Supplementary Material indicate that our new tests generally perform well under heteroskedasticity, even when the null hypothesis entails a non-stationary process different from the IMA(1,1) model, and remain powerful for other forms of regulation. The application to the real exchange rates of a panel of time series of the Eurozone sheds new light upon the PPP debate: the evidence points to a threshold-regulated mean reversion mechanism that takes place in the extreme upper tail of the distribution and, as such, a large sample size is needed before a definitive conclusion can be drawn. This could be the reason behind the failure of many approaches to detect the mean-reversion. Also, the deviations from the random walk regime tend to be compensated within days or weeks and this seems to be consistent with the general expectation. Moreover, the estimated thresholds identify two clusters of countries with different sectorial compositions of trade, thereby prompting further interesting macroeconomic investigations.
\medskip
\subsection*{SUPPLEMENTARY MATERIAL}
The Supplementary material (pdf format) contains all the proofs, further results from the real data analysis, the tabulated quantiles of the null distribution and further Monte Carlo investigations.
\bibliographystyle{plainnat}
\bibliography{Tarma16}
\newpage
\setcounter{table}{0}
\setcounter{section}{0}
\setcounter{page}{1}
\begin{center}