Piecewise autoregression for general integer-valued time series

Abstract

This paper proposes a piecewise autoregression for general integer-valued time series. The conditional mean of the process depends on a parameter which is piecewise constant over time. We derive an inference procedure based on a penalized contrast that is constructed from the Poisson quasi-maximum likelihood of the model. The consistency of the proposed estimator is established. From practical applications, we derive a data-driven procedure based on the slope heuristic to calibrate the penalty term of the contrast; and the implementation is carried out through the dynamic programming algorithm, which leads to a procedure of $\mathcal{O}\left(n^2\right)$ time complexity. Some simulation results are provided, as well as the applications to the US recession data and the number of trades in the stock of Technofirst.

Introduction

We consider a ${\mathbb{N}}_0$-valued (${\mathbb{N}}_0=\mathbb{N}\cup \left\{0\right\}$) process $Y=\{Y_t,t\in \mathbb{Z}\}$ where the conditional mean

$${\lambda }_t={\lambda }_t\left({\theta }_t^{\ast }\right)=\mathbb{E}\left(Y_t\right|{\mathcal{F}}_{t-1})$$

is a function (see below) of the whole information ${\mathcal{F}}_{t-1}$ up to time $t-1$ and of an unknown parameter ${\theta }_t^{\ast }$ belongs to a compact subset $\Theta \subset {\mathbb{R}}^d$ ($d\in \mathbb{N}$). The inference in the cases where ${\theta }_t^{\ast }={\theta }^{\ast }$ is constant or the distribution of $Y_t|{\mathcal{F}}_{t-1}$ is known has been studied by many authors in several directions; see for instance, Fokianos et al. (2009), Fokianos and Tjøstheim, 2011, Fokianos and Tjøstheim, 2012, Davis and Liu (2016), Douc et al. (2017) among others, for some recent works. We consider here a more general setting where ${\theta }_t^{\ast }$ is piecewise constant (multiple change-point problem) and that the distribution of $Y_t|{\mathcal{F}}_{t-1}$ is unknown. We refer to Franke et al. (2012), Kang and Lee (2014), Doukhan and Kengne (2015), Leung et al. (2017) and the references therein for some tests for change-point detection in integer-valued time series.

Let $(Y_1,\dots ,Y_n)$ be a trajectory generated as in model (1.1) and assume that the parameter ${\theta }_t^{\ast }$ is piecewise constant. Also, assume that $\exists K^{\ast }\in \mathbb{N}$, ${\underline{\theta }}^{\ast }=({\theta }_1^{\ast },\dots ,{\theta }_{K^{\ast }}^{\ast })\in {\Theta }^{K^{\ast }}$ and $0<t_1^{\ast }<\cdots <t_{K^{\ast }-1}^{\ast }<n$ such that, $\{Y_t,t_{j-1}^{\ast }<t\le t_j^{\ast }\}$ is generated from the $j$th stationary regime ; i.e., it is a trajectory of the process $\{Y_{t,j},t\in \mathbb{Z}\}$ (which are not actually observed for $j=1,\dots ,K^{\ast }$, see Section 2 for some details) satisfying:

$$\mathbb{E}\left(Y_{t,j}\right|{\mathcal{F}}_{t-1})=f(Y_{t-1,j},Y_{t-2,j},\dots ;{\theta }_j^{\ast }),\forall t_{j-1}^{\ast }<t\le t_j^{\ast }$$

where ${\mathcal{F}}_t=\sigma (Y_{s,j},s\le t,j=1,\dots ,K^{\ast }-1)$ is the $\sigma$-field generated by the whole information up to time $t$ and $f$ is a measurable non-negative function assumed to be known up to the parameter ${\theta }_t^{\ast }$. $K^{\ast }$ is the number of segments (or regimes) of the model; the $j$th segment corresponds to $\{t_{j-1}^{\ast }+1,t_{j-1}^{\ast }+2,\dots ,t_j^{\ast }\}$ and depends on the parameter ${\theta }_j^{\ast }$. $t_1^{\ast },\dots ,t_{K^{\ast }-1}^{\ast }$ are the change-point locations; by convention, $t_0^{\ast }=-\infty$ and $t_{K^{\ast }}^{\ast }=\infty$. To ensure the identifiability of the change-point locations, it is reasonable to assume that ${\theta }_j^{\ast }\ne {\theta }_{j+1}^{\ast }$ for $j=1,\dots ,K^{\ast }-1$. The case $K^{\ast }=1$ corresponds to the model without change. In the sequel, we assume that the random variables $Y_t$, $t\in \mathbb{Z}$ have the same (up to the parameter ${\theta }_t^{\ast }$) distribution $P$ and denote by $P\left(\cdot \right|{\mathcal{F}}_{t-1})$ the distribution of $Y_t|{\mathcal{F}}_{t-1}$. For instance, for an INGARCH$(p^{\ast },q^{\ast })$ representation, we have

