第一章習(xí)題答案 略

　第二章習(xí)題答案 2.1

　 R 命令 x<-1:20 x<-ts(x) plot(x,type="o") acf(x)$acf

　答案：

�。�1）非平穩(wěn)，有典型線性趨勢

�。�2）延遲 1-6 階自相關(guān)系數(shù)如下：

　（3）典型的具有單調(diào)趨勢的時間序列樣本自相關(guān)圖

　2.2 R 命令 #先讀入數(shù)據(jù)文件 co2<-ts(E2_2$co2,start=c(1975,1),frequency = 12) plot(co2,type="o") acf(co2)$acf

�。�1）非平穩(wěn)，時序圖如下

�。�2）1-24 階自相關(guān)系數(shù)如下 [2,]

　0.90750778 [10,]

　0.36433219 [18,]

　0.00135460 [3,]

　0.72171377 [11,]

　0.48471672 [19,] -0.03247916 [4,]

　0.51251814 [12,]

　0.58456166 [20,] -0.02709893 [5,]

　0.34982244 [13,]

　0.60197891 [21,]

　0.01123597 [6,]

　0.24689637 [14,]

　0.51841257 [22,]

　0.08274806 [7,]

　0.20309427 [15,]

　0.36856286 [23,]

　0.17010715 [8,]

　0.21020799 [16,]

　0.20671211 [24,]

　0.24319854 [9,]

　0.26428810 [17,]

　0.08138070 [25,]

　0.25252294

�。�3）自相關(guān)圖呈現(xiàn)典型的長期趨勢與周期并存的特征

　 2.3

　 R R 命令

　#先讀入數(shù)據(jù)文件 rain<-ts(E2_3$rain,start=c(1945,1),frequency = 12) plot(rain,type="o") acf(rain,lag.max =24)$acf for(i in 1:6) print(Box.test(rain,lag=3*i))

　 答案

�。�1）1-24 階自相關(guān)系數(shù) [2,]

　0.012770216 [10,]

　0.013882228 [18,] -0.245494618 [3,]

　0.041600613 [11,]

　0.109450351 [19,]

　0.066461869 [4,] -0.043230426 [12,]

　0.217295088 [20,] -0.139454035 [5,] -0.178692841 [13,]

　0.315872697 [21,] -0.034028373 [6,] -0.251298873 [14,] -0.025053744 [22,]

　0.205723132 [7,] -0.093810241 [15,]

　0.075320665 [23,] -0.009866008 [8,] -0.067777725 [16,] -0.141206897 [24,]

　0.080311638 [9,] -0.071978515 [17,] -0.203589406 [25,]

　0.118056190

�。�2）平穩(wěn)序列

　 （3）非白噪聲序列

　Box-Pierce test data:

　rain X-squared = 0.2709, df = 3, p-value = 0.9654 X-squared = 7.7505, df = 6, p-value = 0.257 X-squared = 8.4681, df = 9, p-value = 0.4877 X-squared = 19.914, df = 12, p-value = 0.06873 X-squared = 21.803, df = 15, p-value = 0.1131 X-squared = 29.445, df = 18, p-value = 0.0432

　2.4 R 命令 Q_test<-function(n,r0){

　 k<-length(ro)

　 Q=0

　 P=0

　 for(i in 1:k) {

　 Q[i]<-n*sum(ro[1:i]^2)

　 P[i]<-1-pchisq(Q[i],df=i)

　 }

　 return(data.frame(Q,P)) } ro<-c(0.02,0.05,0.1,-0.02,0.05,0.01,0.12,-0.06,0.08,-0.05,0.02,-0.05) Q_test(100,ro) 答案：

　我們自定義函數(shù)，計算該序列各階延遲的 Q 統(tǒng)計量及相應(yīng) P 值。由于延遲 1-12 階 Q 統(tǒng)計量的 P 值均顯著大于 0.05，所以該序列為純隨機序列。

　2.5 R 命令 x<-ts(E2_5$x,star=c(2000,1)) par(mfrow=c(1,2)) plot(x,type="o") acf(x) for(i in 1:2) print(Box.test(x,lag=3*i))

　 答案 （1）繪制時序圖與自相關(guān)圖

　 （2）序列時序圖顯示出典型的周期特征，該序列非平穩(wěn) （3）該序列為非白噪聲序列

　Box-Pierce test data:

　x X-squared = 36.592, df = 3, p-value = 5.612e-08 X-squared = 84.84, df = 6, p-value = 3.331e-16

　2.6

　R 命令 x<-ts(E2_6$x) plot(x) acf(x) adf.test(x) for(i in 1:2) print(Box.test(x,lag=3*i)) y<-diff(x) adf.test(y) for(i in 1:2) print(Box.test(y,lag=3*i))

　答案 （1）如果是進行平穩(wěn)性圖識別，該序列自相關(guān)圖呈現(xiàn)一定的趨勢序列特征，可以視為非平穩(wěn)非白噪聲序列。

　如果通過 adf 檢驗進行序列平穩(wěn)性識別，該序列帶漂移項的 0 階滯后 P 值小于 0.05，可以視為平穩(wěn)非白噪聲序列

　Box-Pierce test data:

　x X-squared = 47.99, df = 3, p-value = 2.14e-10 X-squared = 60.084, df = 6, p-value = 4.327e-11

�。�2）差分序列平穩(wěn)，非白噪聲序列

　 Box-Pierce test

　data:

　y X-squared = 22.412, df = 3, p-value = 5.355e-05 X-squared = 27.755, df = 6, p-value = 0.0001045

　2.7 R 命令 x<-ts(E2_7$mortality) par(frown=c(1,2)) plot(x) acf(x,lag.max = 24) adf.test(x) for(i in 1:2) print(Box.test(x,lag=3*i)) adf.test(diff(x)) for(i in 1:2) print(Box.test(diff(x),lag=3*i))

　答案 （1）時序圖和自相關(guān)圖顯示該序列有趨勢特征，所以圖識別為非平穩(wěn)序列。

　 （2）單位根檢驗顯示帶漂移項 0 階延遲的 P 值小于 0.05，所以基于 adf 檢驗可以認(rèn)為該序列平穩(wěn)

�。�3）如果使用 adf 檢驗結(jié)果，認(rèn)為該序列平穩(wěn)，則白噪聲檢驗顯示該序列為非白噪聲序列 Box-Pierce test data:

　x X-squared = 60.252, df = 3, p-value = 5.193e-13 X-squared = 88.061, df = 6, p-value < 2.2e-16 如果使用圖識別認(rèn)為該序列非平穩(wěn)，那么一階差分后序列為平穩(wěn)非白噪聲序列

　Box-Pierce test data:

　diff(x) X-squared = 16.054, df = 3, p-value = 0.001106 X-squared = 20.969, df = 6, p-value = 0.001859

　 2.8 R 命令

　答案 （1）時序圖和自相關(guān)圖都顯示典型的趨勢序列特征

�。�2）單位根檢驗顯示該序列可以認(rèn)為是平穩(wěn)序列（帶漂移項一階滯后 P 值小于 0.05）

　 （3）一階差分后序列平穩(wěn)

　 第三章習(xí)題答案 3.1

　0101 ( ) 01 1 0.7tE x??? ? ?? ?（）

