基于R语言的我国甘蔗产量ARIMA模型建立与预测分析

郝小玲，梁 春，庞新华,2，朱鹏锦,2，严 霖，黄 强,2,龙凌云

(1.广西壮族自治区亚热带作物研究所，广西 南宁 530001；2.广西农垦甘蔗研究所，广西 南宁 530001)



基于R语言的我国甘蔗产量ARIMA模型建立与预测分析

郝小玲1，梁 春1，庞新华1,2，朱鹏锦1,2，严 霖1，黄 强1,2,龙凌云1

(1.广西壮族自治区亚热带作物研究所，广西 南宁 530001；2.广西农垦甘蔗研究所，广西 南宁 530001)

摘要：采用时间序列预测法对我国甘蔗产量预测问题进行了研究。以国家统计局1965～2012年间我国历年甘蔗产量统计数据为基础，依据自回归积分滑动平均(ARIMA)模型理论，运用R语言建立了ARIMA模型对我国未来甘蔗产量进行了预测。数据检验显示：该模型拟合效果较好，预测精度较高，表明运用该模型对我国甘蔗产量的变化趋势进行了分析及预测是可行的，具有较高的实际应用价值。分析了该模型所得预测结果，对我国甘蔗产业未来发展的可预见性趋势与风险提供了应对建议。

关键词：R语言；甘蔗产量；ARIMA模型；预测

1引言

我国是世界第三大产糖国和第一大原糖进口国，据联合国粮农组织(FAO)2013年数据统计，我国产糖量约占全球食糖总产的8 %，原糖进口量约占贸易总量的12 %[1,2]。其中，蔗糖是我国食糖来源的主体。近年来由于我国甘蔗主产区生产技术水平较低，品种退化严重，综合利用不完善，政策扶持力度不够以及农民种蔗积极性降低等因素，我国甘蔗生产呈萎缩态势[3,4]。因此，亟待对甘蔗生产未来发展趋势进行预测，从而制定有效的对策，把握市场供求关系，进而保护相关企业及种植户利益，这将对优化我国甘蔗产业的综合发展产生现实影响。通过研究我国甘蔗生产发展的变动规律并对其产量预测，对甘蔗产业发展的变化趋势进行分析与归纳，为提前防范和应对产业风险提供更有针对性的依据，使甘蔗生产更适应市场经济发展要求，将进一步推动我国甘蔗产业的健康发展。

2我国甘蔗产业发展现状

甘蔗(Saccharum officinarum)是禾本科甘蔗属植物，原产于南亚与东南亚地区，在33°N-33°S之间地区均有分布，22°N-22°S范围内面积较为集中[5]。甘蔗是全球范围内最主要的糖料作物，蔗糖产量约占食糖总产的80 %，主要种植国家包括巴西、印度、中国、泰国以及巴基斯坦等[6,7]。2014年我国甘蔗播种面积为176.05万hm2，仅次于巴西与印度，位居世界第3。据中国糖业协会数据统计，我国2014/15榨季累计产糖量为1055.6万t，其中蔗糖产量981.82万t，占我国食糖总产的93.01 %[8]。

近年来中央一号文件多次涉及“粮棉油糖”四类关系国计民生的重要农产品，足见中央对糖料生产的重视。优化国内糖料作物生产，确保食糖有效供给，对维护社会及经济发展稳定具有重要意义。甘蔗是不仅是我国最主要的糖料作物，也是优良的能源作物之一。此外，甘蔗提供的各类产品还是医疗、化工以及食品加工业等领域的上游产品。由此可见，甘蔗产业已经成为推动我国甘蔗主产地区人民致富、经济发展和税收增长的重要产业，甘蔗生产对农民增收、食品供给安全和国家能源保障以及工业发展等都具有关键意义。

3模型预测理论与实现