$${\lambda }_t={\alpha }_{0,j}^{\ast }+\sum _{i=1}^{q^{\ast }}{\alpha }_{i,j}^{\ast }{Y}_{t-i}+\sum _{i=1}^{p^{\ast }}{\beta }_{i,j}^{\ast }{\lambda }_{t-i},\text{ for all }{t}_{j-1}^{\ast }<t\le t_j^{\ast },$$

where ${\alpha }_{0,j}^{\ast }>0$, ${\alpha }_{1,j}^{\ast },\dots ,{\alpha }_{q^{\ast },j}^{\ast },{\beta }_{1,j}^{\ast },\dots ,{\beta }_{p^{\ast },j}^{\ast }\ge 0$. The parameters vector of the $j$th regime is ${\theta }_j^{\ast }=({\alpha }_{0,j}^{\ast },{\alpha }_{1,j}^{\ast },\dots ,{\alpha }_{q^{\ast },j}^{\ast },{\beta }_{1,j}^{\ast },\dots ,{\beta }_{p^{\ast },j}^{\ast })$. Therefore, $\Theta$ is a compact subset of $(0,\infty )\times {[0,\infty )}^{p^{\ast }+q^{\ast }}$ such that for all $\theta =({\alpha }_0,{\alpha }_1,\dots ,{\alpha }_{q^{\ast }},{\beta }_1,\dots ,{\beta }_{p^{\ast }})\in \Theta$, ${\sum }_{i=1}^{q^{\ast }}{\alpha }_i+{\sum }_{i=1}^{p^{\ast }}{\beta }_i<1$. For all $j=1,\dots ,K^{\ast }$, we assume that ${\theta }_j^{\ast }\in \Theta$; hence, there exists a sequence of non-negative real numbers ${\left({\psi }_k\left({\theta }_j^{\ast }\right)\right)}_{k\ge 0}$ such that ${\lambda }_t={\psi }_0\left({\theta }_j^{\ast }\right)+{\sum }_{k\ge 1}{\psi }_k\left({\theta }_j^{\ast }\right)Y_{t-k}$. Then, $f(y_1,y_2,\dots ;{\theta }_j^{\ast })={\psi }_0\left({\theta }_j^{\ast }\right)+{\sum }_{k\ge 1}{\psi }_k\left({\theta }_j^{\ast }\right)y_k$ for any $(y_1,y_2,\dots )\in {\mathbb{N}}_0^{\infty }$. For instance, if the distribution $P$ is Poisson, negative binomial or binary, then we get respectively a Poisson, negative binomial, binary INGARCH process; see some examples in Section 4.

Our main focus of interest is the estimation of the unknown parameters $\left(K^{\ast },{\left(t_j^{\ast }\right)}_{1\le j\le K^{\ast }-1},{\left({\theta }_j^{\ast }\right)}_{1\le j\le K^{\ast }}\right)$ in the model (1.2). This can be viewed as a classical model selection problem. Assume that the observations $Y_1,\dots ,Y_n$ are generated from (1.2). Let $K_{max}$ be the upper bound of the number of segments (note that $K_{max}<n$). Denote by ${\mathcal{M}}_n$ the set of partitions of $〚1,n〛$ into at most $K_{max}$ contiguous segments. Set $m=\{T_1,\dots ,T_K\}$ a generic element of $K$ segments in ${\mathcal{M}}_n$. Consider the collection $\{{\mathcal{S}}_m,m\in {\mathcal{M}}_n\}$ where, for a given $m\in {\mathcal{M}}_n$, ${\mathcal{S}}_m$ is the families of sequence $\left({\theta }_t\right)$ which are piecewise constant on the partition $m$. Any $\vartheta =\left({\theta }_t\right)\in {\mathcal{S}}_m$ depends on the parameter $\underline{\theta }=({\theta }_1,\dots ,{\theta }_K)$ which is the piecewise values of ${\theta }_t$ on each segment. Set $\mathcal{S}={\cup }_{m\in {\mathcal{M}}_n}{\mathcal{S}}_m$. Denote by $\vartheta$ a generic element of $\mathcal{S}$, with partition $m$ and parameter $\underline{\theta }$. $\left|\underline{\theta }\right|=K$ denotes the number of the piecewise segments, also called the dimension of $\vartheta$. The true model ${\vartheta }^{\ast }$ with dimension $K^{\ast }$, depends on a partition $m^{\ast }$ and the parameter ${\underline{\theta }}^{\ast }$.

