Goodness-of-fit Tests for the Bivariate Poisson Distribution

TL;DR

This paper introduces a new goodness-of-fit test for the bivariate Poisson distribution, capable of detecting various alternatives, and reviews existing tests and extensions, supported by simulation results.

Contribution

It proposes a novel, consistent goodness-of-fit test for the bivariate Poisson distribution that can detect local alternatives and employs bootstrap for null distribution estimation.

Findings

Test is consistent against any fixed alternative.

Bootstrap effectively estimates the null distribution.

Simulation shows good power for finite samples.

Abstract

The bivariate Poisson distribution is commonly used to model bivariate count data. In this paper we study a goodness-of-fit test for this distribution. We also provide a review of the existing tests for the bivariate Poisson distribution, and its multivariate extension. The proposed test is consistent against any fixed alternative. It is also able to detect local alternatives converging to the null at the rate $n^{-\frac{1}{2}}$. The bootstrap can be employed to consistently estimate the null distribution of the test statistic. Through a simulation study we investigated the goodness of the bootstrap approximation and the power for finite sample sizes.

Tables

Table 1. Table I: Simulation results for the probability of type I error, θ 1 = θ 2 = 1 subscript 𝜃 1 subscript 𝜃 2 1 \theta_{1}=\theta_{2}=1 .

		$n=30$			$n=50$			$n=70$
${\theta }_3=\rho$	Test	f05	f10	$KS$	f05	f10	$KS$	f05	f10	$KS$
0.25	$R_{n,(0,0)}$	0.037	0.087	0.863178	0.047	0.100	0.257432	0.044	0.086	0.111356
	$S_{n,(0,0)}$	0.046	0.089	0.934732	0.046	0.089	0.818621	0.045	0.092	0.508494
	$T_{n,(0,0)}$	0.043	0.087	0.329116	0.041	0.089	0.508494	0.047	0.098	0.718379
	$R_{n,(1,0)}$	0.035	0.099	0.329116	0.047	0.103	0.818621	0.038	0.084	0.329116
	$S_{n,(1,0)}$	0.034	0.090	0.902243	0.046	0.097	0.960002	0.041	0.089	0.508494
	$T_{n,(1,0)}$	0.042	0.094	0.329116	0.038	0.095	0.129364	0.046	0.094	0.369615
	$W_n$	0.022	0.056	1.00e-05	0.033	0.078	0.111356	0.038	0.090	0.612128
	$T$	0.011	0.031	2.2e-16	0.046	0.092	0.060937	0.013	0.038	2.2e-16
	$I_B$	0.027	0.061	2.2e-16	0.098	0.144	0.001642	0.022	0.054	2.2e-16
	$N{I}_B$	0.010	0.034	2.2e-16	0.068	0.111	0.003452	0.013	0.033	2.2e-16
0.50	$R_{n,(0,0)}$	0.048	0.112	0.129364	0.044	0.106	0.197933	0.045	0.098	0.559560
	$S_{n,(0,0)}$	0.041	0.094	0.049545	0.049	0.099	0.257432	0.049	0.099	0.413150
	$T_{n,(0,0)}$	0.041	0.098	0.257432	0.046	0.085	0.149677	0.055	0.111	0.413150
	$R_{n,(1,0)}$	0.051	0.101	0.129364	0.044	0.097	0.863178	0.047	0.109	0.172476
	$S_{n,(1,0)}$	0.042	0.099	0.069329	0.050	0.095	0.291736	0.046	0.099	0.197933
	$T_{n,(1,0)}$	0.044	0.095	0.413150	0.043	0.084	0.129364	0.051	0.112	0.459543
	$W_n$	0.022	0.061	0.013476	0.032	0.077	0.111356	0.037	0.081	0.111356
	$T$	0.026	0.049	1.40e-06	0.024	0.039	2.2e-16	0.021	0.053	0.000179
	$I_B$	0.088	0.125	2.2e-16	0.073	0.119	2.2e-16	0.051	0.081	2.2e-16
	$N{I}_B$	0.036	0.074	5.00e-07	0.018	0.049	2.2e-16	0.007	0.035	2.2e-16
0.75	$R_{n,(0,0)}$	0.043	0.089	0.718379	0.060	0.112	0.902243	0.050	0.114	0.508494
	$S_{n,(0,0)}$	0.050	0.092	0.818621	0.062	0.109	0.718379	0.052	0.104	0.612128
	$T_{n,(0,0)}$	0.045	0.084	0.665399	0.053	0.111	0.459543	0.045	0.104	0.226206
	$R_{n,(1,0)}$	0.049	0.090	0.995881	0.060	0.106	0.902243	0.051	0.116	0.612128
	$S_{n,(1,0)}$	0.049	0.088	0.818621	0.062	0.101	0.818621	0.051	0.106	0.459543
	$T_{n,(1,0)}$	0.044	0.084	0.863178	0.056	0.101	0.863178	0.052	0.104	0.459543
	$W_n$	0.029	0.076	0.024117	0.036	0.085	0.111356	0.038	0.088	0.129364
	$T$	0.025	0.049	2.2e-16	0.034	0.065	1.00e-07	0.024	0.058	5.30e-06
	$I_B$	0.116	0.140	2.2e-16	0.141	0.162	2.2e-16	0.129	0.153	2.2e-16
	$N{I}_B$	0.045	0.074	6.10e-06	0.033	0.081	2.2e-16	0.029	0.063	2.2e-16

Table 2. Table II: Simulation results for the probability of type I error, θ 1 = 1.5 , θ 2 = 1 formulae-sequence subscript 𝜃 1 1.5 subscript 𝜃 2 1 \theta_{1}=1.5,\theta_{2}=1 .

		$n=30$			$n=50$			$n=70$
${\theta }_3;\rho$	Test	f05	f10	$KS$	f05	f10	$KS$	f05	f10	$KS$
0.31; 0.25	$R_{n,(0,0)}$	0.054	0.104	0.226206	0.062	0.112	0.559560	0.052	0.111	0.069329
	$S_{n,(0,0)}$	0.054	0.103	0.718379	0.050	0.104	0.291736	0.050	0.095	0.863178
	$T_{n,(0,0)}$	0.056	0.091	0.718379	0.060	0.098	0.989545	0.057	0.107	0.902243
	$R_{n,(1,0)}$	0.050	0.094	0.149677	0.058	0.108	0.960002	0.051	0.107	0.013476
	$S_{n,(1,0)}$	0.047	0.095	0.226206	0.057	0.111	0.413150	0.053	0.097	0.508494
	$T_{n,(1,0)}$	0.043	0.091	0.369615	0.053	0.109	0.559560	0.058	0.115	0.559560
	$R_{n,(0,1)}$	0.050	0.106	0.459543	0.060	0.112	0.902243	0.056	0.109	0.016427
	$S_{n,(0,1)}$	0.051	0.105	0.459543	0.057	0.105	0.413150	0.056	0.099	0.508494
	$T_{n,(0,1)}$	0.052	0.096	0.863178	0.060	0.098	0.863178	0.056	0.106	0.718379
	$W_n$	0.022	0.066	0.041633	0.036	0.076	0.111356	0.037	0.082	0.111356
	$T$	0.018	0.046	1.00e-07	0.021	0.060	0.000318	0.022	0.064	0.009785
	$I_B$	0.031	0.060	2.2e-16	0.013	0.028	2.2e-16	0.007	0.014	2.2e-16
	$N{I}_B$	0.016	0.041	2.2e-16	0.010	0.018	2.2e-16	0.004	0.009	2.2e-16
0.62; 0.51	$R_{n,(0,0)}$	0.047	0.095	0.459543	0.045	0.095	0.863178	0.052	0.114	0.718379
	$S_{n,(0,0)}$	0.048	0.104	0.818621	0.049	0.091	0.818621	0.048	0.093	0.459543
	$T_{n,(0,0)}$	0.042	0.093	0.934732	0.043	0.098	0.934732	0.045	0.099	0.612128
	$R_{n,(1,0)}$	0.045	0.095	0.863178	0.044	0.096	0.978036	0.056	0.101	0.559560
	$S_{n,(1,0)}$	0.051	0.088	0.459543	0.045	0.086	0.718379	0.048	0.100	0.459543
	$T_{n,(1,0)}$	0.037	0.096	0.769894	0.041	0.088	0.291736	0.047	0.109	0.226206
	$R_{n,(0,1)}$	0.049	0.097	0.413150	0.045	0.101	0.718379	0.054	0.104	0.902243
	$S_{n,(0,1)}$	0.052	0.098	0.508494	0.042	0.098	0.612128	0.051	0.091	0.413150
	$T_{n,(0,1)}$	0.043	0.089	0.902243	0.049	0.095	0.508494	0.046	0.088	0.329116
	$W_n$	0.026	0.055	0.003013	0.037	0.071	0.111356	0.039	0.079	0.111356
	$T$	0.056	0.088	0.000526	0.050	0.104	0.011917	0.049	0.096	0.001109
	$I_B$	0.147	0.201	2.2e-16	0.169	0.223	2.2e-16	0.147	0.196	2.2e-16
	$N{I}_B$	0.094	0.152	0.000622	0.082	0.145	0.006666	0.076	0.120	0.078967
