基于ARIMA模型的北京市近十年降水量分析

柳向華++楊碩++劉子康++程一帆++楊毅飛

摘 要：降水量是反映一個地區(qū)氣候狀況的重要指標(biāo)，對林木生長有重要影響，與人類生活息息相關(guān)，因此，研究年降水量的演變趨勢和預(yù)測年降水量具有重大意義。文章對北京地區(qū)2004～2015年的年降水量進(jìn)行時間序列分析，運用R軟件建立預(yù)測模型，并通過模型對北京地區(qū)2016～2020年的年降水量進(jìn)行了預(yù)測分析，從而為北京地區(qū)的各個領(lǐng)域制定相應(yīng)政策措施提供參考資料。

關(guān)鍵詞：降水量；時間序列；白噪聲；ARIMA模型；確定性分析

一、 背景

水資源是人類賴以生存的資源，降水是自然界水循環(huán)的一個環(huán)節(jié)。如能有效干預(yù)和控制降水，對于防洪抗旱，涵養(yǎng)水土，指導(dǎo)生產(chǎn)，均為有利。而干預(yù)和控制降水的前提是能對降水量進(jìn)行有效擬合及預(yù)測。目前有多種方法應(yīng)用于降水量的擬合預(yù)測，如運用神經(jīng)網(wǎng)絡(luò)模型，建立微分方程動態(tài)模型，為考察不確定性對擬合預(yù)測效果的影響，用數(shù)值模擬方法建模，建立降雨的隨機模型進(jìn)行擬合預(yù)測，根據(jù)擬合預(yù)測期長短不同，選用不同的擬合預(yù)測方法。其中，美國統(tǒng)計學(xué)家G.E.P.BOX和英國統(tǒng)計學(xué)家G.M.Jenkins聯(lián)合提出的求和自回歸移動平均（autoregressive integrated moving average，ARIMA）模型，作為經(jīng)典時間序列時域分析方法的核心內(nèi)容，在降水量及其成分的擬合預(yù)測中也有應(yīng)用。為準(zhǔn)確描述降水特征，本文將運用R軟件建立ARIMA月降水量時間序列模型，對北京地區(qū)月降水量時間序列數(shù)據(jù)進(jìn)行模型擬合。

二、 模型介紹

（一） 模型的結(jié)構(gòu)

具有如下結(jié)構(gòu)的模型稱為求和自回歸移動平均（autoregressive integrated moving average）模型，簡記為ARIMA（p，d，q）模型：

Φ（B）

SymbolQC@ dxt=Θ（B）εt

E（εt）=0，Var（εt）=σ2ε，E（εtεs）=0，s≠t

E（xsεt）=0，s

SymbolQC@ d=（1-B）d；Φ（B）=1-φ1B-…-φpBp為平穩(wěn)可逆ARMA（p，q）模型的自回歸系數(shù)多項式；Θ（B）=1-θ1B-…-θqBq 為平穩(wěn)可逆ARMA（p，q）模型的移動平滑系數(shù)多項式。

上式可以簡記為：

SymbolQC@ dxt=Θ（B）Φ（B）εt

{εt} 為零均值白噪聲序列。

（二） 模型的性質(zhì)

1. 平穩(wěn)性

假定{xt}服從ARIMA（p，d，q）模型，記φ（B）=Φ（B）

SymbolQC@ d，φ（B）稱為廣義自回歸系數(shù)多項式。

ARIMA模型的平穩(wěn)性完全由廣義自回歸系數(shù)多項式的根的性質(zhì)決定。可以驗證，自回歸系數(shù)多項式的根序列xt的特征根的倒數(shù)，當(dāng)xt隨t趨于無窮趨向于0時序列平穩(wěn)。即要求xt的特征根都在單位圓內(nèi)。

容易判斷，ARIMA（p，d，q）模型的廣義自回歸系數(shù)多項式共有p+d個根，其中p個根1λ1，…，1λp在單位圓外，d個根在單位圓上。

本例中我們采用圖檢驗的方式判斷序列的平穩(wěn)性。

2. 方差齊性

對于平穩(wěn)的時間序列，有方差齊性的性質(zhì)。故白噪聲序列也是典型的平穩(wěn)序列，具有方差齊性。本例中主要對殘差序列進(jìn)行白噪聲檢驗，要求方差齊性。對于殘差異方差序列，我們可以采取方差齊性變換或者ARCH模型進(jìn)一步提取序列有用信息。

