人大版时间序列分析基于R（第2版）习题答案

　第一章习题答案 略

　第二章习题答案 2.1

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

　答案：

　（1）非平稳，有典型线性趋势

　（2）延迟 1-6 阶自相关系数如下：

　（3）典型的具有单调趋势的时间序列样本自相关图

　2.2 R 命令 #先读入数据文件 co2<-ts(E2_2$co2,start=c(1975,1),frequency = 12) plot(co2,type="o") acf(co2)$acf

　（1）非平稳，时序图如下

　（2）1-24 阶自相关系数如下 [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）自相关图呈现典型的长期趋势与周期并存的特征

　 2.3

　 R R 命令

　#先读入数据文件 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 阶自相关系数 [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）平稳序列

　 （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) 答案：

　我们自定义函数，计算该序列各阶延迟的 Q 统计量及相应 P 值。由于延迟 1-12 阶 Q 统计量的 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）绘制时序图与自相关图

　 （2）序列时序图显示出典型的周期特征，该序列非平稳 （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）如果是进行平稳性图识别，该序列自相关图呈现一定的趋势序列特征，可以视为非平稳非白噪声序列。

　如果通过 adf 检验进行序列平稳性识别，该序列带漂移项的 0 阶滞后 P 值小于 0.05，可以视为平稳非白噪声序列

　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）差分序列平稳，非白噪声序列

　 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）时序图和自相关图显示该序列有趋势特征，所以图识别为非平稳序列。

　 （2）单位根检验显示带漂移项 0 阶延迟的 P 值小于 0.05，所以基于 adf 检验可以认为该序列平稳

　（3）如果使用 adf 检验结果，认为该序列平稳，则白噪声检验显示该序列为非白噪声序列 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 如果使用图识别认为该序列非平稳，那么一阶差分后序列为平稳非白噪声序列

　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）时序图和自相关图都显示典型的趋势序列特征

　（2）单位根检验显示该序列可以认为是平稳序列（带漂移项一阶滞后 P 值小于 0.05）

　 （3）一阶差分后序列平稳

　 第三章习题答案 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             （）

　（ ）模型的平稳条件是 11 21,2 1, 2k k kcc k     （ ）

　 3.5

　证明:

　该序列的特征方程为：3 20 c c        ，解该特征方程得三个特征根：

　11   ，2c   ，3c  

　无论 c 取什么值，该方程都有一个特征根在单位圆上，所以该序列一定是非平稳序列。证毕。

　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          根据待定系数法：

　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）证明：因为对任意常数 C，有 2 2( ) lim(1 )tkVar x kC   

　所以该序列为非平稳序列。

　（2）1 1( 1)t t t t ty x x C        ，则序列  ty 满足如下条件：

　均值、方差为常数， ( ) 0tE y  ,2 2( ) 1 ( 1)tVar y C       自相关系数只与时间间隔长度有关，与起始时间无关

　121, 0, 21 ( 1)kCkC     所以该差分序列为平稳序列。

　3.11

　（1）非平稳，（2）平稳，（3）可逆，（4）不可逆，（5）平稳可逆，（6）不平稳不可逆

　3.12

　该模型的 Green 函数为：

　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  ，根据 (1,1) ARMA 模型 Green 函数的递推公式得：

　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

　 该习题数据文件与 2.7 相同。该题问题设置有问题：是要问如果判断该序列或差分序列是平稳序列，那该平稳序列具有 ARMA 族中哪个模型的特征。

　（1）根据 adf 检验结果可以认为该序列平稳。根据序列的自相关图可以认为是自相关系数拖尾。根据偏自相关图，可以认为是偏自相关 2 阶截尾，所以该序列具有 AR（2）模型的特征。

　 （2）根据图识别也可以认为该序列不平稳，对该序列进行一阶差分。一阶差分后序列可以视为平稳序列，根据差分后序列的自相关图可以认为是自相关系数 1 阶截尾，具有 MA（1）模型的特征。根据偏自相关图，可以认为是偏自相关 3 阶截尾，具有 AR（3）模型的特征。具体哪个模型最适合拟合该序列，下一章介绍。

　3.17

　该习题数据文件与 2.8 相同。该题问题设置有问题：是要问如果判断该序列或差分序列是平

　稳序列，那该平稳序列具有 ARMA 族中哪个模型的特征。

　（1）根据 adf 检验，该序列可以视为平稳序列。自相关图呈现拖尾属性，偏自相关图呈现1 阶截尾特征，所以该序列呈现出 AR(1)模型特征。

