時(shí)間序列整合分析(一)

1、1. 缺省值的補(bǔ)足：2. 時(shí)序圖：（檢驗(yàn)平穩(wěn)性）3. 自相關(guān)函數(shù)：（檢驗(yàn)平穩(wěn)性）4. 計(jì)算標(biāo)準(zhǔn)正態(tài)分布的概率：5. 計(jì)算標(biāo)準(zhǔn)正態(tài)分布的分位數(shù)：6. 計(jì)算標(biāo)準(zhǔn)t分布的概率7. 計(jì)算標(biāo)準(zhǔn)t分布的分位數(shù)8. 計(jì)算標(biāo)準(zhǔn)F分布的概率9. 計(jì)算標(biāo)準(zhǔn)F分布的分位數(shù)10. 計(jì)算標(biāo)準(zhǔn)卡方分布的概率11. 計(jì)算標(biāo)準(zhǔn)卡方分布的分位數(shù)12. 方差的同齊性檢驗(yàn)：將數(shù)據(jù)進(jìn)行適當(dāng)分組，這里將4個(gè)分為一組，一共四組PrF的值大于0.05 故接受H0，認(rèn)為各組方差之間沒(méi)有顯著的差異。13. 方差的同質(zhì)性檢驗(yàn)：將數(shù)據(jù)進(jìn)行適當(dāng)分組，這里將4個(gè)分為一組，一共四組 根據(jù)上結(jié)果列出方差分析表： 方差來(lái)源平方和自由度均方和F值顯著性A誤差

2、878376140372932029279270186.44.17總和228210523 F的p值小于0.05 我們認(rèn)為原始數(shù)據(jù)方差不同質(zhì)。14. 序列的白噪聲檢驗(yàn)（檢驗(yàn)純隨機(jī)性）： 可以看出，LB(6)=95.84，其p值小于0.05； LB(12)=190.40，其p值小于0.05；顯然該序列不是白噪聲序列，即不是純隨機(jī)性序列。（p值都大于0.05時(shí)才是純隨機(jī)序列）15. 平穩(wěn)序列的自相關(guān)函數(shù)和偏自相關(guān)函數(shù)的形式：（沒(méi)有程序的） 模型AR(p)MA(q)ARMA(p,q)ACF自相關(guān)拖尾截尾拖尾PACF偏自相關(guān)截尾拖尾拖尾16. 一個(gè)例子：（利用平穩(wěn)序列建模 進(jìn)行預(yù)測(cè)） 我國(guó)1975-20

3、06年GDP的年增長(zhǎng)率為下表（數(shù)據(jù)略），對(duì)我國(guó)1975-2006年GDP的年增長(zhǎng)率進(jìn)行建模，并對(duì)2007至2011年我國(guó)的GDP增長(zhǎng)率進(jìn)行預(yù)測(cè)。 (1)首先畫(huà)出我國(guó)1975-2006年GDP增長(zhǎng)率的時(shí)序圖。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;從圖中直觀的可以看出有奇異點(diǎn)（2） 將奇異點(diǎn)看成缺省值，利用以下程序來(lái)求缺省點(diǎn)的值：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;結(jié)果可知，缺省點(diǎn)的值為2.4（3） 利用修正后的數(shù)據(jù)再進(jìn)行時(shí)序分析，根據(jù)以下程序：可以看出GDP增長(zhǎng)率修正后的數(shù)據(jù)序列平穩(wěn)。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;從上圖中可以看出，有些參數(shù)不顯著，我們將其去掉，建立最精干的模型。（其中可以看出，AR1,3 AR1,4 AR1,5 的p值遠(yuǎn)遠(yuǎn)大于0.05）所以，將estimate p=5改為estimate p=(1,2)，即程序?yàn)椋篸ata 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;可見(jiàn)所有的p 值都小于0.05 通過(guò)了檢驗(yàn)。所以模型為：最后AR（5）模型的殘差分析為（即模型的殘差白噪聲檢驗(yàn)）：LB(6)=3.45 其p值為0.4851大于0.0

8、5 ，故通過(guò)檢驗(yàn)（其他的也是類(lèi)似），所以該模型的擬合效果很好。（5） 用此模型做預(yù)測(cè)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;可以看出，原始數(shù)據(jù)絕大部分都在預(yù)測(cè)區(qū)域內(nèi)，而且越是近期的數(shù)據(jù)離預(yù)測(cè)曲線(xiàn)越近，表明模型建立的比較合理，預(yù)測(cè)效果比較精準(zhǔn)。17. 線(xiàn)