（三） 模型建立與定階

1. 定階依據(jù)

根據(jù)AR（p）模型的自相關(guān)系數(shù)拖尾性，偏自相關(guān)系數(shù)截尾性；MA（q）模型自相關(guān)系數(shù)截尾性，偏自相關(guān)系數(shù)圖拖尾性；ARMA（p，q）模型自相關(guān)系數(shù)和偏自相關(guān)系數(shù)均拖尾的性質(zhì)進(jìn)行模型的定階。

2. 確定適當(dāng)模型

在分析實際實際的過程中，我們往往不會有原序列就平穩(wěn)的一些數(shù)據(jù)，根據(jù)Creamer定理，任意一個時間序列{xt}都可以分解成兩部分的疊加，其中一部分是由多項式?jīng)Q定的確定性趨勢成分；另一部分是平穩(wěn)的零均值誤差成分。

（1） 對于確定性趨勢的提取，在本例中我們采用了兩種方式，用ARIMA模型擬合或者用確定性分析提取。

ARIMA模型中，由于Cramer定理的理論保證，我們對非平穩(wěn)序列進(jìn)行差分，d階差分一定可以將充分性信息充分提取。

確定性分析中，直接由序列時序圖判斷用幾次多項式擬合原序列。

（2） 對于零均值誤差部分，它的存在即驗證對于數(shù)據(jù)中的確定信息我們是否提取充分，要求殘差序列是白噪聲序列。對此，我們分別用了LB檢驗統(tǒng)計量和DW檢驗統(tǒng)計量。

3. 季節(jié)模型

在實際數(shù)據(jù)中，數(shù)據(jù)往往呈現(xiàn)季節(jié)效應(yīng)。最常用的是加法模型和乘法模型。若差分后序列的季節(jié)效應(yīng)部分振幅根據(jù)水平的變化不發(fā)生變化，則采用加法模型。否則采用乘法模型擬合季節(jié)性趨勢。

（四） 模型優(yōu)化與比較

通過建立不同的ARIMA模型，擬合原序列。根據(jù)原序列或差分后的序列顯示出的特征，不斷進(jìn)行模型的優(yōu)化。

根據(jù)AIC準(zhǔn)則，AIC的值越小，說明模型的擬合效果越好。

（五） 模型預(yù)測

xt+l的真實值為：

xt+l=（εt+l+Ψ1εt+l-1+…+Ψl-1εt+1）+（Ψlεt+Ψl+1εt-1+…）

由于εt+l，εt+l-1，…，εt+1的不可獲得性，所以xt+l的估計值為：

x^t（l）=Ψ*0εt+Ψ*1εt-1+Ψ*2εt-2+…

我們考慮真實值與預(yù)測值之間的均方誤差：

E[xt+l-x^t（l）]2=（1+Ψ21+…+Ψ2l-1）σ2ε+∑∞j=0（Ψl+j-Ψ2j）2σ2ε

在均方誤差最小原則下，l期預(yù)測值為：

x^t（l）=Ψlεt+Ψl+1εt-1+Ψl+2εt-2+…

三、 ARIMA模型

（一） 確定觀察值序列

我們針對北京市04～15年間的年降水量做ARIMA模型的擬合，搜集的數(shù)據(jù)如下：

（二） 平穩(wěn)性檢驗

得到需要的數(shù)據(jù)后，首先將數(shù)據(jù)變成時間序列形式的變量，并做進(jìn)一步分析。

初步畫出它的時序圖、自相關(guān)系數(shù)圖和偏自相關(guān)系數(shù)圖，檢驗序列的平穩(wěn)性。

選用R語言作為基本的分析工具

a=scan（'降水量.txt'）

a_ts=ts（a，start=c（2004，1），frequency=12）

plot（a_ts，type='l'）

acf（a_ts）

由時序圖和自相關(guān)系數(shù)圖判斷，該序列呈現(xiàn)出不平穩(wěn)性。

1. 尋求平穩(wěn)序列

由于Cramer分解定理的支持，我們對該序列進(jìn)行差分運算，可將其確定性信息提取出來。對此我們對序列進(jìn)行一階差分運算，并檢驗其平穩(wěn)性。

a_ts_diff=diff（a_ts）

plot（a_ts_diff，type='l'）

