Regime-Switching Discrete ARMA Models for Categorical Time Series

Abstract

For the modeling of categorical time series, both nominal or ordinal time series, an extension of the basic discrete autoregressive moving-average (ARMA) models is proposed. It uses an observation-driven regime-switching mechanism, leading to the family of RS-DARMA models. After having discussed the stochastic properties of RS-DARMA models in general, we focus on the particular case of the first-order RS-DAR model. This RS-DAR$\left(1\right)$ model constitutes a parsimoniously parameterized type of Markov chain, which has an easy-to-interpret data-generating mechanism and may also handle negative forms of serial dependence. Approaches for model fitting are elaborated on, and they are illustrated by two real-data examples: the modeling of a nominal sequence from biology, and of an ordinal time series regarding cloudiness. For future research, one might use the RS-DAR$\left(1\right)$ model for constructing parsimonious advanced models, and one might adapt techniques for smoother regime transitions.

1. Introduction

Since the pioneering textbook on time series analysis by Box & Jenkins [1], this topic has attracted an immense interest in research and applications. To put it more precisely, real-valued time series (having a range consisting of real numbers or vectors) have been in the limelight of scientists and practitioners since that time. Besides many approaches for analyzing real-valued time series, also enumerable models have been developed, starting with the basic autoregressive moving-average (ARMA) models [1]. ARMA models are characterized by a linear conditional mean and an autocorrelation function (ACF) satisfying the so-called Yule–Walker equations. Their stochastic properties are well understood, but their actual potential for application is limited. Therefore, many alternatives and extensions have been developed, which cover more realistic time series patterns [2]. For example, if being concerned with a time series exhibiting sudden jumps and a piecewise behavior, then it is more appropriate to consider regime-switching models like the self-exciting threshold (SET) AR models proposed by Tong & Lim [3], also see the survey in Tong [4].

During the last decades, discrete-valued time series received more and more attention, see Weiß [5] for a recent survey. Especially count time series, i.e., quantitative time series with a range included in the set ${\mathbb{N}}_0=\{0,1,\dots \}$ of non-negative integers, have been studied intensively. A large number of models have been developed for this type of discrete data, not only integer-valued counterparts to the basic ARMA model, but also to more advanced models like the aforementioned SETAR models. The latter include the proposals by Möller [6], Möller et al. [7], Monteiro et al. [8], Thyregod et al. [9], Wang et al. [10].

In this article, we consider another type of discrete-valued time series, which is somewhat neglected in the time series literature: categorical time series $x_1,\dots ,x_T$ with $T\in \mathbb{N}=\{1,2,\dots \}$ and a qualitative range $\mathcal{S}$ consisting of a finite number of categories [5]. There are fundamental differences between quantitative scales (interval or ratio) and qualitative scales (ordinal or nominal) of measurement, see

Table 1. Some measures of dispersion and serial dependence for categorical time series.

	Nominal Range	Ordinal Range
Dispersion	Index of qualitative variation (Gini index):	Index of ordinal variation:
	$\mathrm{IQV}=\frac{m+1}{m}\left(1-{\sum }_{i=0}^m{p}_i^2\right)$,	$\mathrm{IOV}=\frac{4}{m}{\sum }_{i=0}^{m-1}{f}_i(1-f_i)$,
	Entropy:	Cumulative paired entropy:
	$\mathrm{En}=\frac{-1}{ln(m+1)}{\sum }_{i=0}^m{p}_{i}ln{p}_i$;	$\mathrm{CPE}=\frac{-1}{mln2}{\sum }_{i=0}^{m-1}\left(f_{i}ln{f}_i+(1-f_i)ln(1-f_i)\right)$;
	where $p_i=P(X=s_i)$, $\mathit{p}={(p_0,\dots ,p_m)}^{\top }$	where $f_i=P(X\le s_i)$, $\mathit{f}={(f_0,\dots ,f_{m-1})}^{\top }$
Serial dependence	Cohen’s κ:	Ordinal Cohen’s κ:
	${\displaystyle \kappa \left(h\right)=\frac{{\sum }_{i=0}^m\left(p_{ii}\left(h\right)-p_i^2\right)}{1-{\sum }_{i=0}^m{p}_i^2}}$,	${\displaystyle {\kappa }_{\mathrm{ord}}\left(h\right)=\frac{{\sum }_{i=0}^{m-1}\left(f_{ii}\left(h\right)-f_i^2\right)}{{\sum }_{i=0}^{m-1}{f}_i(1-f_i)}}$,
	where	where
	$p_{ij}\left(h\right)=P(X_t=s_i,X_{t-h}=s_j)$ for time lag $h\in \mathbb{N}$	$f_{ij}\left(h\right)=P(X_t\le s_i,X_{t-h}\le s_j)$ for time lag $h\in \mathbb{N}$

in Stevens [11]. In particular, quantitative scales (such as the aforementioned count data, or real-valued measurements such as temperatures or prices) allow us to use the basic arithmetic operations, which, however, is not permitted for qualitative scales. Qualitative (categorical) ranges $\mathcal{S}$ are further distinguished into the cases that $\mathcal{S}$ exhibits a natural order among the categories (ordinal range), and that the categories in $\mathcal{S}$ are unordered (nominal range). In what follows, we are interested in both types of categorical time series. To simplify notations, we always assume the possible outcomes to be arranged in a certain order (either lexicographical or natural order), i.e., we denote the range as $\mathcal{S}=\{s_0,s_1,\dots ,s_m\}$ with some $m\in \mathbb{N}$. Categorical time series require a tailor-made treatment in any sense. All the moment-based tools developed for quantitative time series cannot be applied due to the inadmissibility of the basic arithmetic operations. Instead, we have to express the location in terms of the mode (for ordinal time series, also the median can be used), and to use one of the customized measures of dispersion and serial dependence from Table 1, see Klein et al. [12], Weiß [13], [14] for further details. Even the basic time series plot can only be done in the ordinal case, whereas we may use the rate evolution graph as a substitute for nominal time series [5].

Finally, the selection of possible models for categorical time series is yet limited, see Weiß [5]. The quite flexible, higher-order Markov models suffer from a huge number of model parameters, which increases exponentially in the model order and polynomially in the number of categories, $m+1$. The latter might be reduced by amalgamating some categories, but this causes a loss of information. The extremely parsimonious discrete ARMA models by Jacobs & Lewis [15], in contrast, have a rather narrow scope of application, see the discussion in Section 2 for further details. Therefore, as a compromise between flexibility and parsimony, we extend the discrete ARMA models for categorical time series by an observation-driven regime-switching mechanism, see Section 3. For ordinal time series, it might be implemented in analogy to the SET approach for quantitative time series. However, the regime-switching can also be applied to nominal time series to capture, e.g., structures or similarities within the categorical range. As an important special case, we obtain a family of parsimonious Markov chain (MC) models, see Section 4. The application potential of the new model family is demonstrated in Section 5 with two real-data applications, where we also illustrate how model fitting might be done. Finally, we conclude in Section 6 and outline issues for future research.

2. About Discrete ARMA Models

Jacobs & Lewis [15] proposed two families of discrete ARMA models. The first of these families (labeled by the acronym “DARMA”) is defined in a nested way and has therefore found less attention in the literature. The second family (labeled as “NDARMA”), in contrast, directly imitates the ordinary ARMA recursion and has been considered in several subsequent works. Both families agree in their boundary cases (i.e., pure AR- and pure MA-type models). We concentrate on the “NDARMA family” in the sequel, and we refer to these models simply as discrete ARMA models. The regime-switching mechanism to be proposed in Section 3 could be applied to the “DARMA family” in an analogous way. According to Weiß & Göb [16], the discrete ARMA models (“NDARMA”) by Jacobs & Lewis [15] might be defined as follows (Jacobs & Lewis [15] provide an equivalent definition based on backshift operators).