0.92; 0.75	$R_{n,(0,0)}$	0.057	0.102	0.612128	0.054	0.097	0.413150	0.046	0.090	0.863178
	$S_{n,(0,0)}$	0.052	0.108	0.413150	0.050	0.091	0.769894	0.044	0.094	0.559560
	$T_{n,(0,0)}$	0.043	0.098	0.934732	0.056	0.102	0.226206	0.041	0.085	0.413150
	$R_{n,(1,0)}$	0.053	0.104	0.559560	0.055	0.103	0.508494	0.043	0.088	0.863178
	$S_{n,(1,0)}$	0.049	0.103	0.769894	0.050	0.093	0.612128	0.045	0.091	0.612128
	$T_{n,(1,0)}$	0.040	0.107	0.665399	0.052	0.110	0.172476	0.037	0.084	0.459543
	$R_{n,(0,1)}$	0.055	0.110	0.459543	0.050	0.094	0.257432	0.044	0.082	0.818621
	$S_{n,(0,1)}$	0.051	0.108	0.665399	0.048	0.087	0.369615	0.045	0.091	0.508494
	$T_{n,(0,1)}$	0.046	0.089	0.769894	0.055	0.109	0.508494	0.043	0.088	0.665399
	$W_n$	0.037	0.081	0.000714	0.042	0.079	0.111356	0.037	0.083	0.149677
	$T$	0.029	0.059	1.70e-06	0.057	0.094	0.008821	0.078	0.109	0.065401
	$I_B$	0.091	0.116	2.2e-16	0.209	0.239	2.2e-16	0.196	0.220	2.2e-16
	$N{I}_B$	0.021	0.051	2.2e-16	0.089	0.152	0.001554	0.094	0.149	0.003483

Table 3. Table III: Simulation results for the power ( n = 50 𝑛 50 n=50 ).

Alternative	$\frac{var\left(X_1\right)}{E\left(X_1\right)}$	$\frac{var\left(X_2\right)}{E\left(X_2\right)}$	$\rho$	$R_{n,(0,0)}$	$S_{n,(0,0)}$	$W_n$	$T_{n,(0,0)}$	$T$	$I_B$	$N{I}_B$
$BB(1;0.41,0.02,0.01)$	0.590	0.980	0.026	0.860	0.871	0.829	0.857	0.103	0.000	0.000
$BB(1;0.41,0.03,0.02)$	0.590	0.970	0.092	0.859	0.879	0.779	0.893	0.122	0.000	0.000
$BB(2;0.42,0.02,0.01)$	0.580	0.980	0.023	0.726	0.677	0.682	0.746	0.251	0.005	0.005
$BB(2;0.51,0.01,0.01)$	0.490	0.990	0.099	0.900	0.862	0.847	0.887	0.656	0.001	0.001
$BB(2;0.61,0.01,0.01)$	0.390	0.990	0.080	0.974	0.946	0.948	0.987	0.938	0.000	0.000
$BNB(4;0.93,0.01,0.01)$	1.930	1.010	0.143	0.793	0.793	0.809	0.851	0.853	0.860	0.853
$BNB(4;0.97,0.01,0.01)$	1.970	1.010	0.141	0.815	0.815	0.802	0.912	0.872	0.880	0.864
$BNB(2;0.97,0.97,0.01)$	1.970	1.970	0.493	0.938	0.908	0.891	0.941	0.895	0.629	0.987
$BNB(4;0.98,0.01,0.01)$	1.980	1.010	0.141	0.832	0.830	0.846	0.925	0.876	0.889	0.873
$BNB(4,0.99,0.01,0.01)$	1.990	1.010	0.140	0.823	0.812	0.817	0.880	0.878	0.881	0.878
$BPP(0.40;(0.2,0.2,0.1);(1.0,0.9,0.1))$	1.226	1.190	0.413	0.956	0.930	0.950	0.989	0.803	0.000	0.000
$BPP(0.40;(0.2,0.3,0.1);(0.9,0.9,0.1))$	1.190	1.131	0.361	0.932	0.895	0.913	0.927	0.747	0.000	0.000
$BPP(0.40;(0.8,0.8,0.1);(0.9,1.0,0.4))$	1.003	1.010	0.322	0.867	0.821	0.834	0.864	0.617	0.000	0.000
$BPP(0.45;(0.8,0.8,0.1);(0.9,0.9,0.2))$	1.003	1.003	0.186	0.873	0.811	0.821	0.898	0.614	0.000	0.000
$BPP(0.7;(0.8,0.8,0.1);(0.9,1.1,0.3))$	1.003	1.021	0.208	0.864	0.809	0.893	0.941	0.600	0.000	0.000
$BNTA(0.15;0.01,0.01,0.97)$	1.990	1.990	0.995	0.800	0.802	0.835	0.898	0.615	0.003	0.682
$BNTA(0.42;0.01,0.01,0.98)$	1.990	1.990	0.995	0.908	0.896	0.907	0.949	0.665	0.003	0.849
$BNTA(0.50;0.01,0.01,0.98)$	1.990	1.990	0.995	0.925	0.919	0.831	0.921	0.684	0.000	0.888
$BNTA(0.70;0.01,0.01,0.98)$	1.990	1.990	0.995	0.937	0.919	0.777	0.923	0.730	0.001	0.899
$BNTA(0.75;0.01,0.01,0.98)$	1.990	1.990	0.995	0.932	0.922	0.796	0.978	0.717	0.001	0.910
$BLS(0.01,0.01,0.07)$	0.156	0.156	0.197	0.876	0.930	0.902	1.000	0.800	0.000	0.000
$BLS(0.01,0.01,0.25)$	0.224	0.224	0.829	0.809	0.895	0.916	0.982	0.749	0.015	0.086
$BLS(0.26,0.01,0.04)$	0.263	0.877	0.054	0.690	0.779	0.863	1.000	0.868	0.001	0.001
$BLS{(3d/7,2d/7,2d/7)}^{\ast }$	1.000	1.000	0.447	0.762	0.876	0.872	0.930	0.198	0.159	0.144
$BLS{(3d/4,d/8,d/8)}^{\ast }$	1.000	1.000	0.267	0.942	1.000	0.981	0.909	0.249	0.205	0.191
${}_{}{}^{\ast }d=1-\exp \left(-1\right)\approx 0.63212$.

Table 4. Table IV: Average CPU time (in seconds).

	$n=30$	$n=50$	$n=70$
$R_{n,(0,0)}$	40,804.73	43,974.45	49,328.93
$S_{n,(0,0)}$	3,040.57	7,375.56	14,502.74
$W_n$	1,452.07	1,807.28	2,142.86
$T_{n,(0,0)}$	252.31	518.03	723.42

Table 5. Table V: Results for the real data set ( n = 80 , 994 ) 𝑛 80 994 (n=80,994) .

	Claims
$(a_1,a_2)$	$R_{n,(a_1,a_2)}$	$S_{n,(a_1,a_2)}$	$T_{n,(a_1,a_2)}$	$W_n$
$(0,0)$	0.001	0.001	0.000	0.020
$(1,0)$	0.005	0.005	0.008
$(0,1)$	0.003	0.006	0.010
${\hat{\theta }}_n$	(0.06702119, 0.08841783, 0.01394778)

Table 6. Table VI: Simulation results for the probability of type I error, θ 1 = θ 2 = θ 3 = 1 subscript 𝜃 1 subscript 𝜃 2 subscript 𝜃 3 1 \theta_{1}=\theta_{2}=\theta_{3}=1 .

		$n=30$			$n=50$			$n=70$
${\theta }_4={\rho }_0$ (∗)	Test	f05	f10	$KS$	f05	f10	$KS$	f05	f10	$KS$
0.25	$R_{3,n,(0,0,0)}$	0.035	0.083	0.329116	0.044	0.089	0.129364	0.047	0.092	0.934732
	$S_{3,n,(0,0,0)}$	0.040	0.087	0.508494	0.043	0.091	0.818621	0.048	0.093	0.863178
	$T_{3,n,(0,0,0)}$	0.041	0.088	0.111356	0.043	0.090	0.508494	0.048	0.096	0.718379
	$R_{3,n,(1,0,0)}$	0.039	0.089	0.329116	0.042	0.091	0.329116	0.045	0.094	0.818621
	$S_{3,n,(1,0,0)}$	0.038	0.090	0.902243	0.043	0.090	0.508494	0.043	0.090	0.960002
	$T_{3,n,(1,0,0)}$	0.040	0.088	0.257432	0.041	0.090	0.329116	0.046	0.094	0.369615
	$W_{3,n}$	0.032	0.075	0.129364	0.038	0.082	0.129364	0.042	0.090	0.612128
0.75	$R_{3,n,(0,0,0)}$	0.038	0.82	0.049545	0.041	0.091	0.197933	0.042	0.093	0.559560
	$S_{3,n,(0,0,0)}$	0.040	0.090	0.069329	0.042	0.092	0.257432	0.047	0.094	0.413150
	$T_{3,n,(0,0,0)}$	0.040	0.088	0.149677	0.042	0.085	0.257432	0.051	0.101	0.413150
	$R_{3,n,(1,0,0)}$	0.041	0.081	0.129364	0.044	0.087	0.172476	0.046	0.102	0.863178
	$S_{3,n,(1,0,0)}$	0.040	0.088	0.129364	0.045	0.090	0.197933	0.046	0.095	0.291736
	$T_{3,n,(1,0,0)}$	0.041	0.090	0.111356	0.043	0.094	0.413150	0.051	0.112	0.459543
	$W_{3,n}$	0.032	0.081	0.013476	0.042	0.087	0.111356	0.043	0.091	0.129364
${}_{}{}^{(\ast )}\rho _{12}^{}={\rho }_{13}={\rho }_{23}={\rho }_0$

Table 7. Table VII: Simulation results for the probability of type I error, θ 1 = θ 2 = θ 3 = 2 subscript 𝜃 1 subscript 𝜃 2 subscript 𝜃 3 2 \theta_{1}=\theta_{2}=\theta_{3}=2 .

		$n=30$			$n=50$			$n=70$
