VMSTA Modern Stochastics: Theory and Applications 2351-6054 2351-6046 2351-6046 VTeXMokslininkų g. 2A, 08412 Vilnius, Lithuania VMSTA130 10.15559/19-VMSTA130 Research Article A copula-based bivariate integer-valued autoregressive process with application ButeikisAndriusandrius.buteikis@mif.vu.lt∗ LeipusRemigijusremigijus.leipus@mif.vu.lt Institute of Applied Mathematics, Faculty of Mathematics and Informatics, Vilnius University, Lithuania Corresponding author. 2019 123201962227249 2182018 12122018 2812019 © 2019 The Author(s). Published by VTeX2019 Open access article under the CC BY license.

A bivariate integer-valued autoregressive process of order 1 (BINAR(1)) with copula-joint innovations is studied. Different parameter estimation methods are analyzed and compared via Monte Carlo simulations with emphasis on estimation of the copula dependence parameter. An empirical application on defaulted and non-defaulted loan data is carried out using different combinations of copula functions and marginal distribution functions covering the cases where both marginal distributions are from the same family, as well as the case where they are from different distribution families.

 Count data BINAR Poisson negative binomial distribution copula FGM copula Frank copula Clayton copula 60G10 62M10 62H12 Introduction

Different financial institutions that issue loans do this following company-specific (and/or country-defined) rules which act as a safeguard against loans issued to people who are known to be insolvent. However, striving for higher profits might motivate some companies to issue loans to higher risk clients. Usually company’s methods for evaluating loan risk are not publicly available. However, one way to evaluate if there aren’t too many knowingly very high-risk loans issued, and if insolvent clients are adequately separated from responsible clients, would be to look at the quantity of defaulted and non-defaulted loans issued each day. The adequacy of company’s rules for issuing loans can be analysed by modelling via copulas the dependence between the number of defaulted loans and the number of non-defaulted loans. The advantage of such approach is that copulas allow to model the marginal distributions (possibly from different distribution families) and their dependence structure (which is described via a copula) separately. Because of this feature, copulas were applied to many different fields, including survival analysis, hydrology, insurance risk analysis as well as finance (for examples of copula applications, see [3] or [4]), which also included the analysis of loans and their default rates.

The dependence of the default rate of loans on different credit risk categories was analysed in [5]. To model the dependence, copulas from ten different families were applied and three model selection tests were carried out. Because of the small sample size (24 observations per risk category) most of the copula families were not rejected and a single best copula model was not selected. To analyse whether dependence is affected by time, Fenech et al. [6] estimated the dependence among four different loan default indexes before the global financial crisis and after. They have found that the dependence was different in these periods. Four copula families were used to estimate the dependence between the default index pairs. While these studies were carried out for continuous data, discrete models created with copulas are less investigated: Genest and Nešlehová [8] discussed the differences and challenges of using copulas for discrete data compared to continuous data. Note that the previously mentioned studies assumed that the data does not depend on its own previous values. By using bivariate integer-valued autoregressive models (BINAR) it is possible to account for both the discreteness and autocorrelation of the data. Furthermore, copulas can be used to model the dependence of innovations in the BINAR(1) models: Karlis and Pedeli [10] used the Frank copula and the normal copula to model the dependence of the innovations of the BINAR(1) model.

In this paper we expand on using copulas in BINAR models by analysing additional copula families for the innovations of the BINAR(1) model and analyse different methods for BINAR(1) model parameter estimation. We also present a two-step method for the parameter estimation of the BINAR(1) model, where we estimate the model parameters separately from the dependence parameter of the copula. These estimation methods (including the one used in [10]) are compared via Monte Carlo simulations. Finally, in order to analyse the presence of autocorrelation and copula dependence in loan data, an empirical application is carried out for empirical weekly loan data.

The paper is organized as follows. Section 2 presents the BINAR(1) process and its main properties, Section 3 presents the main properties of copulas as well as some copula functions. Section 4 compares different estimation methods for the BINAR(1) model and the dependence parameter of copulas via Monte Carlo simulations. In Section 5 an empirical application is carried out using different combinations of copula functions and marginal distribution functions. Conclusions are presented in Section 6.

 The bivariate INAR(1) process

The BINAR(1) process was introduced in [18]. In this section we will provide the definition of the BINAR(1) model and will formulate its properties.

Let Rt=[R1,t,R2,t]′ , t∈Z , be a sequence of independent identically distributed (i.i.d.) nonnegative integer-valued bivariate random variables. A bivariate integer-valued autoregressive process of order 1 (BINAR(1)), Xt=[X1,t,X2,t]′ , t∈Z , is defined as: Xt=A∘Xt−1+Rt=α100α2∘X1,t−1X2,t−1+R1,tR2,t,t∈Z, where αj∈[0,1) , j=1,2 , and the symbol ‘∘’ is the thinning operator which also acts as the matrix multiplication. So the jth ( j=1,2 ) element is defined as an INAR process of order 1 (INAR(1)): Xj,t=αj∘Xj,t−1+Rj,t,t∈Z, where αj∘Xj,t−1:=∑i=1Xj,t−1Yj,t,i and Yj,t,1,Yj,t,2,… is a sequence of i.i.d. Bernoulli random variables with P(Yj,t,i=1)=αj=1−P(Yj,t,i=0) , αj∈[0,1) , such that these sequences are mutually independent and independent of the sequence Rt , t∈Z . For each t, Rt is independent of Xs , s<t .

Properties of the thinning operator are provided in [17] and [19] with proofs for selected few. We present the main properties of the thinning operator which will be used later on in the case of BINAR(1) model. Denote by ‘ =d ’ the equality of distributions.

(Thinning operator properties).

Let X,X1,X2 be nonnegative integer-valued random variables, such that EZ2<∞ , Z∈{X,X1,X2} , α,α1,α2∈[0,1) and let ‘∘’ be the thinning operator. Then the following properties hold:

α1∘(α2∘X)=d(α1α2)∘X ;

α∘(X1+X2)=dα∘X1+α∘X2 ;

E(α∘X)=αE(X) ;

Var(α∘X)=α2Var(X)+α(1−α)E(X) ;

E((α∘X1)X2)=αE(X1X2) ;

Cov(α∘X1,X2)=αCov(X1,X2) ;

E((α1∘X1)(α2∘X2))=α1α2E(X1X2) .

Xj,t , defined in eq. (2), has two random components: the survivors of the elements of the process at time t−1 , each with the probability of survival αj , which are denoted by αj∘Xj,t−1 , and the elements which enter in the system in the interval (t−1,t] , which are called arrival elements and denoted by Rj,t . We can obtain a moving average representation by substitutions and the properties of the thinning operator as in [1] or [11, p. 180]: Xj,t=αj∘Xj,t−1+Rj,t=d ∑k=0∞αjk∘Rj,t−k,j=1,2,t∈Z, where convergence on the right-hand side holds a.s.