Definition 1.

Let ${\left(X_t\right)}_{\mathbb{Z}}$ and ${\left({ϵ}_t\right)}_{\mathbb{Z}}$ be categorical processes with range $\mathcal{S}$, where ${\left({ϵ}_t\right)}_{\mathbb{Z}}$ is independent and identically distributed (i. i. d.) with marginal distribution p (${ϵ}_t\sim \mathit{p}$), and where ${ϵ}_t$ is independent of ${\left(X_s\right)}_{s<t}$. Let

$${\mathit{D}}_t=({\alpha }_{t,1},\dots ,{\alpha }_{t,\mathrm{p}},{\beta }_{t,0},\dots ,{\beta }_{t,\mathrm{q}})\sim \mathrm{Mult}(1;\mathit{\rho })$$

be i. i. d. multinomial random vectors with $\mathit{\rho }=({\varphi }_1,\dots ,{\varphi }_{\mathrm{p}},{\phi }_0,\dots ,{\phi }_{\mathrm{q}})$, which are independent of ${\left({ϵ}_t\right)}_{\mathbb{Z}}$ and of ${\left(X_s\right)}_{s<t}$. So the probabilities in $\mathit{\rho }$, ${\varphi }_1,\dots ,{\phi }_{\mathrm{q}}\in (0;1)$, sum up to one.

Then, ${\left(X_t\right)}_{\mathbb{Z}}$ is said to be a discrete ARMA$(\mathrm{p},\mathrm{q})$ process if it follows the recursion

$$X_t={\alpha }_{t,1}·X_{t-1}+\dots +{\alpha }_{t,\mathrm{p}}·X_{t-\mathrm{p}}+{\beta }_{t,0}·{ϵ}_t+\dots +{\beta }_{t,\mathrm{q}}·{ϵ}_{t-\mathrm{q}}.$$

(1)

(Here, if the range $\mathcal{S}$ is not numerically coded, then we assume $0·s=0$, $1·s=s$ and $s+0=s$ for each $s\in \mathcal{S}$.)

The boundary cases $\mathrm{q}=0$ and $\mathrm{p}=0$ are referred to as a DAR$\left(\mathrm{p}\right)$ process and DMA$\left(\mathrm{q}\right)$ process, respectively. The probability vector p is contained in the $(m+1)$-part unit simplex (recall that the range $\mathcal{S}$ consists of $m+1$ categories),

$${\mathbb{S}}_{m+1}=\left\{\mathit{u}\in {(0;1)}^{m+1}|u_0+\dots +u_m=1\right\},$$

and leads to m model parameters. Analogously, $\mathit{\rho }\in {\mathbb{S}}_{\mathrm{p}+\mathrm{q}+1}$ leads to $\mathrm{p}+\mathrm{q}$ further model parameters. Based on a Markov-chain representation of the discrete ARMA$(\mathrm{p},\mathrm{q})$ process according to Definition 1, Weiß [17] concluded on the ergodicity of the process as well as on the existence of a unique stationary distribution. The initial distribution for achieving stationarity is obtained by solving the invariance equation corresponding to the Markov-chain representation.

Although being denoted in an “ARMA style”, the model recursion of Equation (1) implies that $X_t$ is generated by doing nothing else than simply selecting either one of the past p observations $X_{t-1},\dots ,X_{t-\mathrm{p}}$, or one of the available $\mathrm{q}+1$ innovations ${ϵ}_t,\dots ,{ϵ}_{t-\mathrm{q}}$. As a consequence, $X_t$ and ${ϵ}_t$ have the same stationary marginal distribution, namely p, i.e., $P(X_t=s_i)=p_i=P({ϵ}_t=s_i)$ for all $i=0,\dots ,m$. Furthermore, the random-selection mechanism leads to the following transition probabilities:

$$\begin{array}{c}P(X_t=i_0|X_{t-1}=i_1,\dots X_{t-\mathrm{p}}=i_{\mathrm{p}},{ϵ}_t=j_0,{ϵ}_{t-1}=j_1,\dots ,{ϵ}_{t-\mathrm{q}}=j_{\mathrm{q}})\hfill \\ =\sum _{r=1}^{\mathrm{p}}{\delta }_{i_0{i}_r}{\varphi }_r+{\delta }_{i_0{j}_0}{\phi }_0+\sum _{r=1}^{\mathrm{q}}{\delta }_{i_0{j}_r}{\phi }_r\hfill \\ =:p_{(\mathrm{p},\mathrm{q})}\left(i_0\right|i_1,\dots ,i_{\mathrm{p}},j_0,\dots ,j_{\mathrm{q}};\mathit{\rho }),\hfill \end{array}$$

(2)

where empty sums (in case of $\mathrm{q}=0$ or $\mathrm{p}=0$) are assumed to take the value 0. Here, ${\delta }_{i,j}$ denotes the Kronecker delta, which takes the value 1 (0) iff $i=j$ ($i\ne j$).

Despite their quite extraordinary data-generating mechanism, the discrete ARMA processes according to Definition 1 have an ARMA-like serial dependence structure. If serial dependence is expressed in terms of Cohen’s $\kappa$ from Table 1, then the following Yule–Walker equations hold [16]:

$$\kappa \left(h\right)={\sum }_{j=1}^{\mathrm{p}}{\varphi }_j\kappa \left(\right|h-j\left|\right)+{\sum }_{i=0}^{\mathrm{q}-h}{\phi }_{i+h}r\left(i\right)\mathrm{for}h\ge 1,$$

(3)

where the $r\left(i\right)$ are given by $r\left(i\right)={\sum }_{j=max\{0,i-\mathrm{p}\}}^{i-1}{\varphi }_{i-j}·r\left(j\right)+{\phi }_i\mathbb{1}(0\le i\le \mathrm{q})$. Here, $\mathbb{1}\left(A\right)$ denotes the indicator function, which takes the value 1 (0) iff A is true (false). Furthermore, Weiß [17] showed that the discrete ARMA processes are $\varphi$-mixing with exponentially decreasing weights. As a consequence, one can apply the central limit theorem on p. 200 in Billingsley [18] to establish the asymptotic normality for statistics derived from such a process.

Remark 1.

Equation (3) also applies to the ordinal version of Cohen’s κ in Table 1. The reason for this is given by the fact that the bivariate probabilities at lag h, $p_{ij}\left(h\right)$, satisfy $p_{ij}\left(h\right)=\left(1-\kappa \left(h\right)\right)p_i{p}_j+\kappa \left(h\right){\delta }_{i,j}{p}_j$ for discrete ARMA processes. As a result, one obtains $f_{ij}\left(h\right)=\left(1-\kappa \left(h\right)\right)f_i{f}_j+\kappa \left(h\right)f_{min\{i,j\}}$ in the ordinal case, so $f_{ii}\left(h\right)-f_i^2=\kappa \left(h\right)f_i(1-f_i)$. Thus, for discrete ARMA processes, the identity $\kappa \left(h\right)={\kappa }_{\mathrm{ord}}\left(h\right)$ has to hold, whereas these measures usually differ from each other for other data-generating processes. It should be pointed out that further identities with other measures of serial dependence exist, see Weiß [13], Weiß & Göb [16] for details.

