时间序列：ARIMA模型

1、实验：建立ARIMA模型（综合性实验） 实验题目：某城市连续14年的月度婴儿出生率数据如下表所示：26.66323.59826.93124.74025.80624.36424.47723.90123.17523.22721.67221.87021.43921.08923.70921.66921.75220.76123.47923.82423.10523.11021.75922.07321.93720.03523.59021.67222.22222.12323.95023.50422.23823.14221.05921.57321.54820.00022.42420.61521.76122.87

2、424.10423.74823.26222.90721.51922.02522.60420.89424.67723.67325.32023.58324.67124.45424.12224.25222.08422.99123.28723.04925.07624.03724.43024.66726.45125.61825.01425.11022.96423.98123.79822.27024.77522.64623.98824.73726.27625.81625.21025.19923.16224.70724.36422.64425.56524.06225.43124.63527.00926.60

3、626.26826.46225.24625.18024.65723.30426.98226.19927.21026.12226.70626.87826.15226.37924.71225.68824.99024.23926.72123.47524.76726.21928.36128.59927.91427.78425.69326.88126.21724.21827.91426.97528.52727.13928.98228.16928.05629.13626.29126.98726.58924.84827.54326.89628.87827.39028.06528.14129.04828.48

4、426.63427.73527.13224.92428.96326.58927.93128.00929.22928.75928.40527.94525.91226.61926.07625.28627.66025.95126.39825.56528.86530.00029.26129.01226.99227.897（1）选择适当模型拟和该序列的发展（2）使用拟合模型预测下一年度该城市月度婴儿出生率实验内容：给出实际问题的非平稳时间序列，要求学生利用R统计软件，对该序列进行分析，通过平稳性检验、差分运算、白噪声检验、拟合ARMA模型，建立ARIMA模型，在此基础上进行预测。实验要求：处理数据，掌握

5、非平稳时间序列的ARIMA建模方法，并根据具体的实验题目要求完成实验报告，并及时上传到给定的FTP和课程网站。实验步骤：第一步：编程建立R数据集；第二步：调用plot.ts程序对数据绘制时序图。第三步：从时序图中利用平稳时间序列的定义判断是否平稳？ 第四步：若不满足平稳性，则可利用差分运算是否能使序列平稳？重复第三步步骤第五步：根据Box.test纯随机检验结果，利用LB统计量和白噪声特性检验最后处理的时间序列是否为纯随机序列？第六步：在序列判断为平稳非白噪声序列后，求出该观察值序列的样本自相关系数（ACF）和样本偏自相关系数（PACF）的值，选择阶数适当的ARIMA（p,d,q）模型进行拟合