　2 211 12 ( ) 1.961 1 0.7tVar x?? ? ?? ?（ ）

　2 22 13 =0.7 0.49 ? ? ? ? （）

　121 22221110.49 0.7= 01 1 0.71?? ?????? ??（4）

　 3.2

　 1 1112 222 1 1 2 1 2(2)7= 0.5151 11= 0.3 0.515AR? ???? ?? ? ? ? ? ? ??? ? ?? ?? ? ?? ? ? ?? ? ?? ? ?? ? ? ?? ?? ?模型有：,2115? ?

　 3.3

　 1 201 2(1) (1 0.5 )(1 0.3 ) 0.8 0.15( ) 01t t t t t ttB B x x x xE x? ??? ?? ?? ? ? ? ? ? ?? ?? ?,22 1 2 1 212 ( )(1 )(1 )(1 )1 0.15=(1 0.15)(1 0.8 0.15)(1 0.8 0.15)1.98tVar x?? ? ? ? ???? ? ? ? ??? ? ? ? ??（ ）

　1122 1 1 23 1 2 2 10.83 = 0.701 1 0.150.8 0.7 0.15 0.410.8 0.41 0.15 0.7 0.22???? ? ? ?? ? ? ? ?? ?? ?? ? ? ? ? ?? ? ? ? ? ? ?（ ）

　11 122 2334 0.70.15=0? ?? ??? ?? ? ?（ ）

　3.4

　 1 21 1 11 00 1 1ARc ccc c? ? ? ? ? ??? ? ? ? ?? ?? ? ?? ??（）

�。� ）模型的平穩(wěn)條件是 11 21,2 1, 2k k kcc k?? ? ?? ???????? ? ??（ ）

　 3.5

　證明:

　該序列的特征方程為：3 20 c c ? ? ? ? ? ? ? ，解該特征方程得三個特征根：

　11 ? ? ，2c ? ? ，3c ? ??

　無論 c 取什么值，該方程都有一個特征根在單位圓上，所以該序列一定是非平穩(wěn)序列。證畢。

　3.6

　 （1）錯

　（2）錯 （3）對 （4）錯 （5）對

　3.7

　211 1 1 1 1211=0.4 0.4 0.4 0 21+ 2?? ? ? ? ???? ? ? ? ? ? ? ? ? ? 或者 所以該模型有兩種可能的表達式：11+2t t tx ? ??? 和1+2t t tx ? ??? 。

　3.8

　將1 2 310 0.5 0.8t t t t tx x C ? ? ?? ? ?? ? ? ? ? 等價表達為 2 321 0.810 (1 )1 0.5t t tB cBx aB bBB? ?? ?? ? ? ? ?? 則 ? ?2 3 22 31 0.8 (1 ) 1 0.5=1 ( 0.5) ( 0.5 ) 0.5B cB aB bB Ba B b a B bB? ? ? ? ? ?? ? ? ? ? 根據(jù)待定系數(shù)法：

　0.8= 0.5 0.30 0.5 00.5 0.15a ab bc b a c? ? ? ? ?? ? ? ?? ? ? ? 3.9 1 ( ) 0tE x ? （）

　2 22 ( ) 1 0.7 0.4 1.65tVar x ? ? ? ? （ ）

�。�3）10.7 0.7 0.40.591.65?? ? ?? ? ? ，20.40.241.65? ? ? ， 0, 3kk ? ? ?

　 3.10

�。�1）證明：因為對任意常數(shù) C，有 2 2( ) lim(1 )tkVar x kC????? ? ? ?

　所以該序列為非平穩(wěn)序列。

�。�2）1 1( 1)t t t t ty x x C ? ?? ?? ? ? ? ? ，則序列 ? ?ty 滿足如下條件：

　均值、方差為常數(shù)， ( ) 0tE y ? ,2 2( ) 1 ( 1)tVar y C?? ? ? ? ? ?? ? 自相關(guān)系數(shù)只與時間間隔長度有關(guān)，與起始時間無關(guān)

　121, 0, 21 ( 1)kCkC? ??? ? ?? ? 所以該差分序列為平穩(wěn)序列。

　3.11

　（1）非平穩(wěn)，（2）平穩(wěn)，（3）可逆，（4）不可逆，（5）平穩(wěn)可逆，（6）不平穩(wěn)不可逆

　3.12

　該模型的 Green 函數(shù)為：

　01 G ?

　1 1 0 10.6 0.3 0.3 G G ? ? ? ? ? ? ?

　1 11 1 1 10.3 0.6 , 2k kk kG G G k ? ?? ??? ? ? ? ?

　所以該模型可以等價表示為：100.3 0.6 kt t t kkx ? ??? ??? ? ??

　3.13 21 201 2( ) (1 0.5 ) =0.5 = 0.253( ) 41 1 0.5 0.25tB BE x? ??? ?? ? ? ? ?? ? ?? ? ? ?，

　 3.14

　證明：

　已知112? ? ，114? ? ，根據(jù) (1,1) ARMA 模型 Green 函數(shù)的遞推公式得：

　01 G ? ，21 1 0 1 10.5 0.25 G G ? ? ? ? ? ? ? ? ，1 11 1 1 1 1, 2k kk kG G G k ? ? ?? ??? ? ? ?

　01 ? ?

　52 2 3 211 1 1 12 2 4 50 11 1 1 11 42 42 2( 1)1 1 1120 111 70.271 261 11jj jj jjjj jG GG?? ? ?? ? ? ??? ? ???? ???? ?? ??? ?? ?? ? ?? ? ? ? ? ?? ?? ??? ?? ? ? ?1 1 10 0 01 1 12 2 20 0 0, 2j j k j j k j j kj j jk kj j jj j jG G G G G GkG G G?? ? ? ?? ? ?? ? ? ? ?? ? ?? ? ? ?? ? ?? ? ? ? ?? ? ?? ? ? 證畢。

　3.15

　 （1）成立

　 （2）成立

　 （3）成立

�。�4）成立

　3.16

　 該習(xí)題數(shù)據(jù)文件與 2.7 相同。該題問題設(shè)置有問題：是要問如果判斷該序列或差分序列是平穩(wěn)序列，那該平穩(wěn)序列具有 ARMA 族中哪個模型的特征。

�。�1）根據(jù) adf 檢驗結(jié)果可以認(rèn)為該序列平穩(wěn)。根據(jù)序列的自相關(guān)圖可以認(rèn)為是自相關(guān)系數(shù)拖尾。根據(jù)偏自相關(guān)圖，可以認(rèn)為是偏自相關(guān) 2 階截尾，所以該序列具有 AR（2）模型的特征。

　 （2）根據(jù)圖識別也可以認(rèn)為該序列不平穩(wěn)，對該序列進行一階差分。一階差分后序列可以視為平穩(wěn)序列，根據(jù)差分后序列的自相關(guān)圖可以認(rèn)為是自相關(guān)系數(shù) 1 階截尾，具有 MA（1）模型的特征。根據(jù)偏自相關(guān)圖，可以認(rèn)為是偏自相關(guān) 3 階截尾，具有 AR（3）模型的特征。具體哪個模型最適合擬合該序列，下一章介紹。

　3.17