目前国内利用数学模型对农产品产量进行预测的研究很多，主要有灰色系统预测模型、神经网络、时间序列模型以及多元回归模型等。方孝荣等通过建立灰色-马尔科夫模型对浙江省名优茶年产量进行预测[9]；陈玉佳等运用小波分析和BP神经网络，建立加工番茄产量预测的小波神经网络模型[10]；朱秀红等通过对气候因子和茶叶产量进行相关性分析，建立多元线性回归模型对日照市茶叶产量进行预测[11]。

不同数学模型依据的理论基础、特点以及局限性都各不相同。其中，灰色系统预测模型及神经网络对对象产量的中长期趋势预测较为准确，对短期波动的预测概率不高。多元回归模型主要用来分析多个影响因素对研究对象产量的影响大小。而时间序列模型既考虑了时间序列数据之间的依存性，又考虑到了随机波动的干扰性影响，预测的准确度比较高，是近年应用比较广泛的方法之一[12]。

3.1时间序列预测理论

时间序列预测方法(Time Series Forecasting Method)是通过历史资料数据揭示现象随时间变化的发展过程、方向或趋势，将这种规律进行外推或延伸，从而对该现象的未来做出预测。传统的时间序列分析主要是确定性的时间序列分析方法，对于不平稳的时间序列，通常需要转换成平稳的时间序列[12]。以此为基础，Box和Jenkins于20世纪70年代初提出了自回归积分滑动平均模型(Autoregressive Integrated Moving Average Model，ARIMA模型)，使预测的精确度大大提高[13]。它是在将非平稳时间序列转化为平稳时间序列的前提下，根据对因变量和因变量的滞后值以及随机误差项的现值和滞后值进行回归所建立的模型。

ARIMA模型的基本思想在于，将预测对象随时间推移而形成的数据序列视为一个随机序列，用一定的数学模型来拟合该序列。通过识别这个数学模型，就可以利用时间序列的过去值和当前值来预测未来值[13,14]。

ARIMA模型表示一个数据序列在某一时刻的值不仅与以前的自身值有关，而且还与其历史时刻的扰动项存在一定的依存关系。如果以q代表模型的滑动平均项数，p代表自回归的阶数，d为时间序列成为平稳之前必须做差分的次数，则ARIMA(p，d，q)模型的实质可描述为先对非平稳时间数据进行d次差分处理得到新的平稳的数据序列Xt，将Xt拟合为自回归滑动平均(Auto Regressive Moving Average，ARMA)模型ARMA(p，q)，再将原d次差分还原，从而得到预测数据。由此，可得ARIMA(p，d，q)的一般表达式：

其中，L为滞后算子，φi为模型自回归部分参数，θi为滑动平均部分参数，εt为误差项。

3.2模型的确定

ARIMA 模型的基础是以平稳时间序列为定义的。在某些条件下，需要通过考虑数据之间的相关性来创建更好的预测模型。但是，ARIMA模型包含一个确定的统计模型用于处理时间序列的不规则部分，它也允许不规则部分可以自相关。

(1)根据数据的时序图、自相关函数(Autocorrelation Function，ACF)图和偏自相关函数(Partial Autocorrelation Function，PACF)图，观察其方差、趋势及其季节性变化规律，识别该序列的平稳性。

(2)如果数据序列是非平稳的，则需对数据进行差分或滑动平均法处理，直到所得数据可以通过平稳性检验。

(3)根据时间序列模型的识别规则，建立相应的模型。若数据序列的PACF是截尾的，而ACF是拖尾的，则可判断此序列符合自回归(Auto Regressive，AR)模型；若数据序列的PACF是拖尾的，而ACF是截尾的，则可判断此序列符合滑动平均(Moving Average，MA)模型；若数据序列的PACF和ACF均是拖尾的，则此序列符合ARMA模型。

(4)模型定阶和参数估计。根据时间序列模型的赤池信息量(Akaike Information Criterion，AIC)准则、贝叶斯信息量(Bayesian Information Criterion，BIC)准则等规则，确定模型的阶。然后进一步估计暂定的模型参数，检验是否具有统计意义。

