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

 

     written by MiltonDeng,

     from department of Statistics, XMU

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

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

 

一、数据的读入与时间序列化

 

         时间序列数据是R语言中一种特定形式的数据类型。R语言中有许多专门针对时间序列数据编制的函数。但在运用这些函数前，首先需要对序列进行时间序列化，即运用将一组数据转化为时间序列数据这种数据形式。

 

         ts(data = NA, start = 1, end = numeric(), frequency = 1 )

-          data           ：     要进行时间序列化的向量。如果是矩阵，会按列分别处理。

-          start           ：     起始时点。

-          end             ：     结束时点。

-          frequency ：     频度。取12时自动识别为月，4时自动识别为季度。

 

         首先，应该注意，与多元统计等统计分析不同，单变量时间序列根本上是一个向量，而不是一个矩阵，或者一个表。因此，在读入时间序列数据时，常常会有一些不必要的麻烦。由于数据整理时，数据编排往往不一致，所以应该注意以何种方式读入时间序列数据。

         以txt文件为例，假设这组时间序列为“5,4,2,3,6,1”，在文件中以如下方式出现（人大出版社，王燕老师的《应用时间序列分析》中的数据均以这种形式出现）：

                            5       4       2

                            3       6       1

         如果以     read.table()      读入，读入的是一个2*3的矩阵，如果直接以 ts() 时间序列化，R会认为这是3列数据，将对它们分别生成3组时间序列数据。下面提供两种方法：

         （1）以scan()直接读入为向量。

         （2）以read.table()读入这个矩阵为D，然后

                            D=t(D)                         #因为是以行录入，所以要先转置

                            D=as.vector(D)                  #将矩阵线性化为向量D

 

二、描述时间序列

        

         对于单组时间序列：

         （1）ts.plot(D)        #可以直接对未ts的向量绘图

         （2）ts(D); plot(D)

         如果D是一个矩阵或frame，可以对多个序列画在同一张图里。

 

三、自相关与偏自相关

 

         自相关系数可以通过acf() 实现：

         acf( x, lag.max=NULL, type=c("correlation", "covariance", "partial"), plot=TRUE )

-          x                  ：时间序列数据列。

-          lag.max     ：所计算的最大滞后阶数。

-          type            ：分别为“自相关系数”、“自协方差”、“偏自相关系数”

-          plot             ：是否要绘出acf的图形。

偏自相关可以通过 pacf() 实现，原理相似。也可以acf() 中type=”partial”。

         不过注意 acf 是从0开始的，pacf 是从1阶开始的。

四、平稳性检验

 

         （1）时序图检验

         平稳时间序列均值、方差为常数，因此可以直接用 plot() 检查时序图，平稳时间序列应始终在一个常数值附近随机波动，而且波动范围有界。

         （2）自相关图检验

         平稳序列通常有短期相关性。因此自相关系数会很快地衰减向零。

         （3）单位根检验

         后面再讲。

 

         下面我们观察几种特定时间序列的时序图与自相关图。

         （1）确定趋势时间序列：

 

         首先我们尝试生成一个带有随机扰动的线性趋势序列：

X=seq(0.1,5,0.02)

Y=2+0.8*X

D=jitter(Y,amount=0.2)

D=ts(D)

plot(D)

         这个序列是D=2+0.8*t掺入随机扰动后形成的。

        

         对一个单增序列来说，acf是非常缓慢下降的序列，而pacf=0。

 

（2）周期性时间序列：

X=seq(0.1,15,0.1)

Y=sin(X)

D=jitter(Y,amount=0.2)

 

 

         acf有伪周期的特点。

（3）伪周期时间序列：

X=seq(0.1,20,0.1)

Y=100*sin(X)/X

D=jitter(Y,amount=3)

        

 

 

五、纯随机检验

 

         也称白噪声序列。当时间序列满足（1）随时间变化，均值为常数（2）不同时点互不相关时，序列为纯随机的。它有两个性质（1）纯随机性：任意两时点协方差为0；（2）方差齐性：任意时点方差为定值。因此，从图形上，纯随机序列的 acf 应该均为0，而且时序图中波动范围是恒定的。

         假设检验的检验定理为Barlett 定理，统计量有 Q 统计量和 LB 统计量，两者均渐进服从卡方分布。原假设“该序列为白噪声序列”，即拒绝时认为序列为非白噪声序列，接受原假设时为不能显著拒绝纯随机假定。

        

         Box.test( x, lag=1, type = c(“Box-Pierce”, “Ljung-Box”), fitdf=0 )

-          x         ：待检验序列。

-          lag     ：滞后阶数。

-          type  ：“Box-Pierce”即Q统计量，“Ljung-Box”即LB统计量。

-          fitdf   ：自由度调整。

 

让我们模拟并观察纯随机序列的特性。可以简单生成一组服从正态分布的序列，作为纯随机序列。然后绘制时序图、acf、pacf图。

 

D=rnorm(300,0)

D=ts(D)

par(mfrow=c(3,1))

plot(D)

acf(D)

pacf(D)

 

         可以看到，纯随机序列如前所述的几个特点：

         （1）围绕一个固定水平波动。即均值恒定的特点。

         （2）波动范围基本一致。即方差恒定的特点。

         （3）acf/pacf均基本为0。

 

下面进行纯随机性检验。Box.test() 只能针对一个特定阶数进行检验，我们可以通过下面的一段小程序对各阶均进行检验，并整合到一个表格里。

 

Result=0;LAG=0;LB=0;p=0

 

for(i in 1:12){

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

LAG[i]=i

LB[i]=Btest$statistic

p[i]=Btest$p.value

Result=cbind(LAG,LB,p)

}

 

Result

 

 

      LAG      LB          p

 [1,]   1        0.1813234       0.6702384

 [2,]   2        0.5826921       0.7472571

 [3,]   3        0.5863632       0.8995470

 [4,]   4        0.8876635       0.9263240

 [5,]   5        3.2322386       0.6642308

 [6,]   6        4.9142115       0.5548624

 [7,]   7        6.2369390       0.5123718

 [8,]   8        6.2749378       0.6164647

 [9,]   9        6.4707625       0.6920331

[10,]  10        8.2027704       0.6090380

[11,]  11        8.4882660       0.6690124

[12,]  12        8.8843582       0.7127748

 

         可见，各阶的Bj检验均不能拒绝原假设，因此可以认为该序列是纯随机的。