acf（a_ts_diff）

通過觀察其一階差分的時序圖及自相關(guān)系數(shù)圖，原序列的一階差分序列趨于平穩(wěn)。

同時我們看到原序列的一階差分具有季節(jié)效應(yīng)，即降水量隨月份的變化呈現(xiàn)季節(jié)特征。

觀察一階差分后序列的時序圖，年降水量隨著年份的增加振幅不斷增大，我們采用季節(jié)乘積模型擬合原序列。

2. 模型定階

作出原序列一階差分后的偏自相關(guān)系數(shù)圖

pacf（a_ts_diff）

在上述的自相關(guān)系數(shù)圖和偏自相關(guān)系數(shù)圖中，不考慮季節(jié)因素的影響，我們認(rèn)為原序列滿足ARIMA（11，1，1）

現(xiàn)在考慮季節(jié)效應(yīng)的影響，對一階差分做12步的差分運算，作出其自相關(guān)與偏自相關(guān)系數(shù)圖

a_ts_dif=diff（a_ts，lag=12）

par（mfrow=c（2，1））

acf（a_ts_dif）

pacf（a_ts_dif）

由一階十二步差分的自相關(guān)與偏自相關(guān)圖，認(rèn)為季節(jié)影響可用ARIMA（0，1，1）12模型擬合。

3. 模型建立與參數(shù)估計

根據(jù)選定的模型進(jìn)行參數(shù)估計

estimate=arima（a_ts，order=c（11，1，1），seasonal=list（order=c（0，1，1），period=12））

Call：

arima（x = a_ts，order = c（11，1，1），seasonal = list（order = c（0，1，1），period = 12））

Coefficients：

ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8

-0.0340 0.0099 -0.0534 0.0590 0.0798 -0.0373 0.0379 0.1244

s.e. 0.0882 0.0891 0.0862 0.0855 0.0849 0.0907 0.0891 0.0903

ar9 ar10 ar11 ma1 sma1

0.05000.18100.1961-1.0000-0.8640

s.e.0.08910.08650.08770.07660.1454

sigma^2 estimated as 1054： log likelihood = -653.17， aic = 1334.33

從中我們可以得到原序列的擬合模型的表達(dá)式為

xt=1+B

1+0.034B-0.0099B2+0.0534B3-

0.0590B4-0.0798B5+0.0373B6-

0.0379B7-0.1244B8-0.05B9-

0.181B10-0.1961B11（1-0.1454B12）

4. 模型檢驗

對模型的殘差序列進(jìn)行白噪聲檢驗

tsdiag（estimate）

根據(jù)殘差序列的自相關(guān)系數(shù)圖和p值，p值在0.8到0.1之間，我們不能拒絕殘差序列為白噪聲序列。

5. 模型優(yōu)化

進(jìn)一步我們可以進(jìn)行參數(shù)的顯著性檢驗，將對時間序列影響不明顯的參數(shù)排除，使得模型更精簡。通過偏自相關(guān)系數(shù)圖我們可以看出其2，3，4，5，12階系數(shù)是不明顯的，由此其實我們可以建立疏系數(shù)模型ARIMA（（1，6，7，8，9，10，11），1，1）。在本例中，我們沒有進(jìn)行下一步的參數(shù)顯著性檢驗和疏系數(shù)模型。

四、 確定性分析方法與ARIMA模型的結(jié)合

（一） 對確定性趨勢進(jìn)行擬合

通過原序列的時序圖，我們大致得到年降水量隨著年份的增長呈線性趨勢，故用線性方式以自回歸的方式擬合確定性趨勢

len=length（a_ts）

t=c（1：len）

Tt1=lm（a_ts～t）

得到結(jié)果：

Call：

lm（formula = a_ts ～ t）

Coefficients：

（Intercept） t

36.2128 0.1019endprint

所以趨勢擬合模型為：T=36.2128+0.1019*t

（二） 殘差自相關(guān)檢驗

對原始序列進(jìn)行初步擬合后，我們需要對其進(jìn)行殘差的自相關(guān)檢驗，以確定是否已經(jīng)將序列中的確定性信息提取充分。

殘差自相關(guān)檢驗：dwtest（Tt1）

Durbin-Watson test

data： Tt1

DW = 1.0205，p-value = 9.613e-10

alternative hypothesis： true autocorrelation is greater than 0

