时间序列分析上机指导

1、上机指导第五章5.8.1 拟合ARIMA模型由于ARMA模型是ARIMA模型的一种特例，所以在SAS系统中这两种模型的拟合都放在了ARIMA过程中。我们已经在第3章进行了ARMA模型拟合时介绍了ARIMA过程的基本命令格式。再次以临时数据集example5_1的数据为例介绍ARIMA模型拟合与ARMA模型拟合的不同之处。 data example5_1; input x; difx=dif(x); t=_n_; cards;1.05 -0.84 -1.42 0.20 2.81 6.72 5.40 4.38 5.52 4.46 2.89 -0.43 -4.86 -8.54 -11.54 -16.

2、22-19.41 -21.61 -22.51 -23.51 -24.49 -25.54 -24.06 -23.44-23.41 -24.17 -21.58 -19.00 -14.14 -12.69 -9.48 -10.29-9.88 -8.33 -4.67 -2.97 -2.91 -1.86 -1.91 -0.80 proc gplot; plot x*t difx*t; symbol v=star c=black i=join;run;输出时序图显示这是一个典型的非平稳序列。如图5-49所示图5-49 序列x时序图考虑对该序列进行1阶差分运算，同时考察查分后序列的平稳性，在原程序基础上添加相

3、关命令，程序修改如下： data example5_1; input x; difx=dif(x); t=_n_; cards;1.05 -0.84 -1.42 0.20 2.81 6.72 5.40 4.38 5.52 4.46 2.89 -0.43 -4.86 -8.54 -11.54 -16.22-19.41 -21.61 -22.51 -23.51 -24.49 -25.54 -24.06 -23.44-23.41 -24.17 -21.58 -19.00 -14.14 -12.69 -9.48 -10.29-9.88 -8.33 -4.67 -2.97 -2.91 -1.86 -1.

4、91 -0.80 proc gplot; plot x*t difx*t; symbol v=star c=black i=join; proc arima; identify var=x(1); estimate p=1; forecast lead=5 id=t ; run;语句说明：（1）DATA步中的命令“difx=dif(x);”,这是指令系统对变量x进行1阶差分，差分后的序列值赋值给变量difx。其中dif（）是差分函数，假如要差分的变量名为x，常见的几种差分表示为：1阶差分：dif（x）2阶差分：dif（dif（x）k步差分：difk（x）（2）我们在GPLOT过程中添加绘制了一

5、个时序图“difx*t”，这是为了直观考察1阶差分后序列的平稳性。所得时序图如图5-50所示。图5-50 序列difx时序图时序图显示差分后序列difx没有明显的非平稳特征。（3）“identify var=x（1）；”，使用该命令可以识别查分后序列的平稳性、纯随机性和适当的拟合模型阶数。其中x（1）表示识别变量x的1阶差分后序列。SAS支持多种形式的差分序列识别：var=x（1），表示识别变量x的1阶查分后序列xt；var=x（1,1），表示识别变量x的2阶查分后序列2xt；var=x（k），表示识别变量x的k步差分后序列kxt；var=x（k，s），表示识别变量x的k步差分后，再进行s步查

6、分后序列skxt。识别部分的输出结果显示1阶查分后序列difx为平稳非白噪声序列，且具有显著的自相关系数不截尾、偏自相关系数1截尾的性质。（4）“estimate p=1；”对1阶差分后序列xt拟合AR（1）模型。输出拟合结果显示常数项不显著，添加或修改估计命令如下：estimate p=1 nonit；这是命令系统不要常数项拟合AR(1)模型，拟合结果显示模型显著且参数显著。如图5-51所示。图5-51 序列difx模型拟合结果输出结果显示，序列xt的拟合模型为ARIMA（1,1,0）模型，模型口径为：xt=艾普龙t/1-0.66933B或等阶记为：xt=1.66933xt-1-0.6693

7、3xt-2+艾普龙t（5）“forecast lead=5 id=t；”，利用拟合模型对序列xt作5期预测。一、 建立数据集，绘制时序图data example5_2; input x; lagx=lag(x); t=_n_; cards; 3.03 8.46 10.22 9.80 11.96 2.83 8.43 13.77 16.18 16.84 19.57 13.26 14.78 24.48 28.16 28.27 32.62 18.44 25.25 38.36 43.70 44.46 50.66 33.01 39.97 60.17 68.12 68.84 78.15 49.84 62.2

8、3 91.49 103.20 104.53 118.18 77.8894.75 138.36 155.68 157.46 177.69 117.15 ; proc gplot data=example5_2;plot x*t=1;symbol1 c=black i=join v=star; run; 输出时序图如5-52所示。 图5-52 序列x时序图时序图显示，序列X有一个明显的随时间线性递增的趋势，同时又有一定规律性的波动，所以不妨考虑使用误差自回归模型拟合该序列的发展。2、 因变量关于时间的回归模型proc autoreg data=example5_2;model x=t/ dwpro

