The Table Auto-Regressive Moving-Average Model for (Categorical) Stationary Series: Mathematical Perspectives (Invertibility; Maximum Likelihood Estimation)

Chrysoula Dimitriou-Fakalou

Independent Researcher, Athens, Greece.

DOI: 10.4236/ojs.2022.123025                

Abstract

Once invertibility for a causal TARMA series is defined and accompanied by conditions on the probability parameters of the model, all focus concentrates on the maximum likelihood estimators. Under the coexistence of essential causality and invertibility, the estimators are shown to be convergent to the real values and asymptotically obedient to the Gaussian distribution: their variance matrix identifies with a classic result. Some real-like examples are simulated and simplification attempts include the derivation of the non-parametric chi-square test extension for stationary TAR series.

Keywords

Asymptotic Normality, Consistency, Scores Efficiency, Strict Stationarity

1. Introduction

The scientific progress, from defining a valid time series model to making it useful in practice, depends on the consolidation of inference results. With regard to the non-linear stationary time series, the general TARMA model ( [1]) has recently pledged an unprecedented flexibility: given the mild requirement to categorize the variables’ values, its competence to express any serial order conditional or joint probability dependence is unquestionable. So it will be the aim of this paper to bond the theoretical with the most deserving TARMA sample properties. Before that, a scent of other results and popular models in the field is offered.

The DARMA model undoubtedly triggered a significant number of scientists’ interest for the discrete time series analysis: the reader may skim through [2] to get acquainted with the definition and some asymptotic properties. A recent fix to replication issues of the older model, can be discovered in [3]. Having to do with an additive model though, it should be crystal clear that its (and those pieces’ following it) simplicity and parsimony contribution to the discrete stationary series modelling, are due to its inarguable impairment to manage any better than the marginal and covariance dependence (which is not the case for the TARMA).

A valuable equivalence has been established in the past between the analysis of Bernoulli variables and the strictly stationary series obeying any law ( [4]): this saves the trouble of generalizing the Gaussian ARMA as in [5] for a family of infinitely divisible distributions. Nevertheless, most of the derivations concern the Gaussian stationarity, such as [6] who deals with the alternative (via 0 - 1) estimation of autoregressive parameterizations: those binary data have traditionally benefited from a rich bibliography anyway.

A promotion to the maximum likelihood estimation together with some solid statements can be found in [7], but the dense and short paper uses an approximation to the fixed Markovian or auto-regressive order; the inclusion of a moving-average part only presents itself with the Gaussian ARMA (1, 1) illustration, but it is exactly converting the AR to the ARMA that makes the problem challenging, ever more so when the distribution is not Gaussian. In the best estimators’ properties department, it is worth reading [8] to view the well-known ARMA efficiency result.

Notable attempts have taken place for count series as well; besides modelling as usual, [9] reproduced standard discrete-time renewal process inference results. In the other hand there are the INARMA models, for example, the INAR (1) (with random coefficient) asymptotic behaviour has been studied extensively by [10]. A recent review on the subject of count time series is [11]. One can always incorporate a count variable into a categorical one, encouraging the reader to elect the TARMA model as superior for the series stationarity. For bilinear models (and ‘any’-valued variables) there has been a plethora of inference results too.

Though a distinction wall is often raised between the time series and the Markov chain schemes, the causal TAR model embodies the homogeneous Markov dependence attached to a unique stationary distribution. Hence, [12] were the original contributors for the sake of the Markov chain inference. Later, [13] restricted the derivations to two states basic chain dependence managing explicit efficiency-related results, when (for two parameters only) the inverse of a square matrix can be computed with ease.

In this paper, the main objective is to perform the TARMA parameters inference, filling the gap of estimation for the infinite Markovian dependence: Section 2 summarizes the TARMA definition together with some raw statements concerning the use of the model parameters for a stationary series. Then Section 3 newly defines the invertible TARMA series and Section 3.1 invents a condition that safeguards invertibility, using a parallel with the existing definition and condition for causality. Section 3.2 continues with the invertibility topic, delivering news on the limit of probabilities based on a finite to the countably infinite past. Then all Section 4 is to establish the two grand properties of the maximum likelihood estimators: the weak convergence to the parameters is in Section 4.1, while Section 4.2 uses numerous lemmas and a theorem to conclude with the asymptotic normality, as this is explicitly stated in the beginning of Section 6; there, the normally distributed estimator vector transforms to a chi-square statistic to test a null hypothesis on the specified TARMA model fitting stationarity. The purely theoretical presentations are quickly reminded (Sections 2, 3), established (Sections 3, 4 and 6) or discussed (Section 6). In addition, the simulation Section 5 examines the performance of TARMA estimators (or how to compute the estimates) for small sample sizes, and returns favorable conclusions as well.

2. Reminders

Previously in [1], a strictly stationary time series $\left\{X_t\mathrm{,}t\in \mathbb{Z}\right\}$ was considered: the variables can be of any (bounded or unbounded) range, which for the serial stationarity clothing must be assigned to $\left(k+1\right)$ categories (for fixed $k\in \mathbb{N}$ ) with conventional (or, natural, if that is the case) category codes, say, 0 and $v_1\mathrm{,}\cdots \mathrm{,}{v}_k\ne 0$. The X-variables are built on a multivariate sequence of independent in time and identically distributed ${\left(k+1\right)}^{p+q}\times 1$ random vectors: for each $i=\left(i_1,\cdots ,i_p\right)$, $j=\left(j_1\mathrm{,}\cdots \mathrm{,}{j}_q\right)$, $i_1,\cdots ,i_p,j_1,\cdots ,j_q=0,v_1,\cdots ,v_k$, a univariable $I_t^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}$ ( $I_t^{\left(0_p\mathrm{<mo>\; |\; </mo>}{0}_q\right)}\equiv I_t$ ) with same range as X is considered at time t, and marginally it has to hold $\left\{I_t^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}\right\}~\text{IID}$ with probabilities

${\pi }_x^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}=\mathbb{P}\left(I_t^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}=x\right)\in \left(\mathrm{0,1}\right)\mathrm{,}x=\mathrm{0,}{v}_1\mathrm{,}\cdots \mathrm{,}{v}_k\mathrm{,}{\displaystyle \underset{x=\mathrm{0,}{v}_1\mathrm{,}\cdots \mathrm{,}{v}_k}{\sum }}{\pi }_x^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}=1$

(and ${\pi }_x^{\left(0_p\mathrm{<mo>\; |\; </mo>}{0}_q\right)}\equiv {\pi }_x$ ).

It is reminded that a multivariate IID time series is such that the random vectors are independent in time only, or that any two univariables indexed at different times are independent: the variables’ dependence at the same time within the same vector will be referred to as interdependence with interindependence being a special case; hence the ${\pi }_x^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}$ above are marginal probabilities that do not necessarily suffice to determine the joint dependence at the same timing.

In the general setting, $\left\{I_t^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}=x\right\}$ will simplify to the event that the variable $I_t^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}$ belongs into the category with code $x=0,v_1,\cdots ,v_k$. The Markov chain terminology of a state space, say $\mathcal{S}\mathrm{<mo>\; :\; </mo>}=\left\{\mathrm{0,}{v}_1\mathrm{,}\cdots \mathrm{,}{v}_k\right\}$ and $\cdots \in \mathcal{S}$ will be avoided deliberately. Then the definition introduced in [1] is reminded here.

For fixed $p\mathrm{,}q\in {\mathbb{N}}_0$ (such that $p+q\in \mathbb{N}$ ), $\left\{X_t\mathrm{,}t\in \mathbb{Z}\right\}$ is a TableAuto-Regressive Moving-Average process of order $\left(k\mathrm{,}p\mathrm{,}q\right)$, and it is written $\left\{X_t\right\}~\text{TARMA}\left(k\mathrm{,}p\mathrm{,}q\right)$, if (it is causal and invertible and) it holds that

$X_t\mathrm{<mo>\; :\; </mo>}=I_t^{\left(\left(X_{t-1}\cdots X_{t-p}\right)\mathrm{<mo>\; |\; </mo>}\left(I_{t-1}\cdots I_{t-q}\right)\right)}\mathrm{,}\text{for}\text{every}t\in \mathbb{Z}\mathrm{.}$ (1)

In explanation, $X_t$ is set equal to a $I_t^{\left(\mathrm{...}\right)}$ (at the original range); nevertheless, which I (out of the ${\left(k+1\right)}^{p+q}$ ) contributes its value, depends on the previous $X_{t-1},\cdots ,X_{t-p},I_{t-1},\cdots ,I_{t-q}$, i.e., which category each one of these $\left(p+q\right)$ variable realization falls into.

As an example of a TARMA (1, 1), the real-valued variables may be categorized into $k+1=3$ groups, say $\left(-\infty \mathrm{,}-3\right]$, $\left(-\mathrm{3,3}\right)$ and $\left[\mathrm{3,}\infty \right)$ with codes “-3”, “0” and “3”, respectively; then the (9 × 1) random vectors $I_t^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}$, $i,j=-3,0,3$ are considered with fixed distribution for every $t\in \mathbb{Z}$, and then over time with $I_t^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}$ being independent of $I_{t^{*}}^{\left(i^{*}|j^{*}\right)}$ for $t\ne t^{\mathrm{<mo>\; *\; </mo>}}$. The $X_t$ is defined to be equal to $I_t^{\left(-\mathrm{3\; <mo>\; |\; </mo>}-3\right)}$ if both $X_{t-1}$ and $I_{t-1}^{\left(\mathrm{0\; <mo>\; |\; </mo>0}\right)}$ are less than or equal to −3, or equal to $I_t^{\left(\mathrm{0\; <mo>\; |\; </mo>}-3\right)}$ if $-3<X_{t-1}<3$ and $I_{t-1}^{\left(\mathrm{0\; <mo>\; |\; </mo>0}\right)}\le -3$, or...., or equal to $I_t^{\left(\mathrm{3\; <mo>\; |\; </mo>3}\right)}$ if both $X_{t-1}$ and $I_{t-1}^{\left(\mathrm{0\; <mo>\; |\; </mo>0}\right)}$ are greater than or equal to 3.

Hopefully, the example above clarifies that it is not necessarily a count series that is under study. This is an ‘occurrence-not occurrence’ analysis of categorical variables and, for the serial evolution, it matters not whether $X_t$ is set equal to the exact value or the relevant category code: it is a probability of occurrence analysis and it works subject to the categorization (together with k of course).

Then (1) can be re-arranged as

$\begin{array}{l}{f}_x\left(X_t\right):={\displaystyle \underset{i_1,\cdots ,i_p,j_1,\cdots ,j_q=0,v_1,\cdots ,v_k}{\sum }}{f}_x\left(I_t^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}\right){\displaystyle \underset{l=1}{\overset{p}{\prod }}}{f}_{i_l}\left(X_{t-l}\right){\displaystyle \underset{n=1}{\overset{q}{\prod }}}{f}_{j_n}\left(I_{t-n}\right),\\ \text{for}x=0,v_1,\cdots ,v_k,\end{array}$ (2)

where for $x,y=0,v_1,\cdots ,v_k$, the index functions

$f_x\left(y\right)=\frac{{\displaystyle {\prod }_{x^{*}=0,v_1,\cdots ,v_k,x^{*}\ne x}\left(y-x^{*}\right)}}{{\displaystyle {\prod }_{x^{*}=0,v_1,\cdots ,v_k,x^{*}\ne x}\left(x-x^{*}\right)}}=\{\begin{array}{ll}1,\hfill & \text{if}y=x\hfill \\ 0,\hfill & \text{if}y=0,v_1,\cdots ,v_k,y\ne x\hfill \end{array}$

have been set.

For (2) to hold for every $t\in \mathbb{Z}$, at the very least the assumption of uniqueness of a stationary distribution solution should be applicable; in fact, $\left\{X_t\mathrm{,}t\in \mathbb{Z}\right\}$ should not be considered at all, unless it is a (table) causal process (based on $\left\{I_t\right\}$ ) according to the definition in [1]. The causality element will be revisited here together with the new element of invertibility.

Since the definition of this model has been the focus of past work, the reader is encouraged to look for examples and form a clearer picture from there ( [1]). The understanding and appreciation of the TARMA definition must have already taken place to proceed with this paper’s inference related goals. For the reader that wishes to derive the auto-covariance function of TARMA series, it is highlighted that the second moments stationarity is a special case of the all moments stationarity: the work of [1] has handed over a methodology to compute all joint probabilities, and hence the auto-covariance function is a special case, as it demands the joint probabilities of two variables. Not only are those derivations outside the scope of this paper, but also taming the all- rather than just second-order stationarity is the great TARMA contribution over other models.

From this point on, it will be taken for granted (and without a formal proof) that the TARMA construction may represent any strictly stationary time series, particularly under the convention- if there needs to be one- of tabulating the variables’ values. Together with a solution on identification (of k and p, q) issues, the TARMA theory is a contribution to the non-parametric version of statistical science. When $q=0$ the dependence is pth order Markovian, while for $q\ge 1$ a valid Markovian dependence of infinite order is achieved; when $p=0$, a q-dependent strictly stationary series is concerned, and the joint dependence limits to infinity for $p\ge 1$. Hence it will be considered ( [1]) that, subject to the categorization, the TARMA can be the all-moments stationary analogue to the ARMA model for the covariance stationary series.

Some Basic Requirements

The probabilities $\pi$ will be taking the role of parameters of the model, for which the estimation results will be established. The conditions in [1] are remembered:

(C1): $\left\{I_t^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}\mathrm{,}t\in \mathbb{Z}\right\}$ are jointly interindependent series, i.e., it holds that