6、，并估计模型中未知参数的值。第七步：检验模型的有效性。如果拟合模型通不过检验，转向步骤6，重新选择模型再拟合。第八步：模型优化。如果拟合模型通过检验，仍然转向步骤6，充分考虑各种可能建立多个拟合模型，从所有通过检验的拟合模型中选择最优模型。第九步：利用最优拟合模型，预测下一年度该城市月度婴儿出生率。ex5.2=ts(scan(ex5.2.txt), frequency=4)Read 168 itemsplot.ts(ex5.2)从图中看出序列一开始有下降趋势，后面有明显上升趋势，所以序列不平稳。d12ex5.2 = diff(ex5.2,lag=12)acf(d12ex5.2,48)plot(

7、d12ex5.2)从上面的自相关图中可以看出改做滞后12期差分后为平稳。Box.test(d12ex5.2, lag=17, type=Ljung-Box) Box-Ljung testdata: d12ex5.2 X-squared = 147.9254, df = 17, p-value 2.2e-16P值小于0.05，可以认为是非白噪声序列。par(mfrow=c(2,1); acf(d12ex5.2, 48); pacf(d12ex5.2, 48)ARIMA（0,0,3）、ARIMA（0,0,4）、ARIMA（1,0,3）、ARIMA（1,0,4）四个模型分别进行拟合检验(rec.ol

8、s = arima(d12ex5.2,order=c(0,0,3)Call:arima(x = d12ex5.2, order = c(0, 0, 3)Coefficients: ma1 ma2 ma3 intercept 0.7949 0.4480 0.1156 0.2150s.e. 0.0839 0.0832 0.0885 0.1744sigma2 estimated as 0.8621: log likelihood = -210.12, aic = 430.25rec.pr = predict(rec.ols, n.ahead=5)U = rec.pr$pred + 1.96*rec.

9、pr$seL = rec.pr$pred - 1.96*rec.pr$seminx = min(d12ex5.2,L)maxx = max(d12ex5.2,U)ts.plot(d12ex5.2, rec.pr$pred, ylim=c(minx,maxx)lines(rec.pr$pred, col=red, type=o)lines(U, col=blue, lty=dashed)lines(L, col=blue, lty=dashed)qqnorm(rec.ols$resid)qqline(rec.ols$resid)shapiro.test(rec.ols$resid) Shapir

10、o-Wilk normality testdata: rec.ols$resid W = 0.9777, p-value = 0.0125用shapiro检验，发现p值为0.0125，在5%的显著性水平下显著，所以为ARIMA（0,0,3）模型不合理。(rec.ols = arima(d12ex5.2,order=c(0,0,4)Call:arima(x = d12ex5.2, order = c(0, 0, 4)Coefficients: ma1 ma2 ma3 ma4 intercept 0.8306 0.4943 0.2254 0.2070 0.2041s.e. 0.0902 0.115

11、8 0.0925 0.0889 0.1994sigma2 estimated as 0.828: log likelihood = -207.07, aic = 426.15rec.pr = predict(rec.ols, n.ahead=5)U = rec.pr$pred + 1.96*rec.pr$seL = rec.pr$pred - 1.96*rec.pr$seminx = min(d12ex5.2,L)maxx = max(d12ex5.2,U)ts.plot(d12ex5.2, rec.pr$pred, ylim=c(minx,maxx)lines(rec.pr$pred, co

12、l=red, type=o)lines(U, col=blue, lty=dashed)lines(L, col=blue, lty=dashed)qqnorm(rec.ols$resid)qqline(rec.ols$resid)shapiro.test(rec.ols$resid) Shapiro-Wilk normality testdata: rec.ols$resid W = 0.9689, p-value = 0.001363用shapiro检验，发现p值为0.001363，在5%的显著性水平下显著，所以为ARIMA（0,0,4）模型不合理。(rec.ols = arima(d12

13、ex5.2,order=c(1,0,3)Call:arima(x = d12ex5.2, order = c(1, 0, 3)Coefficients: ar1 ma1 ma2 ma3 intercept 0.9288 -0.1369 -0.2156 -0.1586 0.0240s.e. 0.0669 0.1065 0.0921 0.0879 0.4984sigma2 estimated as 0.7986: log likelihood = -204.35, aic = 420.7rec.pr = predict(rec.ols, n.ahead=5)U = rec.pr$pred + 1.

14、96*rec.pr$seL = rec.pr$pred - 1.96*rec.pr$seminx = min(d12ex5.2,L)maxx = max(d12ex5.2,U)ts.plot(d12ex5.2, rec.pr$pred, ylim=c(minx,maxx)lines(rec.pr$pred, col=red, type=o)lines(U, col=blue, lty=dashed)lines(L, col=blue, lty=dashed)qqnorm(rec.ols$resid)qqline(rec.ols$resid)shapiro.test(rec.ols$resid)

15、 Shapiro-Wilk normality testdata: rec.ols$resid W = 0.9783, p-value = 0.01454(rec.ols = arima(d12ex5.2,order=c(1,0,4)Call:arima(x = d12ex5.2, order = c(1, 0, 4)Coefficients: ar1 ma1 ma2 ma3 ma4 intercept 0.9084 -0.1288 -0.2457 -0.1511 0.1309 0.0493s.e. 0.0725 0.1078 0.1021 0.0778 0.0960 0.4684sigma2 estimated as 0.7891: log likelihood = -203.47, aic = 420.94rec.pr = predict(rec.ols, n.ahead=5)U = rec.pr$pred + 1.96*rec.pr$seL = rec.pr$pred - 1.96*rec.pr$seminx = min(d12ex5.2,L)maxx = max(d12ex5.2,U)ts.plot(d12ex5.2, rec.pr$pred, ylim=c(m