2011年全国大学生统计建模大赛论文二等奖第五名

二重趋势时间序列的组合预测模型

广东外语外贸大学陈文聪、庄文钿、苏玲

摘要

时间序列预测是一个多学科交叉的研究领域，呈现趋势变动性和周期波动性的二重趋势时间序列更是经常出现在各种重要场合，建立有效模型对其准确预测尤为重要。本文基于二重趋势时间序列理论，以我国1992年-2010年各季度的国内生产总值（GDP）为基础，分别建立单一预测模型GM（1，1），BP神经网络和ARIMA，并采取中位移动平均法对数据进行分解，构建基于GM（1，1），BP神经网络和ARIMA的二重趋势时间序列组合预测模型，进而对模型进行检验，确定相对最优模型为二重趋势时间序列的组合预测模型。利用所建立模型对我国2009年-2010年的季度GDP做出预测并与实际情况进行对比分析。

本文总结了前人时间序列预测研究的不足，指出了以组合预测模型代替单一预测模型，基于确定性因素分解以二重趋势组合模型代替非确定性因素分解的一般组合模型，尝试拟合组合模型残差项减少信息丢失，基于非参数Friedman检验科学确认组合模型相对最优的重要性和有效性。

本文在前人研究成果的基础上，针对现实二重趋势时间序列预测模型研究的不足之处，提出作者自己的观点和想法，并将它们付诸于实践，希望可以对二重趋势时间序列组合模型用于实际预测做出一点贡献。

关键词：二重趋势时间序列；组合模型；国内生产总值（GDP）；预测

一、概述......................................................... 3

（一）时间序列研究概况 ....................................... 3 （二）本文研究的意义、亮点 ................................... 3 （三）本文的研究思路及建模准备 ............................... 4

1. 研究思路 ................................................. 4 2. 数据来源 ................................................. 5

二、单一预测模型的介绍及建模过程 ................................ 6

（一）灰色预测模型 ........................................... 6

1. 模型介绍 ................................................. 6 2. 建模过程 ................................................. 6 （二） BP神经网络 ............................................ 12

1. 模型介绍 ................................................ 12 2. 建模过程 ................................................ 12 （三） ARIMA模型 ............................................. 16

1. 模型介绍 ................................................ 16 2. 建模过程 ............................................... 17

三、组合预测模型的引入及建模准备 ............................... 25

（一）若干假设 ................................................. 25 （二）数据的分解 ............................................... 25 四、组合预测模型构建过程 ....................................... 26

（一）长期趋势项模型构建过程 ................................... 26

1. 选用构建长期趋势项模型的单一模型 ........................ 26 2. 长期趋势项的灰色G（1，1）模型构建 ....................... 26 （二）季节波动项模型构建过程 ................................... 29

1. 选用构建季节波动模型的单一模型 .......................... 29 2. 波动项的BP神经网络构建 ................................. 29 3. 波动项的ARIMA模型构建 ................................ 31 4.季节波动模型的最终构建过程 ............................... 32 （三）随机干扰项模型构建过程 ................................... 33 （四）组合预测模型构建最终过结果 ............................ 34 五、组合模型与单一模型预测精度对比分析 ......................... 35

（一）直接结果对比分析 ......................................... 35 （二）非参数Friedman检验 ...................................... 36 六、总结与展望 ................................................. 37

（一）全文总结 ................................................. 37 （二）研究展望 ................................................. 37 参考文献 .......................................................... 39 附录............................................................... 40 致谢............................................................... 56

一、概述

（一）

时间序列研究概况

时间序列分析方法就是从具有先后顺序的信息中提取有用信息,它是数理统计的一个重要分支.其实质是通过对历史数据的处理,寻找出前后数据之间的关系,建立关联模型,然后通过历史数据和所建立的关联模型来预测时间序列的未来值,从而把握事物的未来发展[1]。然而许多商业和经济时间序列并非呈现单一趋势，而是表现出趋势变动性和周期波动性二重特征，该类时间序列的准确预测，对各种业务部门的有效决策起着非常重要的作用[2]。笔者研究发现对于二重趋势时间序列前人研究方法大部分集中在单一模型上面，如最常见的自回归移动平均（ARIMA）模型，灰色GM（1，1）模型已经今年来比较热门的神经网络模型等。也有少部分研究人员开始尝试建立组合模型来进行来提高预测精度，如《利用灰色神经网络预测组合模型在二重趋势性季节型电力负荷的预测》(牛东等,2001)、《回归时间序列组合模型在石油企业费用预算编制中的应用》（胡江波，2006）、《基于GMDH和BP的组合预测模型在工业产品销售收入中的应用》 (王洪涛等,2010)等。这些组合预测方法在具体问题的应用中都表现出了较高的预测精度，使所探讨的实际问题得到了较好的解决。然而更加有针对性的结合二重趋势时间序列特征科学合理选择单一模型建模的组合预测研究少之又少，有必要结合时间序列趋势变动性和周期波动性二重特征选择恰当单一模型建立组合预测模型。

（二）本文研究的意义、亮点

笔者在研究时间序列尤其是二重趋势时间序列建模预测的相关文献过程中得到很多启示，其中宋仙磊等的《二重趋势时间序列的灰色组合预测模型》（2011）对笔者有较大启发，本文基于以上研究提出更加完善的二重趋势时间序列组合预测模型，并给出全部建模过程。研究亮点主要如下：

1. 数据：不同于以往研究只选择年度数据或者月度数据，本文选择季度GDP

数据进行研究，尝试拟补季度GDP预测研究的不足，并且本文研究期数多达76期，使研究结果更有精确；

2. 单一模型：本文对单一模型的选择精雕细琢，结合不同单一模型的预测优

点最终确定以GM（1，1）预测趋势项，以BP神经网络和ARIMA通过最小二乘法加权预测周期项，以自回归模型拟合残差项；

3. 组合模型：为了最大限度减少信息丢失，本文在建立周期项组合模型时采

取最小二乘法加权使信息量提取充分，并且在总得组合模型建立时尝试拟合残差自回归模型，确保较少信息遗漏；

4. 模型检验：对于每个单一模型的建模过程，笔者都科学严谨检验每个过程，

对GM（1，1）进行残差修正，对BP神经网络采取试凑法确定误差最小网络结构，对ARIMA进行严格单位根检验等，并且在组合模型与单一模型对比阶段引进基于非参的Friedman方法，科学确定组合模型是否相对最优。

（三）本文的研究思路及建模准备

1. 研究思路（1）研究大体框架

针对二重趋势时间序列的趋势变动性和周期波动性特征，提出一种组合预测模型。首先对二重趋势时间序列进行分解，得到趋势变动项、周期波动项和随即干扰项，用灰色G（1，1）对趋势项预测，再用BP 神经网络和ARIMA的组合模型对周期项预测，用乘积模型合成这两部分结果，再加上随机干扰项得到二重趋势时间序列的最终预测值。其基本流程如图1所示：

（2）二重趋势时间序列的组合预测模型构造

ˆBP(k)是BP神经网络对周期波动项的预测值，yˆARIMA(k)是ARIMA模型对设y

ˆS(k)为：周期波动项的预测值，则周期项的组合预测值y

y(k)aby

BP

(k)bARIMAy

ARIMA

(k)

注：其中

、

ARIMA

分别为BP神经网络模型和ARIMA模型在周期波动项预测模型中的系数的系数。

结合最优线性组合模型原理，利用样本期的实际值和各单项预测模型的拟合值进行线性回归，然后利用线性回归模型，以原方案的预测值作为外生变量进行外推预测[3]。其中，权重bBP、bARIMA的确定可以通过线性规划法，另外还可有遗传算法、量子搜索算法、粒子群算法等其他的寻优算法。

ˆGM(k)是GM灰色预测模型对趋势项的预测值，把趋势项和组合后的周期波y

ˆ(k)为：动项进行再次组合，则二重趋势时间序列组合预测值y

y(k)y(k)y

GM

(k)(k)

2. 数据来源

本文选取的GDP季度数据均来自国家统计数据库，为实际GDP季度数据。笔者共收集了1992年到2010年19年共76期数据，较大的样本量为我们进行时间序列分析提供了基本条件。本文对1992年至2010年度的76期季度国内生产总值数据进行了分析，只用前68期数据参与建模，再用2009~2010年8期季度的数据检验模型拟合效果（具体数据见附录1）。本文使用多种经济学软件，如MATLAB，Eviews， SAS软件对数据进行处理建模。无论是对数据整理、数据平稳性检验还是模型拟合，效果都很好。

二、单一预测模型的介绍及建模过程

（一）

灰色预测模型

1. 模型介绍

灰色模型（GreyModel，GM模型）是灰色系统理论的基本模型。G（1，1）模型，即单变量一阶灰色模型，通过对原始数据作一次累加生成，并建立微分方程模型，求得拟合曲线，用以对序列进行预测，其建模步骤如下：

（1）给定序列在t时刻的观测值： X={X(1),X(2),···,X(n)}

（2）作一阶累加生成1-AGO，得到新的数列： X={X(1),X(2),···,X(n)}，其中：X (3)构造一阶微分方程：

(1)

(0)

(1)

(k)X(0)(i),k1,2,...n。

i1

dX(1)

aX(1)u dt

利用最小二乘法求得参数a，u：

a

