王超华,李国东
(新疆财经大学 应用数学学院,新疆 乌鲁木齐 830012)
由于影响一个地区降水量的因素多,随机性大,所以有关降水量的研究备受国内外学者的关注。随着对降水量研究的不断深入,在降水量的研究方法上也发生了很多改变,从单一的使用某一种方法,到将多种方法相结合。同时也从更多地关注对降水量的分析上转到了更多地对降水量的分析与预测上。如近两年在预测降水量方面上,黄玺玮等人利用随机森林方法进行预测[1];李吉印利用BP神经网络马尔科夫模型进行预测[2];李亚斌等人利用加权马尔科夫链方法预测[3];Hossein等人利用外生变量改进人工神经网络在干旱和极端干旱气候条件下的降水预报[4];Samra等人利用卡尔曼滤波器对多站点年降水量进行预测[5]。这些方法都是采用多模式集合的方法,不同程度的降低了运算难度、提高了预报灵活性。在2006-2012年利用时间序列相关方法及与其他方法相结合对降水量进行预测的相对较多,刘贤赵等人采用 ARIMA(autoregressive integrated moving average)模型对历年降水量动态数据进行分析[6];王喜华等人采用小波分析与ARMA-GARCH模型相结合的方法预报降水[7];李可柏等人将时间序列频域谱分析和时域ARIMA方法结合对降水量数据进行拟合预测[8];迟道才等人采用自回归移动平均模型(auto regressive moving aver age,ARMA)预测[9]。
时间序列是同一种现象在不同时期相继观察排列而组成的一组数列,其分析方法只需序列本身的历史数据。ARMA模型是时间序列分析中最基本、应用很广的时序模型,它不需要事先假定数据存在一定的结构,只是从数据本身出发,去寻找一种能够较好描述数据的特征的模型,但如果只考虑降水量随时间变化的特征,忽略降水量时间序列的方差随时间变化的特征,往往达不到较为理想的效果。为此,本文选取从数据共享中国气象局网站获得的沙雅县从1981-2010年累年20-20时日降水量共365个数据进行分析。通过考虑降水量时间序列方差随时间变化的特征,利用传统的移动平均模型(moving average,MA)与广义自回归条件异方差模型(generalized auto regressive conditional heteroskedastic,GARCH)相结合的方法,建立MA-GARCH模型分析降水量的波动情况并对降水量进行预测。但在建立GARCH模型后发现序列并不满足GARCH模型的约束条件,最终建立了MA-EGARCH模型,并对降水量进行预测,发现该模型能够较好的描述降水量时间序列的特征,且预测精度较好。
利用时间序列数据建立ARMA模型,要求序列是平稳的,但实际的时间序列并不都表现出平稳性,因此首先利用Eviews8.0软件对沙雅县累年降水量进行数据平稳性检验,仅从时间序列趋势图,如图1,不能确定降水量数据是否平稳,可利用残差自相关图以及单位根检验来确定时间序列的平稳性。
序列的自相关函数可用来判断序列的平稳性,若一个序列是平稳的,则它的自相关函数随滞后阶数的增加而迅速趋于0。从图2降水量序列的自相关函数(auto correlation,AC)随着滞后阶数的增大(滞后10阶),并没有迅速趋于0,只是缓慢的下降,说明该序列不平稳。对其进行一阶差分,差分后序列记作jsl1。由图3可知,其自相关函数在滞后3阶后趋于0,说明降水量序列在经过一阶差分后平稳,降水量一阶差分后结果如图4所示。
图1 沙雅县累年降水量波动图
图2 序列的自相关和偏自相关系数
图3 序列一阶差分后的自相关和偏自相关系数
图4 降水量一阶差分后波动图
常用的单位根检验有ADF检验和DF检验。由于DF检验只适合检验一阶自相关的情况,故本文选用ADF检验。基于AR(P)模型:
将(1)式经过等价变化可得,
ADF检验统计量:
ADF单位根检验方法,如果计算出的统计量值小于临界值,则拒绝原假设,该过程是平稳过程;否则不能拒绝原假设,则该过程是非平稳的。
从表1可知,原降水量序列在5%的显著水平下不能拒绝存在单位根的原假设,说明序列不平稳。从表2可知,一阶差分后降水量序列无常数项(c)、无时间趋势项(t)时,已经拒绝一阶差分后降水量序列存在单位根的原假设,即降水量序列为平稳序列。
在确定序列为平稳序列后,建立均值模型。利用一阶差分后序列的自相关函数(AC)和偏自相关函数(PAC)识别模型阶数。从图4结果可认为降水量序列的自相关函数在滞后一阶后表现为“截尾性”,而偏自相关函数在滞后一阶后显著不为0,表现为“拖尾性”,所以可以考虑MA(1)、MA(2)移动平均降水量模型。MA(1)模型形式如(4)式,MA(2)模型形式如(5)式。
其中,Yt表示降水量,Yt-1,Yt-2分别表示滞后1阶、2阶的降水量,εt表示残差项。
表1 降水量序列的单位根检验
表2 原降水量序列的一阶差分
分别对上述MA(1)、MA(2)降水量模型采用最小二乘法进行拟合,结果如表3所示。从表3的检验结果来看,AIC选择了MA(2),但BIC选择了更为简洁的MA(1),这里不能明显地判断出哪个模型更优,因为εt-2的系数P值为0.2178。但现在无须选择哪个模型更恰当,因为每一个模型都是在假设方差为常数的条件下估计的。若实际的方差是随时间变动的,那么同时估计均值模型和方差模型时,系数的标准差就会发生很大变化[10]。
当模型存在异方差时,仍用最小二乘估计未知参数,参数的估计仍是无偏的,但对于标准差的估计是错误的,所以推断会产生错误。当使用时时间序列数据时,一般假设方差不会随着的变化而变化,即同方差。但当我们关心预测精度时,需要了解误差方差的大小,需要知道随时间变化方差的变化情况,即异方差。波动率建模是对扰动的方差建模,不仅可以修正错误的方差,改进参数估计的有效性,改进预测置信区间的精度还可以预测误差方差的大小。
表3 两模型检验结果比较
-0.064745 0.2178
通过利用Lung-Box Q统计量检验残差平方序列的自相关性,以此来判断均值方程的残差是否存在ARCH效应。
1.假设条件
H0:残差平方序列纯随机
H1:残差平方序列具有自相关性
2.检验统计量
其中计算估计的{ε2t}的样本自相关系数公式为:
其中,T为样本总量,均值模型中误差项方差为:
Lung-Box Q统计量近似服从自由度为n的χ2分布
从MA(1)、MA(2)估计模型得出的残差平方的自相关函数图来看,如图5所示。两模型的Q统计量分别为60.207、64.335,大于在5%显著水平下自由度为17的χ2(1-α)(n)的临界值27.587,所以说明两个均值方程的残差序列存在自相关,即残差序列存在ARCH效应。
当检验到两均值模型的随机误差项存在ARCH效应时,可使用ARCH和广义的GARCH去拟合随机误差项的条件方差。由于ARCH模型描述序列变化所需较大q阶时,是极不方便的,且其经常存在违反系数非负的约束等问题。所以本文选择GARCH模型。在进行GARCH模型估计前本文分析了降水量序列的分布特征,以此来判断该序列服从什么分布,为更准确的描述数据特征做准备。
图5 两模型残差平方图
从图6降水量的统计特征可知,序列JSL的偏度S=2.83476>0,峰度K=13.31173>3,因此,与正态分布相比,沙雅县的降水量呈现“右偏、高瘦”的分布形态。同时Jarque-Bera=2105.976,其概率P值=0.000,说明在5%的显著水平下拒绝接受该序列服从正态分布的原假设。为了更精准的描述这些时间序列分布的“厚尾”特征,需要对误差项的分布进行假设,GARCH模型的随机误差项的分布形式一般有3中假设:正态分布、学生t分布和广义误差分布(GDP)。由于这里分析出降水量序列不服从正态分布,所以使用学生t分布或广义误差分布(GDP)能够更好地描述厚尾特征。这里经过比较后发现使用学生t分布能够较好地描述数据的特征。因此本文选用学生t分布与GARCH模型一起估计参数。
GARCH(p,q)模型的一般形式:
图6 降水量的统计特征
其中,vt为白噪声过程,p是自回归GARCH项的阶数,q是 ARCH 项的阶数。α0、εi、βj>0,i,j=1,2,3……
GARCH(p,q)过程是稳定过程的充要条件为:
利用Eviews8.0,分别对存在异方差的MA(1)、MA(2)进行 GARCH(1,1)模型估计,估计结果如表4-5所示。
根据表4所示结果可知,在条件异方差方程GARCH(1,1)中参数估计得z统计量非常显著,其相应概率值P非常小,说明这些参数估计值都是显著的,且参数值都大于0,这就保证了条件异方差的非负的要求,符合GARCH模型的要求。ARCH项和GARCH项的系数估计值α1,β1分别为0.22778,0.82333,α1+β1=1.05111>1,不能满足GARCH模型参数约束条件,所以不能表明随机误差的条件方差能够收敛到无条件方差σ2。
从表5中也不难看出,虽然GARCH模型中各参数估计值都显著,但是α1+β1=1.02306>1,所以也不能满足GARCH模型参数约束条件,因此表明GARCH(1,1)模型是不平稳的。
由于标准GARCH模型必须保证所有估计的系数都为正且系数之和小于1。但在本文估计的GARCH模型中虽估计的系数都为正,但系数之和大于1。所以考虑考虑GARCH模型的指数GARCH模型,即EGARCH模型。该模型不考虑估计系数非负以及系数之和小于1的限制。条件异方差EGARCH(1,1)的形式如下:
表4 MA(1)GARCH(1,1)参数估计
表5 MA(2)GARCH(1,1)参数估计
其中,条件方差为线性对数形式。同时,EGARCH模型使用标准化的εt-1值(即εt-1/(ht-1)1/2),标准化的εt-1值不受单位的限制,可以更自然的解释冲击的大小和持续性。
对一阶差分后的降水量序列进行EGARCH(1,1)参数估计。由估计结果可知,在MA(1)-EGARCH(1,1)模型中,EGARCH(1,1)模型中的常数项α0的P值为0.1102,在5%显著水平下不显著,系数λ1的P值在5%显著水平下显著。
相比于MA(1)-EGARCH(1,1)模型,MA(2)-EGARCH(1,1)模型中εt-1和εt-2系数分别为-0.6213和-0.0811,εt-1的P值为 0.0692,在5% 的显著水平下不显著。且在EGARCH(1,1)模型中常数项以及系数λ1的P值分别为0.2102、0.0836,在5%的显著水平下都不显著。这时可考虑将MA(2)-EGARCH(1,1)调整为 MA(1)-EGARCH(1,1)。
此外,MA(1)-EGARCH(1,1)模型的对数似然值110.7534大于MA(2)-EGARCH(1,1)模型的对数似然值110.745,MA(1)-EGARCH(1,1)模型的AIC信息准则为-0.7007稍小于MA(2)-EGARCH(1,1)模型的AIC信息准则,其值为-0.6940,MA(1)-EGARCH(1,1)模型的BIC准则为-0.6264也小于MA(2)-EGARCH(1,1)模型的BIC准则,其值为-0.6073。模型的对数似然值越大说明模型越好,AIC与BIC越小越好。综合以上的分析,可以选择MA(1)-EGARCH(1,1)进行进一步的分析。
从上述分析可知,与MA(2)-EHARCH(1,1)相比,Ma(1)-EGARCH(1,1)更能拟合序列的特征,所以对Ma(1)-EGARCH(1,1)模型进行适应性检验。模型恰当程度检验是检验估计的{vt}是否为高斯白噪声,若估计的{vt}是高斯白噪声,即估计的{vt}不存在自相关,则说明均值模型没有问题,否则应适当改变均值模型。若估计的{v2t}不存在自相关,则说明条件异方差模型没有问题,否则应调整条件异方差模型。利用Ljung-Box Q统计量分别对估计的{vt}、估计的{v2t}进行检验。MA(1)在0.05的显著水平下估计的{vt}的Q统计量等于15.002小于χ2(17)的临界值27.587,而估计的{v2t}的Q统计量等于2.0713小于临界值27.587,可知估计的{vt}与估计的{v2t}均不存在自相关,因此所建的均值模型和条件异方差模型都是合适的。
由结果可以得出 MA(1)-EGARCH(1,1)模型的估计结果如下:
降水量序列的均值模型:
Z统计量=-12.32373。其中,假定误差项服从学生t分布。条件异方差方程:
Z统计量=1.5972883.739658 -2.055089389.9581。
利用 MA(1)-EGARCH(1,1)模型分别进行短期预测(第151天到160天,共10天。结果如表6所示)和长期预测(第151天到200天,共50天。结果如图7所示),其中中间线表示预测值,上下两条线为预测置信区间。
利用4种常用的误差评价指标评价预测模型的预测性能,其4种误差评价指标是均方误差(MSE)、平均绝对误差(MAE)、平均绝对百分比误差(MAPE)和均方百分比(MSPE)。
从表7预测模型的误差各项指标可看出,进行短期预测时各个误差评价指标均小于长期预测各项误差评价指标,说明此模型的短期预测效果比长期预测效果好。
本文利用MA-EGARCH模型对沙雅县降水量的动态数据进行了建模与预测。结论如下:
对降水量时间序列进行分析发现,降水量有波动聚集性的特点,其方差随时间的变化而变化。条件异方差模型与均值模型相结合的方法,不仅能反映出数据的部分规律,同时还提高了短期预测区间的精度,对短期预测降水量有一定的借鉴意义。
表6 短期预测结果
图7 长期预测结果
表7 预测模型的误差各项评价指标
从最后的预测结果来看,本文模型的拟合不一定是最佳结果。本文只是利用时间序列分析的方法对类似动态数据进行了尝试。考虑到有诸多因素影响一个地区的降水量,如纬度、气温、气压、植被覆盖率等,可以把诸多影响降水的因素作为模型的自变量,降水量作为因变量进行多变量的时序分析,从而得到更为精准的预测结果。