${\theta }_4=2{\rho }_0$ (∗)	Test	f05	f10	$KS$	f05	f10	$KS$	f05	f10	$KS$
0.5	$R_{3,n,(0,0,0)}$	0.039	0.085	0.111356	0.042	0.091	0.257432	0.043	0.092	0.863178
	$S_{3,n,(0,0,0)}$	0.040	0.082	0.508494	0.042	0.088	0.818621	0.046	0.093	0.934732
	$T_{3,n,(0,0,0)}$	0.041	0.085	0.329116	0.042	0.091	0.508494	0.045	0.094	0.718379
	$R_{3,n,(1,0,0)}$	0.038	0.087	0.329116	0.041	0.093	0.329116	0.048	0.094	0.818621
	$S_{3,n,(1,0,0)}$	0.039	0.089	0.902243	0.043	0.092	0.508494	0.043	0.092	0.960002
	$T_{3,n,(1,0,0)}$	0.040	0.090	0.129364	0.042	0.090	0.329116	0.045	0.095	0.369615
	$W_{3,n}$	0.032	0.086	0.111356	0.043	0.087	0.369615	0.048	0.092	0.612128
1.00	$R_{3,n,(0,0,0)}$	0.038	0.091	0.129364	0.044	0.092	0.197933	0.046	0.094	0.559560
	$S_{3,n,(0,0,0)}$	0.042	0.089	0.049545	0.043	0.093	0.257432	0.046	0.095	0.413150
	$T_{3,n,(0,0,0)}$	0.040	0.090	0.149677	0.043	0.091	0.257432	0.053	0.101	0.413150
	$R_{3,n,(1,0,0)}$	0.051	0.101	0.129364	0.044	0.097	0.172476	0.047	0.109	0.863178
	$S_{3,n,(1,0,0)}$	0.040	0.090	0.069329	0.056	0.109	0.197933	0.052	0.102	0.291736
	$T_{3,n,(1,0,0)}$	0.040	0.089	0.129364	0.044	0.094	0.413150	0.051	0.103	0.459543
	$W_{3,n}$	0.038	0.081	0.111356	0.042	0.087	0.111356	0.047	0.092	0.413150
${}_{}{}^{(\ast )}\rho _{12}^{}={\rho }_{13}={\rho }_{23}={\rho }_0$

Table 8. Table VIII: Simulation results for the probability of type I error, θ 1 = θ 2 = 0.2 subscript 𝜃 1 subscript 𝜃 2 0.2 \theta_{1}=\theta_{2}=0.2 .

		$n=30$			$n=50$			$n=70$
${\theta }_3,{\theta }_4;\rho$	Test	f05	f10	$KS$	f05	f10	$KS$	f05	f10	$KS$
$0.8,0.1;(0.5,0.25,0.25)$	$R_{3,n,(0,0,0)}$	0.046	0.094	0.069329	0.058	0.107	0.559560	0.052	0.101	0.226206
	$S_{3,n,(0,0,0)}$	0.045	0.087	0.291736	0.050	0.104	0.718379	0.051	0.093	0.863178
	$T_{3,n,(0,0,0)}$	0.046	0.090	0.718379	0.056	0.098	0.902243	0.053	0.105	0.989545
	$R_{3,n,(1,0,0)}$	0.048	0.091	0.149677	0.054	0.108	0.013476	0.052	0.107	0.960002
	$S_{3,n,(1,0,0)}$	0.043	0.092	0.226206	0.055	0.111	0.413150	0.053	0.097	0.508494
	$T_{3,n,(1,0,0)}$	0.040	0.091	0.369615	0.053	0.109	0.559560	0.052	0.105	0.559560
	$W_{3,n}$	0.038	0.086	0.041633	0.043	0.087	0.111356	0.047	0.092	0.197933

Table 9. Table IX: Simulation results for the probability of type I error, θ 2 = θ 3 = 0.8 subscript 𝜃 2 subscript 𝜃 3 0.8 \theta_{2}=\theta_{3}=0.8 .

		$n=30$			$n=50$			$n=70$
${\theta }_1,{\theta }_4;\rho$	Test	f05	f10	$KS$	f05	f10	$KS$	f05	f10	$KS$
$1.8,0.6;(0.5,0.5,0.75)$	$R_{3,n,(0,0,0)}$	0.040	0.089	0.069329	0.045	0.090	0.226206	0.045	0.091	0.559560
	$S_{3,n,(0,0,0)}$	0.042	0.086	0.291736	0.045	0.094	0.718379	0.047	0.093	0.902243
	$T_{3,n,(0,0,0)}$	0.044	0.089	0.197933	0.046	0.092	0.863178	0.047	0.101	0.989545
	$R_{3,n,(1,0,0)}$	0.046	0.089	0.013476	0.045	0.094	0.508494	0.046	0.102	0.960002
	$S_{3,n,(1,0,0)}$	0.044	0.088	0.226206	0.045	0.091	0.413150	0.045	0.098	0.559560
	$T_{3,n,(1,0,0)}$	0.041	0.089	0.369615	0.045	0.093	0.508494	0.046	0.101	0.559560
	$W_{3,n}$	0.039	0.087	0.149677	0.043	0.087	0.197933	0.046	0.095	0.291736

Table 10. Table X: Simulation results for the probability of type I error, θ 1 = 8.7 subscript 𝜃 1 8.7 \theta_{1}=8.7 , θ 2 = 8.8 subscript 𝜃 2 8.8 \theta_{2}=8.8 , θ 3 = 8.9 subscript 𝜃 3 8.9 \theta_{3}=8.9 .

		$n=30$			$n=50$			$n=70$
${\theta }_4;\rho$	Test	f05	f10	$KS$	f05	f10	$KS$	f05	f10	$KS$
$2.2;(0.251,0.250,0.249)$	$R_{3,n,(0,0,0)}$	0.036	0.083	0.149677	0.043	0.089	0.226206	0.041	0.089	0.863178
	$S_{3,n,(0,0,0)}$	0.038	0.085	0.291736	0.040	0.089	0.508494	0.043	0.092	0.718379
	$T_{3,n,(0,0,0)}$	0.040	0.086	0.226206	0.042	0.089	0.559560	0.045	0.105	0.902243
	$R_{3,n,(1,0,0)}$	0.036	0.085	0.197933	0.041	0.089	0.413150	0.045	0.104	0.863178
	$S_{3,n,(1,0,0)}$	0.037	0.086	0.369615	0.040	0.090	0.508494	0.044	0.094	0.559560
	$T_{3,n,(1,0,0)}$	0.034	0.084	0.197933	0.042	0.090	0.508494	0.045	0.105	0.559560
	$W_{3,n}$	0.035	0.083	0.069329	0.040	0.087	0.197933	0.046	0.093	0.291736

Table 11. Table XI: Simulation results for the probability of type I error, θ 1 = 9.7 subscript 𝜃 1 9.7 \theta_{1}=9.7 , θ 2 = 9.6 subscript 𝜃 2 9.6 \theta_{2}=9.6 , θ 3 = 9.5 subscript 𝜃 3 9.5 \theta_{3}=9.5 .

		$n=30$			$n=50$			$n=70$
${\theta }_4;\rho$	Test	f05	f10	$KS$	f05	f10	$KS$	f05	f10	$KS$
$2.4;(0.249,0.250,0.251)$	$R_{3,n,(0,0,0)}$	0.035	0.081	0.413150	0.041	0.086	0.508494	0.041	0.088	0.863178
	$S_{3,n,(0,0,0)}$	0.036	0.086	0.291736	0.039	0.087	0.369615	0.044	0.091	0.718379
	$T_{3,n,(0,0,0)}$	0.039	0.086	0.226206	0.040	0.088	0.559560	0.044	0.095	0.508494
	$R_{3,n,(1,0,0)}$	0.037	0.085	0.197933	0.040	0.088	0.149677	0.046	0.094	0.559560
	$S_{3,n,(1,0,0)}$	0.036	0.085	0.226206	0.041	0.089	0.291736	0.045	0.093	0.369615
	$T_{3,n,(1,0,0)}$	0.035	0.084	0.069329	0.040	0.090	0.559560	0.046	0.095	0.508494
	$W_{3,n}$	0.034	0.082	0.197933	0.040	0.088	0.197933	0.045	0.094	0.863178

Table 12. Table XII: Average CPU time (in seconds).

	$n=30$	$n=50$	$n=70$
$R_{3,n,(0,0,0)}$	47,338.59	72,539.85	107,952.20
$S_{3,n,(0,0,0)}$	3,234.27	9,357.65	18,252.43
$T_{3,n,(0,0,0)}$	4,486.30	12,617.74	24,763.87
$W_{3,n}$	1,833.80	2,174.36	2,323.33

Table 13. Table XIII: Results for the real data set ( n = 162 , 019 𝑛 162 019 n=162,019 ).

	Claims
$(a_1,a_2,a_3)$	$R_{3,n,(a_1,a_2,a_3)}$	$S_{3,n,(a_1,a_2,a_3)}$	$T_{3,n,(a_1,a_2,a_3)}$	$W_{3,n}$
$(0,0,0)$	0.001	0.001	0.001	0.0001
$(1,0,0)$	0.002	0.002	0.001
${\hat{\theta }}_n$	(0.249051, 0.03508231, 0.201069, 0.03508218 )

Equations

$$X_1=Y_1+Y_3\text{ and }{X}_2=Y_2+Y_3,$$

$$D^{k_1\cdots k_m}h(u)=\frac{{\partial }^k}{\partial u_1^{k_1}\cdots \partial u_m^{k_m}}h(u),$$

$$H_0:(X_1,X_2)\sim BP({\theta }_1,{\theta }_2,{\theta }_3),\text{for some}({\theta }_1,{\theta }_2,{\theta }_3)\in \Theta ,$$