ˆ(BTB)1BTYN， a

u

1(1)(1)1[X(1)X(1)]2X(0)(2)1(1)(1)(0)X(3)1[X(2)X(3)]

,YN其中：B2。   X(0)(n)1

[X(1)(n1)X(1)(n)]1

2求得微分方程，即可得预测模型：

ˆ(1)(k1)[X(0)(1)u]eaku X

2. 建模过程（1）数据预处理

对原序列做时序图可以看出该序列明显呈现固定周期性变化，即具有“季节”

效应。由于灰色G（1，1）模型仅能较好地拟合时间序列的趋势性部分，而对于周期波动性，其预测精度则明显降低。为了提高模型拟合度，在模型建立前，对数据进行预处理，消除时间序列的“季节”效应。对于存在“季节”效应的时间序列，可以通过原序列除以相应的季节指数予以消除。

表1.运用按季平均求季节指数

如表1求得各季度季节指数，用原序列除以相应的季节指数即可得到消除季

节影响后的时间序列（详见附录2）（2）GM（1,1）模型建立及修正

该模型运用SAS8.0软件进行拟合预测操作，代码详见附录。首先，记原始时间序列为：

X(0)={X(0)(1),X(0)(2),X(0)(3),X(0)(4),······,X(0)(68)}，上标0表示原始时间序列。（SAS中记为xt）

a）构造累加生成列（SAS中记为yt）

X(1)(1)=X(0)(1)= 12672.85 X(1)(2)=X(1)(1)+ X(0)(2)=27729.05 X(1)(3)=X(1)(2)+ X(0)(3)=43475.41 ······

X(1)(68)=X(1)(67)+ X(0)(68)=4793643.35

b）构造数据矩阵B（SAS中记为zt）和数据向量Yn

12[X(1)(1)X(1)

(2)] 11(12672.852772912.05]) (1)(1)1B[X(2)X(3)] 1(27729.0543475.41） 22  12[X(1)(67)X(1)

(68)] 112(4610527.244793643.35）20200.95 135602.23 1 

 4702085.30 115056.20

Yn15746.36

 183116.11

c）计算参数矩阵

ˆ

为待估参数，利用最小二乘法求解，可得 111

  

ˆ(BTB)1BTYn

-0.034202

 

 17253.66

即α=-0.034202，μ=17253.66

d）得出预测模型（预测结果记为xp）

dX(1)

0.034202X(1)17253.66dt

X(0)(1)12672.85 X(0)(1)



-504463.4817



517136.3317

X(1)(k1)517136.3317e0.034202k504463.4817

e）残差检验

表2.GM（1,1）拟合结果残差检验

表2仅截取前10个数据（所以数据详见附录3）。由表2可以看出，相对误差序列中的相对误差较大，所以要对原模型进行残差修正，以提高精度。

f）利用残差对原模型进行修正

将残差记为：

e(0)={e(0)(1’), e(0)(2’),„, e(0)(67’)}={2936.84,2872.72,„,11146.49}

e(0)的累加生成序列为：

e(1)={e(1)(1’),e(1)(2’),„,e(1)(67’)}={29336.84,5809.46,„,346283.25} e(1)可建立相应的GM（1,1）模型，求得参数

e0.016111ˆe2778.27 

eˆ(1)(k1)[e(0)(1)则e

ekek]e175382.3741e0.016111172445.5341 ee

(k1)ˆ(1)(k1)]'2825.585429e0.016111对上述求导，得：[e

这样得到经过残差修正后的模型：

(k1)ˆ(1)(k1)517136.3317e0.034202k504463.4817(k1)2825.585429e0.016111X

1,k2

其中：(k1)

0,k2

表3.GM（1,1）残差模型拟合结果

其中表3为GM（1,1）残差修正模型的预测结果的前10个数据（所有数据详见附录2）。由表3可见，与GM（1,1）相比，修正模型的精度相对提高。

g）关联度检验

记关联系数为n，分析如表4：

表4.关联系数检验

由表4可得，关联系数的均值，即关联度为0.75。根据经验，当ρ=0.5时关联度大于0.6便满意。所以，关联度检验通过。

h）后验差检验

表5.后验差检验

由表5可得，原始序列标准差S1=45773.05，残差标准差S2=4290.74， C= S2/ S1= 0.093739438 S0=0.6745S1= 30873.92223

由于残差的取值范围为[14.5879387,25195.58]，则

ˆ00,20101.43 0

ˆ0S1 Pp00

所以，P=1，C

i）用模型结果预测2009-2010年的数据，预测结果乘以各季度的季节指数得到GDP预测值。

表6.2009-2010年预测结果



由表6预测结果得，平均相对误差较大，达到7.17%进而通过GM（1,1）残差模型进行修正。修正结果关联度检验、后验差检验均可通过，但是残差检验中虽降低了相对误差，但是某些季度预测结果仍不容乐观，可见单一的修正后GM（1，1）残差模型仍存在一定局限性。

（二）

BP神经网络

1. 模型介绍

目前，在众多神经网络中，误差反向传播(ErrorBackPrpoogration)网络由于其良好的逼近能力和成熟的训练方法而得到了最为广泛的应用[4]。BP神经网络模型有良好的适应和自学习能力、较强的抗干扰能力。这使得BP神经网络在时间序列预测中颇受关注。

BP神经网络模型对时间序列预测的求解步骤如下：

（1）网络结构的设计：BP神经网络结构直接影响到模型的精度，但至今尚无统

一的理论方法，本文主要通过经验和多次尝试建立网络结构。

（2）样本的预处理：对输入的训练数据进行归一化处理，使输入数据的取值在

（-1，1）或者（0，1）之间。

（3）BP网络的学习规则：BP网络的学习规则一般包括以下几个步骤：①初始化

各层间节点i 和j 之间的权值wij 和节点j 的阈值θj ；②随机选取输入和目标样本提供给网络；③确定传递函数f (x) ；④ 对第p 样本，计算各层节点的输出；⑤计算各层节点的误差信号；⑥反向传播；⑦计算误差。经过不断的学习训练，使误差趋向于最小。

（4）反归一化：学习结束，对输入和输出样本进行反归一化处理，即可得到预测值。

2. 建模过程

BP网络的设计主要包括输入层,隐层,输出层及各层之间的传输函数几个方面。

（1）输入层：

a) 将数据分为训练样本数据和测试数据

本文采用1992年至2010年各季度全国GDP数据共76期。取其中前68期作为历史数据样本，将前60期数据分为10行6列矩阵作为网络训练输入向量p，将后60期数据分为10行6列矩阵作为网络训练输出向量t，取68期中的第59到68期这10期分为10行1列矩阵作为网络检验用的输入向量h，取第67到76期数据检验训练好的网络的仿真效果，将这10期分为10行1列矩阵作为网络仿真用的输入向量k。 b) 对样本数据进行预处理

数据的预处理和后处理是有效训练神经网络的关键步骤，直接影响到训练后神经网络的性能。常见的方法是将原始数据进行归一化处理，即通过一定的线性变换将输入和输出数据统一限制在（0，1）或（-1，1）区间内。本文选取premnmx函数将数据归一化到-1和1之间。

（2）隐层： a) 确定网络层数

多层的神经网络是非常强的。比方说，一个2层的神经网络，第一层是s型函数，而第2层是线形函数，就可以用来模拟任何函数(必须是连续有界的)[2]。因此，我们采用2层的神经网络。 b) 确定隐层神经元数

神经元数目的确定至关重要，数量太少，网络从样本中获取信息的能力就差，不足以概括和体现样本规律；数量过多，又可能“过度拟合”。本文采取试凑法，确定网络误差最小时对应的神经元，如表8：

表8.隐层神经元试凑过程

注：最终确定30为隐层神经元数目。

（3）训练过程：

a) 确定传递函数和训练函数

一个2层的神经网络，第一层是s型函数，而第2层是线形函数，就可以用来模拟任何函数(必须是连续有界的)。在此我们确定隐层传递函数为“logsin”，对于输出层，如果是Sigmoid函数,那么整个网络的输出就限制在一个较小的范围内( 0～1 之间的连续量) ; 如果是Pureline函数,那么整个网络的输入可以取任意值。为此我们确定输出层传递函数为“purelin”。 b) 训练网络

经过不断训练试验，最终确定神经网络结构为：其中输入输出层各10个神经元，隐层30个神经元，第一层传递函数为logsin，第二层传递函数为pureline。

（4）输出层

a) 网络推广效果检验

经过网络拟合，可以得到经归一化GDP数据拟合结果均方误差（MSE）为0.000101，说明模型预测还是比较精确的，可认为模型具有一定的参考价值，其

实际结果与拟合结果如图3所示，基本拟合了数据的变化趋势：

图3.测试样本拟合情况

b) 预测

将训练样本拟合出来的模型加入测试样本进行预测分析，以得出网络的仿真值。如图所示为仿真结果与目标输出之间的拟合情况，可以看出模型基本拟合了数据的变化趋势：

图4.仿真结果对比图

将预测结果进行反归一化，可以得到模型预测的实际结果值，预测2009到2010年各季度全国GDP数据,将其与真实值进行对比：

表9.单一BP神经网络预测结果对比

