Variable selection in sparse high-dimensional GLARMA models

TL;DR

This paper introduces a new variable selection method for high-dimensional GLARMA models that combines ARMA coefficient estimation with regularized GLM techniques, demonstrating improved accuracy and efficiency.

Contribution

It presents a novel approach integrating ARMA coefficient estimation with regularized GLM methods for sparse high-dimensional models, with proven consistency and low computational load.

Findings

Enhanced variable selection accuracy compared to existing methods

Method shows low computational complexity and high efficiency

Numerical experiments confirm improved performance on synthetic data

Abstract

In this paper, we propose a novel variable selection approach in the framework of sparse high-dimensional GLARMA models. It consists in combining the estimation of the autoregressive moving average (ARMA) coefficients of these models with regularized methods designed for Generalized Linear Models (GLM). The properties of our approach are investigated both from a theoretical and a numerical point of view. More precisely, we establish in a specific case the consistency of the ARMA part coefficient estimators. We explain how to implement our approach and we show that it is very attractive since it benefits from a low computational load. We also assess the performance of our methodology using synthetic data and compare it with alternative approaches. Our numerical experiments show that combining the estimation of the ARMA part coefficients with regularized methods designed for GLM dramatically improves the variable selection performance.

Equations

$$Y_t|{\mathcal{F}}_{t-1}\sim \mathcal{P}\left({\mu }_t^{\star }\right),$$

$${\mu }_t^{\star }=\exp (W_t^{\star })\text{ with }{W}_t^{\star }={\beta }_0^{\star }+\sum _{i=1}^p{\beta }_i^{\star }{x}_{t,i}+Z_t^{\star },$$

$$Z_t^{\star }=\sum _{j=1}^q{\gamma }_j^{\star }{E}_{t-j}^{\star }\text{ with }{E}_t^{\star }=\frac{Y_t-{\mu }_t^{\star }}{{\mu }_t^{\star }}=Y_t\exp (-W_t^{\star })-1,$$

$$L(\mathbf{\delta })=\sum _{t=1}^n\left(Y_t{W}_t(\mathbf{\delta })-\exp (W_t(\mathbf{\delta }))\right),$$

$$W_t(\mathbf{\delta })={\mathbf{\beta }}^{\prime }{x}_t+Z_t(\mathbf{\delta })={\beta }_0+\sum _{i=1}^p{\beta }_i{x}_{t,i}+\sum _{j=1}^q{\gamma }_j{E}_{t-j}(\mathbf{\delta }),$$

$$E_t(\mathbf{\delta })=Y_t\exp (-W_t(\mathbf{\delta }))-1,\text{mbox}ift>0\text{mbox}and{E}_t(\mathbf{\delta })=0\text{mbox},ift\le 0.$$

$$\hat{\mathbf{\delta }}={\text{Argmax}}_{\mathbf{\delta }}L(\mathbf{\delta }),$$

