R语言与时间序列学习笔记.docx

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-wx2=x1+w1for(i in 3:550)xi=xi-1-0.9*xi-2+wix调用AR（W）即可得到。如果对相同的随机数，我们可以发现两个产生的时间序列是一致的。当然对于第二种方法产生的序列需要转换为时间序列格式，用as.ts()处理。类似的，我们给出MA，随机游走的模拟：MA模型：w-rnorm(500)v-filter(w,sides=2,rep(1,3)/3)随机游走：w-rnorm(200)x-cumsum(w)#累计求和，see example：cumsum(1:!0)wd-w+0.2xd-cumsum(wd)可以做出相应的图形：再说一下季节性模型：最简单的季节模型就是一个分段的周期函数。比如说某地区一年的气温就是一个季节性模型。利用TSA包里给出的数据tempdub我们可以发现他就是这样的模型给出验证：library(TSA)data(tempdub)month-season(tempdub)model1|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 standard error: 3.419 on 132 degrees of freedomMultiple R-squared: 0.9712, Adjusted R-squared: 0.9688 F-statistic: 405.1 on 11 and 132 DF, p-value: 2.2e-16 这里2月份系数表明了一月份平均气温与二月份平均气温的差异，以此类推。在介绍一下一个季节模型：余弦趋势1=cos（2pi*f*t+）还是考虑上面气温的例子：验证：har-harmonic(tempdub,1)model2|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 standard error: 3.706 on 141 degrees of freedomMultiple 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:500w-rnorm(500)c x u v sum(x1:9-u)*(x2:10-u)/(9*v) #延迟11 0.7 sum(x1:8-u)*(x3:10-u)/(9*v) #延迟21 0.4121212 sum(x1:7-u)*(x4:10-u)/(9*v) #延迟31 0.1484848在R中也提供了直接计算acf的函数acf（），利用该函数也计算1至3阶的acf，结果如下： a aAutocorrelations of series x, by lag 0 1 2 3 1.000 0.700 0.412 0.148 可以看出，是一样的。利用acf（）可以处理很多阶的acf，以太阳黑子数的数据集做例子： data(sunspots) acf(sunspots) #给出了相应的图形 a aAutocorrelations of series sunspots, by lag0.0000 0.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)/Ex(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 bPartial autocorrelations of series sunspots, by lag0.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 qq qqPartial autocorrelations of series x, by lag1 2 3 4 5 0.532 -0.861 -