中國糧食供需研究及預(yù)測
——以小麥為例

龐碧玉,馮愛芬,曹振雪,胡啟帆,王琰

（河南科技大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院,河南洛陽471023）

糧食問題是關(guān)系民生的問題,也對農(nóng)業(yè)經(jīng)濟發(fā)展有重要的影響.隨著農(nóng)業(yè)經(jīng)濟的發(fā)展,關(guān)于糧食供需的數(shù)據(jù)越來越多,有效地預(yù)測糧食的供應(yīng)和需求可以使農(nóng)業(yè)生產(chǎn)者了解信息,及時調(diào)整生產(chǎn)過程中可能出現(xiàn)的供需矛盾,使糧食產(chǎn)品市場平穩(wěn)運行.

目前國內(nèi)用于糧食產(chǎn)量預(yù)測的方法主要有無偏灰色馬爾可夫模型、BP神經(jīng)網(wǎng)絡(luò)模型、多元線性回歸模型、支持向量機、灰色GM（1,1）模型等.尹邦華等[1]改進了無偏灰色馬爾可夫模型，周永生等[2]采用BP神經(jīng)網(wǎng)絡(luò)模型，王啟平[3]利用多元線性回歸分析方法，趙桂芝等[4]基于混沌-AVM-PSO的模型，楊克磊等[5]應(yīng)用灰色GM（1,1）模型,他們通過建立糧食產(chǎn)量預(yù)測模型,對我國糧食產(chǎn)量進行預(yù)測,但均有一定程度的地域局限性和適應(yīng)條件.國外對產(chǎn)量預(yù)測的研究方法主要采用氣象產(chǎn)量預(yù)測法、遙感技術(shù)、統(tǒng)計動力學(xué)模擬法.這3種方法預(yù)測誤差通常在5%～10%,預(yù)測精度不高[6-8].考慮到中國小麥生產(chǎn)是一個復(fù)雜的動態(tài)系統(tǒng),小麥供需數(shù)據(jù)具有趨勢性、波動性的特點,故從宏觀角度和數(shù)據(jù)本身出發(fā),消除地域限制建立預(yù)測模型.

本文采用ARIMA模型構(gòu)建計量經(jīng)濟學(xué)模型進行多維度分析討論,排除出口量、區(qū)域性差異、季節(jié)等周期性因素的影響,對產(chǎn)量、進口量和需求量進行預(yù)測,并以中國糧食中的主要農(nóng)作物小麥為例進行分析.

1 ARIMA模型

ARMA模型是平穩(wěn)時間序列模型,而糧食產(chǎn)量、進口量和需求量趨勢變化具有波動性、不確定性,是一種非平穩(wěn)隨機過程,所以本文引入了差分算子I（可使原不平穩(wěn)的時間序列經(jīng)d階差分后改進為平穩(wěn)時間序列）,建立糧食產(chǎn)量、進口量和需求量趨勢預(yù)測ARIMA模型.

差分整合移動平均自回歸模型（Auto-regressive integrated movingaveragemodel,ARIMA）是平穩(wěn)時間序列預(yù)測分析模型之一.ARIMA模型考慮研究對象本身的歷史數(shù)據(jù)隨時間的變化而變化的規(guī)律,并以此預(yù)測未來值.即以時間代替各種影響因素,探尋分析事件的變化特征和發(fā)展趨勢[9-10].ARIMA模型由自回歸模型（AR模型）、移動平均模型（MA模型）與差分算子I（d）組成.

1.1 AR模型

AR模型的參數(shù)p代表預(yù)測模型中采用的時序數(shù)據(jù)本身的滯后數(shù).其p階AR模型表達式為

式（1）中：θi（i=1,2,…,p）為AR模型系數(shù),et為隨機干擾項,c1為常數(shù).

定義B為后移算子,B滿足條件

式（1）和式（2）聯(lián)立得到

1.2 MA模型

MA模型參數(shù)q代表預(yù)測模型中采用的預(yù)測誤差的滯后數(shù).其q階MA模型的數(shù)學(xué)表達式為

式（4）中：Wj（j=1,2,…,q）為MA模型系數(shù),et為隨機干擾項,c2為常數(shù).

1.3 差分算子I

差分算子可使不平穩(wěn)的時間序列經(jīng)d階差分后改進為平穩(wěn)時間序列

綜上,ARIMA模型的數(shù)學(xué)表達式為

式（5）中：Θp（B）和Wq（B）是特征多項式,et為隨機干擾項,c為常數(shù).即得到可用于糧食供需預(yù)測的ARIMA時間預(yù)測模型.

2 中國小麥產(chǎn)量ARIMA模型