　該習(xí)題數(shù)據(jù)文件與 2.8 相同。該題問題設(shè)置有問題：是要問如果判斷該序列或差分序列是平

　穩(wěn)序列，那該平穩(wěn)序列具有 ARMA 族中哪個模型的特征。

�。�1）根據(jù) adf 檢驗，該序列可以視為平穩(wěn)序列。自相關(guān)圖呈現(xiàn)拖尾屬性，偏自相關(guān)圖呈現(xiàn)1 階截尾特征，所以該序列呈現(xiàn)出 AR(1)模型特征。

　 （3）如果根據(jù)圖識別，可以認(rèn)為序列蘊含趨勢，視為非平穩(wěn)序列。一階差分后序列平穩(wěn)。一階差分后序列呈現(xiàn)自相關(guān)系數(shù) 2 階截尾，偏自相關(guān)系數(shù) 2 階截尾的特征。所以一階差分后序列具有 MA(2)模型特征，或 AR（2）模型特征。具體哪個模型最適合擬合該序列，下一章介紹。

　第四章習(xí)題答案 4.1 本題 R 代碼 讀入數(shù)據(jù)文件之后 x<-ts(E4_1$x) library(aTSA) plot(x) adf.test(x) Box.test(x) acf(x) pacf(x) fit<-arima(x,order=c(1,0,0),include.mean = F) fit ts.diag(fit) library(forecast) forecast(fit,h=60) 答案 （1）繪制時序圖（略）

�。�2）該序列為平穩(wěn)非白噪聲 （3）自相關(guān)圖拖尾，偏自相關(guān)圖一階截尾 （4）擬合 AR(1)模型

�。�5）五年預(yù)測值見預(yù)測輸出（略）

　4.2 本題 R 代碼 讀入數(shù)據(jù)文件之后 x<-ts(E4_2$x) library(aTSA) plot(x) adf.test(x)

　Box.test(x) acf(x) pacf(x) fit<-arima(x,order=c(0,0,1)) fit ts.diag(fit) library(forecast) forecast(fit,h=60) 答案 （1）繪制時序圖（略）

�。�2）該序列為平穩(wěn)非白噪聲序列 （3）自相關(guān)圖一階截尾，偏自相關(guān)圖拖尾 （4）擬合 MA(1)模型

�。�5）五年預(yù)測值見 R 輸出（略）

　4.3 本題 R 代碼 讀入數(shù)據(jù)文件之后 x<-ts(E4_3$x) library(aTSA) plot(x) adf.test(x) Box.test(x) acf(x) pacf(x) fit<-arima(x,order=c(2,0,0)) fit ts.diag(fit) library(forecast) forecast(fit,h=12)

　答案 （1）繪制時序圖（略）

�。�2）該序列為平穩(wěn)非白噪聲序列 （3）根據(jù)該序列自相關(guān)圖，可以視為：自相關(guān)圖一階截尾，或偏自相關(guān) 2 階截尾 （4）分別擬合 MA(1)模型和 AR(2)模型，兩個模型均參數(shù)顯著非零，殘差為檢驗為白噪聲序列，AIC 和 SBC 的結(jié)果幾乎相等，最后考慮白噪聲檢驗的 P 值，AR(2)模型的白噪聲檢驗 P 值

　更大，說明該模型對序列的相關(guān)信息提取更為充分，所以選擇 AR(2)模型作為最優(yōu)模型。

�。�5）基于 AR(2)模型未來一年預(yù)測值為

　4.4 本題 R 代碼 讀入數(shù)據(jù)文件之后 x<-ts(E4_4$x) library(aTSA) plot(x) adf.test(x) Box.test(x) acf(x) pacf(x) fit<-arima(x,order=c(0,0,1)) fit ts.diag(fit) library(forecast) forecast(fit,h=5) 答案 （1）繪制時序圖（略）

�。�2）該序列為平穩(wěn)非白噪聲序列 （3）根據(jù)該序列自相關(guān)圖，可以視為：自相關(guān)圖一階截尾，或偏自相關(guān) 2 階截尾，或自相關(guān)和偏自相關(guān)均拖尾 （4）分別擬合 MA(1)模型，AR(2)模型和 ARMA(1,1,) 。ARMA(1,1)模型參數(shù)不能拒絕參數(shù)為零的原假設(shè)，所以淘汰。MA(1)模型，AR(2)模型均參數(shù)顯著非零，殘差為檢驗為白噪聲序列，MA(1)模型 SBC 更小一點，所以選擇 MA(1)模型作為最優(yōu)模型。

　 （5）基于 MA(1)模型未來五年預(yù)測值為

　 4.5 （1）

