如有需要，評論私發數據原表（本文利用R語言）

---

文章目錄

• 平穩序列的擬合與預測
• • 1.1DF檢驗
  • 1.2ADF檢驗
  • 1.3模型識別
  • 1.4參數估計
  • 1.5模型檢驗
  • • 1.5.1模型的顯著性檢驗
    • 1.5.3參數的顯著性檢驗
  • 1.6模型優化
  • 1.7序列預測

---

平穩序列的擬合與預測

1.1DF檢驗

類型一：無漂移項自回歸結構
類型二：有漂移項自回歸結構
類型三：帶趨勢回歸結構
【例2-3續(2)】對1915-2004年澳大利亞自殺率序列（每10萬人自殺人口數）進行DF檢驗，判斷該序列的平穩性。

在例2-3中，我們通過圖檢驗方法判斷該序列為非平穩性，但這種判斷帶有很強的個人主觀色彩和經驗主義，現在借助DF統計量，進行序列的平穩性檢驗（$\alpha =0.05$）

R語言中有多個函數可以進行平穩性檢驗，我們利用aTSA包中的adf.test函數，該函數的命令格式為：

adf.test(x,nlag=)

式中：

• x：進行平穩性檢驗的序列名
• nlag：最高延遲階數。
  • 若nlag=1，輸出自回歸0階延遲平穩性檢驗結果；
  • 若nlag=2，輸出自回歸0至1階延遲平穩性檢驗結果；
  • 若用戶不特殊指定延遲階數，系統會自動指定延遲階數。

> #下載并調用aTSA包，使用adf.test函數進行DF檢驗 > library(readxl) > A1_6 <- read_excel("C:/Users/86158/Desktop/A1_6.xlsx") > View(A1_6) > install.packages("aTSA") > library(aTSA) > Suicide<-ts(A1_6$Suicide,start=1915) > adf.test(Suicide,nlag=2)

DF檢驗輸出結果如下：

Augmented Dickey-Fuller Test alternative: stationary Type 1: no drift no trend lag ADF p.value [1,] 0 -1.39 0.179 [2,] 1 -1.32 0.204 Type 2: with drift no trend lag ADF p.value [1,] 0 -1.98 0.337 [2,] 1 -1.31 0.586 Type 3: with drift and trend lag ADF p.value [1,] 0 -2.29 0.449 [2,] 1 -1.65 0.719 ---- Note: in fact, p.value = 0.01 means p.value <= 0.01

檢驗結果顯示，如果序列的結構考慮如上三種類型（六種子類型）的話，$\tau$統計量的P值均顯著大于顯著性水平（$\alpha =0.05$），因此可以判斷，如果序列考慮如上6種結構之一提取確定性信息，則隨機性部分${\xi }_t$都不能實現平穩，也就是說1915-2004年澳大利亞自殺率序列是非平穩序列。

1.2ADF檢驗

【例2-5續(1)】對1900-1998年全球7.0級以上地震發生次數序列進行ADF檢驗，判斷該序列的平穩性。

該序列我們在例2-5通過圖檢驗，判斷該序列平穩，現在基于ADF檢驗，對序列的平穩性進行統計檢驗。

> #讀入序列，并使用adf.test函數進行ADF檢驗 > A1_7 <- read_excel("C:/Users/86158/Desktop/A1_7.xlsx") > View(A1_7) > number<-ts(A1_7$number,start=1900) > library(aTSA) > adf.test(number)

ADF檢驗結果輸出：

Augmented Dickey-Fuller Test alternative: stationary Type 1: no drift no trend lag ADF p.value [1,] 0 -1.585 0.109 [2,] 1 -1.045 0.304 [3,] 2 -0.653 0.445 [4,] 3 -0.509 0.497 Type 2: with drift no trend lag ADF p.value [1,] 0 -5.35 0.0100 [2,] 1 -3.92 0.0100 [3,] 2 -3.18 0.0245 [4,] 3 -2.96 0.0450 Type 3: with drift and trend lag ADF p.value [1,] 0 -5.55 0.0100 [2,] 1 -4.14 0.0100 [3,] 2 -3.51 0.0448 [4,] 3 -3.37 0.0634 ---- Note: in fact, p.value = 0.01 means p.value <= 0.01