Now we present some properties of the BINAR(1) model. They will be used when analysing some of parameter estimation methods. The proofs for these properties can be easily derived and some of them are provided in [17].

(Properties of the BINAR(1) process).

Let Xt=(X1,t,X2,t)′ be a nonnegative integer-valued time series given in Def. 1 and αj∈[0,1) , j=1,2 . Let Rt=(R1,t,R2,t)′ , t∈Z , be nonnegative integer-valued random variables with E(Rj,t)=λj and Var(Rj,t)=σj2<∞ , j=1,2 . Then the following properties hold:

EXj,t=μXj=λj1−αj ;

E(Xj,t|Xj,t−1)=αjXj,t−1+λj ;

Var(Xj,t)=σXj2=σj2+αjλj1−αj2 ;

Cov(Xi,t,Rj,t)=Cov(Ri,t,Rj,t) , i≠j ;

Cov(Xj,t,Xj,t+h)=αjhσXj2,h≥0 ;

Corr(Xj,t,Xj,t+h)=αjh , h≥0 ;

Cov(Xi,t,Xj,t+h)=αjh1−αiαjCov(Ri,t,Rj,t) , i≠j , h≥0 ;

Corr(Xi,t+h,Xj,t)=αih(1−αi2)(1−αj2)(1−αiαj)(σi2+αiλi)(σj2+αjλj)Cov(Ri,t,Rj,t) , i≠j , h≥0 ;

Similarly to (3), we have that Xt=d ∑k=0∞Ak∘Rt−k, where convergence on the right-hand side holds a.s.

Hence, the distributional properties of the BINAR(1) process can be studied in terms of Rt values. Note also, that according to [12], if αj∈[0,1) , j=1,2 , then there exists a unique stationary nonnegative integer-valued sequence Xt , t∈Z , satisfying (1).

From the covariance and correlation (see (g) and (h) in Theorem 2) of the BINAR(1) process we see that the dependence between X1,t and X2,t depends on the joint distribution of the innovations R1,t , R2,t . Pedeli and Karlis [18] analysed BINAR(1) models when the innovations were linked by either a bivariate Poisson or a bivariate negative binomial distribution, where the covariance of the innovations can be easily expressed in terms of their joint distribution parameters. Karlis and Pedeli [10] analysed two cases when the distributions of innovations of a BINAR(1) model are linked by either the Frank copula or a normal copula with either Poisson or negative binomial marginal distributions. We will expand their work by analysing additional copulas for the BINAR(1) model innovation distribution as well as estimation methods for the distribution parameters.

 Copulas

In this section we recall the definition and main properties of bivariate copulas, mainly following [8, 15] and [21] for the continuous and discrete settings.

Copula definition and properties

Copulas are used for modelling the dependence between several random variables. The main advantage of using copulas is that they allow to model the marginal distributions separately from their joint distribution. In this paper we are using two-dimensional copulas which are defined as follows:

A 2-dimensional copula C:[0,1]2→[0,1] is a function with the following properties:

for every u,v∈[0,1] : C(u,0)=C(0,v)=0;

for every u,v∈[0,1] : C(u,1)=u,C(1,v)=v;

for any u1,u2,v1,v2∈[0,1] such that u1≤u2 and v1≤v2 : C(u2,v2)−C(u2,v1)−C(u1,v2)+C(u1,v1)≥0 (this is also called the rectangle inequality).

The theoretical foundation of copulas is given by Sklar’s theorem:

 ([<xref ref-type="bibr" rid="j_vmsta130_ref_020">20</xref>]).

Let H be a joint cumulative distribution function (cdf) with marginal distributions F1,F2 . Then there exists a copula C such that for all (x1,x2)∈[−∞,∞]2 : H(x1,x2)=C(F1(x1),F2(x2)). If Fi is continuous for i=1,2 then C is unique; otherwise C is uniquely determined only on Ran(F1)×Ran(F2) , where Ran(F) denotes the range of the cdf F. Conversely, if C is a copula and F1,F2 are distribution functions, then the function H, defined by equation (7) is a joint cdf with marginal distributions F1,F2 .

If a pair of random variables (X1,X2) has continuous marginal cdfs Fi(x),i=1,2 , then by applying the probability integral transformation one can transform them into random variables (U1,U2)=(F1(X1),F2(X2)) with uniformly distributed marginals which can then be used when modelling their dependence via a copula. More about Copula theory, properties and applications can be found in [15] and [9].

 Copulas with discrete marginal distributions

Since innovations of a BINAR(1) model are nonnegative integer-valued random variables, one needs to consider copulas linking discrete distributions. In this section we will mention some of the key differences when copula marginals are discrete rather than continuous.

Firstly, as mentioned in Theorem 3, if F1 and F2 are discrete marginals then a unique copula representation exists only for values in the range of Ran(F1)×Ran(F2) . However, the lack of uniqueness does not pose a problem in empirical applications because it implies that there may exist more than one copula which describes the distribution of the empirical data. Secondly, regarding concordance and discordance, the discrete case has to allow for ties (i.e. when two variables have the same value), so the concordance measures (Spearman’s rho and Kendal’s tau) are margin-dependent, see [21]. There are several modifications proposed for Spearman’s rho, however, none of them are margin-free. Furthermore, Genest and Nešlehová [8] state that estimators of the dependence parameter θ based on Kendall’s tau or its modified versions are biased, and estimation techniques based on maximum likelihood are recommended. As such, we will not examine estimation methods based on concordance measures. Another difference from the continuous case is the use of the probability mass function (pmf) instead of the probability density function when estimating the model parameters which will be seen in Section 4.

 Some concrete copulas

In this section we will present several bivariate copulas, which will be used later when constructing and evaluating the BINAR(1) model. For all the copulas discussed, the following notation is used: u1:=F1(x1) , u2:=F2(x2) , where F1,F2 are marginal cumulative distribution functions (cdfs) of discrete random variables, and θ is the dependence parameter.

Farlie–Gumbel–Morgenstern copula

The Farlie–Gumbel–Morgenstern (FGM) copula has the following form: C(u1,u2;θ)=u1u2(1+θ(1−u1)(1−u2)). The dependence parameter θ can take values from the interval [−1,1] . If θ=0 , then the FGM copula collapses to independence. Note that the FGM copula can only model weak dependence between two marginals (see [15]). The copula when θ=0 is called a product (or independence) copula: C(u1,u2)=u1u2. Since the product copula corresponds to independence, it is important as a benchmark.

 Frank copula