$\mathbb{P}\left(I_t^{\left(i_1\mathrm{<mo>\; |\; </mo>}{j}_1\right)}=x_1\mathrm{,}\cdots \mathrm{,}{I}_t^{\left(i_n\mathrm{<mo>\; |\; </mo>}{j}_n\right)}=x_n\right)={\displaystyle \underset{m=1}{\overset{n}{\prod }}}\mathbb{P}\left(I_t^{\left(i_m\mathrm{<mo>\; |\; </mo>}{j}_m\right)}=x_m\right)\mathrm{,}t\in \mathbb{Z}$

for any $n=\mathrm{2,}\cdots \mathrm{,}{\left(k+1\right)}^{p+q}$, $\left(i_m\mathrm{,}{j}_m\right)=\left(i_{m\mathrm{,1}}\mathrm{,}\cdots \mathrm{,}{i}_{m\mathrm{,}p}\mathrm{,}{j}_{m\mathrm{,1}}\mathrm{,}\cdots \mathrm{,}{j}_{m\mathrm{,}q}\right)$, $i_{m\mathrm{,1}}\mathrm{,}\cdots \mathrm{,}{i}_{m\mathrm{,}p},j_{m\mathrm{,1}}\mathrm{,}\cdots \mathrm{,}{j}_{m\mathrm{,}q}=\mathrm{0,}{v}_1\mathrm{,}\cdots \mathrm{,}{v}_k$, $m=1,\cdots ,n$, $\left(i_{m_1}\mathrm{,}{j}_{m_1}\right)\ne \left(i_{m_2}\mathrm{,}{j}_{m_2}\right)$, $m_1,m_2=1,\cdots ,n$ ( $m_1\ne m_2$ ) and $x_1\mathrm{,}\cdots \mathrm{,}{x}_n=\mathrm{0,}{v}_1\mathrm{,}\cdots \mathrm{,}{v}_k$.

Condition (C1) enforces the $I_t^{\left(\mathrm{...}\right)}$ variables (at the same time t) to be independent: otherwise, the form of their interdependence should accompany (1), in order to properly define a TARMA process. Under an interdependence scenario, the computation of all joint TARMA probabilities heavily relies on

${\pi }_{x\mathrm{<mo>\; |\; </mo>}y}^{\mathrm{<mo>\; *\; </mo>}\left(i\mathrm{<mo>\; |\; </mo>}j\right)}\mathrm{<mo>\; :\; </mo>}=\mathbb{P}\left(I_t^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}=x\mathrm{<mo>\; |\; </mo>}{I}_t=y\right)\mathrm{,}x\mathrm{,}y=\mathrm{0,}{v}_1\mathrm{,}\cdots \mathrm{,}{v}_k\mathrm{,}\left(i\mathrm{,}j\right)\ne 0_{p+q}$

(as well as $\pi$ ) according to [1]. Hence in this case, ${\pi }^{\mathrm{<mo>\; *\; </mo>}}$ are also parameters for an apposite TARMA definition.

In the absence of (C1), an “h to h + 1” way to plex the interdependence may work by first considering ${\pi }_x$ and then gradually building for $n=1,\cdots ,p+q$ and every $l_1\mathrm{,}\cdots \mathrm{,}{l}_n=\mathrm{1,}\cdots \mathrm{,}p+q$ ( $l_1<\cdots <l_n$ ), $i_{l_1}\mathrm{,}\cdots ,i_{l_n}=v_1\mathrm{,}\cdots \mathrm{,}{v}_k$, the conditional probabilities $\mathbb{P}\left(I_t^{\left(0_{l_1-1}{i}_{l_1}\cdots i_{l_n}{0}_{p+q-l_n}\right)}=x_n\mathrm{<mo>\; |\; </mo>}{I}_t=x\mathrm{,}\cdots \right)$ : “ $\cdots$ ” in the condition sets all $I_t^{\left(0_{l_1^{\mathrm{<mo>\; *\; </mo>}}-1}{i}_{l_1^{\mathrm{<mo>\; *\; </mo>}}}\cdots i_{l_{n^{\mathrm{<mo>\; *\; </mo>}}}^{\mathrm{<mo>\; *\; </mo>}}}{0}_{p+q-l_{n^{\mathrm{<mo>\; *\; </mo>}}}^{\mathrm{<mo>\; *\; </mo>}}}\right)}$, $n^{\mathrm{<mo>\; *\; </mo>}}=\mathrm{1,}\cdots \mathrm{,}n-1$, $l_1^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\cdots \mathrm{,}{l}_{n^{\mathrm{<mo>\; *\; </mo>}}}^{\mathrm{<mo>\; *\; </mo>}}=l_1,\cdots \mathrm{,}{l}_n$ ( $l_1^{*}<\cdots <l_{n^{*}}^{*}$ ). That way one gets a feel of how (C1) can be replaced, especially when other attributes need to be accomplished.

Other conditions in [1] are listed below:

(C2): $\left\{X_t\mathrm{,}t\in \mathbb{Z}\right\}$ is not a deterministic process, i.e.,

$\mathbb{P}\left(X_t=x|X_{t-1}=x_1,X_{t-2}=x_2,\cdots \right)\in \left(0,1\right)$

for all $x,x_n=0,v_1,\cdots ,v_k$, $n\in \mathbb{N}$.

(C3): $\left\{X_t\mathrm{,}t\in \mathbb{Z}\right\}$ is not overparameterized, i.e., in the simplest of cases, it cannot be that ${\pi }_x^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}={\pi }_x^{\left(i^{\prime }\mathrm{<mo>\; |\; </mo>}{j}^{\prime }\right)}$, $x=\mathrm{0,}{v}_1\mathrm{,}\cdots \mathrm{,}{v}_k$ for $\left(i\mathrm{,}j\right)\ne \left(i^{\prime },j^{\prime }\right)$.

It would be misleading that the process $\left\{X_t\right\}$ is of order $\left(p\mathrm{,}q\right)$ with $\left(k+1\right)$ categories, when in fact there are fewer parameters that are active (due to identical distributions), so this is preserved by (C3): the reader may reflect on the extreme example that all distributions are identical and $\left\{X_t\right\}$ is in fact an IID process! Nevertheless, the usual (opposite) notion that the different I-distributions should be as alike as possible, will be fundamental for the process validity according to other conditions (causality/invertibility) as in the next section.

Finally, it is stressed that the assumption of a convoluted multivariate IID time series, which is attached to the definition of the TARMA, is the necessary block for the inference that will be presented in this paper, the same way that the univariate error IID series is for the well-known ARMA.

3. Invertibility

In [1], the process $\left\{X_t\mathrm{,}t\in \mathbb{Z}\right\}$ was called table causal (based on $\left\{I_t\mathrm{,}t\in \mathbb{Z}\right\}$ ) also implying that it is strictly stationary. After setting

$\begin{array}{l}{d}^h{f}_x\left(I_t^{\left(0_{l_1-1}{i}_{l_1}{0}_{l_2-l_1-1}{i}_{l_2}\cdots i_{l_h}{0}_{p+q-l_h}\right)}\right)\\ \mathrm{<mo>\; :\; </mo>}=f_x\left(I_t^{\left(0_{l_1-1}{i}_{l_1}{0}_{l_2-l_1-1}{i}_{l_2}\cdots i_{l_h}{0}_{p+q-l_h}\right)}\right)-{\displaystyle \underset{n=1}{\overset{h}{\sum }}}{f}_x\left(I_t^{\left(0_{l_1-1}\mathrm{,}{i}_{l_1}\mathrm{,}\cdots \mathrm{,}{i}_{l_{n-1}}\mathrm{,}{0}_{l_{n+1}-l_{n-1}-1}\mathrm{,}{i}_{l_{n+1}}\mathrm{,}\cdots \mathrm{,}{i}_{l_h}\mathrm{,}{0}_{p+q-l_h}\right)}\right)\\ +\cdots +{\left(-1\right)}^{h-1}{\displaystyle \underset{n\mathrm{<mo>\; =\; </mo>1}}{\overset{h}{\sum }}}{f}_x\left(I_t^{\left(0_{l_n-1}\mathrm{,}{i}_{l_n}\mathrm{,}{0}_{p+q-l_n}\right)}\right)+{\left(-1\right)}^h{f}_x(\; I\; t\; )\end{array}$

for $h=1,\cdots ,p+q$, $l_1\mathrm{,}\cdots \mathrm{,}{l}_h=\mathrm{1,}\cdots \mathrm{,}p+q$ ( $l_1<\cdots <l_h$ ), $i_{l_1},\cdots ,i_{l_h}=v_1,\cdots ,v_k$, it can be derived from (2), that

$\begin{array}{l}{f}_x\left(I_t\right)=f_x\left(X_t\right)-{\displaystyle \underset{h_p=1}{\overset{p}{\sum }}}{\displaystyle \underset{\begin{array}{c}{l}_1\mathrm{,}\cdots \mathrm{,}{l}_{h_p}=1\\ l_1<\cdots <l_{h_p}\end{array}}{\overset{p}{\sum }}}{\displaystyle \underset{i_{l_1},\cdots ,i_{l_{h_p}}=v_1,\cdots ,v_k}{\sum }}{d}^{h_p}{f}_x\left(I_t^{\left(\left(0_{l_1-1}\mathrm{,}{i}_{l_1}\mathrm{,}\cdots \mathrm{,}{i}_{l_{h_p}}\mathrm{,}{0}_{p-l_{h_p}}\right)\mathrm{<mo>\; |\; </mo>}{0}_q\right)}\right)\\ \cdot \left({\displaystyle \underset{r=1}{\overset{h_p}{\prod }}}{f}_{i_{l_r}}\left(X_{t-l_r}\right)\right)-{\displaystyle \underset{h_q=1}{\overset{q}{\sum }}}{\displaystyle \underset{h_p=0}{\overset{p}{\sum }}}{\displaystyle \underset{\begin{array}{c}{l}_1,\cdots ,l_{h_p}=1\\ l_1<\cdots <l_{h_p}\end{array}}{\overset{p}{\sum }}}{\displaystyle \underset{\begin{array}{c}{n}_1,\cdots ,n_{h_q}=1\\ n_1<\cdots <n_{h_q}\end{array}}{\overset{q}{\sum }}}{\displaystyle \underset{i_{l_1},\cdots ,i_{l_{h_p}},j_{n_1},\cdots ,j_{n_{h_q}}=v_1,\cdots ,v_k}{\sum }}{d}^{h_p+h_q}\\ \cdot f_x\left(I_t^{\left(\left(0_{l_1-1}\mathrm{,}{i}_{l_1}\mathrm{,}\cdots \mathrm{,}{i}_{l_{h_p}}\mathrm{,}{0}_{p-l_{h_p}}\right)\mathrm{<mo>\; |\; </mo>}\left(0_{n_1-1}\mathrm{,}{j}_{n_1}\mathrm{,}\cdots \mathrm{,}{j}_{n_{h_q}}\mathrm{,}{0}_{q-n_{h_q}}\right)\right)}\right)\left({\displaystyle \underset{r=1}{\overset{h_q}{\prod }}}{f}_{j_{n_r}}\left(I_{t-n_r}\right)\right)\left({\displaystyle \underset{r=1}{\overset{h_p}{\prod }}}{f}_{i_{l_r}}\left(X_{t-l_r}\right)\right):\end{array}$ (3)

similarly was done for causality ( [1]) and the representation

$\begin{array}{c}{f}_x\left(X_t\right)=f_x\left(I_t\right)+{\displaystyle \underset{h\in \mathbb{N}}{\sum }}{\displaystyle \underset{\begin{array}{c}{n}_1,\cdots ,n_h\in \mathbb{N}\\ n_1<\cdots <n_h\end{array}}{\sum }}{\displaystyle \underset{j_{n_1},\cdots ,j_{n_h}=v_1,\cdots ,v_k}{\sum }}{f}_x{\Psi }_{t,\left(n_1,\cdots ,n_h\right)}\left(j_{n_1},\cdots ,j_{n_h}\right)\\ \cdot \left({\displaystyle \underset{r=1}{\overset{h}{\prod }}}{f}_{j_{n_r}}\left(I_{t-n_r}\right)\right).\end{array}$ (4)

Here, a new substitution of $\left({\displaystyle {\prod }_{r=1}^{h_q}}{f}_{j_{n_r}}\left(I_{t-n_r}\right)\right)$ and so on, will also lead to a countably infinite representation (of the upcoming form (5)). An adequate definition for invertibility follows.

Definition 3.1: The process $\left\{X_t\mathrm{,}t\in \mathbb{Z}\right\}$ as defined in (1) will be called invertible, in the sense that it can be written

$\begin{array}{c}{f}_x\left(X_t\right)=f_x\left(I_t\right)+{\displaystyle \underset{h\in \mathbb{N}}{\sum }}{\displaystyle \underset{\begin{array}{c}{l}_1,\cdots ,l_h\in \mathbb{N}\\ l_1<\cdots <l_h\end{array}}{\sum }}{\displaystyle \underset{i_{l_1},\cdots ,i_{l_h}=v_1,\cdots ,v_k}{\sum }}{f}_x{\Phi }_{t,\left(l_1,\cdots ,l_h\right)}\left(i_{l_1},\cdots ,i_{l_h}\right)\\ \cdot \left({\displaystyle \underset{r=1}{\overset{h}{\prod }}}{f}_{i_{i_r}}\left(X_{t-l_r}\right)\right)\end{array}$ (5)