檢驗結果顯示，類型二和類型三各州模型的$\tau$統計量的P值均小于顯著性水平（$\alpha =0.05$），所以可以認為該序列顯著平穩。

1.3模型識別

【例4-1】選擇合適的模型擬合1900-1998年全球7級以上地震發生次數序列（數據見表A1-7）

> #繪制該序列自相關圖和偏自相關圖 > acf(number) > pacf(number)

點擊查看自相關圖 和偏自相關圖

從本例的自相關圖可以看出，自相關系數是以一種有規律的方式，按指數函數軌跡衰減的，這說明自相關系數衰減到零不是一個突然截尾的過程，而是一個連續漸變的過程，這是自相關系數拖尾的典型特征，我們可以把拖尾特征形象地描述為“坐著滑梯落水”。

從本例的偏自相關圖可以看出，除了1階偏自相關系數在2倍標準差范圍之外，其他階數的偏自相關系數都在2倍標準差范圍內，這是一個偏自相關系數1階截尾的典型特征，我們可以把這種截尾特征形象地描述為“1階之后高臺跳水”。

本例中，根據自相關系數拖尾，偏自相關系數1階截尾地屬性，我們可以初步確定擬合模型為AR(1)模型。

【例4-2】選擇合適的模型擬合美國科羅拉多州某個加油站連續57天的盈虧(overshort)序列（數據見表A1_8）.

> #讀入數據文件后，繪制時序圖 > library(readxl) > A1_8 <- read_excel("C:/Users/86158/Desktop/A1_8.xlsx") > View(A1_8) > overshort<-ts(A1_8$overshort) > plot(overshort) > #ADF檢驗 > adf.test(overshort)

點擊查看加油站每日盈虧序列時序圖

ADF檢驗結果輸出：

Augmented Dickey-Fuller Test alternative: stationary Type 1: no drift no trend lag ADF p.value [1,] 0 -13.04 0.01 [2,] 1 -7.32 0.01 [3,] 2 -7.01 0.01 [4,] 3 -6.04 0.01 Type 2: with drift no trend lag ADF p.value [1,] 0 -13.10 0.01 [2,] 1 -7.46 0.01 [3,] 2 -7.38 0.01 [4,] 3 -6.61 0.01 Type 3: with drift and trend lag ADF p.value [1,] 0 -12.99 0.01 [2,] 1 -7.39 0.01 [3,] 2 -7.34 0.01 [4,] 3 -6.60 0.01 ---- Note: in fact, p.value = 0.01 means p.value <= 0.01

純隨機性檢驗：

> for(k in 1:2)print(Box.test(overshort,lag=6*k))

純隨機性檢驗結果輸出：

Box-Pierce test data: overshort X-squared = 19.002, df = 6, p-value = 0.00416Box-Pierce test data: overshort X-squared = 28.204, df = 12, p-value = 0.005164

繪制自相關圖和偏自相關圖：

> acf(overshort) > pacf(overshort)

點擊查看加油站每日盈虧序列的自相關圖和偏自相關圖

時序圖顯示該序列沒有明顯的趨勢或周期特征，說明該序列沒有顯著的非平穩特征。進一步進行ADF檢驗，判斷該序列的平穩性。ADF檢驗結果顯示，該序列所有ADF檢驗統計量的P值均小于顯著性水平（$\alpha =0.05$），所以可以確認該序列為平穩序列。再對平穩序列進行純隨機性檢驗，純隨機性檢驗結果顯示，6階和12階延遲的LB統計量的P值都小于顯著性水平（$\alpha =0.05$），所以可以判斷該序列為平穩非白噪聲序列，可以使用ARMA模型擬合該序列。最后考察該序列的樣本自相關圖和偏自相關圖的特征，給ARMA模型定階。

自相關圖顯示除了延遲1階的自相關系數在2倍標準差范圍之外，其他階數的自相關系數都在2倍標準差范圍內波動，且自相關系數衰減沒有顯著的規律性。偏自相關圖顯示出有規律的衰減，這是偏自相關系數拖尾的特征。綜合該序列自相關系數1階截尾和偏自相關系數拖尾的特征，我們將該序列的擬合模型定階為MA(1).

【例4-3】選擇合適的模型擬合1880-1985年全球氣表平均溫度改變值差分序列（全球氣表平均溫度改變值序列見表A1_9)。

> library(readxl) > A1_9 <- read_excel("C:/Users/86158/Desktop/A1_9.xlsx") > View(A1_9) > #讀入數據文件后，對氣表平均溫度改變值序列進行差分運算，對差分后序列繪制時序圖 > dif<-diff(A1_9$change) > dif<-ts(dif,start=1880) > plot(dif) > #ADF檢驗 > adf.test(dif) > #ADF檢驗結果輸出 Augmented Dickey-Fuller Test alternative: stationary Type 1: no drift no trend lag ADF p.value [1,] 0 -13.13 0.01 [2,] 1 -9.60 0.01 [3,] 2 -8.97 0.01 [4,] 3 -7.62 0.01 [5,] 4 -7.68 0.01 Type 2: with drift no trend lag ADF p.value [1,] 0 -13.08 0.01 [2,] 1 -9.58 0.01 [3,] 2 -8.98 0.01 [4,] 3 -7.71 0.01 [5,] 4 -7.85 0.01 Type 3: with drift and trend lag ADF p.value [1,] 0 -13.01 0.01 [2,] 1 -9.53 0.01 [3,] 2 -8.94 0.01 [4,] 3 -7.67 0.01 [5,] 4 -7.82 0.01 ---- Note: in fact, p.value = 0.01 means p.value <= 0.01 > #純隨機性檢驗 > for(k in 1:2)print(Box.test(dif,lag=6*k)) > #純隨機性檢驗結果輸出Box-Pierce test data: dif X-squared = 16.595, df = 6, p-value = 0.01089Box-Pierce test data: dif X-squared = 24.597, df = 12, p-value = 0.01685 > #繪制自相關圖和偏自相關圖 > acf(dif) > pacf(dif)

點擊查看全球氣表平均溫度改變值差分序列時序圖、自相關圖和偏自相關圖

時序圖顯示該序列沒有明顯的趨勢或周期特征。進一步進行ADF檢驗，檢驗結果顯示該序列顯著平穩。接下來檢驗序列的純隨機性，LB檢驗結果顯示該序列為非白噪聲序列。所以全球氣表平均溫度改變值差分序列是平穩非白噪聲序列，可以使用ARMA模型擬合該序列。接下來，序列的自相關圖和偏自相關圖均顯示出不截尾的性質，因此可以嘗試使用ARMA(1,1)模型擬合該序列。

1.4參數估計

R語言中，ARMA模型的參數估計可以通過調用arima函數完成，該函數的命令格式為：

arima(x,order=,include.mean=,method=)

式中：

• x：要進行模型擬合的序列名
• order=c(p,d,q)：指定模型階數
  (1)p為自回歸階數
  (2)d為差分階數，本節不涉及差分問題，故d=0
  (3)q為移動平均階數
• include.mean：模型中要不要包含均值
  (1)include.mean=T，模型中需要擬合均值，這也是系統默認設置
  (2)include.mean=F，模型中不需要擬合均值
• method：指定參數估計方法
  (1)method=“CSS-ML”，默認的是條件最小二乘與極大似然估計混合方法
  (2)method=“ML”，極大似然估計方法
  (3)method="CSS"條件最小二乘估計方法

【例4-1續(1)】使用極大似然估計方法確定1900-1998年全球7級以上地震發生次數序列擬合模型的口徑。

根據該序列自相關圖和偏自相關圖，我們將該序列定階為AR(1)模型，使用極大似然估計確定該模型的口徑，相關命令和輸出結果如下：

> #擬合AR(1)模型 > fit1<-arima(number,order=c(1,0,0),method="ML") > fit1 > #模型擬合結果 Call: arima(x = number, order = c(1, 0, 0), method = "ML") Coefficients:ar1 intercept0.5432 19.8911 s.e. 0.0840 1.3179 sigma^2 estimated as 36.7: log likelihood = -318.98, aic = 643.97

輸出結果分為三部分：

• 說明擬合變量和擬合模型；
• 輸出參數估計結果，第一行輸出的是參數估計值，第二行輸出的是對應參數樣本標準差；
• 輸出殘差序列方差估計值，對數似然函數值和AIC信息量。

根據輸出結果，我們確定該模型口徑為：$$x_t?19.8911=\frac{{\epsilon }_t}{1?0.5432B},Var({\epsilon }_t)=36.7$$
或者對輸出參數進行適當變換$${?}_0=\mu (1?{?}_1)=19.8911\times (1?0.5432)=9.0863$$
該AR(1)模型可以等價表達為：$$x_t=9.0863+0.5432x_{t?1}+{\epsilon }_t,Var({\epsilon }_t)=36.7$$

【例4-2續(1)】確定美國科羅拉多州某個加油站連續57天的盈方序列擬合模型的口徑。

在例4-2中，我們將該序列的擬合模型定價為MA(1)模型，使用條件最小二乘估計方法確定該模型的口徑，相關命令和輸出結果如下：

> #擬合MA(1)模型 > fit2<-arima(overshort,order=c(0,0,1),method="CSS") > fit2 > #模型擬合結果 Call: arima(x = overshort, order = c(0, 0, 1), method = "CSS") Coefficients:ma1 intercept-0.8208 -4.4095 s.e. 0.0996 1.1655 sigma^2 estimated as 2105: part log likelihood = -298.96

根據輸出結果，我們確定該模型口徑為：$$x_t=?4.4095+{\epsilon }_t?0.8208{\epsilon }_{t?1},Var({\epsilon }_t)=2105$$

【例4-3續(1)】確定1880-1985年全球氣表平均溫度改變值差分序列擬合模型的口徑.

在例4-3中，我們將該序列的擬合模型定階為ARMA(1,1)模型，使用條件最小二乘估計確定該模型的口徑，相關命令和輸出結果如下：

> #擬合ARMA(1,1)模型 > fit3<-arima(dif,order=c(1,0,1)) > fit3 > #模型擬合結果 Call: arima(x = dif, order = c(1, 0, 1)) Coefficients:ar1 ma1 intercept0.3926 -0.8867 0.0053 s.e. 0.1180 0.0604 0.0024 sigma^2 estimated as 0.01541: log likelihood = 69.66, aic = -131.32

根據輸出結果，我們確定該模型口徑為：$$x_t?0.0053=\frac{1?0.8867B}{1?0.3926B}{\epsilon }_t,Var({\epsilon }_t)=0.01541$$
或者對輸出參數進行適當變換$${?}_0=\mu (1?{?}_1)=0.0053\times (1?0.3926)=0.0032$$
該ARMA(1,1)模型可以等價表達為：$$x_t=0.0032+0.3926x_{t?1}+{\epsilon }_t?0.8867{\epsilon }_{t?1},Var({\epsilon }_t)=0.01541$$

1.5模型檢驗

1.5.1模型的顯著性檢驗

調用R語言aTSA包中的ts.diag函數，可以進行ARMA模型顯著性檢驗，該函數的命令格式為：

ts.diag(object)

式中：

• object：指ARIMA函數擬合結果，調用ts.diag函數會輸出四個圖：
  (1)左上：殘差序列自相關圖；
  (2)右上：殘差序列偏自相關圖；
  如果該擬合模型不能通過顯著性檢驗，這兩個相關圖可以幫助用戶重新定階。
  (3)左下：殘差序列白噪聲檢驗圖。這個圖的橫軸是延遲階數，縱軸是該延遲階數純隨機性檢驗（Q統計量）的P值，這個圖可以幫助我們直觀判斷擬合模型是否顯著成立。
  如果所有Q統計量的P值都在0.05顯著性參考線之上，可以認為該模型顯著成立。
  (4)右下：殘差序列正態性檢驗QQ圖。我們可以利用QQ圖來對序列進行正態分布假定的檢驗。QQ圖橫軸為正態分布的分位數，縱軸為樣本分位數，如果這兩者構成的點密集地分布在對角線左右，就認為該序列近似服從正態分布。

【例4-1續(2)】檢驗1900-1998年全球7級以上地震發生次數序列擬合模型的顯著性（$\alpha =0.05$）.

我們對該序列擬合了AR(1)模型，擬合模型顯著性檢驗的相關命令和輸出結果如下：

> #擬合模型顯著性檢驗 > library(aTSA) > ts.diag(fit1)

點擊查看全球7級以上地震發生次數序列擬合模型顯著性檢驗圖

考察殘差序列的白噪聲檢驗結果（左下），可以看出各階延遲下白噪聲檢驗統計量的P值都顯著大于0.05（虛線為0.05參考線），我們可以認為這個擬合模型的殘差序列屬于白噪聲序列，即該擬合模型顯著成立。

1.5.3參數的顯著性檢驗

使用arima函數為該序列擬合AR(1)模型，系統會輸出各參數的估計值和系數估計值的樣本標準差，但輸出內容不包括各參數的t統計的值和檢驗P值。如果用戶想獲得參數檢驗統計量的P值，調用t統計量積累分布函數pt即可求得該統計量的P值，pt函數的命令格式為：

pt(t,df=,lower.tail=)

式中：

• t：t統計量的值.
• df：自由度.
• lower.tail：確定計算概率的方向.
  (1)lower.tail=T，計算$Pr(X?x)$.對于參數顯著性檢驗，如果參數估計值為負，選擇lower.tail=T.
  (2)lower.tail=F，計算$Pr(X?x)$.對于參數顯著性檢驗，如果參數估計值為正，選擇lower.tail=F.

【例4-1續(3)】檢驗1900-1998年全球七級以上地震發生次數序列擬合模型參數的顯著性$(\alpha =0.05)$

• 參數顯著性檢驗方法一：因為該序列樣本容量足夠大（序列長度99），所以我們可以基于近似正態分布，非常便捷地得到參數顯著性檢驗結果.

> #調用模型擬合結果 > fit1 > #輸出結果 Call: arima(x = number, order = c(1, 0, 0), method = "ML") Coefficients:ar1 intercept0.5432 19.8911 s.e. 0.0840 1.3179 sigma^2 estimated as 36.7: log likelihood = -318.98, aic = 643.97

因為參數${?}_1$的估計值大于它的兩倍標準差$(0.5432>2\times 0.0840)$,所以參數${?}_1$顯著非零.同理，因為$(19.8911>2\times 1.3179)$，所以參數$\mu$也顯著非零.即我們為1900-1998年全球7級以上地震發生次數序列擬合的AR(1)模型兩參數都顯著非零.

• 參數顯著性檢驗方法二：我們也可以構造參數顯著性檢驗t統計量，調用pt函數求檢驗P值，相關命令和輸出結果如下：

> #構造t統計量，調用pt函數求p值 > t=abs(fit1$coef)/sqrt(diag(fit1$var.coef)) > pt(t,length(number)-length(fit1$coef),lower.tail=F) > #輸出兩參數顯著性檢驗P值ar1 intercept 2.002272e-09 1.687781e-27

因為兩參數t統計量的P值都遠遠小于0.05，所以可以判斷我們為1900-1998年全球7級以上地震發生次數序列擬合的AR(1)模型兩參數都顯著非零.

【例4-2續(2)】對盈虧序列的擬合模型進行檢驗.

> #模型顯著性檢驗 > ts.diag(fit2)

模型顯著性檢驗結果：
點擊查看盈虧序列擬合模型顯著性檢驗圖

> #參數顯著性檢驗(近似方法) > fit2 > #輸出結果 Call: arima(x = overshort, order = c(0, 0, 1), method = "CSS") Coefficients:ma1 intercept-0.8208 -4.4095 s.e. 0.0996 1.1655 sigma^2 estimated as 2105: part log likelihood = -298.96 > #參數顯著性檢驗(精確方法) > t=abs(fit2$coef)/sqrt(diag(fit2$var.coef)) > pt(t,length(overshort)-length(fit2$coef),lower.tail=F) > #輸出t統計量p值ma1 intercept 1.775795e-11 1.919653e-04

模型的顯著性檢驗結果顯示殘差序列可以視為白噪聲序列，所以擬合模型顯著成立.

參數的顯著性檢驗結果顯示，無論使用近似方法還是精確方法，都能得出兩參數顯著非零的結論.

【例4-3續(2)】對1880-1985年全球氣表平均溫度改變值差分序列擬合模型進行檢驗.

> #模型顯著性檢驗 > ts.diag(fit3)

模型顯著性檢驗結果：
點擊查看dif序列擬合模型顯著性檢驗圖

> #參數顯著性檢驗(近似方法) > fit3 > #輸出結果 Call: arima(x = dif, order = c(1, 0, 1)) Coefficients:ar1 ma1 intercept0.3926 -0.8867 0.0053 s.e. 0.1180 0.0604 0.0024 sigma^2 estimated as 0.01541: log likelihood = 69.66, aic = -131.32 > #參數顯著性檢驗(精確方法) > t=abs(fit3$coef)/sqrt(diag(fit3$var.coef)) > pt(t,length(dif)-length(fit3$coef),lower.tail=F) > #輸出t統計量p值ar1 ma1 intercept 6.112701e-04 3.323907e-27 1.555753e-02

模型的顯著性檢驗結果顯示殘差序列可以視為白噪聲序列，所以擬合模型顯著成立.

參數的顯著性檢驗結果顯示，無論使用近似方法還是精確方法，都能得出三參數顯著非零的結論.

1.6模型優化

【例4-7】等時間間隔連續讀取70個某次化學反應的過程數據，構成一時間序列（數據見表A1_10），試對該序列進行擬合$（\alpha =0.05）$.

> library(readxl) > A1_10 <- read_excel("C:/Users/86158/Desktop/A1_10.xlsx") > View(A1_10) > #繪制時序圖 > x<-ts(A1_10$x) > plot(x)

點擊查看化學反應過程時序圖

> #平穩性檢驗 > library(aTSA) > adf.test(x) Augmented Dickey-Fuller Test alternative: stationary Type 1: no drift no trend lag ADF p.value [1,] 0 -1.657 0.093 [2,] 1 -0.872 0.365 [3,] 2 -0.452 0.512 [4,] 3 -0.525 0.489 Type 2: with drift no trend lag ADF p.value [1,] 0 -12.25 0.01 [2,] 1 -5.38 0.01 [3,] 2 -4.36 0.01 [4,] 3 -3.91 0.01 Type 3: with drift and trend lag ADF p.value [1,] 0 -12.61 0.01 [2,] 1 -5.61 0.01 [3,] 2 -4.83 0.01 [4,] 3 -4.41 0.01 ---- Note: in fact, p.value = 0.01 means p.value <= 0.01 > #純隨機性檢驗 > for(k in 1:2)print(Box.test(x,lag=6*k,type="Ljung-Box"))Box-Ljung test data: x X-squared = 21.319, df = 6, p-value = 0.001608Box-Ljung test data: x X-squared = 23.035, df = 12, p-value = 0.02743

平穩性檢驗和純隨機性檢驗顯示該序列為平穩非白噪聲序列，可以使用ARMA模型擬合該序列，下面考察該序列的自相關圖和偏自相關圖的特征，給ARMA模型定階.

> #繪制自相關圖和偏自相關圖 > par(mfrow=c(1,2)) > acf(x) > pacf(x)

點擊查看化學反應過程序列自相關圖和偏自相關圖

根據自相關圖的特征，可能有人會認為自相關系數2階截尾，那么可以對序列擬合MA(2)模型.

根據偏自相關圖的特征，可能有人認為偏自相關系數1階截尾，那么可以對序列擬合AR(1)模型.

下面分別擬合MA(2)模型和AR(1)模型.

> #擬合MA(2)模型 > x.fit1<-arima(x,order=c(0,0,2)) > x.fit1 Call: arima(x = x, order = c(0, 0, 2)) Coefficients:ma1 ma2 intercept-0.3194 0.3019 51.1695 s.e. 0.1160 0.1233 1.2516 sigma^2 estimated as 114.4: log likelihood = -265.35, aic = 538.71 > #擬合AR(1)模型 > x.fit2<-arima(x,order=c(1,0,0)) > x.fit2 Call: arima(x = x, order = c(1, 0, 0)) Coefficients:ar1 intercept-0.4191 51.2658 s.e. 0.1129 0.9137 sigma^2 estimated as 116.6: log likelihood = -265.98, aic = 537.96

MA(2)模型的口徑為：
$$x_t=51.1695+{\epsilon }_t?0.3194{\epsilon }_{t?1}+0.3019{\epsilon }_{t?2},Var({\epsilon }_t)=114.4$$
AR(1)模型的口徑為：
$$x_t?51.2658=\frac{{\epsilon }_t}{1+0.4191B},Var({\epsilon }_t)=116.6$$

對上述兩個擬合模型分別進行模型顯著性檢驗.

> #MA(2)模型顯著性檢驗 > ts.diag(x.fit1)

點擊查看MA(2)模型顯著性檢驗圖

> #AR(1)模型顯著性檢驗 > ts.diag(x.fit2)

點擊查看AR(1)模型顯著性檢驗圖

模型顯著性檢驗結果顯示，這兩個擬合模型的殘差序列都可以視為白噪聲序列，所以這兩個擬合模型均顯著成立.

由于這個化學反應序列的序列長度為70，我們可以基于近似方法判斷參數顯著性，根據模型擬合輸出結果，可以發現無論是MA(2)模型還是AR(1)模型，每個參數估計值的絕對值都大于該參數的2倍標準差，所以這兩個模型的每個參數均顯著非零.

上述分析說明MA(2)模型和AR(1)模型都是這個化學反應序列的有效擬合模型.

同一個序列可以構造兩個甚至更多個擬合模型，多個模型都顯著有效，那么到底該選哪個模型用于統計推斷呢？為了解決這個問題，引進AIC和BIC信息準則的概念，進行模型優化.

【例4-7續(1)】用AIC準則和BIC準則評判例4-7中兩個擬合模型的相對優劣.

> #將兩個模型的AIC和BIC信息量合成一個數據框輸出 > data.frame(AIC(x.fit1),AIC(x.fit2),BIC(x.fit2)) > #輸出結果AIC.x.fit1. AIC.x.fit2. BIC.x.fit2. 1 538.7055 537.9579 544.7033

最小信息量檢驗顯示，無論是使用AIC準則還是使用BIC準則，AR(1)模型都要優于MA(2)模型，所以本例中AR(1)模型是相對最優模型.

AIC準則和BIC準則的提出可以有效彌補根據自相關圖和偏自相關圖定階的主觀性，在有限的階數范圍內幫助我們尋找相對最優擬合模型.

仔細考慮一下，模型定階和模型優化其實是在做同一件事情——為序列選出最優擬合模型，R語言的forecast包提供了auto.arima函數，該函數基于信息量最小原則，可以幫助用戶在一定的范圍內自動識別最優模型的階數，在模型定階時，調用這個函數，可以盡量避免因個人經驗不足導致的模型識別不準確的問題，也可以簡化模型最優化比較的問題.

auto.arima函數的命令格式如下：

auto.arima(x,max.p=,max.q= ic=)

式中

• x：需要定階的序列名
• max.p：自相關系數最高階數，不特殊指定的話，系統默認值為5
• max.q：移動平均系數最高階數，不特殊指定的話，系統默認值為5
• ic：指定信息量準則，有"aic",“bic”,"aicc"三個選項，系統默認AIC準則

