幾種統計軟件建立ARIMA模型的應用比較

錢 俊，王永波

（南方醫科大學 生物醫學工程學院，廣州 510515）

0 引言

20世紀70年代，Box和Jenkins提出了自回歸移動平均模型，即ARIMA模型，也稱Box-Jenkins模型。ARIMA模型被廣泛應用于時間序列資料的分析，成為經典的時間序列分析和預測方法，現在一般提的時間序列模型，指的是該模型或其某種表述形式[1]。ARIMA模型精確度較高但計算復雜，一般需借助軟件實現。目前，常用統計軟件如SAS、SPSS、R和EVIEWS都有ARIMA建模模塊，這給模型的推廣應用提供了強有力的技術支持[1-4]。

在經典的ARIMA模型中，如果時間序列存在趨勢且經過d階逐期差分可以平穩，則可以建立ARIMA(p,d,q)模型[3]。設{Xt}為非平穩序列，d階差分后的平穩序列為{Zt}，即有 Zt=▽dXt(t＞d)，若{Zt}是 ARMA(p,q)序列，則{Xt}稱ARIMA的d階求和序列，并用ARIMA(p,d,q)表示。模型的一般形式：

式（1）中，p、q為自回歸和移動平均階數，d為平穩化過程中差分的階數；?(B)和θ(B)分別為自回歸算子和移動平均算子。若時間序列存在季節性周期波動，則需要加入季節性算子。對于包含有季節和趨勢的非平穩序列，如果可以通過逐期差分和季節差分使序列平穩化，就可運用ARIMA(p，d，q)(P，D，Q)S模型，一般形式：

式（2）中，P、Q為季節性的自回歸和移動平均階數，D為季節差分的階數，s為季節周期。ΦP(Bs)為季節性P階自回歸算子；ΘQ(Bs)為季節性Q階移動回歸算子。式（2）稱為SARIMA模型或乘積季節ARIMA模型，是隨機季節模型與ARIMA模型的結合，用來描述由于季節性變化（包括季度、月度、周度等變化）或其他一些因素引起的周期性變化的序列[3]。SARIMA模型屬于ARIMA建模中較復雜的情形，本文以此情形為例闡述其建模過程。

1 數據來源和建模步驟

本文利用廣東省2009年1月至2015年12月痢疾的月發病例數進行時間序列分析[5]。通過對發病例數的觀察，該數據呈明顯的非平穩性和季節性，并伴隨一定的周期性波動，根據原始數據序列圖以及疾病的發病特點，擬建立SARIMA模型。

將該段觀察時間序列分為2009年1月至2015年6月和2015年7月至12月兩段，前者用于模型的建立，后者用于評價模型預測效能。雖然ARIMA建模的操作并不復雜，但要建立一個好的或較優的模型卻并非易事[6,7]。數學理論已證明，ARIMA模型的形式并不唯一。盡管不少軟件中有自動建模的程序，如SPSS中的專家建模器以及R軟件的auto.arima()命令，但這些都有局限性。任何軟件都不可能編得很復雜，不可能對所有參數和模型形式都做詳盡的測試，因此自動建模得到的結果并不一定是“最優”，理想的模型仍需研究者自行比較和判斷[1,7]。一般ARIMA建模包含3個步驟：模型識別、參數估計、模型檢驗（診斷），通過對這3個步驟的反復進行，篩選出“最優”模型并進行預測[6]。在模型識別階段，根據ACF和PACF的圖形特征，判斷ARIMA模型的階數，這是一種直觀方法。另一種方法就是嘗試階數p和q不同的取值，利用AIC、BIC等準則，選擇“最優”模型[7]。最后得到的模型應具備兩個要素：模型的殘差序列需通過白噪聲檢驗；在模型參數的簡約性和擬合優度指標的優良性方面取得平衡[1]。

2 軟件建模應用及結果分析

下面分別運用SAS 9.2、SPSS 20、EVIEWS 6.0、R3.3.2這4種軟件建立SARIMA模型，并對建模的步驟和計算結果進行分析。

2.1 SAS建模

