传染病异方差性时序数据的建模与拟合

【摘要】 目的 探讨异方差性时间序列模型在传染病疫情数据分析中的应用。方法 分别采用arima和argarch模型对某市淋病发病率月报数据进行建模和拟合。结果 本资料构成的时间序列经检验具有明显异方差性，经模型比较和筛选，ar(1)garch(0,1)模型能够较好的拟合本研究中传染病疫情时序数据。结论 argarch模型适用于传染病疫情数据构成的异方差性时序数据分析。

【关键词】 argarch模型; 时序数据; 异方差; 传染病

modeling and fitting for heteroscedastic timeseries data of infectious diseases zhang ziwu, li xiaosong, quan ting. department of statistics, school of public health, sichuan university, chengdu 610041, china

　　【abstract】 objective to explore the application of heteroscedastic time series model to the analysis of data of infectious diseases. methods arima and argarch models were used to fit the incidence of gonorrhea. results the time series in this study, which was heteroscedastic significantly, finally was well fitted by ar (1)garch (0,1) model through model selecting. conclusions argarch model is suitable for analyzing heteroscedastic timeseries data of infectious diseases.

　　【key words】 argarch model; timeseries data; heteroscedasticity; communicable diseases

　　传染病疫情资料蕴含和反映了传染病传播和流行的相关信息。传染病发病率资料是一组时序数据，运用时间序列分析技术可了解传染病的发展趋势并作出短期预测，有助于及时发现异常并采取措施[1]。目前，传染病疫情数据的时间序列分析方法主要是求和自回归移动平均(autoregressive integrated moving average,arima)模型[1～4]，该模型在拟合非平稳序列时，假定残差序列为零均值白噪声序列，但实际的传染病疫情资料常不满足这一假设，有的残差序列呈现出异方差性，arima模型并不能很好地拟合异方差性时间序列，这一问题在传染病疫情数据的分析和短期预测中并未得到足够关注。因此，有必要采用更为合理的统计模型对这类资料进行拟合。本文采用自回归-广义自回归条件异方差模型(autoregressivegeneralized autoregressive conditional heteroskedastic,argarch模型)对某市淋病发病率月报数据进行拟合，探讨对于具有异方差性的时序数据的建模及拟合方法。

　　1 模型基本原理

