R軟件在“時(shí)間序列分析”教學(xué)中的應(yīng)用
——以ARIMA模型為例

鄧 新 田春雨 葛梅梅 相旭東

1.滁州學(xué)院數(shù)學(xué)與金融學(xué)院 安徽滁州 239000；2.中國(guó)電子科技集團(tuán)公司第五十八研究所 江蘇南京 210000

“時(shí)間序列分析”課程是運(yùn)用時(shí)間序列分析的基本理論與方法，研究隨著時(shí)間的變化，某些事物發(fā)生、發(fā)展的過(guò)程，以尋找其發(fā)展變化的規(guī)律，并預(yù)測(cè)其未來(lái)走勢(shì)的一門(mén)方法論科學(xué)，是統(tǒng)計(jì)學(xué)專(zhuān)業(yè)的必修課程之一。該課程旨在培養(yǎng)學(xué)生運(yùn)用時(shí)間序列分析方法分析、解決實(shí)際問(wèn)題的能力。隨著數(shù)字經(jīng)濟(jì)時(shí)代的到來(lái)，時(shí)間序列分析方法已成為統(tǒng)計(jì)學(xué)、經(jīng)濟(jì)學(xué)等相關(guān)專(zhuān)業(yè)人才必備的數(shù)據(jù)分析方法之一。

為激發(fā)興趣，培養(yǎng)學(xué)生理論聯(lián)系實(shí)際的能力，將統(tǒng)計(jì)軟件引入時(shí)間序列分析教學(xué)中尤其重要。目前，關(guān)于R軟件在統(tǒng)計(jì)學(xué)專(zhuān)業(yè)課程教學(xué)中的應(yīng)用，已有許多學(xué)者做了相關(guān)研究。2011年，閆昭輝通過(guò)聚類(lèi)分析、主成分分析等方面的具體實(shí)際應(yīng)用闡述了該軟件在多元統(tǒng)計(jì)分析教學(xué)中的應(yīng)用；2015年，張毅寧論述了線性、廣義線性、非線性三類(lèi)回歸模型的具體應(yīng)用，以此研究了R軟件在教育統(tǒng)計(jì)中的應(yīng)用；2018年，盧玉桂將R軟件引入抽樣調(diào)查課程，并結(jié)合實(shí)際案例實(shí)現(xiàn)了整群抽樣的樣本抽取和總體參數(shù)估計(jì)；2019年，為培養(yǎng)學(xué)生處理大數(shù)據(jù)的能力，劉君娥以因子分析模型為例，探討了R軟件在統(tǒng)計(jì)學(xué)教學(xué)中的應(yīng)用；2021年，陳超和潘海燕以R軟件輔助醫(yī)學(xué)統(tǒng)計(jì)學(xué)教學(xué)，并以t檢驗(yàn)實(shí)際案例做了進(jìn)一步探究。本研究將以ARIMA模型為例，詳細(xì)說(shuō)明R軟件在時(shí)間序列分析教學(xué)實(shí)踐中的應(yīng)用。

R語(yǔ)言已擁有多個(gè)用于時(shí)間序列分析的程序包。利用這些程序包，能夠完成時(shí)間序列的一系列分析工作：生成時(shí)間序列、繪制圖形、模型識(shí)別、參數(shù)估計(jì)、模型預(yù)測(cè)。除了R軟件，EViews和SAS也常被用于時(shí)間序列分析。相比較于這兩款統(tǒng)計(jì)軟件，R軟件具有如下優(yōu)勢(shì)：第一，與SAS這樣需要昂貴的購(gòu)置費(fèi)及維護(hù)費(fèi)的商業(yè)軟件相比，它是完全免費(fèi)的，通過(guò)官方網(wǎng)站可以免費(fèi)下載安裝最新版本。第二，它占用內(nèi)存非常小，編程代碼簡(jiǎn)單明了，繪圖清晰功能強(qiáng)大。第三，它擁有全球所有優(yōu)秀的統(tǒng)計(jì)軟件包，并向用戶(hù)免費(fèi)開(kāi)放?；谝陨蟽?yōu)點(diǎn)，在實(shí)踐教學(xué)過(guò)程中，不僅可以直接調(diào)用R軟件函數(shù)命令，還可以利用R語(yǔ)言編寫(xiě)自己的程序。這既能滿(mǎn)足教師在教學(xué)過(guò)程中的需求，也能提高學(xué)生的編程能力。

一、ARIMA模型及建模步驟