可以看出所构建的BP神经网络基本上拟合了原序列的发展趋势，训练样本学习效果比较好，然而利用构建好的网络进行推广测试，却发现效果并不稳定，拟合之后部分数据与原序列的误差相当高，说明单一的BP神经网络模型对于二重时间序列的预测精度有限，有必要进行改进。

（三） ARIMA模型

1. 模型介绍

ARIMA模型又称自回归滑动平均模型(Autoregressive Integrated Moving Average Model) ，ARIMA（p,d,q）的一般形式为：

Bdyt(B)t

其中，yt表示原序列，B表示后移算子。Bnytytn，n=1,2,„，d(1B)d

为d阶差分，B11B2B2…pBp为平稳可逆ARMA(p,q)模型的自回归系数多项式，(B)11B2B2…pBpI为平稳可逆ARMA(p,q)模型的移动平滑系数多项式，t~WN(0,2)为零均值的白噪声序列。

然而，当时间序列中具有明显的时间趋势和季节性变化时，单纯使用ARIMA模型进行分析往往会使预测效果不理想，这时只有将ARIMA模型和随机季节模型

（stochastic seasonal model）组合成季节时间序列模型（又称乘积时间序列模型），即SARIMA(seasonal ARIMA model)模型，才能较好地描述该时间序列。

若原时间序列进行d阶逐期差分和D阶长度为s的差分后才平稳，则建立季节时间序列模型ARIMA(p,d,q)(P,D,Q)S，表示为：

p(B)PBSdSytqQ(BS)t

其中，P是季节自回归阶数，Q是季节移动平均阶数，d(1B)d 表示d 阶

SSD

逐期差分，S(1Bs)D表示D阶季节差分，PB和Q(B)为季节P 阶自回

归算子和Q 阶移动平均算子。当P=D=Q=0 时，SARIMA 模型变为一般的ARIMA 模型。当P=D=Q=p=d=q时，SARIMA模型变为白噪声模型。

SARIMA模型的建模过程一般包括：原始数据预处理、模型的识别定阶、参数估计、诊断与检验以及模型的预测几个方面。

2. 建模过程（1）数据预处理

将1992年~2008年的68个GDP季度数据建立时间序列，并命名为GDP序列，运用Eviews6.0 软件对GDP序列进行单位根检验，结果（表10）显示，统计量P值等于0.9998，不能拒绝原假设，该序列未能通过ADF检验，是一个非平稳序列。

表10. GDP序列的单位根检验

绘制GDP序列的时间序列图（图5），观察到1992年~2009年的GDP季度数据时间序列具有很强的非平稳性和明显的指数趋势，因此，我们对这批数据取对数先消除增幅越来越大的现象，得到新的对数序列LOGGDP。取对数之后的时间

序列图（图6）仍存在明显的增长趋势，需要对数据进行一阶逐期差分处理，得到一阶差分序列DLOGGDP。观察一阶差分后的时间序列图，仍存在规律的季节变动，对数据进行差分步长与季节周期一致的季节差分处理，即进行步长为4的季节差分，得到新序列D4LOGGDP。

图5 .1992-2010年我国季度GDP时序图图6.取对数之后的时间序列图

图7.进行一阶差分之后的时间序列图图8.处理后的时间序列图

绘制处理后的时间序列图（图8），同时做单位根检验（表11）。单位根ADF检验结果表明T统计量均小于置信水平为1%、5%、10%下的检验值，P=0.0027，拒绝原假设，所以我们可以认定经过1阶4步差分后的对数序列通过ADF检验，序列的趋势和季节性已经基本消除，成为平稳时间序列，即LOGGDP~I（1，4）。

表11.LOGGDP对数序列1阶4步差分ADF检验

（2）模型的识别定阶

由于对原对数时间序列经过一阶差分后，序列的趋势消除，所以取d=1。又经过一阶四步季节差分后，序列季节性消除，可得D=1，S=4。所以对得到的平稳时间序列初步选定ARIMA(p,d,q)(P,D,Q)4模型中的ARIMA(p,1,q)(P,1,Q)4模型。

SARIMA模型其他阶数的识别与确定可以通过样本的自相关与偏自相关函数的观察获得。LOGGDP一阶四步差分后的自相关图和偏自相关图如图9。

图9.D4LOGGDP的的样本自相关图和偏自相关图

通过观察平稳序列D4LOGGDP的自相关和偏自相关图，差分处理后序列的自相关系数在滞后一期时接近于两倍标准差界限，后呈衰减趋于零，表现为拖尾性，因此阶数q可由显著不为零的偏自相关系数的数目决定，可取q=1；在偏自相关分析图中，滞后一期的偏自相关系数也显著不为零，但之后逐渐衰减趋于零，具

有拖尾性，因此阶数p可由显著不为零的偏自相关系数的数目决定，观察图可以取p=1，季节自回归阶数P可取4，也可取0。

但考虑到对自相关图和偏自相关图的直观判断带有很强的主观色彩，有时不太准确，因此本文通过对非平稳时间序列两种不同模型阶数的各种可能排列进行重复拟合，对明显不适合的排列形式予以删除，对可能性较大的ARIMA1,1,14,1,44，（如表12）。 ARIMA(1,1,4)(0,1,4)4模型分别进行参数比较ARIMA1,1,44,1,44 ARIMA1,1,10,1,44，

模型选定的主要依据是AIC、SC最小准则和修正R²最大原则。AIC（Akaike info criterion）准则兼顾了模型的简洁性和精确性，并排除了建模者的主观因素。BC（Schwarz criterion）准则对AIC 准则做了进一步改善。修正R²除具备R²的特征外，还避免了自变量个数的影响，综合了精度和变量数两个因素，能更有效地判断拟合优度。AIC、SC和修正R²的表达式分别为：

AIC2(lk)/n SC(2lklogn)/n

21

n1

(1R2)

np1

其中，l为对数似然函数，n是观测值数目，k是被估计的参数个数，p为自变量个数。

综合比较表12的各种参数之后，认为模型ARIMA1,1,44,1,44为相对最优模型，ARIMA1之。故选择模型ARIMA1,1,44,1,44和模型,1,44次0,1,4 ARIMA1,1,44进行参数估计。0,1,4

表12.各种拟合模型的参数比较

（3）模型的参数估计

分别对模型ARIMA1,1,44,1,44和模型ARIMA1,1,40,1,44进行参数估计：

表13.模型ARIMA1,1,44,1,44参数估计的结果

表14.模型ARIMA1,1,40,1,44参数估计结果

表13及表14、图10参数估计结果显示，模型ARIMA1,1,4和模型4,1,44

的滞后多项式倒数根均落在单位圆内，满足过程的平稳要求。ARIMA1,1,40,1,44但是模型ARIMA的,1,4参数无法通过显著性检验，故选择模型1,1,444

ARIMA1,1,4

0,更为合适。1,44

图10.ARIMA1,1,44,1,44和ARIMA1,1,40,1,44滞后多项式倒数根的分布图

（4）模型的诊断和体验 a）参数的显著性检验

参数的显著性检验就是要检验每个未知参数是否显著非零,目的在于使模型最精简。如果某个参数不显著, 即表示该参数所对应的那个自变量对因变量的影响不明显, 该自变量就可以从拟合预测模型中删除。最终模型将由一系列参数显著非零的自变量表示。

从表14中ARIMA1,1,40,1,44模型的参数检验结果可以看出，各参数t统计量的P值均小于0.10，故该模型的三个参数均显著。

b）模型的显著性检验

模型的显著性检验即为残差序列的白噪声检验, 若残差序列为白噪声, 则意味着所建立的模型已包含了原始序列的所有趋势, 从而模型适合应用于预测当前时间序列；若残差序列不是白噪声序列, 则认为残差序列还存在有用的信息未被提取，说明模型有改进的必要。

利用Ljung Box Q统计量对拟合预测模型进行残差序列白噪声检验，结果如图3-8所示。

从图11检验结果可以清楚看出，各阶延迟下，检验统计量Q值均大于对应

自由度卡方分布的检验值，且统计量的P值均大于0.05，即拒绝原假设的概率较大，可以认为这个拟合预测模型的残差序列属于白噪声序列，即模型显著有效，可以用于进一步的预测。

图11.对拟合预测模型的残差序列检验结果

模型的拟合效果也可以通过拟合值和真实值的对比看出。运用Eviews软件对ARIMA1,1,40,1,44模型进行拟合，得到拟合值和真实值的时序图（如图12）。该图显示出预测模型的拟合值和真实值之间的拟合结果良好。

图12.真实值、拟合值和残差序列图

因此，1992-2010年我国季度GDP时间序列的预测模型为：

(10.223936B)(1B)(1B4)log(GDP)=(1+0.388649B4)(10.679101B4)t

(5) 模型的预测

利用模型ARIMA1,1,40,1,44对2009-2010年8个季度的GDP数据进行预测。为提高预测精度，本文采用静态滚动预测法，即先利用1992~2008年数据所建立的ARIMA模型对2009年第一季度GDP进行预测，然后将预测值添加到原始序列中，用新的序列建立模型，再用所建立的模型对2009年第二季度GDP进行预测，以此类推，直至预测到2010年第四季度为止。预测结果（表15）显示，2009-2010年各季度预测值与真实值较为接近，相对误差均小于0.1，预测效果良好。预测值与真实值之间存在的偏差，可能是在拟合模型过程中不可避免的误差所致，例如2009年第一、二季度的偏差相对较大，也可能是因为经济危机的残留影响，而在建模过程中未将我国的宏观政策等考虑在内导致的。