　　运用arima模型拟合非平稳序列时，对残差序列{εt｝假定为零均值白噪声序列，即：零均值，e(εt)=0;纯随机性，cov(εt,εt-i=0，i≥1;方差齐性var(εt)=σ2ε。如不满足方差齐性，即残差序列方差不为常数，而随时间变化而变化，表示为时间的函数var(εt)=h(t)，称为异方差。如忽视其存在会导致残差方差严重低估，增加ⅱ型错误，影响模型拟合精度[5]。

　　针对具有异方差性的时间序列，engle于1982年提出条件异方差(autoregressive conditional heteroskedastic,arch)模型，其结构为：

　　xt=f(t,xt-1,xt-2,…)+εt

　　εt=htet

　　ht=ω+∑pj=1λtε2t-j

　　式中，f(t,xt-1,xt-2,…)为{xt｝的自回归模型;et～i,j,dn(0，1);ω、λi为参数。在arch模型基础上，bollerslov在1985年提出广义自回归条件异方差(generalized autoregressive conditional heteroskedustic,garch)模型，增加考虑了异方差函数的p阶自相关性，可有效拟合具有长期记忆性的异方差函数。其结构为：

　　xt=f(t,xt-1,xt-2,…)+εt

　　εt=htet

　　ht=ω+∑pj=1ηtht-i+∑qj=1λiε2t-j

　　式中，f(t,xt-1,xt-2,…)为{xt｝的回归函数;et～i,j,dn(0，1);ω、ηi、λi为参数。

　　argarch模型是garch模型的变体，它考虑了如果回归函数f(t,xt-1,xt-2,…)不能充分提取残差序列{εt｝信息，造成{εt｝不满足纯随机性，应先对{εt｝拟合自回归模型，再考虑自回归残差序列{υt｝的方差齐性，如为异方差，则对它拟合garch模型，称为ar(m)garch(p,q)模型。其结构为：

　　xt=f(t,xt-1,xt-2,…)+εt

　　εt=∑mk=1βkεt-k+υt

　　υt=htet

　　ht=ω+∑pj=1ηtht-i+∑qj=1λiυ2t-j

　　式中，f(t,xt-1,xt-2,…)为{xt｝的回归函数;et～i,j,dn(0，1);β、ω、ηi、λi为参数。

　　2 实例分析

　　根据某市1998-2007年淋病发病人数月报资料及相关人口资料，计算得该市1998-2007年淋病逐月发病率。采用sas 9.0软件编程[6]进行数据分析。首先对该资料拟合arima模型，通过模型筛选，最终选取arima(0,1,1)模型进行拟合，条件最小二乘估计拟合结果为

　　(1-b)xt=(1-0.634 87b)εt

　　var(εt)=67.898 3

　　对残差序列检验结果见表1，残差序列检验p值均大于0.05，可认为残差序列为白噪声序列，拟合系数θ1经参数检验p<0.001，有统计学意义，说明arima(0,1,1)模型可用于拟合本资料。　表1 残差序列检验结果

　　进一步观察原序列1阶差分后的残差平方图(图1)，发现差分后残差序列的方差呈递减趋势，呈异方差性。

　　图1 1阶差分残差平方图

　　figure 1 square of residual after first order differece

　　由时序图可见，数据呈自相关特征，构造自相关函数：xt=α1xt-1+εt，对残差序列{εt｝进行自相关检验，检验结果显示，durbin h=-4.478 0，p=0.000 1，表明残差序列自相关性具有统计学意义。由残差序列5阶延迟自相关图可见其具有1阶显著自相关性，如图2所示，因此对残差序列进一步拟合自回归模型εt=β1εt-1+υt。

　　lagcovariancecorrelation-1987654321012345678910 82.374 7 1.000 000 | |*************|1-22.560 6-0.273 878 | ****| |214.102 40.171 198 | |*** |312.248 40.148 691 | |** |46.976 10.084 687 | |** |57.052 40.085 614 | |** |

　　图2 残差序列自相关图

　　figure 2 autocorrelations of residual

　　表2 异方差检验结果

　　table 2 q and lm tests for arch disturbances

　　orderportmanteau qplmp114.643 5 0.000 1〖〗13.926 00.000 2216.753 80.000 213.934 30.000 9316.866 90.000 814.779 90.002 0416.885 30.002 014.999 80.004 7516.902 80.004 715.092 40.010 0616.906 80.009 615.189 00.018 8729.366 10.000 127.606 40.000 3839.373 2<0.000 127.849 00.000 5943.385 3<0.000 127.933 70.001 01043.453 7<0.000 128.827 40.001 31144.468 8<0.000 130.131 60.001 51246.575 2<0.000 130.273 20.002 5

　　对新产生的残差序列{υt｝进行异方差检验，portmanteau q统计量和拉格朗日乘子(lm)检验统计量均显示{υt｝异方差自相关性具有统计学意义，如表2所示。因此，尝试拟合garch(p,q)模型ht=ω+∑pj=1ηtht-i+∑qj=1λiυ2t-j。

　　使用极大似然估计方法，采用迭代技术对未知参数β、ω、ηi、λi进行估计，按模型参数估计值有意义及施瓦茨贝叶斯(sbc)、赤池信息准则(aic)最小为原则筛选模型，最终确定ar(1)garch(0,1)模型，拟合参数结果如表3所示。

　　表3 ar(1)garch(0,1)模型参数拟合结果

　　table 3 parameters of ar(1)garch(0,1) model

　　parametersdfestimatestandard errortplaga1 1.006 40.009 0113.20<0.000 1ar110.464 30.051 3 9.05<0.000 1arch0126.496 16.054 3 4.38<0.000 1arch110.812 50.229 3 3.54 0.000 4

　　如模型拟合准确，则函数υt/ht呈现标准正态分布。正态性检验tn=5.041 2，p=0.080 4(>0.05)，接受υt/ht服从标准正态分布。最终模型r2=0.954 0，拟合结果为：

　　xt=1.006 4xt-1+εt-1+εt

　　εt=0.464 3εt-k+υt

　　υt=htet，et□i,j,dn(0，73.139 0)

　　ht=26.496 1+0.812 5υ2t-1

　　比较arima(0,1,1)模型和ar(1)garch(0,1)模型的拟合优度，如表4所示，可见ar(1)garch(0,1)模型的aic、sbc均小于arima(0,1,1)模型，前者优于后者;观察拟合效果图，可见ar(1)garch(0,1)模型拟合较优。见图3。

　　表4 拟合优度比较

　　table 4 comparison of goodness of fit

　　modelsaicsbcarima(0,1,1) model833.59836.36ar(1)garch(0,1) model818.03829.11

　　图3 ar(1)garch(0,1)模型拟合效果图

　　figure 3 fitting results of ar(1)garch(0,1) model

　　3 讨论

　　影响传染病传播和流行的因素纷繁复杂，不同病种影响因素差异也很大，而且通常难以获取全面的影响因素的资料，因此对未来传染病发病率进行预测十分困难。正因为如此，时间序列分析被广泛应用于传染病发病率的预测，其过程简单适用，精度较高，适用于仅有发病率资料的传染病短期预测。利用预测数据的可信限在疫情监测中可判断发病率是否在正常范围波动[3]。

　　随着社会经济的发展，近年我国一些传染病发病率数据出现了明显的变化趋势。有的传染病发病率呈迅速增高的趋势，如hiv，梅毒等，这些传染病发病率资料构成的时间序列常具有异方差的特性;此外，对突发性传染病如“非典”，或传染病暴发，因发病率数据快速升高或降低，一般存在一定趋势且波动幅度较大，其构成的时间序列多呈异方差性。对于以上具有异方差性的时间序列单纯采用arima模型拟合并作出预测存在精度较低的问题。

　　argarch模型已广泛应用于计量经济学领域中，并取得良好的拟合和预测效果，但在我国公共卫生领域特别是传染病流行病学研究中应用较少。本文结果显示，应用argarch模型对淋病发病率月报资料拟合较好，通过与相同资料拟合的arima模型比较，argarch模型的信息准则sbc、aic均小于arima模型，拟合优度更佳，说明argarch模型对于具有异方差性的时间序列数据的拟合具有明显优势。

【参考文献】
　　[1] 冯丹,韩晓娜,赵文娟,等. 中国内地法定报告传染病预测和监测的arima模型 [j]. 疾病控制杂志, 2007,11(2):140143.

　　[2] 冯超,白衫. 时间序列模型拟合艾滋病发病趋势预测 [j]. 中国公共卫生, 2005,21(7):893.

　　[3] 吴家兵,叶临湘,尤尔科. arima模型在传染病发病率预测中的应用 [j]. 数理医药学杂志, 2007,20(1):9092.

　　[4] 何庭辉,柳青. 广州市荔湾区恶性肿瘤死亡率时间序列预测模型分析 [j]. 疾病控制杂志, 2006,10(6):630632.

　　[5] 王燕. 应用时间序列分析 [m]. 北京：中国人民大学出版社, 2005.175.