(5)模型检验。检验假设模型残差的ACF和PACF在早期或季节性延迟点(如果是季节性序列)处是否不大于置信区间，同时检验所得模型的残差序列是否为白噪声。

(6)通过上述模型识别和检验，确定最优模型并利用模型进行预测。

3.3实现工具

研究采用的实现工具是R语言。R语言是一种自由软件编程语言与操作环境，主要用于统计分析、绘图以及数据挖掘。R语言最早由新西兰奥克兰大学的Ross Ihaka与Robert Gentleman开发，相比其他统计学或数学专用的编程语言有着更强的面向对象功能。在R语言中有许多扩展包(Packages)可以增强其广泛使用性，从而实现线性和非线性建模，传统统计检验，时间序列分析以及聚类等诸多统计和图形功能。R语言在国外被广泛使用，包括Google、Facebook以及LinkedIn在内的诸多大型企业均使用它进行数据分析。

4数据处理与模型建立

本文根据国家统计局1965～2012年间我国历年甘蔗产量统计数据建立模型。

4.1模型的识别

以原始数据利用R语言ts及plot指令建立时序图(图1)。由ARIMA理论可知，建立该模型的理想化要求为时间序列必须是零均值、正态及平稳的随机过程。根据图形可以辨别数据表现出明显的上升趋势，初步判断数据需要进行差分处理使序列平稳化。

图1　1965～2012年我国甘蔗产量时序

使用diff指令对数据进行差分处理，同时使用plot指令绘制差分后序列图(图2)。根据图形可判断，序列在进行一次及二次差分处理后均较为平稳。根据数据最优化及准确性原则，在一次差分即可使序列平稳化时，通常优先尝试一次差分进行处理。运用R语言tseries扩展包的adf.test指令对一次差分后数据进行ADF单位根检验，得p值为0.01，拒绝5 %显著水平下原假设，提示序列中不存在单位根，从而确认一次差分后的数据已趋于平稳。由此初步确定数据序列需进行一次差分处理。

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

4.2模型选择与参数估计

使用acf及pacf指令绘制一次差分处理后序列的自相关函数(ACF)图和偏自相关函数(PACF)图(图3)。通过观察图像可得，该序列的ACF图截尾；PACF图拖尾，并在五阶后收拢至置信区间。初步判断可建立MA(5)或MA(6)模型。由此使用arima指令建立相应的模型，并根据AIC准则对两个模型进行评估。该指令将自动计算模型的AIC值，非常利于进行模型的比较。由此得到MA(5)模型的AIC值为730.59，MA(6)模型的AIC值为719.74。根据AIC值越小模型越优的理想化原则，可以初步判定MA(6)模型相对于MA(5)模型是更优的预测模型。

图3　一次差分序列ACF(左)与PACF(右)

4.3模型的检验

图4　tsdiag指令检验结果

使用tsdiag指令对选定的MA(6)模型进行检验，得到对应模型的残差时序图、残差ACF图以及Ljung-Box检验p值散点图(图4)。根据时序图可以确认，残差没有表现出明显的周期性，ACF图显示在滞后1～15阶均未超出置信边界，表明模型残差没有明显的自相关性。Ljung-Box检验表明随机干扰项是白噪声。由此可以确定时间序列模型ARIMA(6，1，0)通过检验，最终得到预测模型。

4.4模型评价

使用plot及fitted指令绘制序列实际值与模型拟合值曲线图(图5)，可观察到模型的拟合度较高。

图5　序列实际值与模型拟合值曲线

所得2008～2012年实际值与拟合值对比见表1。

表1　2008～2012年实际值与模型拟合值对比

4.5数据预测

