时间序列ARIMA模型在R中的实现

J I A N G S U U N I V E R S I T Y江苏省餐饮业零售总额分析预测学校：江苏大学学院：财经学院班级：统计1201组员：韩亚琼3120812015马海燕3120812022顾君颖3120812020王培培3120812009陆金龙3120812029白卓3120812028完成时间：2014年12月13日星期六一、摘要二、引言三、数据分析原始数据的获取：本文所有的样本数据均来自《江苏统计年鉴——2014》（/2014nj/nj14.htm）得到的样本数据参见表1：表1 按行业分社会消费品零售总额这里我们仅用到第三列数据，为了方便分析，我们将餐饮业零售总额序列命名为caterts。

第一步序列的平稳性检验为判断一个序列是否平稳，我们主要通过时序图以及自相关图进行检验。

对caterts做时序图，有图形发现有明显的指数趋势，序列非平稳，也可以初步发现江苏省的餐饮业零售总额逐年递增，尤其是在新世纪以后，人们的生活水平逐年提高，对餐饮业的贡献也增大：图1 caterts序列时序图因为原序列有明显的指数趋势，故先对数列进行对数变换得到新的数列logcatets，序列图如下，具有明显的非线性增长趋势：图2 对数化后的时序图对具有明显线性趋势的数列常用的平稳化措施是差分，我们对logcaterts序列进行一阶差分得到新的数列difflogcaterts，时序图如下：图3 对数化和一阶差分后的时序图通过对时序图分析发现数列具有平稳性，为了方便分析，我们对difflogcaterts 序列进行中心化处理，得到新的数列x。

对x进行ADF检验（单位根检验）。

R语言中有专门的fUnitRoots包，里面有urdftest功能，是专门进行序列的ADF单位根检验，通过R软件得到如下结果：Title:Augmented Dickey-Fuller Unit Root TestTest Results:Test regression noneCall:lm(formula = z.diff ~ g.1 - 1 + g)Residuals:Min 1Q Median 3Q Max-0.224894 -0.051073 0.006261 0.043257 0.242110Coefficients:Estimate Std. Error t value Pr(>|t|)g.1 -0.59259 0.19455 -3.046 0.0047 **g 0.02909 0.18067 0.161 0.8731---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 0.08821 on 31 degrees of freedomMultiple R-squared: 0.2853, Adjusted R-squared: 0.2392F-statistic: 6.188 on 2 and 31 DF, p-value: 0.005481Value of test-statistic is: -3.046Critical values for test statistics:1pct 5pct 10pcttau1 -2.62 -1.95 -1.61从结果可以看出ADF统计量为-3.046，在1%，5%，10%的置信水平下均拒绝原假设，认为序列x平稳。

第二步模型的识别与定阶模型的识别与定阶主要是通过对序列的自相关和偏自相关观察得到的。

通过R 软件的到序列x的自相关和偏自相关图：图4 序列x的自相关图图5 序列x的偏自相关图通过自相关图和偏自相关图，可以发现序列x自相关明显拖尾，偏自相关有结尾的性质，在第十阶的时候，PACF超出了两倍的标准差，这可能是由于偶然因素引起的，可以考虑使用AR（1）和AR（11）模型，但是为了数据的完整性，尽量避免滞后期，可以考虑AR（1）模型。

在R语言中有专门的forecast包，包里的auto.arima函数可以对序列模型进行智能识别建模。

通过R语言的智能识别功能得到如下结果：Series: xARIMA(1,0,0) with zero meanCoefficients:ar10.4126s.e. 0.1512sigma^2 estimated as 0.006901: log likelihood=37.33AIC=-70.65 AICc=-70.28 BIC=-67.54通过如剑的输出结果我们发现结果与我们预测的一样，均是AR（1）模型。

得到的模型方程为：X t=0.1426X t−1+a t第三步模型的适应性检验模型的适应性检验实质上就是检验残差序列{a t}是否为白噪声序列。

其中最主要的是{a t}序列的独立性检验。

通过R软件得到模型的残差序列的时序图（偏）自相关图：图6 残差序列的时序、自相关和偏自相关图通过对图形的观察，初步判断残差序列{a t}为白噪声。

R语言中的Box.test函数也可对残差序列进行白噪声检验。

且得到的结果如下：LAG LB p[1,] 1 0.02128392 0.8840081[2,] 2 0.08484764 0.9584635[3,] 3 0.43750962 0.9323894[4,] 4 0.696146050.9518033[5,] 5 2.02913305 0.8450994[6,] 6 2.03204959 0.9167275[7,] 7 2.24882031 0.9447957[8,] 8 2.45473056 0.9638085[9,] 9 8.50017079 0.4846293[10,] 10 10.925501510.3633515[11,] 11 10.939876320.4483139[12,] 12 13.844479050.3107497可以发现无论滞后几期，P值均大于0.05，所以接受原假设，认为残差序列通过纯随机性检验。

故模型最终为ARIMA(1,0,0),即为AR（1）第四步模型参数的估计模型参数的估计主要有矩估计，最小二乘估计（ML估计），极大似然估，这里我们主要用的是最小二乘估计，估计的结果在第二步模型的识别中R语言智能建模中给出了，这里不做重复论述。