for $x=v_1,\cdots ,v_k$ and $t\in \mathbb{Z}$ ; the random variable $f_x{\Phi }_{t,\left(l_1,\cdots ,l_h\right)}\left(i_{l_1},\cdots ,i_{l_h}\right)$, $\left(h\ge 1\right)$ is independent of $I_{t+n}^{\left(\mathrm{...}\right)}$, $n\in \mathbb{N}$, $I_{t-l}^{\left(\mathrm{...}\right)}$, $l\ge l_h$, $l\in \mathbb{N}$ (it is a function of $I_{t-l}^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}$, $0\le l\le l_h-1$ ) and remains unchanged for $t\in \mathbb{Z}$.

Additionally, it must hold that the probabilities

$\mathbb{P}\left(X_t=x|X_{t-l}=i_l,l\in \mathbb{N}\right),i_l=0,v_1,\cdots ,v_k,l\in \mathbb{N}$ (6)

can be uniquely determined from (5).

Remark 1: For $h\mathrm{,}{h}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{h}_L\in \mathbb{N}$, $h_L\le h^{\mathrm{<mo>\; *\; </mo>}}\le h$, $l_1,\cdots ,l_h\in \mathbb{N}$, ( $l_1<\cdots <l_h$ ), consider $\left\{l_1^{*},\cdots ,l_{h_L}^{*}\equiv l_{h^{*}}\right\}\subseteq \left\{l_1,\cdots ,l_h\right\}$, $l_1^{*}<\cdots <l_{h_L}^{*}$ and any $i_{l_1},\cdots ,i_{l_h}=v_1,\cdots ,v_k$, and write the events

${\mathcal{B}}_{\left(t,h\right),\left(\left(l_1,\cdots ,l_h\right),\left(i_{l_1},\cdots ,i_{l_h}\right)\right)}:=\left\{I_{t-l_1}=i_{l_1},\cdots ,I_{t-l_h}=i_{l_h},I_{t-l}=0,l\in \mathbb{N},l\ne l_1,\cdots ,l_h\right\},$

${\mathcal{C}}_{\left(t,h\right),\left(\left(l_1,\cdots ,l_h\right),\left(i_{l_1},\cdots ,i_{l_h}\right)\right)}:=\left\{X_{t-l_1}=i_{l_1},\cdots ,X_{t-l_h}=i_{l_h},X_{t-l}=0,l\in \mathbb{N},l\ne l_1,\cdots ,l_h\right\}$

and the probability of interest

$\mathbb{P}\left(f_x{\Phi }_{t\mathrm{,}\left(l_1^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\cdots \mathrm{,}{l}_{h_L}^{\mathrm{<mo>\; *\; </mo>}}\right)}\left(i_{l_1^{\mathrm{<mo>\; *\; </mo>}}}\mathrm{,}\cdots \mathrm{,}{i}_{l_{h_L}^{\mathrm{<mo>\; *\; </mo>}}}\right)\le y\mathrm{<mo>\; |\; </mo>}{\mathcal{C}}_{\left(t,h\right),\left(\left(l_1,\cdots ,l_h\right),\left(i_{l_1},\cdots ,i_{l_h}\right)\right)}\right)\mathrm{.}$

As opposed to the probability

$\begin{array}{l}\mathbb{P}\left(f_x{\Psi }_{t\mathrm{,}\left(l_1^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\cdots \mathrm{,}{l}_{h_L}^{\mathrm{<mo>\; *\; </mo>}}\right)}\left(i_{l_1^{\mathrm{<mo>\; *\; </mo>}}}\mathrm{,}\cdots \mathrm{,}{i}_{l_{h_L}^{\mathrm{<mo>\; *\; </mo>}}}\right)\le y|{\mathcal{B}}_{\left(t,h\right),\left(\left(l_1,\cdots ,l_h\right),\left(i_{l_1},\cdots ,i_{l_h}\right)\right)}\right)\\ \equiv \mathbb{P}(f_x{\Psi }_{t\mathrm{,}\left(l_1^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\cdots \mathrm{,}{l}_{h_L}^{\mathrm{<mo>\; *\; </mo>}}\right)}\left(i_{l_1^{\mathrm{<mo>\; *\; </mo>}}}\mathrm{,}\cdots \mathrm{,}{i}_{l_{h_L}^{\mathrm{<mo>\; *\; </mo>}}}\right)\le y|I_{t-l_1}=i_{l_1},\cdots ,I_{t-l_{h^{*}}}=i_{l_{h^{*}}},\\ \stackrel{}{}{I}_{t-l}=0,1\le l<l_{h^{*}},l\ne l_1,\cdots ,l_{h^{*}})\end{array}$

from the “causality” topic in [1] when it is used that $\left\{I_t^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}\right\}~\text{IID}$ ( $f_x{\Psi }_{t\mathrm{,}\left(l_1,\cdots ,l_h\right)}$ are functions of $I_{t-l}^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}$, $0\le l<l_h$ ), in the case of “invertibility” $\left\{X_t\right\}$ is a serially dependent series. Nevertheless, a simplification is in order, i.e., it holds that

$\begin{array}{l}\mathbb{P}\left(f_x{\Phi }_{t\mathrm{,}\left(l_1^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\cdots \mathrm{,}{l}_{h_L}^{\mathrm{<mo>\; *\; </mo>}}\right)}\left(i_{l_1^{\mathrm{<mo>\; *\; </mo>}}}\mathrm{,}\cdots \mathrm{,}{i}_{l_{h_L}^{\mathrm{<mo>\; *\; </mo>}}}\right)\le y|{\mathcal{C}}_{\left(t,h\right),\left(\left(l_1,\cdots ,l_h\right),\left(i_{l_1},\cdots ,i_{l_h}\right)\right)}\right)\\ =\mathbb{P}(f_x{\Phi }_{t\mathrm{,}\left(l_1^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\cdots \mathrm{,}{l}_{h_L}^{\mathrm{<mo>\; *\; </mo>}}\right)}\left(i_{l_1^{\mathrm{<mo>\; *\; </mo>}}}\mathrm{,}\cdots \mathrm{,}{i}_{l_{h_L}^{\mathrm{<mo>\; *\; </mo>}}}\right)\le y|{\Phi }_{\left(t-1,h\right),\left(\left(l_1-1,\cdots ,l_h-1\right),\left(i_{l_1},\cdots ,i_{l_h}\right)\right)}^{*}=0,\cdots \\ {\Phi }_{\left(t-l_1,h-1\right),\left(\left(l_2-l_1,\cdots ,l_h-l_1\right),\left(i_{l_2},\cdots ,i_{l_h}\right)\right)}^{*}=i_{l_1},\cdots ,{\Phi }_{\left(t-l_{h-1},1\right),\left(l_h-l_{h-1},i_{l_h}\right)}^{*}=i_{l_{h-1}},\cdots \\ {\Phi }_{\left(t-l_h+1,1\right),\left(1,i_{l_h}\right)}^{*}=0,I_{t-l_h}=i_{l_h})\end{array}$ (7)

as it can be shown (together with the notation ${\Phi }^{\mathrm{<mo>\; *\; </mo>}}$ ).

3.1. Conditions Relating to Causality and Invertibility

To perform the TARMA parameters inference, there will be the need to somehow “squeeze” the random coefficients $f_x{\Phi }_{t\mathrm{,}\left(l_1,\cdots ,l_h\right)}$ to become smaller as $l_h\to \infty$.

To manage convergence results, first set the remainder

$\begin{array}{c}{f}_x{R}_{t,>l}:={\displaystyle \underset{r=l+1}{\overset{\infty }{\sum }}}{\displaystyle \underset{i_r=v_1,\cdots ,v_k}{\sum }}{f}_x{\Phi }_{t,r}\left(i_r\right)\cdot f_{i_r}\left(X_{t-r}\right)\\ +{\displaystyle \underset{r_1=1}{\overset{l}{\sum }}}{\displaystyle \underset{r_2=l+1}{\overset{\infty }{\sum }}}{\displaystyle \underset{i_{r_1},i_{r_2}=v_1,\cdots ,v_k}{\sum }}{f}_x{\Phi }_{t,\left(r_1,r_2\right)}\left(i_{r_1},i_{r_2}\right)\cdot f_{i_{r_1}}\left(X_{t-r_1}\right)f_{i_{r_2}}\left(X_{t-r_2}\right)\\ +{\displaystyle \underset{\begin{array}{c}{r}_1,r_2=l+1\\ r_1<r_2\end{array}}{\overset{\infty }{\sum }}}{\displaystyle \underset{i_{r_1},i_{r_2}=v_1,\cdots ,v_k}{\sum }}{f}_x{\Phi }_{t,\left(r_1,r_2\right)}\left(i_{r_1},i_{r_2}\right)\cdot f_{i_{r_1}}\left(X_{t-r_1}\right)f_{i_{r_2}}\left(X_{t-r_2}\right)+\cdots \end{array}$

from the “invertible” representation (5). Similarly, it has been set by [1] that

$\begin{array}{c}{f}_x{r}_{t,>l}:={\displaystyle \underset{r=l+1}{\overset{\infty }{\sum }}}{\displaystyle \underset{j_r=v_1,\cdots ,v_k}{\sum }}{f}_x{\Psi }_{t,r}\left(j_r\right)\cdot f_{j_r}\left(I_{t-r}\right)\\ +{\displaystyle \underset{r_1=1}{\overset{l}{\sum }}}{\displaystyle \underset{r_2=l+1}{\overset{\infty }{\sum }}}{\displaystyle \underset{j_{r_1},j_{r_2}=v_1,\cdots ,v_k}{\sum }}{f}_x{\Psi }_{t,\left(r_1,r_2\right)}\left(j_{r_1},j_{r_2}\right)\cdot f_{j_{r_1}}\left(I_{t-r_1}\right)f_{j_{r_2}}\left(I_{t-r_2}\right)\\ +{\displaystyle \underset{\begin{array}{c}{r}_1,r_2=l+1\\ r_1<r_2\end{array}}{\overset{\infty }{\sum }}}{\displaystyle \underset{j_{r_1},j_{r_2}=v_1,\cdots ,v_k}{\sum }}{f}_x{\Psi }_{t,\left(r_1,r_2\right)}\left(j_{r_1},j_{r_2}\right)\cdot f_{j_{r_1}}\left(I_{t-r_1}\right)f_{j_{r_2}}\left(I_{t-r_2}\right)+\cdots \end{array}$

from the “causal” representation (4).

Consider generic constants $C>0$ and $\alpha \in \left(\mathrm{0,1}\right)$. The condition is presented below:

(C4): The parameters of the TARMA equation are such that it holds that

$\begin{array}{l}\mathbb{E}\left\{\left|f_x{R}_{t,>l}\right|\mathrm{<mo>\; |\; </mo>}{X}_{t-n}=x_n\mathrm{,}n\in \mathbb{N}\right\}\le C\cdot {\alpha }^{l+1}\mathrm{,}\text{for}x=v_1\mathrm{,}\cdots \mathrm{,}{v}_k\mathrm{,}\\ t\in \mathbb{Z}\mathrm{,}l\in {\mathbb{N}}_0\text{and}\text{any}{x}_n=\mathrm{0,}{v}_1\mathrm{,}\cdots \mathrm{,}{v}_k\mathrm{,}n\in \mathbb{N}\mathrm{.}\end{array}$ (8)

(C5): The parameters of the TARMA equation are such that it holds that

$\mathbb{E}\left|f_x{r}_{t,>l}\right|\le C\cdot {\alpha }^{l+1},\text{for}x=v_1,\cdots ,v_k,t\in \mathbb{Z},l\in {\mathbb{N}}_0.$ (9)

To justify (C4), see that it can be written that

$f_x\left(X_t\right)=f_x\left(I_t\right)+{\displaystyle \underset{l\in {\mathbb{N}}_0}{\sum }}\left(f_x{R}_{t,>l}-f_x{R}_{t,>l+1}\right),$

so that the conditional probability of interest (6) is bounded by

$\mathbb{P}\left(I_t=x|X_{t-n}=i_n,n\in \mathbb{N}\right)+{\displaystyle \underset{l\in {\mathbb{N}}_0}{\sum }}{C}^{*}\cdot {\alpha }^{l+1},$

where a converging geometric series is involved and, under causality, it is $\mathbb{P}\left(I_t=x|X_{t-n}=i_n,n\in \mathbb{N}\right)\equiv {\pi }_x$ ; straight from (8), it holds that $\begin{array}{l}\mathbb{E}\left\{\left|f_x{R}_{t,>l}-f_x{R}_{t,>l+1}\right||X_{t-n}=i_n,n\in \mathbb{N}\right\}\\ \le \mathbb{E}\left\{\left|f_x{R}_{t,>l}\right||X_{t-n}=i_n\right\}+\mathbb{E}\left\{\left|f_x{R}_{t,>l+1}\right||X_{t-n}=i_n\right\}\end{array}$, so that it has been inserted $C^{*}:=C\cdot \left(1+\alpha \right)$.

Regarding (C5) and how it secures that all joint probabilities can be safely bounded, the answer is easy to show.

Since (C4) (as opposed to (C5)) involves conditional expectations, it might be wished for $\left|f_x{R}_{t,>l}\right|$ given “...” to simplify the condition “...” from $\left\{X_{t-n}=x_n,n\in \mathbb{N}\right\}$ to $\left\{X_{t-n}=x_n,n=1,\cdots ,l\right\}$. According to the argument laid out above, it is in fact the absolute value of