　 （3）如果根据图识别，可以认为序列蕴含趋势，视为非平稳序列。一阶差分后序列平稳。一阶差分后序列呈现自相关系数 2 阶截尾，偏自相关系数 2 阶截尾的特征。所以一阶差分后序列具有 MA(2)模型特征，或 AR（2）模型特征。具体哪个模型最适合拟合该序列，下一章介绍。

　第四章习题答案 4.1 本题 R 代码 读入数据文件之后 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）该序列为平稳非白噪声 （3）自相关图拖尾，偏自相关图一阶截尾 （4）拟合 AR(1)模型

　（5）五年预测值见预测输出（略）

　4.2 本题 R 代码 读入数据文件之后 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）该序列为平稳非白噪声序列 （3）自相关图一阶截尾，偏自相关图拖尾 （4）拟合 MA(1)模型

　（5）五年预测值见 R 输出（略）

　4.3 本题 R 代码 读入数据文件之后 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）该序列为平稳非白噪声序列 （3）根据该序列自相关图，可以视为：自相关图一阶截尾，或偏自相关 2 阶截尾 （4）分别拟合 MA(1)模型和 AR(2)模型，两个模型均参数显著非零，残差为检验为白噪声序列，AIC 和 SBC 的结果几乎相等，最后考虑白噪声检验的 P 值，AR(2)模型的白噪声检验 P 值

　更大，说明该模型对序列的相关信息提取更为充分，所以选择 AR(2)模型作为最优模型。

　（5）基于 AR(2)模型未来一年预测值为

　4.4 本题 R 代码 读入数据文件之后 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）该序列为平稳非白噪声序列 （3）根据该序列自相关图，可以视为：自相关图一阶截尾，或偏自相关 2 阶截尾，或自相关和偏自相关均拖尾 （4）分别拟合 MA(1)模型，AR(2)模型和 ARMA(1,1,) 。ARMA(1,1)模型参数不能拒绝参数为零的原假设，所以淘汰。MA(1)模型，AR(2)模型均参数显著非零，残差为检验为白噪声序列，MA(1)模型 SBC 更小一点，所以选择 MA(1)模型作为最优模型。

