1、ARIMI模型介绍
在解决周期性数据的预测问题时，一个较为基本的模型是ARIMA模型。下面我们就来简单的介绍一下这个模型，以及介绍一下使用R语言的ARIMA包求解问题的基本步骤。
ARIMA全称为自回归积分移动平均模型(Autoregressive Integrated Moving Average Model)，简记ARIMA，是由博克思(Box)和詹金斯(Jenkins)于70年代初提出一著名时间序列预测方法，所以又称为box-jenkins模型、博克思-詹金斯法。其中ARIMA（p，d，q）称为差分自回归移动平均模型，AR是自回归， p为自回归项； MA为移动平均，q为移动平均项数，d为时间序列成为平稳时所做的差分次数。ARIMA模型根据原序列是否平稳以及回归中所含部分的不同，包括移动平均过程（MA）、自回归过程（AR）、自回归移动平均过程（ARMA）以及ARIMA过程。季节性的ARIMA模型则加入季节性这个周期性。可以预测含有季节性，趋势性的时间序列。他的形式如下
这里m等于每一季节的周期值。季节项与非季节项的模型非常相近。但是季节项中包含了季节周期性。例如对于ARIMA(1,1,1)(1,1,1)4模型能够写成：
ACF与PACF
ACF是自相关函数（Auto Correlation Function）的简称，它是描述一个时间序列t时刻的Yt值与其周围各个数据(Yt-1,Yt-2…Yt-n)之间线性组合的相关性的；PACF则是偏自相关函数（Partial Auto Correlation Function）的简称，它是描述的是t时刻的Yt值与远端的Yt-k之间的相关性，而不考虑其他的线性组合。对于AR与MA模型的季节项，我们将会在ACF和PACF的lags上看到差异(横坐标lags就是k，也就是相距了多少个时间单位)。例如，ARIMA(0,0,0)(0,0,1)12模型，我们将会在ACF的lag12处看到一个spike(突出的钉子，即非零显著的形象描述)，而在其他地方看不到spike。PACF将会在周期的位置出现指数衰减，即lags 12，24，36…类似的，一个ARIMA(0,0,0)(1,0,0)12模型则会显示出在ACF图的周期性位置显示出指数衰减，而在PACF图中的lag12处看到一个spike。我以中国人民大学统计学院易丹辉教授在时间序列课上给出的数据book19为例进行分析，这里用到的软件是R，课上易丹辉教授用的是EViews6.0。
2、R语言实现:
data=read.csv("~/Downloads/Book19.csv") #读取本地csv文件
data=ts(data[,2],frequency=12,start=c(1990,1)) #[,2]是只读取第2列数据，start是在数据上增加起始时间，frequency是增加周期，ts为转换成ts格式。
library(forecast) #载入程序包forecast，如果提示没有安装则用install.package("forecast")进行安装再载入。
tsdisplay(data,xlab="year",ylab="Retail index") #显示自相关系数ACF,偏自相关系数PACF,图像。
从上图可以看出自相关系数呈现拖尾性，即序列Yt与之前的时刻有强烈的相关性，而且相关项的范围是比较大的，说明序列本身有趋势，为了消除这种趋势我们取每一项与前一项的差。
data_d1=diff(data,1) # 取每一项与前一项的差。
tsdisplay(data_d1,xlab="year",ylab="Retail index")
可以从图中看到趋势基本消除了。序列基本平稳。
data_d1D1=diff(data_diff1,12)
tsdisplay(data_d1D1,xlab="year",ylab="Retail index")
所以对于非季节项，d=1,从PACF看出p=2或3，从PAC看出q=1
对于季节项，D=1，在lag12，24处看是否有spike，从PACF看出P=1或2，从ACF看出Q=1
从而选择模型一共有：
fit1 <- Arima(data, order=c(2,1,1),seasonal=c(1,1,1))
fit1 #显示AICc=1110.09
tsdisplay(residuals(fit1)) #显示残差
fit2 <- Arima(data, order=c(3,1,1),seasonal=c(1,1,1))
fit2 #显示 AICc=1112.3
tsdisplay(residuals(fit1)) #显示残差
plot(forecast(fit1,h=12)) #输出我们最终的预测结果～
5.jpg
还有一个更加傻瓜的做法就是使用auto.arima函数自动定p,d,q,P,D,Q等参数，只需要一行代码。
arima1<-auto.arima(data,trace=T)
如下图：
with drift代表有趋势，最终模型可以加上d=1去除趋势。
最终的模型也可以是：
Arima(2,1,0)(1,1,0)[12]
其AICc=1107.86
plot(forecast(arima(data,order=c(2,1,0),seasonal=c(1,1,0)),h=12))
这个结果是各个模型中AICc最小的。