在SAS/ETS軟件中，ARIMA建模可以通過編程實現，也可以運用菜單過程步。一般而言，SAS編程建模更為靈活實用，也更符合使用習慣，本文在編程環境下說明建模的過程：（1）創建數據集

導入2009年1月至2015年6月廣東省痢疾的月發病例數，通過intnx()函數創建時間，然后輸入發病例數。

input x@@;/*定義月發病例數變量*/

t=intnx('month','01jan2009'd,_n_-1);/*創建時間集*/format t monyy.;cards;353,411…;/*定義時間格式，導入數據*/

（2）模型識別

調用程序PROC GPLOT，繪制序列圖{Xt}，觀察數據的平穩性（圖1）。經觀察發現序列呈現一個略降的長期趨勢和一個周期長度為一年的穩定的季節變動。因此首先消除季節影響，對該序列做12步差分，然后為了消除長期趨勢的影響，運用dif()函數對原始序列進行差分，發現一次差分后序列{??12Xt}基本平穩，則d=1和D=1，s=12。

圖1 廣東省痢疾的月發病病例數時序（SAS）

在SAS中PROC ARIMA程序功能強大，其含有IDENTIFY（模型識別），ESTIMATE（參數估計），FORECAST（短期預測）這三條命令。為了進一步判斷其平穩性，調用PROC ARIMA程序中IDENTIFY命令對序列進行識別，考察差分后序列{??12Xt}的ACF和PACF圖，初步判斷模型中階數p=1,q=1。但為了避免主觀判斷偏差，再運行IDENTIFY命令中自動識階選項MINIC()，計算指定范圍內“最優”模型的階數。運行得最小BIC(0，1)=7.7388，判斷可能最優模型的階數p=0,q=1。再根據圖中ACF和PACF圖形在k=12處顯著，k=24處不顯著，初步判斷季節效應的階數P=1，Q=1。

proc arima; /*建立ARIMA模型*/

identify var=x(1,12)nlag=24 minic p=(0:5)q=(0:5);/*判斷最優模型的階數*/

（3）參數估計

調用PROC ARIMA程序中ESTIMATE命令，對可能的最優模型SARIMA(0,1,1)(1,1,1)12進行參數估計；結果顯示常數項以及參數P=0和Q=0的原假設檢驗P值大于0.05（表1），結果不顯著。說明此模型并非最優模型。

proc arima;estimate p=(12)q=(1)(12);/*擬合帶常數項的SARIMA模型*/

表1 備選模型的參數估計及檢驗（SAS）

（4）模型檢驗

在延遲6階，12階和18階的Ljung-Box統計量(LB統計量)檢驗的P值全部顯著大于0.05，此模型的殘差屬于白噪聲序列（表2）。但由于參數估計得到的系數并不符合“顯著性”要求，需對模型的階數進行調整。在自動識階的過程，根據BIC準則，BIC(0，1)=7.7388最小，BIC(1，0)=7.7498次之，BIC(1，1)=7.7985較小；而季節模型的階數P和Q的取值一般不超過2，對上述參數組合的取值進行試驗，重復步驟（3）和步驟（4）進行比較。根據模型參數估計，擬合效果以及殘差白噪聲檢驗的結果進行綜合判斷篩選最優模型，則SARIMA(1,1,0)(0,1,1)12為最優模型。該模型所有參數都通過了顯著性檢驗（表1）；殘差通過白噪聲檢驗（表2）；擬合優度檢驗的統計量中，除指標BIC略高，其余各項指標值AIC、SBC等均最小（表2），此時得到“最優”模型的數學表達式：

表2 備選模型的殘差白噪聲檢驗和擬合優度檢驗（SAS）

（5）預測

運用SARIMA(1,1,0)(0,1,1)12模型對2015年7月至12月痢疾的發病數進行預測，可以得到具體預測值以及區間估計值（圖2）。將實際值和預測值進行比較，結果顯示，各月實測值都落入了預測值的置信區間內，擬合平均相對誤差MAPE=10.63%，預測平均相對誤差10.33%，該模型具有較好的預測效能（見下頁表3）。

