Next Article in Journal
A New Hybrid Approach for Clustering, Classification, and Prediction of World Development Indicators Combining General Type-2 Fuzzy Systems and Neural Networks
Previous Article in Journal
Certain New Applications of Symmetric q-Calculus for New Subclasses of Multivalent Functions Associated with the Cardioid Domain
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Bivariate Random Coefficient Integer-Valued Autoregressive Model Based on a ρ-Thinning Operator

1
School of Mathematics, Jilin University, Changchun 130012, China
2
School of Mathematics and Statistics, Liaoning University, Shenyang 110136, China
*
Author to whom correspondence should be addressed.
Axioms 2024, 13(6), 367; https://doi.org/10.3390/axioms13060367
Submission received: 17 April 2024 / Revised: 20 May 2024 / Accepted: 28 May 2024 / Published: 29 May 2024
(This article belongs to the Section Mathematical Analysis)

Abstract

:
While overdispersion is a common phenomenon in univariate count time series data, its exploration within bivariate contexts remains limited. To fill this gap, we propose a bivariate integer-valued autoregressive model. The model leverages a modified binomial thinning operator with a dispersion parameter ρ and integrates random coefficients. This approach combines characteristics from both binomial and negative binomial thinning operators, thereby offering a flexible framework capable of generating counting series exhibiting equidispersion, overdispersion, or underdispersion. Notably, our model includes two distinct classes of first-order bivariate geometric integer-valued autoregressive models: one class employs binomial thinning (BVGINAR(1)), and the other adopts negative binomial thinning (BVNGINAR(1)). We establish the stationarity and ergodicity of the model and estimate its parameters using a combination of the Yule–Walker (YW) and conditional maximum likelihood (CML) methods. Furthermore, Monte Carlo simulation experiments are conducted to evaluate the finite sample performances of the proposed estimators across various parameter configurations, and the Anderson-Darling (AD) test is employed to assess the asymptotic normality of the estimators under large sample sizes. Ultimately, we highlight the practical applicability of the examined model by analyzing two real-world datasets on crime counts in New South Wales (NSW) and comparing its performance with other popular overdispersed BINAR(1) models.

1. Introduction

Bivariate count data are prevalent across various scenarios and frequently represent the occurrences of two distinct events, objects, or individuals over a specific time frame. Bivariate integer-valued time series models, in particular, excel at preserving the paired relationship between two correlated count variables observed over certain time intervals. Such data arise in many fields, including guest nights in hotels and cottages [1], tick-by-tick data of two highly traded stocks [2], traffic accidents happening both in daylight and nighttime [3], the occurrence of different offenses in specific areas [4], and the counts of asymptomatic and symptomatic COVID-19 cases within considered regions.
Substantial research interest has been directed toward the analysis of bivariate integer-valued time series, particularly regarding their cross-correlated nature. One direction is to describe the correlation between series by employing different bivariate innovation distributions. For instance, [3] proposed a bivariate diagonal INAR(1) model with bivariate Poisson and bivariate negative binomial innovations. Similarly, [5] found that BINAR(1) models with Poisson–Lindley (PL) innovations outperform other competing INAR(1) models, whether based on diagonal or full coefficient matrices. Furthermore, [6] presented a bivariate full INAR(1) model, assuming a time-dependent innovation vector, where the mean of the innovation vector linearly increases with the previous population size. However, these models are all based on constant coefficients.
To introduce more flexibility into BINAR models, [7] considered inference for a bivariate random coefficient INAR(1) model with different geometric marginals. Thereafter, [8] proposed a more general bivariate diagonal random coefficient INAR(1) process (BRCINAR(1)) with dependent innovations. Moreover, [9] compared the performances of random coefficient BINAR(1) models based on the bivariate negative binomial distributions constructed in different ways with explanatory variables.
While these studies have made significant strides in enhancing BINAR(1) model flexibility through various marginal and innovation distributions, they have primarily focused on binomial thinning operators. The binomial thinning operator, initially proposed by [10], remains the most widely utilized approach in modeling integer-valued time series due to its capability of producing integer-valued results and offering strong interpretability. Denoted as “∘”, this operator is defined as follows:
α X = i = 1 X B i ,
where α ( 0 , 1 ) , X represents a non-negative integer-valued random variable, and { B i } denotes a sequence of independent and identically distributed Bernoulli random variables. Each B i has a probability P ( B i = 1 ) = α and P ( B i = 0 ) = 1 α , independent of X. Essentially, the binomial thinning operator assigns a value of either 0 or 1 to each counting random variable B i , making it suitable for modeling scenarios where random events either survive or vanish after a period of observation.
However, in situations where the observed unit has the potential to generate multiple countable elements or trigger further stochastic occurrences beyond mere survival or disappearance, the Bernoulli random variable may not be the most suitable choice for constructing the counting series. For example, in the context of infectious diseases, an infected individual may not only survive or die but may also contribute to the generation of new cases. To address this limitation, [11] introduced the negative binomial thinning operator, defined as:
α X = i = 1 X G i ,
where α ( 0 , 1 ) , and { G i } represents a sequence of independent and identically distributed geometric random variables with a mean of α , also independent of X.
In fact, while there are many other thinning operators proven useful in univariate integer-valued time series [12,13], the literature on different operators applied to bivariate data is limited. As far as we know, a study by [14] constructed a BINAR(1) model based on the signed thinning operator, which is capable of accommodating data with negative observations. In addition, another study by [15] proposed a new BINAR(1) model that extended the negative binomial thinning operator.
The aim of this paper is to introduce a more sophisticated ρ -binomial thinning operator [12] to bivariate integer-valued time series and explore statistical inference of the proposed model. The motivation behind this endeavor lies in the remarkable versatility exhibited by the counting series derived from this thinning operator. Specifically, this thinning operator enables us to describe equidispersion, overdispersion, or underdispersion characteristics concurrently within both the counting processes. Consequently, the ρ -thinning operator extends the capabilities of the binomial thinning [10] and negative binomial thinning operators [16], offering superior fitting capabilities to paired count data.
The outline of the paper is structured as follows. In Section 2, we introduce a definition of the ρ -BVGINAR(1) model and discuss the basic properties of the thinning operator for bivariate vectors. The properties of the model are further examined in Section 3. Section 4 estimates the model parameters by integrating the Yule–Walker (YW) and maximum likelihood (CML) methods, followed by simulation studies to explore the asymptotic properties of the estimators under various parameter combinations. Moreover, the Anderson-Darling (AD) test is also performed to assess the asymptotic normality of the estimators under large sample sizes. Section 5 illustrates the application of the proposed models to two real-world datasets of crime counts in New South Wales. We examine datasets with varying levels of overdispersion indices and compare the performance of our models with other bivariate integer-valued models. Finally, Section 6 provides concluding remarks and an outlook for future research. Some proofs and figures are provided in the Appendix for reference.

2. Construction of the Model

Ref. [17] introduced a novel variant of the Bernoulli distribution, termed an inflated-parameter Bernoulli (IBe) distribution, designed to model univariate count data exhibiting overdispersion. This distribution incorporates an additional parameter, ρ , allowing for more flexible dispersion indices. The probability mass function (pmf) for the IBe distribution is defined as follows:
Pr ( W = w ) = 1 α , w = 0 , α ( ρ 1 + ρ ) w 1 ( 1 1 + ρ ) , w { 1 , 2 , } ,
where α [ 0 , 1 ] and ρ [ 0 , 1 ) . The distribution can be denoted as W IBe ( α , ρ ) . The mean, variance, and probability generating function (pgf) of the IBe distribution are detailed below:
E ( W ) = α ( 1 + ρ ) , Var ( W ) = α ( 1 + ρ ) [ ( 1 + 2 ρ ) α ( 1 + ρ ) ] , Φ W ( s ) = 1 ( 1 s ) [ α ( 1 + ρ ) ρ ] 1 + ρ ( 1 s ) , | s | < 1 .
Moreover, the dispersion index is expressed as:
I Y : = Var ( W ) E ( W ) = ρ + ( 1 + ρ ) ( 1 α ) .
Interestingly, the IBe distribution presents three dispersion scenarios depending on the values of the parameters:
  • Overdispersion is observed when α 2 α < ρ < 1 .
  • Underdispersion occurs if 0 ρ < α 2 α .
  • Equidispersion is achieved at ρ = α 2 α .
The distribution can also degenerate into two important distributions:
  • The standard Bernoulli distribution, when ρ = 0 and w { 0 , 1 } , with mean α ;
  • The geometric distribution, when α = ρ 1 + ρ , with mean ρ .
Inspired by these properties, [12] formulated a ρ -binomial thinning operator, defined as:
α ρ   X : = i = 1 X W i ,
where X is a non-negative integer-valued random variable, and { W i } i = 1 X is a sequence of inflated-parameter Bernoulli random variables with the pmf given by Equation (1), mutually independent of X. This operator has proven its practical utility in modeling univariate integer-valued time series, known as the ρ -GINAR process.
However, it is crucial to acknowledge that cross-correlations are prevalent in most paired count time series. Hence, this paper aims to enhance the ρ -GINAR process by extending it into a bivariate domain, thereby more effectively capturing the inherent correlations within the data.
We propose the ρ -BVGINAR(1) model, which is a novel bivariate random coefficient INAR(1) process characterized by the following recursive equation:
X t = A t ρ   X t 1 + Z t = U 1 , t U 2 , t V 1 , t V 2 , t   ρ   X 1 , t 1 X 2 , t 1 + Z 1 , t Z 2 , t .
Here,
(i)
A t represents a random coefficient matrix comprising two mutually independent bivariate random vectors, ( U 1 , t , U 2 , t ) and ( V 1 , t , V 2 , t ) , each with independent and identically distributed (i.i.d.) components and with pmf values as:
P ( U 1 , t = α 1 , U 2 , t = 0 ) = 1 P ( U 1 , t = 0 , U 2 , t = α 1 ) = p 1 ,
P ( V 1 , t = α 2 , V 2 , t = 0 ) = 1 P ( V 1 , t = 0 , V 2 , t = α 2 ) = p 2 ,
where α 1 , α 2 ( 0 , 1 ) and p 1 , p 2 [ 0 , 1 ] . The matrix operation A t ρ   replicates matrix multiplication while preserving the properties of random coefficient thinning.
(ii)
The innovation { Z t } t Z is a sequence of i.i.d. bivariate non-negative integer-valued random vectors with mutually independent elements { Z 1 , t } and { Z 2 , t } and independent of X s for s < t .
Remark 1. 
The proposed bivariate INAR(1) process based on the ρ-binomial thinning operator has two sub-models:
  • When ρ = 0 , it corresponds to the bivariate INAR(1) with geometric marginals introduced by [4].
  • When α i = ρ 1 + ρ , for both i = 1 and i = 2 , it aligns with the BVNGINAR(1) proposed by [18].
Next, we further explore the properties of the ρ -thinning operator for vectors.
Lemma 1. 
Consider A t ρ   X t 1 as defined in Equation (4). Then:
(i) 
E ( A t ρ   X t 1 ) = A ˜ E ( X t 1 ) , where
A ˜ = ( 1 + ρ ) A = ( 1 + ρ ) E ( A t ) = ( 1 + ρ ) α 1 p 1 α 1 ( 1 p 1 ) α 2 p 2 α 2 ( 1 p 2 ) .
(ii) 
E ( ( A t ρ   X t 1 ) B ) = A ˜ E ( X t 1 B ) for a random vector B independent of A t .
(iii) 
E ( B ( A t ρ   X t 1 ) ) = E ( B X t 1 ) A ˜ for a random vector B independent of A t .
(iv) 
E ( ( A t ρ   X t 1 ) ( A t ρ   X t 1 ) ) = A ˜ E ( X t 1 X t 1 ) A ˜ + C , where C has elements
c 12 = 0 , c 21 = 0 ,
c 11 = ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) E ( X 1 , t X 2 , t ) 2 + α 1 ( 1 + ρ ) ( 1 + 2 ρ α 1 ( 1 + p 1 ) ) ( p 1 E ( X 1 , t ) + ( 1 p 1 ) E ( X 2 , t ) ) ) ,
c 22 = ( 1 + ρ ) 2 α 2 2 p 2 ( 1 p 2 ) E ( X 1 , t X 2 , t ) 2 + α 2 ( 1 + ρ ) ( 1 + 2 ρ α 2 ( 1 + p 2 ) ) ( p 2 E ( X 1 , t ) + ( 1 p 2 ) E ( X 2 , t ) ) ) .
Lemma 2. 
If we assume α 1 ( 1 + ρ ) , α 2 ( 1 + ρ ) ( 0 , 1 ) , then all eigenvalues of matrix A ˜ lie within the unit circle.
Proof. 
Similar to [19], we outline the key steps. Let λ 1 and λ 2 denote the eigenvalues, with λ 2 > λ 1 assumed without loss of generality.
(i)
We have λ 1 + λ 2 = ( 1 + ρ ) α 1 p 1 + ( 1 + ρ ) α 2 ( 1 p 2 ) > 0 , indicating λ 2 > 0 .
(ii)
Furthermore, we calculate λ 1 λ 2 = ( 1 + ρ ) 2 α 1 p 1 α 2 ( 1 q 1 ) ( 1 + ρ ) 2 α 1 ( 1 p 1 ) α 2 p 2 ; then, we have ( 1 λ 1 ) ( 1 λ 2 ) = 1 ( λ 1 + λ 2 ) + λ 1 λ 2 = ( 1 ( 1 + ρ ) α 1 p 1 ) ( 1 ( 1 + ρ ) α 2 ( 1 p 2 ) ) ( 1 + ρ ) 2 α 1 ( 1 p 1 ) α 2 p 2 ( 1 p 1 ) ( 1 ( 1 p 2 ) ) ( 1 p 1 ) p 2 = 0 . Hence, it necessitates that either both λ 1 , λ 2 < 1 or both λ 1 , λ 2 > 1 . Since λ 1 + λ 2 = ( 1 + ρ ) α 1 p 1 + ( 1 + ρ ) α 2 ( 1 p 2 ) < 2 , we deduce both λ 1 , λ 2 < 1 .
(iii)
Similarly, evaluating ( 1 + λ 1 ) ( 1 + λ 2 ) = 1 + ( λ 1 + λ 2 ) + λ 1 λ 2 = ( 1 + ( 1 + ρ ) α 1 p 1 ) ( 1 + ( 1 + ρ ) α 2 ( 1 p 2 ) ) ( 1 + ρ ) 2 α 1 ( 1 p 1 ) α 2 p 2 > 0 . As λ 2 > 0 , then λ 1 < 1 .
Consequently, we conclude 1 < λ 1 < 1 and 0 < λ 2 < 1 under the conditions α 1 ( 1 + ρ ) , α 2 ( 1 + ρ ) ( 0 , 1 ) . Moreover, this condition ensures the stationarity of both the first- and second-order moments of the process. □
Proposition 1. 
If α 1 ( 1 + ρ ) , α 2 ( 1 + ρ ) ( 0 , 1 ) and p 1 , p 2 [ 0 , 1 ] , a strictly stationary bivariate integer-valued time series { X t } t Z satisfying Equation (4) exists uniquely. Moreover, the process is ergodic.
The proof of Proposition 1 is provided in Appendix B. Now let us derive the moments and conditional moments of the ρ -BVGINAR(1) process.
Proposition 2. 
Suppose the bivariate time series { X t } { t Z } is a stationary process defined by Equation (4); then, for t 1 , i = 1 , 2 , we have
(i) 
E ( X t ) = ( I A ˜ ) 1 E ( Z ) .
(ii) 
Σ X t + 1 = A ˜ Σ X t A ˜ + C + Σ Z , where C has elements
c 12 = 0 , c 21 = 0 ,
c 11 = ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) E ( X 1 , t X 2 , t ) 2 + α 1 ( 1 + ρ ) ( 1 + 2 ρ α 1 ( 1 + p 1 ) ) ( p 1 E ( X 1 , t ) + ( 1 p 1 ) E ( X 2 , t ) ) ) ,
c 22 = ( 1 + ρ ) 2 α 2 2 p 2 ( 1 p 2 ) E ( X 1 , t X 2 , t ) 2 + α 2 ( 1 + ρ ) ( 1 + 2 ρ α 2 ( 1 + p 2 ) ) ( p 2 E ( X 1 , t ) + ( 1 p 2 ) E ( X 2 , t ) ) ) .
(iii) 
E ( X i , t + 1 | X 1 , t , X 2 , t ) = ( 1 + ρ ) α i [ p i X 1 , t + ( 1 p i ) X 2 , t ] + μ Z i , t .
(iv) 
Var ( X i , t + 1 | X 1 , t , X 2 , t ) = ( 1 + ρ ) 2 α i 2 p i ( 1 p i ) ( X 1 , t 2 + X 2 , t 2 ) + ( 1 + ρ ) [ 1 + 2 ρ ( 1 + ρ ) α i ] α i [ p i X 1 , t + ( 1 p i ) X 2 , t ] 2 ( 1 + ρ ) 2 α i 2 p i ( 1 p i ) X 1 , t X 2 , t + σ Z i , t 2 .
where μ Z i , t and σ Z i , t 2 are the mean and variance, respectively, of { Z i , t } , for i = 1 , 2 .
We define a stationary bivariate time series ( X 1 , t , X 2 , t ) t Z according to Equation (4). By specifying appropriate marginal distributions for X 1 , t and X 2 , t , we can deduce the respective marginal distributions of the innovation Z 1 , t and Z 2 , t . This process is clarified by the theorem below.
Theorem 1. 
Let α 1 ( 1 + ρ ) , α 2 ( 1 + ρ ) ( 0 , 1 ) , p 1 , p 2 [ 0 , 1 ] if { X 1 , t } = d { X 2 , t } = d G e o m ( μ 1 + μ ) , and μ > 0 ; then the distributions of the innovation processes { Z 1 , t } and { Z 2 , t } are as follows:
Z 1 , t = d G e o m ( ρ 1 + ρ ) , w . p . α 1 μ ( 1 + ρ ) μ ρ , G e o m ( μ 1 + μ ) , w . p . 1 α 1 μ ( 1 + ρ ) μ ρ ,
Z 2 , t = d G e o m ( ρ 1 + ρ ) , w . p . α 2 μ ( 1 + ρ ) μ ρ , G e o m ( μ 1 + μ ) , w . p . 1 α 2 μ ( 1 + ρ ) μ ρ .
Proof. 
For intuitive purposes, the bivariate time series model { ( X 1 , t , X 2 , t ) } t Z can be represented as
X 1 , t = α 1 ρ   X 1 , t 1 + Z 1 , t , w . p . p 1 , α 1 ρ   X 2 , t 1 + Z 1 , t , w . p . 1 p 1 ,
X 2 , t = α 2 ρ   X 1 , t 1 + Z 2 , t , w . p . p 2 , α 2 ρ   X 2 , t 1 + Z 2 , t , w . p . 1 p 2 .
Given { X 1 , t } = d { X 2 , t } = d G e o m ( μ 1 + μ ) and leveraging the properties of the thinning operator and the stationarity of the process, we have
Φ α 1 ρ   X 1 , t ( s ) = E ( s α 1 ρ   X 1 , t ) = E [ E ( s j = 1 X 1 , t W j ( α 1 , ρ ) | X 1 , t ) ] = E [ E ( s W j ( α 1 , ρ ) | X 1 , t ) X 1 , t ] = E [ ( Φ W ( s ) ) X 1 , t ] = Φ X 1 , t ( Φ W ( s ) ) = Φ X 1 , t 1 ( 1 s ) [ α 1 ( 1 + ρ ) ρ ] 1 + ρ ( 1 s ) .
Then the pgf of Z 1 , t can be obtained by
Φ Z 1 , t ( s ) = Φ X 1 , t ( s ) p 1 Φ X 1 , t ( Φ W ( s ) ) + ( 1 p 1 ) Φ X 1 , t ( Φ W ( s ) ) = 1 1 + μ μ s p 1 · Φ X 1 , t ( 1 ( 1 s ) [ α 1 ( 1 + ρ ) ρ ] 1 + ρ ( 1 s ) ) + ( 1 p 1 ) · Φ X 2 , t ( 1 ( 1 s ) [ α 1 ( 1 + ρ ) ρ ] 1 + ρ ( 1 s ) ) = 1 + ρ ( 1 s ) + α 1 μ ( 1 + ρ ) ( 1 s ) [ 1 + μ ( 1 s ) ] [ 1 + ρ ( 1 s ) ] = [ 1 + ρ ( 1 s ) + α 1 μ ( 1 + ρ ) ( 1 s ) ] ( μ ρ ) [ 1 + μ ( 1 s ) ] [ 1 + ρ ( 1 s ) ] ( μ ρ ) = ( μ ρ ) α 1 μ ( 1 + ρ ) ( μ ρ ) [ 1 + μ ( 1 s ) ] + α 1 μ ( 1 + ρ ) ( μ ρ ) [ 1 + ρ ( 1 s ) ] = 1 α 1 μ ( 1 + ρ ) μ ρ 1 1 + μ ( 1 s ) + α 1 μ ( 1 + ρ ) μ ρ 1 1 + ρ ( 1 s ) .
The innovation Z 1 , t clearly consists of a combination of two geometrically distributed random variables. Similar derivation holds for Z 2 , t . Thus, the distributions of innovation process can be expressed as Equations (5) and (6). Notably, it is emphasized that 0 < ρ < min ( μ ( 1 α 1 ) 1 + α 1 μ , μ ( 1 α 2 ) 1 + α 2 μ ) is necessary to ensure the non-negativity of all probabilities for Z 1 , t and Z 2 , t . □

