R语言与时间序列学习笔记（1）

       今天分享的是R语言中时间序列的有关内容。主要有：时间序列的创建，ARMA模型的建立与自相关和偏自相关函数。

一、          时间序列的创建

时间序列的创建函数为：ts().函数的参数列表如下：

ts(data = NA, start = 1, end = numeric(),frequency = 1,
  deltat = 1, ts.eps = getOption("ts.eps"), class = , names = )
参数说明：data：这个必须是一个矩阵，或者向量，再或者数据框frame

          Frequency：这个是时间观测频率数，也就是每个时间单位的数据数目

          Start：时间序列开始值，允许第一个个时间单位出现数据缺失

举例：ts(matrix(c(NA,NA,NA,1:31,NA),byrow=T,5,7),frequency=7,names=c("Sun","Mon ","Tue", "Wen" ,"Thu","Fri"," Sat"))

运行上面的代码就可以得到一个日历：

Sun  Mon Tue Wen Thu  Fri  Sat

  NA   NA  NA   1  2    3    4

   5    6   7   8  9   10   11

  12   13  14  15 16   17   18

  19   20  21  22 23   24   25

  26    27 28  29  30   31   NA

在R语言中本身也有不少数据集，比如统计包中的sunspots，你可以通过函数data(sunspots)来调用它们。

二、       一些时间序列模型

这里主要介绍AR，MA，随机游走，余弦曲线趋势，季节趋势等

首先介绍一下AR模型：AR模型，即自回归（AutoRegressive,AR）模型，数学表达式为：　　AR :y(t)=a1y(t-1)+...any(t-n)+e(t)

其中，e(t)为均值为0，方差为某值的白噪声信号。

那么产生AR模型的数据，我们就有两种方法：1、调用R中的函数filter（线性滤波器）去产生AR模型；2、根据AR模型的定义自己编写函数

先说第一种方法：调用R中的函数filter（线性滤波器）去产生AR模型

介绍函数filter的用法如下：

filter(x, filter, method = c("convolution", "recursive"),

       sides = 2, circular = FALSE, init)

对于AR（2）模型x（t）=x(t-1)--0.9x(t-2)+e(t)

w<-rnorm(550)#我们假定白噪声的分布是正态的。

x<-filter(w,filter=c(1,-0.9),"recursive")

#方法：无论是“卷积”或“递归”（可以缩写）。如果使用移动平均选择“卷积”：如果“递归”便是选择了自回归。

再说第二种方法：依据定义自己编程产生AR模型，还是以AR（2）模型x（t）=x(t-1)--0.9x(t-2) +e(t)为例，可编写函数如下：

w<-rnorm(550)

AR<-function(w){

x<-w

x[2]=x[1]+w[1]

for(i in 3:550)

x[i]=x[i-1]-0.9*x[i-2]+w[i]

x

}

调用AR（W）即可得到。如果对相同的随机数，我们可以发现两个产生的时间序列是一致的。当然对于第二种方法产生的序列需要转换为时间序列格式，用as.ts()处理。

类似的，我们给出MA，随机游走的模拟：

MA模型：

w<-rnorm(500)

v<-filter(w,sides=2,rep(1,3)/3)

随机游走：

w<-rnorm(200)

x<-cumsum(w)#累计求和，seeexample：cumsum(1:!0)

wd<-w+0.2

xd<-cumsum(wd)

可以做出相应的图形：

 




再说一下季节性模型：

最简单的季节模型就是一个分段的周期函数。比如说某地区一年的气温就是一个季节性模型。利用TSA包里给出的数据tempdub我们可以发现他就是这样的模型

给出验证：

library(TSA)

data(tempdub)

month<-season(tempdub)

model1<-lm(tempdub~month)

summary(model1)

根据R输出的结果：

Call:

lm(formula = tempdub ~month)

 

Residuals:

    Min     1Q  Median      3Q    Max

-8.2750 -2.2479  0.1125 1.8896  9.8250

 

Coefficients:

               Estimate Std. Error t valuePr(>|t|)   

(Intercept)      16.608      0.987 16.828  < 2e-16 ***

monthFebruary     4.042     1.396   2.896  0.00443 **

monthMarch       15.867     1.396  11.368  < 2e-16 ***

monthApril       29.917      1.396 21.434  < 2e-16 ***

monthMay         41.483      1.396 29.721  < 2e-16 ***

monthJune        50.892      1.396 36.461  < 2e-16 ***

monthJuly        55.108      1.396 39.482  < 2e-16 ***

monthAugust      52.725      1.396 37.775  < 2e-16 ***

monthSeptember   44.417     1.396  31.822  < 2e-16 ***

monthOctober     34.367     1.396  24.622  < 2e-16 ***

monthNovember    20.042     1.396  14.359  < 2e-16 ***

monthDecember     7.033     1.396   5.039 1.51e-06 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’1

 

Residual standarderror: 3.419 on 132 degrees of freedom

Multiple R-squared:0.9712,     Adjusted R-squared: 0.9688

F-statistic: 405.1 on11 and 132 DF,  p-value: < 2.2e-16

 

这里2月份系数表明了一月份平均气温与二月份平均气温的差异，以此类推。

在介绍一下一个季节模型：余弦趋势μ1=βcos（2pi*f*t+φ）

还是考虑上面气温的例子：

验证：

har<-harmonic(tempdub,1)