$$H_1:(X_1,X_2)\nsim BP({\theta }_1,{\theta }_2,{\theta }_3),\forall ({\theta }_1,{\theta }_2,{\theta }_3)\in \Theta ,$$

$$g(u;\theta)=E_{\theta}(u_{1}^{X_{1}}u_{2}^{X_{2}})=\exp\bigl{\{}\theta_{1}(u_{1}-1)+\theta_{2}(u_{2}-1)+\theta_{3}(u_{1}-1)(u_{2}-1)\bigr{\}}.$$

$$T=\frac{n}{2}\frac{{\bar{X}}_2^2{\left(S_{X_1}^2-{\bar{X}}_1\right)}^2-2S_{X_1{X}_2}^2\left(S_{X_1}^2-{\bar{X}}_1\right)\left(S_{X_2}^2-{\bar{X}}_2\right)+{\bar{X}}_1^2{\left(S_{X_2}^2-{\bar{X}}_2\right)}^2}{{\bar{X}}_1^2{\bar{X}}_2^2-S_{X_1{X}_2}^4},T\ge {\chi }_{2,\alpha }^2.$$

$$I_B=\frac{n({\bar{X}}_2{S}_{X_1}^2-2S_{X_1{X}_2}^2+{\bar{X}}_1{S}_{X_2}^2)}{{\bar{X}}_1{\bar{X}}_2-S_{X_1{X}_2}^2},I_B\ge {\chi }_{2n-3,\alpha }^2.$$

$$N{I}_B=\frac{n}{1-r^2}\left(\frac{S_{X_1}^2}{{\bar{X}}_1}-2r^2\sqrt{\frac{S_{X_1}^2{S}_{X_2}^2}{{\bar{X}}_1{\bar{X}}_2}}+\frac{S_{X_2}^2}{{\bar{X}}_2}\right).$$

$$R_{n,w}={\int }_0^1{\int }_0^1{G}_n^2(u;{\hat{\theta }}_n)w(u)du,$$

$$S_{n,w}=n{\int }_0^1{\int }_0^1\left\{B_{1n}^2(u;{\hat{\theta }}_n)+B_{2n}^2(u;{\hat{\theta }}_n)\right\}w(u)du,$$

$$\begin{array}[]{ll}B_{1n}(u;\widehat{\theta}_{n})&=\displaystyle\frac{\partial}{\partial u_{1}}g_{n}(u_{1},u_{2})-\left\{\widehat{\theta}_{1n}+\widehat{\theta}_{3n}(u_{2}-1)\right\}g_{n}(u_{1},u_{2}),\\[8.5359pt] B_{2n}(u;\widehat{\theta}_{n})&=\displaystyle\frac{\partial}{\partial u_{2}}g_{n}(u_{1},u_{2})-\left\{\widehat{\theta}_{2n}+\widehat{\theta}_{3n}(u_{1}-1)\right\}g_{n}(u_{1},u_{2}),\end{array}$$

$$W_n=\sum _{r_1\ge 0}\sum _{r_2\ge 0}\{b_1^2(r_1,r_2;{\hat{\theta }}_n)+b_2^2(r_1,r_2;{\hat{\theta }}_n)\}=\sum _{r_1,r_2=0}^M\{b_1^2(r_1,r_2;{\hat{\theta }}_n)+b_2^2(r_1,r_2;{\hat{\theta }}_n)\},$$

$$b_1(r_1,r_2;{\hat{\theta }}_n)$$

$$b_2(r_1,r_2;{\hat{\theta }}_n)$$

$$X_1=Y_1+Y_{m+1},X_2=Y_2+Y_{m+1},\dots ,X_m=Y_m+Y_{m+1},$$

$$g(u;\theta )=\exp \left\{\sum _{i=1}^m{\theta }_i\left(u_i-1\right)+{\theta }_{m+1}\left(\prod _{i=1}^m{u}_i-\sum _{i=1}^m{u}_i+m-1\right)\right\},\forall u\in {\mathbb{R}}^m.$$

$$g_n(u)=\frac{1}{n}\sum _{i=1}^n{u}_1^{X_{i1}}\cdots u_m^{X_{im}},u=(u_1,\dots ,u_m).$$

$$H_{0m}:(X_{1},X_{2},\ldots,X_{m})\mbox{ has a $d$-variate Poisson distribution}.$$

$$S_{m,n,w}=n{\int }_{[0,1{]}^m}\left\{B_{1n}^2(u;{\hat{\theta }}_n)+\cdots +B_{mn}^2(u;{\hat{\theta }}_n)\right\}w(u)du,$$

$$B_{in}(u;{\hat{\theta }}_n)=\frac{\partial }{\partial u_i}{g}_n(u)-\left\{{\hat{\theta }}_{i,n}+{\hat{\theta }}_{m+1,n}\left(\prod _{j\ne i}{u}_j-1\right)\right\}{g}_n(u),1\le i\le m,g_n(u)\text{ as in (refefgp-multivariante)}.$$

$$W_{m,n}=\sum _{r_1,r_2,\dots ,r_m\ge 0}\left\{\sum _{j=1}^m{b}_j^2(r_1,r_2,\dots ,r_m;{\hat{\theta }}_n)\right\}=\sum _{r_1,r_2,\dots ,r_m=0}^M\left\{\sum _{j=1}^m{b}_j^2(r_1,r_2,\dots ,r_m;{\hat{\theta }}_n)\right\},$$

$$b_j(r_1,\dots ,r_m;{\hat{\theta }}_n)$$

$$-{\hat{\theta }}_{m+1,n}{p}_n(r_1-1,\dots ,r_{j-1}-1,r_j,r_{j+1}-1,\dots ,r_m-1),1\le j\le m,$$

$$\left.\begin{array}[]{rcl}D_{1}(u;\theta)&=&\displaystyle\frac{\partial g(u_{1},1)}{\partial u_{1}}-\theta_{1}\,g(u_{1},1)=0,\\[8.5359pt] D_{2}(u;\theta)&=&\displaystyle\frac{\partial g(1,u_{2})}{\partial u_{2}}-\theta_{2}\,g(1,u_{2})=0,\\[8.5359pt] D_{3}(u;\theta)&=&\displaystyle\frac{\partial^{2}g(u_{1},u_{2})}{\partial u_{1}\,\partial u_{2}}-f(u_{1},u_{2};\theta)\,g(u_{1},u_{2})=0,\end{array}\right\}$$

$$\begin{array}[]{rcl}D_{1n}\bigl{(}u;\widehat{\theta}_{n}\bigr{)}&=&\displaystyle\frac{\partial g_{n}(u_{1},1)}{\partial u_{1}}-\widehat{\theta}_{1n}\ g_{n}(u_{1},1),\\[8.5359pt] D_{2n}\bigl{(}u;\widehat{\theta}_{n}\bigr{)}&=&\displaystyle\frac{\partial g_{n}(1,u_{2})}{\partial u_{2}}-\widehat{\theta}_{2n}\ g_{n}(1,u_{2}),\\[8.5359pt] D_{3n}\bigl{(}u;\widehat{\theta}_{n}\bigr{)}&=&\displaystyle\frac{\partial^{2}g_{n}(u_{1},u_{2})}{\partial u_{1}\,\partial u_{2}}-f(u_{1},u_{2};\widehat{\theta}_{n})\,g_{n}(u_{1},u_{2}),\end{array}$$

$$g_n(u_1,u_2)=\frac{1}{n}\sum _{i=1}^n{u}_1^{X_{i1}}{u}_2^{X_{i2}}.$$

$$T_{n,w}=n\int_{0}^{1}\int_{0}^{1}\left\{D^{2}_{1n}\bigl{(}u;\widehat{\theta}_{n}\bigr{)}+D^{2}_{2n}\bigl{(}u;\widehat{\theta}_{n}\bigr{)}+D^{2}_{3n}\bigl{(}u;\widehat{\theta}_{n}\bigr{)}\right\}w(u)\ du,$$

$$\frac{T_{n,w}}{n}\stackrel{a.s.}{⟶}{\int }_0^1{\int }_0^1\left\{D_1^2(u;\theta )+D_2^2(u;\theta )+D_3^2(u;\theta )\right\}w(u)du=\eta (g;\theta )\ge 0.$$

$$\mathcal{H}=\{\phi :[0,1{]}^2\to \mathbb{R},\text{with}\Vert \phi {\Vert }_{{}_{\mathcal{H}}}^2={\int }_0^1{\int }_0^1{\phi }^2(u)w(u)du<\infty \}.$$

$$Z_{kn}(u)=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}R_{k}\bigl{(}\boldsymbol{X}_{i};\widehat{\theta}_{n};u\bigr{)},\ \ k=1,2,3,$$

$$\displaystyle R_{1}\bigl{(}\boldsymbol{X}_{i};\widehat{\theta}_{n};u\bigr{)}$$

Full text

GOODNESS-OF-FIT TESTS FOR THE BIVARIATE POISSON DISTRIBUTION

F. Novoa-Muñoz

Departamento de Estadística

Universidad del Bío-Bío

fnovoa@ubiobio.cl

Key Words: bivariate Poisson distribution; goodness-of-fit; empirical probability generating function; consistency against fixed alternatives; bootstrap distribution estimator.

ABSTRACT

The bivariate Poisson distribution is commonly used to model bivariate count data. In this paper we study a goodness-of-fit test for this distribution. We also provide a review of the existing tests for the bivariate Poisson distribution, and its multivariate extension. The proposed test is consistent against any fixed alternative. It is also able to detect local alternatives converging to the null at the rate $n^{-\frac{1}{2}}$. The bootstrap can be employed to consistently estimate the null distribution of the test statistic. Through a simulation study we investigated the goodness of the bootstrap approximation and the power for finite sample sizes.