9、b;run;语句说明：（1) “proc autoreg data=example5_2;”指令SAS系统对临时数据集example5_2进行回归程序分析。（2) “model x=t/ dwprob;”指令SAS系统以变量t作为自变量，变量x作为因变量，建立线性模型：并给出残差序列 DW检验统计量的分为点。本例中，序列x关于变量t的线性回归模型最小二乘估计输出结果如图5-53所示。图5-53 序列x关于变量t的线性回归模型最小二乘估计结果本例输出结果显示，DW统计量的值等于0.7628，输出概率显示残差序列显著正相关。所以应该考虑对残差序列拟合自相关模型，修改AUTOREG程序如下：proc

10、 autoreg data=example5_2; model x=t/nlag=5 backstep method=ml; run;Model语句是指令系统对线性回归模型 的残差序列显示延迟5阶的自相关图，并拟合延迟5阶自相关模型，特别注意，SAS输出的自回归模型结构为：即输出的自相关回归参数值与我们习惯定义的自回归参数值相差一个负号。 由于自相关延迟阶数的确定是由我们尝试选择的，所以nlag得阶数通常会指得大一些。这就导致残差自回归模型中可能有部分参数不显著，因而添加逐步回归选项backstep,指令系统使用逐步回归的方法筛选出显著自相关因子，并使用极大似然的方法进行参数估计。输出如下四方

11、面的结果：1因变量说明如图5-54所示：图5-54 因变量说明2普通最小二乘估计结果该部分输出信息包括差平方和（SSE)、自由度（DFE)、均方误差（MSE)、根号均方误差（Root MSE)、SBC信息量。回归部分相关系数平方（regress rsquare)、总的相关系数平方（totel rsquare),DW统计量（durbin watson)及所有待参数的自由度、估计值、标准差、t值和统计量的P值。如图5-55所示。图5-55 普通最小二乘估计结果3回归误差分析该部分共输出四方面的信息：残差序列自相关图、逐步回归消除的不显著项报告、初步均方误差（MSE)、自回归参数估计值。本例该部分输

12、出结果如图5-56所示。图5-56 自回归误差分析输出结果本例输出的残差序列自相关图显示残差序列有非常显著的1阶正相关性。逐步回归消除报告显示除了1阶的序列值显著自相关外，延迟其他阶数的序列值均不具有显著的自相关性，因此延迟2阶-5阶的自相关项被剔除。初步均方误差为4、最终拟合模型该部分包括三方面的汇总信息：收敛状况、极大似然估计结果和回归系数估计。本例该部分输出结果如图5-57所示。图5-57 最终拟合模型输出结果本例得到最终拟合模型为：为了得到直观的拟合效果，我们可以利用OUTPUT命令将拟合结果存入SAS数据集中，并对输出结果作图，相关命令如下：proc autoreg data=exa

13、mple5_2; model x=t/nlag=5 backstep method=ml; output out=out p=xp pm=trend;proc autoreg data=example5_2; model x=t/nlag=5 backstep method=ml noint; output out=out p=xp pm=trend; proc gplot data=out; plot x*t=2 xp*t=3 trend*t=4 / overlay ; symbol2 v=star i=none c=black; symbol3 v=none i=join c=red w=

14、2 l=3; symbol4 v=none i=join c=green w=2;run; 语句说明：“output out=out p=xp pm=trend;”，该命令是指令系统将部分结果输入临时数据集OUT,选择输出的第一个信息为整体模型的拟合值（P值），该拟合变量取名为XP；选择输出的第二个信息为线性趋势拟合值（PM值），还可以选择R选项输出拟合残差项，本例不要求输出此项。输出图像如图5-58所示。图5-58 拟合效果图3、 延迟因变量回归模型proc autoreg data=example5_2; model x=lagx/lagdep=lagx; run；语句说明：（1） 首先在

15、DATA步中添加命令“lagx=lgax(x);”，该语句指令系统使用延迟函数生成序列x的1阶延迟序列，并将该序列赋值给变量lagx，即（2） “model x=lagx/lagdep=lagx; ”指令 系统建立带有延迟变量的回归模型并通过LAGDEP选项指定被延迟的因变量名。本例输出结果如图5-59所示。图5-59带延迟因变量回归分析结果由于带有延迟因变量，所以这种场合在回归模型估计结果中输出的是durbin h统计量。本例中durbin h统计量的分布函数达到0.3853，这表示残差序列不存在显著的相关性，不需要考虑对残差序列继续拟合自回归模型。在注意参数检验结果，如图5-60所示。图5