10、性擬合：序列呈現(xiàn)出顯著的線(xiàn)性特征我國(guó)1995-2006年期刊種類(lèi)數(shù)據(jù)如表（數(shù)據(jù)略），試建模并進(jìn)行預(yù)測(cè)。（1） 畫(huà)出改時(shí)間序列的散點(diǎn)圖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;發(fā)現(xiàn)該序列有顯著的線(xiàn)性遞增趨勢(shì)，于是考慮使用線(xiàn)性模型：來(lái)擬合該趨勢(shì)的發(fā)展。（2） 擬合模型，進(jìn)行參數(shù)估計(jì)和檢驗(yàn)。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檢驗(yàn)中，p值都小于0.05 拒絕H0 接受H1，認(rèn)為參數(shù)都顯著非0。在F檢驗(yàn)中，p值小于0.05 認(rèn)為方程顯著。得到趨勢(shì)模型：（3） 對(duì)殘差進(jìn)行白噪聲檢驗(yàn)。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，遠(yuǎn)遠(yuǎn)大于0.05，故通過(guò)參差檢驗(yàn)，殘差為白噪聲序列。表明模型信息提取充分。我們將利用前面得到的趨勢(shì)模型進(jìn)行預(yù)測(cè)，將t=13 t=14分別代入模型中，可得：2007年我國(guó)期刊種類(lèi)預(yù)測(cè)值：7423.40909+188.01399*13=9867.592008年我國(guó)期刊種類(lèi)預(yù)測(cè)值：7423.40909+188.01399*14=10055.6018. 曲線(xiàn)擬合：序列在發(fā)展過(guò)程中呈現(xiàn)出以不同的速率增長(zhǎng)或下降，或者由逐漸增長(zhǎng)到逐漸衰退等各種不同的性態(tài)，即非線(xiàn)性特征對(duì)1980-2007年GDP的數(shù)據(jù)進(jìn)行建模和預(yù)測(cè)（數(shù)據(jù)略）（1）畫(huà)出1980

13、-2007年全國(guó)GDP散點(diǎn)圖：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;由時(shí)序圖可以看出序列呈現(xiàn)指數(shù)增長(zhǎng)趨勢(shì)，因此，我們提出初步模型（2） 用最小二乘法進(jìn)行參數(shù)估計(jì)：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;得到非線(xiàn)性模型的擬合圖：非線(xiàn)性模型的參數(shù)估計(jì)通過(guò)了檢驗(yàn)。（3） 對(duì)殘差進(jìn)行白噪聲檢驗(yàn)：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，沒(méi)有通過(guò)檢驗(yàn)。 LB(12)=60.12其p值小于0.05，沒(méi)有通過(guò)檢驗(yàn)。表明模型參差不是白噪聲序列，信息提取不充分，我們還需要對(duì)模型進(jìn)行進(jìn)一步改進(jìn)。（4） 由（1）中的散點(diǎn)圖可以看出，序列具有一定的指數(shù)趨勢(shì)，首先對(duì)其取對(duì)數(shù)：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;取對(duì)數(shù)后的序列呈現(xiàn)線(xiàn)性增長(zhǎng)

19、趨勢(shì)，對(duì)取對(duì)數(shù)后的數(shù)據(jù)進(jìn)行一階差分：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;取對(duì)數(shù)后為平穩(wěn)序列，接下來(lái)按照16中平穩(wěn)序列的處理方法進(jìn)行建模和預(yù)測(cè)。19. 季節(jié)指數(shù)水平法（適用于無(wú)明顯上升或下降變動(dòng)趨勢(shì)，主要受季節(jié)變動(dòng)和不規(guī)則變動(dòng)影響的 時(shí)間序列）下面是臨海大藥房季節(jié)藥博利康尼2004-2007年各季度的銷(xiāo)售量。 季度2004200520062007190110103982605666583253532364120155123115（1） 先做出時(shí)序圖：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;由上圖可以看出，臨海大藥房季節(jié)藥博利康尼為季節(jié)性藥品，在不同的季度其銷(xiāo)售量的波動(dòng)很明顯，出現(xiàn)明顯的銷(xiāo)售旺季和銷(xiāo)售淡季，改時(shí)間序列為非平穩(wěn)的。（2） 得到藥品銷(xiāo)售各平均值和季節(jié)指數(shù)如下： 總平均=80.125 第一季度平均=100.25 第二季度平均=60.00 第三季度平均=32.00 第四季度平均=128.25季節(jié)指數(shù)S1=100.25/80.125=1.2511

22、7 S2=0.748830 S3=0.399376 S4=1.600624可以看出季節(jié)指數(shù)差別很大，說(shuō)明藥品銷(xiāo)售量有明顯的季節(jié)性。（3） 對(duì)2008年進(jìn)行預(yù)測(cè)：進(jìn)行預(yù)測(cè)事，一般選擇最近年份的平均值，這里，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. 季節(jié)指數(shù)趨勢(shì)法：如果序列既存在季節(jié)變動(dòng)，同時(shí)各年