1 Introduction

The univariate Poisson distribution (UPD) has helped to model many real life situations. For a survey of statistical issues, problems and applications associated with the UPD the reader is referred to the text of Haight (1967) and Johnson and Kotz (1969). For the other hand, the bivariate Poisson distribution (BPD) is appropriate for modelling paired count data exhibiting positive correlation.

Several definitions for the BPD have been given (see, e.g. Kocherlakota and Kocherlakota, 1992). In this paper we will work with the following one, because it has received the most attention in the statistical literature (see, e.g. Holgate, 1964; Johnson, Kotz and Balakrishnan, 1997). Let

[TABLE]

where $Y_{1},Y_{2}$ and $Y_{3}$ are independent Poisson random variables with means $\theta^{\prime}_{1}=\theta_{1}-\theta_{3}>0$, $\theta^{\prime}_{2}=\theta_{2}-\theta_{3}>0$ and $\theta_{3}>0$, respectively. The joint distribution of the vector $(X_{1},X_{2})$ is called BPD with parameter $\theta=(\theta_{1},\theta_{2},\theta_{3})$, $(X_{1},X_{2})\sim BP(\theta)$ for short.

In the statistical literature on goodness-of-fit (gof) tests for the BPD, which is rather sparse in comparison with the univariate case, we found the following: the tests given by Crockett (1979), Loukas and Kemp (1986), Rayner and Best (1995) -these three tests are not consistent against each fixed alternative- and, more recently, the tests in Novoa-Muñoz and Jiménez-Gamero (2014), and Novoa-Muñoz and Jiménez-Gamero (2016) (hereafter abbreviated to NJ (2014) and NJ (2016), respectively).

The tests in NJ (2014) and NJ (2016) are consistent against each fixed alternative. The results in Janssen (2000) assert that the global power function of any nonparametric test is flat on balls of alternatives except for alternatives coming from a finite-dimensional subspace. Therefore, it is interesting to propose new gof tests able to detect different sets of alternatives.

The present work proposes a new consistent gof test for the BPD. To derive it we first show that the probability generating function (pgf) of the BPD is the only pgf that satisfies a certain system of partial differential equations. Therefore, under the null hypothesis, the empirical probability generating function (epgf), which is a consistent estimator of the pgf (see, e.g. NJ, 2014), should approximately satisfy such system. The proposed test statistic can be seen as a bivariate extension of the one in Baringhaus and Henze (1992) designed for testing gof to the univariate Poisson distribution.

The asymptotic behavior of the proposed test under alternatives is shared with the ones in NJ (2014) or NJ (2016). An advantage of the test proposed in this paper over those in NJ (2014) and NJ (2016) is its speed for the delivery of results.

In order to consistently approximate the null distribution of the test statistic, we propose to use a parametric bootstrap estimator. The finite-sample size performance of the test is numerically evaluated through a simulation study. The power of the test is compared with the tests mentioned above. There is no test yielding the highest power against each considered alternative, as expected from the results in Janssen (2000). In most cases, the power of the proposed test is quite close to the highest one; in other cases, the proposed test is the most powerful. In addition, from a computational point of view, the test proposed in this paper is more efficient than its competitors.

Next we show the notation used in this work: all vectors are row vectors and $v^{\top}$ is the transpose of the row vector $v$; for any vector $v,\,v_{k}$ denotes its $k$th coordinate, and $\|v\|$ its Euclidean norm. We put $\,\mathbb{N}_{0}=\{0,1,2,3,\ldots\}$ and write $I_{A}\,$ for the indicator function of the set $A$; $P_{\theta}$ denotes the probability law of the BPD with parameter $\theta$; $P$ denotes the probability law of the data; $E_{\theta}$ denotes expectation regarding the probability function $P_{\theta}$; $E$ denotes expectation with respect to the true probability function of the data; $P_{*}$ denote the probability law, given the data; all limits in this work are taken as n\rightarrow\infty;\mathop{\longrightarrow}\limits^{L}\ denotes convergence in distribution; \mathop{\longrightarrow}\limits^{a.s.}\ denotes almost sure (a.s.) convergence. For any function $h:S\subset\mathbb{R}^{m}\to\mathbb{R}$, for some fixed $m\in\mathbb{N}$, we will denote

[TABLE]

for each choice of nonnegative integers $k_{1},\ldots,k_{m}$ such that $k=k_{1}+\cdots+k_{m}$.

2 Review of the existing tests for the BPD, and their multivariate extension

2.1 Tests for the BPD

Let $\boldsymbol{X}_{1}=(X_{11},X_{12}),\boldsymbol{X}_{2}=(X_{21},X_{22}),\ldots,\boldsymbol{X}_{n}=(X_{n1},X_{n2})$ be independent identically distributed (iid) from a random vector $\boldsymbol{X}=(X_{1},X_{2})$ taking values in $\mathbb{N}_{0}^{2}$. Based on the sample $\boldsymbol{X}_{1},\boldsymbol{X}_{2},\ldots,\boldsymbol{X}_{n}$, the objective is to test the hypothesis

[TABLE]

against the alternative

[TABLE]

where $\Theta=\left\{\left(\theta_{1},\theta_{2},\theta_{3}\right)\in\mathbb{R}^{3}\ :\ \theta_{1}>\theta_{3},\ \theta_{2}>\theta_{3},\ \theta_{3}>0\right\}$. From NJ (2014) the distribution of $\boldsymbol{X}=(X_{1},X_{2})$ is determined by its pgf $g(u)=E\left(u_{1}^{X_{1}}u_{2}^{X_{2}}\right),\,u=(u_{1},u_{2})\in[0,1]^{2}$, and the joint pgf of a random vector $\boldsymbol{X}\sim BP(\theta)$ is

[TABLE]

The empirical counterpart of pgf is epgf of the data given by $g_{n}(u)=\frac{1}{n}\sum_{i=1}^{n}u_{1}^{X_{i1}}u_{2}^{X_{i2}}$.

Next we will briefly expose three non-consistent tests that we found in the statistical literature, where $\bar{X}_{1},\bar{X}_{2},S^{2}_{X_{1}}$ and $S^{2}_{X_{2}}$ are the sample means and variances, respectively, $S^{2}_{X_{1}X_{2}}$ is the sample covariance, $r$ is the sample correlation coefficient and $\chi^{2}_{k,\alpha}$, for $0<\alpha<1$ and $k\in\mathbb{N}$, denotes the upper $\alpha$-percentile of the $\chi^{2}$ distribution with $k$ degrees of freedom.

2.1.1 Crockett test $\boldsymbol{T}$

The statistic $T$ (say) proposed by Crockett (1979) is based on a quadratic form in $Z_{X_{1}}=S^{2}_{X_{1}}-\bar{X}_{1}$ and $Z_{X_{2}}=S^{2}_{X_{2}}-\bar{X}_{2}$. He shows that, under $H_{0}$, $T=ZV^{-1}Z^{\top}\ \mathop{\longrightarrow}\limits^{L}\ Y\sim\chi^{2}_{2}$, where $V$ denotes the matrix of variances and covariances of $Z=(Z_{X_{1}},Z_{X_{2}})$. Thus, the statistic and its critical region are given by

[TABLE]

2.1.2 Test $\boldsymbol{I_{B}}$ of Loukas and Kemp

Loukas and Kemp (1986) developed a test based on what they called the bivariate dispersion index, $I_{B}=\frac{1}{1-\rho^{2}}\sum_{i=1}^{n}\left(W_{i1}^{2}-2\rho W_{i1}W_{i2}+W_{i2}^{2}\right),$ where $W_{ik}=\frac{X_{ik}-\theta_{k}}{\sqrt{\theta_{k}}},k=1,2,\,i=1,2,\ldots,n$ and $\ \rho=\frac{\theta_{3}}{\sqrt{\theta_{1}\theta_{2}}}$. If $\theta_{1},\theta_{2}$ and $\theta_{3}$ are known, these authors show that $I_{B}$ is distributed approximately as a variable $\chi^{2}_{2n}$. If $\theta_{1},\theta_{2}$ and $\theta_{3}$ are unknown, the statistic and its critical region are given by

[TABLE]

2.1.3 Test $\boldsymbol{NI_{B}}$ of Rayner and Best

Rayner and Best (1995) expressed the statistics of Loukas and Kemp (1986) as $I_{B}=\frac{n}{1-\widehat{\rho}^{\,2}}\!\left(\frac{S_{X_{1}}^{2}}{\bar{X}_{1}}-2\,\frac{S_{X_{1}X_{2}}^{2}}{\bar{X}_{1}\bar{X}_{2}}+\frac{S_{X_{2}}^{2}}{\bar{X}_{2}}\right),$ where $\widehat{\rho}=\frac{S_{X_{1}X_{2}}}{\sqrt{\bar{X}_{1}\bar{X}_{2}}}$ is an estimator of $\rho$. If $\widehat{\rho}^{\,2}>\frac{1}{2}\left(\frac{S_{X_{1}}^{2}}{\bar{X}_{1}}+\frac{S_{X_{2}}^{2}}{\bar{X}_{2}}\right)$, then $I_{B}<0$, and its distribution is not well approximated by a $\chi^{2}$. They proposed the statistic

[TABLE]

Under $H_{0}$, $NI_{B}$ is approximately distributed as $\chi^{2}_{2n-3}$ if $n$ is large. Therefore, $H_{0}$ is rejected if $NI_{B}\geq\chi^{2}_{2n-3,\alpha}.$

Note that the statistical tests $T,I_{B}$ and $NI_{B}$ are not consistent, because they are based on the moments, specifically based on the fact that the first two population moments are equal. In contrast, the tests presented below are consistent.