$\begin{array}{l}{f}_x{R}_{t,>l}-f_x{R}_{t,>l+1}\equiv {\displaystyle \underset{i_{l+1}=v_1,\cdots ,v_k}{\sum }}{f}_x{\Phi }_{t,l+1}\left(i_{l+1}\right)\cdot f_{i_{l+1}}\left(X_{t-\left(l+1\right)}\right)\\ +{\displaystyle \underset{r_1=1}{\overset{l}{\sum }}}{\displaystyle \underset{i_{r_1},i_{l+1}=v_1,\cdots ,v_k}{\sum }}{f}_x{\Phi }_{t,\left(r_1,l+1\right)}\left(i_{r_1},i_{r_2}\right)\cdot f_{i_{r_1}}\left(X_{t-r_1}\right)f_{i_{l+1}}\left(X_{t-\left(l+1\right)}\right)+\cdots \end{array}$

(instead of that of $f_x{R}_{t,>l}$ ) that is needed, which is “contained” within the random coefficients $f_x{\Phi }_{t\mathrm{,}\left(l_1,\cdots ,l_h\right)}$, $l_1<\cdots <l_h\le l+1$ (that are functions of $I_{t-n}^{\left(\mathrm{...}\right)}$, $n=0,1,\cdots ,l$ ). In the case that $X_{t-n}=0$, $n\ge l+1$, it is clear from Remark 1 that (under causality) a simplification is possible, though not concerning the variables X in the condition: otherwise, it’d better not be attempted. The interested reader may look at Section 3.2 to verify what happens in the general case given the requirement in (C4) as it is.

For causality, due to the difficulties in determining the random coefficients $f_x{\Psi }_t$ (as functions of $I_{t-l}^{\left(\mathrm{...}\right)}$, $l\in {\mathbb{N}}_0$ ), [1] resorted to an alternative (to (4)) representation: then based on the new (rather than (4)) form, salvaged a condition relating to causality. The equivalent representation for invertibility may be demonstrated with a Proposition 1 and Propositions 2 and 3 may lead to an invertibility-related condition. Nevertheless, it might be stressed here that Proposition 2 (and, consequently, Proposition 3) is established using the prerequisite of causality; this is because an “invertible” representation relies on the index variables $f\left(X_{t-l}\right)$ (not $f\left(I_{t-l}\right)$ ) and, under causality, it can be certified that $X_{t-l}$, $l\in \mathbb{N}$ is independent of $I_t^{\left(\mathrm{...}\right)}$.

Furthermore, remember that the causality consideration is attached to the definition of a TARMA process, as it guarantees probability stationarity, which is what this theory is about: so this is before the results for the inference are seeked. In the contrary, the contribution of invertibility will shine, when the weak convergence (consistency) of the maximum likelihood estimators for the relevant TARMA probabilities is established. Nevertheless, it is not wished to undermine the value of invertibility as compared to that of causality: after all, both type coefficients $f_x{\Psi }_t$ and $f_x{\Phi }_t$ are built as functions of the $I^{\left(\mathrm{...}\right)}$ variables from present and past, using a similar mechanism. The conditions obtained seem to complete rather than contradict each other too, so it is clear that causality and invertibility need to work together for the TARMA body to stand straight.

The reader needs to acquaint themselves with the following notation: it is set for $x\mathrm{<mo>\; =\; </mo>}{v}_1\mathrm{,}\dots \mathrm{,}{v}_k$, that

${\gamma }_x^{\left(\nu \right)}\mathrm{<mo>\; :\; </mo>}=\{\begin{array}{ll}\max \left\{\mathbb{E}\left|d^{\nu }{f}_x\left(I_t^{\left(0_p\mathrm{<mo>\; |\; </mo>}\left(0_{n_1-1}{j}_{n_1}\cdots j_{n_{\nu }}{0}_{q-n_{\nu }}\right)\right)}\right)\right|\right\}\mathrm{,}\hfill & \text{if}\nu =\mathrm{1,}\cdots \mathrm{,}q\hfill \\ \underset{i\mathrm{,}j}{\max }\left\{{\pi }_x^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}\right\}\mathrm{,}\hfill & \text{if}\nu =0\hfill \end{array}\mathrm{,}$

as well as for each $h_p=1,\cdots ,p$, $\nu =0,1,\cdots ,q$, that

${\gamma }_{p\mathrm{,}x}^{\left(\nu \right)}\left(h_p\right)\mathrm{<mo>\; :\; </mo>}=\max \left\{\mathbb{E}\left|d^{h_p+\nu }{f}_x\left(I_t^{\left(\left(0_{l_1-1}{i}_{l_1}\cdots i_{l_{h_p}}{0}_{p-l_{h_p}}\right)\mathrm{<mo>\; |\; </mo>}\left(0_{n_1-1}{j}_{n_1}\cdots j_{n_{\nu }}{0}_{q-n_{\nu }}\right)\right)}\right)\right|\right\}\mathrm{,}$

where both maximums take place for any $n_1,\cdots ,n_{\nu }=1,\cdots ,q$ ( $n_1<\cdots <n_{\nu }$ ), $j_{n_1},\cdots ,j_{n_{\nu }}=v_1,\cdots ,v_k$ ; the second also taking place for any $l_1,\cdots ,l_{h_p}=1,\cdots ,p$ ( $l_1<\cdots <l_{h_p}$ ), $i_{l_1},\cdots ,i_{l_{h_p}}=v_1,\cdots ,v_k$. It is defined for $x=v_1,\cdots ,v_k$

${\gamma }_x:=\underset{\nu =1,\cdots ,q}{\max }\left\{{\gamma }_x^{\left(\nu \right)}+{\displaystyle \underset{h_p=1}{\overset{p}{\sum }}}\left(\begin{array}{c}p\\ h_p\end{array}\right){\gamma }_{p,x}^{\left(\nu \right)}\left(h_p\right)\right\}$ and

${\gamma }^{*}:=\underset{\nu =0,\cdots ,q}{\max }\left\{\left({\displaystyle \underset{x=v_1,\cdots ,v_k}{\sum }}{\gamma }_x^{\left(\nu \right)}\right)+{\displaystyle \underset{h_p=1}{\overset{p}{\sum }}}\left(\begin{array}{c}p\\ h_p\end{array}\right)\left({\displaystyle \underset{x=v_1,\cdots ,v_k}{\sum }}{\gamma }_{p,x}^{\left(\nu \right)}\left(h_p\right)\right)\right\}.$

By relaxing the sum of maximums into a maximum of sums, ${\displaystyle {\sum }_{x=v_1,\cdots ,v_k}}{\gamma }_x^{\left(0\right)}$ is replaced by ${\max }_{i\mathrm{,}j}\left\{1-{\pi }_0^{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}\right\}$. More research is welcome on the subject of the equivalence with (C4), which was demonstrated to be sufficient for invertibility.

Remark 2: (i) Under (C4) and (C5), all joint probabilities and all probabilities $\mathbb{P}\left(X_t=x|X_{t-n}=x_n,n\in \mathbb{N}\right)$, $x,x_n=0,v_1,\cdots ,v_k$ are contained (away from infinity): it can be concluded that any joint probability is away from zero (hence away from one as well), because if any such probability, say $\mathbb{P}\left({\mathcal{X}}_t^{\left(a\right)}\right)$, was 0 it would have to be that any conditional probability as above- with ${\mathcal{X}}_{t^{\mathrm{<mo>\; *\; </mo>}}}^{\left(a\right)}\subseteq \left\{X_{t-n}=x_n\mathrm{,}n\in \mathbb{N}\right\}$ for some $t^{\mathrm{<mo>\; *\; </mo>}}$ - would not be bounded and properly defined; then it can be concluded (by division of a joint non zero probability over a joint non infinity probability) that all conditional probabilities are strictly larger than zero (hence smaller than one as well). Consequently, it will be taken that (C4) and (C5) can suffice for (C2).

(ii) For the $\pi$ that yield a causal and invertible TARMA series $\left\{X_t\right\}$, it can be arranged under (8) in (C4) (or, (9) in (C5)) that

$\left|\frac{\partial }{\partial \pi }\mathbb{E}\left(f_x{R}_{t,>l}|X_{t-n}=x_n,n\in \mathbb{N}\right)\right|\le C\cdot {\alpha }^{l+1},\left(\text{or}\left|\frac{\partial }{\partial \pi }\mathbb{E}\left(f_x{r}_{t,>l}\right)\right|\le C\cdot {\alpha }^{l+1}\right)$ (10)

for some $0<C<\infty$ and $\alpha \in \left(\mathrm{0,1}\right)$, and

$\left|\frac{{\partial }^2}{\partial \pi \partial \pi }\mathbb{E}\left(f_x{R}_{t,>l}|X_{t-n}=x_n,n\in \mathbb{N}\right)\right|\le C\cdot {\alpha }^{l+1}\left(\text{or}\left|\frac{{\partial }^2}{\partial \pi \partial \pi }\mathbb{E}\left(f_x{r}_{t,>l}\right)\right|\le C\cdot {\alpha }^{l+1}\right)$ (11)

under (8) (or (9)): this can be shown.

It is being taken for granted that the $\frac{\partial }{\partial \pi }\mathbb{P}\left(X_t=x|X_{t-n}=x_n,n\in \mathbb{N},\pi \right)$, $\frac{{\partial }^2}{\partial \pi \partial \pi }\mathbb{P}\left(X_t=x|X_{t-n}=x_n,n\in \mathbb{N},\pi \right)$ derivatives exist.

3.2. Random Coefficient Modelling Based on the Past

For $l\mathrm{,}n\in {\mathbb{N}}_0$, $n\le l$ and a set $\left\{l_1^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\cdots \mathrm{,}{l}_n^{\mathrm{<mo>\; *\; </mo>}}\right\}\subseteq \left\{\mathrm{1,}\cdots \mathrm{,}l\right\}$, it is written (for fixed $i_{l_1^{*}},\cdots ,i_{l_n^{*}}=v_1,\cdots ,v_k$ but this is omitted from the symbol $\mathcal{N}$ ), that

${\mathcal{N}}_{t,n\le l}:=\left\{X_{t-l_1^{*}}=i_{l_1^{*}},\cdots ,X_{t-l_n^{*}}=i_{l_n^{*}},X_{t-l^{*}}=0,l^{*}=1,\cdots ,l,l^{*}\ne l_1^{*},\cdots ,l_n^{*}\right\}.$

Eventually it can be shown that

$\underset{l\to \infty }{\lim }\mathbb{P}\left(\left|\mathbb{P}\left(X_t=x\mathrm{<mo>\; |\; </mo>}{\mathcal{N}}_{t\mathrm{,}n\le l}\mathrm{,}{X}_{t-l-l^{\mathrm{<mo>\; *\; </mo>}}}\mathrm{,}{l}^{\mathrm{<mo>\; *\; </mo>}}\in \mathbb{N}\right)-\mathbb{P}\left(X_t=x\mathrm{<mo>\; |\; </mo>}{\mathcal{N}}_{t\mathrm{,}n\le l}\right)\right|\le C\cdot {\alpha }^{l+1}\right)=1$ (12)

( $n,l_1^{\mathrm{<mo>\; *\; </mo>}},\cdots ,l_n^{\mathrm{<mo>\; *\; </mo>}},i_{l_1^{\mathrm{<mo>\; *\; </mo>}}},\cdots ,i_{l_n^{\mathrm{<mo>\; *\; </mo>}}}$ remain unchanged as $l\to \infty$ ).

Remark 3: From (10) and (11) and in the same manner, it will be considered occasionally that

$\underset{l\to \infty }{\lim }\mathbb{P}\left(\left|\frac{\partial }{\partial \pi }\mathbb{P}\left(X_t=x|{\mathcal{N}}_{t,n\le l}\right)-\frac{\partial }{\partial \pi }\mathbb{P}\left(X_t=x|{\mathcal{N}}_{t,n\le l},X_{t-l^{*}},l^{*}\ge l+1\right)\right|\le C\cdot {\alpha }^{l+1}\right)=1$ (13)

and

$\underset{l\to \infty }{\lim }\mathbb{P}\left(\left|\frac{{\partial }^2}{\partial \pi \partial \pi }\mathbb{P}\left(X_t=x|{\mathcal{N}}_{t,n\le l}\right)-\frac{{\partial }^2}{\partial \pi \partial \pi }\mathbb{P}\left(X_t=x|{\mathcal{N}}_{t,n\le l},X_{t-l^{*}},l^{*}\ge l+1\right)\right|\le C\cdot {\alpha }^{l+1}\right)=1,$ (14)

respectively. A justification may be offered.

4. Maximum Likelihood Estimation

First, the parameter space, say $\Theta$, must be considered; the fixed $k\in \mathbb{N}$, categories $v_1\mathrm{,}\cdots \mathrm{,}{v}_k$ and order $p\mathrm{,}q\in {\mathbb{N}}_0$, $p+q\in \mathbb{N}$ are attached to it:

The parameter space $\Theta$ is a set that includes all candidate parameter vectors $\pi \in \Theta$ that model a process of interest $\left\{X_t\mathrm{,}t\in \mathbb{Z}\right\}$ using a TARMA (k, p, q) equation with a predetermined form of interdependence, such that conditions (C3), (C4) and (C5) are satisfied (plus any extra requirements added in this section).

Now, suppose that $\left\{X_1\mathrm{,}\cdots \mathrm{,}{X}_T\right\}$ have been made available from (1), and the target is to estimate the true parameters of the model, say ${\pi }_0\in \Theta$, based on these observations. The likelihood expressed as a function of $\pi$, takes the form

