ARMA模型擬合指令：

arima(x,order=c(p,d,q),method=c("ML","CSS"),include.mean)

其中:

• x—帶估計序列；
• Order—指定模型階數(shù)；其中，P為自回歸階數(shù)；d為差分階數(shù)，若為平穩(wěn)序列，則不需要差分，d=0；q為移動平均階數(shù)。
• method—估計方法。method =CSS-ML，系統(tǒng)默認(rèn)的是條件最小二乘估計和極大似然估計的混合方法；method =ML，極大似然估計；method =CSS，條件最小二乘估計；
• include.mean—是否包含均值項，即估計結(jié)果中的intercept項。系統(tǒng)默認(rèn)含有常數(shù)項。

案例：選擇合適的模型擬合1950-2008年我國油路及農(nóng)村投遞線路每年新增里程數(shù)序列。
（1）讀取數(shù)據(jù)

file.choose #彈出對話框，可以讓我們選擇文件的位置。R的運行結(jié)果： #"C:\\Users\\w\\Documents\\ licheng.txt"d=read.table("C:\\Users\\w\\Documents\\licheng.txt") #讀入數(shù)據(jù)。 # 注意：讀取txt格式文件時用“read.table”指令，header默認(rèn)為FALSE，也即認(rèn)為第一行就是數(shù)據(jù)；所以如果第一行不是數(shù)據(jù)而是標(biāo)簽，要改寫為header=T/TURE；讀取csv格式文件時用“read.csv”，header默認(rèn)為Ture，也即認(rèn)為第一行是標(biāo)題行；licheng=ts(d,start=1950,end=2008,freq=1) #轉(zhuǎn)化為時間序列數(shù)據(jù)。年度數(shù)據(jù)則freq=1; 半年度數(shù)據(jù)則freq=2; 季度數(shù)據(jù)則freq=4; 月度數(shù)據(jù)則freq=12; 周度數(shù)據(jù)則freq=52; 日度數(shù)據(jù)則freq=365;

（2）序列的平穩(wěn)性檢驗：
時序圖，自相關(guān)圖，ADF檢驗

plot(licheng) #作時序圖。大致圍繞0值上下波動，且波動有界，該序列大致為平穩(wěn)時間序列。

圖1 時序圖

acf(licheng,14) #作自相關(guān)圖，自相關(guān)系數(shù)具有短期相關(guān)特性，因此進一步認(rèn)為該序列為平穩(wěn)時間序列。且具有拖尾特性。 pacf(licheng,14) #作偏自相關(guān)圖，2階截尾

圖2 自相關(guān)圖

圖3 偏自相關(guān)圖

單位根檢驗的做法：

library(fBasics) library(fUnitRoots) #下載并安裝“fUnitRoots”包，并用library函數(shù)進行調(diào)用。該程序包中的adfTest函數(shù)可以很便捷的進行單位根檢驗。for(i in 1:3)print(adfTest(licheng,lag=i,type="c"))

對“l(fā)icheng”序列進行檢驗類型為類型2（無趨勢，有常數(shù)項），延遲階數(shù)分別取1,2,3的單位根檢驗。
R的運行結(jié)果如下：

Title:（延遲階數(shù)為1時的ADF檢驗）Augmented Dickey-Fuller Test Test Results:PARAMETER:Lag Order: 1STATISTIC:Dickey-Fuller: -6.9674P VALUE:0.01 Title:（延遲階數(shù)為2階時的ADF檢驗）Augmented Dickey-Fuller Test Test Results:PARAMETER:Lag Order: 2STATISTIC: Dickey-Fuller: -5.4821P VALUE:0.01 Title: （延遲階數(shù)為3時的ADF檢驗）Augmented Dickey-Fuller Test Test Results:PARAMETER:Lag Order: 3STATISTIC: Dickey-Fuller: -4.2721P VALUE:0.01

adfTest檢驗，表明無論延遲階數(shù)取1,2或3， 統(tǒng)計量對應(yīng)的P值都為0.01，都拒絕原假設(shè)：序列非平穩(wěn)。
表明序列為平穩(wěn)序列。結(jié)合前文序列平穩(wěn)性的時序圖檢驗和自相關(guān)圖檢驗，認(rèn)為該序列為平穩(wěn)序列。

