时间序列整合分析(一)汇总

1、1. 缺省值的补足：2. 时序图：（检验平稳性）3. 自相关函数： （检验平稳性）4. 计算标准正态分布的概率：5. 计算标准正态分布的分位数：6. 计算标准 t 分布的概率7. 计算标准 t 分布的分位数8. 计算标准 F 分布的概率9. 计算标准 F 分布的分位数10. 计算标准卡方分布的概率11. 计算标准卡方分布的分位数12. 方差的同齐性检验： 将数据进行适当分组，这里将 4 个分为一组，一共四组Pr>F 的值大于 0.05 故接受 H0，认为各组方差之间没有显著的差异。13.方差的同质性检验：将数据进行适当分组，这里将 4 个分为一组，一共四组根据上结果列出方差分析表：方差来

2、源平方和自由度均方和F值显著性A87837632927924.17误差14037292070186.4总和228210523F 的 p 值小于 0.05 我们认为原始数据方差不同质。14. 序列的白噪声检验（检验纯随机性） ：可以看出， LB（6）=95.84 ，其 p 值小于 0.05 ； LB（12）=190.40 ，其 p 值小于 0.05 ；显然该序列不是白噪声序列，即不是纯随机性序列。（ p 值都大于 0.05 时才是纯随机序列）15. 平稳序列的自相关函数和偏自相关函数的形式：（没有程序的）模型AR(p)MA(q)ARMA(p,q)ACF 自相关拖尾截尾拖尾PACF 偏自相关截尾拖

3、尾拖尾16. 一个例子： (利用平稳序列建模 进行预测)我国 1975-2006 年 GDP的年增长率为下表(数据略) ，对我国 1975-2006 年 GDP 的 年增长率进行建模，并对 2007 至 2011 年我国的 GDP 增长率进行预测。(1) 首先画出我国 1975-2006 年 GDP 增长率的时序图。data ex; input x;t=_n_; cards ;8.7 -1.6 7.6 11.7 7.6 7.8 5.2 9.1 10.9 15.2 13.5 8.8 11.6 11.3 4.13.8 9.2 14.2 14 13.1 10.9 10 9.3 7.8 7.6 8.4

4、 8.3 9.1 10 10.1 10.4 10.7 ; proc gplot ; symbol i =jiontv =dot; plot x*t; run ;从图中直观的可以看出有奇异点(2) 将奇异点看成缺省值，利用以下程序来求缺省点的值： data ex; input x;time=intnx( 'month' , '01jan1975'd ,_n_- 1); format time data;cards ;8.7 . 7.6 11.7 7.6 7.8 5.2 9.1 10.9 15.2 13.5 8.8 11.6 11.3 4.1 3.89.2 14.2

5、 14 13.1 10.9 10 9.3 7.8 7.6 8.4 8.3 9.1 10 10.1 10.4 10.7proc expand data=ex out=ex1;id time; proc print data =ex1;run ;结果可知，缺省点的值为 2.4(3 )利用修正后的数据再进行时序分析，根据以下程序：可以看出 GDP 增长率修正后的数据序列平稳。BIC (5,0)=-0.24488 的值最小，考虑建立 AR(5)模型。(4 )模型的建立data ex; input x;time=intnx( 'month' , '01jan1975'd

6、,_n_- 1); format time year4. ;cards ;arima ; identifyestimate8.7 2.4 7.6 11.7 7.6 7.8 5.2 9.1 10.9 15.2 13.5 8.8 11.6 11.3 4.1 3.8 9.2 14.2 14 13.1 10.9 10 9.3 7.8 7.6 8.4 8.3 9.1 10 10.1 10.4 10.7var =x nlag =12 minic p=( 0: 5) q=( 0: 5);p=5; run从上图中可以看出， 有些参数不显著，我们将其去掉，建立最精干的模型。(其中可以看出，AR1,3 AR1,4

7、 AR1,5 的 p 值远远大于 0.05 )所以， 将 estimatep=5 改为 estimatep=( 1, 2) ，即程序为：data ex; input x;time=intnx( 'month' , '01jan1975'd ,_n_- 1);format time year4. ;cards ;8.7 2.4 7.6 11.7 7.6 7.8 5.2 9.1 10.9 15.2 13.5 8.8 11.6 11.3 4.1 3.89.2 14.2 14 13.1 10.9 10 9.3 7.8 7.6 8.4 8.3 9.1 10 10.1 10