3. Properties

Lemma 3. 
Let 0 < ρ < min ( μ ( 1 α 1 ) 1 + α 1 μ , μ ( 1 α 2 ) 1 + α 2 μ ) , and p 1 , p 2 [ 0 , 1 ] . The correlation coefficient between { X 1 , t } and { X 2 , t } lies in [ 0 , 1 ) and is expressed as:
γ = ( 1 + ρ ) 2 α 1 α 2 ( p 1 p 2 + ( 1 p 1 ) ( 1 p 2 ) ) 1 ( 1 + ρ ) 2 α 1 α 2 ( p 1 ( 1 p 2 ) + ( 1 p 1 ) p 2 ) .
It is evident that the correlation coefficient γ lies in the interval [ 0 , 1 ) .
From Lemma 1, we find that the covariance matrix between the random vectors X t and X 0 is given by
Cov ( X t , X 0 ) = A ˜ Cov ( X t 1 , X 0 ) = = A ˜ t Var ( X 0 ) 0 , t .
The one-step-ahead conditional expectation can be derived as follows:
E ( X t + 1 | X t ) = A ˜ X t + E ( Z t + 1 ) .
Typically, the k-step ahead conditional expectation of X t is
E ( X t | X t k ) = E ( A t ρ   X t 1 + Z t | X t k ) = E ( A t k ρ   X t k + j = 1 k 1 A t j ρ   Z t j | X t k ) = A ˜ k X t k + j = 1 k 1 [ A ˜ ] j E ( Z t ) = A ˜ k X t k + [ I A ˜ k ] [ I A ˜ ] 1 E ( Z t ) .
Meanwhile, we observe that E ( Z t ) = [ I A ˜ ] E ( X t ) , thereby implying that
E ( X t | X t k ) = A ˜ k X t k + [ I A ˜ k ] E ( X t ) E ( X t ) , k .
This finding validates the characteristic of autoregressive processes, whereby the conditional expectation converges to the unconditional expectation as the number of steps approaches infinity.
Due to the conditional independence of the random variables X 1 , t and X 2 , t conditioned on X 1 , t 1 and X 2 , t 1 , respectively, the conditional probability function can be represented as the product of individual conditional probabilities. Therefore, we can derive the conditional probability function of the random vector ( X 1 , t , X 2 , t ) as follows:
P ( X 1 , t = x , X 2 , t = y | X 1 , t 1 = u , X 2 , t 1 = v ) = P ( X 1 , t = x | X 1 , t 1 = u , X 2 , t 1 = v ) · P ( X 2 , t = y | X 1 , t 1 = u , X 2 , t 1 = v ) ,
where
P ( X 1 , t = x | X 1 , t 1 = u , X 2 , t 1 = v ) = p 1 ( P ( α 1 ρ   X 1 , t 1 = 0 | X 1 , t 1 = u ) P ( Z 1 , t = x ) + k = 1 x P ( α 1 ρ   X 1 , t 1 = k | X 1 , t 1 = u ) P ( Z 1 , t = x k ) ) + ( 1 p 1 ) ( P ( α 1 ρ   X 2 , t 1 = 0 | X 2 , t 1 = v ) P ( Z 1 , t = x ) + k = 1 x P ( α 1 ρ   X 2 , t 1 = k | X 2 , t 1 = v ) P ( Z 1 , t = x k ) ) ,
and
P ( X 2 , t = y | X 1 , t 1 = u , X 2 , t 1 = v ) = p 2 ( P ( α 2 ρ   X 2 , t 1 = 0 | X 1 , t 1 = u ) P ( Z 2 , t = y ) + k = 1 y P ( α 2 ρ   X 2 , t 1 = k | X 1 , t 1 = u ) P ( Z 2 , t = y k ) ) + ( 1 p 2 ) ( P ( α 2 ρ   X 2 , t 1 = 0 | X 1 , t 1 = v ) P ( Z 2 , t = y ) + k = 1 y P ( α 2 ρ   X 2 , t 1 = k | X 2 , t 1 = v ) P ( Z 2 , t = y k ) ) .
Here, the distributions of the random variables Z 1 , t and Z 2 , t are defined by Theorem 1, so their probability mass functions are given by
P ( Z 1 , t = z 1 ) = μ z 1 ( 1 + μ ) z 1 + 1 ( 1 α 1 μ ( 1 + ρ ) μ ρ ) + ρ z 1 ( 1 + ρ ) z 1 + 1 ( α 1 μ ( 1 + ρ ) μ ρ )
and
P ( Z 2 , t = z 2 ) = μ z 2 ( 1 + μ ) z 2 + 1 ( 1 α 2 μ ( 1 + ρ ) μ ρ ) + ρ z 2 ( 1 + ρ ) z 2 + 1 ( α 2 μ ( 1 + ρ ) μ ρ ) .

4. Estimation Procedure

In this section, we consider X t as a strictly stationary and ergodic solution of the ρ -BVGINAR(1) process, with { X t } t Z representing a series of observations generated from this process. We discuss the estimation of the model parameters, comprising six parameters: one for thinning the distribution ( ρ ), two for the autocorrelation coefficients ( α 1 , α 2 ), two for specifying the dependence between processes X 1 , t and X 2 , t ( p 1 , p 2 ), and one for the marginal distributions ( μ ). Considering the unique characteristics of these parameters, we integrate two estimation approaches: the Yule–Walker (YW) and the conditional maximum likelihood (CML) methods.
The sample mean is commonly employed for estimating model parameters in time series analysis. Since the model assumption is that the marginal distribution { X i , t } = d G e o m ( μ 1 + μ ) , then μ = E ( X 1 , t ) = E ( X 2 , t ) . Thus, the reasonable estimate would be:
μ ^ Y W = 1 2 n j = 1 n ( X 1 , j + X 2 , j ) .
Theorem 2. 
The estimator μ ^ Y W is strongly consistent.
Proof. 
Proposition 1 proved that process { X t } is stationary and ergodic. Then, processes { X 1 , t } and { X 2 , t } are jointly stationary, which implies that 1 2 ( X 1 , t + X 2 , t ) is also stationary and ergodic. We have
1 2 n j = 1 n ( X 1 , j + X 2 , j ) a . s . 1 2 E ( X 1 , t + X 2 , t ) .
Therefore, μ ^ Y W a . s . μ . The proof of Theorem 2 is complete. □
Theorem 3. 
The estimator μ ^ Y W is asymptotically normally distributed with parameters ( μ , t 1 v s ) , where v s = k = C o v ( S t , S t k ) and S t = 1 2 ( X 1 , t + X 2 , t ) .
In addition, the conditional maximum likelihood (CML) stands out as one of the most commonly employed techniques for parameter estimation. The CML estimator of parameter vector θ = ( α 1 , α 2 , p 1 , p 2 , ρ ) is the value θ ^ = ( α ^ 1 , α ^ 2 , p ^ 1 , p ^ 2 , ρ ^ ) that maximizes the conditional log-likelihood function L ( θ ) . Suppose that X 1 is fixed. The conditional log-likelihood function is given by
L ( θ ) = t = 2 T log { P ( X 1 , t = x , X 2 , t = y | X 1 , t 1 = u , X 2 , t 1 = v , θ ) } .
Under the given conditions, the conditional probability mass functions of processes X 1 , t and X 2 , t can be represented as the products of their respective conditional probabilities. These probabilities result from the convolution of the ρ -binomial distribution and the probability mass function of the corresponding innovation processes. Specifically,
L ( v ) = t = 2 n log [ P ( X 1 , t = x | X 1 , t 1 = u , X 2 , t 1 = v ) ] + t = 2 n log [ P ( X 2 , t = y | X 1 , t 1 = u , X 2 , t 1 = v ) ] = t = 2 n log [ p 1 ( 1 α 1 ) u + ( 1 p 1 ) ( 1 α 1 ) v ( μ x ( 1 + μ ) x + 1 ( 1 α 1 μ ( 1 + ρ ) μ ρ ) + ρ x ( 1 + ρ ) x + 1 × α 1 μ ( 1 + ρ ) μ ρ ) + ( p 1 k = 1 x ( 1 α 1 ) u ( ρ 1 + ρ ) k i = 1 min ( k , u ) ( u i ) ( k 1 i 1 ) ( α 1 ρ ( 1 α 1 ) ) i + ( 1 p 1 ) × k = 1 x ( 1 α 1 ) v ( ρ 1 + ρ ) k i = 1 min ( k , v ) ( v i ) ( k 1 i 1 ) ( α 1 ρ ( 1 α 1 ) ) i ) ( μ x k ( 1 + μ ) x + 1 k 1 α 1 μ ( 1 + ρ ) μ ρ ) + ρ x k ( 1 + ρ ) x + 1 k ( α 1 μ ( 1 + ρ ) μ ρ ) ) ] + t = 2 n log [ p 2 ( 1 α 2 ) u + ( 1 p 2 ) ( 1 α 2 ) v × μ y ( 1 + μ ) y + 1 ( 1 α 2 μ ( 1 + ρ ) μ ρ ) + ρ y ( 1 + ρ ) y + 1 α 2 μ ( 1 + ρ ) μ ρ + ( p 2 k = 1 y ( 1 α 2 ) u ( ρ 1 + ρ ) k × i = 1 min ( k , u ) ( u i ) ( k 1 i 1 ) ( α 2 ρ ( 1 α 2 ) ) i + ( 1 p 2 ) k = 1 y ( 1 α 2 ) v ( ρ 1 + ρ ) k i = 1 min ( k , v ) ( v i ) ( k 1 i 1 ) × ( α 2 ρ ( 1 α 2 ) ) i ) μ y k ( 1 + μ ) y + 1 k ( 1 α 2 μ ( 1 + ρ ) μ ρ ) + ρ y k ( 1 + ρ ) y + 1 k ( α 2 μ ( 1 + ρ ) μ ρ ) ] .
The likelihood function reveals that the parameters α 1 , α 2 , ρ , and μ often interact multiplicatively, posing significant challenges for the optimization process. To address potential issues with parameter identifiability and leverage the specific characteristics of these parameters, we have implemented a stepwise optimization strategy. Initially, we estimate μ using the Yule–Walker method. This estimate μ ^ Y W is then incorporated back into the likelihood function to facilitate the optimization of the remaining parameters. This tailored approach effectively mitigates the risk of converging to local optima: a prevalent concern with non-convex objective functions.
For numerical maximization, we employ the “nlm” function from R programming software. All computational experiments were performed using R version 4.0.3 on a system equipped with an Intel Xeon Gold 6154 processor (Intel Corporation, Santa Clara, CA, USA) and 256 GB of RAM.
Next, we present the detailed simulation study design and results. We generated ρ -BVGINAR(1) samples with various model parameterizations and sample sizes n = 150 , 300 , 600 , 1200 , 3000 , where n = 300 is close to the length of the crime counts that will be analyzed in Section 5. We considered the following parameter configurations:
  • Model (A): ( α 1 , α 2 , p 1 , p 2 , ρ , μ ) = (0.3, 0.25, 0.2, 0.15, 0.1, 5);
  • Model (B): ( α 1 , α 2 , p 1 , p 2 , ρ , μ ) = (0.3, 0.25, 0.2, 0.15, 0.3, 5);
  • Model (C): ( α 1 , α 2 , p 1 , p 2 , ρ , μ ) = (0.4, 0.4, 0.2, 0.15, 0.25, 5);
  • Model (D): ( α 1 , α 2 , p 1 , p 2 , ρ , μ ) = (0.6, 0.4, 0.3, 0.7, 0.3, 3);
  • Model (E): ( α 1 , α 2 , p 1 , p 2 , ρ , μ ) = (0.6, 0.4, 0.7, 0.3, 0.3, 3).