23、水平或同月（同季）水平呈現(xiàn)上升或下降趨勢(shì)，這是應(yīng)該采用季節(jié)指數(shù)趨勢(shì)法。K期的移動(dòng)平均：K期中心移動(dòng)平均： 最后計(jì)算：例：下表為某市某品牌桶裝水銷(xiāo)售記錄數(shù)據(jù)，我們據(jù)此數(shù)據(jù)做出2006年該桶裝水的銷(xiāo)售記錄的預(yù)測(cè)。年份季度銷(xiāo)售量Xt/萬(wàn)桶4期中心化移動(dòng)平均季節(jié)比率消除季節(jié)影響的銷(xiāo)售量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） 先做出序列的散點(diǎn)圖：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;觀察該市桶裝水銷(xiāo)售量序列由長(zhǎng)期趨勢(shì)和一年為周期的季節(jié)波動(dòng)共同影響，因此要考慮剔除趨勢(shì)再分析季節(jié)特征。（2） 剔除趨勢(shì)：NO1：以一年的季度數(shù)4為k，對(duì)時(shí)間序列觀測(cè)值進(jìn)行k項(xiàng)中心移動(dòng)平均，形成新的序列，如表中第四列所示，以此為長(zhǎng)期趨勢(shì)Tt。在消除了時(shí)間序列的季節(jié)和不規(guī)則性之后，可以看出，具

26、有很明顯的趨勢(shì)性。 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：將各期觀察值除以同期移動(dòng)均值作為季節(jié)比率，即表中倒數(shù)第二列。NO3：各年同季的季節(jié)比率平均，季節(jié)平均比率可以消除不規(guī)則變動(dòng)，即得其季節(jié)指數(shù)St。

27、做法與19中（2）中的做法是一樣的。 得到：季度季節(jié)不規(guī)則成分的數(shù)值（StIt）季節(jié)指數(shù)（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季節(jié)指數(shù)模型需要平均季節(jié)指數(shù)等于1.00，在表中，季節(jié)指數(shù)總和為4.009084，平均值基本等于1，所以不用調(diào)整。對(duì)于不等于1的，要進(jìn)行一些調(diào)整：用每一個(gè)季節(jié)指數(shù)乘以季度總和再除調(diào)整之前的季節(jié)指數(shù)之和

28、。NO4：用時(shí)間序列的每個(gè)觀察值處以相應(yīng)的季節(jié)指數(shù)，消除時(shí)間序列季節(jié)影響，得到趨勢(shì)序列，即，表中最后一列。（3） 利用消除季節(jié)影響后的線(xiàn)性趨勢(shì)預(yù)測(cè)值建立回歸模型，進(jìn)行回歸分析與白噪聲檢驗(yàn)。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;趨勢(shì)模型的參數(shù)和模型顯著，都通過(guò)了檢驗(yàn)。其表達(dá)式為：進(jìn)行白噪聲檢驗(yàn)：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;可以看出，白噪聲檢驗(yàn)通過(guò)。趨勢(shì)模型是適應(yīng)的。（4） 進(jìn)行預(yù)測(cè)和季節(jié)調(diào)整。對(duì)同時(shí)有趨勢(shì)和季節(jié)成分的時(shí)間序列，進(jìn)行預(yù)測(cè)的最后是用季節(jié)指數(shù)調(diào)整趨勢(shì)預(yù)測(cè)值。2006.1：Tt=57.05063+1.64266*17=91.54649 調(diào)整后：91.54649*0.77364=70.82362006.2：Tt=57.05063+1.64266*18=93.18915 調(diào)整后：93.18915*1.08315=100.9同理2006.3 2006.4。21. X11方法（變量的第二類(lèi)調(diào)查統(tǒng)計(jì)方法季節(jié)調(diào)整方案）：比季節(jié)指數(shù)法預(yù)測(cè)更精確例：某市醫(yī)院三

31、年中各個(gè)季度接受腫瘤治療的人次如表所示，是進(jìn)行季節(jié)調(diào)整，計(jì)算出季節(jié)指數(shù)，并進(jìn)行預(yù)報(bào)。季度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列出季節(jié)指數(shù)d11列出季節(jié)調(diào)整后的序列 d12趨勢(shì)擬合 d13列出最后的不規(guī)則波動(dòng)值*/run;圖1 生成的原始數(shù)據(jù)表圖2 時(shí)間序列的季節(jié)因子St（%）圖3 季節(jié)調(diào)整后的序