季节性ARIMA模型【R语言】

季节性的ARIMA模型可以预测含有季节性，趋势性的时间序列。他的形式如下

这里m是每一季节的周期值。季节项与非季节项的模型非常相近。但是季节项中包含了季节周期性。例如对于ARIMA(1,1,1)(1,1,1)4模型能够写成：

ACF与PACF
对于AR与MA模型的季节项，我们将会在ACF和PACF的lags上看到差异。例如，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。
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”)进行安装，再载入forecast包。
tsdisplay(data,xlab=“year”,ylab=“Retail index”) #显示自相关系数ACF,偏自相关系数PACF,与数据原始图像。

data_d1=diff(data,1)
从上图可以看出自相关系数呈现拖尾性，即序列Yt与之前的时刻有强烈的相关性，而且相关项的范围是比较大的，说明序列本身有趋势，为了消除这种趋势我们取每一项与前一项的差。
tsdisplay(data_d1,xlab=“year”,ylab=“Retail index”)

可以看出来，趋势基本消除了。序列较平稳。而PACF在季节性周期lags12处有spike，因此进一步做季节差分。
data_d1D1=diff(data_diff1,12)
tsdisplay(data_d1D1,xlab=“year”,ylab=“Retail index”)

所以对于非季节项，我们只做了一阶非季节差分，故d=1,从PACF看出p=2或3，从ACF看出q=1
对于季节项，我们也只做了一阶季节差分，故D=1，再看在lag12，24处看是否有spike，从PACF看出P=1，从ACF看出Q=1
从而选择模型一共有：Arima(2,1,1)(1,1,1)[12]和Arima(3,1,1)(1,1,1)[12]

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)) #选择AICc值小的，输出我们最终的预测结果～

还有一个更加傻瓜的做法就是使用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最小的。