A Novel Unconstrained Geometric BINAR(1) Model ^†

Abstract

Modelling the non-stationary unconstrained bivariate integer-valued autoregressive of order 1 (NSUBINAR(1)) model is challenging due to the complex cross-correlation relationship between the counting series. Hence, this paper introduces a novel non-stationary unconstrained BINAR(1) with geometric marginals (NSUBINAR(1)GEOM) based on the assumption that the counting series are both influenced by the same time-dependent explanatory variables. The generalized quasi-likelihood (GQL) estimation method is used to estimate the regression and dependence parameters. Monte Carlo simulations and an application to a real-life accident series data are presented.

1. Introduction

In the literature, several researchers have developed first-order bivariate integer-valued autoregressive (BINAR(1)) models to analyse bivariate time series of counts. Originally, Pedeli and Karlis [1,2] developed two constrained BINAR(1) models with Poisson (CBINAR(1)P) and negative binomial (NB) (CBINAR(1)NB) innovations by extending the classical INAR(1) model of McKenzie [3] based on the binomial thinning mechanism [4]. These two models were developed under stationary moment assumptions only and the cross-correlation between the bivariate series was induced by the correlated Poisson and NB innovations, hence implying a constrained relationship. By the same token, Pedeli and Karlis [5] extended the CBINAR(1)P model to an unconstrained BINAR(1) model with Poisson innovations (UBINAR(1)P) under the same condition of stationarity. In this latter model, the cross-correlation between the series was induced by the correlated Poisson innovation terms and the relationship between the observations of each counting series with previous-lagged observations of the other series.

Likewise, Ristic, Nastic, Jayakumar and Bakouch [6] and Nastic, Ristic and Popovic [7] developed a stationary UBINAR(1) model with geometric marginals (UBINAR(1)GEOM) with independent mixed geometric innovations. Hence, the cross-correlation relationship was induced only by the relation of the current observations with previous-lagged observations of the other series via the negative binomial (NB) thinning operator. Interestingly, Nastic et al. [7] showed in their paper that the UBINAR(1)GEOM yields better AICs than the above over-dispersed BINAR(1) models. However, it is worth mentioning that the UBINAR(1)GEOM model was developed only for stationary time series and hence, cannot be used to analyse non-stationary real-life over-dispersed series.

As for non-stationary time series of counts, few BINAR(1) models have been developed. Mamodekhan, Sunecher and Jowaheer [8] developed a CBINAR(1) model with Poisson innovations under non-stationarity assumption (NSCBINAR(1)P) induced by time-dependent explanatory variables. In a similar context, Sunecher, Mamodekhan and Jowaheer [9] developed a non-stationary CBINAR(1) model with NB innovations for over-dispersed time series of counts (NSCBINAR(1)NB). However, an unconstrained BINAR(1) model under the non-stationarity context that can model over-dispersed counting series has not yet been developed in the literature. Hence, in this paper we propose to develop a non-stationary unconstrained BINAR(1) model with geometric innovations (NSUBINAR(1)GEOM) similar to Pedeli and Karlis [5]. However, such model development poses some computational challenges in estimating the unknown model parameters as it is rather difficult to specify the joint generating function [10].

In the paper by Pedeli and Karlis [2,5], the conditional maximum likelihood (CML) approach was compared with the method of moments (MoM), where the authors concluded that CML yields far better estimates than MoM, but at huge computational costs. A similar conclusion was drawn by Nastic et al. [7] who used the least-square (LS) technique, as an alternative to CML, which omits the likelihood function. In the same way, due to the computational challenges of the CML, Mamodekhan et al. [8] and Sunecher et al. [9] developed the generalized quasi-likelihood (GQL) approach in the non-stationary bivariate context. Mamodekhan et al. [8] compared the GQL with CML, where it was shown that GQL yields asymptotically equally efficient estimates as CML. Hence, based on the above findings, the unknown parameters of the NSUBINAR(1)GEOM will be estimated using the GQL.

The organization of the paper is as follows: In the next section, the NSUBINAR(1) model is developed. In Section 3, the GQL approach is developed under the non-stationarity bivariate context. Section 4 focuses on the simulation part where BINAR(1) data with geometric marginals are generated and the GQL approach is used to estimate the model parameters. In Section 5, the model is applied to the accident data in Mauritius. The conclusion is provided in the final section.

2. The Non-Stationary Unconstrained BINAR(1) with Geometric Marginals (NSUBINAR(1)GEOM)

The UBINAR(1) model is specified as:

$$Y_t^{\left[1\right]}={\rho }_{11}\ast Y_{t-1}^{\left[1\right]}+{\rho }_{12}\ast Y_{t-1}^{\left[2\right]}+R_t^{\left[1\right]}$$

(1)

$$Y_t^{\left[2\right]}={\rho }_{21}\ast Y_{t-1}^{\left[1\right]}+{\rho }_{22}\ast Y_{t-1}^{\left[2\right]}+R_t^{\left[2\right]}$$

(2)

based on the following assumptions:

(a)

