時(shí)間序列：ARIMA模型

1、實(shí)驗(yàn)：建立ARIMA模型（綜合性實(shí)驗(yàn)） 實(shí)驗(yàn)題目：某城市連續(xù)14年的月度嬰兒出生率數(shù)據(jù)如下表所示：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）選擇適當(dāng)模型擬和該序列的發(fā)展（2）使用擬合模型預(yù)測(cè)下一年度該城市月度嬰兒出生率實(shí)驗(yàn)內(nèi)容：給出實(shí)際問(wèn)題的非平穩(wěn)時(shí)間序列，要求學(xué)生利用R統(tǒng)計(jì)軟件，對(duì)該序列進(jìn)行分析，通過(guò)平穩(wěn)性檢驗(yàn)、差分運(yùn)算、白噪聲檢驗(yàn)、擬合ARMA模型，建立ARIMA模型，在此基礎(chǔ)上進(jìn)行預(yù)測(cè)。實(shí)驗(yàn)要求：處理數(shù)據(jù)，掌握

5、非平穩(wěn)時(shí)間序列的ARIMA建模方法，并根據(jù)具體的實(shí)驗(yàn)題目要求完成實(shí)驗(yàn)報(bào)告，并及時(shí)上傳到給定的FTP和課程網(wǎng)站。實(shí)驗(yàn)步驟：第一步：編程建立R數(shù)據(jù)集；第二步：調(diào)用plot.ts程序?qū)?shù)據(jù)繪制時(shí)序圖。第三步：從時(shí)序圖中利用平穩(wěn)時(shí)間序列的定義判斷是否平穩(wěn)？ 第四步：若不滿(mǎn)足平穩(wěn)性，則可利用差分運(yùn)算是否能使序列平穩(wěn)？重復(fù)第三步步驟第五步：根據(jù)Box.test純隨機(jī)檢驗(yàn)結(jié)果，利用LB統(tǒng)計(jì)量和白噪聲特性檢驗(yàn)最后處理的時(shí)間序列是否為純隨機(jī)序列？第六步：在序列判斷為平穩(wěn)非白噪聲序列后，求出該觀察值序列的樣本自相關(guān)系數(shù)（ACF）和樣本偏自相關(guān)系數(shù)（PACF）的值，選擇階數(shù)適當(dāng)?shù)腁RIMA（p,d,q）模型進(jìn)行擬合

6、，并估計(jì)模型中未知參數(shù)的值。第七步：檢驗(yàn)?zāi)Ｐ偷挠行?。如果擬合模型通不過(guò)檢驗(yàn)，轉(zhuǎn)向步驟6，重新選擇模型再擬合。第八步：模型優(yōu)化。如果擬合模型通過(guò)檢驗(yàn)，仍然轉(zhuǎn)向步驟6，充分考慮各種可能建立多個(gè)擬合模型，從所有通過(guò)檢驗(yàn)的擬合模型中選擇最優(yōu)模型。第九步：利用最優(yōu)擬合模型，預(yù)測(cè)下一年度該城市月度嬰兒出生率。ex5.2=ts(scan(ex5.2.txt), frequency=4)Read 168 itemsplot.ts(ex5.2)從圖中看出序列一開(kāi)始有下降趨勢(shì)，后面有明顯上升趨勢(shì)，所以序列不平穩(wěn)。d12ex5.2 = diff(ex5.2,lag=12)acf(d12ex5.2,48)plot(

7、d12ex5.2)從上面的自相關(guān)圖中可以看出改做滯后12期差分后為平穩(wěn)。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，可以認(rèn)為是非白噪聲序列。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）四個(gè)模型分別進(jìn)行擬合檢驗(yàn)(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檢驗(yàn)，發(fā)現(xiàn)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檢驗(yàn)，發(fā)現(xiàn)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