时间序列整合分析(一)

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

2、878376140372932029279270186.44.17总和228210523 F的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偏自相关截尾拖尾拖尾16. 一个例子：（利用平稳序列建模 进行预测） 我国1975-20

3、06年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.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.7;proc gplot;symbol i=jiont v=dot

4、;plot x*t;run;从图中直观的可以看出有奇异点（2） 将奇异点看成缺省值，利用以下程序来求缺省点的值：data ex;input x;time=intnx(month,01jan1975d,_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.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.7;proc expand data=ex out=ex1;id time

5、;proc print data=ex1;run;结果可知，缺省点的值为2.4（3） 利用修正后的数据再进行时序分析，根据以下程序：可以看出GDP增长率修正后的数据序列平稳。BIC（5,0）=-0.24488的值最小，考虑建立AR(5)模型。（4） 模型的建立data ex;input x;time=intnx(month,01jan1975d,_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.8 9.2 14.2 14 13.1 10.9 10

6、 9.3 7.8 7.6 8.4 8.3 9.1 10 10.1 10.4 10.7;proc arima;identify var=x nlag=12 minic p=(0:5) q=(0:5);estimate p=5;run;从上图中可以看出，有些参数不显著，我们将其去掉，建立最精干的模型。（其中可以看出，AR1,3 AR1,4 AR1,5 的p值远远大于0.05）所以，将estimate p=5改为estimate p=(1,2)，即程序为：data ex;input x;time=intnx(month,01jan1975d,_n_-1);format time year4.;car

7、ds;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.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.7;proc arima;identify var=x nlag=12 minic p=(0:5) q=(0:5);estimate p=(1,2);run;可见所有的p 值都小于0.05 通过了检验。所以模型为：最后AR（5）模型的残差分析为（即模型的残差白噪声检验）：LB(6)=3.45 其p值为0.4851大于0.0

8、5 ，故通过检验（其他的也是类似），所以该模型的拟合效果很好。（5） 用此模型做预测data ex;input x;t=intnx(year,01jan1975d,_n_-1);format t 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.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.7;proc arima;identify var=x nlag=12 minic p=(0 5) q=(0

9、 5);estimate p=(1,2) method=cls;forecast lead=5 id=t out=results;proc gplot data=results;plot x*t=1 forecast*t=2 l95*t=3/overlay;symbol1 c=blue i=jion v=star;symbol2 c=red i=jion v=none l=1 w=1;symbol3 c=green i=jion v=none l=2 w=2;run;可以看出，原始数据绝大部分都在预测区域内，而且越是近期的数据离预测曲线越近，表明模型建立的比较合理，预测效果比较精准。17. 线

10、性拟合：序列呈现出显著的线性特征我国1995-2006年期刊种类数据如表（数据略），试建模并进行预测。（1） 画出改时间序列的散点图data ex;input x;t=_n_;cards;7583 7916 7918 7999 8187 8725 8889 9029 9074 9490 9468 9468;proc gplot;symbol i=jiont v=dot;plot x*t;run;发现该序列有显著的线性递增趋势，于是考虑使用线性模型：来拟合该趋势的发展。（2） 拟合模型，进行参数估计和检验。data ex;input x;t=_n_;cards;7583 7916 7918 79

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

12、lag=12;run;LB(6)=5.88 其p值为0.4367，远远大于0.05，故通过参差检验，残差为白噪声序列。表明模型信息提取充分。我们将利用前面得到的趋势模型进行预测，将t=13 t=14分别代入模型中，可得：2007年我国期刊种类预测值：7423.40909+188.01399*13=9867.592008年我国期刊种类预测值：7423.40909+188.01399*14=10055.6018. 曲线拟合：序列在发展过程中呈现出以不同的速率增长或下降，或者由逐渐增长到逐渐衰退等各种不同的性态，即非线性特征对1980-2007年GDP的数据进行建模和预测（数据略）（1）画出1980

13、-2007年全国GDP散点图：data ex;input GDP;t=intnx(year,01jan1980d,_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.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

14、;proc gplot;plot GDP*t;symbol c=blue i=jion v=dot;run;由时序图可以看出序列呈现指数增长趋势，因此，我们提出初步模型（2） 用最小二乘法进行参数估计：data ex;input x_t;t=_n_;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 1180

15、20.7 135822.8 159878.3 183217.4 211923 249530;proc 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(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=bl

16、ue i=jion v=dot l=2 w=2;Run;得到非线性模型的拟合图：非线性模型的参数估计通过了检验。（3） 对残差进行白噪声检验：data ex;input x_t;t=_n_;epsilon=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.297

17、314.8 118020.7 135822.8 159878.3 183217.4 211923 249530;proc arima;identify var=epsilon nlag=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,01jan1980d,_n_-1);f

18、ormat 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.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 y*t;symbol c=blue i=jion v=dot;run;取对数后的序列呈现线性增长

19、趋势，对取对数后的数据进行一阶差分：data ex;input GDP;y=log(GDP);y1=dif(y);t=intnx(year,01jan1980d,_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.446759.4 58478.1 67884.6 74462.6 78345.2 82067.5 89442.297314.8 118020.7 135822.8 159878.

20、3 183217.4 211923 249530;proc gplot;plot y1*t;symbol c=blue i=jion v=dot;run;取对数后为平稳序列，接下来按照16中平稳序列的处理方法进行建模和预测。19. 季节指数水平法（适用于无明显上升或下降变动趋势，主要受季节变动和不规则变动影响的 时间序列）下面是临海大药房季节药博利康尼2004-2007年各季度的销售量。 季度2004200520062007190110103982605666583253532364120155123115（1） 先做出时序图：data ex;input x;t=_n_;format t ye

21、ar4.;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=dot;run;由上图可以看出，临海大药房季节药博利康尼为季节性药品，在不同的季度其销售量的波动很明显，出现明显的销售旺季和销售淡季，改时间序列为非平稳的。（2） 得到药品销售各平均值和季节指数如下： 总平均=80.125 第一季度平均=100.25 第二季度平均=60.00 第三季度平均=32.00 第四季度平均=128.25季节指数S1=100.25/80.125=1.2511

22、7 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=30720. 季节指数趋势法：如果序列既存在季节变动，同时各年

23、水平或同月（同季）水平呈现上升或下降趋势，这是应该采用季节指数趋势法。K期的移动平均：K期中心移动平均： 最后计算：例：下表为某市某品牌桶装水销售记录数据，我们据此数据做出2006年该桶装水的销售记录的预测。年份季度销售量Xt/万桶4期中心化移动平均季节比率消除季节影响的销售量Xt/St200114659.4595126358.1637738862.51.40864.6961245163.8750.7984364.38623200215065.1250.7677564.629927066.251.056664.6264139168.251.333366.9016745770.50.808517

24、1.96108200316072.50.8275977.5558827874.251.0505172.0122939974.6251.3266372.7831346375.6250.8330679.53593200415778.3750.7272773.6780928979.3751.1212682.16787311080.1251.3728580.8701546082.3750.7283875.74851200516685.50.7719385.3114729888.751.1042390.47698312692.6330847088.37326（1） 先做出序列的散点图：data ex;i

25、nput x;t=intnx(quarter,01jan2001d,_n_-1);format t yyq4.;cards;466388515070915760789963578911060669812670;proc gplot;plot x*t;symbol c=blue i=jion v=dot;run;观察该市桶装水销售量序列由长期趋势和一年为周期的季节波动共同影响，因此要考虑剔除趋势再分析季节特征。（2） 剔除趋势：NO1：以一年的季度数4为k，对时间序列观测值进行k项中心移动平均，形成新的序列，如表中第四列所示，以此为长期趋势Tt。在消除了时间序列的季节和不规则性之后，可以看出，具

26、有很明显的趋势性。 data ex;input x; t=intnx(quarter,01jan2001d,_n_-1);format t yyq4.; cards; 62.563.87565.12566.2568.2570.572.574.2574.625 75.62578.37579.37580.12582.37585.588.75 ; proc gplot;plot x*t;symbol c=blue i=jion v=dot;run;NO2：将各期观察值除以同期移动均值作为季节比率，即表中倒数第二列。NO3：各年同季的季节比率平均，季节平均比率可以消除不规则变动，即得其季节指数St。

27、做法与19中（2）中的做法是一样的。 得到：季度季节不规则成分的数值（StIt）季节指数（St）10.76775 0.82759 0.72727 0.771930.7736421.0566 1.05051 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的，要进行一些调整：用每一个季节指数乘以季度总和再除调整之前的季节指数之和

28、。NO4：用时间序列的每个观察值处以相应的季节指数，消除时间序列季节影响，得到趋势序列，即，表中最后一列。（3） 利用消除季节影响后的线性趋势预测值建立回归模型，进行回归分析与白噪声检验。data ex;input T_t;t=_n_;format t yyq4.;cards;59.4595158.1637764.6961264.3862364.629964.6264166.9016771.9610877.5558872.0122972.7831379.5359373.6780982.1678780.8701575.7485185.3114790.4769892.6330888.37326;p

29、roc reg;model T_t=t;run;趋势模型的参数和模型显著，都通过了检验。其表达式为：进行白噪声检验：data ex;input T_t;t=_n_;format t yyq4.;res=T_t-57.05063-1.64266*t;cards;59.4595158.1637764.6961264.3862364.629964.6264166.9016771.9610877.5558872.0122972.7831379.5359373.6780982.1678780.8701575.7485185.3114790.4769892.6330888.37326;proc arima

30、;identify var=res nlag=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. X11方法（变量的第二类调查统计方法季节调整方案）：比季节指数法预测更精确例：某市医院三

31、年中各个季度接受肿瘤治疗的人次如表所示，是进行季节调整，计算出季节指数，并进行预报。季度1234199920002001142451461014976143351470115066144261479215066145181488415249（1）data ex;input x;t=intnx(quarter,01jan1999d,_n_-1);format t yyq.4;cards;142451433514426145181461014701147921488414976150661506615249;proc 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 季节调整后的序