# 时间序列分析之ARIMA模型预测__R篇

from http://www.cnblogs.com/bicoffee/p/3838049.html

1. 处理数据

1.1. 导入forecast包

forecast包是一个封装的ARIMA统计软件包，在默认情况下，R没有预装forecast包，因此需要先安装该包

`> install.packages("forecast')`

```> library("zoo")
> library("forecast")```

1.2. 导入数据

```> airline <- read.table("airline.txt")
> airline```

V1 V2
1 1 112
2 2 118
3 3 132
4 4 129
5 5 121
6 6 135
7 7 148
8 8 148
9 9 136
10 10 119

(……..)

1.3. 将数据转化为时间序列格式(ts)

```> airline2 <- ariline[2]
> airts <- ts(airline2,start=1,frequency=12)```

2. 识别模型

2.1. 查看趋势图

`> plot.ts(airts)`

```> airlog <- log(airts)
> airdiff <- diff(airlog, differences=1)
> plot.ts(airdiff)```

2.2. 查看ACF和PACF

```> acf(airdff, lag.max=30)
> acf(airdff, lag.max=30,plot=FALSE)
Autocorrelations of series ‘airdiff’, by lag
0.0000 0.0833 0.1667 0.2500 0.3333 0.4167 0.5000 0.5833 0.6667 0.7500 0.8333
1.000 0.188 -0.127 -0.154 -0.326 -0.066 0.041 -0.098 -0.343 -0.109 -0.120
0.9167 1.0000 1.0833 1.1667 1.2500 1.3333 1.4167 1.5000 1.5833 1.6667 1.7500
0.199 0.833 0.198 -0.143 -0.110 -0.288 -0.046 0.036 -0.104 -0.313 -0.106
1.8333 1.9167 2.0000 2.0833 2.1667 2.2500 2.3333 2.4167 2.5000
-0.085 0.185 0.714 0.175 -0.126 -0.077 -0.214 -0.046 0.029```

```> pacf(airdff, lag.max=30)
> pacf(airdff, lag.max=30,plot=FALSE)
Partial autocorrelations of series ‘airdiff’, by lag
0.0833 0.1667 0.2500 0.3333 0.4167 0.5000 0.5833 0.6667 0.7500 0.8333 0.9167
0.188 -0.169 -0.101 -0.317 0.018 -0.072 -0.199 -0.509 -0.171 -0.553 -0.300
1.0000 1.0833 1.1667 1.2500 1.3333 1.4167 1.5000 1.5833 1.6667 1.7500 1.8333
0.551 0.010 -0.200 0.164 -0.052 -0.037 -0.108 0.094 0.005 -0.095 -0.001
1.9167 2.0000 2.0833 2.1667 2.2500 2.3333 2.4167 2.5000
0.057 -0.074 -0.048 0.024 0.073 0.047 0.010 0.033```

2.3. 利用auto.arima

```> auto.arima(airlog,trace=T)

ARIMA(2,1,2)(1,1,1)[12]                    : -354.4719
ARIMA(0,1,0)(0,1,0)[12]                    : -316.8213
ARIMA(1,1,0)(1,1,0)[12]                    : -356.4353
ARIMA(0,1,1)(0,1,1)[12]                    : -359.7679
ARIMA(0,1,1)(1,1,1)[12]                    : -354.9069
ARIMA(0,1,1)(0,1,0)[12]                    : -327.5759
ARIMA(0,1,1)(0,1,2)[12]                    : -357.6861
ARIMA(0,1,1)(1,1,2)[12]                    : -363.2418
ARIMA(1,1,1)(1,1,2)[12]                    : -359.6535
ARIMA(0,1,0)(1,1,2)[12]                    : -346.1537
ARIMA(0,1,2)(1,1,2)[12]                    : -361.1765
ARIMA(1,1,2)(1,1,2)[12]                    : 1e+20
ARIMA(0,1,1)(1,1,2)[12]                    : -363.2418
ARIMA(0,1,1)(2,1,2)[12]                    : -368.8244
ARIMA(0,1,1)(2,1,1)[12]                    : -368.1761
ARIMA(1,1,1)(2,1,2)[12]                    : -367.0903
ARIMA(0,1,0)(2,1,2)[12]                    : -363.7024
ARIMA(0,1,2)(2,1,2)[12]                    : -366.6877
ARIMA(1,1,2)(2,1,2)[12]                    : 1e+20
ARIMA(0,1,1)(2,1,2)[12]                    : -368.8244

Best model: ARIMA(0,1,1)(2,1,2)[12]

Series: airlog
ARIMA(0,1,1)(2,1,2)[12]

Coefficients:
ma1     sar1     sar2     sma1     sma2
-0.2710  -0.4764  -0.1066  -0.0098  -0.1987
s.e.   0.0995   0.1432   0.1087   0.1567   0.1130

sigma^2 estimated as 0.001188:  log likelihood=231.88
AIC=-369.57   AICc=-368.82   BIC=-352.9```

auto.arima提供的最佳模型为ARIMA(0,1,1)(2,1,2)[12]，我们可以同时测试两个模型，看看哪个更适合。

3. 参数估计

```> airarima1 <- arima(airlog,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12),method="ML")
> airarima1
Series: airlog
ARIMA(0,1,1)(0,1,1)[12]

Coefficients:
ma1     sma1
-0.3484  -0.5622
s.e.   0.0943   0.0774

sigma^2 estimated as 0.001313:  log likelihood=223.63
AIC=-441.26   AICc=-441.05   BIC=-432.92

> airarima2 <- arima(airlog,order=c(0,1,1),seasonal=list(order=c(2,1,2),period=12),method="ML")
> airarima2
Series: airlog
ARIMA(0,1,1)(2,1,2)[12]

Coefficients:
ma1    sar1     sar2     sma1    sma2
-0.3546  1.0614  -0.1211  -1.9130  0.9962
s.e.   0.0995  0.1094   0.1844   0.3887  0.3812

sigma^2 estimated as 0.0009811:  log likelihood=225.56
AIC=-439.12   AICc=-438.37   BIC=-422.44```

ARIMA(0,1,1)(0,1,1)[12] ：t(ma1)=-39.1791, t(sma1)=-93.8445

ARIMA(0,1,1)(2,1,2)[12] : t(ma1)=-35.8173,t(sar1)=88.68383,t(sar2)=-3.56141,t(sma1)=-12.6615,t(sma2)= 6.855526

4. 预测

预测五年后航空公司的销售额：

```> airforecast <- forecast.Arima(airarima1,h=5,level=c(99.5))
> airforecast
Point Forecast Lo 99.5 Hi 99.5
Jan 12 6.038649 5.936951 6.140348
Feb 12 5.988762 5.867380 6.110143
Mar 12 6.145428 6.007137 6.283719
Apr 12 6.118993 5.965646 6.272340
May 12 6.159657 5.992605 6.326709

> plot.forecast(airforecast)```