基于R語言的我國甘蔗產量ARIMA模型建立與預測分析

郝小玲，梁 春，龐新華,2，朱鵬錦,2，嚴 霖，黃 強,2,龍凌云

(1.廣西壯族自治區(qū)亞熱帶作物研究所，廣西 南寧 530001；2.廣西農墾甘蔗研究所，廣西 南寧 530001)

?

基于R語言的我國甘蔗產量ARIMA模型建立與預測分析

郝小玲1，梁 春1，龐新華1,2，朱鵬錦1,2，嚴 霖1，黃 強1,2,龍凌云1

(1.廣西壯族自治區(qū)亞熱帶作物研究所，廣西 南寧 530001；2.廣西農墾甘蔗研究所，廣西 南寧 530001)

摘要：采用時間序列預測法對我國甘蔗產量預測問題進行了研究。以國家統(tǒng)計局1965～2012年間我國歷年甘蔗產量統(tǒng)計數(shù)據為基礎，依據自回歸積分滑動平均(ARIMA)模型理論，運用R語言建立了ARIMA模型對我國未來甘蔗產量進行了預測。數(shù)據檢驗顯示：該模型擬合效果較好，預測精度較高，表明運用該模型對我國甘蔗產量的變化趨勢進行了分析及預測是可行的，具有較高的實際應用價值。分析了該模型所得預測結果，對我國甘蔗產業(yè)未來發(fā)展的可預見性趨勢與風險提供了應對建議。

關鍵詞：R語言；甘蔗產量；ARIMA模型；預測

1引言

我國是世界第三大產糖國和第一大原糖進口國，據聯(lián)合國糧農組織(FAO)2013年數(shù)據統(tǒng)計，我國產糖量約占全球食糖總產的8 %，原糖進口量約占貿易總量的12 %[1,2]。其中，蔗糖是我國食糖來源的主體。近年來由于我國甘蔗主產區(qū)生產技術水平較低，品種退化嚴重，綜合利用不完善，政策扶持力度不夠以及農民種蔗積極性降低等因素，我國甘蔗生產呈萎縮態(tài)勢[3,4]。因此，亟待對甘蔗生產未來發(fā)展趨勢進行預測，從而制定有效的對策，把握市場供求關系，進而保護相關企業(yè)及種植戶利益，這將對優(yōu)化我國甘蔗產業(yè)的綜合發(fā)展產生現(xiàn)實影響。通過研究我國甘蔗生產發(fā)展的變動規(guī)律并對其產量預測，對甘蔗產業(yè)發(fā)展的變化趨勢進行分析與歸納，為提前防范和應對產業(yè)風險提供更有針對性的依據，使甘蔗生產更適應市場經濟發(fā)展要求，將進一步推動我國甘蔗產業(yè)的健康發(fā)展。

2我國甘蔗產業(yè)發(fā)展現(xiàn)狀

甘蔗(Saccharum officinarum)是禾本科甘蔗屬植物，原產于南亞與東南亞地區(qū)，在33°N-33°S之間地區(qū)均有分布，22°N-22°S范圍內面積較為集中[5]。甘蔗是全球范圍內最主要的糖料作物，蔗糖產量約占食糖總產的80 %，主要種植國家包括巴西、印度、中國、泰國以及巴基斯坦等[6,7]。2014年我國甘蔗播種面積為176.05萬hm2，僅次于巴西與印度，位居世界第3。據中國糖業(yè)協(xié)會數(shù)據統(tǒng)計，我國2014/15榨季累計產糖量為1055.6萬t，其中蔗糖產量981.82萬t，占我國食糖總產的93.01 %[8]。

近年來中央一號文件多次涉及“糧棉油糖”四類關系國計民生的重要農產品，足見中央對糖料生產的重視。優(yōu)化國內糖料作物生產，確保食糖有效供給，對維護社會及經濟發(fā)展穩(wěn)定具有重要意義。甘蔗是不僅是我國最主要的糖料作物，也是優(yōu)良的能源作物之一。此外，甘蔗提供的各類產品還是醫(yī)療、化工以及食品加工業(yè)等領域的上游產品。由此可見，甘蔗產業(yè)已經成為推動我國甘蔗主產地區(qū)人民致富、經濟發(fā)展和稅收增長的重要產業(yè)，甘蔗生產對農民增收、食品供給安全和國家能源保障以及工業(yè)發(fā)展等都具有關鍵意義。