While the discrete ARMA models are very attractive in view of parameter parsimony (only $m+\mathrm{p}+\mathrm{q}$ parameters) and some of its model properties (e.g., Yule–Walker equations for $\kappa \left(h\right)$), they suffer from the fact that only positive forms of serial dependence are possible (in the sense that always $\kappa \left(h\right)\ge 0$). Furthermore, because of the simple selection mechanism in Equation (1), the sample paths generated by discrete ARMA models are characterized by long constant segments being finished by abrupt changes, which will often be inappropriate for real applications. For quantitative time series, a possible remedy is to use additional variation operators, see Möller & Weiß [19]. However, this solution cannot be applied to qualitative time series. Therefore, in Section 3, the novel regime-switching discrete ARMA models are proposed to offer solutions to both of the above drawbacks, the limitation to positive dependence and the piecewise constant sample paths.

3. Regime-Switching Discrete ARMA Models

Let ${\left(X_t\right)}_{\mathbb{Z}}$ be a categorical process with range $\mathcal{S}=\{s_0,s_1,\dots ,s_m\}$. If the range $\mathcal{S}$ is ordinal, then the states are strictly ordered imposing a distinct structure on $\mathcal{S}$. Even for a nominal range, the states are not necessarily free of any relations, but similarities might exist within it. An example is given by biological sequences such as deoxyribonucleic acid (DNA) and protein sequences. The four DNA bases ‘a’, ‘c’, ‘g’, and ‘t’ (adenine, cytosine, guanine, and thymine, respectively) are divided into the group of pyrimidines (c, t) and purines (a, g). In this sense, c is more similar to t than to a or g. The twenty different amino acids exhibit an even more refined similarity structure, which is visualized by the Venn diagram in Figure 1, see Taylor [20] for further details. If developing a stochastic model for a biological sequence (also see Section 5.1 below), it is reasonable to try to account for the apparent structure within the range.

As a possible solution for this task, let us now introduce a novel regime-switching (RS) extension of the discrete ARMA model. It is defined with respect to a partition of the range $\mathcal{S}=\{s_0,s_1,\dots ,s_m\}$ into K non-empty subsets, where $1\le K\le m+1$. Here, ${\mathcal{S}}_1,\dots ,{\mathcal{S}}_K$ constitute a partition of $\mathcal{S}$ iff these sets are pairwise disjoint and satisfy $\mathcal{S}={\mathcal{S}}_1\cup \dots \cup {\mathcal{S}}_K$. The K regimes shall be used as a means to account for a structure within the categorical range, such as the grouping of the DNA bases in pyrimidines and purines, or the ordering of the categories in the case of an ordinal time series. The current regime is determined by the last (or even earlier) observation: if $X_{t-1}\in {\mathcal{S}}_k$, then the process is in the kth regime at time t, and the upcoming observation $X_t$ is generated according to a regime-specific model. So we are concerned with an observation-driven (“self-exciting”) RS-mechanism, which is in contrast to, e.g., the Hidden-Markov model [21], where the regimes are defined by a latent process. As an example, for the ordinal cloudiness time series to be discussed in Section 5.2, we consider (among others) a two-regime model (so $K=2$), where the “lower regime” refers to a sky with at most scattered clouds, and the “upper regime” to broken clouds or an even overcast sky. Being in the lower regime, the upcoming cloudiness state follows a different model as if being in the upper regime.

Definition 2.

Like in Definition 1, let ${\left(X_t\right)}_{\mathbb{Z}}$ and ${\left({ϵ}_t\right)}_{\mathbb{Z}}$ be categorical processes, where the range $\mathcal{S}$ is partitioned into ${\mathcal{S}}_1,\dots ,{\mathcal{S}}_K$, and let $D_t=({\alpha }_{t,1},\dots ,{\beta }_{t,\mathrm{q}})$ denote the multinomial mixture vectors.

Let ${\mathit{p}}_{ϵ}^{\left(1\right)},\dots ,{\mathit{p}}_{ϵ}^{\left(K\right)}\in {\mathbb{S}}_{m+1}$ and let ${\mathit{\rho }}^{\left(1\right)},\dots ,{\mathit{\rho }}^{\left(K\right)}\in {\mathbb{S}}_{\mathrm{p}+\mathrm{q}+1}$ be K state-dependent probability vectors. Then, the regime-switching discrete ARMA$(\mathrm{p},\mathrm{q})$ process (“RS-DARMA”) is defined by the recursive scheme

$$\begin{array}{cc}\hfill X_t=& {\alpha }_{t,1}·X_{t-1}+\dots +{\alpha }_{t,\mathrm{p}}·X_{t-\mathrm{p}}+{\beta }_{t,0}·{ϵ}_t+\dots +{\beta }_{t,\mathrm{q}}·{ϵ}_{t-\mathrm{q}}\hfill \\ & with{\mathit{D}}_t\sim \mathrm{Mult}(1,{\mathit{\rho }}_t)\mathit{and}{ϵ}_t\sim {\mathit{p}}_t,\hfill \\ & where{\mathit{\rho }}_t={\mathit{\rho }}^{\left(k\right)}\mathit{and}{\mathit{p}}_t={\mathit{p}}_{ϵ}^{\left(k\right)}iff{X}_{t-1}\in {\mathcal{S}}_k.\hfill \end{array}$$

(4)

Note that the boundary case $K=1$ leads to the ordinary discrete ARMA model. In Definition 2, we stated the most basic type of RS-condition, where the last observation $X_{t-1}$ determines the current regime. Certainly, one may use other conditions as well, e.g., based on more delayed observations $X_{t-d}$ with $d>1$.

If $\mathcal{S}={\mathcal{S}}_1\cup \dots \cup {\mathcal{S}}_K$ denotes a partition, then each $s\in \mathcal{S}$ belongs to exactly one of the ${\mathcal{S}}_1,\dots ,{\mathcal{S}}_K$. To simplify notations, we introduce the mapping $\pi :\mathcal{S}\to \{1,\dots ,K\}$ defined by $\pi \left(s\right)=k$ iff $s\in {\mathcal{S}}_k$. So the element $s\in \mathcal{S}$ belongs to the $\pi \left(s\right)$-th subset. Then, the RS-condition in Equation (4) can be rewritten as ${\mathit{p}}_t={\mathit{p}}_{ϵ}^{\left(\pi \left(X_{t-1}\right)\right)}$ and ${\mathit{\rho }}_t={\mathit{\rho }}^{\left(\pi \left(X_{t-1}\right)\right)}$.

Note that the number of possible partitions of a set of size $m+1$ into K subsets is equal to $S(m+1,K)$, a Stirling number of the second kind [22], p. 5. These numbers can be computed recursively according to $S(m+1,K)=K·S(m,K)+S(m,K-1)$ with $S(0,0)=1$ and $S(m,0)=0=S(0,m)$. The total number of all partitions of $\mathcal{S}$ into non-empty subsets, $B_{m+1}={\sum }_{K=1}^{m+1}S(m+1,K)$, is the $(m+1)$-th Bell number [22], p. 5.

Example 1.

Recall the DNA example mentioned in the beginning of Section 3. There, the range consists of $m+1=4$ different states, $\mathcal{S}=\{\mathtt{a},\mathtt{c},\mathtt{g},\mathtt{t}\}$. Since $m=3$, there are $B_4=15$ possible partitions, namely $1,7,6,1$ different partitions into $K=1,2,3,4$ subsets, respectively. One of them splits into the purines and pyrimidines, i.e., $\mathcal{S}$ is partitioned into ${\mathcal{S}}_1=\{\mathtt{a},\mathtt{g}\}$ and ${\mathcal{S}}_2=\{\mathtt{c},\mathtt{t}\}$. In this particular case, π maps $\mathtt{a}\mapsto 1$, $\mathtt{c}\mapsto 2$, $\mathtt{g}\mapsto 1$, and $\mathtt{t}\mapsto 2$.