Recalling the properties of IBe random variables discussed in Section 2, we explored diverse parameter combinations of α 1 , α 2 , and ρ in our simulation study. Models (A), (B), and (C) represent scenarios of underdispersion ( ρ < α 1 2 α 1 , α 2 2 α 2 ), overdispersion ( ρ > α 1 2 α 1 , α 2 2 α 2 ), and equidispersion ( ρ = α 1 2 α 1 , α 2 2 α 2 ) of the IBe random variables W i ( α 1 , ρ ) and W i ( α 2 , ρ ) under ρ -thinning. Conversely, Models (D) and (E) are characterized by distinct dispersion patterns of W i ( α 1 , ρ ) and W i ( α 2 , ρ ) . Specifically, Model (D) sets underdispersion of W i ( α 1 , ρ ) and overdispersion of W i ( α 2 , ρ ) ( ρ < α 1 2 α 1 , ρ > α 2 2 α 2 ), while Model (E) shows overdispersion of W i ( α 1 , ρ ) and underdispersion of W i ( α 2 , ρ ) ( ρ > α 1 2 α 1 , ρ < α 2 2 α 2 ).
To assess model performance, we employed two widely recognized criteria: mean absolute error (MAE) and root mean squared error (RMSE), based on M = 500 replications for each model parametrization. MAE is preferred for its robustness against outliers, while RMSE provides a more detailed measure of errors and is particularly sensitive to larger deviations. The use of both MAE and RMSE allows for a comprehensive evaluation of estimation accuracy. They are defined as follows:
MAE = 1 M m = 1 M | θ ^ m θ m | , RMSE = 1 M m = 1 M ( θ ^ m θ m ) 2 ,
where θ ^ m is the estimate of θ at the m-th replication.
In addition, to demonstrate the asymptotic normality of the estimators, we conducted a goodness-of-fit test for normality. The Anderson-Darling (AD) test, proposed by [20], is a statistical test used to assess whether data come from a specific distribution: typically, a normal distribution. Unlike other normality tests, the AD test gives more weight to the tails of the distribution, making it more powerful for detecting non-normality. Therefore, we selected the AD test to examine the asymptotic normality of the estimators by using the procedure from the R package “nortest” authored by [21].
Table 1, Table 2 and Table 3 report the estimates, biases, MAEs, and RMSEs for Models (A)–(E) across various sample sizes, as well as the AD test statistics (denoted AD in the tables) and corresponding p-values for n = 3000 . From Table 1, we observe that the biases, MAEs, and RMSEs of the estimates for Models (A) and (B) decrease as the sample size n increases, as expected. Figure 1 and Figure 2 also illustrate the notable downward trends in MAEs and RMSEs of the estimates, implying the consistency of the proposed estimators with increasing values of n. A similar conclusion can be drawn from Table 2 and Table 3, along with the corresponding visual curves in Figure 3, Figure 4 and Figure 5. Based on the above discussions, we conclude that the estimation method can produce reliable parameter estimators.
Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15 display the Gaussian QQ plots of the proposed estimators for Models (A)–(E) across various sample sizes. For small sample sizes, particularly when n = 150 , the data points are concentrated along a non-45-degree diagonal. However, as n increases, more data points align closely with the 45-degree line, indicating a good match between the estimator and normal distributions. Furthermore, all p-values of the AD normality test for the estimates when n = 3000 for Models (A)–(E) are greater than the significance level of 0.05, as shown in Table 1, Table 2 and Table 3, further confirming the asymptotic normality of the estimators for large sample sizes. Based on the above facts, we conclude that the proposed estimation method is trustworthy for the models under consideration and can yield estimators with asymptotic normality.

5. Real Data Examples

In this section, we present two applications to demonstrate the effectiveness of our proposed ρ -BVGINAR(1) model in capturing overdispersion and zero inflation phenomena. We compare our model with existing bivariate INAR(1) models utilizing standard binomial and negative binomial thinning operators. Specifically, we consider the following models:
  • BVNGINAR(1) model ([18]);
  • BVPOINAR(1) model ([4]);
  • BVMIXINAR(1) model ([22]).
Note that both the binomial and negative binomial thinning operators are special cases of the ρ -binomial thinning operator: the three models above are considered. For fairness, we estimate the mean parameter μ using the Yule–Walker method across all models, while employing maximum likelihood estimation for other parameters. Model performance was evaluated using the Akaike information criterion (AIC) and Bayesian information criterion (BIC).
The dataset used in this analysis was obtained from the NSW Bureau of Crime Statistics and Research (BOCSAR). It is organized by local government area (LGA), offense category (including subcategories), and month. Covering the period from January 1995 to December 2023, the dataset includes 348 monthly crime counts in New South Wales (NSW), Australia. The data can be downloaded from the following website: https://www.bocsar.nsw.gov.au/Pages/bocsar_datasets/Offence.aspx, accessed on 13 May 2024.
To illustrate the characteristics of the observations, we present descriptive statistics for each series. These statistics include the minimum (Min), maximum (Max), median, mean, variance (Var), dispersion index ( I d ), and zero inflation index ( z 0 ). The I d is defined in Equation (2). The z 0 index, introduced by [23], is employed to assess the excess occurrence of zeros in count data and is formulated as:
z 0 = 1 + ln ( p 0 ) λ ,
where p 0 denotes the proportion of zeros, and λ represents the mean. A z 0 value greater than 0 indicates zero inflation, while a value less than 0 suggests zero deflation.

5.1. Crime Data: Disorderly Conduct Counts in Carrathool

In the first application, we analyze disorderly conduct counts in Carrathool, including three subcategories: “offensive conduct” (OCND), “offensive language” (OLNG), and “criminal intent”. Evidently, OCND and OLNG often co-occur due to similar contexts, likely indicating a significant degree of mutual association between their counts. Therefore, we applied the ρ -BVGINAR(1) model to fit the counts of OCND and OLNG.
The time series, autocorrelation function (ACF), and cross-correlation function (CCF) plots for the OCND and OLNG series are depicted in Figure 16 and Figure 17. The ACF plots show autocorrelation in both series, while the values in the CCF plot surpass the confidence interval, suggesting non-independence between the two series. Table 4 presents descriptive statistics for both series. The empirical mean values for OCND and OLNG are relatively close, at 0.3448 and 0.2960, respectively. Both sequences exhibit dispersion indices I d of 1.6098 and 1.3487 , slightly exceeding 1, indicating marginal overdispersion. Moreover, the z 0 values for both series exceed 0, indicating zero inflation characteristics in the data. Further insight is provided by their histograms in Figure 18, which highlight a notable proportion of zeros in each series.
The fitted results of the proposed ρ -BVGINAR(1), BVNGINAR(1), BVPOINAR(1), and BVMIXINAR(1) models are summarized in Table 5. Despite incorporating a mixture of binomial and negative binomial thinning, the BVMIXINAR(1) model exhibits the poorest performance, evidenced by its minimum log-likelihood value ( 526.14 ), maximum AIC (1062.29) and BIC (1081.55) values. The BVNGINAR(1) and BVPOINAR(1) models perform comparably based on their log-likelihood values, appearing only suboptimal compared to the ρ -BVGINAR(1) model. However, the ρ -BVGINAR(1) model achieves the lowest AIC (991.27) and BIC (1008.53) values, indicating superior data fitting.

5.2. Crime Data: Theft Counts in Narrandera

As suggested by a referee, to demonstrate the flexibility and applicability of our model, we selected a real-world example with higher levels of overdispersion for analysis. We focused on the theft counts in Narrandera, encompassing five subcategories: “break and enter dwelling”, “break and enter non-dwelling”, “receiving or handling stolen goods”, “motor vehicle theft”, and “steal from motor vehicle”. Notably, “break and enter dwelling” and “break and enter non-dwelling” exhibit a correlation due to their similar modus operandi and motivations, likely originating from the actions of the same group of offenders. Therefore, we chose to examine the counts of “break and enter thefts into dwelling” (BETD) and “break and enter theft into non-dwelling” (BETND) for further investigation.
Figure 19 and Figure 20 display the time series, ACF, and CCF plots, revealing significant autocorrelation within each series and cross-correlation between the BETD and BETND series. Table 6 presents the dispersion indices I d values of 4.6408 and 3.9099 , markedly exceeding 1, indicating a higher degree of overdispersion compared to the I d values of 1.6098 and 1.3487 observed in the OCND and OLNG counts. Moreover, the z 0 values for both series indicate notable zero inflation, with values of 0.7587 and 0.7681 , respectively. Their histograms are depicted in Figure 21.
Table 7 presents the fitted results. We observe that the BVMIXINAR(1) model exhibits the highest AIC and BIC values, indicating that the model fails to capture the overdispersion and zero inflation characteristic of the dataset. We also notice that the BIC value of the fitted BVPOINAR(1) models is also large, indicating that the model is unsuitable to fit this dataset. While the BVNGINAR(1) model yields satisfactory results, the ρ -BVGINAR(1) model outperforms it in terms of AIC and BIC values, indicating its greater suitability for fitting the BETD and BETND counts.
In conclusion, the ρ -BVGINAR(1) model outperforms the other three models and more effectively captures the overdispersion and zero inflation features in count time series data. It demonstrates superiority in model fitting by striking a balance between flexibility and complexity.

6. Conclusions

This paper introduces a more flexible ρ -BVGINAR(1) model tailored to analyzing bivariate integer-valued time series data with overdispersion characteristic. It extends the ρ -GINAR(1) model [12] to the two-dimensional case. Meanwhile, it is also a generalization of the BVGNAR(1) model or the BVNGINAR(1) model [18], offering enhanced capability for handling excess zeros and overdispersed data. Furthermore, the paper derives the innovation structure of the proposed model, discusses its essential properties, and describes the methodologies for YW and CML estimation. A comprehensive simulation study is conducted to evaluate the finite sample performances of the estimators and their asymptotic properties under various parameter combinations. Two real applications showcase the effectiveness of the proposed model relative to existing ones, demonstrating its utility in practical settings.
Moving forward, there are several promising avenues for future research in the field of bivariate INAR-type models. One promising direction involves exploring the application of the zero-modified geometric distribution as the marginal distribution. This distribution offers the capability to effectively model features such as zero inflation, zero deflation, overdispersion, and underdispersion, as discussed in detail in [24]. In addition, there is potential for further investigation into the modification of various marginal parameters and thinning parameters. Previous works, such as those by [7,19], have demonstrated the effectiveness of such modifications for analyzing bivariate time series data. These approaches hold promise for enhancing the flexibility and applicability of bivariate models and warrant thorough exploration in future research projects.

Author Contributions

Conceptualization and methodology, C.L.; validation and review, D.W. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Social Science Planning Foundation of Liaoning Province (No. L22ZD065) and the National Natural Science Foundation of China (Nos. 12271231, 12001229, 11901053).

Data Availability Statement

Publicly available data sets were analyzed in this study. These data can be found here: https://www.bocsar.nsw.gov.au/Pages/bocsar_datasets/Offence.aspx (accessed on 13 May 2024).

Conflicts of Interest

The authors declare no conflicts of interest.

Appendix A. Proof in Lemma 1

