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服從正態分布相吻合，所以可以認為該模型擬合成功。最