2.1.4 Test $\boldsymbol{R_{n,w}}$ of Novoa-Muñoz and Jiménez-Gamero (2014)

From NJ (2014) the distribution of $\boldsymbol{X}$ is uniquely determined by its pgf, $g(u)$, $u\in[0,1]^{2}$, a reasonable test for testing $H_{0}$ should reject the null hypothesis for large values of

[TABLE]

where $G_{n}(u;\vartheta)=\sqrt{n}\{g_{n}(u)-g(u;\vartheta)\}$, $w(u)=u_{1}^{a_{1}}u_{2}^{a_{2}}$ is a measurable weight function $\forall u\in[0,1]^{2}$, $a_{1},a_{2}\in(-1,\infty)$, and $\widehat{\theta}_{n}=\widehat{\theta}_{n}(\boldsymbol{X}_{1},\ldots,\boldsymbol{X}_{n})=(\widehat{\theta}_{1n},\widehat{\theta}_{2n},\widehat{\theta}_{3n})$ is a consistent estimator of $\theta$.

2.1.5 Test $\boldsymbol{S_{n,w}}$ of Novoa-Muñoz and Jiménez-Gamero (2014)

Since the pgf $g(x)$ of the univariate Poisson distribution, with parameter $\lambda$, is the only pgf satisfying the differential equation $g^{\prime}(x)=\lambda g(x)$, Baringhaus and Henze (1992) proposed a test statistic which is based on an empirical counterpart of this equation. With the aim of extending this result to the bivariate case, NJ (2014) proposed to reject $H_{0}$ for large values of

[TABLE]

where $w(u)=u_{1}^{a_{1}}u_{2}^{a_{2}}$, $a_{1},a_{2}\in(-1,\infty)$, $\widehat{\theta}_{n}=(\widehat{\theta}_{1n},\widehat{\theta}_{2n},\widehat{\theta}_{3n})$ is a consistent estimator of $\theta$, and

[TABLE]

should be close to 0 when $H_{0}$ is true. These functions are the empirical counterpart of the system of partial differential equations of Proposition 2 in NJ (2014).

2.1.6 Test $\boldsymbol{W_{n}}$ of Novoa-Muñoz and Jiménez-Gamero (2016)

When $H_{0}$ is true, NJ (2016) presented another interpretation of the fact that $S_{n,w}=n\int_{0}^{1}\int_{0}^{1}\{B^{2}_{1n}(u;\widehat{\theta}_{n})+B^{2}_{2n}(u;\widehat{\theta}_{n})\}w(u)\ du\approx 0$. Reasoning as Nakamura and Pérez-Abreu (1993) for the univariate case and noting that $B_{kn}(u;\widehat{\theta}_{n})=\sum_{r_{1}\geq 0}\sum_{r_{2}\geq 0}b_{k}(r_{1},r_{2};\widehat{\theta}_{n})u_{1}^{r_{1}}u_{2}^{r_{2}}$, $k=1,2$. NJ (2016) proposed to reject $H_{0}$ for large values of

[TABLE]

where $M=\max\{X_{(n)1},X_{(n)2}\}$, $X_{(n)k}=\max_{1\leq i\leq n}X_{ik}$, $k=1,2$,

[TABLE]

and $p_{n}(r_{1},r_{2})=\frac{1}{n}\sum_{i=1}^{n}I_{\{X_{i1}=r_{1},X_{i2}=r_{2}\}}$ is the relative frequency of the pair $(r_{1},r_{2})$,

2.2 The general $m-$variate case

For the multivariate case, for each integer $m>2$, let

[TABLE]

where $Y_{1},Y_{2},\ldots,Y_{m+1}$ are independent Poisson random variables with means $\theta^{\prime}_{1}=\theta_{1}-\theta_{m+1}>0,\ldots,\theta^{\prime}_{m}=\theta_{m}-\theta_{m+1}>0$ and $\theta_{m+1}>0$, respectively. The joint distribution of the vector $(X_{1},X_{2},\ldots,X_{m})$ is called a $m$-variate Poisson distribution with parameter $\theta=(\theta_{1},\theta_{2},\ldots,\theta_{m},\theta_{m+1})$ (see Johnson, Kotz and Balakrishnan, 1997). The joint pgf of $(X_{1},X_{2},\ldots,X_{m})$ is

[TABLE]

The empirical counterpart of pgf is epgf of the data given by

[TABLE]

Now, the objective is to test the hypothesis

[TABLE]

The tests proposed by Crockett (1979), Loukas and Kemp (1986), and Rayner and Best (1995) do not have a multivariate extension. However, NJ (2014) and NJ (2016) proposed a natural extension of their tests, which will be presented below.

2.2.1 Test $\boldsymbol{R_{m,n,w}}$ of Novoa-Muñoz and Jiménez-Gamero (2014)

NJ (2014) affirmed that the extension of the test $R_{n,w}$ is direct, it is enough to consider pgf $g(u;\theta)$ as in (2), epgf $g_{n}(u;\widehat{\theta}_{n})$ as in (3) and $w(u)$ is a measurable nonnegative weight function with finite integral over $[0,1]^{m}$.

2.2.2 Test $\boldsymbol{S_{m,n,w}}$ of Novoa-Muñoz and Jiménez-Gamero (2014)

To test $H_{0m}$, NJ (2014) considered the test statistic

[TABLE]

where $w(u)$ is a measurable nonnegative weight function with finite integral over $[0,1]^{m}$, and

[TABLE]

2.2.3 Test $\boldsymbol{W_{m,n}}$ of Novoa-Muñoz and Jiménez-Gamero (2016)

With the aim of extending $W_{n}$ to the multivariate case, NJ (2016) proposed the following statistic for testing $H_{0m}$,

[TABLE]

where $M=\max\{X_{(n)1},X_{(n)2},\ldots,X_{(n)m}\}$, $X_{(n)k}=\max_{1\leq i\leq n}X_{ik}$, $1\leq k\leq m$, and

[TABLE]

and $p_{n}(r_{1},\ldots,r_{m})=\frac{1}{n}\sum_{i=1}^{n}I_{\{X_{i1}=r_{1},\ldots,X_{im}=r_{m}\}}$ is the relative frequency of $(r_{1},\ldots,r_{m})$.

3 A new characterization of the BPD

In order to obtain a new test to test the hypothesis $H_{0}$ against the alternative $H_{1}$ and based on the fact that the distribution of $\boldsymbol{X}=(X_{1},X_{2})$ is determined by its pgf, we give a different characterization for the BPD.

Proposition 1

Let $g(u_{1},u_{2};\theta)\,$ be as defined in (1). Then $g(u_{1},u_{2};\theta)\,$ is the only pgf satisfying the following system of partial differential equations

[TABLE]

where $\,f(u_{1},u_{2};\theta)=\theta_{3}+\{\theta_{2}+\theta_{3}(u_{1}-1)\}\{\theta_{1}+\theta_{3}(u_{2}-1)\}$.

The system of equations (4) has the following nice interpretation: first and the second equation characterize the marginal distributions, i. e., they are equivalent to saying that the marginal distributions are univariate Poisson; the last equation characterizes the dependence structure.

By Proposition 1 in NJ (2014), $g(u)$ and its derivatives can be consistently estimated by the epgf and the derivatives of the epgf, respectively. Thus, if $H_{0}$ is true, then the functions

[TABLE]

should be close to [math], $\forall(u_{1},u_{2})\in[0,1]^{2}$, where $\widehat{\theta}_{n}=(\widehat{\theta}_{1n},\widehat{\theta}_{2n},\widehat{\theta}_{3n})$ is a consistent estimator of $\theta$ and $g_{n}(u_{1},u_{2})$ is the epgf associated with the data, i. e.,

[TABLE]

Thus, to test $H_{0}$ we consider the following test statistic

[TABLE]

where $w(u)$ is a non-negative function on $[0,1]^{2}$.

In order to give a sound justification of $\,T_{n,w}$ as a test statistic for testing $H_{0}$ we next derive its almost sure limit.

Theorem 1

Let $\boldsymbol{X}_{1},\boldsymbol{X}_{2},\ldots,\boldsymbol{X}_{n}$ be iid from $\boldsymbol{X}=(X,Y)\in\mathbb{N}_{0}^{2}$ with pgf $g(u)$ such that $\frac{\partial g(u_{1},1)}{\partial u_{1}}$, $\frac{\partial g(1,u_{2})}{\partial u_{2}}$ and $\frac{\partial^{2}g(u_{1},u_{2})}{\partial u_{1}\,\partial u_{2}}$, exist and are continuous functions on a region containing $[0,1]^{2}$. If $\widehat{\theta}_{n}\mathop{\longrightarrow}\limits^{a.s.}\theta$, for some $\,\theta\in\mathbb{R}^{3}$, then

[TABLE]

Note that if $w>0$ almost everywhere (a.e.) on $[0,1]^{2}$, then $\eta(g;\theta)=0$ if and only if $H_{0}$ is true. Therefore, a reasonable test for testing $H_{0}$ should reject the null hypothesis for large values of $T_{n,w}$. Now, to determine what are large values of $T_{n,w}$, we must calculate its null distribution, or at least an approximation to it. Clearly, the null distribution of $T_{n,w}$ is unknown. A classical way of approximating the null distribution of a test statistic is through its asymptotic null distribution. The next section studies this issue.

4 A bootstrap estimator of the null distribution

In order to derive the asymptotic null distribution of the test statistic $T_{n,w}$ we will assume that the estimator $\widehat{\theta}_{n}$ is asymptotically linear, as expressed in Assumption 1 in NJ (2014) and we will consider the separable Hilbert space