16、-60 参数估计结果在显著性水平默认为0.05的条件下，截距项不显著（P值大于0.05），所以可以考虑在模型拟合命令中增加NOINT选项。最后输出拟合结果，并绘制拟合时序图相关命令如下：proc autoreg data=example5_2; model x=lagx/lagdep=lagx; model x=lagx/lagdep=lagx noint; output out=out p=xp; proc gplot data=out; plot x*t=2 xp*t=3 / overlay;symbol2 v=star i=none c=black; symbol3 v=none i=j

17、oin c=red w=2 l=3; run;模型拟合部分输出的信息表示最终的拟合模型为：输出的拟合效果图如图5-61所示。图5-61带有延迟因变量的回归模型拟合效果图本例中由于没有残差自回归项，最终拟合值就是趋势值，所以只需绘制其中之一即可。介绍GARCH模型拟合，相关程序如下：Data example5_3;Input x;T=_n_;Cards;10.77 13.30 16.64 19.54 18,97 10.52 24.3623.51 27.16 30.80 31.84 31.63 32.68 34.9033.85 33.09 35.46 35.32 39.94 37.47 35.24

18、33.03 32.67 35.20 32.36 32.34 38.45 38.1732.14 39.70 49.42 47.86 48.34 63.50 63.5667.61 64.59 66.17 67.50 76.12 79.31 78.8581.34 87.06 86.41 93.20 82.95 72.96 61.1061.27 71.58 88.34 98.70 97.31 97.17 91.1780.20 85.12 81.40 70.87 57.75 52.35 67.5087.95 85.46 84.55 98.16 102.42 113.02 119.95122.37 126

19、.96 122.79 127.96 139.20 141.05 140.87137.08 145.53 145.59 134.36 122.54 106.92 97.23110.39 132.40 152.30 154.91 152.69 162.67 160.31142.57 146.54 153.83 141.81 157.83 162.79 142.07139.43 181.49 166.84 154.28 150.12 165.17 170.32；proc gplot data= example5_3:plot x*t=1;symboll c=black i=jion v=star;p

20、roc autoreg data=example5_3;model x=t/nalg=5 dwprob archtest;model x=t/nalg=2 noint garch= (p=1,q=1);output out=out p=p residual=residual lcl=lcl ucl=ucl cev=cev;data out;set out;l95=-1.96*sqrt(51.42515);u95=1.96*sqrt(51.42515);Lcl_GARCH=-1.96*sqrt(cev);Ucl_GARCH=1.96*sqrt(cev);Lcl_p=p-1.96*sqrt(cev

21、);Ucl_p=p+1.96*sqrt(cev);proc gplot data=out;polt residual *t=2l95*t=3LCL-garch*t=4u95*t=UCL-garch*t=4/ovctlaypolt X *t=5 lcl * t =3LCL-p * t = 3UCL-p * t = 4/overlay;symbol2 c = green i = needle v = none;symbol3 c = black i =join v = none w 2 l =2;symbol4 c = red i =join v none;symbol5 c = green i

22、= join v = none;run;该时序输出图如图5-62所示图5-62 序列时序图时序图显示序列具有显著线性递增趋势，且波动幅度随时间递增，所以考虑使用AUTOREG过程建立序列xt,关于时间t的线性回归模型,并检验残差序列的自相关和易方差性，如果检验结果显示残差序列具有显著之相关性，建立参差回归模型；如果残差序列具有显著的异方差性，则需建立条件异方差模型，则要建立条件异方差模型。语句说明：（1）“model x =t/nlag=5 dwprob archtest；”该命令指令系统建立序列xt，关于时间t的线性回归模型，并检验残差系列5阶延迟的自相关性，并输出DW检验的p值，同时对残差

23、序列进行异方差检验。DW检验结果显示残差序列具有显著的正自相关性，如图5-63所示。图5-63 普通最小二乘输出结果图5-64 残差序列自相关图参数估计结果显示回归模型常数截距项不显著，如图5-65所示。图5-65 线性回归模型参数估计结果异方差检验结果显示，残差序列具有显著的异方差性，且具有显著的长期相关性，如图5-66所示。图5-66 异方差检验结果（2）“modelx=t/nlag=2 noint garch=（p=1，q=1）；”综合考虑残差序列自相关性和异方差性检验结果，尝试拟合无回归常数项的AR(2)-GARCH(1，1)模型。模型最终拟合结果如图5-67所示。参数检验结果显示除GARCH(1，1)模型中的常数项不显著外，其他模型均显著，整个模型的R2高达0.9954，且正态性检验不显著（p值为0.3106），这与假定GARCH的残差函数t/ht服从正态分布相吻合，所以可以认为该模型拟合成功。最