幾種統(tǒng)計軟件建立ARIMA模型的應(yīng)用比較

錢 俊，王永波

（南方醫(yī)科大學(xué) 生物醫(yī)學(xué)工程學(xué)院，廣州 510515）

0 引言

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

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

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

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

1 數(shù)據(jù)來源和建模步驟

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

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

2 軟件建模應(yīng)用及結(jié)果分析

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

2.1 SAS建模

在SAS/ETS軟件中，ARIMA建?？梢酝ㄟ^編程實現(xiàn)，也可以運用菜單過程步。一般而言，SAS編程建模更為靈活實用，也更符合使用習(xí)慣，本文在編程環(huán)境下說明建模的過程：（1）創(chuàng)建數(shù)據(jù)集

導(dǎo)入2009年1月至2015年6月廣東省痢疾的月發(fā)病例數(shù)，通過intnx()函數(shù)創(chuàng)建時間，然后輸入發(fā)病例數(shù)。

input x@@;/*定義月發(fā)病例數(shù)變量*/

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

（2）模型識別

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

圖1 廣東省痢疾的月發(fā)病病例數(shù)時序（SAS）

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

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

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

（3）參數(shù)估計

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

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

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

（4）模型檢驗

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

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

（5）預(yù)測

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

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

forecast lead=6 id=t out=out;/*預(yù)測后6個月的數(shù)據(jù)*/

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

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

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

2.2 SPSS建模

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

2.2.1 專家建模器建模

（1）導(dǎo)入原始數(shù)據(jù)并創(chuàng)建時間序列：將數(shù)據(jù)導(dǎo)入數(shù)據(jù)編輯器，設(shè)置痢疾病例數(shù)變量“x”（發(fā)病數(shù)）為數(shù)值變量，將時間變量“t”（日期）定義為日期和時間變量，設(shè)定為“年、月”，則建立時間序列｛Xt｝。

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

表4 備選模型的參數(shù)估計結(jié)果（SPSS）

2.2.2 模型識別與參數(shù)估計

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

（1）模型識別

“分析”→“預(yù)測”→“序列圖”，畫出時間序列｛Xt｝以及差分、季節(jié)差分后的時序圖，判斷平穩(wěn)性；“自相關(guān)”畫出時間序列｛Xt｝的ACF和PACF圖形，初步判斷模型的階數(shù)。

（2）參數(shù)估計和模型檢驗

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

表5 備選模型的殘差白噪聲檢驗和擬合優(yōu)度檢驗的結(jié)果（SPSS）

（3）模型預(yù)測

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

2.3 R軟件

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

（1）創(chuàng)建數(shù)據(jù)集

生成時間序列，定義為月度數(shù)據(jù)。

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

（2）模型識別

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

plot(Xt)#繪制序列圖

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

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

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

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

（3）參數(shù)估計

本例對階數(shù)p≤2，q≤1，季節(jié)模型階數(shù)P和Q不超過2的情形進行參數(shù)組合，建立SARIMA模型，通過（3）和（4）這兩個步驟反復(fù)進行，比較備選模型的各項指標值，篩選得“最優(yōu)”模型 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階數(shù)建立的ARIMA模型的BIC值（R軟件）

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

表6 備選模型的參數(shù)估計和模型檢驗結(jié)果（R軟件）

（4）模型檢驗

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

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

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

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

（5）模型預(yù)測

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

圖5 備選模型的預(yù)測值和置信區(qū)間（R軟件）

2.4 EVIEWS建模

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

（1）導(dǎo)入數(shù)據(jù)

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

（2）模型識別

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

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

（3）參數(shù)估計

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

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

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

（4）模型檢驗

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

（5）模型預(yù)測

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

圖8 ARIMA模型的預(yù)測和擬合結(jié)果（EVIEWS）

3 比較討論

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

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

總之，不同軟件在ARIMA建模各有特色和優(yōu)點，應(yīng)用者可依據(jù)自身專業(yè)背景和建模需求選擇不同的軟件建立ARIMA模型。