　 （5）基于 MA(1)模型未来五年预测值为

　 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   的 的置信区间等于 ，即( ）

　（2）更新数据后 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   的 的置信区间等于 ，即( ）

　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）平稳非白噪声序列 （2）自相关系数拖尾，偏自相关系数 1 阶截尾，拟合 AR（1）模型

　（3）未来 5 年的降雪量

　 4.7 本题 R 指令 读入数据文件之后 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）平稳非白噪声序列 （2）自相关系数和偏自相关系数都拖尾，拟合 ARMA(1,1)模型

　（3）未来 5 年的谷物产量预测：

　 4.8 本题 R 指令 读入数据文件之后 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）平稳非白噪声序列 （2）自相关系数 1 阶截尾，偏自相关系数拖尾，拟合 MA(1)模型

　 （3）下一刻的 95%置信区间为：（80.43481 90.92585）

　4.9 本题 R 指令 读入数据文件之后 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）平稳非白噪声序列 （2）自相关系数和偏自相关系数都拖尾，拟合 ARMA(4,1)模型（该模型有部分系数不能显著非零）。所以最好是拟合疏系数 ARMA((3,4),1) 模型

　（3）未来 5 年的预测

　 第 五 章习题答案 5.1 本题 R 指令 读入数据文件之后 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 阶差分后平稳，根据差分后序列的自相关图拖尾和偏自相关图 1 阶截尾特征，对该序列拟合 ARIMA(1,1,0)模型，模型参数如下：

　根据该模型，得到 2009-2013 年我国铁路货运量的预测值为：

　 5.2 本题 R 指令 读入数据文件之后 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 检验平稳，根据自相关图和偏自相关图特征，可以识别为偏自相关图 1 阶截尾，拟合 AR（1）模型。参数输出结果如下：

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

　 5.3 本题 R 指令 读入数据文件之后 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）根据时序图和 adf 检验结果可知，原序列非平稳，一阶差分后序列平稳。

　（2）一阶差分后序列自相关图拖尾，偏自相关图 3 阶截尾，所以对原序列拟合模型ARIMA(3,1,0)，参数输出结果显示 2 阶自相关系数与均值都不显著非零。所以最后拟合不带均值项的疏系数模型 ARIMA((1,3),1,0)，最后参数估计值输出结果如下

　 （3）基于模型对 1939—1945 年英国绵羊的数量的预测为：

　5_4 本题 R 指令 读入数据文件之后 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）我国人口出生率、死亡率和自然增长率序列均有显著线性趋势，为典型非平稳序列。

　（2）根据图检验和 adf 检验结果显示，我国人口出生率、死亡率和自然增长率序列的一阶差分后序列均平稳。

　（3）我国人口出生率可以拟合 ARIMA(0,1,0)模型，也可以拟合 ARIMA(0,1,1)模型。其中ARIMA（0,1,1）模型 AIC,SBC 更小一点。得到的拟合模型参数为

　基于 ARIMA（0,1,1）模型，预测 1998 年之后人口出生率未来十年的预测值为：

　 （4）我国人口死亡率可以拟合随机游走模型 ARIMA(0,1,0)。基于 ARIMA（0,1,0）模型，预测 1998 年之后人口死亡率未来十年的预测值为：

　 （5）（5）我国人口自然增长率可以拟合 ARIMA(0,1,0)模型，ARIMA(1,1,0)模型，也可以拟合 ARIMA(0,1,1)模型。其中 ARIMA（0,1,1）模型 AIC,SBC 更小一点。得到的拟合模型参数为

　基于 ARIMA（0,1,1）模型，预测 1998 年之后人口自然增长率未来十年的预测值为：

　 5_5 本题 R 指令 读入数据文件之后

　答案 （1）玉米价格、玉米产量、工人薪水、生猪价格、生猪产量序列都不平稳， （2）玉米价格、玉米产量、工人薪水、生猪价格、生猪产量序列一阶差分都实现平稳。

　（3）-（4）分序列分析结果：

　① 对玉米价格序列拟合疏系数模型 ARIMA(0,1,(7))，参数估计值如下：

　对玉米价格序列进行 10 期预测，预测值如下：

　对玉米价格序列的拟合与预测效果图如下：

　 ② 对玉米产量序列拟合模型 ARIMA(0,1,1)，参数估计值如下：

　 对玉米产量序列进行 10 期预测，预测值如下：

　对玉米产量列的拟合与预测效果图如下：

　③ 对工人薪水序列拟合疏系数模型 ARIMA(（1,13）,1,0)，参数估计值如下：

　对工人薪水序列进行 10 期预测，预测值如下：

　对工人薪水序列的拟合与预测效果图如下：

　④ 对生猪价格序列拟合模型 ARIMA(3,1,0)，参数估计值如下：

　 对生猪价格序列进行 10 期预测，预测值如下：

　对生猪价格序列的拟合与预测效果图如下：

　⑤ 对生猪产量序列拟合疏系数模型 ARIMA((2,3,5),1,0)，参数估计值如下：

　 对生猪产量序列进行 10 期预测，预测值如下：

　对生猪产量序列的拟合与预测效果图如下：

　 第六章习题答案 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）该序列受到趋势、季节、随机波动等三个因素的影响。由于季节性没有随趋势变化而发生显著变化，可以选择加法模型分解各因素。

　（2）提取该序列的趋势效应（趋势数据见 fit$trend，具体数据略）， 趋势效应图如下：

　（3）提取该序列的季节效应（季节效应数据见 fit$figure，具体数据略））。季节效应图如下：

　（4）使用 Holt-Winters 加法模型拟合该序列的发展，幵做 2 年期预测。相关参数如下：

　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 年预测值如下：

　 预测效果图：

　 （5）用 ARIMA 季节模型拟合幵预测该序列的发展。

　拟合模型：ARIMA(1,1,0)×(0,1,1) 12, 模型参数估计输出如下：

　拟合不预测效果图：

　 (6)比较分析上面使用过的三种模型的拟合精度。

　由于因素分解方法和 Holt-Winters 三参数指数平滑法丌提供极大似然函数值，所以无法计算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 ）该序列受到季节、趋势和随机波动的影响。可以拟合加法因素分解模型。

　（2）提取该序列的趋势效应（数据略）， 趋势效应图如下

　（3）提取该序列的季节效应（数据略），季节效应图如下

　 （4）用指数平滑法对该序列做 2 年期预测（预测值略，预测效果图如下）

　 （5）用 ARIMA 季节模型拟合幵预测该序列的发展。拟合模型：ARIMA(1,1,0)×(0,1,1) 12, 模型参数估计输出如下：

　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 中心移动平均做预测,则 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 期预测值 ˆ (2)tx中3 tx前面的系数是532，4 tx前面的系数是164。

　6_4 解下面的方程组，得到 0.4  

　5.25 5(1 )5.26 5.5 (1 )tt          6_5 根据指数平滑的定义有（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 年我国人口总数序列为显著的线性递增序列，可以使用 holt 两参数指数平滑法进行趋势拟合和预测，或使用 ARIMA（1,1,0）模型进行拟合和预测。

　(2) ARIMA（1,1,0）模型拟合效果最优的模型相对最优,该模型参数估计结果如下， Call: arima(x = x, order = c(1, 1, 0))

　Coefficients:

　ar1

　 0.9473 s.e.

　0.0345 使用该模型对 2009—2016 年中国人口总数进行预测（预测值略），预测效果图如下

　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 指数平滑模型拟合也可以用 2阶差分的 ARIMA 模型拟合。时序图如下

　(2) 拟合效果是疏系数 ARIMA(0,

　2，(1, 16))相对较优。该模型的参数估计情况为

　（3）5 年的预测（预测值略），拟合与预测效果图如下：

　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) 时序图显示该序列时序图有季节效应、趋势（循环）效应和随机波动。

　 (2) 使用因素分解加法模型对该序列进行因素分解

　 (3) 可以选择 Holt-Winters 模型对该序列进行拟合与预测，也可以选择 ARIMA 模型对该序列进行拟合与预测。

　方法一：使用 Holt-Winters 模型对该序列进行预测，相关参数如下

　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

　预测效果图

　方法二：使用 ARIMA(2,1,0)  (0,1,1) 12 对该序列进行拟合与预测，相关参数如下 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

　预测效果图

　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) 考察该序列时序图，时序图显示该序列具有显著的线性递增趋势和年度季节效应。

