人大版時(shí)間序列分析基于R（第2版）習(xí)題答案

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

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

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

　答案：

　（1）非平穩(wěn)，有典型線性趨勢(shì)

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

　（3）典型的具有單調(diào)趨勢(shì)的時(shí)間序列樣本自相關(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)，時(shí)序圖如下

　（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)典型的長(zhǎng)期趨勢(shì)與周期并存的特征

　 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ù)，計(jì)算該序列各階延遲的 Q 統(tǒng)計(jì)量及相應(yīng) P 值。由于延遲 1-12 階 Q 統(tǒng)計(jì)量的 P 值均顯著大于 0.05，所以該序列為純隨機(jī)序列。

　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）繪制時(shí)序圖與自相關(guān)圖

　 （2）序列時(shí)序圖顯示出典型的周期特征，該序列非平穩(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）如果是進(jìn)行平穩(wěn)性圖識(shí)別，該序列自相關(guān)圖呈現(xiàn)一定的趨勢(shì)序列特征，可以視為非平穩(wěn)非白噪聲序列。

　如果通過(guò) adf 檢驗(yàn)進(jìn)行序列平穩(wěn)性識(shí)別，該序列帶漂移項(xiàng)的 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）時(shí)序圖和自相關(guān)圖顯示該序列有趨勢(shì)特征，所以圖識(shí)別為非平穩(wěn)序列。

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

�。�3）如果使用 adf 檢驗(yàn)結(jié)果，認(rèn)為該序列平穩(wěn)，則白噪聲檢驗(yà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 如果使用圖識(shí)別認(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）時(shí)序圖和自相關(guān)圖都顯示典型的趨勢(shì)序列特征

�。�2）單位根檢驗(yàn)顯示該序列可以認(rèn)為是平穩(wěn)序列（帶漂移項(xiàng)一階滯后 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 ? ? ? ? ? ? ? ，解該特征方程得三個(gè)特征根：

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

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

　3.6

　 （1）錯(cuò)

�。�2）錯(cuò) （3）對(duì) （4）錯(cuò) （5）對(duì)

　3.7

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

　3.8

　將1 2 310 0.5 0.8t t t t tx x C ? ? ?? ? ?? ? ? ? ? 等價(jià)表達(dá)為 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）證明：因?yàn)閷?duì)任意常數(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ù)只與時(shí)間間隔長(zhǎng)度有關(guān)，與起始時(shí)間無(wú)關(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 ? ?? ??? ? ? ? ?

　所以該模型可以等價(jià)表示為：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 相同。該題問(wèn)題設(shè)置有問(wèn)題：是要問(wèn)如果判斷該序列或差分序列是平穩(wěn)序列，那該平穩(wěn)序列具有 ARMA 族中哪個(gè)模型的特征。

�。�1）根據(jù) adf 檢驗(yàn)結(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ù)圖識(shí)別也可以認(rèn)為該序列不平穩(wěn)，對(duì)該序列進(jìn)行一階差分。一階差分后序列可以視為平穩(wěn)序列，根據(jù)差分后序列的自相關(guān)圖可以認(rèn)為是自相關(guān)系數(shù) 1 階截尾，具有 MA（1）模型的特征。根據(jù)偏自相關(guān)圖，可以認(rèn)為是偏自相關(guān) 3 階截尾，具有 AR（3）模型的特征。具體哪個(gè)模型最適合擬合該序列，下一章介紹。

　3.17

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

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

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

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

　第四章習(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）繪制時(shí)序圖（略）

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

�。�5）五年預(yù)測(cè)值見(jiàn)預(yù)測(cè)輸出（略）

　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）繪制時(shí)序圖（略）

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

�。�5）五年預(yù)測(cè)值見(jiàn) 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）繪制時(shí)序圖（略）

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

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

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

　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）繪制時(shí)序圖（略）

�。�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ù)顯著非零，殘差為檢驗(yàn)為白噪聲序列，MA(1)模型 SBC 更小一點(diǎn)，所以選擇 MA(1)模型作為最優(yōu)模型。

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

　 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）未來(lái) 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）未來(lái) 5 年的谷物產(chǎn)量預(yù)測(cè)：

　 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）未來(lái) 5 年的預(yù)測(cè)

　 第 五 章習(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 階截尾特征，對(duì)該序列擬合 ARIMA(1,1,0)模型，模型參數(shù)如下：

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

　 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 檢驗(yàn)平穩(wěn)，根據(jù)自相關(guān)圖和偏自相關(guān)圖特征，可以識(shí)別為偏自相關(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ù)時(shí)序圖和 adf 檢驗(yàn)結(jié)果可知，原序列非平穩(wěn)，一階差分后序列平穩(wěn)。

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

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

　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）我國(guó)人口出生率、死亡率和自然增長(zhǎng)率序列均有顯著線性趨勢(shì)，為典型非平穩(wěn)序列。

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

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

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

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

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

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

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

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

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

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

　對(duì)玉米價(jià)格序列進(jìn)行 10 期預(yù)測(cè)，預(yù)測(cè)值如下：

　對(duì)玉米價(jià)格序列的擬合與預(yù)測(cè)效果圖如下：

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

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

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

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

　對(duì)工人薪水序列進(jìn)行 10 期預(yù)測(cè)，預(yù)測(cè)值如下：

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

　④ 對(duì)生豬價(jià)格序列擬合模型 ARIMA(3,1,0)，參數(shù)估計(jì)值如下：

　 對(duì)生豬價(jià)格序列進(jìn)行 10 期預(yù)測(cè)，預(yù)測(cè)值如下：

　對(duì)生豬價(jià)格序列的擬合與預(yù)測(cè)效果圖如下：

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

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

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

　 第六章習(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）該序列受到趨勢(shì)、季節(jié)、隨機(jī)波動(dòng)等三個(gè)因素的影響。由于季節(jié)性沒(méi)有隨趨勢(shì)變化而發(fā)生顯著變化，可以選擇加法模型分解各因素。

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

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

�。�4）使用 Holt-Winters 加法模型擬合該序列的發(fā)展，幵做 2 年期預(yù)測(cè)。相關(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ù)測(cè)值如下：

　 預(yù)測(cè)效果圖：

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

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

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

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

　由于因素分解方法和 Holt-Winters 三參數(shù)指數(shù)平滑法丌提供極大似然函數(shù)值，所以無(wú)法計(jì)算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é)、趨勢(shì)和隨機(jī)波動(dòng)的影響�？梢詳M合加法因素分解模型。

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

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

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

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

　Coefficients:

　 ar1

　 sma1

　 -0.3339

　-0.5960 s.e.

　 0.1224

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

　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 中心移動(dòng)平均做預(yù)測(cè),則 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ù)測(cè)值 ˆ (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）式等號(hào)兩邊同乘 (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 年我國(guó)人口總數(shù)序列為顯著的線性遞增序列，可以使用 holt 兩參數(shù)指數(shù)平滑法進(jìn)行趨勢(shì)擬合和預(yù)測(cè)，或使用 ARIMA（1,1,0）模型進(jìn)行擬合和預(yù)測(cè)。

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

　Coefficients:

　ar1

　 0.9473 s.e.

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

　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) 該序列時(shí)序圖顯示典型的非線性趨勢(shì)，可以用 Holt-Winters 指數(shù)平滑模型擬合也可以用 2階差分的 ARIMA 模型擬合。時(shí)序圖如下

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

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

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

　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) 時(shí)序圖顯示該序列時(shí)序圖有季節(jié)效應(yīng)、趨勢(shì)（循環(huán)）效應(yīng)和隨機(jī)波動(dòng)。

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

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

　方法一：使用 Holt-Winters 模型對(duì)該序列進(jìn)行預(yù)測(cè)，相關(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ù)測(cè)效果圖

　方法二：使用 ARIMA(2,1,0) ? (0,1,1) 12 對(duì)該序列進(jìn)行擬合與預(yù)測(cè)，相關(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ù)測(cè)效果圖

　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) 考察該序列時(shí)序圖，時(shí)序圖顯示該序列具有顯著的線性遞增趨勢(shì)和年度季節(jié)效應(yīng)。

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

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

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

　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ù)測(cè)效果圖如下

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

　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ù)測(cè)效果圖如下