（3）序列的白噪聲檢驗：LB統(tǒng)計量

Box.test(licheng, type="Ljung-Box",lag=6) Box.test(licheng, type="Ljung-Box",lag=12)

純隨機性檢驗,p值都小于5%，序列為非白噪聲序列。

（4）平穩(wěn)非白噪聲序列的擬合：ARMA（p,q）模型
針對平穩(wěn)非白噪聲序列，可以運用ARMA類模型進行擬合。
根據(jù)ACF拖尾，PACF2階截尾，建議采用AR（2）模型。

arima(licheng, order = c(2,0,0),method="ML")

用AR(2)模型擬合，這里選擇參數(shù)估計方法method=“ML”。
若估計方法為條件最小二乘法CSS時則不顯示AIC。

R運行結(jié)果：

Call: arima(x = licheng, order = c(2, 0, 0), method = "ML") Coefficients:ar1 ar2 intercept0.7185 -0.5294 11.0223 s.e. 0.1083 0.1067 3.0906 sigma^2 estimated as 365.2: log likelihood = -258.23, aic = 524.46

#注意：若要用中心化的AR(2)模型擬合，即假定不含截距項，則要補充語句include.mean = F，即指令為：
arima(licheng, order = c(2,0,0),method="ML", include.mean = F)
運行該指令，得到不含截距項的AR（2）模型的AIC=532.03。
所以選擇非中心化的模型。

根據(jù)非中心化模型的參數(shù)估計結(jié)果，得到均值為
=11.0223，而
=0.7185，=-0.5294，計算截距項
=8.9380。即此時的AR（2）模型為：

（5）平穩(wěn)非白噪聲序列擬合模型的檢驗
模型的檢驗1：模型的顯著性檢驗。

a= arima(licheng, order = c(2,0,0),method="ML") #把所有的擬合信息都存在一個叫“a”的文件（倉庫）里。 r=a$residuals #用r來保存殘差 Box.test(r,type="Ljung-Box",lag=6,fitdf=2)

對殘差序列進行純隨機性檢驗。
注意：fitdf的值等于（p+q），number of degrees of freedom to be subtracted if x is a series of residuals。當(dāng)檢驗的序列是殘差序列時，需要加上命令fitdf，表示減去的自由度。

R運行結(jié)果：

data: r X-squared = 2.3416, df = 4, p-value = 0.6732 表明當(dāng)檢驗階數(shù)為6階時，殘差序列是白噪聲序列。 Box.test(r,type="Ljung-Box",lag=12,fitdf=2) R運行結(jié)果：Box-Ljung test data: r X-squared = 3.2655, df = 10, p-value = 0.9745

表明當(dāng)檢驗階數(shù)為12階時，殘差序列依然是白噪聲序列。
綜合檢驗階數(shù)為6階的檢驗結(jié)果，模型顯著有效。

模型的檢驗2：參數(shù)的顯著性檢驗。
注意：在時間序列中對參數(shù)顯著性的要求與回歸模型不同，我們更多的是考察模型整體的好壞，而不是參數(shù)。所以，R中的arima擬合結(jié)果中沒有給出參數(shù)的t統(tǒng)計量值和p值。
如需計算參數(shù)的t統(tǒng)計量值和p值，利用下面的公式。
（1）參數(shù) 的t統(tǒng)計量值=

如在本例中，ar1即參數(shù)的估計值為0.7185，其估計值的標(biāo)準(zhǔn)誤差為0.1083，則其t統(tǒng)計量=
。
其他的參數(shù)的t統(tǒng)計量的求法也是一樣的。

（2）P值的求法。
以參數(shù)的t檢驗統(tǒng)計量的P值的計算為例。

p=pt(6.5420,df=56,lower.tail = F)*2