For any $\vartheta \in \mathcal{S}$, set ${\lambda }_t^{\vartheta }={\sum }_{k=1}^K{\lambda }_t\left({\theta }_k\right)1_{t\in T_k}$ and denote by $P\left(\cdot \right|{\mathcal{F}}_{t-1},\vartheta )$ the distribution of $Y_t|{\mathcal{F}}_{t-1},\vartheta$; let $p\left(\cdot \right|{\mathcal{F}}_{t-1},\vartheta )=p(\cdot ;{\lambda }_t^{\vartheta })$ be the probability density function of this distribution. For $\vartheta \in \mathcal{S}$, let $P_{n,\vartheta }$ be the conditional distribution of $(Y_1,\dots ,Y_n)|{\mathcal{F}}_{n-1},\vartheta$. We consider the log-likelihood contrast conditioned to $Y_0,Y_{-1},\dots$: $\forall \vartheta \in \mathcal{S}$,

$${\gamma }_n\left(\vartheta \right)≔{\gamma }_n\left(P_{n,\vartheta }\right)=-log{P}_{n,\vartheta }(Y_1,\dots ,Y_n)=-\sum _{t=1}^{n}logp\left(Y_t\right|{\mathcal{F}}_{t-1},\vartheta )=-\sum _{t=1}^{n}logp(Y_t;{\lambda }_t^{\vartheta }).$$

Thus, the minimal contrast estimator ${\widehat{\vartheta }}_m$ of ${\vartheta }^{\ast }$ on the collection ${\mathcal{S}}_m$ is obtained by minimizing the contrast ${\gamma }_n\left(\vartheta \right)$ over $\vartheta \in {\mathcal{S}}_m$; that is, ${\widehat{\vartheta }}_m=\underset{\vartheta \in {\mathcal{S}}_m}{\text{argmin}}{\gamma }_n\left(\vartheta \right)$. The main approaches of the model selection procedures take into account the model complexity and select the estimator ${\widehat{\vartheta }}_{m_n}$ such that, $m_n$ minimizes the penalized criterion

$${\text{crit}}_n\left(m\right)={\gamma }_n\left({\widehat{\vartheta }}_m\right)+{\text{pen}}_n\left(m\right),\text{for all}m\in {\mathcal{M}}_n$$

where ${\text{pen}}_n:{\mathcal{M}}_n\to {\mathbb{R}}_{+}$ is a penalty function, possibly data-dependent. We now address the following issues.

(i) Semi-parametric setting. Kashikar et al. (2013) have carried out structural breaks in Poisson INAR process from the MCMC and Gibbs sampling approach. Cleynen and Lebarbier, 2014, Cleynen and Lebarbier, 2017 have recently considered the change-point type problem (1.2) with i.i.d. observations; in their works, the distribution $P$ is assumed to be known and could be Poisson, Negative binomial or belongs to the exponential family distribution. From the practical viewpoint, we consider the case where $P$ is unknown and deal with the Poisson quasi-likelihood (see for instance, Ahmad and Francq, 2016). So in the sequel, ${\gamma }_n$ is the Poisson quasi-likelihood contrast and ${\widehat{\vartheta }}_m$ is the Poisson quasi-maximum likelihood estimator (PQMLE).

(ii) Multiple change-point problem from a non-asymptotic point of view. This question is tacked by model selection approach. Numerous works have been devoted to this issue; see among others, Lebarbier (2005), Arlot and Massart (2009), Cleynen and Lebarbier, 2014, Cleynen and Lebarbier, 2017 and Arlot et al. (2016).

In this (quasi)log-likelihood framework, it is more usual to consider the Kullback–Leibler risk. For any $\vartheta \in \mathcal{S}$, the Kullback–Leibler divergence between $P_{n,{\vartheta }^{\ast }}$ and $P_{n,\vartheta }$ is

$$KL({\vartheta }^{\ast },\vartheta )≔KL(P_{n,{\vartheta }^{\ast }},P_{n,\vartheta })=\mathbb{E}\left[log\frac{P_{n,{\vartheta }^{\ast }}(Y_1,\dots ,Y_n)}{P_{n,\vartheta }(Y_1,\dots ,Y_n)}\right]=\sum _{t=1}^n\mathbb{E}\left[log\frac{p\left(Y_t\right|{\mathcal{F}}_{t-1},{\vartheta }^{\ast })}{p\left(Y_t\right|{\mathcal{F}}_{t-1},\vartheta )}\right]=\sum _{t=1}^n\mathbb{E}\left[logp(Y_t;{\lambda }_t^{{\vartheta }^{\ast }})\right]-\sum _{t=1}^n\mathbb{E}\left[logp(Y_t;{\lambda }_t^{\vartheta })\right],$$