　(2)根据该序列呈现的规律,可以使用 Holt-Winters 模型对该序列进行拟合与预测，也可以选择 ARIMA 模型对该序列进行拟合与预测。

　(3)选择拟合效果相对最好的模型 ,预测该序列未来 3 年的旅馆入住情况（没有唯一解，给出可能参考答案）

　方法一：可以基于 Holt-Winters 模型对该序列进行拟合与预测 参数估计如下：

　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 拟合与预测效果图如下

　方法二：基于 ARIMA(1,1,1)  (0,1,1) 12 模型对该序列进行拟合与预测，参数估计如下：

　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 拟合与预测效果图如下

　 第七章习题答案 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)两个序列皆平稳。

　(2)两个序列皆为白噪声序列。

　(3)因为这两个变量都具有白噪声属性，所以不考虑建立协整方程。

　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）掠食者和被掠食者数量都呈现出显著的周期特征，两个序列均为非平稳序列。但是掠食者和被掠食者延迟2阶序列具有协整关系。协整方程参数如下：

　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  一周的预测值为

　 预测效果图如下

　掠食者的拟合模型为：

　22166.9407 0.11431 1.1892 0.5688t t ty xB B    一周的预测值为：

　预测效果图如下

　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) 进口总额和出口总额序列均丌平稳，均为 2 阶单整序列

　(2)

　出口总额序列拟合模型 ARIMA(0,2,0)：

　2t tx   

　进口总额序列拟合模型：ARIMA((2),2,0):

　(3) 这两个序列具有协整关系

　 (4)协整模型如下：

　 （5）构造该协整模型的误差修正模型：

　 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 年预测值为：

　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 国内生产总值序列拟合 ARIMA((2),2,0),5 年预测值为：

　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）这两个序列之间不存在协整关系。

　（3）没有协整方程,所以无法写出协整方程 （4）因为这两个变量不存在协整关系，所以无法直接分析如果中国国内社会消费品零售总额增长 1%,对国内生产总值有什么影响。需要引入其他变量，让整个系统构成协整关系之后，才能准确分析这个问题。

　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) 安全带强制法律的执行对司机伤亡数据有显著的干预作用。

　(2) 司机伤亡数据不汽油价格及安全带强制法律的执行之间具有协整关系，具体协整方程如下：

　(3) 安全带强制法律的执行对前座乘客伤亡数据有显著的干预作用。

　(4) 前座乘客伤亡数据汽油价格及安全带强制法律的执行之间具有协整关系，具体协整方程如下

　(5)安全带强制法律的执行对后座乘客伤亡数据没有显著的干预作用。

　(6)后座乘客伤亡数据不行驶里程、汽油价格之间具有协整关系，具体协整方程如下：

　(7) 司机伤亡数据、前座乘客伤亡数据和后座乘客伤亡数据之间具有协整关系

　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)农场工人薪水为果，玉米价格、玉米产量、生猪价格、生猪产量为因。

　(3)分析农场工人平均工资、农场作物、家畜的价格不供应链之间具有协整关系。拟合协整模型如下（参数含义解释略）

　误差修正模型如下（参数含义解释略）

　 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     工人薪水 玉米价格 生猪价格 生猪产量3=4.0423+0.0486 +0.144 +0.243+1 1.19 0.21t t t ttB B 工人薪水 玉米价格 生猪价格 生猪产量

　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 因果关系：loggnp 为果，其他三个序列为因。

　(3)协整模型为

　 （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  