表15.2009-2010年我国GDP真实值与预测值的比较

三、组合预测模型的引入及建模准备

（一）若干假设

1.假设一：时间序列是二重趋势时间序列，具备长期趋势性和周期波动性，在分

析时可以对二重趋势序列进行分解，分解为趋势波动变量和周期波动变量。

2.假设二：趋势波动变量和周期波动变量相互独立，即长期趋势和周期变动之间

互不影响。

3.假设三：时间序列是等时距时间序列，即时间序列相邻两个数据之间的时间间

隔相等。

（二）数据的分解

分析应用二重趋势时间序列时，通常可分解为趋势变动项和周期波动项，典型的分解方法是乘积模型：

X(t)=P(t)*T(t)

其中：X(t)为季度观测值，P(t)为观测值的趋势波动项，T(t)为观测值的周期波动项，t与年份i（i=1,2，···，n；n为观测数据年份个数），季度j满足关系式：

t=(i-1)*4+j

设序列X(1),X(2),···，X(h)，其中h为序列的长度，则观测数据序列季度个数h及季度周期数4的关系为h=4i。用中位移动平均法，可以提取不含周期波动的趋势项：

111

P(t)x(ti)[x(t2)x(t2)]，t=3,4，···，h-2

412

式中是以4为周期，以t为中心，2阶对称滑动平均数字滤波，P(t)经过滤波后则不再含有周期波动项，而周期波动项为：

T(t)

X(t)

P(t)

经过以上分解，就把观测数据的时间序列分解为趋势项和周期项，可分别对它们进行分析、建模，进行预测[6]。

四、组合预测模型构建过程

（一）长期趋势项模型构建过程

1. 选用构建长期趋势项模型的单一模型

由于灰色G（1，1）模型仅能较好地拟合时间序列的趋势性部分，而对于周期波动与随机波动，其预测精度则明显降低[5]。所以运用灰色G（1，1）模型对长期趋势项进行拟合预测。

2. 长期趋势项的灰色G（1，1）模型构建

运用1992年第三季度至2008年第四季度的GDP数据进行模型拟合，拟合得出的模型对2009-2010年数据进行预测，把预测结果与真实值进行对比。（1）构造累加生成列（SAS中记为yt）

X(1)(1)=X(0)(1)= 15611.12 X(1)(2)=X(1)(1)+ X(0)(2)=31814.44 X(1)(3)=X(1)(2)+ X(0)(3)=49061.69 ······

X(1)(66)=X(1)(66)+ X(0)(65)=4779756.25

（2）构造数据矩阵B（SAS中记为zt）和数据向量Yn

1(1)1(1)

