R软件在“时间序列分析”教学中的应用——以ARIMA模型为例

邓 新 田春雨 葛梅梅 相旭东

1.滁州学院数学与金融学院 安徽滁州 239000；2.中国电子科技集团公司第五十八研究所 江苏南京 210000

“时间序列分析”课程是运用时间序列分析的基本理论与方法，研究随着时间的变化，某些事物发生、发展的过程，以寻找其发展变化的规律，并预测其未来走势的一门方法论科学，是统计学专业的必修课程之一。该课程旨在培养学生运用时间序列分析方法分析、解决实际问题的能力。随着数字经济时代的到来，时间序列分析方法已成为统计学、经济学等相关专业人才必备的数据分析方法之一。

为激发兴趣，培养学生理论联系实际的能力，将统计软件引入时间序列分析教学中尤其重要。目前，关于R软件在统计学专业课程教学中的应用，已有许多学者做了相关研究。2011年，闫昭辉通过聚类分析、主成分分析等方面的具体实际应用阐述了该软件在多元统计分析教学中的应用；2015年，张毅宁论述了线性、广义线性、非线性三类回归模型的具体应用，以此研究了R软件在教育统计中的应用；2018年，卢玉桂将R软件引入抽样调查课程，并结合实际案例实现了整群抽样的样本抽取和总体参数估计；2019年，为培养学生处理大数据的能力，刘君娥以因子分析模型为例，探讨了R软件在统计学教学中的应用；2021年，陈超和潘海燕以R软件辅助医学统计学教学，并以t检验实际案例做了进一步探究。本研究将以ARIMA模型为例，详细说明R软件在时间序列分析教学实践中的应用。

R语言已拥有多个用于时间序列分析的程序包。利用这些程序包，能够完成时间序列的一系列分析工作：生成时间序列、绘制图形、模型识别、参数估计、模型预测。除了R软件，EViews和SAS也常被用于时间序列分析。相比较于这两款统计软件，R软件具有如下优势：第一，与SAS这样需要昂贵的购置费及维护费的商业软件相比，它是完全免费的，通过官方网站可以免费下载安装最新版本。第二，它占用内存非常小，编程代码简单明了，绘图清晰功能强大。第三，它拥有全球所有优秀的统计软件包，并向用户免费开放。基于以上优点，在实践教学过程中，不仅可以直接调用R软件函数命令，还可以利用R语言编写自己的程序。这既能满足教师在教学过程中的需求，也能提高学生的编程能力。

一、ARIMA模型及建模步骤

ARIMA(Auto Regressive Integrated Moving Average)模型，全称是差分自回归移动平均模型，是由博克斯-詹金斯提出的一种专门用于非平稳时间序列分析和预测的方法。ARIMA(p，d，q)模型基本结构如下：

其中∇=(1-)，Φ()=1--…-为平稳可逆的ARMA(p，q)模型的自回归系数多项式，Θ()=1--…-为平稳可逆的ARMA(p，q)模型的移动平均系数多项式。

使用ARIMA模型建模的具体过程：首先，对获取的观察值序列进行平稳性检验，若通过检验，进入下一步，否则，利用差分运算实现平稳；其次，进行白噪声检验，若检验通过，分析结束，否则，拟合ARIMA模型并开展模型检验；最后，模型检验通过，分析结束，否则，重新拟合模型。

二、案例分析——滁州市GDP预测分析

GDP，即国内生产总值，是对一个国家、一个地区的总体经济状况进行衡量的重要指标。滁州市地处安徽省的东部，是南京与合肥都市圈的重要成员。《滁州市统计年鉴》指出，在2019年，滁州市以2909.06亿元的GDP总量，居于安徽省第三名，同时滁州市的GDP以9.7%的增速排名安徽省第一。滁州经济如此惊人的成绩和增速，使其备受关注。通过滁州市2020年统计年鉴和《滁州市2020年国民经济和社会发展统计公报》，获取了1990—2020年滁州市GDP数据(单位：亿元)：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。将该数据分为1990—2018年和2019—2020年两段，前者用于建立模型，后者用于检验模型预测效果。

实验目的：利用ARIMA模型对滁州市GDP进行分析及预测。

实验环境：Window 10，R 4.0.3。

实验步骤：

(一)创建数据集

生成时间序列(对原始数据进行自然对数处理)

a<-read.table("滁州市GDP.csv"，sep="，"，header=T)#读入数据

lnGDP<-ts(log(a$GDP)，start=1990)#建立时间序列{lnGDP}

(二)平稳性检验

对lnGDP序列进行单位根检验：adf.test(lnGDP)，此时需要先下载安装tseries程序包，并用library(tseries)载入。单位根检验结果如下。

Augmented Dickey-Fuller Test

data：lnGDP

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

alternative hypothesis：stationary

根据该检验的P值大于显著性水平=005，说明该序列是非平稳的。

(三)平稳化处理(一阶差分)

对序列进行一阶差分：lnGDP.dif<-diff(lnGDP)，绘制一阶差分序列图：plot(lnGDP.dif)并进行单位根检验：adf.test(lnGDP.dif)。综合一阶差分时序图和单位根检验结果(P值为0.04664)，一阶差分序列是平稳的。

(四)白噪声检验

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

运行结果如下。

Box-Ljung test

data：lnGDP.dif

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

检验结果显示，LB统计量的P值小于显著性水平=005，说明该一阶差分序列属于非白噪声序列。

(五)模型识别

根据序列平稳化过程可知，=1。为确定p和q的值，分别利用函数acf和pacf绘制一阶差分序列的自相关系数(ACF)图和偏自相关系数(PACF)图。根据图1，自相关系数和偏自相关系数都为1阶截尾，所以初步选定模型阶数为(0，1)和(1，0)。

图1 一阶差分自相关系数图和偏自相关系数图

(六)模型优化和参数估计

分别利用ARIMA(1，1，0)和ARIMA(0，1，1)拟合1990—2018年滁州市ln(GDP)序列，其代码及运行结果如下。

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

根据AIC准则选择ARIMA(1，1，0)，确定模型口径为：

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

等价为=18559-1-08559-2+，～(0，000655)。

(七)模型检验

调用aTSA包中的ts.diag函数，对建立的ARIMA(1，1，0)模型进行显著性检验。根据白噪声检验结果(左下图)，各阶延迟下白噪声检验统计量的P值都显著大于显著性水平，因此认为该拟合模型的残差序列属于白噪声序列，即该拟合模型显著成立。

图2 拟合模型显著性检验图

(八)模型预测

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))

调用forecast函数，运行上述代码预测2019—2023年滁州市GDP(单位：亿元)：2894.83，3179.79，3445.85，3691.17，3914.96。经计算，ARIMA(1，1，0)模型的平均相对误差为5.36%，较小，说明该模型预测精度较高。

结语

ARIMA模型是处理非平稳时间序列的一种分析方法。利用R软件实现时间序列分析，程序简单灵活，输出结果清晰。通过R软件开展案例教学，一方面能进一步巩固和加深学生对理论知识的理解和认识，更好地促进学生对理论知识的学习；另一方面，有利于提高学生实践动手能力，培养创新精神。