The Frank copula has the following form: C(u1,u2;θ)=−1θlog(1+(exp(−θu1)−1)(exp(−θu2)−1)exp(−θ)−1). The dependence parameter θ can take values from (−∞,∞)∖{0} . The Frank copula allows for both positive and negative dependence between the marginals.

 Clayton copula

The Clayton copula has the following form: C(u1,u2;θ)=max{u1−θ+u2−θ−1,0}−1θ, with the dependence parameter θ∈[−1,∞)∖{0} . The marginals become independent when θ→0 . It can be used when the correlation between two random variables exhibits a strong left tail dependence – if smaller values are strongly correlated and hight values are less correlated. The Clayton copula can also account for negative dependence when θ∈[−1,0) . For more properties of this copula, see the recent paper by Manstavičius and Leipus [14].

 Parameter estimation of the copula-based BINAR(1) model

In this section we examine different BINAR(1) model parameter estimation methods and provide a two-step method for separate estimation of the copula dependence parameter. Estimation methods are compared via Monte Carlo simulations. Let Xt=(X1,t,X2,t)′ be a non-negative integer-valued time series given in Def. 1, where the joint distribution of (R1,t,R2,t)′ , with marginals F1,F2 , is linked by a copula C(·,·) : P(R1,t≤x1,R2,t≤x2)=C(F1(x1),F2(x2)) and let C(u1,u2)=C(u1,u2;θ) , where θ is a dependence parameter.

Conditional least squares estimation

The Conditional least squares (CLS) estimator minimizes the squared distance between Xt and its conditional expectation. Similarly to the method in [19] for the INAR(1) model, we construct the CLS estimator in the case of the BINAR(1) model.

Using Theorem 1 we can write the vector of conditional means as μt|t−1:=E(X1,t|X1,t−1)E(X2,t|X2,t−1)=α1X1,t−1+λ1α2X2,t−1+λ2, where λj:=ERj,t , j=1,2 . In order to calculate the CLS estimators of (α1,α2,λ1,λ2) we define the vector of residuals as the difference between the observations and their conditional expectation: Xt−μt|t−1=X1,t−α1X1,t−1−λ1X2,t−α2X2,t−1−λ2. Then, given a sample of N observations, X1,…,XN , the CLS estimators of αj,λj , j=1,2 , are found by minimizing the sum Qj(αj,λj):= ∑t=2N(Xj,t−αjXj,t−1−λj)2⟶minαj,λj,j=1,2. By taking the derivatives with respect to αj and λj , j=1,2 , and equating them to zero we get: αˆjCLS=∑t=2N(Xj,t−X¯j)(Xj,t−1−X¯j)∑t=2N(Xj,t−1−X¯j)2 and λˆjCLS=1N−1( ∑t=2NXj,t−αˆjCLS ∑t=2NXj,t−1). The asymptotic properties of the CLS estimators for the INAR(1) model case are provided in [13, 19, 2] and can be applied to the BINAR(1) parameter estimates, specified via equations (12) and (13). By the fact that the j-th component of the BINAR(1) process is an INAR(1) itself, we can formulate the following theorem for the marginal parameter vector distributions (see [2]):

Let Xt=(X1,t,X2,t)′ be defined in Def. 1 and let the parameter vector of (2) be (αj,λj)′ . Assume that αˆjCLS and λˆjCLS are the CLS estimators of αj and λj , j=1,2 . Then: NαˆjCLS−αjλˆjCLS−λj⟶dN(02,Bj), where Bj=EXj,t2EXj,tEXj,t1−1AjEXj,t2EXj,tEXj,t1−1,Aj=αj(1−αj)EXj,t3EXj,t2EXj,t2EXj,t+σj2EXj,t2EXj,tEXj,t1,j=1,2. Here, according to BINAR(1) properties in Theorem 2, EXj,t=λj1−αj,EXj,t2=σj2+αjλj1−αj2+λj2(1−αj)2,EXj,t3=ERj,t3−3σj2(1+λj)−λj3+2λj1−αj3+3σj2+αjλj1−αj2−2λj1−αj+3λj(σj2+αjλj)(1−αj)(1−αj2)+λj3(1−αj)3.

For the Poisson marginal distribution case the asymptotic variance matrix can be expressed as (see [7]) Bj=αj(1−αj)2λj+1−αj2−(1+αj)λj−(1+αj)λjλj+1+αj1−αjλj2,j=1,2. Furthermore, for a more general case, [12] proved that the CLS estimators of a multivariate generalized integer-valued autoregressive process (GINAR) are asymptotically normally distributed.

Note that E(X1,t−α1X1,t−1−λ1)(X2,t−α2X2,t−1−λ2)=Cov(R1,t,R2,t), which follows from E(X1,t−α1X1,t−1−λ1)(X2,t−α2X2,t−1−λ2)=E(α1∘X1,t−1−α1X1,t−1)(α2∘X2,t−1−α2X2,t−1)+E(α1∘X1,t−1−α1X1,t−1)(R2,t−λ2)+E(α2∘X2,t−1−α2X2,t−1)(R1,t−λ1)+E(R1,t−λ1)(R2,t−λ2) since the first three summands are zeros.

Assume that the joint pmf of (R1,t,R2,t) is given by bivariate Poisson distribution: P(R1,t=k,R2,t=l)= ∑i=0min{k,l}(λ1−λ)k−i(λ2−λ)l−iλi(k−i)!(l−i)!i!e−(λ1+λ2−λ), where k,l=0,1,... , λj>0 }0$]]>, j=1,2 , 0≤λ<min{λ1,λ2} . Then, for each j=1,2 , the marginal distribution of Rj,t is Poisson with parameter λj and Cov(R1,t,R2,t)=λ . If λ=0 then the two variables are independent.

Assume that the joint pmf of (R1,t,R2,t) is bivariate negative binomial distribution given by P(R1,t=k,R2,t=l)=Γ(β+k+l)Γ(β)k!l!(λ1λ1+λ2+β)k(λ2λ1+λ2+β)l×(βλ1+λ2+β)β, where k,l=0,1,... , λj>0 }0$]]>, j=1,2 , β>0 }0$]]>. Then, for each j=1,2 , the marginal distribution of Rj,t is negative binomial with parameters β and pj=β/(λj+β) and ERj,t=λj , Var(Rj,t)=λj(1+β−1λj) , Cov(R1,t,R2,t)=β−1λ1λ2 . Thus, bivariate negative binomial distribution is more flexible than bivariate Poisson due to overdispersion parameter β.