　0 1123(1 ) 10 (1 0.3) 7ˆ 7 0.3 9.6 9.88ˆ 7 0.3 9.88 9.964ˆ 7 0.3 9.964 9.9892tttxxx? ? ????? ? ? ? ? ?? ? ? ?? ? ? ?? ? ? ? 2 2 2 2 2 43 0 1 2ˆ ( ) ( ) (1 0.3 +0.3 9=9.8829tVar x G G G???? ? ? ? ? ? ）

　所以3ˆ 95% 9.9892 1.96 9.8829 3.83,16.15tx ? ? 的 的置信區(qū)間等于 ，即( ）

�。�2）更新數(shù)據(jù)后 23ˆ 7 0.3 10.5 10.15ˆ 7 0.3 10.15 10.045ttxx??? ? ? ?? ? ? ? 2 2 2 23 0 1ˆ ( ) ( ) (1 0.3 9=9.81tVar x G G???? ? ? ? ? ）

　所以3ˆ 95% 10.045 1.96 9.819 3.91,16.18tx ? ? 的 的置信區(qū)間等于 ，即( ）

　4.6 本題 R 指令 x<-ts(E4_6$x) library(aTSA) plot(x) adf.test(x)

　Box.test(x) acf(x) pacf(x) fit<-arima(x,order=c(1,0,0)) fit ts.diag(fit) library(forecast) forecast(fit,h=5)

　答案 （1）平穩(wěn)非白噪聲序列 （2）自相關(guān)系數(shù)拖尾，偏自相關(guān)系數(shù) 1 階截尾，擬合 AR（1）模型

�。�3）未來 5 年的降雪量

　 4.7 本題 R 指令 讀入數(shù)據(jù)文件之后 x<-ts(E4_7$x) library(aTSA) plot(x) adf.test(x) Box.test(x) acf(x) pacf(x) fit<-arima(x,order=c(1,0,1)) fit ts.diag(fit) library(forecast) forecast(fit,h=5)

　答案 （1）平穩(wěn)非白噪聲序列 （2）自相關(guān)系數(shù)和偏自相關(guān)系數(shù)都拖尾，擬合 ARMA(1,1)模型

�。�3）未來 5 年的谷物產(chǎn)量預(yù)測：

　 4.8 本題 R 指令 讀入數(shù)據(jù)文件之后 x<-ts(E4_8$x) library(aTSA) plot(x) adf.test(x) Box.test(x) acf(x) pacf(x) fit<-arima(x,order=c(0,0,1)) fit ts.diag(fit) library(forecast) forecast(fit,h=1)

　答案 （1）平穩(wěn)非白噪聲序列 （2）自相關(guān)系數(shù) 1 階截尾，偏自相關(guān)系數(shù)拖尾，擬合 MA(1)模型

　 （3）下一刻的 95%置信區(qū)間為：（80.43481 90.92585）

　4.9 本題 R 指令 讀入數(shù)據(jù)文件之后 x<-ts(E4_9$x) library(aTSA) plot(x) adf.test(x) Box.test(x) acf(x) pacf(x) fit<-arima(x,order=c(4,0,1),transform.pars = F,fixed=c(0,0,NA,NA,NA,NA)) fit ts.diag(fit) library(forecast) forecast(fit,h=5)

　答案 （1）平穩(wěn)非白噪聲序列 （2）自相關(guān)系數(shù)和偏自相關(guān)系數(shù)都拖尾，擬合 ARMA(4,1)模型（該模型有部分系數(shù)不能顯著非零）。所以最好是擬合疏系數(shù) ARMA((3,4),1) 模型

�。�3）未來 5 年的預(yù)測

　 第 五 章習(xí)題答案 5.1 本題 R 指令 讀入數(shù)據(jù)文件之后 x<-ts(E5_1$volume,start=1949) par(mfrow=c(1,2)) plot(x) plot(diff(x)) adf.test(diff(x)) Box.test(diff(x)) acf(diff(x)) pacf(diff(x)) fit<-Arima(x,order=c(1,1,0),include.drift=T) fit tsdiag(fit) fore<-forecast(fit,h=5) fore par(mfrow=c(1,1)) plot(fore) lines(fore$fitted,col=2)

　答案：

　該序列 1 階差分后平穩(wěn)，根據(jù)差分后序列的自相關(guān)圖拖尾和偏自相關(guān)圖 1 階截尾特征，對該序列擬合 ARIMA(1,1,0)模型，模型參數(shù)如下：

　根據(jù)該模型，得到 2009-2013 年我國鐵路貨運量的預(yù)測值為：

　 5.2 本題 R 指令 讀入數(shù)據(jù)文件之后 x<-ts(E5_2$birth_rate,start=1750) plot(x) adf.test(x) Box.test(x) par(mfrow=c(1,2)) acf(x) pacf(x) fit<-arima(x,order=c(1,0,0)) fit

　答案 該序列 adf 檢驗平穩(wěn)，根據(jù)自相關(guān)圖和偏自相關(guān)圖特征，可以識別為偏自相關(guān)圖 1 階截尾，擬合 AR（1）模型。參數(shù)輸出結(jié)果如下：

　Forecasts from ARIMA(1,1,0) with drift1950 1960 1970 1980 1990 2000 201060000 100000 140000

　 5.3 本題 R 指令 讀入數(shù)據(jù)文件之后 x<-ts(E5_3$number,start=1867) par(mfrow=c(1,2)) plot(x) plot(diff(x)) adf.test(diff(x)) Box.test(diff(x)) acf(diff(x)) pacf(diff(x)) fit1<-Arima(x,order=c(3,1,0),include.drift=T) fit1 fit2<-arima(x,order=c(3,1,0),transform.pars = F,fixed=c(NA,0,NA)) fit2 tsdiag(fit2) fore<-forecast(fit2,h=7) fore par(mfrow=c(1,1)) plot(fore) lines(fore$fitted,col=2)

　 答案 （1）根據(jù)時序圖和 adf 檢驗結(jié)果可知，原序列非平穩(wěn)，一階差分后序列平穩(wěn)。

�。�2）一階差分后序列自相關(guān)圖拖尾，偏自相關(guān)圖 3 階截尾，所以對原序列擬合模型ARIMA(3,1,0)，參數(shù)輸出結(jié)果顯示 2 階自相關(guān)系數(shù)與均值都不顯著非零。所以最后擬合不帶均值項的疏系數(shù)模型 ARIMA((1,3),1,0)，最后參數(shù)估計值輸出結(jié)果如下

　 （3）基于模型對 1939—1945 年英國綿羊的數(shù)量的預(yù)測為：

　5_4 本題 R 指令 讀入數(shù)據(jù)文件之后 x1<-ts(E5_4$birth_rate,start=1980) x2<-ts(E5_4$mortality,start=1980) x3<-ts(E5_4$ngr,start=1980)

　par(mfrow=c(1,2)) plot(x1) plot(diff(x1)) adf.test(diff(x1)) for(i in 1:6) print(Box.test(diff(x1),i)) acf(diff(x1)) pacf(diff(x1)) fit11<-arima(x1,order=c(0,1,0)) fit12<-arima(x1,order=c(0,1,1)) fit13<-arima(x1,order=c(1,1,0)) c(AIC(fit11),AIC(fit12),AIC(fit13),BIC(fit11),BIC(fit12),BIC(fit13)) fit12 fore1<-forecast(fit12,h=10) fore1

　plot(x2) plot(diff(x2)) adf.test(diff(x2)) for(i in 1:6) print(Box.test(diff(x2),i))

　acf(diff(x2)) pacf(diff(x2)) fit21<-arima(x2,order=c(0,1,0)) fit22<-arima(x2,order=c(0,1,1)) fit23<-arima(x2,order=c(1,1,0)) c(AIC(fit21),AIC(fit22),AIC(fit23),BIC(fit21),BIC(fit22),BIC(fit23)) fore2<-forecast(fit21,h=10) fore2

　plot(x3) plot(diff(x3)) adf.test(diff(x3)) for(i in 1:6) print(Box.test(diff(x3),i)) acf(diff(x3)) pacf(diff(x3)) fit31<-arima(x3,order=c(0,1,0)) fit32<-arima(x3,order=c(0,1,1)) fit33<-arima(x3,order=c(1,1,0)) c(AIC(fit31),AIC(fit32),AIC(fit33),BIC(fit31),BIC(fit32),BIC(fit33))

　fit32 fore3<-forecast(fit32,h=10) fore3

　答案 （1）我國人口出生率、死亡率和自然增長率序列均有顯著線性趨勢，為典型非平穩(wěn)序列。

�。�2）根據(jù)圖檢驗和 adf 檢驗結(jié)果顯示，我國人口出生率、死亡率和自然增長率序列的一階差分后序列均平穩(wěn)。

�。�3）我國人口出生率可以擬合 ARIMA(0,1,0)模型，也可以擬合 ARIMA(0,1,1)模型。其中ARIMA（0,1,1）模型 AIC,SBC 更小一點。得到的擬合模型參數(shù)為

　基于 ARIMA（0,1,1）模型，預(yù)測 1998 年之后人口出生率未來十年的預(yù)測值為：

　 （4）我國人口死亡率可以擬合隨機游走模型 ARIMA(0,1,0)�；�于 ARIMA（0,1,0）模型，預(yù)測 1998 年之后人口死亡率未來十年的預(yù)測值為：

　 （5）（5）我國人口自然增長率可以擬合 ARIMA(0,1,0)模型，ARIMA(1,1,0)模型，也可以擬合 ARIMA(0,1,1)模型。其中 ARIMA（0,1,1）模型 AIC,SBC 更小一點。得到的擬合模型參數(shù)為

　基于 ARIMA（0,1,1）模型，預(yù)測 1998 年之后人口自然增長率未來十年的預(yù)測值為：

　 5_5 本題 R 指令 讀入數(shù)據(jù)文件之后

　答案 （1）玉米價格、玉米產(chǎn)量、工人薪水、生豬價格、生豬產(chǎn)量序列都不平穩(wěn)， （2）玉米價格、玉米產(chǎn)量、工人薪水、生豬價格、生豬產(chǎn)量序列一階差分都實現(xiàn)平穩(wěn)。

�。�3）-（4）分序列分析結(jié)果：

�、� 對玉米價格序列擬合疏系數(shù)模型 ARIMA(0,1,(7))，參數(shù)估計值如下：

　對玉米價格序列進行 10 期預(yù)測，預(yù)測值如下：

　對玉米價格序列的擬合與預(yù)測效果圖如下：

　 ② 對玉米產(chǎn)量序列擬合模型 ARIMA(0,1,1)，參數(shù)估計值如下：

　 對玉米產(chǎn)量序列進行 10 期預(yù)測，預(yù)測值如下：

　對玉米產(chǎn)量列的擬合與預(yù)測效果圖如下：

�、� 對工人薪水序列擬合疏系數(shù)模型 ARIMA(（1,13）,1,0)，參數(shù)估計值如下：

　對工人薪水序列進行 10 期預(yù)測，預(yù)測值如下：

　對工人薪水序列的擬合與預(yù)測效果圖如下：

�、� 對生豬價格序列擬合模型 ARIMA(3,1,0)，參數(shù)估計值如下：

　 對生豬價格序列進行 10 期預(yù)測，預(yù)測值如下：

　對生豬價格序列的擬合與預(yù)測效果圖如下：

�、� 對生豬產(chǎn)量序列擬合疏系數(shù)模型 ARIMA((2,3,5),1,0)，參數(shù)估計值如下：

　 對生豬產(chǎn)量序列進行 10 期預(yù)測，預(yù)測值如下：

　對生豬產(chǎn)量序列的擬合與預(yù)測效果圖如下：

　 第六章習(xí)題答案 6_1 R 命令 x<-ts(E6_1$x,start=c(1962,1),frequency = 12) plot(x) fit1<-decompose(x) fit1$trend plot(fit1$trend) fit1$figure plot(1:12,fit$figure,type="o",xlab="Month") fit2<-HoltWinters(x) fit2 library(forecast) fore<-forecast(fit2,h=24) fore plot(fore) lines(fore$fitted,col=2) acf(diff(x)) pacf(diff(x)) fit3<-arima(x,order=c(1,1,0),seasonal = list(order=c(0,1,1),period=12)) fit3 plot(x) lines(fitted(fit3),col=2)

　c(mean(na.exclude(fit1$random)^2),mean((x-fit2$fit[,1])^2),mean(fit3$residuals^2))

　答案

�。�1）該序列受到趨勢、季節(jié)、隨機波動等三個因素的影響。由于季節(jié)性沒有隨趨勢變化而發(fā)生顯著變化，可以選擇加法模型分解各因素。

�。�2）提取該序列的趨勢效應(yīng)（趨勢數(shù)據(jù)見 fit$trend，具體數(shù)據(jù)略）， 趨勢效應(yīng)圖如下：

�。�3）提取該序列的季節(jié)效應(yīng)（季節(jié)效應(yīng)數(shù)據(jù)見 fit$figure，具體數(shù)據(jù)略））。季節(jié)效應(yīng)圖如下：

�。�4）使用 Holt-Winters 加法模型擬合該序列的發(fā)展，幵做 2 年期預(yù)測。相關(guān)參數(shù)如下：

　Smoothing parameters:

　alpha: 0.68933

　beta : 0

　gamma: 0.8362592

　Coefficients:

　 [,1] a

　 885.775547 b

　 1.278118 s1

　-16.743296 s2

　-59.730034 s3

　 47.492731 s4

　 56.203890

　s5

　115.537545 s6

　 84.554817 s7

　 39.580306 s8

　 -4.702033 s9

　-54.554684 s10 -51.582594 s11 -85.953466 s12 -42.907363 2 年預(yù)測值如下：

　 預(yù)測效果圖：

　 （5）用 ARIMA 季節(jié)模型擬合幵預(yù)測該序列的發(fā)展。

　擬合模型：ARIMA(1,1,0)×(0,1,1) 12, 模型參數(shù)估計輸出如下：

　擬合不預(yù)測效果圖：

　 (6)比較分析上面使用過的三種模型的擬合精度。

　由于因素分解方法和 Holt-Winters 三參數(shù)指數(shù)平滑法丌提供極大似然函數(shù)值，所以無法計算AIC 和 BIC 信息量的值。如果單純考察方差，本題是因素分解方法擬合的模型方差最小。

　c(mean(na.exclude(fit1$random)^2),mean((x-fit2$fit[,1])^2),mean(fit3$residuals^2)) [1] 45.04238 67.52938 48.64248 如果基于信息量最小法則，通常情況下是 ARIMA , 模型的擬合精度最高。

　6_2

　R 命令 x<-ts(E6_2$x,start=c(1973,1),frequency = 12) plot(x) fit1<-decompose(x) fit1$trend plot(fit1$trend) fit1$figure plot(1:12,fit$figure,type="o",xlab="Month") fit2<-HoltWinters(x) fit2 library(forecast) fore<-forecast(fit2,h=24) fore plot(fore) lines(fore$fitted,col=2) acf(diff(x)) pacf(diff(x)) fit3<-arima(x,order=c(1,1,0),seasonal = list(order=c(0,1,1),period=12)) fit3 plot(x) lines(fitted(fit3),col=2) c(mean(na.exclude(fit1$random)^2),mean((x-fit2$fit[,1])^2),mean(fit3$residuals^2))

　答案 （1 ）該序列受到季節(jié)、趨勢和隨機波動的影響�？梢詳M合加法因素分解模型。

�。�2）提取該序列的趨勢效應(yīng)（數(shù)據(jù)略）， 趨勢效應(yīng)圖如下

�。�3）提取該序列的季節(jié)效應(yīng)（數(shù)據(jù)略），季節(jié)效應(yīng)圖如下

　 （4）用指數(shù)平滑法對該序列做 2 年期預(yù)測（預(yù)測值略，預(yù)測效果圖如下）

　 （5）用 ARIMA 季節(jié)模型擬合幵預(yù)測該序列的發(fā)展。擬合模型：ARIMA(1,1,0)×(0,1,1) 12, 模型參數(shù)估計輸出如下：

　Coefficients:

　 ar1

　 sma1

　 -0.3339

　-0.5960 s.e.

　 0.1224

　 0.1817 （6）比較分析上面使用過的三種模型的擬合精度。如果單純考察方差，本題是因素分解方法擬合的模型方差最小。

　c(mean(na.exclude(fit1$random)^2),mean((x-fit2$fit[,1])^2),mean(fit3$residuals^2)) [1]

　48265.31 144573.60

　83713.63 如果基于信息量最小法則，通常情況下是 ARIMA , 模型的擬合精度最高。

　6_3 使用 M 2 × 4 中心移動平均做預(yù)測,則 4 3 2 13 2 14 3 2 11 1 1 1 1ˆ (1)8 4 4 4 81 1 1 1 1ˆ ˆ (2) (1)8 4 4 4 81 5 9 9 1764 32 32 32 64t t t t t tt t t t t tt t t t tx x x x x xx x x x x xx x x x x? ? ? ?? ? ?? ? ? ?? ? ? ? ?? ? ? ? ?? ? ? ? ? 求在 2 期預(yù)測值 ˆ (2)tx中3 tx?前面的系數(shù)是532，4 tx?前面的系數(shù)是164。

　6_4 解下面的方程組，得到 0.4 ? ?

　5.25 5(1 )5.26 5.5 (1 )tt? ? ?? ? ?? ? ? ??? ? ?? 6_5 根據(jù)指數(shù)平滑的定義有（1）式成立，（1）式等號兩邊同乘 (1 ) ? ? 有（2）式成立

　2 32 3( 1) (1 ) ( 2) (1 ) ( 2) (1 ) (1)(1 ) (1 ) ( 1) (1 ) ( 2) (1 ) (2)ttx t t t tx t t t? ? ? ? ? ? ?? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? （1）-（2）得 22(1 ) (1 )(1 ) (1 )1ttx tx tt? ? ? ? ? ?? ???? ? ? ? ? ?? ? ? ? ? ??? ? 則1lim lim 1tt ttxt t???? ??? ? ??? ?? ?? ?? ?? ?。

　6_6 R 命令 x<-ts(E6_6$x,start=c(1973,1),frequency = 12) plot(x) fit1<-HoltWinters(x,gamma=F) fit1 acf(diff(x)) pacf(diff(x)) fit2<-arima(x,order=c(1,1,0)) fit2 c(mean((x-fit1$fit[,1])^2),mean(fit2$residuals^2)) library(forecast) fore<-forecast(fit2,h=24) plot(fore) lines(fore$fitted,col=2)

　答案 (1)1949—2008 年我國人口總數(shù)序列為顯著的線性遞增序列，可以使用 holt 兩參數(shù)指數(shù)平滑法進行趨勢擬合和預(yù)測，或使用 ARIMA（1,1,0）模型進行擬合和預(yù)測。

　(2) ARIMA（1,1,0）模型擬合效果最優(yōu)的模型相對最優(yōu),該模型參數(shù)估計結(jié)果如下， Call: arima(x = x, order = c(1, 1, 0))

　Coefficients:

　ar1

　 0.9473 s.e.

　0.0345 使用該模型對 2009—2016 年中國人口總數(shù)進行預(yù)測（預(yù)測值略），預(yù)測效果圖如下

　6_7 R 命令 x<-ts(E6_7$x,start=c(1948,1),frequency = 4) plot(x) fit1<-HoltWinters(x,gamma=F) fit1 library(aTSA) adf.test(diff(x)) adf.test(diff(diff(x))) acf(diff(diff(x))) pacf(diff(diff(x))) fit2<-Arima(x,order=c(0,2,16),transform.pars=F, fixed=c(NA,0,0,0,0,0,0,0,0,0,0,0,0,0,0,NA)) fit2

　c(mean((x-fit1$fit[,1])^2),mean(fit2$residuals^2)) library(forecast) fore<-forecast::forecast(fit2,h =20) plot(fore) lines(fore$fitted,col=2)

　答案 (1) 該序列時序圖顯示典型的非線性趨勢，可以用 Holt-Winters 指數(shù)平滑模型擬合也可以用 2階差分的 ARIMA 模型擬合。時序圖如下

　(2) 擬合效果是疏系數(shù) ARIMA(0,

　2，(1, 16))相對較優(yōu)。該模型的參數(shù)估計情況為

　（3）5 年的預(yù)測（預(yù)測值略），擬合與預(yù)測效果圖如下：

　6_8 R 命令 x<-ts(E6_8$x,start=c(1980,1),frequency = 12) plot(x) fit1<-decompose(x) plot(fit1) fit2<-HoltWinters(x) fit2 library(forecast) fore1<-forecast::forecast(fit2,h =60) plot(fore1) lines(fore1$fitted,col=2) library(aTSA) adf.test(diff(diff(x,12))) fit3<-Arima(x,order=c(2,1,0),seasonal = list(order=c(0,1,1),period=12)) fit3 fore2<-forecast::forecast(fit3,h =60) plot(fore2) lines(fore2$fitted,col=2)

　答案 (1) 時序圖顯示該序列時序圖有季節(jié)效應(yīng)、趨勢（循環(huán)）效應(yīng)和隨機波動。

　 (2) 使用因素分解加法模型對該序列進行因素分解

　 (3) 可以選擇 Holt-Winters 模型對該序列進行擬合與預(yù)測，也可以選擇 ARIMA 模型對該序列進行擬合與預(yù)測。

　方法一：使用 Holt-Winters 模型對該序列進行預(yù)測，相關(guān)參數(shù)如下

　Holt-Winters exponential smoothing with trend and additive seasonal component.

　Call: HoltWinters(x = x)

　Smoothing parameters:

　alpha: 0.3074109

　beta : 0.003657678

　gamma: 0.4375144

　Coefficients:

　 [,1] a

　 100572.16590 b

　211.71046 s1

　 486.50198 s2

　 -2040.98232 s3

　 -4953.90627 s4

　1931.00406 s5

　-14281.72507 s6

　-11711.73657 s7

　3014.95841 s8

　 -7888.68734 s9

　4979.71026 s10

　-88.60042 s11

　-5717.59270 s12

　 -254.48340

　預(yù)測效果圖

　方法二：使用 ARIMA(2,1,0) ? (0,1,1) 12 對該序列進行擬合與預(yù)測，相關(guān)參數(shù)如下 Series: x

　ARIMA(2,1,0)(0,1,1)[12]

　Coefficients:

　 ar1

　ar2

　 sma1

　 -0.6801

　-0.5520

　-0.7165 s.e.

　 0.0668

　 0.0704

　 0.1065

　預(yù)測效果圖

　6_9 R 命令 x<-ts(E6_9$x,start=c(1963,1),frequency = 12) plot(x) fit1<-HoltWinters(x,seasonal ="multi") fit1 library(forecast) fore1<-forecast::forecast(fit1,h =60) plot(fore1) lines(fore1$fitted,col=2) library(aTSA) adf.test(diff(diff(x,12))) fit2<-Arima(x,order=c(1,1,1),seasonal = list(order=c(0,1,1),period=12)) fit2 fore2<-forecast::forecast(fit2,h =60) plot(fore2) lines(fore2$fitted,col=2)

　答案

　(1) 考察該序列時序圖，時序圖顯示該序列具有顯著的線性遞增趨勢和年度季節(jié)效應(yīng)。

　(2)根據(jù)該序列呈現(xiàn)的規(guī)律,可以使用 Holt-Winters 模型對該序列進行擬合與預(yù)測，也可以選擇 ARIMA 模型對該序列進行擬合與預(yù)測。

　(3)選擇擬合效果相對最好的模型 ,預(yù)測該序列未來 3 年的旅館入住情況（沒有唯一解，給出可能參考答案）

　方法一：可以基于 Holt-Winters 模型對該序列進行擬合與預(yù)測 參數(shù)估計如下：

　Holt-Winters exponential smoothing with trend and multiplicative seasonal component. Call: HoltWinters(x = x, seasonal = "multi")

　Smoothing parameters:

　alpha: 0.01567674

　beta : 0.008068444

　gamma: 0.4392594

　Coefficients:

　[,1] a

　 875.5123349 b

　 1.9568538 s1

　0.9301067 s2

　0.8613974 s3

　0.8741871 s4

　0.9798006 s5

　0.9624445

　s6

　1.0964284 s7

　1.2846348 s8

　1.3110818 s9

　1.0004965 s10

　 0.9987417 s11

　 0.8628932 s12

　 0.9793806 擬合與預(yù)測效果圖如下

　方法二：基于 ARIMA(1,1,1) ? (0,1,1) 12 模型對該序列進行擬合與預(yù)測，參數(shù)估計如下：

　Series: x

　ARIMA(1,1,1)(0,1,1)[12]

　 Coefficients:

　ar1

　ma1

　 sma1

　 0.1873

　-1.0000

　-0.4046 s.e.

　0.0822

　 0.0181

　 0.0705 擬合與預(yù)測效果圖如下

　 第七章習(xí)題答案 7_1 R 程序 x<-ts(E7_1$x) y<-ts(E7_1$y) par(mfrow=c(1,2)) plot(x) plot(y,col=2) library(aTSA) adf.test(x) for(i in 1:6) print(Box.test(x,lag=3*i)) adf.test(y) for(i in 1:6) print(Box.test(y,lag=3*i)) 答案 (1)兩個序列皆平穩(wěn)。

　(2)兩個序列皆為白噪聲序列。

　(3)因為這兩個變量都具有白噪聲屬性，所以不考慮建立協(xié)整方程。

　7_2 R 程序 library(forecast) library(aTSA) x<-ts(E7_2$x) y<-ts(E7_2$y) plot(x) lines(y,col=2) x1<-x[-c(1,2)] y1<-y[-c(37,38)] fit<-Arima(y1,xreg=x1 ,order=c(2,0,0)) fit tsdiag(fit) acf(diff(x,10)) pacf(diff(x,10)) fitx<-arima(x,order=c(2,0,0),seasonal = list(order=c(0,1,1),period=10)) fitx ts.diag(fitx) forex<-forecast::forecast(fitx,h=14) forex plot(forex) lines(forex$fitted,col=2) forey<-forecast::forecast(fit,xreg=forex$mean,h=14) forey plot(forey) lines(forey$fitted,col=2) 答案 （1）掠食者和被掠食者數(shù)量都呈現(xiàn)出顯著的周期特征，兩個序列均為非平穩(wěn)序列。但是掠食者和被掠食者延遲2階序列具有協(xié)整關(guān)系。協(xié)整方程參數(shù)如下：

　22166.9407 0.11431 1.1892 0.5688t t ty xB B??? ? ?? ? （2）被掠食者擬合乘積模型：8(2,0,0) (0,1,1) ARIMA ? ,模型口徑為: 101021+=1 1.109 +0.4791t tBxB B? ?? 一周的預(yù)測值為

　 預(yù)測效果圖如下

　掠食者的擬合模型為：

　22166.9407 0.11431 1.1892 0.5688t t ty xB B??? ? ?? ? 一周的預(yù)測值為：

　預(yù)測效果圖如下

　8_3 R 程序 library(forecast) library(aTSA) x<-ts(E7_3$x,start=1950) y<-ts(E7_3$y,start=1950) plot(x) lines(y,col=2) acf(diff(diff(x))) pacf(diff(diff(x))) fitx<-arima(x,order=c(0,2,0)) fitx ts.diag(fitx) acf(diff(diff(y))) pacf(diff(diff(y))) fity<-arima(y,order=c(2,2,0),transform.pars = F,fixed=c(0,NA)) fity ts.diag(fity) coint.test(y,x) fit<-Arima(y,xreg=x,order=c(2,0,0)) fit tsdiag(fit) ecm(y,x) 答案 (1) 進口總額和出口總額序列均丌平穩(wěn)，均為 2 階單整序列

　(2)

　出口總額序列擬合模型 ARIMA(0,2,0)：

　2t tx ? ? ?

　進口總額序列擬合模型：ARIMA((2),2,0):

　(3) 這兩個序列具有協(xié)整關(guān)系

　 (4)協(xié)整模型如下：

　 （5）構(gòu)造該協(xié)整模型的誤差修正模型：

　 8_4 R 程序 library(forecast) library(aTSA) x<-ts(E7_4$x,start=1979) y<-ts(E7_4$y,start=1979) plot(x) plot(y) acf(diff(diff(x))) pacf(diff(diff(x))) fitx<-arima(x,order=c(3,2,0),transform.pars = F,fixed=c(0,0,NA)) fitx ts.diag(fitx) forex<-forecast::forecast(fitx,h=5) forex acf(diff(diff(y))) pacf(diff(diff(y))) fity<-arima(y,order=c((2),2,0),transform.pars = F,fixed=c(0,NA)) fity ts.diag(fity) forey<-forecast::forecast(fity,h=5) forey coint.test(y,x) 答案 （1）消費品零售總額序列擬合 ARIMA((3),2,0)模型，5 年預(yù)測值為：

　2211 0.403t txB? ? ??2269.5097+0.83581 1.2583 0.7701tt ty xB B?? ?? ?10.70185 0.37337t t ty x ECM?? ? ? ?

　 Point Forecast

　Lo 80

　Hi 80

　Lo 95

　Hi 95 2015

　 300674.5 297983.1 303366.0 296558.3 304790.7 2016

　 332630.6 326612.3 338648.8 323426.4 341834.7 2017

　 362786.5 352716.0 372856.9 347385.0 378187.9 2018

　 392800.3 377014.0 408586.6 368657.2 416943.4 2019

　 424456.7 401808.6 447104.7 389819.5 459093.9 國內(nèi)生產(chǎn)總值序列擬合 ARIMA((2),2,0),5 年預(yù)測值為：

　Point Forecast

　Lo 80

　Hi 80

　Lo 95

　Hi 95 2015

　 682666.3 671838.8 693493.8 666107.0 699225.5 2016

　 731554.3 707343.3 755765.4 694526.7 768581.9 2017

　 781093.1 744033.6 818152.6 724415.5 837770.7 2018

　 829667.2 779061.8 880272.7 752272.9 907061.5 2019

　 877975.4 811780.8 944170.1 776739.4 979211.4 （2）這兩個序列之間不存在協(xié)整關(guān)系。

�。�3）沒有協(xié)整方程,所以無法寫出協(xié)整方程 （4）因為這兩個變量不存在協(xié)整關(guān)系，所以無法直接分析如果中國國內(nèi)社會消費品零售總額增長 1%,對國內(nèi)生產(chǎn)總值有什么影響。需要引入其他變量，讓整個系統(tǒng)構(gòu)成協(xié)整關(guān)系之后，才能準(zhǔn)確分析這個問題。

　8_5 R 程序 library(forecast) library(aTSA) Seatbelts drivers<-Seatbelts[,2] front<-Seatbelts[,3] rear<-Seatbelts[,4] kms<-Seatbelts[,5] petroprice<-Seatbelts[,6] law<-Seatbelts[,8] X<-matrix(Seatbelts[,c(5,6,8)],ncol=3) coint.test(drivers,X) fit1<-arima(drivers,xreg=data.frame(petroprice,law),order=c(2,0,0),

　 seasonal =list(order=c(1,1,0),period=12) ) fit1 ts.diag(fit1) coint.test(front,X) fit2<-arima(front,xreg=data.frame(petroprice,law),order=c(2,0,0),

　 seasonal =list(order=c(1,1,0),period=12) ) fit2 ts.diag(fit2)

　coint.test(rear,X) fit3<-arima(rear,xreg=data.frame(kms,petroprice),order=c(1,0,0),

　 seasonal =list(order=c(1,1,0),period=12) ) fit3 ts.diag(fit3) Y<-matrix(Seatbelts[,c(3,4)],ncol=2) coint.test(drivers,Y) 答案 (1) 安全帶強制法律的執(zhí)行對司機傷亡數(shù)據(jù)有顯著的干預(yù)作用。

　(2) 司機傷亡數(shù)據(jù)不汽油價格及安全帶強制法律的執(zhí)行之間具有協(xié)整關(guān)系，具體協(xié)整方程如下：

　(3) 安全帶強制法律的執(zhí)行對前座乘客傷亡數(shù)據(jù)有顯著的干預(yù)作用。

　(4) 前座乘客傷亡數(shù)據(jù)汽油價格及安全帶強制法律的執(zhí)行之間具有協(xié)整關(guān)系，具體協(xié)整方程如下

　(5)安全帶強制法律的執(zhí)行對后座乘客傷亡數(shù)據(jù)沒有顯著的干預(yù)作用。

　(6)后座乘客傷亡數(shù)據(jù)不行駛里程、汽油價格之間具有協(xié)整關(guān)系，具體協(xié)整方程如下：

　(7) 司機傷亡數(shù)據(jù)、前座乘客傷亡數(shù)據(jù)和后座乘客傷亡數(shù)據(jù)之間具有協(xié)整關(guān)系

　8_6 R 程序 library(forecast) library(aTSA) x1<-ts(E7_6$maize_price,start=1967); x2<-ts(E7_6$maize_yield,,start=1967); x3<-ts(E7_6$pig_price,start=1967); x4<-ts(E7_6$pig_yield,start=1967); x5<-ts(E7_6$salary,start=1967) plot(x1) ? ?? ?122 12313 6053.91 0.25 0.23 1 0.45tt t tdrivers law petropriceB B B?? ?? ? ?? ? ?? ? ? ?12120.014 710.21 0.15 1 0.36tt t trear kms petropriceB B?? ?? ? ?? ?? ?? ?122 12217.5 3441.21 0.27 0.15 1 0.34tt t tfront law petropriceB B B?? ?? ? ?? ? ?

　adf.test(x1) adf.test(diff(x1)) plot(x2) adf.test(x2) adf.test(diff(x2)) plot(x3) adf.test(x3) adf.test(diff(x3)) plot(x4) adf.test(x4) adf.test(diff(x5)) plot(x5) adf.test(x5) adf.test(diff(x5)) X<-data.frame(x1,x2,x3,x4) X<-as.matrix(X) coint.test(x5,X,d=1) fit<-arima(x5,xreg=data.frame(x1,x3,x4),order=c(3,0,0),transform.pars = T

　,fixed=c(NA,0,NA,NA,NA,NA,NA)) fit ts.diag(fit) ecm(x5,as.matrix(data.frame(x1,x3,x4))) 答案 (1)這五個變量均為 1 階單整序列 (2)農(nóng)場工人薪水為果，玉米價格、玉米產(chǎn)量、生豬價格、生豬產(chǎn)量為因。

　(3)分析農(nóng)場工人平均工資、農(nóng)場作物、家畜的價格不供應(yīng)鏈之間具有協(xié)整關(guān)系。擬合協(xié)整模型如下（參數(shù)含義解釋略）

　誤差修正模型如下（參數(shù)含義解釋略）

　 8_7 library(forecast) library(aTSA) logM1<-ts(E7_7$logM1,start=c(1954,4),frequency=4); logGDP<-ts(E7_7$logGNP,start=c(1954,4),frequency=4); 1=0.04 +0.21 +0.31 0.14t t t t tECM?? ? ? ? ? 工人薪水 玉米價格 生豬價格 生豬產(chǎn)量3=4.0423+0.0486 +0.144 +0.243+1 1.19 0.21t t t ttB B?? ?工人薪水 玉米價格 生豬價格 生豬產(chǎn)量

　sr<--ts(E7_7$sr,start=c(1954,4),frequency=4); lr<-ts(E7_7$lr,start=c(1954,4),frequency=4); plot(logM1) adf.test(logM1) adf.test(diff(logM1)) acf(diff(logM1)) pacf(diff(logM1)) fit1<-arima(logM1,order=c(1,1,0)) fit1 ts.diag(fit1) plot(logGDP) adf.test(logGDP) adf.test(diff(logGDP)) acf(diff(logGDP)) pacf(diff(logGDP)) fit2<-arima(logGDP,order=c(2,1,1)) fit2 ts.diag(fit2) plot(sr) adf.test(sr) adf.test(diff(sr)) acf(diff(sr)) pacf(diff(sr)) fit3<-arima(sr,order=c(7,1,1),transform.pars = F, fixed=c(0,NA,0,0,0,0,NA,NA)) fit3 ts.diag(fit3) plot(lr) adf.test(lr) adf.test(diff(lr) acf(diff(lr)) pacf(diff(lr)) fit4<-arima(lr,order=c(14,1,0),transform.pars = F, fixed=c(NA,0,0,0,0,0,0,0,0,0,0,0,0,NA)) fit4 ts.diag(fit4) lagGNP<-lag(logGDP) y<-logGDP[1:134] X<-data.frame(logM1[1:134],sr[1:134],lr[1:134],lagGNP[2:135]) X<-as.matrix(X) coint.test(y,X) fit<-arima(y,xreg=data.frame(sr[1:134],lr[1:134],lagGNP[2:135])) fit ts.diag(fit) ecm(y,as.matrix(data.frame(sr[1:134],lr[1:134],lagGNP[2:135]))) 答案

　(1) 分別繪制這四個序列的時序圖（圖略）,這四個序列都有趨勢。它們均為 1 階單整序列。分別擬合單變量 ARIMA 模型如下：

　Logm1 序列：

　 Loggnp 序列：

　 短期利率序列：

　 長期利率序列：

　 (2)考察這四個變量的 Granger 因果關(guān)系：loggnp 為果，其他三個序列為因。

　(3)協(xié)整模型為

　 （4）誤差修正模型為 -1 1log 0.50876 0.1372 0.97062 log 0.8944t t tGNP sr lr GNP ECM?? ?? ? ? ? ? ? ? log 11 0.58061ttMB?? ??21 0.994log1 1.318 0.319t tBGNPB B??? ?? ?2 71+0.3271+0.253 0.3t tBsrB B? ? ??141l1 0.234 0.253t trB B? ? ?? ?1log 0.0683 0.9893log 0.3257 0.1650t t t tGNP GNP sr lr?? ? ? ?1logM11 0.5811t tB? ? ??

相關(guān)熱詞搜索：習(xí)題 序列 人大