(i)
E ( A t ρ   X t 1 ) = A ˜ E ( X t 1 ) , where A ˜ = ( 1 + ρ ) A = ( 1 + ρ ) E ( A t ) = ( 1 + ρ ) α 1 p 1 α 1 ( 1 p 1 ) α 2 p 2 α 2 ( 1 p 2 ) .
Proof. 
According to the definition of the model, we have
E ( A t ρ   X t 1 ) = E U 1 , t ρ   X 1 , t 1 + U 2 , t ρ   X 2 , t 1 V 1 , t ρ   X 1 , t 1 + V 2 , t ρ   X 2 , t 1 .
We decompose the computation of each element of the matrix.
E ( U 1 , t ρ   X 1 , t 1 + U 2 , t ρ   X 2 , t 1 ) = E { E [ ( U 1 , t ρ   X 1 , t 1 + U 2 , t ρ   X 2 , t 1 ) | ( X 1 , t 1 , X 2 , t 1 ) ] } = E E [ ( i = 1 X 1 , t 1 W i ( U 1 , t , ρ ) + j = 1 X 2 , t 1 W j ( U 2 , t , ρ ) | ( X 1 , t 1 , X 2 , t 1 ) ] = E E [ i = 1 X 1 , t 1 W i ( U 1 , t , ρ ) | X 1 , t 1 ] + E [ j = 1 X 2 , t 1 W j ( U 2 , t , ρ ) | X 2 , t 1 ] = E { X 1 , t 1 · E [ W i ( U 1 , t , ρ ) ] + X 2 , t 1 · E [ W j ( U 2 , t , ρ ) ] } = E [ W i ( U 1 , t , ρ ) ] · E ( X 1 , t 1 ) + E [ W j ( U 2 , t , ρ ) ] · E ( X 2 , t 1 ) = ( 1 + ρ ) · E ( U 1 , t ) · E ( X 1 , t 1 ) + ( 1 + ρ ) · E ( U 2 , t ) · E ( X 2 , t 1 ) = ( 1 + ρ ) [ E ( U 1 , t ) · E ( X 1 , t 1 ) + E ( U 2 , t ) · E ( X 2 , t 1 ) ] ,
where
E [ W i ( U 1 , t , ρ ) ] = E { E [ W i ( U 1 , t , ρ ) | U 1 , t ] } = E [ U 1 , t ( 1 + ρ ) ] = ( 1 + ρ ) E ( U 1 , t ) = ( 1 + ρ ) α 1 p 1 . E [ W j ( U 2 , t , ρ ) ] = E { E [ W j ( U 2 , t , ρ ) | U 2 , t ] } = ( 1 + ρ ) E ( U 2 , t ) = ( 1 + ρ ) α 1 ( 1 p 1 ) .
By the same token,
E ( V 1 , t ρ   X 1 , t 1 + V 2 , t ρ   X 2 , t 1 ) = ( 1 + ρ ) [ E ( V 1 , t ) · E ( X 1 , t 1 ) + E ( V 2 , t ) · E ( X 2 , t 1 ) ] .
Therefore,
E ( A t ρ   X t 1 ) = ( 1 + ρ ) [ E ( U 1 , t ) · E ( X 1 , t 1 ) + E ( U 2 , t ) · E ( X 2 , t 1 ) ] ( 1 + ρ ) [ E ( V 1 , t ) · E ( X 1 , t 1 ) + E ( V 2 , t ) · E ( X 2 , t 1 ) ] = ( 1 + ρ ) · E U 1 , t U 2 , t V 1 , t V 2 , t · E X 1 , t 1 X 2 , t 1 = ( 1 + ρ ) · E ( A t ) · E ( X t 1 ) = A ˜ · E ( X t 1 ) .
(ii)
E [ ( A t ρ   X t 1 ) B ] = A ˜ E ( X t 1 B ) for a random vector B independent of A t .
Proof. 
From the definition of the model, it follows that
E [ ( A t ρ   X t 1 ) B ] = E U 1 , t ρ   X 1 , t 1 + U 2 , t ρ   X 2 , t 1 V 1 , t ρ   X 1 , t 1 + V 2 , t ρ   X 2 , t 1 · ( B 1 , B 2 ) = E ( U 1 , t ρ   X 1 , t 1 + U 2 , t ρ   X 2 , t 1 ) · B 1 ( U 1 , t ρ   X 1 , t 1 + U 2 , t ρ   X 2 , t 1 ) · B 2 ( V 1 , t ρ   X 1 , t 1 + V 2 , t ρ   X 2 , t 1 ) · B 1 ( V 1 , t ρ   X 1 , t 1 + V 2 , t ρ   X 2 , t 1 ) · B 2 .
Decomposing the computation of each element of the matrix, we have
E [ ( U 1 , t ρ   X 1 , t 1 + U 2 , t ρ   X 2 , t 1 ) · B 1 ] = E { E [ ( U 1 , t ρ   X 1 , t 1 ) · B 1 + ( U 2 , t ρ   X 2 , t 1 ) · B 1 | ( X 1 , t 1 , X 2 , t 1 ) ] } = E E [ i = 1 X 1 , t 1 W i ( U 1 , t , ρ ) · B 1 | X 1 , t 1 ] + E [ j = 1 X 2 , t 1 W j ( U 2 , t , ρ ) · B 1 | X 2 , t 1 ] = E { X 1 , t 1 · E [ W i ( U 1 , t , ρ ) · B 1 | X 1 , t 1 ] + X 2 , t 1 · E [ W j ( U 2 , t , ρ ) · B 1 | X 2 , t 1 ] } = E { X 1 , t 1 · E [ W i ( U 1 , t , ρ ) ] · E [ B 1 | X 1 , t 1 ] + X 2 , t 1 · E [ W j ( U 2 , t , ρ ) ] · E [ B 1 | X 2 , t 1 ] } = E { E [ W i ( U 1 , t , ρ ) ] · X 1 , t 1 · E [ B 1 | X 1 , t 1 ] } + E [ W j ( U 2 , t , ρ ) ] · X 2 , t 1 · E [ B 1 | X 2 , t 1 ] } = ( 1 + ρ ) · E ( U 1 , t ) · E { E [ X 1 , t 1 · B 1 | X 1 , t 1 ] } + ( 1 + ρ ) · E ( U 2 , t ) · E { E [ X 2 , t 1 · B 1 | X 2 , t 1 ] } = ( 1 + ρ ) [ E ( U 1 , t ) E ( X 1 , t 1 B 1 ) + E ( U 2 , t ) E ( X 2 , t 1 B 1 ) ] .
Similarly,
E [ ( U 1 , t ρ   X 1 , t 1 + U 2 , t [ ρ ]   X 2 , t 1 ) · B 2 ] = ( 1 + ρ ) [ E ( U 1 , t ) E ( X 1 , t 1 B 2 ) + E ( U 2 , t ) E ( X 2 , t 1 B 2 ) ] . E [ ( V 1 , t ρ   X 1 , t 1 + V 2 , t [ ρ ]   X 2 , t 1 ) · B 1 ] = ( 1 + ρ ) [ E ( V 1 , t ) E ( X 1 , t 1 B 1 ) + E ( V 2 , t ) E ( X 2 , t 1 B 1 ) ] . E [ ( V 1 , t ρ   X 1 , t 1 + V 2 , t [ ρ ]   X 2 , t 1 ) · B 2 ] = ( 1 + ρ ) [ E ( V 1 , t ) E ( X 1 , t 1 B 2 ) + E ( V 2 , t ) E ( X 2 , t 1 B 2 ) ] .
Therefore,
E [ ( A t ρ   X t 1 ) B ] = ( 1 + ρ ) E ( U 1 , t ) E ( U 2 , t ) E ( V 1 , t ) E ( V 2 , t ) · E ( X 1 , t 1 B 1 ) E ( X 1 , t 1 B 2 ) E ( X 2 , t 1 B 1 ) E ( X 2 , t 1 B 2 ) = A ˜ E ( X t 1 B ) .
(iii)
E [ B ( A t ρ   X t 1 ) ] = E ( B X t 1 ) A ˜ for a random vector B independent of A t .
Proof. 
The proof process for this is similar to (ii), and we only provide key expressions.
E [ B 1 · ( U 1 , t ρ   X 1 , t 1 + U 2 , t ρ   X 2 , t 1 ) ] = E { E [ B 1 · ( U 1 , t ρ   X 1 , t 1 ) | X 1 , t 1 ] + E [ B 1 · ( U 2 , t ρ   X 2 , t 1 ) | X 2 , t 1 ] } = E E [ B 1 · i = 1 X 1 , t 1 W i ( U 1 , t , ρ ) | X 1 , t 1 ] + E [ B 1 · j = 1 X 2 , t 1 W j ( U 2 , t , ρ ) | X 2 , t 1 ] = E { X 1 , t 1 · E [ B 1 · W i ( U 1 , t , ρ ) | X 1 , t 1 ] + X 2 , t 1 · E [ B 1 · W j ( U 2 , t , ρ ) | X 2 , t 1 ] } = E { X 1 , t 1 · E [ W i ( U 1 , t , ρ ) ] · E [ B 1 | X 1 , t 1 ] + X 2 , t 1 · E [ W j ( U 2 , t , ρ ) ] · E [ B 1 | X 2 , t 1 ] } = E [ W i ( U 1 , t , ρ ) ] · E { X 1 , t 1 · E [ B 1 | X 1 , t 1 ] } + E [ W j ( U 2 , t , ρ ) ] · E { X 2 , t 1 · E [ B 1 | X 2 , t 1 ] } = E [ ( 1 + ρ ) · U 1 , t ] · E { E [ X 1 , t 1 · B 1 | X 1 , t 1 ] } + E [ ( 1 + ρ ) · U 2 , t ] · E { E [ X 2 , t 1 · B 1 | X 2 , t 1 ] } = ( 1 + ρ ) · E ( U 1 , t ) · E ( X 1 , t 1 · B 1 ) + ( 1 + ρ ) · E ( U 2 , t ) · E ( X 2 , t 1 · B 1 ) = ( 1 + ρ ) [ E ( U 1 , t ) · E ( X 1 , t 1 · B 1 ) + E ( U 2 , t ) · E ( X 2 , t 1 · B 1 ) ] ,
and
E [ B 1 · ( V 1 , t ρ   X 1 , t 1 + V 2 , t ρ   X 2 , t 1 ) ] = ( 1 + ρ ) [ E ( V 1 , t ) E ( X 1 , t 1 B 1 ) + E ( V 2 , t ) E ( X 2 , t 1 B 1 ) ] . E [ B 2 · ( U 1 , t ρ   X 1 , t 1 + U 2 , t ρ   X 2 , t 1 ) ] = ( 1 + ρ ) [ E ( U 1 , t ) E ( X 1 , t 1 B 2 ) + E ( U 2 , t ) E ( X 2 , t 1 B 2 ) ] . E [ B 2 · ( V 1 , t ρ   X 1 , t 1 + V 2 , t ρ   X 2 , t 1 ) ] = ( 1 + ρ ) [ E ( V 1 , t ) E ( X 1 , t 1 B 2 ) + E ( V 2 , t ) E ( X 2 , t 1 B 2 ) ] .
Therefore,
E [ B ( A t ρ   X t 1 ) ] = E B 1 B 2 · ( X 1 , t 1 , X 2 , t 1 ) · E ( 1 + ρ ) U 1 , t ( 1 + ρ ) V 1 , t ( 1 + ρ ) U 2 , t ( 1 + ρ ) V 2 , t = E ( B X t 1 ) A ˜ .
(iv)
E [ ( A t ρ   X t 1 ) ( A t ρ   X t 1 ) ] = A ˜ E ( X t 1 X t 1 ) A ˜ + C , where c 12 = 0 , c 21 = 0 ,
c 11 = ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) E ( X 1 , t 1 X 2 , t 1 ) 2 + α 1 ( 1 + ρ ) ( 1 + 2 ρ ( 1 + ρ ) α 1 ) [ p 1 E ( X 1 , t 1 ) + ( 1 p 1 ) E ( X 2 , t 1 ) ] c 22 = ( 1 + ρ ) 2 α 2 2 p 2 ( 1 p 2 ) E ( X 1 , t 1 X 2 , t 1 ) 2 + α 2 ( 1 + ρ ) ( 1 + 2 ρ ( 1 + ρ ) α 2 ) [ p 2 E ( X 1 , t 1 ) + ( 1 p 2 ) E ( X 2 , t 1 ) ] .
Proof. 
Based on the model definition, we expand the expectations.
E [ ( A t ρ   X t 1 ) ( A t ρ   X t 1 ) ] = E U 1 , t ρ   X 1 , t 1 + U 2 , t ρ   X 2 , t 1 V 1 , t ρ   X 1 , t 1 + V 2 , t ρ   X 2 , t 1 · U 1 , t ρ   X 1 , t 1 + U 2 , t ρ   X 2 , t 1 V 1 , t ρ   X 1 , t 1 + V 2 , t ρ   X 2 , t 1 e 11 e 12 e 21 e 22 .
For representational convenience, we denote the corresponding elements of the aforementioned matrix as e i , j , i , j = 1 , 2 . We then proceed to compute these elements individually.
e 11 = E [ ( U 1 , t ρ   X 1 , t 1 + U 2 , t ρ   X 2 , t 1 ) ( U 1 , t ρ   X 1 , t 1 + U 2 , t ρ   X 2 , t 1 ) ] = E [ ( U 1 , t ρ   X 1 , t 1 ) 2 + 2 ( U 1 , t ρ   X 1 , t 1 ) ( U 2 , t ρ   X 2 , t 1 ) + ( U 2 , t ρ   X 2 , t 1 ) 2 ] = E { i = 1 X 1 , t 1 W i ( U 1 , t , ρ ) 2 + 2 i = 1 X 1 , t 1 W i ( U 1 , t , ρ ) j = 1 X 2 , t 1 W j ( U 2 , t , ρ ) + j = 1 X 2 , t 1 W j ( U 2 , t , ρ ) 2 } = E i = 1 X 1 , t 1 W i 2 ( U 1 , t , ρ ) + i m X 1 , t 1 W i ( U 1 , t , ρ ) · W m ( U 1 , t , ρ ) + E { 2 i = 1 X 1 , t 1 j = 1 X 2 , t 1 W i ( U 1 , t , ρ ) × W j ( U 2 , t , ρ ) } + E j = 1 X 2 , t 1 W j 2 ( U 2 , t , ρ ) + j n X 2 , t 1 W j ( U 2 , t , ρ ) · W n ( U 2 , t , ρ ) s 1 + s 2 + s 3 .
s 1 = E i = 1 X 1 , t 1 W i 2 ( U 1 , t , ρ ) + i m X 1 , t 1 W i ( U 1 , t , ρ ) · W m ( U 1 , t , ρ ) = E E i = 1 X 1 , t 1 W i 2 ( U 1 , t , ρ ) + i m X 1 , t 1 W i ( U 1 , t , ρ ) · W m ( U 1 , t , ρ ) | X 1 , t 1 = E { X 1 , t 1 E [ W i 2 ( U 1 , t , ρ ) | X 1 , t 1 ] + X 1 , t 1 ( X 1 , t 1 1 ) E [ W i ( U 1 , t , ρ ) · W m ( U 1 , t , ρ ) | X 1 , t 1 ] } = E { X 1 , t 1 · E [ W i 2 ( U 1 , t , ρ ) ] + X 1 , t 1 ( X 1 , t 1 1 ) · E [ W i ( U 1 , t , ρ ) · W m ( U 1 , t , ρ ) ] } = E { X 1 , t 1 · ( 1 + ρ ) ( 1 + 2 ρ ) E ( U 1 , t ) + X 1 , t 1 ( X 1 , t 1 1 ) · ( 1 + ρ ) 2 E ( U 1 , t 2 ) } = E ( X 1 , t 1 ) · ( 1 + ρ ) ( 1 + 2 ρ ) E ( U 1 , t ) + E ( X 1 , t 1 2 X 1 , t 1 ) · ( 1 + ρ ) 2 E ( U 1 , t 2 ) = ( 1 + ρ ) ( 1 + 2 ρ ) α 1 p 1 · E ( X 1 , t 1 ) + ( 1 + ρ ) 2 α 1 2 p 1 2 · E ( X 1 , t 1 2 ) ( 1 + ρ ) 2 α 1 2 p 1 · E ( X 1 , t 1 ) = ( ( 1 + ρ ) ( 1 + 2 ρ ) α 1 p 1 ( 1 + ρ ) 2 α 1 2 p 1 ) E ( X 1 , t 1 ) + ( 1 + ρ ) 2 α 1 2 p 1 · E ( X 1 , t 1 2 ) ,
where
E [ W i 2 ( U 1 , t , ρ ) ] = E [ E [ W i 2 ( U 1 , t , ρ ) | U 1 , t ] ] = E [ V a r ( W i ( U 1 , t , ρ ) | U 1 , t ) + ( E ( W i ( U 1 , t , ρ ) | U 1 , t ) ) 2 ] = E [ U 1 , t ( 1 + ρ ) [ ρ + ( 1 + ρ ) ( 1 U 1 , t ) ] + ( U 1 , t ( 1 + ρ ) ) 2 ] = E [ U 1 , t ( 1 + ρ ) ( 1 + 2 ρ ) ] = ( 1 + ρ ) ( 1 + 2 ρ ) E ( U 1 , t ) = ( 1 + ρ ) ( 1 + 2 ρ ) α 1 p 1 .
E [ W i ( U 1 , t , ρ ) W j ( U 1 , t , ρ ) ] = E [ E [ W i ( U 1 , t , ρ ) W j ( U 1 , t , ρ ) | U 1 , t ] ] = E [ E [ W i ( U 1 , t , ρ ) | U 1 , t ] E [ W j ( U 1 , t , ρ ) | U 1 , t ] ] = E [ U 1 , t ( 1 + ρ ) · U 1 , t ( 1 + ρ ) ] = ( 1 + ρ ) 2 E ( U 1 , t 2 ) = ( 1 + ρ ) 2 α 1 2 p .
E [ W i ( U 1 , t , ρ ) W j ( U 2 , t , ρ ) ] = E [ E [ W i ( U 1 , t , ρ ) W j ( U 2 , t , ρ ) | ( U 1 , t , U 2 , t ) ] ] = E [ E [ W i ( U 1 , t , ρ ) | U 1 , t ] E [ W j ( U 2 , t , ρ ) | U 2 , t ] ] = E [ U 1 , t ( 1 + ρ ) · U 2 , t ( 1 + ρ ) ] = ( 1 + ρ ) 2 E ( U 1 , t U 2 , t ) = 0 .
E [ W i ( U 1 , t , ρ ) W n ( V 1 , t , ρ ) ] = E [ E [ W i ( U 1 , t , ρ ) W j ( V 1 , t , ρ ) | ( U 1 , t , V 1 , t ) ] ] = E [ E [ W i ( U 1 , t , ρ ) | U 1 , t ] E [ W j ( V 2 , t , ρ ) | V 1 , t ] ] = E [ U 1 , t ( 1 + ρ ) · V 1 , t ( 1 + ρ ) ] = ( 1 + ρ ) 2 E ( U 1 , t ) E ( V 1 , t ) = ( 1 + ρ ) 2 α 1 p 1 α 2 p 2 .
s 2 = E 2 i = 1 X 1 , t 1 j = 1 X 2 , t 1 W i ( U 1 , t , ρ ) W j ( U 2 , t , ρ ) = E E 2 i = 1 X 1 , t 1 j = 1 X 2 , t 1 W i ( U 1 , t , ρ ) W j ( U 2 , t , ρ ) | ( X 1 , t 1 , X 2 , t 1 ) = E { 2 X 1 , t 1 X 2 , t 1 · E [ W i ( U 1 , t , ρ ) W j ( U 2 , t , ρ ) | ( X 1 , t 1 , X 2 , t 1 ) ] } = 2 E ( X 1 , t 1 X 2 , t 1 ) · E [ W i ( U 1 , t , ρ ) W j ( U 2 , t , ρ ) ] = 2 E ( X 1 , t 1 X 2 , t 1 ) · ( 1 + ρ ) 2 E ( U 1 , t U 2 , t ) = 0 . s 3 = E j = 1 X 2 , t 1 W j 2 ( U 2 , t , ρ ) + j n X 2 , t 1 W j ( U 2 , t , ρ ) · W n ( U 2 , t , ρ ) = E { X 2 , t 1 · ( 1 + ρ ) ( 1 + 2 ρ ) E ( U 2 , t ) + X 2 , t 1 ( X 2 , t 1 1 ) · ( 1 + ρ ) 2 E ( U 2 , t 2 ) } = ( ( 1 + ρ ) ( 1 + 2 ρ ) α 1 ( 1 p 1 ) ( 1 + ρ ) 2 α 1 2 ( 1 p 1 ) ) E ( X 2 , t 1 ) + ( 1 + ρ ) 2 α 1 2 ( 1 p 1 ) E ( X 2 , t 1 2 ) .
e 11 = s 1 + s 2 + s 3 = ( ( 1 + ρ ) ( 1 + 2 ρ ) α 1 p 1 ( 1 + ρ ) 2 α 1 2 p 1 ) E ( X 1 , t 1 ) + ( 1 + ρ ) 2 α 1 2 p 1 · E ( X 1 , t 1 2 ) + ( ( 1 + ρ ) × ( 1 + 2 ρ ) α 1 ( 1 p 1 ) ( 1 + ρ ) 2 α 1 2 ( 1 p 1 ) ) E ( X 2 , t 1 ) + ( 1 + ρ ) 2 α 1 2 ( 1 p 1 ) · E ( X 2 , t 1 2 ) = [ ( 1 + ρ ) ( 1 + 2 ρ ) · α 1 p 1 ( 1 + ρ ) 2 α 1 2 p 1 ] E ( X 1 , t 1 ) + ( 1 + ρ ) 2 α 1 2 p 1 E ( X 1 , t 1 2 ) + [ ( 1 + ρ ) × ( 1 + 2 ρ ) α 1 ( 1 p 1 ) ( 1 + ρ ) 2 α 1 2 ( 1 p 1 ) ] E ( X 2 , t 1 ) + ( 1 + ρ ) 2 α 1 2 ( 1 p 1 ) E ( X 2 , t 1 2 ) = ( 1 + ρ ) α 1 [ ( 1 + 2 ρ ) ( 1 + ρ ) α 1 ] [ p 1 E ( X 1 , t 1 ) + ( 1 p 1 ) E ( X 2 , t 1 ) ] + ( 1 + ρ ) 2 α 1 2 [ p 1 × E ( X 1 , t 1 2 ) + ( 1 p 1 ) E ( X 2 , t 1 2 ) ] = ( 1 + ρ ) α 1 [ ( 1 + 2 ρ ) ( 1 + ρ ) α 1 ] [ p 1 E ( X 1 , t 1 ) + ( 1 p 1 ) E ( X 2 , t 1 ) ] + ( 1 + ρ ) 2 α 1 2 { p 1 [ p 1 × E ( X 1 , t 1 2 ) + ( 1 p 1 ) E ( X 2 , t 1 2 ) ] + ( 1 p 1 ) [ p 1 E ( X 1 , t 1 2 ) + ( 1 p 1 ) E ( X 2 , t 1 2 ) ] } = ( 1 + ρ ) α 1 [ ( 1 + 2 ρ ) ( 1 + ρ ) α 1 ] [ p 1 E ( X 1 , t 1 ) + ( 1 p 1 ) E ( X 2 , t 1 ) ] + ( 1 + ρ ) 2 α 1 2 [ p 1 2 × E ( X 1 , t 1 2 ) + ( 1 p 1 ) 2 E ( X 2 , t 1 2 ) ] + ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) [ E ( X 1 , t 1 2 ) + E ( X 2 , t 1 2 ) ] = ( 1 + ρ ) α 1 [ ( 1 + 2 ρ ) ( 1 + ρ ) α 1 ] [ p 1 E ( X 1 , t 1 ) + ( 1 p 1 ) E ( X 2 , t 1 ) ] + ( 1 + ρ ) 2 α 1 2 [ p 1 2 × E ( X 1 , t 1 2 ) + ( 1 p 1 ) 2 E ( X 2 , t 1 2 ) ] + ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) { [ E ( X 1 , t 1 ) E ( X 2 , t 1 ) ] 2 + 2 E ( X 1 , t 1 X 2 , t 1 ) } = ( 1 + ρ ) 2 [ α 1 2 p 2 E ( X 1 , t 1 2 ) + α 1 2 ( 1 p 1 ) 2 E ( X 2 , t 1 2 ) ] + 2 ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) × E ( X 1 , t 1 X 2 , t 1 ) + ( 1 + ρ ) α 1 [ ( 1 + 2 ρ ) ( 1 + ρ ) α 1 ] [ p 1 E ( X 1 , t 1 ) + ( 1 p 1 ) × E ( X 2 , t 1 ) ] + ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) [ E ( X 1 , t 1 ) E ( X 2 , t 1 ) ] 2 .
e 12 = E [ ( U 1 , t ρ   X 1 , t 1 + U 2 , t ρ   X 2 , t 1 ) ( V 1 , t ρ   X 1 , t 1 + V 2 , t ρ   X 2 , t 1 ) ] = E [ i = 1 X 1 , t 1 W i ( U 1 , t , ρ ) · m = 1 X 1 , t 1 W m ( V 1 , t , ρ ) + i = 1 X 1 , t 1 W i ( U 1 , t , ρ ) · n = 1 X 2 , t 1 W n ( V 2 , t , ρ ) + j = 1 X 2 , t 1 W j ( U 2 , t , ρ ) · m = 1 X 1 , t 1 W m ( V 1 , t , ρ ) + j = 1 X 2 , t 1 W j ( U 2 , t , ρ ) · n = 1 X 2 , t 1 W n ( V 2 , t , ρ ) ] s ˜ 1 + s ˜ 2 + s ˜ 3 + s ˜ 4 .
s ˜ 1 = E i = 1 X 1 , t 1 W i ( U 1 , t , ρ ) · m = 1 X 1 , t 1 W m ( V 1 , t , ρ ) = E E i = 1 X 1 , t 1 W i ( U 1 , t , ρ ) · m = 1 X 1 , t 1 W m ( V 1 , t , ρ ) | ( X 1 , t 1 , X 2 , t 1 ) = E { X 1 , t 1 2 E [ W i ( U 1 , t , ρ ) · W m ( V 1 , t , ρ ) ] | ( X 1 , t 1 , X 2 , t 1 ) } = E { X 1 , t 1 2 E [ W i ( U 1 , t , ρ ) · W m ( V 1 , t , ρ ) ] } = E ( X 1 , t 1 2 ) · ( 1 + ρ ) 2 E ( U 1 , t ) E ( V 1 , t ) = ( 1 + ρ ) 2 α 1 p 1 α 2 p 2 E ( X 1 , t 1 2 ) .
s ˜ 2 = E i = 1 X 1 , t 1 W i ( U 1 , t , ρ ) · n = 1 X 2 , t 1 W n ( V 2 , t , ρ ) = ( 1 + ρ ) 2 · E ( U 1 , t ) E ( V 2 , t ) · E ( X 1 , t 1 X 2 , t 1 ) = ( 1 + ρ ) 2 α 1 p 1 α 2 ( 1 p 2 ) E ( X 1 , t 1 X 2 , t 1 ) .
s ˜ 3 = E j = 1 X 2 , t 1 W j ( U 2 , t , ρ ) · m = 1 X 1 , t 1 W m ( V 1 , t , ρ ) = ( 1 + ρ ) 2 · E ( U 2 , t ) E ( V 1 , t ) · E ( X 1 , t 1 X 2 , t 1 ) = ( 1 + ρ ) 2 α 1 ( 1 p 1 ) α 2 p 2 E ( X 1 , t 1 X 2 , t 1 ) .
s ˜ 4 = E j = 1 X 2 , t 1 W j ( U 2 , t , ρ ) · n = 1 X 2 , t 1 W n ( V 2 , t , ρ ) = ( 1 + ρ ) 2 · E ( U 2 , t ) E ( V 2 , t ) · E ( X 2 , t 1 2 ) = ( 1 + ρ ) 2 α 1 ( 1 p 1 ) α 2 ( 1 p 2 ) E ( X 2 , t 1 2 ) .
e 12 = s ˜ 1 + s ˜ 2 + s ˜ 3 + s ˜ 4 = ( 1 + ρ ) 2 · α 1 p 1 α 2 p 2 · E ( X 1 , t 1 2 ) + ( 1 + ρ ) 2 · α 1 p 1 α 2 ( 1 p 2 ) · E ( X 1 , t 1 X 2 , t 1 ) + ( 1 + ρ ) 2 α 1 ( 1 p 1 ) α 2 p 2 E ( X 1 , t 1 X 2 , t 1 ) + ( 1 + ρ ) 2 α 1 ( 1 p 1 ) α 2 ( 1 p 2 ) E ( X 2 , t 1 2 ) . e 21 = E [ ( V 1 , t ρ   X 1 , t 1 + V 2 , t ρ   X 2 , t 1 ) ( U 1 , t ρ   X 1 , t 1 + U 2 , t ρ   X 2 , t 1 ) ] .
Similarly,
e 22 = ( 1 + ρ ) 2 [ α 2 2 q 2 E ( X 1 , t 1 2 ) + α 2 2 ( 1 p 2 ) 2 E ( X 2 , t 1 2 ) ] + 2 ( 1 + ρ ) 2 α 2 2 p 2 ( 1 p 2 ) E ( X 1 , t 1 X 2 , t 1 ) + ( 1 + ρ ) α 2 [ ( 1 + 2 ρ ) ( 1 + ρ ) α 2 ] [ p 2 E ( X 1 , t 1 ) + ( 1 p 2 ) E ( X 2 , t 1 ) ] + ( 1 + ρ ) 2 α 2 2 p 2 × ( 1 p 2 ) [ E ( X 1 , t 1 ) E ( X 2 , t 1 ) ] 2 .
Therefore,
E [ ( A t ρ   X t 1 ) ( A t ρ   X t 1 ) ] = e 11 e 12 e 21 e 22 = ( 1 + ρ ) 2 α 1 p 1 α 1 ( 1 p 1 ) α 2 p 2 α 2 ( 1 p 2 ) · E X 1 , t 1 2 X 1 , t 1 X 2 , t 1 X 1 , t 1 X 2 , t 1 X 2 , t 1 2 · α 1 p 1 α 2 p 2 α 1 ( 1 p 1 ) α 2 ( 1 p 2 ) + c 11 0 0 c 22
where c 12 = 0 , c 21 = 0 ,
c 11 = ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) E ( X 1 , t 1 X 2 , t 1 ) 2 + α 1 ( 1 + ρ ) ( 1 + 2 ρ ( 1 + ρ ) α 1 ) [ p 1 E ( X 1 , t 1 ) + ( 1 p 1 ) E ( X 2 , t 1 ) ] c 22 = ( 1 + ρ ) 2 α 2 2 p 2 ( 1 p 2 ) E ( X 1 , t 1 X 2 , t 1 ) 2 + α 2 ( 1 + ρ ) ( 1 + 2 ρ ( 1 + ρ ) α 2 ) [ p 2 E ( X 1 , t 1 ) + ( 1 p 2 ) E ( X 2 , t 1 ) ] .