Assume now that the Poisson innovations R1,t and R2,t with parameters λ1 and λ2 , respectively, are linked by a copula with the dependence parameter θ. Taking into account equality (14), we can estimate θ by minimizing the sum of squared differences S= ∑t=2N(R1,tCLSR2,tCLS−γ(λˆ1CLS,λˆ2CLS;θ))2, where Rj,tCLS:=Xj,t−αˆjCLSXj,t−1−λˆjCLS,j=1,2,γ(λ1,λ2;θ):=Cov(R1,t,R2,t)= ∑k,l=1∞klc(F1(k;λ1),F2(l;λ2);θ)−λ1λ2. Here, c(F1(k;λ1),F2(s;λ2);θ) is the joint pmf: c(F1(k;λ1),F2(l;λ2);θ)=P(R1,t=k,R2,t=l)=C(F1(k;λ1),F2(s;λ2);θ)−C(F1(k−1;λ1),F2(l;λ2);θ)−C(F1(k;λ1),F2(l−1;λ2);θ)+C(F1(k−1;λ1),F2(l−1;λ2);θ),k≥1,l≥1.

Our estimation method is based on the approximation of covariance γ(λˆ1CLS,λˆ2CLS;θ) by γ(M1,M2)(λˆ1CLS,λˆ2CLS;θ)= ∑k=1M1 ∑l=1M2klc(F1(k;λˆ1CLS),F2(l;λˆ2CLS);θ)−λˆ1CLSλˆ2CLS. For example, if the marginals are Poisson with parameters λ1=λ2=1 and their joint distribution is given by the FGM copula in (8), then the covariance γ(M1,M2)(1,1;θ) stops changing significantly after setting M1=M2=M=8 , regardless of the selected dependence parameter θ. We used this approximation methodology when carrying out a Monte Carlo simulation in Section 4.4.

For the FGM copula, if we take the derivative of the sum S(M1,M2)= ∑t=2N(R1,tCLSR2,tCLS−γ(M1,M2)(λˆ1CLS,λˆ2CLS;θ))2, equate it to zero and use equation (17), we get θˆFGM=∑t=2N(X1,t−αˆ1CLSX1,t−1−λˆ1CLS)(X2,t−αˆ2CLSX2,t−1−λˆ2CLS)(N−1)∑k=1M1k(F1,kF‾1,k−F1,k−1F‾1,k−1)∑l=1M2l(F2,lF‾2,l−F2,l−1F‾2,l−1), where Fj,k:=Fj(k;λˆjCLS) , F‾j,k:=1−Fj,k , j=1,2 . The derivation of equation (19) is straightforward and thus omitted.

Depending on the selected copula family, calculation of (16) to get the analytical expression of the estimator θˆ may be difficult. However, we can use the function optim in the R statistical software to minimize (15). For other cases, where the marginal distribution has parameters other than expected value λj , equation (15) would need to be minimized by those additional parameters. For example, in the case of negative binomial marginals with corresponding mean λj and variance σj2 , i.e. when P(Rj,t=k)=Γ(k+λj2σj2−λj)Γ(λj2σj2−λj)k!(λjσj2)λj2σj2−λj(σj2−λjσj2)k,k=0,1,…,j=1,2, the additional parameters are σ12,σ22 , and the minimization problem becomes S(M1,M2)⟶minσ12,σ22,θ.

Conditional maximum likelihood estimation

BINAR(1) models can be estimated via conditional maximum likelihood (CML) (see [18] and [10]). The conditional distribution of the BINAR(1) process is: P(X1,t=x1,t,X2,t=x2,t|X1,t−1=x1,t−1,X2,t−1=x2,t−1)=P(α1∘x1,t−1+R1,t=x1,t,α2∘x2,t−1+R2,t=x2,t)= ∑k=0x1,t ∑l=0x2,tP(α1∘x1,t−1=k)P(α2∘x2,t−1=l)P(R1,t=x1,t−k,R2,t=x2,t−l). Here, αj∘x is the sum of x independent Bernoulli trials. Hence, P(αj∘xj,t−1=k)=xj,t−1kαjk(1−αj)xj,t−1−k,k=0,…,xj,t−1,j=1,2.

In the case of copula-based BINAR(1) model with Poisson marginals, P(R1,t=x1,t−k,R2,t=x2,t−l)=c(F1(x1,t−k,λ1),F2(x2,t−l,λ2);θ). Thus, we obtain P(X1,t=x1,t,X2,t=x2,t|X1,t−1=x1,t−1,X2,t−1=x2,t−1)= ∑k=0x1,t ∑l=0x2,tx1,t−1kα1k(1−α1)x1,t−1−kx2,t−1lα2l(1−α2)x2,t−1−l×c(F1(x1,t−k,λ1),F2(x2,t−l,λ2);θ) and the log conditional likelihood function, for estimating the marginal distribution parameters λ1,λ2 , the probabilities of the Bernoulli trial successes α1,α2 and the dependence parameter θ, is ℓ(α1,α2,λ1,λ2,θ)= ∑t=2NlogP(X1,t=x1,t,X2,t=x2,t|X1,t−1=x1,t−1,X2,t−1=x2,t−1) for some initial values x1,1 and x2,1 . In order to estimate the unknown parameters we maximize the log conditional likelihood: ℓ(α1,α2,λ1,λ2,θ)⟶maxα1,α2,λ1,λ2,θ. Numerical maximization is straightforward with the optim function in the R statistical software.

As for the CLS estimator, in other cases, where the marginal distribution has parameters other than λj , equation (20) would need to be maximized by those additional parameters. The CML estimator is asymptotically normally distributed under standard regularity conditions and its variance matrix is the inverse of the Fisher information matrix [18].

 Two-step estimation based on CLS and CML

Depending on the range of attainable values of the parameters and the sample size, CML maximization might take some time to compute. On the other hand, since CLS estimators of αj and λj are easily derived (compared to the CLS estimator of θ, which depends on the copula pmf form and needs to be numerically maximized), we can substitute the parameters of the marginal distributions in eq. (20) with CLS estimates from equations (12) and (13). Then we will only need to maximize ℓ with respect to a single dependence parameter θ for the Poisson marginal distribution case.

Summarizing, the two-step approach to estimating unknown parameters is to find (αˆjCLS,λˆjCLS)=argminQj(αj,λj),j=1,2, and to take these values as given in the second step: θˆCML=argmaxℓ(αˆ1CLS,αˆ2CLS,λˆ1CLS,λˆ2CLS,θ). For other cases of marginal distribution, any additional parameters, other than αj and λj would be estimated in the second step.

 Comparison of estimation methods via Monte Carlo simulation