model2<-lm(tempdub~har)

summary(model2)

看看结果：

Call:

lm(formula = tempdub ~har)

 

Residuals:

     Min      1Q   Median       3Q     Max

-11.1580  -2.2756 -0.1457   2.3754  11.2671

 

Coefficients:

               Estimate Std. Error t valuePr(>|t|)   

(Intercept)     46.2660    0.3088 149.816  < 2e-16 ***

harcos(2*pi*t)-26.7079     0.4367 -61.154  < 2e-16 ***

harsin(2*pi*t)  -2.1697    0.4367  -4.968 1.93e-06 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’1

 

Residual standarderror: 3.706 on 141 degrees of freedom

Multiple R-squared:0.9639,     Adjusted R-squared: 0.9634

F-statistic:  1882 on 2 and 141 DF,  p-value: < 2.2e-16

我们可以作图来看拟合效果：


 

 

顺便指出季节模型也可以模拟：比如μ1=βcos（2pi*f*t+φ）模型可以模拟如下：

t<-1:500

w<-rnorm(500)

c<-2*cos(2*pi*t/50+0.6*pi+w)

 

三、       自相关与偏自相关

给出自相关的定义：在信息分析中，通常将自相关函数称之为自协方差方程。 用来描述信息在不同时间的，信息函数值的相关性。详情可参见wiki：http://zh.wikipedia.org/wiki/%E8%87%AA%E7%9B%B8%E5%85%B3

我们可以根据定义给出自相关系数（ACF）的算法：

例如数据：
> x<-1:10
>u<-mean(x)
>v<-var(x)
>sum((x[1:9]-u)*(x[2:10]-u))/(9*v) #延迟1
[1] 0.7
>sum((x[1:8]-u)*(x[3:10]-u))/(9*v)  #延迟2
[1]0.4121212
>sum((x[1:7]-u)*(x[4:10]-u))/(9*v)  #延迟3
[1]0.1484848
在R中也提供了直接计算acf的函数acf（），利用该函数也计算1至3阶的acf，结果如下：
>a<-acf(x,3)
> a
 
Autocorrelationsof series ‘x’, by lag
 
    0    1     2     3
1.0000.700 0.412 0.148
可以看出，是一样的。
利用acf（）可以处理很多阶的acf，以太阳黑子数的数据集做例子：
>data(sunspots)
>acf(sunspots)  #给出了相应的图形
>a<-acf(sunspots,6)  #为下面做估计做铺垫，列出前6阶的acf
> a
 
Autocorrelationsof series ‘sunspots’, by lag
 
0.00000.0833 0.1667 0.2500 0.3333 0.4167 0.5000
 1.000 0.922  0.890  0.875 0.864  0.850  0.836
偏自相关：
对于一个平稳AR(p)模型，求出滞后k自相关系数p(k)时，实际上得到并不是x(t)与x(t-k)之间单纯的相关关系。因为x(t)同时还会受到中间k-1个随机变量x(t-1)、x(t-2)、……、x(t-k+1)的影响，而这k-1个随机变量又都和x(t-k)具有相关关系，所以自相关系数p(k)里实际掺杂了其他变量对x(t)与x(t-k)的影响。 
为了能单纯测度x(t-k)对x(t)的影响，引进偏自相关系数的概念。 
对于平稳时间序列{x(t)}，用数学语言描述就是： 
　　p[(x(t),x(t-k)]|(x(t-1),……，x(t-k+1)={E[(x(t)-Ex(t)][x(t-k)-Ex(t-k)]}/E{[x(t-k)-Ex(t-k)]^2} 
这就是滞后k偏自相关系数的定义。
总之，偏自相关就是在试图解释在剔除了中间k-1个随机变量x(t-1)、x(t-2)、……、x(t-k+1)的干扰之后，x(t-k)对x(t)影响的相关程度。
在R语言中，使用函数PACF（）可求解
还是使用太阳黑子数的例子：
> b<-pacf(sunspots,6)
> b
 
Partial autocorrelations of series ‘sunspots’, bylag
 
0.0833 0.1667 0.2500 0.3333 0.4167 0.5000
 0.922  0.272 0.189  0.135  0.064 0.044
最后，我们利用这两个函数来看看AR（p）,MA(q)的自相关函数与偏自相关函数的截尾性与拖尾性。
利用二中所介绍的方法生成AR（2），MA（2）的数据。
AR（2）模型：

w<-rnorm(550)#我们假定白噪声的分布是正态的。

x<-filter(w,filter=c(1,-0.9),"recursive")

MA(3)模型：

w<-rnorm(500)

v<-filter(w,sides=2,rep(1,3)/3)

> qq<-pacf(x,5)
> qq
 
Partial autocorrelations of series ‘x’, by lag
1   2    3   4    5    
 0.532-0.861 -0.082  0.000
可以看出AR（2）模型的偏自相关函数是截尾的（但由于这个是数据，所以出现pacf只能看出趋势，而不是在2步后直接变为0）
对于MA（3）模型的自相关函数，由于v的第一项与最后一项缺失，不妨截取v的一部分数据，命名为a,有：
> y<-acf(a,5)
> y
 
Autocorrelations of series ‘a’, by lag
 
    0     1    2     3       4     5
1.000   0.652   0.397   0.059  0.067  0.035
也可以看出趋势。
关于给出模型后的参数估计，我们将在下一篇博文中讨论。