有時選擇不同的信息量準則會導致模型的識別階數不同，另外auto.arima函數沒有考慮系數不顯著需要剔除的問題，有時它指定的最優模型階數會高于真實階數，因此，auto.arima函數是給用戶作為定階的參考信息，而不是最優模型的準確信息，用戶需要自己結合自相關圖、偏自相關圖、參數顯著性檢驗等多方面的信息綜合分析，合理定階.

【例4-7續(2)】使用auto.arima函數對化學反應序列進行定階.

> #下載forecast包 > install.packages("forecast") > #調用forecast包 > library("forecast") > #調用auto.arima函數 > auto.arima(x) > #輸出結果 Series: x ARIMA(2,0,0) with non-zero mean Coefficients:ar1 ar2 mean-0.3407 0.1873 51.2263 s.e. 0.1218 0.1223 1.1007 sigma^2 = 117.8: log likelihood = -264.83 AIC=537.66 AICc=538.27 BIC=546.65

輸出結果顯示，基于AIC準則，系統推薦的最優模型為AR(2)模型，但是考慮到該序列偏自相關圖1階截尾，再看上面輸出結果中顯示的ar2的參數估計值為0.1873，它的樣本標準差為0.1223，參數估計值小于2倍的樣本標準差，所以ar2的參數不能認為顯著非零，綜上考慮，基于AIC準則，系統推薦的最優模型階數偏高了，我們應該給該序列擬合AR(1)模型.

如果修改命令，指定使用BIC準則，則系統推薦的最優模型為AR(1)模型，相關命令和輸出結果如下：

> #基于BIC準則，選擇最優擬合模型 > auto.arima(x,ic="bic") > #輸出結果 Series: x ARIMA(1,0,0) with non-zero mean Coefficients:ar1 mean-0.4191 51.2658 s.e. 0.1129 0.9137 sigma^2 = 120: log likelihood = -265.98 AIC=537.96 AICc=538.32 BIC=544.7

1.7序列預測

R語言中有多個函數可以進行ARMA模型預測，我們介紹最常用的一個函數，forecast包中的forecast函數的命令格式如下：

forecast(object,h=,level=)

式中：

• object：擬合模型對象名.
• h：預測期數.
• level：置信區間的置信水平，不特殊指定的話，系統會自動給出置信水平分別為80%和95%的雙層置信區間.
• 若aTSA包同時打開，需要使用forecast::forecast格式調用forecast包中的forecast函數.

【例4-1續(4)】根據1900-1998年全球7級以上地震發生次數的觀察值，預測1999-2008年全球7級以上地震發生次數.

> #調用forecast函數，做10期預測 > library("forecast") > fore1<-forecast::forecast(fit1,h=10) > fore1 > #輸出預測結果Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 1999 17.77725 10.013811 25.54069 5.904094 29.65041 2000 18.74274 9.907694 27.57778 5.230705 32.25477 2001 19.26723 10.139955 28.39451 5.308265 33.22620 2002 19.55217 10.340414 28.76392 5.464007 33.64032 2003 19.70695 10.470420 28.94349 5.580895 33.83301 2004 19.79104 10.547207 29.03487 5.653817 33.92826 2005 19.83672 10.590734 29.08271 5.696204 33.97724 2006 19.86154 10.614915 29.10816 5.720048 34.00303 2007 19.87502 10.628208 29.12183 5.733243 34.01679 2008 19.88234 10.635476 29.12921 5.740482 34.02420 > #輸出預測圖 > plot(fore1,lty=2) > lines(fore1$fitted,col=4)

點擊查看全球7級以上地震發生次數序列擬合與預測效果圖

圖中，虛線為觀察值，實線為擬合值，深色陰影部分為置信水平為80%的預測值置信區間，淺色陰影部分為置信水平為95%的預測值置信區間.

---

如有需要，評論私發數據.

總結

以上是生活随笔為你收集整理的平稳与非平稳序列的拟合及预测的全部內容，希望文章能夠幫你解決所遇到的問題。

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