NGINAR(1)模型的Bayes估计及预测

杨艳秋,王德辉

(1.吉林大学 数学学院,长春 130012;2.吉林师范大学 数学学院,吉林 四平 136000)

现实生活中的许多计数过程,如某医院某时刻住院的病患人数、某时刻某地区的生物种群数量、某地区某阶段的犯罪数量等,这些数据通常取非负整数值,用一般的时间序列模型拟合这些数据通常会产生异常预测,因此需要引入整值自回归模型.目前对整值时间序列数据的建模主要有如下两种方法: 借助于潜过程的状态空间建模过程[1]；借助于稀疏算子的建模过程,这类方法是建模整值时间序列数据的主要方法.二项稀疏算子“∘”[2]是整值时间序列发展的基础,定义为

其中:α∈(0,1);X为非负整值随机变量;{Bi}为独立同分布(i.i.d.)的Bernoulli随机变量序列且与X互相独立,满足

P(Bi=1)=1-P(Bi=0)=α.

但二项稀疏算子有一个局限,即求和序列{Bi}为i.i.d.的Bernoulli随机变量序列,只能取0或1值,因此α∘X的取值总是小于或等于X的值.但在实际问题中,每一个事件都可能关联更多相关的计数事件,因此用几何随机变量来描述这些事件更合适.文献[3]引入了负二项稀疏算子“*”,定义为

k∈0.由于负二项稀疏算子中求和序列{Wj}的取值是非负整数,使得β*X的取值可能大于也可能小于X的取值,从而很好地突破了二项稀疏算子的局限.本文给出基于负二项稀疏算子的一阶整值自回归模型参数的Bayes估计,先进行数值模拟,再与条件最小二乘估计和Yule-Walker估计进行均方误差的比较,最后对新西兰牛皮肤病数据进行实例分析及模型预测.

1 模型简介

基于二项稀疏算子“∘”的一阶整值自回归模型[4-5]为

Xt=α∘Xt-1+Zt.

若Xt表示t时刻住院的患者数,α∘Xt-1表示上个月以概率α仍继续住院的患者数,Zt为t时刻新住院的患者数,则Xt为α∘Xt-1与Zt之和.但对于过度分散的计数过程,负二项稀疏算子“*”则更合适.文献[6]基于负二项稀疏算子“*”,提出了一阶整数值自回归模型(new geometric first-order integer-valued autoregressive process,NGINAR(1))如下:

Xt=α*Xt-1+εt,

其中: 负二项稀疏算子“*”定义为

2 Bayes估计

其中

i=1,2,…,T.

考虑α的先验分布取(0,1)上的均匀分布,即π(α)=1,0<α<1.根据Bayes原理,参数α的后验分布为

综上可得：

定理1设样本x0,x1,…,xT来自于NGINAR(1)模型,在二次损失函数和均匀先验分布下,参数α的Bayes估计形如式(1).

3 数值模拟

下面通过数值模拟给出NGINAR(1)模型参数Bayes估计的优良性,将NGINAR(1)模型参数的Bayes估计与条件最小二乘(CLS)估计和Yule-Walker估计(Y-W)进行均方误差的比较.先给出Bayes估计的算法:

1) 选择迭代初值α(0),并令i=1;

5) 令i=i+1,返回步骤2),直到算法达到事先约定的收敛标准.

下面进行数值模拟,样本容量分别取T=100,500.取μ=0.5,5,10,表1列出了Bayes估计值的偏差(Bias)和均方误差(MSE),表2列出了参数的Bayes估计与条件最小二乘估计和Yule-Walker估计的均方误差比.模拟运行过程中先进行2 000次预迭代,以确保参数的收敛性,然后再进行500次迭代,得到模拟结果.

表1 Bayes估计的偏差和均方误差

表2 3种估计方法的均方误差比

由表1可见,Bayes估计的偏差和均方误差都比较小.以均方误差为准则,表2中3种估计方法的均方误差比可见,Bayes估计优于条件最小二乘估计和Yule-Walker估计.

4 实例分析及模型预测

图1 牛皮肤病数据样本路径Fig.1 Sample path of skin-lesions data

