时间序列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》（https://www.360docs.net/doc/6c6207039.html,/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 Test

Test Results:

Test regression none

Call:

lm(formula = z.diff ~ https://www.360docs.net/doc/6c6207039.html,g.1 - 1 + https://www.360docs.net/doc/6c6207039.html,g)

Residuals:

Min 1Q Median 3Q Max

-0.224894 -0.051073 0.006261 0.043257 0.242110

Coefficients:

Estimate Std. Error t value Pr(>|t|)

https://www.360docs.net/doc/6c6207039.html,g.1 -0.59259 0.19455 -3.046 0.0047 **

https://www.360docs.net/doc/6c6207039.html,g 0.02909 0.18067 0.161 0.8731

---

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.08821 on 31 degrees of freedom

Multiple R-squared: 0.2853, Adjusted R-squared: 0.2392

F-statistic: 6.188 on 2 and 31 DF, p-value: 0.005481

Value of test-statistic is: -3.046

Critical values for test statistics:

1pct 5pct 10pct

tau1 -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: x

ARIMA(1,0,0) with zero mean

Coefficients:

ar1

0.4126

s.e. 0.1512

sigma^2 estimated as 0.006901: log likelihood=37.33

AIC=-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步预测，得到点预测，区间预测。结果如下：

$pred

Time Series:

Start = 2014

End = 2018

Frequency = 1

[1] -0.025******* -0.010******* -0.0043260335 -0.0017851186 -0.0007366213 $se

Time Series:

Start = 2014

End = 2018

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

fit

resid<- fit$residual #生成残差序列

resid

plotForecastErrors(fit$residuals)

# 白噪声检验

Result=0

LAG=0

LB=0

p=0

for(i in 1:12){

Btest=Box.test(resid,type="Ljung-Box",lag=i)

LAG[i]=i

LB[i]=Btest$statistic

p[i]=Btest$p.value

Result=cbind(LAG,LB,p)

}

Result #接受原假设，认为resid序列为白噪声

#预测

pre=predict(fit,n.ahead=5)

pre

U=pre$pred + 1.96*pre$se

L=pre$pred - 1.96*pre$se

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

forecast1

forecast2

forecast3

forecast4

forecast5

y <- c(forecast1,forecast2,forecast3,forecast3,forecast4,forecast5) y <- t(t(y))

x_fit<- fitted(fit)

x_fit<- as.matrix(x_fit,ncol=1)

i<- 1

z <- 0

for (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")

```