第28章如何用SAS实现时间序列分析所谓时间序列，就是将某一指标在不同时间上的不同数值，按照时间先后次序排列而成的数列，这种数列由于受到各
种偶然因素的影响，往往表现出某种随机性，彼此之间存在统计上的依赖关系。因此，可以通过对时间序列的研究来认识所研究系统的结构特征（如
波动的周期、振幅、趋势的种类），揭示其运行规律，进而用以预测、控制未来行为，修正和重新设计系统。时间序列分析是一种重要的现代统计学
方法，主要有确定性时间序列分析和随机时间序列分析方法。另外，在实际问题中会遇到这样的情况，一个时间序列目前的表现，不仅受过去行为的
影响，而且与另一个时间序列相关。某地区经济增长的情况，不仅与过去有关，还受到投资、政策等因素的影响，进行多个因素对结果变量的影响要
进行多重时间序列分析。28.1求和自回归滑动平均模型（integratedautoregressivemovingaver
agemodel,ARIMA）原理概述在SAS软件中，采用ARIMA过程进行分析和预测等间隔的时间序列。ARIMA过程提供了一个
综合的工具包来进行模型的识别、参数估计及预测。ARIMA模型通过其自身的过去值、过去误差、其他时间序列的当前值和过去值的线性组合来
预测响应时间序列。其中，差分具有强大的确定性信息提取能力，许多非平稳序列差分后会显示出平稳序列的性质，称该非平稳序列为差分平稳序列
。对该种序列常用的方法就是本章介绍的齐次非平稳序列，简记为ARIMA(p,d,q)模型。ARIMA(p,d,q)模型的结构为：对d
阶齐次非平稳序列而言，{}是一个平稳序列，设其适合ARIMA（p，q）模型，即或表达为其ARIMA模型的构建由3个阶段组成：（
1）模型的识别阶段：在识别阶段，可通过identify语句识别差分数、计算自相关、偏自相关、逆相关、互相关系数。还可进行平稳性检验
和模型阶数的识别。另外，还可同时写多个identify语句，用以寻找模型的适合形式。另外，需要指出的是，时间序列的每个时间点只有一
个观测值，是无法采用数理统计学的方法进行推断的。因此，如果时间序列不平稳，就其均值来说，就是每个都不相等，而每个都只能用一个观测去
估计，这样做效果肯定不佳。但是，若时间序列是平稳的，就可以将整个时间序列看做相同的总体，也就是将它们的观测值看作从同一个总体中抽出
的样本，从而采用经典数理统计学的方法进行处理。因此，在时间序列建立的过程中，一般首先将序列转换为平稳时间序列。在实际数据分析中，可
首先绘制时序图、自相关图，还可进行逆序检验或游程检验法（最简单的判断随机性的方法），观察序列是否是平稳序列，游程检验亦称“连贯检验
”，是根据样本标志表现排列所形成的游程的多少进行判断的检验方法。游程检验的原则：如果序列为真随机序列，那么游程的总数应该不太多也不
太少。如果游程的总数极少，就说明样本缺乏独立性，内部存在一定的趋势或者结构，这可能由于观察值间不独立，或者来自不同的总体。如果样本
间存在大量游程，则可能有系统的短周期波动影响观察结果。同样认为序列非随机。所以，可以用游程的个数来检验样本的随机性，或总体的分布特
征。游程检验的本质：首先，变量的类型必须为二分变量，例如性别变量，只有二个数组成的变量。然后，游程检验的分析目的是用于判断观察值的
顺序是否随机。这一点非常重要，因为，许多遇到的实际问题中并不只是使研究者关心分布的位置或者形状，也包括样本的随机性。如果样本不是从
总体中随机抽取的，则所做的任何推断都将没有价值。游程检验是最简单的判断随机性的方法。（2）模型的估计和检验阶段：可通过estima
te语句进行模型参数的估计及其显著性检验。并且，estimate语句还可进行模型拟合优度检验、残差的白噪声检验等。对于模型的识别方
法有很多，对于Box-Jenkins的模型识别方法是根据样本自相关、偏相关函数的截尾和拖尾性初步判断时间序列所适合的模型，这种方法
既可以对AR（p）模型、MA（q）模型进行初步识别，且方法简单，但精确度欠缺，尤其是当样本序列未达到足够长度时，其精度不够理想。另
外，还可采用最佳准则函数定阶法（如采用AIC准则、SBC准则选取阶数）以及残差方差图定阶法。表28-4平稳时间序列的自相关函数
和偏自相关函数的形式模型AR(p)MA(q）ARIMA(p,q)ACF（自相关函数）拖尾截尾拖尾PACF（偏自相关函数）截尾拖尾拖
尾注：当d=0时，ARIMA(p,d,q)模型就是ARIMA模型另外，研究中的数据属于随时间推移而收集的数据，不同时期的观察
往往是有关联的，或者称为自相关的。并且，序列之间也可能存在相关性，称为互相关。因此，对于传递函数时间序列数据，传递函数模型的识别和
干扰模式的识别，模型识别采用的是计算机输入和输出序列的互相关函数，根据互相关函数显著不为零的期数确定p和q值。平稳的序列的自相关图
和偏相关图不是拖尾就是截尾。截尾就是在某阶之后，系数都为0，怎么理解呢，看上面偏相关的图，当阶数为1的时候，系数值还是很大，0.9
14二阶长的时候突然就变成了0.050后面的值都很小，认为是趋于0，这种状况就是截尾。再就是拖尾，拖尾就是有一个衰减的趋势，但是不
都为0。自相关图既不是拖尾也不是截尾。以上的图的自相关是一个三角对称的形式，这种趋势是单调趋势的典型图形。下面是通过自相关的其他功
能如果自相关是拖尾，偏相关截尾，则用AR算法如果自相关截尾，偏相关拖尾，则用MA算法如果自相关和偏相关都是拖尾，则用ARMA算法，
ARIMA是ARMA算法的扩展版，用法类似（3）模型的预测阶段：可通过forecast语句进行时间序列的预测及置信区间的计算。值得
注意的是：对于多重时间序列，需要有对应的的数值，才能对进行预报。28.2数据与分析【分析与解答】研究中的数据属于时间序列数据，对
该类型数据进行建模和预测，建立时间序列模型，需要几个步骤：模型的识别、模型的定阶、模型参数的估计。【SAS程序】ODSHTML;
DATAMXWTTJXS28_1;INPUTx@@;y=dif12(x);/对原始数据进行12次差分/z=dif(y);
/对原始数据进行1次差分/t=intnx(''month'',''01jan2000''d,_n_-1);/INTNX(interv
al,from,n)计算从from开始经过n个in间隔后的SAS日期。_n_是data步的自动变量，_n_表示观测的序/FOR
MATtmonyy.;/format用来设定输出数据的格式，将时间t用monyy.格式输出/CARDS;2962.9280
4.92626.62571.52636.92645.22596.92636.33136.93347.33107.83680333
2.83047.12876.12820.92929.62908.72851.42889.42854.33029.33421.740
33.33324.43596.13114.83052.23202.13158.83096.63143.73422.43661.9
3733.14404.43907.43706.43494.83406.93463.33576.93562.13609.63971.
84204.44202.74735.74569.44211.44049.84001.84166.14250.74209.2426
2.74717.74983.24965.65562.55300.95012.24799.14663.34899.24935493
4.95040.85495.25846.659096850.46641.66001.95796.75774.66175.66057
.86012.26077.46553.66997.76821.77499.27488.37013.76685.86672.571
57.570266998.27116.67668.482638104.79015.3;PROCGPLOT;/对数据集绘
图/PLOTxt;/指出以变量名t为横坐标，变量名x为纵坐标绘图/SYMBOLi=joinv=dot;/i=joi
n表示用直线连接各点，v=dot表示用圆点表示点/PLOTzt;SYMBOLi=joinv=dot;RUN;PROC
ARIMA;/表示调用时间序列分析的过程/IDENTIFYVAR=XNLAG=12;/var=X表示对x变量进行识别，输
出均值，标准差等描述性统计量；输出自相关，逆自相关和偏自相关系数及图（平稳性检验和参数定阶基础）；输出白噪声检验结果//nla
g=12表示指明计算自相关系数和互相关系数过程中需要考虑的延迟阶数。如果不特别指定nlag的阶数，计算机默认的输出阶数是min(2
4,n/4)/RUN;PROCARIMA;IDENTIFYVAR=X(1,12)NLAG=12;/表示对原始数据X进行
1次差分和12步差分后变量进行识别（通常用于有趋势和有周期的序列识别），并求数据自相关函数和偏相关函数nlag步长为12/RUN
;PROCARIMA;IDENTIFYVAR=X(1,12)NLAG=12MINICp=(0:5)q=(0:5);/
MINICp=(0:5)q=(0:5)在自相关阶数跑遍0-5，移动平均阶数跑遍0-5的范围内，寻找AIC最小的模型阶数（一种
傻瓜型定阶方法）/ESTIMATEp=4q=2method=CLS;/p=自相关阶数，q=移动平均阶数/ESTIMA
TEp=1q=1method=CLS;/CLS表示条件最小二乘估计/FORECASTLEAD=5ID=tOUT=r
esults;/forecast给出估计值及置信区间；给出指定期数的预测值及置信区间；选择是否输出预测结果//命令格式：FO
RECASTlead=预测期数id=时间标示interval=单位时间间隔alpha=aout=输出结果数据集/RUN
;ODSHTMLCLOSE;【主要输出结果及解释】第一步，GPLOT过程用于绘制时序图。如图所示，该序列具有明显的趋势性和季节
性，为非平稳序列。第二步，为了进一步确定其非平稳性，采用第一个ARIMA过程来分析确定其非平稳性。生成序列的自相关图，见下图，从自
相关图可以看出自相关函数出现非常缓慢的衰减，表明该序列是不平稳的。为了使得序列稳定，需要消除序列的季节性和趋势性，故对原始序列进行
季节差分和1次差分。即SAS中第二个ARIMA过程进行该项工作。结果见下图。可见，经过12步和1次差分后，序列基本平稳。绘制出{}
的自相关函数图和偏相关函数图，如图28-4所示。可知{}序列的自相关和偏相关函数均拖尾，因此对序列建立ARMA(p,q)模型。第四
步，利用第三个ARIMA过程，通过计算不同阶数的自回归和移动平均所对应的BIC，找到最小的BIC所对应的p，q值，并将其作为最后建
模的候选阶数。结果如下图所示，可见，p=4,q=2时BIC(4,2)=9.776787最小，故最后选择ARMA(4,2)模型。Mi
nimumTableValue:BIC(4,2)=9.776787第五步，对ARMA(4,2)模型进行参数估计和假设检验
。通过estimate语句实现。如下图所示：初始模型的参数估计和检验的结果由上图可知，部分参数未通过检验，所以要去掉没有意义的变量
，尝试重新拟合。最后采用第二个estimate语句进行最终的估计，结果如下图所示。记作，可见最终模型的表达为：。第六步，对最终模型
进行残差检验，结果如图28-8所示。可见，残差序列是白噪声序列，说明该模型是合适的。白噪声序列，是指白噪声过程的样本实称，简称白噪
声。随机变量X（t）（t=1，2，3……），如果是由一个不相关的随机变量的序列构成的，即对于所有S不等于T，随机变量Xt和Xs的协
方差为零，则称其为纯随机过程。对于一个纯随机过程来说，若其期望为0，方差为常数，则称之为白噪声过程。之所以称为白噪声，是因为它和白
光的特性类似，白光的光谱在各个频率上有相同的强度，白噪声的谱密度在各个频率上的值相同。残差白噪声检验结果第七步，利用模型进行预测，
是通过forecast命令实现的。预测的结果存入了数据集results。【分析与解答】研究中的数据属于多重时间序列数据[因变量为Y
(obsend)，除隐含的自变量t以外，还有另外4个以显变量形式呈现的自变量z1~z4（即dhd、dhdlag、wind、xwee
kend）]。程序名为mxwttjxs28_02.sas，具体如下：【SAS程序】odshtml;DataMXWTTJXS28
_2;inputobsenddhddhdlagwindxweekend@@;time=_n_;cards;227323
01212363132812283031802523430802382834120195242880261392416040662
39161428626216136649621604266249160402536280301385312024730388031
04630161401594614136350591603434650204486046120450576014035245578
03174645813764746121348384716029837388034245378032837452002893937
16025938391212883538161263323580258313214029937311603304237120290
35422002253035161378503018144759501203864259220416544216046961541
20439586114037848581813494248181363454212043156451204245256160384
47528025233478019425331611842325161192252312021330251202333130802
12303112019226301602133226121296413216133346418026633468028038331
8038652382204155752180;procgplot;plotobsendtime;symboli=join
v=dotc=blue;run;procarima;identifyvar=obsendnlag=15;run;proc
arimadata=MXWTTJXS28_2;/调用ARIMA过程/identifyvar=obsendcrosscor
r=(dhddhdlagwindxweekend);/crosscorr是计算各原因变量与结果变量之间的互相关函数//
estimatemethod=mlinput=(dhddhdlagwindxweekend)plot;/estim
atep=(17)noconstantmethod=mlinput=(dhddhdlagwindxweekend
)plot;/noconstant是不输出常数项/run;quit;odshtmlclose;【主要输出结果及解释】对序
列进行识别。首先，判定序列的稳定性，由各时点天然气的输送量趋势图可以看出，序列基本稳定。图28-9不同时间点天然气输送趋势图经过
尝试，在模型中，不仅要加进当前的DHD值，而且要加进DHD的一阶滞后值，首先不考虑天然气输送量的自相关，建立回归模型，通过绘制该模
型残差的ACF和PACF图，用以判断自相关和偏自相关系数。由上可知，残差的滞后一期和七期自相关系数分别为,有统计学意义，可见在构造
模型时，要考虑当天与前一周的数值有关。图28-11是模型估计的结果，可见NUM1~NUM4分别指dhd、dhdlag、wind、x
weekend.需要指出的是，xweekend经统计学检验P=0.0896，但考虑到本例样本量比较小（n=36），故在模型中保留，
亦可将其去掉，重新拟合模型。图28-11参数估计结果残差白噪声检验的结果见图28-12。图28-12残差白噪声检验的结果建立的
回归模型为：天然气输送量。【分析与解答】为估计呼叫次数减少的数量，并预报将来助手的呼叫次数，在建立干预模型之前，先建立单变量时间序
列模型。【SAS程序】odshtml;DATAMXWTTJXS28_3;inputx@@;time+1;cards;3503
39351364369331331340346341357398381367383375353361375371373366382
42940640342942542740940240941940442946342844944446747446343245346
24564745144894754925255275335275225265135645995725875996016116205
79582592581630663638631645682601595521521516496538575537534542538
54754052654855554559464365261664062563763462164165464966269967270
47007117157186526646957047337727167127327557617487487507447317828
10777816840868872811810762634626649697657549162177175162161165171
72178186178178189205202185193200196204206227225217219236253213205
210216218235241;procgplot;plotxtime;symboli=joinv=dotc=blue
;run;procarima;identifyvar=xnlag=12;run;dataMXWTTJXS28_3;set
MXWTTJXS28_3;y=dif(dif12(x));/y=(1-B)(1-B12)Xt/run;procgplot;p
lotytime;symboli=joinv=dotc=blue;run;procarima;identifyvar
=x(1,12)nlag=12;run;dataa;dotime=1to192;output;end;run;data
MXWTTJXS28_3;mergeMXWTTJXS28_3a;iftime<147thens=0;elses=1;
run;procarimadata=MXWTTJXS28_3;identifyvar=s(1,12);estimateq=
1;identifyvar=x(1,12)crosscorr=(s(1,12));estimateq=(12)input=(0$s)noconstantaltparmmaxit=30backlim=-3plot;forecastlead=12interval=month;run;quit;odshtmlclose;/backlim是指设定参数收敛时的方法，altparm是定义总体离散参数估计的方法。需要指出的是，在该程序中，如要对结果变量进行预报，需要事先知道原因变量，才能对结果变量进行预报。/【主要输出结果】不同时间点呼叫次数{}散点图如下：可见{}不平稳。由下图可见，{}经1步和12步差分后，序列基本平稳，但在time=147处出现异常情况。首先，不考虑异常点，由程序得到序列的ACF图和PACF图。因此，根据上图，提出模型：由于1974年3月政策出台，使得在147时间点以后出现了大幅度变化，也就是干预点的出现，使呼叫次数显著下降。所以，考虑干预点的影响，令得到干预模型为：将干预措施引入模型后，拟合结果如下：图28-16干预模型的参数估计结果图28-17残差白噪声检验的结果建立的干预模型为：可见，采取该政策后，平均每个月减少的呼叫次数大约为40000次。由程序得到将来12个月的月平均呼叫次数如下图所示： |
|