　 第七章習(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)兩個(gè)序列皆平穩(wěn)。

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

　(3)因?yàn)檫@兩個(gè)變量都具有白噪聲屬性，所以不考慮建立協(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)出顯著的周期特征，兩個(gè)序列均為非平穩(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ù)測(cè)值為

　 預(yù)測(cè)效果圖如下

　掠食者的擬合模型為：

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

　預(yù)測(cè)效果圖如下

　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) 進(jìn)口總額和出口總額序列均丌平穩(wěn)，均為 2 階單整序列

　(2)

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

　2t tx ? ? ?

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

　(3) 這兩個(gè)序列具有協(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）消費(fèi)品零售總額序列擬合 ARIMA((3),2,0)模型，5 年預(yù)測(cè)值為：

　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 國(guó)內(nèi)生產(chǎn)總值序列擬合 ARIMA((2),2,0),5 年預(yù)測(cè)值為：

　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）這兩個(gè)序列之間不存在協(xié)整關(guān)系。

　（3）沒(méi)有協(xié)整方程,所以無(wú)法寫(xiě)出協(xié)整方程 （4）因?yàn)檫@兩個(gè)變量不存在協(xié)整關(guān)系，所以無(wú)法直接分析如果中國(guó)國(guó)內(nèi)社會(huì)消費(fèi)品零售總額增長(zhǎng) 1%,對(duì)國(guó)內(nèi)生產(chǎn)總值有什么影響。需要引入其他變量，讓整個(gè)系統(tǒng)構(gòu)成協(xié)整關(guān)系之后，才能準(zhǔn)確分析這個(gè)問(wè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) 安全帶強(qiáng)制法律的執(zhí)行對(duì)司機(jī)傷亡數(shù)據(jù)有顯著的干預(yù)作用。

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

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

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

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

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

　(7) 司機(jī)傷亡數(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)這五個(gè)變量均為 1 階單整序列 (2)農(nóng)場(chǎng)工人薪水為果，玉米價(jià)格、玉米產(chǎn)量、生豬價(jià)格、生豬產(chǎn)量為因。

　(3)分析農(nóng)場(chǎng)工人平均工資、農(nóng)場(chǎng)作物、家畜的價(jià)格不供應(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?? ? ? ? ? 工人薪水 玉米價(jià)格 生豬價(jià)格 生豬產(chǎn)量3=4.0423+0.0486 +0.144 +0.243+1 1.19 0.21t t t ttB B?? ?工人薪水 玉米價(jià)格 生豬價(jià)格 生豬產(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) 分別繪制這四個(gè)序列的時(shí)序圖（圖略）,這四個(gè)序列都有趨勢(shì)。它們均為 1 階單整序列。分別擬合單變量 ARIMA 模型如下：

　Logm1 序列：

　 Loggnp 序列：

　 短期利率序列：

　 長(zhǎng)期利率序列：

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

　(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í)題 序列 人大