[TABLE]

In this framework, $T_{n,w}$ can be expressed as $T_{n,w}=\|Z_{1n}\|_{{}_{\mathcal{H}}}^{\,2}+\|Z_{2n}\|_{{}_{\mathcal{H}}}^{\,2}+\|Z_{3n}\|_{{}_{\mathcal{H}}}^{\,2}$, with

[TABLE]

where, for $1\leq i\leq n$,

[TABLE]

The next result gives the asymptotic null distribution of $T_{n,w}$.

Theorem 2

Let $\,\boldsymbol{X}_{1},\boldsymbol{X}_{2},\ldots,\boldsymbol{X}_{n}\,$ be iid from $\,\boldsymbol{X}=(X_{1},X_{2})\sim BP(\theta)$. Suppose that Assumption 1 in NJ (2014) holds and that $\widehat{\theta}_{n}\mathop{\longrightarrow}\limits^{a.s.}\theta$. Then

[TABLE]

where $\ P_{\theta}(|r_{n}|>\varepsilon)\to 0$, $\forall\varepsilon>0$,

[TABLE]

[TABLE]

$1\leq i\leq n$, $B(u;\theta)=\!\left(b_{1}(u;\theta),b_{2}(u;\theta),b_{3}(u;\theta)\right),$ where $\,b_{1}(u;\theta)=\theta_{2}+\theta_{3}(u_{1}-1)$, $b_{2}(u;\theta)=\theta_{1}+\theta_{3}(u_{2}-1)\,$ and $\ b_{3}(u;\theta)=1+\theta_{1}(u_{1}-1)+\theta_{2}(u_{2}-1)+2\theta_{3}(u_{1}-1)(u_{2}-1)$. Moreover,

[TABLE]

where $\chi^{2}_{11},\chi^{2}_{12},\ldots$ are independent $\chi^{2}$ variates with one degree of freedom and the set $\{\lambda_{j}\}$ are the non-null eigenvalues of the operator $C(\theta)$ defined on the function space $\{\tau:\mathbb{N}_{0}^{2}\to\mathbb{R},\ \text{such that}\ E_{\theta}\!\left\{\tau^{2}(\boldsymbol{X})\right\}<\infty,\forall\theta\in\Theta\}$, as follows

[TABLE]

with

[TABLE]

The asymptotic null distribution of $T_{n,w}$ does not provide a useful approximation to its null distribution since it depends on the unknown true value of $\theta$. This could be overcome by replacing $\theta$ by $\widehat{\theta}_{n}$. But the greatest difficulty is to determine the set $\{\lambda_{j}\}$, since, in general, calculating the eigenvalues of an operator is not an easy task and in our case we must also obtain expression (6), which is not easy to derive. So, we next consider another way of approximating the null distribution of the test statistic, the bootstrap.

The following result proves that the bootstrap method consistently approximates the null distribution of $T_{n,w}$, for which we require the Assumption 2 in NJ (2014) and the previous explanations for that assumption.

Theorem 3

Let $\,\boldsymbol{X}_{1},\ldots,\boldsymbol{X}_{n}\,$ be iid random vectors from $\boldsymbol{X}=(X_{1},X_{2})\in\mathbb{N}_{0}^{2}$. Suppose that Assumption 2 in NJ (2014) holds, $\widehat{\theta}_{n}\mathop{\longrightarrow}\limits^{a.s.}\theta$, for some $\theta\in\Theta$. Then

[TABLE]

It is important to note that analogous comments follow those given after Theorem 2 in NJ (2014) and the test function for our case is presented below.

Let t^{*}_{n,w,\alpha}=\inf\{x:P_{*}\bigl{(}T^{*}_{n,w}\geq x\bigr{)}\leq\alpha\}. From Theorem 3, the test function

[TABLE]

or equivalently, the test that rejects $H_{0}$ when $p^{*}=P_{*}\!\left(T^{*}_{n,w}\geq T_{obs}\right)\leq\alpha$, is asymptotically correct, in the sense that the type I error is asymptotically equal to the nominal value $\alpha$, where $T_{obs}$ is the observed value of the test statistic $T_{n,w}$.

5 Behaviour against alternatives

As an immediate consequence of Theorems 1, 2 and 3, the next result gives the asymptotic power of the test $\Psi^{*}$ against fixed alternatives.

Corollary 1

Let $\boldsymbol{X}_{1},\boldsymbol{X}_{2},\ldots,\boldsymbol{X}_{n}$ be iid from $\boldsymbol{X}\in\mathbb{N}_{0}^{2}$ with pgf $g(u)$. Suppose that assumptions in Theorems 1 and 3 hold. If $\eta(g;\theta)>0$, then P\bigl{(}\Psi^{*}=1\bigr{)}\to 1.

As commented after Theorem 1, a simple way to ensure that $\eta(g;\theta)>0$, $\forall(X_{1},X_{2})\nsim BP(\theta)$, $\forall\,(\theta)\in\Theta$, and thus the consistency against any fixed alternative, is by choosing the weight function $w$ positive a.e. on $[0,1]^{2}$.

For the local power, the next result ensures that the test $\Psi^{*}$ is able to detect alternatives as defined in (11) in NJ (2014), which converge to the BPD at the rate $n^{-1/2}$. With this aim, let $\{\phi_{j}\}$ be the set of orthonormal eigenfunctions corresponding to the eigenvalues $\{\lambda_{j}\}$ of the operator $C(\theta)$ given in (5).

Theorem 4

Let $\boldsymbol{X}_{1},\boldsymbol{X}_{2},\ldots,\boldsymbol{X}_{n}$ be iid from $\boldsymbol{X}\in\mathbb{N}_{0}^{2}$, with pmf $P_{n}(x,y)$ as defined in (11) in NJ (2014). Suppose that Assumptions 1 and 3 in NJ (2014) hold. Then

[TABLE]

where $\,c_{k}=\mathop{\sum}\limits_{x,\,y}b(x,y)\,\phi_{k}(x,y)\,$ and $Z_{1},Z_{2},\ldots\,$ are independent standard normal variates.

6 Some computational issues

6.1 On the calculation of the test statistic

Using the weight function (11) in NJ (2014) we obtained the following expression of our statistic.

[TABLE]

[TABLE]

where $\,B_{rs}=\{X_{rs}\geq 1\},\,\ 1\leq r\leq n,\ s=1,2$.

6.2 On the calculation of the null bootstrap distribution estimator

In practice, the exact bootstrap estimator of the null distribution of $T_{n,w}$ cannot be calculated, we will approximate it by simulation following the parametric bootstrap procedure (PB algorithm) given in section 4.1 in NJ (2016).

7 Numerical results

The properties studied so far describe the behavior of the proposed test for very large samples. We carried a simulation experiment in order to study the goodness of the bootstrap approximation as well as to compare the power of the proposed test with other tests for finite sample sizes. We briefly describe it in this section and display a summary of the results obtained. All computations were performed by using programs written in the R language.

7.1 Simulated data

In addition to the test proposed in this paper, $T_{n,a}$, we also considered the tests given in Crockett (1979) (denoted by $T$, see subsection 2.1.1), Loukas and Kemp (1986) (denoted by $I_{B}$, see subsection 2.1.2), Rayner and Best (1995) (denoted by $NI_{B}$, see subsection 2.1.3), NJ (2014) (denoted by $R_{n,a}$ and $S_{n,a}$, see subsections 2.1.4 and 2.1.5, respectively) and NJ (2016) (denoted by $W_{n}$, see subsection 2.1.6).

We studied the goodness of the proposed bootstrap approximations to the null distribution of the test statistic for finite sample sizes. With this aim, we generated 1,000 samples of size $n=30(20)70$ from $BP(\theta_{1},\theta_{2},\theta_{3})$, with $\theta_{1}=\theta_{2}=1$ and $\theta_{3}$ such that the correlation coefficient, $\rho={\theta_{3}}/{\sqrt{\theta_{1}\,\theta_{2}}}$, equals 0.25, 0.5 and 0.75. To estimate $\theta$ we employed the maximum likelihood method. Then we approximate the $p$-values bootstrap of the proposed tests with weight function (11) in NJ (2014) for $a=(a_{1},a_{2})\in\{(0,0),(1,0)\}$ and 500 bootstrap samples, as well as the (asymptotic) $p$-values associated with the test statistics $T$, $I_{B}$ and $NI_{B}$.

We repeated the above experiment for $\theta_{1}=1.5$, $\theta_{2}=1$ and $\theta_{3}$ such that the correlation coefficient (approximately) equals 0.25, 0.5 and 0.75. In this case, since $\theta_{1}\neq\theta_{2}$, we considered $(a_{1},a_{2})\in\{(0,0),(1,0),(0,1)\}$ for $R_{n,a}$, $S_{n,a}$ and $T_{n,a}$ in order to examine the effect of giving different weight to each component when they have different means.

Tables I and II display the fraction of estimated $p$-values less than or equal to $0.05$ and $0.10$, which are the estimated type I error probabilities for $\alpha=0.05$ and $0.10$ (denoted as f05 and f10 in the tables), respectively.

To measure the performance of the considered approximations, we calculated the $p$-value of the Kolmogorov-Smirnov test statistic of uniformity (KS) for each set of 1,000 values obtained for each test statistic. These values were rounded to 2 decimal places.

Looking at these tables we conclude that the asymptotic approximation to the $p$-values works better for $T$ than for $I_{B}$ and $NI_{B}$. Nevertheless, none of them give satisfactory results even for $n=70$. By contrast, the bootstrap provides an accurate approximation of the null distribution of $T_{n,a}$ in all tried cases. As for the choice of $a_{1}$ and $a_{2}$, we observe that there is no gain in performance when $a_{1}\neq a_{2}$.