Let us now discuss the stochastic properties of the RS-DARMA$(\mathrm{p},\mathrm{q})$ process according to Definition 2. The ordinary discrete ARMA’s transition probabilities from Equation (2) change to

$$\begin{array}{c}P(X_t=i_0|X_{t-1}=i_1,\dots X_{t-\mathrm{p}}=i_{\mathrm{p}},{ϵ}_t=j_0,{ϵ}_{t-1}=j_1,\dots ,{ϵ}_{t-\mathrm{q}}=j_{\mathrm{q}})\hfill \\ =p_{(\mathrm{p},\mathrm{q})}\left(i_0|i_1,\dots ,i_{\mathrm{p}},j_0,\dots ,j_{\mathrm{q}};{\mathit{\rho }}^{\left(\pi \left(i_1\right)\right)}\right),\hfill \end{array}$$

(5)

i.e., depending on the outcome $X_{t-1}=i_1$, another set of dependence parameters ${\mathit{\rho }}^{\left(\pi \left(i_1\right)\right)}$ is plugged-in into Equation (2), otherwise, we just have the ordinary discrete-ARMA structure. Therefore, the MC-representation of the discrete ARMA process as well as the resulting proofs for existence, ergodicity, and mixing properties ($\varphi$-mixing with exponentially decreasing weights), as provided by Section 2.2 in Weiß [17], can be adapted to the RS-DARMA process, see Appendix A for details.

Example 2.

Let us consider the case $\mathrm{q}=0$, i.e., the purely autoregressive RS-DARMA model according to Definition 2. It can be understood as the direct counterpart to the popular SETAR model. The RS-DAR$\left(\mathrm{p}\right)$ process constitutes a pth-order Markov process, where the transition probabilities compute as

$$\begin{array}{c}P(X_t=s_i|X_{t-1}=s_{i_1},\dots X_{t-\mathrm{p}}=s_{i_{\mathrm{p}}})\hfill \\ =\sum _{j=0}^{m}P(X_t=s_i|X_{t-1}=s_{i_1},\dots X_{t-\mathrm{p}}=s_{i_{\mathrm{p}}},{ϵ}_t=s_j)P({ϵ}_t=s_j|X_{t-1}=s_{i_1})\hfill \\ \stackrel{\left(5\right)}{=}\sum _{j=0}^m{p}_{(\mathrm{p},\mathrm{q})}\left(s_i|s_{i_1},\dots ,s_{i_{\mathrm{p}}},s_j;{\mathit{\rho }}^{\left(\pi \left(s_{i_1}\right)\right)}\right)p_{ϵ,j}^{\left(\pi \left(s_{i_1}\right)\right)}\hfill \\ \stackrel{\left(2\right)}{=}\sum _{j=0}^m\left(\sum _{r=1}^{\mathrm{p}}{\delta }_{s_i{s}_{i_r}}{\varphi }_r^{\left(\pi \left(s_{i_1}\right)\right)}+{\delta }_{s_i{s}_j}{\phi }_0^{\left(\pi \left(s_{i_1}\right)\right)}\right)p_{ϵ,j}^{\left(\pi \left(s_{i_1}\right)\right)}\hfill \\ =\sum _{r=1}^{\mathrm{p}}{\delta }_{s_i{s}_{i_r}}{\varphi }_r^{\left(\pi \left(s_{i_1}\right)\right)}+{\phi }_0^{\left(\pi \left(s_{i_1}\right)\right)}{p}_{ϵ,i}^{\left(\pi \left(s_{i_1}\right)\right)}.\hfill \end{array}$$

These can now be used for likelihood computations or forecasting purposes.

4. A Class of Parsimonious Markov Chains

A particularly important special case of the RS-DARMA family is obtained by setting $(\mathrm{p},\mathrm{q})=(1,0)$, because such a RS-DAR$\left(1\right)$ process constitutes a parsimoniously parameterized Markov chain (MC), also see Example 2. MCs, in turn, are of great relevance, because (1) such a first-order memory is often sufficient in practice, and because (2) MCs may constitute the starting point for defining more complex time series models [5]. Well-known examples regarding (2) are the so-called mixture transition distribution (MTD) model proposed by Raftery [23], which extends an underlying MC to a higher-order Markov model with only one additional parameter for each increment of the model order, or the hidden-Markov model (HMM), where the observable process is controlled by a latent MC [21]. Recall that an HMM can be interpreted as a parameter-driven RS-model, whereas we consider observation-driven RS-models in this article.

For any of the extensions in the sense of (2), it would be of relevance to start with a maximally parsimonious MC model to keep the overall number of model parameters at a feasible level. A full MC model on $\mathcal{S}$ has $m(m+1)$ model parameters, which increases quadratically in m. The lower bound of model parameters is determined through the i. i. d.-case, where only the marginal distribution p has to be specified (which requires m parameters). So any non-i. i. d. model with unspecified marginal distribution must have $\ge m+1$ parameters. The ordinary DAR$\left(1\right)$ model with its $m+1$ parameters reaches this lower bound (and also the so-called “Negative Markov model”, see, e.g., Weiß [13] for details). It may sometimes be too simplistic for practice, recall the discussion in Section 2. In fact, an MTD$\left(\mathrm{p}\right)$ model relying on a DAR$\left(1\right)$-MC just leads to a DAR$\left(\mathrm{p}\right)$ model. Thus, using a (true) RS-DAR$\left(1\right)$ model as a base for defining a HMM or MTD model, respectively, might turn out as a reasonable compromise between model flexibility and parameter parsimony.

The model recursion of an ordinary DAR$\left(1\right)$ process can be denoted as

$$X_t={\alpha }_t{X}_{t-1}+(1-{\alpha }_t){ϵ}_t\mathrm{with}{\alpha }_t\sim \mathrm{Bin}(1,\varphi ),{ϵ}_t\sim {\mathit{p}}_{ϵ},$$

(6)

where $\varphi ={\varphi }_1$ and ${\mathit{p}}_{ϵ}=\mathit{p}$ according to Definition 1. Its transition probabilities equal $p_{i|j}=P(X_t=i|X_{t-1}=j)=(1-\varphi )p_{ϵ,i}+\varphi {\delta }_{i,j}$. Let $\mathcal{S}={\mathcal{S}}_1\cup \dots \cup {\mathcal{S}}_K$ be a partition, then the RS-DAR$\left(1\right)$ model is generally defined by

$$X_t={\alpha }_t{X}_{t-1}+(1-{\alpha }_t){ϵ}_t\mathrm{with}{\alpha }_t\sim \mathrm{Bin}(1,{\varphi }^{\left(\pi \left(X_{t-1}\right)\right)}),{ϵ}_t\sim {\mathit{p}}_t^{\left(\pi \left(X_{t-1}\right)\right)},$$

(7)

see Example 2. For applications, however, it might be better to impose further restrictions such that the resulting model is better interpretable. Therefore, we shall now propose two special cases of Equation (7), where the regimes affect either the marginal distribution or the dependence parameter.

4.1. Marginal Regimes

A RS-DAR$\left(1\right)$ model with respect to the marginals is defined by

$$\begin{array}{cc}\hfill X_t=& {\alpha }_t{X}_{t-1}+(1-{\alpha }_t){ϵ}_t\mathrm{with}{\alpha }_t\sim \mathrm{Bin}(1,\varphi ),{ϵ}_t\sim {\mathit{p}}_t,\hfill \\ & \mathrm{where}{\mathit{p}}_t={\mathit{p}}_{ϵ}^{\left(k\right)}\mathrm{iff}{X}_{t-1}\in {\mathcal{S}}_k,\mathrm{i}.\mathrm{e}.,{\mathit{p}}_t={\mathit{p}}_{ϵ}^{\left(\pi \left(X_{t-1}\right)\right)}.\hfill \end{array}$$