We carried out a Monte Carlo simulation 1000 times to test the estimation methods with sample size 50 and 500. The generated model was a BINAR(1) with innovations joined by either the FGM, Frank or Clayton copula with Poisson marginal distributions, as well as with marginal distributions from different families: one is a Poisson distribution and the other is a negative binomial one. Note that for the two-step method only the estimates of θ and σ22 are included because estimated values of α1CLS,α2CLS,λ1CLS,λ2CLS are used in order to estimate the remaining parameters via CML.

Monte Carlo simulation results for a BINAR(1) model with Poisson innovations linked by the FGM, Frank or Clayton copula

Copula Sample size Parameter True value CLS CML Two-Step
MSE Bias MSE Bias MSE Bias
FGM N=50 α1 0.6 0.01874 −0.05823 0.00887 −0.01789 – –
α2 0.4 0.02033 −0.05223 0.01639 −0.02751 – –
λ1 1 0.12983 0.13325 0.06514 0.03366 – –
λ2 2 0.25625 0.16029 0.19939 0.07597 – –
θ −0.5 0.29789 0.12568 0.33840 0.07568 0.3311 0.0876
N=500 α1 0.6 0.00147 −0.00432 0.00073 −0.00122 – –
α2 0.4 0.00184 −0.00505 0.00129 −0.00157 – –
λ1 1 0.01012 0.00968 0.00556 0.00215 – –
λ2 2 0.02413 0.01843 0.01763 0.00678 – –
θ −0.5 0.04679 0.00668 0.04271 −0.00700 0.04265 −0.00443
Frank N=50 α1 0.6 0.02023 −0.06039 0.00950 −0.01965 – –
α2 0.4 0.02005 −0.05251 0.01630 −0.02858 – –
λ1 1 0.13562 0.13536 0.06740 0.03625 – –
λ2 2 0.25687 0.16392 0.19975 0.08291 – –
θ −1 1.83454 0.12394 2.05786 0.00860 1.97515 0.04216
N=500 α1 0.6 0.00153 −0.00595 0.00075 −0.00249 – –
α2 0.4 0.00181 −0.00582 0.00129 −0.00132 – –
λ1 1 0.01033 0.01269 0.00550 0.00421 – –
λ2 2 0.02442 0.02129 0.01785 0.00629 – –
θ −1 0.22084 0.01746 0.20138 −0.01779 0.20070 −0.01342
Clayton N=50 α1 0.6 0.01826 −0.05489 0.00799 −0.013295 – –
α2 0.4 0.01976 −0.05057 0.01585 −0.02427 – –
λ1 1 0.12679 0.12104 0.06080 0.01743 – –
λ2 2 0.25725 0.15704 0.19934 0.06499 – –
θ 1 0.71845 0.02621 0.72581 0.22628 0.62372 0.13283
N=500 α1 0.6 0.00146 −0.00518 0.00070 0.00016 – –
α2 0.4 0.00189 −0.00350 0.00120 −0.00049 – –
λ1 1 0.00973 0.01137 0.00513 −0.00150 – –
λ2 2 0.02447 0.01113 0.01707 0.00065 – –
θ 1 0.11578 0.03556 0.05864 0.04250 0.03199 −0.01342

The results for the Poisson marginal distribution case are provided in Table 1. The results for the case when one innovation follows a Poisson distribution and the other follows a negative binomial one are provided in Table 2. The lowest MSE values of θˆ are highlighted in bold. It is worth noting that CML estimation via numerical maximization depends heavily on the initial parameter values. If the initial values are selected too low or too high from the actual value, then the global maximum may not be found. In order to overcome this, we have selected the starting values equal to the CLS parameter estimates.

As can be seen in Table 1, the estimated values of αj and λj , j=1,2 , have a smaller bias and MSE when parameters are estimated via CML. On the other hand, estimation of θ via CLS exhibits a smaller MSE in the Frank copula case for smaller samples. For larger samples, the estimates of θ via the Two-step estimation method are very close to the CML estimates in terms of MSE and bias, and are closer to the true parameter values than the CLS estimates. Furthermore, since in the Two-step estimation numerical maximization is only carried out via a single parameter θ, the initial parameter values have less effect on the numerical maximization.

Monte Carlo simulation results for a BINAR(1) model with one innovation following a Poisson distribution and the other – a negative binomial one, where both innovations are linked by the FGM, Frank or Clayton copula

Copula Sample size Parameter True value CLS CML Two-Step
MSE Bias MSE Bias MSE Bias
FGM N=50 α1 0.6 0.01895 −0.05858 0.00845 −0.01513 – –
α2 0.4 0.01936 −0.04902 0.00767 −0.01953 – –
λ1 1 0.12940 0.12812 0.05424 0.01879 – –
λ2 2 0.39724 0.15151 0.24138 0.04833 – –
θ −0.5 0.31467 0.14070 0.29415 0.06674 0.29949 0.09693
σ22 9 27.87327 1.15731 15.12863 −0.14888 21.68229 0.72326
N=500 α1 0.6 0.00156 −0.00695 0.00076 −0.00153 – –
α2 0.4 0.00194 −0.00373 0.00053 0.00016 – –
λ1 1 0.01041 0.01201 0.00543 0.00290 – –
λ2 2 0.03882 0.01843 0.02362 −0.00057 – –
θ −0.5 0.06670 −0.02014 0.04298 −0.00268 0.04313 0.00562
σ22 9 6.24237 −1.99232 1.81265 0.00611 1.85222 −0.03506
Frank N=50 α1 0.6 0.02049 −0.06064 0.00912 −0.01594 – –
α2 0.4 0.01951 −0.04936 0.00772 −0.02070 – –
λ1 1 0.13769 0.13467 0.05748 0.02280 – –
λ2 2 0.40626 0.15408 0.23717 0.05534 – –
θ −1 1.81788 0.12516 1.75638 −0.01239 1.68019 0.06211
σ22 9 25.10400 0.49423 14.86812 −0.10034 21.92090 0.74026
N=500 α1 0.6 0.00161 −0.00702 0.00075 −0.00239 – –
α2 0.4 0.00187 −0.00364 0.00050 −0.00046 – –
λ1 1 0.01093 0.01652 0.00562 0.00501 – –
λ2 2 0.03728 0.01217 0.02335 0.00203 – –
θ −1 0.31942 −0.05593 0.18960 −0.01481 0.1902 −0.0079
σ22 9 4.82620 −1.75765 1.83082 0.02144 1.85852 −0.02690
Clayton N=50 α1 0.6 0.01987 −0.06159 0.00903 −0.01671 – –
α2 0.4 0.01879 −0.04928 0.00632 −0.01644 – –
λ1 1 0.13479 0.14072 0.06096 0.03052 – –
λ2 2 0.40675 0.14807 0.23171 0.02871 – –
θ 1 0.78497 0.07464 0.67837 0.21235 0.57454 0.10972
σ22 9 24.40051 0.17321 15.29879 −0.08379 23.73506 0.73754
N=500 α1 0.6 0.00153 −0.00722 0.00075 −0.00197 – –
α2 0.4 0.00196 −0.00385 0.00047 −0.00083 – –
λ1 1 0.01036 0.01745 0.00517 0.00409 – –
λ2 2 0.03999 0.01227 0.02304 0.00110 – –
θ 1 0.09927 0.04408 0.05557 0.03556 0.05559 0.02310
σ22 9 2.95995 −0.68733 1.79836 0.01348 1.87740 −0.02407