Appendix B. Proof in Proposition 1

Proof. 
We begin by introducing a sequence of bivariate random series { X t ( m ) = ( X 1 , t ( m ) , X 2 , t ( m ) ) } m Z , defined as follows:
X t ( m ) = 0 , m < 0 , Z t , m = 0 , A t ρ   X t 1 ( m 1 ) + Z t , m > 0 ,
where, for any m, Z t is independent of A t ρ   X t 1 ( m 1 ) and X s ( m ) , s < t .
For any given m, we define a random vector X = ( X 1 , X 2 ) and the space L 2 ( Ω , F , P ) as the set of all random vectors X satisfying E ( X X ) < , where the measure between two random vectors X and Y is denoted as d ( X , Y ) = E ( X Y ) . lt is easy to obtain that L 2 ( Ω , F , P ) is a Hilbert space. Our goal is to establish the convergence of X t ( m ) X t as m . To achieve this, we aim to demonstrate that { X t ( m ) } forms a Cauchy sequence in the Hilbert space L 2 ( Ω , F , P ) . Now, let us proceed to the proof.
Existence.
We will prove the existence of the model through three key points.
1.
X t ( m ) is non-decreasing for all t.
To substantiate this claim, we must demonstrate that for all t Z and n 1 , X t ( m ) X t ( m 1 ) . For m = 1 , we have that
X t ( 1 ) = A t ρ   X t 1 ( 0 ) + Z t = A t ρ   Z t + Z t Z t = X t ( 0 ) .
Now suppose that X t ( k ) X t ( k 1 ) for k m and all t Z ; we will demonstrate that X t ( m + 1 ) X t ( m ) 0 . We consider the components X 1 , t ( m + 1 ) X 1 , t ( m ) and X 2 , t ( m + 1 ) X 2 , t ( m ) :
X 1 , t ( m + 1 ) X 1 , t ( m ) = U 1 , t ρ   X 1 , t 1 ( m ) U 1 , t ρ   X 1 , t 1 ( m 1 ) + U 2 , t ρ   X 2 , t 1 ( m ) U 2 , t ρ   X 2 , t 1 ( m 1 ) = i = 1 X 1 , t 1 ( m ) X 1 , t 1 ( m 1 ) W i ( U 1 , t , ρ ) + j = 1 X 2 , t 1 ( m ) X 2 , t 1 ( m 1 ) W j ( U 2 , t , ρ ) = U 1 , t ρ   ( X 1 , t 1 ( m ) X 1 , t 1 ( m 1 ) ) + U 2 , t ρ   ( X 2 , t 1 ( m ) X 2 , t 1 ( m 1 ) ) 0 .
Similarly,
X 2 , t ( m + 1 ) X 2 , t ( m ) = V 1 , t ρ   ( X 1 , t 1 ( m ) X 1 , t 1 ( m 1 ) ) + V 2 , t ρ   ( X 2 , t 1 ( m ) X 2 , t 1 ( m 1 ) ) 0 .
Therefore, by mathematical induction, we establish that { X t ( m ) } m Z is non-decreasing for all t.
2.
X t ( m ) L 2 ( Ω , F , P ) for m > 0 .
To demonstrate this, let us denote μ Z = E ( Z t ) . From Lemma 1, it follows that:
E ( X t ( m ) ) = A ˜ E ( X t 1 ( m 1 ) ) + μ Z = A ˜ [ A ˜ E ( X t 2 ( m 2 ) ) + μ Z ] + μ Z = A ˜ 2 E ( X t 2 ( m 2 ) ) + A ˜ μ Z + μ Z = A ˜ m E ( X t m 0 ) + A ˜ m 1 μ Z + + A ˜ μ Z + μ Z = ( A ˜ m + A ˜ m 1 + + A ˜ 0 ) · μ Z = ( j = 0 m A ˜ j ) · μ Z .
Assuming 0 < ( 1 + ρ ) α 1 , ( 1 + ρ ) α 2 < 1 , we find:
det ( I A ˜ ) = 1 ( 1 + ρ ) α 1 p 1 [ 1 ( 1 + ρ ) α 2 ] ( 1 + ρ ) α 2 [ 1 ( 1 ( 1 + ρ ) α 1 ) p 2 ] > 1 [ 1 ( 1 + ρ ) α 2 ] ( 1 + ρ ) α 2 = 0 ,
where I represents the 2 × 2 identity matrix. Hence, the matrix I A ˜ is invertible, and lim m j = 0 m A ˜ j = ( I A ˜ ) 1 . Consequently, E ( X t ( m ) ) becomes independent of t as E ( X t ( m ) ) < and tends to ( I A ˜ ) 1 μ Z as m .
Let us consider E [ X t ( m ) ( X t ( m ) ) ] . Leveraging the findings from Lemma 1, we derive:
E [ X t ( m ) ( X t ( m ) ) ] = E [ ( A t ρ   X t 1 ( m 1 ) + Z t ) ( A t ρ   X t 1 ( m 1 ) + Z t ) ] = E [ ( A t ρ   X t 1 ( m 1 ) ) ( A t ρ   X t 1 ( m 1 ) ) ] + E [ ( A t ρ   X t 1 ( m 1 ) ) Z t ] + E [ Z t ( A t ρ   X t 1 ( m 1 ) ) ] + E [ Z t Z t ] = A ˜ E ( X t 1 ( m 1 ) ( X t 1 ( m 1 ) ) ) A ˜ + C + A ˜ E ( X t 1 ( m 1 ) ) μ Z + μ Z E ( X t 1 ( m 1 ) ) A ˜ + E ( Z t Z t ) ,
where C = diag ( c 11 , c 22 ) , with
c 11 c 22 = ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) E ( X 1 , t 1 2 + X 2 , t 1 2 2 X 1 , t 1 X 2 , t 1 ) ( 1 + ρ ) 2 α 2 2 p 2 ( 1 p 2 ) E ( X 1 , t 1 2 + X 2 , t 1 2 2 X 1 , t 1 X 2 , t 1 ) + ( 1 + ρ ) ( 2 ρ + 1 ( 1 + ρ ) α 1 ) α 1 [ p 1 E ( X 1 , t 1 ) + ( 1 p 1 ) E ( X 2 , t 1 ) ] ( 1 + ρ ) ( 2 ρ + 1 ( 1 + ρ ) α 2 ) α 2 [ p 2 E ( X 1 , t 1 ) + ( 1 p 2 ) E ( X 2 , t 1 ) ] = ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) ( 1 + ρ ) 2 α 2 2 p 2 ( 1 p 2 ) ( 1 + ρ ) 2 α 2 2 p 2 ( 1 p 2 ) · E X 1 , t 1 2 X 2 , t 1 2 + ( 1 + ρ ) ( 2 ρ + 1 ( 1 + ρ ) α 1 ) α 1 p 1 ( 1 + ρ ) ( 2 ρ + 1 ( 1 + ρ ) α 1 ) α 1 ( 1 p 1 ) ( 1 + ρ ) ( 2 ρ + 1 ( 1 + ρ ) α 2 ) α 2 p 2 ( 1 + ρ ) ( 2 ρ + 1 ( 1 + ρ ) α 2 ) α 2 ( 1 p 2 ) · E X 1 , t 1 X 2 , t 1 2 ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) 2 ( 1 + ρ ) 2 α 2 2 p 2 ( 1 p 2 ) · E X 1 , t 1 X 2 , t 1 X 1 , t 1 X 2 , t 1 ( 1 + ρ ) α 1 p 1 ( 1 p 1 ) · 1 ( 1 + ρ ) α 1 ( 1 p 1 ) · 1 ( 1 + ρ ) α 2 p 2 ( 1 p 2 ) · 1 ( 1 + ρ ) α 2 ( 1 p 2 ) · 1 · E [ ( X t 1 ( m ) ) 2 ] + ( 1 + ρ ) ( 2 ρ + 1 ) α 1 p 1 ( 1 + ρ ) ( 2 ρ + 1 ) α 1 ( 1 p 1 ) ( 1 + ρ ) ( 2 ρ + 1 ) α 2 p 2 ( 1 + ρ ) ( 2 ρ + 1 ) α 2 ( 1 p 2 ) · E [ X t 1 ( m ) ] = A ˜ · E [ ( X t 1 ( m ) ) 2 ] + ( 2 ρ + 1 ) A ˜ · E [ X t 1 ( m ) ] .
If 0 < ( 1 + ρ ) α 1 , ( 1 + ρ ) α 2 < 1 , all entries of the matrix A ˜ fall within the interval ( 0 , 1 ) . Through iterative recursion repeated m times, we establish its independence from t. Consequently, E [ X t ( m ) ( X t ( m ) ) ] < . Hence, X t ( m ) L 2 ( Ω , F , P ) for m > 0 .
3.
X t ( m ) is a Cauchy sequence.
Let D m , t , k = X t ( m ) X t ( m k ) for all t Z , k = 1 , 2 , , m . From Equation (A1), it is straightforward to get
X t ( m ) X t ( m k ) = d A t ρ   ( X t 1 ( m 1 ) X t 1 ( m 1 k ) ) = A t ρ   D m 1 , t 1 , k .
Then, we have
E ( D m , t , k ) = A ˜ E ( D m 1 , t 1 , k ) = A ˜ m E ( Z t ) 0 , m .
Next, we aim to prove that
E ( D m , t , k · D m , t , k ) 0 , m .
Proof. 
Let H m , t , k [ ( X 1 , t ( m ) X 1 , t ( m k ) ) 2 , ( X 2 , t ( m ) X 2 , t ( m k ) ) 2 ] , D m , t , k = [ ( X 1 , t ( m ) X 1 , t ( m k ) ) , ( X 2 , t ( m ) X 2 , t ( m k ) ) ] . According to Lemma 1, it follows that
E ( D m , t , k · D m , t , k ) = E [ ( X t ( m ) X t ) · ( X t ( m ) X t ) ] = E [ A t ρ   ( X t 1 ( m 1 ) X t 1 ) · ( A t ρ   ( X t 1 ( m 1 ) X t 1 ) ) ] = A ˜ E [ ( X t 1 ( m 1 ) X t 1 ) · ( X t 1 ( m 1 ) X t 1 ) ] A ˜ + M ,
where M = diag ( m 1 , m 2 ) and
m 1 m 2 = ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) α 1 2 p ( 1 p 1 ) α 2 2 p 2 ( 1 p 2 ) α 2 2 p 2 ( 1 p 2 ) E ( H m 1 , t 1 , k ) + ( 1 + ρ ) α 1 p 1 ( ( 1 + 2 ρ ) α 1 ( 1 + ρ ) ) α 1 ( 1 p 1 ) ( ( 1 + 2 ρ ) α 1 ( 1 + ρ ) ) α 2 p 2 ( ( 1 + 2 ρ ) α 2 ( 1 + ρ ) ) α 2 ( 1 p 2 ) ( ( 1 + 2 ρ ) α 2 ( 1 + ρ ) ) E ( D m 1 , t 1 , k ) 2 ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) · E [ ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) ] α 2 2 p 2 ( 1 p 2 ) · E [ ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) ] .
Let us start by analyzing M . On the one hand, considering the non-negativity of the random variables X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) and X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) , we find that E ( ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) ) 0 . On the other hand, we proceed to derive E ( H m , t , k ) as follows:
E ( X 1 , t ( m ) X 1 , t ( m k ) ) 2 ( X 2 , t ( m ) X 2 , t ( m k ) ) 2 = E [ U 1 , t ρ   ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) + U 2 , t ρ   ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) ] 2 E [ V 1 , t ρ   ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) + V 2 , t ρ   ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) ] 2 = E [ U 1 , t ρ   ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) ] 2 E [ V 1 , t ρ   ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) ] 2 + E [ U 2 , t ρ   ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) ] 2 E [ V 2 , t ρ   ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) ] 2 + 2 E [ U 1 , t ρ   ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) · U 2 , t ρ   ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) ] 2 E [ V 1 , t ρ   ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) · V 2 , t ρ   ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) ] s 11 s 21 + s 12 s 22 + s 13 s 23 .
s 11 = E [ U 1 , t ρ   ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) ] 2 = E i = 1 X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) W i ( U 1 , t , ρ ) 2 = E i = 1 X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) W i 2 ( U 1 , t , ρ ) + i n X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) W i ( U 1 , t , ρ ) · W n ( U 1 , t , ρ ) = E E i = 1 X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) W i 2 ( U 1 , t , ρ ) | X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) + E E i n X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) W i ( U 1 , t , ρ ) · W n ( U 1 , t , ρ ) | X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) = E ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) · E [ W i 2 ( U 1 , t , ρ ) ] W i ( U 1 , t , ρ ) + E { ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) · ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) 1 ) · E [ W i ( U 1 , t , ρ ) · W n ( U 1 , t , ρ ) ] } = E ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) · E { U 1 , t ( 1 + ρ ) ( 1 + 2 ρ ) } + E [ ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) · ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) 1 ) ] · E [ U 1 , t 2 ( 1 + ρ ) 2 ] = E ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) · ( 1 + ρ ) ( 1 + 2 ρ ) · α 1 p 1 + E [ ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) · ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) 1 ) ] · ( 1 + ρ ) 2 · α 1 2 p 1 = α 1 p 1 ( 1 + ρ ) [ ( 1 + 2 ρ ) α 1 ( 1 + ρ ) ] · E ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) + α 1 2 p 1 ( 1 + ρ ) 2 · E ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) 2 .
Similarly,
s 12 = E [ U 2 , t ρ   ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) ] 2 = α 1 ( 1 p 1 ) ( 1 + ρ ) [ ( 1 + 2 ρ ) α 1 ( 1 + ρ ) ] · E ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) + α 1 2 ( 1 p 1 ) ( 1 + ρ ) 2 · E ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) 2 .
s 13 = 2 E [ U 1 , t ρ   ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) · U 2 , t ρ   ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) ] = 2 E { E [ U 1 , t ρ   ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) · U 2 , t ρ   ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) | ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) , X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) ] } = 2 E { E [ i = 1 X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) W i ( U 1 , t , ρ ) · j = 1 X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) W j ( U 2 , t , ρ ) | ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) , X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) ] } = 2 E { ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) · ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) E [ W i ( U 1 , t , ρ ) · W j ( U 2 , t , ρ ) ] } = 2 E { ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) } · E { E [ W i ( U 1 , t , ρ ) · W j ( U 2 , t , ρ ) ] } = 2 E [ ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) ] · E [ U 1 , t ( 1 + ρ ) · U 2 , t ( 1 + ρ ) ] = 2 ( 1 + ρ ) 2 E ( U 1 , t U 2 , t ) · E [ ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) ] = 0 = s 23 .
s 21 = E [ V 1 , t ρ   ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) ] 2 = α 2 p 2 ( 1 + ρ ) [ ( 1 + 2 ρ ) α 2 ( 1 + ρ ) ] · E ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) + α 2 2 p 2 ( 1 + ρ ) 2 · E ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) 2 .
s 22 = E [ V 2 , t ρ   ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) ] 2 = α 2 ( 1 p 2 ) ( 1 + ρ ) [ ( 1 + 2 ρ ) α 2 ( 1 + ρ ) ] · E ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) + α 2 2 ( 1 p 2 ) ( 1 + ρ ) 2 · E ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) 2 .
Then,
E ( H m , t , k ) = ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) α 1 2 p ( 1 p 1 ) α 2 2 p 2 ( 1 p 2 ) α 2 2 p 2 ( 1 p 2 ) · E ( X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) ) 2 ( X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ) 2 + ( 1 + ρ ) α 1 p 1 ( ( 1 + 2 ρ ) α 1 ( 1 + ρ ) ) α 1 ( 1 p 1 ) ( ( 1 + 2 ρ ) α 1 ( 1 + ρ ) ) α 2 p 2 ( ( 1 + 2 ρ ) α 2 ( 1 + ρ ) ) α 2 ( 1 p 2 ) ( ( 1 + 2 ρ ) α 2 ( 1 + ρ ) ) · E X 1 , t 1 ( m 1 ) X 1 , t 1 ( m k 1 ) X 2 , t 1 ( m 1 ) X 2 , t 1 ( m k 1 ) ( 1 + ρ ) α 1 p 1 · 1 α 1 ( 1 p 1 ) · 1 α 2 p 2 · 1 α 2 ( 1 p 2 ) · 1 · E ( H m 1 , t 1 , k ) + ( 1 + ρ ) α 1 p 1 ( 1 + 2 ρ ) α 1 ( 1 p 1 ) ( 1 + 2 ρ ) α 2 p 2 ( 1 + 2 ρ ) α 2 ( 1 p 2 ) ( 1 + 2 ρ ) E ( D m 1 , t 1 , k ) = A ˜ · E ( H m 1 , t 1 , k ) + ( 1 + 2 ρ ) A ˜ · E ( D m 1 , t 1 , k ) .
Therefore,
m 1 m 2 ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) α 1 2 p ( 1 p 1 ) α 2 2 p 2 ( 1 p 2 ) α 2 2 p 2 ( 1 p 2 ) · E ( H m 1 , t 1 , k ) + ( 1 + ρ ) α 1 p 1 ( 1 + 2 ρ α 1 ( 1 + ρ ) ) α 1 ( 1 p 1 ) ( 1 + 2 ρ α 1 ( 1 + ρ ) ) α 2 p 2 ( 1 + 2 ρ α 1 ( 1 + ρ ) ) α 2 ( 1 p 2 ) ( 1 + 2 ρ α 1 ( 1 + ρ ) ) · E ( D m 1 , t 1 , k ) < ( 1 + ρ ) α 1 p 1 · 1 α 1 ( 1 p 1 ) · 1 α 2 p 2 · 1 α 2 ( 1 p 2 ) · 1 · E ( H m 1 , t 1 , k ) + ( 1 + ρ ) α 1 p 1 · ( 1 + 2 ρ ) α 1 ( 1 p 1 ) · ( 1 + 2 ρ ) α 2 p 2 · ( 1 + 2 ρ ) α 2 ( 1 p 2 ) · ( 1 + 2 ρ ) · E ( D m 1 , t 1 , k ) = A ˜ · E ( H m 1 , t 1 , k ) + ( 1 + 2 ρ ) A ˜ · E ( D m 1 , t 1 , k ) A ˜ · A ˜ · E ( H m 2 , t 2 , k ) + ( 1 + 2 ρ ) A ˜ · E ( D m 2 , t 22 , k ) + ( 1 + 2 ρ ) A ˜ · E ( D m 1 , t 1 , k ) A ˜ m · E ( Z t 2 ) + ( 1 + 2 ρ ) j = 1 m A ˜ j E ( D m j , t j , k ) A ˜ m · E ( Z t 2 ) + ( 1 + 2 ρ ) m A ˜ m E ( Z t ) .
If 0 < ( 1 + ρ ) α 1 , ( 1 + ρ ) α 2 < 1 , it is straightforward to demonstrate that the eigenvalues (denoted as λ i , for i = 1 , 2 ) of the matrix A ˜ lie within the unit circle. Consequently, we have m λ i m 0 , i = 1 , 2 as m , which implies that Equation (A2) converges to 0 as m . Furthermore, M 0 as m .
Next, we examine
E ( D m , t , k · D m , t , k ) = A ˜ E ( D m 1 , t 1 , k · D m 1 , t 1 , k ) A ˜ + M A ˜ m E ( Z t Z t ) A ˜ m + i = 0 m 1 A ˜ i M A ˜ i .
As m , we have:
vec [ A ˜ m E ( Z t Z t ) A ˜ m ] = ( A ˜ m A ˜ m ) vec [ E ( Z t Z t ) ] = ( A ˜ A ˜ ) m vec [ E ( Z t Z t ) ] 0 ,
and
i = 0 m 1 vec ( A ˜ i M A ˜ i ) = i = 0 m 1 ( A ˜ A ˜ ) i vec ( M ) = I { A ˜ A ˜ } m I A ˜ A ˜ vec ( M ) { I A ˜ A ˜ } 1 .
The validity of Equations (A4) and (A5) relies on the fact that the eigenvalues of the matrix A ˜ A ˜ lie within the unit circle. Hence, E ( D m , t , k · D m , t , k ) 0 when m . □
We have proved E [ ( X t ( m ) X t ( m k ) ) · ( X t ( m ) X t ( m k ) ) ] 0 . This implies that the sequence { X t ( m ) } m Z is a Cauchy sequence, and thus, lim m X t ( m ) = X t L 2 ( Ω , F , P ) . Finally, by taking limits on both sides of Equation (A1) and letting m , we obtain X t = A t ρ   X t 1 + Z t , where Z t is independent of A t ρ   X t 1 and X s for s < t .
Uniqueness
Let us delve into uniqueness. Assume there exists another series { Y t } t Z satisfying Equation (4). Then, we can express the difference between X t and Y t as
X t Y t = A t ρ   X t A t ρ   Y t .
Define
B i , m = { ω : X i , t ( m ) ( ω ) Y i , t ( ω ) > 0 } , i = 1 , 2 , m 1 , B = i = 1 2 { ω : X i , t ( ω ) Y i , t ( ω ) > 0 } = i = 1 2 B i , = i = 1 2 n = 1 m = n B i , m .
We then establish:
P ( B i , m ) = P { ω : X i , t ( n ) ( ω ) Y i , t ( ω ) > 0 } k = 1 P { ω : X i , t ( m ) ( ω ) Y i , t ( ω ) = k } k = 1 k P { ω : X i , t ( m ) ( ω ) Y i , t ( ω ) = k } = E X i , t ( m ) Y i , t .
We introduce new notations:
L t ( 0 ) = Z t , L t ( m ) = ( L 1 , t ( m ) , L 2 , t ( m ) ) = ( | X 1 , t ( m ) Y 1 , t | , | X 2 , t ( m ) Y 2 , t | ) , for m 1 .
According to Lemma 1, we derive:
E [ L t ( m ) ] A ˜ E [ L t 1 ( m 1 ) ] A ˜ m E [ L t m ( 0 ) ] = A ˜ m μ Z t .
Consequently,
m = 1 Pr B i , m m = 1 A ˜ m μ Z t ( I A ˜ ) 1 μ Z t < .
By the Borel–Cantelli lemma, we conclude that Pr ( B i , ) = 0 . Thus, Pr ( B ) = 0 , i.e., X t = Y t almost surely.
Strictly stationary.
We will employ mathematical induction to demonstrate that for all h and k, ( X h + 1 ( m ) , , X h + k ( m ) ) and ( X 1 ( m ) , , X k ( m ) ) are identically distributed. Firstly, when m = 0 , we have
X h + 1 ( 0 ) X h + k ( 0 ) = Z h + 1 Z h + k
and
X 1 ( 0 ) X k ( 0 ) = Z 1 Z k .
Since ( Z h + 1 , , Z h + k ) = d ( Z 1 , , Z k ) (here, = d stands for having the same distribution), then ( X h + 1 ( m ) , , X h + k ( m ) ) and ( X 1 ( m ) , , X k ( m ) ) have identical distributions. Consequently, { X t ( 0 ) } t Z is strictly stationary.
Next, suppose { X t ( m 1 ) } t Z is strictly stationary; then we have
X 1 ( m ) X k ( m ) = Z 1 Z k + A 1 0 0 A k   ρ   X 1 ( m 1 ) X k ( m 1 ) ,
X h + 1 ( m ) X h + k ( m ) = Z h + 1 Z h + k + A h + 1 0 0 A h + k   ρ   X h + 1 ( m 1 ) X h + k ( m 1 ) .
Likewise, since { X t ( m 1 ) } t Z is strictly stationary, we have
( X 1 ( m 1 ) , , X k ( m 1 ) ) = d ( X h + 1 ( m 1 ) , , X h + k ( m 1 ) ) .
Then ( X 1 ( m ) , , X k ( m ) ) = d ( X h + 1 ( m ) , , X h + k ( m ) ) . Thus, { X t ( m ) } t Z forms a strictly stationary process. Furthermore, since lim m X t ( m ) = X t , i.e., X t ( m ) L 2 X t , then X t ( m ) P X t . Therefore, { X t } t Z is also a strictly stationary process.
Ergodicity.
At time t, the random matrical operation A t   ρ involves two random coefficient-thinning operations, i.e., “ U 1 , t   ρ ” or “ U 2 , t   ρ ” and “ V 1 , t   ρ ” or “ V 2 , t   ρ ”. Let W ( t ) denote all counting series involved in the matrix operation. Obviously, W ( t ) is a 2-dimensional series. Let σ ( X ) represent the σ -algebra rendering the vector X measurable. According to Equation (4), for any t, we have
σ ( X t , X t + 1 , ) σ ( Z t , A t , W ( t ) , Z t + 1 , A t + 1 , W ( t + 1 ) , ) ,
and consequently,
t = 1 σ ( X t , X t + 1 , ) D = t = 1 σ ( Z t , A t , W ( t ) , Z t + 1 , A t + 1 , W ( t + 1 ) , ) .
Given that the sequence { Z t , A t , W ( t ) } is an i.i.d. sequence of random vectors, then { Z t , A t , W ( t ) } is ergodic. According to Kolmogorov’s 0 1 law, for any event B within D , the probability P ( B ) is either 0 or 1. This means that the tail of the σ field of { X t } contains only the measure sets with probability 0 or 1. Consistent with findings akin to those in [8], the sequence X t is considered ergodic. □