根據(jù)p值，我們判斷殘差存在自相關(guān)性，說明信息提取不夠充分，需要對殘差進(jìn)行再次擬合

我們?nèi)圆扇∏懊鍭RIMA模型的方式，對于殘差序列進(jìn)行擬合

res1=Tt1$res

par（mfrow=c（2，1））

acf（res1）

pacf（res1）

觀察自相關(guān)和偏自相關(guān)系數(shù)圖，偏自相關(guān)具有截尾性，并且自相關(guān)3階后遞減趨于零。對于殘差序列運用模型AR（3）進(jìn)行擬合。

rfit1=arima（res1，order=c（2，0，0），include.mean=FALSE）

Call：

arima（x = res1，order = c（2，0，0），include.mean = FALSE）

Coefficients：

ar1 ar2

0.5071 -0.0427

s.e.0.0831 0.0830

sigma^2 estimated as 2308： log likelihood=-762.04， aic = 1530.07

（三） 模型比較

此時我們選用1階延遲序列值為變量再次對殘差序列進(jìn)行擬合

Tt2=lm（a_ts[-1]～a_ts[-len]-1）

可得模型為xt=0.6875x（t-1）

res2=Tt2$res

acf（res2）

pacf（res2）

rfit2=arima（res2，order=c（6，0，0），include.mean=FALSE）

Call：

arima（x = res2，order = c（6，0，0），include.mean = FALSE）

Coefficients：

ar1 ar2 ar3 ar4 ar5 ar6

-0.0455 0.1148 -0.0339 -0.0214 -0.0059 -0.0965

s.e. 0.0829 0.0828 0.0835 0.0830 0.0826 0.0849

sigma^2 estimated as 2557： log likelihood = -763.99， aic = 1541.97

根據(jù)AIC準(zhǔn)則，在上述三個模型中，我們認(rèn)為ARIMA（12，1，1）×ARMA（0，1，1）12模型能更好地擬合北京市近十年的降水量序列。

在用確定性分析方法提取趨勢時，Tt1模型能夠更有效擬合北京市近十年降水量的序列分布

（四） 擬合效果圖

對于上述兩種較好的擬合模型做出其擬合效果圖，并與原序列進(jìn)行比較

我們更進(jìn)一步的認(rèn)定了ARIMA（12，1，1）×ARMA（0，1，1）12對于模型的擬合效果是最好的

（五） 最優(yōu)模型結(jié)構(gòu)

經(jīng)過上述步驟，我們認(rèn)為建立的考慮季節(jié)因素的ARIMA模型能夠有效擬合北京市近十年降水量的序列分布

其模型的具體結(jié)構(gòu)為

xt=1+B

1+0.034B-0.0099B2+0.0534B3-

0.0590B4-0.0798B5+0.0373B6-

0.0379B7-0.1244B8-0.05B9-

0.181B10-0.1961B11（1-0.1454B12）

五、 基于ARIMA預(yù)測

現(xiàn)在利用我們擬合的ARIMA（12，1，1）×ARMA（0，1，1）12模型對未來五年北京市的降水量做預(yù)測

fore=predict（estimate，n.ahead=60）

藍(lán)色部分即為對序列的預(yù)測。

六、 模型評價

經(jīng)過上述步驟的建模我們可以觀測到，ARIMA（12，1，1）×ARMA（0，1，1）12對于原始序列的擬合仍然存在一些偏差，對此我們可以利用確定性分析與Auto-Regreesive 模型結(jié)合進(jìn)行擬合。更多的還可以利用Holt 三參數(shù)模型對于原始序列建模進(jìn)行擬合。通過我們所建立的模型進(jìn)行預(yù)測可以得到，如上圖所示，北京市未來五年的降水量在0～150ml左右波動，不會出現(xiàn)極端天氣如暴雨或干旱等災(zāi)害。

注：本論文只是基于模型進(jìn)行預(yù)測，不保證預(yù)測的準(zhǔn)確性。

參考文獻(xiàn)：

[1] 中國統(tǒng)計年鑒2005～2017.

[2] 王燕.應(yīng)用時間序列分析（第四版）.中國人民大學(xué)出版社.

作者簡介：

柳向華，楊碩，劉子康，程一帆，楊毅飛，北京市，中國礦業(yè)大學(xué)（北京）理學(xué)院2014級。endprint