ARIMA(Auto Regressive Integrated Moving Average)模型，全稱(chēng)是差分自回歸移動(dòng)平均模型，是由博克斯-詹金斯提出的一種專(zhuān)門(mén)用于非平穩(wěn)時(shí)間序列分析和預(yù)測(cè)的方法。ARIMA(p，d，q)模型基本結(jié)構(gòu)如下：

其中?=(1-)，Φ()=1--…-為平穩(wěn)可逆的ARMA(p，q)模型的自回歸系數(shù)多項(xiàng)式，Θ()=1--…-為平穩(wěn)可逆的ARMA(p，q)模型的移動(dòng)平均系數(shù)多項(xiàng)式。

使用ARIMA模型建模的具體過(guò)程：首先，對(duì)獲取的觀察值序列進(jìn)行平穩(wěn)性檢驗(yàn)，若通過(guò)檢驗(yàn)，進(jìn)入下一步，否則，利用差分運(yùn)算實(shí)現(xiàn)平穩(wěn)；其次，進(jìn)行白噪聲檢驗(yàn)，若檢驗(yàn)通過(guò)，分析結(jié)束，否則，擬合ARIMA模型并開(kāi)展模型檢驗(yàn)；最后，模型檢驗(yàn)通過(guò)，分析結(jié)束，否則，重新擬合模型。

二、案例分析——滁州市GDP預(yù)測(cè)分析

GDP，即國(guó)內(nèi)生產(chǎn)總值，是對(duì)一個(gè)國(guó)家、一個(gè)地區(qū)的總體經(jīng)濟(jì)狀況進(jìn)行衡量的重要指標(biāo)。滁州市地處安徽省的東部，是南京與合肥都市圈的重要成員?！冻菔薪y(tǒng)計(jì)年鑒》指出，在2019年，滁州市以2909.06億元的GDP總量，居于安徽省第三名，同時(shí)滁州市的GDP以9.7%的增速排名安徽省第一。滁州經(jīng)濟(jì)如此驚人的成績(jī)和增速，使其備受關(guān)注。通過(guò)滁州市2020年統(tǒng)計(jì)年鑒和《滁州市2020年國(guó)民經(jīng)濟(jì)和社會(huì)發(fā)展統(tǒng)計(jì)公報(bào)》，獲取了1990—2020年滁州市GDP數(shù)據(jù)(單位：億元)：53.97，49.38，62.90，87.13，116.55，148.85，188.21，204.50，210.28，214.28，232.91，259.81，273.08，288.62，330.33，350.26，395.51，492.35，582.47，679.95，837.80，1042.56，1212.45，1418.89，1596.45，1812.28，2035.94，2282.01，2594.07，2909.06，3032.1。將該數(shù)據(jù)分為1990—2018年和2019—2020年兩段，前者用于建立模型，后者用于檢驗(yàn)?zāi)Ｐ皖A(yù)測(cè)效果。

實(shí)驗(yàn)?zāi)康模豪肁RIMA模型對(duì)滁州市GDP進(jìn)行分析及預(yù)測(cè)。

實(shí)驗(yàn)環(huán)境：Window 10，R 4.0.3。

實(shí)驗(yàn)步驟：

(一)創(chuàng)建數(shù)據(jù)集

生成時(shí)間序列(對(duì)原始數(shù)據(jù)進(jìn)行自然對(duì)數(shù)處理)

a<-read.table("滁州市GDP.csv"，sep="，"，header=T)#讀入數(shù)據(jù)

lnGDP<-ts(log(a$GDP)，start=1990)#建立時(shí)間序列{lnGDP}

(二)平穩(wěn)性檢驗(yàn)

對(duì)lnGDP序列進(jìn)行單位根檢驗(yàn)：adf.test(lnGDP)，此時(shí)需要先下載安裝tseries程序包，并用library(tseries)載入。單位根檢驗(yàn)結(jié)果如下。

Augmented Dickey-Fuller Test

data：lnGDP

Dickey-Fuller =-2.8888，Lag order=3，p-value=0.2319

alternative hypothesis：stationary

根據(jù)該檢驗(yàn)的P值大于顯著性水平=005，說(shuō)明該序列是非平穩(wěn)的。

(三)平穩(wěn)化處理(一階差分)

對(duì)序列進(jìn)行一階差分：lnGDP.dif<-diff(lnGDP)，繪制一階差分序列圖：plot(lnGDP.dif)并進(jìn)行單位根檢驗(yàn)：adf.test(lnGDP.dif)。綜合一階差分時(shí)序圖和單位根檢驗(yàn)結(jié)果(P值為0.04664)，一階差分序列是平穩(wěn)的。