pt()為求t分布的p值的函數(shù)；
6.5420為t統(tǒng)計量的絕對值；
df為自由度=數(shù)據(jù)個數(shù)-參數(shù)個數(shù)，此處，數(shù)據(jù)有59個，參數(shù)3個，所以自由度為56；
lower.tail = F表示所求p值為P[T > t]即右側(cè)概率，如不加入這個參數(shù)表示所求p值為P[T <=t]；
乘2表示p值是雙側(cè)的。

R運行的結(jié)果為：

p=pt(6.5420,df=56,lower.tail = F)*2 p [1] 6.94276e-09 由于該P值小于0.05，表明參數(shù)顯著。 類似的，對ar(2)參數(shù)和常數(shù)項的顯著性進行檢驗。表明模型通過了參數(shù)的顯著性檢驗。

（6）基于擬合模型的平穩(wěn)非白噪聲序列的預(yù)測
接下來，基于該AR（2）模型進行未來5期的預(yù)測：
預(yù)測方法1：直接用predict函數(shù)進行預(yù)測。

predict(arima(licheng, order = c(2,0,0)), n.ahead =5) #預(yù)測未來5期

R運行結(jié)果如下：

$pred #pred為未來5期的預(yù)測值 Time Series: Start = 2009 End = 2013 Frequency = 1 [1] 9.465515 6.215194 8.392736 11.678108 12.885903$se #se表示未來5期的預(yù)測標(biāo)準(zhǔn)差 Time Series: Start = 2009 End = 2013 Frequency = 1 [1] 19.10953 23.53092 23.53226 24.68326 25.22913

置信區(qū)間的求法：以95%的置信區(qū)間的求解為例

licheng.forecast= predict(arima(licheng, order = c(2,0,0)), n.ahead =5) #將未來5期預(yù)測過程中的產(chǎn)生的結(jié)果保存在licheng.forecast變量中 U = licheng.forecast$pred + 1.96* licheng.forecast$se L = licheng.forecast$pred - 1.96* licheng.forecast$se #算出未來5期預(yù)測值的95%置信區(qū)間

R的運行結(jié)果：

U Time Series: Start = 2009 End = 2013 Frequency = 1 [1] 46.92020 52.33580 54.51596 60.05729 62.33500L Time Series: Start = 2009 End = 2013 Frequency = 1 [1] -27.98917 -39.90541 -37.73049 -36.70108 -36.56319

所以，未來5期即2009-2013年的新增里程的預(yù)測結(jié)果為：

預(yù)測年份 預(yù)測值 95％置信區(qū)間 2009 9.465515 （-27.98917,46.92020） 2010 6.215194 （-39.90541,52.33580） 2011 8.392736 （-37.73049,54.51596） 2012 11.678108 （-36.70108,60.05729） 2013 12.885903 （-36.56319,62.33500） ts.plot(licheng,licheng.forecast$pred,col=1:2) #作時序圖，含預(yù)測。 lines(U, col="blue", lty="dashed") lines(L, col="blue", lty="dashed") #在時序圖中作出95%置信區(qū)間，并指定線型為短劃線。lty:連線的線型（1.實線，2.虛線，3.點線，4.點虛線，5.長虛線，6.雙虛線）

預(yù)測方法2：下載并調(diào)用forecast程序包后，用forecast函數(shù)完成預(yù)測。
forecast命令的介紹：

forecast(object,h=,level=)

其中
object：擬合信息文件名；
h：預(yù)測期數(shù)；
level:置信區(qū)間的置信水平。不特殊指定的話，系統(tǒng)默認(rèn)給出置信水平分別為80%和90%的雙層置信區(qū)間。

本例中：

library(forecast) licheng.forecast=forecast(a,5) licheng.forecast Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 2009 9.465302 -15.02455 33.95516 -27.98870 46.91930 2010 6.214789 -23.94131 36.37089 -39.90499 52.33456 2011 8.392250 -21.76556 38.55006 -37.73015 54.51465 2012 11.677647 -19.95516 43.31046 -36.70056 60.05586 2013 12.885518 -19.44684 45.21788 -36.56256 62.33360

畫圖程序如上。

總結(jié)

以上是生活随笔為你收集整理的平稳时间序列分析：ARMA模型的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。