3模型預測理論與實現(xiàn)

目前國內利用數(shù)學模型對農產品產量進行預測的研究很多，主要有灰色系統(tǒng)預測模型、神經網絡、時間序列模型以及多元回歸模型等。方孝榮等通過建立灰色-馬爾科夫模型對浙江省名優(yōu)茶年產量進行預測[9]；陳玉佳等運用小波分析和BP神經網絡，建立加工番茄產量預測的小波神經網絡模型[10]；朱秀紅等通過對氣候因子和茶葉產量進行相關性分析，建立多元線性回歸模型對日照市茶葉產量進行預測[11]。

不同數(shù)學模型依據的理論基礎、特點以及局限性都各不相同。其中，灰色系統(tǒng)預測模型及神經網絡對對象產量的中長期趨勢預測較為準確，對短期波動的預測概率不高。多元回歸模型主要用來分析多個影響因素對研究對象產量的影響大小。而時間序列模型既考慮了時間序列數(shù)據之間的依存性，又考慮到了隨機波動的干擾性影響，預測的準確度比較高，是近年應用比較廣泛的方法之一[12]。

3.1時間序列預測理論

時間序列預測方法(Time Series Forecasting Method)是通過歷史資料數(shù)據揭示現(xiàn)象隨時間變化的發(fā)展過程、方向或趨勢，將這種規(guī)律進行外推或延伸，從而對該現(xiàn)象的未來做出預測。傳統(tǒng)的時間序列分析主要是確定性的時間序列分析方法，對于不平穩(wěn)的時間序列，通常需要轉換成平穩(wěn)的時間序列[12]。以此為基礎，Box和Jenkins于20世紀70年代初提出了自回歸積分滑動平均模型(Autoregressive Integrated Moving Average Model，ARIMA模型)，使預測的精確度大大提高[13]。它是在將非平穩(wěn)時間序列轉化為平穩(wěn)時間序列的前提下，根據對因變量和因變量的滯后值以及隨機誤差項的現(xiàn)值和滯后值進行回歸所建立的模型。

ARIMA模型的基本思想在于，將預測對象隨時間推移而形成的數(shù)據序列視為一個隨機序列，用一定的數(shù)學模型來擬合該序列。通過識別這個數(shù)學模型，就可以利用時間序列的過去值和當前值來預測未來值[13,14]。

ARIMA模型表示一個數(shù)據序列在某一時刻的值不僅與以前的自身值有關，而且還與其歷史時刻的擾動項存在一定的依存關系。如果以q代表模型的滑動平均項數(shù)，p代表自回歸的階數(shù)，d為時間序列成為平穩(wěn)之前必須做差分的次數(shù)，則ARIMA(p，d，q)模型的實質可描述為先對非平穩(wěn)時間數(shù)據進行d次差分處理得到新的平穩(wěn)的數(shù)據序列Xt，將Xt擬合為自回歸滑動平均(Auto Regressive Moving Average，ARMA)模型ARMA(p，q)，再將原d次差分還原，從而得到預測數(shù)據。由此，可得ARIMA(p，d，q)的一般表達式：

其中，L為滯后算子，φi為模型自回歸部分參數(shù)，θi為滑動平均部分參數(shù)，εt為誤差項。

3.2模型的確定

ARIMA 模型的基礎是以平穩(wěn)時間序列為定義的。在某些條件下，需要通過考慮數(shù)據之間的相關性來創(chuàng)建更好的預測模型。但是，ARIMA模型包含一個確定的統(tǒng)計模型用于處理時間序列的不規(guī)則部分，它也允許不規(guī)則部分可以自相關。

(1)根據數(shù)據的時序圖、自相關函數(shù)(Autocorrelation Function，ACF)圖和偏自相關函數(shù)(Partial Autocorrelation Function，PACF)圖，觀察其方差、趨勢及其季節(jié)性變化規(guī)律，識別該序列的平穩(wěn)性。

