R语言在时间序列分析中的应用（二）

 （复制之后排版太乱，并且没有图片，提供一个下载版

http://wenku.baidu.com/view/6efb18c8d5bbfd0a7956734c.html?st=1）

         written by MiltonDeng,

     from department of Statistics, XMU

六、ARMA模型模拟与识别

 

         ARMA模型基本理论不再赘述。基本的识别过程如下表：

 

 

         ARMA模型需要用到和arima相关的一族函数。这里暂时不涉及差分阶数I。

         下面将模拟几个特殊的ARMA模型。需要用到的函数是：

        

                   arima.sim(model, n)

-          model        ：生成的模型。具体看下面的例子。

-          n                  ：序列的长度。

 

注意这个函数不能生成非平稳序列。

 

（1）AR(1) :   X(t) = 0.8*X(t-1) + E(t)

D = arima.sim ( list( order = c(1, 0, 0), ar=0.8 ), n = 200 )

plot(D) ; 

par(mfrow=c(1,2)) ; acf(D) ; pacf(D)

         这里，order的三个参数分别对应AR,I,MA的阶数。

         如果用Box.test验证纯随机性的话，可以很容易拒绝它不是白噪声序列。

 

 

 

         可以明显看出acf的拖尾性，与pacf的一阶截尾性。

（2）AR(2) :   X(t) = 0.8*X(t-1) + 0.1X(t-1) + E(t)

D = arima.sim ( list( order = c(2, 0, 0), ar=c(0.8,0.1) ), n = 200 )

 

 

 

         可以明显看出acf的拖尾性，与pacf的二阶截尾性。

 

（2）MA(2) :   X(t) = E(t) – 0.7*E(t-1) – 0.2*E(0.2)   

##这里符号我还没验证是指正还是指负。

         D = arima.sim ( list( order = c(0, 0, 2), ma=c(0.7,0.2) ), n = 200 )

 

 

可以明显看出pacf的拖尾性，与acf的二阶截尾性。

        

（2）MA(2) :   X(t) = E(t) – 0.7*E(t-1) – 0.2*E(0.2)   

D = arima.sim ( list( order = c(2, 0, 2), ar=c(0.3,0.4),ma=c(0.5,0.4) ), n = 200 )

 

 

         都是拖尾的。但在实际情况中，也经常遇到难以判断的情况。

 

         下面我们尝试生成几组非平稳序列。

 

 

七、ARMA模型的估计

 

对于一个给定的序列，如何用R求出相关的参数呢？

                   arima( x, order=c(0,0,0), method=c(“CSS-ML”,”ML”,”CSS”),include.mean)

-          x                            ：待估计序列。

-          order                             ：指定阶数。

-          method               ：估计方法。ML为极大似然估计，CSS为条件最小二乘估计。

-          include.mean    ：是否包含均值项(intercept)。

 

另一个可用的函数是auto.arima()，这个函数在forecast包里，它的估计依据是AIC/BIC标准。尤其对ARMA模型定阶困难时，可以考虑用auto.ARIMA。

现在模拟生成一个ARMA(2,2)的时间序列，然后分别运用两个函数进行估计。

         D = arima.sim ( list( order = c(2, 0, 2), ar=c(-0.5,0.4),ma=c(0.5,0.4) ), n = 1000 )

         ARIMA=arima(D,order=c(2,0,2),method="ML"); summary(ARIMA)

library(forecast);  ARIMA.auto=auto.arima(D);  summary(ARIMA.auto)

 

arima()的估计结果：

 

Coefficients:

          ar1     ar2     ma1     ma2  intercept

      -0.5038  0.3334  0.5354  0.4566    -0.0358

s.e.   0.0480  0.0476  0.0458  0.0402     0.0566

 

arima.auto()的估计结果：

 

Coefficients:

         ar1     ar2     ma1     ma2

      -0.503  0.3342  0.5350  0.4564

s.e.   0.048  0.0476  0.0458  0.0403

 

可以看到两者的结果是一样的，也基本接近模拟时设定的参数。

但是还是存在一些困难难以克服。在我们难以定阶时，尤其是ARMA(p,q)模型中，难以确定arima() 函数中的阶数，因为我们并不知道p,q真正是多少。而在试验过程中，我测试了多组时序数据，在样本量200时auto.arima()确定的阶数也常常和模拟的参数不一样，这就是说auto.arima的结果（即AIC准则结果）经常不准确。

不过，在评价模型结果时，可以采取这样的标准：只要残差序列为纯随机序列，就可以认为模型已经建立成功。可以这样获取残差序列：

resid=ARIMA$residual

         再重复前面的随机性检验即可。不过严格说，对残差序列的随机性检验要扣除一些自由度，要用到前面Box.test的dffit，令其等于p+q。

 

八、ARMA模型的预测

 

         建立好ARMA模型后，我们可以进行预测。

 

                  predict(model,n.ahead)

-          model        ：这里需要填入估计好的ARIMA模型。

-          n.ahead     ：需要预测几期。

 

predict() 将会返回两个时间序列数据。一个是$pred，即点预测，另一个是$se即对应的估计误差。

下面给出一个较为完整的例子：

 

D = arima.sim ( list( order = c(2, 0, 2), ar=c(0.5,-0.2),ma=c(0.5,0.2) ), n = 100 )

plot(D) ; 

par(mfrow=c(1,2)) ; acf(D) ; pacf(D)

 

ARIMA=arima(D,order=c(2,0,2),method="ML",include.mean=FALSE)

summary(ARIMA)

 

library(forecast)

ARIMA.auto=auto.arima(D)

summary(ARIMA.auto)

 

resid=ARIMA $residual

Box.test(resid,lag=6)

Box.test(resid,lag=12)

 

pre=predict(ARIMA,n.ahead=10)

U=pre$pred + 1.96*pre$se

L=pre$pred - 1.96*pre$se

ts.plot(D,pre$pred,col=1:2)

lines(U,col="blue",lty="dashed")

lines(L,col="blue",lty="dashed")

        

 

         最终可以得到一个这样的图形，最后的红蓝线即预测结果。

         另外提供一个forecast.Arima() 函数，具体用法不再详述，不过使用更简便一些。

pre2=forecast.Arima(ARIMA,h=10,level=c(0.95))

pre2

plot.forecast(pre2)

 

        

         两者的结果是完全一样的。