$L\left(\pi \right)={\displaystyle \underset{t=1}{\overset{T}{\prod }}}\left({\tilde{p}}_{t,t-1}\left(\pi \right)\right),$

where

${\tilde{p}}_{t,0}\left(\pi \right):=\mathbb{P}\left(X_t=x;\pi \right),\text{if}{X}_t=x,\text{where}x=0,v_1,\cdots ,v_k,$

and for the general $i\in \mathbb{N}$, it is written

${\tilde{p}}_{t,i}\left(\pi \right):=\mathbb{P}\left(X_t=x|X_{t-1},\cdots ,X_{t-i};\pi \right),\text{if}{X}_t=x,\text{where}x=0,v_1,\cdots ,v_k;$

${\tilde{p}}_{t\mathrm{,}i}$ are functions of $X_{t-1},\cdots ,X_{t-i}$ as well as $X_t$ and it is considered that these (conditional) probabilities result from $\left\{X_t\mathrm{,}t\in \mathbb{Z}\right\}$ as if those have been generated by $\pi \in \Theta$. Of course, the natural logarithm of the likelihood may be taken

$l\left(\pi \right)={\displaystyle \underset{t=1}{\overset{T}{\sum }}}\ln \left({\tilde{p}}_{t\mathrm{,}t-1}(\; \pi \; )\right)$

as the maximum likelihood estimators $\hat{\pi }$ then satisfy the equations

${\frac{\partial }{\partial \pi }l\left(\pi \right)|}_{\pi =\hat{\pi }}=0\text{or}{\displaystyle \underset{t=1}{\overset{T}{\sum }}}\frac{1}{{\tilde{p}}_{t,t-1}\left(\hat{\pi }\right)}{\frac{\partial }{\partial \pi }{\tilde{p}}_{t,t-1}\left(\pi \right)|}_{\pi =\hat{\pi }}=0.$ (15)

4.1. Weak Convergence of Estimators

For $\pi \in \Theta$ and any $t\in \mathbb{Z}$, it is written

$p_t\left(\pi \right):=\mathbb{P}\left(X_t=x|X_{t-i},i\in \mathbb{N};\pi \right),\text{if}{X}_t=x,\text{where}x=0,v_1,\cdots ,v_k,$

which is a function of $X_{t-i}$, $i\in {\mathbb{N}}_0$ and the (conditional) distribution law is generated by $\pi \in \Theta$.

According to (12) (thanks to invertibility) and for any $\pi \in \Theta$ and $x_1,\cdots ,x_T=0,v_1,\cdots ,v_k$, it holds as $T\to \infty$ that

$\begin{array}{l}(\frac{\mathbb{P}\left(X_1=x_1;\pi \right)}{\mathbb{P}\left(X_1=x_1|X_{1-n},n\in \mathbb{N};\pi \right)}\frac{\mathbb{P}\left(X_2=x_2|X_1=x_1;\pi \right)}{\mathbb{P}\left(X_2=x_2|X_1=x_1,X_{1-n},n\in \mathbb{N};\pi \right)}\cdots \\ {\frac{\mathbb{P}\left(X_T=x_T|X_{T-1}=x_{T-1},\dots ,X_1=x_1;\pi \right)}{\mathbb{P}\left(X_T=x_T|X_{T-1}=x_{T-1},\cdots ,X_1=x_1,X_{1-n},n\in \mathbb{N};\pi \right)})}^{1/T}\stackrel{P}{\to }1,\end{array}$ (16)

because the geometric mean weighs with (rather than just multiplies) the previous values: as the new values get closer to 1, so must the mean itself. The geometric mean from (T + 1) observations is a weighted product mean of the geometric mean of the T observations and the (T + 1)th observation (with a $C_A$ adjusted constant for the ratio of probabilities to one, by allowing the exponential rate for the geometric mean of ratios, this becomes apparent with the bound ${\left\{C_A^T\cdot {\alpha }^{\frac{\left(T+1\right)T}{2}}\right\}}^{1/T}\to 0$ as $T\to \infty$ ).

Once the fixated sample series $x_1\mathrm{,}\cdots \mathrm{,}{x}_T$ has been collected (this has been generated by ${\pi }_0$ ), remember that the maximum likelihood estimate $\hat{\pi }$ is set that way, such that it is true for any $\pi \in \Theta$, that

$\begin{array}{l}{\displaystyle \underset{t=1}{\overset{T}{\prod }}}\mathbb{P}\left(X_t=x_t\mathrm{<mo>\; |\; </mo>}{X}_{t-1}=x_{t-1}\mathrm{,}\cdots \mathrm{,}{X}_1=x_1\mathrm{<mo>\; ;\; </mo>}\hat{\pi }\right)\\ \ge {\displaystyle \underset{t=1}{\overset{T}{\prod }}}\mathbb{P}\left(X_t=x_t\mathrm{<mo>\; |\; </mo>}{X}_{t-1}=x_{t-1}\mathrm{,}\cdots \mathrm{,}{X}_1=x_1\mathrm{<mo>\; ;\; </mo>}\pi \right)\end{array}$

where the (conditional) probabilities (the distribution law of $X_1,\cdots ,X_T$ ) are calculated as if $\pi$ are the real parameters that generated this realization. Thanks to (16) (and invertibility), this becomes

${\left({\displaystyle \underset{t=1}{\overset{T}{\prod }}}\frac{\mathbb{P}\left(X_t=x_t|X_{t-1}=x_{t-1},\cdots ,X_1=x_1,X_{1-n},n\in \mathbb{N};\hat{\pi }\right)}{\mathbb{P}\left(X_t=x_t|X_{t-1}=x_{t-1},\cdots ,X_1=x_1,X_{1-n},n\in \mathbb{N};\pi \right)}\right)}^{1/T}\ge 1$

with probability that tends to one as $T\to \infty$. Thinking about the random variables (rather than the realizations) now, it can be written

$\begin{array}{l}{\left({\displaystyle \underset{t=1}{\overset{T}{\prod }}}\frac{p_t\left(\hat{\pi }\right)}{p_t\left(\pi \right)}\right)}^{1/T}\\ ={\displaystyle \underset{x_1,\cdots ,x_T=0,v_1,\cdots ,v_k}{\prod }}{\left\{{\left({\displaystyle \underset{t=1}{\overset{T}{\prod }}}\frac{\mathbb{P}\left(X_t=x_t|X_{t-1}=x_{t-1},\cdots ,X_1=x_1,X_{1-n},n\in \mathbb{N};\hat{p}\right)}{\mathbb{P}\left(X_t=x_t|X_{t-1}=x_{t-1},\cdots ,X_1=x_1,X_{1-n},n\in \mathbb{N};p\right)}\right)}^{1/T}\right\}}^{f_{x_1}\left(X_1\right)\cdots f_{x_T}\left(X_T\right)},\end{array}$

where $X_1,\cdots ,X_T$ are the same random variables for both numerator/denominator with realizations as in $f_{x_1}\left(X_1\right)\cdots f_{x_T}\left(X_T\right)$ (those have been generated by the real ${\pi }_0$ ), and it is verified that

$\underset{T\to \infty }{\lim }\mathbb{P}\left({\left({\displaystyle \underset{t=1}{\overset{T}{\prod }}}\frac{p_t\left(\hat{\pi }\right)}{p_t\left(\pi \right)}\right)}^{1/T}\ge 1\right)=1\text{for}\text{any}\pi \in \Theta \mathrm{.}$ (17)

In the other hand, it is true for any $\pi \in \Theta$, that

$\mathbb{E}\left({\displaystyle \underset{t=1}{\overset{T}{\prod }}}\frac{p_t\left(\pi \right)}{p_t\left({\pi }_0\right)}\right)=\mathbb{E}\left(\mathbb{E}\left({\displaystyle \underset{t=1}{\overset{T}{\prod }}}\frac{p_t\left(\pi \right)}{p_t\left({\pi }_0\right)}|X_{1-n},n\in \mathbb{N};{\pi }_0\right)\right)$

and it is essential that it is the “real” ${\pi }_0$ that generates the random variables and governs the (conditional) expectation, i.e.,

$\begin{array}{l}\mathbb{E}\left({\displaystyle \underset{t=1}{\overset{T}{\prod }}}\frac{p_t\left(\pi \right)}{p_t\left({\pi }_0\right)}|X_{1-n},n\in \mathbb{N};{\pi }_0\right)\\ ={\displaystyle \underset{x_1,\cdots ,x_T=0,v_1,\cdots ,v_k}{\sum }}\mathbb{P}\left(X_1=x_1,\cdots ,X_T=x_T|X_{1-n},n\in \mathbb{N};{\pi }_0\right)\\ \cdot {\displaystyle \underset{t=1}{\overset{T}{\prod }}}\frac{\mathbb{P}\left(X_t=x_t|X_{t-1}=x_{t-1},\cdots ,X_1=x_1,X_{1-n},n\in \mathbb{N};\pi \right)}{\mathbb{P}\left(X_t=x_t|X_{t-1}=x_{t-1},\cdots ,X_1=x_1,X_{1-n},n\in \mathbb{N};{\pi }_0\right)}\\ \equiv {\displaystyle \underset{x_1,\cdots ,x_T=0,v_1,\cdots ,v_k}{\sum }}\mathbb{P}\left(X_1=x_1,\cdots ,X_T=x_T|X_{1-n},n\in \mathbb{N};\pi \right)=1:\end{array}$

the last statement is true since these are probabilities adding to one (regardless of the $\pi \in \Theta$ ). Then it holds that

$\mathbb{E}\left({\displaystyle \underset{t=1}{\overset{T}{\prod }}}\frac{p_t\left(\pi \right)}{p_t\left({\pi }_0\right)}\right)=1\text{for}\text{any}\pi \in \Theta .$

As in Jensen’s inequality (use the function $g\left(y\right)=y^{1/T}$, $y>0$ with $g^{″}\left(y\right)=\left(1/T\right)\left(\left(1-T\right)/T\right)y^{\left(1/T\right)-2}<0$ for $T>1$ ), this last statement can be transformed to

$\mathbb{E}\left\{{\left({\displaystyle \underset{t=1}{\overset{T}{\prod }}}\frac{p_t\left(\pi \right)}{p_t\left({\pi }_0\right)}\right)}^{1/T}\right\}\le 1\text{for}\text{any}\pi \in \Theta ,$ (18)

with the equality holding if and only if $\pi ={\pi }_0$ (it is considered here that $p_t\left({\pi }_1\right)=p_t\left({\pi }_2\right)$, ${\pi }_1,{\pi }_2\in \Theta$ can only result from ${\pi }_1={\pi }_2$ ).

In order to conclude, one can combine (17) and write

$\underset{T\to \infty }{\lim }\mathbb{P}\left({\left({\displaystyle \underset{t=1}{\overset{T}{\prod }}}\frac{p_t\left(\hat{\pi }\right)}{p_t\left({\pi }_0\right)}\right)}^{1/T}\ge 1\right)=1$

with (18) that writes

$\mathbb{E}\left\{{\left({\displaystyle \underset{t=1}{\overset{T}{\prod }}}\frac{p_t\left(\hat{\pi }\right)}{p_t\left({\pi }_0\right)}\right)}^{1/T}\right\}\le 1$

and they are both satisfied when

${\left({\displaystyle \underset{t=1}{\overset{T}{\prod }}}{p}_t\left(\hat{\pi }\right)\right)}^{1/T}-{\left({\displaystyle \underset{t=1}{\overset{T}{\prod }}}{p}_t\left({\pi }_0\right)\right)}^{1/T}\stackrel{P}{\to }0\text{as}T\to \infty$

or, in other words, when

$\hat{\pi }\stackrel{P}{\to }{\pi }_0\text{as}T\to \infty \mathrm{.}$ (19)

4.2. Asymptotic Distribution

Thanks to a Taylor’s expansion, (15) can be turned to

$\begin{array}{l}{\displaystyle \underset{t=1}{\overset{T}{\sum }}}\frac{1}{{\tilde{p}}_{t,t-1}\left({\pi }_0\right)}{\frac{\partial }{\partial {\pi }_i}{\tilde{p}}_{t,t-1}\left(\pi \right)|}_{\pi ={\pi }_0}+{\displaystyle \underset{t=1}{\overset{T}{\sum }}}{\displaystyle \underset{j}{\sum }}\{\frac{1}{{\tilde{p}}_{t,t-1}\left({\pi }_0\right)}{\frac{{\partial }^2}{\partial {\pi }_j\partial {\pi }_i}{\tilde{p}}_{t,t-1}\left(\pi \right)|}_{\pi ={\pi }_0}\\ -\frac{1}{{\tilde{p}}_{t,t-1}^2\left({\pi }_0\right)}\left({\frac{\partial }{\partial {\pi }_j}{\tilde{p}}_{t,t-1}\left(\pi \right)|}_{\pi ={\pi }_0}\right)\left({\frac{\partial }{\partial {\pi }_i}{\tilde{p}}_{t,t-1}\left(\pi \right)|}_{\pi ={\pi }_0}\right)\}\left({\hat{\pi }}_j-{\pi }_{j,0}\right)\\ +{\displaystyle \underset{t=1}{\overset{T}{\sum }}}{E}_{t,i}\left(\hat{\pi }\right)=0,\end{array}$

where i (and j) are indexes that refer to all different scalar $\pi$ in $\pi \in \Theta$ and $E_{t\mathrm{,}i}\left(\hat{\pi }\right)$ are scalars; to understand their role better, consider the expansion for fixed t and define the function