8、.4 10.7arima ; identifyvar =x nlag =12 minic p=( 0: 5) q=( 0: 5);estimate p=( 1, 2); run可见所有的 p 值都小于 0.05 通过了检验。所以模型为：Xt9.58 0.64Xt 1 0.38Xt 2 t最后 AR（ 5）模型的残差分析为（即模型的残差白噪声检验）LB（6）=3.45 其 p 值为 0.4851 大于 0.05 ，故通过检验（其他的也是类似） ，所以该模型的拟合效果很好。（5 ）用此模型做预测data ex; input x;t=intnx('year', '01jan1

9、975'd,_n_- 1);format tyear4. ;cards ;8.7 2.47.6 11.7 7.6 7.85.2 9.110.9 15.2 13.58.8 11.6 11.3 4.1 3.89.2 14.2 14 13.1 10.9 10 9.3 7.8 7.6 8.4 8.3 9.1 10 10.1 10.4 10.7proc arima ; identify var =x nlag =12 minic p=( 0 5) q=( 0 5); estimate p=( 1, 2) method =cls;forecast lead =5 id =t out =result

10、s;proc gplot data =results;plot x*t= 1 forecast*t= 2 l95*t= symbol1 c =blue i =jion v =star; symbol2 c=red i =jion v =none l symbol3 c =green i =jion v =none=1 w=1;l =2 w=2;3/ overlay ;run ;可以看出， 原始数据绝大部分都在预测区域内， 而且越是近期的数据离预测曲线越近， 表明 模型建立的比较合理，预测效果比较精准。17. 线性拟合：序列呈现出显著的线性特征我国 1995-2006 年期刊种类数据如表（数据略

11、） ，试建模并进行预测。（1）画出改时间序列的散点图data ex; input x;t=_n_; cards ;7583 7916 7918 7999 8187 8725 8889 9029 9074 9490 9468 9468; proc gplot ; symbol i =jiontv =dot; plot x*t; run ;发现该序列有显著的线性递增趋势，于是考虑使用线性模型：X t a bt t tt2E( t) 0,Var ( t)22）拟合模型，进行参数估计和检验。data ex; input x;t=_n_; cards7583 7916 7918 7999 8187 87

12、25 8889 9029 9074 9490 9468 9468 ; proc reg ; model x=t; run ;在 t 检验中， p 值都小于 0.05 拒绝 H0 接受 H1 ，认为参数都显著非 0。 在 F 检验中， p 值小于 0.05 认为方程显著。得到趋势模型： Tt 7423.40909 188.01399t7423.40909 - 188.01399 *t;(3)对残差 t Xt Tt 进行白噪声检验。 data ex; input x;t=_n_;epsilon=x- cards ;7583 7916 7918 7999 8187 8725 8889 9029 90

13、74 9490 9468 9468 ; proc arima ; identify var =epsilon nlag =12 ; runLB（6）=5.88 其 p 值为 0.4367 ，远远大于 0.05 ，故通过参差检验，残差为白噪声序列。 表明模型信息提取充分。 我们将利用前面得到的趋势模型 Tt 7423.40909 188.01399t 进 行预测，将 t=13 t=14 分别代入模型中，可得：2007 年我国期刊种类预测值： 7423.40909+188.01399*13=9867.592008 年我国期刊种类预测值： 7423.40909+188.01399*14=10055.

14、6018. 曲线拟合：序列在发展过程中呈现出以不同的速率增长或下降，或者由逐渐增长到逐渐 衰退等各种不同的性态，即非线性特征对 1980-2007 年 GDP 的数据进行建模和预测（数据略）（1）画出 1980-2007 年全国 GDP 散点图：data ex; input GDP;t=intnx（ 'year' , '01jan1980'd ,_n_- 1）;format t year4. ;cards ;4517.8 4862.4 5294.7 5934.5 7171 8964.4 10202.211962.5 14928.3 16909.2 18547.9

15、 21617.8 26638.1 34634.446759.4 58478.1 67884.6 74462.6 78345.2 82067.5 89442.297314.8 118020.7 135822.8 159878.3 183217.4 211923 249530 proc gplot ; plot GDP*t; symbol c =blue i =jion v =dot; run ;由时序图可以看出序列呈现指数增长趋势，因此，我们提出初步模型Xt a bect t（2）用最小二乘法进行参数估计：data ex; input x_t;t=_n_;cards ;4517.8 4862.4

16、 5294.7 5934.5 7171 8964.4 10202.211962.5 14928.3 16909.2 18547.9 21617.8 26638.1 34634.446759.4 58478.1 67884.6 74462.6 78345.2 82067.5 89442.297314.8 118020.7 135822.8 159878.3 183217.4 211923 249530proc nlin method =gauss; model x_t=a+b*exp(c*t);parameters a= 1 b= 0.5 c= 0.5 ;der.a= 1 ;der.b=exp(

17、c*t);der.c=t*b*exp(c*t);output predicted =x_that out =result;proc gplot data =result; plot x_t*t= 1 x_that*t= 2/ overlay symbol1 c =black i =none v=dot l =1 w=1;symbol2 c =blue i =jion v=dot l =2 w=2;Run;得到非线性模型的拟合图：非线性模型的参数估计通过了检验。(3) 对残差进行白噪声检验：t Xt 3989.1 7252.7e0.1305tdata ex; input x_t;t=_n_; e

18、psilon=x_t+ 3989.1 - 7252.7 *exp( 0.1305 *t); cards ;4517.8 4862.4 5294.7 5934.5 7171 8964.4 10202.211962.5 14928.3 16909.2 18547.9 21617.8 26638.1 34634.446759.4 58478.1 67884.6 74462.6 78345.2 82067.5 89442.297314.8 118020.7 135822.8 159878.3 183217.4 211923 249530procarima ; identifyvar =epsilonn

19、lag =12;run ;可以看出， LB(6)=54.08 其 p 值小于 0.05 ，没有通过检验。 LB(12)=60.12 其 p 值小 于 0.05 ，没有通过检验。表明模型参差不是白噪声序列，信息提取不充分，我们还需要对模型进行进一步改进。(4 )由( 1)中的散点图可以看出，序列具有一定的指数趋势，首先对其取对数： data ex; input GDP;y=log(GDP);t=intnx( 'year' , '01jan1980'd ,_n_- 1);format t year4. ;cards ;4517.8 4862.4 5294.7 593

20、4.5 7171 8964.4 10202.211962.5 14928.3 16909.2 18547.9 21617.8 26638.1 34634.446759.4 58478.1 67884.6 74462.6 78345.2 82067.5 89442.297314.8 118020.7 135822.8 159878.3 183217.4 211923 249530proc gplot ; plot y*t; symbol c =blue i =jion v =dot;run ;取对数后的序列呈现线性增长趋势，对取对数后的数据进行一阶差分： data ex; input GDP;y

21、=log(GDP);y1=dif(y);t=intnx( 'year' , '01jan1980'd ,_n_- 1);format t year4. ;cards ;4517.8 4862.4 5294.7 5934.5 7171 8964.4 10202.211962.5 14928.3 16909.2 18547.9 21617.8 26638.1 34634.4proc gplot ; plot y1*t; symbol c =blue i =jion v =dot; run ;取对数后为平稳序列，接下来按照 16 中平稳序列的处理方法进行建模和预测。面

22、是临海大药房季节药博利康尼19. 季节指数水平法（适用于无明显上升或下降变动趋势，主要受季节变动和不规则变动影 响的 时间序列）季度20042005200620071901101039826056665832535323641201551231152004-2007 年各季度的销售量。1）先做出时序图：data ex; input x;t=_n_;formatt year4. ;cards ;90 60 25 120 110 56 35 155 103 66 32 123 98 58 36 115 proc gplot ; plot x*t; symbol c =blue i =jion v

23、=dot; run ;由上图可以看出， 临海大药房季节药博利康尼为季节性药品， 在不同的季度其销售量的K 期中心移动平均：先计算：Xt 12k为偶数：在计算：X1k1 (X k Xt k 1 L Xt K Xt k 1)221X t 1 (X k K Xt K X k X k) 2 k t k 1 t t k 1 t k 2 2 2波动很明显，出现明显的销售旺季和销售淡季，改时间序列为非平稳的。（2）得到药品销售各平均值和季节指数如下：总平均 =80.125第一季度平均 =100.25第二季度平均 =60.00第三季度平均 =32.00第四季度平均 =128.25季节指数 S1=100.25/

24、80.125=1.25117 S2=0.748830 S3=0.399376 S4=1.600624 可以看出季节指数差别很大，说明药品销售量有明显的季节性。（3）对 2008 年进行预测： 进行预测事，一般选择最近年份的平均值，这里， 2007 年的每季度的平均值为 （98+58+36+115 ） /4=76.75第一季度： 76.75*1.25117=95.9375第二季度： 76.75*0.74883=57.5625第三季度： 76.75*0.399376=30.7第四季度： 76.75*1.600624=122.8 全年： 95.9375+57.5625+30.7+122.8=3072

25、0. 季节指数趋势法：如果序列既存在季节变动，同时各年水平或同月（同季）水平呈现上 升或下降趋势，这是应该采用季节指数趋势法。K 期的移动平均： Xt Xt k 1 X t k 2 K Xt 1 X tk1Xt k 1 1Xt k 1)22k为奇数：X t 1(X k 1 Xk t 21最后计算： Xt(Xt 12 Xt 21 )2 2 2例：下表为某市某品牌桶装水销售记录数据，我们据此数据做出2006 年该桶装水的销售记录的预测。年份季度销售量 Xt/ 万桶4 期中心化移 动平均季节比率消除季节影响的销售 量 Xt/St200114659.4595126358.1637738862.51.4

26、0864.6961245163.8750.7984364.38623200215065.1250.7677564.629927066.251.056664.6264139168.251.333366.9016745770.50.8085171.96108200316072.50.8275977.5558827874.251.0505172.0122939974.6251.3266372.7831346375.6250.8330679.53593200415778.3750.7272773.6780928979.3751.1212682.16787311080.1251.3728580.8701

27、546082.3750.7283875.74851200516685.50.7719385.3114729888.751.1042390.47698312692.6330847088.37326（1）先做出序列的散点图： data ex; input x;t=intnx('quarter', '01jan2001'd,_n_- 1); format tyyq4. ;cards46 6388 51 5070 91 5760 7899 63 57 89110 60 66 98 126 70procgplot ; plotx*t; symbolc =bluei =ji

28、on v =dot;run ;观察该市桶装水销售量序列由长期趋势和一年为周期的季节波动共同影响， 因此要考虑剔除 趋势再分析季节特征。（2）剔除趋势：NO1 ：以一年的季度数 4为 k，对时间序列观测值进行 k项中心移动平均，形成新的序列， 如表中第四列所示，以此为长期趋势 Tt 。在消除了时间序列的季节和不规则性之后，可以 看出，具有很明显的趋势性。data ex; input x;t=intnx('quarter' , '01jan2001'd,_n_- 1); format tyyq4. ;cards ;62.563.875 65.125 66.2568.

29、2570.572.574.2574.62575.62578.375 79.375 80.12582.37585.588.75proc gplot ; plot x*t; symbol c =blue i =jionv =dot;run ;NO2 ：将各期观察值除以同期移动均值作为季节比率，即表中倒数第二列。NO3 ：各年同季的季节比率平均，季节平均比率可以消除不规则变动，即得其季节指数St。做法与 19 中（ 2）中的做法是一样的。得到：季度季节不规则成分的数值（ StIt ）季节指数（ St）10.76775 0.82759 0.72727 0.771930.7736421.0566 1.0

30、5051 1.12126 1.104231.0831531.408 1.3333 1.32663 1.372851.3602140.79843 0.80851 0.83306 0.728380.79209季节指数模型需要平均季节指数等于1.00，在表中， 季节指数总和为 4.009084，平均值基本等于 1，所以不用调整。对于不等于 1 的，要进行一些调整：用每一个季节指数乘以季度总 和再除调整之前的季节指数之和。NO4 ：用时间序列的每个观察值处以相应的季节指数， 消除时间序列季节影响， 得到趋势序 X列，即 t Tt ，表中最后一列。St（3）利用消除季节影响后的线性趋势预测值建立回归模型

31、Tt a bt ，进行回归分析与白噪声检验。data ex; input T_t;t=_n_; format t yyq4. ;cards ;59.45951 58.16377 64.69612 64.38623 64.6299 64.6264166.9016771.9610877.5558872.0122972.7831379.5359373.6780982.1678780.8701575.7485185.3114790.4769892.63308 88.37326; proc reg ; modelT_t=t;run ;趋势模型的参数和模型显著，都通过了检验。其表达式为：Tt 57.050

32、63 1.64266t进行白噪声检验：data ex; inputT_t;t=_n_; format t yyq4. ;res=T_t- 57.05063 - 1.64266 *t;cards ;59.45951 58.16377 64.69612 64.38623 64.6299 64.6264166.9016771.9610877.5558872.0122972.7831379.5359373.6780982.1678780.8701575.7485185.3114790.4769892.63308 88.37326proc arima ; identify var =res nlag =

33、12; run ;可以看出，白噪声检验通过。趋势模型是适应的。（ 4 ）进行预测和季节调整。对同时有趋势和季节成分的时间序列，进行预测的最后是用季节指数调整趋势预测值。2006.1：Tt=57.05063+1.64266*17=91.54649调整后：91.54649*0.77364=70.82362006.2： Tt=57.05063+1.64266*18=93.18915调整后：93.18915*1.08315=100.9同理 2006.3 2006.4。21.X 11方法（变量的第二类调查统计方法季节调整方案）：比季节指数法预测更精确例：某市医院三年中各个季度接受肿瘤治疗的人次如表所示，

34、是进行季节调整，计算出季 节指数，并进行预报。季度1234199914245143351442614518200014610147011479214884200114976150661506615249(1)data ex; input x;t=intnx( 'quarter' , '01jan1999'd ,_n_- 1); format t yyq.4 ;cards ;14245 14335 14426 1451814610 14701 14792 1488414976 15066 15066 15249proc x11 data= ex; quarterly date=t;output out=out b1=x d10=season d11=adjusted d12=trend d13=irr;/*b1 为原序列 d10 列出季节指数 d11 列出季节调整后的序列 d12 趋势拟合 d13 列出最后 的不规则波动值 */run ;图 1 生成的原始数据表图2 时间序列的季节因子 St（ %）图 3 季节