[X(1)X(2)] 1(15611.1231814.44]) 1221(1)1(1)[X(2)X(3)] 1(31814.4449061.69） 12B2  11

(1)(1)

[X(65)X(66)] 1(4593448.954779756.25） 1 2223712.78 140438.07 1 

 4686602.60 116203.32

17247.25Yn 

186307.30

（3）计算参数矩阵

ˆ为待估参数，利用最小二乘法求解，可得 -0.033804T1T

ˆ(BB)BYn18804.21

即α=-0.033804，μ=18804.21

（4）得出预测模型（预测结果记为xp）

dX(1)

0.033804X(1)18804.21dt

X(0)(1)15611.12 X(0)(1)



-556272



556271.7

X(1)(k1)556271.7e0.033804k-556272

（5）残差检验

表16.GM（1,1）模型拟合趋势项

其中xp为预测值，error为绝对误差序列，rerror为相对误差序列。由于

篇幅有限，仅截取后10个数据（所有数据见附录4）。由拟合结果表16可以看出，相对误差序列中的相对误差较大，所以要对原模型进行残差修正，以提高精度。

（7）关联度检验

记关联系数为n，分析如下表17

表17.关联度检验结果

由表可得，关联系数的均值，即关联度为0.72。根据经验，当ρ=0.5时关联度大于0.6便满意。所以，关联度检验通过。

（8）后验差检验

表18.后验差检验结果

由表18可得，原始序列标准差S1=45064.35，残差标准差S2=4174.53， C= S2/ S1= 0.09 S0=0.6745S1= 30395.9

由于残差的取值范围为[0,19505.51]，则

ˆ00,14722.79 0

ˆ0S1 Pp00

所以，P=1，C

（8）预测



表19.GM（1,1）模型趋势项预测结果

运用GM（1,1）模型对2009年第一季度至2010年第二季度的GDP进行预测，然后与实际值进行对比，所得结果的相对误差在0.07左右。

（二）季节波动项模型构建过程

1. 选用构建季节波动模型的单一模型

周期波动项具有非常复杂的非线性结构，对其准确预测是又一难点。BP神经网络模型是目前比较成熟的算法，在函数逼近和数据拟合方面具有很大优越性，但它常会忽略一些大的噪音数据[4-5]。而自回归滑动平均（ARIMA）模型对噪音数据具有很强的预测能力。由此，结合这两个模型的优势，分别用BP 神经网络模型和ARIMA模型对周期项波动预测，根据权重优化模型，建立周期波动项的组合预测模型。

2. 波动项的BP神经网络构建

（1）采用1992年第三季度至2010年第二季度各季度全国GDP周期项数据共72

期。取其中前68期作为历史数据样本，将前60期数据分为10行6列矩阵作为网络训练输入向量p，将后60期数据分为10行6列矩阵作为网络训练输出向量t，取68期中的第59到68期这10期分为10行1列矩阵作为网络检验用的输入向量h，取第63到72期数据检验训练好的网络的仿真效果，将这10期分为10行1列矩阵作为网络仿真用的输入向量k。

（2）经过不断训练试验，最终确定神经网络结构为：其中输入输出层各10个神

经元，隐层25个神经元，第一层传递函数为logsin，第二层传递函数为

pureline。

（3）拟合情况：均方误差（MSE）为：

9.24753e-028

图13.测试样本拟合情况

（4）仿真结果：

图14.仿真结果对比图

将预测结果进行反归一化，可以得到模型预测的实际结果值，预测2009到2010年第二季度各季度全国GDP数据波动项,将其与真实值进行对比：

表20.BP神经网络对波动项的预测值

1. 波动项的ARIMA模型构建

用ARIMA模型拟合季节波动项，通过绘制时序图和观察自相关图和偏自相关图，发现模型具有步长为4的周期波动，对模型进行4步差分。再对差分后的数据进行识别，经过一系列的检验认为选择模型ARIMA(20,0,20)4 较为适合，模型的参数估计结果如下：

表21.拟合季节波动的ARIMA(20,0,20)4模型参数估计结果

由表21可得，拟合季节项的模型表达式为：

(10.530140B20)(1B4)yARIMA=(1+0.853537B20)t

根据所拟合的ARIMA(20,0,20)4模型对2009-2010年8个季度的季节波动项进行预测，结果如表22所示：

表22.2009-2010年的季节波动项预测值

4.季节波动模型的最终构建过程（1）确定权重

采用最优线性组合模型原理利用样本期实际值和各单项预测模型的拟合值进行线性回归，可以确定各单项预测模型的权重。

由于不同模型建模时造成部分时期数据缺失，如ARIMA模型差分导致信息丢失等。为了建立更加接近2009年和2010年真实情况的模型，分别选取BP神经网络模型和ARIMA模型2006年第二季度到2008年第四季度共10期波动项的

ˆBP(k)、yˆARIMA(k)作为自变量，选取对应真实值y(k)作为因变量，最终建拟合值y

立线性回归模型：

ˆy(k)=1.038y

ARIMA

ˆBP(k) (k)-0.035y

表23.波动项的拟合值与真实值

（2）周期波动项组合预测模型

最终建立周期波动项组合预测模型为：

ˆy(k)=1.038y



ARIMA

ˆBP(k) (k)-0.035y

（三）随机干扰项模型构建过程

由前面数据的分解过程可知，X(t)=P(t)*T(t)（X(t)为季度观测值，P(t)为观测值的趋势波动项，T(t)为观测值的周期波动项，t为时期数），因此有：

y(k)

y(k)y

(k)

为了改进组合模型，提高模型预测的精度，解决传统确定性因素分解对残差信息完全浪费的缺点，本文尝试对模型随即干扰项进一步进行拟合，以便验证是否存在有效信息没有提取，避免信息浪费，由组合模型：y(k)式：

y(k)y

GM

(k)(k)可得随机干扰项模型表达

(k)y(k)y(k)y

GM

(k)

注：其中，

ˆy(k)=1.038y



ARIMA

ˆBP(k) (k)-0.035y

因此，根据数据可得随机干扰项数值：

表24.随即干扰项的计算

由残差项(k)的散点图可知其具有一定的线性趋势，为此对(k)进行残差自回归拟

合检验，在SAS8.0拟合结果如下：

表25.残差自回归检验结果

由于DW值为1.9412接近2，且p值为0.3013大于0.05，所有有95%的把握说明以上所得组合预测模型对时间序列的信息提取比较充分，无需对残差项进行拟合，分析结束。

（一）组合预测模型构建最终过结果

经过不断的检验与尝试，最终二重时间序列组合预测模型拟合口径如下：



GM

y(k)y(k)y

(k)

ˆy(k)=1.038y



ARIMA

ˆBP(k)y(k) (k)-0.035y



(10.530140B)(1B)yARIMA=(1+0.853537B)tdX(1)

0.033804X(1)18804.21dt

X(0)(1)15611.12 X(0)(1)

y(k)





-556272

GM

(k)



556271.7

X(1)(k1)556271.7e0.033804k-556272

其中，y

GM

ˆBP(k)为周期波动项单一拟合模，y(k)为趋势项拟合模型GM（1，1）

ˆARIMA(k)为周期波动项单一拟合模型型BP神经网络（网络结构图如图15），yARIMA，y(k)为周期波动项组合拟合模型，y(k)为二重时间序列最终组合预测模型。

s



注：第一层传递函数为logsin，第二层传递函数为pureline。

五、组合模型与单一模型预测精度对比分析

（一）直接结果对比分析

将组合模型的预测值和三个单一模型的预测值分别与真实值对比发现，组合模型平均绝对误差百分比最低，仅为0.060618，GM（1，1）模型次之，神经网络模型预测精度最低。为了进一步验证组合模型相对最优，笔者引进基于非参数检验的Friedman检验。

表26.组合模型与单一模型拟合结果对比

（二）非参数Friedman检验

Friedman检验是利用秩实现对多个总体分布是否存在显著差异的非参数检验方法，其原假设是：多个配对样本来自的多个总体分布无显著差异。多配对样本的Friedman检验时，首先以行为单位将数据按升序排序，并求得各变量值在各自行中的秩；然后，分别计算各组样本下的秩总和与平均秩。[7]

经过Friedman检验，卡方值为14.000，Sig.值为0.003，结合各种模型的平均绝对误差百分比，有99%的把握确定所建立的组合预测模型相对最优，采取组合预测法建立的模型预测精度更佳。

至此，模型建立结束。

六、总结与展望

（一）全文总结

现实中许多时间序列都具备趋势波动和周期波动的二重性质，运用单一模型对二重趋势时间序列的拟合预测分析具有一定的不足之处，于是提出了组合模型。本文针对1992-2010年国内生产总值（GDP）季度累加数据，分别运用单一模型和组合模型进行拟合预测进行分析。

运用单一GM （1,1）灰色预测模型对消除季节影响后的GDP时间序列进行拟合，并用GM （1,1）残差模型进行修正，修正后相对误差较大。单一BP神经网络对GDP的周期波动项进行拟合，仅对训练样本拟合效果好，对于测试样本的拟合结果误差相当高。而单一ARIMA模型和随机季节模型结合为SARIMA对季度GDP1992年-2008年数据进行建模，对2009-2010年进行预测并与真实值进行对比，相对误差也不容乐观。

由于单一模型对二重趋势时间序列拟合预测的局限性，引入SARIMA模型、GM（1,1）模型和BP神经网络的组合模型。首先对原始序列进行分解得到趋势波动项和周期波动项，运用BP神经网络和SARIMA模型对周期项进行预测，预测结果通过最小二乘法得到周期项的组合预测值；而对长期趋势项则运用GM（1,1）模型进行预测，结果与周期项组合预测值相乘对比组合模型的预测值得到残差项。对残差项进行自回归拟合发现信息提取充分，没必要拟合残差，由前面相乘值即为组合模型预测值。将组合模型预测结果与单一模型对比，相对误差明显降低，经过Friedman检验确定组合模型优于单一模型。

对于二重趋势时间序列，组合模型消除了单一模型对二重趋势片面拟合的不足之处，把GM（1,1）模型、SARIMA模型和BP神经网络各自的优势予以结合，取长补短，得到较为精确的预测结果。

（二）研究展望

非线性时间序列预测问题是当今统计分析领域研究的难点。在宏观经济领域或者实际商业运作中，许多数据都具有二重趋势的时间序列，存在着大量的不确定性，因此本文虽然只研究分析国内生产总值GDP的预测问题，但是该组合模

型实际上具有普遍的意义，可以应用到更广泛地应用到其他宏观经济或者企业生产情况或工业设备管理的预测问题。

另外，由于技术的原因和笔者自身学识有限，本文的研究尚有许多不足之处。二重趋势时间序列组合模型虽然在预测精度上优于单一模型，但在一定程度上仍受到组合模型各个部分的预测准确性影响。而在周期波动项呈现出太多的干扰数据时，BP神经网络模型和ARIMA模型的预测准确率都有待改善。因此，未来的研究方向可以建立在对周期波动项噪音数据的消除上，如引入小波变换或小波包分解的方法对这类时间序列的组合预测模型进行完善。

[1]张迎春，肖冬荣，赵远东.基于时间序列神经网络的气象预测研究. 武汉理工大学学报.2003，4（29）

[2] 宋仙磊，刘业政，陈思凤，许波.二重趋势时间序列的灰色组合预测模型.计算机工程与应用. 2011，47（8）

[3]徐国祥主编：《统计预测与决策》，上海财经大学出版社2010年版

[4] 宋玉强.人工神经网络在时间序列预测中的应用研究. 西安建筑科技大学硕士论文.2005

[5]戴丹. BP神经网络用于市场预测的研究. 武汉理工大学硕士论文.2006

[6] 汪世汉, 张国忠, 黄克明. 山东省发电量预测模型的研究.理论分析与研究.2002（01） [7]王静龙，梁小筠编著：《非参数统计分析》，高等教育出版社2010年版

[8] 颜杰，相丽驰，方积乾.灰色预测模型及SAS实现[J].中国卫生统计，2006年，第23卷第1期：75-85页

[9]蔡建平，朱秀萍等，SAS社会统计实用教程[M].清华大学出版社，2006年，105-107页 [10] 张晓峒.EViews使用指南与案例 [M].北京：机械工业出版社，2007.2.

[11] 牛东晓，乞建勋，刑棉.二重趋势性季节型电力负荷预测组合灰色神经网络模型[J].中国管理科学，2001.12（第9卷第6期）.

[12] 胡江波，刘金兰.回归_时间序列组合模型在石油企业费用预算编制中的应用[N].沈阳理工大学学报，2006.2(第25卷第1期).

[13] 王洪涛，柯跃.基于GMDH和BP的组合预测模型在工业产品销售收入中的应用[J].技术与市场，2010（第17卷第2期）.

[14]赵喜仓，周作杰. 基于SARIMA 模型的我国季度GDP 时间序列分析与预测[J]. 统计与决策，2010,(22).

附录1原始数据

附录3.GM（1,1）拟合结果残差检验（消除季节影响后）

附录4 GM（1,1）模型拟合趋势项

附录5.单一GM（1，1）模型SAS代码

data a1; input xt@@; t=_n_; yt+xt; index=1;

zt=-(yt+LAG(yt))/2; datalines;

12672.85081 15056.19753 15746.36017 15698.7564 16560.99578 19322.99267 20132.84194 20602.78793 23093.76977 26771.54194 27817.76417 28103.59229 30211.23111 34501.29932 35418.00683 35448.09473 36332.60614 41003.89927 41417.79147 41502.21739 41416.3521

46441.36032 46170.59515 46048.23002

44587.22996 49455.23588 49151.39902 49213.95758 47869.54894 52553.82147 52410.49038 52289.61565 52601.24142 58125.21977 58105.42191 57850.81728 59359.0815

65037.60352 64702.74268 63938.61548

64648.40121 70870.4977

70880.04594 70164.54864

73529.78403 79543.64758 79646.60569 79196.6465 85144.01389 93543.42734 93845.62058 93223.12019 99657.47022 108831.5801 108088.2439 107834.781 115448.8281 126789.2173 125740.0883 126130.3336 139498.9052 154119.6785 153697.0844 154990.7995 168867.9673 186643.0921 185208.6345 183116.1081 ;

proc gplot data=a1; plot xt*t;

symbol c=black v=star i=join; run; proc iml; use a1;

read all var{zt index} into B where (zt^=.); read all var{xt} into Yn where (zt^=.); ahat=inv(B`*B)*B`*Yn; ahatt=ahat`;na={a u};

create a2 from ahatt[colname=na]; append from ahatt; quit; data a3; set a2;index=1; data a4;

set a1;if _n_=1;xt0=xt; keep xt0 index; data a5;

merge a1 a3 a4;by index; if _n_=1 then xp=xt; else do

yt1=(xt0-u/a)*exp(-a*(t-1))+u/a; yt0=(xt0-u/a)*exp(-a*(t-2))+u/a; xp=yt1-yt0; end;

error=abs(xp-xt); rerror=error/xt*100;

drop yt index zt a u xt0 yt1 yt0; proc means data=a5; var xt; var error; output var=val; run;

ODS rtf file='f:/aa.rtf'; proc print data=a5; run;

ods rtf close;

data a6; /*建立GM（1,1）残差模型*/ set a5; et+error;

kt=-(et+lag(et))/2; proc iml;

use a6;

read all var{kt index} into C where(kt^=.); read all var{error} into Zn where(kt^=.); bhat=inv(C`*C)*

C`*Zn;

bhatt=bhat`;na2={ae ue};

create a7 from bhatt[colname=na2]; append from bhatt; quit; data a8; set a7;index=1; data a9;

set a6;if _n_=1;error0=error; keep error0 index; data a10;

merge a1 a3 a6 a8 a9;by index; if _n_=1 then ep=error; else do;

if t>=3 then d=1; else d=0;

yte1=(xt0-u/a)*exp(-a*(t-1))+u/a+d*(-ae)*(2936.84-ue/ae)*exp(-a*(t-2)); yte0=(xt0-u/a)*exp(-a*(t-2))+u/a+d*(-ae)*(2936.84-ue/ae)*exp(-a*(t-3)); xpe=yte1-yte0; end;

error2=abs(xpe-xt); rerror2=error2/xt*100;

drop yt index zt a u xt0 yt1 yt0 et kt ae ue error0 ep d yte1 yte0; ODS rtf file='f:/aa.rtf'; proc print data=a10; run;

ods rtf close;

proc means data=a10; /*残差检验*/ var xt; var error2; output var=val2; run;

data a11; /*关联度检验和后验差检验*/ set a10;

n=(14.5879387+0.5*25195.58)/(error2+0.5*25195.58); proc means data=a11; var n;

output mean=avgn; run;

data p; /*预测*/ input k@@;

cards; 69 70 71 72 73 74 75 76 ;

data pre;

set p;

xps=517136.3317*exp(0.034202*(k-1))-504463.4817+2825.585429*exp(0.016111*(k-2)); abc=lag1(xps);

do;

if k=69 then

xpre=xps-(517136.3317*exp(0.034202*67)-504463.4817+2825.585429*exp(0.016111*66)); else xpre=xps-abc;

end;

drop xps abc;

ODS rtf file='f:/aa.rtf';

proc print data=pre;

run;

ods rtf close;

附录6. GM（1，1）趋势项拟合SAS代码

data a1;

input xt@@; t=_n_; yt+xt;

index=1;

zt=-(yt+LAG(yt))/2;

datalines; 15611.12 16203.32 17247.25 18941.06 20312.9 21334.2 23160.62 25894.25 27851.46 28927.9 30768.37 33456.1 35330.92 36243.04 37733.63 39910.3 41457.6 42218.59 43426.32 45097.04 46227.18 46666.3 47386.47 48501.73 49341.43 49794 50562.89 51699.61 52591.12 53347.45 54705.77 56732.12 58255.88 59237.78 60854.45 63125.86 64690.46 65498.75 66952.33 69191.84 70962.29 72214.04 74314.1 77534.44 80040.55 81927.53 85324.44 90411.16 94130.2 96280.64 99805.15 105023.7 108930.9 111395.2 115670.2 122177.9 127280 131031.3 137697.6 147979.6 155607.5 160108.4 167783.9 178428.9 184899.9 186307.3

;

二重趋势时间序列的组合预测模型

广东外语外贸大学陈文聪、庄文钿、苏玲

摘要

关键词：二重趋势时间序列；组合模型；国内生产总值（GDP）；预测

一、概述......................................................... 3

1. 研究思路 ................................................. 4 2. 数据来源 ................................................. 5

二、单一预测模型的介绍及建模过程 ................................ 6

（一）灰色预测模型 ........................................... 6

1. 模型介绍 ................................................ 16 2. 建模过程 ............................................... 17

三、组合预测模型的引入及建模准备 ............................... 25

（一）长期趋势项模型构建过程 ................................... 26

一、概述

（一）

时间序列研究概况

（二）本文研究的意义、亮点

1. 数据：不同于以往研究只选择年度数据或者月度数据，本文选择季度GDP

数据进行研究，尝试拟补季度GDP预测研究的不足，并且本文研究期数多达76期，使研究结果更有精确；

2. 单一模型：本文对单一模型的选择精雕细琢，结合不同单一模型的预测优

点最终确定以GM（1，1）预测趋势项，以BP神经网络和ARIMA通过最小二乘法加权预测周期项，以自回归模型拟合残差项；

3. 组合模型：为了最大限度减少信息丢失，本文在建立周期项组合模型时采

取最小二乘法加权使信息量提取充分，并且在总得组合模型建立时尝试拟合残差自回归模型，确保较少信息遗漏；

4. 模型检验：对于每个单一模型的建模过程，笔者都科学严谨检验每个过程，

（三）本文的研究思路及建模准备

1. 研究思路（1）研究大体框架

（2）二重趋势时间序列的组合预测模型构造

ˆBP(k)是BP神经网络对周期波动项的预测值，yˆARIMA(k)是ARIMA模型对设y

ˆS(k)为：周期波动项的预测值，则周期项的组合预测值y

y(k)aby

BP

(k)bARIMAy

ARIMA

(k)

注：其中

、

ARIMA

分别为BP神经网络模型和ARIMA模型在周期波动项预测模型中的系数的系数。

ˆGM(k)是GM灰色预测模型对趋势项的预测值，把趋势项和组合后的周期波y

ˆ(k)为：动项进行再次组合，则二重趋势时间序列组合预测值y

y(k)y(k)y

GM

(k)(k)

2. 数据来源

二、单一预测模型的介绍及建模过程

（一）

灰色预测模型

1. 模型介绍

（1）给定序列在t时刻的观测值： X={X(1),X(2),···,X(n)}

（2）作一阶累加生成1-AGO，得到新的数列： X={X(1),X(2),···,X(n)}，其中：X (3)构造一阶微分方程：

(1)

(0)

(1)

(k)X(0)(i),k1,2,...n。

i1

dX(1)

aX(1)u dt

利用最小二乘法求得参数a，u：

a

ˆ(BTB)1BTYN， a

u

1(1)(1)1[X(1)X(1)]2X(0)(2)1(1)(1)(0)X(3)1[X(2)X(3)]

,YN其中：B2。   X(0)(n)1

[X(1)(n1)X(1)(n)]1

2求得微分方程，即可得预测模型：

ˆ(1)(k1)[X(0)(1)u]eaku X

2. 建模过程（1）数据预处理

对原序列做时序图可以看出该序列明显呈现固定周期性变化，即具有“季节”

表1.运用按季平均求季节指数

如表1求得各季度季节指数，用原序列除以相应的季节指数即可得到消除季

节影响后的时间序列（详见附录2）（2）GM（1,1）模型建立及修正

该模型运用SAS8.0软件进行拟合预测操作，代码详见附录。首先，记原始时间序列为：

X(0)={X(0)(1),X(0)(2),X(0)(3),X(0)(4),······,X(0)(68)}，上标0表示原始时间序列。（SAS中记为xt）

a）构造累加生成列（SAS中记为yt）

X(1)(1)=X(0)(1)= 12672.85 X(1)(2)=X(1)(1)+ X(0)(2)=27729.05 X(1)(3)=X(1)(2)+ X(0)(3)=43475.41 ······

X(1)(68)=X(1)(67)+ X(0)(68)=4793643.35

b）构造数据矩阵B（SAS中记为zt）和数据向量Yn

12[X(1)(1)X(1)

(2)] 11(12672.852772912.05]) (1)(1)1B[X(2)X(3)] 1(27729.0543475.41） 22  12[X(1)(67)X(1)

(68)] 112(4610527.244793643.35）20200.95 135602.23 1 

 4702085.30 115056.20

Yn15746.36

 183116.11

c）计算参数矩阵

ˆ

为待估参数，利用最小二乘法求解，可得 111

  

ˆ(BTB)1BTYn

-0.034202

 

 17253.66

即α=-0.034202，μ=17253.66

d）得出预测模型（预测结果记为xp）

dX(1)

0.034202X(1)17253.66dt

X(0)(1)12672.85 X(0)(1)



-504463.4817



517136.3317

X(1)(k1)517136.3317e0.034202k504463.4817

e）残差检验

表2.GM（1,1）拟合结果残差检验

表2仅截取前10个数据（所以数据详见附录3）。由表2可以看出，相对误差序列中的相对误差较大，所以要对原模型进行残差修正，以提高精度。

f）利用残差对原模型进行修正

将残差记为：

e(0)={e(0)(1’), e(0)(2’),„, e(0)(67’)}={2936.84,2872.72,„,11146.49}

e(0)的累加生成序列为：

e(1)={e(1)(1’),e(1)(2’),„,e(1)(67’)}={29336.84,5809.46,„,346283.25} e(1)可建立相应的GM（1,1）模型，求得参数

e0.016111ˆe2778.27 

eˆ(1)(k1)[e(0)(1)则e

ekek]e175382.3741e0.016111172445.5341 ee

(k1)ˆ(1)(k1)]'2825.585429e0.016111对上述求导，得：[e

这样得到经过残差修正后的模型：

(k1)ˆ(1)(k1)517136.3317e0.034202k504463.4817(k1)2825.585429e0.016111X

1,k2

其中：(k1)

0,k2

表3.GM（1,1）残差模型拟合结果

其中表3为GM（1,1）残差修正模型的预测结果的前10个数据（所有数据详见附录2）。由表3可见，与GM（1,1）相比，修正模型的精度相对提高。

g）关联度检验

记关联系数为n，分析如表4：

表4.关联系数检验

由表4可得，关联系数的均值，即关联度为0.75。根据经验，当ρ=0.5时关联度大于0.6便满意。所以，关联度检验通过。

h）后验差检验

表5.后验差检验

由表5可得，原始序列标准差S1=45773.05，残差标准差S2=4290.74， C= S2/ S1= 0.093739438 S0=0.6745S1= 30873.92223

由于残差的取值范围为[14.5879387,25195.58]，则

ˆ00,20101.43 0

ˆ0S1 Pp00

所以，P=1，C

i）用模型结果预测2009-2010年的数据，预测结果乘以各季度的季节指数得到GDP预测值。

表6.2009-2010年预测结果



（二）

BP神经网络

1. 模型介绍

BP神经网络模型对时间序列预测的求解步骤如下：

（1）网络结构的设计：BP神经网络结构直接影响到模型的精度，但至今尚无统

一的理论方法，本文主要通过经验和多次尝试建立网络结构。

（2）样本的预处理：对输入的训练数据进行归一化处理，使输入数据的取值在

（-1，1）或者（0，1）之间。

（3）BP网络的学习规则：BP网络的学习规则一般包括以下几个步骤：①初始化

（4）反归一化：学习结束，对输入和输出样本进行反归一化处理，即可得到预测值。

2. 建模过程

BP网络的设计主要包括输入层,隐层,输出层及各层之间的传输函数几个方面。

（1）输入层：

a) 将数据分为训练样本数据和测试数据

（2）隐层： a) 确定网络层数