Appendix C. Proof in Proposition 2

(i)
E ( X t ) = ( I A ˜ ) 1 E ( Z ) .
Proof. 
E ( X t ) = E ( A t   ρ   X t 1 + Z t ) = A ˜ E ( X t 1 ) + E ( Z t ) . □
(ii)
Σ X 1 , t = A ˜ Σ X t 1 A ˜ + C + Σ Z t , where c 12 = 0 , c 21 = 0 ,
c 11 = ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) E ( X 1 , t 1 X 2 , t 1 ) 2 + α 1 ( 1 + ρ ) ( 1 + 2 ρ ( 1 + ρ ) α 1 ) [ p 1 E ( X 1 , t 1 ) + ( 1 p 1 ) E ( X 2 , t 1 ) ] , c 22 = ( 1 + ρ ) 2 α 2 2 p 2 ( 1 p 2 ) E ( X 1 , t 1 X 2 , t 1 ) 2 + α 2 ( 1 + ρ ) ( 1 + 2 ρ ( 1 + ρ ) α 2 ) [ p 2 E ( X 1 , t 1 ) + ( 1 p 2 ) E ( X 2 , t 1 ) ] .
Proof. 
Var ( X 1 , t ) = Var ( U 1 , t   ρ   X 1 , t 1 + U 2 , t   ρ   X 2 , t 1 + Z 1 , t ) = Var ( U 1 , t   ρ   X 1 , t 1 ) + Var ( U 2 , t   ρ   X 2 , t 1 ) + 2 Cov ( U 1 , t   ρ   X 1 , t 1 , U 2 , t   ρ   X 2 , t 1 ) + σ Z 1 2 = ( 1 + ρ ) 2 α 1 2 [ p 1 2 Var ( X 1 , t 1 ) + ( 1 p 1 ) 2 Var ( X 2 , t 1 ) ] + ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) [ E ( X 1 , t 1 2 ) + E ( X 2 , t 1 2 ) 2 E ( X 1 , t 1 ) E ( X 2 , t 1 ) ] + [ ( 1 + ρ ) ( 1 + 2 ρ ) α 1 ( 1 + ρ ) 2 α 1 2 ] [ p 1 × E ( X 1 , t 1 ) + ( 1 p 1 ) E ( X 2 , t 1 ) ] + σ Z 1 2 ,
where
Var ( U 1 , t   ρ   X 1 , t 1 ) = E [ ( U 1 , t   ρ   X 1 , t 1 ) 2 ] [ E ( U 1 , t   ρ   X 1 , t 1 ) ] 2 = E i = 1 X 1 , t 1 W i ( U 1 , t , ρ ) 2 E ( i = 1 X 1 , t 1 W i ( U 1 , t , ρ ) ) 2 = E E i = 1 X 1 , t 1 W i 2 ( U 1 , t , ρ ) + i m X 1 , t 1 W i ( U 1 , t , ρ ) · W m ( U 1 , t , ρ ) | X 1 , t 1 E E ( i = 1 X 1 , t 1 W i ( U 1 , t , ρ ) | X 1 , t 1 ) 2 = E X 1 , t 1 ( 1 + ρ ) ( 1 + 2 ρ ) E ( U 1 , t ) + X 1 , t 1 ( X 1 , t 1 1 ) ( 1 + ρ ) 2 E ( U 1 , t 2 ) E X 1 , t 1 ( 1 + ρ ) E ( U 1 , t ) 2 = ( 1 + ρ ) 2 α 1 2 p 2 Var ( X 1 , t 1 ) + ( 1 + æ ) 2 ff 1 2 p ( 1 p 1 ) E ( X 1 , t 1 2 ) ( 1 + æ ) 2 ff 1 2 p 1 × E ( X 1 , t 1 ) 1 + ( 1 + ρ ) ( 1 + 2 ρ ) α 1 p 1 E ( X 1 , t 1 ) .
Var ( U 2 , t   ρ   X 1 , t 1 ) = ( 1 + ρ ) 2 α 1 2 ( 1 p 1 ) 2 Var ( X 2 , t 1 ) + ( 1 + æ ) 2 ff 1 2 p ( 1 p 1 ) E ( X 2 , t 1 2 ) ( 1 + ρ ) 2 α 1 2 ( 1 p 1 ) E ( X 2 , t 1 ) + ( 1 + ρ ) ( 1 + 2 ρ ) α 1 ( 1 p 1 ) E ( X 2 , t 1 ) .
Var ( X 2 , t ) = Var ( V 1 , t   ρ   X 1 , t 1 + V 2 , t   ρ   X 2 , t 1 + Z 2 , t ) = Var ( V 1 , t   ρ   X 1 , t 1 ) + Var ( V 2 , t   ρ   X 2 , t 1 ) + 2 Cov ( V 1 , t   æ   X 1 , t 1 , V 2 , t   æ   X 2 , t 1 ) + σ Z 2 , t 2 = ( 1 + ρ ) 2 α 2 2 [ p 2 2 Var ( X 1 , t 1 ) + ( 1 p 2 ) 2 Var ( X 2 , t 1 ) ] + ( 1 + æ ) 2 ff 2 2 p 2 ( 1 p 2 ) × [ E ( X 1 , t 1 2 ) + E ( X 2 , t 1 2 ) 2 E ( X 1 , t 1 ) E ( X 2 , t 1 ) ] + [ ( 1 + ρ ) ( 1 + 2 ρ ) α 2 ( 1 + ρ ) 2 α 2 2 ] [ p 2 E ( X 1 , t 1 ) + ( 1 p 2 ) E ( X 2 , t 1 ) ] + σ Z 2 , t 2 .
Cov ( X 1 , t , X 2 , t ) = ( 1 + ρ ) 2 α 1 p 1 α 2 p 2 Var ( X 1 , t 1 ) + ( 1 + æ ) 2 ff 1 ( 1 p 1 ) ff 2 ( 1 p 2 ) Var ( X 2 , t 1 ) + ( 1 + ρ ) 2 α 1 p 1 α 2 ( 1 p 2 ) Cov ( X 1 , t 1 , X 2 , t 1 ) + ( 1 + æ ) 2 ff 1 ( 1 p 1 ) ff 2 p 2 × Cov ( X 1 , t 1 , X 2 , t 1 ) .
Cov ( U 1 , t   ρ   X 1 , t 1 , U 2 , t   ρ   X 2 , t 1 ) = 2 ( 1 + ρ ) 2 α 1 p 1 α 1 ( 1 p 1 ) E ( X 1 , t 1 ) E ( X 2 , t 1 ) .
Therefore,
Σ X 1 , t = Var ( X 1 , t ) Cov ( X 1 , t , X 2 , t ) Cov ( X 2 , t , X 1 , t ) Var ( X 2 , t ) = ( 1 + ρ ) 2 α 1 p 1 α 1 ( 1 p 1 ) α 2 p 2 α 2 ( 1 p 2 ) Var ( X 1 , t 1 ) Cov ( X 1 , t , X 2 , t 1 ) Cov ( X 2 , t 1 , X 1 , t ) Var ( X 2 , t 1 ) · α 1 p 1 α 2 p 2 α 1 ( 1 p 1 ) α 2 ( 1 p 2 ) + C + Σ Z t = A ˜ Σ X t 1 A ˜ + diag ( c 11 , c 22 ) + Σ Z t ,
where
c 11 = ( 1 + ρ ) 2 α 1 2 p ( 1 p 1 ) E ( X 1 , t 1 X 2 , t 1 ) 2 + α 1 ( 1 + ρ ) ( 1 + 2 ρ ( 1 + ρ ) α 1 ) [ p 1 E ( X 1 , t 1 ) + ( 1 p 1 ) E ( X 2 , t 1 ) ] , c 22 = ( 1 + ρ ) 2 α 2 2 p 2 ( 1 p 2 ) E ( X 1 , t 1 X 2 , t 1 ) 2 + α 2 ( 1 + ρ ) ( 1 + 2 ρ ( 1 + ρ ) α 2 ) [ p 2 E ( X 1 , t 1 ) + ( 1 p 2 ) E ( X 2 , t 1 ) ] .