模型的估计方程如下：X t=0.1426X t−1+a t第五步模型的预测将模型的拟合值与模型的真实值进行画图对比，得到如下图：图7 拟合值—真实值发现拟合值与真实值之间有一定差别，这是由于残差序列的存在，将拟合值与残差序列相加得到的结果与X序列作图，结果如下：图8 加上残差的拟合值—真实值观察图8，可以发现拟合的特别好，所有的点都落在了真实值上。

我们用拟合的有效模型进行短期预测，比如我们预测2014年、2015年、2016、2017年和2018年的餐饮业零售总额。

先预测2014年、2015年、2016年、2017年和2018年的x，再预测餐饮业零售总额。

在R语言中可以使用predict函数对序列进行预测，这里我们进行向前5步预测，得到点预测，区间预测。

结果如下：$predTime Series:Start = 2014End = 2018Frequency = 1[1] -0.025******* -0.010******* -0.0043260335 -0.0017851186 -0.0007366213 $seTime Series:Start = 2014End = 2018Frequency = 1[1] 0.08307163 0.08986632 0.09097274 0.09115980 0.09119162利用公式X t±Zα×S.E.2可以得到预测值的区间估计，点估计和区间估计具体结果如下：通过还原零均值、差分和对数变换，就可以得到江苏省餐饮业零售总额（caterts）2014到2015年的预测值：利用公式：logcaterts2014=X2014+logcaterts2013+mean（difflogcaterts）caterts2014=exp⁡(logcaterts2014)得到江苏省2014年餐饮业的零售手总额预测值，以此类推，分别得到2015,2016,2017,2018年餐饮业零售总额的预测值。

值得一体的是时间序列预测效果越好，则向前预测的阶数就越少。

所以在该模型中预测最好的是2014年。

值得一体的是时间序列预测效果越好，则向前预测的阶数就越少。

所以在该模型中预测最好的是2014年。

四、参考文献[1] 王振龙.应用时间序列分析[M].北京：中国统计出版社，2010[2] Paul Teetor .R语言经典实例[M]. 北京：机械工业出版社，2013[3] Robert I.Kabacoff.R语言实战[M].北京：人民邮电出版社，2013五、附录R语言代码：```{r}catering<- read.table("E:\\catering.txt",header=F)caterts<- ts(catering,start=1978,frequency=1) #数据时间序列化plot.ts(caterts,main=("CATERING"),ylab="CATERING") #序列图logcaterts<- log(caterts)difflogcaterts<- diff(log(caterts)) #对数差分plot(difflogcaterts, main="logged and diffed")x <- scale(difflogcaterts,T,F) #零均值化acf(x,10) #自相关pacf(x,10) #偏自相关corr<- acf(x,10)cov<- acf(x,10,type = "covariance")library(fUnitRoots) #ADF单位根检验urdfTest(x)library(forecast) #建模fit<- auto.arima(x)fitresid<- fit$residual #生成残差序列residplotForecastErrors(fit$residuals)# 白噪声检验Result=0LAG=0LB=0p=0for(i in 1:12){Btest=Box.test(resid,type="Ljung-Box",lag=i)LAG[i]=iLB[i]=Btest$statisticp[i]=Btest$p.valueResult=cbind(LAG,LB,p)}Result #接受原假设，认为resid序列为白噪声#预测pre=predict(fit,n.ahead=5)preU=pre$pred + 1.96*pre$seL=pre$pred - 1.96*pre$sets.plot(x,pre$pred,col=1:2)lines(U,col="blue",lty="dashed")lines(L,col="blue",lty="dashed")plot(x,type="o")lines(fitted(fit),col="2",pch=7)points(fitted(fit),col="2",pch=7)plot(x,type="o")points(fitted(fit)+resid,col="blue",pch=20)# 预测2014~2015年x2014<- -0.025*******+logcaterts[36]+mean(difflogcaterts)x2015<- -0.010*******+x2014+mean(difflogcaterts)x2016<- -0.0043260335+x2015+mean(difflogcaterts)x2017<- -0.0017851186+x2016+mean(difflogcaterts)x2018<- -0.0007366213+x2017+mean(difflogcaterts)forecast1 <- exp(x2014)forecast2 <- exp(x2015)forecast3 <- exp(x2016)forecast4 <- exp(x2017)forecast5 <- exp(x2018)forecast1forecast2forecast3forecast4forecast5y <- c(forecast1,forecast2,forecast3,forecast3,forecast4,forecast5) y <- t(t(y))x_fit<- fitted(fit)x_fit<- as.matrix(x_fit,ncol=1)i<- 1z <- 0for (i in 1:35) {z[i] <- x_fit[i]+logcaterts[i]+mean(difflogcaterts)} x_fit_p<- as.matrix(z,ncol=1)s <- c(x2014,x2015,x2016,x2017,x2018)s <- t(t(s))x_fit_p<- rbind(x_fit_p,s)x_fit_pts<- ts(x_fit_p,start=1979,frequency=1)#真实值与预测值图plot.ts(caterts,main=("CATERING"),ylab="CATERING",ylim=c(0,4500),xlim=c(197 8,2018),type="l",col="black")points(exp(x_fit_pts),pch=8,col="red")```。