表8.隐层神经元试凑过程

注：最终确定30为隐层神经元数目。

（3）训练过程：

a) 确定传递函数和训练函数

经过不断训练试验，最终确定神经网络结构为：其中输入输出层各10个神经元，隐层30个神经元，第一层传递函数为logsin，第二层传递函数为pureline。

（4）输出层

a) 网络推广效果检验

经过网络拟合，可以得到经归一化GDP数据拟合结果均方误差（MSE）为0.000101，说明模型预测还是比较精确的，可认为模型具有一定的参考价值，其

实际结果与拟合结果如图3所示，基本拟合了数据的变化趋势：

图3.测试样本拟合情况

b) 预测

图4.仿真结果对比图

将预测结果进行反归一化，可以得到模型预测的实际结果值，预测2009到2010年各季度全国GDP数据,将其与真实值进行对比：

表9.单一BP神经网络预测结果对比

（三） ARIMA模型

1. 模型介绍

ARIMA模型又称自回归滑动平均模型(Autoregressive Integrated Moving Average Model) ，ARIMA（p,d,q）的一般形式为：

Bdyt(B)t

其中，yt表示原序列，B表示后移算子。Bnytytn，n=1,2,„，d(1B)d

然而，当时间序列中具有明显的时间趋势和季节性变化时，单纯使用ARIMA模型进行分析往往会使预测效果不理想，这时只有将ARIMA模型和随机季节模型