To study the power we repeated the above experiment for samples with size $n=50$ and we use the same alternative distributions used in NJ (2014), some of which have also been taken as alternatives by other researchers (see, e.g. Loukas and Kemp, 1986; Rayner and Best, 1995, and NJ, 2016).

The parameters of these alternatives were chosen for the same reason given by NJ (2014). We took $a_{1}=a_{2}=0$ because, as observed from the results in the previous experiment, there is no gain in performance when $a_{1}\neq a_{2}$ when approximating the probability of type I error. In addition, taking $a_{1}=a_{2}=0$ is less time consuming.

Table III displays the alternatives considered and the estimated power for nominal significance level $\alpha=0.05$. The results presented in this table allow us to conclude that the new test proposed in this paper is able to detect all the alternatives treated and with a power as good or better than the other tests based on the bootstrap method, while the non-consistent tests are not able to detect most of these alternatives, especially tests $I_{B}$ and $NI_{B}$.

As we stated, the test we propose is faster than its competitors, the Table IV presents the results obtained.

7.2 Real data sets

To end this section, $T_{n,a}$ is applied to a real data set. This data set was analyzed in Bermúdez (2009), who used two variables, the number of claims for third-party liability ($X_{1}$) and the number of claims for the rest of guarantees ($X_{2}$). The original sample comprised a ten percent sample of the automobile portfolio of a major insurance company operating in Spain in 1995. The author assumed that $(X_{1},X_{2})$ has a BPD, but according to the report shown in Table V, the data set is not well modeled by a BPD. The blanks are due to the fact that $W_{n}$ does not depend on the value of $(a_{1},a_{2})$.

7.3 Case $\boldsymbol{\theta_{3}=0}$

This case has been excluded from $H_{0}$ because it is a boundary point. This situation occurs when the variables $X_{1}$ and $X_{2}$ are independent and is analyzed in NJ (2016) where different ways of approaching it are given, besides references are cited for a detailed treatment, even it is a subject for a future research.

8 Extension of $\boldsymbol{T_{n,w}}$

In principle, the approach can be generalized to the case $m\geq 3$ and we refer to a manuscript that is uploaded to arXive math.

To illustrate this situation we will present the case for $m=3$, in which we need to satisfy $7=\sum_{i=1}^{3}\binom{3}{i}$ equations to obtain a characterization of the respective Poisson distribution. It can be seen that the number of equations grows following a sum of combinatorial numbers due to differential equations of different order that must be verified, which range from order 1 to order $m$. These equations arise due to the philosophy of the method to characterize the respective Poisson distribution.

8.1 Trivariate case

For this particular case, from section 2.2, for $m=3$, let

[TABLE]

where $Y_{1},Y_{2},Y_{3},Y_{4}$ are independent Poisson random variables with means $\theta^{\prime}_{1}=\theta_{1}-\theta_{4}>0,\theta^{\prime}_{2}=\theta_{2}-\theta_{4}>0,\theta^{\prime}_{3}=\theta_{3}-\theta_{4}>0$ and $\theta_{4}>0$, respectively. The joint distribution of the vector $(X_{1},X_{2},X_{3})$ is called a trivariate Poisson distribution (TPD) with parameter $\theta=(\theta_{1},\theta_{2},\theta_{3},\theta_{4})$ (see, e.g. Johnson, Kotz and Balakrishnan, 1997; Loukas and Papageorgiou, 1991). The joint pgf of $(X_{1},X_{2},X_{3})$ is

[TABLE]

The empirical counterpart of pgf is epgf of the data given by

[TABLE]

Now, the objective is to test the hypothesis

[TABLE]

To achieve this new objective, we give a characterization for the TPD.

Proposition 2

Let $g(u_{1},u_{2},u_{3};\theta)\,$ be as defined in (7). Then $g(u_{1},u_{2},u_{3};\theta)\,$ is the only pgf satisfying the following system of partial differential equations

[TABLE]

where $\,h(u;\theta)=\prod_{i=1}^{3}\left\{\theta_{i}+\theta_{4}\left(\prod_{j\neq i}u_{j}-1\right)\right\}+\theta_{4}\left(1+\sum_{k=1}^{3}u_{k}\left\{\theta_{k}+\theta_{4}\left(\prod_{j\neq k}u_{j}-1\right)\right\}\right)$.

By Proposition 1 in NJ (2014), $g(u)$ and its derivatives can be consistently estimated by the epgf and the derivatives of the epgf, respectively. Thus, if $H_{03}$ is true, then the functions

[TABLE]

should be close to [math], $\forall(u_{1},u_{2},u_{3})\in[0,1]^{3}$, where $\widehat{\theta}_{n}=(\widehat{\theta}_{1n},\widehat{\theta}_{2n},\widehat{\theta}_{3n})$ is a consistent estimator of $\theta$ and $g_{n}(u_{1},u_{2},u_{3})$ is the epgf associated with the data, i. e.,

[TABLE]

Thus, to test $H_{03}$ we consider the following test statistic

[TABLE]

where $w(u)$ is a measurable non-negative function with finite integral over $[0,1]^{3}$. Similar results to those stated in Sections 3, 4, and 5 for the bivariate case can be established for $T_{3,n,w}$.

Remark 1

So far we have not managed to obtain numerical results for the case $m=3$ due to the large number of calculations involved in $T_{3,n,w}$. We can assure that this new test is not recommended for $m\geq 3$ and it is preferable to use the $W_{n}$ statistic.

8.2 Simulated data for the trivariate case

To simulate type I error we follow a procedure similar to that described for the bivariate case, but we do not have competitors. We consider three situations: a) $\theta_{1}=\theta_{2}=\theta_{3}$, b) $\theta_{1}=\theta_{2}\neq\theta_{3}$, $\theta_{2}=\theta_{3}\neq\theta_{1}$ and c) $\theta_{1}\neq\theta_{2}$ and $\theta_{1}\neq\theta_{3}$ and $\theta_{2}\neq\theta_{3}$. In each of these cases $\theta_{1}$, $\theta_{2}$, $\theta_{3}>\theta_{4}$. In addition, $\theta_{4}$ was chosen in such a way that the correlation coefficients, $\rho=(\rho_{12},\rho_{13},\rho_{23})$, were equal or very close to 0.25, 0.5, 0.75 and 1.00, where $\rho_{ij}=\frac{Cov(X_{i},X_{j})}{\sqrt{Var(X_{i})Var(X_{j})}}$.

Tables VI and VII display the fraction of estimated $p$-values less than or equal to $0.05$ and $0.10$, which are the estimated type I error probabilities for $\alpha=0.05$ and $0.10$ (denoted as f05 and f10 in the tables), respectively.

As we had anticipated in Remark 1, the proposed new test $T_{3,n}$ is not faster than some of its competitors, as can be seen in the Table XII.

8.3 Real data set for trivariate case

The data set was analyzed in Catalina Bolancé & Raluca Vernic (2017), the data come from the Spanish insurance market and consist of a random sample of 162,019 policyholders who had had one or more auto and home policies during the decade 2006-2015. Catalina Bolancé & Raluca Vernic (2017) used three dependent variables: the number of claims in auto insurance at fault involving only property damage ($X_{1}$); the number of claims in auto insurance at fault with bodily injury ($X_{2}$); and, the number of claims in home insurance at fault ($X_{3}$).

Table XIII shows the p-values obtained by applying the test we propose. It is concluded that the data do not come from a trivariate Poisson distribution, this is in agreement with the researchers who used this data set to model a trivariate Sarmanov distribution.

Acknowledgements

The author would like to thank the Departamento de Investigación de la Universidad del Bío-Bío and the Grupo de Investigación Matemática Aplicada GI 172409/C de la Universidad del Bío-Bío, Chile. He also thanks the anonymous reviewers and the editor of this journal for their valuable time and their careful comments and suggestions with which the quality of this paper has been improved. Thanks to Florencia Osorio for her valuable comments in the revision of the manuscript.

Appendix

Proofs

The proofs of Theorems 1 to 4 are quite similar to those of Theorems 3, 1, 2 and 4 in NJ (2014), respectively. Here we give a sketch of the proof of Proposition 1. The proof of Proposition 2 follows the steps of the proof of Proposition 1 by occupying the recurrence relationships for the probabilities and their respective partial derivatives given in Loukas and Papageorgiou (1991). A detailed derivation of the results can be obtained from the authors upon request.

Proof of Proposition 1   Let $(X_{1},X_{2})$ be a random vector and let $g(u_{1},u_{2})=\sum_{i,j\geq 0}P_{ij}\,u_{1}^{i}u_{2}^{j}$ be its pgf, where $P_{ij}=P(X_{1}=i,X_{2}=j)$. Let $f(u_{1},u_{2};\theta)=c_{0}+c_{1}u_{1}+c_{2}u_{2}+c_{3}u_{1}u_{2},\,$ with $\,c_{0}=\theta_{3}+(\theta_{1}-\theta_{3})(\theta_{2}-\theta_{3}),\ c_{1}=\theta_{3}(\theta_{1}-\theta_{3}),$ $c_{2}=\theta_{3}(\theta_{2}-\theta_{3})\,$ and $\,c_{3}=\theta_{3}^{2}.\,$ Then

[TABLE]

From the first equation in (4), $D_{1}(u;\theta)=0$, then by matching coefficients, we obtain

[TABLE]

With enough algebraic work we can demonstrate that equations (8) satisfy (1) or (2) and (5) in Kawamura (1985). Moreover, the last two equations in (4) satisfy (3) and (4) in Kawamura (1985). Therefore, the result is obtained by applying Theorem 3 in Kawamura (1985).