(2)如果數(shù)據序列是非平穩(wěn)的，則需對數(shù)據進行差分或滑動平均法處理，直到所得數(shù)據可以通過平穩(wěn)性檢驗。

(3)根據時間序列模型的識別規(guī)則，建立相應的模型。若數(shù)據序列的PACF是截尾的，而ACF是拖尾的，則可判斷此序列符合自回歸(Auto Regressive，AR)模型；若數(shù)據序列的PACF是拖尾的，而ACF是截尾的，則可判斷此序列符合滑動平均(Moving Average，MA)模型；若數(shù)據序列的PACF和ACF均是拖尾的，則此序列符合ARMA模型。

(4)模型定階和參數(shù)估計。根據時間序列模型的赤池信息量(Akaike Information Criterion，AIC)準則、貝葉斯信息量(Bayesian Information Criterion，BIC)準則等規(guī)則，確定模型的階。然后進一步估計暫定的模型參數(shù)，檢驗是否具有統(tǒng)計意義。

(5)模型檢驗。檢驗假設模型殘差的ACF和PACF在早期或季節(jié)性延遲點(如果是季節(jié)性序列)處是否不大于置信區(qū)間，同時檢驗所得模型的殘差序列是否為白噪聲。

(6)通過上述模型識別和檢驗，確定最優(yōu)模型并利用模型進行預測。

3.3實現(xiàn)工具

研究采用的實現(xiàn)工具是R語言。R語言是一種自由軟件編程語言與操作環(huán)境，主要用于統(tǒng)計分析、繪圖以及數(shù)據挖掘。R語言最早由新西蘭奧克蘭大學的Ross Ihaka與Robert Gentleman開發(fā)，相比其他統(tǒng)計學或數(shù)學專用的編程語言有著更強的面向對象功能。在R語言中有許多擴展包(Packages)可以增強其廣泛使用性，從而實現(xiàn)線性和非線性建模，傳統(tǒng)統(tǒng)計檢驗，時間序列分析以及聚類等諸多統(tǒng)計和圖形功能。R語言在國外被廣泛使用，包括Google、Facebook以及LinkedIn在內的諸多大型企業(yè)均使用它進行數(shù)據分析。

4數(shù)據處理與模型建立

本文根據國家統(tǒng)計局1965～2012年間我國歷年甘蔗產量統(tǒng)計數(shù)據建立模型。

4.1模型的識別

以原始數(shù)據利用R語言ts及plot指令建立時序圖(圖1)。由ARIMA理論可知，建立該模型的理想化要求為時間序列必須是零均值、正態(tài)及平穩(wěn)的隨機過程。根據圖形可以辨別數(shù)據表現(xiàn)出明顯的上升趨勢，初步判斷數(shù)據需要進行差分處理使序列平穩(wěn)化。

圖1　1965～2012年我國甘蔗產量時序

使用diff指令對數(shù)據進行差分處理，同時使用plot指令繪制差分后序列圖(圖2)。根據圖形可判斷，序列在進行一次及二次差分處理后均較為平穩(wěn)。根據數(shù)據最優(yōu)化及準確性原則，在一次差分即可使序列平穩(wěn)化時，通常優(yōu)先嘗試一次差分進行處理。運用R語言tseries擴展包的adf.test指令對一次差分后數(shù)據進行ADF單位根檢驗，得p值為0.01，拒絕5 %顯著水平下原假設，提示序列中不存在單位根，從而確認一次差分后的數(shù)據已趨于平穩(wěn)。由此初步確定數(shù)據序列需進行一次差分處理。

圖2　一次差分及二次差分序列

4.2模型選擇與參數(shù)估計

使用acf及pacf指令繪制一次差分處理后序列的自相關函數(shù)(ACF)圖和偏自相關函數(shù)(PACF)圖(圖3)。通過觀察圖像可得，該序列的ACF圖截尾；PACF圖拖尾，并在五階后收攏至置信區(qū)間。初步判斷可建立MA(5)或MA(6)模型。由此使用arima指令建立相應的模型，并根據AIC準則對兩個模型進行評估。該指令將自動計算模型的AIC值，非常利于進行模型的比較。由此得到MA(5)模型的AIC值為730.59，MA(6)模型的AIC值為719.74。根據AIC值越小模型越優(yōu)的理想化原則，可以初步判定MA(6)模型相對于MA(5)模型是更優(yōu)的預測模型。