(四)白噪聲檢驗(yàn)

Box.test(lnGDP.dif，type="Ljung-Box")

運(yùn)行結(jié)果如下。

Box-Ljung test

data：lnGDP.dif

X-squared=5.9843，df=1，p-value=0.01443

檢驗(yàn)結(jié)果顯示，LB統(tǒng)計(jì)量的P值小于顯著性水平=005，說(shuō)明該一階差分序列屬于非白噪聲序列。

(五)模型識(shí)別

根據(jù)序列平穩(wěn)化過(guò)程可知，=1。為確定p和q的值，分別利用函數(shù)acf和pacf繪制一階差分序列的自相關(guān)系數(shù)(ACF)圖和偏自相關(guān)系數(shù)(PACF)圖。根據(jù)圖1，自相關(guān)系數(shù)和偏自相關(guān)系數(shù)都為1階截尾，所以初步選定模型階數(shù)為(0，1)和(1，0)。

圖1 一階差分自相關(guān)系數(shù)圖和偏自相關(guān)系數(shù)圖

(六)模型優(yōu)化和參數(shù)估計(jì)

分別利用ARIMA(1，1，0)和ARIMA(0，1，1)擬合1990—2018年滁州市l(wèi)n(GDP)序列，其代碼及運(yùn)行結(jié)果如下。

lnGDP.fit<-arima(lnGDP，order=c(1，1，0))

lnGDP.fit

Call：

arima(x=lnGDP，order=c(1，1，0))

Coefficients：

ar1

0.8559

s.e.0.0861

sigma^2 estimated as 0.00655：log likelihood=30.01，aic=-56.01

lnGDP.fit1<-arima(lnGDP，order=c(0，1，1))

lnGDP.fit1

Call：

arima(x=lnGDP，order=c(0，1，1))

Coefficients：

ma1

0.8559

s.e.0.1161

sigma^2 estimated as 0.01062：log likelihood=23.24，aic=-42.49

根據(jù)AIC準(zhǔn)則選擇ARIMA(1，1，0)，確定模型口徑為：

(1-08559)(1-)=，～(0，000655)

等價(jià)為=18559-1-08559-2+，～(0，000655)。

(七)模型檢驗(yàn)

調(diào)用aTSA包中的ts.diag函數(shù)，對(duì)建立的ARIMA(1，1，0)模型進(jìn)行顯著性檢驗(yàn)。根據(jù)白噪聲檢驗(yàn)結(jié)果(左下圖)，各階延遲下白噪聲檢驗(yàn)統(tǒng)計(jì)量的P值都顯著大于顯著性水平，因此認(rèn)為該擬合模型的殘差序列屬于白噪聲序列，即該擬合模型顯著成立。

圖2 擬合模型顯著性檢驗(yàn)圖

(八)模型預(yù)測(cè)

GDP.fore<-forecast(GDP.fit，h=7)

exp(GDP.fore$mean)

GDP<-ts(a$GDP，start=1990)

plot(exp(lnGDP.fore$fitted)，col=4)

lines(GDP，lty=2，col=6)

legend("top"，legend=c("ARIMA模型擬合曲線"，"1990—2018滁州市GDP")，ncol=2，cex=0.8，bty="n"，col=c(4，2)，lty=c(1，2))

調(diào)用forecast函數(shù)，運(yùn)行上述代碼預(yù)測(cè)2019—2023年滁州市GDP(單位：億元)：2894.83，3179.79，3445.85，3691.17，3914.96。經(jīng)計(jì)算，ARIMA(1，1，0)模型的平均相對(duì)誤差為5.36%，較小，說(shuō)明該模型預(yù)測(cè)精度較高。

結(jié)語(yǔ)

ARIMA模型是處理非平穩(wěn)時(shí)間序列的一種分析方法。利用R軟件實(shí)現(xiàn)時(shí)間序列分析，程序簡(jiǎn)單靈活，輸出結(jié)果清晰。通過(guò)R軟件開(kāi)展案例教學(xué)，一方面能進(jìn)一步鞏固和加深學(xué)生對(duì)理論知識(shí)的理解和認(rèn)識(shí)，更好地促進(jìn)學(xué)生對(duì)理論知識(shí)的學(xué)習(xí)；另一方面，有利于提高學(xué)生實(shí)踐動(dòng)手能力，培養(yǎng)創(chuàng)新精神。