2.1 數(shù)據(jù)準(zhǔn)備與處理

本文提取《中國統(tǒng)計年鑒》及“中華人民共和國農(nóng)業(yè)農(nóng)村部”網(wǎng)站1979—2018年中國小麥產(chǎn)量,以萬t為產(chǎn)量單位,以年為時間單位觀察數(shù)據(jù).

對樣本數(shù)據(jù)進行處理,剔除異常數(shù)據(jù),使預(yù)測值盡量準(zhǔn)確.數(shù)據(jù)顯示,1949—1978這30 a中,數(shù)據(jù)缺失現(xiàn)象嚴(yán)重,在時間序列模型中,這段時間樣本不具有分析意義,刪除樣本數(shù)據(jù).

2.2 數(shù)據(jù)平穩(wěn)性檢驗

根據(jù)上述數(shù)據(jù)準(zhǔn)備與處理原則,做出原始產(chǎn)量時間序列圖,見圖1.

圖1 原始產(chǎn)量時間序列Fig.1 Original yield timeseriesdiagram

由圖1可知,原始數(shù)據(jù)趨勢圖不具有平穩(wěn)性,對原始數(shù)據(jù)進行一階差分處理,處理后產(chǎn)量時間序列見圖2.

圖2 一階差分處理后產(chǎn)量時間序列Fig.2 Yield timeseries diagramafter first orderdifference treatment

由圖2可知,數(shù)據(jù)在“0”上下浮動,初步判斷一階差分以后數(shù)據(jù)平穩(wěn).

2.3 模型確定

對一階差分后的數(shù)據(jù)進行自相關(guān)與偏自相關(guān)檢驗,創(chuàng)建自相關(guān)圖與偏自相關(guān)圖,根據(jù)圖片的截尾、拖尾情況估計p、q的數(shù)值.一階差分后自相關(guān)與偏自相關(guān)圖見圖3.

圖3 一階差分后自相關(guān)圖與偏相關(guān)圖Fig.3 Autocorrelation and partial correlation after the first order difference

圖3的ACF（自相關(guān)圖）顯示時間序列所有自相關(guān)函數(shù)都在置信區(qū)間內(nèi)且數(shù)值趨向于0,因此一階差分后其具有隨機性與平穩(wěn)性特征.進一步得知,圖ACF、PACF均截尾,可得到備選模型ARIMA（1,1,1）、ARIMA（1,1,0）.

對備選模型 ARIMA（1,1,1）、ARIMA（1,1,0）做參數(shù)p、d、q估計與假設(shè)檢驗,見表 1.

表1 參數(shù)估計與假設(shè)檢驗Tab.1 Parameter estimation and hypothesistesting

表 1顯示,ARIMA（1,1,1）的P值 ＜0.05,ARIMA（1,1,0）的P值 ＞0.05,且由 BIC最優(yōu)準(zhǔn)則得,ARIMA（1,1,1）為最優(yōu)模型.由表1得到產(chǎn)量預(yù)測模型

利用式（6）對2019和2020年的產(chǎn)量進行預(yù)測,得到預(yù)測產(chǎn)量分別為13 152.20萬t、13 160.79萬t.

2.4 模型檢驗與誤差分析

將數(shù)據(jù)分成實驗組1979—2016年與檢驗組2017—2018年,首先將實驗組1979—2016年的數(shù)據(jù)帶入ARIMA預(yù)測模型,得到實驗組1979—2016年產(chǎn)量實際值與擬合值的誤差,結(jié)果見表2.

表2 實驗組1979—2016年產(chǎn)量實際值與擬合值誤差表Tab.2 Error tablebetween actual and fitting values of annual output in experimental group

由表2計算得小麥實驗組1979—2016年實際值與擬合值平均誤差為5.29%,誤差較小,擬合程度較好,模型合適.觀察發(fā)現(xiàn),表2中1980年、1983年、1998年、2000年小麥產(chǎn)量實際值與擬合值差值較大,通過查閱資料并分析原因如下：1983年小麥產(chǎn)量實際值大于擬合值,該年種植面積無特別增減,但畝產(chǎn)量有明顯提升,所以該年小麥總產(chǎn)量增加.而1980年、1998年、2000年小麥產(chǎn)量實際值小于擬合值,其中1980年與1998年因秋冬連續(xù)干旱,11月份突然降溫,使小麥抗寒能力降低,而早春3—4月份的低溫光照不足,又使小麥抽穗期延遲,灌漿期縮短,所以這兩年產(chǎn)量減少；2000年干旱少雨,水資源緊張,種植面積減少,而且黃淮冬麥區(qū)小麥在收獲期間遭受暴雨,無法及時收割及晾曬,小麥質(zhì)量嚴(yán)重受損,所以2000年小麥產(chǎn)量也有所減少.