（stochastic seasonal model）组合成季节时间序列模型（又称乘积时间序列模型），即SARIMA(seasonal ARIMA model)模型，才能较好地描述该时间序列。

若原时间序列进行d阶逐期差分和D阶长度为s的差分后才平稳，则建立季节时间序列模型ARIMA(p,d,q)(P,D,Q)S，表示为：

p(B)PBSdSytqQ(BS)t

其中，P是季节自回归阶数，Q是季节移动平均阶数，d(1B)d 表示d 阶

SSD

逐期差分，S(1Bs)D表示D阶季节差分，PB和Q(B)为季节P 阶自回

归算子和Q 阶移动平均算子。当P=D=Q=0 时，SARIMA 模型变为一般的ARIMA 模型。当P=D=Q=p=d=q时，SARIMA模型变为白噪声模型。

SARIMA模型的建模过程一般包括：原始数据预处理、模型的识别定阶、参数估计、诊断与检验以及模型的预测几个方面。

2. 建模过程（1）数据预处理

表10. GDP序列的单位根检验

图5 .1992-2010年我国季度GDP时序图图6.取对数之后的时间序列图

图7.进行一阶差分之后的时间序列图图8.处理后的时间序列图

表11.LOGGDP对数序列1阶4步差分ADF检验

（2）模型的识别定阶

SARIMA模型其他阶数的识别与确定可以通过样本的自相关与偏自相关函数的观察获得。LOGGDP一阶四步差分后的自相关图和偏自相关图如图9。

图9.D4LOGGDP的的样本自相关图和偏自相关图

有拖尾性，因此阶数p可由显著不为零的偏自相关系数的数目决定，观察图可以取p=1，季节自回归阶数P可取4，也可取0。

AIC2(lk)/n SC(2lklogn)/n

21

n1

(1R2)

np1

其中，l为对数似然函数，n是观测值数目，k是被估计的参数个数，p为自变量个数。

表12.各种拟合模型的参数比较

（3）模型的参数估计

分别对模型ARIMA1,1,44,1,44和模型ARIMA1,1,40,1,44进行参数估计：

表13.模型ARIMA1,1,44,1,44参数估计的结果

表14.模型ARIMA1,1,40,1,44参数估计结果

表13及表14、图10参数估计结果显示，模型ARIMA1,1,4和模型4,1,44

ARIMA1,1,4

0,更为合适。1,44

图10.ARIMA1,1,44,1,44和ARIMA1,1,40,1,44滞后多项式倒数根的分布图

（4）模型的诊断和体验 a）参数的显著性检验

从表14中ARIMA1,1,40,1,44模型的参数检验结果可以看出，各参数t统计量的P值均小于0.10，故该模型的三个参数均显著。

b）模型的显著性检验

利用Ljung Box Q统计量对拟合预测模型进行残差序列白噪声检验，结果如图3-8所示。

从图11检验结果可以清楚看出，各阶延迟下，检验统计量Q值均大于对应

图11.对拟合预测模型的残差序列检验结果

图12.真实值、拟合值和残差序列图

因此，1992-2010年我国季度GDP时间序列的预测模型为：

(10.223936B)(1B)(1B4)log(GDP)=(1+0.388649B4)(10.679101B4)t

(5) 模型的预测

表15.2009-2010年我国GDP真实值与预测值的比较

三、组合预测模型的引入及建模准备

（一）若干假设

1.假设一：时间序列是二重趋势时间序列，具备长期趋势性和周期波动性，在分

析时可以对二重趋势序列进行分解，分解为趋势波动变量和周期波动变量。

2.假设二：趋势波动变量和周期波动变量相互独立，即长期趋势和周期变动之间

互不影响。

3.假设三：时间序列是等时距时间序列，即时间序列相邻两个数据之间的时间间

隔相等。

（二）数据的分解

分析应用二重趋势时间序列时，通常可分解为趋势变动项和周期波动项，典型的分解方法是乘积模型：

X(t)=P(t)*T(t)

t=(i-1)*4+j

111

P(t)x(ti)[x(t2)x(t2)]，t=3,4，···，h-2

412

式中是以4为周期，以t为中心，2阶对称滑动平均数字滤波，P(t)经过滤波后则不再含有周期波动项，而周期波动项为：

T(t)

X(t)

P(t)

经过以上分解，就把观测数据的时间序列分解为趋势项和周期项，可分别对它们进行分析、建模，进行预测[6]。

四、组合预测模型构建过程

（一）长期趋势项模型构建过程

1. 选用构建长期趋势项模型的单一模型

2. 长期趋势项的灰色G（1，1）模型构建

X(1)(1)=X(0)(1)= 15611.12 X(1)(2)=X(1)(1)+ X(0)(2)=31814.44 X(1)(3)=X(1)(2)+ X(0)(3)=49061.69 ······

X(1)(66)=X(1)(66)+ X(0)(65)=4779756.25

（2）构造数据矩阵B（SAS中记为zt）和数据向量Yn

1(1)1(1)

(1)(1)