(8)

Here, ${\mathit{p}}_{ϵ}^{\left(1\right)},\dots ,{\mathit{p}}_{ϵ}^{\left(K\right)}\in {\mathbb{S}}_{m+1}$ are the K state-dependent probability vectors for the innovations ${ϵ}_t$, implying $Km+1$ parameters for model Equation (8). $K=m+1$ would lead to a full MC model, but then $\varphi$ would not be identifiable anymore. So model Equation (8) requires $K\le m$. The transition probabilities follow as $p_{i|j}=(1-\varphi )p_{ϵ,i}^{\left(\pi \right(j\left)\right)}+\varphi {\delta }_{i,j}$, see Example 2. In contrast to the ordinary DAR$\left(1\right)$ model, the RS-DAR$\left(1\right)$ model Equation (8) also allows for negative serial dependence. This can be obtained by choosing the ${\mathit{p}}_{ϵ}^{\left(k\right)}$ such that $p_{ϵ,i}^{\left(\pi \right(i\left)\right)}\to 0$ for all $i=0,\dots ,m$.

Example 3.

Let $m=3$, and define the partition ${\mathcal{S}}_1=\{s_0,s_1\}$ and ${\mathcal{S}}_2=\{s_2,s_3\}$. Furthermore, let us consider the boundary case

$${\mathit{p}}_{ϵ}^{\left(1\right)}=\left(0,0,p^{\left(1\right)},1-p^{\left(1\right)}\right){}^{\top },{\mathit{p}}_{ϵ}^{\left(2\right)}=\left(p^{\left(2\right)},1-p^{\left(2\right)},0,0\right){}^{\top },$$

where $p^{\left(1\right)},p^{\left(2\right)}\in (0;1)$. Then $p_{ϵ,i}^{\left(\pi \right(i\left)\right)}=0$ for all $i=0,\dots ,m$. The transition matrix $\mathbf{P}={\left(p_{i|j}\right)}_{i,j=0,\dots ,m}$ equals

$$\mathbf{P}=\left(\begin{array}{cccc}\varphi & 0& (1-\varphi )p^{\left(2\right)}& (1-\varphi )p^{\left(2\right)}\\ 0& \varphi & (1-\varphi )(1-p^{\left(2\right)})& (1-\varphi )(1-p^{\left(2\right)})\\ (1-\varphi )p^{\left(1\right)}& (1-\varphi )p^{\left(1\right)}& \varphi & 0\\ (1-\varphi )(1-p^{\left(1\right)})& (1-\varphi )(1-p^{\left(1\right)})& 0& \varphi \end{array}\right).$$

Solving the invariance equation $\mathbf{P}\mathit{p}=\mathit{p}$, the stationary marginal distribution, $X_t\sim \mathit{p}$, equals $\mathit{p}=\frac{1}{2}\left(p^{\left(2\right)},1-p^{\left(2\right)},p^{\left(1\right)},1-p^{\left(1\right)}\right){}^{\top }$. From the diagonal of $\mathbf{P}$, it becomes clear that the bivariate probabilities $p_{ii}\left(1\right)=\varphi p_i$. Furthermore, $2{\sum }_{i=0}^m{p}_i^2=1-p^{\left(1\right)}(1-p^{\left(1\right)})-p^{\left(2\right)}(1-p^{\left(2\right)})$, so one computes Cohen’s κ at lag 1, see Table 1, as 

$$\kappa \left(1\right)=1-\frac{2(1-\varphi )}{1+p^{\left(1\right)}(1-p^{\left(1\right)})+p^{\left(2\right)}(1-p^{\left(2\right)})}.$$

This expression might also become negative, which is illustrated in Figure 2a, where $\kappa \left(1\right)$ is plotted against ϕ and p with $p^{\left(1\right)}=p^{\left(2\right)}=p$.

Note that if the range is assumed ordinal, then ${\kappa }_{\mathrm{ord}}\left(1\right)$ differs from $\kappa \left(1\right)$ (in contrast to the case of ordinary DAR$\left(1\right)$ models). With analogous computations as before, one obtains

$${\kappa }_{\mathrm{ord}}\left(1\right)=1-\frac{2(1-\varphi )(2+p^{\left(2\right)}-p^{\left(1\right)})}{2+p^{\left(2\right)}-p^{\left(1\right)}+p^{\left(1\right)}(1-p^{\left(1\right)})+p^{\left(2\right)}(1-p^{\left(2\right)})}.$$

But like for $\kappa \left(1\right)$, also ${\kappa }_{\mathrm{ord}}\left(1\right)$ might take negative values, see Figure 2b.

4.2. Dependence Regimes

As a possible alternative to Equation (8), a RS-DAR$\left(1\right)$ model with respect to the dependence parameter is defined by

$$\begin{array}{cc}\hfill X_t=& {\alpha }_t{X}_{t-1}+(1-{\alpha }_t){ϵ}_t\mathrm{with}{\alpha }_t\sim \mathrm{Bin}(1,{\varphi }_t),{ϵ}_t\sim {\mathit{p}}_{ϵ},\hfill \\ & \mathrm{where}{\varphi }_t={\varphi }^{\left(k\right)}\mathrm{iff}{X}_{t-1}\in {\mathcal{S}}_k,\mathrm{i}.\mathrm{e}.,{\varphi }_t={\varphi }^{\left(\pi \left(X_{t-1}\right)\right)}.\hfill \end{array}$$

(9)

Here, ${\varphi }^{\left(1\right)},\dots ,{\varphi }^{\left(K\right)}\in (0;1)$ are the K state-dependent dependence parameters, so model Equation (9) has altogether $m+K$ parameters. The transition probabilities compute as $p_{i|j}=(1-{\varphi }^{\left(\pi \right(j\left)\right)})p_{ϵ,i}+{\varphi }^{\left(\pi \right(j\left)\right)}{\delta }_{i,j}$, see Example 2.

Example 4.

As a possible application of model Equation (9), consider an ordinal range $\mathcal{S}$ with $m\ge 2$, and let us assume individual dependence parameters for each state. More precisely, let us assume the partition $\mathcal{S}={\bigcup }_{k=0}^m{\mathcal{S}}_k$ with ${\mathcal{S}}_k=\left\{s_k\right\}$, and the corresponding $m+1$ state-dependent dependence parameters ${\varphi }^{\left(0\right)},\dots ,{\varphi }^{\left(m\right)}\in (0;1)$. So we have $2m+1$ model parameters, whereas a full MC would have $m(m+1)$ parameters. The transition probabilities equal $p_{i|j}=(1-{\varphi }^{\left(j\right)})p_{ϵ,i}+{\varphi }^{\left(j\right)}{\delta }_{i,j}$. Then $X_t\sim \mathit{p}$ with

$$p_j=\frac{\frac{p_{ϵ,j}}{1-{\varphi }^{\left(j\right)}}}{{\sum }_{l=0}^m\frac{p_{ϵ,l}}{1-{\varphi }^{\left(l\right)}}},and{p}_{ij}=p_{i|j}{p}_j=\frac{p_{ϵ,i}{p}_{ϵ,j}+{\delta }_{i,j}\frac{{\varphi }^{\left(i\right)}}{1-{\varphi }^{\left(i\right)}}{p}_{ϵ,i}}{{\sum }_{l=0}^m\frac{p_{ϵ,l}}{1-{\varphi }^{\left(l\right)}}}.$$

The formula for $p_i$ is easily verified by computing