Table 2 demonstrates the estimation results when one innovation has a Poisson distribution and the other has a negative binomial one. With the inclusion of an additional variance parameter, the CLS estimation methods exhibit larger MSE and bias than the CML and Two-step estimation methods, for both the dependence and variance parameter estimates. Furthermore, the MSE of σˆ22 is smallest when the CML estimation method is used. On the other hand, both the Two-step and CML estimation methods produce similar estimates of θ in terms of MSE, regardless of sample size and copula function.

We can conclude that it is possible to accurately estimate the dependence parameter via CML using the CLS estimates of αˆj and λˆj . The resulting θˆ will be closer to the actual value of θ than θˆCLS and will not differ much from θˆCML . Additional inference on the bias of the estimates can be found in Appendix A.

 Application to loan default data

In this section we estimate a BINAR(1) model with the joint innovation distribution modelled by a copula cdf for empirical data. The data set consists of loan data which includes loans that have defaulted and loans that were repaid without missing any payments (non-defaulted loans). We will analyse and model the dependence between defaulted and non-defaulted loans as well as the presence of autocorrelation.

Loan default data

The data sample used is from Bondora, the Estonian peer-to-peer lending company. In November of 2014 Bondora introduced a loan rating system which assigns loans to different groups, based on their risk level. There are 8 groups ranging from the lowest risk group, ‘AA’, to the highest risk group, ‘HR’. However, the loan rating system could not be applied to most older loans due to a lack of data needed for Bondora’s rating model. Although Bondora issues loans in 4 different countries: Estonia, Finland, Slovakia and Spain, we will only focus on the loans issued in Spain. Since a new rating model indicates new rules for accepting or rejecting loans, we have selected the data sample from 21 October 2013, because from that date forward all loans had a rating assigned to them, to 1 January 2016. The time series are displayed in Figure 1. We are analysing data consisting of 115 weekly records.

‘CompletedLoans’ – the amount of non-defaulted loans issued per week which are repaid and have never defaulted (a loan that is 60 or more days overdue is considered defaulted);

‘DefaultedLoans’ – the amount of defaulted loans issued per week.

 The loan statistics are provided in Table 3:

Summary statistics of the weekly data of defaulted and non-defaulted loans issued in Spain

min max mean variance
DefaultedLoans 1.00 60.00 22.60 158.66
CompletedLoans 0.00 15.00 5.30 11.67

Bondora loan data: non-defaulted and defaulted loans by their issue date

AC function and PAC function plots of Bondora loan data

The mean, minimum, maximum and variance is higher for defaulted loans than for non-defaulted loans. As can be seen from Figure 2, the numbers of defaulted and non-defaulted loans might be correlated since they both exhibit increase and decrease periods at the same times.

The correlation between the two time series is 0.6684. We also note that the mean and variance are lower in the beginning of the time series. This feature could be due to various reasons: the effect of the new loan rating system, which was officially implemented in December of 2014, the effect of advertising or the fact that the amount of loans, issued to people living outside of Estonia, increased. The analysis of the significance of these effects is left for future research.

The sample autocorrelation (AC) function and the partial autocorrelation (PAC) function are displayed in Figure 2. We can see that the AC function is decaying over time and the PAC function has a significant first lag which indicates that the non-negative integer-valued time series could be autocorrelated.

In order to analyse if the amount of defaulted loans depends on the amount of non-defaulted loans on the same week, we will consider a BINAR(1) model with different copulas for the innovations. For the marginal distributions of the innovations we will consider the Poisson distribution as well as the negative binomial one. Our focus is the estimation of the dependence parameter, and we will use the Two-step estimation method, based on the Monte Carlo simulation results presented in Section 4.

 Estimated models

We estimated a number of BINAR(1) models with different distributions of innovations which include combinations of:

different copula functions: FGM, Frank or Clayton;

different combinations of the Poisson and negative binomial distributions: both marginals are Poisson, both marginals are negative binomial, or a mix of both.

In the first step of the Two-step method, we estimated αˆ1 and λˆ1 for non-defaulted loans, and αˆ2 and λˆ2 for defaulted loans via CLS. The results are provided in Table 4 with standard errors for the Poisson case in parenthesis:

Parameter estimates for BINAR(1) model via the Two-step estimation method: parameter CLS estimates from the first step with standard errors for the Poisson marginal distribution case in parenthesis

αˆ1 αˆ2 λˆ1 λˆ2
0.53134 0.75581 2.52174 5.58940
(0.08151) (0.06163) (0.45012) (1.41490)

Because the CLS estimation of parameters αj and λj , j=1,2 , does not depend on the selected copula and the marginal distribution family, these parameters will remain the same for each of the different distribution combinations for innovations. We can see that defaulted loans exhibit a higher degree of autocorrelation than non-defaulted loans do, due to a larger value of αˆ2 . The innovation mean parameter for defaulted loans is also higher, what indicates that random shocks have a larger effect on the number of defaulted loans.

The parameter estimation results from the second-step are provided in Table 5 with standard errors in parenthesis. σˆ12 is the innovation variance estimate of non-defaulted loans and σˆ22 is the innovation variance estimate of defaulted loans. According to [16], the observed Fisher information is the negative Hessian matrix, evaluated at the maximum likelihood estimator (MLE). The asymptotic standard errors reported in Table 5 are derived under the assumption that αj and λj , j=1,2 , are known, ignoring that the true values are substituted in the second step with their CLS estimates.

From the results in Table 5 we see that, according to the Akaike information criterion (AIC) and log-likelihood values, in most cases the FGM copula most accurately describes the relationship between the innovations of defaulted and non-defaulted loans, with the Frank copula being very close in terms of the AIC value. The Clayton copula is the least accurate in describing the innovation joint distribution, when compared to the FGM and Frank copula cases, which indicates that defaulted and non-defaulted loans do not exhibit strong left tail dependence.