$\begin{array}{c}{e}_{t,i}\left(\pi \right):=\frac{1}{{\tilde{p}}_{t,t-1}\left(\pi \right)}\frac{\partial }{\partial {\pi }_i}{\tilde{p}}_{t,t-1}\left(\pi \right)\\ -{\displaystyle \underset{x=0,1}{\sum }}\frac{{\partial }^x}{\partial {\pi }^x}{\left(\frac{1}{{\tilde{p}}_{t,t-1}\left(\pi \right)}\frac{\partial }{\partial {\pi }_i}{\tilde{p}}_{t,t-1}\left(\pi \right)\right)|}_{\pi ={\pi }_0}\frac{{\left(\pi -{\pi }_0\right)}^x}{x!}\end{array}$

with “ $x=0$ ” applying no derivative at all to the function ( ${\left(\pi -{\pi }_0\right)}^0\equiv 1$ ) and “ $x=1$ ” being the usual row vector of first derivatives (times a column vector): then it is clear that $e\left({\pi }_0\right)=0$. Due to (19) and the continuity of the function e at ${\pi }_0$, it can be concluded that $e\left(\hat{\pi }\right)\stackrel{P}{\to }0$ as $T\to \infty$ (without worrying about the $\left(t\mathrm{,}t-1\right)$ label, i.e., the convergence to zero takes place anyway for $T\to \infty$ ). It will be taken for granted then that

$\frac{{\displaystyle {\sum }_{t=1}^T{E}_{t\mathrm{,}i}\left(\hat{\pi }\right)}}{\mathrm{...}}\stackrel{P}{\to }0\text{as}T\to \infty \mathrm{,}$ (20)

where “…” is the convenient divisor T or even T^1/2. Note that (20) has been justified without the presumption of existence of a higher than second derivative of the conditional probabilities under study.

After omitting the extra terms, all the equations are stacked to come up with

${\displaystyle \underset{t=1}{\overset{T}{\sum }}}\left(\frac{1}{{\tilde{p}}_{t,t-1}^2\left({\pi }_0\right)}{\tilde{\delta }}_t{\tilde{\delta }}_t^{\tau }-\frac{1}{{\tilde{p}}_{t,t-1}\left({\pi }_0\right)}{\tilde{D}}_t\right)\left(\hat{\pi }-{\pi }_0\right)={\displaystyle \underset{t=1}{\overset{T}{\sum }}}\frac{1}{{\tilde{p}}_{t,t-1}\left({\pi }_0\right)}{\tilde{\delta }}_t,$

where ${\tilde{\delta }}_t$ is the column vector of ${\frac{\partial }{\partial {\pi }_i}{\tilde{p}}_{t\mathrm{,}t-1}\left(\pi \right)|}_{\pi ={\pi }_0}$, $i=\mathrm{1,2,}\cdots$ and ${\tilde{D}}_t$ is the matrix with elements ${\frac{{\partial }^2}{\partial {\pi }_i\partial {\pi }_j}{\tilde{p}}_{t\mathrm{,}t-1}\left(\pi \right)|}_{\pi ={\pi }_0}$, $i,j=1,2,\cdots$ and “ $\tau$ ” stands for the transpose operator. It is re-written that

$\begin{array}{l}{T}^{1/2}\left(\hat{\pi }-{\pi }_0\right)\\ ={\left\{\frac{1}{T}{\displaystyle \underset{t=1}{\overset{T}{\sum }}}\left(\frac{1}{{\tilde{p}}_{t\mathrm{,}t-1}^2\left({\pi }_0\right)}{\tilde{\delta }}_t{\tilde{\delta }}_t^{\tau }-\frac{1}{{\tilde{p}}_{t\mathrm{,}t-1}\left({\pi }_0\right)}{\tilde{D}}_t\right)\right\}}^{-1}{T}^{-1/2}{\displaystyle \underset{t=1}{\overset{T}{\sum }}}\frac{1}{{\tilde{p}}_{t\mathrm{,}t-1}\left({\pi }_0\right)}{\tilde{\delta }}_t\mathrm{.}\end{array}$ (21)

Next, consider ${\delta }_t$ to be the column vector of ${\frac{\partial }{\partial {\pi }_i}{p}_t\left(\pi \right)|}_{\pi ={\pi }_0}$, $i=\mathrm{1,2,}\cdots$.

Lemma 4.1: It holds that

$T^{-1/2}{\displaystyle \underset{t=1}{\overset{T}{\sum }}}\left(\frac{1}{p_t\left({\pi }_0\right)}{\delta }_t-\frac{1}{{\tilde{p}}_{t\mathrm{,}t-1}\left({\pi }_0\right)}{\tilde{\delta }}_t\right)\stackrel{P}{\to }0\text{as}T\to \infty \mathrm{.}$

Proof: This may be presented.

Once Lemma 4.1 has been established, Equation (22) should be replaced by

$T^{1/2}\left(\hat{\pi }-{\pi }_0\right)={\left\{\frac{1}{T}{\displaystyle \underset{t=1}{\overset{T}{\sum }}}\left(\frac{1}{{\tilde{p}}_{t\mathrm{,}t-1}^2\left({\pi }_0\right)}{\tilde{\delta }}_t{\tilde{\delta }}_t^{\tau }-\frac{1}{{\tilde{p}}_{t\mathrm{,}t-1}\left({\pi }_0\right)}{\tilde{D}}_t\right)\right\}}^{-1}{T}^{-1/2}{\displaystyle \underset{t=1}{\overset{T}{\sum }}}\frac{1}{p_t\left({\pi }_0\right)}{\delta }_t\mathrm{.}$ (22)

To proceed further, consider $D_t$ to be the matrix with elements ${\frac{{\partial }^2}{\partial {\pi }_i\partial {\pi }_j}{p}_t\left(\pi \right)|}_{\pi ={\pi }_0}$, $i\mathrm{,}j=\mathrm{1,}\cdots$.

Lemma 4.2: It holds that

$\frac{1}{T}{\displaystyle \underset{t=1}{\overset{T}{\sum }}}\left(\left(\frac{1}{{\tilde{p}}_{t\mathrm{,}t-1}\left({\pi }_0\right)}{\tilde{D}}_t-\frac{1}{{\tilde{p}}_{t\mathrm{,}t-1}^2\left({\pi }_0\right)}{\tilde{\delta }}_t{\tilde{\delta }}_t^{\tau }\right)-\left(\frac{1}{p_t\left({\pi }_0\right)}{D}_t-\frac{1}{p_t^2\left({\pi }_0\right)}{\delta }_t{\delta }_t^{\tau }\right)\right)\stackrel{P}{\to }O$

as $T\to \infty$, where $O$ (use bold for emphasis) is the square matrix of zeros.

Proof: This may be presented.

Once Lemma 4.2 has been established, Equation (22) should be replaced by

$T^{1/2}\left(\hat{\pi }-{\pi }_0\right)={\left\{\frac{1}{T}{\displaystyle \underset{t=1}{\overset{T}{\sum }}}\left(\frac{1}{p_t^2\left({\pi }_0\right)}{\delta }_t{\delta }_t^{\tau }-\frac{1}{p_t\left({\pi }_0\right)}{D}_t\right)\right\}}^{-1}{T}^{-1/2}{\displaystyle \underset{t=1}{\overset{T}{\sum }}}\frac{1}{p_t\left({\pi }_0\right)}{\delta }_t\mathrm{.}$ (23)

Theorem 4.3: It holds that

$T^{-1/2}{\displaystyle \underset{t=1}{\overset{T}{\sum }}}\frac{1}{p_t\left({\pi }_0\right)}{\delta }_t\stackrel{D}{\to }N\left(0\mathrm{,}\text{Var}\left(\frac{1}{p_t\left({\pi }_0\right)}{\delta }_t\right)\right)\text{as}T\to \infty \mathrm{.}$

Proof: This may be presented.

Lemma 4.4: It holds that

$\frac{1}{T}{\displaystyle \underset{t=1}{\overset{T}{\sum }}}\frac{1}{p_t^2\left({\pi }_0\right)}{\delta }_t{\delta }_t^{\tau }\stackrel{P}{\to }\text{Var}\left(\frac{1}{p_t\left({\pi }_0\right)}{\delta }_t\right)\text{as}T\to \infty \mathrm{.}$

Proof: This may be presented.

Lemma 4.5: It holds that

$\frac{1}{T}{\displaystyle \underset{t=1}{\overset{T}{\sum }}}\frac{1}{p_t\left({\pi }_0\right)}{D}_t\stackrel{P}{\to }O\text{as}T\to \infty \mathrm{.}$

Proof: This may be presented.

5. Empirical Illustrations

This section serves in practice, when for special cases of TARMA models the maximum likelihood estimation takes place, so that (i) it is refreshed how to compute joint probabilities to insert in the likelihood (and estimate) and, more importantly, (ii) it is examined how well the estimators perform for moderate sample sizes: for both (i), (ii) particular interest lies in $q\ge 1$, when “naïve” moments estimates cannot be computed directly from the data, as opposed to the more traditional TAR cases. At this point, it is reminded that the TAR is a special case of the TARMA model with AR parameters only and explicit results that have been obtained for the Markov chains. The TARMA with the inclusion of MA parts is a parsimonious way to succeed the infinite order of a TAR model and where the paper interest lies.

Hence the groundwork equation

$\begin{array}{c}{X}_t=I_t\left(1-X_{t-1}\right)\left(1-I_{t-1}\right)+I_t^{\left(\mathrm{1\; <mo>\; |\; </mo>0}\right)}{X}_{t-1}\left(1-I_{t-1}\right)\\ +I_t^{\left(\mathrm{0\; <mo>\; |\; </mo>1}\right)}\left(1-X_{t-1}\right)I_{t-1}+I_t^{\left(\mathrm{1\; <mo>\; |\; </mo>1}\right)}{X}_{t-1}{I}_{t-1}\end{array}$ (24)

will, in general, govern a series $\left\{X_t\right\}~\text{TARMA}\left(\mathrm{1,1,1}\right)$ (with realizations in $\left\{\mathrm{0,1}\right\}$ ). By writing ${\psi }_{\left(i|j\right)}:=\mathbb{P}\left(X_t=i,I_t=j\right)$, $i,j=0,1$, then from the parameters $\pi =\mathbb{P}\left(I_t=1\right)$ and ${\pi }_{|v}^{\ast \left(i\mathrm{<mo>\; |\; </mo>}j\right)}=\mathbb{P}\left(I_t^{\left(i|j\right)}=1|I_t=v\right)$, $\left(i\mathrm{,}j\right)\ne \left(\mathrm{0,0}\right)$, $v=0,1$ ( ${\pi }^{\left(i|j\right)}=\mathbb{P}\left(I_t^{\left(i|j\right)}=1\right)\equiv \pi {\pi }_{|1}^{\ast \left(i|j\right)}+\left(1-\pi \right){\pi }_{|0}^{\ast \left(i|j\right)}$ ), [1] has contributed a methodology that eventually (under causality) computes the $\psi$ as in

$\begin{array}{l}{\psi }_{\left(1|1\right)}=\left\{\pi \left[\left(1-\left(1-\pi \right){\pi }_{|0}^{\ast \left(1|0\right)}\right)\left(1-\pi \left(1-{\pi }_{|1}^{\ast \left(0|1\right)}\right)\right)-\left(1-\pi \right){\pi }_{|0}^{\ast \left(0|1\right)}\pi \left(1-{\pi }_{|1}^{\ast \left(1|0\right)}\right)\right]\right\}\\ /\{\left(1+\pi -\pi {\pi }_{|1}^{\ast \left(1|1\right)}\right)\left[\left(1-\left(1-\pi \right){\pi }_{|0}^{\ast \left(1|0\right)}\right)\left(1-\pi \left(1-{\pi }_{|1}^{\ast \left(0|1\right)}\right)\right)-\left(1-\pi \right){\pi }_{|0}^{\ast \left(0|1\right)}\pi \left(1-{\pi }_{|1}^{\ast \left(1|0\right)}\right)\right]\\ +\left(\pi -\pi {\pi }_{|1}^{\ast \left(1|0\right)}\right)\left[\left(1-\pi \left(1-{\pi }_{|1}^{\ast \left(0|1\right)}\right)\right)\left(1-\pi \right){\pi }_{|0}^{\ast \left(1|1\right)}+\left(1-\pi \right){\pi }_{|0}^{\ast \left(0|1\right)}\pi \left(1-{\pi }_{|1}^{\ast \left(1|1\right)}\right)\right]\\ +\left(\pi -\pi {\pi }_{|1}^{\ast \left(0|1\right)}\right)\left[\pi \left(1-{\pi }_{|1}^{\ast \left(1|0\right)}\right)\left(1-\pi \right){\pi }_{|0}^{\ast \left(1|1\right)}+\left(1-\left(1-\pi \right){\pi }_{|0}^{\ast \left(1|0\right)}\right)\pi \left(1-{\pi }_{|1}^{\ast \left(1|1\right)}\right)\right]\},\end{array}$ (25)

followed by