$$\sum _{j=0}^m{p}_{ij}=\frac{p_{ϵ,i}+\frac{{\varphi }^{\left(i\right)}}{1-{\varphi }^{\left(i\right)}}{p}_{ϵ,i}}{{\sum }_{l=0}^m\frac{p_{ϵ,l}}{1-{\varphi }^{\left(l\right)}}}=\frac{\frac{p_{ϵ,i}}{1-{\varphi }^{\left(i\right)}}}{{\sum }_{l=0}^m\frac{p_{ϵ,l}}{1-{\varphi }^{\left(l\right)}}}.$$

The cumulative probabilities f of $X_t$ are given by $f_i={\sum }_{j=0}^i\frac{p_{ϵ,j}}{1-{\varphi }^{\left(j\right)}}/{\sum }_{l=0}^m\frac{p_{ϵ,l}}{1-{\varphi }^{\left(l\right)}}$. It follows that

$$p_{ii}=\frac{\frac{p_{ϵ,i}}{1-{\varphi }^{\left(i\right)}}-p_{ϵ,i}(1-p_{ϵ,i})}{{\sum }_{l=0}^m\frac{p_{ϵ,l}}{1-{\varphi }^{\left(l\right)}}},f_{ii}=\sum _{r,s=0}^i{p}_{rs}=\frac{{\sum }_{j=0}^i\frac{p_{ϵ,j}}{1-{\varphi }^{\left(j\right)}}-f_{ϵ,i}(1-f_{ϵ,i})}{{\sum }_{l=0}^m\frac{p_{ϵ,l}}{1-{\varphi }^{\left(l\right)}}},$$

which can be used to compute the dependence measures $\kappa \left(1\right)$ and ${\kappa }_{\mathrm{ord}}\left(1\right)$ from Table 1.

4.3. Statistical Inference

Let $\mathit{\theta }$ denote the vector of all model parameters, i.e., for the RS-DAR$\left(1\right)$ model Equation (8), we have $\mathit{\theta }=\left(\varphi ,p_{ϵ,1}^{\left(1\right)},\dots ,p_{ϵ,m}^{\left(K\right)}\right)\in {(0;1)}^{Km+1}$, whereas $\mathit{\theta }=\left({\varphi }^{\left(1\right)},\dots ,{\varphi }^{\left(K\right)},p_{ϵ,1},\dots ,p_{ϵ,m}\right)\in {(0;1)}^{K+m}$ for model Equation (9). To estimate $\mathit{\theta }$ from a given time series $x_1,\dots ,x_T$, we use the maximum likelihood (ML) approach. The (conditional) ML estimate $\widehat{\mathit{\theta }}$ is obtained by numerically maximizing the log-likelihood function,

$$\ell \left(\mathit{\theta }\right)={\sum }_{t=2}^{T}ln{p}_{x_t|x_{t-1}}\left(\mathit{\theta }\right).$$

(10)

We denote the maximized log-likelihood by ${\ell }_{max}=\frac{T}{T-1}\ell \left(\widehat{\mathit{\theta }}\right)$, where the factor $\frac{T}{T-1}$ corrects for the conditioning on $x_1$ [5], p. 236. The existence, consistency, and asymptotic normality are easily established by proving that Condition 5.1 in Billingsley [24] holds. This condition requires that

• the set $D=\left\{(k,l)|p_{k|l}\left(\mathit{\theta }\right)>0\right\}$ does not dependent on $\mathit{\theta }$;
• each $p_{k|l}\left(\mathit{\theta }\right)$ has continuous partial derivatives in $\theta$;
• the Jacobian matrix of ${\left(\dots ,\partial p_{k|l}\left(\mathit{\theta }\right),\dots \right)}_{(k,l)\in D}$ has full rank, i.e., the rank $n_{\mathrm{model}}=\dim \left(\mathit{\theta }\right)$;
• the MC is irreducible.

Since the transition probabilities are quadratic polynomials in the model parameters, part 2 is always satisfied. Part 1 holds by restricting the model parameters to the open interval $(0;1)$, then all $p_{k|l}\left(\mathit{\theta }\right)>0$. This also implies the irreducibility of the transition matrix. Part 3 is ensured by an appropriate design of the model (identifiability of parameters).

Finally, if the model design is not fixed by the considered application scenario, then one will commonly try out multiple types of partitioning (between the boundary cases of an ordinary DAR$\left(1\right)$ model and a full MC). In this case, the model selection might be done based on a certain type of information criterion, such as Akaike’s information criterion (AIC) or the Bayesian information criterion (BIC), see Burnham & Anderson [25]. These are given by

$$\mathrm{AIC}=-2{\ell }_{max}+2n_{\mathrm{model}},\mathrm{BIC}=-2{\ell }_{max}+n_{\mathrm{model}}lnT,$$

(11)

respectively. The performance of AIC and BIC if selecting among general Markov models was investigated by Katz [26]. It was shown that only the BIC is consistent while the AIC tends to overfitting.

5. Real-Data Applications

In what follows, we apply the RS-DAR$\left(1\right)$ models to two data examples. The first one refers to a DNA sequence, which constitutes a nominal time series (Section 5.1). The second example, in contrast, is about an ordinal time series of cloudiness states (Section 5.2).

5.1. DNA Sequence Modeling

Let us pick up the discussion in the beginning of Section 3 as well as in Example 1, where we considered a nominal DNA sequence having the range $\mathcal{S}=\{\mathtt{a},\mathtt{c},\mathtt{g},\mathtt{t}\}$ (so $m=3$). Although such a sequence does not constitute a “time” series in the original sense, it is common practice to use models for stochastic processes as a tool for summarizing its main properties, see Churchill [27], Dehnert et al. [28]. In what follows, we consider the DNA sequence of the Bovine leukemia virus (length $T=8419$), which is published by the National Center for Biotechnology Information at https://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/nuccore/NC_001414?%3Fdb=nucleotide. Its rate evolution graph (a time series plot is not possible for nominal data, see Weiß [5]) and its sample Cohen’s $\kappa$ (recall Table 1) are plotted in Figure 3.

The Bovine sequence was comprehensively analyzed in Section 7.2 of Weiß [5], where several types of models were fitted to the data. The most parsimonious model used was the ordinary DAR$\left(1\right)$ model (four parameters), but the best fit was actually obtained by a full MC model (12 parameters). Weiß [5] also tried out a 2-state HMM (eight parameters), which improved over the DAR$\left(1\right)$ model but was inferior to the full MC. Nevertheless, it is worth noting that the two hidden states that were constructed during model fitting closely matched the purines (a, g) and pyrimidines (c, t), respectively. This indicates a RS-approach might turn out to be appropriate for the data.

Therefore, we now fit several types of RS-DAR$\left(1\right)$ model to these data. As the criterion for model selection, we use the BIC, which is very suitable for time series of large length T. First, we try out the approach of Section 4.2, where the DAR$\left(1\right)$’s dependence parameter $\varphi$ depends on the current regime. The obtained BIC results are summarized in

Table 2. Bovine DNA data: Bayesian information criterion (BIC) of RS-DAR$\left(1\right)$ models with respect to dependence parameter $\varphi$, compared to those of ordinary DAR$\left(1\right)$ and Markov chain (MC) model, respectively.

Model	DAR(1)	${\mathcal{S}}_1=\{\mathtt{a},\mathtt{g}\},$	${\mathcal{S}}_1=\left\{\mathtt{a}\right\},{\mathcal{S}}_2=\left\{\mathtt{g}\right\},$	Full MC
		${\mathcal{S}}_{\mathbf{2}}=\{\mathtt{c},\mathtt{t}\}$	${\mathcal{S}}_{\mathbf{3}}=\left\{\mathtt{c}\right\},{\mathcal{S}}_{\mathbf{4}}=\left\{\mathtt{t}\right\}$