Appendix D. Proof in Lemma 3

Proof. 
Since
Cov ( X 1 , t , X 2 , t ) = ( 1 + ρ ) 2 α 1 α 2 [ p q Var ( X 1 , t 1 ) + ( 1 p ) ( 1 q ) Var ( X 2 , t 1 ) ] + ( 1 + ρ ) 2 α 1 α 2 × [ p 1 ( 1 p 2 ) + ( 1 p 1 ) p 2 ] Cov ( X 1 , t 1 , X 2 , t 1 ) ,
under the conditions of the lemma, we have that
( 1 + ρ ) 2 α 1 α 2 [ p 1 ( 1 p 2 ) + ( 1 p 1 ) p 2 ) ] < ( 1 + ρ ) 2 α 1 α 2 [ p 1 + ( 1 p 1 ) ] = 1 .
Building on the equation and the fact that the time series model { X 1 , t , X 2 , t } is stationary, we obtain that the covariance between the random variables X 1 , t and X 2 , t is given by
cov ( X 1 , t , X 2 , t ) = ( 1 + ρ ) 2 α 1 α 2 [ p 1 p 2 Var ( X 1 , t 1 ) + ( 1 p 1 ) ( 1 p 2 ) Var ( X 2 , t 1 ) ] 1 ( 1 + ρ ) 2 α 1 α 2 [ p 1 ( 1 p 2 ) + ( 1 p 1 ) p 2 ] .
Since { X 1 , t } = d { X 2 , t } = d Geom ( μ 1 + μ ) , then we obtain that the correlation coefficient between X 1 , t and X 2 , t is
γ = cov ( X 1 , t , X 2 , t ) D ( X 1 , t ) · D ( X 2 , t ) = ( 1 + ρ ) 2 α 1 α 2 [ p 1 p 2 + ( 1 p 1 ) ( 1 p 2 ) ] 1 ( 1 + ρ ) 2 α 1 α 2 [ p 1 ( 1 p 2 ) + ( 1 p 1 ) p 2 ] = ( 1 + ρ ) 2 α 1 α 2 ( 1 + ρ ) 2 α 1 α 2 ( p 1 + p 2 2 p 1 p 2 ) 1 ( 1 + ρ ) 2 α 1 α 2 ( p 1 + p 2 2 p 1 p 2 ) .

Appendix E. Proof in Theorem 3

Proof. 
From the proof of Theorem 2, S t = 1 2 ( X 1 , t + X 2 , t ) is also stationary and ergodic. Then, the process { S t μ } is a zero-mean, stationary, ergodic process. According to the Wold decomposition theorem (see [25] Section 2.6), this process can be represented as
S t μ = k = 0 d k ξ t k ,
where d 0 = 1 , k = 0 ( d k ) 2 < , and { ξ t } is white noise with parameters ( 0 , σ 2 ) . Then, the processes X i , t , i = 1 , 2 can be decomposed as
X 1 , t = a + k = d k ξ t k ,
where d k = 0 for k < 0 , and k = | d k | < . Also,
C o v ( S t , S t k ) = 1 4 C o v ( X 1 , t , X 1 , t k ) + C o v ( X 1 , t , X 2 , t k ) + C o v ( X 1 , t k , X 2 , t ) + C o v ( X 2 , t , X 2 , t k ) ) .
These terms represent the components of the matrix C o v ( X t , X t k ) = A ˜ k V a r ( X t ) . According to the correlation structure of the model and the properties of the matrix A ˜ , all terms in the equation are nonnegative, with the first and last terms being strictly positive. Equation (A6) indicates that k = C o v ( S t , S t k ) = σ 2 ( k = d k ) 2 , implying k = d k 0 . We can now apply the theorem presented in [26], thereby completing the proof. □

References

  1. Brannas, K.; Nordstrom, J. A Bivariate Integer Valued Allocation Model for Guest Nights in Hotels and Cottages. Umea Economic Studies Working Paper No. 547. 2001. Available online: https://ssrn.com/abstract=255292 (accessed on 20 May 2024). [CrossRef]
  2. Quoreshi, A.S. Bivariate time series modeling of financial count data. Commun. Stat. Theory Methods 2006, 35, 1343–1358. [Google Scholar] [CrossRef]
  3. Pedeli, X.; Karlis, D. A bivariate INAR (1) process with application. Stat. Model. 2011, 11, 325–349. [Google Scholar] [CrossRef]
  4. Nastić, A.S.; Ristić, M.M.; Popović, P.M. Estimation in a bivariate integer-valued autoregressive process. Commun. Stat. Theory Methods 2016, 45, 5660–5678. [Google Scholar] [CrossRef]
  5. Khan, N.M.; Oncel Cekim, H.; Ozel, G. The family of the bivariate integer-valued autoregressive process (BINAR (1)) with Poisson–Lindley (PL) innovations. J. Stat. Comput. Simul. 2020, 90, 624–637. [Google Scholar] [CrossRef]
  6. Chen, H.; Zhu, F.; Liu, X. A new bivariate INAR (1) model with time-dependent innovation vectors. Stats 2022, 5, 819–840. [Google Scholar] [CrossRef]
  7. Popović, P.M.; Ristić, M.M.; Nastić, A.S. A geometric bivariate time series with different marginal parameters. Stat. Pap. 2016, 57, 731–753. [Google Scholar] [CrossRef]
  8. Yu, M.; Wang, D.; Yang, K.; Liu, Y. Bivariate first-order random coefficient integer-valued autoregressive processes. J. Stat. Plan. Inference 2020, 204, 153–176. [Google Scholar] [CrossRef]
  9. Su, B.; Zhu, F. Comparison of BINAR (1) models with bivariate negative binomial innovations and explanatory variables. J. Stat. Comput. Simul. 2021, 91, 1616–1634. [Google Scholar] [CrossRef]
  10. Steutel, F.W.; van Harn, K. Discrete analogues of self-decomposability and stability. Ann. Probab. 1979, 7, 893–899. [Google Scholar] [CrossRef]
  11. Al-Osh, M.A.; Aly, E.-E.A. First order autoregressive time series with negative binomial and geometric marginals. Commun. Stat. Theory Methods 1992, 21, 2483–2492. [Google Scholar] [CrossRef]
  12. Borges, P.; Molinares, F.F.; Bourguignon, M. A geometric time series model with inflated-parameter Bernoulli counting series. Stat. Probab. Lett. 2016, 119, 264–272. [Google Scholar] [CrossRef]
  13. Kachour, M.; Truquet, L. A p-order signed integer-valued autoregressive (SINAR (p)) model. J. Time Ser. Anal. 2011, 32, 223–236. [Google Scholar] [CrossRef]
  14. Bulla, J.; Chesneau, C.; Kachour, M. A bivariate first-order signed integer-valued autoregressive process. Commun. Stat. Theory Methods 2017, 46, 6590–6604. [Google Scholar] [CrossRef]
  15. Zhang, Q.; Wang, D.; Fan, X. A negative binomial thinning-based bivariate INAR (1) process. Stat. Neerl. 2020, 74, 517–537. [Google Scholar] [CrossRef]
  16. Ristić, M.M.; Bakouch, H.S.; Nastić, A.S. A new geometric first-order integer-valued autoregressive (NGINAR (1)) process. J. Stat. Plan. Inference 2009, 139, 2218–2226. [Google Scholar] [CrossRef]
  17. Kolev, N.; Minkova, L.; Neytchev, P. Inflated-parameter family of generalized power series distributions and their application in analysis of overdispersed insurance data. ARCH Res. Clear. House 2000, 2, 295–320. [Google Scholar]
  18. Ristić, M.M.; Nastić, A.S.; Jayakumar, K.; Bakouch, H.S. A bivariate INAR (1) time series model with geometric marginals. Appl. Math. Lett. 2012, 25, 481–485. [Google Scholar] [CrossRef]
  19. Popović, P.M. A bivariate INAR (1) model with different thinning parameters. Stat. Pap. 2016, 57, 517–538. [Google Scholar] [CrossRef]
  20. Anderson, T.W.; Darling, D.A. A test of goodness of fit. J. Am. Stat. Assoc. 1954, 49, 765–769. [Google Scholar] [CrossRef]
  21. Gross, L. Tests for Normality, R Package Version 1.0-2. 2013. Available online: http://CRAN.R-project.org/package=nortest (accessed on 20 May 2024).
  22. Popović, P.M.; Nastić, A.S.; Ristić, M.M. Residual analysis with bivariate INAR (1) models. REVSTAT-Stat. J. 2018, 16, 349–363. [Google Scholar]
  23. Weiss, C.H.; Homburg, A.; Puig, P. Testing for zero inflation and overdispersion in inar (1) models. Stat. Pap. 2019, 60, 823–848. [Google Scholar] [CrossRef]
  24. Kang, Y.; Zhu, F.; Wang, D.; Wang, S. A zero-modified geometric INAR (1) model for analyzing count time series with multiple features. Can. J. Stat. 2023. [Google Scholar] [CrossRef]
  25. Brockwell, P.J.; Davis, R.A. Introduction to Time Series and Forecasting; Springer: Berlin/Heidelberg, Germany, 2002. [Google Scholar]
  26. Brockwell, P.J.; Davis, R.A. Time Series: Theory and Methods; Springer Science & Business Media: Berlin, Germany, 1991. [Google Scholar]