proc arima;estimate p=1 q=(12)noconstant；/*擬合不帶常數項的模型*/

forecast lead=6 id=t out=out;/*預測后6個月的數據*/

圖2 SARIMA(1,1,0)(0,1,1)12模型的預測值和置信度為95%置信區間（SAS）

表3 SAS、SPSS、EVIEWS、R軟件操作ARIMA模型的比較

SAS編程建模的過程中，每個步驟都能得到比較詳盡的結果，建模方法靈活，功能強大。比如，模型識別時，可以自動識階；參數估計時可以選擇條件最小二乘法、無約束最小二乘法和最大似然法，nonconstant選項則模型不帶常數項[8,9]。

2.2 SPSS建模

SPSS中的ARIMA建模操作簡單，只需進行幾個參數的設置即可運行。特別地，SPSS軟件對時間序列分析有強大的自動建模功能，即“專家建模器”，可以由軟件自動生成模型，提高建模速度。

2.2.1 專家建模器建模

（1）導入原始數據并創建時間序列：將數據導入數據編輯器，設置痢疾病例數變量“x”（發病數）為數值變量，將時間變量“t”（日期）定義為日期和時間變量，設定為“年、月”，則建立時間序列｛Xt｝。

（2）專家自動建模：打開菜單項“分析”→“預測”→“創建模型”，在“時間序列建模器”選項卡中確定因變量x（發病數），在選項卡的“方法”選項中選擇“專家建模器（僅限ARIMA模型）”，不需要設置具體參數。軟件將自動計算，建立模型SARIMA(0,0,2)(1,0,0)12（表4）。

表4 備選模型的參數估計結果（SPSS）

2.2.2 模型識別與參數估計

建模的思路和操作步驟和SAS軟件大致相同，只是調用菜單項來實現。

（1）模型識別

“分析”→“預測”→“序列圖”，畫出時間序列｛Xt｝以及差分、季節差分后的時序圖，判斷平穩性；“自相關”畫出時間序列｛Xt｝的ACF和PACF圖形，初步判斷模型的階數。

（2）參數估計和模型檢驗

“分析”→“預測”→“創建模型”，打開“時間序列建模器”選項卡，在選項卡的“方法”選項中選擇“ARIMA模型”，輸入自回歸、差分和移動平均數的階數（包括季節因子的階數）。SPSS中沒有自動識階的功能，因此需要把所有可能階數的組合都計算一遍，建立多個模型，對比它們的參數估計、擬合優度檢驗和殘差白噪聲檢驗等計算結果，選擇“最優”模型。本例通過上述計算，篩選得到模型SARIMA(1,1,0)(0,1,1)12，其參數估計的各項結果符合“顯著性”要求（P＜0.05）（表4）；模型擬合結果的標準化BIC值和MAPE最小，決定系數R2最大（表5）。殘差白噪聲檢驗中LB統計量P值都大于0.05，說明殘差序列為白噪聲序列。但專家建模結果SARIMA(0,0,2)(1,0,0)12的LB統計量值為57.36(P＜0.001)，顯示不是白噪聲序列，模型并不理想。

表5 備選模型的殘差白噪聲檢驗和擬合優度檢驗的結果（SPSS）

（3）模型預測

在上述的“時間序列建模器”選項卡，選擇“擬合值”以及定義預測階段，則可以計算預測值的點估計和區間估計，輸出圖形。此例中，兩個模型各月的實測值都落在預測值95%置信區間內，但專家建模結果SARIMA(0,0,2)(1,0,0)12模型預測值的相對誤差14.61%，預測精度差（表3）。

2.3 R軟件

R軟件提供了彈性、互動的環境分析和數據處理功能。它可以輕松地加載以庫或者程序包的形式存在的補充工具，里面含有各種數學和統計計算的函數，以實現一些復雜的建模功能。在R中建立ARIMA模型，需先加載程序包FORECAST、TSERIES、TSA[4,7]，編程建模的思路和SAS基本相同：

（1）創建數據集

生成時間序列，定義為月度數據。