BIC	22927.4	22926.3	22928.9	22824.6

. The RS-DAR$\left(1\right)$ model with two regimes (purines ${\mathcal{S}}_1$ and pyrimidines ${\mathcal{S}}_2$) and, thus, five parameters, leads to a slight improvement compared to the DAR$\left(1\right)$ model, but the four-regime model (seven parameters) performs even worse in terms of the BIC. Any of these models is much worse than the full MC, so a regime-dependent dependence structure does not seem to be appropriate for the data.

Hence, let us now use the RS-DAR$\left(1\right)$ model with regime-dependent innovations’ distribution, see Section 4.1. Starting again with the two-regime model (purines vs. pyrimidines; seven parameters), we achieve a considerable improvement over the DAR$\left(1\right)$ model, but still do not reach the full MC’s BIC, see

Table 3. Bovine DNA data: BICs of RS-DAR$\left(1\right)$ models with respect to marginals ${\mathit{p}}_{ϵ}$, compared to those of ordinary DAR$\left(1\right)$ and MC model, respectively.

Model	DAR(1)		${\mathcal{S}}_1=\left\{\mathtt{a}\right\},$	${\mathcal{S}}_1=\{\mathtt{a},\mathtt{g}\},$	Full MC
		${\mathcal{S}}_{\mathbf{1}}=\{\mathtt{a},\mathtt{g}\},$	${\mathcal{S}}_{\mathbf{2}}=\left\{\mathtt{g}\right\},$	${\mathcal{S}}_{\mathbf{2}}=\left\{\mathtt{c}\right\},$
		${\mathcal{S}}_{\mathbf{2}}=\{\mathtt{c},\mathtt{t}\}$	${\mathcal{S}}_{\mathbf{3}}=\{\mathtt{c},\mathtt{t}\}$	${\mathcal{S}}_{\mathbf{3}}=\left\{\mathtt{t}\right\}$
BIC	22927.4	22869.6	22875.0	22822.0	22824.6

. Thus, we try splitting either the purine or the pyrimidine regime, leading to a three-regime model with ten parameters. If we split the purines into ${\mathcal{S}}_1=\left\{\mathtt{a}\right\},{\mathcal{S}}_2=\left\{\mathtt{g}\right\}$, then the BIC deteriorates again. However, if splitting the pyrimidines into ${\mathcal{S}}_2=\left\{\mathtt{c}\right\},{\mathcal{S}}_3=\left\{\mathtt{t}\right\}$, then the BIC improves and even becomes better than that of the full MC. So this type of RS-DAR$\left(1\right)$ model is best among all candidate models.

Thus, let us analyze this model fit in some more detail. The dependence parameter is estimated as $\widehat{\varphi }\approx 0.061$, and the three regime-dependent innovations’ distributions are

$$\begin{array}{c}{\widehat{\mathit{p}}}_{ϵ}^{\left(1\right)}\approx {(0.244,0.299,0.245,0.211)}^{\top },{\widehat{\mathit{p}}}_{ϵ}^{\left(2\right)}\approx {(0.216,0.352,0.143,0.289)}^{\top },\\ {\widehat{\mathit{p}}}_{ϵ}^{\left(3\right)}\approx {(0.181,0.360,0.240,0.219)}^{\top }.\end{array}$$

This implies, for example, that if $X_{t-1}\in {\mathcal{S}}_3=\left\{\mathtt{t}\right\}$, then the probability for ${ϵ}_t\in {\mathcal{S}}_2=\left\{\mathtt{c}\right\}$ is quite large, whereas $X_{t-1}\in {\mathcal{S}}_1=\{\mathtt{a},\mathtt{g}\}$ also leads to a rather large probability for ${ϵ}_t\in {\mathcal{S}}_1$. These “transition rules” can also be seen from the resulting transition matrix

$${\mathbf{P}}_{\mathrm{fit}}\approx \left(\begin{array}{cccc}0.291& 0.203& 0.229& 0.170\\ 0.281& 0.392& 0.281& 0.338\\ 0.230& 0.134& 0.292& 0.225\\ 0.198& 0.271& 0.198& 0.267\end{array}\right),$$

which implies the following stationary marginal distribution for $X_t$:

$${\mathit{p}}_{\mathrm{fit}}\approx {(0.220,0.331,0.210,0.239)}^{\top }.$$

For the given three-digit rounding, this perfectly agrees with the vector of relative frequencies computed from the data. This excellent agreement between fitted and observed marginal distribution carries over to the measures of dispersion given in Table 1, with an IQV value of $\approx 0.988$ and an entropy of $\approx 0.987$. Both values are very close to 1, because the marginal distribution is quite close to a uniform distribution, which is considered as the maximally dispersed scenario for nominal data [13].

Finally, let us do a diagnostic check of the serial dependence structure. On p. 146 in Weiß [5], it was pointed out that the sample value of Cohen’s $\kappa$ at lag 1, $\widehat{\kappa }\left(1\right)\approx 0.080$, deviates from the corresponding sample value of Cramer’s v, $\widehat{v}\left(1\right)\approx 0.113$. Here, Cramer’s v is defined by $v\left(k\right)={\left(\frac{1}{m}{\sum }_{i,j\in \mathcal{S}}{\left(p_{ij}\left(k\right)-p_i{p}_j\right)}^2/\left(p_i{p}_j\right)\right)}^{1/2}$, constituting a so-called “unsigned” measure of serial dependence [5]. This discrepancy between $\widehat{\kappa }\left(1\right)$ and $\widehat{v}\left(1\right)$ contradicts a DAR$\left(1\right)$ model, where $\kappa \left(h\right)$ and $v\left(h\right)$ exactly agree, but it was reproduced by the fitted full MC. For the fitted three-regime model, we obtain ${\kappa }_{\mathrm{fit}}\left(1\right)\approx 0.080$ and $v_{\mathrm{fit}}\left(1\right)\approx 0.111$. The fitted RS-DAR$\left(1\right)$ model reproduces the discrepancy between $\kappa$ and v, confirming its adequacy for the Bovine data.

5.2. Cloudiness Time Series

The amount of cloud coverage is measured in “okta”, i.e., in eighths of the sky being covered by clouds. Then, a common classification of cloudiness consists of the following five ordinal states (so $m=4$), ordered from lowest to highest: ‘SKC’ (sky clear, 0 oktas), ‘FEW’ (few, 1–2 oktas), ‘SCT’ (scattered, 3–4 oktas), ‘BKN’ (broken, 5–7 oktas), and  ‘OVC’ (overcast, 8 oktas). We considered a time series obtained from the “DWD Climate Data Center” offered by the Deutscher Wetterdienst (German Weather Service) at https://cdc.dwd.de/portal/201912031600/mapview. The data refer to the hourly observations of cloudiness at the weather station in Schleswig (a town in the north of Germany) in May 2011. So the time series is of total length $T=744$. A plot of the data as well as the corresponding sample ordinal Cohen’s $\kappa$ (recall Table 1) are shown in Figure 4. We have a rather strong degree of serial dependence ($\widehat{\kappa }\left(1\right)\approx 0.741$). The marginal cumulative frequencies are given by $\widehat{\mathit{f}}\approx {(0.050,0.253,0.433,0.825)}^{\top }$, leading to the dispersion values 0.626 for the sample IOV and 0.689 for the CPE (recall Table 1). Note that in the ordinal case, maximal dispersion does not go along with a uniform distribution, but with an extreme two-point distribution, i.e., $\mathit{f}={(0.5,\dots ,0.5)}^{\top }$.