實驗組數(shù)據(jù)的誤差結(jié)果顯示ARIMA（1,1,1）模型合適,故進一步分析檢驗組2017—2018年產(chǎn)量實際值與預(yù)測值的誤差結(jié)果,見表3.

表3 檢驗組2017-2018年產(chǎn)量實際值與預(yù)測值誤差Tab.3 Error between actual annual output and predicted value of test group

由表3可知,檢驗組2017—2018年產(chǎn)量實際值與預(yù)測值平均誤差為1.65%.

綜上所述,ARIMA（1,1,1）模型的實驗組與檢驗組的結(jié)果顯示,擬合與預(yù)測的誤差較小,擬合效果較好,模型合適且支持應(yīng)用于糧食產(chǎn)量預(yù)測.

2.5 白噪聲檢驗

若殘差序列通過白噪聲檢驗[11],說明產(chǎn)量數(shù)據(jù)中的有用信息已經(jīng)被提取完畢,此時建?？梢越K止,ARIMA最優(yōu)模型可用來預(yù)測產(chǎn)量.ARIMA（1,1,1）模型殘差序列的自相關(guān)與偏自相關(guān)圖見圖4.

圖4 殘差序列的自相關(guān)與偏自相關(guān)圖Fig.4 Autocorrelation and partical autocorrelation of residual sequence

由圖4可知,殘差序列的自相關(guān)函數(shù)和偏自相關(guān)函數(shù)基本都在95%的置信區(qū)間內(nèi),說明殘差序列不存在自相關(guān).該模型通過白噪聲檢驗,可以用來預(yù)測短期中國小麥的產(chǎn)量.

3 中國小麥供需預(yù)測

3.1 進口量預(yù)測

在多個樣本數(shù)據(jù)中,進口量在供給量占一定比重,具有分析意義,所以我們考慮進口量的影響,并進行預(yù)測.從《中國統(tǒng)計年鑒》中選取1999—2018年中國小麥進口量、年需求量與人口數(shù),見表4.

表4 1999—2018年中國年小麥進口量、年需求量與人口Tab.4 Annual wheat import,annual consumption and population in China,1999—2018

利用上述同樣ARIMA方法,代入表4中進口量數(shù)據(jù),建立ARIMA（3,2,3）進口量預(yù)測模型,得到關(guān)系式

通過式（7）對2019和2020年的進口量進行預(yù)測,得到進口量實際值與擬合值趨勢圖,見圖5.

圖5 進口量實際值與擬合值趨勢Fig.5 Trend chart of actual and fitting values of import volume

由圖5可知,2019年和2020年中國小麥進口量預(yù)計分別為210.98萬t和121.15萬t.

供給量=產(chǎn)量+進口量

計算可得,2019年和2020年小麥供給量分別為13 363.18萬t和13 281.94萬t.

3.2 需求量預(yù)測

利用同樣ARIMA方法,代入表4中的需求量,建立ARIMA（3,2,2）需求預(yù)測模型,得到關(guān)系式

通過式（8）對2019和2020年的需求量進行預(yù)測,得到1999—2019需求量的實際值與擬合值,見表5.

表5 1999—2019需求量的實際值與擬合值Tab.5 Actual and fittingvalues of thedemand in 1999—2019

由表5可知,2019年和2020年需求量預(yù)測結(jié)果分別為12 683.65萬t和13 030.17萬t.

4 小結(jié)

本文建立了中國糧食供應(yīng)和需求預(yù)測模型,得到的預(yù)測值和觀測值在95%的置信區(qū)間內(nèi)擬合很好.然而,ARIMA模型也具有一定的局限性,它是依靠歷史的統(tǒng)計數(shù)據(jù)建立的數(shù)學(xué)模型,并未考慮國家重大政策改變和調(diào)整、環(huán)境變化、區(qū)域性差異等外部因素的影響.因此ARIMA模型僅適用于短期結(jié)果預(yù)測,預(yù)測未來2～3年的趨勢是可行的.預(yù)測時間過長會增大預(yù)測誤差,影響預(yù)測精度.因此需要及時了解市場行情、重大氣候變化,更新觀測值,重新擬合,確保預(yù)測的精度,或在此模型基礎(chǔ)上進行優(yōu)化,繼而得到適應(yīng)新情況的模型.

從中國小麥的供需數(shù)據(jù)與預(yù)測值來看,我國小麥供需接近平衡.長遠來看,應(yīng)該提高小麥畝產(chǎn)量,增加本國糧食生產(chǎn)總量,實現(xiàn)糧食安全,維護社會穩(wěn)定.