$Y_t^{\left[k\right]}$ is geometric such that $Y_t^{\left[k\right]}\sim \mathrm{Geom}\left(\frac{{\mu }_t^{\left[k\right]}}{1+{\mu }_t^{\left[k\right]}}\right)$. Hence, E($Y_t^{\left[k\right]})={\mu }_t^{\left[k\right]}$ and Var($Y_t^{\left[k\right]})={\mu }_t^{\left[k\right]}(1+{\mu }_t^{\left[k\right]})$, where ${\mu }_t^{\left[k\right]}=exp\left({\mathit{x}}_t^{{}^{\prime }}{\mathit{\beta }}^{\left[k\right]}\right)$ with ${\mathit{x}}_t={[{x_t}_1,{x_t}_2,\dots ,{x_t}_j,\dots ,{x_t}_p]}^{{}^{\prime }}$ is a $(p\times 1)$ vector of covariates influencing both $Y_t^{\left[1\right]}$ and $Y_t^{\left[2\right]}$, with corresponding regression coefficients ${\mathit{\beta }}^{\left[k\right]}={[{\beta }_1^{\left[k\right]},{\beta }_2^{\left[k\right]},\dots ,{\beta }_j^{\left[k\right]},\dots ,{\beta }_p^{\left[k\right]}]}^{{}^{\prime }}$ for $t=1,2,\dots ,T$ and $k\in \{1,2\}$.

(b)

∗ is the binomial thinning operator [4] such that ${\rho }_{ij}\ast Y_{t-1}^{\left[j\right]}={\sum }_{m=1}^{Y_{t-1}^{\left[j\right]}}{Z}_m$ with $Z_m\sim \mathrm{Geom}\left(\frac{{\rho }_{ij}}{1+{\rho }_{ij}}\right)$. Hence, $E({\rho }_{ij}\ast Y_{t-1}^{\left[j\right]})={\rho }_{ij}E\left(Y_{t-1}^{\left[j\right]}\right)$ and $\mathrm{Var}({\rho }_{ij}\ast Y_{t-1}^{\left[j\right]})=$${\rho }_{ij}(1+{\rho }_{ij})E\left(Y_{t-1}^{\left[j\right]}\right)+{\rho }_{ij}^2\mathrm{Var}\left(Y_{t-1}^{\left[j\right]}\right)$.

(c)