$$\frac{\partial L}{\partial \mathbf{\delta }}(\mathbf{\delta })=\sum _{t=1}^n(Y_t-\exp (W_t(\mathbf{\delta }))\frac{\partial W_t}{\partial \mathbf{\delta }}(\mathbf{\delta }),$$

$$\frac{\partial W_t}{\partial \mathbf{\delta }}(\mathbf{\delta })=\frac{\partial {\mathbf{\beta }}^{\prime }{x}_t}{\partial \mathbf{\delta }}+\frac{\partial Z_t}{\partial \mathbf{\delta }}(\mathbf{\delta }),$$

$${\mathbf{\delta }}^{(r)}={\mathbf{\delta }}^{(r-1)}-\frac{{\partial }^{2}L}{\partial {\mathbf{\delta }}^{\prime }\partial \mathbf{\delta }}({\mathbf{\delta }}^{(r-1)}{)}^{-1}\frac{\partial L}{\partial \mathbf{\delta }}({\mathbf{\delta }}^{(r-1)}),$$

$$\frac{{\partial }^{2}L}{\partial {\mathbf{\delta }}^{\prime }\partial \mathbf{\delta }}(\mathbf{\delta })=\sum _{t=1}^n(Y_t-\exp (W_t(\mathbf{\delta }))\frac{{\partial }^2{W}_t}{\partial {\mathbf{\delta }}^{\prime }\partial \mathbf{\delta }}(\mathbf{\delta })-\sum _{t=1}^n\exp (W_t(\mathbf{\delta }))\frac{\partial W_t}{\partial {\mathbf{\delta }}^{\prime }}(\mathbf{\delta })\frac{\partial W_t}{\partial \mathbf{\delta }}(\mathbf{\delta }).$$

$$\tilde{L}(\mathbf{\beta }):=L({\beta }_0,\dots ,{\beta }_p,\hat{\gamma })=\tilde{L}(\tilde{\mathbf{\beta }})+\frac{\partial L}{\partial \mathbf{\beta }}(\tilde{\mathbf{\beta }},\hat{\mathbf{\gamma }})(\mathbf{\beta }-\tilde{\mathbf{\beta }})+\frac{1}{2}(\mathbf{\beta }-\tilde{\mathbf{\beta }}{)}^{\prime }\frac{{\partial }^{2}L}{\partial \mathbf{\beta }\partial {\mathbf{\beta }}^{\prime }}(\tilde{\mathbf{\beta }},\hat{\mathbf{\gamma }})(\mathbf{\beta }-\tilde{\mathbf{\beta }}),$$

$$\frac{\partial L}{\partial \mathbf{\beta }}=\left(\frac{\partial L}{\partial {\beta }_0},\dots ,\frac{\partial L}{\partial {\beta }_p}\right)\text{ and }\frac{{\partial }^{2}L}{\partial \mathbf{\beta }\partial {\mathbf{\beta }}^{\prime }}={\left(\frac{{\partial }^{2}L}{\partial {\beta }_j\partial {\beta }_k}\right)}_{0\le j,k\le p}.$$

$$\tilde{L}(\mathbf{\beta })=\tilde{L}(\tilde{\mathbf{\beta }})+\frac{\partial L}{\partial \mathbf{\beta }}(\tilde{\mathbf{\beta }},\hat{\mathbf{\gamma }})U(\mathbf{\nu }-\tilde{\mathbf{\nu }})-\frac{1}{2}(\mathbf{\nu }-\tilde{\mathbf{\nu }}{)}^{\prime }\Lambda (\mathbf{\nu }-\tilde{\mathbf{\nu }}),$$

$$\hat{\mathbf{\beta }}(\lambda )={\text{Argmin}}_{\mathbf{\beta }}\left\{-{\tilde{L}}_Q(\mathbf{\beta })+\lambda \Vert \mathbf{\beta }{\Vert }_1\right\},$$

$$-{\tilde{L}}_Q(\mathbf{\beta })=\frac{1}{2}\Vert \mathcal{Y}-\mathcal{X}\mathbf{\beta }{\Vert }_2^2,$$

$$\mathcal{Y}={\Lambda }^{1/2}{U}^{\prime }\tilde{\mathbf{\beta }}-{\Lambda }^{-1/2}{U}^{\prime }{\left(\frac{\partial L}{\partial \mathbf{\beta }}(\tilde{\mathbf{\beta }},\hat{\mathbf{\gamma }})\right)}^{\prime },\mathcal{X}={\Lambda }^{1/2}{U}^{\prime }$$

$${\mu }_t^{\star }=\exp (W_t^{\star })\text{ with }{W}_t^{\star }={\beta }_0^{\star }+Z_t^{\star }.$$

$${\hat{\gamma }}_1={\text{Argmax}}_{{\gamma }_1\in \Gamma }L({\beta }_0^{\star },{\gamma }_1),$$

$$L({\beta }_0^{\star },{\gamma }_1)=\sum _{t=1}^n\left(Y_t{W}_t({\beta }_0^{\star },{\gamma }_1)-\exp (W_t({\beta }_0^{\star },{\gamma }_1)\right),$$

$$W_t({\beta }_0^{\star },{\gamma }_1)={\beta }_0^{\star }+Z_t({\gamma }_1)={\beta }_0^{\star }+{\gamma }_1{E}_{t-1}({\gamma }_1),$$

$$E_{t-1}({\gamma }_1)=Y_{t-1}\exp (-W_{t-1}({\beta }_0^{\star },{\gamma }_1))-1,\text{ if }t>1\text{ and }{E}_{t-1}({\gamma }_1)=0,\text{ if }t\le 1.$$

$$\frac{1}{n}L({\beta }_0^{\star },{\gamma }_1)\stackrel{p}{⟶}\mathcal{L}({\gamma }_1):=\mathbb{E}\left[Y_3{W}_3({\beta }_0^{\star },{\gamma }_1)-\exp (W_3({\beta }_0^{\star },{\gamma }_1)\right],\text{ as n tends to infinity.}$$

$$\underset{{\gamma }_1\in \Gamma }{\sup }\left|\frac{L({\beta }_0^{\star },{\gamma }_1)}{n}-\mathcal{L}({\gamma }_1)\right|\stackrel{p}{⟶}0,\text{ as n tends to infinity,}$$

$$\frac{\partial W_t}{\partial \mathbf{\delta }}(\mathbf{\delta })=\frac{\partial {\mathbf{\beta }}^{\prime }{x}_t}{\partial \mathbf{\delta }}+\frac{\partial Z_t}{\partial \mathbf{\delta }}(\mathbf{\delta }),$$

$$\frac{\partial W_t}{\partial {\beta }_k}$$

$$=x_{t,k}-\sum _{j=1}^{q\wedge (t-1)}{\gamma }_j{Y}_{t-j}\frac{\partial W_{t-j}}{\partial {\beta }_k}\exp (-W_{t-j})=x_{t,k}-\sum _{j=1}^{q\wedge (t-1)}{\gamma }_j(1+E_{t-j})\frac{\partial W_{t-j}}{\partial {\beta }_k},$$

$$\frac{\partial W_t}{\partial {\gamma }_{\ell }}$$

$$=E_{t-\ell }-\sum _{j=1}^{q\wedge (t-1)}{\gamma }_j{Y}_{t-j}\frac{\partial W_{t-j}}{\partial {\gamma }_{\ell }}\exp (-W_{t-j})=E_{t-\ell }-\sum _{j=1}^{q\wedge (t-1)}{\gamma }_j(1+E_{t-j})\frac{\partial W_{t-j}}{\partial {\gamma }_{\ell }},$$

$$\frac{\partial W_1}{\partial {\beta }_k}$$

$$\frac{\partial W_2}{\partial {\beta }_k}$$

$$W_1={\mathbf{\beta }}^{\prime }{x}_1\text{ and }{E}_1=Y_1\exp (-W_1)-1.$$

Full text

1. Introduction

Discrete-valued time series arise in a wide variety of fields ranging from finance to molecular biology and public health. For instance, we can mention the number of transactions in stocks in the finance field, see Brännäs and Quoreshi (2010). In the field of molecular biology, modeling RNA-Seq kinetics data is a challenging issue, see Thorne (2018) and in the public health context, there is an interest in the modeling of daily asthma presentations in a given hospital, see Souza et al. (2014).

The literature on modeling discrete-valued time series is becoming increasingly abundant, see Davis et al. (2016) for a review. Different classes of models have been proposed such as the Integer Autoregressive Moving Average (INARMA) models and the generalized state space models.

The Integer Autoregressive process of order 1 (INAR(1)) was first introduced by McKenzie (1985) and the Integer-valued Moving Average (INMA) process is described in Al-Osh and Alzaid (1988). One of the attractive features of INARMA processes is that their autocorrelation structure is similar to the one of autoregressive moving average (ARMA) models. However, it has to be noticed that statistical inference in these models is generally complicated and requires to develop intensive computational approaches such as the efficient MCMC algorithm devised by Neal and Subba Rao (2007) for INARMA processes of known AR and MA orders. This strategy was extended to unknown AR and MA orders by Enciso-Mora et al. (2009). For further references on INARMA models, we refer the reader to Weiss (2018).

The other important class of models for discrete-valued time series is the one of generalized state space models which can have a parameter-driven and an observation-driven version, see Davis et al. (1999) for a review.

The main difference between these two versions is that in parameter-driven models, the state vector evolves independently of the past history of the observations whereas the state vector depends on the past observations in observation-driven models. More precisely, in parameter-driven models, let $(\nu_{t})$ be a stationary process, the observations $Y_{t}$ are thus modeled as follows: conditionally on $(\nu_{t})$, $Y_{t}$ has a Poisson distribution of parameter $\exp(\beta_{0}^{\star}+\sum_{i=1}^{p}\beta_{i}^{\star}x_{t,i}+\nu_{t})$, where the $x_{t,i}$’s are the $p$ regressor variables (or covariates). Estimating the parameters in such models has a very high computational load, see Jung and Liesenfeld (2001).

Observation-driven models initially proposed by Cox et al. (1981) and further studied in Zeger and Qaqish (1988) do not have this computational drawback and are thus considered as a promising alternative to parameter-driven models. Different kinds of observation-driven models can be found in the literature: the Generalized Linear Autoregressive Moving Average (GLARMA) models introduced by Davis et al. (1999) and further studied in Davis et al. (2003, 2005); Dunsmuir (2015) and the (log-)linear Poisson autoregressive models introduced in Fokianos et al. (2009); Fokianos and Tjøstheim (2011, 2012). Note that GLARMA models cannot be seen as a particular case of the log-linear Poisson autoregressive models.

In the following, we shall consider the GLARMA model introduced in Davis et al. (2005) with additional covariates. More precisely, given the past history $\mathcal{F}_{t-1}=\sigma(Y_{s},s\leq t-1)$, we assume that

[TABLE]

where $\mathcal{P}(\mu)$ denotes the Poisson distribution with mean $\mu$. In (1),

[TABLE]

where the $x_{t,i}$’s are the $p$ regressor variables ($p\geq 1$),

[TABLE]

where $1\leq q\leq\infty$ and $E_{t}^{\star}=0$ for all $t\leq 0$. Here, the $E_{t}^{\star}$’s correspond to the working residuals in classical Generalized Linear Models (GLM), which means that we limit ourselves to the case $\lambda=1$ in the more general definition: $E_{t}^{\star}=(Y_{t}-\mu_{t}^{\star}){\mu_{t}^{\star}}^{-\lambda}$. Note that in the case where $q=\infty$, $(Z_{t}^{\star})$ satisfies the ARMA-like recursions given in Equation (4) of Davis et al. (2005). The model defined by (1), (2) and (3) is thus referred as a GLARMA model.

The main goal of this paper is to introduce a novel variable selection approach in the deterministic part (covariates) of high-dimensional sparse GLARMA models that is in (1) and (2) where $p$ is large and when the vector of the $\beta_{i}$’s is sparse. The novel approach that we propose consists in combining a procedure for estimating the ARMA part coefficients with regularized methods designed for GLM.

The paper is organized as follows. We firstly describe an estimation procedure for the ARMA part of the GLARMA model defined in (1), (2) and (3), see Section 2.1, and establish its consistency in a specific case, see Section 2.4. Secondly, we propose a novel variable selection approach in the regression part of the sparse high-dimensional model (1) and explain how to combine it with the estimation procedure of the ARMA part coefficients, see Section 2.2. The practical implementation of our approach is given in Section 2.3. Thirdly, in Section 3, some numerical experiments are provided to illustrate our method and to compare its performance to alternative approaches on finite sample size data. The proofs of the theoretical results are given in Section 4.

2. Statistical inference

2.1. Estimation procedure

For estimating the parameter $\boldsymbol{\delta}^{\star}=(\boldsymbol{\beta}^{\star\prime},\boldsymbol{\gamma}^{\star\prime})$ where $\boldsymbol{\beta}^{\star}=(\beta_{0}^{\star},\beta_{1}^{\star},\dots,\beta_{p}^{\star})^{\prime}$ is the vector of regressor coefficients defined in (2) and $\boldsymbol{\gamma}^{\star}=(\gamma_{1}^{\star},\dots,\gamma_{q}^{\star})^{\prime}$ is the vector of the ARMA part coefficients defined in (3), we maximize the following criterion, based on the conditional log-likelihood, with respect to $\boldsymbol{\delta}=(\boldsymbol{\beta}^{\prime},\boldsymbol{\gamma}^{\prime})$, with $\boldsymbol{\beta}=(\beta_{0},\beta_{1},\dots,\beta_{p})^{\prime}$ and $\boldsymbol{\gamma}=(\gamma_{1},\dots,\gamma_{q})^{\prime}$:

[TABLE]

where

[TABLE]

with $x_{t}=(x_{t,0},x_{t,1},\dots,x_{t,p})^{\prime}$, $x_{t,0}=1$ for all $t$ and

[TABLE]

For further details on the choice of this criterion, we refer the reader to Davis et al. (2005). To obtain $\widehat{\boldsymbol{\delta}}$ defined by

[TABLE]

we consider the first derivatives of $L$:

[TABLE]

where

[TABLE]

$\boldsymbol{\beta}$, $x_{t}$ and $Z_{t}$ being defined in (5). The computations of the first derivatives of $W_{t}$ are detailed in Section 4.1.1.

Since the first derivatives of $W_{t}$ are recursively defined, it is not possible to obtain a closed-form formula for $\widehat{\boldsymbol{\delta}}$. Thus, in order to compute $\widehat{\boldsymbol{\delta}}$, we shall use the following Newton-Raphson algorithm. More precisely, we start from an initial value for $\boldsymbol{\delta}$ denoted $\boldsymbol{\delta}^{(0)}$. Then, we use the following recursion for $r\geq 1$:

[TABLE]

where $\frac{\partial^{2}L}{\partial\boldsymbol{\delta}^{\prime}\partial\boldsymbol{\delta}}$ corresponds to the Hessian matrix of $L$ and is defined in (9). Hence, it requires the computation of the first and second derivatives of $L$. We already explained how to compute the first derivatives of $L$. As for the second derivatives of $L$, it can be obtained as follows:

[TABLE]

The computations of the second derivatives of $W_{t}$ are detailed in Section 4.1.2.

Further details on the choice of $\boldsymbol{\delta}^{(0)}$ and the number of iterations to use will be given in Section 3.

2.2. Variable selection

To perform variable selection in the $\beta_{i}^{\star}$ of Model (2) that is to obtain a sparse estimator of $\beta_{i}^{\star}$, we shall use a methodology inspired by Friedman et al. (2010) for fitting generalized linear models with $\ell_{1}$ penalties. It consists in penalizing a quadratic approximation to the log-likelihood obtained by a Taylor expansion. Hence, denoting $\widetilde{\boldsymbol{\beta}}=(\widetilde{\beta}_{0},\dots,\widetilde{\beta}_{p})^{\prime}$ the current estimate of the parameter $\boldsymbol{\beta}^{\star}=({\beta_{0}}^{\star},\dots,{\beta_{p}}^{\star})^{\prime}$, we obtain the following quadratic approximation where $\widehat{\boldsymbol{\gamma}}=(\widehat{\gamma}_{1},\dots,\widehat{\gamma}_{q})^{\prime}$ is the estimate of $\boldsymbol{\gamma}^{\star}=({\gamma_{1}}^{\star},\dots,{\gamma_{q}}^{\star})^{\prime}$ obtained in Section 2.1:

[TABLE]

where

[TABLE]

Thus,

[TABLE]

where $U\Lambda U^{\prime}$ is the singular value decomposition of the positive semidefinite symmetric matrix $-\frac{\partial^{2}L}{\partial\boldsymbol{\beta}\partial\boldsymbol{\beta}^{\prime}}(\widetilde{\boldsymbol{\beta}},\widehat{\boldsymbol{\gamma}})$ and $\boldsymbol{\nu}-\widetilde{\boldsymbol{\nu}}=U^{\prime}(\boldsymbol{\beta}-\widetilde{\boldsymbol{\beta}})$.

In order to obtain a sparse estimator of $\boldsymbol{\beta}^{\star}$, we propose using $\widehat{\boldsymbol{\beta}}(\lambda)$ defined by

[TABLE]

for a positive $\lambda$, where $\|\boldsymbol{\beta}\|_{1}=\sum_{k=0}^{p}|\beta_{k}|$ and $\widetilde{L}_{Q}(\boldsymbol{\beta})$ denotes the quadratic approximation of the log-likelihood. This quadratic approximation is defined by

[TABLE]

with

[TABLE]

and $\|\cdot\|_{2}$ denoting the $\ell_{2}$ norm in $\mathbb{R}^{p+1}$. Computational details for obtaining the expression (12) of $\widetilde{L}_{Q}(\boldsymbol{\beta})$ appearing in Criterion (11) are provided in Section 4.2.

The parameter $\lambda$ is tuned thanks to the stability selection approach devised by Meinshausen and Bühlmann (2010). For further details, we refer the reader to Section 3.

2.3. Practical implementation

We summarize hereafter the different steps of our methodology.

• •

First step: Initialization. We take for $\boldsymbol{\beta}^{(0)}$ the estimator of $\boldsymbol{\beta}^{\star}$ obtained by fitting a GLM to the observations $Y_{1},\dots,Y_{n}$ thus ignoring the ARMA part of the model. For $\boldsymbol{\gamma}^{(0)}$, we take the null vector.

• •

Second step: Newton-Raphson algorithm. We use the recursion defined in (8) with the initialization $\boldsymbol{\delta}^{(0)}=(\boldsymbol{\beta}^{(0)},\boldsymbol{\gamma}^{(0)})$ obtained in the first step and we stop at the iteration $R$ such that $\|\boldsymbol{\delta}^{(R)}-\boldsymbol{\delta}^{(R-1)}\|_{\infty}<10^{-6}$.

• •

Third step: Variable selection. To obtain a sparse estimator of $\boldsymbol{\beta}^{\star}$, we use the criterion (11) where $\widetilde{\beta}$ and $\widehat{\boldsymbol{\gamma}}$ appearing in (22) are replaced by $\boldsymbol{\beta}^{(0)}$ and $\boldsymbol{\gamma}^{(R)}$ obtained in the first and second steps, respectively. This step provides $\widehat{\boldsymbol{\beta}}(\lambda)$ for different values of $\lambda$.

• •

Fourth step: Choice of $\lambda$. To choose the value of $\lambda$ and thus the final estimator $\widehat{\boldsymbol{\beta}}$ of $\boldsymbol{\beta}^{\star}$, we use the stability selection approach devised by Meinshausen and Bühlmann (2010).

2.4. Consistency results

In this section, we shall establish in the case where $q=1$ the consistency of the parameter $\gamma_{1}^{\star}$ from $Y_{1},\dots,Y_{n}$ defined in (1) and (3) where (2) is replaced by

[TABLE]

Note that some theoretical results have already been obtained in this framework (no covariates and $q=1$) by Davis et al. (2003) and Davis et al. (2005). However, here, we provide, on the one hand, a more detailed version of the proof of these results and on the other hand, a proof of the consistency of $\gamma_{1}^{\star}$ based on a stochastic equicontinuity result. We limit ourselves to this framework since the more general one is a framework where the consistency is much more tricky to handle and is beyond the scope of this paper.

Theorem 1**.**

Assume that $Y_{1},\dots,Y_{n}$ satisfy the model defined by (1), (14) and (3) with $q=1$ and $\gamma_{1}^{\star}\in\Gamma$ where $\Gamma$ is a compact set of $\mathbb{R}$ which does not contain 0. Assume also that $(W_{t}^{\star})$ is started with its stationary invariant distribution. Let $\widehat{\gamma}_{1}$ be defined by:

[TABLE]

where

[TABLE]

with

[TABLE]

[TABLE]

Then $\widehat{\gamma}_{1}\stackrel{{\scriptstyle p}}{{\longrightarrow}}\gamma_{1}^{\star}$, as $n$ tends to infinity, where $\stackrel{{\scriptstyle p}}{{\longrightarrow}}$ denotes the convergence in probability.

The proof of Theorem 1 is based on the following propositions which are proved in Section 4. These propositions are the classical arguments for establishing consistency results of maximum likelihood estimators. Note that we shall explain in the proof of Proposition 1 why a stationary invariant distribution for $(W_{t}^{\star})$ does exist. The main tools used for proving Propositions 1 and 3 are the Markov property and the ergodicity of $(W_{t}^{\star})$.

Proposition 1**.**

For all fixed $\gamma_{1}$, under the assumptions of Theorem 1,

[TABLE]

Proposition 2**.**

The function $\mathcal{L}$ defined in (17) has a unique maximum at the true parameter $\gamma_{1}=\gamma_{1}^{\star}$.

Proposition 3**.**

Under the assumptions of Theorem 1

[TABLE]

where $\mathcal{L}(\gamma_{1})$ is defined in (17).

3. Numerical experiments

3.1. Statistical performance

The goal of this section is to investigate the performance of our method both from a statistical and a numerical point of view.

3.1.1. Estimation of the ARMA part coefficients when $p=0$

In this section, we investigate the statistical performance of our methodology in the case where $Y_{1},\dots,Y_{n}$ satisfy the model defined by (1), (2) and (3) for $n$ in $\{50,100,250,500,1000\}$ in the case where $p=0$, namely when there are no covariates and for $q$ in $\{1,2,3\}$. The results are displayed in Figures 1, 2 and 3. We can see from these figures that the accuracy of the parameters estimations is improved when $n$ increases.

Moreover, it has to be noticed that in this particular context where there are no covariates ($p=0$), the performance of our approach in terms of parameters estimation is similar to the one of the package glarma described in Dunsmuir and Scott (2015).

3.1.2. Sparse estimation of the ${\beta}_{i}^{\star}$

In this section, we assess the performance of our methodology in the case where $Y_{1},\dots,Y_{n}$ satisfy the model defined by (1), (2) and (3) for $n=1000$, $q\in\{1,2,3\}$ and $p=100$. We shall moreover assume that the sparsity in the $\beta_{i}^{\star}$ is very high, namely all the $\beta_{i}^{\star}$ are assumed to be equal to zero except for five of them which are equal to 1.739, 0.387, 0.295, -0.644 and -0.135. The corresponding results are displayed in Figures 4, 5 and 6.

The ROC curves of Figures 4, 5 and 6 display the True Positive Rate (TPR) with respect to the False Positive Rate (FPR). On the one hand, we can see from these figures that the performance of our methodology when $\boldsymbol{\gamma}^{\star}$ is known is on a par with the one of our methodology when $\boldsymbol{\gamma}^{\star}$ is unknown. On the other hand, our methodology outperforms the variable selection approach described in Friedman et al. (2010) which assumes that the observations are the realizations of a Poisson distribution but does not take into account the dependence between the observations.

We can also observe from these figures that the performance of our methodology is not altered by the underestimation of $\boldsymbol{\gamma}^{\star}$ in the different situations: $q=1$, 2 or 3.

3.1.3. Choice of $\lambda$

In order to improve our methodology, we propose hereafter a strategy for tuning the parameter $\lambda$ appearing in (11).

We first take the smallest $\lambda$ provided by the glmnet package for computing (11). This $\lambda$ denoted $\lambda_{\textrm{min}}$ is then used in the stability selection procedure proposed by Meinshausen and Bühlmann (2010) which guarantees the robustness of the selected variables. This latter approach can be described as follows. The vector $\mathcal{Y}$ defined in (22) is randomly split into several subsamples of size $(p+1)/2$, which corresponds to the half of the length of $\mathcal{Y}$. For each subsample, the LASSO criterion is applied with $\lambda=\lambda_{\textrm{min}}$ and the indices $i$ of the non null $\widehat{\beta}_{i}$ are stored. Then, for a given threshold, we keep in the final set of selected variables only the variables appearing a number of times larger than this threshold. In practice, we generated $1000$ subsamples of $\mathcal{Y}$.

Figure 7 displays the results obtained when applying this strategy to observations $Y_{1},\dots,Y_{n}$ satisfying the model defined by (1), (2) and (3) for $n=1000$, $q=1$, $p=100$ and when only five coefficients $\beta_{i}^{\star}$ are not null. We can see from this figure that the positions of the non null coefficients are well retrieved for most of the thresholds and that the number of false positive is higher when the threshold is too low. Based on this figure, taking a threshold equal to 0.9 seems to achieve an interesting trade-off between false and true positives.

This choice is also confirmed by the results of Figure 8 which gives the means of selection frequencies for each position.

3.2. Numerical performance

Figure 9 displays the means and standard errors of the computational times for our variable selection method. We can see from this figure that it takes only around 30 seconds to process observations $Y_{1},\dots,Y_{n}$ satisfying (1), (2) and (3) when $n=1000$, $p=100$, $q=3$ and when the number of replications used in the stability selection step described in Section 3.1.3 is equal to 1000.

4. Proofs

4.1. Computation of the first and second derivatives of $W_{t}$ defined in (5)

The computations given below are similar to those provided in Davis et al. (2005) but are specific to the parametrization $\boldsymbol{\delta}=(\boldsymbol{\beta}^{\prime},\boldsymbol{\gamma}^{\prime})$ considered in this paper.

4.1.1. Computation of the first derivatives of $W_{t}$

By the definition of $W_{t}$ given in (5), we get

[TABLE]

where $\boldsymbol{\beta}$, $x_{t}$ and $Z_{t}$ are defined in (5). More precisely, for all $k\in\{0,\dots,p\}$, $\ell\in\{1,\dots,q\}$ and $t\in\{1,\dots,n\}$, by (6),

[TABLE]

where we used that $E_{t}=0,\;\forall t\leq 0$.

The first derivatives of $W_{t}$ are thus obtained from the following recursive expressions. For all $k\in\{0,\dots,p\}$

[TABLE]

where

[TABLE]

Moreover,

[TABLE]

where

[TABLE]

and so on. In the same way, for all $\ell\in\{1,\dots,q\}$

[TABLE]

and so on, where $E_{t}=0,\;\forall t\leq 0$ and $E_{1}$, $E_{2}$ are defined in (20) and (21), respectively.

4.1.2. Computation of the second derivatives of $W_{t}$

Using (4.1.1) and (19), we get that for all $j,k\in\{0,\dots,p\}$, $\ell,m\in\{1,\dots,q\}$ and $t\in\{1,\dots,n\}$,

[TABLE]

To compute the second derivatives of $W_{t}$, we shall use the following recursive expressions for all $j,k\in\{0,\dots,p\}$

[TABLE]

where $E_{1}$ is defined in (20) and so on. Moreover, for all $k\in\{0,\dots,p\}$ and $\ell\in\{1,\dots,q\}$

[TABLE]

where $E_{t}=0$ for all $t\leq 0$ and the first derivatives of $W_{t}$ are computed in (4.1.1). Note also that

[TABLE]

and so on.

4.2. Computational details for obtaining Criterion (11)

By (10),

[TABLE]

where $\boldsymbol{\nu}-\widetilde{\boldsymbol{\nu}}=U^{\prime}(\boldsymbol{\beta}-\widetilde{\boldsymbol{\beta}})$. Hence,

[TABLE]

where the $\lambda_{k}$’s are the diagonal terms of $\Lambda$.

Since the only term depending on $\boldsymbol{\beta}$ is the second one in the last expression of $\widetilde{L}(\boldsymbol{\beta})$, we define $\widetilde{L}_{Q}(\boldsymbol{\beta})$ appearing in Criterion (11) as follows:

[TABLE]

where

[TABLE]

4.3. Proofs of Propositions 1, 2 and 3 and of Lemma 1

This section contains the proofs of Propositions 1, 2 and 3.

4.3.1. Proof of Proposition 1

We first establish the following lemma for proving Proposition 1.

Lemma 1**.**

$(W_{t}^{\star})$* is an aperiodic Markov process satisfying Doeblin’s condition.*

Proof of Lemma 1.

By (14) and (3), we observe that:

[TABLE]

Thus, $\mathcal{F}_{t-2}=\mathcal{F}_{t-1}^{W^{\star}}:=\sigma(W_{s},s\leq t-1)$. By (1), the distribution of $Y_{t-1}$ conditionally to $\mathcal{F}_{t-2}$ is $\mathcal{P}(\exp(W_{t-1}^{\star}))$. Hence, the distribution of $W_{t}^{\star}$ conditionally to $\mathcal{F}_{t-1}^{W^{\star}}$ is the same as distribution of $W_{t}^{\star}$ conditionally to $W_{t-1}^{\star}$, which means that $(W_{t}^{\star})$ has the Markov property.

Let us now prove that $(W_{t}^{\star})$ is strongly aperiodic which implies that it is aperiodic.

[TABLE]

where the first equality comes from (23) and the last equality comes from (1) since $\mathcal{F}_{t-2}=\mathcal{F}_{t-1}^{W^{\star}}$.

To prove that $(W_{t}^{\star})$ satisfies Doeblin’s condition namely that there exists a probability measure $\nu$ with the property that, for some $m\geq 1$, $\varepsilon>0$ and $\delta>0$,

[TABLE]

for all $x$ in the state space $X$ of $W_{t}^{\star}$ and $B$ in the Borel sets of $X$, we refer the reader to the proof of Proposition 2 in Davis et al. (2003).

∎

Proof of Proposition 1.

For proving Proposition 1, we shall use Theorems 1.3.3 and 1.3.5 of Taniguchi and Kakizawa (2012). In order to apply these theorems it is enough to prove that $(W_{t}^{\star})$ is a strictly stationary and ergodic process since $Y_{t}W_{t}(\beta_{0}^{\star},\gamma_{1})-\exp(W_{t}(\beta_{0}^{\star},\gamma_{1}))$ is a measurable function of $W_{t+1}^{\star},W_{t}^{\star},\dots,W_{2}^{\star}$. Note that the latter fact comes from (14) and (3) for $Y_{t}$ and from (5) with $q=1$ and $p=0$ for $W_{t}$.

In order to prove that $(W_{t}^{\star})$ is a strictly stationary and ergodic process, we have first to prove that $(W_{t}^{\star})$ is an aperiodic Markov process satisfying Doeblin’s condition, see Lemma 1.

The statement of Lemma 1 corresponds to Assertion (iv) of Theorem 16.0.2 of Meyn and Tweedie (1993) which is equivalent to Assertion (i) of this theorem, and implies that $(W_{t}^{\star})$ is uniformly ergodic.

Hence, by Definition (16.6) of uniform ergodicity given in Meyn and Tweedie (1993), there exists a unique stationary invariant measure for $(W_{t}^{\star})$, see also the paragraph below Equation (1.3) of Sandrić (2017) for an additional justification. Combining that existence of a unique stationary invariant measure for $(W_{t}^{\star})$ with the following arguments shows that $(W_{t}^{\star})$ is a strictly stationary process and also an ergodic Markov process.

By Theorem 3.6.3, Corollary 3.6.1 and Definition 3.6.6 of Stout (1974), if the process $(W_{t}^{\star})$ is started with its unique stationary invariant distribution, $(W_{t}^{\star})$ is a strictly stationary process.

By Definition 3.6.8 of Stout (1974), the existence of a unique stationary invariant measure for $(W_{t}^{\star})$ means that $(W_{t}^{\star})$ is an ergodic Markov process, see also the paragraph below (b) (Sandrić, 2017, p. 717).

Finally, by Theorem 3.6.5 of Stout (1974), since $(W_{t}^{\star})$ is an ergodic Markov process and a strictly stationary process, $(W_{t}^{\star})$ is an ergodic and strictly stationary process in the sense of the assumption of Theorem 1.3.5 of Taniguchi and Kakizawa (2012). ∎

4.3.2. Proof of Proposition 2

Note that for all $\gamma_{1}$,

[TABLE]

where the inequality comes from the following inequality $x-\exp(x)\leq-1$, for all $x\in\mathbb{R}$. This inequality is an equality only when $x=0$ which means that $\gamma_{1}=\gamma_{1}^{\star}$.

4.3.3. Proof of Proposition 3

The proof of this proposition comes from Proposition 1 and the stochastic equicontinuity of $n^{-1}L(\beta_{0}^{\star},\gamma_{1})$. Thus, it is enough to prove that there exists a positive $\delta$ such that

[TABLE]

Observe that, by (15),

[TABLE]

Let us first focus on bounding the following expression for $t\geq 2$ (since $W_{1}(\beta_{0}^{\star},\gamma)=\beta_{0}^{\star}$, for all $\gamma$). By (16)

[TABLE]

where we used in the last inequality that for all $x$ and $y$ in $\mathbb{R}$,

[TABLE]

Observing that

[TABLE]

and $|Z_{2}(\gamma_{1})-Z_{2}(\gamma_{2})|\leq\delta[Y_{1}\textrm{e}^{-\beta_{0}^{\star}}+1]$ we get, for $\gamma_{1}$ and $\gamma_{2}$ such that $|\gamma_{1}-\gamma_{2}|\leq\delta$, that

[TABLE]

where $F$ is a measurable function. By (25),

[TABLE]

where the last inequality comes from (27), (26) and (16) and where $G$ is a measurable function. Thus, we get that

[TABLE]

which gives the result by using similar arguments as those given in the proof of Proposition 1 namely that $(Y_{t})$ is strictly stationary and ergodic. By Theorem 1.3.3 of Taniguchi and Kakizawa (2012), $H(Y_{t},Y_{t-1},\dots,Y_{1})$ is strictly stationary and ergodic since $(Y_{t})$ has these properties. Thus, $\mathbb{E}[|H(Y_{t},Y_{t-1},\dots,Y_{1})|]<\infty$, which concludes the proof by Theorem 1.3.5 of Taniguchi and Kakizawa (2012).