SAS在时间序列分析中的应用

时间序列分析我们查询中国统计局网站，得到2008-2012年全国旅客运输量（单位：万人）。数据如下：datenumberdatenumberdatenumberJan-08186900Sep-092231999May-111442798Feb-08386465Oct-092502024Jun-111723719Mar-08583536Nov-092744358Jul-112014344Apr-08776252Dec-092977374Aug-112311493May-08976847Jan-10254630Sep-112611184Jun-081171441Feb-10539287Oct-112924965Jul-081371459Mar-10805047Nov-113220005Aug-081573119Apr-101067262Dec-113517042Sep-081774113May-101333855Jan-12320017Oct-08214452Jun-101593075Feb-12637323Nov-082186100Jul-101862716Mar-12939672Dec-082396031Aug-102137100Apr-121245066Jan-09213183Sep-102419615May-121557178Feb-09424989Oct-102711737Jun-121865719Mar-09631200Nov-102981391Jul-122176606Apr-09705047Dec-103278775Aug-122494741May-091234301Jan-11284044Sep-122815605Jun-091476702Feb-11585065Oct-123151315Jul-091726529Mar-11868311Nov-123469603Aug-091977861Apr-111151719Dec-123789855表1:2008-2012全国旅客运输量第一步：平稳性的判断作出该序列的时序图。SAS代码如下：Procgplotdata=gyf.yunshu1;/*对数据集销售绘制连线图*/plotnumber*date;/*指出以变量data为横坐标，变量number为横坐标*/symboli=splinev=starci=greencv=red;/*指出用样条插值来连线*/run;/*星号标记点，点的颜色为红色，连线的颜色为绿色*/该序列的时序图如图1所示：从序列图中，我们很容易发现该序列式非平稳序列，具有趋势性和周期性。因此，为了消除趋势性与周期性的影响，下面我们会对原始数据分别进行1步和12步差分。图1：2008-2012全国旅客运输量时序图第二步：平稳化序列对原始数据分别进行1次差分和12次差分，并画出差分后的序列图，SAS代码如下：datagyf.yunshu2;/*建立差分后的数据集yunshu2*/setgyf.yunshu1;/*读入数据集yunshu1中的观测*/y=dif12(number);/*对原始数据进行12次差分*/z=dif(y);/*对原始数据进行1次差分*/run;/**/procgplotdata=gyf.yunshu2;/*对数据集销售绘制连线图*/plotz*date;/*指出以变量date为横坐标，变量number为横坐标*/symboli=splinev=starci=greencv=red;/*指出用样条插值来连线*/run;/*星号标记点，点的颜色为红色，连线的颜色为绿色*/差分后的序列图如图2所示：图2：原数据差分后的序列图从图2中我们可以看到差分后的序列基本平稳，且均值在0附近。第三步：模型定阶根据差分后序列的自相关函数和偏相关函数，我们可以确定模型的阶数。编写SAS程序如下：Procarimadata=gyf.yunshu1;Identifyvar=number(1,12)nlag=12;/*对原数据number进行1步12步差分后*/run;/*的数据求自相关函数和偏自相关函数nlag步长为12*/自相关函数如图3所示：图3：旅客运输量差分后自相关图偏自相关函数如图4所示：图4：旅客运输量差分后偏自相关图从图3中我们可以看出，差分后的数据的自相关函数呈现1步延迟之后都减小到两倍的标准误以内，且趋近于0，具有一阶截尾的特点。从图4中可以看到偏自相关函数在第一步之后数据均落在两倍的标准误以内，可以认为具有截尾性，然后数据并没有最终趋近于0，而表现出拖尾性。根据自相关函数和偏自相关函数的判断，我们分别选择AR(1)模型，MA(1)模型，以及ARMA(1，1)模型进行拟合。编写SAS程序如下：Procarimadata=gyf.yunshu1;Identifyvar=number(1,12)nlag=12;Estimatep=1noconstant;/*拟合AR(1)模型，且不对常数项进行拟合*/Estimateq=1noconstant;/*拟合MA(1)模型，且不对常数项进行拟合*/estimateq=1p=1noconstant;/*拟合ARMA(1)模型，且不对常数项进行拟合*/run;下面对各型进行比较，见表2：模型AIC准则值BIC准则值AR(1)1333.4011335.251MA(1)1330.7731332.623ARMA(11331.1481334.848表2准则函数值比较从表2中，我们可以看到MA(1)模型的AIC准则值与BIC准则值都最小。因此我们可以认为MA(1)模型要优于AR(1)模型与ARMA(1，1)模型。因此选择用MA(1)模型进行拟合。第四步：参数估计与模型检验第三步中我们确定了使用MA（1）模型进行拟合。求得参数如图5所示：图5MA（1）模型的参数估计和拟合统计量从图5中，可以看到参数通过检验，若原序列为tX，且记12ttZX，得到MA（1）模型为：10.64133tttZ(1)残差的自相关检验如图6所示：图6残差的自相关检验从图6中可以看出，对MA（1）模型的残差序列进行自相关检验的p值均明显大于显著性水平=0.05，因而接受残差无自相关性的假设。这说明MA（1）模型从原数据序列中提取出足够多的信息，建模后的残差不存在自相关性，因此这是一个合理的模型。（1）式写为原始表达式为：121(1)(1)0.64133tttBBX（2）（2）式还可以写为：1121310.64133ttttttXXXX（3）第五步：模型预测在第四步我们已经确定了MA(1)模型的参数，最终确定模型为（3）式。下面利用所见模型进行预测，并与实际数据进行比较。编写SAS代码如下：Procarimadata=gyf.yunshu1;Identifyvar=number(1,12)nlag=12;Estimateq=1noconstant;forecastlead=12id=dateout=resultinterval=month;/*计算12步预测值*/run;/*以数据集yunshu1中变量date作为标识变量，将预测值输出到数据集result中，时间间隔为1个月*/预测结果见图7：图7预测2013年全国旅客运输量我们查阅中国统计局网站，得到2013年全国全年旅客运输量，实际值与预测值见表3。利用预测值的输出数据集，我们绘制出了样本观测值与预测值的曲线图。编写代码如下：procgplotdata=result;/*对观测集result画图*/plotnumber*date=1forecast*date=2/overlay;/*分别以number为纵坐标，date为横坐标，以forecast为纵坐标，date为横坐标，并且将两个图形重叠*/symbol1c=bluei=splinev=star;/*指出第一个图中以插值方式连线，点的符号位型号，点和线的颜色为蓝色*/symbol2c=redi=splinev=circle;/*指出第二个图中以插值方式连线，点的符号位圆圈，点和线的颜色为红色*/run;提交上述程序后，得到曲线图如图8所示：表3预测值与真实值比较时间真实值预测值Jan-13332321553603.2Feb-13672346870909.2Mar-139966971173258Apr-1313202631478652May-1316474261790764Jun-1319721392099305Jul-1323025642410192Aug-1326400772728327Sep-1329797203049191Oct-1333324393384901Nov-1336671393703189Dec-1340190614023441图8样本观测值与预测值的曲线图图8中，用星号代表样本观测值，并用蓝色曲线连接，用圆圈表示预测值，并用红色曲线连接。结合表3，从图8中，我们可以看出，利用模型作出的预测值与样本的原始数据非常的接近。这也说明可拟合时间序列的效果较好。