圖3　一次差分序列ACF(左)與PACF(右)

4.3模型的檢驗

圖4　tsdiag指令檢驗結果

使用tsdiag指令對選定的MA(6)模型進行檢驗，得到對應模型的殘差時序圖、殘差ACF圖以及Ljung-Box檢驗p值散點圖(圖4)。根據時序圖可以確認，殘差沒有表現(xiàn)出明顯的周期性，ACF圖顯示在滯后1～15階均未超出置信邊界，表明模型殘差沒有明顯的自相關性。Ljung-Box檢驗表明隨機干擾項是白噪聲。由此可以確定時間序列模型ARIMA(6，1，0)通過檢驗，最終得到預測模型。

4.4模型評價

使用plot及fitted指令繪制序列實際值與模型擬合值曲線圖(圖5)，可觀察到模型的擬合度較高。

圖5　序列實際值與模型擬合值曲線

所得2008～2012年實際值與擬合值對比見表1。

表1　2008～2012年實際值與模型擬合值對比

4.5數(shù)據預測

使用R語言forecast擴展包的forecast.Arima指令，運用該模型對我國2013～2014年甘蔗產量進行估測，得到兩年產量分別為13158.92萬t及13451.02萬t，與實際值誤差分別為2.64 %及7.08 %。考慮到近兩年我國甘蔗主產區(qū)產量受行情變動及自然災害影響波動較大，該誤差處于合理區(qū)間，可以認為該模型用于預測未來產量是相對準確的。運用該模型對我國2016～2036年甘蔗產量進行99.5 %置信度預測，同時繪制相關圖表見圖6。

圖6　2016～2016年甘蔗產量預測

5結語

通過對模型的驗證表明，運用ARIMA模型對我國甘蔗產量進行預測是較為準確的，對于時間序列預測方法而言，由于各種因素的相對穩(wěn)定性，在一個較短時期內，可以大致認為相關因素對預測對象的影響及其自身的變化趨勢是規(guī)律性的，利用歷史數(shù)據進行預測就可以保證一定的預測精度[13]。因此，將該模型用于我國甘蔗產量的中短期實際預測是可行的。

由預測數(shù)據可得，未來20年間我國甘蔗產量將會持續(xù)波動上升，但整體上升態(tài)勢將會趨緩，預測年均增長率約為1.2 %。其中，2025年甘蔗產量預計將達到14972.19萬t，2035年將達到15871.17萬t。這是符合我國甘蔗生產周期以及經濟政策發(fā)展趨勢的。我國是世界食糖消費大國，然而當前人均食糖消費量約為10～12kg，遠低于24kg的世界平均水平[1,14]。隨著未來我國人口增長、居民消費水平提高以及食品加工等相關產業(yè)的發(fā)展，從中長期來看，我國食糖需求量將會進一步提高。在我國食糖供給以蔗糖為主的前提下，放緩的甘蔗生產趨勢勢必會進一步拉大現(xiàn)已存在的食糖供需缺口，對我國的食糖自給率造成沖擊[3,4]。

為維護我國食品供給安全與經濟穩(wěn)定，促進市場健康發(fā)展，需要政府決策者與相關企業(yè)從統(tǒng)籌經濟效益、社會效益與生態(tài)效益的角度出發(fā)，在穩(wěn)定價格與保障蔗農利益的前提下，通過政策引導、良種推廣、機械化與規(guī)模化生產等多渠道措施的整合，最終推動我國甘蔗產業(yè)的可持續(xù)發(fā)展。

參考文獻：

[1]廖平偉，張華，羅俊，等.我國甘蔗生產現(xiàn)狀及競爭力分析[J].中國糖料，2010(4):44～45.

[2]聯(lián)合國糧食及農業(yè)組織. FAOSTAT[DB/OL].[2016-04-01]http://faostat3.fao.org.

[3]韋小蕾.廣西甘蔗產業(yè)化現(xiàn)狀研究[J].中國市場，2014(46)：44～46.