下面将NGINAR(1)模型应用到一组牛皮肤病患病数据集中[9],并进行分析及预测.该数据集来源于新西兰农林部,记录了新西兰某地区动物卫生实验室记录的2003—2009年间每月患皮肤病的牛数量.将该数据集分为两部分：2003-01—2009-08的数据用于估计参数值,2009-09—2009-12的数据作为样本外待预测值.

将数据集的前80个观测数据记作X1,X2,…,X80,统计结果表明,牛皮肤病数据的样本均值为1.5,样本方差为3.417 72,方差比均值为Id=2.278 5.牛皮肤病数据样本路径如图1所示,自相关函数图像及偏自相关函数图像分别如图2和图3所示.由图2和图3可见,该组数据为一阶相关,因此可以用NGINAR(1)模型对其进行拟合.

图2 牛皮肤病数据自相关函数图像Fig.2 Autocorrelation function plot of skin-lesions data

图3 牛皮肤病数据偏自相关函数图像Fig.3 Partial autocorrelation function plot of skin-lesions data

下面通过序列{Xt}的近似h步(h∈0)预测条件分布方法对NGINAR(1)模型进行预测(简称条件分布预测).条件分布预测方法较传统条件期望预测方法更适用于整数值时间序列.虽然条件期望预测方法可以使预测值的均方误差最小,但当观察值和待预测值为整数值时,利用条件期望预测方法得到的预测值却很少取到整数值点.为了解决上述问题,文献[10]提出通过条件分布预测方法对整数值模型进行预测,用这种预测方法得到的预测值和整数值时间序列本身的状态空间一致,而且利用条件分布预测方法来计算条件中位数、条件均值及条件众数等点的预测,甚至预测值的置信区间都比较容易,能得到较理想的预测值.

由于NGINAR(1)过程具有Markov性,在给定Xn的条件下,Xn+h的条件分布(即Xn的条件预测分布)为

P(Xn+h=xn+h|Xn=xn)=[Ph]xn+h,xn,

其中转移概率为

下面考察NGINAR(1)模型的预测效果.利用牛皮肤病数据集对模型进行预测.将条件分布预测和条件期望预测方法进行比较,结果列于表3.由表3可见,当h=1,2,3,4时,使用条件分布预测方法得到的预测值均为0,与实际值相符,而条件期望预测方法的预测结果分别为1.412 7,1.470 4,1.478 5,1.479 7,与实际值0有一定的偏差,而且条件分布预测方法的预测均值绝对偏差为0,条件期望预测方法的均值绝对偏差为1.460 3,因此,条件分布预测的方法更适用于整数值时间序列.

表3 牛皮肤病数据的条件期望与条件分布预测结果比较

图4为牛皮肤病数据的h步条件预测分布图像.用NGINAR(1)模型、ZIPINAR(1)模型、PINAR(1)模型来拟合该组数据,并用AIC(Akaike信息准则)、BIC(Bayes信息准则)、均方根值(RMS)和方差比均值对上述3个模型进行比较,结果列于表4.由表4可见,NGINAR(1)模型具有最小的AIC值和BIC值,3个模型的均方根值相差不大,而NGINAR(1)模型的方差比均值为2.479 9,更接近于数据集自身的方差比均值2.278 5,这些数据均表明用NGINAR(1)模型拟合该组数据集较合适.

图4 牛皮肤病数据的h步条件预测分布Fig.4 h-Step conditional prediction distribution of skin-lesions data

模型估计AICBIC均方根值IdZIPINAR(1)^α=0.164 1276.875 6284.021 71.806 81.674 0^λ=2.031 1^ρ=0.386 3PINAR(1)^α=0.157 3293.339298.103 11.807 51.000 0^λ=1.256 7NGINAR(1)^α=0.140 0271.103275.8671.809 72.479 9^μ=1.479 9

通过上述模拟结果可知,NGINAR(1)模型参数的Bayes估计效果优于条件最小二乘估计和Yule-Walker估计,且条件分布预测方法比条件期望预测方法更适用于整数值时间序列.