Xt=ts(c(353,411,…),start=c(2009/01),frequency=12)#建立時間序列｛Xt｝

（2）模型識別

繪制序列圖，分析時序特性；根據ACF和PACF等結果對序列進行識別、定階。R軟件的程序包TSA中，armasubsets()函數有自動識階功能，它根據最小BIC準則來挑選“最優”模型，結果具有參考價值。本例中自動識階結果如下頁圖3，提示模型的階數p=2，q=1。根據ACF和PACF圖形特征，結合自動識階的結果，可以設定階數的取值范圍，建立多個備選模型進行比較。

plot(Xt)#繪制序列圖

acf(as.vector(Xt),lag.max=24)#序列自相關圖ACF

pacf(as.vector(Xt),lag.max=24)#序列偏相關圖PACF

resbic=armasubsets(y=Xt,nar=7,nma=7,ar.method='ols')#ARIMA模型自動識階

plot(resbic)#不同ARIMA模型的BIC值

（3）參數估計

本例對階數p≤2，q≤1，季節模型階數P和Q不超過2的情形進行參數組合，建立SARIMA模型，通過（3）和（4）這兩個步驟反復進行，比較備選模型的各項指標值，篩選得“最優”模型 SARIMA(1,1,0)(0,1,1)12：

m.Xt=arima(Xt,order=c(1,1,0),seasonal=list(order=c(0,1,1),period=12))#建立模型accuracy(m.Xt)#計算模型的各項擬合指標值

圖3 不同p,q階數建立的ARIMA模型的BIC值（R軟件）

R軟件的FORECAST程序包還提供了auto.arima（）函數，具有自動建模功能。本例中調用函數auto.arima(Xt)，得到最優模型SARIMA(1,1,0)(1,0,0)12及參數估計的結果（表6）。但R軟件參數估計只提供系數和標準誤，顯著性檢驗需自行計算統計量t值（系數除以標準誤的平方）來推斷。通過計算，表6中兩個備選模型的系數都符合“顯著性”要求（P＜0.05）。

表6 備選模型的參數估計和模型檢驗結果（R軟件）

（4）模型檢驗

R軟件可以計算各延遲階數的LB統計量及P值，并通過圖形直觀表達。從圖4看出，上述兩個模型的殘差白噪聲檢驗，P值均大于0.05，說明所建立模型的殘差通過白噪聲檢驗。各延遲階數的LB統計量具體值也可以計算，如lag=18時，SARIMA(1,1,0)(0,1,1)12模型LB統計量16.405(P=0.5643)；自動建模SARIMA(1,1,0)(1,0,0)12模型LB統計量20.83(P=0.2882)。

BOX.test（m.Xt$residuals,lag=18,type=”Ljung-Box”）#計算LB統計量

tsdiag(m.Xt,gof=24,omit.initial=F)#LB統計量對應的P值圖

圖4 備選模型的殘差白噪聲檢驗（R軟件）

（5）模型預測

預測指定時間范圍序列｛Xt｝的值和置信區間，預測圖形的繪制比較靈活、輸出美觀。圖5是圖形的輸出，從具體預測值計算，擬合效果SARIMA(1,1,0)(0,1,1)12模型較好，MAPE為8.67%，自動建模SARIMA(1,1,0)(1,0,0)12模型預測精度較高，平均相對誤差小，為6.64%（表3）。

圖5 備選模型的預測值和置信區間（R軟件）

2.4 EVIEWS建模

EVIEWS是廣泛應用的計量經濟學軟件，能進行傳統的時間序列分析，建立各種時序模型[3,10]。使用該軟件可以采用編程建模和菜單過程步建模兩種方式，本文以菜單功能實現說明ARIMA建模的步驟：

（1）導入數據

建立一個新的 Workfile，“File”→“Workfile Create”，在對話框輸入起始日期與結束日期，新建時間序列，將案例中的數據導入，即可建立時間序列｛Xt｝。

（2）模型識別

自相關和偏自相關函數是識別模型的主要工具。先通過菜單項對時間序列｛Xt｝的時序特性進行識別，“View”→“Graph Option”，設置繪制序列對話框，繪制序列圖、序列差分圖等。再運行“Quick”→“Series Statistic”→“Correclogram”，得出序列的ACF和PACF等結果，可根據這些信息對序列模型進行識別、定階，根據圖6結果，初步判定p=1和q=1，P=1和Q=1。

圖6 差分后序列{??12Xt}ACF和PACF結果（EVIEWS）

（3）參數估計

EVIEWS中參數估計采用非線性算法。建模過程為“Quick”→“Estimate Equation”，打開方程定義對話框，若擬建立的模型為SARIMA(1,1,0)(0,1,1)12，就輸入相關的命令行：diffXt ar(1)ma(12)，就可得到結果輸出(diffXt代表差分后序列{??12Xt})。

圖7 SARIMA(1,1,0)(0,1,1)12模型的參數估計結果（EVIEWS）

圖7上半部分是參數估計系數及顯著性結果，下半部分是模型擬合結果：R2=0.46，AIC=10.26，SC=10.32。

（4）模型檢驗

檢驗殘差序列是否為白噪聲，在方程輸出窗口選擇“views”→“Residual Tests”→“Correlogram-Q-Statistics”，在彈出的對話框中輸入最大的滯后期lag=18，得LB統計量Q=18.940（P=0.27），其他延遲階數LB統計量對應的P值都大于0.05，模型通過白噪聲檢驗。

（5）模型預測

在方程輸出窗口選擇“forecast”，預測方法可以選擇靜態預測或追溯預測（Dynamic forecast），動態（向前多步）預測。得到擬合結果的MAPE為8.86%，Theil系數為0.047等結果（圖8），預測值的平均相對誤差為14.28%（表3）。

圖8 ARIMA模型的預測和擬合結果（EVIEWS）

3 比較討論

4種軟件建模方法和應用中：SPSS操作過程最為簡便，適合非專業統計建模的需求。其“專家建模器”操作簡單，可以快速建模；運用菜單項建模也是“按部就班”，不需要復雜操作，有利于初學者使用；但算法呆板，靈活性差也是其不足之處。EVIEWS軟件易學易用，輸出整齊美觀，是計量經濟學常用的軟件，但算法和功能不夠全面有時會降低建模的精度。SAS軟件提供很多命令和選項，建模功能強大，用戶可以根據自己的需求靈活建模。但SAS輸出不夠美觀，SAS編程有一定難度，需要應用者有一定編程基礎。R軟件是免費軟件，建模靈活多變，目前的3.3.2版本已具有自動識階、自動建模等多種功能；但R軟件中會不斷更新程序包，一方面使得建模功能越來越強大，另一方面也需應用者不斷學習探索；也有一定的學習難度。

以上運用4種軟件進行ARIMA建模，不同方法得到了不同的“最優”模型和預測結果。在本實例分析中，SARIMA(1,1,0)(0,1,1)12模型是編程得到的“最優”模型，擬合精度若以指標MAPE比較，R軟件結果最佳；預測精度以預測平均相對誤差比較，EVIEWS最差。這是因為建模過程中，不同軟件提供或默認的算法不同，“最優”模型盡管形式相同，但參數估計的系數不同，擬合和預測的結果就略有差異。在軟件自動建模方面，SPSS專家建模器得到的SARIMA(0,0,2)(1,0,0)12模型雖然參數估計結果均顯著不為0，但殘差白噪聲檢驗未通過，因此擬合和預測效果都不理想。R軟件自動建模得到的SARIMA(1,1,0)(1,0,0)12模型，在模型檢驗中各項指標AIC、BIC、MAPE等不是最優，但在此實例中預測精度最高（預測平均相對誤差6.64%）（表3）。無論是SPSS還是R軟件，自動建模的結果僅可做為參考。在編程建模時，結合ACF、PACF圖和自動識階結果，篩選“最優”模型，避免出現主觀偏差。

總之，不同軟件在ARIMA建模各有特色和優點，應用者可依據自身專業背景和建模需求選擇不同的軟件建立ARIMA模型。