Since the summary statistics of the data sample showed that the variance of the data is larger than the mean, a negative binomial marginal distribution may provide a better fit. Additionally, because copulas can link different marginal distributions, it is interesting to see if copulas with different discrete marginal distributions would also improve the model fit. BINAR(1) models where non-defaulted loan innovations are modelled with negative binomial distributions and defaulted loan innovations are modelled with Poisson marginal distributions, and vice versa, were estimated. In general, changing one of the marginal distributions to a negative binomial provides a better fit in terms of AIC than the Poisson marginal distribution case. However, the smallest AIC value is achieved when both marginal distributions are modelled with negative binomial distributions, linked via the FGM copula. Furthermore, the estimated innovation variance, σˆ22 , is much larger for defaulted loans, and this is similar to what we observed from the defaulted loan data summary statistics.

Parameter estimates for BINAR(1) model via Two-step estimation method: parameter CML estimates from the second-step for different innovation marginal and joint distribution combinations with standard errors in parenthesis, derived under the assumption that the values λˆj and αˆj , j=1,2 , from the first step are true

Marginals Copula θˆ σˆ12 σˆ22 AIC Log-likelihood
Both Poisson FGM 0.89270 – – 1763.48096 −880.74048
(0.18671)
Frank 2.38484 – – 1760.15692 −879.07846
(0.53367)
Clayton 0.39357 – – 1761.12369 −879.56185
(0.11697)
Negative binomial and Poisson FGM 1.00000 6.46907 – 1731.57339 −863.78670
(0.22914) (1.01114)
Frank 2.14329 6.10242 – 1731.95241 −863.97620
(0.45100) (1.15914)
Clayton 0.34540 5.73731 – 1736.47641 −866.23821
(0.12859) (0.52831)
Poisson and negative binomial FGM 1.00000 – 44.83107 1498.29563 −747.14782
(0.26357) (7.37423)
Frank 2.01486 – 44.10555 1498.81039 −747.40519
(0.61734) (7.33169)
Clayton 0.38310 – 43.42739 1503.55388 −749.77694
(0.17376) (7.29842)
Both negative binomial FGM 1.00000 6.55810 45.36834 1466.15418 −730.07709
(0.31675) (1.24032) (7.55217)
Frank 2.21356 6.58754 45.42601 1466.97947 −730.48973
(0.68192) (1.26126) (7.57743)
Clayton 0.55939 6.64478 45.78307 1470.73515 −732.36758
(0.24652) (1.25833) (7.66324)

Overall, both Frank and FGM copulas provide similar fit in terms of log-likelihood, regardless of the selected marginal distributions. We note, however, that for some FGM copula cases, the estimated value of parameter θ is equal to the maximal attainable value 1. Based on copula descriptions from Section 3, the FGM copula is used to model weak dependence. Given a larger sample size, the Frank copula might be more appropriate because it can capture a stronger dependence than the FGM copula can do. The negative binomial marginal distribution case θˆ≈2.21356 for the Frank copula indicates that there is a positive dependence between defaulted and non-defaulted loans, just as in the FGM copula case.

 Conclusions

The analysis via Monte Carlo simulations of different estimation methods shows that, although the estimates of BINAR(1) parameters via CML has the smallest MSE and bias, estimates of the dependence parameter has smaller differences of MSE and bias than for other estimation methods, indicating that estimations of the dependence parameter via different methods do not exhibit large differences. While CML estimates exhibit the smallest MSE, their calculation via numerical optimization relies on the selection of the initial parameter values. These values can be selected via CLS estimation.

An empirical application of BINAR models for loan data shows that, regardless of the selected marginal distributions, the FGM copula provides the best model fit in almost all cases. Models with the Frank copula are similar to FGM copula models in terms of AIC values. For some of these cases, the estimated FGM copula dependence parameter value was equal to the maximum that can be attained by an FGM copula. In such cases, a larger sample size could help to determine whether the FGM or Frank copula is more appropriate to model the dependence between amounts of defaulted and non-defaulted loans.

Although selecting marginal distributions from different families (Poisson or negative binomial) provided better models than those with only Poisson marginal distributions, the models with both marginal distributions modelled via negative binomial distributions provide the smallest AIC values which reflects overdispersion in amounts of both defaulted and non-defaulted loans. The FGM copula, which provides the best model fit, models variables which exhibit weak dependence. Furthermore, the estimated copula dependence parameter indicates that the dependence between amounts of defaulted and non-defaulted loans is positive.

Finally, one can apply some other copulas in order to analyse whether the loan data exhibits different forms of dependence from the ones discussed in this paper. Lastly, the approach can be extended by analysing the presence of structural changes within the data, or checking the presence of seasonality as well as extending the BINAR(1) model with copula joined innovations to account for the past values of other time series rather than only itself.

 Appendix

Standard errors of the bias of the estimated parameters from the Monte Carlo simulation

Copula Sample size Parameter True value CLS CML Two-Step
P-P P-NB P-P P-NB P-P P-NB
FGM N=50 α1 0.6 0.12396 0.12465 0.09252 0.09073 – –
α2 0.4 0.13274 0.13029 0.12510 0.08541 – –
λ1 1 0.33494 0.33631 0.25311 0.23225 – –
λ2 2 0.48040 0.61210 0.44024 0.48916 – –
θ −0.5 0.53139 0.54330 0.57707 0.53850 0.56899 0.53887
σ22 9 – 5.15368 – 3.88865 – 4.60221
N=500 α1 0.6 0.03813 0.03893 0.02706 0.02745 – –
α2 0.4 0.04258 0.04392 0.03585 0.02306 – –
λ1 1 0.10018 0.10076 0.07455 0.07367 – –
λ2 2 0.15433 0.19676 0.13266 0.15377 – –
θ −0.5 0.21631 0.25760 0.20666 0.20741 0.20657 0.20770
σ22 9 – 1.50841 – 1.34701 – 1.36119
Frank N=50 α1 0.6 0.12882 0.12975 0.09552 0.09420 – –
α2 0.4 0.13158 0.13073 0.12448 0.08543 – –
λ1 1 0.34266 0.34594 0.25719 0.23879 – –
λ2 2 0.47982 0.61879 0.43939 0.48409 – –
θ −1 1.34944 1.34314 1.43522 1.32589 1.40547 1.29538
σ22 9 – 4.98845 – 3.85654 – 4.62540
N=500 α1 0.6 0.03862 0.03951 0.02734 0.02727 – –
α2 0.4 0.04212 0.04312 0.03591 0.02240 – –
λ1 1 0.10091 0.10329 0.07409 0.07481 – –
λ2 2 0.15490 0.19278 0.13351 0.15287 – –
θ −1 0.46985 0.56268 0.44862 0.43540 0.44802 0.43627
σ22 9 – 1.31856 – 1.35359 – 1.36369
Clayton N=50 α1 0.6 0.12352 0.12684 0.08846 0.09360 – –
α2 0.4 0.13123 0.12798 0.12361 0.07779 – –
λ1 1 0.33505 0.33926 0.24609 0.24514 – –
λ2 2 0.48252 0.62066 0.44194 0.48075 – –
θ 1 0.84763 0.88328 0.82176 0.79618 0.77890 0.75037
σ22 9 – 4.93912 – 3.91243 – 4.81812
N=500 α1 0.6 0.03782 0.03850 0.02641 0.02742 – –
α2 0.4 0.04337 0.04410 0.03468 0.02176 – –
λ1 1 0.09804 0.10033 0.07162 0.07180 – –
λ2 2 0.15612 0.19969 0.13071 0.15185 – –
θ 1 0.33857 0.31212 0.23852 0.23316 0.23717 0.23476
σ22 9 – 1.57798 – 1.34163 – 1.37066