使用R语言forecast扩展包的forecast.Arima指令，运用该模型对我国2013～2014年甘蔗产量进行估测，得到两年产量分别为13158.92万t及13451.02万t，与实际值误差分别为2.64 %及7.08 %。考虑到近两年我国甘蔗主产区产量受行情变动及自然灾害影响波动较大，该误差处于合理区间，可以认为该模型用于预测未来产量是相对准确的。运用该模型对我国2016～2036年甘蔗产量进行99.5 %置信度预测，同时绘制相关图表见图6。

图6　2016～2016年甘蔗产量预测

5结语

通过对模型的验证表明，运用ARIMA模型对我国甘蔗产量进行预测是较为准确的，对于时间序列预测方法而言，由于各种因素的相对稳定性，在一个较短时期内，可以大致认为相关因素对预测对象的影响及其自身的变化趋势是规律性的，利用历史数据进行预测就可以保证一定的预测精度[13]。因此，将该模型用于我国甘蔗产量的中短期实际预测是可行的。

由预测数据可得，未来20年间我国甘蔗产量将会持续波动上升，但整体上升态势将会趋缓，预测年均增长率约为1.2 %。其中，2025年甘蔗产量预计将达到14972.19万t，2035年将达到15871.17万t。这是符合我国甘蔗生产周期以及经济政策发展趋势的。我国是世界食糖消费大国，然而当前人均食糖消费量约为10～12kg，远低于24kg的世界平均水平[1,14]。随着未来我国人口增长、居民消费水平提高以及食品加工等相关产业的发展，从中长期来看，我国食糖需求量将会进一步提高。在我国食糖供给以蔗糖为主的前提下，放缓的甘蔗生产趋势势必会进一步拉大现已存在的食糖供需缺口，对我国的食糖自给率造成冲击[3,4]。

为维护我国食品供给安全与经济稳定，促进市场健康发展，需要政府决策者与相关企业从统筹经济效益、社会效益与生态效益的角度出发，在稳定价格与保障蔗农利益的前提下，通过政策引导、良种推广、机械化与规模化生产等多渠道措施的整合，最终推动我国甘蔗产业的可持续发展。

参考文献：

[1]廖平伟，张华，罗俊，等.我国甘蔗生产现状及竞争力分析[J].中国糖料，2010(4):44～45.

[2]联合国粮食及农业组织. FAOSTAT[DB/OL].[2016-04-01]http://faostat3.fao.org.

[3]韦小蕾.广西甘蔗产业化现状研究[J].中国市场，2014(46)：44～46.

[4]赵建屹. 甘蔗种业发展现状与对策研究[D].福州：福建农林大学,2014.

[5] 黄中艳. 中国甘蔗气候类型和特点的客观分析[J]. 作物杂志, 2009(02): 21～25.

[6]李杨瑞，杨丽涛.20世纪90年代以来我国甘蔗产业和科技的新发展[J].西南农业学报，2009(5)：1469～1476.

[7]邓军,蔡晓琳,付思明,等. 中国蔗糖产业布局及发展对策[J]. 甘蔗糖业,2011(1):57～60.

[8]刘晓雪,张宸,宋杰. 2014/15榨季国内外食糖市场回顾与2015/16榨季展望[J]. 农业展望,2015(10):12～17.

[9]方孝荣,丁希斌,李晓丽. 基于灰色马尔柯夫链模型的浙江省名优茶产量预测[J]. 农机化研究,2014(7):18～21.

[10]陈玉佳,姜波. 基于小波神经网络的加工番茄产量预测模型[J]. 深圳大学学报(理工版),2015(5):546～550.

[11] 朱秀红,郑美琴,姚文军,于怀征. 基于SPSS的日照市茶叶产量预测模型的建立[J]. 河南农业科学,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]. 农业展望,2013(2):71～75.

收稿日期：2016-04-28

基金项目：广西科学研究与技术开发计划项目(编号：桂科攻1598006-1-1B)；南宁市科学研究与技术开发计划项目(编号：科技攻关20152063)；广西农垦科学研究与技术开发计划项目(编号：桂垦攻科201502)；广西区直公益性科研院所基本科研业务经费专项资助(编号：桂热研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