ARMA模型擬合指令：

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

其中:

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

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

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

（2）序列的平穩性檢驗：
時序圖，自相關圖，ADF檢驗

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

圖1 時序圖

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

圖2 自相關圖

圖3 偏自相關圖

單位根檢驗的做法：

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

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

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

adfTest檢驗，表明無論延遲階數取1,2或3， 統計量對應的P值都為0.01，都拒絕原假設：序列非平穩。
表明序列為平穩序列。結合前文序列平穩性的時序圖檢驗和自相關圖檢驗，認為該序列為平穩序列。

（3）序列的白噪聲檢驗：LB統計量

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

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

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

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

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

R運行結果：

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。
所以選擇非中心化的模型。

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

（5）平穩非白噪聲序列擬合模型的檢驗
模型的檢驗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。當檢驗的序列是殘差序列時，需要加上命令fitdf，表示減去的自由度。

R運行結果：

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

表明當檢驗階數為12階時，殘差序列依然是白噪聲序列。
綜合檢驗階數為6階的檢驗結果，模型顯著有效。

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

如在本例中，ar1即參數的估計值為0.7185，其估計值的標準誤差為0.1083，則其t統計量=
。
其他的參數的t統計量的求法也是一樣的。

（2）P值的求法。
以參數的t檢驗統計量的P值的計算為例。

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

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

R運行的結果為：

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

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

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

R運行結果如下：

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

置信區間的求法：以95%的置信區間的求解為例

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

R的運行結果：

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年的新增里程的預測結果為：

預測年份 預測值 95％置信區間 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) #作時序圖，含預測。 lines(U, col="blue", lty="dashed") lines(L, col="blue", lty="dashed") #在時序圖中作出95%置信區間，并指定線型為短劃線。lty:連線的線型（1.實線，2.虛線，3.點線，4.點虛線，5.長虛線，6.雙虛線）

預測方法2：下載并調用forecast程序包后，用forecast函數完成預測。
forecast命令的介紹：

forecast(object,h=,level=)

其中
object：擬合信息文件名；
h：預測期數；
level:置信區間的置信水平。不特殊指定的話，系統默認給出置信水平分別為80%和90%的雙層置信區間。

本例中：

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

畫圖程序如上。

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。