In analogy to Section 5.1, we now fit types of MC model to the cloudiness data. The most parsimonious model is again the DAR$\left(1\right)$ model, but this does not account for the ordinal structure of the range. Therefore, besides the full MC model, we also considered RS-DAR$\left(1\right)$ models (with respect to the marginal distribution), which are designed such that they account for the natural order within the range (the RS-DAR$\left(1\right)$ model with $m+1$ dependence regimes, see Example 4, performs clearly worse and is, thus, not reported here). This is achieved by imitating the threshold approach for quantitative time series: we considered the partitioning ${\mathcal{S}}_1=\{\mathtt{SKC},\mathtt{FEW},\mathtt{SCT}\}$, ${\mathcal{S}}_2=\{\mathtt{BKN},\mathtt{OVC}\}$ corresponding to a threshold at $\mathtt{SCT}$, and the refinement ${\mathcal{S}}_1=\{\mathtt{SKC},\mathtt{FEW}\}$, ${\mathcal{S}}_2=\left\{\mathtt{SCT}\right\}$, ${\mathcal{S}}_3=\{\mathtt{BKN},\mathtt{OVC}\}$ corresponding to two thresholds at $\mathtt{FEW},\mathtt{SCT}$. The BICs of the four candidate models are summarized in

Table 4. Cloudiness data: BICs of RS-DAR$\left(1\right)$ models with respect to marginals ${\mathit{p}}_{ϵ}$, compared to those of ordinary DAR$\left(1\right)$ and MC model, respectively.

Model	DAR(1)		${\mathcal{S}}_1=\{\mathtt{SKC},\mathtt{FEW}\},$	Full MC
		${\mathcal{S}}_{\mathbf{1}}=\{\mathtt{SKC},\mathtt{FEW},\mathtt{SCT}\},$	${\mathcal{S}}_{\mathbf{2}}=\left\{\mathtt{SCT}\right\},$
		${\mathcal{S}}_{\mathbf{2}}=\{\mathtt{BKN},\mathtt{OVC}\}$	${\mathcal{S}}_{\mathbf{3}}=\{\mathtt{BKN},\mathtt{OVC}\}$
BIC	1423.4	1345.5	1350.1	1392.6

.

According to the BIC, we select the two-regime model, where the dependence parameter is estimated as $\widehat{\varphi }\approx 0.547$, and where the two regime-dependent innovations’ distributions are

$$\begin{array}{cc}\hfill {\widehat{\mathit{p}}}_{ϵ}^{\left(1\right)}\approx & {(0.091,0.395,0.238,0.269,0.006)}^{\top },\hfill \\ \hfill {\widehat{\mathit{p}}}_{ϵ}^{\left(2\right)}\approx & {(0.000,0.050,0.195,0.560,0.194)}^{\top }.\hfill \end{array}$$

According to this fitted model, we mainly produce innovations from the lower regime if staying in the lower regime, and vice versa. This type of “inertia” can also be seen from the right part of Figure 5, where a time series was simulated according to the fitted two-regime model. This sample path looks much more similar to the original time series in Figure 4 than the simulated DAR$\left(1\right)$ path in the left part of Figure 5, which does not exhibit a piecewise behavior. The (cumulative) stationary marginal distribution of the fitted two-regime model results as

$${\mathit{f}}_{\mathrm{fit}}\approx {(0.043,0.256,0.471,0.894)}^{\top },$$

which is reasonably close to $\widehat{\mathit{f}}$. In fact, looking at the corresponding dispersion measures, we get ${\mathrm{IOV}}_{\mathrm{fit}}\approx 0.575$ and ${\mathrm{CPE}}_{\mathrm{fit}}\approx 0.640$, which is only slightly below the above sample values. The dependence structure is captured quite well, with ${\kappa }_{\mathrm{ord},\mathrm{fit}}\left(1\right)\approx 0.697$.

Remark 2.

Since $T=744$ is not that large, one might also think of using the AIC for model selection. In this case, the three-regime model ($\mathrm{AIC}\approx 1290.1$) would be preferred over the two-regime model ($\mathrm{AIC}\approx 1304.0$). The parameter estimates are $\widehat{\varphi }\approx 0.592$ and

$$\begin{array}{cc}\hfill {\widehat{\mathit{p}}}_{ϵ}^{\left(1\right)}\approx & {(0.154,0.414,0.302,0.130,0.000)}^{\top },\hfill \\ \hfill {\widehat{\mathit{p}}}_{ϵ}^{\left(2\right)}\approx & {(0.031,0.486,0.001,0.464,0.018)}^{\top },\hfill \\ \hfill {\widehat{\mathit{p}}}_{ϵ}^{\left(3\right)}\approx & {(0.000,0.055,0.204,0.361,0.380)}^{\top }.\hfill \end{array}$$

It is interesting to note that the innovations’ distribution of the central regime, ${\mathcal{S}}_2=\left\{\mathtt{SCT}\right\}$, will hardly produce the value $\mathtt{SCT}$ itself. Instead, there is a nearly fifty-fifty chance of falling either above or below this value. The four additional model parameters are outweighed by a closer fit of the marginal distribution, now

$${\mathit{f}}_{\mathrm{fit}}\approx {(0.052,0.297,0.490,0.802)}^{\top }$$

with stronger dispersion (${\mathrm{IOV}}_{\mathrm{fit}}\approx 0.666$, ${\mathrm{CPE}}_{\mathrm{fit}}\approx 0.722$), and by more serial dependence ${\kappa }_{\mathrm{ord},\mathrm{fit}}\left(1\right)\approx 0.739$.

The fitted two-regime model can now be applied to forecasting the cloudiness. The one-step-ahead conditional mode or median is always equal to the given state, i.e., the point forecast for the next hour is equal to the current cloudiness state. But if using 95% prediction intervals for weather forecasting, then more complex rules of the form “$X_t=a\Rightarrow X_{t+1}\in B$” are obtained:

$$\begin{array}{c}\mathtt{SKC}\Rightarrow \{\mathtt{SKC},\mathtt{FEW},\mathtt{SCT},\mathtt{BKN}\},\mathtt{FEW}\mathrm{or}\mathtt{SCT}\Rightarrow \{\mathtt{FEW},\mathtt{SCT},\mathtt{BKN}\},\hfill \\ \mathtt{BKN}\mathrm{or}\mathtt{OVC}\Rightarrow \{\mathtt{SCT},\mathtt{BKN},\mathtt{OVC}\}.\hfill \end{array}$$

6. Conclusions

We extended the basic discrete ARMA model of Jacobs & Lewis [15] by an observation-driven regime-switching mechanism, leading to the family of RS-DARMA models. Particular attention was given to two instances of the RS-DAR$\left(1\right)$ model, because they constitute an easy-to-interpret type of parsimoniously parameterized MC model. Furthermore, in contrast to the ordinary DAR$\left(1\right)$ model, the RS-DAR$\left(1\right)$ model may even handle negative forms of serial dependence. Model fitting was illustrated by two real-data examples: a nominal DNA sequence, and an ordinal time series of cloudiness states. Besides such an immediate application of the novel models, it was also pointed out that types of RS-DAR$\left(1\right)$ model might serve as the base for constructing parsimonious advanced models, such as MTD models or HMMs. This direction deserves further attention by future research. Furthermore, the special case of regime-switching ordinal processes, which has various applications in practice, should be further elaborated. Due to the close connection to bounded-counts time series [14], one may try to adapt regime-switching techniques known from count processes, e.g., hysteresis zones for enabling smoother regime transitions.