Figure 1. Variation of MAE and RMSE for Model (A) estimates across various sample sizes.
Figure 1. Variation of MAE and RMSE for Model (A) estimates across various sample sizes.
Axioms 13 00367 g001
Figure 2. Variation of MAE and RMSE for Model (B) estimates across various sample sizes.
Figure 2. Variation of MAE and RMSE for Model (B) estimates across various sample sizes.
Axioms 13 00367 g002
Figure 3. Variation of MAE and RMSE for Model (C) estimates across various sample sizes.
Figure 3. Variation of MAE and RMSE for Model (C) estimates across various sample sizes.
Axioms 13 00367 g003
Figure 4. Variation of MAE and RMSE for Model (D) estimates across various sample sizes.
Figure 4. Variation of MAE and RMSE for Model (D) estimates across various sample sizes.
Axioms 13 00367 g004
Figure 5. Variation of MAE and RMSE for Model (E) estimates across various sample sizes.
Figure 5. Variation of MAE and RMSE for Model (E) estimates across various sample sizes.
Axioms 13 00367 g005
Figure 6. Gaussian QQ plots of the estimates of α 1 , α 2 , and p 1 for Model (A) across various sample sizes.
Figure 6. Gaussian QQ plots of the estimates of α 1 , α 2 , and p 1 for Model (A) across various sample sizes.
Axioms 13 00367 g006
Figure 7. Gaussian QQ plots of the estimates of p 2 , ρ , and μ for Model (A) across various sample sizes.
Figure 7. Gaussian QQ plots of the estimates of p 2 , ρ , and μ for Model (A) across various sample sizes.
Axioms 13 00367 g007
Figure 8. Gaussian QQ plots of the estimates of α 1 , α 2 , and p 1 for Model (B) across various sample sizes.
Figure 8. Gaussian QQ plots of the estimates of α 1 , α 2 , and p 1 for Model (B) across various sample sizes.
Axioms 13 00367 g008
Figure 9. Gaussian QQ plots of the estimates of p 2 , ρ , and μ for Model (B) across various sample sizes.
Figure 9. Gaussian QQ plots of the estimates of p 2 , ρ , and μ for Model (B) across various sample sizes.
Axioms 13 00367 g009
Figure 10. Gaussian QQ plots of the estimates of α 1 , α 2 , and p 1 for Model (C) across various sample sizes.
Figure 10. Gaussian QQ plots of the estimates of α 1 , α 2 , and p 1 for Model (C) across various sample sizes.
Axioms 13 00367 g010
Figure 11. Gaussian QQ plots of the estimates of p 2 , ρ , and μ for Model (C) across various sample sizes.
Figure 11. Gaussian QQ plots of the estimates of p 2 , ρ , and μ for Model (C) across various sample sizes.
Axioms 13 00367 g011
Figure 12. Gaussian QQ plots of the estimates of α 1 , α 2 , and p 1 for Model (D) across various sample sizes.
Figure 12. Gaussian QQ plots of the estimates of α 1 , α 2 , and p 1 for Model (D) across various sample sizes.
Axioms 13 00367 g012
Figure 13. Gaussian QQ plots of the estimates of p 2 , ρ , and μ for Model (D) across various sample sizes.
Figure 13. Gaussian QQ plots of the estimates of p 2 , ρ , and μ for Model (D) across various sample sizes.
Axioms 13 00367 g013
Figure 14. Gaussian QQ plots of the estimates of α 1 , α 2 , and p 1 for Model (E) across various sample sizes.
Figure 14. Gaussian QQ plots of the estimates of α 1 , α 2 , and p 1 for Model (E) across various sample sizes.
Axioms 13 00367 g014
Figure 15. Gaussian QQ plots of the estimates of p 2 , ρ , and μ for Model (E) across various sample sizes.
Figure 15. Gaussian QQ plots of the estimates of p 2 , ρ , and μ for Model (E) across various sample sizes.
Axioms 13 00367 g015
Figure 16. Sample paths of OCND and OLNG series.
Figure 16. Sample paths of OCND and OLNG series.
Axioms 13 00367 g016
Figure 17. The autocorrelation function (ACF) and cross-correlation (CCF) plots of OCND and OLNG series.
Figure 17. The autocorrelation function (ACF) and cross-correlation (CCF) plots of OCND and OLNG series.
Axioms 13 00367 g017
Figure 18. Histograms of OCND and OLNG counts.
Figure 18. Histograms of OCND and OLNG counts.
Axioms 13 00367 g018
Figure 19. Sample paths of BETD and BETND series.
Figure 19. Sample paths of BETD and BETND series.
Axioms 13 00367 g019
Figure 20. The autocorrelation function (ACF) and cross-correlation (CCF) plots of BETD and BETND series.
Figure 20. The autocorrelation function (ACF) and cross-correlation (CCF) plots of BETD and BETND series.
Axioms 13 00367 g020
Figure 21. Histograms of BETD and BETND counts.
Figure 21. Histograms of BETD and BETND counts.
Axioms 13 00367 g021
Table 1. Simulation results for Model (A) and (B) under different sample sizes.
Table 1. Simulation results for Model (A) and (B) under different sample sizes.
Model (A)Model (B)
α ^ 1 α ^ 2 p ^ 1 p ^ 2 ρ ^ μ ^ α ^ 1 α ^ 2 p ^ 1 p ^ 2 ρ ^ μ ^
SizeMetrics0.300.250.200.150.105.000.300.250.200.150.305.00
150Est.0.29130.24100.16650.10280.13884.99290.29660.24130.17820.08230.31044.9649
Bias−0.0087−0.0090−0.0335−0.04720.0388−0.0071−0.0034−0.0087−0.0218−0.06770.0104−0.0351
MAE0.05000.04650.10360.13410.12840.32990.04490.04650.10470.16210.17760.3656
RMSE0.05550.05270.15590.20080.34830.41130.05470.05820.14760.44200.23440.4522
300Est.0.29330.24310.18870.13510.10715.00620.30450.25210.19830.14400.29215.0100
Bias−0.0067−0.0069−0.0113−0.01490.00710.00620.00450.0021−0.0017−0.0060−0.00790.0100
MAE0.03270.03260.06380.08020.07980.25060.02980.03120.06320.07460.11490.2724
RMSE0.03440.03540.08220.10590.10040.30800.03730.03870.08210.09800.14470.3337
600Est.0.29640.24570.19850.13980.10124.99580.30090.24970.20250.14560.29755.0169
Bias−0.0036−0.0043−0.0015−0.01020.0012−0.00420.0009−0.00030.0025−0.0044−0.00250.0169
MAE0.02340.02240.04630.05350.05790.16700.02090.01910.04200.04940.07890.1719
RMSE0.02590.02520.05680.06690.07060.21250.02590.02390.05400.06410.09900.2194
1200Est.0.29890.24920.19920.14490.10335.00210.30110.24950.19800.14750.29915.0012
Bias−0.0011−0.0008−0.0008−0.00510.00330.00210.0011−0.0005−0.0020−0.0025−0.00090.0012
MAE0.01460.01430.03030.03800.04230.12520.01520.01430.03090.03530.05710.1424
RMSE0.01730.01730.03770.04780.05200.15630.01920.01780.03910.04450.07240.1784
3000Est.0.29980.25050.19900.14920.10015.00580.29990.24970.19920.14850.29974.9988
Bias−0.00020.0005−0.0010−0.00080.00010.0058−0.0001−0.0003−0.0008−0.0015−0.0003−0.0012
MAE0.01180.01110.02420.02860.03350.09640.00970.00970.01820.02200.03430.0893
RMSE0.00950.00870.01930.02280.02660.07740.01220.01200.02290.02770.04210.1130
AD0.40140.44460.64560.46470.23350.16600.29790.48770.30910.45790.73240.4977
p-value0.35850.28320.09160.25340.79570.93950.58720.22270.55690.26330.05590.2106
Table 2. Simulation results of Model (C) under different sample sizes.
Table 2. Simulation results of Model (C) under different sample sizes.
Model (C)
α ^ 1 α ^ 2 p ^ 1 p ^ 2 ρ ^ μ ^
SizeMetrics0.400.400.200.150.255.00
150Est.0.39900.39620.18600.13200.25465.0048
Bias−0.0010−0.0038−0.0140−0.01800.00460.0048
MAE0.03830.03820.07670.08330.09850.4698
RMSE0.04690.04760.10180.11780.12340.5860
300Est.0.39970.40050.19570.14920.24855.0276
Bias−0.00030.0005−0.0043−0.0008−0.00150.0276
MAE0.02620.02700.05280.04970.06510.3519
RMSE0.03330.03390.06750.06420.08110.4382
600Est.0.40020.40070.19560.14920.25235.0087
Bias0.00020.0007−0.0044−0.00080.00230.0087
MAE0.01880.01900.03460.03510.04800.2383
RMSE0.02340.02390.04370.04330.06200.3015
1200Est.0.40040.40010.19780.15050.24794.9909
Bias0.00040.0001−0.00220.0005−0.0021−0.0091
MAE0.01320.01200.02440.02510.03330.1745
RMSE0.01660.01510.03090.03130.04150.2153
3000Est.0.40050.40040.20080.14930.24905.0044
Bias0.00050.00040.0008−0.0007−0.00100.0044
MAE0.00790.00790.01550.01500.02080.1091
RMSE0.01010.00980.01950.01860.02580.1347
AD0.27180.20610.43860.25950.48900.5042
p-value0.67030.86960.29280.71180.22110.2029
Table 3. Simulation results of Model (D) and (E) under different sample sizes.
Table 3. Simulation results of Model (D) and (E) under different sample sizes.
Model (D)Model (E)
α ^ 1 α ^ 2 p ^ 1 p ^ 2 ρ ^ μ ^ α ^ 1 α ^ 2 p ^ 1 p ^ 2 ρ ^ μ ^
SizeMetrics0.600.400.300.700.303.000.600.400.700.300.303.00
150Est.0.59890.39430.30250.69470.31062.96140.59610.39310.69870.28610.31112.9720
Bias−0.0011−0.00570.0025−0.00530.0106−0.0386−0.0039−0.0069−0.0013−0.01390.0111−0.0280
MAE0.03970.04410.05810.08560.08150.36660.03790.04580.05880.08560.08070.3892
RMSE0.04880.05470.07460.11030.11190.46510.04790.05830.07440.11110.11790.4905
300Est.0.60130.39910.29500.70520.29492.98460.59890.39860.70140.29710.29992.9953
Bias0.0013−0.0009−0.00500.0052−0.0051−0.0154−0.0011−0.00140.0014−0.0029−0.0001−0.0047
MAE0.02550.03020.04170.05540.04560.25890.02560.02770.03950.05820.04950.2759
RMSE0.03220.03780.05140.07130.05750.33110.03220.03570.04970.07250.06470.3487
600Est.0.59940.39750.29950.70070.30122.98960.59890.40070.70000.30070.29922.9992
Bias−0.0006−0.0025−0.00050.00070.0012−0.0104−0.00110.00070.00000.0007−0.0008−0.0008
MAE0.01730.02110.02850.04330.03440.17380.01760.02050.02890.03930.03500.1881
RMSE0.02190.02690.03550.05350.04300.21700.02190.02560.03630.04990.04390.2390
1200Est.0.60040.40130.29990.69810.29843.00390.59990.39790.69940.30010.29993.0026
Bias0.00040.0013−0.0001−0.0019−0.00160.0039−0.0001−0.0021−0.00060.0001−0.00010.0026
MAE0.01220.01440.02080.02730.02370.12700.01260.01560.02040.02900.02430.1447
RMSE0.01560.01790.02620.03540.02960.16020.01600.01970.02570.03650.03010.1791
3000Est.0.60000.40050.29980.69820.29983.00300.60070.39970.70000.30020.29823.0000
Bias0.00000.0005−0.0002−0.0018−0.00020.00300.0007−0.00030.00000.0002−0.00180.0000
MAE0.00810.00930.01260.01830.01520.08360.00790.00930.01220.01810.01450.0856
RMSE0.01020.01180.01550.02320.01920.10540.00980.01170.01530.02260.01800.1074
AD0.11860.39620.22910.13780.58150.70150.70480.42530.27180.24140.27660.5041
p-value0.98980.36880.80880.97630.12970.06670.06540.31500.67050.77130.65420.2031
Table 4. Descriptive statistics of OCND and OLNG series.
Table 4. Descriptive statistics of OCND and OLNG series.
CrimeMinMaxMedianMeanVar I d z 0
Offensive Conduct (OCND)0500.34480.55511.60980.6616
Offensive Language (OLNG)0400.29600.39921.34870.6438
Table 5. Fitting results of the monthly OCND and OLNG counts across different models.
Table 5. Fitting results of the monthly OCND and OLNG counts across different models.
Estimate ρ -BVGINAR(1)BVNGINAR(1)BVPOINAR(1)BVMIXINAR(1)
α ^ 1 0.12780.23640.07690.3306
α ^ 2 0.34350.16950.53400.2339
p ^ 1 0.57240.22890.69020.3071
p ^ 2 0.62630.25720.20780.0961
ρ ^ 0.0777
μ ^ 0.34480.34480.34480.3448
LogLik−489.64−492.15−492.25−526.14
AIC991.27994.31994.511062.29
BIC1008.531013.571013.771081.55
Table 6. Descriptive statistics of BETD and BETND series.
Table 6. Descriptive statistics of BETD and BETND series.
CrimeMinMaxMedianMeanVar I d z 0
Break and Enter Thefts into Dwellings (BETD)04045.543125.72444.64080.7587
Break and Enter Thefts into Non-Dwellings (BETND)02234.051715.84173.90990.7681
Table 7. Fitting results of the monthly BETD and BETND counts across different models.
Table 7. Fitting results of the monthly BETD and BETND counts across different models.
Estimate ρ -BVGINAR(1)BVNGINAR(1)BVPOINAR(1)BVMIXINAR(1)
α ^ 1 0.24340.70280.34830.6987
α ^ 2 0.26860.68550.56140.7911
p ^ 1 0.74510.39110.97410.4548
p ^ 2 0.24280.24310.32230.4693
ρ ^ 1.5649
μ ^ 5.54315.54315.54315.5431
LogLik−1742.96−1761.20−1987.37−2216.58
AIC3497.913532.403984.754443.16
BIC3515.173551.664004.014462.42
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Liu, C.; Wang, D. Bivariate Random Coefficient Integer-Valued Autoregressive Model Based on a ρ-Thinning Operator. Axioms 2024, 13, 367. https://doi.org/10.3390/axioms13060367

AMA Style

Liu C, Wang D. Bivariate Random Coefficient Integer-Valued Autoregressive Model Based on a ρ-Thinning Operator. Axioms. 2024; 13(6):367. https://doi.org/10.3390/axioms13060367

Chicago/Turabian Style

Liu, Chang, and Dehui Wang. 2024. "Bivariate Random Coefficient Integer-Valued Autoregressive Model Based on a ρ-Thinning Operator" Axioms 13, no. 6: 367. https://doi.org/10.3390/axioms13060367

APA Style

Liu, C., & Wang, D. (2024). Bivariate Random Coefficient Integer-Valued Autoregressive Model Based on a ρ-Thinning Operator. Axioms, 13(6), 367. https://doi.org/10.3390/axioms13060367

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop