GPIGINAR(1)模型的贝叶斯估计

Bayesian Estimation of the GPIGINAR(1) Model

摘要: 当前研究表明，整数值时间序列通常呈现出重尾特征。然而，能够处理这类具有重尾特征数据的模型相对较少。本文聚焦于广义泊松逆高斯INAR(1) (GPIGINAR(1))模型，该模型能够拟合重尾特征。文中采用自适应马尔可夫链蒙特卡洛(MCMC)算法，对该模型中所关注的参数进行贝叶斯估计，并将其结果与极大似然估计结果进行对比，以此证明所提方法的有效性。最后，通过具体实例进一步说明该模型在处理此类数据时的有效性，以及该算法在解决这类问题时的可行性。

Abstract: Current research shows that integer-valued time series usually exhibit heavy-tailed characteristics. However, relatively few models are capable of handling data with such heavy-tailed characteristics. This paper focuses on the Generalized Poisson Inverse Gaussian INAR(1) (GPIGINAR(1)) model, which is able to fit heavy-tailed features. In this paper, the Adaptive Markov Chain Monte Carlo (MCMC) algorithm is adopted to perform Bayesian estimation on the parameters of interest in this model. The results are then compared with those of maximum likelihood estimation to prove the effectiveness of the proposed method. Finally, specific examples are used to further illustrate the effectiveness of this model in processing such data and the feasibility of the algorithm in solving such problems.

1. 引言

对于重尾型的整数值时间序列，还很少被考虑，在时间序列中，经常可以观察到数据呈现重尾的特征，这也就意味着尾部概率是不可忽略的。如果忽略尾部，那么可能会导致有用信息的丢失。对于具有重尾特征的数据，Maiti等(2018) [1]提出了一类新的INAR(1)过程，用于模拟存在大量零和一的低计数的重尾时间序列。泊松逆高斯(PIG)分布，可以建模比负二项分布和泊松分布尾部要厚的数据。该分布已被许多研究者研究[2]-[5]，Barreto-Souza (2019) [6]介绍了边际分布为PIG分布的INAR模型。钱连勇(2021) [7]研究了若干的重尾整数值时间序列模型，其中基于GPIG分布为新息项的INAR模型为本文所讨论的模型。

随着贝叶斯统计影响的扩大，学者们对贝叶斯方法进行了深入研究。Bolstad (2007) [8]指出，使用先验信息的贝叶斯分析通常会比忽略先验信息的频率论者提供更好的推断。通常情况下，后验分布的计算需要进行积分，但大部分情况是该积分没有显式解，因此可以采用MCMC的方法来从复杂的后验分布中获取样本，从而实现对参数的贝叶斯推断。Hay和Pettitt (2001) [9]率先对广义线性混合模型的参数展开了全面的贝叶斯分析，并深入探究其在传染病控制方面的应用，为后续贝叶斯方法在各类模型中的运用奠定了基础。从自适应MCMC算法的历程来看，Chen和Lee (1995) [10]创新性地提出自适应MCMC算法，用于估计阈值自回归(TAR)模型和双线性模型，开启了自适应MCMC算法在模型参数估计领域的应用篇章。进入21世纪，相关研究不断拓展深化。Chen和So (2006) [11]运用MCMC自适应算法对阈值异方差模型进行参数估计，进一步拓展MCMC自适应算法的应用范围。Chen等(2019) [12]针对气象变量存在的滞后依赖、过分散、连续零、非线性动力学及时变系数等复杂问题，在贝叶斯框架下成功构建马尔可夫切换泊松INGARCH模型，有效解决了气象数据处理中的难题。与此同时，整数值时间序列模型的研究也取得显著进展。Ausín和Galeano (2007) [13]研究高斯混合广义自回归条件异方差(GARCH)模型的贝叶斯估计问题，Bu和McCabe (2008) [14]运用MCMC方法深入探讨INAR ($p$ )模型的模型选择、估计及预测问题，朱复康和李琦(2009) [15]聚焦INGARCH模型的贝叶斯估计问题，Neal和Rao (2007) [16]探究整数值ARMA过程的贝叶斯估计，张哲等(2010) [17]借助贝叶斯方法对INAR(1)模型的参数估计展开研究，Drovandi等(2016) [18]则考虑了低维整数值时间序列模型的贝叶斯推理。这些研究逐步完善了整数值时间序列模型的贝叶斯分析体系。在特定模型的参数估计与推断方面，众多学者成果丰硕。Truong等(2017) [19]提出在MCMC采样方法下对泊松INGARCH模型进行参数估计，并采用贝叶斯信息准则进行模型比较，Chen和Lee (2017) [20]提出基于泊松、负二项和对数线性泊松INGARCH模型的贝叶斯因果检验方法，并应用于气候和犯罪数据，Chu和Yu (2023) [21]对对数BNB-INGARCH模型应用MCMC自适应算法进行贝叶斯推断。此外，王佩(2019) [22]对门限ARMA模型的贝叶斯估计展开研究，Li等(2019) [23]提出带有解释变量的随机波动模型，并思考该模型的贝叶斯估计问题，Garay等(2020) [24]通过设置潜在变量的形式，利用M-H算法对ZINAR($p$ )模型进行贝叶斯推断，黄菊(2021) [25]研究空间自回归条件异方差模型的贝叶斯估计，Yang等(2022) [26]采用MCMC算法对阈值一阶自回归模型进行贝叶斯推断，于鑫洋(2022) [27]研究自激励整数值门限自回归模型的贝叶斯估计及应用。Kang Y等(2024) [28]对改进的混合二项自回归模型进行了贝叶斯分析，并将其应用于雨天数和空气质量水平数据研究。因此本文将基于MCMC算法讨论GPIGINAR(1)模型的贝叶斯估计问题。

2. 模型定义及其性质

2.1. GPIG分布

在本节中介绍GPIG分布的定义及其基本性质，Zhu (2009) [29]中给出了定义：

定义2.1. 若随机变量$Y$ 有如下所示的概率生成函数的形式：

$G_y\left(h;a,b,c\right)=\exp \left\{b\left[{\left(1-c\right)}^a-{\left(1-ch\right)}^a\right]\right\}$ (1)

其中$0<a\le 1$ ，$b>0$ ，$0\le c\le 1$ 。则称该随机变量$Y$ 服从GPIG分布，记为$GPIG\left(a;b;c\right)$ 。

注2.1． 接下来讨论边界与特殊情况：

(1) 当$a=1$ ，$c\ne 0$ 时，其概率生成函数为$G_y\left(h;a,b,c\right)=\exp \left\{bc\left(h-1\right)\right\}$ ，是均值为$bc$ 的泊松分布。

(2) 当$c=1$ 且$0<a\le 1$ 时，其概率生成函数为$G_y\left(h;a,b,c\right)=\exp \left\{-b{\left(1-h\right)}^a\right\}$ ，是稳定的离散分布。

(3) 当$a=0$ 或$c=0$ 时，其概率生成函数为$G_y\left(h;a,b,c\right)=1$ ，是零膨胀分布。

(4) 当$a=\frac{1}{2}$ ，$b=\lambda {\theta }^{-1}\sqrt{1+2{\theta }^2/\lambda }$ 且$c=2{\theta }^2{\lambda }^{-1}/\left(1+2{\theta }^2/\lambda \right)$ 时，其概率生成函数为 $G_y\left(h;a,b,c\right)=\exp \left\{\lambda {\theta }^{-1}\left(1-{\left[1+2{\theta }^2{\lambda }^{-1}\left(1-h\right)\right]}^{\frac{1}{2}}\right)\right\}$ ，

其中$\lambda >0$ ，$\theta >0$ 时，为PIG分布。

注2.2. 假设$Y~GPIG\left(a,b,c\right)$ ，有

${\left(1-ch\right)}^{1-a}{G}^{\prime }\left(h\right)=abcG\left(h\right)$ (2)

当$0\le c<1$ 时，其均值$\mu$ ，方差${\sigma }^2$ ，离散指数$I$ 以及偏度$\gamma$ 分别为：

$\mu =abc/{\left(1-c\right)}^{1-a}$ ，${\sigma }^2=\frac{abc\left(1-ac\right)}{{\left(1-c\right)}^{2-a}}$ ，$I=\frac{{\sigma }^2}{\mu }=\frac{1-ac}{1-c}$ ，$\gamma =\frac{{\left(1-c\right)}^{a-3}abc\left[1+c-3ac+a^2{c}^2\right]}{{\sigma }^3}$ 。

概率密度函数$\left\{m_k\right\}$ 的具体形式通过递推算法得到，对于GPIG分布来说，$m_0$ 与$m_1$ 很容易被计算，因为它们能够通过计算概率生成函数在$h=0$ 处的一阶导数和二阶导数的值。对式(2)进行泰勒展开，存在实数序列$\left\{l_i:i=1,2,\cdot \cdot \cdot \right\}$ ，满足以下等式：

$\begin{array}{c}{\left(1-ch\right)}^{1-a}=1-\left(1-ah\right)ch-\frac{\left(1-a\right)a}{2!}{c}^2{h}^2-\frac{\left(1-a\right)a\left(1+a\right)}{3!}{c}^3{h}^3-\cdot \cdot \cdot \\ =1-l_{1}h-l_2{h}^2-l_3{h}^3-\cdot \cdot \cdot \end{array}$

$l_1=\left(1-a\right)c$ ，$l_{i+1}=\left(\frac{i-1+a}{i+1}\right)c{l}_i$ ，$i=1,2,\cdot \cdot \cdot$ 。

如果$G_Y\left(h;a,b,c\right)={\displaystyle \sum _i{m}_i}{h}^i$ 且有$m_i=m_i\left(a,b,c\right)$ ，对式(2)进行泰勒展开并进行再一次整理，可以得到：

$\begin{array}{c}abc\left[m_0+m_{1}h+m_2{h}^2+m_3{h}^3+\cdots \right]\\ =\left[1-l_{1}h-l_2{h}^2-l_3{h}^3-\cdots \right]\times \left[m_1+2m_{2}h+3m_3{h}^2+4m_4{h}^3+\cdots \right]\\ =m_1+\left(2m_2-m_1{l}_1\right)h+\left(3m_3-2m_2{l}_1-m_1{l}_2\right)h^2+\left(4m_4-3m_3{l}_1-m{p}_2{l}_2-m_1{l}_3\right)h^3\end{array}$

从而可以得到概率密度函数$\left\{m_k\right\}$ 的递归关系有：

$m_0=e^{b\left[{\left(1-c\right)}^a-1\right]}$ ，$m_1=abc{m}_0$ ，$m_{r+1}=\frac{1}{1+r}\left(abc{m}_r+{{\displaystyle \sum }}_{i=1}^{r}i{l}_{r+1-i}{m}_i\right)$ ，$r=1,2,\cdot \cdot \cdot$ 。

2.2. GPIGINAR(1)模型

定义2.2. 基于二项稀疏算子“$\circ$ ”，以GPIG分布为新息项的一阶整值自回归INAR(1)过程，即有如下表示形式：

$Y_t=\alpha \circ Y_{t-1}+{\epsilon }_{t}t\in Z$ (3)

其中$0\le \alpha <1$ ，${\epsilon }_t$ 是$i.i.d$ 的$GPIG\left(a,b,c\right)$ 随机变量，与$X_s\left(s<t\right)$ 独立。这个过程简称为GPIGINAR(1)过程。

由模型的定义及概率论的基础知识可以容易得到该模型如下的一些基本性质：

(1) $X_t$ 的条件期望为：$E\left(X_t|X_{t-1}\right)=\alpha X_{t-1}+E\left({\epsilon }_t\right)$ ，

(2) $X_t$ 的条件方差为：$Var\left(X_t|X_{t-1}\right)=\alpha \left(1-\alpha \right)X_{t-1}+Var\left({\epsilon }_t\right)$ ，

(3) $X_t$ 的期望为：$E\left(X_t\right)=\frac{{\mu }_{ϵ}}{1-\alpha }=\frac{abc}{\left(1-\alpha \right){\left(1-c\right)}^{1-a}}$ ，

(4) $X_t$ 的方差为：$Var\left(X_t\right)=\frac{a{\mu }_{ϵ}+{\sigma }_{ϵ}^2}{1-{\alpha }^2}=\frac{\alpha abc}{{\left(1-\alpha \right)}^2{\left(1-c\right)}^{1-a}}+\frac{abc\left(1-ac\right)}{{\left(1-\alpha \right)}^2{\left(1-c\right)}^{2-a}}$ ，

(5) 自协方差函数与自相关函数为：${\gamma }_k=\mathrm{Cov}\left(X_t,X_{t-k}\right)={\alpha }^k{\gamma }_0$ ，${\rho }_k={\alpha }^k$ 。

该模型的转移概率为：

$p_{mj}=p\left(X_t=j|X_{t-1}=m\right)={\displaystyle \sum _{l=0}^{min(j,m)}\left(\begin{array}{c}m\\ l\end{array}\right)}{\alpha }^l{\left(1-\alpha \right)}^{m-l}p\left({\epsilon }_t=j-l\right)$ (4)

其中：

$p\left({ϵ}_t=j-l\right)=\{\begin{array}{c}{p}_0=e^{b\left[{\left(1-c\right)}^a-1\right]}，j-l=0,\\ p_1=abc{p}_0，j-l=1,\\ p_{k+1}\frac{1}{1+k}\left(abc{p}_k+{{\displaystyle \sum }}_{i=1}^{k}i{r}_{k+1-i}{p}_i\right)，j-l=k+1，k\ge 1.\end{array}$

3. 参数估计

设序列$X={\left(X_1,X_2,\cdot \cdot \cdot ,X_n\right)}^T$ 是来自GPIGINAR(1)模型的一组观测数据，其中参数$\theta ={\left(a,b,c,\alpha \right)}^T$ 为未知参数，要进行贝叶斯估计，首先写出似然函数为：

$L\left(X|\theta \right)={\displaystyle \sum _{t=2}^{n}l}\left(X|\theta \right)={\displaystyle \sum _{t=2}^{n}l}ogP\left(X_t=x_t|X_{t-1}=x_{t-1}\right)$ (5)

令$S_t=\alpha \circ X_{t-1}$ ，$t=1,\cdot \cdot \cdot ,n$ ，那么

${\epsilon }_t=X_t-S_t$ (6)

在给定$\left(\alpha ,a,b,c,X_{t-1}\right)$ 的情况下，$\left(S_t,{\epsilon }_t\right)$ 是相互独立的。因此，对于$t\ge 1$ 有

$P_{S_t,{\epsilon }_t}\left(S_t,{\epsilon }_t|X_{t-1},\alpha ,a,b,c\right)={\displaystyle \sum _{t=1}^n\left(\begin{array}{c}{X}_{t-1}\\ S_t\end{array}\right)}{\alpha }^{S_t}{\left(1-\alpha \right)}^{X_{t-1}-S_t}P\left({\epsilon }_t=X_t-S_t\right)$ (7)

其中$S_t\le X_{t-1}$ ，新的观测值由$X_t=S_t+{\epsilon }_t$ 给出。令$X=\left(X_0,X_1,\cdot \cdot \cdot ,X_n\right)$ ，$S_t=\left(S_1,S_2,\cdot \cdot \cdot ,S_n\right)$ ， ${\epsilon }_t=\left({\epsilon }_1,{\epsilon }_2,\cdot \cdot \cdot ,{\epsilon }_n\right)$ ，那么，模型(3)的对数似然函数可以改写为

$L_{S_t,{\epsilon }_t}\left(S_t,{\epsilon }_t|\alpha ,a,b,c\right)=\log P_{S_t,{\epsilon }_t}\left(S_t,{\epsilon }_t|X_{t-1},\alpha ,a,b,c\right)$ (8)

3.1. 先验

为了完成GPIGINAR(1)过程的贝叶斯分析，需要定义模型中未知参数的先验分布，假设参数$\theta$ 是相互独立的，有

$P(\theta )=P(a)P(b)P(c)P(\alpha )$ (9)

本文为参数选取先验分布，考虑以下先验设置：

$P\left(a\right)=U\left(0,1\right)$ ，$P\left(\alpha \right)=U\left(0,1\right)$ ，$P\left(c\right)=U\left(0,1\right)$ ，$P\left(b\right)=Gamma\left(m,n\right)\propto b^{m-1}{e}^{-nb}$ 。

即为参数$a$ ，$c$ ，$\alpha$ 选取无信息先验，为参数$b$ 选取共轭伽马先验，其中参数$m$ 与$n$ 都为超参数，且它们都大于0。结合似然和先验，得到如下的联合后验密度函数：

$\begin{array}{r}f\left(S,\epsilon ,\alpha ,a,b,c|X\right)={\displaystyle \sum _{t=2}^{n}l}og{P}_{S_t,{\epsilon }_t}\left(S_t,{\epsilon }_t|X_{t-1},\alpha ,a,b,c\right)P\left(a\right)P\left(b\right)P\left(c\right)P\left(\alpha \right)\\ \propto {\displaystyle \sum _{S_t=0}^{min\left(x_t,x_{t-1}\right)}\left(\begin{array}{c}{x}_{t-1}\\ S_t\end{array}\right)}{\alpha }^{S_t}{\left(1-\alpha \right)}^{x_{t-1}-S_t}p\left({\epsilon }_t=x_t-S_t\right)b^{m-1}{e}^{-nb}\end{array}$ (10)

在理论上，可以对式(10)中的$\theta$ 进行积分，以获得参数的边际后验密度。然而这一过程复杂没有显示解。因此，本文采用Gibbs抽样、随机游走链与M-H算法相结合的抽样方法来进行贝叶斯估计。为此推导出参数$\alpha$ 与$S_t$ 的满条件分布，其形式如下：

(1) 参数$\alpha$ 的满条件分布

$\begin{array}{r}P(\alpha |S,X)\propto {\displaystyle \prod _{t=1}^n{\alpha }^{S_t}{\left(1-\alpha \right)}^{X_{t-1}-S_t}}\\ \propto {\alpha }^{{\displaystyle \sum _{t=1}^n{S}_t}}{\left(1-\alpha \right)}^{{\displaystyle \sum _{t=1}^n{X}_{t-1}-S_t}}\end{array}$ (11)

由式(11)可知，$\alpha$ 的满条件分布为$Beta\left(1+{\displaystyle \sum _{t=1}^n{S}_t},1+{\displaystyle \sum _{t=1}^n{X}_{t-1}-S_t}\right)$ 。

(2) 潜变量$S_t$ 的满条件分布，$t=1,\cdot \cdot \cdot ,n$

$P(S_t|\alpha ,X)\propto \left(\begin{array}{c}{X}_{t-1}\\ S_t\end{array}\right){\alpha }^{S_t}{\left(1-\alpha \right)}^{x_{t-1}-S_t}$ (12)

由式(12)可知，$S_t$ 的满条件分布为$Bin\left(X_{t-1},\alpha \right)$ ，$t=1,\cdot \cdot \cdot ,n$ 。

3.2. 抽样方案

(1) 更新参数$\alpha$ ：从$Beta\left(1+{\displaystyle \sum _{t=1}^n{S}_t},1+{\displaystyle \sum _{t=1}^n{X}_{t-1}-S_t}\right)$ 中抽取$\alpha$ 。

(2) 更新$\left(S,\epsilon \right)$ ，采用M-H算法：

步骤1：从$Bin\left(X_{t-1},\alpha \right)$ 中抽取$S^{\prime }$ ，

步骤2：若$X_t<{{\epsilon }^{\prime }}_t$ ，则重复步骤1，

步骤3：令${{\epsilon }^{\prime }}_t=X_t-S_t$ ，

步骤4：计算接受概率：$P=\min \left\{1,\frac{f\left({\epsilon }_t{}^{\prime }\right)}{f\left({\epsilon }_t\right)}\right\}$ ，其中${{\epsilon }^{\prime }}_t$ 是${\epsilon }_t$ 的后一次抽样值。

(3) 由于${\epsilon }_t$ 的后验分布没有明确的表达式，故不能直接得出参数$\left(a,b,c\right)$ 的估计值，本文采用了Chen和So (2006) [11]提出的自适应MCMC算法。该过程一共进行$N$ 次迭代，其中包括在前$M$ 次迭代中采用随机游走M-H算法来抽取样本，为燃烧期。之后的$N-M$ 次迭代中采用独立核M-H算法来抽取参数样本。具体步骤如下：

步骤1：设置初始值${\theta }^{\left(0\right)}={\left(a^{\left(0\right)},b^{\left(0\right)},c^{\left(0\right)}\right)}^T$ 。

步骤2：当$1\le h\le M$ 时，对于参数${\theta }^{\left(j\right)}$ 采用随机游走M-H算法，具体过程如下所示：

步骤2.1：选取候选样本${\theta }^{\ast }=\left\{{\theta }_j^{\ast }\right\}$ ，$j=1,2,\cdot \cdot \cdot ,4$ 。其中${\theta }_j^{\ast }={\theta }_j{}^{\left(h-1\right)}+m\epsilon$ ，$\epsilon ~N\left(0,{\sigma }_j^2\right)$ 。其中${\theta }_j{}^{\left(h-1\right)}$ 表示第$h-1$ 次迭代，$m$ 表示步长。

步骤2.2：选取步骤2.1中${\theta }^{\ast }$ 满足$0<{\theta }_1<1$ ，${\theta }_2>0$ ，$0<{\theta }_3<1$ 条件的候选值。否则返回步骤2.1。

步骤2.3：计算接受率

$P\left({\theta }^{\ast },{\theta }_j{}^{\left(h-1\right)}\right)=\min \left\{1,\frac{f\left({\theta }^{\ast }\right)}{f\left({\theta }_j{}^{\left(h-1\right)}\right)}\right\}$ (13)

其中$f(\cdot )$ 表示目标后验分布。

步骤2.4：从均匀分布中随机生成$U$ ，$U\in \left[0,1\right]$ 。若$U<P\left({\theta }^{\ast },{\theta }_j{}^{\left(h-1\right)}\right)$ ，则接受新的候选值，即${\theta }_j^{\ast }={\theta }_j{}^{\left(h\right)}$ 。否则，设置${\theta }_j^{\left(h\right)}={\theta }_j{}^{\left(h\right)}$ 。

步骤3：当$h\ge M+1$ ，对于样本${\theta }_j^{\left(h\right)}$ 采用独立核M-H算法：

步骤3.1：从${\theta }_j^{\ast }={\mu }_{\theta }+\epsilon$ 中生成候选值${\theta }_j^{\ast }$ ，其中$\epsilon ~N\left(0,{\Omega }_{\theta }\right)$ 。${\mu }_{\theta }$ 和${\Omega }_{\theta }$ 分别为样本均值和样本协方差矩阵，其计算采用的是燃烧期前$M$ 次迭代的样本。

步骤3.2：选取步骤3.1中${\theta }_j{}^{\ast }$ 满足$0<{\theta }_1^{\ast }\le 1$ ，${\theta }_2^{\ast }>0$ ，$0\le {\theta }_3^{\ast }<1$ 条件的候选值。否则返回步骤3.1。

步骤3.3：计算接受率

$P({\theta }_j^{\ast },{\theta }_j^{\left(h-1\right)})=\min \left\{1,\frac{f\left({\theta }_j^{\ast }\right)k\left({\theta }_j^{\left(h-1\right)}\right)}{f\left({\theta }_j^{\left(h-1\right)}\right)k\left({\theta }_j^{\ast }\right)}\right\}$ (14)

其中${\theta }_j^{\left(h\right)}$ 表示${\theta }_j$ 的第$h$ 次迭代，$k(\theta )$ 表示是多元正态分布，样本均值${\mu }_{\theta }$ 和样本协方差${\Omega }_{\theta }$ 采用前$M$ 次迭代${\theta }_j$ 的样本。

步骤3.4：从均匀分布中随机生成$U$ ，$U\in \left[0,1\right]$ 。若$U<P\left({\theta }^{\ast },{\theta }_j^{\left(h-1\right)}\right)$ ，则接受新的候选值，即${\theta }_j^{\ast }={\theta }_j^{\left(h\right)}$ 。否则，设置${\theta }_j{}^{(h)}={\theta }_j{}^{(h-1)}$ 。

步骤4：如果链没有达到收敛，则重复进行直至链收敛。

4. 数值模拟

为了检验MCMC算法对GPIGINAR(1)模型的性能，从GPIGINAR(1)模型中生成不同样本量大小的样本，在此数据下使用提出的MCMC算法对每个数据集进行拟合，获得后验样本，从而估计模型的参数。先验分布为3.1节中所提到的，计数时间序列$X_t$ 从GPIGINAR(1)模型中生成，在接下来的数值模拟中，对不同的样本量的数据执行30000次迭代，前15000次迭代为燃烧期，重复模拟1000次，获得结果。

在本节中，考虑以下方案设置：

方案1：$\left(a,b,c,\alpha \right)=\left(0.5,3,0.55,0.1\right)$ ；

方案2：$\left(a,b,c,\alpha \right)=\left(0.55,2,0.6,0.4\right)$ 。

如图1所示，展示了方案1在样本量为300时序列的样本路径图、直方图、ACF图及PACF图。可以看出，序列为平稳序列并且存在短期自相关。可以看出，序列为平稳序列并且存在短期自相关。

Figure 1. Sequence, histogram, ACF and PACF plots with $n=300$

图1. 样本量$n=300$ 时序列的样本路径图、直方图、ACF图及PACF图

Figure 2. Convergence plots, histogram plots and ergodic mean plots for each parameter with sample size of 300

图2. 每个参数在样本量为300时的链轨迹图、后验分布直方图以及遍历均值图

在接下来的数值模拟中通过对不同样本量的样本分别进行30,000次迭代，丢弃前15,000次，重复模拟进行1000次得到其结果。为了展示3.2节中抽样方案的收敛情况，画出每个参数在样本量为300时的链轨迹图、后验分布直方图以及遍历均值图，如图2所示，可以很明显地看到其后验样本的路径没有明显的波动，后验样本的遍历均值轨迹在15,000次迭代后趋于稳定，表明后验样本的平稳性，后验样本的直方图清楚地说明各参数后验样本的分布情况。总而言之，图2提供证据表明MCMC算法的收敛性，从而也进一步表明所提出的贝叶斯估计的可行性。

为了展现两种估计效果的好坏，本节用Bias和MSE两种度量指标来衡量。为了与此方法所得出的估计结果进行比较，采用极大似然估计方法进行模拟，其结果如表1与表2所示，表1与表2中分别为方案1与方案2中重复模拟不同样本量的样本1000次后，得到的贝叶斯方法与极大似然方法的比较结果，其中包括贝叶斯方法所得到的后验均值(MC Mean)、后验中位数(MC Median)、后验标准差(MC SD)、95%置信区间、Bias、MSE，极大似然方法所得到的Bias、MSE。

根据表1及表2的结果可知：

(1) 随着样本量$n$ 从80增大到700，贝叶斯估计中多数参数的均值、中位数与真实值的Bias逐渐减小，标准差也在降低，表明估计值更集中于真实值。极大似然估计下部分参数表现类似，但参数$b$ 的变化规律较复杂，其估计值与真实值的偏差并非随样本量单调变化。

(2) 贝叶斯估计下，多数参数的Bias绝对值在不同样本量时相对较小；极大似然估计中，部分参数如$b$ 的Bias绝对值较大。多数情况下，贝叶斯估计的MSE小于极大似然估计，说明贝叶斯估计的准确性更高。例如参数$a$ ，$c$ ，$\alpha$ 在各样本量下，贝叶斯估计的MSE普遍更低。贝叶斯估计的$Q_{0.025}$ 和$Q_{0.975}$ 所表示的置信区间，以及标准差，在多数参数和样本量下相对更优，显示出估计值的离散程度和不确定性控制得更好。

(3) 综合来看，采用贝叶斯方法估计该模型更精确一点。

Table 1. Simulation results

表1. 模拟结果

参数	真实值	贝叶斯估计							极大似然估计
		MC Mean	MC Median	MC SD	$Q_{0.025}$	$Q_{0.975}$	Bias	MSE	Bias	MSE
$n=80$
$a$	0.5	0.4723	0.4742	0.2458	0.0465	0.9006	−0.0277	0.0092	−0.1252	0.1645
$b$	3	3.0130	3.0178	0.5926	1.9375	4.1455	0.0130	1.3423	1.2168	12.2101
$c$	0.55	0.5575	0.5452	0.1124	0.3797	0.7983	0.0075	0.0112	−0.0261	0.0488
$\alpha$	0.1	0.1383	0.1326	0.0700	0.0235	0.2860	0.0383	0.0051	0.0008	0.0037
$n=150$
$a$	0.5	0.4768	0.4794	0.2192	0.0921	0.8589	−0.0232	0.0220	−0.1264	0.1432
$b$	3	2.6291	2.6290	0.4246	1.8974	3.4110	−0.3709	0.6300	0.6017	2.7939
$c$	0.55	0.5128	0.5003	0.0982	0.3614	0.7222	−0.0372	0.0150	−0.0247	0.0397
$\alpha$	0.1	0.1160	0.1121	0.0556	0.0226	0.2327	0.0160	0.0026	−0.0052	0.0045
$n=300$
$a$	0.5	0.4804	0.4929	0.2214	0.0595	0.8574	−0.0196	0.0069	−0.0335	0.1101

续表

$b$	3	3.2591	3.2464	0.3390	2.6545	3.8989	0.2591	0.6603	0.4759	2.6815
$c$	0.55	0.5734	0.5632	0.0868	0.4293	0.7587	0.0234	0.0158	−0.0225	0.0287
$\alpha$	0.1	0.1044	0.1024	0.0438	0.0267	0.1942	0.0044	0.0016	−0.0040	0.0025
$n=700$
$a$	0.5	0.4877	0.4965	0.1890	0.1062	0.8184	−0.0123	0.0255	−0.0882	0.0820
$b$	3	3.1331	3.1470	0.4017	2.4321	3.8330	0.1331	0.2805	0.2661	0.6620
$c$	0.55	0.5169	0.5075	0.0682	0.4112	0.6718	−0.0331	0.0035	−0.0173	0.0155
$\alpha$	0.1	0.0988	0.0984	0.0323	0.0379	0.1626	−0.0012	0.0008	−0.0021	0.0010

Table 2. Simulation results

表2. 模拟结果

参数	真实值	贝叶斯估计							极大似然估计
		MC Mean	MC Median	MC SD	$Q_{0.025}$	$Q_{0.975}$	Bias	MSE	Bias	MSE
$n=80$
$a$	0.55	0.4459	0.4458	0.2084	0.0879	0.8133	−0.1041	0.0288	−0.1749	0.1457
$b$	2	2.6114	2.5424	0.6506	1.5788	3.9098	0.6114	0.8749	0.6475	2.5182
$c$	0.6	0.5098	0.4919	0.1341	0.3075	0.7940	−0.0902	0.0191	−0.0709	0.0533
$\alpha$	0.4	0.3992	0.3990	0.0552	0.2937	0.5075	−0.0008	0.0037	−0.0079	0.0036
$n=150$
$a$	0.55	0.5057	0.5140	0.2357	0.0881	0.9013	−0.0443	0.0348	−0.1602	0.1297
$b$	2	2.5319	2.4617	0.6610	1.5058	3.9007	0.5319	0.6970	0.5649	0.5069
$c$	0.6	0.5195	0.5049	0.1110	0.3509	0.7512	−0.0805	0.0149	−0.0606	0.0389
$\alpha$	0.4	0.3952	0.3942	0.0947	0.2157	0.5833	−0.0048	0.0110	−0.0240	0.0117
$n=300$
$a$	0.55	0.5231	0.5342	0.2298	0.1078	0.8998	−0.0269	0.0293	−0.1473	0.1140
$b$	2	2.5127	2.4588	0.5293	1.6697	3.6001	0.5127	0.6531	0.4822	1.1736
$c$	0.6	0.5350	0.5238	0.1015	0.3789	0.7410	−0.0650	0.0128	−0.0349	0.0274
$\alpha$	0.4	0.3996	0.3994	0.0364	0.3298	0.4706	−0.0004	0.0014	−0.0041	0.0020
$n=700$
$a$	0.55	0.5352	0.5061	0.0949	0.3861	0.7249	−0.0148	0.0184	−0.1197	0.1031
$b$	2	2.4024	2.3495	0.4620	1.6729	3.3589	0.4024	0.3889	0.4058	0.9847
$c$	0.6	0.5353	0.5261	0.0949	0.3861	0.7249	−0.0647	0.0117	−0.0350	0.0231
$\alpha$	0.4	0.3999	0.3996	0.4291	0.3176	0.4837	−0.0004	0.0020	−0.0034	0.0014

5. 实例分析

通过对1877年12月1日至1879年4月15日期间每周太阳黑子群的数量数据应用说明了所提出的方法，采用雅典的SCHMIDT独立观测的数据，其太阳黑子群的数量数据由91个观测值组成。样本均值和样本方差分别为1.154和4.651，表现出相当大的过分散。图3中描述太阳黑子群的数量数据的序列图、直方图、ACF图和PACF图。其中ACF图与PACF图表明数据存在显著的自相关。且从计数序列图与序列直方图中可以看出数据具有拖尾性质。数据可由美国国家地理数据中心 (http://www.ncei.noaa.gov/)获得。

Figure 3. Sequence, histogram, ACF and PACF plots of the sunspot subgroups data

图3. 太阳黑子群的数量数据的序列图、直方图、ACF图与PACF图

本文采用以下三个模型拟合该数据，其拟合结果进行比较：

(1) PINAR(1)模型[30]；

(2) PIGINAR(1)模型[6]；

(3) GPIGINAR(1)模型[31]。

本文运用MCMC算法对三个模型进行贝叶斯估计，经30,000次迭代后，舍弃前15,000次，进而得到最终估计结果。表3中为GPIGINAR(1)模型与其他模型的比较结果，包括参数估计结果、期望、方差、$I_d$ 、AIC与BIC的值。从表3结果可知，PINAR(1)模型表现欠佳，离散度最低。PINAR(1)模型和GPIGINAR(1)模型的均值与样本均值相近，然而，仅GPIGINAR(1)模型的方差与样本方差接近，且其离散度也最接近样本离散度。在所有模型中，GPIGINAR(1)模型的AIC和BIC值最小。综合各项指标，GPIGINAR(1)模型最适宜模拟该数据。

为了检验所讨论拟合的模型是否真的适合分析的数据。检验模型充分性最常见的方法是比较拟合模型的均值和离散特性，采用标准的皮尔逊残差检验(Jung等(2006) [32])：

$Z_t=\frac{Y_t-E\left(Y_t|Y_{t-1}\right)}{\sqrt{Var\left(Y_t|Y_{t-1}\right)}}$ (15)

根据GPIGINAR(1)模型的条件期望与条件方差公式，计算了在贝叶斯估计方法下残差的均值和方差，认为拟合残差的均值和方差越接近0和1，模型的拟合效果越好。在本节的实际数据中，所拟合的贝叶斯方法GPIGINAR(1)模型残差的均值和方差分别为−0.0013和1.0111，该结果也再次说明应用GPIGINAR(1)模型拟合该数据是合适的。

Table 3. Comparison results

表3. 比较结果

模型	参数	估计值	期望	方差	$I_d$	AIC	BIC
PINAR (1)			1.0608	1.0608	1	308.8153	313.6287
	$\widehat{\alpha }$	0.210
	$\widehat{\lambda }$	0.838
PIGINAR (1)			0.951	11.8474	12.4578	299.0965	306.3167
	$\widehat{\alpha }$	0.109
	$\widehat{\mu }$	0.951
	$\widehat{\varphi }$	0.083
GPIGINAR (1)			1.0597	5.4371	5.1308	194.9537	203.609
	$\widehat{\alpha }$	0.1064
	$\widehat{a}$	0.1841
	$\widehat{b}$	0.8633
	$\widehat{c}$	0.1122

为了形象的展示残差的性质，基于拟合的GPIGINAR (1)模型绘制了残差序列图、直方图、ACF及PACF图，如图4所示。为了进一步验证残差序列是一个白噪声序列，采用Ljung-Box (LB)检验的方法对残差$Z_t$ 进行检验，对其延迟从1到12，对于每一个延迟，计算LB统计量和其$p$ 值，其结果在表4中，可以看到其$p$ 值均大于0.05，从而说明该残差序列是一个白噪声序列。进一步说明GPIGINAR(1)模型拟合这组数据是合适的。

Figure 4. Residual sequence, histograms, ACF and PACF plots

图4. 残差序列图、直方图、ACF及PACF图

Table 4. LB test results

表4. LB检验结果

	$\left\{{\widehat{Z}}_t\right\}$		$\left\{{\widehat{Z}}_t^2\right\}$
延迟阶数	LB统计量	$p$ 值	LB统计量	$p$ 值
1	0.5057	0.4769	3.8224	0.0506
2	0.9343	0.6268	4.8255	0.0896
3	1.2958	0.7301	5.7183	0.1262
4	1.4084	0.8427	5.7979	0.2148
5	1.4601	0.9176	6.6746	0.2460
6	2.4681	0.8720	7.6963	0.2612
7	3.0047	0.8846	8.6061	0.2822
8	3.0646	0.9302	8.6208	0.3753
9	3.4702	0.9427	8.8396	0.4522
10	4.2350	0.9361	9.1885	0.5143
11	4.2904	0.9606	9.1936	0.6040
12	4.3733	0.9757	9.4592	0.6633

为了说明GPIGINAR(1)模型的马尔可夫链收敛，如图5所示，通过拟合模型的每个参数的轨迹图、直方图以及累积均值图，研究了马尔可夫链的收敛诊断。通过观察轨迹图，没有明显的上升或下降趋势，也不存在周期性等规律变化，从整体上可以看到每个参数随着时间的推移逐渐收敛到稳定状态，认为马尔可夫链已经收敛。通过观察直方图，可以知道每个参数都已经收敛到目标分布，认为其马尔可夫链已经收敛。通过观察累积均值图，可以看到曲线在某个值附近波动较小，说明参数的均值已经稳定，暗示MCMC算法已经收敛。

Figure 5. Fitting data convergence plots, histogram plots and ergodic mean plots

图5. 拟合数据链收敛图、后验密度直方图、累计均值图

6. 结论

在对GPIGINAR(1)模型进行研究时，基于MCMC算法开展贝叶斯估计工作。具体而言，通过在似然函数中巧妙引入潜变量，成功得到完全数据似然，进而在此基础上推导出参数的满条件分布。对于那些具有复杂后验分布的参数，采用随机游走加独立核的方法进行采样操作，以此获得参数的后验样本。随后，通过数值模拟的方式，将该方法与极大似然估计方法进行对比，结果充分验证了贝叶斯方法的优越性。此外，以雅典SCHMIDT观测的太阳黑子群数据作为实例，详细介绍了该模型在实际中的应用过程。

基金项目

辽宁科技大学博士启动资金(6003000310)。

NOTES

^*通讯作者。