[X(65)X(66)] 1(4593448.954779756.25） 1 2223712.78 140438.07 1 

 4686602.60 116203.32

17247.25Yn 

186307.30

（3）计算参数矩阵

ˆ为待估参数，利用最小二乘法求解，可得 -0.033804T1T

ˆ(BB)BYn18804.21

即α=-0.033804，μ=18804.21

（4）得出预测模型（预测结果记为xp）

dX(1)

0.033804X(1)18804.21dt

X(0)(1)15611.12 X(0)(1)



-556272



556271.7

X(1)(k1)556271.7e0.033804k-556272

（5）残差检验

表16.GM（1,1）模型拟合趋势项

其中xp为预测值，error为绝对误差序列，rerror为相对误差序列。由于

（7）关联度检验

记关联系数为n，分析如下表17

表17.关联度检验结果

由表可得，关联系数的均值，即关联度为0.72。根据经验，当ρ=0.5时关联度大于0.6便满意。所以，关联度检验通过。

（8）后验差检验

表18.后验差检验结果

由表18可得，原始序列标准差S1=45064.35，残差标准差S2=4174.53， C= S2/ S1= 0.09 S0=0.6745S1= 30395.9

由于残差的取值范围为[0,19505.51]，则

ˆ00,14722.79 0

ˆ0S1 Pp00

所以，P=1，C

（8）预测



表19.GM（1,1）模型趋势项预测结果

运用GM（1,1）模型对2009年第一季度至2010年第二季度的GDP进行预测，然后与实际值进行对比，所得结果的相对误差在0.07左右。

（二）季节波动项模型构建过程

1. 选用构建季节波动模型的单一模型

2. 波动项的BP神经网络构建

（1）采用1992年第三季度至2010年第二季度各季度全国GDP周期项数据共72

（2）经过不断训练试验，最终确定神经网络结构为：其中输入输出层各10个神

经元，隐层25个神经元，第一层传递函数为logsin，第二层传递函数为

pureline。

（3）拟合情况：均方误差（MSE）为：

9.24753e-028

图13.测试样本拟合情况

（4）仿真结果：

图14.仿真结果对比图

将预测结果进行反归一化，可以得到模型预测的实际结果值，预测2009到2010年第二季度各季度全国GDP数据波动项,将其与真实值进行对比：

表20.BP神经网络对波动项的预测值

1. 波动项的ARIMA模型构建

表21.拟合季节波动的ARIMA(20,0,20)4模型参数估计结果

由表21可得，拟合季节项的模型表达式为：

(10.530140B20)(1B4)yARIMA=(1+0.853537B20)t

根据所拟合的ARIMA(20,0,20)4模型对2009-2010年8个季度的季节波动项进行预测，结果如表22所示：

表22.2009-2010年的季节波动项预测值

4.季节波动模型的最终构建过程（1）确定权重

采用最优线性组合模型原理利用样本期实际值和各单项预测模型的拟合值进行线性回归，可以确定各单项预测模型的权重。

ˆBP(k)、yˆARIMA(k)作为自变量，选取对应真实值y(k)作为因变量，最终建拟合值y

立线性回归模型：

ˆy(k)=1.038y

ARIMA

ˆBP(k) (k)-0.035y

表23.波动项的拟合值与真实值

（2）周期波动项组合预测模型

最终建立周期波动项组合预测模型为：

ˆy(k)=1.038y



ARIMA

ˆBP(k) (k)-0.035y

（三）随机干扰项模型构建过程

由前面数据的分解过程可知，X(t)=P(t)*T(t)（X(t)为季度观测值，P(t)为观测值的趋势波动项，T(t)为观测值的周期波动项，t为时期数），因此有：

y(k)

y(k)y

(k)

y(k)y

GM

(k)(k)可得随机干扰项模型表达

(k)y(k)y(k)y

GM

(k)

注：其中，

ˆy(k)=1.038y



ARIMA

ˆBP(k) (k)-0.035y

因此，根据数据可得随机干扰项数值：

表24.随即干扰项的计算

由残差项(k)的散点图可知其具有一定的线性趋势，为此对(k)进行残差自回归拟

合检验，在SAS8.0拟合结果如下：

表25.残差自回归检验结果

（一）组合预测模型构建最终过结果

经过不断的检验与尝试，最终二重时间序列组合预测模型拟合口径如下：



GM

y(k)y(k)y

(k)

ˆy(k)=1.038y



ARIMA

ˆBP(k)y(k) (k)-0.035y



(10.530140B)(1B)yARIMA=(1+0.853537B)tdX(1)

0.033804X(1)18804.21dt

X(0)(1)15611.12 X(0)(1)

y(k)





-556272

GM

(k)



556271.7

X(1)(k1)556271.7e0.033804k-556272

其中，y

GM

ˆBP(k)为周期波动项单一拟合模，y(k)为趋势项拟合模型GM（1，1）

ˆARIMA(k)为周期波动项单一拟合模型型BP神经网络（网络结构图如图15），yARIMA，y(k)为周期波动项组合拟合模型，y(k)为二重时间序列最终组合预测模型。

s



注：第一层传递函数为logsin，第二层传递函数为pureline。

五、组合模型与单一模型预测精度对比分析

（一）直接结果对比分析

表26.组合模型与单一模型拟合结果对比

（二）非参数Friedman检验

至此，模型建立结束。

六、总结与展望

（一）全文总结

（二）研究展望

型实际上具有普遍的意义，可以应用到更广泛地应用到其他宏观经济或者企业生产情况或工业设备管理的预测问题。

[1]张迎春，肖冬荣，赵远东.基于时间序列神经网络的气象预测研究. 武汉理工大学学报.2003，4（29）

[2] 宋仙磊，刘业政，陈思凤，许波.二重趋势时间序列的灰色组合预测模型.计算机工程与应用. 2011，47（8）

[3]徐国祥主编：《统计预测与决策》，上海财经大学出版社2010年版

[4] 宋玉强.人工神经网络在时间序列预测中的应用研究. 西安建筑科技大学硕士论文.2005

[5]戴丹. BP神经网络用于市场预测的研究. 武汉理工大学硕士论文.2006

[8] 颜杰，相丽驰，方积乾.灰色预测模型及SAS实现[J].中国卫生统计，2006年，第23卷第1期：75-85页

[9]蔡建平，朱秀萍等，SAS社会统计实用教程[M].清华大学出版社，2006年，105-107页 [10] 张晓峒.EViews使用指南与案例 [M].北京：机械工业出版社，2007.2.

[11] 牛东晓，乞建勋，刑棉.二重趋势性季节型电力负荷预测组合灰色神经网络模型[J].中国管理科学，2001.12（第9卷第6期）.

[12] 胡江波，刘金兰.回归_时间序列组合模型在石油企业费用预算编制中的应用[N].沈阳理工大学学报，2006.2(第25卷第1期).

[13] 王洪涛，柯跃.基于GMDH和BP的组合预测模型在工业产品销售收入中的应用[J].技术与市场，2010（第17卷第2期）.

[14]赵喜仓，周作杰. 基于SARIMA 模型的我国季度GDP 时间序列分析与预测[J]. 统计与决策，2010,(22).

附录1原始数据

附录3.GM（1,1）拟合结果残差检验（消除季节影响后）

附录4 GM（1,1）模型拟合趋势项

附录5.单一GM（1，1）模型SAS代码

data a1; input xt@@; t=_n_; yt+xt; index=1;

zt=-(yt+LAG(yt))/2; datalines;

46441.36032 46170.59515 46048.23002

44587.22996 49455.23588 49151.39902 49213.95758 47869.54894 52553.82147 52410.49038 52289.61565 52601.24142 58125.21977 58105.42191 57850.81728 59359.0815

65037.60352 64702.74268 63938.61548

64648.40121 70870.4977

70880.04594 70164.54864

proc gplot data=a1; plot xt*t;

symbol c=black v=star i=join; run; proc iml; use a1;

read all var{zt index} into B where (zt^=.); read all var{xt} into Yn where (zt^=.); ahat=inv(B`*B)*B`*Yn; ahatt=ahat`;na={a u};

create a2 from ahatt[colname=na]; append from ahatt; quit; data a3; set a2;index=1; data a4;

set a1;if _n_=1;xt0=xt; keep xt0 index; data a5;

merge a1 a3 a4;by index; if _n_=1 then xp=xt; else do

yt1=(xt0-u/a)*exp(-a*(t-1))+u/a; yt0=(xt0-u/a)*exp(-a*(t-2))+u/a; xp=yt1-yt0; end;

error=abs(xp-xt); rerror=error/xt*100;

drop yt index zt a u xt0 yt1 yt0; proc means data=a5; var xt; var error; output var=val; run;

ODS rtf file='f:/aa.rtf'; proc print data=a5; run;

ods rtf close;

data a6; /*建立GM（1,1）残差模型*/ set a5; et+error;

kt=-(et+lag(et))/2; proc iml;

use a6;

read all var{kt index} into C where(kt^=.); read all var{error} into Zn where(kt^=.); bhat=inv(C`*C)*

C`*Zn;

bhatt=bhat`;na2={ae ue};

create a7 from bhatt[colname=na2]; append from bhatt; quit; data a8; set a7;index=1; data a9;

set a6;if _n_=1;error0=error; keep error0 index; data a10;

merge a1 a3 a6 a8 a9;by index; if _n_=1 then ep=error; else do;

if t>=3 then d=1; else d=0;

yte1=(xt0-u/a)*exp(-a*(t-1))+u/a+d*(-ae)*(2936.84-ue/ae)*exp(-a*(t-2)); yte0=(xt0-u/a)*exp(-a*(t-2))+u/a+d*(-ae)*(2936.84-ue/ae)*exp(-a*(t-3)); xpe=yte1-yte0; end;

error2=abs(xpe-xt); rerror2=error2/xt*100;

drop yt index zt a u xt0 yt1 yt0 et kt ae ue error0 ep d yte1 yte0; ODS rtf file='f:/aa.rtf'; proc print data=a10; run;

ods rtf close;

proc means data=a10; /*残差检验*/ var xt; var error2; output var=val2; run;

data a11; /*关联度检验和后验差检验*/ set a10;

n=(14.5879387+0.5*25195.58)/(error2+0.5*25195.58); proc means data=a11; var n;

output mean=avgn; run;

data p; /*预测*/ input k@@;

cards; 69 70 71 72 73 74 75 76 ;

data pre;

set p;

xps=517136.3317*exp(0.034202*(k-1))-504463.4817+2825.585429*exp(0.016111*(k-2)); abc=lag1(xps);

do;

if k=69 then

xpre=xps-(517136.3317*exp(0.034202*67)-504463.4817+2825.585429*exp(0.016111*66)); else xpre=xps-abc;

end;

drop xps abc;

ODS rtf file='f:/aa.rtf';

proc print data=pre;

run;

ods rtf close;

附录6. GM（1，1）趋势项拟合SAS代码

data a1;

input xt@@; t=_n_; yt+xt;

index=1;

zt=-(yt+LAG(yt))/2;

;

2011年全国大学生统计建模大赛论文二等奖第五名

相关内容

热门内容

标签