$$\mathrm{Corr}(R_t^{\left[1\right]},R_{t^{{}^{\prime }}}^{\left[2\right]})=\left\{\begin{array}{cc}{\kappa }_{12,t}\hfill & t=t^{{}^{\prime }},\hfill \\ 0\hfill & t\ne t^{{}^{\prime }}.\hfill \end{array}\right.$$

$$\begin{array}{cc}\hfill E\left(Y_t^{\left[1\right]}\right)& =E({\rho }_{11}\ast Y_{t-1}^{\left[1\right]})+E({\rho }_{12}\ast Y_{t-1}^{\left[2\right]})+E\left(R_t^{\left[1\right]}\right)\hfill \\ \hfill {\mu }_t^{\left[1\right]}& ={\rho }_{11}{\mu }_{t-1}^{\left[1\right]}+{\rho }_{12}{\mu }_{t-1}^{\left[2\right]}+E\left(R_t^{\left[1\right]}\right).\hfill \end{array}$$

(3)

Re-arranging Equation (3), we have

$$E\left(R_t^{\left[1\right]}\right)={\lambda }_t^{\left[1\right]}={\mu }_t^{\left[1\right]}-{\rho }_{11}{\mu }_{t-1}^{\left[1\right]}-{\rho }_{12}{\mu }_{t-1}^{\left[2\right]}.$$

(4)

Similarly,

$$E\left(R_t^{\left[2\right]}\right)={\lambda }_t^{\left[2\right]}={\mu }_t^{\left[2\right]}-{\rho }_{21}{\mu }_{t-1}^{\left[1\right]}-{\rho }_{22}{\mu }_{t-1}^{\left[2\right]}.$$

(5)

$$\begin{array}{cc}\hfill \mathrm{Var}\left(Y_t^{\left[1\right]}\right)& =\mathrm{Var}({\rho }_{11}\ast Y_{t-1}^{\left[1\right]}+{\rho }_{12}\ast Y_{t-1}^{\left[2\right]}+R_t^{\left[1\right]})\hfill \\ \hfill \mathrm{Var}\left(Y_t^{\left[1\right]}\right)& ={\rho }_{11}(1+{\rho }_{11})E\left(Y_{t-1}^{\left[1\right]}\right)+{\rho }_{11}^2\mathrm{Var}\left(Y_{t-1}^{\left[1\right]}\right)+{\rho }_{12}(1+{\rho }_{12})E\left(Y_{t-1}^{\left[2\right]}\right)+{\rho }_{12}^2\mathrm{Var}\left(Y_{t-1}^{\left[2\right]}\right)\hfill \\ \hfill & +2{\rho }_{11}{\rho }_{12}\mathrm{Cov}(Y_{t-1}^{\left[1\right]},Y_{t-1}^{\left[2\right]})+\mathrm{Var}\left(R_t^{\left[1\right]}\right)\hfill \\ \hfill {\mu }_t^{\left[1\right]}+{{\mu }_t^{\left[1\right]}}^2& ={\rho }_{11}(1+{\rho }_{11}){\mu }_{t-1}^{\left[1\right]}+{\rho }_{11}^2({\mu }_{t-1}^{\left[1\right]}+{{\mu }_{t-1}^{\left[1\right]}}^2)+{\rho }_{12}(1+{\rho }_{12}){\mu }_{t-1}^{\left[2\right]}+{\rho }_{12}^2({\mu }_{t-1}^{\left[2\right]}+{{\mu }_{t-1}^{\left[2\right]}}^2)\hfill \\ \hfill & +2{\rho }_{11}{\rho }_{12}\mathrm{Cov}(Y_{t-1}^{\left[1\right]},Y_{t-1}^{\left[2\right]})+\mathrm{Var}\left(R_t^{\left[1\right]}\right)\hfill \end{array}$$

(6)

Re-arranging Equation (6), we have

$$\begin{array}{cc}\hfill \mathrm{Var}\left(R_t^{\left[1\right]}\right)& ={\mu }_t^{\left[1\right]}(1+{\mu }_t^{\left[1\right]})-{\rho }_{11}(1+{\rho }_{11}){\mu }_{t-1}^{\left[1\right]}-{\rho }_{11}^2{\mu }_{t-1}^{\left[1\right]}(1+{\mu }_{t-1}^{\left[1\right]})-{\rho }_{12}(1+{\rho }_{12}){\mu }_{t-1}^{\left[2\right]}\hfill \\ \hfill & -{\rho }_{12}^2{\mu }_{t-1}^{\left[2\right]}(1+{\mu }_{t-1}^{\left[2\right]})-2{\rho }_{11}{\rho }_{12}\mathrm{Cov}(Y_{t-1}^{\left[1\right]},Y_{t-1}^{\left[2\right]})\hfill \end{array}$$

(7)

and similarly,

$$\begin{array}{cc}\hfill \mathrm{Var}\left(R_t^{\left[2\right]}\right)& ={\mu }_t^{\left[2\right]}(1+{\mu }_t^{\left[2\right]})-{\rho }_{21}(1+{\rho }_{21}){\mu }_{t-1}^{\left[1\right]}-{\rho }_{21}^2{\mu }_{t-1}^{\left[1\right]}(1+{\mu }_{t-1}^{\left[1\right]})-{\rho }_{22}(1+{\rho }_{22}){\mu }_{t-1}^{\left[2\right]}\hfill \\ \hfill & -{\rho }_{22}^2{\mu }_{t-1}^{\left[2\right]}(1+{\mu }_{t-1}^{\left[2\right]})-2{\rho }_{21}{\rho }_{22}\mathrm{Cov}(Y_{t-1}^{\left[1\right]},Y_{t-1}^{\left[2\right]}),\hfill \end{array}$$

(8)

The above moments clearly indicate that the marginal distribution of $R_t^{\left[k\right]}$ is rather complex to derive. To facilitate the derivation of the cross-covariances, we write Equations (1) and (2) in vector form as follows:

$${\mathit{Y}}_t=\mathit{A}\ast {\mathit{Y}}_{t-1}+{\mathit{R}}_t$$

(9)

with ${\mathit{Y}}_{\mathit{t}}={[Y_t^{\left[1\right]},Y_t^{\left[2\right]}]}^{{}^{\prime }}$, ${\mathit{R}}_{\mathit{t}}={[R_t^{\left[1\right]},R_t^{\left[2\right]}]}^{{}^{\prime }}$, $\mathit{A}=\left(\begin{array}{cc}{\rho }_{11}& {\rho }_{12}\\ {\rho }_{21}& {\rho }_{22}\end{array}\right)$.

Assuming ${\mathsf{\Sigma }}_{h,t}=\left[\begin{array}{cc}\mathrm{Cov}(Y_t^{\left[1\right]},Y_{t+h}^{\left[1\right]})& \mathrm{Cov}(Y_{t+h}^{\left[1\right]},Y_t^{\left[2\right]})\\ \mathrm{Cov}(Y_t^{\left[1\right]},Y_{t+h}^{\left[2\right]})& \mathrm{Cov}(Y_t^{\left[2\right]},Y_{t+h}^{\left[2\right]})\end{array}\right]$, from Pedeli and Karlis [5] and Ristic et al. [6], it was shown that ${\mathsf{\Sigma }}_{h,t}={\mathit{A}}^h{\mathsf{\Sigma }}_{0,t}$ and hence

$$\begin{array}{cc}\hfill \mathrm{Cov}(Y_t^{\left[1\right]},Y_t^{\left[2\right]})& =\mathrm{Cov}({\rho }_{11}\ast Y_{t-1}^{\left[1\right]}+{\rho }_{12}\ast Y_{t-1}^{\left[2\right]}+R_t^{\left[1\right]},{\rho }_{21}\ast Y_{t-1}^{\left[1\right]}+{\rho }_{22}\ast Y_{t-1}^{\left[2\right]}+R_t^{\left[2\right]})\hfill \\ \hfill & =({\rho }_{11}{\rho }_{22}+{\rho }_{12}{\rho }_{21})\mathrm{Cov}(Y_{t-1}^{\left[1\right]},Y_{t-1}^{\left[2\right]})+{\rho }_{11}{\rho }_{21}({\mu }_{t-1}^{\left[1\right]}+{{\mu }_{t-1}^{\left[1\right]}}^2)\hfill \\ \hfill & +{\rho }_{22}{\rho }_{12}({\mu }_{t-1}^{\left[2\right]}+{{\mu }_{t-1}^{\left[2\right]}}^2)+\left[{\kappa }_{12,t}\sqrt{\mathrm{Var}\left(R_t^{\left[1\right]}\right)}\sqrt{\mathrm{Var}\left(R_t^{\left[2\right]}\right)}\right].\hfill \end{array}$$

(10)

Note that if $\alpha ={\rho }_{11}={\rho }_{12}$, $\gamma ={\rho }_{21}={\rho }_{22}$, $\mu ={\mu }_{t-1}^{\left[1\right]}={\mu }_{t-1}^{\left[2\right]}$ and ${\kappa }_{12,t}=0$, Equation (10) simply reduces to $\mathrm{Cov}(Y_t^{\left[1\right]},Y_t^{\left[2\right]})=\frac{2\alpha \gamma }{1-2\alpha \gamma }\mu (1+\mu )$, which is the same as in Ristic et al. [6].

3. Estimation Method

The GQL equation to estimate the regression parameters is specified as:

$${\mathit{D}}_{\mathit{\beta }}^{{}^{\prime }}{\mathsf{\Sigma }}_{\mathit{\beta }}^{-1}(\mathit{f}-\mathit{\mu })=0$$

(11)

with score vector $\mathit{f}=[{\mathit{f}}_1,{\mathit{f}}_2,\dots ,{\mathit{f}}_t,\dots ,{\mathit{f}}_{t+h},\dots ,{\mathit{f}}_T]$ with ${\mathit{f}}_t={[Y_t^{\left[1\right]},Y_t^{\left[2\right]}]}^{{}^{\prime }}$ and $\mathit{\mu }=[{\mathit{\mu }}_1,{\mathit{\mu }}_2,\dots ,{\mathit{\mu }}_t,\dots ,{\mathit{\mu }}_T]$ with corresponding mean ${\mathit{\mu }}_t={[{\mu }_t^{\left[1\right]},{\mu }_t^{\left[2\right]}]}^{{}^{\prime }}$ for $t=1,2,\dots ,T$.

The covariance matrix ${\mathsf{\Sigma }}_{\mathit{\beta }}$ is a $(2T\times 2T)$. The derivative matrix ${\mathit{D}}_{\mathit{\beta }}$ is denoted by ${\mathit{D}}_{\mathit{\beta }}={[{\mathit{D}}_1,{\mathit{D}}_2,\dots ,{\mathit{D}}_t,\dots ,{\mathit{D}}_T]}^{\prime }$ with

$${\mathit{D}}_t={\left(\begin{array}{cc}\frac{\partial {\mu }_t^{\left[1\right]}}{\partial {\mathit{\beta }}^{\left[1\right]}}& 0\\ 0& \frac{\partial {\mu }_t^{\left[2\right]}}{\partial {\mathit{\beta }}^{\left[2\right]}}\end{array}\right)}_{2p\times 2}$$

where $\frac{\partial {\mu }_t^{\left[k\right]}}{\partial {\beta }_j^{\left[k\right]}}={\mu }_t^{\left[k\right]}{x}_{tj}^{\prime }$.

The Newton–Raphson iterative technique is used to estimate the regression parameters as follows:

$$\left(\begin{array}{c}{\widehat{\mathit{\beta }}}_{r+1}^{\left[1\right]}\\ {\widehat{\mathit{\beta }}}_{r+1}^{\left[2\right]}\end{array}\right)=\left(\begin{array}{c}{\widehat{\mathit{\beta }}}_r^{\left[1\right]}\\ {\widehat{\mathit{\beta }}}_r^{\left[2\right]}\end{array}\right)+{\left[{\mathit{D}}_{\mathit{\beta }}^{{}^{\prime }}{\mathsf{\Sigma }}_{\mathit{\beta }}^{-1}{\mathit{D}}_{\mathit{\beta }}\right]}_r^{-1}{\left[{\mathit{D}}_{\mathit{\beta }}^{{}^{\prime }}{\mathsf{\Sigma }}_{\mathit{\beta }}^{-1}(\mathit{f}-\mathit{\mu })\right]}_r$$

(12)

where ${\widehat{\mathit{\beta }}}_r^{\left[k\right]}$ are the estimates at the rth iteration and ${[.]}_r$ are the values of the expression at the rth iteration.

For an initial value of $[{\widehat{\rho }}_{11},{\widehat{\rho }}_{12},{\widehat{\rho }}_{21},{\widehat{\rho }}_{22}$, ${\widehat{\kappa }}_{12,t},\mathit{\beta }]$, we solve the iterative Equation (12) until convergence. These estimates are consistent and under mild regulatory conditions, ${(\widehat{\mathit{\beta }}-\mathit{\beta })}^{\prime }$ is asymptotically normal with a mean of 0 and a covariance matrix of ${\left[{\mathit{D}}_{\mathit{\beta }}^{{}^{\prime }}{\mathsf{\Sigma }}_{\mathit{\beta }}^{-1}{\mathit{D}}_{\mathit{\beta }}\right]}^{-1}$ as shown in [8,9,11].

A second GQL is specified to estimate the dependence parameter $\mathit{\psi }=[{\rho }_{11},{\rho }_{12},{\rho }_{21},{\rho }_{22}]$ as follows:

$${\mathit{D}}_{\mathit{\psi }}^{{}^{\prime }}{\mathsf{\Sigma }}_{\mathit{\psi }}^{-1}({\mathit{Y}}_{\mathit{\psi }}-{\mu }_{\mathit{\psi }})=0,$$

(13)

with ${\mathit{Y}}_{\mathit{\psi }}=[Y_1^{\left[1\right]}{Y}_1^{\left[2\right]}|Y_0^{\left[1\right]},Y_0^{\left[2\right]},Y_2^{\left[1\right]}{Y}_2^{\left[2\right]}|Y_1^{\left[1\right]},Y_1^{\left[2\right]},\dots ,Y_t^{\left[1\right]}{Y}_t^{\left[2\right]}|Y_{t-1}^{\left[1\right]},Y_{t-1}^{\left[2\right]},\dots ,Y_T^{\left[1\right]}{Y}_T^{\left[2\right]}{|Y_{T-1}^{\left[1\right]},Y_{T-1}^{\left[2\right]}]}_{T\times 1}^{{}^{\prime }}$ and ${\mathit{\mu }}_{\mathit{\psi }}=E\left({\mathit{Y}}_{\mathit{\psi }}\right)$.

The $(T\times T)$ covariance matrix ${\mathsf{\Sigma }}_{\mathit{\psi }}$ comprises of Var$(Y_t^{\left[1\right]}{Y}_t^{\left[2\right]}\mid Y_{t-1}^{\left[1\right]},Y_{t-1}^{\left[2\right]})$ along the diagonal and Cov$(Y_t^{\left[1\right]}{Y}_t^{\left[2\right]}{Y}_{t+h}^{\left[1\right]}{Y}_{t+h}^{\left[2\right]}\mid Y_{t-1}^{\left[1\right]},Y_{t-1}^{\left[2\right]},Y_{t+h-1}^{\left[1\right]},Y_{t+h-1}^{\left[2\right]})$ in the off-diagonal entries. All the entries are of higher-order moments and hence, the ’working’ multivariate normality assumption structure is used to compute these entries as in [9,12].

As for the $(T\times 4)$ derivative matrix ${\mathit{D}}_{\mathit{\psi }}$,

$E(Y_t^{\left[1\right]}{Y}_t^{\left[2\right]}\mid Y_{t-1}^{\left[1\right]},Y_{t-1}^{\left[2\right]})$

$$\begin{array}{cc}\hfill & =[{\kappa }_{12,t}[{\mu }_t^{\left[1\right]}(1+{\mu }_t^{\left[1\right]})-{\rho }_{11}(1+{\rho }_{11}){\mu }_{t-1}^{\left[1\right]}-{\rho }_{11}^2{\mu }_{t-1}^{\left[1\right]}(1+{\mu }_{t-1}^{\left[1\right]})-{\rho }_{12}(1+{\rho }_{12}){\mu }_{t-1}^{\left[2\right]}\hfill \\ \hfill & -{\rho }_{12}^2{\mu }_{t-1}^{\left[2\right]}(1+{\mu }_{t-1}^{\left[2\right]}){]}^{\frac{1}{2}}[{\mu }_t^{\left[2\right]}(1+{\mu }_t^{\left[2\right]})-{\rho }_{21}(1+{\rho }_{21}){\mu }_{t-1}^{\left[1\right]}-{\rho }_{21}^2{\mu }_{t-1}^{\left[1\right]}(1+{\mu }_{t-1}^{\left[1\right]})\hfill \\ \hfill & -{\rho }_{22}(1+{\rho }_{22}){\mu }_{t-1}^{\left[2\right]}-{\rho }_{22}^2{\mu }_{t-1}^{\left[2\right]}(1+{\mu }_{t-1}^{\left[2\right]}){]}^{\frac{1}{2}}+({\rho }_{11}{Y}_{t-1}^{\left[1\right]}+{\rho }_{12}{Y}_{t-1}^{\left[2\right]}+{\mu }_t^{\left[1\right]}-{\rho }_{11}{\mu }_{t-1}^{\left[1\right]}\hfill \\ \hfill & -{\rho }_{12}{\mu }_{t-1}^{\left[2\right]})\times ({\rho }_{21}{Y}_{t-1}^{\left[1\right]}+{\rho }_{22}{Y}_{t-1}^{\left[2\right]}+{\mu }_t^{\left[2\right]}-{\rho }_{21}{\mu }_{t-1}^{\left[1\right]}-{\rho }_{22}{\mu }_{t-1}^{\left[2\right]})]\hfill \end{array}$$

(14)

As for the estimates of ${\kappa }_{12,t}$, they are estimated using the method of moments as follows:

$${\kappa }_{12,t}=\frac{\widetilde{\mathrm{Cov}}(Y_t^{\left[1\right]},Y_t^{\left[2\right]})-({\widehat{\rho }}_{11}{\widehat{\rho }}_{22}+{\widehat{\rho }}_{12}{\widehat{\rho }}_{21})\widetilde{\mathrm{Cov}}(Y_{t-1}^{\left[1\right]},Y_{t-1}^{\left[2\right]})-{\widehat{\rho }}_{11}{\widehat{\rho }}_{21}({\widehat{\mu }}_{t-1}^{\left[1\right]}+{\widehat{{\mu }_{t-1}^{\left[1\right]}}}^2)-{\widehat{\rho }}_{22}{\widehat{\rho }}_{12}({\widehat{\mu }}_{t-1}^{\left[2\right]}+{\widehat{{\mu }_{t-1}^{\left[2\right]}}}^2)}{\sqrt{\mathrm{Var}\left(R_t^{\left[1\right]}\right)}\sqrt{\mathrm{Var}\left(R_t^{\left[2\right]}\right)}}$$

(15)

where ${\widehat{\mu }}_0^{\left[k\right]}={\widehat{\mu }}_1^{\left[k\right]}$, $\widetilde{\mathrm{Cov}}(Y_t^{\left[1\right]},Y_t^{\left[2\right]})$=$\frac{1}{T}{\sum }_{t=1}^T(y_t^{\left[1\right]}-{\widehat{\mu }}_t^{\left[1\right]})(y_t^{\left[2\right]}-{\widehat{\mu }}_t^{\left[2\right]})$ and

$\widetilde{\mathrm{Cov}}(Y_{t-1}^{\left[1\right]},Y_{t-1}^{\left[2\right]})$=$\frac{1}{T-1}{\sum }_{t=2}^T(y_{t-1}^{\left[1\right]}-{\widehat{\mu }}_{t-1}^{\left[1\right]})(y_{t-1}^{\left[2\right]}-{\widehat{\mu }}_{t-1}^{\left[2\right]})$.

The Newton–Raphson iteration for the second GQL yields

$$\left(\begin{array}{c}{\widehat{\mathit{\psi }}}_{r+1}\end{array}\right)=\left(\begin{array}{c}{\widehat{\mathit{\psi }}}_r\end{array}\right)+{\left[{\mathit{D}}_{\mathit{\psi }}^{{}^{\prime }}{{\widehat{\mathsf{\Sigma }}}_{\mathit{\psi }}}^{-1}{\mathit{D}}_{\mathit{\psi }}\right]}_r^{-1}{\left[{\mathit{D}}_{\mathit{\psi }}^{{}^{\prime }}{\widehat{\mathsf{\Sigma }}}_{\mathit{\psi }}^{-1}({\mathit{Y}}_{\mathit{\psi }}-{\mathit{\mu }}_{\mathit{\psi }})\right]}_r$$

(16)

where ${\widehat{\mathit{\psi }}}_r$ are the estimates at the rth iteration and ${[.]}_r$ are the values of the expression at the rth iteration.

The regression estimates $\widehat{\mathit{\beta }}$ obtained from Equation (12) are used to solve Equation (16) until convergence. These updated values of $\widehat{\mathit{\psi }}$ are in turn replaced in Equation (12) to obtain a new set of regression parameters, which are again used to obtain a new set of dependence parameters and this cycle continues until convergence of the two sets of parameters. These estimates are consistent and under mild regulatory conditions, ${(\widehat{\mathit{\psi }}-\mathit{\psi })}^{{}^{\prime }}$ is asymptotically normal with a mean of 0 and a covariance matrix ${\left[{\mathit{D}}_{\mathit{\psi }}^{{}^{\prime }}{\mathsf{\Sigma }}_{\mathit{\psi }}^{-1}{\mathit{D}}_{\mathit{\psi }}\right]}^{-1}$.

Forecasting Equations

The forecasting equations are derived as follows:

Given $Y_t^{\left[k\right]}$, the forecasting function is expressed as

$$E\left(Y_{t+1}^{\left[1\right]}\right|Y_t^{\left[1\right]},Y_t^{\left[2\right]})={\widehat{\mu }}_{t+1}^{\left[1\right]}+{\widehat{\rho }}_{11}(Y_t^{\left[1\right]}-{\widehat{\mu }}_t^{\left[1\right]})+{\widehat{\rho }}_{12}(Y_t^{\left[2\right]}-{\widehat{\mu }}_t^{\left[2\right]}),$$

(17)

$$E\left(Y_{t+1}^{\left[2\right]}\right|Y_t^{\left[1\right]},Y_t^{\left[2\right]})={\widehat{\mu }}_{t+1}^{\left[2\right]}+{\widehat{\rho }}_{21}(Y_t^{\left[1\right]}-{\widehat{\mu }}_t^{\left[1\right]})+{\widehat{\rho }}_{22}(Y_t^{\left[2\right]}-{\widehat{\mu }}_t^{\left[2\right]})$$

(18)

4. Simulation Study

In this section we generate BINAR(1) time series data with geometric marginals under the following time-varying covariate design:

$$x_{t1}=\left\{\begin{array}{cc}-\cos \left(2\pi t\right)+0.01\hfill & (t=1,\dots ,T/4)\hfill \\ \sin \left(2\pi t\right)+0.05\hfill & (t=(T/4)+1,\dots ,3T/4)\hfill \\ \cos \left(2\pi t\right)+0.10\hfill & (t=(3T/4)+1,\dots ,T)\hfill \end{array}\right.$$

$$x_{t2}=\left\{\begin{array}{cc}(1/t)\hfill & (t=1,\dots ,T/4)\hfill \\ (-1/t)\hfill & (t=(T/4)+1,\dots ,3T/4)\hfill \\ t\hfill & (t=(3T/4)+1,\dots ,T)\hfill \end{array}\right.$$

where ${\mu }_t^{\left[k\right]}=exp(x_{t1}{\beta }_1^{\left[k\right]}+x_{t2}{\beta }_2^{\left[k\right]})$. Assuming $[{\rho }_{12},{\rho }_{21}]=[0.5,0.5]$, $[{\rho }_{11},{\rho }_{22}]=$$[0.9,0.9],[0.3,0.9],[0.3,0.3]$, ${\beta }^{\left[1\right]}=0.5$ and ${\beta }^{\left[2\right]}=0.9$ for $t=1,2,\dots ,T=100,500,1000$, we generate $R_t^{\left[k\right]}$ using the inverse transformation method as in [13]. A total of 5000 Monte Carlo replications are made under the above combinations and the simulated mean estimates are shown below:

From

Table 1. GQL estimates of the parameters and standard errors under the non-stationary geometric BINAR(1) model.

${\mathit{\rho }}_{11}$	${\mathit{\rho }}_{22}$	T	Methods	${\widehat{\mathit{\beta }}}_1^{\left[1\right]}$	${\widehat{\mathit{\beta }}}_2^{\left[1\right]}$	${\widehat{\mathit{\beta }}}_1^{\left[2\right]}$	${\widehat{\mathit{\beta }}}_2^{\left[2\right]}$	${\widehat{\mathit{\rho }}}_{11}$	${\widehat{\mathit{\rho }}}_{22}$	${\widehat{\mathit{\rho }}}_{12}$	${\widehat{\mathit{\rho }}}_{21}$	${\widehat{\mathit{\kappa }}}_{12,1}$
0.9	0.9	100	GQL	0.4823	0.4870	0.8876	0.8840	0.8819	0.8847	0.4859	0.4847	0.9815
				(0.0910)	(0.0931)	(0.0977)	(0.0946)	(0.1163)	(0.1115)	(0.1125)	(0.1153)
		500	GQL	0.4923	0.4914	0.8950	0.8917	0.8942	0.8920	0.4945	0.4940	0.9909
				(0.0517)	(0.0512)	(0.0547)	(0.0585)	(0.0649)	(0.0625)	(0.0623)	(0.0630)
		1000	GQL	0.4991	0.4988	0.8960	0.8981	0.8984	0.8994	0.5004	0.5006	0.9979
				(0.0122)	(0.0171)	(0.0194)	(0.0146)	(0.0271)	(0.0224)	(0.0209)	(0.0237)
0.3	0.9	100	GQL	0.4894	0.4891	0.8870	0.8840	0.2856	0.8868	0.4899	0.4888	0.9826
				(0.0969)	(0.0977)	(0.0915)	(0.0959)	(0.1133)	(0.1172)	(0.1128)	(0.1119)
		500	GQL	0.4926	0.4927	0.8944	0.8959	0.2915	0.8928	0.4965	0.4940	0.9918
				(0.0526)	(0.0556)	(0.0507)	(0.0512)	(0.0681)	(0.0694)	(0.0671)	(0.0631)
		1000	GQL	0.4995	0.4975	0.8969	0.8994	0.2988	0.8963	0.5008	0.5011	0.9995
				(0.0137)	(0.0175)	(0.0111)	(0.0141)	(0.0213)	(0.0293)	(0.0251)	(0.0221)
0.3	0.3	100	GQL	0.4823	0.4870	0.8804	0.8896	0.2812	0.2834	0.4854	0.4835	0.9819
				(0.0981)	(0.0911)	(0.0931)	(0.0928)	(0.1169)	(0.1160)	(0.1135)	(0.1197)
		500	GQL	0.4929	0.4942	0.8910	0.8935	0.2964	0.2931	0.4920	0.4915	0.9913
				(0.0594)	(0.0589)	(0.0562)	(0.0580)	(0.0614)	(0.0621)	(0.0677)	(0.0681)
		1000	GQL	0.4956	0.4992	0.8987	0.8990	0.2988	0.2980	0.5004	0.5001	0.9966
				(0.0152)	(0.0125)	(0.0135)	(0.0133)	(0.0241)	(0.0208)	(0.0219)	(0.0211)

, we observe that the GQL estimates are consistent and that the cross-correlation parameter ${\kappa }_{12,1}$ is close to unity. In addition, as the time points increase, we notice a decrease in the standard errors, with GQL yielding low standard errors as also demonstrated in [14]. Some details on the number of non-convergent simulations include: for ${\rho }_{11}={\rho }_{22}=0.9$ under GQL, around 360 simulations failed for $T=100$, 300 for $T=500$ and 220 for $T=1000$. For ${\rho }_{11}=0.3$ and ${\rho }_{22}=0.9$, around 340 GQL simulations failed for $T=100$, 275 for $T=500$ and 190 for $T=1000$. However, when ${\rho }_{11}={\rho }_{22}=0.3$, the GQL algorithms failed in 315 simulations for $T=100$, 215 for $T=500$ and 170 for $T=1000$. The failures were mainly due to either an ill-conditioned covariance matrix or the Hessian structure in Equation (12). Hence, it is concluded from this section that GQL yields far superior estimates than GLS and GMM, and constitutes of a slightly better non-convergent computational problem.

5. Analysing the Time Series of Day and Night Road Accidents in Mauritius

In this section we analyse the monthly day and night accident series data in Mauritius collected from January 2011 to January 2020 that connects the capital city of Port-Louis and the tourist zone Grand-Bay, Mauritius, totalling 109 bivariate time series data. With a sample cross-correlation of 0.3267, it is rationale to believe that there exists a cross-correlation between the two series since both sets of data were collected on the same route. The summary statistics illustrate that day and night accidents have means (variance) of 8.6422 (27.5652) and 4.2110 (13.1310), respectively. Given the significant over-dispersion, the NSUBINAR(1)GEOM in Section 2 is applied to analyse the time series. The covariates we consider are: number of speed cameras (SC) in this area, the number of police officers deployed on street patrol in the different police stations in this area (PO), the number of times the streets in the area have been re-maintained during the years (NS) and number of roundabouts (RA) from Port-Louis to Grand-Bay.

Table 2. Monthly day and night accidents: GQL estimates of the regression parameters.

Model	Time Series	Intercept	NS	SC	PO	RA	$\widehat{\mathit{c}}$
	Day Accidents	2.5353	−0.0815	−0.0942	−0.0934	0.0760	1.8475
	s.e	(0.0633)	(0.0396)	(0.0456)	(0.0397)	(0.0289)	(0.0948)
NSCBINAR(1)NB	Night Accidents	0.9272	−0.0790	−0.0671	−0.0824	0.0943	0.8965
	s.e	(0.0742)	(0.0245)	(0.0213)	(0.0329)	(0.0330)	(0.0980)
	Day Accidents	2.4445	−0.0824	−0.0952	−0.0955	0.0668
	s.e	(0.0560)	(0.0273)	(0.0314)	(0.0215)	(0.0164)
NSUBINAR(1)GEOM	Night Accidents	0.9106	−0.0714	−0.0572	−0.0747	0.0817
	s.e	(0.0710)	(0.0180)	(0.0157)	(0.0260)	(0.0205)
	Day Accidents	2.413	−0.0929	−0.0866	−0.0961	0.0852
	s.e	(0.0958)	(0.0353)	(0.0385)	(0.0419)	(0.0342)
NSCBINAR(1)GEOM	Night Accidents	0.9623	−0.0876	−0.0894	−0.0785	0.0899
	s.e	(0.0952)	(0.0367)	(0.0387)	(0.0375)	(0.0345)

and

Table 3. Monthly day and night accidents: GQL estimates of the dependence parameters.

Model	Time Series	$\widehat{\mathit{\rho }}-\mathbf{Serial}$	$\widehat{\mathit{\rho }}-\mathbf{Cross}$	${\widehat{\mathit{\kappa }}}_{12,1}$
	Day Accidents	0.2620		0.0065
	s.e
NSCBINAR(1)NB	Night Accidents	0.2941
	s.e
	Day Accidents	0.2748	0.0728	0.0025
	s.e	(0.0346)	(0.0234)
NSUBINAR(1)GEOM	Night Accidents	0.2438	0.0563
	s.e	(0.0388)	(0.0191)
	Day Accidents	0.2145		0.0034
	s.e
NSCBINAR(1)GEOM	Night Accidents	0.2442
	s.e

show the regression estimates and the serial and dependence estimates of the in-sample accident data from January 2011 to August 2019, totalling 104 paired observations, while the out-sample data from September 2019 to January 2020 were used to validate the model. Note, the in-sample data were also analysed using the NSCBINAR(1)NB from [9] and the estimates were compared with the NSCBINAR(1)GEOM and NSUBINAR(1)GEOM.

From Table 2, using the NSUBINAR(1)GEOM, we can note an expected decrease in the number of day and night accidents by 9 and 5%, respectively, if there is an installation of an additional speed camera along the motorway M2. As more police patrols are re-enforced in the area, the number of accidents is expected to decrease by 9% during the day and 7% during the night. Similarly, better road maintenance contributes to a decrease in day accidents by 7% and 6% in night accidents. Roundabout construction must be carefully monitored as this factor leads to an expected increase in the number of day accidents by 6% and 8% during the night. In comparison with the NSCBINAR(1)NB, the signs estimated effects are the same as in NSUBINAR(1)GEOM with little fluctuation in the estimates and their corresponding standard errors, but far better than NSCBINAR(1)GEOM. Using the corresponding forecasting in Equations (17) and (18) and in [9], and based on the out-sample observations from September 2019 to January 2020,

Table 4. RMSE and MAD for the one step-ahead forecast for the number of monthly day and night accidents.

Model	RMSE	RMSE	MAD	MAD
	${\mathit{Y}}_{\mathit{t}}^{\left[\mathbf{1}\right]}$	${\mathit{Y}}_{\mathit{t}}^{\left[\mathbf{2}\right]}$	${\mathit{Y}}_{\mathit{t}}^{\left[\mathbf{1}\right]}$	${\mathit{Y}}_{\mathit{t}}^{\left[\mathbf{2}\right]}$
NSCBINAR(1)NB	0.132	0.141	0.109	0.120
NSUBINAR(1)GEOM	0.120	0.129	0.098	0.104
NSCBINAR(1)GEOM	0.196	0.189	0.155	0.144

displays the RMSEs and mean absolute deviation (MAD) under the NSCBINAR(1)NB, NSCBINAR(1)GEOM and NSUBINAR(1)GEOM.

The standard errors, RMSE and MAD, illustrate that NSUBINAR(1)GEOM yields almost the same measures as NSCBINAR(1)NB, but better estimates than the NSCBINAR(1)GEOM.

6. Conclusions

In this paper, the unconstrained non-stationary BINAR(1) model with geometric marginals was modelled. However, during the development of this model, it is observed that the joint probability function of the innovation series is rather complex to derive and this limits the construction of a conditional likelihood function to estimate the unknown model parameters. Hence, this paper proposes an alternative GQL estimation approach that only requires the correct specification of the score and moment vectors. As for the derivation of the higher-order entries of the auto-covariance matrix, the multivariate normality assumption was used. With regard to the simulation study, it is shown that the GQL approach provides consistent parameter estimates and statistically more efficient estimates. In the analysis of the accident data, reliable estimates of the different covariates were obtained and comparable to NSCBINAR(1)NB and NSCBINAR(1)GEOM. The RMSE and MAD show that NSUBINAR(1)GEOM yields better forecasts than the other two competing models, a similar conclusion illustrated in [7]. This model is commendable to analyse over-dispersed series under non-stationary setups characterized by time-dependent effects.