${\psi }_{\left(1|0\right)}=\frac{\left(1-\pi \left(1-{\pi }_{|1}^{\ast \left(0|1\right)}\right)\right)\left(1-\pi \right){\pi }_{|0}^{\ast \left(1|1\right)}+\left(1-\pi \right){\pi }_{|0}^{\ast \left(0|1\right)}\pi \left(1-{\pi }_{|1}^{\ast \left(1|1\right)}\right)}{\left(1-\left(1-\pi \right){\pi }_{|0}^{\ast \left(1|0\right)}\right)\left(1-\pi \left(1-{\pi }_{|1}^{\ast \left(0|1\right)}\right)\right)-\left(1-\pi \right){\pi }_{|0}^{\ast \left(0|1\right)}\pi \left(1-{\pi }_{|1}^{\ast \left(1|0\right)}\right)}{\psi }_{\left(1|1\right)}$ (26)

and

${\psi }_{\left(0|1\right)}=\frac{\pi \left(1-{\pi }_{|1}^{\ast \left(1|0\right)}\right)\left(1-\pi \right){\pi }_{|0}^{\ast \left(1|1\right)}+\left(1-\left(1-\pi \right){\pi }_{|0}^{\ast \left(1|0\right)}\right)\pi \left(1-{\pi }_{|1}^{\ast \left(1|1\right)}\right)}{\left(1-\left(1-\pi \right){\pi }_{|0}^{\ast \left(1|0\right)}\right)\left(1-\pi \left(1-{\pi }_{|1}^{\ast \left(0|1\right)}\right)\right)-\left(1-\pi \right){\pi }_{|0}^{\ast \left(0|1\right)}\pi \left(1-{\pi }_{|1}^{\ast \left(1|0\right)}\right)}{\psi }_{\left(1|1\right)}$

(and ${\psi }_{\left(0|0\right)}=1-{\psi }_{\left(1|0\right)}-{\psi }_{\left(0|1\right)}-{\psi }_{\left(1|1\right)}$ ). (27)

Additionally, after writing for $h\in \mathbb{N}$ that ${\psi }_{\left(\left(i_1,\cdots ,i_{h+1}\right)|j\right)}\left(h\right):=\mathbb{P}\left(X_t=i_1,\cdots ,X_{t-h}=i_{h+1},I_t=j\right)$ (and ${\psi }_{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}\left(0\right)\equiv {\psi }_{\left(i\mathrm{<mo>\; |\; </mo>}j\right)}$ ), the methodology also contributes the recursive formulae

${\psi }_{\left(\left(1,0,i_2,\cdots ,i_{h+1}\right)|1\right)}\left(h+1\right)=\pi {\psi }_{\left(\left(0,i_2,\cdots ,i_{h+1}\right)|0\right)}\left(h\right)+\pi {\pi }_{|1}^{\ast \left(0|1\right)}{\psi }_{\left(\left(0,i_2,\cdots ,i_{h+1}\right)|1\right)}\left(h\right),$ (28)

${\psi }_{\left(\left(1,1,i_2,\cdots ,i_{h+1}\right)|1\right)}\left(h+1\right)={\displaystyle \underset{j=0,1}{\sum }}\pi {\pi }_{|1}^{\ast \left(1|j\right)}{\psi }_{\left(\left(1,i_2,\cdots ,i_{h+1}\right)|j\right)}\left(h\right),$

${\psi }_{\left(\left(1,0,i_2,\cdots ,i_{h+1}\right)|0\right)}\left(h+1\right)=\left(1-\pi \right){\pi }_{|0}^{\ast \left(0|1\right)}{\psi }_{\left(\left(0,i_2,\cdots ,i_{h+1}\right)|1\right)}\left(h\right),$

${\psi }_{\left(\left(1,1,i_2,\cdots ,i_{h+1}\right)|0\right)}\left(h+1\right)={\displaystyle \underset{j=0,1}{\sum }}\left(1-\pi \right){\pi }_{|0}^{\ast \left(1|j\right)}{\psi }_{\left(\left(1,i_2,\cdots ,i_{h+1}\right)|j\right)}\left(h\right),$

and the recursive formulae

${\psi }_{\left(\left(0,0,i_2,\cdots ,i_{h+1}\right)|1\right)}\left(h+1\right)=\pi \left(1-{\pi }_{|1}^{\ast \left(0|1\right)}\right){\psi }_{\left(\left(0,i_2,\cdots ,i_{h+1}\right)|1\right)}\left(h\right),$

${\psi }_{\left(\left(0,1,i_2,\cdots ,i_{h+1}\right)|1\right)}\left(h+1\right)={\displaystyle \underset{j=0,1}{\sum }}\pi \left(1-{\pi }_{|1}^{\ast \left(1|j\right)}\right){\psi }_{\left(\left(1,i_2,\cdots ,i_{h+1}\right)|j\right)}\left(h\right),$

$\begin{array}{l}{\psi }_{\left(\left(0,0,i_2,\cdots ,i_{h+1}\right)|0\right)}\left(h+1\right)\\ =\left(1-\pi \right){\psi }_{((0,i_2,\dots ,i_{h+1})|0)}\left(h\right)+\left(1-\pi \right)\left(1-{\pi }_{|0}^{\ast \left(0|1\right)}\right){\psi }_{\left(\left(0,i_2,\cdots ,i_{h+1}\right)|1\right)}\left(h\right),\end{array}$

${\psi }_{\left(\left(0,1,i_2,\cdots ,i_{h+1}\right)|0\right)}\left(h+1\right)={\displaystyle \underset{j=0,1}{\sum }}\left(1-\pi \right)\left(1-{\pi }_{|0}^{\ast \left(1|j\right)}\right){\psi }_{\left(\left(1,i_2,\cdots ,i_{h+1}\right)|j\right)}\left(h\right),$ (29)

from which the probabilities

$\mathbb{P}\left(X_t=i_0,X_{t-1}=i_1,\cdots ,X_{t-\left(h+1\right)}=i_{h+1}\right)={\displaystyle \underset{j=0,1}{\sum }}{\psi }_{\left(\left(i_0,i_1,\dots ,i_{h+1}\right)|j\right)}\left(h+1\right)$

can be computed for $i_0,i_1,\cdots ,i_{h+1}=0,1$.

To generate $X_t$, $t=1,\cdots ,T$, it has been used that $X_t=Y_t+{\displaystyle {\sum }_{i\in \mathbb{N}}}\left({\displaystyle {\prod }_{j=0}^{i-1}}{\left(DY\right)}_{t-j}\right)Y_{t-i}$, $Y_t:=I_t^{\left(0|1\right)}{I}_{t-1}+I_t\left(1-I_{t-1}\right)$, ${\left(DY\right)}_t\mathrm{<mo>\; :\; </mo>}=\left(I_t^{\left(\mathrm{1\; <mo>\; |\; </mo>1}\right)}-I_t^{\left(\mathrm{0\; <mo>\; |\; </mo>1}\right)}\right)I_{t-1}+\left(I_t^{\left(\mathrm{1\; <mo>\; |\; </mo>0}\right)}-I_t\right)\left(1-I_{t-1}\right)$, with the sum being held up to a finite constant $c=20$ for $t=1$, and so on. The two presented conditions for causality and invertibility write together that

$\pi +\mathbb{E}\left|d{f}_1\left(I_t^{\left(\mathrm{0\; <mo>\; |\; </mo>1}\right)}\right)\right|\mathrm{,}\mathbb{E}\left|d{f}_1\left(I_t^{\left(\mathrm{1\; <mo>\; |\; </mo>0}\right)}\right)\right|+\mathbb{E}\left|d^2{f}_1\left(I_t^{\left(\mathrm{1\; <mo>\; |\; </mo>1}\right)}\right)\right|\mathrm{,}$

$\underset{i,j=0,1}{\max }\left\{{\pi }^{\left(i|j\right)}\right\}+\mathbb{E}\left|d{f}_1\left(I_t^{\left(1|0\right)}\right)\right|,\mathbb{E}\left|d{f}_1\left(I_t^{\left(0|1\right)}\right)\right|+\mathbb{E}\left|d^2{f}_1\left(I_t^{\left(1|1\right)}\right)\right|<1$

(the $d^h{f}_1\left(\mathrm{..}\right)$, $h=1,2$ notation was re-introduced in the very beginning of Section 3). Hence by imposing the strong requirement ${\max }_{i,j=0,1}\left\{{\pi }^{\left(i|j\right)}\right\}<1/6$, the series is secured to be where it should be (without worrying about the interdependence scenario).

Firstly the case $\pi =0.165$, ${\pi }^{\left(1|0\right)}=0.152$, ${\pi }^{\left(0|1\right)}=0.125$ and ${\pi }^{\left(1|1\right)}=0.108$ (and under interindependence) has been studied, for a number of observations $T=6,10,14$ (the maximum sample size 14 has been picked due to a purely technical reason of restoring in one vector $2^{T+1}$ probabilities for the last recursion, or otherwise a modification must take place in the code). A search of the four parameters $0<{\pi }^{\left(\mathrm{...}\right)}<0.1667$, ${\pi }^{\left(\mathrm{...}\right)}={\pi }^{\left(\mathrm{...}\right)}+0.01$ has materialized and for each attempted value the probabilities (25), (26), (27) and recursions (28)...(29) must be computed: the likelihood is none other than one of these final probabilities indicated by the position of zero-one(s) in the sample series (this is the same position regardless of the parameter values under search). Then the attempted parameter values that offered the biggest of those probabilities are taken to be the maximum likelihood estimates.

Meanwhile, an opportunity of a comparison with the TAR model (of order $p=2$ ) with same parameters ( $\pi =0.165$, ${\pi }^{\left(10\right)}=0.152$, ${\pi }^{\left(01\right)}=0.125$ and ${\pi }^{\left(11\right)}=0.108$ ) has arisen, i.e., $\left\{X_1\mathrm{,}\cdots \mathrm{,}{X}_T\right\}$ have now been born from

$\begin{array}{c}{X}_t=I_t\left(1-X_{t-1}\right)\left(1-X_{t-2}\right)+I_t^{\left(10\right)}{X}_{t-1}\left(1-X_{t-2}\right)\\ +I_t^{\left(01\right)}\left(1-X_{t-1}\right)X_{t-2}+I_t^{\left(11\right)}{X}_{t-1}{X}_{t-2}.\end{array}$

Note that it can be derived $X_t=I_t+{\displaystyle {\sum }_{n\in \mathbb{N}}}{C}_{t,n}^{\left(n-1\right)}{I}_{t-n}$ with $C_{t,1}^{\left(0\right)}:=D{I}_t^{\left(10\right)}\equiv I_t^{\left(10\right)}-I_t$, $C_{t,2}^{\left(0\right)}:=D{I}_t^{\left(01\right)}\equiv I_t^{\left(01\right)}-I_t$, $C_{t,(1,2)}^{\left(0\right)}:=D^2{I}_t^{\left(11\right)}\equiv I_t^{\left(11\right)}-I_t^{\left(10\right)}-I_t^{\left(01\right)}+I_t$ and for $n\in \mathbb{N}$

$\begin{array}{l}{C}_{t,n+1}^{\left(n\right)}:=C_{t,n+1}^{\left(n-1\right)}+C_{t,n}^{\left(n-1\right)}D{I}_{t-n}^{\left(10\right)}+C_{t,(n,n+1)}^{\left(n-1\right)}{I}_{t-n}^{\left(10\right)},C_{t,n+2}^{\left(n\right)}:=C_{t,n}^{\left(n-1\right)}D{I}_{t-n}^{\left(01\right)},\text{and}\\ C_{t,\left(n+1,n+2\right)}^{\left(n\right)}:=C_{t,n}^{\left(n-1\right)}{D}^2{I}_{t-n}^{\left(11\right)}+C_{t,\left(n,n+1\right)}^{\left(n-1\right)}{D}^{*}{I}_{t-n}^{\left(11\right)},D^{*}{I}_t^{\left(11\right)}:=I_t^{\left(11\right)}-I_t^{\left(10\right)}.\end{array}$

Computing the probabilities $\mathbb{P}\left(X_t=i,X_{t-1}=j\right)$, $i,j=0,1$ (from the same methodology and under causality) is straightforward, while the condition for causality now requires $\left[{\left(1+{\beta }^{*}{\displaystyle {\sum }_{m=0}^1}{\displaystyle {\sum }_{m_1=0}^m}1\right)}^2-1\right]<1$ or ${\beta }^{*}<\frac{\sqrt{2}-1}{3}\approx 0.138$, where

${\beta }^{*}:=\max \left\{\pi ,\max \left\{\mathbb{E}\left|D{I}_t^{\left(10\right)}\right|,\mathbb{E}\left|D{I}_t^{\left(01\right)}\right|\right\},\mathbb{E}\left|D^2{I}_t^{\left(11\right)}\right|\right\}:$

so it might be that it is violated (the search is no different than for the TARMA above).

According to Table 1, it is indeed remarkable how the performances for the two models almost coincide. That is such an encouraging sign for the TARMA, as it tempts the researcher to dare include moving-average parts and benefit from the rare flexibility, if they are willing to take the extra computational burden of the recursions: exact likelihoods can be computed as functions of the parameters and, even for small sample sizes, the estimators seem to perform quite satisfactorily. It is highlighted that the precision here is as returned by the code with no formal justification attached to it (for a smaller $R=100$ repetitions, the differences to the table were still minimal). From Table 1 it looks like there is a better bias (and MSE) performance of the $\pi$ estimator over the other three ones, so this is investigated further.