where $\mathbb{E}$ denotes the expectation with respect to the true distribution of the observations. In the case where ${\gamma }_n$ is the likelihood contrast, we get $KL({\vartheta }^{\ast },\vartheta )=\mathbb{E}[{\gamma }_n\left(\vartheta \right)-{\gamma }_n\left({\vartheta }^{\ast }\right)]$. The “ideal” partition $m\left({\vartheta }^{\ast }\right)$ (the one whose estimator is closest to ${\vartheta }^{\ast }$ according to the Kullback–Leibler risk) satisfying:

$$m\left({\vartheta }^{\ast }\right)=\underset{m\in {\mathcal{M}}_n}{\text{argmin}}\mathbb{E}\left[KL({\vartheta }^{\ast },{\widehat{\vartheta }}_m)\right].$$

The corresponding estimator ${\widehat{\vartheta }}_{m\left({\vartheta }^{\ast }\right)}$, called the oracle, depends on the true sample distribution, and cannot be computed in practice. The goal is to calibrate the penalty term, such that the segmentation $\widehat{m}$ provides an estimator ${\widehat{\vartheta }}_{\widehat{m}}$ where the risk of ${\widehat{\vartheta }}_{\widehat{m}}$ is close as possible to the risk of the oracle, namely such that

$$\mathbb{E}\left[KL({\vartheta }^{\ast },{\widehat{\vartheta }}_{\widehat{m}})\right]\le C\mathbb{E}\left[KL({\vartheta }^{\ast },{\widehat{\vartheta }}_{m\left({\vartheta }^{\ast }\right)})\right]$$

for a non-negative constant $C$, expected close to 1. This issue is addressed in the above mentioned papers, and the results obtained are heavily relied on the independence of the observations. In our setting here, it seems to be a more difficult task. But, we believe that the coupling method can be used as in Lerasle (2011) to overcome this difficulty. We leave this question as the topic of a different research project.

(iii) Multiple change-point problem from an asymptotic point of view. The aim here is to consistently estimate the parameters of the change-point model. This issue has been addressed by several authors using the classical contrast/criteria optimization or binary/sequential segmentation/estimation; see for instance, Bai and Perron (1998), Davis et al. (2008), Harchaoui and Lévy-Leduc (2010), Bardet et al. (2012), Davis and Yau (2013), Davis et al. (2016), Ma and Yau (2016), Yau and Zhao (2016), Inclán and Tiao (1994), Bai (1997), Fryzlewicz and Subba Rao (2014), Fryzlewicz (2014), among others, for some advanced towards this issue. These works and many other papers in the literature on the asymptotic study of multiple change-point problem are often focused on continuous valued time series; moreover, the case of a large class of semi-parametric model for discrete-valued time series (such as those discussed earlier) have not yet addressed.

We consider (1.2) and derive a penalized contrast of type (1.3). We assume that there exists a partition ${\underline{\tau }}^{\ast }$ of $[0,1]$ such that $\left[{\underline{\tau }}^{\ast }n\right]=m^{\ast }$, where $\left[{\underline{\tau }}^{\ast }n\right]$ is the corresponding partition of $〚1,n〛$ obtained from ${\underline{\tau }}^{\ast }$. We provide sufficient conditions on the penalty ${\text{pen}}_n$, for which the estimators $\widehat{m}$ and ${\widehat{\vartheta }}_{\widehat{m}}$ are consistent; that is:

$$\left(\left|\widehat{m}\right|,\frac{\widehat{m}}{n},{\widehat{\vartheta }}_{\widehat{m}}\right)\begin{array}{c}\hfill \stackrel{\mathcal{P}}{⟶}\hfill \\ \hfill n\to \infty \hfill \end{array}\left(K^{\ast },{\underline{\tau }}^{\ast },{\vartheta }^{\ast }\right)$$

where $\frac{\widehat{m}}{n}$ is the corresponding partition of $[0,1]$ obtained from $\widehat{m}$.

The paper is organized as follows. In Section 2, we set some notations, assumptions and define the Poisson QMLE. In Section 3, we derive the estimation procedure and provide the main results. Some simulation results are displayed in Section 4 whereas Section 5 focuses on applications on the US recession data and the daily number of trades in the stock of Technofirst. Section 6 is devoted to a summary and conclusion. The Supporting Information provides the proofs of the main results.