[4]趙建屹. 甘蔗種業(yè)發(fā)展現(xiàn)狀與對策研究[D].福州：福建農林大學,2014.

[5] 黃中艷. 中國甘蔗氣候類型和特點的客觀分析[J]. 作物雜志, 2009(02): 21～25.

[6]李楊瑞，楊麗濤.20世紀90年代以來我國甘蔗產業(yè)和科技的新發(fā)展[J].西南農業(yè)學報，2009(5)：1469～1476.

[7]鄧軍,蔡曉琳,付思明,等. 中國蔗糖產業(yè)布局及發(fā)展對策[J]. 甘蔗糖業(yè),2011(1):57～60.

[8]劉曉雪,張宸,宋杰. 2014/15榨季國內外食糖市場回顧與2015/16榨季展望[J]. 農業(yè)展望,2015(10):12～17.

[9]方孝榮,丁希斌,李曉麗. 基于灰色馬爾柯夫鏈模型的浙江省名優(yōu)茶產量預測[J]. 農機化研究,2014(7):18～21.

[10]陳玉佳,姜波. 基于小波神經網絡的加工番茄產量預測模型[J]. 深圳大學學報(理工版),2015(5):546～550.

[11] 朱秀紅,鄭美琴,姚文軍,于懷征. 基于SPSS的日照市茶葉產量預測模型的建立[J]. 河南農業(yè)科學,2010(7):31～33.

[12]張伶燕,葛翔. 時間序列模型在我國牛肉產量預測中的應用[J]. 中國畜牧雜志,2008(7):42～45.

[13]George E. P. Box, Gwilym M. Jenkins, Gregory C. Reinsel. Time Series Analysis: Forecasting and Control [M]. Wiley, 2015.

[14]溫嵐,陳基權,戴志剛,龔友才,劉倩,李楠,粟建光. 長蒴黃麻產葉量的多元回歸與偏相關的R語言分析[J]. 作物雜志,2013(1):49～53.

[15]劉曉雪,王沈南,鄭傳芳. 2015—2030年中國食糖消費量預測和供需缺口分析[J]. 農業(yè)展望,2013(2):71～75.

收稿日期：2016-04-28

基金項目：廣西科學研究與技術開發(fā)計劃項目(編號：桂科攻1598006-1-1B)；南寧市科學研究與技術開發(fā)計劃項目(編號：科技攻關20152063)；廣西農墾科學研究與技術開發(fā)計劃項目(編號：桂墾攻科201502)；廣西區(qū)直公益性科研院所基本科研業(yè)務經費專項資助(編號：桂熱研201403)

作者簡介：郝小玲(1988—)，女，碩士，主要從事熱帶亞熱帶作物栽培與育種方面的研究工作。 通訊作者：龐新華(1968—)，男，碩士，高級農藝師，主要從事熱帶亞熱帶作物栽培與育種方面的研究工作。

中圖分類號：C32

文獻標識碼：A

文章編號：1674-9944(2016)12-0257-04

Establishment of ARIMA Model and Forecast for Sugar Cane Production in China Based on R Language

Hao Xiaoling1, Liang Chun1,Pang Xinhua1,2, Zhu Pengjin1,2, Yan Lin1,2, Huang Qiang1,2，Long Lingyun1

(1.GuangxiInstituteofSubtropicalCrops,Naning530001;2.SugarcaneResearchInstitute,GuangxiAgriculturalReclamationBureau,Naning530001)

Abstract:This paper studied the prediction of the yield of sugar cane in China with time series forecasting method. Based on statistical data of sugar cane annual production during 1965-2012 years provided by National Bureau of Statistics, this paper established an autoregressive integrated moving average (ARIMA) model with R Language to predict the future of sugar cane production in China, on the basis of ARIMA model theory. Data test showed that the model fitted the data well, with good precision. The test result also indicated that utilizing the model to forecast and analyze the trend for sugar cane production of China is feasible, with a high practical value. Through the analysis of the prediction result of the model, this paper provided advice on trend forecasting and risk response for the future development of sugar cane industry of China.

Key words:R Language；sugarcane yields；ARIMA modelestablishment；forcast