Table 2 picture still favors the $\pi$ estimator with a low absolute bias: the two ${\pi }^{\left(0|1\right)}$, ${\pi }^{\left(1|1\right)}$ estimators score the highest bias, naturally resulting from the real values being quite high. Similarly, for both numerical cases the ${\pi }^{\left(1|0\right)}$ estimator, with the lowest real value, escapes the high bias values without outperforming the $\pi$ in the same department though. The conclusions for the variance can be reversed. At this stage, it is highlighted that the easily obtainable (conditional

Table 1. Approximate (from $R=1000$ replications) bias, variance and Mean Squared Error of the Maximum Likelihood Estimators based on T consecutive observations from the TARMA (1, 1) and TAR (2) models ( $k=1$, Bernoulli variables) with real values $\pi =0.165$, ${\pi }^{\left(1|0\right)}=0.152$, ${\pi }^{\left(0|1\right)}=0.125$ and ${\pi }^{\left(1|1\right)}=0.108$ (under interindependence).

Table 2. Approximate (from $R=10000$ replications) bias, variance and Mean Squared Error of the Maximum Likelihood Estimators based on $T=6$ consecutive observations from the TARMA (1, 1) model ( $k=1$, Bernoulli variables) with real values, for the “Case 1.1”, $\pi =0.105$, ${\pi }^{\left(1|0\right)}=0.081$, ${\pi }^{\left(0|1\right)}=0.148$, ${\pi }^{\left(1|1\right)}=0.152$, or for the “Case 1.2”, $\pi =0.165$, ${\pi }^{\left(1|0\right)}=0.081$, ${\pi }^{\left(0|1\right)}=0.152$ and ${\pi }^{\left(1|1\right)}=0.148$ (both under interindependence).

likelihood) asymptotic result (for the estimators ${\hat{\pi }}^{\left(1\right)}$, $\hat{\pi }$ of ${\pi }^{\left(1\right)}:=\mathbb{P}\left(X_t=1|X_{t-1}=1\right)$, $\pi :=\mathbb{P}\left(X_t=1|X_{t-1}=0\right)$, respectively) regarding the $\left\{X_t\right\}~\text{TAR}\left(1\right)$ model, i.e.,

$T^{1/2}\left(\begin{array}{c}{\hat{\pi }}^{\left(1\right)}-{\pi }^{\left(1\right)}\\ \hat{\pi }-\pi \end{array}\right)\stackrel{D}{\to }N\left(\left(\begin{array}{c}0\\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{cc}{\pi }^{\left(1\right)}\left(1-{\pi }^{\left(1\right)}\right)/\mathbb{E}\left\{X_t\right\}& 0\\ 0& \pi \left(1-\pi \right)/\mathbb{E}\left\{1-X_t\right\}\end{array}\right)\right)$

has also been approved by simulation indications (for sample sizes $T=10,20,30$ and $R=100$ ), but those will not be presented here.

The investigation is not over yet, as it is desirable to smooth out the differences in the bias of estimators for the TARMA (1, 1) model, and the solution to that could be an “agreement” between the different $I^{\left(\mathrm{...}\right)}$ variables (at the same time t) via the presence of interdependence. Indeed according to Table 3 and for the only small sample size ( $T=6$ ), “Case 2.1” exhibits low biases in all but the -0.11726 approximation from the 1000 estimates of ${\pi }_{\mathrm{<mo>\; |\; </mo>1}}^{\ast \left(0|1\right)}$ ; similarly, for the “Case 2.2” and excluding the highest (absolute) bias for $\pi$ (its true value has been set high), almost all other estimates are approximated to be near their real value. The variance results for both cases are very impressive too. Note that the conditions for causality-invertibility might not have been followed faithfully in Table 3 (hence the real value $\pi =0.25$ in “Case 2.2”) and the “agreement” search has been set to ≥0.85 for ${\pi }_{\mathrm{<mo>\; |\; </mo>1}}^{\ast \left(\mathrm{...}\right)}$, or ≤0.15 for ${\pi }_{\mathrm{<mo>\; |\; </mo>0}}^{\ast \left(\mathrm{...}\right)}$ (together with a search $\pi \le 0.1667$ that justifies the -0.11655 bias in “Case 2.2”).

Again it is reminded that the results obtained from fewer simulations (Table 2: $R=100,1000$, Table 3: $R=100$ ) are hardly any different. The simulation indications of this section favor the addition of moving-average parts for the modelling of strict stationarity, if one is willing to slightly inconvenience themselves with a more sophisticated code for the likelihood computation. There

Table 3. Approximate (from $R=1000$ replications) bias, variance and Mean Squared Error of the Maximum Likelihood Estimators based on $T=6$ consecutive observations from the TARMA (1, 1) model ( $k=1$, Bernoulli variables) with real values, for the “Case 2.1”, $\pi =0.165$, ${\pi }_{\mathrm{<mo>\; |\; </mo>1}}^{\ast \left(1|0\right)}=0.95$, ${\pi }_{\mathrm{<mo>\; |\; </mo>0}}^{\ast \left(1|0\right)}=0.03$, ${\pi }_{\mathrm{<mo>\; |\; </mo>1}}^{\ast \left(0|1\right)}=0.98$, ${\pi }_{\mathrm{<mo>\; |\; </mo>0}}^{\ast \left(0|1\right)}=0.05$, ${\pi }_{\mathrm{<mo>\; |\; </mo>1}}^{\ast \left(1|1\right)}=0.966$, ${\pi }_{\mathrm{<mo>\; |\; </mo>0}}^{\ast \left(1|1\right)}=0.044$, or for the “Case 2.2”, $\pi =0.25$, ${\pi }_{\mathrm{<mo>\; |\; </mo>1}}^{\ast \left(1|0\right)}=0.95$, ${\pi }_{\mathrm{<mo>\; |\; </mo>0}}^{\ast \left(1|0\right)}=0.03$, ${\pi }_{\mathrm{<mo>\; |\; </mo>1}}^{\ast \left(0|1\right)}=0.98$, ${\pi }_{\mathrm{<mo>\; |\; </mo>0}}^{\ast \left(0|1\right)}=0.05$ and ${\pi }_{\mathrm{<mo>\; |\; </mo>1}}^{\ast \left(1|1\right)}=0.966$, ${\pi }_{\mathrm{<mo>\; |\; </mo>0}}^{\ast \left(1|1\right)}=0.044$ (the interdependence is implied for both).

have been no signs of distinction between the auto-regressive and moving-average estimators’ merits: not only does this strengthen the views of Sections 4.1 and 4.2, but also it brings to mind the classic Gaussian ARMA gem of [8], that the moving-average estimation transforms in theory to an auto-regression like situation.

6. Conclusions and Extending to the χ² Test for Stationarity

Straight from (23), Theorem 4.3 and Lemmas 4.4, 4.5, there is the important derivation

$T^{1/2}\left(\hat{\pi }-{\pi }_0\right)\stackrel{D}{\to }N\left(0\mathrm{,}\text{Var}{\left(\frac{1}{p_t\left({\pi }_0\right)}{\delta }_t\right)}^{-1}\right)\text{as}T\to \infty$ (30)

for the Maximum Likelihood Estimators of the parameters of a process that assigns its appearance to a causal and invertible TARMA model. The normal distribution convergence (30) implies the chi-square distribution convergence

$T{\left(\hat{\pi }-{\pi }_0\right)}^{\tau }\text{Var}\left(\frac{1}{p_t\left({\pi }_0\right)}{\delta }_t\right)\left(\hat{\pi }-{\pi }_0\right)\stackrel{D}{\to }{\chi }_{df}^2\mathrm{,}$ (31)

where $df$ equals the size of the parameter vector.

The asymptotic result (31) can be a real asset when testing whether a sample series has been generated by a stationary ( $\mathbb{Z}$ -indexed) process: the null hypothesis will be of the form

H₀: {X_t} is stationary with “distribution”...

where “distribution” is being used in the wide sense of a hyper-model that might specify as much information is allowed for the researcher (for example, marginal distribution, pairwise distribution,... or marginal and conditional distributions,... with or without knowledge of parameters etc). Under H₀, a proper identification of k (this might not be necessary if the variables are discrete with “Bernoulli” being the best of cases), as well as p and q (some faint suggestions for setting $p=0$ or $q=0$ might be found in [1]) must take place to proceed. Identification issues must still be resolved in the case that the researcher is not testing an assumption, but opts for the alternative task of point and interval estimation. For either inference route taken, further research is mandatory, in order to approximate the variance matrix $\text{Var}\left(p_t^{-1}{\delta }_t\right)$. The reader might wish to look for the special case of a causal TAR model, which due to its simplicity, is believed here that it should be supported as the most appropriate test for the non-parametric time series stationarity.

In a more structured context, an equivalent to the [14] ARMA methodology of identification, estimation, validation and forecasting can now become reality for the strictly stationary processes: any distribution (heavy-tailed, with extreme values, asymmetric or with levels of skewness) can be tamed by a proper categorization and then p and q need to be selected as well. This paper has contributed the complete guide for the second step with the best of all methods of estimation, by putting on the table all necessary properties of the TARMA maximum likelihood estimators.

In order to use the TARMA equation, the sample series needs to be produced by a stationary one. Any sign of persistent tendency, such as trends, seasonality or cycles forbid the use of TARMA: joint distributions (of equal “windows”) must look alike in time. For the second-order stationarity [15] give answers on how the series can be fixed to ARMA; or, there is a whole philosophy of modelling those tendencies directly. For the strict stationarity, one can still resort to (31), testing and computing whether the sample series is close enough to a hypothesis of hyper- dependence: a large value of the statistic, as compared to a chi-square value, will suggest that the hypothesis collapses. Otherwise, there are parametric/theoretical approaches that work under specific distributional assumptions.

To wrap it up, it is strongly believed that a precious link has been solidified between the “multiplicative linear” (as called in [1])—strictly stationary- times series and the non-parametric inferential statistics. The well-known ${\chi }^2$ test for the distribution of a variable that has generated a random sample, can be extended now to the ${\chi }^2$ test for the stationary-principled hyper-model, dominating a $\mathbb{Z}$ -indexed process that has generated a sample series: as in the former case, when extra (than the prespecified) parameters need to be estimated, it might be worth investigating the reduction of degrees of freedom in the statistic ${\chi }^2$ distribution. Regardless of various such minor issues that might be studied further, the main achievement of this paper is that it has dealt quite satisfactorily with a complex problem: the inference for any stationary time series that can be clothed by a causal and invertible TARMA equation.

Acknowledgments

This is for the author’s beloved father (deceased).

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Dimitriou-Fakalou, C. (2019) The Table Auto-Regressive Moving-Average Model for (Categorical) Stationary Series: Statistical Properties (Causality; from the All Random to the Conditional Random). Journal of Nonparametric Statistics, 31, 31-63. https://doi.org/10.1080/10485252.2018.1527912
[2]	Jacobs, P.A. and Lewis, A.W. (1978) Discrete Time Series Generated by Mixtures II: Asymptotic Properties. Journal of the Royal Statistical Society: Series B, 40, 222-228. https://doi.org/10.1111/j.2517-6161.1978.tb01667.x
[3]	Möller, T.A. and Wei, C.H. (2020) Generalized Discrete Autoregressive Moving-Average Models. Applied Stochastic Models in Business and Industry, 36, 641-659. https://doi.org/10.1002/asmb.2520
[4]	Lomnicki, Z.A. and Zaremba, S.K. (1955) Some Applications of Zero-One Processes. Journal of the Royal Statistical Society: Series B, 17, 243-255. https://doi.org/10.1111/j.2517-6161.1955.tb00198.x
[5]	Joe, H. (1996) Time Series Models with Univariate Margins in the Convolution-Closed Infinitely Divisible Class. Journal of Applied Probability, 33, 664-677. https://doi.org/10.1017/S0021900200100105
[6]	Kedem, B. (1980) Estimation of the Parameters in Stationary Autoregressive Processes after Hard Limiting. Journal of the American Statistical Association, 75, 146-153. https://doi.org/10.1080/01621459.1980.10477445
[7]	Azzalini, A. (1983) Maximum Likelihood Estimation of Order m for Stationary Stochastic Processes. Biometrika, 70, 381-387. https://doi.org/10.1093/biomet/70.2.381
[8]	Hannan, E.J. (1973) The Asymptotic Theory of Linear Time Series Models. Journal of Applied Probability, 10, 130-145. https://doi.org/10.1017/S0021900200042145
[9]	Cui, Y. and Lund, R.B. (2009) A New Look at Time Series of Counts. Biometrika, 96, 781-792. https://doi.org/10.1093/biomet/asp057
[10]	Roitershtein, A. and Zhong, Z. (2013) On Random Coefficient INAR(1) Processes. Science China Mathematics, 56, 177-200. https://doi.org/10.1007/s11425-012-4547-z
[11]	Davis, R.A., Fokianos, K., Holan, S.H., Joe, H., Livsey, J., Lund, R., Pipiras, V. and Ravishanker, N. (2021) Count Time Series: A Methodological Review. Journal of the American Statistical Association, 116, 1533-1547. https://doi.org/10.1080/01621459.2021.1904957
[12]	Anderson, T.W. and Goodman, L.A. (1957) Statistical Inference about Markov Chains. The Annals of Mathematical Statistics, 28, 89-110. https://doi.org/10.1214/aoms/1177707039
[13]	Klotz, J. (1973) Statistical Inference in Bernoulli Trials with Dependence. The Annals of Statistics, 1, 373-379. https://doi.org/10.1214/aos/1176342377
[14]	Box, G.E.P. and Jenkins, G.M. (1970) Time Series Analysis: Forecasting and Control. Holden-Day, San Francisco.
[15]	Brockwell, P.J. and Davis, R.A. (1991) Time Series: Theory and Methods. 2nd Edition, Springer-Verlag, New York. https://doi.org/10.1007/978-1-4419-0320-4