Let our Monte Carlo simulation data be Xj,1(i),…,Xj,N(i) for simulated sample i=1,…,M and j=1,2 . Let η∈{α1,α2,λ1,λ2,θ,σ22} and let ηˆ(i) be either a CLS, CML or Two-step estimate of the true parameter value η for the simulated sample i.

The mean squared error and the bias are calculated as follows: MSE(ηˆ)=1M ∑i=1M(ηˆ(i)−η)2,Bias(ηˆ)=1M ∑i=1M(ηˆ(i)−η). Calculating the per-sample bias for each simulated sample i would also allow us to calculate the sample variance of biases Bias(ηˆ(i))=ηˆ(i)−η , i=1,…,M : Varˆ(Bias(ηˆ))=1M−1 ∑i=1M[Bias(ηˆ(i))−Bias(ηˆ)]2, which we can use to calculate the standard error of the bias.

Kernel density estimate for the bias of the dependence parameter estimates in the Monte Carlo simulation

The standard errors of the parameter bias of the Monte Carlo simulation are presented in Table 6. The columns labelled ‘P-P’ indicate the cases where both innovations have Poisson marginal distributions, while columns labelled ‘P-NB’ is for the cases where one innovation component follows the Poisson distribution and the other follows a negative binomial one. The kernel density estimate for the bias of the dependence parameter estimate, θˆ , is presented in Figure 3 for the Monte Carlo simulation cases, where the sample size was 500.

The results in Table 6 are in line with the conclusions presented in Section 4.4 – for αˆj , λˆj , j=1,2 , and σˆ22 the standard error of the bias is smaller for CML than for CLS. On the other hand, θˆ has a similar standard error of the bias for CML and Two-step estimation methods. From Figure 3 we see that the CML and Two-step estimates of the dependence parameter θ are similar to each other and have a lower standard error of the bias than the CLS estimate.

 Acknowledgement

The authors would like to thank the anonymous referee for his/her feedback and constructive insights, which helped to improve this paper.

 References Al-Osh, M., Alzaid, A.: First-order integer-valued autoregressive (INAR(1)) process. J. Time Ser. Anal. 8, 261–275 (1987) MR0903755. https://doi.org/10.1111/j.1467-9892.1987.tb00438.x Barczy, M., Ispány, M., Pap, G., Scotto, M., Silva, M.E.: Innovational outliers in INAR(1) models. Commun. Stat., Theory Methods 39(18), 3343–3362 (2010) MR2747588. https://doi.org/10.1080/03610920903259831 Brigo, D., Pallavicini, A., Torresetti, R.: Credit Models and the Crisis: A Journey Into CDOs, Copulas, Correlations and Dynamic Models. Wiley, United Kingdom (2010) Cherubini, U., Mulinacci, S., Gobbi, F., Romagnoli, S.: Dynamic Copula Methods in Finance. Wiley, United Kingdom (2011) Crook, J., Moreira, F.: Checking for asymmetric default dependence in a credit card portfolio: A copula approach. J. Empir. Finance 18, 728–742 (2011) Fenech, J.P., Vosgha, H., Shafik, S.: Loan default correlation using an Archimedean copula approach: A case for recalibration. Econ. Model. 47, 340–354 (2015) Freeland, R.K., McCabe, B.: Asymptotic properties of CLS estimators in the Poisson AR(1) model. Stat. Probab. Lett. 73(2), 147–153 (2005) MR2159250. https://doi.org/10.1016/j.spl.2005.03.006 Genest, C., Nešlehová, J.: A primer on copulas for count data. ASTIN Bull. 37(2), 475–515 (2007) MR2422797. https://doi.org/10.2143/AST.37.2.2024077 Joe, H.: Dependence Modeling with Copulas. Chapman & Hall/CRC Monographs on Statistics and Applied probability 134 (2015) MR3328438 Karlis, D., Pedeli, X.: Flexible bivariate INAR(1) processes using copulas. Commun. Stat., Theory Methods 42, 723–740 (2013) MR3211946. https://doi.org/10.1080/03610926.2012.754466 Kedem, B., Fokianos, K.: Regression Models for Time Series Analysis. Wiley-Interscience, New Jersey (2002) MR1933755. https://doi.org/10.1002/0471266981 Latour, A.: The multivariate GINAR(p) process. Adv. Appl. Probab. 29(1), 228–248 (1997) MR1432938. https://doi.org/10.2307/1427868 Latour, A.: Existence and stochastic structure of a non-negative integer-valued autoregressive process. J. Time Ser. Anal. 19(4), 439–455 (1998) MR1652193. https://doi.org/10.1111/1467-9892.00102 Manstavičius, M., Leipus, R.: Bounds for the Clayton copula. Nonlinear Anal., Model. Control 22, 248–260 (2017) MR3608075. https://doi.org/10.15388/na.2017.2.7 Nelsen, R.: An Introduction to Copulas, 2nd Edition. Springer (2006) MR2197664. https://doi.org/10.1007/s11229-005-3715-x Pawitan, Y.: In All Likelihood: Statistical Modelling and Inference Using Likelihood. Oxford University Press, New York (2001) MR3668697. https://doi.org/10.1080/00031305.2016.1202140 Pedeli, X.: Modelling multivariate time series for count data. PhD thesis, Athens University of Economics and Business (2011) Pedeli, X., Karlis, D.: A bivariate INAR(1) process with application. Stat. Model., Int. J. 11(4), 325–349 (2011) MR2906704. https://doi.org/10.1177/1471082X1001100403 Silva, I.M.M.: Contributions to the analysis of discrete-valued time series. PhD thesis, University of Porto (2005) Sklar, M.: Fonctions de répartition à n dimensions et leurs marges. Publ. Inst. Stat. Univ. Paris 8, 229–231 (1959) MR0125600 Trivedi, P.K., Zimmer, D.M.: Copula modelling: An introduction for practitioners. Found. Trends Econom. 1(1), 1–111 (2007)